医疗方面的图像处理

2019-06-18 15:55:17 weixin_37619439 阅读数 973

卷积神经网络基础(CNN)

为了理解CNN的基础,首先我们需要了解什么是卷积 convolution。

什么是卷积convolution ?

维基百科将卷积定义为“对两个函数(f和g)的数学运算; 它产生第三个函数,通常被视为原始函数之一的修改版本,给出两个函数的逐点乘法的积分作为原始函数之一的转换量的函数。 可以将其简单理解为应用于矩阵的滑动窗口函数。

Convolution with 3×3 Filter. Source: deeplearning.stanford.edu/wiki/index.php/Feature_extraction_using_convolution

上图显示了以绿色应用于矩阵的滑动窗口,其中滑动窗口矩阵为红色。 输出是卷积特征矩阵。 下图显示了两个方形脉冲(蓝色和红色)的卷积和结果。

Source: Wikipedia.

Jeremy Howard在他的MOOC中使用excel表解释了卷积,这是理解基本原理的好方法。 考虑2个矩阵f和g。 f和g的卷积输出是由2个矩阵的点积给出的第三个矩阵“Conv层1”。 2个矩阵的点积是标量,如下所示。 这里可以找到一个很好的数学函数源。

Dot product of 2 matrices.

按照Jeremy的建议,使用excel,我们的输入矩阵是函数f(),滑动窗口矩阵是滤波函数g()。 点积是excel中2个矩阵的和积,如下所示。

Convolution of 2 matrices.

让我们将其扩展为字母“A”的图像。 我们知道任何图像都是由像素组成的。 所以我们的输入矩阵f是“A”。 我们选择滑动窗函数为随机矩阵g。 然后,该矩阵的点积的回旋输出如下所示。

什么是卷积神经网络CNN?

Source: cs231n.github.io/convolutional-networks/

在我看来,简单的卷积神经网络(CNN)是一系列层。 每一层都有一些特定的功能。 每个卷积层都是3维的,因此我们使用volume作为度量。 此外,CNN的每一层通过可微分函数将一个体积的激活转换为另一个体积。 这种功能称为激活或传递功能。

CNN的不同类型的实体是:输入,过滤器(或内核),卷积层,激活层,池化层和批量标准化层。 这些层以不同的排列组合,当然有些规则为我们提供了不同的深度学习架构。

输入层:CNN的常用输入是n维数组。 对于图像,我们输入了3个维度 - 长度,宽度和深度(这是颜色通道)

Source: http://xrds.acm.org/blog/2016/06/convolutional-neural-networks-cnns-illustrated-explanation/

过滤器或内核:如下图l所示,过滤器或内核滑动到图像的每个位置,并计算一个新像素作为这些像素的加权和。 在我们上面的excel示例中,我们的滤波器是g,在输入矩阵f上移动。

source: intellabs.github.io/RiverTrail/tutorial/

卷积层:输入矩阵和内核的点积产生一个新的矩阵,称为卷积矩阵或层。

Source: https://docs.gimp.org/en/plug-in-convmatrix.html

可以在下面找到一个非常好的视觉图表,了解如何填充,跨步和转置。

Source: https://github.com/vdumoulin/conv_arithmetic

激活层:激活功能可以分为两类 - 饱和和不饱和。

饱和激活函数有sigmoid和tanh,而非饱和激活函数是ReLU及其变体。使用非饱和激活函数的优点在于两个方面:

首先是解决所谓的“爆炸/消失梯度”。

第二是加快收敛速度​​。

Sigmoid:获取实值输入并将其压缩到[0,1]之间的范围

σ(x)= 1 /(1 + exp(-x))

tanh:获取实值输入并将其压缩到范围[-1,1]

tanh(x)=2σ(2x) - 1

RELU

ReLU全程叫整流线性单元。它是输入为x的最大函数(x,0),例如来自卷积图像的矩阵。然后,ReLU将矩阵x中的所有负值设置为零,并且所有其他值保持不变.ReLU在卷积之后计算,因此是非线性激活函数,如tanh或sigmoid。这是Geoff Hinton在他的nature报纸中首次讨论的。

ELUs

指数线性单位试图使平均激活接近于零,加速学习速度。 ELU同样通过正值的标识避免消失的梯度。从下图可以看出,ELU比ReLU获得更高的分类精度。

Source: http://image-net.org/challenges/posters/JKU_EN_RGB_Schwarz_poster.pdf [15 layer CNN with stacks of (1×96×6, 3×512×3, 5×768×3, 3×1024×3, 2×4096×FC, 1×1000×FC) layers×units×receptive fields or fully-connected (FC). 2×2 max-pooling with a stride of 2 after each stack, spatia

Source: Wikipedia.

Leaky ReLUs

与ReLU相反,relu完全舍掉负数部分,leaky ReLU则为其指定非零斜率。 泄漏整流线性激活首先在声学模型中引入(Maas等,2013)。 在数学上,我们有:

Source: Empirical Evaluation of Rectified Activations in Convolution Network.

其中ai是范围内的固定参数(1,+无穷大)。

参数整流线性单元(PReLU)

PReLU可以被视为Leaky ReLU的变种。 在PReLU中,负面部分的斜率是从数据中学习而不是预定义的。 作者声称,PReLU是超越ImageNet分类(Russakovsky等,2015)任务的人类表现的关键因素。 除了ai是通过反向传播在训练中学习,PReLU与ReLU相同。

随机泄漏整流线性单元(RReLU)

随机整流线性单元(RReLU)也是Leaky ReLU的变体。 在RReLU中,负部分的斜率在训练中的给定范围内随机化,然后在测试中固定。 RReLU的亮点在于,在训练过程中,aji是从均匀分布U(l,u)中采样的随机数。 在形式上:


ReLU,Leaky ReLU,PReLU和RReLU之间的比较如下所示。

Source :https://arxiv.org/pdf/1505.00853.pdf ReLU, Leaky ReLU, PReLU and RReLU. For PReLU, ai is learned and for Leaky ReLU ai is fixed. For RReLU, aji is a random variable keeps sampling in a given range, and remains fixed in testing.

噪声激活函数

这些激活函数,扩展到包括高斯噪声。 在这里可以找到关于噪声如何帮助的良好理解。

Source: Wikipedia.

池化层(Pooling layer)

