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