ในการพัฒนาระบบ AI สำหรับวินิจฉัยโรค ชนิดของข้อมูล Dataset ที่เราจะพบบ่อย ๆ หนึ่งในนั้นคือไฟล์ DICOM ซึ่งเป็นไฟล์ภาพถ่ายทางการแพทย์ Mdical Imaging สำหรับฟิล์ม X-Ray, MRI, CT-Scan, Ultrasound ใน ep นี้เราจะมาทำความเข้าใจ ไฟล์ DICOM กันให้มากขึ้น

DICOM คืออะไร

DICOM Logo Digital Imaging and Communications in Medicine
DICOM Logo Digital Imaging and Communications in Medicine

DICOM ย่อมาจาก Digital Imaging and Communications in Medicine เป็นมาตรฐานกลางในการสื่อสาร จัดเก็บ เรียกดู ประมวลผล จัดพิมพ์ และแสดงผลข้อมูลภาพถ่ายทางการแพทย์ Mdical Imaging กำหนดโดย องค์กรชื่อ National Electrical Manufacturers Association (NEMA)

มาตรฐาน DICOM ทำให้เราสามารถใช้ข้อมูลภาพถ่ายทางการแพทย์ Mdical Imaging ต่าง ๆ ร่วมกันได้ข้ามระบบ ระหว่างอุปกรณ์ต่าง ๆ ตั้งแต่ อุปกรณ์ถ่ายภาพ, PACS (Picture Archiving and Communication System), Workstation, VNAs และเครื่องพิมพ์ จากหลากหลายยี่ห้อ

Whole Slide Imaging Scanner. Credit https://commons.wikimedia.org/wiki/File:WSI_Scanner.jpg
Whole Slide Imaging Scanner. Credit https://commons.wikimedia.org/wiki/File:WSI_Scanner.jpg

ภายในไฟล์จะประกอบด้วยทั้งส่วน Header ที่เป็น Metada และส่วนของข้อมูลรูปภาพ ซึ่งมีได้มากกว่า 1 รูปต่อ 1 ไฟล์ เรามาดูตัวอย่างไฟล์ และรูปแบบข้อมูลในไฟล์กันดีกว่า

Search Results Web result with site links  SIIM-ACR Pneumothorax Classification Dataset
Search Results Web result with site links SIIM-ACR Pneumothorax Classification Dataset

Weighted Cross Entropy Loss

ในเคสนี้จำนวนข้อมูลตัวอย่าง ในแต่ละ Class แตกต่างกันมาก เรียกว่า Class Imbalance เราจะใช้ Loss Function แบบใหม่ ชื่อว่า Weighted Cross Entropy Loss

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

Open In Colab

ในการพัฒนาระบบ AI สำหรับวินิจฉัยโรค ชนิดของข้อมูล Dataset ที่เราจะพบบ่อย ๆ หนึ่งในนั้นคือไฟล์ DICOM ซึ่งเป็นไฟล์สำหรับฟิล์มเ X-Ray, MRI, CT-Scan, Ultrasound ใน ep นี้เราจะมาทำความเข้าใจ ไฟล์ DICOM กันให้มากขึ้น

DICOM คืออะไร

DICOM ย่อมาจาก Digital Imaging and Communications in Medicine เป็นมาตรฐานกลางในการสื่อสาร จัดเก็บ เรียกดู ประมวลผล จัดพิมพ์ และแสดงผลข้อมูลรูปภาพทางการแพทย์ กำหนดโดย องค์กรชื่อ National Electrical Manufacturers Association (NEMA)

มาตรฐาน DICOM ทำให้เราสามารถใช้ข้อมูลต่าง ๆ ร่วมกันได้ข้ามระบบ ระหว่างอุปกรณ์ต่าง ๆ ตั้งแต่ อุปกรณ์ถ่ายภาพ, PACS (Picture Archiving and Communication System), Workstation, VNAs และเครื่องพิมพ์ จากหลากหลายยี่ห้อ

