乳腺癌图像数据预处理--mias-database1

深度学习小白，最近在做乳腺癌良恶性分类的内容，试着复现论文：
乳腺钼靶图像中的肿块良恶性分类 - 中国知网

这里面的图像预处理花了不少时间（因为是小白嘛）

而且我用的数据集是mias-database1和论文中的不一样，所以后面效果方面和论文相比会有出入

相关数据集下载参看：DDSM database、INbreast database、MIAS等乳腺MG数据获取方式

（事先说好哈，水平有限，有问题那那那…就是水平有限，勿怪勿怪）

大致流程如下：

先采用OTSU（大津法）方法将乳腺部分和背景区域分割开，之后再使用截断归一化对乳房的最小外接矩形区域做归一化处理，之后采用两种不同的阈值对矩形区域进行自适应直方图均衡化操作，最后将以上三幅图像拼接成最终的三通道增强图像。
​​
具体流程：

1.先将图像高斯滤波一次，再用大津法处理，得到一张二值图。

可以看到上图除了有乳房图像之外还有医生铭牌信息，而且背景占比太多。

大津法实现可参看：
CV笔记5：图像分割之最大类间方差法、自适应阈值分割(基于python-opencv实现)_my_kun的博客-CSDN博客

2.计算二值图中的连通区域，选取最大的一个（除背景外，最大的区域是乳腺），其他的舍去。

然后进行裁剪

其中第一张是原图，第二张是根据二值图计算连通区域，并保存乳房区域舍去其他区域之后的结果

第三张图是用opencv计算图像最小外接矩形之后剪切得到。第一二步代码如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 这个后面用于对得到的坐标进行排序，因为得到的坐标
# 有的时候第三个表示左下角有的时候是右上角，排序之后才好进行计算
def judge_box(box_):
    box = [[0,0],[0,0],[0,0],[0,0]]
    for i in range(4):
        for j in range(2):
           box[i][j] = box_[i][j]
    a = [0,0,0,0]
    b = [sum(box[0]), sum(box[1]), sum(box[2]), sum(box[3])]
    max_ = b.index(max(b))
    min_ = b.index(min(b))
    a[3] = box[max_]
    a[0] = box[min_]
    if max_ >= min_:
        box.pop(min_)
        box.pop(max_-1)
    else:
        box.pop(max_)
        box.pop(min_-1)
    if box[0][0] >= box[1][0]:
        a[1] = box[0]
        a[2] = box[1]
    else:
        a[1] = box[1]
        a[2] = box[0]
    return a

# 以下为主函数（mias-database1数据集中一共322张图像）
for n in range(1,323):
    # 以下用于得到只含一个区域的图像
    src = cv2.imread('./mias-database1/mdb {}.pgm'.format(n),0)  # 这个0不要忘了
    blur = cv2.GaussianBlur(src, (5, 5), 0)  # 高斯滤波
    ret3, th3 = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)  # 大津法 （得到的th3是二值图）
    thu = cv2.connectedComponentsWithStats(th3)  # 得到连通区域数据,thu是一个数组，含四个数据，第二个是全图的label值
    list1 = np.reshape(thu[1], (thu[1].shape[0]*thu[1].shape[1]))
    list2 = np.bincount(list1)
    label1 = np.argmax(list2[1:]) + 1   # 得到除背景外面积最大的区域的label
    dst = src.copy()
    # 根据得到的标签将图中的其他如医生铭牌图像去掉，只保留乳房图像
    for i in range(src.shape[0]):
        for j in range(src.shape[1]):
            if thu[1][i,j] != label1:
                dst[i,j] = 0
    
    # 得到新的二值图
    th4 = np.zeros_like(dst)
    for i in range(th4.shape[0]):
        for j in range(th4.shape[1]):
            if dst[i,j] != 0:
                th4[i,j] = 1

    (cnts, _) = cv2.findContours(th4, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)  # 得到物体轮廓信息
    c = sorted(cnts, key = cv2.contourArea, reverse=True)[0]
    rect = cv2.minAreaRect(c)         # 得到物体轮廓最小外接矩形的四点坐标
    box = np.int0(cv2.boxPoints(rect))  # 处理得到的数据
    box = judge_box(box)  # 将坐标进行排序

    # 下面操作是为了得到较为精确的图片，因为得到的物体轮廓可能并不和图片边界平行，直接裁剪会遗漏信息。
    points1 = np.array([box[0], box[1], box[2], box[3]], dtype='float32') 
    row = ((box[0][0]-box[1][0])**2+(box[0][1]-box[1][1])**2)**0.5
    col = ((box[0][0]-box[3][0])**2+(box[0][1]-box[3][1])**2)**0.5          # 根据得到的矩形坐标，计算矩形的长宽
    points2 = np.array([[0, 0], [row, 0], [0, col], [row, col]], dtype='float32') # 根据原有矩阵得到一个新的矩阵
    
    M = cv2.getPerspectiveTransform(points1, points2)  # 得到新矩阵和旧矩阵之间的映射关系，为后面的函数提供参数

    outpic = cv2.warpPerspective(dst, M, (int(row), int(col)))  # 根据映射关系，将原图经过变换后得到新的图片，数据类型一定要注意
    plt.subplot(1,3,1), plt.imshow(src, cmap='gray')   # 原图
    plt.subplot(1,3,2), plt.imshow(dst, cmap='gray')   # 只保留乳房信息
    plt.subplot(1,3,3), plt.imshow(outpic, cmap='gray') # 裁剪后的图片
    plt.show()
    
    cv2.imwrite('./min_box_image2/{}.pgm'.format(n), outpic)

