ITK  5.4.0
Insight Toolkit
SphinxExamples/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py
1 #!/usr/bin/env python
2 
3 # Copyright NumFOCUS
4 #
5 # Licensed under the Apache License, Version 2.0 (the "License");
6 # you may not use this file except in compliance with the License.
7 # You may obtain a copy of the License at
8 #
9 # https://www.apache.org/licenses/LICENSE-2.0.txt
10 #
11 # Unless required by applicable law or agreed to in writing, software
12 # distributed under the License is distributed on an "AS IS" BASIS,
13 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 # See the License for the specific language governing permissions and
15 # limitations under the License.
16 
17 import itk
18 import argparse
19 
20 parser = argparse.ArgumentParser(
21  description="Segment With Geodesic Active Contour Level Set."
22 )
23 parser.add_argument("input_image")
24 parser.add_argument("output_image")
25 parser.add_argument("seed_x", type=int)
26 parser.add_argument("seed_y", type=int)
27 parser.add_argument("initial_distance", type=float)
28 parser.add_argument("sigma", type=float)
29 parser.add_argument("sigmoid_alpha", type=float)
30 parser.add_argument("sigmoid_beta", type=float)
31 parser.add_argument("propagation_scaling", type=float)
32 parser.add_argument("number_of_iterations", type=int)
33 args = parser.parse_args()
34 
35 seedValue = -args.initial_distance
36 
37 Dimension = 2
38 
39 InputPixelType = itk.F
40 OutputPixelType = itk.UC
41 
42 InputImageType = itk.Image[InputPixelType, Dimension]
43 OutputImageType = itk.Image[OutputPixelType, Dimension]
44 
45 ReaderType = itk.ImageFileReader[InputImageType]
46 WriterType = itk.ImageFileWriter[OutputImageType]
47 
48 reader = ReaderType.New()
49 reader.SetFileName(args.input_image)
50 
52  InputImageType, InputImageType
53 ]
54 smoothing = SmoothingFilterType.New()
55 smoothing.SetTimeStep(0.125)
56 smoothing.SetNumberOfIterations(5)
57 smoothing.SetConductanceParameter(9.0)
58 smoothing.SetInput(reader.GetOutput())
59 
61  InputImageType, InputImageType
62 ]
63 gradientMagnitude = GradientFilterType.New()
64 gradientMagnitude.SetSigma(args.sigma)
65 gradientMagnitude.SetInput(smoothing.GetOutput())
66 
67 SigmoidFilterType = itk.SigmoidImageFilter[InputImageType, InputImageType]
68 sigmoid = SigmoidFilterType.New()
69 sigmoid.SetOutputMinimum(0.0)
70 sigmoid.SetOutputMaximum(1.0)
71 sigmoid.SetAlpha(args.sigmoid_alpha)
72 sigmoid.SetBeta(args.sigmoid_beta)
73 sigmoid.SetInput(gradientMagnitude.GetOutput())
74 
75 FastMarchingFilterType = itk.FastMarchingImageFilter[InputImageType, InputImageType]
76 fastMarching = FastMarchingFilterType.New()
77 
78 GeoActiveContourFilterType = itk.GeodesicActiveContourLevelSetImageFilter[
79  InputImageType, InputImageType, InputPixelType
80 ]
81 geodesicActiveContour = GeoActiveContourFilterType.New()
82 geodesicActiveContour.SetPropagationScaling(args.propagation_scaling)
83 geodesicActiveContour.SetCurvatureScaling(1.0)
84 geodesicActiveContour.SetAdvectionScaling(1.0)
85 geodesicActiveContour.SetMaximumRMSError(0.02)
86 geodesicActiveContour.SetNumberOfIterations(args.number_of_iterations)
87 geodesicActiveContour.SetInput(fastMarching.GetOutput())
88 geodesicActiveContour.SetFeatureImage(sigmoid.GetOutput())
89 
90 ThresholdingFilterType = itk.BinaryThresholdImageFilter[InputImageType, OutputImageType]
91 thresholder = ThresholdingFilterType.New()
92 thresholder.SetLowerThreshold(-1000.0)
93 thresholder.SetUpperThreshold(0.0)
94 thresholder.SetOutsideValue(itk.NumericTraits[OutputPixelType].min())
95 thresholder.SetInsideValue(itk.NumericTraits[OutputPixelType].max())
96 thresholder.SetInput(geodesicActiveContour.GetOutput())
97 
98 seedPosition = itk.Index[Dimension]()
99 seedPosition[0] = args.seed_x
100 seedPosition[1] = args.seed_y
101 
102 node = itk.LevelSetNode[InputPixelType, Dimension]()
103 node.SetValue(seedValue)
104 node.SetIndex(seedPosition)
105 
106 seeds = itk.VectorContainer[itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New()
107 seeds.Initialize()
108 seeds.InsertElement(0, node)
109 
110 fastMarching.SetTrialPoints(seeds)
111 fastMarching.SetSpeedConstant(1.0)
112 
113 CastFilterType = itk.RescaleIntensityImageFilter[InputImageType, OutputImageType]
114 
115 caster1 = CastFilterType.New()
116 caster2 = CastFilterType.New()
117 caster3 = CastFilterType.New()
118 caster4 = CastFilterType.New()
119 
120 writer1 = WriterType.New()
121 writer2 = WriterType.New()
122 writer3 = WriterType.New()
123 writer4 = WriterType.New()
124 
125 caster1.SetInput(smoothing.GetOutput())
126 writer1.SetInput(caster1.GetOutput())
127 writer1.SetFileName("GeodesicActiveContourImageFilterOutput1.png")
128 caster1.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
129 caster1.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
130 writer1.Update()
131 
132 caster2.SetInput(gradientMagnitude.GetOutput())
133 writer2.SetInput(caster2.GetOutput())
134 writer2.SetFileName("GeodesicActiveContourImageFilterOutput2.png")
135 caster2.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
136 caster2.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
137 writer2.Update()
138 
139 caster3.SetInput(sigmoid.GetOutput())
140 writer3.SetInput(caster3.GetOutput())
141 writer3.SetFileName("GeodesicActiveContourImageFilterOutput3.png")
142 caster3.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
143 caster3.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
144 writer3.Update()
145 
146 caster4.SetInput(fastMarching.GetOutput())
147 writer4.SetInput(caster4.GetOutput())
148 writer4.SetFileName("GeodesicActiveContourImageFilterOutput4.png")
149 caster4.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
150 caster4.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
151 
152 fastMarching.SetOutputSize(reader.GetOutput().GetBufferedRegion().GetSize())
153 
154 writer = WriterType.New()
155 writer.SetFileName(args.output_image)
156 writer.SetInput(thresholder.GetOutput())
157 writer.Update()
158 
159 print(
160  "Max. no. iterations: " + str(geodesicActiveContour.GetNumberOfIterations()) + "\n"
161 )
162 print("Max. RMS error: " + str(geodesicActiveContour.GetMaximumRMSError()) + "\n")
163 print(
164  "No. elpased iterations: "
165  + str(geodesicActiveContour.GetElapsedIterations())
166  + "\n"
167 )
168 print("RMS change: " + str(geodesicActiveContour.GetRMSChange()) + "\n")
169 
170 writer4.Update()
171 
172 InternalWriterType = itk.ImageFileWriter[InputImageType]
173 
174 mapWriter = InternalWriterType.New()
175 mapWriter.SetInput(fastMarching.GetOutput())
176 mapWriter.SetFileName("GeodesicActiveContourImageFilterOutput4.mha")
177 mapWriter.Update()
178 
179 speedWriter = InternalWriterType.New()
180 speedWriter.SetInput(sigmoid.GetOutput())
181 speedWriter.SetFileName("GeodesicActiveContourImageFilterOutput3.mha")
182 speedWriter.Update()
183 
184 gradientWriter = InternalWriterType.New()
185 gradientWriter.SetInput(gradientMagnitude.GetOutput())
186 gradientWriter.SetFileName("GeodesicActiveContourImageFilterOutput2.mha")
187 gradientWriter.Update()
itk::Index
Represent a n-dimensional index in a n-dimensional image.
Definition: itkIndex.h:70
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
itk::SigmoidImageFilter
Computes the sigmoid function pixel-wise.
Definition: itkSigmoidImageFilter.h:144
itk::LevelSetNode
Represent a node in a level set.
Definition: itkLevelSetNode.h:45
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::CurvatureAnisotropicDiffusionImageFilter
This filter performs anisotropic diffusion on a scalar itk::Image using the modified curvature diffus...
Definition: itkCurvatureAnisotropicDiffusionImageFilter.h:58
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:88
itk::FastMarchingImageFilter
Solve an Eikonal equation using Fast Marching.
Definition: itkFastMarchingImageFilter.h:135
itk::GradientMagnitudeRecursiveGaussianImageFilter
Computes the Magnitude of the Gradient of an image by convolution with the first derivative of a Gaus...
Definition: itkGradientMagnitudeRecursiveGaussianImageFilter.h:50
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:59
itk::RescaleIntensityImageFilter
Applies a linear transformation to the intensity levels of the input Image.
Definition: itkRescaleIntensityImageFilter.h:133
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
New
static Pointer New()
itk::VectorContainer
Define a front-end to the STL "vector" container that conforms to the IndexedContainerInterface.
Definition: itkVectorContainer.h:48