SimpleITK  1.3.0.dev466
Python/GeodesicActiceContourSegmentation.py
1 #=========================================================================
2 #
3 # Copyright Insight Software Consortium
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 #=========================================================================
18 
19 from __future__ import print_function
20 
21 import SimpleITK as sitk
22 import os
23 import sys
24 
25 if len ( sys.argv ) < 8:
26  print( "Usage:", sys.argv[0], " <inputImage> <outputImage> <sigma> <InitialDistance> <PropagationScaling> <seedX> <seedY> <?seedZ>" )
27  sys.exit ( 1 )
28 
29 inputImage = sys.argv[1]
30 outputImage = sys.argv[2]
31 sigma = float(sys.argv[3])
32 initialDistance = float(sys.argv[4])
33 propagationScaling = float(sys.argv[5])
34 seed = [float(sys.argv[6]), float(sys.argv[7]) ]
35 
36 if len( sys.argv ) > 8:
37  seed.append( float(sys.argv[8]) )
38 
39 
40 reader = sitk.ImageFileReader()
41 reader.SetFileName ( inputImage )
42 image = reader.Execute()
43 
45 gradientMagnitude.SetSigma( sigma )
46 
47 featureImage = sitk.BoundedReciprocal( gradientMagnitude.Execute( image ) )
48 
49 seedImage = sitk.Image( image.GetSize()[0], image.GetSize()[1], sitk.sitkUInt8 )
50 seedImage.SetSpacing( image.GetSpacing() )
51 seedImage.SetOrigin( image.GetOrigin() )
52 seedImage.SetDirection( image.GetDirection() )
53 seedImage[ seedImage.TransformPhysicalPointToIndex(seed) ] = 1
54 
56 distance.InsideIsPositiveOff()
57 distance.UseImageSpacingOn()
58 
59 initialImage = sitk.BinaryThreshold( distance.Execute( seedImage ), -1000, 10 )
60 initialImage = sitk.Cast( initialImage, featureImage.GetPixelID() ) * -1 + 0.5
61 
62 
63 geodesicActiveContour = sitk.GeodesicActiveContourLevelSetImageFilter()
64 geodesicActiveContour.SetPropagationScaling( propagationScaling )
65 geodesicActiveContour.SetCurvatureScaling( .5 )
66 geodesicActiveContour.SetAdvectionScaling( 1.0 )
67 geodesicActiveContour.SetMaximumRMSError( 0.01 )
68 geodesicActiveContour.SetNumberOfIterations( 1000 )
69 
70 levelset = geodesicActiveContour.Execute( initialImage, featureImage )
71 
72 print( "RMS Change: ", geodesicActiveContour.GetRMSChange() )
73 print( "Elapsed Iterations: ", geodesicActiveContour.GetElapsedIterations() )
74 
75 contour = sitk.BinaryContour( sitk.BinaryThreshold( levelset, -1000, 0 ) )
76 
77 if ( not "SITK_NOSHOW" in os.environ ):
78  sitk.Show( sitk.LabelOverlay( image, contour ), "Levelset Countour" )