[Insight-users] 2D deformable registration

[ping chen]

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);
towarpreader->SetFileName(towarpimage);

try
    {
    fieldreader1->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during inputimage1
reading " << std::endl;
    std::cerr << e << std::endl;
    return -1;
    }

	try
    {
    towarpreader->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during towarp
reading " << std::endl;
    std::cerr << e << std::endl;
    return -1;
    }
 
  WarperType::Pointer warper = WarperType::New();
  InterpolatorType::Pointer interpolator =
InterpolatorType::New();
  warper->SetInterpolator( interpolator );
  warper->SetOutputSpacing(
towarpreader->GetOutput()->GetSpacing() );
  warper->SetOutputOrigin(
towarpreader->GetOutput()->GetOrigin() );
 
warper->SetDeformationField(fieldreader1->GetOutput());
  warper->SetInput( towarpreader->GetOutput() );
  	try
    {
       warper->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during warper
process " << std::endl;
    std::cerr << e << std::endl;
    return -1;
    }

	
	
  itk::ImageFileWriter<OutputImageType>::Pointer
writer;
  writer =
itk::ImageFileWriter<OutputImageType>::New();
  writer->SetFileName(warpedimage);
  writer->SetInput(warper->GetOutput()); 
  
   	try
    {
         writer->Write();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during warped
writing " << std::endl;
    std::cerr << e << std::endl;
    return -1;
    }


	
  return 0;
}
//  Software Guide : EndLatex 





--- Luis Ibanez <luis.ibanez@kitware.com> wrote:
> 
> Hi Ping,
> 
> 
> Please do the following:
> 
> 
> 1) Go to
> 
>         http://www.paraview.org
> 
>     download and intall the binary
>     version of ParaView corresponding
>     to your system.
> 
> 
> 2) Run ParaView and load your deformation
>     field. (It must be in MetaImage format).
> 
> 
> 3) Use the Probes in ParaView in order
>     to check if your deformation field
>     has any values different from zero.
> 
> 
> 
> It looks like your deformation field
> is being generated with zero values
> or at least, very low values. 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.
> 
> 
> 
> Please let us know what you find.
> 
> 
> 
>    Thanks.
> 
> 
> 
>      Luis
> 
> 
> 
> 
> ---------------
> ping chen wrote:
> 
> > Hi Luis,
> > 
> > I just CVS Insight yesterday and reinstall the
> > Insight, I Compile the code with the newly
> installed
> > ITK, the warped image is still the same as input
> > image. do you have any suggestions? Thanks.
> > 
> > -Ping 
> > 
> > 
> > 
> > 
> > --- Luis Ibanez <luis.ibanez@kitware.com> wrote:
> > 
> >>Hi Ping,
> >>
> >>A quick question before we go into details...
> >>
> >>Did you tried this process with a CVS updated
> >>*this week*  ?
> >>
> >>A bug was recently found in the process of
> >>reading Vector images. The effect is that
> >>when reading an image of vectors, such as a
> >>deformation field, all the pixels turned out
> >>to be null. This bug was fixed at the end
> >>of last week.
> >>
> >>If you are using a version with this bug,
> >>your deformation field is being set as null
> >>vectors in memory and therefore is not
> >>applying *any* deformation to your  data.
> >>
> >>Please try updating your CVS checkout and
> >>let us know if you continue experiencing
> >>any problems.
> >>
> >>
> >>   Thanks,
> >>
> >>
> >>      Luis
> >>
> >>
> >>
> >>-------------------
> >>ping chen wrote:
> >>
> >>>Hi Luis, 
> >>>
> >>>I try to use the deformation field to warp a
> >>
> >>image. i
> >>
> >>>use the deformation field generated by
> >>>DeformableRegistration1.cxx, in which
> displacement
> >>
> >>for
> >>
> >>>x and displacement for y are both generated, also
> >>
> >>a
> >>
> >>>vector displacement image
> >>
> >>VectorDeformationField.mhd
> >>
> >>>is generated. 
> >>>
> >>>the 1st way for me to use the generated x y
> >>>displacement, i use itkCompose2DVectorImageFilter
> >>
> >>to
> >>
> >>>combine x displacement and y displacement into a
> >>>vector image. and then use this vector image to
> >>
> >>set
> >>
> >>>the deformation field for the warper. 
> >>>
> >>>the problem with this way is that the warped
> >>
> >>result
> >>
> >>>from this way is different from the warped image
> >>
> >>from
> >>
> >>>the DeformableRegistration1.cxx, though similiar.
> 
> >>>
> >>>the 2nd way, i just generated the imagereader to
> >>>directly read VectorDeformationField.mhd, and
> then
> >>
> >>use
> >>
> >>>this reader output to set the deformation field
> of
> >>
> >>the
> >>
> >>>warper.
> >>>
> >>>the problem with this way is that this method
> >>
> >>didnt
> >>
> >>>work. the output of warper is exactly the same as
> >>
> >>the
> >>
> >>>moving image. 
> >>>
> >>>below is part of the code: the 1st way
> >>>typedef itk::Compose2DVectorImageFilter<
> >>>              
> >>
> >>InputImageType,DeformationFieldType>
> >>
> >>>FilterType;
> >>> FilterType::Pointer filter = FilterType::New();
> >>> filter->SetInput1(inputreader1->GetOutput());
> >>> filter->SetInput2(inputreader1->GetOutput());
> >>> filter->Update();
> >>>  warper->SetInterpolator( interpolator );
> >>>  warper->SetOutputSpacing(
> >>>towarpreader->GetOutput()->GetSpacing() );
> >>>  warper->SetOutputOrigin(
> >>>towarpreader->GetOutput()->GetOrigin() );
> >>>  warper->SetDeformationField(
> filter->GetOutput()
> >>
> >>);
> >>
> >>>  warper->SetInput( towarpreader->GetOutput() );
> >>>  	try
> >>>    {
> >>>       warper->Update();
> >>>    }
> >>>
> >>>the 2nd way:
> >>>typedef itk::Vector< float, 2 >   
> >>
> >>VectorPixelType;
> 
=== message truncated ===


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