Chuyển đến nội dung chính

Hướng dẫn Super-Resolution ảnh với Khung đa Decoder

Trong hướng dẫn này, bạn sẽ triển khai một bài báo về ảnh y khoa dùng deep learning với Python trong Keras.
· 15 phút đọc

Khám phá cùng AI

Trong phần Training Script, bạn sẽ làm việc với bài toán siêu phân giải ảnh sử dụng một kiến trúc deep learning mới. Nhiệm vụ là ánh xạ phi tuyến từ ảnh MRI não 3 Tesla (từ trường thấp) sang ảnh MRI não 7 Tesla (từ trường cao). Trong phần Testing Script này, bạn sẽ dùng các trọng số đã được huấn luyện ở Phần 1 và dự đoán trên dữ liệu chưa thấy. Bạn cũng sẽ học cách cuối cùng lưu các ảnh 2D thành một khối thể tích kết hợp.

Hướng dẫn được chia thành hai phần, phần đầu tiên sẽ dẫn bạn qua quy trình huấn luyện và phần thứ hai sẽ trình bày quy trình kiểm thử.

Lưu ý: Có thể bạn sẽ muốn đọc bài báo; bạn có thể tìm bài viết tại đây.

Training Script

Tóm lược, trong Phần 1 của hướng dẫn bạn sẽ đề cập các chủ đề sau:

  • Bắt đầu bằng cách import các module cần thiết để huấn luyện mô hình deep learning,
  • Sau đó giới thiệu ngắn về bộ dữ liệu MRI 3T và 7T,
  • Tiếp theo bạn sẽ định nghĩa các initializer và tải bộ dữ liệu 3T và 7T, đồng thời trong quá trình tải sẽ thay đổi kích thước ảnh ngay lập tức,
  • Kế đến, bạn sẽ tiền xử lý dữ liệu đã tải: chuyển danh sách train và test thành ma trận numpy, chuyển kiểu dữ liệu ma trận thành float32, chuẩn tỉ lệ ma trận theo chiến lược max-min, reshape mảng và cuối cùng chia dữ liệu thành 80% huấn luyện và 20% còn lại làm tập validation,
  • Sau đó, bạn sẽ tạo Kiến trúc 1-Encoder-3-Decoder: gồm các kết nối ghép (merge) và nhiều decoder,
  • Tiếp theo, bạn sẽ định nghĩa hàm loss, định nghĩa ba mô hình khác nhau và cuối cùng compile chúng,
  • Cuối cùng, đến lúc bạn sẽ huấn luyện mô hình merge và multi-decoder, kiểm thử trên dữ liệu validation và tính toán kết quả định lượng.

Các phụ thuộc module Python

Trước khi làm theo hướng dẫn, hãy đảm bảo bạn có đúng phiên bản các module như dưới đây:

  • Keras==2.0.4
  • tensorflow==1.8.0
  • scipy==0.19.0
  • numpy==1.14.5
  • Pillow==4.1.1
  • nibabel==2.1.0
  • scikit_learn==0.18.1

Lưu ý: Trước khi bắt đầu, lưu ý rằng mô hình sẽ được huấn luyện trên hệ thống với GPU Nvidia 1080 Ti, CPU Xeon e5 GeForce và RAM 32GB. Nếu bạn dùng Jupyter Notebook, bạn sẽ cần thêm ba dòng mã để chỉ định thứ tự thiết bị CUDA và các thiết bị CUDA khả dụng bằng module os.

Trong đoạn mã dưới, bạn đặt biến môi trường trong notebook bằng os.environ. Nên làm điều này trước khi khởi tạo Keras để giới hạn Keras backend TensorFlow sử dụng GPU đầu tiên. Nếu máy bạn dùng để huấn luyện có GPU ở chỉ số 0, hãy dùng 0 thay vì 1. Bạn có thể kiểm tra bằng lệnh đơn giản trên terminal: ví dụ, nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU

Import các module

Trước tiên, bạn import tất cả module cần thiết như tensorflow, numpy và quan trọng nhất là keras cùng các hàm hoặc layer như Input, Conv2D,MaxPooling2D... vì bạn sẽ dùng chúng để huấn luyện mô hình!
Để đọc ảnh định dạng nifti, bạn cũng cần import module nibabel.

import os
import numpy as np
import math
import tensorflow as tf
import nibabel as nib
import numpy as np
from keras.layers import Input,Dense,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import RMSprop
from keras import backend as K
import scipy.misc
from sklearn.utils import shuffle
from sklearn.cross_validation import train_test_split
import matplotlib.pyplot as plt
from keras.models import model_from_json

Using TensorFlow backend.
/usr/local/lib/python3.5/dist-packages/sklearn/cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also, note that the interface of the new CV iterators is different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)

Tìm hiểu bộ dữ liệu MRI não 3T và 7T

