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

[Ming Chao]

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 --------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20060317/6654d3f8/attachment.html