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

[Ming Chao]

Hi Luis,
Thanks a lot for looking into the problem. You are right that the problem
was related to the numbr (7426, too large) of landmark points I used to
deform the image. Anyhow, you could find the image I wanted to deform and
the polydata input file in the following locations:
http://www.stanford.edu/~mingchao/image/phase1.vtk (image)
and
http://www.stanford.edu/~mingchao/image/LTLung.vtk (polydata).
Is there an easy way to reduce the number of landmark points in my polydata
input files?

Cheers,
Ming

On 3/20/06, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
>
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20060320/e0532546/attachment.html