Bộ dữ liệu MRI não 3T và 7T gồm các thể tích 3D, mỗi thể tích có tổng cộng 207 lát/cảnh MRI não chụp ở các lát cắt khác nhau của não. Mỗi lát có kích thước 173 x 173. Ảnh là ảnh xám một kênh. Có tổng cộng 39 đối tượng, mỗi đối tượng chứa ảnh MRI của một bệnh nhân. Định dạng ảnh không phải jpeg, png, v.v. mà là định dạng nifti. Ở phần sau bạn sẽ thấy cách đọc ảnh định dạng nifti.

Bộ dữ liệu gồm ảnh MR thể thức T1; chuỗi T1 truyền thống được xem là tốt cho việc đánh giá cấu trúc giải phẫu. Bộ dữ liệu bạn sẽ làm việc hôm nay gồm MRI não 3T và 7T.

Bộ dữ liệu là công khai và có thể tải tại nguồn này.

28 đối tượng được dùng để huấn luyện, 11 đối tượng còn lại dùng để kiểm thử.

Định nghĩa các initializer

Trước hết, hãy định nghĩa kích thước dữ liệu. Bạn sẽ thay đổi kích thước ảnh từ 173x173 lên 176x176 trong phần đọc dữ liệu. Ở đây bạn cũng định nghĩa thư mục dữ liệu, batch size để huấn luyện mô hình, số kênh, một lớp Input(), ma trận train và test dưới dạng danh sách, và cuối cùng để chuẩn tỉ lệ bạn sẽ tải tệp văn bản chứa giá trị nhỏ nhất và lớn nhất của bộ dữ liệu MRI.

Lưu ý: Với việc chuẩn tỉ lệ, bạn cũng có thể lấy giá trị lớn nhất và nhỏ nhất của bộ dữ liệu vì bạn sẽ không có tệp maxANDmin.txt.

x,y = 173,173
full_z = 207
resizeTo=176
batch_size = 32
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))
inp = "ground3T"
out = "ground7T"     
train_matrix = []
test_matrix = []
min_max = np.loadtxt('maxANDmin.txt')

Tải dữ liệu

Tiếp theo, bạn tải dữ liệu mri bằng thư viện nibabel và thay đổi kích thước ảnh từ 173 x 173 lên 176 x 176 bằng cách đệm số 0 theo chiều x và y.

Lưu ý khi bạn tải một thể tích định dạng Nifti, Nibabel chưa tải mảng ảnh ngay. Nó chờ đến khi bạn yêu cầu dữ liệu mảng. Cách chuẩn để yêu cầu dữ liệu mảng là gọi phương thức get_data().

Vì bạn muốn các lát 2D thay vì 3D, bạn sẽ dùng các danh sách train và test đã khởi tạo trước; mỗi lần đọc một thể tích, bạn sẽ lặp qua toàn bộ 207 lát của thể tích 3D và thêm từng lát một vào danh sách.

folder = os.listdir(inp)

Hãy tải dữ liệu 3T trước.

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load(inp + f)
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        train_matrix.append(a[:,j,:])

Sau đó, tải dữ liệu 7T, và bạn dùng cùng biến folder để tải dữ liệu 7T vì số lượng thể tích MR 3T và 7T bằng nhau.

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    b = nib.load(out + f)
    b = b.get_data()
    temp[3:,:,3:] = b
    b = temp
    for j in range(full_z):
        test_matrix.append(b[:,j,:])

Tiền xử lý dữ liệu

Vì train và test matrices là danh sách, bạn sẽ dùng numpy để chuyển danh sách thành mảng numpy.

Tiếp đó, bạn chuyển type của mảng numpy thành float32 và chuẩn tỉ lệ cả input và ground truth.

train_matrix = np.asarray(train_matrix)
train_matrix = train_matrix.astype('float32')
m = min_max[0]
mi = min_max[1]
train_matrix = (train_matrix - mi) / (m - mi)

test_matrix = np.asarray(test_matrix)
test_matrix = test_matrix.astype('float32')
test_matrix = (test_matrix - mi) / (m - mi)

Hãy in nhanh kích thước của train_matrix (3T) và test_matrix (7T). Chúng sẽ có tổng 28 x 207 = 5796 ảnh, mỗi ảnh kích thước 176 x 176.

train_matrix.shape

(5796, 176, 176)

test_matrix.shape

(5796, 176, 176)

Tiếp theo, bạn sẽ tạo hai biến mới augmented_images(3T/input) và Haugmented_images (7T/ground truth) có cùng hình dạng với train và test matrix. Đây sẽ là ma trận 4 chiều trong đó chiều thứ nhất là tổng số ảnh, chiều hai và ba là kích thước mỗi ảnh và chiều cuối là số kênh, trong trường hợp này là một.

augmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])
Haugmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])

Sau đó bạn sẽ lặp qua tất cả ảnh, mỗi lần reshape ma trận train và test về 176 x 176 và thêm vào augmented_images (3T/input) và Haugmented_images (7T/ground truth) tương ứng.

