<div>Hi,</div>
<div> </div>
<div>I encountered a problem when trying to deform a 3D image. The way I deformed the image is as the following: I created a 3D contour for an organ, say, lung for this image. Along the contours I have a bunch of points. Then I generated a vector field for these points and used these vectors to deform the image. I didn't have trouble deforming any 2D image. But for 3D image I got an error. I traced down to the point where the problem occured: in the try-catch block with deformer->UpdateLargestPossibleRegion(); debugging reports a problem on line 417 of vnl_matrix.txx file. For reference purpose I copy and paste that part of code here:
</div>
<div>
<p>//: Sets all elements of matrix to specified value. O(m*n).</p>
<p>template<class T><br>void vnl_matrix<T>::fill (T const& value)<br>{<br> for (unsigned int i = 0; i < this->num_rows; i++)<br> for (unsigned int j = 0; j < this->num_cols; j++)<br> this->data[i][j] = value; <========================================= debugging reports problem !
<br>}</p> </div>
<div> </div>
<div> </div>
<div>Also, I am attaching my code of deforming 3D images below. I would greatly appreciate it if anybody can help look into it.</div>
<div> </div>
<div>Best regards,</div>
<div>Ming</div>
<div> </div>
<div>
<p>// ---------------- starting of the code -------------------------<br>#include "itkVector.h"<br>#include "itkImage.h"<br>#include "itkImageRegionIteratorWithIndex.h"<br>#include "itkDeformationFieldSource.h
"<br>#include "itkImageFileReader.h"<br>#include "itkImageFileWriter.h"<br>#include "itkWarpImageFilter.h"<br>#include "itkLinearInterpolateImageFunction.h"</p>
<p>#include "itkNormalVariateGenerator.h"<br>#include "vtkPolyDataReader.h"<br>#include "vtkPolyDataWriter.h"<br>#include "vtkFloatArray.h"<br>#include "vtkPolyData.h"<br>
#include "vtkPointData.h"<br>#include "vtkPoints.h"</p>
<p>#include <fstream></p>
<p>int main( int argc, char * argv[] )<br>{</p>
<p> if( argc < 2 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br> std::cerr << "Usage: " << argv[0];<br> std::cerr << " landmarksFile fixedImage ";
<br> std::cerr << "movingImage deformedMovingImage" << std::endl;<br> return 1;<br> }</p>
<p> const unsigned int Dimension = 3;<br> typedef float VectorComponentType;</p>
<p> typedef itk::Vector< VectorComponentType, Dimension > VectorType;<br> typedef itk::Image< VectorType, Dimension > DeformationFieldType;<br> typedef float PixelType;<br> typedef itk::Image< PixelType, Dimension > ImageType ;
<br> typedef itk::ImageFileReader< ImageType > ReaderType ;<br> typedef itk::ImageFileWriter< ImageType > WriterType ;</p>
<p> ReaderType::Pointer Reader = ReaderType::New();<br> Reader->SetFileName( argv[1] );</p>
<p> try<br> {<br> Reader->Update();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << "Exception thrown " << std::endl;<br> std::cerr << excp << std::endl;
<br> return EXIT_FAILURE;<br> }</p>
<p> WriterType::Pointer Writer = WriterType::New();<br> Writer->SetFileName( argv[2] );</p>
<p> ImageType::ConstPointer fixedImage = Reader->GetOutput();</p>
<p><br> typedef itk::DeformationFieldSource< DeformationFieldType > DeformationSourceType;<br> DeformationSourceType::Pointer deformer = DeformationSourceType::New();<br> deformer->SetOutputSpacing( fixedImage->GetSpacing() );
<br> deformer->SetOutputOrigin( fixedImage->GetOrigin() );<br> deformer->SetOutputRegion( fixedImage->GetLargestPossibleRegion() );</p>
<p> typedef DeformationSourceType::LandmarkContainer LandmarkContainerType;<br> typedef DeformationSourceType::LandmarkPointType LandmarkPointType;<br> LandmarkContainerType::Pointer sourceLandmarks = LandmarkContainerType::New();
<br> LandmarkContainerType::Pointer targetLandmarks = LandmarkContainerType::New();</p>
<p> vtkPolyDataReader *reader = vtkPolyDataReader::New();<br> reader->SetFileName( argv[2] );<br> reader->Update();<br> vtkPolyData *poly = vtkPolyData::New();<br> poly->DeepCopy( reader->GetOutput());<br> reader->Delete();
<br> vtkPoints *points = poly->GetPoints();<br> int NumberOfTotalPoints = points->GetNumberOfPoints();<br></p>
<p> unsigned int pointId = 0;<br> for ( int i = 0; i< NumberOfTotalPoints; i++)<br> {<br> double pos[3]; points->GetPoint(i,pos);</p>
<p> double x = 0;<br> double y = 0;<br> double z = 0;</p>
<p> if ( i < (int)NumberOfTotalPoints/2 ) {<br> x = 0.002*i;<br> y = 0.003*i;<br> z = 0.0001*i;<br> } else {<br> x = 0.002*(NumberOfTotalPoints-1-i);<br> y = 0.003*(NumberOfTotalPoints-1-i);<br> z = 0.0001*
(NumberOfTotalPoints-1-i);<br> }</p>
<p> std::cout << "generator 1 and 2 :: " << i << ": \t" << x << "\t" << y << "\t" << z << std::endl;</p>
<p> // Create source and target landmarks.<br> // <br> LandmarkPointType sourcePoint;<br> LandmarkPointType targetPoint;</p>
<p> sourcePoint[0] = pos[0] ;<br> sourcePoint[1] = pos[1] ;<br> sourcePoint[2] = pos[2] ;</p>
<p> targetPoint[0] = sourcePoint[0] + x ;<br> targetPoint[1] = sourcePoint[1] + y ;<br> targetPoint[2] = sourcePoint[2] + z ;</p>
<p> sourceLandmarks->InsertElement( pointId, sourcePoint );<br> targetLandmarks->InsertElement( pointId, targetPoint );<br> pointId++;</p>
<p> }</p>
<p> deformer->SetSourceLandmarks( sourceLandmarks.GetPointer() );<br> deformer->SetTargetLandmarks( targetLandmarks.GetPointer() );<br> <br> std::cout << " Up to this point we are still okay...." << std::endl;
</p>
<p> try<br> {<br> deformer->UpdateLargestPossibleRegion();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << "Exception thrown " << std::endl;<br> std::cerr << excp << std::endl;
<br> return EXIT_FAILURE;<br> }</p>
<p> std::cout << "the deformer is done at this point .... " << std::endl;</p>
<p> DeformationFieldType::ConstPointer deformationField = deformer->GetOutput();<br> typedef itk::WarpImageFilter< ImageType, ImageType, <br> DeformationFieldType > FilterType;</p>
<p> FilterType::Pointer warper = FilterType::New();<br> typedef itk::LinearInterpolateImageFunction< <br> ImageType, double > InterpolatorType;</p>
<p> InterpolatorType::Pointer interpolator = InterpolatorType::New();<br> warper->SetInterpolator( interpolator );<br> warper->SetOutputSpacing( deformationField->GetSpacing() );<br> warper->SetOutputOrigin( deformationField->GetOrigin() );
<br> warper->SetDeformationField( deformationField );<br> warper->SetInput( Reader->GetOutput() );<br> Writer->SetInput( warper->GetOutput() );</p>
<p><br> <br> try<br> {<br> Writer->Update();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << "Exception thrown " << std::endl;<br> std::cerr << excp << std::endl;
<br> return EXIT_FAILURE;<br> }</p>
<p><br> return EXIT_SUCCESS;</p>
<p>}</p>
<p>// ---------------- end --------------------</p></div>