池化层的目标是逐渐减小矩阵的空间大小,以减少网络中的参数和计算量,从而也控制过度拟合。 Pooling Layer在输入的每个深度切片上独立运行,并使用MAX或Average操作在空间上调整其大小。 最常见的形式是一个池化层,其过滤器大小为2*2,步长为2,沿着宽,高,输入为2,对每个深度切片进行降低采样,丢弃75%的激活。 在这种情况下,每个MAX操作将采用最多超过4个数字(在某个深度切片中的小2x2区域)。 深度维度保持不变。 简单来说,池化层:

Source: http://cs231n.github.io/convolutional-networks/#pool

Source: https://ujjwalkarn.me/2016/08/11/intuitive-explanation-convnets/

注意:这里我们以步长为2,滑动2*2的窗口并且去该区域最大值。

批量标准化层 Batch Normalization layer

批量标准化是标准化每个中间层的有效方式,包括权重和激活函数。 batchnorm有两个主要的好处:

将batchnorm添加到模型中可以使训练速度提高10倍或更多

由于归一化极大地降低了少量外围输入过度影响训练的能力,因此也倾向于减少过度拟合。

全连接层 fully connected layer

完全连接层是传统的多层感知器,在输出层使用softmax激活函数。 术语“完全连接”意味着前一层中的每个神经元都连接到下一层的每个神经元。 softmax函数是泛华函数的推广,其将任意实数值的K维向量“压缩”到范围(0,1)中的实数值的K维向量,其加起来为1。

Source: Wikipedia

Softmax激活通常用于最终的完全连接层,以获得概率,它的值为0到1之间。

现在,我们了解CNN中的不同层。 有了这些知识,我们将在下一篇文章中使用Keras开发肺癌检测所需的深度学习架构。

 

参考:

  1. Jeremy Howard’s MOOC (course.fast.ai)
  2. http://www.wildml.com/2015/11/understanding-convolutional-neural-networks-for-nlp/
  3. https://medium.com/towards-data-science/linear-algebra-cheat-sheet-for-deep-learning-cd67aba4526c
  4. https://ujjwalkarn.me/2016/08/11/intuitive-explanation-convnets/
  5. https://medium.com/technologymadeeasy/the-best-explanation-of-convolutional-neural-networks-on-the-internet-fbb8b1ad5df8
  6. http://image-net.org/challenges/posters/JKU_EN_RGB_Schwarz_poster.pdf
2017-06-09 15:02:13 u013635029 阅读数 11934

常见医疗扫描图像处理步骤


一 下载必要包

skimage需要更新到最新0.13版本,否则会报错,ImportError: cannot import name label。

sudo pip install scikit-image  -U -i  https://pypi.tuna.tsinghua.edu.cn/simple

二 数据格式

2.1 dicom

当前医疗竞赛官方数据集格式 data-science-bowl-2017。数据列表如下:
这里写图片描述
后缀为 .dcm。
每个病人的一次扫描CT(scan)可能有几十到一百多个dcm数据文件(slices)。可以使用 python的dicom包读取,读取示例代码如下:

dicom.read_file('/data/lung_competition/stage1/7050f8141e92fa42fd9c471a8b2f50ce/498d16aa2222d76cae1da144ddc59a13.dcm')

一般使用其pixl_array数据

slices = [dicom.read_file(os.path.join(folder_name,filename)) for filename in os.listdir(folder_name)]
slices = np.stack([s.pixel_array for s in slices])

2.2 mhd格式

mhd格式是另外一种数据格式,来源于(LUNA2016)[https://luna16.grand-challenge.org/data/]。每个病人一个mhd文件和一个同名的raw文件。如下:

这里写图片描述

一个mhd通常有几百兆,对应的raw文件只有1kb。mhd文件需要借助python的SimpleITK包来处理。SimpleITK 示例代码如下:

import SimpleITK as sitk
itk_img = sitk.ReadImage(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)

需要注意的是,SimpleITK的img_array的数组不是直接的像素值,而是相对于CT扫描中原点位置的差值,需要做进一步转换。转换步骤参考 SimpleITK图像转换

2.3 查看CT扫描文件软件

一个开源免费的查看软件 mango
这里写图片描述

三 dicom格式数据处理过程

3.1 处理思路

首先,需要明白的是医学扫描图像(scan)其实是三维图像,使用代码读取之后开源查看不同的切面的切片(slices),可以从不同轴切割

这里写图片描述

如下图展示了一个病人CT扫描图中,其中部分切片slices

这里写图片描述

其次,CT扫描图是包含了所有组织的,如果直接去看,看不到任何有用信息。需要做一些预处理,预处理中一个重要的概念是放射剂量,衡量单位为HU(Hounsfield Unit),下表是不同放射剂量对应的组织器官

这里写图片描述

Hounsfield Unit = pixel_value * rescale_slope + rescale_intercept

一般情况rescale slope = 1, intercept = -1024。
灰度值是pixel value经过重重LUT转换得到的用来进行显示的值,而这个转换过程是不可逆的,也就是说,灰度值无法转换为ct值。只能根据窗宽窗位得到一个大概的范围。 pixel value经过modality lut得到Hu,但是怀疑pixelvalue的读取出了问题。dicom文件中存在(0028,0106)(0028,0107)两个tag,分别是最大最小pixel value,可以用来检验你读取的pixel value 矩阵是否正确。
LUT全称look up table,实际上就是一张像素灰度值的映射表,它将实际采样到的像素灰度值经过一定的变换如阈值、反转、二值化、对比度调整、线性变换等,变成了另外一 个与之对应的灰度值,这样可以起到突出图像的有用信息,增强图像的光对比度的作用。
首先去除超过 -2000的pixl_array,CT扫描边界之外的像素值固定为-2000(dicom和mhd都是这个值)。第一步是设定这些值为0,当前对应为空气(值为0)

slices[slices == -2000] = 0

上表中肺部组织的HU数值为-500,但通常是大于这个值,比如-320、-400。挑选出这些区域,然后做其他变换抽取出肺部像素点。抽取出这些特征,最后将其存储为ndarray的npy格式,供给卷积神经网络。

3.2 先载入必要的包

# -*- coding:utf-8 -*-
'''
this script is used for basic process of lung 2017 in Data Science Bowl
'''

import glob
import os
import pandas as pd
import SimpleITK as sitk

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import skimage, os
from skimage.morphology import ball, disk, dilation, binary_erosion, remove_small_objects, erosion, closing, reconstruction, binary_closing
from skimage.measure import label,regionprops, perimeter
from skimage.morphology import binary_dilation, binary_opening
from skimage.filters import roberts, sobel
from skimage import measure, feature
from skimage.segmentation import clear_border
from skimage import data
from scipy import ndimage as ndi
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import dicom
import scipy.misc
import numpy as np

DICOM是医学图像中标准文件,这些文件包含了诸多的
元数据信息(比如像素尺寸,每个维度的一像素代表真实世界里的长度)。如下代码是载入一个扫描面,包含了多个切片(slices),我们仅简化的将其存储为python列表。数据集中每个目录都是一个扫描面集(一个病人)。有个元数据域丢失,即Z轴方向上的像素尺寸,也即切片的厚度 。所幸,我们可以用其他值推测出来,并加入到元数据中。

# Load the scans in given folder path
def load_scan(path):
    slices = [dicom.read_file(path + '/' + s) for s in os.listdir(path)]
    slices.sort(key = lambda x: int(x.ImagePositionPatient[2]))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)

    for s in slices:
        s.SliceThickness = slice_thickness

    return slices

