สมมติว่าอยู่ดี ๆ เราก็หายใจลำบาก หอบตัวโยน โดยไม่มีสาเหตุ ไม่มีอาการล่วงหน้าใด ๆ หรือว่าเราจะเป็น ภาวะปอดรั่ว ใน ep นี้ เราจะมาใช้ Machine Learning และ Deep Neural Network พัฒนาโปรแกรม AI การแพทย์ ช่วยวินิจฉัยภาวะปอดรั่ว หรือ Pneumothorax นี้กัน

Pneumothorax หรือ ภาวะปอดรั่ว คืออะไร

Illustration depicting a collapsed lung or pneumothorax. Credit https://en.wikipedia.org/wiki/File:Blausen_0742_Pneumothorax.png
Illustration depicting a collapsed lung or pneumothorax. Credit https://en.wikipedia.org/wiki/File:Blausen_0742_Pneumothorax.png

Pneumothorax หรือ ภาวะปอดรั่ว คือการที่มีอากาศรั่วไหลเข้าไปอยู่ในช่องระหว่างปอด และ ผนังอกด้านใน ทำให้เบียดเนื้อปอดให้ขยายตัวได้ไม่เต็มที่ ปอดทำงานได้ไม่ดี ส่งผลต่อการหายใจ ภาวะปอดรั่วนี้ต้องได้รับการดูแลจากแพทย์โดยด่วน หากปล่อยไว้จะเป็นอันตรายถึงแก่ชีวิตได้

Pneumothorax หรือ ภาวะปอดรั่ว เกิดจากอะไร

Right-sided pneumothorax (right side of image) on CT scan of the chest with chest tube in place. Credit https://en.wikipedia.org/wiki/File:Pneumothorax_CT.jpg
Right-sided pneumothorax (right side of image) on CT scan of the chest with chest tube in place. Credit https://en.wikipedia.org/wiki/File:Pneumothorax_CT.jpg

Pneumothorax หรือ ภาวะปอดรั่ว เป็นภาวะที่พบได้ในอุบัติเหตุที่มีการกระแทกบริเวณหน้าอก มีลมซึมเข้ามาจากภายนอกทรวงอก หรือเกิดขึ้นเองจากการติดเชื้อภายใน โรคเกี่ยวกับปอด หอบหืด มะเร็งปอด etc.

การวินิจฉัย Pneumothorax หรือ ภาวะปอดรั่ว

Pneumothorax Overlay. Credit https://www.kaggle.com/c/siim-acr-pneumothorax-segmentation/
Pneumothorax Overlay. Credit https://www.kaggle.com/c/siim-acr-pneumothorax-segmentation/

Pneumothorax หรือ ภาวะปอดรั่ว มักจะถูกวินิจฉัยโดย Radiologist (รังสีแพทย์ / แพทย์รังสีวิทยา) เป็นผู้อ่านฟิล์ม X-Ray เป็นหลัก แต่บางครั้งก็เป็นการยากที่จะยืนยันว่าผู้ป่วยเป็นภาวะปอดรั่วจริง ๆ จากการเอกซเรย์ทรวงอกเพียงอย่างเดียว บางทีอาจต้องใช้วิธีการตรวจอื่น ๆ ร่วมด้วย เช่น CT Scan, ฟังเสียงหายใจ และ ตรวจเลือดเพื่อวัดระดับออกซิเจนในเลือด

AI การแพทย์ ช่วยวิจิฉัยภาวะปอดรั่ว (Pneumothorax)

AI Diagnose Pneumothorax. Credit https://www.kaggle.com/c/siim-acr-pneumothorax-segmentation/
AI Diagnose Pneumothorax. Credit https://www.kaggle.com/c/siim-acr-pneumothorax-segmentation/

AI การแพทย์เฉพาะทาง ช่วยวิจิฉัยภาวะปอดรั่ว (Pneumothorax) อย่างแม่นยำ จะเป็นประโยชน์อย่างมาก ในการช่วยคัดแยกผู้ป่วย ในสถานพยาบาล ช่วย Radiologist (รังสีแพทย์ / แพทย์รังสีวิทยา) จัดลำดับความสำคัญตามอาการมากน้อย

หรือเป็นเครื่องมือช่วยในการวินิจฉัยเบื้องต้น สำหรับบุคลากรทางการแพทย์ที่ไม่ใช่ Radiologist (รังสีแพทย์ / แพทย์รังสีวิทยา) กรณีบุคลากรไม่เพียงพอ เพื่อเพิ่มอัตราการรอดชีวิตของผู้ป่วย

Model Architecture

