ITK  4.13.0
Insight Segmentation and Registration Toolkit
SphinxExamples/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py
1 #!/usr/bin/env python
2 
3 #==========================================================================
4 #
5 # Copyright Insight Software Consortium
6 #
7 # Licensed under the Apache License, Version 2.0 (the "License");
8 # you may not use this file except in compliance with the License.
9 # You may obtain a copy of the License at
10 #
11 # http://www.apache.org/licenses/LICENSE-2.0.txt
12 #
13 # Unless required by applicable law or agreed to in writing, software
14 # distributed under the License is distributed on an "AS IS" BASIS,
15 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 # See the License for the specific language governing permissions and
17 # limitations under the License.
18 #
19 #==========================================================================*/
20 
21 import sys
22 import itk
23 
24 if len(sys.argv) != 11:
25  print(
26  "Usage: " + sys.argv[0] + "<InputFileName> <OutputFileName>"
27  " <seedX> <seedY> <InitialDistance> <Sigma> <SigmoidAlpha> "
28  "<SigmoidBeta> <PropagationScaling> <NumberOfIterations>")
29  sys.exit(1)
30 
31 inputFileName = sys.argv[1]
32 outputFileName = sys.argv[2]
33 seedPosX = int(sys.argv[3])
34 seedPosY = int(sys.argv[4])
35 
36 initialDistance = float(sys.argv[5])
37 sigma = float(sys.argv[6])
38 alpha = float(sys.argv[7])
39 beta = float(sys.argv[8])
40 propagationScaling = float(sys.argv[9])
41 numberOfIterations = int(sys.argv[10])
42 seedValue = -initialDistance
43 
44 Dimension = 2
45 
46 InputPixelType = itk.F
47 OutputPixelType = itk.UC
48 
49 InputImageType = itk.Image[InputPixelType, Dimension]
50 OutputImageType = itk.Image[OutputPixelType, Dimension]
51 
52 ReaderType = itk.ImageFileReader[InputImageType]
53 WriterType = itk.ImageFileWriter[OutputImageType]
54 
55 reader = ReaderType.New()
56 reader.SetFileName(inputFileName)
57 
58 SmoothingFilterType = itk.CurvatureAnisotropicDiffusionImageFilter[
59  InputImageType, InputImageType]
60 smoothing = SmoothingFilterType.New()
61 smoothing.SetTimeStep(0.125)
62 smoothing.SetNumberOfIterations(5)
63 smoothing.SetConductanceParameter(9.0)
64 smoothing.SetInput(reader.GetOutput())
65 
66 GradientFilterType = itk.GradientMagnitudeRecursiveGaussianImageFilter[
67  InputImageType, InputImageType]
68 gradientMagnitude = GradientFilterType.New()
69 gradientMagnitude.SetSigma(sigma)
70 gradientMagnitude.SetInput(smoothing.GetOutput())
71 
72 SigmoidFilterType = itk.SigmoidImageFilter[InputImageType, InputImageType]
73 sigmoid = SigmoidFilterType.New()
74 sigmoid.SetOutputMinimum(0.0)
75 sigmoid.SetOutputMaximum(1.0)
76 sigmoid.SetAlpha(alpha)
77 sigmoid.SetBeta(beta)
78 sigmoid.SetInput(gradientMagnitude.GetOutput())
79 
80 FastMarchingFilterType = itk.FastMarchingImageFilter[
81  InputImageType, InputImageType]
82 fastMarching = FastMarchingFilterType.New()
83 
84 GeoActiveContourFilterType = itk.GeodesicActiveContourLevelSetImageFilter[
85  InputImageType, InputImageType, InputPixelType]
86 geodesicActiveContour = GeoActiveContourFilterType.New()
87 geodesicActiveContour.SetPropagationScaling(propagationScaling)
88 geodesicActiveContour.SetCurvatureScaling(1.0)
89 geodesicActiveContour.SetAdvectionScaling(1.0)
90 geodesicActiveContour.SetMaximumRMSError(0.02)
91 geodesicActiveContour.SetNumberOfIterations(numberOfIterations)
92 geodesicActiveContour.SetInput(fastMarching.GetOutput())
93 geodesicActiveContour.SetFeatureImage(sigmoid.GetOutput())
94 
95 ThresholdingFilterType = itk.BinaryThresholdImageFilter[
96  InputImageType, OutputImageType]
97 thresholder = ThresholdingFilterType.New()
98 thresholder.SetLowerThreshold(-1000.0)
99 thresholder.SetUpperThreshold(0.0)
100 thresholder.SetOutsideValue(itk.NumericTraits[OutputPixelType].min())
101 thresholder.SetInsideValue(itk.NumericTraits[OutputPixelType].max())
102 thresholder.SetInput(geodesicActiveContour.GetOutput())
103 
104 seedPosition = itk.Index[Dimension]()
105 seedPosition[0] = seedPosX
106 seedPosition[1] = seedPosY
107 
108 node = itk.LevelSetNode[InputPixelType, Dimension]()
109 node.SetValue(seedValue)
110 node.SetIndex(seedPosition)
111 
112 seeds = itk.VectorContainer[
113  itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New()
114 seeds.Initialize()
115 seeds.InsertElement(0, node)
116 
117 fastMarching.SetTrialPoints(seeds)
118 fastMarching.SetSpeedConstant(1.0)
119 
120 CastFilterType = itk.RescaleIntensityImageFilter[
121  InputImageType, OutputImageType]
122 
123 caster1 = CastFilterType.New()
124 caster2 = CastFilterType.New()
125 caster3 = CastFilterType.New()
126 caster4 = CastFilterType.New()
127 
128 writer1 = WriterType.New()
129 writer2 = WriterType.New()
130 writer3 = WriterType.New()
131 writer4 = WriterType.New()
132 
133 caster1.SetInput(smoothing.GetOutput())
134 writer1.SetInput(caster1.GetOutput())
135 writer1.SetFileName("GeodesicActiveContourImageFilterOutput1.png")
136 caster1.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
137 caster1.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
138 writer1.Update()
139 
140 caster2.SetInput(gradientMagnitude.GetOutput())
141 writer2.SetInput(caster2.GetOutput())
142 writer2.SetFileName("GeodesicActiveContourImageFilterOutput2.png")
143 caster2.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
144 caster2.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
145 writer2.Update()
146 
147 caster3.SetInput(sigmoid.GetOutput())
148 writer3.SetInput(caster3.GetOutput())
149 writer3.SetFileName("GeodesicActiveContourImageFilterOutput3.png")
150 caster3.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
151 caster3.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
152 writer3.Update()
153 
154 caster4.SetInput(fastMarching.GetOutput())
155 writer4.SetInput(caster4.GetOutput())
156 writer4.SetFileName("GeodesicActiveContourImageFilterOutput4.png")
157 caster4.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
158 caster4.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
159 
160 fastMarching.SetOutputSize(
161  reader.GetOutput().GetBufferedRegion().GetSize())
162 
163 writer = WriterType.New()
164 writer.SetFileName(outputFileName)
165 writer.SetInput(thresholder.GetOutput())
166 writer.Update()
167 
168 print(
169  "Max. no. iterations: " +
170  str(geodesicActiveContour.GetNumberOfIterations()) + "\n")
171 print(
172  "Max. RMS error: " +
173  str(geodesicActiveContour.GetMaximumRMSError()) + "\n")
174 print(
175  "No. elpased iterations: " +
176  str(geodesicActiveContour.GetElapsedIterations()) + "\n")
177 print("RMS change: " + str(geodesicActiveContour.GetRMSChange()) + "\n")
178 
179 writer4.Update()
180 
181 InternalWriterType = itk.ImageFileWriter[InputImageType]
182 
183 mapWriter = InternalWriterType.New()
184 mapWriter.SetInput(fastMarching.GetOutput())
185 mapWriter.SetFileName("GeodesicActiveContourImageFilterOutput4.mha")
186 mapWriter.Update()
187 
188 speedWriter = InternalWriterType.New()
189 speedWriter.SetInput(sigmoid.GetOutput())
190 speedWriter.SetFileName("GeodesicActiveContourImageFilterOutput3.mha")
191 speedWriter.Update()
192 
193 gradientWriter = InternalWriterType.New()
194 gradientWriter.SetInput(gradientMagnitude.GetOutput())
195 gradientWriter.SetFileName("GeodesicActiveContourImageFilterOutput2.mha")
196 gradientWriter.Update()
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Tue Dec 19 2017 03:56:59 for ITK by   doxygen 1.8.5