连通区域计算可参看：cv2.connectedComponentsWithStats()_zhangjipinggom的博客-CSDN博客

得到物体轮廓的最小外接矩形可参看：使用OpenCV和Python检测图像中的物体并裁剪下来_余小晴的博客-CSDN博客

后续的坐标排序和变换矩阵可参看（只需要看实现思路和实现用到的函数，尽管所用计算机语言不一样但函数也差不多，自己多去试几次就可以发现了）：利用OpenCV提取图像中的矩形区域（PPT屏幕等）_weixin_34023982的博客-CSDN博客

具体的其中用到的函数cv2.getPerspectiveTransform和函数cv2.warpPerspective可参看：使用python opencv透视变换getperspectiveTransform获得图像的鸟瞰图_Abc2149的博客-CSDN博客

3.对得到的图像进行截断归一化处理：


其中第一张图像是第二步得到的结果图，第二张图像是截断归一化后的。后面两张图是前两张图像的灰度分布 。

公式如下：

代码如下：

for i in range(1,323):
    img = cv2.imread('./min_box_image2/{}.pgm'.format(i), 0)  # 读取图片

    list_ = []
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if img[i,j] != 0:
                list_.append(img[i,j])   # 将数据中的0全部删去

    list_.sort()         # 将数据按照大小进行排序
    max_ = list_[-int(len(list_)*0.01)]  # 得到靠近大端1%位置的像素值--作为最大像素
    min_ = list_[int(len(list_)*0.05)]   # 得到靠近小端5%位置的像素值--作为最小像素
    # print(max_, min_)

    img2 = np.zeros_like(img, np.float64)
    # 根据得到的最大最小像素值对图像进行处理
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            if img[i,j] != 0:
                if img[i,j] <= min_:
                    img2[i,j] = min_
                elif img[i,j] >= max_:
                    img2[i,j] = max_
                else:
                    img2[i,j] = img[i,j]
    
    # 再次处理图像
    for i in range(img2.shape[0]):
        for j in range(img2.shape[1]):
            if img2[i,j] != 0:
                img2[i,j] = (float(img2[i,j]) - float(min_))/(float(max_) - float(min_))
    
    # plt.subplot(1,2,1), plt.imshow(img, cmap='gray')
    # plt.subplot(1,2,2), plt.imshow(img2, cmap='gray')
    # plt.show()

    # plt.subplot(1,2,1), plt.hist(img.ravel(), 255)
    # plt.subplot(1,2,2), plt.hist(img2.ravel(), 255)
    # plt.show()

    img2 = np.uint8(img2*255.)   # 转换数据格式，CV2保存灰度图要uint8数据类型
    cv2.imwrite('./normalization/{}.pgm'.format(n), img2)

4.限制对比度自适应直方图均衡(CLAHE)：

