ITK  5.2.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 # http://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 sys
18 import itk
19 
20 if len(sys.argv) != 11:
21  print(
22  "Usage: " + sys.argv[0] + "<InputFileName> <OutputFileName>"
23  " <seedX> <seedY> <InitialDistance> <Sigma> <SigmoidAlpha> "
24  "<SigmoidBeta> <PropagationScaling> <NumberOfIterations>")
25  sys.exit(1)
26 
27 inputFileName = sys.argv[1]
28 outputFileName = sys.argv[2]
29 seedPosX = int(sys.argv[3])
30 seedPosY = int(sys.argv[4])
31 
32 initialDistance = float(sys.argv[5])
33 sigma = float(sys.argv[6])
34 alpha = float(sys.argv[7])
35 beta = float(sys.argv[8])
36 propagationScaling = float(sys.argv[9])
37 numberOfIterations = int(sys.argv[10])
38 seedValue = -initialDistance
39 
40 Dimension = 2
41 
42 InputPixelType = itk.F
43 OutputPixelType = itk.UC
44 
45 InputImageType = itk.Image[InputPixelType, Dimension]
46 OutputImageType = itk.Image[OutputPixelType, Dimension]
47 
48 ReaderType = itk.ImageFileReader[InputImageType]
49 WriterType = itk.ImageFileWriter[OutputImageType]
50 
51 reader = ReaderType.New()
52 reader.SetFileName(inputFileName)
53 
55  InputImageType, InputImageType]
56 smoothing = SmoothingFilterType.New()
57 smoothing.SetTimeStep(0.125)
58 smoothing.SetNumberOfIterations(5)
59 smoothing.SetConductanceParameter(9.0)
60 smoothing.SetInput(reader.GetOutput())
61 
63  InputImageType, InputImageType]
64 gradientMagnitude = GradientFilterType.New()
65 gradientMagnitude.SetSigma(sigma)
66 gradientMagnitude.SetInput(smoothing.GetOutput())
67 
68 SigmoidFilterType = itk.SigmoidImageFilter[InputImageType, InputImageType]
69 sigmoid = SigmoidFilterType.New()
70 sigmoid.SetOutputMinimum(0.0)
71 sigmoid.SetOutputMaximum(1.0)
72 sigmoid.SetAlpha(alpha)
73 sigmoid.SetBeta(beta)
74 sigmoid.SetInput(gradientMagnitude.GetOutput())
75 
76 FastMarchingFilterType = itk.FastMarchingImageFilter[
77  InputImageType, InputImageType]
78 fastMarching = FastMarchingFilterType.New()
79 
80 GeoActiveContourFilterType = itk.GeodesicActiveContourLevelSetImageFilter[
81  InputImageType, InputImageType, InputPixelType]
82 geodesicActiveContour = GeoActiveContourFilterType.New()
83 geodesicActiveContour.SetPropagationScaling(propagationScaling)
84 geodesicActiveContour.SetCurvatureScaling(1.0)
85 geodesicActiveContour.SetAdvectionScaling(1.0)
86 geodesicActiveContour.SetMaximumRMSError(0.02)
87 geodesicActiveContour.SetNumberOfIterations(numberOfIterations)
88 geodesicActiveContour.SetInput(fastMarching.GetOutput())
89 geodesicActiveContour.SetFeatureImage(sigmoid.GetOutput())
90 
91 ThresholdingFilterType = itk.BinaryThresholdImageFilter[
92  InputImageType, OutputImageType]
93 thresholder = ThresholdingFilterType.New()
94 thresholder.SetLowerThreshold(-1000.0)
95 thresholder.SetUpperThreshold(0.0)
96 thresholder.SetOutsideValue(itk.NumericTraits[OutputPixelType].min())
97 thresholder.SetInsideValue(itk.NumericTraits[OutputPixelType].max())
98 thresholder.SetInput(geodesicActiveContour.GetOutput())
99 
100 seedPosition = itk.Index[Dimension]()
101 seedPosition[0] = seedPosX
102 seedPosition[1] = seedPosY
103 
104 node = itk.LevelSetNode[InputPixelType, Dimension]()
105 node.SetValue(seedValue)
106 node.SetIndex(seedPosition)
107 
108 seeds = itk.VectorContainer[
109  itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New()
110 seeds.Initialize()
111 seeds.InsertElement(0, node)
112 
113 fastMarching.SetTrialPoints(seeds)
114 fastMarching.SetSpeedConstant(1.0)
115 
116 CastFilterType = itk.RescaleIntensityImageFilter[
117  InputImageType, OutputImageType]
118 
119 caster1 = CastFilterType.New()
120 caster2 = CastFilterType.New()
121 caster3 = CastFilterType.New()
122 caster4 = CastFilterType.New()
123 
124 writer1 = WriterType.New()
125 writer2 = WriterType.New()
126 writer3 = WriterType.New()
127 writer4 = WriterType.New()
128 
129 caster1.SetInput(smoothing.GetOutput())
130 writer1.SetInput(caster1.GetOutput())
131 writer1.SetFileName("GeodesicActiveContourImageFilterOutput1.png")
132 caster1.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
133 caster1.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
134 writer1.Update()
135 
136 caster2.SetInput(gradientMagnitude.GetOutput())
137 writer2.SetInput(caster2.GetOutput())
138 writer2.SetFileName("GeodesicActiveContourImageFilterOutput2.png")
139 caster2.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
140 caster2.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
141 writer2.Update()
142 
143 caster3.SetInput(sigmoid.GetOutput())
144 writer3.SetInput(caster3.GetOutput())
145 writer3.SetFileName("GeodesicActiveContourImageFilterOutput3.png")
146 caster3.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
147 caster3.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
148 writer3.Update()
149 
150 caster4.SetInput(fastMarching.GetOutput())
151 writer4.SetInput(caster4.GetOutput())
152 writer4.SetFileName("GeodesicActiveContourImageFilterOutput4.png")
153 caster4.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
154 caster4.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
155 
156 fastMarching.SetOutputSize(
157  reader.GetOutput().GetBufferedRegion().GetSize())
158 
159 writer = WriterType.New()
160 writer.SetFileName(outputFileName)
161 writer.SetInput(thresholder.GetOutput())
162 writer.Update()
163 
164 print(
165  "Max. no. iterations: " +
166  str(geodesicActiveContour.GetNumberOfIterations()) + "\n")
167 print(
168  "Max. RMS error: " +
169  str(geodesicActiveContour.GetMaximumRMSError()) + "\n")
170 print(
171  "No. elpased iterations: " +
172  str(geodesicActiveContour.GetElapsedIterations()) + "\n")
173 print("RMS change: " + str(geodesicActiveContour.GetRMSChange()) + "\n")
174 
175 writer4.Update()
176 
177 InternalWriterType = itk.ImageFileWriter[InputImageType]
178 
179 mapWriter = InternalWriterType.New()
180 mapWriter.SetInput(fastMarching.GetOutput())
181 mapWriter.SetFileName("GeodesicActiveContourImageFilterOutput4.mha")
182 mapWriter.Update()
183 
184 speedWriter = InternalWriterType.New()
185 speedWriter.SetInput(sigmoid.GetOutput())
186 speedWriter.SetFileName("GeodesicActiveContourImageFilterOutput3.mha")
187 speedWriter.Update()
188 
189 gradientWriter = InternalWriterType.New()
190 gradientWriter.SetInput(gradientMagnitude.GetOutput())
191 gradientWriter.SetFileName("GeodesicActiveContourImageFilterOutput2.mha")
192 gradientWriter.Update()
itk::Index
Represent a n-dimensional index in a n-dimensional image.
Definition: itkIndex.h:66
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:139
itk::SigmoidImageFilter
Computes the sigmoid function pixel-wise.
Definition: itkSigmoidImageFilter.h:150
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:87
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
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:58
itk::RescaleIntensityImageFilter
Applies a linear transformation to the intensity levels of the input Image.
Definition: itkRescaleIntensityImageFilter.h:154
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::VectorContainer
Define a front-end to the STL "vector" container that conforms to the IndexedContainerInterface.
Definition: itkVectorContainer.h:48