Trong hướng dẫn này, tôi sẽ hướng dẫn từng bước về cách áp dụng các phương pháp phân cụm thống kê, thuật toán đồ họa máy tính và kỹ thuật xử lý hình ảnh cho hình ảnh y tế để giúp hiểu và trực quan hóa dữ liệu ở cả 2D và 3D [tất cả mã . ]
Bước 1. Cài đặt các gói cần thiết
Hai cách phổ biến để cài đặt các gói là sử dụng conda hoặc pip. Cá nhân tôi nghĩ. tại sao không phải cả hai?
- Trăn 3. 6 [Tôi tìm thấy các phiên bản lớn hơn 3. 6 gặp sự cố tương thích với pydicom]
- pydicom > để đọc tệp DICOM
- numpy > cho các ứng dụng điện toán
- scipy > cũng dành cho các ứng dụng điện toán
- Skipage > cho các chức năng xử lý hình ảnh
- matplotlib -> để trực quan hóa
- plotly > để hiển thị tương tác [xem ghi chú bên dưới trước khi cài đặt]
- ipywidget > cho các lô tương tác
*GHI CHÚ. Để sử dụng cốt truyện, bạn sẽ cần phải. [1. ] Cài đặt tất cả các gói/tiện ích mở rộng cần thiết theo các hướng dẫn sau. [2. ] đăng ký tài khoản miễn phí và kích hoạt khóa của bạn theo các hướng dẫn này [chỉ đọc phần 'Khởi tạo cho âm mưu trực tuyến']
Bước 2. Đang tải dữ liệu DICOMchúc mừng. Bạn đã vượt qua phần tồi tệ nhất… cài đặt các gói. Đối với phần còn lại của hướng dẫn, tôi sẽ sử dụng máy quét CT ngực và sẽ chạy các tập lệnh của mình trên máy tính xách tay jupyter
Trước tiên, hãy nhập các gói cần thiết
# common packages
import numpy as np
import os
import copy
from math import *
import matplotlib.pyplot as plt
from functools import reduce# reading in dicom files
import pydicom# skimage image processing packages
from skimage import measure, morphology
from skimage.morphology import ball, binary_closing
from skimage.measure import label, regionprops# scipy linear algebra functions
from scipy.linalg import norm
import scipy.ndimage# ipywidgets for some interactive plots
from ipywidgets.widgets import *
import ipywidgets as widgets# plotly 3D interactive graphs
import plotly
from plotly.graph_objs import *
import chart_studio.plotly as py# set plotly credentials here
# this allows you to send results to your account plotly.tools.set_credentials_file[username=your_username, api_key=your_key]
Tiếp theo, hãy tải hình ảnh DICOM của bạn
def load_scan[path]:
slices = [pydicom.dcmread[path + ‘/’ + s] for s in
os.listdir[path]]
slices = [s for s in slices if ‘SliceLocation’ in s]
slices.sort[key = lambda x: int[x.InstanceNumber]]
try:
slice_thickness = np.abs[slices[0].ImagePositionPatient[2] —
slices[1].ImagePositionPatient[2]]
except:
slice_thickness = np.abs[slices[0].SliceLocation —
slices[1].SliceLocation]
for s in slices:
s.SliceThickness = slice_thickness
return slicesdef get_pixels_hu[scans]:
image = np.stack[[s.pixel_array for s in scans]]
image = image.astype[np.int16] # Set outside-of-scan pixels to 0
# The intercept is usually -1024, so air is approximately 0
image[image == -2000] = 0
# Convert to Hounsfield units [HU]
intercept = scans[0].RescaleIntercept
slope = scans[0].RescaleSlope
if slope != 1:
image = slope * image.astype[np.float64]
image = image.astype[np.int16]
image += np.int16[intercept]
return np.array[image, dtype=np.int16]
Chạy đoạn mã sau để trích xuất các pixel DICOM cho từng vị trí lát cắt và hiển thị một lát cắt duy nhất
# set path and load files
path = ‘path/to/your/dicom/directory/’
patient_dicom = load_scan[path]
patient_pixels = get_pixels_hu[patient_dicom]#sanity check
plt.imshow[patient_pixels[326], cmap=plt.cm.bone]
Hãy sử dụng một số thao tác phân ngưỡng và hình thái học để chỉ phân chia phổi từ lồng ngực
def largest_label_volume[im, bg=-1]:
vals, counts = np.unique[im, return_counts=True] counts = counts[vals != bg]
vals = vals[vals != bg] if len[counts] > 0:
return vals[np.argmax[counts]]
else:
return Nonedef segment_lung_mask[image, fill_lung_structures=True]:
# not actually binary, but 1 and 2.
# 0 is treated as background, which we do not want
binary_image = np.array[image >= -700, dtype=np.int8]+1
labels = measure.label[binary_image]
# Pick the pixel in the very corner to determine which label is
air.
# Improvement: Pick multiple background labels from around the
patient
# More resistant to “trays” on which the patient lays cutting
the air around the person in half
background_label = labels[0,0,0]
# Fill the air around the person
binary_image[background_label == labels] = 2
# Method of filling the lung structures [that is superior to
# something like morphological closing]
if fill_lung_structures:
# For every slice we determine the largest solid structure
for i, axial_slice in enumerate[binary_image]:
axial_slice = axial_slice — 1
labeling = measure.label[axial_slice]
l_max = largest_label_volume[labeling, bg=0]
if l_max is not None: #This slice contains some lung
binary_image[i][labeling != l_max] = 1
binary_image -= 1 #Make the image actual binary
binary_image = 1-binary_image # Invert it, lungs are now 1
# Remove other air pockets inside body
labels = measure.label[binary_image, background=0]
l_max = largest_label_volume[labels, bg=0]
if l_max is not None: # There are air pockets
binary_image[labels != l_max] = 0
return binary_image
Bằng cách chạy mã bên dưới, bạn đang sử dụng các chức năng đọc lướt từ phía trên để tạo mặt nạ che phủ phổi. Chúng tôi sẽ sử dụng cả fill_lung_structures=True và fill_lung_structures=False, để cô lập phổi và các cấu trúc bên trong. Hãy chạy mã bên dưới và hiển thị một ví dụ về cách ly phổi khỏi ngực
# get masks
segmented_lungs = segment_lung_mask[patient_pixels,
fill_lung_structures=False]
segmented_lungs_fill = segment_lung_mask[patient_pixels,
fill_lung_structures=True]
internal_structures = segmented_lungs_fill - segmented_lungs# isolate lung from chest
copied_pixels = copy.deepcopy[patient_pixels]
for i, mask in enumerate[segmented_lungs_fill]:
get_high_vals = mask == 0
copied_pixels[i][get_high_vals] = 0
seg_lung_pixels = copied_pixels# sanity check
plt.imshow[seg_lung_pixels[326], cmap=plt.cm.bone]
Segmented Lung [vẫn duy trì các giá trị pixel HU ban đầu]
Nếu nó trông như thế này, thì hoàn hảo. Bạn vừa tách phổi ra khỏi phần còn lại của quá trình quét. Đừng lo lắng về mặt nạ, tôi sẽ cho bạn thấy một số hình dung sau
Phương pháp mạnh mẽ để trích xuất các tính năng bằng cách sử dụng GK Clustering
Tại sao chỉ dừng lại ở lá phổi? . Tùy thuộc vào nhiệm vụ nhất định của bạn, bạn có thể muốn phân tích các khu vực quét cụ thể hơn. Một cách khả thi để làm điều này là sử dụng GK Clustering
class GK:
def __init__[self, n_clusters=4, max_iter=100, m=2, error=1e-6]:
super[].__init__[]
self.u, self.centers, self.f = None, None, None
self.clusters_count = n_clusters
self.max_iter = max_iter
self.m = m
self.error = errordef fit[self, z]:
N = z.shape[0]
C = self.clusters_count
centers = [] u = np.random.dirichlet[np.ones[N], size=C] iteration = 0
while iteration < self.max_iter:
u2 = u.copy[] centers = self.next_centers[z, u]
f = self._covariance[z, centers, u]
dist = self._distance[z, centers, f]
u = self.next_u[dist]
iteration += 1 # Stopping rule
if norm[u - u2] < self.error:
break self.f = f
self.u = u
self.centers = centers
return centersdef gk_segment[img, clusters=5, smooth=False]: # expand dims of binary image [1 channel in z axis]
new_img = np.expand_dims[img, axis=2] # reshape
x, y, z = new_img.shape
new_img = new_img.reshape[x * y, z] # segment using GK clustering
algorithm = GK[n_clusters=clusters]
cluster_centers = algorithm.fit[new_img]
output = algorithm.predict[new_img]
segments = cluster_centers[output].astype[np.int32].reshape[x,
y] # get cluster that takes up least space [nodules / airway]
min_label = min_label_volume[segments]
segments[np.where[segments != min_label]] = 0
segments[np.where[segments == min_label]] = 1 return segments
Bây giờ, hãy chỉ phân đoạn một vài lát cắt [5 trong số đó] bằng cách sử dụng lớp GK Clustering và hiển thị một trong số chúng
# cluster slices [#324 - #328]
dist = 2
selected_slices = seg_lung_pixels[326-dist:326+[dist+1]]
gk_clustered_imgs = np.array[[gk_segment[x] for x in
selected_slices]]# display middle slice [slice #326]
plt.imshow[gk_clustered_imgs[2], cmap=plt.cm.bone]
Cấu trúc bên trong của phổi
Chỉ các cấu trúc bên trong. Hơi ngầu nhỉ?
Không tương tác
Khi trực quan hóa dữ liệu, tôi thấy rất hữu ích khi trực quan hóa từng quy trình trong tập lệnh của bạn. Điều này không chỉ giúp bạn hiểu từng bước mã của mình mà còn tạo ra một bản trình bày rất đẹp, rõ ràng
Về bản chất, đoạn mã dưới đây sẽ truyền tải một cách trực quan câu chuyện về cách bạn phân đoạn dữ liệu của mình. [1. ] ảnh gốc > [2. ] mặt nạ nhị phân bao phủ phổi > [3. ] làm nổi bật các cấu trúc bên trong phổi bằng cách sử dụng một trong các mặt nạ > [4. ] trích xuất các cấu trúc bên trong bằng phân cụm GK
# pick random slice
slice_id = 326plt.figure[1]
plt.title['Original Dicom']
plt.imshow[patient_pixels[slice_id], cmap=plt.cm.bone]plt.figure[2]
plt.title['Lung Mask']
plt.imshow[segmented_lungs_fill[slice_id], cmap=plt.cm.bone]plt.figure[3]
plt.title['Parenchyma Mask on Segmented Lung']
plt.imshow[seg_lung_pixels[slice_id], cmap=plt.cm.bone]
plt.imshow[internal_structures[slice_id], cmap='jet', alpha=0.7]plt.figure[4]
plt.title['Internal Structures Using GK Clustering']
plt.imshow[gk_segment_lung[2], cmap=plt.cm.bone]
tương tác
Thực sự chỉ có một số cách để hiển thị nhiều hình ảnh trên máy tính xách tay jupyter - vẽ từng hình ảnh theo cách thủ công hoặc tạo một biểu đồ có nhiều cột và hàng để hiển thị trong một hình. Thật không may, điều này không hữu ích cho những người muốn quét qua từng hình ảnh và hiểu rõ hơn về dữ liệu. Với mã bên dưới, bạn tạo một thanh trượt tương tác cho phép bạn cuộn qua các hình ảnh
________số 8Bước 5. Kỹ thuật hình ảnh 3D
Trong hình ảnh y tế, người ta thường sử dụng thuật toán đồ họa máy tính được gọi là Marching Cubes để trích xuất một bề mặt từ dữ liệu ba chiều. Trong trường hợp của chúng tôi, chúng tôi sẽ sử dụng biện pháp chức năng tích hợp của skigage. hành quân_cubes_lewiner
không tương tác
Cho dù bạn đang áp dụng mạng thần kinh tích chập 2D hay 3D cho tập dữ liệu của mình, thì vẫn có những lợi ích khi xem dữ liệu DICOM của bạn ở chế độ 3D. Nói một cách triết học, nếu bạn muốn thuật toán của mình hiểu rõ về dữ liệu để thực hiện theo cách bạn muốn — thì bạn cũng nên hiểu rõ
def plot_3d[image]:
# Position the scan upright,
# so the head of the patient would be at the top facing the
# camera
p = image.transpose[2,1,0]
verts, faces, _, _ = measure.marching_cubes_lewiner[p] fig = plt.figure[figsize=[10, 10]]
ax = fig.add_subplot[111, projection='3d'] # Fancy indexing: `verts[faces]` to generate a collection of
# triangles
mesh = Poly3DCollection[verts[faces], alpha=0.70]
face_color = [0.45, 0.45, 0.75]
mesh.set_facecolor[face_color]
ax.add_collection3d[mesh] ax.set_xlim[0, p.shape[0]]
ax.set_ylim[0, p.shape[1]]
ax.set_zlim[0, p.shape[2]] plt.show[]# run visualization
plot_3d[internal_structures_mask]
Khá tuyệt phải không?
tương tác
Vì việc tạo biểu đồ 3D tương tác có thể mất một chút thời gian [nếu bạn cung cấp nhiều dữ liệu], tốt nhất nên sử dụng điều này khi bạn đã tách riêng một vùng cụ thể trong quá trình quét của mình. Ví dụ: có thể bạn đã tìm thấy một tính năng bất thường trong hình ảnh của mình và bạn muốn tương tác với nó ở chế độ 3D
Mã bên dưới hiển thị một ví dụ về cắt xén một vùng cụ thể, thu được các lát lân cận [để chúng tôi có thể thu được nhiều lần quét 2D để tạo thành một khối 3D] và vẽ đồ thị tương tác trong 3D
def load_scan[path]:0
slices = [pydicom.dcmread[path + ‘/’ + s] for s in
os.listdir[path]]
slices = [s for s in slices if ‘SliceLocation’ in s]
slices.sort[key = lambda x: int[x.InstanceNumber]]
try:
slice_thickness = np.abs[slices[0].ImagePositionPatient[2] —
slices[1].ImagePositionPatient[2]]
except:
slice_thickness = np.abs[slices[0].SliceLocation —
slices[1].SliceLocation]
for s in slices:
s.SliceThickness = slice_thickness
return slicesdef get_pixels_hu[scans]:
image = np.stack[[s.pixel_array for s in scans]]
image = image.astype[np.int16] # Set outside-of-scan pixels to 0
# The intercept is usually -1024, so air is approximately 0
image[image == -2000] = 0
# Convert to Hounsfield units [HU]
intercept = scans[0].RescaleIntercept
slope = scans[0].RescaleSlope
if slope != 1:
image = slope * image.astype[np.float64]
image = image.astype[np.int16]
image += np.int16[intercept]
return np.array[image, dtype=np.int16]
Bây giờ, hãy sử dụng các chức năng của trình trợ giúp, cùng với gk_clustered_imgs để cắt bỏ một số vùng và hiển thị nó trông như thế nào trong 3D tương tác [bên dưới chỉ là một hình ảnh, nhưng trong sổ ghi chép của bạn, bạn có thể nhấp và di chuyển nó xung quanh, lol]
def load_scan[path]:1
slices = [pydicom.dcmread[path + ‘/’ + s] for s in
os.listdir[path]]
slices = [s for s in slices if ‘SliceLocation’ in s]
slices.sort[key = lambda x: int[x.InstanceNumber]]
try:
slice_thickness = np.abs[slices[0].ImagePositionPatient[2] —
slices[1].ImagePositionPatient[2]]
except:
slice_thickness = np.abs[slices[0].SliceLocation —
slices[1].SliceLocation]
for s in slices:
s.SliceThickness = slice_thickness
return slicesdef get_pixels_hu[scans]:
image = np.stack[[s.pixel_array for s in scans]]
image = image.astype[np.int16] # Set outside-of-scan pixels to 0
# The intercept is usually -1024, so air is approximately 0
image[image == -2000] = 0
# Convert to Hounsfield units [HU]
intercept = scans[0].RescaleIntercept
slope = scans[0].RescaleSlope
if slope != 1:
image = slope * image.astype[np.float64]
image = image.astype[np.int16]
image += np.int16[intercept]
return np.array[image, dtype=np.int16]
Sự kết luận
chết tiệt. Bạn vừa lập trình mã không-không-không-ghê-tồi nhất liên quan đến các phương pháp phân cụm thống kê, thuật toán đồ họa máy tính và xử lý hình ảnh. Chúng tôi thậm chí không cần sử dụng những từ buzz lớn để làm cho nó nghe hay hơn;]