U-Net architecture. Credit https://arxiv.org/abs/1505.04597
U-Net architecture. Credit https://arxiv.org/abs/1505.04597

Model Architecture เหมือน Image Segmentation ep.1 คือ U-Net โดยใช้ Pretrained ResNet34 เป็น Encoder และเทรนด้วย Progressive Resizing

Dice Coefficient

Metric ที่เราจะใช้ใน ep นี้ คือ Dice Score หรือ Sørensen–Dice coefficient เป็นการคำนวนสถิติที่ใช้วัดค่าความเหมือนกัน ของ 2 ข้อมูลตัวอย่าง ดังสมการด้านล่าง Dice Coefficient มีอีกชื่อหนึ่ง คือ F1 Score

รายละเอียดเพิ่มเติมใน ep Metrics / Confusion Matrix

\( DSC = \frac{2 |X \cap Y|}{|X|+ |Y|} = \frac{2 TP}{2 TP + FP + FN} \)

Dice Metric ค่อนข้าง ใกล้เคียงกับ Intersection-Over-Union (IoU, Jaccard Index) จะอธิบายต่อไป

เรามาเริ่มกันเลยดีกว่า

Open In Colab

สมมติว่าอยู่ดี ๆ เราก็หายใจลำบาก หอบตัวโยน โดยไม่มีสาเหตุ ไม่มีอาการล่วงหน้าใด ๆ หรือว่าเราจะเป็น ภาวะปอดรั่ว เราจะมาใช้ AI ช่วยวินิจฉัยเรื่องนี้กัน

Pneumothorax หรือ ภาวะปอดรั่ว คือการที่มีอากาศรั่วไหลเข้าไปอยู่ในช่องระหว่างปอด และ ผนังอกด้านใน ทำให้เบียดเนื้อปอดให้ขยายตัวได้ไม่เต็มที่ ปอดทำงานได้ไม่ดี ส่งผลต่อการหายใจ ภาวะปอดรั่วนี้ต้องได้รับการดูแลจากแพทย์โดยด่วน หากปล่อยไว้จะเป็นอันตรายถึงแก่ชีวิตได้

Pneumothorax หรือ ภาวะปอดรั่ว เป็นภาวะที่พบได้ในอุบัติเหตุที่มีการกระแทกบริเวณหน้าอก มีลมซึมเข้ามาจากภายนอกทรวงอก หรือเกิดขึ้นเองจากการติดเชื้อภายใน โรคเกี่ยวกับปอด หอบหืด มะเร็งปอด etc.

Pneumothorax หรือ ภาวะปอดรั่ว มักจะถูกวินิจฉัยโดย Radiologist (แพทย์รังสีวิทยา) อ่านฟิล์ม X-Ray แต่บางครั้งก็เป็นการยากที่จะยืนยันภาวะปอดรั่ว

AI วิจิฉัยภาวะปอดรั่ว (Pneumothorax) อย่างแม่นยำ จะเป็นประโยชน์อย่างมาก ในการช่วยคัดแยกผู้ป่วย ในสถานพยาบาล จัดลำดับความสำคัญตามอาการมากน้อย หรือช่วยในการวินิจฉัยเบื้องต้น สำหรับบุคลากรทางการแพทย์ที่ไม่ใช่ Radiologist (แพทย์รังสีวิทยา)

In [0]:
# Github URL
# https://github.com/gnoparus/bualabs/blob/master/nbs/03b-image-segmentation-siim-acr-pneumothorax-segmentation.ipynb

0. Install

เราจะต้อง Install kaggle เพื่อ Download Dataset

In [0]:
## Colab
# ! curl -s https://course.fast.ai/setup/colab | bash

# ! pip install kaggle --upgrade
# ! pip install pydicom

เช็ค GPU Memory ใน Colab

In [0]:
# # memory footprint support libraries/code
# !ln -sf /opt/bin/nvidia-smi /usr/bin/nvidia-smi
# !pip install gputil
# !pip install psutil
# !pip install humanize
# import psutil
# import humanize
# import os
# import GPUtil as GPU

# GPUs = GPU.getGPUs()
# # XXX: only one GPU on Colab and isn’t guaranteed
# gpu = GPUs[0]
# def printm():
#  process = psutil.Process(os.getpid())
#  print("Gen RAM Free: " + humanize.naturalsize( psutil.virtual_memory().available ), " | Proc size: " + humanize.naturalsize( process.memory_info().rss))
#  print("GPU RAM Free: {0:.0f}MB | Used: {1:.0f}MB | Util {2:3.0f}% | Total {3:.0f}MB".format(gpu.memoryFree, gpu.memoryUsed, gpu.memoryUtil*100, gpu.memoryTotal))
# printm() 

