ITK  5.3.0
Insight Toolkit
SphinxExamples/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
int
main(int argc, char * argv[])
{
if (argc != 11)
{
std::cerr << "Usage: " << argv[0];
std::cerr << " <InputFileName> <OutputFileName>";
std::cerr << " <seedX> <seedY> <InitialDistance>";
std::cerr << " <Sigma> <SigmoidAlpha> <SigmoidBeta>";
std::cerr << " <PropagationScaling> <NumberOfIterations>" << std::endl;
return EXIT_FAILURE;
}
const char * inputFileName = argv[1];
const char * outputFileName = argv[2];
const int seedPosX = std::stoi(argv[3]);
const int seedPosY = std::stoi(argv[4]);
const double initialDistance = std::stod(argv[5]);
const double sigma = std::stod(argv[6]);
const double alpha = std::stod(argv[7]);
const double beta = std::stod(argv[8]);
const double propagationScaling = std::stod(argv[9]);
const double numberOfIterations = std::stoi(argv[10]);
const double seedValue = -initialDistance;
constexpr unsigned int Dimension = 2;
using InputPixelType = float;
using InputImageType = itk::Image<InputPixelType, Dimension>;
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
const auto input = itk::ReadImage<InputImageType>(inputFileName);
SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
smoothing->SetTimeStep(0.125);
smoothing->SetNumberOfIterations(5);
smoothing->SetConductanceParameter(9.0);
smoothing->SetInput(input);
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
gradientMagnitude->SetSigma(sigma);
gradientMagnitude->SetInput(smoothing->GetOutput());
SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
sigmoid->SetOutputMinimum(0.0);
sigmoid->SetOutputMaximum(1.0);
sigmoid->SetAlpha(alpha);
sigmoid->SetBeta(beta);
sigmoid->SetInput(gradientMagnitude->GetOutput());
FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();
geodesicActiveContour->SetPropagationScaling(propagationScaling);
geodesicActiveContour->SetCurvatureScaling(1.0);
geodesicActiveContour->SetAdvectionScaling(1.0);
geodesicActiveContour->SetMaximumRMSError(0.02);
geodesicActiveContour->SetNumberOfIterations(numberOfIterations);
geodesicActiveContour->SetInput(fastMarching->GetOutput());
geodesicActiveContour->SetFeatureImage(sigmoid->GetOutput());
ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
thresholder->SetLowerThreshold(-1000.0);
thresholder->SetUpperThreshold(0.0);
thresholder->SetOutsideValue(itk::NumericTraits<OutputPixelType>::min());
thresholder->SetInsideValue(itk::NumericTraits<OutputPixelType>::max());
thresholder->SetInput(geodesicActiveContour->GetOutput());
using NodeContainer = FastMarchingFilterType::NodeContainer;
using NodeType = FastMarchingFilterType::NodeType;
seedPosition[0] = seedPosX;
seedPosition[1] = seedPosY;
NodeContainer::Pointer seeds = NodeContainer::New();
NodeType node;
node.SetValue(seedValue);
node.SetIndex(seedPosition);
seeds->Initialize();
seeds->InsertElement(0, node);
fastMarching->SetTrialPoints(seeds);
fastMarching->SetSpeedConstant(1.0);
CastFilterType::Pointer caster1 = CastFilterType::New();
CastFilterType::Pointer caster2 = CastFilterType::New();
CastFilterType::Pointer caster3 = CastFilterType::New();
CastFilterType::Pointer caster4 = CastFilterType::New();
caster1->SetInput(smoothing->GetOutput());
caster1->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
caster1->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
itk::WriteImage(caster1->GetOutput(), "GeodesicActiveContourImageFilterOutput1.png");
caster2->SetInput(gradientMagnitude->GetOutput());
caster2->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
caster2->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
itk::WriteImage(caster2->GetOutput(), "GeodesicActiveContourImageFilterOutput2.png");
caster3->SetInput(sigmoid->GetOutput());
caster3->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
caster3->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
itk::WriteImage(caster3->GetOutput(), "GeodesicActiveContourImageFilterOutput3.png");
caster4->SetInput(fastMarching->GetOutput());
caster4->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
caster4->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
fastMarching->SetOutputDirection(input->GetDirection());
fastMarching->SetOutputOrigin(input->GetOrigin());
fastMarching->SetOutputRegion(input->GetBufferedRegion());
fastMarching->SetOutputSpacing(input->GetSpacing());
try
{
itk::WriteImage(thresholder->GetOutput(), outputFileName);
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
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;
try
{
itk::WriteImage(caster4->GetOutput(), "GeodesicActiveContourImageFilterOutput4.png");
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
try
{
itk::WriteImage(fastMarching->GetOutput(), "GeodesicActiveContourImageFilterOutput4.mha");
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
try
{
itk::WriteImage(sigmoid->GetOutput(), "GeodesicActiveContourImageFilterOutput3.mha");
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
try
{
itk::WriteImage(gradientMagnitude->GetOutput(), "GeodesicActiveContourImageFilterOutput2.mha");
}
catch (itk::ExceptionObject & error)
{
std::cerr << "Error: " << error << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
itkSigmoidImageFilter.h
itk::GeodesicActiveContourLevelSetImageFilter
Segments structures in images based on a user supplied edge potential map.
Definition: itkGeodesicActiveContourLevelSetImageFilter.h:105
itk::BinaryThresholdImageFilter
Binarize an input image by thresholding.
Definition: itkBinaryThresholdImageFilter.h:132
itkImageFileReader.h
itk::SigmoidImageFilter
Computes the sigmoid function pixel-wise.
Definition: itkSigmoidImageFilter.h:144
itkFastMarchingImageFilter.h
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itk::CurvatureAnisotropicDiffusionImageFilter
This filter performs anisotropic diffusion on a scalar itk::Image using the modified curvature diffus...
Definition: itkCurvatureAnisotropicDiffusionImageFilter.h:58
itkCurvatureAnisotropicDiffusionImageFilter.h
itk::FastMarchingImageFilter
Solve an Eikonal equation using Fast Marching.
Definition: itkFastMarchingImageFilter.h:136
itk::GradientMagnitudeRecursiveGaussianImageFilter
Computes the Magnitude of the Gradient of an image by convolution with the first derivative of a Gaus...
Definition: itkGradientMagnitudeRecursiveGaussianImageFilter.h:50
itkRescaleIntensityImageFilter.h
itkImageFileWriter.h
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:58
itkGeodesicActiveContourLevelSetImageFilter.h
itk::RescaleIntensityImageFilter
Applies a linear transformation to the intensity levels of the input Image.
Definition: itkRescaleIntensityImageFilter.h:147
itkBinaryThresholdImageFilter.h
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itk::WriteImage
ITK_TEMPLATE_EXPORT void WriteImage(TImagePointer &&image, const std::string &filename, bool compress=false)
Definition: itkImageFileWriter.h:254
itkGradientMagnitudeRecursiveGaussianImageFilter.h