碼上歡樂
相關內容    簡體    繁體

[轉]常見醫療掃描圖像處理步驟

本文轉載自   查看原文   2018-12-06 18:38   1684    DICOM

DICOM是醫學圖像中標准文件，這些文件包含了諸多的元數據信息（比如像素尺寸，每個維度的一像素代表真實世界里的長度）。此處以kaggle Data Science Bowl 數據集為例。

data-science-bowl-2017。數據列表如下:

dicom格式的圖像

后綴為 .dcm。

每個病人的一次掃描CT(scan)可能有幾十到一百多個dcm數據文件(slices)。可以使用 python的dicom包讀取，讀取示例代碼如下:

  1. dicom.read_file('/data/lung_competition/stage1/7050f8141e92fa42fd9c471a8b2f50ce/498d16aa2222d76cae1da144ddc59a13.dcm')

其pixl_array包含了真實數據。

  1. slices = [dicom.read_file(os.path.join(folder_name,filename)) for filename in os.listdir(folder_name)]
  2. slices = np.stack([s.pixel_array for s in slices])

mhd格式是另外一種數據格式，來源於(LUNA2016)[https://luna16.grand-challenge.org/data/]。每個病人一個mhd文件和一個同名的raw文件。如下:

dicom格式的圖像

一個mhd通常有幾百兆，對應的raw文件只有1kb。mhd文件需要借助python的SimpleITK包來處理。SimpleITK 示例代碼如下:

  1. import SimpleITK as sitk
  2. itk_img = sitk.ReadImage(img_file)
  3. img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
  4. num_z, height, width = img_array.shape #heightXwidth constitute the transverse plane
  5. origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)
  6. spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)

需要注意的是，SimpleITK的img_array的數組不是直接的像素值，而是相對於CT掃描中原點位置的差值，需要做進一步轉換。轉換步驟參考 SimpleITK圖像轉換

一個開源免費的查看軟件 mango

dicom格式的圖像

首先，需要明白的是醫學掃描圖像(scan)其實是三維圖像，使用代碼讀取之后開源查看不同的切面的切片(slices)，可以從不同軸切割

dicom格式的圖像

如下圖展示了一個病人CT掃描圖中，其中部分切片slices

dicom格式的圖像

其次，CT掃描圖是包含了所有組織的，如果直接去看，看不到任何有用信息。需要做一些預處理,預處理中一個重要的概念是放射劑量，衡量單位為HU(Hounsfield Unit)，下表是不同放射劑量對應的組織器官

substance HU
空氣 -1000
肺 -500
脂肪 -100到-50
水 0
CSF 15
腎 30
血液 +30到+45
肌肉 +10到+40
灰質 +37到+45
白質 +20到+30
Liver +40到+60
軟組織,constrast +100到+300
骨頭 +700(軟質骨)到+3000(皮質骨)
  1. Hounsfield Unit = pixel_value * rescale_slope + rescale_intercept

一般情況rescale slope = 1, intercept = -1024。

上表中肺部組織的HU數值為-500,但通常是大於這個值，比如-320、-400。挑選出這些區域，然后做其他變換抽取出肺部像素點。

  1. # -*- coding:utf-8 -*-
  2. '''
  3. this script is used for basic process of lung 2017 in Data Science Bowl
  4. '''
  6. import glob
  7. import os
  8. import pandas as pd
  9. import SimpleITK as sitk
  11. import numpy as np # linear algebra
  12. import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
  13. import skimage, os
  14. from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
  15. from skimage.measure import label,regionprops, perimeter
  16. from skimage.morphology import binary_dilation, binary_opening
  17. from skimage.filters import roberts, sobel
  18. from skimage import measure, feature
  19. from skimage.segmentation import clear_border
  20. from skimage import data
  21. from scipy import ndimage as ndi
  22. import matplotlib
  23. #matplotlib.use('Agg')
  24. import matplotlib.pyplot as plt
  25. from mpl_toolkits.mplot3d.art3d import Poly3DCollection
  26. import dicom
  27. import scipy.misc
  28. import numpy as np

如下代碼是載入一個掃描面，包含了多個切片(slices)，我們僅簡化的將其存儲為python列表。數據集中每個目錄都是一個掃描面集（一個病人）。有個元數據域丟失，即Z軸方向上的像素尺寸，也即切片的厚度 。所幸，我們可以用其他值推測出來，並加入到元數據中。

  1. # Load the scans in given folder path
  2. def load_scan(path):
  3. slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
  4. slices.sort(key = lambda x: int(x.ImagePositionPatient[2]))
  5. try:
  6. slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
  7. except:
  8. slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
  10. for s in slices:
  11. s.SliceThickness = slice_thickness
  13. return slices