สำหรับ Clear Memory

In [0]:
# !kill -9 -1

1. Import Library

Import pydicom สำหรับเปิดไฟล์ภาพเอ็กซ์เรย์ DICOM นามสกุล dcm และ Libray อื่น ๆ ที่เราต้องการใช้

In [0]:
import math
import numpy as np
import pandas as pd

import os
import glob

import pydicom

from pathlib import Path
from matplotlib import cm
from matplotlib import pyplot as plt

import fastai
from fastai.vision import *
In [0]:
fastai.__version__
Out[0]:
'1.0.59'

2. เตรียม Path สำหรับดาวน์โหลดข้อมูล

กำหนด path ของ Config File และ Dataset ว่าจะอยู่ใน Google Drive ถ้าเราใช้ Google Colab หรือ อยู่ใน HOME ถ้าเราใช้ VM ธรรมดา และกำหนด Environment Variable ไปยังโฟลเดอร์ที่เก็บ kaggle.json

ในกรณีใช้ Colab ให้ Mount Google Drive เพื่อดึง Config File มาจาก Google Drive ส่วนตัวของเรา เมื่อเรารัน Cell ด้านล่างจะมีลิงค์ปรากฎขึ้นมาให้เรา Login กด Approve แล้ว Copy Authorization Code มาใส่ในช่องด้านล่าง แล้วกด Enter

ในเคสนี้เราจะดึงข้อมูลจากหลาย Dataset เนื่องจากเราจะใช้การเทรนแบบ Progressive Resizing ใช้รูปหลายขนาดในการเทรน

In [0]:
dataset = 'siim-acr-pneumothorax-segmentation'
dataset2 = 'jesperdramsch/siim-acr-pneumothorax-segmentation-data'
dataset3 = 'iafoss/siimacr-pneumothorax-segmentation-data-128'
dataset4 = 'iafoss/siimacr-pneumothorax-segmentation-data-256'
dataset5 = 'iafoss/siimacr-pneumothorax-segmentation-data-512'
dataset6 = 'iafoss/siimacr-pneumothorax-segmentation-data-1024'
dataset_test = 'ivanzhovannik/siim_stage2_png'

# # Google Colab
# config_path = Path('/content/drive')
# data_path_base = Path('/content/datasets/')

# data_path = data_path_base/dataset
# data_path2 = data_path_base/dataset2
# data_path3 = data_path_base/dataset3
# data_path4 = data_path_base/dataset4
# data_path5 = data_path_base/dataset5
# data_path6 = data_path_base/dataset6
# data_path_test = data_path_base/dataset_test

# from google.colab import drive
# drive.mount(str(config_path))
# os.environ['KAGGLE_CONFIG_DIR'] = f"{config_path}/My Drive/.kaggle"

# VM
config_path = Path(os.getenv("HOME"))
data_path = config_path/"datasets"/dataset
data_path2 = config_path/"datasets"/dataset2
data_path3 = config_path/"datasets"/dataset3
data_path4 = config_path/"datasets"/dataset4
data_path5 = config_path/"datasets"/dataset5
data_path6 = config_path/"datasets"/dataset6
data_path_test = config_path/"datasets"/dataset_test

data_path.mkdir(parents=True, exist_ok=True)
os.environ['KAGGLE_CONFIG_DIR'] = f"{config_path}/.kaggle"

3. Dataset

ในเคสนี้ เราจะ Download ข้อมูล Dataset ที่เกี่ยวข้องทั้งหมดมาเก็บไว้ แบ่งเป็นรูปฟิล์ม X-Ray ขนาดต่าง ๆ ข้อมูลต้นฉบับ ข้อมูล Test Set, etc.

Dataset เราจะดึงจาก Kaggle วิธี Download kaggle.json ให้ดูจาก ep ที่แล้ว

เมื่อได้ kaggle.json มาแล้ว ในกรณีใช้ Google Colab ให้นำมาใส่ไว้ในโฟลเดอร์ My Drive/.kaggle ใน Google Drive ของเรา เป็น My Drive/.kaggle/kaggle.json ถ้าใช้ VM ให้ใส่ใน HOME/.kaggle/

สั่งดาวน์โหลด Dataset จาก Kaggle พร้อมทั้ง unzip ไว้ใน data_path

In [0]:
# !kaggle competitions download -c {dataset} -p "{data_path}"
# !kaggle datasets download {dataset2} -p "{data_path2}" --unzip
# !kaggle datasets download {dataset3} -p "{data_path3}" --unzip
# !kaggle datasets download {dataset4} -p "{data_path4}" --unzip
# !kaggle datasets download {dataset5} -p "{data_path5}" --unzip
# !kaggle datasets download {dataset6} -p "{data_path6}" --unzip
# !kaggle datasets download {dataset_test} -p "{data_path_test}" --unzip

