初步肺部分割---记录一个星期的苦逼生活

初步肺部分割---记录一个星期的苦逼生活

肺部分割结果不是很好，望大神们指教，不慎感激

一：一开始的想法：

由于最近在学习计算机视觉，使用OpenCV开发程序，所以第一个想法是使用OpenCV里自带的一个分割函数Canny进行分割。我写了一个简单的程序，输入我是从给的肺的资料中随便挑了一个，如下：

 

然后运行canny边缘检测，通过调整canny函数中的小阈值和大阈值，得到如下结果：

二：自己的算法思想：

然后我了解到canny边缘检测是根据梯度变化检测边缘的，只要梯度出现明显变化就认为是边缘，在上图中边缘明显出来了，但是我们只要肺部的边缘，而不需要其他边缘，所以使用canny边缘检测是不行的。

然后我看图像想出了一个方法可以检测肺部的区域，具体的思想如下：

首先读入图像，然后对图像一行一行的扫描。首先从左往右扫描，如果出现了像素值明显变化，即从白变成黑，那么它一定是左肺的起始点，然后再往后扫描，如果出现像素值从黑变到白，那么这个点就是左肺的末点。此时停止扫描，转化到从右开始往左扫描，利用同前面的思想可以找到右肺的部分。

伪代码如下：

对于图像的每一行像素：

从左往右扫描此行：

Flag=0;//flag是一个标志，0表示未找到起始点，1表示找到起始点，2表示找到终止点

如果flag == 0{

如果当前像素大于某个阈值的时候（此处的阈值表示白色的范围），且这个像素的后一个像素小于某个阈值的时候（此处的阈值表示黑色的范围），那么可以断定是左肺的起始点（此行的左肺的起始点）。

Flag = 1;

}

如果flag ==1{

如果当前像素小于某个阈值的时候（此处的阈值表示黑色的范围），且这个像素的后一个像素大于某个阈值的时候（此处的阈值表示白色的范围），那么断定是左肺的终止点（此行的左肺的终止点）。

Flag = 2；

}