for i in range(train_matrix.shape[0]):
    augmented_images[i,:,:,0] = train_matrix[i,:,:].reshape(resizeTo,resizeTo)
    Haugmented_images[i,:,:,0] = test_matrix[i,:,:].reshape(resizeTo,resizeTo)

Sau tất cả, việc phân chia dữ liệu là quan trọng. Để mô hình tổng quát tốt, bạn chia dữ liệu thành hai phần: tập huấn luyện và tập validation. Bạn sẽ huấn luyện mô hình trên 80% dữ liệu và xác thực trên 20% còn lại.

Điều này cũng giúp giảm khả năng overfitting, vì bạn sẽ xác thực mô hình trên dữ liệu mà mô hình không thấy trong giai đoạn huấn luyện.

Bạn có thể dùng module train_test_split của scikit-learn đã định nghĩa lúc đầu để chia dữ liệu đúng cách:

data,Label = shuffle(augmented_images,Haugmented_images, random_state=2)
X_train, X_test, y_train, y_test = train_test_split(data, Label, test_size=0.2, random_state=2)

X_test = np.array(X_test)
y_test = np.array(y_test)
X_test = X_test.astype('float32')
y_test = y_test.astype('float32')

X_train = np.array(X_train)
y_train = np.array(y_train)
X_train = X_train.astype('float32')
y_train = y_train.astype('float32')

Mô hình: 1-Encoder-3-Decoders!

model
Hình: Lan truyền ngược Auto-Encoder chọn lọc
Ảnh lấy từ bài báo này.

Kết nối Merge

Tiếp theo, bạn định nghĩa kiến trúc đề xuất với các khối lớp lọc liên tiếp theo sau là một lớp max pooling ở phần encoder như thể hiện trong cell tiếp theo. Để tái tạo kích thước gốc của ảnh ở đầu ra, một lớp upsampling được thêm vào mỗi khối của các decoder. Trong quá trình upsampling, có thể xuất hiện một số artifact do thiếu chi tiết trong đầu vào đã bị giảm mẫu của decoder. Do vậy, bạn nối (concatenate) đầu vào của decoder với phiên bản đã được phóng to từ encoder nhằm cung cấp bản chất chi tiết phóng to cho việc tái tạo tốt hơn khi upsampling ở decoder. Thực tế, thêm các kết nối merge sẽ đem lại cải thiện PSNR đáng kể (khoảng 5 dB). Thiết lập kiến trúc này lấy cảm hứng từ bài báo này.

Nhiều Decoder

Cách tiếp cận đề xuất sử dụng kiến trúc một encoder và nhiều decoder với đầu vào một kênh. Ba lớp tích chập được dùng trong mỗi khối của encoder và cả ba decoder, theo sau là một lớp batch normalization để giữ ổn định số học.

Khối tích chập đầu tiên trong encoder có 32 bộ lọc, và số bộ lọc tăng gấp đôi sau mỗi khối tích chập. Trong tất cả decoder, khối tích chập đầu tiên có 256 bộ lọc, và số bộ lọc giảm một nửa sau mỗi khối. Bạn sẽ dùng kích thước bộ lọc 3x3 trong tất cả các khối tích chập.

Hàm kích hoạt ReLU sẽ được dùng ở tất cả các lớp ngoại trừ lớp cuối. Vì dữ liệu được chuẩn hóa trong khoảng 0 đến 1, hàm kích hoạt Sigmoid được dùng ở lớp cuối.

Đã biết rằng chi tiết ảnh cục bộ ở nhiều thang đo đóng vai trò quan trọng trong tái tạo ảnh. Kiến trúc đề xuất xét ảnh ở các thang đo khác nhau bằng cách dùng các lớp phân cấp cho giảm mẫu (maxpooling) và tăng mẫu (mỗi lần hệ số 2) ở encoder và decoder tương ứng. Biểu diễn được mã hóa sau ba lần giảm mẫu đưa dữ liệu từ đầu vào chiều cao về biểu diễn không gian tiềm ẩn.

Encoder

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1

Decoder 1

def decoder_1(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_1 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_1

Decoder 2

def decoder_2(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_2 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_2

Decoder 3

def decoder_3(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_3 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_3

Trong 4 cell trên, bạn đã định nghĩa bốn hàm. Một cho encoder và ba hàm còn lại cho decoder. Vì là hàm, bạn có thể định nghĩa decoder() rồi gọi ba lần, nhưng để dễ hiểu hơn, bạn sẽ định nghĩa ba lần và tránh vấn đề ngẫu nhiên giữa ba decoder.

Hàm Loss

Tiếp theo, bạn sẽ dùng mean square error trong đó bạn loại trừ các giá trị (pixel) của cả ground truth y_t và y_p bằng không.

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    where = tf.not_equal(y_t, 0)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))

Định nghĩa và biên dịch mô hình

Trước tiên, bạn sẽ gọi hàm encoder bằng cách truyền input vào. Vì bạn dùng kết nối merge trong kiến trúc, hàm encoder sẽ trả về đầu ra của năm lớp convolution mà bạn sẽ merge với đầu ra của cả ba decoder riêng rẽ.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)

