Kaggle lung cancer detection – Phác thảo hệ thống (Prototype)

compressed-images

compressed-images

Một Data scientist cần có khả năng prototype nhanh mô hình dự đoán của mình bằng cách sử dụng mẫu dữ liệu nhỏ có thể lưu trữ ngay trên laptop. Khi mô hình đã được kiểm chứng và hoạt động, ta có thể tiến hành làm việc với các thành viên khác trong nhóm để tích hợp các tác vụ rút trích đặc trưng (feature extraction), quản lý và mở rộng prototype được viết bởi nhiều thành viên khác trong nhóm. Sau khi hoàn tất phần prototype, ta có thể làm việc với engineers/developers để hiện thực hóa sản phẩm thông qua mô hình đã huấn luyện.

Khi prototype, ta hoàn toàn có quyền viết scripts không trật tự và rõ ràng để hoàn tất công việc, nhưng cần đảm bảo code được viết càng đơn giản càng tốt để sau này có thể đọc hiểu và phát triển bởi các thành viên khác.

Tiếp tục với cuộc thi Kaggle lung cancer detection, trong bài viết này, ta sẽ cùng nhau prototype hệ thống chẩn đoán ung thư phổi đơn giản. Đầu vào là ma trận đặc trưng (sử dụng ngay ảnh raw, chưa áp dụng các kỹ thuật rút trích đặc trưng). Đầu ra là kết quả đánh giá và so sánh giữa các mô hình dự đoán.

Github sources:

Chuyển đổi dữ liệu ảnh

Đầu tiên, ta cần chuyển đổi chuỗi ảnh các lát cắt thành một ảnh duy nhất. Ý tưởng chính

  • Duyệt từng thư mục các bệnh nhân
    • Khởi tạo biến images để lưu danh sách ảnh
    • Đối với mỗi thư mục bệnh nhân, đọc ra danh sách các file .dcm
      • Chuẩn hóa giá trị ảnh vào khoảng [0, 255]
      • Áp dụng bộ lọc meadian để làm trơn ảnh
      • Phân đoạn ảnh
      • Thêm vào images ảnh đã qua xử lý
    • Tính trung bình ảnh (cộng giá trị của các ảnh và lấy trung bình)
    • Lưu lại ảnh với định dạng .png

Sau khi hoàn tất, ta sẽ có danh sách ảnh phổi cho từng bệnh nhân như bên dưới

processed-images

processed-images

Để tiến hành cài đặt, bạn cần các thư viện bên dưới

import argparse
import datetime
import glob
import os
import sys
import time

import numpy as np
# fix lỗi không nhận ra module cv2
sys.path.append("/usr/local/lib/python2.7/site-packages")
import cv2
import dicom as dicomio

Trong đó, thư viện xử lý ảnh quen thuộc OpenCV, các bạn có thể tham khảo cách cài đặt tại đây. Ta thiết lập để nhận tham số từ dòng lệnh, việc làm này khá hữu ích, nó giúp cho script của bạn không bị viết cứng (hard coding).

# construct the argument parse and parse the arguments
ap = argparse.ArgumentParser()
ap.add_argument("-d", "--dataset", required=True, help="path to input dataset")
ap.add_argument("-s", "--saveto", required=True, help="path to saved processed data")
args = vars(ap.parse_args())

Giả sử bạn tạo một script tên là image_segmentation.py, ta thực hiện chạy trên dòng lệnh như sau:

python image_segmentation.py --dataset=<path_to_folder_dataset> --saveto=<path_to_save_dataset>

Tiếp đến, ta hiện thực hóa ý tưởng như đoạn script bên dưới:

list_dir = os.listdir(args["dataset"])
for (i, dir) in enumerate(list_dir):
    if os.path.isfile(dir) is False:
        basePath = args["dataset"] + "/" + dir
        os.chdir(basePath)
        images = []
        for f in glob.glob("*.dcm"):
            # read dcm file
            ds = dicomio.read_file(f)
            img = ds.pixel_array

            # normalize image values to [0, 255]
            cv2.normalize(img, img, 0, 255, cv2.NORM_MINMAX)
            img = cv2.medianBlur(img.astype(np.uint8), 5)

            # image segmentation
            thresh = cv2.adaptiveThreshold(img,
                                           255,
                                           cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
                                           cv2.THRESH_BINARY,
                                           11,
                                           2)
            images.append(thresh)

        data = np.array(images)
        mean_img = np.mean(data, axis=0)
        save_name = args["saveto"] + "/" + dir + ".png"
        cv2.imwrite(save_name, mean_img)
        print "Saved processed image:", save_name