Unzip ไฟล์ที่ดาวน์โหลดจาก Kaggle Competition

In [0]:
# ! unzip -q {data_path}/siim-acr-pneumothorax-segmentation.zip -d {data_path}
# ! unzip -q {data_path}/stage_2_images.zip -d {data_path}/stage_2_images
# ! unzip -q {data_path}/stage_2_train.csv.zip -d {data_path}

Import ฟังก์ชันสำหรับแปลงข้อมูล Label ที่อยู่ในรูปแบบ Run-length encoding (RLE) to Mask ที่ให้มากับ Dataset

In [0]:
import sys
sys.path.insert(0, str(data_path))

from mask_functions import *

4. Data

4.1 DICOM Image

เปิดดูข้อมูลฟิล์มเอ็กซ์เรย์

In [0]:
def show_dcm_info(dataset):
    print("Filename.........:", file_path)
    print("Storage type.....:", dataset.SOPClassUID)
    print()

    pat_name = dataset.PatientName
    display_name = pat_name.family_name + ", " + pat_name.given_name
    print("Patient's name......:", display_name)
    print("Patient id..........:", dataset.PatientID)
    print("Patient's Age.......:", dataset.PatientAge)
    print("Patient's Sex.......:", dataset.PatientSex)
    print("Modality............:", dataset.Modality)
    print("Body Part Examined..:", dataset.BodyPartExamined)
    print("View Position.......:", dataset.ViewPosition)
    
    if 'PixelData' in dataset:
        rows = int(dataset.Rows)
        cols = int(dataset.Columns)
        print("Image size.......: {rows:d} x {cols:d}, {size:d} bytes".format(
            rows=rows, cols=cols, size=len(dataset.PixelData)))
        if 'PixelSpacing' in dataset:
            print("Pixel spacing....:", dataset.PixelSpacing)

def plot_pixel_array(dataset, figsize=(10,10)):
    plt.figure(figsize=figsize)
    plt.imshow(dataset.pixel_array, cmap=plt.cm.bone)
    plt.show()

วนลูปเปิดดูไฟล์ DICOM ที่จะมีทั้งรูป และ Metadata ข้อมูลประกอบ อยู่ในไฟล์เดียวกัน

In [0]:
for file_path in glob.glob(str(data_path2/'pneumothorax/dicom-images-train/*/*/*.dcm')):
    dataset = pydicom.dcmread(file_path)
    show_dcm_info(dataset)
    plot_pixel_array(dataset)
    break # Comment this out to see all
Filename.........: /home/jupyter/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/pneumothorax/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.4344.1517875182.377197/1.2.276.0.7230010.3.1.3.8323329.4344.1517875182.377196/1.2.276.0.7230010.3.1.4.8323329.4344.1517875182.377198.dcm
Storage type.....: 1.2.840.10008.5.1.4.1.1.7

Patient's name......: 30bf74c6-0ba7-4d68-bfc5-892841547503, 
Patient id..........: 30bf74c6-0ba7-4d68-bfc5-892841547503
Patient's Age.......: 35
Patient's Sex.......: F
Modality............: CR
Body Part Examined..: CHEST
View Position.......: PA
Image size.......: 1024 x 1024, 132184 bytes
Pixel spacing....: ['0.171', '0.171']

ลองดูทีละหลาย ๆ ไฟล์

In [0]:
start = 22   # Starting index of images
num_img = 3 # Total number of images to show

fig, ax = plt.subplots(nrows=1, ncols=num_img, sharey=True, figsize=(num_img*5,5))
for q, file_path in enumerate(glob.glob(str(data_path2/'pneumothorax/dicom-images-train/*/*/*.dcm'))[start:start+num_img]):
    dataset = pydicom.dcmread(file_path)
    #show_dcm_info(dataset)
    
    ax[q].imshow(dataset.pixel_array, cmap=plt.cm.bone)

ดู Label ที่มาในรูปแบบ Run-length encoding (RLE) ดังที่เห็นใน Column EncodedPixels

-1 หมายถึง ปอดปกติ

