簡介 3D 電腦斷層影像判讀新冠肺炎、心電圖異常訊號偵測

相信不管 Run 到內、外、婦、兒、急、或其他各個科別，心電圖都是非常常見的檢查；而電腦斷層更是幫助臨床醫師掌握病人外觀看不見的體腔內所發生疾病問題的利器。這次，就讓我們來瞧瞧 CNN 如何運用在立體的胸部電腦斷層與一維的心電圖上做分類吧！

🫁 3D Image：胸部電腦斷層

之前在做 Deep Learning 單張 2D 影像辨識時，我們把圖片輸入 Model 前，會先掌握它的寬度 (Width)、高度 (Height)、深度 (Depth) 或通道數 (Channel)，而 Channel 會取決於資料的形式，像黑白圖片是 1、彩色圖片是 3、MRI 有 T1 / T1 weighted / T2 weighted 3 個參數也有 3 個通道，也就是每個 Pixel 我們用幾個數值來代表它，(H, W, C) 會讓我們定義的 Kernel 在圖片上做掃瞄卷積得到 Feature Maps，Shape 也同樣有 (H, W, C)，這裡 Output 的 C 則取決於該卷積層用了幾個 Filters。

那麼，3D 的立體影像，譬如 3D CT、MRI 的影像辨識要怎麼做呢？其實3D 影像和 2D 只相差了一個維度：拍攝的張數 (Slice)，Image Volume 可用 (Slice, H, W, C) 來表示；影片也可以想像成是連續照片的拼接，由許多 2D 圖像組合成 (Time, H, W, C)，因此實作上也會把影片當成 3D 影像來做。

實作上，如果各位朋友對於先前所分享的 2D Convolution 有一定的認識，套用到 3D 影像就變得相對簡單許多。3D 對於 Input Image 也有一定的規範，只是比 2D 多了一個 Slice 的方向，3D 的 Kernel 在移動時會有張數、寬度、高度 3 個方向，產生的 Feature Maps 也是 (Slice, H, W, C) 的 Shape。卷積層參數的設定上，與 2D 類似，都要設定 Filter 個數、Kernel 大小、Activation Function，Output 的 Channel 取決於 Filter 的數量。

無論是分類、切割、物件偵測的模型，2D 都能轉為 3D 的形式，只是多了 Slice 的維度，Layers 則是把 2D 改成 3D，例如 Conv3D、Conv3DTranspose, GlobalAveragePooling3D、MaxPooling3D、UpSampling3D… 把原始 2D 的模型改為 3D 基本上就把影像辨識網路建立完成，除了 Shape 多了一個方向要拼接與疊加、運算資源的要求比較高之外，其他的差距不會太大。

我們用立體的 Chest CT 胸部電腦斷層影像分類來舉例，一般在醫院 PACs 系統會把立體影像的不同張數用 2D 圖片顯示。

使用的套件和 2D 不會差異太多，但會多引用一個 scipy 的套件，來針對 2D 以上的立體體積做操作 (像是 Resize、Reshape，針對目標大小來內插設定每個方向要放大縮小幾倍)，並用 nibabel 來讀取 Nifti 格式的影像。

import os
import zipfile
import numpy as np
import tensorflow as tf
from tqdm.auto import tqdm
import matplotlib.pyplot as plt

from tensorflow import keras
from tensorflow.keras import layers

from scipy import ndimage

大家知道 CT Scan 原始的 Voxel Intensity 是用 Hounsfield units (HU) 來表示，在我們的資料集 HU 從 -1024 到 2000，超過 400 大概就是骨頭，所以我們的 Threshold 先設在 -1000 到 400 之間，並會 Normalize HU 的值到 0–1 之間的數值。此外，資料前處理的部分，我們也會旋轉 90 度、Resize 寬度、高度、深度成 128x128x64，再把資料及分成 Training 和 Validation 的子集。

def read_nifti_file(filepath):
    # Read file and get raw data
    return nib.load(filepath).get_fdata()


def normalize(volume):
    # Normalize the volume -1000 ~ 400 -> 0~1
    min = -1000
    max = 400
    # 2D: pixel -> 3D: voxel
    volume[volume < min] = min
    volume[volume > max] = max
    volume = (volume - min) / (max - min)
    volume = volume.astype("float32")
    return volume


def resize_volume(img):
    """Resize across z-axis"""
    # Set the desired depth
    desired_depth = 64
    desired_width = 128
    desired_height = 128
    # Get current depth
    current_depth = img.shape[-1]
    current_width = img.shape[0]
    current_height = img.shape[1]
    # Compute depth factor
    depth = current_depth / desired_depth
    width = current_width / desired_width
    height = current_height / desired_height
    depth_factor = 1 / depth
    width_factor = 1 / width
    height_factor = 1 / height
    # Rotate
    # img = ndimage.rotate(img, 90, reshape=False)
    # Resize across z-axis
    img = ndimage.zoom(img, (width_factor, height_factor, depth_factor), order=1)
    return img


