[Insight-users] 2D deformable registration
Luis Ibanez
luis.ibanez@kitware.com
Thu May 20 00:31:15 EDT 2004
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);
> 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
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list