首先去除灰度值為 -2000的pixl_array，CT掃描邊界之外的灰度值固定為-2000(dicom和mhd都是這個值)。第一步是設定這些值為0，當前對應為空氣（值為0）

回到HU單元，乘以rescale比率並加上intercept(存儲在掃描面的元數據中)。

  1. def get_pixels_hu(slices):
  2. image = np.stack([s.pixel_array for s in slices])
  3. # Convert to int16 (from sometimes int16),
  4. # should be possible as values should always be low enough (<32k)
  5. image = image.astype(np.int16)
  6. # Set outside-of-scan pixels to 0
  7. # The intercept is usually -1024, so air is approximately 0
  8. image[image == -2000] = 0
  10. # Convert to Hounsfield units (HU)
  11. for slice_number in range(len(slices)):
  13. intercept = slices[slice_number].RescaleIntercept
  14. slope = slices[slice_number].RescaleSlope
  16. if slope != 1:
  17. image[slice_number] = slope * image[slice_number].astype(np.float64)
  18. image[slice_number] = image[slice_number].astype(np.int16)
  20. image[slice_number] += np.int16(intercept)
  22. return np.array(image, dtype=np.int16)

可以查看病人的掃描HU值分布情況

  1. first_patient = load_scan(INPUT_FOLDER + patients[0])
  2. first_patient_pixels = get_pixels_hu(first_patient)
  3. plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
  4. plt.xlabel("Hounsfield Units (HU)")
  5. plt.ylabel("Frequency")
  6. plt.show()

dicom格式的圖像

不同掃描面的像素尺寸、粗細粒度是不同的。這不利於我們進行CNN任務，我們可以使用同構采樣。

一個掃描面的像素區間可能是[2.5,0.5,0.5],即切片之間的距離為2.5mm。可能另外一個掃描面的范圍是[1.5,0.725,0.725]。這可能不利於自動分析。

常見的處理方法是從全數據集中以固定的同構分辨率重新采樣，將所有的東西采樣為1mmx1mmx1mm像素。

  1. def resample(image, scan, new_spacing=[1,1,1]):
  2. # Determine current pixel spacing
  3. spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))
  4. spacing = np.array(list(spacing))
  5. resize_factor = spacing / new_spacing
  6. new_real_shape = image.shape * resize_factor
  7. new_shape = np.round(new_real_shape)
  8. real_resize_factor = new_shape / image.shape
  9. new_spacing = spacing / real_resize_factor
  11. image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
  13. return image, new_spacing

現在重新取樣病人的像素，將其映射到一個同構分辨率 1mm x1mm x1mm。

  1. pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])

使用matplotlib輸出肺部掃描的3D圖像方法。可能需要一兩分鍾。

  1. def plot_3d(image, threshold=-300):
  3. # Position the scan upright,
  4. # so the head of the patient would be at the top facing the camera
  5. p = image.transpose(2,1,0)
  7. verts, faces = measure.marching_cubes(p, threshold)
  8. fig = plt.figure(figsize=(10, 10))
  9. ax = fig.add_subplot(111, projection='3d')
  10. # Fancy indexing: `verts[faces]` to generate a collection of triangles
  11. mesh = Poly3DCollection(verts[faces], alpha=0.1)
  12. face_color = [0.5, 0.5, 1]
  13. mesh.set_facecolor(face_color)
  14. ax.add_collection3d(mesh)
  15. ax.set_xlim(0, p.shape[0])
  16. ax.set_ylim(0, p.shape[1])
  17. ax.set_zlim(0, p.shape[2])
  18. plt.show()

打印函數有個閾值參數，來打印特定的結構，比如tissue或者骨頭。400是一個僅僅打印骨頭的閾值(HU對照表)

  1. plot_3d(pix_resampled, 400)

dicom格式的圖像

  1. def plot_ct_scan(scan):
  2. '''
  3. plot a few more images of the slices
  4. :param scan:
  5. :return:
  6. '''
  7. f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(50, 50))
  8. for i in range(0, scan.shape[0], 5):
  9. plots[int(i / 20), int((i % 20) / 5)].axis('off')
  10. plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)

此方法的效果示例如下:

dicom格式的圖像