ภายในไฟล์จะประกอบด้วยทั้งส่วน Header ที่เป็น Metada และส่วนของข้อมูลรูปภาพ ซึ่งมีได้มากกว่า 1 รูปต่อ 1 ไฟล์ เรามาดูตัวอย่างไฟล์ และรูปแบบข้อมูลในไฟล์กันดีกว่า

เช็ค GPU

In [1]:
! nvidia-smi
Wed Jun 24 04:12:18 2020       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 450.36.06    Driver Version: 418.67       CUDA Version: 10.1     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  Tesla P100-PCIE...  Off  | 00000000:00:04.0 Off |                    0 |
| N/A   36C    P0    26W / 250W |      0MiB / 16280MiB |      0%      Default |
|                               |                      |                 ERR! |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|  No running processes found                                                 |
+-----------------------------------------------------------------------------+

0. Magic Command

In [2]:
%reload_ext autoreload
%autoreload 2
%matplotlib inline

1. Install Library

Install Library ที่จำเป็น ในที่นี้เราจะใช้ โมดูล medical.imaging ของ fastai2

In [9]:
# ## Colab
# ! pip install fastai2 pydicom kornia -q
In [10]:
# ## Colab

# !pip install torch-geometric \
#   torch-sparse==latest+cu101 \
#   torch-scatter==latest+cu101 \
#   torch-cluster==latest+cu101 \
#   -f https://pytorch-geometric.com/whl/torch-1.4.0.html

2. Import Library

Import Library ที่จำเป็น

In [5]:
from fastai2.basics import *
from fastai2.callback.all import *
from fastai2.metrics import *
from fastai2.vision.all import *
from fastai2.medical.imaging import *

import pydicom
import kornia

import pandas as pd
In [6]:
set_seed(123456)

3. Dataset

ในเคสนี้เราจะใช้ Dataset ฟิล์ม X-Ray จาก Kaggle เหมือนใน ep ก่อน พัฒนาโปรแกรม AI การแพทย์ วินิจฉัยภาวะปอดรั่ว Pneumothorax

เราจะ Mount Drive ไปยัง Google Drive ที่เก็บ Token File ไว้

In [7]:
dataset = 'jesperdramsch/siim-acr-pneumothorax-segmentation-data'

# Google Colab
config_path = Path('/content/drive')
learner_path = config_path/"My Drive"
data_path_base = Path('/content/datasets/')

path = data_path_base/dataset

from google.colab import drive

drive.mount(str(config_path))
os.environ['KAGGLE_CONFIG_DIR'] = f"{config_path}/My Drive/.kaggle"
Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&response_type=code&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly

Enter your authorization code:
··········
Mounted at /content/drive

ในการจะ Download ข้อมูลจาก Kaggle ต้องใช้ Token ดูวิธีได้ใน ep ก่อน

In [8]:
# !kaggle datasets download {dataset} -p "{path}" --unzip
Downloading siim-acr-pneumothorax-segmentation-data.zip to /content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data
100% 3.03G/3.03G [00:50<00:00, 22.4MB/s]
100% 3.03G/3.03G [00:50<00:00, 63.9MB/s]

ls ดูว่าได้ Folder อะไรมาบ้าง