autoencoder_1 = Model(input_img, decoder_1(conv5,conv4,conv3,conv2,conv1))
autoencoder_1.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_2 = Model(input_img, decoder_2(conv5,conv4,conv3,conv2,conv1))
autoencoder_2.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_3 = Model(input_img, decoder_3(conv5,conv4,conv3,conv2,conv1))
autoencoder_3.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:2: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.

/usr/local/lib/python3.5/dist-packages/keras/legacy/layers.py:460: UserWarning: The `Merge` layer is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  name=name)
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:10: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  # Remove the CWD from sys.path while we load stuff.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:18: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:26: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.


WARNING:tensorflow:From /usr/local/lib/python3.5/dist-packages/keras/backend/tensorflow_backend.py:1257: calling reduce_mean (from tensorflow.python.ops.math_ops) with keep_dims is deprecated and will be removed in a future version.
Instructions for updating:
keep_dims is deprecated, use keepdims instead

Huấn luyện mô hình

Bạn sẽ chỉ lưu trọng số khi peak signal-to-noise ratio trên dữ liệu validation được cải thiện. Vì vậy, bạn sẽ định nghĩa một danh sách psnr_gray_channel trong đó thêm giá trị mặc định là 1; bạn thêm 1 vì cần khởi tạo danh sách bằng một số ngẫu nhiên và PSNR sẽ cao hơn nhiều so với giá trị này ngay từ giai đoạn đầu huấn luyện, nên đây chỉ là con số giả lập.

Bạn khởi tạo learning rate là 1e-3 và dùng chiến lược suy giảm learning rate, giảm 10% so với giá trị hiện tại sau mỗi 20 epoch!

psnr_gray_channel = []
psnr_gray_channel.append(1)
learning_rate = 0.001
j=0

Lưu ý: Đoạn mã đầy đủ tiếp theo cần chạy trong một cell hoàn chỉnh, nhưng để dễ hiểu toàn bộ quy trình huấn luyện, nó sẽ được chia thành các cell nhỏ.

Mô hình được huấn luyện trong 500 epoch. Learning rate ban đầu là 1e-3 như đã định nghĩa ở cell trước. Bạn cũng sẽ lưu các giá trị PSNR, MSE sau mỗi epoch giữa 7T và dự đoán. Bạn sẽ dùng hàm K.set_value để thay đổi learning rate của cả ba mô hình giảm 10% sau mỗi 20 epoch.

for jj in range(500):
    myfile_valid_psnr_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_psnr7T_1encoder_3decoders.txt', 'a')
    myfile_valid_mse_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_mse7T_1encoder_3decoders.txt', 'a')

    K.set_value(autoencoder_1.optimizer.lr, learning_rate)
    K.set_value(autoencoder_2.optimizer.lr, learning_rate)
    K.set_value(autoencoder_3.optimizer.lr, learning_rate)

Tiếp theo, bạn xáo trộn ảnh đầu vào 3T và ảnh ground truth 7T để tránh overfitting vì nếu không xáo trộn, mô hình sẽ thấy mẫu theo cùng thứ tự sau mỗi epoch. Sau đó, bạn tính số batch sẽ được tạo dựa trên batch_size đã định nghĩa trước. Cuối cùng, bạn bắt đầu lặp qua num_batches.

train_X,train_Y = shuffle(X_train,y_train)
print ("Epoch is: %d\n" % j)
print ("Number of batches: %d\n" % int(len(train_X)/batch_size))
num_batches = int(len(train_X)/batch_size)
for batch in range(num_batches):

Ngoài việc lưu giá trị PSNR bạn cũng sẽ lưu loss của cả ba autoencoder và ba decoder tương ứng.

myfile_ae1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae1_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae2_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae3_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec1_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec2_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec3_train_loss_1encoder_3decoders.txt', 'a')

Vì ở mỗi batch bạn muốn mô hình thấy 32 mẫu tiếp theo (batch_size), cell tiếp theo sẽ làm điều đó cho bạn!

batch_train_X = train_X[batch*batch_size:min((batch+1)*batch_size,len(train_X)),:]
batch_train_Y = train_Y[batch*batch_size:min((batch+1)*batch_size,len(train_Y)),:]

Tiếp theo, để chiến lược Minimum autoencoder hoạt động, sau mỗi epoch bạn kiểm thử trên dữ liệu huấn luyện bằng hàm test_on_batch của Keras để nhận ba loss tương ứng, và cuối cùng in chúng ra.

loss_1 = autoencoder_1.test_on_batch(batch_train_X,batch_train_Y)
loss_2 = autoencoder_2.test_on_batch(batch_train_X,batch_train_Y)
loss_3 = autoencoder_3.test_on_batch(batch_train_X,batch_train_Y)
print ('epoch_num: %d batch_num: %d Test_loss_1: %f\n' % (j,batch,loss_1))
print ('epoch_num: %d batch_num: %d Test_loss_2: %f\n' % (j,batch,loss_2))
print ('epoch_num: %d batch_num: %d Test_loss_3: %f\n' % (j,batch,loss_3))

