医学图像处理_医学图像处理案例（二十）——医学图像处理案例代码详解
发布时间
阅读量：
阅读量
在前面分享的医学图像处理案例中，给出了很多具体案例，但有些读者还是渴望可以深入分享案例代码详解。那么今天我将从骨骼分割，气管分割，肺组织分割，血管分割这四个具体案例来详细讲解如何来实现。 1、骨骼分割代码详解 1.1、采用固定阈值分割，一般骨骼的HU值在300以上，但是发现心脏和主动脉的HU值与一些骨骼的HU比较接近，所以我们设置固定阈值为100，把骨骼和心脏及主动脉都分割出来。
1.2、 首先采用形态学开操作，将骨骼和心脏和主动脉连接的部分断开，然后再取最大连通域可以得到粗略的心脏和主动脉图像。
1.3、 再将步骤1.1的结果与步骤1.2的结果相减。
1.4、 对步骤1.3的结果求取最大连通域得到骨骼Mask区域。
1.5、将得到的Mask图像与原始图像进行逻辑与操作得到最后的气管分割结果图像。
代码：
def skeletonSegment(pathDicom): sitk_src = dicomseriesReader(pathDicom) # 1 sitk_seg = BinaryThreshold(sitk_src, lowervalue=100, uppervalue=3000) sitk.WriteImage(sitk_seg, 'step1.mha') # 2 sitk_open = MorphologicalOperation(sitk_seg, kernelsize=2, name='open') sitk_open = GetLargestConnectedCompont(sitk_open) sitk.WriteImage(sitk_open, 'step2.mha') # 3 array_open = sitk.GetArrayFromImage(sitk_open) array_seg = sitk.GetArrayFromImage(sitk_seg) array_mask = array_seg - array_open sitk_mask = sitk.GetImageFromArray(array_mask) sitk_mask.SetDirection(sitk_seg.GetDirection()) sitk_mask.SetSpacing(sitk_seg.GetSpacing()) sitk_mask.SetOrigin(sitk_seg.GetOrigin()) sitk.WriteImage(sitk_mask, 'step3.mha') # 4 skeleton_mask = GetLargestConnectedCompont(sitk_mask) sitk.WriteImage(skeleton_mask, 'step4.mha') # 5 sitk_skeleton = GetMaskImage(sitk_src, skeleton_mask, replacevalue=-1500) sitk.WriteImage(sitk_skeleton, 'step5.mha') return sitk_skeleton2、气管分割代码详解 2.1、采用区域生长方法来提取气管，一般空气的HU在-1024到-900左右。
2.2、 将得到的Mask图像与原始图像进行逻辑与操作得到最后的气管分割结果图像。
代码：
def tracheaSegment(pathDicom): sitk_src = dicomseriesReader(pathDicom) lstSeeds = [] seed1 = [246, 268, 835] seed2 = [241, 287, 798] seed3 = [272, 277, 798] lstSeeds.append(seed1) lstSeeds.append(seed2) lstSeeds.append(seed3) sitk_mask = RegionGrowThreshold(sitk_src, lstSeeds, -1024, -900) sitk.WriteImage(sitk_mask, 'tracheamask.mha') sitk_trachea = GetMaskImage(sitk_src, sitk_mask, replacevalue=-1500) return sitk_trachea3、肺组织分割代码详解 3.1、 首先采用固定阈值区分肺组织和身体其它组织，设置固定阈值为-300。
3.2、 然后采用水漫金山法将肺组织和人体其它组织填充成一个区域。
3.3、 然后对固定阈值结果和水漫金山法结果先进行逻辑异或操作，然后再对其结果进行取反。
3.4、 对肺部内部区域进行填小洞操作。
3.5、采用连通域分析，去除小于最大物体一定比例的小目标物体。
3.6、采用区域生长得到气管分割结果。
3.7、 将第3.5步结果减去第3.6步结果。
3.8、 首先采用形态学开操作去除多余的边界区域，然后再采用形态学闭操作平滑肺组织轮廓，并填充肺部中中型空洞区域。
3.9、将得到的Mask图像与原始图像进行逻辑与操作得到最后的肺部分割结果图像。
代码：
def lungSegment(pathDicom): sitk_src = dicomseriesReader(pathDicom) # 1 sitk_seg = BinaryThreshold(sitk_src, lowervalue=-300, uppervalue=2000) sitk.WriteImage(sitk_seg, 'step1.mha') # 2 sitk_floodfilled = FloodFilled(sitk_seg) sitk.WriteImage(sitk_floodfilled, 'step2.mha') # 3 sitk_xorop = sitk.XorImageFilter() sitk_mask1 = sitk_xorop.Execute(sitk_seg, sitk_floodfilled) sitk_notop = sitk.NotImageFilter() sitk_mask2 = sitk_notop.Execute(sitk_mask1) sitk.WriteImage(sitk_mask2, 'step3.mha') # 4 sitk_mask3 = FillHole(sitk_mask2) sitk.WriteImage(sitk_mask3, 'step4.mha') # 5 sitk_mask4 = RemoveSmallConnectedCompont(sitk_mask3, 0.2) sitk.WriteImage(sitk_mask4, 'step5.mha') # 6 segtrachea lstSeeds = [] seed1 = [246, 268, 835] seed2 = [241, 287, 798] seed3 = [272, 277, 798] lstSeeds.append(seed1) lstSeeds.append(seed2) lstSeeds.append(seed3) sitk_tracheamask = RegionGrowThreshold(sitk_src, lstSeeds, -1024, -900) sitk.WriteImage(sitk_tracheamask, 'step6.mha') # 7 lung reduce trachea array_tracheamask = sitk.GetArrayFromImage(sitk_tracheamask) array_mask4 = sitk.GetArrayFromImage(sitk_mask4) array_mask4 = array_mask4 - array_tracheamask sitk_mask4 = sitk.GetImageFromArray(array_mask4) sitk_mask4.SetDirection(sitk_tracheamask.GetDirection()) sitk_mask4.SetSpacing(sitk_tracheamask.GetSpacing()) sitk_mask4.SetOrigin(sitk_tracheamask.GetOrigin()) sitk.WriteImage(sitk_mask4, 'step7.mha') # 8 sitk_mask4 = MorphologicalOperation(sitk_mask4, kernelsize=3, name='open') sitk_mask5 = MorphologicalOperation(sitk_mask4, kernelsize=9, name='close') sitk.WriteImage(sitk_mask5, 'step8.mha') # 9 sitk_lung = GetMaskImage(sitk_src, sitk_mask5, replacevalue=-1500) sitk.WriteImage(sitk_lung, 'step9.mha') return sitk_lung4、血管分割代码详解
4.1、采用itk的多尺度海森矩阵血管增强
4.2、 将增强后的血管归一化到0-255范围
4.3、 对归一化的血管进行固定阈值分割
代码：
def vessleSegment(pathDicom): sigma_minimum = 0.2 sigma_maximum = 3. number_of_sigma_steps = 8 lowerThreshold = 40 output_image = 'vessel.mha' input_image = itk.imread(pathDicom, itk.F) # 1 ImageType = type(input_image) Dimension = input_image.GetImageDimension() HessianPixelType = itk.SymmetricSecondRankTensor[itk.D, Dimension] HessianImageType = itk.Image[HessianPixelType, Dimension] objectness_filter = itk.HessianToObjectnessMeasureImageFilter[HessianImageType, ImageType].New() objectness_filter.SetBrightObject(True) objectness_filter.SetScaleObjectnessMeasure(True) objectness_filter.SetAlpha(0.5) objectness_filter.SetBeta(1.0) objectness_filter.SetGamma(5.0) multi_scale_filter = itk.MultiScaleHessianBasedMeasureImageFilter[ImageType, HessianImageType, ImageType].New() multi_scale_filter.SetInput(input_image) multi_scale_filter.SetHessianToMeasureFilter(objectness_filter) multi_scale_filter.SetSigmaStepMethodToLogarithmic() multi_scale_filter.SetSigmaMinimum(sigma_minimum) multi_scale_filter.SetSigmaMaximum(sigma_maximum) multi_scale_filter.SetNumberOfSigmaSteps(number_of_sigma_steps) itk.imwrite(multi_scale_filter.GetOutput(), "step1.mha") # 2 OutputPixelType = itk.UC OutputImageType = itk.Image[OutputPixelType, Dimension] rescale_filter = itk.RescaleIntensityImageFilter[ImageType, OutputImageType].New() rescale_filter.SetInput(multi_scale_filter) itk.imwrite(rescale_filter.GetOutput(), "step2.mha") # 3 thresholdFilter = itk.BinaryThresholdImageFilter[OutputImageType, OutputImageType].New() thresholdFilter.SetInput(rescale_filter.GetOutput()) thresholdFilter.SetLowerThreshold(lowerThreshold) thresholdFilter.SetUpperThreshold(255) thresholdFilter.SetOutsideValue(0) thresholdFilter.SetInsideValue(255) itk.imwrite(thresholdFilter.GetOutput(), "step3.mha")为了方便大家更高效地学习，我将代码进行了整理并更新到github上，点击原文链接即可访问。测试数据集可以在百度云盘上下载： https://pan.baidu.com/s/1caOdrgmtBkbSCQvM7Y0t6Q 密码： j6s6 。 如果大家觉得这个项目还不错，希望大家给个Star并Fork，可以让更多的人学习。如果有任何问题，随时给我留言我会及时回复的。
全部评论 (0)
还没有任何评论哟~