[Insight-users] 2D deformable registration

[ping chen]

hi Luis,


i add your code in, but it runs into such errors, 

Building dependencies cmake.check_depends...
Building object file fieldreader.o...
/Users/Ping/Desktop/field2/fieldreader.cxx: In
function `int main(int, 
   char**)':
/Users/Ping/Desktop/field2/fieldreader.cxx:86: error:
no match for `* 
   itk::ImageRegionIterator<DeformationFieldType>&'
operator
/Users/Ping/Insight/Utilities/vxl/core/vnl/vnl_transpose.h:69:
error: candidates
   are: vnl_matrix<double> operator*(const
vnl_matrix<double>&, const 
   vnl_transpose&)
/Users/Ping/Desktop/field2/fieldreader.cxx:87: error:
no match for `* 
   itk::ImageRegionIterator<DeformationFieldType>&'
operator
/Users/Ping/Insight/Utilities/vxl/core/vnl/vnl_transpose.h:69:
error: candidates
   are: vnl_matrix<double> operator*(const
vnl_matrix<double>&, const 
   vnl_transpose&)


the error happens where we add 
  const double x = fabs( (*it)[0] );
  const double y = fabs( (*it)[1] );

i dont know how to use iterator to get specific voxel
values but that doesnt matter, instead i just use
it.Get( ) to get all the values and save in a logfile 
so below is the code i add to the original code

DeformationFieldType::Pointer  fieldimage2=
DeformationFieldType::New();
 fieldimage2 = fieldreader1->GetOutput();

double maxx = 0;
double maxy = 0;
itk::ImageRegionIterator<
           DeformationFieldType >  it( fieldimage2,
               fieldimage2->GetBufferedRegion() );

it.GoToBegin();

while( !it.IsAtEnd() )
   {
//   const double x = fabs( (*it)[0] );
 //  const double y = fabs( (*it)[1] );
//   if( x > maxx ) maxx = x ;
//   if( y > maxy ) maxy = y ;
  std::cout << "vector values = " << it.Get( ) <<
std::endl; 
   ++it;
   }

std::cout << "Max fabs X = " << maxx << std::endl;
std::cout << "Max fabs Y = " << maxy << std::endl;

and run the program as 
./fieldreader VectorDeformationField.mhd
BrainProtonDensitySliceBorder20.png try.img >> logfile

when checking the logfile, we can see the vector
vaules are all zeros. 


vector values = [0, 0]
vector values = [0, 0]
vector values = [0, 0]
vector values = [0, 0]
vector values = [0, 0]
vector values = [0, 0]


Thanks,
Ping 


--- Luis Ibanez <luis.ibanez@kitware.com> wrote:
> 
> Hi Ping,
> 
> Your code looks fine.
> 
> 
> Note that you are applying a deformable registration
> methods to a pair of images that are only different
> in their translation.
> 
> The ranges that you obtain from ParaView seems to
> be ok with the expected translations: 13 in X and 17
> in Y.
> 
> Please try the following:
> 
> After updating the reader of the deformation field,
> 
> add a while loop and visit all the pixels in the
> vector field, and compute the max abs value of
> the components.
> 
> The fact is that if you still are observing that
> the warped image looks exactly the same as the
> input image... that means that the vector field
> is null... and there may be still some problem
> when reading vector images.
> 
> Your experiment with the while loop will help
> to decide the question.
> 
> The loop will look like
> 
> double maxx = 0;
> double maxy = 0;
> itk::ImageRegionIterator<
>            DeformationFieldType >  it( fieldimage,
>                fieldimage->GetBufferedRegion() );
> 
> it.GoToBegin();
> 
> while( !it.IsAtEnd() )
>    {
>    const double x = fabs( (*it)[0] );
>    const double y = fabs( (*it)[1] );
>    if( x > maxx ) maxx = x ;
>    if( y > maxy ) maxy = y ;
>    ++it;
>    }
> 
> std::cout << "Max fabs X = " << maxx << std::endl;
> std::cout << "Max fabs Y = " << maxy << std::endl;
> 
> 
> 
> 
> Please let us know what you find.
> 
> 
>     Thanks
> 
> 
> 
>         Luis
> 
> 
> 
> 
> -----------------
> ping chen wrote:
> 
> > Hi Luis,
> > 
> > I dont understand when you say 
> > 
> >  Note that the values of the deformation
> > 
> >>field vectors are in physical units,
> >>not pixels, therefore you have to
> >>interpret them in the context of the
> >>pixel spacing in your image.
> > 
> > 
> > 
> > in DeformableRegistration1, the program will
> produce a
> > deformationation field which is
> > VectorDeformationField.mhd     
> > VectorDeformationField.raw 
> > which warps BrainProtonDensitySliceBorder20.png as
> > moving image to
> > BrainProtonDensitySliceShifted13x17y.png as fixed
> > image. 
> > so for this 2d deformableregistration, i didnt
> make
> > any change, i suppose the deformationfield should
> be
> > all right. i have viewed the deformationfield in
> the
> > paraview, it is nonzeros. 
> > 
> > bounds:
> > x range 0 to 12.9
> > y range 0 to 15.9 
> > z range 0 to 0. 
> >  
> > 
> > then in my code, i just create the
> VectorImageReader,
> > and read in this VectorDeformationField.mhd file,
> and
> > then update, and then GetOutput to set the
> deformation
> > field in the warper. thank you very much for your
> > suggestions. 
> > 
> > -Ping 
> > 
> > 
> > 
> > below is my code 
> > 
> > //  Software Guide : BeginLatex
> > //
> > //  The following code is an implementation of a
> small
> > Insight
> > //  program. It tests including header files and
> > linking with ITK
> > //  libraries.
> > //
> > //  Software Guide : EndLatex 
> > 
> > // Software Guide : BeginCodeSnippet
> > #include "itkImage.h"
> > #include "itkImageFileReader.h" 
> > #include "itkImageFileWriter.h" 
> > #include "itkWarpImageFilter.h"
> > #include "itkLinearInterpolateImageFunction.h"
> > #include <iostream>
> > 
> > typedef itk::Image<unsigned char, 2>
> InputImageType;
> > typedef itk::Image<float, 2>	OutputImageType;
> > 
> > typedef itk::ImageFileReader< InputImageType >
> > InputReaderType;
> > 
> > typedef itk::Vector< float, 2 >   
> VectorPixelType;
> > typedef itk::Image<  VectorPixelType, 2 >
> > DeformationFieldType;
> > 
> > typedef itk::ImageFileReader< DeformationFieldType
> >
> > FieldReaderType;
> > 
> > typedef itk::WarpImageFilter<
> >                           InputImageType, 
> >                           OutputImageType,
> >                           DeformationFieldType  > 
>   
> > WarperType;
> > typedef itk::LinearInterpolateImageFunction<
> >                                    InputImageType,
> >                                    double         
> > 
> > InterpolatorType;
> > 					
> > 
> > int main(int argc, char *argv[])
> > {
> >   char *fieldimage;
> >   char *towarpimage;
> >   char *warpedimage;
> >   
> >   if ( argc < 4 )
> >     {
> >     std::cout << "Parameter file name missing" <<
> > std::endl;
> > 	std::cerr << " fieldimage towarpimagename
> > warpedimagename";
> >     std::cout << "Usage: " << argv[0] << "
> param.file"
> > << std::endl;
> > 	
> >     return -1;
> >     } 
> >   else 
> >     { 
> >     fieldimage=argv[1];
> > 	towarpimage = argv[2];
> > 	warpedimage = argv[3];
> >     }
> > 	
> > FieldReaderType::Pointer fieldreader1 =
> > FieldReaderType::New();
> > InputReaderType::Pointer towarpreader =
> > InputReaderType::New();
> > 
> > fieldreader1->SetFileName(fieldimage);
> 
=== message truncated ===


	
		
__________________________________
Do you Yahoo!?
Yahoo! Domains – Claim yours for only $14.70/year
http://smallbusiness.promotions.yahoo.com/offer