Do dung lượng dữ liệu khá lớn, nên thời gian xử lý ảnh khá lâu. Ta có thể thêm vào vài dòng để theo dõi thời gian xử lý.

def time_diff_str(t1, t2):
    """
    Calculates time durations.
    """
    diff = t2 - t1
    mins = int(diff / 60)
    secs = round(diff % 60, 2)
    return str(mins) + " mins and " + str(secs) + " seconds"

t_start = time.time()

# process here ...

print "[INFO]", datetime.datetime.now(), "* DONE After *", time_diff_str(t_start, time.time())

Sau khi viết lại, ta sẽ được đoạn script như bên dưới

t_start = time.time()

# construct the argument parse and parse the arguments
ap = argparse.ArgumentParser()
ap.add_argument("-d", "--dataset", required=True, help="path to input dataset")
ap.add_argument("-s", "--saveto", required=True, help="path to saved processed data")
args = vars(ap.parse_args())

list_dir = os.listdir(args["dataset"])
for (i, dir) in enumerate(list_dir):
    if os.path.isfile(dir) is False:
        basePath = args["dataset"] + "/" + dir
        os.chdir(basePath)
        images = []
        for f in glob.glob("*.dcm"):
            # read dcm file
            ds = dicomio.read_file(f)
            img = ds.pixel_array

            # normalize image values to [0, 255]
            cv2.normalize(img, img, 0, 255, cv2.NORM_MINMAX)
            img = cv2.medianBlur(img.astype(np.uint8), 5)

            # image segmentation
            thresh = cv2.adaptiveThreshold(img,
                                           255,
                                           cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
                                           cv2.THRESH_BINARY,
                                           11,
                                           2)
            images.append(thresh)

        data = np.array(images)
        mean_img = np.mean(data, axis=0)
        save_name = args["saveto"] + "/" + dir + ".png"
        cv2.imwrite(save_name, mean_img)
        print "Saved processed image:", save_name

    # show an update every 10 patients
    if i > 0 and i % 10 == 0:
        print "[INFO] processed {}/{} patients".format(i, len(list_dir))
        print "[INFO] time passed", time_diff_str(t_start, time.time())

print "[INFO]", datetime.datetime.now(), "* DONE After *", time_diff_str(t_start, time.time())

Đối với tập dữ liệu gồm 1595 bệnh nhân, thời gian xử lý của chúng ta khoảng 1 giờ đồng hồ. Tập dữ liệu đã qua xử lý có thể download ở đây.

[INFO] processed 1590/1595 patients
[INFO] time passed 56 mins and 42.61 seconds
Saved processed image: /Users/hongong/Downloads/processed/fe5c37e82b412833b8ad0abb57978377.png
Saved processed image: /Users/hongong/Downloads/processed/ff5d8e90500cf324e7b04a2f07cf0399.png
Saved processed image: /Users/hongong/Downloads/processed/ff8599dd7c1139be3bad5a0351ab749a.png
Saved processed image: /Users/hongong/Downloads/processed/ffe02fe7d2223743f7fb455dfaff3842.png
0
[INFO] 2017-02-09 16:13:22.629087 * DONE After * 57 mins and 12.51 seconds

Chuẩn bị đầu vào

raw-features

raw-features

Ta cần chuẩn bị ma trận N x (M + 1). Trong đó, N là số lượng quan sát/dòng dữ liệu/bệnh nhân (trong bài toán này). M là số features/thuộc tính/cột dữ liệu/giá trị màu của ảnh,… cộng thêm 1 giá trị là nhãn phân loại của quan sát đó (có ung thu: 1, không có ung thư: 0).

Trong 2 files stage1_labels.csv và stage1_sample_submission.csv có nội dung như bên dưới dùng để xác định bệnh nhân đang khảo sát có nguy cơ bị ung thư hay không:

id,cancer
0015ceb851d7251b8f399e39779d1e7d,1
0030a160d58723ff36d73f41b170ec21,0
003f41c78e6acfa92430a057ac0b306e,0
006b96310a37b36cccb2ab48d10b49a3,1
008464bb8521d09a42985dd8add3d0d2,1
0092c13f9e00a3717fdc940641f00015,0
00986bebc45e12038ef0ce3e9962b51a,0
00cba091fa4ad62cc3200a657aeb957e,0
00edff4f51a893d80dae2d42a7f45ad1,1
0121c2845f2b7df060945b072b2515d7,0
013395589c01aa01f8df81d80fb0e2b8,0
01de8323fa065a8963533c4a86f2f6c1,0
01e349d34c06410e1da273add27be25c,0
01f1140c8e951e2a921b61c9a7e782c2,0
024efb7a1e67dc820eb61cbdaa090166,0
...