3.3 像素转换为HU单元

有些扫描面有圆柱形扫描边界,但是输出图像是正方形。边界之外的像素值固定为-2000,。第一步是设定这些值为0,当前对应为空气(值为0)。然后回到HU单元,乘以rescale比率并加上intercept(存储在扫描面的元数据中)。

def get_pixels_hu(slices):
    image = np.stack([s.pixel_array for s in slices])
    # Convert to int16 (from sometimes int16), 
    # should be possible as values should always be low enough (<32k)
    image = image.astype(np.int16)
    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    image[image == -2000] = 0

# Convert to Hounsfield units (HU)
    for slice_number in range(len(slices)):

        intercept = slices[slice_number].RescaleIntercept
        slope = slices[slice_number].RescaleSlope

        if slope != 1:
            image[slice_number] = slope * image[slice_number].astype(np.float64)
            image[slice_number] = image[slice_number].astype(np.int16)

        image[slice_number] += np.int16(intercept)

    return np.array(image, dtype=np.int16)

可以查看病人的扫描HU值分布情况

first_patient = load_scan(INPUT_FOLDER + patients[0])
first_patient_pixels = get_pixels_hu(first_patient)
plt.hist(first_patient_pixels.flatten(), bins=80, color='c')
plt.xlabel("Hounsfield Units (HU)")
plt.ylabel("Frequency")
plt.show()

这里写图片描述

3.4重新采样

不同扫描面的像素尺寸、粗细粒度是不同的。这不利于我们进行CNN任务,我们可以使用同构采样。
一个扫描面的像素区间可能是[2.5,0.5,0.5],即切片之间的距离为2.5mm。可能另外一个扫描面的范围是[1.5,0.725,0.725]。这可能不利于自动分析。
常见的处理方法是从全数据集中以固定的同构分辨率重新采样。如果我们选择,将所有的东西采样为1mmx1mmx1mm像素,我们可以使用3D卷积网络

def resample(image, scan, new_spacing=[1,1,1]):
    # Determine current pixel spacing
    spacing = map(float, ([scan[0].SliceThickness] + scan[0].PixelSpacing))
    spacing = np.array(list(spacing))
    resize_factor = spacing / new_spacing
    new_real_shape = image.shape * resize_factor
    new_shape = np.round(new_real_shape)
    real_resize_factor = new_shape / image.shape
    new_spacing = spacing / real_resize_factor

    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')

    return image, new_spacing

现在重新取样病人的像素,将其映射到一个同构分辨率 1mm x1mm x1mm。

pix_resampled, spacing = resample(first_patient_pixels, first_patient, [1,1,1])

输出肺部扫描的3D图像方法

def plot_3d(image, threshold=-300):

    # Position the scan upright, 
    # so the head of the patient would be at the top facing the camera
    p = image.transpose(2,1,0)

    verts, faces = measure.marching_cubes(p, threshold)
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    # Fancy indexing: `verts[faces]` to generate a collection of triangles
    mesh = Poly3DCollection(verts[faces], alpha=0.1)
    face_color = [0.5, 0.5, 1]
    mesh.set_facecolor(face_color)
    ax.add_collection3d(mesh)
    ax.set_xlim(0, p.shape[0])
    ax.set_ylim(0, p.shape[1])
    ax.set_zlim(0, p.shape[2])
    plt.show()

打印函数有个阈值参数,来打印特定的结构,比如tissue或者骨头。400是一个仅仅打印骨头的阈值(HU对照表)

plot_3d(pix_resampled, 400)

这里写图片描述

3.5 输出一个病人scans中所有切面slices

def plot_ct_scan(scan):
    '''
            plot a few more images of the slices
    :param scan:
    :return:
    '''
    f, plots = plt.subplots(int(scan.shape[0] / 20) + 1, 4, figsize=(50, 50))
    for i in range(0, scan.shape[0], 5):
        plots[int(i / 20), int((i % 20) / 5)].axis('off')
        plots[int(i / 20), int((i % 20) / 5)].imshow(scan[i], cmap=plt.cm.bone)

此方法的效果示例如下:

这里写图片描述

3.6 定义分割出CT切面里面肺部组织的函数