In [11]:
path.ls()
Out[11]:
(#4) [Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/pneumothorax'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/train-rle.csv'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-test')]

ดูข้อมูลใน Training Folder

In [12]:
(path/'dicom-images-train').ls()
Out[12]:
(#10712) [Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.1243.1517875166.938628'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.32404.1517875160.439331'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.12884.1517875242.357985'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.3891.1517875179.987983'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.32176.1517875158.889266'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.13369.1517875244.873269'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.12298.1517875238.439592'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.10206.1517875222.624432'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.440.1517875162.950677'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.4008.1517875180.621634')...]

ใช้ฟังก์ชัน get_dicom_files ดึงไฟล์ทั้งหมดมาใส่ List ไว้ก่อน

In [13]:
items = get_dicom_files(path/'dicom-images-train')
items
Out[13]:
(#10712) [Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.1243.1517875166.938628/1.2.276.0.7230010.3.1.3.8323329.1243.1517875166.938627/1.2.276.0.7230010.3.1.4.8323329.1243.1517875166.938629.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.32404.1517875160.439331/1.2.276.0.7230010.3.1.3.8323329.32404.1517875160.439330/1.2.276.0.7230010.3.1.4.8323329.32404.1517875160.439332.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.12884.1517875242.357985/1.2.276.0.7230010.3.1.3.8323329.12884.1517875242.357984/1.2.276.0.7230010.3.1.4.8323329.12884.1517875242.357986.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.3891.1517875179.987983/1.2.276.0.7230010.3.1.3.8323329.3891.1517875179.987982/1.2.276.0.7230010.3.1.4.8323329.3891.1517875179.987984.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.32176.1517875158.889266/1.2.276.0.7230010.3.1.3.8323329.32176.1517875158.889265/1.2.276.0.7230010.3.1.4.8323329.32176.1517875158.889267.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.13369.1517875244.873269/1.2.276.0.7230010.3.1.3.8323329.13369.1517875244.873268/1.2.276.0.7230010.3.1.4.8323329.13369.1517875244.873270.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.12298.1517875238.439592/1.2.276.0.7230010.3.1.3.8323329.12298.1517875238.439591/1.2.276.0.7230010.3.1.4.8323329.12298.1517875238.439593.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.10206.1517875222.624432/1.2.276.0.7230010.3.1.3.8323329.10206.1517875222.624431/1.2.276.0.7230010.3.1.4.8323329.10206.1517875222.624433.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.440.1517875162.950677/1.2.276.0.7230010.3.1.3.8323329.440.1517875162.950676/1.2.276.0.7230010.3.1.4.8323329.440.1517875162.950678.dcm'),Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.4008.1517875180.621634/1.2.276.0.7230010.3.1.3.8323329.4008.1517875180.621633/1.2.276.0.7230010.3.1.4.8323329.4008.1517875180.621635.dcm')...]

4. DICOM File

เลือกคนไข้ขึ้นมา 1 คน แล้วใช้ฟังก์ชัน dcmread เปิดไฟล์ขึ้นมาแสดงผลดูรูปภาพภายใน ไฟล์ DICOM

In [14]:
patient = 123
items[patient]
Out[14]:
Path('/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.2615.1517875173.725632/1.2.276.0.7230010.3.1.3.8323329.2615.1517875173.725631/1.2.276.0.7230010.3.1.4.8323329.2615.1517875173.725633.dcm')

4.1 Image

ดูรูปภาพ

In [15]:
xray_sample = dcmread(items[patient])
xray_sample.show(figsize=(12, 12))

4.2 Metadata

ดู Metadata ในไฟล์ DICOM

In [16]:
xray_sample.as_dict()
Out[16]:
{'AccessionNumber': '',
 'BitsAllocated': 8,
 'BitsStored': 8,
 'BodyPartExamined': 'CHEST',
 'Columns': 1024,
 'ConversionType': 'WSD',
 'HighBit': 7,
 'InstanceNumber': 1,
 'LossyImageCompression': '01',
 'LossyImageCompressionMethod': 'ISO_10918_1',
 'Modality': 'CR',
 'MultiPixelSpacing': 1,
 'PatientAge': '68',
 'PatientBirthDate': '',
 'PatientID': '55362559-978b-489f-b851-e1279bd40f61',
 'PatientName': '55362559-978b-489f-b851-e1279bd40f61',
 'PatientOrientation': '',
 'PatientSex': 'F',
 'PhotometricInterpretation': 'MONOCHROME2',
 'PixelRepresentation': 0,
 'PixelSpacing': 0.139,
 'PixelSpacing1': 0.139,
 'ReferringPhysicianName': '',
 'Rows': 1024,
 'SOPClassUID': '1.2.840.10008.5.1.4.1.1.7',
 'SOPInstanceUID': '1.2.276.0.7230010.3.1.4.8323329.2615.1517875173.725633',
 'SamplesPerPixel': 1,
 'SeriesDescription': 'view: AP',
 'SeriesInstanceUID': '1.2.276.0.7230010.3.1.3.8323329.2615.1517875173.725631',
 'SeriesNumber': 1,
 'SpecificCharacterSet': 'ISO_IR 100',
 'StudyDate': '19010101',
 'StudyID': '',
 'StudyInstanceUID': '1.2.276.0.7230010.3.1.2.8323329.2615.1517875173.725632',
 'StudyTime': '000000.00',
 'ViewPosition': 'AP',
 'fname': '/content/datasets/jesperdramsch/siim-acr-pneumothorax-segmentation-data/dicom-images-train/1.2.276.0.7230010.3.1.2.8323329.2615.1517875173.725632/1.2.276.0.7230010.3.1.3.8323329.2615.1517875173.725631/1.2.276.0.7230010.3.1.4.8323329.2615.1517875173.725633.dcm',
 'img_max': 255,
 'img_mean': 126.84897518157959,
 'img_min': 0,
 'img_pct_window': 0.22364425659179688,
 'img_std': 64.46480392944859}

จะเห็นว่า มีหลายตัวที่น่าสนใจ เช่น

  • BitsStored เป็นตัวบอกว่ารูปภาพกี่ Bit อาจจะเป็น 8, 14, 16 Bit การแปลงเป็น 8 Bit x 3 Channel RGB เพื่อแสดงผลบนจอภาพ อาจจะทำให้สูญเสียข้อมูลได้
  • PixelRepresentation บอกว่าข้อมูลเป็น Signed หรือ Unsigned
  • img_max, img_mean, img_min บอกสถิติของข้อมูล Pixel

4.3 Pixels Data

ในเคสนี้ ขนาดของรูปภาพคือ 1024 x 1024 Pixel

In [17]:
xray_sample.pixels.shape
Out[17]:
torch.Size([1024, 1024])

เป็นรูปภาพขาวดำ Grayscale 0-255

In [18]:
xray_sample.pixels.min(), xray_sample.pixels.max(), xray_sample.scaled_px.min(), xray_sample.scaled_px.max()
Out[18]:
(tensor(0.), tensor(255.), tensor(0.), tensor(255.))
In [19]:
xray_sample.pixels[300:310, 400:410]
Out[19]:
tensor([[ 90.,  91.,  89.,  86.,  85.,  85.,  85.,  84.,  83.,  70.],
        [ 88.,  91.,  92.,  88.,  85.,  84.,  85.,  84.,  81.,  72.],
        [ 85.,  90.,  92.,  89.,  86.,  85.,  84.,  83.,  79.,  74.],
        [ 84.,  89.,  92.,  90.,  88.,  87.,  84.,  81.,  76.,  75.],
        [ 86.,  85.,  87.,  91.,  94.,  95.,  89.,  79.,  77.,  79.],
        [ 87.,  87.,  91.,  94.,  95.,  94.,  89.,  81.,  79.,  83.],
        [ 86.,  86.,  91.,  94.,  94.,  93.,  89.,  83.,  83.,  87.],
        [ 86.,  84.,  86.,  91.,  94.,  96.,  93.,  86.,  87.,  89.],
        [ 91.,  84.,  83.,  89.,  96., 100.,  96.,  88.,  87.,  88.],
        [ 96.,  87.,  84.,  89.,  96.,  98.,  93.,  85.,  85.,  85.]])
In [20]:
xray_sample.scaled_px[300:310, 400:410]
Out[20]:
tensor([[ 90.,  91.,  89.,  86.,  85.,  85.,  85.,  84.,  83.,  70.],
        [ 88.,  91.,  92.,  88.,  85.,  84.,  85.,  84.,  81.,  72.],
        [ 85.,  90.,  92.,  89.,  86.,  85.,  84.,  83.,  79.,  74.],
        [ 84.,  89.,  92.,  90.,  88.,  87.,  84.,  81.,  76.,  75.],
        [ 86.,  85.,  87.,  91.,  94.,  95.,  89.,  79.,  77.,  79.],
        [ 87.,  87.,  91.,  94.,  95.,  94.,  89.,  81.,  79.,  83.],
        [ 86.,  86.,  91.,  94.,  94.,  93.,  89.,  83.,  83.,  87.],
        [ 86.,  84.,  86.,  91.,  94.,  96.,  93.,  86.,  87.,  89.],
        [ 91.,  84.,  83.,  89.,  96., 100.,  96.,  88.,  87.,  88.],
        [ 96.,  87.,  84.,  89.,  96.,  98.,  93.,  85.,  85.,  85.]])

4.4 Histogram

ลองพล็อต Histogram จะเห็นว่า ข้อมูลมีบางระดับที่มากน้อยต่างกัน บางค่าก็น้อยมาก

In [21]:
hist = plt.hist(xray_sample.pixels.reshape(-1).numpy(), bins=255)

4.5 Normalize

เราสามารถแบ่งข้อมูลออกเป็น Bins ตามความถี่ให้ข้อมูลกระจายออกไปทั่ว ๆ เป็นการ Normalize

In [22]:
bins = xray_sample.pixels.freqhist_bins()
bins.shape, bins
Out[22]:
(torch.Size([99]),
 tensor([  0.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.,  20.,  30.,  38.,
          46.,  53.,  58.,  62.,  66.,  69.,  73.,  76.,  79.,  82.,  85.,  87.,
          90.,  92.,  94.,  96.,  99., 101., 103., 105., 107., 109., 111., 113.,
         115., 117., 119., 120., 122., 123., 125., 126., 128., 129., 130., 132.,
         133., 134., 135., 137., 138., 139., 141., 142., 143., 145., 146., 148.,
         149., 151., 153., 155., 157., 160., 162., 165., 167., 170., 173., 176.,
         180., 183., 186., 189., 191., 194., 196., 198., 200., 202., 203., 205.,
         206., 208., 210., 211., 212., 213., 215., 216., 217., 218., 219., 221.,
         223., 227., 231.]))

แสดง Histogram อีกครั้ง จาก Bins ด้านบน จะเห็นว่ากระจายดีขึ้น โมเดลจะได้ทำงานง่ายขึ้น

In [23]:
hist = plt.hist(xray_sample.pixels.reshape(-1).numpy(), bins=bins)

4.6 DICOM Windows

ในการวิเคราะห์ วินิจฉัยโรค โดยแพทย์ผู้ชำนาญการ เนื่องจากสายตาของมนุษย์ โดยทั่วไปสามารถแยกแยะได้แค่ประมาณ 100 ระดับสีเท่านั้น ในการแสดงผล จึงมีการแยกเฉพาะ Window ในช่วงข้อมูลต่าง ๆ มาแสดงผล ตามความหนาแน่นทางรังสี Radiodensity ตามช่วง Hounsfield scale ของอวัยวะที่เราสนใจ

In [24]:
dicom_windows
Out[24]:
namespace(abdomen_soft=(400, 50), brain=(80, 40), brain_bone=(2800, 600), brain_soft=(375, 40), liver=(150, 30), lungs=(1500, -600), mediastinum=(350, 50), spine_bone=(1800, 400), spine_soft=(250, 50), stroke=(8, 32), subdural=(254, 100))

แสดงข้อมูล ใน Windows ต่าง ๆ เช่น เน้นปอด เน้นกระดูกสันหลัง

In [25]:
scales = False, True, dicom_windows.spine_bone, dicom_windows.lungs
titles = 'raw','normalized','spine_bone windowed','lungs windowed'
for s,a,t in zip(scales, subplots(2,2,imsize=5)[1].flat, titles):
    xray_sample.show(scale=s, ax=a, title=t)