[ITK-users] Chan and Vese Segmentation for 3D data
Param Rajpura
Param.Rajpura at LntTechservices.com
Tue Feb 10 07:27:21 EST 2015
Hi all,
I am using ITK for segmentation from 3D MR Images data.
Initially applying ChanandVeseFilter on 2D images I could segment successfully but 3D data generates a problem.
I am using the same code provided in the example: "Single Phase Chan and Vese" the only difference is I am loading fast marching output instead of generating here.
I am unable to catch any exceptions here, but while debugging my code, first chance exception is thrown at:
1. Inside levelSetFilter->Update();
2. In ImportImageContainer.hxx
3. Line 181, data = new TElement[size];
Below is my code.
typedef itk::ImageFileReader< MainType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName("E:\\Intermediate Results\\Sample001_MRI.mhd");
try
{
reader->Update();
}
catch (itk::ExceptionObject & excep)
{
std::cerr << "Exception caught in read !" << std::endl;
std::cerr << excep << std::endl;
return -1;
}
typedef itk::CastImageFilter<MainType, InternalImageType> FloatConversionFilter;
FloatConversionFilter::Pointer ImgTypeConverter = FloatConversionFilter::New();
ImgTypeConverter->SetInput(reader->GetOutput());
ImgTypeConverter->Update();
typedef itk::ImageFileWriter< InternalImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("E:\\Intermediate Results\\Out.mhd");
typedef itk::ImageFileReader< InternalImageType > ReaderType1;
ReaderType1::Pointer reader1 = ReaderType1::New();
reader1->SetFileName("E:\\Intermediate Results\\FastMarching.mhd");
reader1->Update();
InternalImageType::Pointer Image = InternalImageType::New();
Image = ImgTypeConverter->GetOutput();
Image->Print(std::cout);
typedef itk::ChangeInformationImageFilter<InternalImageType> ChangeFilter;
ChangeFilter::Pointer change = ChangeFilter::New();
change->ChangeAll();
change->UseReferenceImageOn();
change->SetReferenceImage(Image);
change->SetInput(reader1->GetOutput());
try
{
change->Update();
}
catch (itk::ExceptionObject & excep)
{
std::cerr << "Exception caught in change !" << std::endl;
std::cerr << excep << std::endl;
return -1;
}
Image = change->GetOutput();
Image->Print(std::cout);
unsigned int nb_iteration = 1;
double rms = 0.;
double epsilon = 2.;
double curvature_weight = 0.;
double area_weight = 0.;
double reinitialization_weight = 0.;
double volume_weight = 0.5;
double volume = 0.5;
double l1 = 1.;
double l2 = 1.;
double overlap_weight = 100.0;
//
// We now define the image type using a particular pixel type and
// dimension. In this case the \code{float} type is used for the pixels
// due to the requirements of the smoothing filter.
//
const unsigned int Dimension = 3;
typedef float ScalarPixelType;
typedef itk::Image< ScalarPixelType, Dimension > InternalImageType;
typedef itk::Image< unsigned short, Dimension > MainType;
typedef itk::ScalarChanAndVeseLevelSetFunctionData< InternalImageType,
InternalImageType > DataHelperType;
typedef itk::ConstrainedRegionBasedLevelSetFunctionSharedData<
InternalImageType, InternalImageType, DataHelperType > SharedDataHelperType;
typedef itk::ScalarChanAndVeseLevelSetFunction< InternalImageType,
InternalImageType, SharedDataHelperType > LevelSetFunctionType;
// We declare now the type of the numerically discretized Step and Delta functions that
// will be used in the level-set computations for foreground and background regions
//
typedef itk::AtanRegularizedHeavisideStepFunction< ScalarPixelType,
ScalarPixelType > DomainFunctionType;
DomainFunctionType::Pointer domainFunction = DomainFunctionType::New();
domainFunction->SetEpsilon(epsilon);
typedef itk::ScalarChanAndVeseSparseLevelSetImageFilter< InternalImageType,
InternalImageType, InternalImageType, LevelSetFunctionType,
SharedDataHelperType > MultiLevelSetType;
MultiLevelSetType::Pointer levelSetFilter = MultiLevelSetType::New();
// We set the function count to 1 since a single level-set is being evolved.
//
levelSetFilter->SetFunctionCount(1);
// Set the feature image and initial level-set image as output of the
// fast marching image filter.
//
levelSetFilter->SetFeatureImage(ImgTypeConverter->GetOutput());
levelSetFilter->SetLevelSet(0, change->GetOutput());
//levelSetFilter->SetLevelSet(1, fastMarching1->GetOutput());
// Once activiated the level set evolution will stop if the convergence
// criteria or if the maximum number of iterations is reached. The
// convergence criteria is defined in terms of the root mean squared (RMS)
// change in the level set function. The evolution is said to have
// converged if the RMS change is below a user specified threshold. In a
// real application is desirable to couple the evolution of the zero set
// to a visualization module allowing the user to follow the evolution of
// the zero set. With this feedback, the user may decide when to stop the
// algorithm before the zero set leaks through the regions of low gradient
// in the contour of the anatomical structure to be segmented.
//
levelSetFilter->SetNumberOfIterations(nb_iteration);
levelSetFilter->SetMaximumRMSError(rms);
// Often, in real applications, images have different pixel resolutions. In such
// cases, it is best to use the native spacings to compute derivatives etc rather
// than sampling the images.
//
//levelSetFilter->SetUseImageSpacing(1);
// For large images, we may want to compute the level-set over the initial supplied
// level-set image. This saves a lot of memory.
//
levelSetFilter->SetInPlace(false);
// For the level set with phase 0, set different parameters and weights. This may
// to be set in a loop for the case of multiple level-sets evolving simultaneously.
//
for (unsigned int i = 0; i < 1; i++)
{
levelSetFilter->GetDifferenceFunction(i)->SetDomainFunction(domainFunction);
levelSetFilter->GetDifferenceFunction(i)->SetCurvatureWeight(curvature_weight);
levelSetFilter->GetDifferenceFunction(i)->SetAreaWeight(area_weight);
//levelSetFilter->GetDifferenceFunction(i)->SetOverlapPenaltyWeight(overlap_weight);
levelSetFilter->GetDifferenceFunction(i)->SetVolumeMatchingWeight(volume_weight);
levelSetFilter->GetDifferenceFunction(i)->SetVolume(volume);
levelSetFilter->GetDifferenceFunction(i)->SetLambda1(l1);
levelSetFilter->GetDifferenceFunction(i)->SetLambda2(l2);
}
try
{
levelSetFilter->Update();
}
catch (itk::ExceptionObject & excep)
{
std::cerr << "Exception caught !" << std::endl;
std::cerr << excep << std::endl;
return -1;
}
catch (...)
{
std::cerr << "Exception Caught!!!!!" << std::endl;
return -1;
}
Image = levelSetFilter->GetOutput();
Image->Print(std::cout);
writer->SetInput(levelSetFilter->GetOutput());
writer->Update();
return EXIT_SUCCESS;
}
Also I would like to know what difference it makes when I use 3D data for fast marching. By looking at the distance map it doesn't look like a shere centered at the given seed point in 3D.
And what data should ideally be passed to ChanVese from FastMarching for 3D dataset.
Pls help!
Param
L&T Technology Services Ltd
www.LntTechservices.com<http://www.lnttechservices.com/>
This Email may contain confidential or privileged information for the intended recipient (s). If you are not the intended recipient, please do not use or disseminate the information, notify the sender and delete it from your system.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20150210/38dad64e/attachment.html>
More information about the Insight-users
mailing list