使用腾讯云 GPU 学习深度学习系列之四：深度学习的特征工程

from __future__ import print_function, division

import numpy as np
import os
import csv
from glob import glob
import pandas as pd
import numpy as np
import SimpleITK as sitk

from skimage import measure,morphology
from sklearn.cluster import KMeans
from skimage.transform import resize
import matplotlib.pyplot as plt
import seaborn as sns

from glob import glob

try:
    from tqdm import tqdm # long waits are not fun
except:
    print('TQDM does make much nicer wait bars...')
    tqdm = lambda x: x

def make_mask(center,diam,z,width,height,spacing,origin):
    '''
    Center : centers of circles px -- list of coordinates x,y,z
    diam : diameters of circles px -- diameter
    widthXheight : pixel dim of image
    spacing = mm/px conversion rate np array x,y,z
    origin = x,y,z mm np.array
    z = z position of slice in world coordinates mm
    '''
    mask = np.zeros([height,width]) # 0's everywhere except nodule swapping x,y to match img
    #convert to nodule space from world coordinates

    # Defining the voxel range in which the nodule falls
    v_center = (center-origin)/spacing
    v_diam = int(diam/spacing[0]+5)
    v_xmin = np.max([0,int(v_center[0]-v_diam)-5])
    v_xmax = np.min([width-1,int(v_center[0]+v_diam)+5])
    v_ymin = np.max([0,int(v_center[1]-v_diam)-5])
    v_ymax = np.min([height-1,int(v_center[1]+v_diam)+5])

    v_xrange = range(v_xmin,v_xmax+1)
    v_yrange = range(v_ymin,v_ymax+1)

    # Convert back to world coordinates for distance calculation
    x_data = [x*spacing[0]+origin[0] for x in range(width)]
    y_data = [x*spacing[1]+origin[1] for x in range(height)]

    # Fill in 1 within sphere around nodule
    for v_x in v_xrange:
        for v_y in v_yrange:
            p_x = spacing[0]*v_x + origin[0]
            p_y = spacing[1]*v_y + origin[1]
            if np.linalg.norm(center-np.array([p_x,p_y,z]))<=diam:
                mask[int((p_y-origin[1])/spacing[1]),int((p_x-origin[0])/spacing[0])] = 1.0
    return(mask)

def matrix2int16(matrix):
    '''
    matrix must be a numpy array NXN
    Returns uint16 version
    '''
    m_min= np.min(matrix)
    m_max= np.max(matrix)
    matrix = matrix-m_min
    return(np.array(np.rint( (matrix-m_min)/float(m_max-m_min) * 65535.0),dtype=np.uint16))

df_node = pd.read_csv('./annotation.csv')
for fcount, img_file in enumerate(tqdm(df_node['seriesuid'])):
    mini_df = df_node[df_node["seriesuid"]==img_file] #get all nodules associate with file
    if mini_df.shape[0]>0: # some files may not have a nodule--skipping those
        # load the data once
        itk_img = sitk.ReadImage("%s.mhd" % img_file)
        img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
        num_z, height, width = img_array.shape        #heightXwidth constitute the transverse plane
        origin = np.array(itk_img.GetOrigin())      # x,y,z  Origin in world coordinates (mm)
        spacing = np.array(itk_img.GetSpacing())    # spacing of voxels in world coor. (mm)
        # go through all nodes (why just the biggest?)
        for node_idx, cur_row in mini_df.iterrows():       
            node_x = cur_row["coordX"]
            node_y = cur_row["coordY"]
            node_z = cur_row["coordZ"]
            diam = cur_row["diameter_mm"]
            # just keep 3 slices
            imgs = np.ndarray([3,height,width],dtype=np.float32)
            masks = np.ndarray([3,height,width],dtype=np.uint8)
            center = np.array([node_x, node_y, node_z])   # nodule center
            v_center = np.rint((center-origin)/spacing)  # nodule center in voxel space (still x,y,z ordering)
            for i, i_z in enumerate(np.arange(int(v_center[2])-1,
                             int(v_center[2])+2).clip(0, num_z-1)): # clip prevents going out of bounds in Z
                mask = make_mask(center, diam, i_z*spacing[2]+origin[2],
                                     width, height, spacing, origin)
                masks[i] = mask
                imgs[i] = img_array[i_z]
            np.save(os.path.join("./images_%04d_%04d.npy" % (fcount,node_idx)),imgs)
            np.save(os.path.join("./masks_%04d_%04d.npy"  % (fcount,node_idx)),masks)

img_file = "./images_0000_0000.npy"
imgs_to_process = np.load(img_file).astype(np.float64)

fig = plt.figure(figsize=(12,4))
for i in range(3):
    ax = fig.add_subplot(1,3,i+1)
    ax.imshow(imgs_to_process[i,:,:], 'bone')
    ax.set_axis_off()

i = 0
img = imgs_to_process[i]

#Standardize the pixel values
mean = np.mean(img)
std = np.std(img)
img = img-mean
img = img/std

# Find the average pixel value near the lungs
#         to renormalize washed out images
middle = img[100:400,100:400]
mean = np.mean(middle)  
max = np.max(img)
min = np.min(img)

# To improve threshold finding, I'm moving the
#         underflow and overflow on the pixel spectrum
img[img==max]=mean
img[img==min]=mean


# Using Kmeans to separate foreground (radio-opaque tissue)
#     and background (radio transparent tissue ie lungs)
# Doing this only on the center of the image to avoid
#     the non-tissue parts of the image as much as possible
kmeans = KMeans(n_clusters=2).fit(np.reshape(middle,[np.prod(middle.shape),1]))
centers = sorted(kmeans.cluster_centers_.flatten())
threshold = np.mean(centers)
thresh_img = np.where(img<threshold,1.0,0.0)  # threshold the image

# 对一张图中所有像素点的亮度做概率密度分布, 用竖线标注阈值所在
fig = plt.figure(figsize=(8,12))
ax1 = fig.add_subplot(321)
sns.distplot(middle.ravel(), ax=ax1)
ax1.vlines(x=threshold, ymax=10, ymin=0)
ax1.set_title('Threshold: %1.2F' % threshold)
ax1.set_xticklabels([])

# 展示阈值对图像切割的结果。小于阈值的点标注为1，白色。大于阈值的点标注为0，黑色。
ax2 = fig.add_subplot(322)
ax2.imshow(thresh_img, "gray")
ax2.set_axis_off()
ax2.set_title('Step1, using threshold as cutoff')

# 增大黑色部分（非ROI）的区域，使之尽可能的连在一起
eroded = morphology.erosion(thresh_img,np.ones([4,4]))
ax3 = fig.add_subplot(323)
ax3.imshow(eroded, "gray")
ax3.set_axis_off()
ax3.set_title('Step2,  erosion shrinks bright\nregions and enlarges dark regions.')

# 增大白色部分（ROI）的区域，尽可能的消除面积较小的黑色区域
dilation = morphology.dilation(eroded,np.ones([10,10]))
ax4 = fig.add_subplot(324)
ax4.imshow(dilation, "gray")
ax4.set_axis_off()
ax4.set_title('Step3,  dilation shrinks dark\nregions and enlarges bright regions.')

# 上一张图中共有三片连续区域，即最外层的体外区域，内部的肺部区域，以及二者之间的身体轮廓区域。这里将其分别标出
labels = measure.label(dilation)
ax5 = fig.add_subplot(325)
ax5.imshow(labels, "gray")
#ax5.set_axis_off()
ax5.set_title('Step4, label connected regions\n of an integer array.')


# 提取regions 信息，这张图片的 region的 bbox位置分别在 [[0,0,512,512],[141, 86, 396, 404]],
#   分别对应 体外+轮廓 以及 肺部区域的左上角、右下角坐标。
#   于是这里通过区域的宽度 B[2]-B[0]、高度 B[3]-B[1]
#   以及距离图片上下的距离  B[0]>40 and B[2]<472,
#   最终保留需要的区域。
regions = measure.regionprops(labels)
good_labels = []
for prop in regions:
    B = prop.bbox
    if B[2]-B[0]<475 and B[3]-B[1]<475 and B[0]>40 and B[2]<472:
        good_labels.append(prop.label)

mask = np.zeros_like(labels)
for N in good_labels:
    mask = mask + np.where(labels==N,1,0)

mask = morphology.dilation(mask,np.ones([10,10])) # one last dilation

ax6 = fig.add_subplot(326)
ax6.imshow(mask, "gray")
ax6.set_axis_off()
ax6.set_title('Step5, remain the region of interests.')

fig = plt.figure(figsize=(8,4))

ax1 = fig.add_subplot(1,2,1)
ax1.imshow(imgs_to_process[0,:,:], 'bone')
ax1.set_axis_off()
ax1.set_title("Raw CT image")

ax2 = fig.add_subplot(1,2,2)
ax2.imshow(imgs_to_process[0,:,:]*mask, 'bone')
ax2.set_axis_off()
ax2.set_title("Pre-processed Images with\nregion of interest in lung")

wget https://github.com/daviddao/spatial-transformer-tensorflow/raw/master/data/mnist_sequence1_sample_5distortions5x5.npz

import tensorflow as tf
# https://github.com/tensorflow/models/tree/master/transformer
from spatial_transformer import transformer  
import numpy as np
from tf_utils import weight_variable, bias_variable, dense_to_one_hot
import matplotlib.pyplot as plt
from keras.backend.tensorflow_backend import set_session

np.random.seed(0)
tf.set_random_seed(0)

config = tf.ConfigProto()
config.gpu_options.allow_growth=True
set_session(tf.Session(config=config))
%matplotlib inline

mnist_cluttered = np.load('./mnist_sequence1_sample_5distortions5x5.npz')

X_train = mnist_cluttered['X_train']
y_train = mnist_cluttered['y_train']
X_valid = mnist_cluttered['X_valid']
y_valid = mnist_cluttered['y_valid']
X_test = mnist_cluttered['X_test']
y_test = mnist_cluttered['y_test']

Y_train = dense_to_one_hot(y_train, n_classes=10)
Y_valid = dense_to_one_hot(y_valid, n_classes=10)
Y_test = dense_to_one_hot(y_test, n_classes=10)

x = tf.placeholder(tf.float32, [None, 1600])
keep_prob = tf.placeholder(tf.float32)

iter_per_epoch = 100
n_epochs = 500
train_size = 10000
indices = np.linspace(0, 10000 - 1, iter_per_epoch)
indices = indices.astype('int')
iter_i = 0

batch_xs = X_train[indices[iter_i]:indices[iter_i+1]]
x_tensor = tf.reshape(x, [-1, 40, 40, 1])

sess = tf.Session()
sess.run(tf.global_variables_initializer())
xout = sess.run(x_tensor, feed_dict={x: batch_xs})

x = tf.placeholder(tf.float32, [None, 1600])
y = tf.placeholder(tf.float32, [None, 10])

x_tensor = tf.reshape(x, [-1, 40, 40, 1])
W_fc_loc1 = weight_variable([1600, 20])
b_fc_loc1 = bias_variable([20])

W_fc_loc2 = weight_variable([20, 6])
initial = np.array([[1., 0, 0], [0, 1., 0]])
initial = initial.astype('float32')
initial = initial.flatten()
b_fc_loc2 = tf.Variable(initial_value=initial, name='b_fc_loc2')

# %% Define the two layer localisation network
h_fc_loc1 = tf.nn.tanh(tf.matmul(x, W_fc_loc1) + b_fc_loc1)

# %% We can add dropout for regularizing and to reduce overfitting like so:
keep_prob = tf.placeholder(tf.float32)
h_fc_loc1_drop = tf.nn.dropout(h_fc_loc1, keep_prob)

# %% Second layer
h_fc_loc2 = tf.nn.tanh(tf.matmul(h_fc_loc1_drop, W_fc_loc2) + b_fc_loc2)

# %% We'll create a spatial transformer module to identify discriminative
# %% patches
out_size = (40, 40)
h_trans = transformer(x_tensor, h_fc_loc2, out_size)

iter_i = 0
batch_xs = X_train[0:101]
batch_ys = Y_train[0:101]

sess = tf.Session()
sess.run(tf.global_variables_initializer())

xtransOut = sess.run(h_trans,
                            feed_dict={
                                x: batch_xs,
                                y: batch_ys,
                                keep_prob: 1.0
})

fig = plt.figure(figsize=(10,2))

for idx in range(10):
    ax1 = fig.add_subplot(2,10,idx+1)
    ax2 = fig.add_subplot(2,10,idx+11)

    ax1.imshow(xout[idx,:,:,0], "gray")
    ax2.imshow(xtransOut[idx,:,:,0], "gray")
    ax1.set_axis_off()
    ax2.set_axis_off()
    ax1.set_title(np.argmax(batch_ys, axis=1)[idx])
    ax2.set_title(np.argmax(batch_ys, axis=1)[idx])

x = tf.placeholder(tf.float32, [None, 1600])
y = tf.placeholder(tf.float32, [None, 10])
keep_prob = tf.placeholder(tf.float32)

def Networks(x,keep_prob, SpatialTrans=True):
    x_tensor = tf.reshape(x, [-1, 40, 40, 1])
    W_fc_loc1 = weight_variable([1600, 20])
    b_fc_loc1 = bias_variable([20])

    W_fc_loc2 = weight_variable([20, 6])
    initial = np.array([[1., 0, 0], [0, 1., 0]])
    initial = initial.astype('float32')
    initial = initial.flatten()
    b_fc_loc2 = tf.Variable(initial_value=initial, name='b_fc_loc2')

    # %% Define the two layer localisation network
    h_fc_loc1 = tf.nn.tanh(tf.matmul(x, W_fc_loc1) + b_fc_loc1)
    # %% We can add dropout for regularizing and to reduce overfitting like so:

    h_fc_loc1_drop = tf.nn.dropout(h_fc_loc1, keep_prob)
    # %% Second layer
    h_fc_loc2 = tf.nn.tanh(tf.matmul(h_fc_loc1_drop, W_fc_loc2) + b_fc_loc2)

    # %% We'll create a spatial transformer module to identify discriminative
    # %% patches
    out_size = (40, 40)
    h_trans = transformer(x_tensor, h_fc_loc2, out_size)

    # %% We'll setup the first convolutional layer
    # Weight matrix is [height x width x input_channels x output_channels]
    filter_size = 3
    n_filters_1 = 16
    W_conv1 = weight_variable([filter_size, filter_size, 1, n_filters_1])

    # %% Bias is [output_channels]
    b_conv1 = bias_variable([n_filters_1])

    # %% Now we can build a graph which does the first layer of convolution:
    # we define our stride as batch x height x width x channels
    # instead of pooling, we use strides of 2 and more layers
    # with smaller filters.
    if SpatialTrans:
        h_conv1 = tf.nn.relu(
            tf.nn.conv2d(input=h_trans,
                         filter=W_conv1,
                         strides=[1, 2, 2, 1],
                         padding='SAME') +
            b_conv1)
    else:
        h_conv1 = tf.nn.relu(
            tf.nn.conv2d(input=x_tensor,
                         filter=W_conv1,
                         strides=[1, 2, 2, 1],
                         padding='SAME') +
            b_conv1)

    # %% And just like the first layer, add additional layers to create
    # a deep net
    n_filters_2 = 16
    W_conv2 = weight_variable([filter_size, filter_size, n_filters_1, n_filters_2])
    b_conv2 = bias_variable([n_filters_2])
    h_conv2 = tf.nn.relu(
        tf.nn.conv2d(input=h_conv1,
                     filter=W_conv2,
                     strides=[1, 2, 2, 1],
                     padding='SAME') +
        b_conv2)

    # %% We'll now reshape so we can connect to a fully-connected layer:
    h_conv2_flat = tf.reshape(h_conv2, [-1, 10 * 10 * n_filters_2])

    # %% Create a fully-connected layer:
    n_fc = 1024
    W_fc1 = weight_variable([10 * 10 * n_filters_2, n_fc])
    b_fc1 = bias_variable([n_fc])
    h_fc1 = tf.nn.relu(tf.matmul(h_conv2_flat, W_fc1) + b_fc1)

    h_fc1_drop = tf.nn.dropout(h_fc1, keep_prob)

    # %% And finally our softmax layer:
    W_fc2 = weight_variable([n_fc, 10])
    b_fc2 = bias_variable([10])
    y_logits = tf.matmul(h_fc1_drop, W_fc2) + b_fc2
    return y_logits


# %% We'll now train in minibatches and report accuracy, loss:
iter_per_epoch = 100
n_epochs = 100
train_size = 10000

indices = np.linspace(0, 10000 - 1, iter_per_epoch)
indices = indices.astype('int')

y_logits_F = Networks(x, keep_prob, SpatialTrans=False)
cross_entropy_F = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=y_logits_F, labels=y))
opt = tf.train.AdamOptimizer()
optimizer_F = opt.minimize(cross_entropy_F)
grads_F = opt.compute_gradients(cross_entropy_F, [b_fc_loc2])
correct_prediction_F = tf.equal(tf.argmax(y_logits_F, 1), tf.argmax(y, 1))
accuracy_F = tf.reduce_mean(tf.cast(correct_prediction_F, 'float'))

sessF = tf.Session()
sessF.run(tf.global_variables_initializer())

l_acc_F = []
for epoch_i in range(n_epochs):
    for iter_i in range(iter_per_epoch - 1):
        batch_xs = X_train[indices[iter_i]:indices[iter_i+1]]
        batch_ys = Y_train[indices[iter_i]:indices[iter_i+1]]
        sessF.run(optimizer_F, feed_dict={
            x: batch_xs, y: batch_ys, keep_prob: 0.8})

    acc = sessF.run(accuracy_F,
                                                     feed_dict={
                                                         x: X_valid,
                                                         y: Y_valid,
                                                         keep_prob: 1.0
                                                     })
    l_acc_F.append(acc)
    if epoch_i % 10 == 0:
        print('Accuracy (%d): ' % epoch_i + str(acc))

Accuracy (0): 0.151
Accuracy (10): 0.813
Accuracy (20): 0.832
Accuracy (30): 0.825
Accuracy (40): 0.833
Accuracy (50): 0.837
Accuracy (60): 0.832
Accuracy (70): 0.837
Accuracy (80): 0.833
Accuracy (90): 0.843

y_logits = Networks(x, keep_prob)
cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=y_logits, labels=y))
opt = tf.train.AdamOptimizer()
optimizer = opt.minimize(cross_entropy)
grads = opt.compute_gradients(cross_entropy, [b_fc_loc2])
correct_prediction = tf.equal(tf.argmax(y_logits, 1), tf.argmax(y, 1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, 'float'))

sess = tf.Session()
sess.run(tf.global_variables_initializer())

l_acc = []
for epoch_i in range(n_epochs):
    for iter_i in range(iter_per_epoch - 1):
        batch_xs = X_train[indices[iter_i]:indices[iter_i+1]]
        batch_ys = Y_train[indices[iter_i]:indices[iter_i+1]]
        sess.run(optimizer, feed_dict={
            x: batch_xs, y: batch_ys, keep_prob: 0.8})

    acc = sess.run(accuracy,
                                                     feed_dict={
                                                         x: X_valid,
                                                         y: Y_valid,
                                                         keep_prob: 1.0
                                                     })
    l_acc.append(acc)
    if epoch_i % 10 == 0:
        print('Accuracy (%d): ' % epoch_i + str(acc))

Accuracy (0): 0.25
Accuracy (10): 0.92
Accuracy (20): 0.94
Accuracy (30): 0.955
Accuracy (40): 0.943
Accuracy (50): 0.944
Accuracy (60): 0.957
Accuracy (70): 0.948
Accuracy (80): 0.941
Accuracy (90): 0.948

plt.plot(l_acc, label="Using Spatial Transform")
plt.plot(l_acc_F, label="Raw input")
plt.legend(loc=8)