Hi,<div>I'm trying to use GeodesicActiveContoursSegmentation with a 3D volume. I want to use a PET gradient image as the feature image and the CT image as the image to be segmented.</div><div>This is my code:</div><div>
<br></div><div><div>#include "itkImage.h"</div><div>#include "itkGeodesicActiveContourLevelSetImageFilter.h"</div><div><br></div><div>#include "itkCurvatureAnisotropicDiffusionImageFilter.h"</div>
<div>#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"</div><div>#include "itkSigmoidImageFilter.h"</div><div>#include "itkFastMarchingImageFilter.h"</div><div>#include "itkRescaleIntensityImageFilter.h"</div>
<div>#include "itkBinaryThresholdImageFilter.h"</div><div>#include "itkImageFileReader.h"</div><div>#include "itkImageFileWriter.h"</div><div><br></div><div><br></div><div>int main( int argc, char *argv[] )</div>
<div>{</div><div> if( argc < 7 )</div><div> {</div><div> std::cerr << "Missing Parameters " << std::endl;</div><div> std::cerr << "Usage: " << argv[0];</div><div> std::cerr << " inputImageCT inputImagePT outputImage";</div>
<div> std::cerr << " seedX seedY seedZ";</div><div> return 1;</div><div> }</div><div><br></div><div><br></div><div> typedef float InternalPixelType;</div><div> const unsigned int Dimension = 3;</div>
<div> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;</div><div><br></div><div> typedef unsigned char OutputPixelType;</div><div> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;</div>
<div> typedef itk::BinaryThresholdImageFilter< </div><div> InternalImageType, </div><div> OutputImageType > ThresholdingFilterType;</div><div> </div><div> ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();</div>
<div> </div><div> thresholder->SetLowerThreshold( -1000.0 );</div><div> thresholder->SetUpperThreshold( 0.0 );</div><div><br></div><div> thresholder->SetOutsideValue( 0 );</div><div>
thresholder->SetInsideValue( 255 );</div><div><br></div><div><br></div><div> typedef itk::ImageFileReader< InternalImageType > ReaderType;</div><div> typedef itk::ImageFileWriter< OutputImageType > WriterType;</div>
<div><br></div><div> ReaderType::Pointer reader = ReaderType::New();</div><div> ReaderType::Pointer readerPT = ReaderType::New();</div><div> WriterType::Pointer writer = WriterType::New();</div><div><br></div><div> reader->SetFileName( argv[1] );</div>
<div> readerPT->SetFileName( argv[2] );</div><div> writer->SetFileName( argv[3] );</div><div><br></div><div><br></div><div> typedef itk::RescaleIntensityImageFilter< </div><div> InternalImageType, </div>
<div> OutputImageType > CastFilterType;</div><div><br></div><div>///filter in order to obtain the feature image///</div><div> typedef itk::CurvatureAnisotropicDiffusionImageFilter< </div>
<div> InternalImageType, </div><div> InternalImageType > SmoothingFilterType;</div><div> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< </div>
<div> InternalImageType, </div><div> InternalImageType > GradientFilterType;</div><div> typedef itk::SigmoidImageFilter<</div><div> InternalImageType, </div>
<div> InternalImageType > SigmoidFilterType;</div><div><br></div><div> SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();</div><div> GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();</div>
<div> SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();</div><div><br></div><div> smoothing->SetTimeStep( 0.0625 );</div><div> smoothing->SetNumberOfIterations( 5 );</div><div> smoothing->SetConductanceParameter( 9.0 );</div>
<div><br></div><div> const double sigma = 1.0;</div><div> gradientMagnitude->SetSigma( sigma );</div><div><br></div><div> const double alpha = -0.5;</div><div> const double beta = 3.0;</div><div> sigmoid->SetAlpha( alpha );</div>
<div> sigmoid->SetBeta( beta );</div><div> sigmoid->SetOutputMinimum( 0.0 );</div><div> sigmoid->SetOutputMaximum( 1.0 );</div><div>///filter in order to obtain the feature image///</div><div><br></div><div>
<br></div><div>///level set///</div><div> typedef itk::FastMarchingImageFilter< </div><div> InternalImageType, </div><div> InternalImageType > FastMarchingFilterType; </div>
<div> typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, </div><div> InternalImageType > GeodesicActiveContourFilterType;</div><div><br></div><div> FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();</div>
<div> GeodesicActiveContourFilterType::Pointer geodesicActiveContour = </div><div> GeodesicActiveContourFilterType::New();</div><div><br></div><div> const double propagationScaling = 2.0;</div>
<div><br></div><div> geodesicActiveContour->SetPropagationScaling( propagationScaling );</div><div> geodesicActiveContour->SetCurvatureScaling( 1.0 );</div><div> geodesicActiveContour->SetAdvectionScaling( 1.0 );</div>
<div><br></div><div> geodesicActiveContour->SetMaximumRMSError( 0.02 );</div><div> geodesicActiveContour->SetNumberOfIterations( 800 );</div><div>///level set///</div><div><br></div><div><br></div><div>///utility for the fast marching///</div>
<div> typedef FastMarchingFilterType::NodeContainer NodeContainer;</div><div> typedef FastMarchingFilterType::NodeType NodeType;</div><div><br></div><div> NodeContainer::Pointer seeds = NodeContainer::New();</div>
<div><br></div><div> InternalImageType::IndexType seedPosition;</div><div> </div><div> seedPosition[0] = atoi( argv[4] );</div><div> seedPosition[1] = atoi( argv[5] );</div><div> seedPosition[1] = atoi( argv[6] );</div>
<div><br></div><div> const double initialDistance = 5.0;</div><div><br></div><div> NodeType node;</div><div><br></div><div> const double seedValue = - initialDistance;</div><div> </div><div> node.SetValue( seedValue );</div>
<div> node.SetIndex( seedPosition );</div><div><br></div><div> seeds->Initialize();</div><div> seeds->InsertElement( 0, node );</div><div><br></div><div> fastMarching->SetTrialPoints( seeds );</div><div><br>
</div><div> fastMarching->SetSpeedConstant( 1.0 );</div><div>///utility for the fast marching///</div><div><br></div><div><br></div><div> //<span class="Apple-tab-span" style="white-space:pre"> </span>1° step</div><div>
smoothing->SetInput( readerPT->GetOutput() );</div><div> gradientMagnitude->SetInput( smoothing->GetOutput() );</div><div> sigmoid->SetInput( gradientMagnitude->GetOutput() );</div><div><br></div><div>
//<span class="Apple-tab-span" style="white-space:pre"> </span>2° step</div><div> geodesicActiveContour->SetInput( fastMarching->GetOutput() );</div><div> geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );</div>
<div><br></div><div> //<span class="Apple-tab-span" style="white-space:pre"> </span>3° step</div><div> thresholder->SetInput( geodesicActiveContour->GetOutput() );</div><div> writer->SetInput( thresholder->GetOutput() );</div>
<div> </div><div><br></div><div> try</div><div> {</div><div> writer->Update();</div><div> }</div><div> catch( itk::ExceptionObject & excep )</div><div> {</div><div> std::cerr << "Exception caught !" << std::endl;</div>
<div> std::cerr << excep << std::endl;</div><div> }</div><div>}</div><div><br></div><div><br></div><div>unfortunatly, when i try to run it, i obtain this output message:</div><div><br></div><div><div>Exception caught !</div>
<div><br></div><div>itk::ExceptionObject (009EEB20)</div><div>Location: "__thiscall itk::ImageConstIterator<class itk::Image<class itk::FixedArray<float,3>,3> >:</div><div>ImageConstIterator(const class itk::Image<class itk::FixedArray<float,3>,3> *,const class itk::Imag</div>
<div>Region<3> &)"</div><div>File: c:\insighttoolkit-3.16.0\code\common\itkImageConstIterator.h</div><div>Line: 182</div><div>Description: itk::ERROR: Region ImageRegion (009EE8DC)</div><div> Dimension: 3</div>
<div> Index: [0, 0, 0]</div><div> Size: [128, 128, 287]</div><div> is outside of buffered region ImageRegion (024035E0)</div><div> Dimension: 3</div><div> Index: [0, 0, 0]</div><div> Size: [16, 16, 16]</div><div><br>
</div><div><br></div><div>i don't know why. Please help me, thanks, bye,</div><div><br></div><div>Marco</div></div></div>