Ta sẽ load nội dung 2 files này vào một biến kiểu dữ liệu Dataframe (key là patient id -> value là nhãn phân loại 0|1).

def load_csv(file_path):
    """Get data, from local csv."""
    if os.path.exists(file_path):
        print "[INFO] load", file_path, "file..."
        df = pd.read_csv(file_path)

    return df

# load train/test labels
stage1_labels = load_csv("../data/stage1_labels.csv")
stage1_sample_submission = load_csv("../data/stage1_sample_submission.csv")

Tiếp theo, ta sẽ xây dựng ma trận N x (M + 1) hay #patients x (#pixels + label value).

def image_to_feature_vector(image, size=(32, 32)):
    # resize the image to a fixed size, then flatten the image into
    # a list of raw pixel intensities
    return cv2.resize(image, size).flatten()


def get_simple_feature_labels(patient_df, img_paths):
    features = []
    labels = []

    patient_ids = patient_df["id"].tolist()

    # loop over the input images
    for (i, img_path) in enumerate(img_paths):
        # get only training labels
        base = os.path.basename(img_path)
        patient_id = os.path.splitext(base)[0]
        if patient_id in patient_ids:
            label = patient_df[patient_df["id"] == patient_id].iloc[0]["cancer"]
            labels.append(label)
        else:
            continue

        # load the image
        image = cv2.imread(img_path)
        feat = image_to_feature_vector(image)

        # update features
        features.append(feat)

        # show an update every 100 images
        if i > 0 and i % 100 == 0:
            print("[INFO] processed {}/{}".format(i, len(img_paths)))

    return features, labels

img_paths = list(paths.list_images(args["dataset"]))
train_features, train_labels = get_simple_feature_labels(stage1_labels, img_paths)
test_features, test_labels = get_simple_feature_labels(stage1_sample_submission, img_paths)

Ngoài ra, sau khi hoàn tất, ta có thể khảo sát một xíu về ma trận thu được về kích thước hay hình ảnh trực quan như thế nào.

train_labels = np.array(train_labels)
train_features = np.array(train_features)
print "[INFO] labels vector shape:", train_labels.shape
print "[INFO] features matrix shape:", train_features.shape
print("[INFO] features matrix size: {:.2f}MB".format(train_features.nbytes / (1024 * 1000.0)))

plt.imshow(train_features)
plt.show()

Cuối cùng, ta sẽ có kết quả quan sát như sau:

[INFO] load ../data/stage1_labels.csv file...
[INFO] load ../data/stage1_sample_submission.csv file...
[INFO] processed 100/1595
[INFO] processed 200/1595
[INFO] processed 300/1595
[INFO] processed 400/1595
[INFO] processed 500/1595
[INFO] processed 600/1595
[INFO] processed 700/1595
[INFO] processed 800/1595
[INFO] processed 900/1595
[INFO] processed 1000/1595
[INFO] processed 1100/1595
[INFO] processed 1200/1595
[INFO] processed 1300/1595
[INFO] processed 1400/1595
[INFO] processed 1500/1595
[INFO] labels vector shape: (1397,)
[INFO] features matrix shape: (1397, 3072)
[INFO] features matrix size: 4.19MB
features-labels-matrix

features-labels-matrix

Tiến hành huấn luyện

Nếu bạn đã đọc qua các bài viết liên quan đến bài toán classification thì ở bước này khá đơn giản. Đầu tiên, ta cần nạp các thư viện huấn luyện cần thiết.

import argparse
import csv
import datetime
import os
import sys
import time
from operator import itemgetter

import numpy as np
import pandas as pd
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.ensemble import AdaBoostClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier

Thêm vào đó, bạn có thể quyết định bao nhiêu core bộ xử lý tham gia vào quá trình huấn luyện giúp tăng cường tốc độ tính toán ở một số mô hình.

ap.add_argument("-j", "--jobs", type=int, default=-1, help="# of jobs (-1 uses all available cores)")

print "---------------------------"
print "Training"
print "---------------------------"

classifiers = {
    "Nearest Neighbors": KNeighborsClassifier(3, n_jobs=args["jobs"]),
    "Linear SVM": SVC(kernel="linear", C=0.025),
    "RBF SVM": SVC(gamma=2, C=1),
    "Gaussian Process": GaussianProcessClassifier(1.0 * RBF(1.0), warm_start=True, n_jobs=args["jobs"]),
    "Decision Tree": DecisionTreeClassifier(max_depth=5),
    "Random Forest": RandomForestClassifier(max_depth=5, n_estimators=10, max_features=1, n_jobs=args["jobs"]),
    "Neural Net": MLPClassifier(alpha=1),
    "AdaBoost": AdaBoostClassifier(),
    "Naive Bayes": GaussianNB(),
    "QDA": QuadraticDiscriminantAnalysis()
}

