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

• 批量输入模块
• 各种深度学习零件搭建的深度神经网络
• 凸优化模块
• 模型的训练与评估

## 1. 结合传统数据处理方法的特征提取

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

'''
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:

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))

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
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)
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
width, height, spacing, origin)
imgs[i] = img_array[i_z]
np.save(os.path.join("./images_%04d_%04d.npy" % (fcount,node_idx)),imgs)


x

seriesuid

coordX

coordY

coordZ

diameter_mm

0

1.3.6.1.4.1.14519.5.2.1.6279.6001.179049373636...

56.208405

86.343413

-115.867579

23.350644

img_file = "./images_0000_0000.npy"

fig = plt.figure(figsize=(12,4))
for i in range(3):
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))
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.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.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.imshow(dilation, "gray")
ax4.set_axis_off()
ax4.set_title('Step3,  dilation shrinks dark\nregions and enlarges bright regions.')

# 上一张图中共有三片连续区域，即最外层的体外区域，内部的肺部区域，以及二者之间的身体轮廓区域。这里将其分别标出
labels = measure.label(dilation)
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)

for N in good_labels:

ax6.set_axis_off()
ax6.set_title('Step5, remain the region of interests.')


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

ax1.imshow(imgs_to_process[0,:,:], 'bone')
ax1.set_axis_off()
ax1.set_title("Raw CT image")

ax2.set_axis_off()
ax2.set_title("Pre-processed Images with\nregion of interest in lung")


## 2. 结合深度学习技术的特征提取增强

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

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.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],
b_conv1)
else:
h_conv1 = tf.nn.relu(
tf.nn.conv2d(input=x_tensor,
filter=W_conv1,
strides=[1, 2, 2, 1],
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],
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))
optimizer_F = opt.minimize(cross_entropy_F)
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))
optimizer = opt.minimize(cross_entropy)
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)


0 条评论

• ### 使用腾讯云 GPU 学习深度学习系列之二：Tensorflow 简明原理

本文介绍了如何通过简单的 Python 代码，重点实现深度学习框架的计算图模型，以及自动求导过程。

• ### 使用腾讯云 GPU 学习深度学习系列之六：物体的识别与定位

本文以如何识别马路上的行人、车辆为主题，介绍了基于 Tensorflow 的 SSD 模型如何应用在物体识别定位项目中。

• ### 使用腾讯云 GPU 学习深度学习系列之三：搭建深度神经网络

本文进一步详细介绍了 Tensorflow 中 Keras 工具包提供的几种深度神经网络模块，包括其功能以及用途。

• ### CV 新手避坑指南：计算机视觉常见的8个错误

人类并不是完美的，我们经常在编写软件的时候犯错误。有时这些错误很容易找到：你的代码根本不工作，你的应用程序会崩溃。但有些 bug 是隐藏的，很难发现，这使它们更...

• ### 4. 经典卷积网络之AlexNet

原文：《ImageNet Classification with Deep Convolutional Neural Networks》 我没有读原文，这个已...

• ### 【leetcode刷题】T16-乘积最大子序列

今天分享leetcode第16篇文章，也是leetcode第152题—乘积最大子序列（Maximum Product Subarray），地址是：https:/...

• ### 【整理】上架4.3被拒，我做了这些，正在等结果

刚刚开始触发到4.3被拒这个苹果爸爸安排的隐藏剧情，很懵，不知道它是个啥。经过看官方文档和网上的资料，知道了自己的程序被认定为马甲包（几乎就是一套代码一套UI换...