[Insight-users] Copy filter output
Silvano Agliozzo
agliozzo at isiosf.isi.it
Mon, 19 Jan 2004 19:17:55 +0100 (CET)
Hi Luis,
I update the new version of the ImageDuplicator class.
It works even though it produces the following error each time I
calculated a derivative.B
WARNING: In /usr/local/itk/Insight/Code/Common/itkImageDuplicator.txx,
line 52
Self (0x8099fe0): Input image has not been modified
I suppose this is due to the modification you did.
In the documentation of the RecursiveGaussianImageFilter, the derivative
order of this filter is set by the method SetOrder(OrderEnumType _arg), but in
the example
Insight/Examples/Filtering/SecondDerivativeRecursiveGaussianImageFilter.cxx
the order is set differently by the SetZeroOrder(), SetFirstOrder(),
SetSecondOrder() methods. I checked both way ot setting and I did not
check any changes, but I would like to know the difference,
Moreover I would like to point out that there could beB a mistake in the
same example.
For calculating the mixed second order derivative, e.g. Ixy,
for a 3D image, one should set the filters filterX, filterY, filterZ
in the following way:
filterX->SetFirstOrder();
filterY->SetFirstOrder();
filterZ->SetZeroOrder();
Conversely, in the example this derivative is calculated by setting:
filterX->SetSecondOrder();
filterY->SetSecondOrder();
filterZ->SetZeroOrder();
Theoretically, if you set this derivative orders, you get a mixed forth
order derivative. I did not get insight the implementation of the code,
and may-be you have an explanation for that.
Concerning the calculated values, they are close to the expected
ones (I calculated a spherical and cilinder intensity fields), a part the
second order derivative about z, i.e. Izz, which is always smaller than
Ixx and Iyy.
a+
Silvano
On Thu, 15 Jan 2004, Luis Ibanez wrote:
>
> Hi Silvano,
>
> It seems that you found an interesting Bug.
>
> For some reason the pipeline is not updating
> the modified time of the output image in the
> RecursiveGaussianImageFilter.
>
> This has been entered as Bug#517.
> http://www.itk.org/Bug/bug.php?op=show&bugid=517&pos=11
>
> In the meantime, the following actions have
> been taken to make your life easier:
>
>
> 1) The check for the modified time has been
> temporarily removed from the ImageDuplicator
> class. So it will always duplicate the image
> regardless of the time stamp. This file is in
>
> Insight/Code/Common/itkImageDuplicator.txx
>
>
> 2) The full code example of the Second derivative
> computation has been commited under
>
>
> Insight/Examples/Filtering
> SecondDerivativeRecursiveGaussianImageFilter.cxx
>
>
>
> You may simply update these two files in your CVS
> checkout and give them a try.
>
>
>
> Please let us know if you find any further problems,
>
>
> Thanks
>
>
> Luis
>
>
>
> --------------------------
> Silvano Agliozzo wrote:
>
> > Hi Luis,
> >
> > thank you for your answer. However I still have a problem.
> >
> > I've used the duplicator as you suggested, because I have to process
> > large 3D images (CT scans). Nevertheless I still get the same results
> > for all the different images even though I set the
> > RecursiveGaussianImageFilters to calculate different derivatives.
> > Moreover the values of the filter outputs do not correspond to the
> > expected one.
> > If you don't mind, I wite hereafter some line of the code:
> >
> > ....
> >
> > 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 = data to be processedA
> > vol->SetRegions(region);
> > vol->Allocate();
> > vol->FillBuffer( 0.0 );
> > vol->SetSpacing( spacing );
> > vol->SetOrigin( origin );
> >
> > pdxx->SetRegions(region);// second order derivative
> > ...
> > ...
> >
> > the same for pdyy, pdzz, pdxy, pdxz, pdyz
> >
> > 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());
> >
> > filterX->SetSigma(1.0);
> > filterY->SetSigma(1.0);
> > filterZ->SetSigma(1.0);
> >
> > typedef itk::ImageDuplicator< ImageType > ImageDuplicatorType;
> > ImageDuplicatorType::Pointer duplicator = ImageDuplicatorType::New();
> >
> > std::cout<<"Getting pdxx ...";
> > filterX->SetOrder(FilterType::SecondOrder);
> > filterY->SetOrder(FilterType::ZeroOrder);
> > filterZ->SetOrder(FilterType::ZeroOrder);
> >
> > std::cout<<"filter updating ...";
> > filterZ->Update();
> >
> > duplicator->SetInputImage( filterZ->GetOutput() );
> > std::cout<<"duplicator updating ...";
> > duplicator->Update();
> >
> > pdxx = duplicator->GetOutput();
> >
> > pdyy ...
> >
> > std::cout<<"Getting pdzz ...";
> > filterX->SetOrder(FilterType::ZeroOrder);
> > filterY->SetOrder(FilterType::SecondOrder);
> > filterZ->SetOrder(FilterType::ZeroOrder);
> >
> > std::cout<<"filter updating ...";
> > filterZ->Update();
> >
> > std::cout<<"duplicator updating ...";
> > duplicator->Update();
> >
> > pdzz = duplicator->GetOutput();
> >
> > std::cout<<"Getting pdxy ...";
> > filterX->SetOrder(FilterType::FirstOrder);
> > filterY->SetOrder(FilterType::FirstOrder);
> > filterZ->SetOrder(FilterType::ZeroOrder);
> >
> > std::cout<<"filter updating ...";
> > filterZ->Update();
> >
> > std::cout<<"duplicator updating ...";
> > duplicator->Update();
> >
> > pdxy = duplicator->GetOutput();
> >
> > ....and so on for the rest ...(pdxz, pdyz)
> >
> >
> > but pdxx, pdyy, ..., pdyz are the same.
> > I also tried to insert the method Modified() before updating the filter
> > without get any changes.
> >
> >
> > thank you in advance,
> >
> > a+
> >
> > Silvano
> > On Wed, 14 Jan 2004, Luis Ibanez wrote:
> >
> >
> >>Hi Silvano,
> >>
> >>Welcome to ITK !
> >>
> >> From your message it looks like you are trying to
> >>create one RecursiveGaussianImageFilter and
> >>reuse it in order to compute multiple second
> >>derivatives.
> >>
> >>Reusing filters is a bad practice in a data
> >>pipepline system as ITK or VTK.
> >>
> >>Don't be affraid to be generous in creating ITK
> >>filters (unless your images are really large).
> >>
> >>
> >>
> >>I assume that you are trying to compute the full
> >>set of second derivatives in a 3D image, which
> >>is equivalent to compute the Hessian Matrix for
> >>every pixel
> >>
> >>
> >> Ixx Ixy Ixz
> >> Iyx Iyy Iyz
> >> Izx Izy Izz
> >>
> >>
> >>given that this is a symmetric matrix, you only
> >>need to compute the 6 images:
> >>
> >> { Ixx Iyy Izz Ixy Ixz Iyz }
> >>
> >>
> >>Let's take Izz first
> >>
> >>
> >>You need 3 RecursiveGaussianImageFilters in order
> >>to compute the Izz image.
> >>
> >>FilterType::Pointer ga = FilterType::New();
> >>FilterType::Pointer gb = FilterType::New();
> >>FilterType::Pointer gc = FilterType::New();
> >>
> >>ga->SetDirection( 0 );
> >>gb->SetDirection( 1 );
> >>gc->SetDirection( 2 );
> >>
> >>ga->SetZeroOrder();
> >>gb->SetZeroOrder();
> >>gc->SetSecondOrder();
> >>
> >>ga->SetInput( inputImage );
> >>gb->SetInput( ga->GetOutput() );
> >>gc->SetInput( gb->GetOutput() );
> >>
> >>gc->Update();
> >>
> >>At this point you have Izz as the output of the
> >>filter "gc".
> >>
> >>You can copy this output out of the pipeline
> >>by using the ImageDuplicator class
> >>http://www.itk.org/Insight/Doxygen/html/classitk_1_1ImageDuplicator.html
> >>
> >>DuplicatorType::Pointer duplicator = DuplicatorType::New();
> >>
> >>duplicator->SetInput( gc->GetOutput() );
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Izz = duplicator->GetOutput();
> >>
> >>
> >>Note the "duplicator" is NOT an ITK filter.
> >>It doesn't to make part of the pipeline,
> >>despite the fact that its API look like
> >>a filter.
> >>
> >>
> >>
> >>
> >>Now, you can readjust the pipeline in order
> >>to compute Iyy, just do
> >>
> >>gc->SetDirection( 1 ); // gc now works along Y
> >>gb->SetDirection( 2 ); // gb now works along Z
> >>
> >>gc->Update();
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Iyy = duplicator->GetOutput();
> >>
> >>
> >>
> >>
> >>In order to get Ixx you do
> >>
> >>gc->SetDirection( 0 ); // gc now works along X
> >>ga->SetDirection( 1 ); // ga now works along Y
> >>
> >>gc->Update();
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Ixx = duplicator->GetOutput();
> >>
> >>
> >>
> >>
> >>
> >>
> >>For the cross derivatives you do
> >>
> >>ga->SetDirection( 0 );
> >>gb->SetDirection( 1 );
> >>gc->SetDirection( 2 );
> >>
> >>ga->SetZeroOrder();
> >>gb->SetSecondOrder();
> >>gc->SetSecondOrder();
> >>
> >>gc->Update();
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Iyz = duplicator->GetOutput();
> >>
> >>
> >>ga->SetDirection( 1 );
> >>gb->SetDirection( 0 );
> >>gc->SetDirection( 2 );
> >>
> >>ga->SetZeroOrder();
> >>gb->SetSecondOrder();
> >>gc->SetSecondOrder();
> >>
> >>gc->Update();
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Ixz = duplicator->GetOutput();
> >>
> >>
> >>ga->SetDirection( 2 );
> >>gb->SetDirection( 0 );
> >>gc->SetDirection( 1 );
> >>
> >>ga->SetZeroOrder();
> >>gb->SetSecondOrder();
> >>gc->SetSecondOrder();
> >>
> >>gc->Update();
> >>duplicator->Update();
> >>
> >>SecondDerivativeImageType::Pointer Ixy = duplicator->GetOutput();
> >>
> >>
> >>
> >>
> >>
> >>Regards,
> >>
> >>
> >> Luis
> >>
> >>
> >>------------------------
> >>Silvano Agliozzo wrote:
> >>
> >>
> >>>Dear all,
> >>>
> >>>I'am a newcomer to the Itk libraries, so I could ask a trivial and I
> >>>apologize in advance.
> >>>
> >>>I need to use the RecursiveGuassianImageFilter in a 3D image in order to
> >>>obtain second order derivatives, but I need to store the output of the filter
> >>>to different 3D images. I do not want to get the pointer of the filter
> >>>output since each time the pipeline of the filter is updated, the voxel
> >>>values of all the images previously calculated, are updated too.
> >>>
> >>>I would like to know how can I perform this task ?
> >>>
> >>>thank you very much,
> >>>
> >>>Silvano
> >>>
> >>>_______________________________________________
> >>>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
> >
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>