Có sáu điều kiện có thể xảy ra trong mạng của bạn:
loss_1: Loss của Autoencoder 1
loss_2: Loss của Autoencoder 2
loss_3: Loss của Autoencoder 3

    • loss_1 có thể lớn hơn loss_2 và loss_3. Nếu đúng, bạn chỉ huấn luyện Autoencoder 1. Bạn đặt phần Encoder của Autoencoder 2 và Autoencoder 3 thành False và chỉ huấn luyện các decoder của chúng, và cuối cùng ghi tất cả loss vào tệp văn bản.
    • loss_2 có thể lớn hơn loss_1 và loss_3. Nếu đúng, bạn chỉ huấn luyện Autoencoder 2. Bạn đặt phần Encoder của Autoencoder 1 và Autoencoder 3 thành False và chỉ huấn luyện các decoder của chúng, và cuối cùng ghi tất cả loss vào tệp văn bản.
    • loss_3 có thể lớn hơn loss_1 và loss_2. Nếu đúng, bạn chỉ huấn luyện Autoencoder 3. Bạn đặt phần Encoder của Autoencoder 1 và Autoencoder 2 thành False và chỉ huấn luyện các decoder của chúng, và cuối cùng ghi tất cả loss vào tệp văn bản.
    • loss_1 có thể bằng loss_2. Nếu đúng, bạn huấn luyện một trong Autoencoder 1 hoặc 2. Bạn đặt phần Encoder của Autoencoder 3 thành False cùng với autoencoder không được chọn giữa 1 và 2, và chỉ huấn luyện các decoder của chúng, sau cùng ghi tất cả loss vào tệp văn bản.
    • loss_2 có thể bằng loss_3. Nếu đúng, bạn huấn luyện một trong Autoencoder 2 hoặc 3. Bạn đặt phần Encoder của Autoencoder 1 thành False cùng với autoencoder không được chọn giữa 2 và 3, và chỉ huấn luyện các decoder của chúng, sau cùng ghi tất cả loss vào tệp văn bản.
    • loss_3 có thể bằng loss_1. Nếu đúng, bạn huấn luyện một trong Autoencoder 3 hoặc 1. Bạn đặt phần Encoder của Autoencoder 2 thành False cùng với autoencoder không được chọn giữa 3 và 1, và chỉ huấn luyện các decoder của chúng, sau cùng ghi tất cả loss vào tệp văn bản.
model
Hình: Kiến trúc mô hình
Ảnh lấy từ bài báo này.
if loss_1 < loss_2 and loss_1 < loss_3:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
elif loss_2 < loss_1 and loss_2 < loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
elif loss_3 < loss_1 and loss_3 < loss_2:
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_2.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" %(train_1))
    myfile_dec2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
elif loss_1 == loss_2:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_2 == loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False

    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_3 == loss_1:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))


    myfile_ae1_loss.close()
    myfile_ae2_loss.close()
    myfile_ae3_loss.close()
    myfile_dec1_loss.close()
    myfile_dec2_loss.close()
    myfile_dec3_loss.close()

Đây là bước thiết yếu vì bạn đã đặt một số layer thành False trong các điều kiện trên; bạn cần đặt lại các layer đó thành True để tất cả layer của autoencoder có thể được dùng cho hàm test_on_batch và bạn không muốn chúng bị giữ False trong suốt quá trình huấn luyện.

for layer in autoencoder_1.layers[:34]:
            layer.trainable = True
for layer in autoencoder_2.layers[:34]:
            layer.trainable = True
for layer in autoencoder_3.layers[:34]:
            layer.trainable = True

Bạn sẽ lưu trọng số theo ba cách, dưới đây là hai trong số đó. Một là lưu trọng số sau mỗi 100 epoch, và hai là lưu trọng số sau mỗi epoch. Về cơ bản, chúng sẽ bị ghi đè sau mỗi epoch.

if jj % 100 ==0:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")


autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")
autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")
autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")

Kiểm thử trên dữ liệu validation

Trước tiên bạn xáo trộn dữ liệu validation rồi dùng cùng 6 điều kiện có thể xảy ra như đã định nghĩa ở trên. Điều kiện nào đúng thì bạn sẽ test mô hình bằng autoencoder tương ứng.

Sau đó bạn tính hai chỉ số MSE và PSNR giữa decoded_imgs (dự đoán) và ground truth. Cuối cùng, bạn lưu chúng vào các tệp văn bản.

X_test,y_test = shuffle(X_test,y_test)
if loss_1 < loss_2 and loss_1 < loss_3:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)


    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_2 < loss_1 and loss_2 < loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_3 < loss_2 and loss_3 < loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))
    #print (check)

elif loss_1 == loss_2:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_2 == loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_3 == loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

