[Insight-users] RecursiveGaussianImageFilter

[Silvano Agliozzo]

Hi Luis and Mathieu,

I'm working with double values. The min of the 3D images is around -3000 
housfield, whereas the max varies between 500-1000.

On Mon, 9 Feb 2004, Mathieu Malaterre wrote:

> Luis Ibanez wrote:
> > 
> > Hi Silvano,
> > 
> > You are right, the results of the derivative shouldn't
> > change when you add a constant value to the input image.
> 
> Except if you are overflow. Are you working on char / unsigned char 
> image ? Then you have great chance to be overflow even if you only add a 
> small constant. BTW if you are unsigned don't add any negative value or 
> check first the scalar range (=min/max) of your image.
> 
> 
> Mathieu
> 
> 

I'm going to joint hereafter a code that generate a spherical intensity 
field with the following equation x^2 + y^2 + z^2. So the derivatives 
should be egual to :
pdx = 2*x, pdy = 2*y, pdx = 2*z, 
pdxy = pdxz = pdyz = 0,
pdxx = pdyy = pdzz = 2.

the values are almost close to the expected ones execpt for the pdxx, 
pdyy, pdzz, but with the right choice of the sigma of the different 
filters, one can obtain better results. Nevertheless if you add a constant, 
like reported in the code you'll get different results, at least this is what 
I obtained. I have to leave now, but I'm looking forward to read your 
answer tomorrow morning. Thank you,

#include <iostream>
#include <vector>

#include "itkRecursiveGaussianImageFilter.h"
#include "itkImportImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageDuplicator.h"
#include "itkImage.h"
#include <string.h>
    
struct PartialDerivative
{
    double x, y, z, xy, xz, yz, xx, yy, zz;
};
    
