[Insight-users] problem of deforming a 3D image

[Luis Ibanez]

Hi Ming,

Thanks for your detailed report and for posting your source code.

Your code seems to be very similar to the example:


       Insight/
          Examples/
              Registration/
                  DeformationFieldInitialization.cxx


The problem seems to be related to the number of landmark points
that you are passing for initializing the deformation field.

Could you please let us know how many points you are using ?

Could you send us the input file with the vtkPolyData that you
are using for feeding this program ?
That data could help us to figure out what may be going wrong
with your code.


   Please let us know,


     Thanks


       Luis


--------------
Ming Chao wrote:
> Hi,
>  
> 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:
> 
> //: Sets all elements of matrix to specified value. O(m*n).
> 
> template<class T>
> void vnl_matrix<T>::fill (T const& value)
> {
>   for (unsigned int i = 0; i < this->num_rows; i++)
>     for (unsigned int j = 0; j < this->num_cols; j++)
>       this->data[i][j] = value;   
> <========================================= debugging reports problem !
> }
> 
>  
>  
>  
> Also, I am attaching my code of deforming 3D images below. I would 
> greatly appreciate it if anybody can help look into it.
>  
> Best regards,
> Ming
>  
> 
> // ---------------- starting of the code  -------------------------
> #include "itkVector.h"
> #include "itkImage.h"
> #include "itkImageRegionIteratorWithIndex.h"
> #include "itkDeformationFieldSource.h "
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkWarpImageFilter.h"
> #include "itkLinearInterpolateImageFunction.h"
> 
> #include "itkNormalVariateGenerator.h"
> #include "vtkPolyDataReader.h"
> #include "vtkPolyDataWriter.h"
> #include "vtkFloatArray.h"
> #include "vtkPolyData.h"
> #include "vtkPointData.h"
> #include "vtkPoints.h"
> 
> #include <fstream>
> 
> int main( int argc, char * argv[] )
> {
> 
>  if( argc < 2 )
>  {
>   std::cerr << "Missing Parameters " << std::endl;
>   std::cerr << "Usage: " << argv[0];
>   std::cerr << " landmarksFile fixedImage ";
>   std::cerr << "movingImage deformedMovingImage" << std::endl;
>   return 1;
>     }
> 
>  const     unsigned int   Dimension = 3;
>  typedef   float          VectorComponentType;
> 
>  typedef   itk::Vector< VectorComponentType, Dimension >    VectorType;
>  typedef   itk::Image< VectorType,  Dimension >   DeformationFieldType;
>  typedef   float  PixelType;
>  typedef   itk::Image< PixelType, Dimension >  ImageType  ;
>  typedef   itk::ImageFileReader< ImageType  >  ReaderType ;
>  typedef   itk::ImageFileWriter< ImageType  >  WriterType ;
> 
>  ReaderType::Pointer Reader = ReaderType::New();
>  Reader->SetFileName( argv[1] );
> 
>  try
>     {
>   Reader->Update();
>     }
>  catch( itk::ExceptionObject & excp )
>     {
>   std::cerr << "Exception thrown " << std::endl;
>   std::cerr << excp << std::endl;
>   return EXIT_FAILURE;
>     }
> 
>  WriterType::Pointer Writer = WriterType::New();
>  Writer->SetFileName( argv[2] );
> 
>  ImageType::ConstPointer fixedImage = Reader->GetOutput();
> 
> 
>  typedef itk::DeformationFieldSource< DeformationFieldType >  
> DeformationSourceType;
>  DeformationSourceType::Pointer deformer = DeformationSourceType::New();
>  deformer->SetOutputSpacing( fixedImage->GetSpacing() );
>  deformer->SetOutputOrigin(  fixedImage->GetOrigin() );
>  deformer->SetOutputRegion(  fixedImage->GetLargestPossibleRegion() );
> 
>  typedef DeformationSourceType::LandmarkContainer          
> LandmarkContainerType;
>  typedef DeformationSourceType::LandmarkPointType          
> LandmarkPointType;
>  LandmarkContainerType::Pointer sourceLandmarks = 
> LandmarkContainerType::New();
>  LandmarkContainerType::Pointer targetLandmarks = 
> LandmarkContainerType::New();
> 
>  vtkPolyDataReader *reader = vtkPolyDataReader::New();
>  reader->SetFileName( argv[2] );
>  reader->Update();
>  vtkPolyData *poly = vtkPolyData::New();
>  poly->DeepCopy( reader->GetOutput());
>  reader->Delete();
>  vtkPoints *points = poly->GetPoints();
>  int NumberOfTotalPoints = points->GetNumberOfPoints();
> 
>  unsigned int pointId = 0;
>  for ( int i = 0; i< NumberOfTotalPoints; i++)
>  {
>   double pos[3]; points->GetPoint(i,pos);
> 
>   double x = 0;
>   double y = 0;
>   double z = 0;
> 
>   if ( i < (int)NumberOfTotalPoints/2 ) {
>    x = 0.002*i;
>    y = 0.003*i;
>    z = 0.0001*i;
>   } else {
>    x = 0.002*(NumberOfTotalPoints-1-i);
>    y = 0.003*(NumberOfTotalPoints-1-i);
>    z = 0.0001* (NumberOfTotalPoints-1-i);
>   }
> 
>   std::cout << "generator 1 and 2 :: " << i << ": \t" << x << "\t" << y 
> << "\t" << z << std::endl;
> 
>   //  Create source and target landmarks.
>   // 
>   LandmarkPointType sourcePoint;
>   LandmarkPointType targetPoint;
> 
>   sourcePoint[0] = pos[0] ;
>   sourcePoint[1] = pos[1] ;
>   sourcePoint[2] = pos[2] ;
> 
>   targetPoint[0] = sourcePoint[0] + x ;
>   targetPoint[1] = sourcePoint[1] + y ;
>    targetPoint[2] = sourcePoint[2] + z ;
> 
>   sourceLandmarks->InsertElement( pointId, sourcePoint );
>   targetLandmarks->InsertElement( pointId, targetPoint );
>   pointId++;
> 
>  }
> 
>  deformer->SetSourceLandmarks( sourceLandmarks.GetPointer() );
>  deformer->SetTargetLandmarks( targetLandmarks.GetPointer() );
>  
>  std::cout << " Up to this point we are still okay...." << std::endl;
> 
>  try
>     {
>   deformer->UpdateLargestPossibleRegion();
>     }
>  catch( itk::ExceptionObject & excp )
>     {
>   std::cerr << "Exception thrown " << std::endl;
>   std::cerr << excp << std::endl;
>   return EXIT_FAILURE;
>     }
> 
>  std::cout << "the deformer is done at this point .... " << std::endl;
> 
>  DeformationFieldType::ConstPointer deformationField = 
> deformer->GetOutput();
>  typedef itk::WarpImageFilter< ImageType, ImageType,
>                          DeformationFieldType  >  FilterType;
> 
>  FilterType::Pointer warper = FilterType::New();
>  typedef itk::LinearInterpolateImageFunction<
>                        ImageType, double >  InterpolatorType;
> 
>  InterpolatorType::Pointer interpolator = InterpolatorType::New();
>  warper->SetInterpolator( interpolator );
>  warper->SetOutputSpacing( deformationField->GetSpacing() );
>  warper->SetOutputOrigin(  deformationField->GetOrigin() );
>  warper->SetDeformationField( deformationField );
>  warper->SetInput( Reader->GetOutput() );
>  Writer->SetInput( warper->GetOutput() );
> 
> 
>  
>  try
>     {
>   Writer->Update();
>     }
>  catch( itk::ExceptionObject & excp )
>     {
>   std::cerr << "Exception thrown " << std::endl;
>   std::cerr << excp << std::endl;
>   return EXIT_FAILURE;
>     }
> 
> 
>  return EXIT_SUCCESS;
> 
> }
> 
> // ---------------- end --------------------
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users