def process_scan(path):
    """Read and resize volume"""
    # Read scan
    volume = read_nifti_file(path)
    # Normalize
    volume = normalize(volume)
    # Resize width, height and depth
    volume = resize_volume(volume)
    # add channel dim
    volume = np.expand_dims(volume, axis=-1) # (H, W, Slice, 1)
    return volume

接著，會透過在訓練時隨機旋轉圖片做 Data Augmentation，3D 圖像比較少套件在實作，每次針對立體影像中不同 Slice 做完 Augmentation 再拼接起來；因為 Input 的體積多了 Slice 的維度，Batch Size 通常不會開到太大（Colab 的話頂多 3–4），GPU 運氣不好抽到比較小的記憶體就可能會爆掉，無法處理太龐大的模型。DataGenerator 來讀取資料，一次讀一個 Batch，也可以先視覺化分好的資料是不是符合我們預期的樣子、答案對得上 Input 的圖片。

# Define data loaders.

class DataGenerator(tf.keras.utils.Sequence):
    def __init__(self, x, y, batch_size, shuffle=True, aug=False):
        self.x = x
        self.y = y
        self.batch_size = batch_size
        self.shuffle = shuffle
        self.aug = aug
        self.indexes = np.arange(len(x))

    def __len__(self):
        return int(np.ceil(len(self.x)/self.batch_size))

    def __getitem__(self, index):
        idxs = self.indexes[index*self.batch_size: (index+1)*self.batch_size]
        x = self.x[idxs]
        y = self.y[idxs]
        if self.aug:
            x = self.augmentation(x)
        return x, y
    
    def on_epoch_end(self):
        if self.shuffle: np.random.shuffle(self.indexes)
    
    def augmentation(self, volumes):
        for i, volume in enumerate(volumes):
            volumes[i] = rotate(volume)
        return volumes
        
batch_size=4
train_gen = DataGenerator(x_train, y_train, batch_size=batch_size, shuffle=True, aug=True)
val_gen = DataGenerator(x_val, y_val, batch_size=batch_size, shuffle=False, aug=False)

資料都處理好之後，我們就要建模囉！這裡引用的是偵測 Tuberculosis 的 3D CNNs Resnet 模型架構，基本上由 Conv3D、BatchNormalization、ReLU Function 構成，Filter 大小是 3x3x3，有加上 Padding，會再接上 MaxPooling3D，最後 GlobalAveragePooling3D 並通過 Sigmoid Function 來做二元分類。

def conv_bn(x, filters, kernel_size=3, strides=1, activation=True):
    x = layers.Conv3D(filters, kernel_size, padding='same', strides=strides)(x)
    x = layers.BatchNormalization()(x)
    if activation:
        x = layers.ReLU()(x)
    return x

def res_block(input_tensor, out_dims, strides=1):
    x = conv_bn(input_tensor, out_dims, strides=strides)
    x = conv_bn(x, out_dims, activation=False)

    if strides == 1:
        short_cut = input_tensor
    else:
        short_cut = conv_bn(input_tensor, out_dims, kernel_size=1, strides=strides, activation=False)

    x = layers.add([x, short_cut])
    x = layers.ReLU()(x)
    return x

def get_model(width=128, height=128, depth=64):
    """Build a 3D convolutional neural network model."""

    inputs = keras.Input((width, height, depth, 1))
    x = conv_bn(inputs, 64)
    x = layers.MaxPooling3D()(x)

    for i in range(2): x = res_block(x, 64)
    for i in range(2): x = res_block(x, 128, strides=2)
    for i in range(2): x = res_block(x, 256, strides=2)

    x = layers.GlobalAveragePooling3D()(x)
    outputs = layers.Dense(units=1, activation="sigmoid")(x)

    # Define the model.
    model = keras.Model(inputs, outputs, name="Sigrid3Dcnn")
    return model


# Build model.
model = get_model(width=128, height=128, depth=64)
model.summary()

Metric 先用 Accuracy（當然自己寫 Dice Coefficient 也是可以，放在 Compile 裡就會自動運算），Loss Function 用二元分類的 binary_crossentropy，開始來做訓練；和 2D 相比，不論是 Resize 還是讀取資料量，花費的時間都會多上許多。

# Compile model.
model.compile(
    loss="binary_crossentropy",
    optimizer=keras.optimizers.Adam(),
    metrics=["accuracy"],
)

# Define callbacks.
checkpoint_cb = keras.callbacks.ModelCheckpoint("3d_classification.h5", save_best_only=True)
early_stopping_cb = keras.callbacks.EarlyStopping(monitor="val_loss", patience=15)

# Train the model, doing validation at the end of each epoch
epochs = 100
model.fit(
    train_gen,
    validation_data=val_gen,
    epochs=epochs,
    callbacks=[checkpoint_cb, early_stopping_cb],
)

訓練好之後，我們來瞧瞧 Training 和 Validation 資料集的 Accuracy 和 Loss 👇

🫀 1D Image：心電圖