In [0]:
df = pd.read_csv(data_path2/'train-rle.csv', index_col=0)
df.head()
Out[0]:
EncodedPixels
ImageId
1.2.276.0.7230010.3.1.4.8323329.5597.1517875188.959090-1
1.2.276.0.7230010.3.1.4.8323329.12515.1517875239.501137-1
1.2.276.0.7230010.3.1.4.8323329.4904.1517875185.355709175349 7 1013 12 1009 17 1005 19 1003 20 1002...
1.2.276.0.7230010.3.1.4.8323329.32579.1517875161.299312407576 2 1021 7 1015 10 1013 12 1011 14 1008 ...
1.2.276.0.7230010.3.1.4.8323329.32579.1517875161.299312252069 1 1021 3 1020 4 1018 5 1018 6 1016 7 1...
In [0]:
df.shape
Out[0]:
(11582, 1)

4.2 Mask

ใช้ฟังก์ชันที่ได้มา แปลง RLE เป็น Mask ขนาด 1024 x 1024

In [0]:
mask = rle2mask(df.loc['1.2.276.0.7230010.3.1.4.8323329.4904.1517875185.355709'][0].strip(), 1024, 1024).T
# a = mask.sum(axis=1)
# a = mask.sum(axis=0)
# a.argmax()

mask[105:115, 300:310]
Out[0]:
array([[  0.,   0.,   0.,   0., ...,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0., ...,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0., ...,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0., ...,   0.,   0.,   0.,   0.],
       ...,
       [  0.,   0.,   0.,   0., ..., 255., 255., 255., 255.],
       [  0.,   0.,   0., 255., ..., 255., 255., 255., 255.],
       [  0., 255., 255., 255., ..., 255., 255., 255., 255.],
       [255., 255., 255., 255., ..., 255., 255., 255., 255.]])

เพื่อให้เห็นชัด ๆ เราจะแปลงเป็น Transparent ไว้สำหรับ Overlay กับภาพฟิล์ม X-Ray

In [0]:
mask_rgba = mask/255
mask_rgba = np.repeat(mask_rgba[:, :, np.newaxis], 4, axis=2)
mask_rgba[:, :, [1,2]] = 0
mask_rgba[:, :, 3] = mask_rgba[:, :, 3]*0.7

mask_rgba[105:115, 300:310, 0]
Out[0]:
array([[0., 0., 0., 0., ..., 0., 0., 0., 0.],
       [0., 0., 0., 0., ..., 0., 0., 0., 0.],
       [0., 0., 0., 0., ..., 0., 0., 0., 0.],
       [0., 0., 0., 0., ..., 0., 0., 0., 0.],
       ...,
       [0., 0., 0., 0., ..., 1., 1., 1., 1.],
       [0., 0., 0., 1., ..., 1., 1., 1., 1.],
       [0., 1., 1., 1., ..., 1., 1., 1., 1.],
       [1., 1., 1., 1., ..., 1., 1., 1., 1.]])

ลองแสดง Mask เป็นรูป

In [0]:
plt.imshow(mask_rgba)
Out[0]:
<matplotlib.image.AxesImage at 0x7f0f6231b7f0>

4.3 DICOM Image and Mask

แสดงฟิล์ม X-Ray พร้อมระบายสี Pneumothorax ด้วย Mask สีแดง

In [0]:
start = 44   # Starting index of images
num_img = 3 # Total number of images to show
In [0]:
fig, ax = plt.subplots(nrows=1, ncols=num_img, sharey=True, figsize=(num_img*5,5))
for q, file_path in enumerate(glob.glob(str(data_path2/'pneumothorax/dicom-images-train/*/*/*.dcm'))[start:start+num_img]):
    dataset = pydicom.dcmread(file_path)
    print(file_path.split('/')[-1])
    ax[q].imshow(dataset.pixel_array, cmap=plt.cm.bone)
    fn = file_path.split('/')[-1][:-4]
    rle = df.loc[fn][0].strip()
    # print(rle)
    if rle != '-1':
        print('See Marker')
        mask = rle2mask(rle, 1024, 1024).T
        ax[q].set_title('See Marker')
#         # Normal Mask
#         ax[q].imshow(mask, alpha=0.3, cmap="rainbow")
        
        # Make alpha mask
        mask_rgba = mask/255
        mask_rgba = np.repeat(mask_rgba[:, :, np.newaxis], 4, axis=2)
        mask_rgba[:, :, [1,2]] = 0
        mask_rgba[:, :, 3] = mask_rgba[:, :, 3]*0.4

        ax[q].imshow(mask_rgba)
    else:
        print('Nothing to see')
        ax[q].set_title('Nothing to see')
1.2.276.0.7230010.3.1.4.8323329.10359.1517875223.342952.dcm
See Marker
1.2.276.0.7230010.3.1.4.8323329.321.1517875162.374614.dcm
Nothing to see
1.2.276.0.7230010.3.1.4.8323329.3666.1517875178.910556.dcm
See Marker