[Insight-users] 2D deformable registration

[Luis Ibanez]

Hi Ping,

You are right, the Get() method is what
should be used to retrieve the pixel
values.

It was my mistake to put (*it) in the
previous email.

Sorry about that.

---

Your test shows that the reading
of vector images may still have some
problems. We may have to reopen that
bug in the database...



    Luis



----------------
ping chen wrote:

> 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 
>