[Insight-users] problem of deforming a 3D image
Luis Ibanez
luis.ibanez at kitware.com
Mon Mar 20 17:47:11 EST 2006
Hi Ming,
When you use Kernel Splines (which is what the DeformationFieldSource
class does), you don't need many landmarks, unless you have a lot of
local deformations.
In general something between 20 to 100 well located points should be
enough for driving the deformation of a 3D image, specially for the
case of biological structures, since they don't tend to have many
local deformations.
In order to reduce the number of points in your vtkPolyData you could
use the decimate filters available in VTK.
For example:
http://www.vtk.org/doc/nightly/html/classvtkDecimatePro.html
Regards,
Luis
==================
Ming Chao wrote:
> 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
> <mailto: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 <mailto:Insight-users at itk.org>
> > http://www.itk.org/mailman/listinfo/insight-users
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list