下面的代碼使用了pythonde 的圖像形態學操作。具體可以參考python高級形態學操作

  1. def get_segmented_lungs(im, plot=False):
  3. '''
  4. This funtion segments the lungs from the given 2D slice.
  5. '''
  6. if plot == True:
  7. f, plots = plt.subplots(8, 1, figsize=(5, 40))
  8. '''
  9. Step 1: Convert into a binary image.
  10. '''
  11. binary = im < 604
  12. if plot == True:
  13. plots[0].axis('off')
  14. plots[0].set_title('binary image')
  15. plots[0].imshow(binary, cmap=plt.cm.bone)
  17. '''
  18. Step 2: Remove the blobs connected to the border of the image.
  19. '''
  20. cleared = clear_border(binary)
  21. if plot == True:
  22. plots[1].axis('off')
  23. plots[1].set_title('after clear border')
  24. plots[1].imshow(cleared, cmap=plt.cm.bone)
  26. '''
  27. Step 3: Label the image.
  28. '''
  29. label_image = label(cleared)
  30. if plot == True:
  31. plots[2].axis('off')
  32. plots[2].set_title('found all connective graph')
  33. plots[2].imshow(label_image, cmap=plt.cm.bone)
  34. '''
  35. Step 4: Keep the labels with 2 largest areas.
  36. '''
  37. areas = [r.area for r in regionprops(label_image)]
  38. areas.sort()
  39. if len(areas) > 2:
  40. for region in regionprops(label_image):
  41. if region.area < areas[-2]:
  42. for coordinates in region.coords:
  43. label_image[coordinates[0], coordinates[1]] = 0
  44. binary = label_image > 0
  45. if plot == True:
  46. plots[3].axis('off')
  47. plots[3].set_title(' Keep the labels with 2 largest areas')
  48. plots[3].imshow(binary, cmap=plt.cm.bone)
  49. '''
  50. Step 5: Erosion operation with a disk of radius 2. This operation is
  51. seperate the lung nodules attached to the blood vessels.
  52. '''
  53. selem = disk(2)
  54. binary = binary_erosion(binary, selem)
  55. if plot == True:
  56. plots[4].axis('off')
  57. plots[4].set_title('seperate the lung nodules attached to the blood vessels')
  58. plots[4].imshow(binary, cmap=plt.cm.bone)
  59. '''
  60. Step 6: Closure operation with a disk of radius 10. This operation is
  61. to keep nodules attached to the lung wall.
  62. '''
  63. selem = disk(10)
  64. binary = binary_closing(binary, selem)
  65. if plot == True:
  66. plots[5].axis('off')
  67. plots[5].set_title('keep nodules attached to the lung wall')
  68. plots[5].imshow(binary, cmap=plt.cm.bone)
  69. '''
  70. Step 7: Fill in the small holes inside the binary mask of lungs.
  71. '''
  72. edges = roberts(binary)
  73. binary = ndi.binary_fill_holes(edges)
  74. if plot == True:
  75. plots[6].axis('off')
  76. plots[6].set_title('Fill in the small holes inside the binary mask of lungs')
  77. plots[6].imshow(binary, cmap=plt.cm.bone)
  78. '''
  79. Step 8: Superimpose the binary mask on the input image.
  80. '''
  81. get_high_vals = binary == 0
  82. im[get_high_vals] = 0
  83. if plot == True:
  84. plots[7].axis('off')
  85. plots[7].set_title('Superimpose the binary mask on the input image')
  86. plots[7].imshow(im, cmap=plt.cm.bone)
  88. return im

此方法每個步驟對圖像做不同的處理，依次為二值化、清除邊界、連通區域標記、腐蝕操作、閉合運算、孔洞填充、效果如下:

dicom格式的圖像

為了減少有問題的空間，我們可以分割肺部圖像（有時候是附近的組織）。這包含一些步驟，包括區域增長和形態運算，此時，我們只分析相連組件。

