~ubuntu-branches/debian/sid/simpleitk/sid

#=========================================================================

#

#  Copyright Insight Software Consortium

#

#  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.

#

#=========================================================================

 

from __future__ import print_function

 

import SimpleITK as sitk

import os

import sys

 

if len ( sys.argv ) < 8:

    print( "Usage:", sys.argv[0], " <inputImage> <outputImage> <sigma> <InitialDistance> <PropagationScaling> <seedX> <seedY> <?seedZ>" )

    sys.exit ( 1 )

 

inputImage = sys.argv[1]

outputImage = sys.argv[2]

sigma = float(sys.argv[3])

initialDistance = float(sys.argv[4])

propagationScaling = float(sys.argv[5])

seed = [float(sys.argv[6]), float(sys.argv[7]) ]

 

if len( sys.argv ) > 8:

    seed.append( float(sys.argv[8]) )

 

 

reader = sitk.ImageFileReader()

reader.SetFileName ( inputImage )

image = reader.Execute()

 

gradientMagnitude = sitk.GradientMagnitudeRecursiveGaussianImageFilter()

gradientMagnitude.SetSigma( sigma )

 

featureImage = sitk.BoundedReciprocal( gradientMagnitude.Execute( image ) )

 

seedImage = sitk.Image( image.GetSize()[0], image.GetSize()[1], sitk.sitkUInt8 )

seedImage.SetSpacing( image.GetSpacing() )

seedImage.SetOrigin( image.GetOrigin() )

seedImage.SetDirection( image.GetDirection() )

seedImage[ seedImage.TransformPhysicalPointToIndex(seed) ] = 1

 

distance = sitk.SignedMaurerDistanceMapImageFilter()

distance.InsideIsPositiveOff()

distance.UseImageSpacingOn()

 

initialImage = sitk.BinaryThreshold( distance.Execute( seedImage ), -1000, 10 )

initialImage = sitk.Cast( initialImage, featureImage.GetPixelID() ) * -1 + 0.5

 

 

geodesicActiveContour = sitk.GeodesicActiveContourLevelSetImageFilter()

geodesicActiveContour.SetPropagationScaling( propagationScaling )

geodesicActiveContour.SetCurvatureScaling( .5 )

geodesicActiveContour.SetAdvectionScaling( 1.0 )

geodesicActiveContour.SetMaximumRMSError( 0.01 )

geodesicActiveContour.SetNumberOfIterations( 1000 )

 

levelset = geodesicActiveContour.Execute( initialImage, featureImage )

 

print( "RMS Change: ", geodesicActiveContour.GetRMSChange() )

print( "Elapsed Iterations: ", geodesicActiveContour.GetElapsedIterations() )

 

contour = sitk.BinaryContour( sitk.BinaryThreshold( levelset, -1000, 0 ) )

 

if ( not "SITK_NOSHOW" in os.environ ):

    sitk.Show( sitk.LabelOverlay( image, contour ), "Levelset Countour" )