Ở đây bạn chỉ lưu trọng số khi PSNR giữa 7T dự đoán và ground truth (7T) là lớn nhất so với các giá trị PSNR trước đó đã lưu trong danh sách psnr_gray_channel.

if max(psnr_gray_channel) < psnr_7T:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")

    psnr_gray_channel.append(psnr_7T)

Lưu Input, Ground Truth và ảnh đã giải mã: Kết quả định lượng

Bạn sẽ định nghĩa một ma trận numpy temp kích thước 176 x 528 vì sẽ có 3 ảnh lưu trên một hàng, mỗi ảnh kích thước 176 x 176. Bạn sẽ lưu một trong các ảnh từ toàn bộ tập validation 3T, 7T và dự đoán, đồng thời nhân ma trận với 255 vì ảnh đã được chuẩn tỉ lệ về khoảng 0 đến 1.

Cuối cùng dùng hàm của scipy, bạn sẽ lưu ảnh.

temp = np.zeros([resizeTo,resizeTo*3])
temp[:resizeTo,:resizeTo] = X_test[0,:,:,0]
temp[:resizeTo,resizeTo:resizeTo*2] = y_test[0,:,:,0]
temp[:resizeTo,2*resizeTo:] = decoded_imgs[0,:,:,0]
temp = temp*255
scipy.misc.imsave('../Results/1_encoder_3_decoders_complete_slices_single_channel/' + str(j) + '.jpg', temp)
j +=1

Hãy đóng các tệp PSNR và MSE ở cuối.

myfile_valid_psnr_7T.close()
myfile_valid_mse_7T.close()

Cuối cùng, ở cuối, bạn sẽ giảm learning rate đi 10% so với giá trị hiện tại sau mỗi 20 epoch.

if jj % 20 ==0:
        learning_rate = learning_rate - learning_rate * 0.10

Testing Script

Trong phần 2 của hướng dẫn này, bạn sẽ đề cập các chủ đề sau:

  • Bắt đầu bằng cách import các module cần thiết để huấn luyện mô hình deep learning,
  • Sau đó giới thiệu ngắn về bộ dữ liệu MRI 3T và 7T,
  • Tiếp theo bạn sẽ định nghĩa các initializer và tải bộ dữ liệu kiểm thử 3T và 7T, đồng thời trong quá trình tải sẽ thay đổi kích thước ảnh ngay lập tức,
  • Kế đến bạn sẽ tiền xử lý dữ liệu đã tải: chuyển danh sách train và test thành ma trận numpy, chuyển kiểu dữ liệu ma trận thành float32, chuẩn tỉ lệ ma trận theo chiến lược max-min, reshape mảng và cuối cùng chia dữ liệu thành 80% huấn luyện và 20% còn lại làm tập validation,
  • Sau đó, bạn sẽ tạo Kiến trúc 1-Encoder-3-Decoder: gồm các kết nối merge và nhiều decoder,
  • Tiếp theo, bạn sẽ định nghĩa: hàm loss, ba mô hình khác nhau và cuối cùng tải các trọng số đã huấn luyện,
  • Cuối cùng, đến lúc bạn sẽ dự đoán bằng mô hình merge và multi-decoder trên dữ liệu chưa thấy và lưu cả kết quả định lượng lẫn định tính. Bạn cũng sẽ học cách lưu các ảnh 2D thành một khối thể tích kết hợp bằng nibabel.

Các phụ thuộc module Python

Trước khi làm theo hướng dẫn, hãy đảm bảo bạn có đúng phiên bản các module như dưới đây:

  • Keras==2.0.4
  • tensorflow==1.8.0
  • scipy==0.19.0
  • numpy==1.14.5
  • Pillow==4.1.1
  • nibabel==2.1.0
  • scikit_learn==0.18.1

Lưu ý: Trước khi bắt đầu, lưu ý rằng mô hình sẽ được huấn luyện trên hệ thống với GPU Nvidia 1080 Ti, CPU Xeon e5 GeForce và RAM 32GB. Nếu bạn dùng Jupyter Notebook, bạn sẽ cần thêm ba dòng mã để chỉ định thứ tự thiết bị CUDA và các thiết bị CUDA khả dụng bằng module os.

Trong đoạn mã dưới, về cơ bản bạn đặt biến môi trường trong notebook bằng os.environ. Nên làm điều này trước khi khởi tạo Keras để giới hạn Keras backend TensorFlow sử dụng GPU đầu tiên. Nếu máy bạn dùng để huấn luyện có GPU ở chỉ số 0, hãy dùng 0 thay vì 1. Bạn có thể kiểm tra bằng lệnh đơn giản trên terminal: ví dụ, nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU

Import các module

Trước tiên, bạn import tất cả module cần thiết như tensorflow, numpy và quan trọng nhất là keras cùng các hàm hoặc layer như Input, Conv2D,MaxPooling2D,... vì bạn sẽ dùng chúng để huấn luyện mô hình!
Để đọc ảnh định dạng nifti, bạn cũng cần import module nibabel.

