[ITK Community] [Insight-users] Making a deformation field denser in order to deform points
Yago Diez
yagodiezdonoso at gmail.com
Sat Dec 7 13:23:35 EST 2013
Mat, thank you very much for your answer, my code seems to be working and
it looks like this:
I had some additional problems with the itk::DisplacementFieldTransform
class when trying to use it with floats but things seem to be working with
doubles. I do not include yet the inversion of the deformation fields that
Brian suggested but I am pretty confident I can add it without many more
problems.
I am including the code in case anyone is interested, please note that it
is in no way tested (I will update it if I find major problems with it).
Thanks for your help.
Yago
CODE TO DEFORM A POINT WITH AN ITK FORMATTED DEFORMATION FIELD (ITK4)
#include "itkCompositeTransform.h"
#include "itkDisplacementFieldTransform.h"
#include "itkImageFileReader.h"
#include<iostream>
#include <itkVector.h>
#include <itkImage.h>
using namespace std;
int main( int argc, char * argv [] )
{
if( argc < 3 )
{
std::cerr << "Missing command line arguments" << std::endl;
std::cerr << "Usage : DeformationFieldFile pointX pointY pointZ";
return -1;
}
typedef double PixelType;
typedef double VectorComponentType;
const unsigned int Dimension = 3;
typedef itk::Vector< PixelType, Dimension > VectorPixelType;
typedef itk::Image< VectorPixelType, Dimension > DeformationFieldType;
typedef itk::ImageFileReader< DeformationFieldType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(argv[1]);
reader->Update();
DeformationFieldType::Pointer deffield = reader->GetOutput();
typedef itk::DisplacementFieldTransform<VectorComponentType, 3>
DisplacementFieldTransformType;
DisplacementFieldTransformType::Pointer displacementFieldTransform =
DisplacementFieldTransformType::New();
displacementFieldTransform->SetDisplacementField( reader->GetOutput() );
// read the point to be transformed
itk::Point<double,3> p0,p1;
p0[0] = atof(argv[2]);
p0[1] = atof(argv[3]);
p0[2] = atof(argv[4]);
typedef itk::CompositeTransform<double, Dimension>
CompositeTransformType;
typename CompositeTransformType::InputPointType point_in;
typename CompositeTransformType::OutputPointType point_out;
typename CompositeTransformType::Pointer compositeTransform =
CompositeTransformType::New();
compositeTransform->AddTransform( displacementFieldTransform);
p1=compositeTransform->TransformPoint( p0 );
//cout<<"initial point "<<"p0 "<<p0<<"p1
"<<compositeTransform->TransformPoint( p0 )<<endl;
cout<<p1[0]<<" "<<p1[1]<<" "<<p1[2]<<endl;
return 0;
}
Yago Diez Donoso (PhD)
Computer Vision and Robotics Group http://vicorob.udg.es/
E-mail: yago at eia.udg.es; yagodiezdonoso at gmail.com
Phone: (int. code) 34 972418013
University of Girona
On Sat, Dec 7, 2013 at 6:23 PM, Matt McCormick
<matt.mccormick at kitware.com>wrote:
> Hi Yago,
>
> An Image that defines the displacement field vectors can be set on a
> DisplacementFieldTransfrom with SetDisplacementFieldTransform [1]. It
> is also possible to save the displacement field as a transform
> directly, instead of an image, with the TransformFileWriter [2].
>
> Hope this helps,
> Matt
>
> [1]
> http://www.itk.org/Doxygen/html/classitk_1_1DisplacementFieldTransform.html#ab0b5ac8b6792ed0acddbd0a6666b39d2
> [2]
> http://www.itk.org/Doxygen/html/classitk_1_1TransformFileWriterTemplate.html
>
> On Sat, Dec 7, 2013 at 10:48 AM, Yago Diez <yagodiezdonoso at gmail.com>
> wrote:
> > First of all, thanks a lot for your answer and sorry I couldn't solve the
> > problem on my own. I have been fighting for a while with this but so far
> I
> > have not had much success. At least I am now a bit more aware of how many
> > things I don't know.
> >
> > First of all, I realise what I want to do boils down to :
> >
> > typename CompositeTransformType::Pointer compositeTransform =
> > CompositeTransformType::New();
> > compositeTransform->AddTransform( deffield);
> > compositeTransform->TransformPoint( p0 );
> >
> > The problem now is that my deformation field is not recognised as a valid
> > transformation.
> >
> > My deformation field which is the output of demons registration and it is
> > stored in a file in mha format. I read it like this (as a vector image as
> > used to be done in itk3):
> >
> >> typedef float PixelType;
> >> typedef float VectorComponentType;
> >>
> >> const unsigned int Dimension = 3;
> >>
> >> typedef itk::Vector< PixelType, Dimension > VectorPixelType;
> >> typedef itk::Image< VectorPixelType, Dimension > DeformationFieldType;
> >> typedef itk::ImageFileReader< DeformationFieldType > ReaderType;
> >>
> >> ReaderType::Pointer reader = ReaderType::New();
> >> reader->SetFileName(argv[1]);
> >> reader->Update();
> >> DeformationFieldType::Pointer deffield = reader->GetOutput();
> >
> >
> >
> > When I try to add it to the composite transform, however, I get the
> error:
> >
> > /media/disc4/code/newTransformPoints.cxx:112:47: error: no matching
> function
> > for call to ‘itk::CompositeTransform<float,
> > 3u>::AddTransform(itk::Image<itk::Vector<float, 3u>, 3u>::Pointer&)’
> > /media/disc4/code/newTransformPoints.cxx:112:47: note: candidate is:
> > /usr/local/include/ITK-4.4/itkMultiTransform.h:145:16: note: void
> > itk::MultiTransform<TScalar, NDimensions,
> > NSubDimensions>::AddTransform(itk::MultiTransform<TScalar, NDimensions,
> > NSubDimensions>::TransformType*) [with TScalar = float, unsigned int
> > NDimensions = 3u, unsigned int NSubDimensions = 3u,
> > itk::MultiTransform<TScalar, NDimensions, NSubDimensions>::TransformType
> =
> > itk::Transform<float, 3u, 3u>]
> > /usr/local/include/ITK-4.4/itkMultiTransform.h:145:16: note: no known
> > conversion for argument 1 from ‘itk::Image<itk::Vector<float, 3u>,
> > 3u>::Pointer {aka itk::SmartPointer<itk::Image<itk::Vector<float, 3u>,
> 3u>
> >>}’ to ‘itk::MultiTransform<float, 3u, 3u>::TransformType* {aka
> > itk::Transform<float, 3u, 3u>*}’
> >
> >
> > From what I understand, instead of reading my deformation field as a
> vector
> > image, I should probably be using the itk4 "displacementfieldtransform".
> >
> >> typedef itk::DisplacementFieldTransform<VectorComponentType, 3>
> >> DeformationFieldTransformType;
> >
> >
> > My main problem now is that I cannot seem to "connect" my mha file
> > containing the deformation field with the
> displacementFieldTransformType. I
> > probably need some type of basic "read DisplacementFieldTransform from
> file"
> > function, but so far I have not been able to make one.
> >
> > Any ideas that might help? I realise that the bigger problem here is my
> own
> > lack of knowledge on the details of itk4, so please do not hesitate to
> point
> > me to as basic a source as necessary.
> > Gratefully
> > Yago Diez Donoso
> >
> >
> >
> > Yago Diez Donoso (PhD)
> > Computer Vision and Robotics Group http://vicorob.udg.es/
> > E-mail: yago at eia.udg.es; yagodiezdonoso at gmail.com
> > Phone: (int. code) 34 972418013
> > University of Girona
> >
> >
> > On Wed, Dec 4, 2013 at 8:48 PM, brian avants <stnava at gmail.com> wrote:
> >>
> >> see line 142 here:
> >>
> >>
> >>
> https://github.com/stnava/ANTs/blob/d99d9c9e23ae17c116377a248cd370ddbd8c2599/Examples/antsApplyTransformsToPoints.cxx
> >>
> >> you can do the same thing in your application but note that the
> >> deformation field , when applied to a point , should be the inverse of
> the
> >> deformation field applied to an image
> >>
> >> see my previous email regarding SyN and inverse warps.
> >>
> >>
> >> brian
> >>
> >>
> >>
> >>
> >> On Wed, Dec 4, 2013 at 2:33 PM, Yago Diez <yagodiezdonoso at gmail.com>
> >> wrote:
> >>>
> >>> Hi everyone,
> >>>
> >>> I have a problem regarding itk deformation fields:
> >>>
> >>> Specifically, I want to know where a particular point ends up after
> >>> registration. Said registration has produced a deformation field. I
> have
> >>> already used this deformation field to deform the source image and
> produce
> >>> the output image for registration (this was done using
> itkWarpImageFilter).
> >>>
> >>> My problem is that at this moment I only want to transform one point.
> The
> >>> obvious solution was to read the corresponding voxel in the deformation
> >>> field and add it to the point, but the problem is that the deformation
> field
> >>> is not "dense enough", so in most cases, the read deformation for a
> point
> >>> ends up being zero (see attached image for a screen capture of a 2D
> slice of
> >>> my 3D deformation field).
> >>>
> >>> I have noticed how when I deform an image with itkWarpImageFilter,
> there
> >>> is an "extra" interpolation step which would avoid the problem I am
> having,
> >>> but so far I have not been able to interpolate a deformation field. Any
> >>> ideas on how to do that?
> >>>
> >>> Thank you
> >>>
> >>>
> >>>
> >>>
> >>> Yago Diez Donoso (PhD)
> >>> Computer Vision and Robotics Group http://vicorob.udg.es/
> >>> E-mail: yago at eia.udg.es; yagodiezdonoso at gmail.com
> >>> Phone: (int. code) 34 972418013
> >>> University of Girona
> >>>
> >>> _____________________________________
> >>> Powered by www.kitware.com
> >>>
> >>> Visit other Kitware open-source projects at
> >>> http://www.kitware.com/opensource/opensource.html
> >>>
> >>> Kitware offers ITK Training Courses, for more information visit:
> >>> http://www.kitware.com/products/protraining.php
> >>>
> >>> Please keep messages on-topic and check the ITK FAQ at:
> >>> http://www.itk.org/Wiki/ITK_FAQ
> >>>
> >>> Follow this link to subscribe/unsubscribe:
> >>> http://www.itk.org/mailman/listinfo/insight-users
> >>>
> >>
> >
> >
> > _____________________________________
> > Powered by www.kitware.com
> >
> > Visit other Kitware open-source projects at
> > http://www.kitware.com/opensource/opensource.html
> >
> > Kitware offers ITK Training Courses, for more information visit:
> > http://www.kitware.com/products/protraining.php
> >
> > Please keep messages on-topic and check the ITK FAQ at:
> > http://www.itk.org/Wiki/ITK_FAQ
> >
> > Follow this link to subscribe/unsubscribe:
> > http://www.itk.org/mailman/listinfo/insight-users
> >
> > _______________________________________________
> > Community mailing list
> > Community at itk.org
> > http://public.kitware.com/cgi-bin/mailman/listinfo/community
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20131207/8e0230e7/attachment-0001.htm>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users
More information about the Community
mailing list