<div class="gmail_quote"><div><div class="Wj3C7c">Hi,<br><div class="gmail_quote"><br>I've been trying to implement the DeformableMesh3DFilter and run into some problems. I'm basically using a variation of the DeformableModel1.cxx
algorithm in which I read in a vesselness image (which I get using itk's vesselness filter) and input his image to the gradient filter (I bypass the gradient magnitude calculation since this image is giving a maximum value at maximum vesselness measure). Instead of creating a mesh from a binary image, I read in points in 3D and create a polyline mesh from them (this is my initial mesh). I am using an open mesh as an initial guess that should be deformed to the vesselness image (which in my application corresponds to the coronary arteries). The program compiles fine but when I run it, it crashes when it attempts to do the deformable mesh filter. Is there another way to tackle this problem or solve the issue in this application? Anyway, below is a copy of my code.
<br><br>I really appreciate your help and time in advance.<br><br><br>Best Regards,<br clear="all"><br>-- <br>Federico J. Lopez Bertoni<br>Graduate Research Assistant<br>Duke University<br>Department of Biomedical Engineering
<br>Office: 1373 CIEMAS<br>Phone: 919-660-5125<br><a href="mailto:fjlb@duke.edu" target="_blank">fjlb@duke.edu</a><br><br>#if defined(_MSC_VER)<br>#pragma warning ( disable : 4786 )<br>#endif<br>#include "itkMesh.h
"<br>#include "
itkVertexCell.h"<br>#include "itkLineCell.h"<br>#include <cstdlib><br>#include <iostream><br>#include <fstream><br>#include <stdio.h><br>#include "itkImageFileReader.h"<br>
#include "itkImageFileWriter.h"<br>#include "itkCastImageFilter.h"<br>#include "itkRescaleIntensityImageFilter.h"<br>#include "itkGradientRecursiveGaussianImageFilter.h"<br>//#include "
itkGradientMagnitudeRecursiveGaussianImageFilter.h"<br>#include "itkDeformableMesh3DFilter.h"<br>#include "itkImage.h"<br>#include "itkCovariantVector.h"<br>#include "itkPointSetToImageFilter.h
"<br>using namespace std;<br><br><br><br>int main( int argc, char *argv[] )<br>{<br> if( argc < 4 )<br> {<br> std::cerr << "Missing Parameters " << std::endl;<br> std::cerr << "Usage: VesselnessImage RCApoints ImageOut" << std::endl;
<br> return EXIT_FAILURE;<br> }<br><br> //Read RCA, LCA points in.<br> ifstream infile;<br> infile.open( argv[2] , ios::in );<br> if (!infile.good())<br> {<br> cout << "Error Opening file " << endl;
<br> exit (1);<br> }<br> int X[200], Y[200], Z[200], rad[200], i=0;<br> while(!infile.eof())<br> {<br> infile >> X[i] >> Y[i] >> Z[i] >> rad[i];<br> if (X[i] == -1 )
<br> break;<br> else<br> //cout << i << "\t" << X[i] << "\t" << Y[i] << "\t" << Z[i] << "\t" << rad[i] << endl;
<br> i++;<br> }<br><br> /*ifstream infile1;<br> infile1.open( argv[3] , ios::in );<br> if (!infile1.good())<br> {<br> cout << "Error Opening file "<< endl;<br> exit (1);
<br> }<br> int X1[200], Y1[200], Z1[200], rad1[200], i1=0;<br> while(!infile1.eof())<br> {<br> infile1 >> X1[i1] >> Y1[i1] >> Z1[i1] >> rad1[i1];<br> if (X1[i1] == -1 )<br>
break;<br> else<br> cout << i1 << "\t" << X1[i1] << "\t" << Y1[i1] << "\t" << Z1[i1] << "\t" << rad1[i1] << endl;
<br> i1++;<br> }*/<br><br> //typedef unsigned short PixelType;<br> typedef double PixelType;<br> const unsigned int Dimension = 3;<br> typedef itk::Mesh< PixelType, Dimension > MeshType;
<br> typedef MeshType::CellType CellType;<br> typedef itk::VertexCell< CellType > VertexType;<br> typedef itk::LineCell< CellType > LineType;<br> typedef itk::Image< PixelType, Dimension > ImageType;
<br> <br> MeshType::Pointer mesh = MeshType::New();<br> MeshType::PointType points;<br> //Create the mesh<br> for(unsigned int j=0; j<92; j++)<br> {<br> points[0] = X[j];<br> points[1] = Y[j];
<br> points[2] = Z[j];<br> mesh->SetPoint( j, points );<br> }<br> /*for(unsigned int j=0; j<92; j++)<br> {<br> points[0][0] = X[j];<br> points[0][1] = Y[j];<br> points[0][2] = Z[j];
<br> mesh->SetPoint( j, points );<br> }*/<br><br> CellType::CellAutoPointer cellpointer;<br> for(unsigned int j=0; j<91; j++)<br> {<br> cellpointer.TakeOwnership( new LineType );<br> cellpointer->SetPointId( 0, j );
<br> cellpointer->SetPointId( 0, j+1 );<br> mesh->SetCell( j, cellpointer );<br><br> cellpointer.TakeOwnership( new VertexType );<br> cellpointer->SetPointId( 0, j );<br> mesh->SetCell( j + 90, cellpointer );
<br> }<br><br> // Print out the number of points and the number of cells.<br> std::cout << "# Points= " << mesh->GetNumberOfPoints() << std::endl;<br> std::cout << "# Cell = " << mesh->GetNumberOfCells() << std::endl;
<br><br> typedef MeshType::PointsContainer::ConstIterator PointIterator;<br> PointIterator pointIterator = mesh->GetPoints()->Begin();<br> PointIterator pointEnd = mesh->GetPoints()->End();<br>
<br> while( pointIterator != pointEnd ) <br> {<br> //std::cout << pointIterator.Value() << std::endl;<br> ++pointIterator;<br> }<br><br> typedef MeshType::CellsContainer::ConstIterator CellIterator;
<br><br> CellIterator cellIterator = mesh->GetCells()->Begin();<br> CellIterator cellEnd = mesh->GetCells()->End();<br> <br> while( cellIterator != cellEnd ) <br> {<br> CellType * cell =
cellIterator.Value();<br> //std::cout << cell->GetNumberOfPoints() << std::endl;<br> ++cellIterator;<br> }<br><br> cellIterator = mesh->GetCells()->Begin();<br> cellEnd = mesh->GetCells()->End();
<br> <br> while( cellIterator != cellEnd ) <br> {<br> CellType * cell = cellIterator.Value();<br> //std::cout << "cell with " << cell->GetNumberOfPoints();<br> //std::cout << " points " << std::endl;
<br><br> // Software Guide : BeginCodeSnippet<br> typedef CellType::PointIdIterator PointIdIterator;<br><br> PointIdIterator pointIditer = cell->PointIdsBegin();<br> PointIdIterator pointIdend = cell->PointIdsEnd();
<br><br> while( pointIditer != pointIdend )<br> {<br> //std::cout << *pointIditer << std::endl;<br> ++pointIditer;<br> }<br> // Software Guide : EndCodeSnippet<br><br> ++cellIterator;
<br> }<br><br> //Mesh creation done at this point. <br> std::cout << "PolyLine mesh created..." << std::endl;<br> std::cout << "Read in Input Image..." << std::endl;
<br> //Now read in Vesselness image.<br> //Read in Vesselness (Cost) Image<br> typedef itk::ImageFileReader< ImageType > ReaderType;<br> ReaderType::Pointer reader = ReaderType::New();<br> reader->SetFileName( argv[1] );
<br> //Set image to variable<br> ImageType::Pointer image = ImageType::New();<br> image = reader->GetOutput();<br> //Define Gradient Filter Type<br> typedef itk::CovariantVector< double, Dimension > GradientPixelType;
<br> typedef itk::Image< GradientPixelType, Dimension > GradientImageType;<br> typedef itk::GradientRecursiveGaussianImageFilter< ImageType, GradientImageType > GradientFilterType;<br> //typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< ImageType,ImageType > GradientMagnitudeFilterType;
<br> GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New(); <br> //GradientMagnitudeFilterType::Pointer gradientMapFilter = GradientMagnitudeFilterType::New();<br> //Set input of gradient filter to be vesselnes image
<br> gradientMapFilter->SetInput( reader->GetOutput() );<br> gradientMapFilter->SetSigma( 1.0 );<br> //trigger execution<br> try<br> {<br> gradientMapFilter->Update();<br> }<br> catch( itk::ExceptionObject & e )
<br> {<br> std::cerr << "Exception caught when updating gradientMapFilter " << std::endl;<br> std::cerr << e << std::endl;<br> return -1;<br> }<br><br> std::cout << "The gradient Map was created successfully" << std::endl;
<br><br> // Mesh Source filter<br> //MeshSourceType::Pointer meshSource = MeshSourceType::New();<br><br> //Deformable model filter<br> typedef itk::DeformableMesh3DFilter< MeshType, MeshType > DeformableFilterType;
<br> DeformableFilterType::Pointer deformableModelFilter = DeformableFilterType::New();<br> //Set input to deformable model filter<br> deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );
<br>
deformableModelFilter->SetInput( mesh ); <br> // Set parameters for deformable mesh<br> typedef itk::CovariantVector< double, 2 > double2DVector;<br> typedef itk::CovariantVector< double, 3 > double3DVector;
<br><br> double2DVector stiffness;<br> stiffness[0] = 0.01;<br> stiffness[1] = 0.1;<br><br> double3DVector scale;<br> scale[0] = 1.0;<br> scale[1] = 1.0;<br> scale[2] = 1.0;<br><br> deformableModelFilter->SetStiffness( stiffness );
<br> deformableModelFilter->SetScale( scale );<br> deformableModelFilter->SetGradientMagnitude( 0.8 );<br> deformableModelFilter->SetTimeStep( 0.01 );<br> deformableModelFilter->SetStepThreshold( 60 );
<br><br> std::cout << "Deformable mesh fitting..." << std::endl;<br><br> //Trigger execution of deformable mesh filter<br> try<br> {<br> deformableModelFilter->Update();<br> }
<br> catch( itk::ExceptionObject & excep )<br> {<br> std::cerr << "Exception caught while using the deformable mesh filter" << std::endl;<br> std::cerr << excep << std::endl;
<br> }<br><br> std::cout << "Mesh Source: " << std::endl;<br><br> typedef itk::PointSetToImageFilter< MeshType, ImageType > MeshFilterType;<br> MeshFilterType::Pointer meshFilter = MeshFilterType::New();
<br> meshFilter->SetOrigin( image->GetOrigin() );<br> meshFilter->SetSize( image->GetLargestPossibleRegion().GetSize() );<br> meshFilter->SetSpacing( image->GetSpacing() );<br> meshFilter->SetInput( mesh );
<br><br> try<br> {<br> meshFilter->Update();<br> }<br> catch( itk::ExceptionObject & excep )<br> {<br> std::cerr << "Exception Caught!" << std::endl;<br> std::cerr << excep << std::endl;
<br> }<br><br> typedef itk::ImageFileWriter< ImageType > WriterType;<br> WriterType::Pointer writer = WriterType::New();<br> writer->SetInput( meshFilter->GetOutput() );<br> writer->SetFileName( argv[3] );
<br> writer->Update();<br><br> return EXIT_SUCCESS; <br>}<br><br>
</div><br><br clear="all"><br>-- <br>Federico J. Lopez Bertoni<br>Graduate Research Assistant<br>Duke University<br>Department of Biomedical Engineering<br>Office: 1373 CIEMAS<br>Phone: 919-660-5125<br><a href="mailto:fjlb@duke.edu" target="_blank">
fjlb@duke.edu</a>
</div></div></div><br><br clear="all"><br>-- <br>Federico J. Lopez Bertoni<br>Graduate Research Assistant<br>Duke University<br>Department of Biomedical Engineering<br>Office: 1373 CIEMAS<br>Phone: 919-660-5125<br><a href="mailto:fjlb@duke.edu">
fjlb@duke.edu</a>