步驟如下：

  • 閾值圖像（-320HU是個極佳的閾值，但是此方法中不是必要）

  • 處理相連的組件，以決定當前患者的空氣的標簽，以1填充這些二值圖像

  • 可選：當前掃描的每個軸上的切片，選定最大固態連接的組織（當前患者的肉體和空氣），並且其他的為0。以掩碼的方式填充肺部結構。

  • 只保留最大的氣袋（人類軀體內到處都有氣袋）

  1. def largest_label_volume(im, bg=-1):
  2. vals, counts = np.unique(im, return_counts=True)
  3. counts = counts[vals != bg]
  4. vals = vals[vals != bg]
  5. if len(counts) > 0:
  6. return vals[np.argmax(counts)]
  7. else:
  8. return None
  9. def segment_lung_mask(image, fill_lung_structures=True):
  11. # not actually binary, but 1 and 2.
  12. # 0 is treated as background, which we do not want
  13. binary_image = np.array(image > -320, dtype=np.int8)+1
  14. labels = measure.label(binary_image)
  16. # Pick the pixel in the very corner to determine which label is air.
  17. # Improvement: Pick multiple background labels from around the patient
  18. # More resistant to "trays" on which the patient lays cutting the air
  19. # around the person in half
  20. background_label = labels[0,0,0]
  22. #Fill the air around the person
  23. binary_image[background_label == labels] = 2
  26. # Method of filling the lung structures (that is superior to something like
  27. # morphological closing)
  28. if fill_lung_structures:
  29. # For every slice we determine the largest solid structure
  30. for i, axial_slice in enumerate(binary_image):
  31. axial_slice = axial_slice - 1
  32. labeling = measure.label(axial_slice)
  33. l_max = largest_label_volume(labeling, bg=0)
  35. if l_max is not None: #This slice contains some lung
  36. binary_image[i][labeling != l_max] = 1
  38. binary_image -= 1 #Make the image actual binary
  39. binary_image = 1-binary_image # Invert it, lungs are now 1
  41. # Remove other air pockets insided body
  42. labels = measure.label(binary_image, background=0)
  43. l_max = largest_label_volume(labels, bg=0)
  44. if l_max is not None: # There are air pockets
  45. binary_image[labels != l_max] = 0
  47. return binary_image

查看切割效果

  1. segmented_lungs = segment_lung_mask(pix_resampled, False)
  2. segmented_lungs_fill = segment_lung_mask(pix_resampled, True)
  3. plot_3d(segmented_lungs, 0)

dicom格式的圖像

我們可以將肺內的結構也包含進來（結節是固體），不僅僅只是肺部內的空氣

  1. plot_3d(segmented_lungs_fill, 0)

dicom格式的圖像

使用mask時，要注意首先進行形態擴充(python的skimage的skimage.morphology)操作（即使用圓形kernel，結節是球體），參考 python形態操作。這會在所有方向（維度）上擴充mask。僅僅肺部的空氣+結構將不會包含所有結節，事實上有可能遺漏黏在肺部一側的結節（這會經常出現，所以建議最好是擴充mask）。

歸一化處理

當前的值范圍是[-1024,2000]。而任意大於400的值並不是處理肺結節需要考慮，因為它們都是不同反射密度下的骨頭。LUNA16競賽中常用來做歸一化處理的閾值集是-1000和400.以下代碼

歸一化

  1. MIN_BOUND = -1000.0
  2. MAX_BOUND = 400.0
  4. def normalize(image):
  5. image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
  6. image[image>1] = 1.
  7. image[image<0] = 0.
  8. return image

0值中心化

簡單來說就是所有像素值減去均值。LUNA16競賽中的均值大約是0.25.

不要對每一張圖像做零值中心化（此處像是在kernel中完成的）CT掃描器返回的是校准后的精確HU計量。不會出現普通圖像中會出現某些圖像低對比度和明亮度的情況

  1. PIXEL_MEAN = 0.25
  3. def zero_center(image):
  4. image = image - PIXEL_MEAN
  5. return image

歸一化和零值中心化的操作主要是為了后續訓練網絡，零值中心化是網絡收斂的關鍵。

mhd的數據只是格式與dicom不一樣，其實質包含的都是病人掃描。處理MHD需要借助SimpleIKT這個包，處理思路詳情可以參考Data Science Bowl2017的toturail Data Science Bowl 2017。需要注意的是MHD格式的數據沒有HU值，它的值域范圍與dicom很不同。

我們以LUNA2016年的數據處理流程為例。參考代碼為 LUNA2016數據切割

  1. import SimpleITK as sitk
  2. import numpy as np
  3. import csv
  4. from glob import glob
  5. import pandas as pd
  6. file_list=glob(luna_subset_path+"*.mhd")
  7. #####################
  8. #
  9. # Helper function to get rows in data frame associated
  10. # with each file
  11. def get_filename(case):
  12. global file_list
  13. for f in file_list:
  14. if case in f:
  15. return(f)
  16. #
  17. # The locations of the nodes
  18. df_node = pd.read_csv(luna_path+"annotations.csv")
  19. df_node["file"] = df_node["seriesuid"].apply(get_filename)
  20. df_node = df_node.dropna()
  21. #####
  22. #
  23. # Looping over the image files
  24. #
  25. fcount = 0
  26. for img_file in file_list:
  27. print "Getting mask for image file %s" % img_file.replace(luna_subset_path,"")
  28. mini_df = df_node[df_node["file"]==img_file] #get all nodules associate with file
  29. if len(mini_df)>0: # some files may not have a nodule--skipping those
  30. biggest_node = np.argsort(mini_df["diameter_mm"].values)[-1] # just using the biggest node
  31. node_x = mini_df["coordX"].values[biggest_node]
  32. node_y = mini_df["coordY"].values[biggest_node]
  33. node_z = mini_df["coordZ"].values[biggest_node]
  34. diam = mini_df["diameter_mm"].values[biggest_node]