import os
from keras.layers import Input,Dense,Flatten,Dropout,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adadelta, RMSprop,SGD,Adam
from keras import regularizers
from keras import backend as K
import numpy as np
import scipy.misc
import numpy.random as rng
from sklearn.utils import shuffle
import nibabel as nib
from sklearn.cross_validation import train_test_split
import math

Tìm hiểu bộ dữ liệu MRI não 3T và 7T

Bộ dữ liệu MRI não 3T và 7T gồm các thể tích 3D, mỗi thể tích có tổng cộng 207 lát/cảnh MRI não được chụp tại các lát cắt khác nhau của não. Mỗi lát có kích thước 173 x 173. Các ảnh là ảnh thang độ xám một kênh. Tổng cộng có 39 đối tượng, mỗi đối tượng chứa ảnh quét MRI của một bệnh nhân. Định dạng ảnh không phải jpeg, png, v.v. mà là định dạng nifti. Ở phần sau, bạn sẽ thấy cách đọc ảnh định dạng nifti.

Bộ dữ liệu gồm ảnh MR theo chuỗi T1; các chuỗi T1 truyền thống được xem là tốt cho việc đánh giá cấu trúc giải phẫu. Bộ dữ liệu mà bạn sẽ làm việc hôm nay gồm MRI não 3T và 7T.

Bộ dữ liệu là công khai và có thể tải về tại nguồn này.

Định nghĩa các bộ khởi tạo

Trước hết, hãy xác định kích thước của dữ liệu. Bạn sẽ đổi kích thước ảnh từ 173x173 lên 176x176 trong phần đọc dữ liệu. Ở đây, chúng ta cũng xác định thư mục dữ liệu, kích thước lô dùng để huấn luyện mô hình, số kênh, một tầng Input(), các ma trận train và test dưới dạng danh sách, và cuối cùng để tái tỷ lệ, chúng ta nạp tệp văn bản chứa giá trị nhỏ nhất và lớn nhất của bộ dữ liệu MRI.

x,y = 173,173
full_z = 207
resizeTo=176
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))

train_matrix = []
test_matrix = []
ff = os.listdir("../test_crossval1")
save = "../Result_nii_crossval1"
folder_ground = os.listdir("../test_g_crossval1")
ToPredict_images=[]
predict_matrix=[]

ground_images=[]
ground_matrix=[]
min_max=np.loadtxt('../maxANDmin.txt')

Tải các thể tích kiểm thử

Tiếp theo, chúng ta tải dữ liệu mri bằng thư viện nibabel và đổi kích thước ảnh từ 173 x 173 lên 176 x 176 bằng cách đệm các số 0 theo chiều x và y.

Lưu ý rằng khi bạn tải một thể tích định dạng Nifti, Nibabel không nạp mảng ảnh ngay. Nó chờ cho đến khi bạn yêu cầu dữ liệu mảng. Cách tiêu chuẩn để yêu cầu dữ liệu mảng là gọi phương thức get_data().

Vì bạn muốn các lát 2D thay vì 3D, bạn sẽ dùng các danh sách train và test đã khởi tạo trước; mỗi lần đọc một thể tích, bạn sẽ lặp qua toàn bộ 207 lát của thể tích 3D và thêm từng lát một vào danh sách.

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        predict_matrix.append(a[:,j,:])

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_g_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        ground_matrix.append(a[:,j,:])

Tiền xử lý dữ liệu

Vì ma trận kiểm thử 3T và 7T là danh sách, bạn sẽ dùng module numpy để chuyển danh sách thành một mảng numpy.

Tiếp nữa, bạn sẽ chuyển type của mảng numpy thành float32 và cũng tái tỷ lệ cả input và ground truth.

ToPredict_images = np.asarray(predict_matrix)
ToPredict_images = ToPredict_images.astype('float32')
mx = min_max[0]
mn = min_max[1]
ToPredict_images[:,:,:,0] = (ToPredict_images[:,:,:,0] - mn ) / (mx - mn)

ground_images = np.asarray(ground_matrix)
ground_images = ground_images.astype('float32')
ground_images[:,:,:,0] = (ground_images[:,:,:,0] - mn ) / (mx - mn)

Tiếp theo, bạn sẽ tạo hai biến mới ToPredict_images (3T Test/input) và ground_images (7T Test/ground truth) có cùng dạng với ma trận train và test. Đây sẽ là ma trận 4D, trong đó chiều thứ nhất là tổng số ảnh, chiều thứ hai và thứ ba là kích thước của mỗi ảnh, và chiều cuối là số kênh, trong trường hợp này là một.

ToPredict_images=np.zeros(shape=[(ToPredict_images.shape[0]),(ToPredict_images.shape[1]),(ToPredict_images.shape[2]),(1)])
ground_images=np.zeros(shape=[(ground_images.shape[0]),(ground_images.shape[1]),(ground_images.shape[2]),(1)])

Sau đó, bạn sẽ lặp qua toàn bộ ảnh từng cái một, mỗi lần sẽ reshape ma trận train và test về 176 x 176 và thêm vào ToPredict_images (3T Test/input) và ground_images (7T Test/ground truth) tương ứng.

for i in range(ToPredict_images.shape[0]):
    ToPredict_images[i,:,:,0] = ToPredict_images[i,:,:].reshape(resizeTo,resizeTo)

for i in range(ground_images.shape[0]):
    ground_images[i,:,:,0] = ground_images[i,:,:].reshape(resizeTo,resizeTo)

Mô hình!

model
Hình: Kiến trúc mô hình
Hình lấy từ bài báo này.

Encoder

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1

Decoder

def decoder(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded

Hàm mất mát

Tiếp theo, bạn sẽ dùng sai số bình phương trung bình, trong đó bạn loại trừ các giá trị (pixel) của cả y_t ground truth và y_p bằng không.

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    zero = tf.constant(0, dtype=tf.float32)
    where = tf.not_equal(y_t, zero)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))

Định nghĩa mô hình và nạp trọng số cho cả ba Autoencoder!

Trước tiên, hãy gọi hàm encoder bằng cách truyền input vào đó. Vì bạn sử dụng kết nối merge trong kiến trúc, hàm encoder sẽ trả về đầu ra của năm tầng convolution mà sau đó bạn sẽ merge với đầu ra của cả ba decoder riêng biệt.

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)

Bây giờ, hãy tạo ba mô hình khác nhau và nạp ba bộ trọng số đã huấn luyện vào chúng.

autoencoder_1 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_1.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")

autoencoder_2 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_2.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")

autoencoder_3 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_3.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")

Dự đoán mô hình trên các thể tích kiểm thử: Kết quả định lượng và định tính!

Hãy nhanh chóng khởi tạo hai mảng numpy, mỗi mảng kích thước 11 x 3 x 1, trong đó chiều thứ nhất biểu diễn số thể tích bạn sẽ dùng để kiểm thử mô hình. Chiều thứ hai biểu diễn MSE và PSNR giữa: Ảnh MR 7T đầu ra dự đoán và ảnh 7T ground truth, đầu ra dự đoán và ảnh MR 3T đầu vào, ảnh MR 3T đầu vào và ảnh MR 7T ground truth; chiều thứ ba biểu diễn số kênh bạn sẽ đưa vào mô hình!

mse= np.zeros([11,3,1])
psnr= np.zeros([11,3,1])
i=0 #for iterating over the slices of all the 11 volumes

Ở phần tiếp theo của ô, bạn sẽ lần lượt lặp qua cả 11 thể tích. Mỗi vòng lặp, bạn sẽ dự đoán từng thể tích bằng cả ba autoencoder và cuối cùng lấy trung bình các dự đoán.

Trong mỗi thể tích, bạn sẽ lặp qua số kênh của thể tích đó và tính PSNR và MSE cho ba trường hợp như đã nêu ở trên.

Sau đó, dùng thư viện nibabel, bạn sẽ lưu đầu ra dự đoán, đầu vào (3T) và ground truth (7T) dưới dạng tệp .nii: mỗi trong số 11 thể tích gồm 207 lát.

Cuối cùng, bạn sẽ lưu ma trận PSNR vào một tệp văn bản bằng thư viện numpy.

Như đã nêu trong bài báo, việc lấy trung bình các đầu ra dự đoán giúp giảm nhiễu nhưng vẫn giữ được các đặc trưng cục bộ trong ảnh tái tạo, nhờ đó PSNR được cải thiện so với đầu ra từng decoder riêng lẻ.

for j in range(11):
    decoded_imgs_1 = autoencoder_1.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_2 = autoencoder_2.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_3 = autoencoder_3.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs = np.mean( np.array([ decoded_imgs_1, decoded_imgs_2,decoded_imgs_3 ]), axis=0 )
    for channel in range(1):
        mse[j,0,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        psnr[j,0,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,0,channel]))
        mse[j,1,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - ToPredict_images[i:i+207,:,:,channel])** 2)
        psnr[j,1,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,1,channel]))
        mse[j,2,channel] =  np.mean((ToPredict_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        checklt = math.sqrt(mse[j,2,channel])
        psnr[j,2,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,2,channel]))
    obj = nib.Nifti1Image(decoded_imgs, affine)
    string =str(j)+'_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ground_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ground_images_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ToPredict_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ToPredict_images_crossval1.nii'
    nib.save(obj, save + string)
    i=i+207


np.savetxt('psnr_all_slices.txt',psnr[:,:,0])

Nếu bạn muốn học thêm về Python, hãy tham gia khóa Giới thiệu Trực quan hóa Dữ liệu với Matplotlib của DataCamp.

Hãy xem Hướng dẫn Keras: Deep Learning với Python.

Chủ đề

PythonMachine LearningHọc Sâu

Liên quan

Xem thêm →