三維區(qū)域生長from operator import eq import numpy as np import pydicom as dicom from skimage.measure import label from skimage.measure import regionprops from skimage.morphology import thin import os import SimpleITK as sptk import matplotlib.pyplot as plt import matplotlib.image as mpimg import cv2def regiongrowth(orig_image, seed_matrix, threshold):"""3-D區(qū)域生長算法# Arguments：orig_image：待分割圖像seed_matrix：種子矩陣threshold：生長條件，即待分割圖中的像素值與目標像素值之間的可承受的差值# Returnsseg：最終分割后的masknum_region：分割區(qū)域的個數(shù)# Exampleregiongrow(orig_image, seed_matrix, 50)"""orig_image = np.array(orig_image, dtype=np.float32) # 將圖像數(shù)值格式轉為floatmarkers = seed_matrix#thin(seed_matrix) # 將種子所在區(qū)域縮小為一個點，后面限定區(qū)域時使用??????coordinate_nonzero = markers.nonzero() # 當種子為矩陣時，獲取種子點的坐標值，返回兩個數(shù)組，分別是非零元素所在行和所在列coordinate_nonzero = np.transpose(coordinate_nonzero) # 獲取種子點的坐標值的轉置seed_value = []for coordinate_i in coordinate_nonzero:seed_value.append(orig_image[coordinate_i[0], coordinate_i[1], coordinate_i[2]]) # 獲取種子點的灰度值size_x, size_y, size_z = orig_image.shapeseg_image = np.array(np.zeros((size_x, size_y, size_z)), dtype=bool) # 初始化一個bool格式的矩陣儲存分割結果for i, each_seed_value in enumerate(seed_value):all_satisfied = abs(orig_image - each_seed_value) <= threshold # 所有滿足單個種子點閾值條件的像素點seg_image = seg_image | all_satisfied # 所有滿足多有種子點閾值條件的像素點label_seg_image = label(seg_image) # 將所有連通區(qū)域打label？？？？？？？？area_region = regionprops(label_seg_image) # 統(tǒng)計被標記的區(qū)域的面積分布，返回值為顯示區(qū)域總數(shù)seg = np.array(np.zeros((size_x, size_y, size_z)), dtype=np.float32) # 初始化一個float格式的矩陣儲存分割結果for masker in coordinate_nonzero: # 提取maskers所在的連通區(qū)域for each_area_region in area_region:for coord in each_area_region.coords:if sum(eq(coord, masker)) == 2: # 若連通區(qū)域中包含masker，則保留此連通區(qū)域for idx in each_area_region.coords:seg[idx[0], idx[1], idx[2]] = 1breakseg = label(seg) # 將有效的label進行標記num_region = seg.max() # 有效連通label的個數(shù)return seg, num_regionif __name__ == "__main__":#one_slices = dicom.dcmread(r'D:\LIDC\dicm\LIDC-IDRI\LIDC-IDRI-0001\01-01-2000-30178\3000566-03192\000100.dcm', force=True)#orig_image = one_slices.pixel_array#orig_image= np.array(cv2.imread("D:\LIDC\cup.jpg"))PathDicom = "D:/LIDC/dicm/LIDC-IDRI/LIDC-IDRI-0001/01-01-2000-30178/3000566-03192" # 與python文件同一個目錄下的文件夾lstFilesDCM = []for dirName, subdirList, fileList in os.walk(PathDicom):for filename in fileList:if ".dcm" in filename.lower(): # 判斷文件是否為dicom文件lstFilesDCM.append(os.path.join(dirName, filename)) # 加入到列表中array3D = np.zeros((1, 512, 512))for dicmImage in lstFilesDCM:image = sptk.ReadImage(dicmImage)image_array2D = np.squeeze(sptk.GetArrayFromImage(image))image_array = np.append(array3D, image_array2D)dim = array3D.shapearray3D = image_array.reshape(dim[0] + 1, dim[1], dim[2])array3D = np.delete(array3D, 0, 0)array3D=array3D.transpose((1,2,0))#print(orig_image.shape)plt.figure()plt.imshow(array3D[:,:,2])plt.show()#orig_image = (orig_image - np.min(orig_image))/(np.max(orig_image) - np.min(orig_image)) * 255#orig_image = np.array(orig_image, dtype=np.uint8)#plt.figure()#plt.imshow(orig_image)#plt.show()# dst_image = np.zeros((orig_image.shape[0], orig_image.shape[1]))# dst_image = cv2.adaptiveThreshold(orig_image, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 511, 2)# dst_image = dst_image / 255 * (np.max(orig_image) - np.min(orig_image)) + np.min(orig_image)# plt.figure()# plt.subplot(1, 2, 1)# plt.imshow(orig_image)# plt.subplot(1, 2, 2)# plt.imshow(dst_image)# plt.show()size_x, size_y, size_z = array3D.shapeorig_image = np.array(array3D, dtype=np.float32)orig_image_min = orig_image.min()orig_image = orig_image - orig_image_minorig_image_max = orig_image.max()orig_image = orig_image / orig_image_max * 255orig_image = np.array(orig_image, dtype=np.uint8)seed_matrix = np.zeros((size_x, size_y, size_z))# 指定想要分割的像素的坐標，按照dicom的坐標系coordinate1_x = 260coordinate1_y = 158coordinate1_z = 2coordinate2_x = 260coordinate2_y = 370coordinate2_z = 2#coordinate1_x = 138#coordinate1_y = 252#coordinate2_x = 346#coordinate2_y = 252#coordinate1_x = 226#coordinate1_y = 290#coordinate2_x = 308#coordinate2_y = 290# 矩陣按行號、列號，與dicom正好相反seed_matrix[coordinate1_x, coordinate1_y, coordinate1_z] = orig_image[coordinate1_x, coordinate1_y, coordinate1_z]seed_matrix[coordinate2_x, coordinate2_y, coordinate2_z] = orig_image[coordinate2_x, coordinate2_y, coordinate2_z]seg_image, num_region = regiongrowth(orig_image, seed_matrix, 3) # 需指定閾值print("返回的矩陣維數(shù)：", seg_image.shape)print("連通區(qū)域數(shù)目：", num_region)plt.figure()fig1=plt.subplot(1, 2, 1)fig2 = plt.subplot(1, 2, 2)plt.sca(fig1)plt.imshow(orig_image[:,:,2])plt.sca(fig2)plt.imshow(seg_image[:,:,2])plt.show()

總結

以上是生活随笔為你收集整理的regionGrowth3D的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。