<div>Hi Luis,</div>
<div>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:
</div>
<div><a href="http://www.stanford.edu/~mingchao/image/phase1.vtk">http://www.stanford.edu/~mingchao/image/phase1.vtk</a> (image)</div>
<div>and</div>
<div><a href="http://www.stanford.edu/~mingchao/image/LTLung.vtk">http://www.stanford.edu/~mingchao/image/LTLung.vtk</a> (polydata).</div>
<div>Is there an easy way to reduce the number of landmark points in my polydata input files? <br> </div>
<div>Cheers,</div>
<div>Ming<br> </div>
<div><span class="gmail_quote">On 3/20/06, <b class="gmail_sendername">Luis Ibanez</b> <<a href="mailto:luis.ibanez@kitware.com">luis.ibanez@kitware.com</a>> wrote:</span>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid"><br>Hi Ming,<br><br>Thanks for your detailed report and for posting your source code.<br><br>Your code seems to be very similar to the example:
<br><br><br> Insight/<br> Examples/<br> Registration/<br> DeformationFieldInitialization.cxx<br><br><br>The problem seems to be related to the number of landmark points<br>that you are passing for initializing the deformation field.
<br><br>Could you please let us know how many points you are using ?<br><br>Could you send us the input file with the vtkPolyData that you<br>are using for feeding this program ?<br>That data could help us to figure out what may be going wrong
<br>with your code.<br><br><br> Please let us know,<br><br><br> Thanks<br><br><br> Luis<br><br><br>--------------<br>Ming Chao wrote:<br>> Hi,<br>><br>> I encountered a problem when trying to deform a 3D image. The way I
<br>> deformed the image is as the following: I created a 3D contour for an<br>> organ, say, lung for this image. Along the contours I have a bunch of<br>> points. Then I generated a vector field for these points and used these
<br>> vectors to deform the image. I didn't have trouble deforming any 2D<br>> image. But for 3D image I got an error. I traced down to the point where<br>> the problem occured: in the try-catch block with<br>> deformer->UpdateLargestPossibleRegion(); debugging reports a problem on
<br>> line 417 of vnl_matrix.txx file. For reference purpose I copy and paste<br>> that part of code here:<br>><br>> //: Sets all elements of matrix to specified value. O(m*n).<br>><br>> 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;
<br>> <========================================= debugging reports problem !<br>> }<br>><br>><br>><br>><br>> Also, I am attaching my code of deforming 3D images below. I would<br>> greatly appreciate it if anybody can help look into it.
<br>><br>> Best regards,<br>> Ming<br>><br>><br>> // ---------------- 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"<br>><br>> #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"<br>><br>> #include <fstream>
<br>><br>> int main( int argc, char * argv[] )<br>> {<br>><br>> 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>> }<br>><br>> const unsigned int Dimension = 3;
<br>> typedef float VectorComponentType;<br>><br>> 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 ;
<br>><br>> ReaderType::Pointer Reader = ReaderType::New();<br>> Reader->SetFileName( argv[1] );<br>><br>> 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>> }<br>><br>> WriterType::Pointer Writer = WriterType::New();
<br>> Writer->SetFileName( argv[2] );<br>><br>> ImageType::ConstPointer fixedImage = Reader->GetOutput();<br>><br>><br>> typedef itk::DeformationFieldSource< DeformationFieldType ><br>> DeformationSourceType;
<br>> DeformationSourceType::Pointer deformer = DeformationSourceType::New();<br>> deformer->SetOutputSpacing( fixedImage->GetSpacing() );<br>> deformer->SetOutputOrigin( fixedImage->GetOrigin() );
<br>> deformer->SetOutputRegion( fixedImage->GetLargestPossibleRegion() );<br>><br>> typedef DeformationSourceType::LandmarkContainer<br>> LandmarkContainerType;<br>> typedef DeformationSourceType::LandmarkPointType
<br>> LandmarkPointType;<br>> LandmarkContainerType::Pointer sourceLandmarks =<br>> LandmarkContainerType::New();<br>> LandmarkContainerType::Pointer targetLandmarks =<br>> LandmarkContainerType::New();<br>
><br>> 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>><br>> unsigned int pointId = 0;<br>> for ( int i = 0; i< NumberOfTotalPoints; i++)
<br>> {<br>> double pos[3]; points->GetPoint(i,pos);<br>><br>> double x = 0;<br>> double y = 0;<br>> double z = 0;<br>><br>> 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>> }<br>><br>> std::cout << "generator 1 and 2 :: " << i << ": \t" << x << "\t" << y<br>> << "\t" << z << std::endl;
<br>><br>> // Create source and target landmarks.<br>> //<br>> LandmarkPointType sourcePoint;<br>> LandmarkPointType targetPoint;<br>><br>> sourcePoint[0] = pos[0] ;<br>> sourcePoint[1] = pos[1] ;
<br>> sourcePoint[2] = pos[2] ;<br>><br>> targetPoint[0] = sourcePoint[0] + x ;<br>> targetPoint[1] = sourcePoint[1] + y ;<br>> targetPoint[2] = sourcePoint[2] + z ;<br>><br>> sourceLandmarks->InsertElement( pointId, sourcePoint );
<br>> targetLandmarks->InsertElement( pointId, targetPoint );<br>> pointId++;<br>><br>> }<br>><br>> deformer->SetSourceLandmarks( sourceLandmarks.GetPointer() );<br>> deformer->SetTargetLandmarks(
targetLandmarks.GetPointer() );<br>><br>> std::cout << " Up to this point we are still okay...." << std::endl;<br>><br>> 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>> }<br>><br>> std::cout << "the deformer is done at this point .... " << std::endl;<br>><br>> DeformationFieldType::ConstPointer deformationField =<br>> deformer->GetOutput();
<br>> typedef itk::WarpImageFilter< ImageType, ImageType,<br>> DeformationFieldType > FilterType;<br>><br>> FilterType::Pointer warper = FilterType::New();<br>> typedef itk::LinearInterpolateImageFunction<
<br>> ImageType, double > InterpolatorType;<br>><br>> 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() );
<br>><br>><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>> }<br>><br>><br>> return EXIT_SUCCESS;<br>><br>> }<br>><br>> // ---------------- end --------------------
<br>><br>><br>> ------------------------------------------------------------------------<br>><br>> _______________________________________________<br>> Insight-users mailing list<br>> <a href="mailto:Insight-users@itk.org">
Insight-users@itk.org</a><br>> <a href="http://www.itk.org/mailman/listinfo/insight-users">http://www.itk.org/mailman/listinfo/insight-users</a><br><br></blockquote></div><br>