Ta tiến hành phân chia tập dữ liệu ban đầu thành tập train và dev (do tập dữ liệu test dùng để submit lên Kaggle ta sẽ gọi tập dữ liệu dùng để đánh giá mô hình là tập dev). Cuối cùng, ta sẽ lưu giá trị đánh giá của từng mô hình, sắp xếp và hiển thị kết quả đánh giá.

# train/dev split
# X_train, X_dev, Y_train, Y_dev
(for_train_features, dev_features, for_train_labels, dev_labels) = train_test_split(train_features,
                                                                                    train_labels,
                                                                                    test_size=0.25,
                                                                                    random_state=42)

# iterate over classifiers
results = {}

for name in classifiers:
    print "[INFO]" + name + " classifier..."
    clf = classifiers[name]
    clf.fit(for_train_features, for_train_labels)
    score = clf.score(dev_features, dev_labels)
    results[name] = score

print "---------------------------"
print "Evaluation results"
print "---------------------------"

# sorting results and print out
sorted(results.items(), key=itemgetter(1))
for name in results:
    print "[INFO]", name, "accuracy: %0.3f" % results[name]

Output:

---------------------------
Training
---------------------------
[INFO]Neural Net classifier...
[INFO]Gaussian Process classifier...
[INFO]Decision Tree classifier...
[INFO]Random Forest classifier...
[INFO]AdaBoost classifier...
[INFO]QDA classifier...
[INFO]Naive Bayes classifier...
[INFO]Linear SVM classifier...
[INFO]Nearest Neighbors classifier...
[INFO]RBF SVM classifier...
---------------------------
Evaluation results
---------------------------
[INFO] Gaussian Process accuracy: 0.760
[INFO] QDA accuracy: 0.531
[INFO] Decision Tree accuracy: 0.734
[INFO] Naive Bayes accuracy: 0.657
[INFO] Linear SVM accuracy: 0.589
[INFO] Neural Net accuracy: 0.711
[INFO] RBF SVM accuracy: 0.760
[INFO] AdaBoost accuracy: 0.683
[INFO] Random Forest accuracy: 0.760
[INFO] Nearest Neighbors accuracy: 0.689

Đối với feature và các mô hình dự đoán đã sử dụng, ta có thể thấy mô hình Gausian cho kết quả cao nhất khoảng 76%. Ta sẽ sử dụng mô hình này để huấn luyện và dự đoán kết quả cho tập test. Sau đó, ta sẽ upload kết quả này lên Kaggle.

print "---------------------------"
print "Training for submission"
print "---------------------------"

name = list(results)[0]
clf = classifiers[name]
print "[INFO]" + name + " classifier..."
clf.fit(train_features, train_labels)
predict_submission = clf.predict(test_features)

# update submission
submission = {}
patient_ids = stage1_sample_submission["id"].tolist()
for (i, patient_id) in enumerate(patient_ids):
    submission[patient_id] = predict_submission[i]

with open(args["out"], "wb") as f:
    writer = csv.writer(f, delimiter=',')
    writer.writerow(["id", "cancer"])
    for key, value in submission.items():
        writer.writerow([key, value])

---------------------------
Training for submission
---------------------------
[INFO]Gaussian Process classifier...
[INFO] 2017-02-12 16:35:57.880718 * DONE After * 1 mins and 16.71 seconds

Sau khi submit kết quả dự đoán lên Kaggle, ta không quá ngạc nhiên khi kết quả dự đoán không cao. Tuy nhiên, ta đã đạt được mục tiêu ban đầu đó là phác thảo được một hệ thống phân loại đơn giản để có thể tiếp tục cải tiến mô hình về sau.

kaggle-submission-01

kaggle-submission-01

Advertisements

2 thoughts on “Kaggle lung cancer detection – Phác thảo hệ thống (Prototype)

  1. Cám ơn bạn, bài chia sẽ rất chi tiết.
    Mình chỉ comment 1 chút đoạn này:

    # train/dev split
    (for_train_features, dev_features, for_train_labels, dev_labels) = train_test_split(…)

    Mình nghĩ nên dùng cross-validation sẽ evaluate model chính xác hơn (K-Fold hoặc Stratified fold)

    Liked by 1 person

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s