如果flag == 2{

Break;

从右往左扫描此行：

Flag=0;//flag是一个标志，0表示未找到起始点，1表示找到起始点，2表示找到终止点

如果flag == 0{

如果当前像素大于某个阈值的时候（此处的阈值表示白色的范围），且这个像素的前一个像素小于某个阈值的时候（此处的阈值表示黑色的范围），那么可以断定是右肺的起始点（此行的右肺的起始点）。

Flag = 1;

}

如果flag ==1{

如果当前像素小于某个阈值的时候（此处的阈值表示黑色的范围），且这个像素的前一个像素大于某个阈值的时候（此处的阈值表示白色的范围），那么断定是右肺的终止点（此行的右肺的终止点）。

Flag = 2；

}

如果flag == 2{

Break;

//算法结束

注释：

1.程序会用另一个同输入图像一样大小的图像显示结果，这个图像是黑白图像，肺内部区域使用白色标记，非肺部内部区域使用黑色标记

2.仔细观察图像，如果单纯应用上面的算法会将一处特殊的地方表示为肺内部，这处如下图红色部分所圈出来（所以要进行特殊的处理，可以使用记录左肺和右肺的起始和终止点位置，然后判断是否重合来判断是否为特殊的不符合区域）：

三：图像的位深度问题

然后我就开始写程序实现自己的想法。我一直认为我们的测试图像是灰度图像，即只有一个通道的图像。我写完程序运行，通过调整两个阈值发现，无论你怎么调整，都不会显示出肺部，即输出图像全部是黑色的。我就纳闷了好长时间，怎么会全黑的呢。最后我怀疑测试图像不是灰度图像，我查看了一下测试图像的属性，发现确实不是灰度图像（完全被表象迷惑了），测试图像的属性如下图：

注意上图中黄色部分，其显示位深度是24位，而灰度图像的位深度应该是8位，24位位深度图像说明测试图像不是灰度图像。

四：图像的色彩空间

既然不是灰度图像，且看到24位位深度，我第一反应就是这个测试图像是RGB色彩空间图像，RGB色彩空间正好使用24位来表示每一个像素。其中R、G、B分别用8位来表示。我又认为测试图像的R=G=B，这样即使是真彩色图像也会表现出灰度图像的样子（可以迷惑观察者认为是灰度图像）。然后我又写了一个程序，我只考虑R通道的灰度值。通过调整阈值发现运行结果仍然不对。然后我又测试了G和B通道，运行结果和R通道一样都是不正确的。运行结果如下图：

这个问题我也纳闷了很长时间。我开始怀疑测试图像的色彩空间不是RGB色彩空间，于是我网上搜索了一下jpg图像的介绍，发现果然jpg图像不是使用的RGB色彩空间，而是使用的YCrCb色彩空间（即YUV色彩空间）（jpg百度百科如下）：

五：色彩空间的转化

既然jpg格式的图像不是灰度图像且不是使用的RGB色彩空间，我的第一个反应就是将其转化成灰度图像。然后我写了一个小程序，将测试图像转化成灰度图像，这个可以直接使用OpenCV中的转化函数就可以实现了，转化结果如下：

上面的图像才是真正意义上的灰度图像，具体属性如下：

上面黄色区域显示此时的图像是8位位深度的图像，可以确保是灰度图像。

六：灰度图像的观察

仔细观察上面转化后的灰度图像，然后和原有的测试图像进行对比，发现转化后的灰度图像显示效果明显好于原来图像，自己原来的算法能更好的运行。

七：灰度图像的灰度值

然后我很高兴，认为基本可以搞定了，只要我的算法对转化后的灰度图像进行处理就好了。于是我又实现了一下自己的算法。通过调整阈值，发现效果不好，只有很少的阈值对效果稍微好点，而且这些阈值对中的两个阈值很靠近。我就又纳闷了，转化后的图像黑白很分明，为什么要设置很靠近的阈值才能有稍微好一点的显示效果呢？难道我看上去是黑的不是那么黑，看上去是白的不是那么白？

于是我又写了一个小程序，将转化后的灰度图像的每一个像素值都输出到文件中，然后我去查看灰度图像的灰度值，部分灰度值截图如下（我是图像的每一行数据输出为一行，但是测试图像是1200*1116大小的，即每行有1116个像素值，虽然我程序是将这1116个数据显示在一行中的，但是当使用记事本打开文件时这1116个数据是不会显示在一行中的，原因是记事本一行显示的个数是有限的，所以这1116个数据在记事本中要分成多行显示）：

                                                                        第一行的数据

从第一行的数据可以看出，虽然第一行看上去是黑的，但是实际上它的灰度值并不低，为79。

我从文件中找到特殊部分需要另处理的一行数据如下（黄色区域表示起始和终止部分）：

                                                                      需要另处理的一行数据

从上图中可以发现两点：第一，灰度图像中的白色的像素值大约是178。第二，灰度图像的的从白到黑以及从黑到白之间的变化不是想象中的那么明显（至少在灰度值是这样的，感觉上变化明显但是灰度值上却不是），注意到上图黄色标注的地方，它们就是变化的地方，可以看出它们变化不是很剧烈，上图中显示的变化区域还算可以，继续观察灰度图像的数据可以发现，有些行的变化区域的相邻两个间距是5。

从灰度图像的灰度值可以看书，如果程序中两个阈值设置的过大，那么会将本来是存在肺内部的一行判断成没有肺内部，即全部是黑色的。

八：最后的运行结果

九：接下来的工作

1、程序的运行结果还是不够理想，我仔细看了一下灰度图像的每行的起始部分灰度值，发现都是79，但是从输出可以看书存在大于85，小于100的像素值（开始部分是白色表示存在这样的灰度值），这个不知道是什么回事？至少我觉得每行的一开始不会有错误的，但是为什么有些行一开始就出现错误，而有些就没有，这个研究了很长时间还是不知道错在哪里。

2、初始图像在转化成RGB图像的过程中像素值我觉得都向中间靠拢，使得黑白像素差距没想象中的那么大，是否可以直接对原始的dicom格式操作，需要了解一下dicom格式，这样是否黑白像素值差距会大一点

3、特殊情况（如上面所述），在程序中我已经考虑了，可是仍然输出，表示很奇怪呀

4、使用其他的测试图像跑一下程序，看看这里设置的85、100这两个阈值是不是适合其他的测试图像

十：源代码

#include<iostream>

#include<fstream>

#include<string>

#include<opencv\highgui.h>

#include<opencv\cv.h>

using namespace std;

int main()

{

IplImage * orignalImage;

IplImage * grayImage;

IplImage * resultImage;

IplImage * tmpImage;

uchar * ptr_gray;

uchar * ptr_result;

int position1;

int position2;

int position3;

int position4;

int flag1;

int flag2;

orignalImage = cvLoadImage("D:\\opencv program\\fengge5\\1.jpg");

cvNamedWindow("Orignal Image",CV_WINDOW_AUTOSIZE);

cvNamedWindow("Gray Image",CV_WINDOW_AUTOSIZE);

cvNamedWindow("Result Image",CV_WINDOW_AUTOSIZE);

tmpImage = cvCreateImage(cvSize(orignalImage->width,orignalImage->height),IPL_DEPTH_8U,3);

grayImage = cvCreateImage(cvSize(orignalImage->width,orignalImage->height),IPL_DEPTH_8U,1);

resultImage = cvCreateImage(cvSize(orignalImage->width,orignalImage->height),IPL_DEPTH_8U,1);

cvCvtColor(orignalImage,tmpImage,CV_YCrCb2RGB);

cvCvtColor(tmpImage,grayImage,CV_RGB2GRAY);

for(int i = 0 ; i < grayImage->height ; ++i)

{

ptr_gray = (uchar *)(grayImage->imageData + i * grayImage->widthStep);

ptr_result = (uchar *)(resultImage->imageData + i * resultImage->widthStep);

flag1 = 0;

for(int j = 0 ; j < grayImage->width - 1; ++j)

{

if(flag1 == 0)

{

if(*(ptr_gray + j) > 100 && *(ptr_gray + j + 1) < 85)

{

*ptr_result = 255;

position1 = j;

flag1 = 1;

}

else

{

*ptr_result = 0;

}

}

else

{

if(*(ptr_gray + j) < 85 && *(ptr_gray + j + 1) > 100)

{

*ptr_result = 255;

position2 = j;

break;

}

else

{

*ptr_result = 255;

}

}

++ptr_result;

}

if(flag1 == 1)

{

ptr_gray = (uchar *)(grayImage->imageData + i * grayImage->widthStep + grayImage->width - 1);

ptr_result = (uchar *)(resultImage->imageData + i * resultImage->widthStep + resultImage->width - 1);

flag2 = 0;

for(int j = 0 ; j < orignalImage->width - 1 ; ++j)

{

if(flag2 == 0)

{

if(*(ptr_gray - j) > 100 && *(ptr_gray - j - 1) < 85)

{

*ptr_result = 255;

flag2 = 1;

position3 = orignalImage->width - j - 1;

}

else

{

*ptr_result = 0;

}

}

else

{

if(*(ptr_gray - j) < 85  && *(ptr_gray - j - 1) > 100)

{

*ptr_result = 255;

position4 = orignalImage->width - j - 1;

if(!(position3 == position2 && position4 == position1))

{

ptr_result = (uchar *)(resultImage->imageData + i * resultImage->widthStep);

for(int k = position2 ; k < position4 ; ++k)

{

*(ptr_result + k) = 0;

}

}

break;

}

else

{

*ptr_result = 255;

}

}

--ptr_result;

}

//--ptr_result;

}

else

{

*ptr_result = 0;

}

}

/*

string filename;

filename = "data.txt";

ofstream fout(filename.c_str());

for(int i = 0 ; i < resultImage->height ; ++i)

{

uchar * ptr = (uchar *)(resultImage->imageData + i * resultImage->widthStep);

for(int j = 0 ; j < resultImage->width ; ++j)

{

fout<<(int)(*(ptr+j))<<" ";

}

fout<<endl;

}

*/

cvSaveImage("3.jpg",resultImage);

cvShowImage("Orignal Image",orignalImage);

cvShowImage("Gray Image",grayImage);

cvShowImage("Result Image",resultImage);

cvWaitKey(0);

cvReleaseImage(&orignalImage);

cvReleaseImage(&grayImage);

cvReleaseImage(&resultImage);

cvDestroyWindow("Orignal Image");

cvDestroyWindow("Gray Image");

cvDestroyWindow("Result Image");

return 0;

}

十一：总结

找出问题所在实在是痛苦的事情，这周我终于领悟到一个方法：可以将问题的每个前提条件或者假设先列出来，然后一个一个进行怀疑，逐个去验证是否真的正确