Python SimpleITK 水平集_栖迟於一丘-CSDN博客
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 osimport 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对原图进行处理 -> 得到featureImagegradientMagnitude = 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对初始轮廓处理 -> 得到initialImagedistance = 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")12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
Reference
“Geodesic Active Contours”, V. Caselles, R. Kimmel and G. Sapiro. International Journal on Computer Vision, Vol 22, No. 1, pp 61-97, 1997
赞 (0)