[Insight-users] RecursiveGaussianImageFilter (fwd)
Silvano Agliozzo
agliozzo at isiosf.isi.it
Tue, 10 Feb 2004 20:32:49 +0100 (CET)
Hi Luis,
I think too that this could be also related to the smoothing effect.
Actually, I had to spend some time in order to find the right sigma values
of the filters. However I wonder how much the output of this filter
depends on the image values which the filter is applied to. In this case
its setting is very delicate. thank you anyway,
a+
Silvano
Hi Silvano,
Thanks a lot for posting your
compact example.
The good news:
We were able to reproduce the effect that you
reported. That is, the values of the derivatives
change when you add a constant to the input image.
The bad news:
We still have to figure out why ...
and... how to fix it.
My guess is that this effect is related to the
fact that this filter is not just doing derivatives
but smoothing followed by derivatives.
This issue has been entered in the Bug database
under ID # 584
http://www.itk.org/Bug/bug.php?op=show&bugid=584&pos=0
If you want to get updates on the status of this bug, you
are welcom to create an account on the bug tracker at
http://www.itk.org/Bug/
You will just need to provide an email address and
make up your own password.
Regards,
Luis
-------------------------
Silvano Agliozzo wrote:
> 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;
> }
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users