看完 3D Image 的影像辨識後，1D Data 也是依樣畫葫蘆：3D 比 2D 多一個 Slice 的維度，而 1D 則是比 2D 少一個方向。1D Data 通常和時間序列有關，腦波、血壓、呼吸、心電圖等訊號，資料格式通常是 (Time, Channels)，根據不同時間點及通道數量決定。心電圖最簡單的只有 1 個 Channel，而在醫院我們當 Intern / PGY 時推推車去幫病人做的 12-Lead EKG，則會根據 Leads 之間的電壓差產出 12 個 Channels 的訊號。

1D Time-Series Data 用 CNN 處理時，會有時間及通道兩個方向，Kernel 沿著某個連續時間維度的資訊做 Convolution，不像 2D 圖片是沿著高度、寬度，比 2D 還要單純些。卷積層參數的設定上，與 2D 和 3D 很像，也都要設定 Filter 個數、Kernel 大小、Activation Function，Output 的 Channel 取決於 Filter 的數量。從這些例子我們可以知道，「X」D 指的是 Kernel 移動的方向有幾個，2D 沿著寬、高做卷積，3D 還多了 Slice，1D 則只有時間的方向。

這個例子中，我們要把心電圖 1D Signal 分成 5 個類別，把資料轉成 numpy array 就可以來做訓練，這個資料集共有 1459 筆資料、每筆資料有 188 個欄位，最後一個 Column 代表資料的類別，其他 Column 裡的數值就是每個時間點的訊號。Training 和 Validation Datasets 用 8:2 的比例來分，前處理時把 Label 和訊號區分出來，並把 Label 做 One-Hot Encoding 的處理；讓我們來畫一張 EKG 看看長什麼樣子，像底下這張的 Label 是第 0 類，也就是正常的。

和 2D 轉 3D 模型一樣，我們建模時把 2D 轉成 1D 即可，像 Conv1D、MaxPooling1D、GlobalAveragePooling1D… 抽取特徵時我們用 1D 的卷積，再接全連接層做分類，因為是多元分類所以最後的 Activation Function 用 Softmax。Compile 時用的 Optimizer 是 Adam、Loss Function 用多元分類的 categorical_crossentropy，Metric 以 Accuracy 來評估。因為參數用得比 2D、3D 還少，相較之下是比較輕盈的模型，像這裡我只花 58 秒就 Train 好了，但 3D 的要花將近 50 分鐘。

def build_model():
    inputs = layers.Input(shape=x_train.shape[1:])
    conv1_1 = layers.Conv1D(64, (6), activation='relu')(inputs)
    conv1_1 = layers.BatchNormalization()(conv1_1)
    pool1 = layers.Conv1D(64, 3, strides=(2), padding="same")(conv1_1)
    conv2_1 = layers.Conv1D(128, (3), activation='relu')(pool1)
    conv2_1 = layers.BatchNormalization()(conv2_1)
    pool2 = layers.Conv1D(128, 3, strides=(2), padding="same")(conv2_1)
    conv3_1 = layers.Conv1D(256, (3), activation='relu')(pool2)
    conv3_1 = layers.BatchNormalization()(conv3_1)
    pool3 = layers.Conv1D(256, 3, strides=(2), padding="same")(conv3_1)
    pool3 = layers.MaxPooling1D()(pool3)
    flatten = layers.GlobalAveragePooling1D()(pool3)
    dense_end1 = layers.Dense(64, activation='relu')(flatten)
    dense_end2 = layers.Dense(32, activation='relu')(dense_end1)
    main_output = layers.Dense(num_classes, activation='softmax', name='main_output')(dense_end2)
    
    model = models.Model(inputs, main_output)
    model.compile(optimizer='adam', 
                  loss='categorical_crossentropy',
                  metrics = ['accuracy'])
    return model

分類報告說整體 Accuracy 是 94%，但要注意是不是每個類別都做得很好？細看之下，第 0 個類別 f1-score 是0.97，但第 1、3 個類別則分別只有 0.63 和 0.67，這是因為資料量極度不平均，我們可以看到 Training 資料裡，第 0 類有 239 筆，第 1、3 個類別則分別只有 9 和 4 筆，醫療資料蠻常碰到這類 Imbalance 的問題，將來可能可以針對資料量比較少的類別來做修正。

實務上，每間醫院機器的訊號輸出格式都不太一樣，需要視個別資料的狀況來做轉換；而訊號量是否足夠、取樣頻率要設定多少、分析時要看連續的時間序列還是分段做分析，這些都是根據我們想問的臨床問題來定義的。但 1D 和 3D CNN 的小麻煩點，就在於可以直接拿來做 Transfer Learning 的 Model 很少，常常要自己經過挑選和修改來做成適合手邊資料的模型，各位朋友可以自己試試看把知名 2D Model 拿來做調整，再來分析看看結果的差異是多少。

基因序列在某種邏輯程度上也是 1D Data，Input 是基因的長度，每個位置都有 A、T、C、G 四種可能 (Channel = 4)，就可以用 1D CNN 針對 Gene Sequence 來做分類。

總算把 1D、2D、3D 醫學影像辨識的 Computer Vision 簡單地梳理完畢，相信以後還是會持續有新的 Model、影像擷取及處理的技術出現，讓我們一起期待未來機器學習與醫學碰撞出的火花吧 🤩

#CC