[Insight-users] FastMarching-output is empty using GeodesicActiveContourLevelSetImageFilter for 3D data sets

[Oliver Gloger]

Dear Insight-users,


I want to apply the "itkGeodesicActiveContourLevelSetImageFilter" for 3D datasets. For 2D images I had success in applying this filter!
However, for 3D datasets (I tried .gipl and .vtk-formats by using the corresponding IO-classes) I have problems using this filter!

I store the outputs from the "itkGradientMagnitudeImageFilter" and "itkSigmoidImageFilter" and they look fine! Just how I expected, they show the 3D edges (and inverse edges after applying the sigmoid filter) from the 3D data set.

However, the Output of the "itkFastMarchingImageFilter" is empty <=> contrarily to the 2D-case. In the 2D-case there is always a signed distance image!
In 3D there is a gipl-OutputImage, but there is nothing to see!!! (I provide the fastMarchingfilter 1 seepoint!!!)

Consequently, the final binary result of the "itkGeodesicActiveContourLevelSetImageFilter" is also empty! Does the "itkGeodesicActiveContourLevelSetImageFilter" or the "itkFastMarchingImageFilter" not function for 3D datasets? Or could anybody of you give me a hint what I made wrong?I send the code below!

I have also another theoretical question: in SETHIAN's book the fast Marching approach is used for the boundary value problem by solving the Eikonal equation. However, this approach can only expand or shrink! Hence, the geodesic active contours could also move in both directions at certain locations due to curvature or edge attraction force! So, why is it justified to use the fast Marching approach in combination with the geodesic active contours?

Thank you very much for your help!

Olli

---

#include "itkImage.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkGDCMImageIO.h"
#include "itkGiplImageIO.h"
#include "itkVTKImageIO.h"
#include "DicomSeriesReadImageWrite2.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"
#include "DicomSeriesReadImageWrite2.h"
using namespace std;

//reading DicomSeries and storing volume data in gipl-file
void DicomSeriesReadImageWrite2(char* DicomDirectory, char* outputFileName, char* seriesName)
//This function is okay => I do not send it here!

int main( int argc, char *argv[] )
{
if( argc < 11 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " seedX seedY seedZ InitialDistance";
std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
std::cerr << " PropagationScaling CurvatureScaling AdvectionScaling"<< std::endl;
return 1;
}

//char *path = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\";
// String path = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\";

char *src_dicom = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated";
char *dst_3D = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\dreiD.gipl";
char *serName ="";

DicomSeriesReadImageWrite2(src_dicom,dst_3D,serName);

typedef float InternalPixelType;
const unsigned int Dimension = 3; //here all should be in 3D !!!
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;

typedef itk::BinaryThresholdImageFilter< 
InternalImageType, 
OutputImageType > ThresholdingFilterType;

ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();

thresholder->SetLowerThreshold( -1000.0 );
thresholder->SetUpperThreshold( 0.0 );

thresholder->SetOutsideValue( 0 );
thresholder->SetInsideValue( 255 );

typedef itk::ImageFileReader< InternalImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;

ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();


typedef itk::GDCMImageIO ImageIOType;
ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
typedef itk::VTKImageIO VTKIOType;
VTKIOType::Pointer vtkIO = VTKIOType::New();
typedef itk::GiplImageIO GiplIOType;
GiplIOType::Pointer giplIO = GiplIOType::New(); 

reader->SetImageIO(giplIO);
reader->SetFileName(dst_3D);

char *name = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\volume_dataset.gipl";
writer->SetFileName(name);
writer->SetImageIO(giplIO);
typedef itk::RescaleIntensityImageFilter< 
InternalImageType, 
OutputImageType > CastFilterType;

typedef itk::CurvatureAnisotropicDiffusionImageFilter< 
InternalImageType, 
InternalImageType > SmoothingFilterType;

SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();

typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< 
InternalImageType, 
InternalImageType > GradientFilterType;

typedef itk::SigmoidImageFilter< 
InternalImageType, 
InternalImageType > SigmoidFilterType;

GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();

SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

sigmoid->SetOutputMinimum( 0.0 );
sigmoid->SetOutputMaximum( 1.0 );

typedef itk::FastMarchingImageFilter< 
InternalImageType, 
InternalImageType > FastMarchingFilterType;


FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();

typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, 
InternalImageType > GeodesicActiveContourFilterType;
GeodesicActiveContourFilterType::Pointer geodesicActiveContour = 
GeodesicActiveContourFilterType::New();

const double propagationScaling = atof( argv[8] );

geodesicActiveContour->SetPropagationScaling( propagationScaling );
geodesicActiveContour->SetCurvatureScaling( atof(argv[9]) );
geodesicActiveContour->SetAdvectionScaling( atof(argv[10]) );
geodesicActiveContour->SetMaximumRMSError( 0.02 );
geodesicActiveContour->SetNumberOfIterations( 5000 );

smoothing->SetInput( reader->GetOutput() );
gradientMagnitude->SetInput( smoothing->GetOutput() );
sigmoid->SetInput( gradientMagnitude->GetOutput() );

geodesicActiveContour->SetInput( fastMarching->GetOutput() );
geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );

thresholder->SetInput( geodesicActiveContour->GetOutput() );
writer->SetInput( thresholder->GetOutput() );

smoothing->SetTimeStep( 0.0625 );
smoothing->SetNumberOfIterations( 10 );
smoothing->SetConductanceParameter( 1.0 );

const double sigma = atof( argv[5] );
gradientMagnitude->SetSigma( sigma );

const double alpha = atof( argv[6] );
const double beta = atof( argv[7] );

sigmoid->SetAlpha( alpha );
sigmoid->SetBeta( beta );

typedef FastMarchingFilterType::NodeContainer NodeContainer;
typedef FastMarchingFilterType::NodeType NodeType;

NodeContainer::Pointer seeds = NodeContainer::New();

InternalImageType::IndexType seedPosition;

seedPosition[0] = atoi( argv[1] );
seedPosition[1] = atoi( argv[2] );
seedPosition[2] = atoi( argv[3] );

const double initialDistance = atof( argv[4] );

NodeType node;

const double seedValue = - initialDistance;

node.SetValue( seedValue );
node.SetIndex( seedPosition );

seeds->Initialize();
seeds->InsertElement( 0, node );

fastMarching->SetTrialPoints( seeds );
fastMarching->SetSpeedConstant( 1.0 );
fastMarching->DebugOn();

CastFilterType::Pointer caster1 = CastFilterType::New();
CastFilterType::Pointer caster2 = CastFilterType::New();
CastFilterType::Pointer caster3 = CastFilterType::New();
CastFilterType::Pointer caster4 = CastFilterType::New();

WriterType::Pointer writer1 = WriterType::New();
WriterType::Pointer writer2 = WriterType::New();
WriterType::Pointer writer3 = WriterType::New();
WriterType::Pointer writer4 = WriterType::New();

caster1->SetInput( smoothing->GetOutput() );
writer1->SetInput( caster1->GetOutput() );
char *name1 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\smoothed_dataset.gipl";
writer1->SetFileName(name1);
caster1->SetOutputMinimum( 0 );
caster1->SetOutputMaximum( 255 );
writer1->SetImageIO(giplIO);
try{
writer1->Update();
}
catch( itk::ExceptionObject & excep ){
std::cerr << "Exception caught at writer1!" << std::endl;
std::cerr << excep << std::endl;
}

caster2->SetInput( gradientMagnitude->GetOutput() );
writer2->SetInput( caster2->GetOutput() );
char *name2 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\3DKanten.gipl";
writer2->SetFileName(name2);
caster2->SetOutputMinimum( 0 );
caster2->SetOutputMaximum( 255 );
writer2->SetImageIO(giplIO);
try{
writer2->Update();
}
catch( itk::ExceptionObject & excep ){
std::cerr << "Exception caught at writer2!" << std::endl;
std::cerr << excep << std::endl;
}

caster3->SetInput( sigmoid->GetOutput() );
writer3->SetInput( caster3->GetOutput() );
char *name3 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\inverse3DKanten.gipl";
writer3->SetFileName(name3);
caster3->SetOutputMinimum( 0 );
caster3->SetOutputMaximum( 255 );
writer3->SetImageIO(giplIO);
try{
writer3->Update();
}
catch( itk::ExceptionObject & excep ){
std::cerr << "Exception caught at writer3 !" << std::endl;
std::cerr << excep << std::endl;
}

caster4->SetInput( fastMarching->GetOutput() );
writer4->SetInput( caster4->GetOutput() );
char *name4 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\FFM.gipl";
writer4->SetFileName(name4);
caster4->SetOutputMinimum( 0 );
caster4->SetOutputMaximum( 255 );
writer4->SetImageIO(giplIO);
try{
writer4->Update();
std::cout<< "FMIF: GetOutputSize: "<<fastMarching->GetOutputSize()<<std::endl;
}
catch( itk::ExceptionObject & excep ){
std::cerr << "Exception caught at writer4!" << std::endl;
std::cerr << excep << std::endl;
}

fastMarching->SetOutputSize( 
reader->GetOutput()->GetBufferedRegion().GetSize() );
try
{
writer->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception caught at main writer!" << std::endl;
std::cerr << excep << std::endl;
}
std::cout << std::endl;
std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
std::cout << std::endl;
std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;

return 0;
}


I call the filter with those parameters (for instance, but I also changed the values in a small range!!!)
sigma=1, sigmiodAlpha = -4, SigmoidBeta=0, PrpagationScaling=3, CurvatureScaling=2, AdvectionScaling=3



______________________________________________________
GRATIS für alle WEB.DE-Nutzer: Die maxdome Movie-FLAT!
Jetzt freischalten unter http://movieflat.web.de