第一张是截断归一化处理后的图像，第二张是阈值为5.0的CLAHE处理后的，第三张是阈值为10.0的（具体的参数论文中没有给，只能自己试着来）

后三张是前三张图像的灰度分布

代码：

for i in range(1, 323):
    src = cv2.imread('./normalization/{}.pgm'.format(i), 0)
    clahe = cv2.createCLAHE(clipLimit=5.0,tileGridSize=(8,8))  # 实例化均衡直方图参数
                                                               # 第一个参数表示阈值大小， 
                                                               # 第二个表示像素均衡化网格大小
    dst = clahe.apply(src)        # 进行均衡化操作
    cv2.imwrite('./limit/{}.pgm'.format(i), dst)

cv2.createCLAHE(clipLimit=5.0,tileGridSize=(8,8)) 函数详情见链接：机器学习进阶-直方图与傅里叶变化-直方图均衡化 1.cv2.equalizeHist(进行直方图均衡化) 2. cv2.createCLAHA(用于生成自适应均衡化图像) - python我的最爱 - 博客园

接下来换一个阈值（10.0），即更换上面函数中的clipLimit参数，得到新的图像（代码同上）（关于阈值的后面总结部分还有点内容，大致就是…我这阈值可能应该大概准确来说不太合理）

5.将前面得到的截断归一化后的图像和经过两种不同阈值的CLAHE处理得到的图像组成三通道图像：

将三张图拼接很简单，但是不知道为什么画出的图却没有论文中的效果（差距还挺大的，不过后面看用python的matplotlib.pyplot画出的图和计算机自带软件展示的图有差别，实际的图会好一点）：

两张原图不一样的，因为我的数据集和论文里的不一样，大家看个大概就行。

展示的目的就是想说我直接三通道处理（即三张图拼接）然后展示，就展示不出那种效果（虽然计算机自带软件展示好了不少）。我觉得可能是数据集不一样的问题（这个感觉概率不大，因为都是差不多的灰度图），也可能是其他问题，不管了，大致特征还在就行。

还有一个小小的问题：

本来第二步开始我没有计算连通区域然后去掉铭牌信息的，因为当时看不懂那个计算连通区域的算法，后面也没找到已有的函数（之后在导师帮助下找到了）。

然后我就直接将大津法得到的二值图输入cv2.findContours函数中找到物体轮廓的最小外接矩形，后面一对比，发现这个和计算连通区域之后的效果相差不大，虽然还是有铭牌信息，但是黑色背景什么的都去除的差不多。

（这个是对比图，第二张图没有进行其他处理，只是在第一张图的基础上进行了变换）

但是一些其他的可能因为铭牌信息所以计算最小外接矩阵的时候矩阵稍微大了点，然后得到的图黑色背景多了点

（虽然感觉效果不差，但咱还是不省那点代码，总共也不用多久）

------------------------分割线-------------------------
以下内容与图像预处理无关，而且为个人浅薄观点，如有问题，就当我是瞎写的（本人水平本就不高（′Д`））

后续我将处理好的图用VGG16网络训练后效果不好，过拟合，验证精度没有变化，应该是没学到特征。

我总结了下（就我直观感受）：
第一，应该是我对这个数据集分类不恰当，因为mias-database1数据集中包含多类，比如正常、钙化、肿块等，而在那些有问题的图像中又分了良性和恶性。而我只是简单地将正常的和所有的良性分为一类，所有的恶性单独为一类，简化为二分类问题。这是很不对的，因为论文中只是针对肿块的良恶性进行分类，目的性强，而且他的制作的数据集只包含这两类（且每张都有肿块，应该是每张都有）。不像我这，什么类型的数据都有，感觉很散，分类又不合理，自然很难学到好的特征。（而且在我经过翻转对称操作之后数据量增加了三倍，但还是没用，哈哈还是分类不够合理）
第二，前面提到的两种不同阈值的CLAHE处理，感觉这个阈值设置的不够好，一些图片高亮部分连成一片，啥细节都没了，好失败 o(≧口≦)o
见图：
第三，网络用的可能不对（VGG16不是论文中用到的，但是这篇博客是用来介绍我复现的图像预处理过程的，训练什么的…小白一枚，慢慢摸索吧━┳━　━┳━）