def get_segmented_lungs(im, plot=False):

    '''
    This funtion segments the lungs from the given 2D slice.
    '''
    if plot == True:
        f, plots = plt.subplots(8, 1, figsize=(5, 40))
    '''
    Step 1: Convert into a binary image.
    '''
    binary = im < 604
    if plot == True:
        plots[0].axis('off')
        plots[0].set_title('binary image')
        plots[0].imshow(binary, cmap=plt.cm.bone)

    '''
    Step 2: Remove the blobs connected to the border of the image.
    '''
    cleared = clear_border(binary)
    if plot == True:
        plots[1].axis('off')
        plots[1].set_title('after clear border')
        plots[1].imshow(cleared, cmap=plt.cm.bone)

    '''
    Step 3: Label the image.
    '''
    label_image = label(cleared)
    if plot == True:
        plots[2].axis('off')
        plots[2].set_title('found all connective graph')
        plots[2].imshow(label_image, cmap=plt.cm.bone)
    '''
    Step 4: Keep the labels with 2 largest areas.
    '''
    areas = [r.area for r in regionprops(label_image)]
    areas.sort()
    if len(areas) > 2:
        for region in regionprops(label_image):
            if region.area < areas[-2]:
                for coordinates in region.coords:
                       label_image[coordinates[0], coordinates[1]] = 0
    binary = label_image > 0
    if plot == True:
        plots[3].axis('off')
        plots[3].set_title(' Keep the labels with 2 largest areas')
        plots[3].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 5: Erosion operation with a disk of radius 2. This operation is
    seperate the lung nodules attached to the blood vessels.
    '''
    selem = disk(2)
    binary = binary_erosion(binary, selem)
    if plot == True:
        plots[4].axis('off')
        plots[4].set_title('seperate the lung nodules attached to the blood vessels')
        plots[4].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 6: Closure operation with a disk of radius 10. This operation is
    to keep nodules attached to the lung wall.
    '''
    selem = disk(10)
    binary = binary_closing(binary, selem)
    if plot == True:
        plots[5].axis('off')
        plots[5].set_title('keep nodules attached to the lung wall')
        plots[5].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 7: Fill in the small holes inside the binary mask of lungs.
    '''
    edges = roberts(binary)
    binary = ndi.binary_fill_holes(edges)
    if plot == True:
        plots[6].axis('off')
        plots[6].set_title('Fill in the small holes inside the binary mask of lungs')
        plots[6].imshow(binary, cmap=plt.cm.bone)
    '''
    Step 8: Superimpose the binary mask on the input image.
    '''
    get_high_vals = binary == 0
    im[get_high_vals] = 0
    if plot == True:
        plots[7].axis('off')
        plots[7].set_title('Superimpose the binary mask on the input image')
        plots[7].imshow(im, cmap=plt.cm.bone)

    return im

此方法每个步骤对图像做不同的处理,效果如下:

这里写图片描述

3.7 肺部图像分割

为了减少有问题的空间,我们可以分割肺部图像(有时候是附近的组织)。这包含一些步骤,包括区域增长和形态运算,此时,我们只分析相连组件。
步骤如下:

  • 阈值图像(-320HU是个极佳的阈值,但是此方法中不是必要)
  • 处理相连的组件,以决定当前患者的空气的标签,以1填充这些二值图像
  • 可选:当前扫描的每个轴上的切片,选定最大固态连接的组织(当前患者的肉体和空气),并且其他的为0。以掩码的方式填充肺部结构。
  • 只保留最大的气袋(人类躯体内到处都有气袋)
def largest_label_volume(im, bg=-1):
    vals, counts = np.unique(im, return_counts=True)
    counts = counts[vals != bg]
    vals = vals[vals != bg]
    if len(counts) > 0:
        return vals[np.argmax(counts)]
    else:
        return None
def segment_lung_mask(image, fill_lung_structures=True):

    # not actually binary, but 1 and 2. 
    # 0 is treated as background, which we do not want
    binary_image = np.array(image > -320, dtype=np.int8)+1
    labels = measure.label(binary_image)

    # Pick the pixel in the very corner to determine which label is air.
    #   Improvement: Pick multiple background labels from around the patient
    #   More resistant to "trays" on which the patient lays cutting the air 
    #   around the person in half
    background_label = labels[0,0,0]

    #Fill the air around the person
    binary_image[background_label == labels] = 2


    # Method of filling the lung structures (that is superior to something like 
    # morphological closing)
    if fill_lung_structures:
        # For every slice we determine the largest solid structure
        for i, axial_slice in enumerate(binary_image):
            axial_slice = axial_slice - 1
            labeling = measure.label(axial_slice)
            l_max = largest_label_volume(labeling, bg=0)

            if l_max is not None: #This slice contains some lung
                binary_image[i][labeling != l_max] = 1

    binary_image -= 1 #Make the image actual binary
    binary_image = 1-binary_image # Invert it, lungs are now 1

    # Remove other air pockets insided body
    labels = measure.label(binary_image, background=0)
    l_max = largest_label_volume(labels, bg=0)
    if l_max is not None: # There are air pockets
        binary_image[labels != l_max] = 0

    return binary_image

查看切割效果

segmented_lungs = segment_lung_mask(pix_resampled, False)
segmented_lungs_fill = segment_lung_mask(pix_resampled, True)
plot_3d(segmented_lungs, 0)

这里写图片描述

我们可以将肺内的结构也包含进来(结节是固体),不仅仅只是肺部内的空气

plot_3d(segmented_lungs_fill, 0)

这里写图片描述

使用掩码时,要注意首先进行形态扩充(python的skimage的skimage.morphology)操作(即使用圆形kernel,结节是球体),参考 python形态操作。这会在所有方向(维度)上扩充掩码。仅仅肺部的空气+结构将不会包含所有结节,事实上有可能遗漏黏在肺部一侧的结节(这会经常出现,所以建议最好是扩充掩码)。

3.8 数据预处理

归一化处理
当前的值范围是[-1024,2000]。任意大于400的值都是我们所感兴趣的,因为它们都是不同反射密度下的骨头。LUNA16竞赛中常用来做归一化处理的阈值集是-1000和400.以下代码

#归一化
MIN_BOUND = -1000.0
MAX_BOUND = 400.0

def normalize(image):
    image = (image - MIN_BOUND) / (MAX_BOUND - MIN_BOUND)
    image[image>1] = 1.
    image[image<0] = 0.
    return image

0值中心化
简单来说就是所有像素值减去均值。LUNA16竞赛中的均值大约是0.25.
不要对每一张图像做零值中心化(此处像是在kernel中完成的)CT扫描器返回的是校准后的精确HU计量。不会出现普通图像中会出现某些图像低对比度和明亮度的情况

PIXEL_MEAN = 0.25

def zero_center(image):
    image = image - PIXEL_MEAN
    return image

3.9 存储每个病人scan的所有slice肺部特征

下述代码为main函数中整体过程

 lung_head_dir = '/data/lung_competition/stage1/'
    img_save_head = '/data/lung_competition/roi_images/'
    patient_scans_list = [x for x in os.listdir(lung_head_dir) if os.path.isdir(os.path.join(lung_head_dir,x))]
    #patient_labels_file = '/data/lung_competition/stage1_labels.csv'
    #patient_labels = pd.read_csv(patient_labels_file)
    selem = ball(2)
    #print(patient_scans_list)
    for scan in patient_scans_list:
        scan_files = os.path.join(lung_head_dir,scan)
        ct_scan = read_ct_scan(scan_files)
        save_npy_path = os.path.join(img_save_head, scan)
        if os.path.exists(save_npy_path):
            os.mkdir(save_npy_path)

存储病人scan的所有slice肺部特征的关键代码为

segmented_ct_scan = segment_lung_from_ct_scan(ct_scan)
save_npy = save_npy_path+'.npy'
np.save(save_npy,segmented_ct_scan)
print('file %s saved ..'%save_npy)
plot_ct_scan(segmented_ct_scan)

四 mhd格式数据处理过程

4.1 处理思路

mhd的数据只是格式与dicom不一样,其实质包含的都是病人扫描。处理MHD需要借助SimpleIKT这个包,处理思路详情可以参考Data Science Bowl2017的toturail Data Science Bowl 2017。需要注意的是MHD格式的数据没有HU值,它的值域范围与dicom很不同。
我们以LUNA2016年的数据处理流程为例。参考代码为 LUNA2016数据切割

4.2 载入必要的包

import SimpleITK as sitk
import numpy as np
import csv
from glob import glob
import pandas as pd
file_list=glob(luna_subset_path+"*.mhd")
#####################
#
# Helper function to get rows in data frame associated 
# with each file
def get_filename(case):
    global file_list
    for f in file_list:
        if case in f:
            return(f)
#
# The locations of the nodes
df_node = pd.read_csv(luna_path+"annotations.csv")
df_node["file"] = df_node["seriesuid"].apply(get_filename)
df_node = df_node.dropna()
#####
#
# Looping over the image files
#
fcount = 0
for img_file in file_list:
    print "Getting mask for image file %s" % img_file.replace(luna_subset_path,"")
    mini_df = df_node[df_node["file"]==img_file] #get all nodules associate with file
    if len(mini_df)>0:       # some files may not have a nodule--skipping those 
        biggest_node = np.argsort(mini_df["diameter_mm"].values)[-1]   # just using the biggest node
        node_x = mini_df["coordX"].values[biggest_node]
        node_y = mini_df["coordY"].values[biggest_node]
        node_z = mini_df["coordZ"].values[biggest_node]
        diam = mini_df["diameter_mm"].values[biggest_node]

4.3 LUNA16的MHD格式数据的值

一直在寻找MHD格式数据的处理方法,对于dicom格式的CT有很多论文根据其HU值域可以轻易地分割肺、骨头、血液等,但是对于MHD没有这样的参考。从LUNA16论坛得到的解释是,LUNA16的MHD数据已经转换为HU值了,不需要再使用slope和intercept来做rescale变换了。此论坛主题下,有人提出MHD格式没有提供pixel spacing(mm) 和 slice thickness(mm) ,而标准文件annotation.csv文件中结节的半径和坐标都是mm单位,最后确认的是MHD格式文件中只保留了体素尺寸以及坐标原点位置,没有保存slice thickness。如此来说,dicom才是原始数据格式。

4.4 坐标体系变换

MHD值得坐标体系是体素,以mm为单位。结节的位置是CT scanner坐标轴里面相对原点的mm值,需要将其转换到真实坐标轴位置,可以使用SimpleITK包中的 GetOrigin() GetSpacing()。图像数据是以512x512数组的形式给出的。
坐标变换如

这里写图片描述

相应的代码处理如下:

k_img = sitk.ReadImage(img_file) 
img_array = sitk.GetArrayFromImage(itk_img) # indexes are z,y,x (notice the ordering)
center = np.array([node_x,node_y,node_z])   # nodule center
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)
v_center =np.rint((center-origin)/spacing)  # nodule center in voxel space (still x,y,z ordering)

在LUNA16的标注CSV文件中标注了结节中心的X,Y,Z轴坐标,但是实际取值的时候取的是Z轴最后三层的数组(img_array)。
下述代码只提取了包含结节的最后三个slice的数据

i=0
for i_z in range(int(v_center[2])-1,int(v_center[2])+2):
    mask = make_mask(center,diam,i_z*spacing[2]+origin[2],width,height,spacing,origin)
    masks[i] = mask
    imgs[i] = matrix2int16(img_array[i_z])
    i+=1
np.save(output_path+"images_%d.npy" % (fcount) ,imgs)
np.save(output_path+"masks_%d.npy" % (fcount) ,masks)

4.5 查看结节

以下代码用于查看原始CT和结mask

import matplotlib.pyplot as plt
imgs = np.load(output_path+'images_0.npy')
masks = np.load(output_path+'masks_0.npy')
for i in range(len(imgs)):
    print "image %d" % i
    fig,ax = plt.subplots(2,2,figsize=[8,8])
    ax[0,0].imshow(imgs[i],cmap='gray')
    ax[0,1].imshow(masks[i],cmap='gray')
    ax[1,0].imshow(imgs[i]*masks[i],cmap='gray')
    plt.show()
    raw_input("hit enter to cont : ")

这里写图片描述

4.6 包含结节位置信息的mhd格式数据特征

参考一个开源代码实现 UNET训练肺组织结节分割

2019-10-02 16:03:57 qq_41985559 阅读数 789

这些方面形成此领域的三个主要过程——图像形成、图像计算和图像管理。医疗图像处理是一个非常复杂的跨学科领域,涵盖从数学、计算

机科学到物理学和医学的众多科学学科。
在这里插入图片描述
分析法的典型示例包括广泛用于断层扫描的滤波反投影(FBP);在MRI中尤为重要的傅里叶变换(FT);以及延时叠加(DAS)波束成

型,这是超声检查中一种不可或缺的技术。这些算法在所需的处理能力和计算时间方面精巧而高效。

核磁共振扫描仪(MRI)是使用非常强的磁场和无线电波,这些磁场和无线电波与组织中的质子相互作用,产生一个信号,然后经过处理,

形成人体图像。质子(氢原子)可以被认为是微小的条形磁铁,有北极和南极,绕轴旋转——就像行星一样。正常情况下,质子是随机排列

的,但当施加强磁场时,质子磁场方向会与这个磁场方向对齐。

由于所涉及的强磁场,该设备不能用于心脏起搏器可能会被破坏的患者,或金属植入物弹片可能会在手术过程中被磁铁吸引和

移动的患者。此外,铁磁性物体会被磁铁强烈吸引,并对抛射物造成严重的危险。因为这些原因,这些物体被禁止靠近核磁共振设

备。

-滤波反投影算法:

  1. 计算每一个投影的一维傅里叶变换;
    
  2. 用滤波函数 乘以每一个傅里叶变换;
    
  3. 得到每一个滤波后的变换的一维反傅里叶变换;
    
  4. 对步骤3得到的所以一维反变换积分(求和)
    

在这里插入图片描述
X光特别适合用于骨骼外伤,如果怀疑伤到了骨头,优选X光,检查结果快速易得;也适用于胸部的粗略检查,X光胸片可粗略检查心脏、主

动脉、肺、胸膜、肋骨等,可检查有无肺纹理增多、肺部钙化等。

在这里插入图片描述

CT是用X射线束对人体某部位一定厚度地扫描,根据人体不同组织对X线的吸收与透过率的不同,应用灵敏度极高的仪器对人体进行测量,

然后将测量所获取的数据输入电子计算机,电子计算机对数据进行处理后,就可摄下人体被检查部位的断面或立体的图像,发现体内任何部

位的细小病变。它是通过数据计算得到的重建图像。

在这里插入图片描述
B超的原理是用超声波穿透人体,当声波遇到人体组织时会产生反射波,通过计算反射波成像。就像挑西瓜一样,边敲边看显示病灶情况。B

超的原理是用超声波穿透人体,当声波遇到人体组织时会产生反射波,通过计算反射波成像。就像挑西瓜一样,边敲边看显示病灶情况。

在这里插入图片描述

迭代法有很多种,包括最大似然期望最大化(MLEM)、最大后验(MAP)、代数重建(ARC)技术以及许多其他目前广泛应用于医疗成像

模式的方法。

微信公众号“计算机基础学”关注我哟

2018-11-24 17:06:23 weixin_43789163 阅读数 3085

电子科技大学 格拉斯哥学院 2017级 徐冰砚

浅谈图像处理在医学方面的应用

1、背景:

在上学期的新生研讨课中,曾兵院长介绍了图像处理的相关原理和应用。图像处理(image processing)是一种用计算机对图像进行分析,以达到所需结果的技术。在获得图像之后,需要用专门的设备将其数字化,即通过取样和量化过程将一个以自然形式存在的图像变换为适合计算机处理的数字形式。图像在计算机内部被表示为一个数字矩阵,矩阵中每一元素称为像素。
图像处理有各种应用途径,卫星图像处理、面孔识别特征识别、显微图像处理等等,给我印象最为深刻的是图像处理在医学方面的应用。目前临床广泛使用的医学成像模式主要分为四类:X- 射线成像、核磁共振成像(MRI)、核医学成像(NMI)、超声波成像(Ultrasonic Imaging)。图像分析可以将医学模拟图像转化为数字图像,开展了计算机辅助诊断(computer aided diagnosis)的初步研究,在一定程度上可以辅助医生分析医学图像,从而排除人为主观因素,提高诊断准确性和效率。

2、医学图像处理技术:

(1) 图像分割:
由于人体的组织器官不均匀、器官蠕动等造成医学图像一般具有噪声、病变组织边缘模糊等特点, 医学图像分割技术的目的就是将图像中感兴趣的区域清楚的提取出来, 这样就能为后续的定量、定性分析提供图像基础,同时它也是三维可视化的基础。现在有的图像分割方法有如下几种:基于阈值的分割方法、基于区域的分割方法、基于边缘的分割方法以及基于特定理论的分割方法等。

(2) 图像配准和图像融合:
医学图像配准是指对于一幅医学图像通过一种或一系列的空间变换,使它与另一幅医学图像上的对应点达到空间上的一致。配准的结果应使两幅图像上所有的解剖点,或至少是所有具有诊断意义的点及手术感兴趣的点都达到匹配,配准处理一般可以分为图像变换和图像定位两种。
医学图像在空间域配准之后,就可以进行图像融合,融合图像的创建又分为图像数据的融合与融合图像的显示两部分来完成。图像融合的目的是通过综合处理应用这些成像设备所得信息以获得新的有助于临床诊断的信息。利用可视化软件,对多种模态的图像进行图像融合,可以准确地确定病变体的空间位置、大小、几何形状及它与周围生物组织之间的空间关系,从而及时高效地诊断疾病,也可以用在手术计划的制定、病理变化的跟踪、治疗效果的评价等方面。

(3) 伪彩色处理技术:

伪彩色图像处理技术是将黑白图像经过处理变为彩色图像, 可以充分发挥人眼对彩色的视觉能力, 从而使观察者能从图像中取得更多的信息。经过伪彩色处理技术, 提高了对图像特征的识别。临床研究对CT、MRI、B 超和电镜等图片均进行了伪彩色技术的尝试, 取得了良好的效果, 部分图片经过处理后可以显现隐性病灶。

3、总结:

随着医疗技术的蓬勃发展,对医学图像处理提出的要求也越来越高。医学图像处理技术发展至今,仍然还有很多亟待解决的问题。有效地提高医学图像处理技术的水平、与多学科理论的交叉融合、医务人员和计算机理论技术人员之间的交流就显得越来越重要。总之,医学图像作为现代医疗诊断的重要依据,必将在医药信息研究领域和计算机图像处理领域受到更多的关注。

4、参考文献:
[1]王新成.高级图像处理技术[M].北京:中国科学技术出版社,2001.
[2]丁莹.图像配准技术在医学图像处理中的应用研究[M].长春理工大学,2006.12.
[3]田捷.医学影像处理与分析[M], 电子工业出版社, 2003.
[4]田娅, 饶妮妮, 蒲立新.国内医学图像处理技术的最新动态[J].电子科技大学学报, 2002(5): 485- 489.
图片来源:(https://baike.baidu.com/item/医学图像分析/3939451#2)

2016-02-26 17:48:13 baimafujinji 阅读数 70342

什么是数字图像处理?历史、以及它所研究的内容。

 

说起图像处理,你会想到什么?你是否真的了解这个领域所研究的内容。纵向来说,数字图像处理研究的历史相当悠久;横向来说,数字图像处理研究的话题相当广泛。

数字图像处理的历史可以追溯到近百年以前,大约在1920年的时候,图像首次通过海底电缆从英国伦敦传送到美国纽约。图像处理的首次应用是为了改善伦敦和纽约之间海底电缆发送的图片质量,那时就应用了图像编码,被编码后的图像通过海底电缆传送至目的地,再通过特殊设备进行输出。这是一次历史性的进步,传送一幅图片的时间从原来的一个多星期减少到了3小时。

1950年,美国的麻省理工学院制造出了第一台配有图形显示器的电子计算机——旋风I号(Whirlwind I)。旋风I号的显示器使用一个类似于示波器的阴极射线管(Cathode Ray Tube,CRT)来显示一些简单的图形。1958年美国Calcomp公司研制出了滚筒式绘图仪,GerBer公司把数控机床发展成为平板式绘图仪。在这一时期,电子计算机都主要应用于科学计算,而为这些计算机配置的图形设备也仅仅是作为一种简单的输出设备。

随着计算机技术的进步,数字图像处理技术也得到了很大的发展。1962年,当时还在麻省理工学院攻读博士学位的伊凡·苏泽兰(Ivan Sutherland)成功开发了具有划时代意义的“画板”(Sketchpad)程式。而这正是有史以来第一个交互式绘图系统,同时这也是交互式电脑绘图的开端。从此计算机和图形图像被更加紧密地联系到了一起。鉴于伊凡·苏泽兰为计算机图形学创立所做出的杰出贡献,他于1988年被授予计算机领域最高奖——图灵奖。

1964年,美国加利福尼亚的喷气推进实验室用计算机对“旅行者七号”太空船发回的大批月球照片进行处理,以校正航天器上摄影机中各种类型的图像畸变,收到了明显的效果。在后来的宇航空间技术中,数字图像处理技术都发挥了巨大的作用。

到了20世纪60年代末期,数字图像处理已经形成了比较完善的学科体系,这套理论在20世纪70年代发展得十分迅速,并开始应用于医学影像和天文学等领域。1972年,美国物理学家阿伦·马克利奥德·柯麦科(Allan MacLeodCormack)和英国电机工程师戈弗雷·纽博尔德·豪恩斯弗尔德(Godfrey Newbold Housfield)发明了轴向断层术,并将其用于头颅诊断。世界第一台X射线计算机轴向断层摄影装置由EMI公司研制成功,这也就是人们通常所说的CT(Computer Tomograph)。CT可通过一些算法用感知到的数据去重建通过物体的“切片”图像。这些图像组成了物体内部的再现图像,也就是根据人的头部截面的投影,经计算机处理来进行图像重建。鉴于CT对于医学诊断技术的发展所起到的巨大推动作用,柯麦科和豪恩斯弗尔德于1979年获得了诺贝尔生理或医学奖。

随后在2003年,诺贝尔生理或医学奖的殊荣再次授予了两位在医疗影像设备研究方面做出杰出贡献的科学家——美国化学家保罗·劳特伯尔(Paul Lauterbur)和英国物理学家彼得·曼斯菲尔(Peter Mansfield)。两位获奖者在利用磁共振成像(Magnetic Resonance Imaging,MRI)显示不同结构方面分别取得了开创性成就。瑞典卡罗林斯卡医学院称,这两位科学家在MRI领域的开创性工作,代表了医学诊疗和研究的重大突破。而事实上,核磁共振的成功同样也离不开数字图像处理方面的发展。即使在今天,诸如MRI图像降噪等问题依然是数字图像处理领域的热门研究方向。

说到数字图像的发展历程,还有一项至关重要的成果不得不提,那就是电荷耦合元件(Charge-coupled Device,CCD)。CCD最初是由美国贝尔实验室的科学家维拉德·波义耳(Willard Sterling Boyle)和乔治·史密斯(George Elwood Smith)于1969年发明的。CCD的作用就像胶片一样,它能够把光学影像转化为数字信号。今天人们所广泛使用的数码照相机、数码摄影机和扫描仪都是以CCD为基础发展而来的。换句话说,我们现在所研究的数字图像主要也都是通过CCD设备获取的。由于波义耳和史密斯在CCD研发上所做出的巨大贡献,他们两人共同荣获了2009年度的诺贝尔物理学奖。

数字图像处理在今天是非常热门的技术之一,生活中无处不存在着它的影子,可以说它是一种每时每刻都在改变着人类生活的技术。但长久以来,很多人对数字图像处理存在着较大的曲解,人们总是不自觉地将图像处理和Photoshop联系在一起。大名鼎鼎的Photoshop无疑是当前使用最为广泛的图像处理工具。类似的软件还有Corel公司生产的CorelDRAW等软件。

尽管Photoshop是一款非常优秀的图像处理软件,但它的存在并不代表数字图像处理的全部理论与方法。它所具有的功能仅仅是数字图像处理中的一部分。总的来说,数字图像处理研究的内容主要包括如下几个方面:

  • 1)图像获取和输出
  • 2)图像编码和压缩
  • 3)图像增强与复原
  • 4)图像的频域变换
  • 5)图像的信息安全
  • 6)图像的区域分割
  • 7)图像目标的识别
  • 8)图像的几何变换

但图像处理的研究内容,又不仅限于上述内容!所以说图像处理的研究话题是相当宽泛的。那现在图像处理都应用在哪些领域呢?或许我们可能熟知的例子有(当然,你应该还能举出更多例子):

  • 1)一些专业图像处理软件:Photoshop、CorelDRAW……
  • 2)一些手机APP应用:美图秀秀、玩图……
  • 3)一些医学图像处理应用:MRI、彩超图像处理……
  • 4)一些制造业上的应用:元器件检测、瑕疵检测……
  • 5)一些摄像头、相机上的应用:夜间照片的质量改善……
  • 6)一些电影工业上是应用:换背景、电影特技……

 

什么样的人会去学(或者需要学)图像处理?

 

1)如果你是我上述那些应用领域的从业者,你当然需要掌握图像方面的理论和技术;2)相关专业的研究人员、大专院校的博士生、研究生。

所谓相关专业又是指什么呢?这个答案也可能相当宽泛,例如(但不仅限于此):Computer Science, Software Engineering, Electronic Engineering, Biomedical Engineering, Automation, Control, Applied Mathematics……

 

如何学好图像处理——我的一些箴言

 

1)对于初级入门者

 

一个扎实的基础和对于图像处理理论的完整的、系统的整体认识对于后续的深入研究和实践应用具有非常非常重要的意义。

我经常喜欢拿武侠小说《天龙八部》中的一段情节来向读者说明此中的道理,相信读者对这部曾经被多次搬上银幕的金庸作品已经耳熟能详了。书中讲到有个名叫鸠摩智的番僧一心想练就绝世武学,而且他也算是个相当勤奋的人了。但是,他错就错在太过于急功近利,甚至使用道家的小无相功来催动少林绝技。看上去威力无比,而且可以在短时间内“速成”,但实则后患无穷。最终鸠摩智走火入魔,前功尽废,方才大彻大悟。这个故事其实就告诉我们打牢基础是非常重要的,特别是要取得更长足的发展,就更是要对基本原理刨根问底,力求甚解,从而做到庖丁解牛,游刃有余。

一些看似高深的算法往往是许多基础算法的组合提升。例如,令很多人望而却步的SIFT特征构建过程中,就用到了图像金字塔、直方图、高斯滤波这些非常非常基础的内容。但是,它所涉及的基础技术显然有好几个,如果缺乏对图像处理理论的系统认识,你可能会感觉事倍功半。因为所有的地方好像都是沟沟坎坎。

关于课程——

在这个阶段其实对于数学的要求并不高,你甚至可以从一些感性的角度去形象化的理解图像处理中很多内容(但不包括频域处理方面的内容)。具体到学习的建议,如果有条件(例如你还在高校里读书)你最好能选一门图像处理方面的课程,系统地完整的地去学习一下。这显然是入门的最好办法。如此一来,在建立一个完整的、系统的认知上相当有帮助。如果你没办法在学校里上一门这样的课,网上的一些公开课也可以试试。但现在中文MOOC上还没有这方面的优质课程推荐。英文的课程则有很多,例如美国加州伦斯勒理工学院Rich教授的数字图像处理公开课——https://www.youtube.com/channel/UCaiJlKxXamoODQtlx486qJA?spfreload=10。

关于教材——

显然,只听课其实还不太够,如果能一并读一本书就最好了。其实不用参考很多书,只要一本,你能从头读到尾就很好了。如果你没有条件去上一门课,那读一本来完整的自学一下就更有必要了。这个阶段,去网上到处找博客、看帖子是不行的。因为你特别需要在这个阶段对这门学问建立一个系统的完整的知识体系。东一块、西一块的胡拼乱凑无疑是坑你自己,你的知识体系就像一个气泡,可能看起来很大,但是又脆弱的不堪一击。

现在很多学校采用冈萨雷斯的《数字图像处理》一书作为教材。这是一本非常非常经典的著作。但是我必须要提醒读者:

1)这是一本专门为Electronic Engineering专业学生所写的书。它需要有信号与系统、数字信号处理这两门课作为基础。如果你没有这两门课的基础,你读这本书要么是看热闹,要么就是看不懂。

下面是冈书中的一张插图。对于EE的学生来说,这当然不是问题。但是如果没有我说的那两门课的基础,其实你很难把握其中的精髓。H和h,一个大小一个小写,冈书中有的地方用H,有的地方用h,这都是有很深刻用意的。原作者并没有特别说明它们二者的区别,因为他已经默认你应该知道二者是不同的。事实上,它们一个表示频域信号,一个表示时域信号,这也导致有时候运算是卷积,有时候运算是乘法(当然这跟卷积定理有关)。所以我并不太建议那些没有这方面基础的学生在自学的时候读这本书。

 

2)冈萨雷斯教授的《数字图像处理》第一版是在1977年出版的,到现在已经快40年了;现在国内广泛使用的第二版是2002年出版的(第三版是2007年但是其实二者差异并不大),到现在也有20年左右的时间了。事实上,冈萨雷斯教授退休也有快30年了。所以这本书的内容已经偏于陈旧。数字图像处理这个领域的发展绝对是日新月异,突飞猛进的。特别在最近二三十年里,很多新思路,新方法不断涌现。如果你看了我前面推荐的Rich教授的公开课(这也是当前美国大学正在教学的内容),你一下子就会发现,原来我们的教育还停留在改革开放之前外国的水平上。这其实特别可怕。所以我觉得冈萨雷斯教授的《数字图像处理》作为学习过程中的一个补充还是不错的,但是如果把它作为主参考,那真的就是:国外都洋枪洋炮了,我们还在大刀长矛。

 

那么现在问题来了,对于图像处理学习者而言到底看什么书好呢?我的意见是你可以选择下面两本书中的任何一本《数字图像处理原理与实践(Matlab版)》,以及《数字图像处理:技术详解与Visual C++实践》,当然选择的标准之一就是到底你更擅长使用MATLAB还是C++。

   

 

 

 

2)对于中级水平者

 

纸上得来终觉浅,绝知此事要躬行。对于一个具有一定基础的,想更进一步的中级水平的人来说,这个阶段最重要的就是增强动手实践的能力。

还是说《天龙八部》里面的一个角色——口述武功、叹为观止的王语嫣。王语嫣的脑袋里都是武功秘籍,但问题是她从来都没练过一招一式。结果是,然并卵。所以光说不练肯定不灵啊。特别是,如果你将来想从事这个行业,结果一点代码都不会写,那几乎是不可想象的。学习阶段,最常被用来进行算法开发的工具是Matlab和OpenCV。你可以把这两个东西都理解为一个相当完善的库。当然,在工业中C++用得更多,所以Matlab的应用还是很有限的。前面我们讲到,图像处理研究内容其实包括:图像的获取和编解码,但使用Matlab和OpenCV就会掩盖这部分内容的细节。你当然永远不会知道,JPEG文件到底是如何被解码的。

如果你的应用永远都不会涉及这些话题,那么你一直用Matlab和OpenCV当然无所谓。例如你的研究领域是SIFT、SURF这种特征匹配,可以不必理会编解码方面的内容。但是如果你的研究话题是降噪或者压缩,可能你就绕不开这些内容。最开始学的时候,如果能把这部分内容也自己写写,可能会加深你的理解。以后做高级应用开发时,再调用那些库。所以具体用什么,要不要自己写,是要视你所处的阶段和自己的实际情况而定的。以我个人的经验,在我自学的时候,我就动手写了Magic House,我觉得这个过程为我奠定了一个非常夯实的基础,对于我后续的深入研究很有帮助。

 

下面这个文中,我给出了一些这方面的资源,代码多多,很值得参考学习:图像处理与机器视觉网络资源收罗

http://blog.csdn.net/baimafujinji/article/details/32332079

 

3)对于高级进阶者

 

到了这个程度的读者,编程实现之类的基本功应该不在话下。但是要往深,往高去学习、研究和开发图像处理应用,你最需要的内容就变成了数学。这个是拦在很多处于这个阶段的人面前的一大难题。如果你的专业是应用数学,当然你不会感觉有问题。但如果是其他专业背景的人就会越发感觉痛苦。

如果你的图像处理是不涉及机器学习内容的,例如用Poisson方程来做图像融合,那你就要有PDE数值解方面的知识;如果你要研究KAZE特征,你就必须要知道AOS方面的内容。如果你研究TV降噪,你又要知道泛函分析中的BV空间内容……这些词你可能很多都没听过。总的来说,这块需要的内容包括:复变函数、泛函分析、偏微分方程、变分法、数学物理方法……

如果你要涉足机器视觉方法的内容,一些机器学习和数据挖掘方法的内容就不可或缺。而这部分内容同样需要很强大的数学基础,例如最大似然方法、梯度下降法、欧拉-拉格朗日方程、最小二乘估计、凸函数与詹森不等式……

当然,走到这一步,你也已经脱胎换骨,从小白到大神啦!路漫漫其修远兮,吾将上下而求索。

 

(全文完)

 

 

图像处理的应用

阅读数 3030