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

这是《使用腾讯云GPU学习深度学习》系列文章的第四篇，主要举例介绍了深度学习计算过程中的一些数据预处理方法。本系列文章主要介绍如何使用 腾讯云GPU服务器 进行深度学习运算，前面主要介绍原理部分，后期则以实践为主。

上一节，我们基于Keras设计了一个用于 CIFAR-10 数据集的深度学习网络。我们的代码主要包括以下部分：

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

我们注意到，批量输入模块中，实际上就是运用了一个生成器，用来批量读取图片文件，保存成矩阵，直接用于深度神经网络的训练。由于在训练的过程中，图片的特征，是由卷积神经网络自动获取的，因此深度学习通常被认为是一种 端对端（End to end） 的训练方式，期间不需要人为的过多干预。

但是，在实际的运用过程中，这一条并不总是成立。深度神经网络在某些特定的情况下，是需要使用某些特定方法，在批量输入模块之前，对输入数据进行预处理，并且处理结果会极大的改善。

本讲将主要介绍几种数据预处理方法，并且通过这些预处理方法，进行特征提取，提升模型的准确性。

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

这一部分我们举医学影像学的一个例子，以 Kaggle 社区第三届数据科学杯比赛的肺部 CT 扫描结节数据为例，来说明如何进行数据的前处理。以下代码改编自该 kaggle 比赛的官方指导教程,主要是特异性的提取 CT 影像图片在肺部的区域的扫描结果，屏蔽无关区域，进而对屏蔽其他区域后的结果，使用深度学习方法进行进一步分析。

屏蔽的程序本身其实并未用到深度学习相关内容，这里主要使用了skimage库。下面我们详细介绍一下具体方法。

第一步，读取医学影像图像。这里以 LUNA16数据集 中的 1.3.6.1.4.1.14519.5.2.1.6279.6001.179049373636438705059720603192 这张CT 影像数据为例，这张片子可以在这里下载，然后解压缩，用下面的代码分析。其他片子请在 LUNA16 数据集)下载：

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)

简单解释下，首先，CT 影像是一个三维的图像，以三维矩阵的形式保存在 1.3.6.1.4.1.14519.5.2.1.6279.6001.179049373636438705059720603192.raw 这个文件中，.mhd文件则保存了影像文件的基本信息。具体而言，annotation.csv 文件中，图像中结节的坐标是：

这里结节坐标的 coordX~Z 都是物理坐标， .mhd文件保存的，就是从这些物理坐标到 .raw文件中矩阵坐标的映射。于是上面整个函数，其实就是在从 CT 影像仪器的原始文件读取信息，转换物理坐标为矩阵坐标，并且 将结节附近的CT 切片存储成对应的 python 矩阵，用来进行进一步的分析。

然后，我们看一下读取的结果。可见输入文件中标注的结节就在右下方。

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
						 技术沙龙 教程文章 热点综合