Python SimpleITK 水平集

GeodesicActiveContourLevelSetImageFilter是一个水平集图像分割方法，根据已有的初始轮廓，向内/外延伸并找到分割边缘。

输入

• initial image：初始水平集，即第0层的图像。把初始轮廓代入signed distance function（比如SignedMaurerDistanceMapImageFilter）。这个初始轮廓可以和图像的分割边缘重合或相交。
• feature image：待分割的原图的边缘图（edge map），通常用gradient计算。在边缘图中，边缘处的值一般接近于0，在分割区域内的值接近于1。

输出

图像矩阵，以正负值区分分割区域的不同区域，0代表边缘/传播面。

# =========================================================================
 #
 #  Copyright NumFOCUS
 #
 #  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 os
import sys
  
import SimpleITK as sitk
  
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()

# 使用GradientMagnitudeRecursiveGaussianImageFilter对原图进行处理 -> 得到featureImage
gradientMagnitude = sitk.GradientMagnitudeRecursiveGaussianImageFilter()
gradientMagnitude.SetSigma(sigma)
featureImage = sitk.BoundedReciprocal(gradientMagnitude.Execute(image))

# 得到初始轮廓（mask）
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

# 使用SignedMaurerDistanceMapImageFilter对初始轮廓处理 -> 得到initialImage
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) # 控制圆滑/smooth程度
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 ("SITK_NOSHOW" not in os.environ):
  sitk.Show(sitk.LabelOverlay(image, contour), "Levelset Countour")

Reference

“Geodesic Active Contours”, V. Caselles, R. Kimmel and G. Sapiro. International Journal on Computer Vision, Vol 22, No. 1, pp 61-97, 1997