一直在尋找MHD格式數據的處理方法，對於dicom格式的CT有很多論文根據其HU值域可以輕易地分割肺、骨頭、血液等，但是對於MHD沒有這樣的參考。從LUNA16論壇得到的解釋是，LUNA16的MHD數據已經轉換為HU值了，不需要再使用slope和intercept來做rescale變換了。此論壇主題下，有人提出MHD格式沒有提供pixel spacing(mm) 和 slice thickness(mm) ，而標准文件annotation.csv文件中結節的半徑和坐標都是mm單位，最后確認的是MHD格式文件中只保留了體素尺寸以及坐標原點位置，沒有保存slice thickness。即，dicom才是原始數據格式。

MHD值的坐標體系是體素，以mm為單位（dicom的值是GV灰度值）。結節的位置是CT scanner坐標軸里面相對原點的mm值，需要將其轉換到真實坐標軸位置，可以使用SimpleITK包中的 GetOrigin() ` GetSpacing()`。圖像數據是以512x512數組的形式給出的。

坐標變換如下：

dicom格式的圖像

相應的代碼處理如下:

  1. itk_img = sitk.ReadImage(img_file)
  2. img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
  3. center = np.array([node_x,node_y,node_z]) # nodule center
  4. origin = np.array(itk_img.GetOrigin()) # x,y,z Origin in world coordinates (mm)
  5. spacing = np.array(itk_img.GetSpacing()) # spacing of voxels in world coor. (mm)
  6. v_center =np.rint((center-origin)/spacing) # nodule center in voxel space (still x,y,z ordering)

在LUNA16的標注CSV文件中標注了結節中心的X,Y,Z軸坐標，但是實際取值的時候取的是Z軸最后三層的數組(img_array)。

下述代碼只提取了包含結節的最后三個slice的數據，代碼參考自 LUNA_mask_extraction.py

  1. i = 0
  2. for i_z in range(int(v_center[2])-1,int(v_center[2])+2):
  3. mask = make_mask(center,diam,i_z*spacing[2]+origin[2],width,height,spacing,origin)
  4. masks[i] = mask
  5. imgs[i] = matrix2int16(img_array[i_z])
  6. i+=1
  7. np.save(output_path+"images_%d.npy" % (fcount) ,imgs)
  8. np.save(output_path+"masks_%d.npy" % (fcount) ,masks)

以下代碼用於查看原始CT和結節mask。其實就是用matplotlib打印上一步存儲的npy文件。

  1. import matplotlib.pyplot as plt
  2. imgs = np.load(output_path+'images_0.npy')
  3. masks = np.load(output_path+'masks_0.npy')
  4. for i in range(len(imgs)):
  5. print "image %d" % i
  6. fig,ax = plt.subplots(2,2,figsize=[8,8])
  7. ax[0,0].imshow(imgs[i],cmap='gray')
  8. ax[0,1].imshow(masks[i],cmap='gray')
  9. ax[1,0].imshow(imgs[i]*masks[i],cmap='gray')
  10. plt.show()
  11. raw_input("hit enter to cont : ")

示例結節和mask

接下來的處理和DICOM格式數據差不多，腐蝕膨脹、連通區域標記等。

參考信息：

灰度值是pixel value經過重重LUT轉換得到的用來進行顯示的值，而這個轉換過程是不可逆的，也就是說，灰度值無法轉換為ct值。只能根據窗寬窗位得到一個大概的范圍。 pixel value經過modality lut得到Hu，但是懷疑pixelvalue的讀取出了問題。dicom文件中存在（0028，0106）（0028，0107）兩個tag，分別是最大最小pixel value，可以用來檢驗你讀取的pixel value 矩陣是否正確。

LUT全稱look up table，實際上就是一張像素灰度值的映射表，它將實際采樣到的像素灰度值經過一定的變換如閾值、反轉、二值化、對比度調整、線性變換等，變成了另外一 個與之對應的灰度值，這樣可以起到突出圖像的有用信息，增強圖像的光對比度的作用。


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 掃描儀掃描圖像 從 TWAIN 設備中掃描圖像 matlab中圖像處理常見用法 了解GDAL圖像處理（轉） python數字圖像處理（一）圖像的常見操作 圖像處理的基本知識（轉） [轉]用Python做圖像處理 【轉】OpenCV圖像處理 圖像的點運算 ( 灰度直方圖 ) opencv圖像處理之常見濾波器 圖像處理常見面試題
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM