使用腾讯云 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
<em style="font-style:italic"># Defining the voxel range in which the nodule falls</em>
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)

<em style="font-style:italic"># Convert back to world coordinates for distance calculation</em>
x_data = [x*spacing[0]+origin[0] for x in range(width)]
y_data = [x*spacing[1]+origin[1] for x in range(height)]

<em style="font-style:italic"># Fill in 1 within sphere around nodule</em>
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]))&lt;=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')

<em style="font-style:italic"># %% Define the two layer localisation network</em>
h_fc_loc1 = tf.nn.tanh(tf.matmul(x, W_fc_loc1) + b_fc_loc1)
<em style="font-style:italic"># %% We can add dropout for regularizing and to reduce overfitting like so:</em>

h_fc_loc1_drop = tf.nn.dropout(h_fc_loc1, keep_prob)
<em style="font-style:italic"># %% Second layer</em>
h_fc_loc2 = tf.nn.tanh(tf.matmul(h_fc_loc1_drop, W_fc_loc2) + b_fc_loc2)

<em style="font-style:italic"># %% We&#39;ll create a spatial transformer module to identify discriminative</em>
<em style="font-style:italic"># %% patches</em>
out_size = (40, 40)
h_trans = transformer(x_tensor, h_fc_loc2, out_size)

<em style="font-style:italic"># %% We&#39;ll setup the first convolutional layer</em>
<em style="font-style:italic"># Weight matrix is [height x width x input_channels x output_channels]</em>
filter_size = 3
n_filters_1 = 16
W_conv1 = weight_variable([filter_size, filter_size, 1, n_filters_1])

<em style="font-style:italic"># %% Bias is [output_channels]</em>
b_conv1 = bias_variable([n_filters_1])

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

<em style="font-style:italic"># %% And just like the first layer, add additional layers to create</em>
<em style="font-style:italic"># a deep net</em>
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=&#39;SAME&#39;) +
    b_conv2)

<em style="font-style:italic"># %% We&#39;ll now reshape so we can connect to a fully-connected layer:</em>
h_conv2_flat = tf.reshape(h_conv2, [-1, 10 * 10 * n_filters_2])

<em style="font-style:italic"># %% Create a fully-connected layer:</em>
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)

<em style="font-style:italic"># %% And finally our softmax layer:</em>
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)