int main()
{   
    int i, j, k;
    unsigned short dimx = 100;
    unsigned short dimy = 100;
    unsigned short dimz = 100;
    unsigned long numberOfElements = dimx * dimy *dimz;

    int half = 49;

    double *volume = new double[numberOfElements];

    typedef double pixelType;

    typedef itk::Image<pixelType, 3> ImageType;

    typedef itk::ImageRegionConstIterator< ImageType > 
RegionConstIteratorType;

    typedef itk::ImageDuplicator< ImageType > ImageDuplicatorType;
        
    ImageType::Pointer vol = ImageType::New();

    ImageType::IndexType begin;
    
    begin[0] = 0;
    begin[1] = 0;
    begin[2] = 0;

    ImageType::SizeType size;
    
    size[0] = dimx;
    size[1] = dimy;
    size[2] = dimz;

    ImageType::RegionType region;

    region.SetIndex(begin);
    region.SetSize(size);
    double spacing[ ImageType::ImageDimension ];
    spacing[0] = 1.0; // along X direction
    spacing[1] = 1.0; // along Y direction
    spacing[2] = 1.0; // along Z direction

    double origin[ ImageType::ImageDimension ];
    origin[0] = 0.0; // X coordinate
    origin[1] = 0.0; // Y coordinate
    origin[2] = 0.0; // Z coordinate

    vol->SetRegions(region);
    vol->Allocate();
    vol->SetSpacing( spacing );
    vol->SetOrigin( origin );

    ImageType::IndexType pixelIndex;
    ImageType::PixelType value;

    std::cout<<"Getting Volume ...";

    //here you can set the value of the constant 
    //to be added to the values of the 3D image
    ImageType::PixelType constant = 0;// +/- 1000 or even much more 

    for(k = 0; k < dimz ; k++)
        for(j = 0; j < dimy; j++)
            for(i = 0; i < dimx; i++)
            {
                pixelIndex[0] = i;
                pixelIndex[1] = j;
                pixelIndex[2] = k;

                // getting a spherical intensity field
                value = (static_cast<double>(i-half)*static_cast<double>(i-half)+static_cast<double>(j-half)*static_cast<double>(j-half)+static_cast<double>(k-half)*static_cast<double>(k-half));
                value += constant;
                vol->SetPixel(pixelIndex, value);
            }

    std::cout<<"done"<<std::endl;

    std::cout<<"Getting a surface ...";

    std::vector<ImageType::IndexType> s ;

    RegionConstIteratorType volIt( vol, vol->GetRequestedRegion() );

    for(volIt.GoToBegin(); !volIt.IsAtEnd(); ++volIt)
    {
            if(volIt.Value() == 324+constant)//324 = 18^2
                s.push_back(volIt.GetIndex());
    }

    std::cout<<"done"<<std::endl;
    typedef itk::RecursiveGaussianImageFilter< ImageType, ImageType > 
FilterType;

    FilterType::Pointer filterX = FilterType::New();
    FilterType::Pointer filterY = FilterType::New();
    FilterType::Pointer filterZ = FilterType::New();

    filterX->SetNormalizeAcrossScale(false);
    filterY->SetNormalizeAcrossScale(false);
    filterZ->SetNormalizeAcrossScale(false);

    filterX->SetDirection(0);
    filterY->SetDirection(1);
    filterZ->SetDirection(2);

    filterX->SetInput(vol);
    filterY->SetInput(filterX->GetOutput());
    filterZ->SetInput(filterY->GetOutput());

    ImageDuplicatorType::Pointer duplicator = ImageDuplicatorType::New();

    duplicator->SetInputImage( filterZ->GetOutput() );

    std::cout<<"Getting pdx ...";
    filterX->SetOrder(FilterType::FirstOrder);
    filterY->SetOrder(FilterType::ZeroOrder);
    filterZ->SetOrder(FilterType::ZeroOrder);

    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdx = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

    //pdy

    std::cout<<"Getting pdy ...";
    filterX->SetOrder(FilterType::ZeroOrder);
    filterY->SetOrder(FilterType::FirstOrder);
    filterZ->SetOrder(FilterType::ZeroOrder);
    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdy = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

    //pdz

    std::cout<<"Getting pdz ...";

    filterX->SetOrder(FilterType::ZeroOrder);
    filterY->SetOrder(FilterType::ZeroOrder);
    filterZ->SetOrder(FilterType::FirstOrder);
    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdz =duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

    //pdxy

    std::cout<<"Getting pdxy ...";

    filterX->SetOrder(FilterType::FirstOrder);
    filterY->SetOrder(FilterType::FirstOrder);
    filterZ->SetOrder(FilterType::ZeroOrder);

    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdxy = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

        //pdxz

    std::cout<<"Getting pdxz ...";
    filterX->SetOrder(FilterType::FirstOrder);
    filterY->SetOrder(FilterType::ZeroOrder);
    filterZ->SetOrder(FilterType::FirstOrder);

    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();
    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdxz = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

        //pdyz

    std::cout<<"Getting pdyz ...";

    filterX->SetOrder(FilterType::ZeroOrder);
    filterY->SetOrder(FilterType::FirstOrder);
    filterZ->SetOrder(FilterType::FirstOrder);
    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();
    std::cout<<"getting output ...";

    ImageType::Pointer pdyz = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

        //pdxx

    std::cout<<"Getting pdxx ...";

    filterX->SetOrder(FilterType::SecondOrder);
    filterY->SetOrder(FilterType::ZeroOrder);
    filterZ->SetOrder(FilterType::ZeroOrder);
    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdxx = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

        //pdyy

    std::cout<<"Getting pdyy ...";

    filterX->SetOrder(FilterType::ZeroOrder);
    filterY->SetOrder(FilterType::SecondOrder);
    filterZ->SetOrder(FilterType::ZeroOrder);

    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdyy = duplicator->GetOutput();
    std::cout<<"done"<<std::endl;

    //pdzz

    std::cout<<"Getting pdzz ...";

    filterX->SetOrder(FilterType::ZeroOrder);
    filterY->SetOrder(FilterType::ZeroOrder);
    filterZ->SetOrder(FilterType::SecondOrder);

    filterX->SetSigma(1.0);
    filterY->SetSigma(1.0);
    filterZ->SetSigma(1.0);

    std::cout<<"filter updating ...";

    filterZ->Update();

    duplicator->Update();

    std::cout<<"getting output ...";

    ImageType::Pointer pdzz = duplicator->GetOutput();

    std::cout<<"done"<<std::endl;

    PartialDerivative PD;

    unsigned long pos;

    for( i = 0; i < s.size(); ++i )
    {
        pixelIndex = s.at(i);

        //first order derivatives
        PD.x = pdx->GetPixel(pixelIndex);
        PD.y = pdy->GetPixel(pixelIndex);
        PD.z = pdz->GetPixel(pixelIndex);

        //second order derivatives
        PD.xy = pdxy->GetPixel(pixelIndex);
        PD.xz = pdxz->GetPixel(pixelIndex);
        PD.yz = pdyz->GetPixel(pixelIndex);

        PD.xx = pdxx->GetPixel(pixelIndex);
        PD.yy = pdyy->GetPixel(pixelIndex);
        PD.zz = pdzz->GetPixel(pixelIndex);

        std::cout<<"x "<<pixelIndex[0]<<" y "<<pixelIndex[1]<<" z 
"<<pixelIndex[2]<<std::endl;
        std::cout<<"pdx "<<PD.x<<" pdy "<<PD.y<<" pdz "<<PD.z<<std::endl;
        std::cout<<"pdxy "<<PD.xy<<" pdxz "<<PD.xz<<" pdyz 
"<<PD.yz<<std::endl;
        std::cout<<"pdxx "<<PD.xx<<" pdyy "<<PD.yy<<" pdzz 
"<<PD.zz<<std::endl;

    }

    return 0;
}