边缘检测、阈值处理和区域生长

1、医学图像处理实验报告实验十三：边缘检测、阈值处理和区域生长日期： 2014年05月27日摘要本次实验的目的是： 了解边缘检测原理，用梯度阈值法，使用Sobel算子结合平滑处理和阈值处理提取边缘，； 了解阈值处理的计算方法，进行全局阈值和Otsu（大律法）阈值处理。 了解区域生长原理。对图像做区域生长提取图像特征。本次实验的内容是： 阈值处理 全局阈值处理 OTSU（大律法）阈值处理 梯度法检测边缘 区域生长一、技术讨论1.1实验原理1.1.1图像的边缘检测边缘：是指图像局部特征的不连续性。灰度或结构信息的突变。边缘检测：一种定位二维或三维图像（特别是医学图像）中的对象的边缘的系统。通过输入端（

2、310）接收表示该图像的各元素值的数据元素集。该数据集被存储在存储装置（320）中。处理器（340）确定该图像中的对象的边缘。该处理器计算所述数据元素的至少一阶和/或二阶导数，并且计算该图像的等照度线曲率，所述曲率由标识。该处理器还确定校正因数，该校正因数对于由对象的曲率和/或所述数据的模糊造成的边缘错位进行校正。该校正因数取决于所述等照度线曲率。然后，该处理器确定取决于所计算出的导数和所述等照度线曲率的算子的过零点。该系统的输出端（330）提供对于该图像中的边缘位置的指示。1.1.2图像的Sobel梯度算子主要用作边缘检测。在技术上，它是一离散性差分算子，用来运算图像亮度函数的梯度之近似值。

3、该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像，其公式如下:图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小。公式如下：计算梯度方向公式如下：在以上例子中，如果以上的角度等于零，即代表图像该处拥有纵向边缘，左方较右方暗。1.1.3全局阈值法全局阈值是指整幅图像使用同一个阈值做分割处理，适用于背景和前景有明显对比的图像。它是根据整幅图像确定的：T=T(f)。但是这种方法只考虑像素本身的灰度值，一般不考虑空间特征，因而对噪声很敏感。常用的全局阈值

4、选取方法有利用图像灰度直方图的峰谷法、最小误差法、最大类间方差法、最大熵自动阈值法以及其它一些方法。1.1.4最佳阈值法（Otsu算法）该方法是根据图像的灰度特征，将图像划分为目标和背景两种类型。目标类和背景类之间的类间方差越大，说明图像的目标和背景之间的差别越大。当背景类信息被错误的分为目标类信息，或目标类信息被错误的分为背景类信息时，都会导致目标和背景之间的差别变小。正是由于这一特性，只要使类间方差最大或类内方差最小，这样的分割就意味着最小的误差，并且根据该准则选取最佳阈值使得目标和背景之间的分离性最好。1.1.5区域生长区域生长的基本思想是将具有相似性质的像素集合起来构成区域。具体先对每

5、个需要分割的区域找一个种子像素作为生长的起点，然后将种子像素周围邻域中与种子像素有相同或相似性质的像素(根据某种事先确定的生长或相似准则来判定)合并到种子像素所在的区域中。将这些新像素当作新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来。这样一个区域就长成了。1.2实验函数cvThreshold阈值处理.void cvThreshold( const CvArr* src, CvArr* dst,double threshold, double max_value, int threshold_type )src：输入图像；dst：输出图像；threshold：阈值；cvC

6、reateTrackbar创建trackbar并将它添加到指定的窗口。CV_EXTERN_C_FUNCPTR( void (*CvTrackbarCallback)(int pos) );int cvCreateTrackbar( const char* trackbar_name, const char* window_name, int* value, int count, CvTrackbarCallback on_change );trackbar_name被创建的trackbar名字。window_name窗口名字value整数指针，它的值将反映滑块的位置。这个变量指定创建时的滑块位

7、置。count滑块位置的最大值。最小值一直是0。on_change每次滑块位置被改变的时候，被调用函数的指针 二、结果与讨论2.1实验结果（运行每个程序，根据以下实验要求得到结果）1. 实验要求（与实验指导的实验要求相同）1、10-02（2）：可以随时自选阈值，大家可以先用此程序尝试出最佳阈值 （a） （b）图2.1 程序10-02（2）结果，图为大致最佳阈值2、10-02：本程序使用迭代法，比书本上的程序更加精确(完全使用书本上的公式做的，但不知道为什么就是出来这个结果)，但还不如Otsu，能直接输出阈值，因为是直接计算，因此没有参数要求，大家可以多换几幅图片试试。 （a）histogram

8、 （b）original pic （c）result图2.2 10-02迭代法实验结果4、10-01：Sobel算子检验边缘，用blur()函数做了一个平滑处理，还用了一个阈值处理，因为选择的图像特征比较复杂，因此可以看到提取效果不是很好，自己找出平滑的模板大小为多少时，对应阈值为多少时效果最好（尽量保留血管部分，除去其他部分）。由于程序会保存图像，因此可以用10-03的程序找出最佳阈值。大家先熟悉10-03，找起来会很快Ø 先利用10-03求出最佳阈值，输出结果为：125 init done.Ø 设置不同大小的平滑模板，并用不同的最佳阈值处理，结果如下：1）模板为5x5大

9、小时，55 init done （a）original pic （b）gradient （c）result2）模板为8x8大小时，30 init done （d）original pic （e）gradient （f）result3） 模板为15x15大小时，27 init done （g）original pic （h）gradient （i）result4）模板为20x20大小时，25 init done （j）original pic （k）gradient （l）result图2.3 Sobel算子检验边缘实验结果5、10-04：区域生长，目的是从图片中提取出两个大圈的轮廓，程序中从阈

10、值化图像中提取出种子（可以改变阈值因子S，在195行），并使用种子做区域生长（可以改变阈值因子Q，在225行），得到最佳边界。（a）取不同的S和Q处理结果如下：1）S=0.95,Q=35 (a) 种子 (b)阈值处理 (c)区域生长的结果2）S=0.9,Q=35 （d）种子 (e)阈值处理 (f)区域生长的结果 3) S=0.85,Q=35 （g）种子 (h)阈值处理 (i)区域生长的结果4) S=0.95,Q=25 （j）种子 (k)阈值处理 (l)区域生长的结果图2.4 区域生长实验结果Ø 2.2实验讨论开始10-03：Otsu阈值处理，也叫大律法阈值处理，要求写出程序流程图读入

11、图像遍历每个像素得到整幅图像的均值、前景均值、背景均值计算最佳阈值利用阈值处理图像输出结果像结束10-01：Sobel算子检验边缘，用blur()函数做了一个平滑处理，还用了一个阈值处理，因为选择的图像特征比较复杂，因此可以看到提取效果不是很好，自己找出平滑的模板大小为多少时，对应阈值为多少时效果最好（尽量保留血管部分，除去其他部分）。由于程序会保存图像，因此可以用10-03的程序找出最佳阈值。大家先熟悉10-03，找起来会很快，讨论一下平滑模板和阈值分别对结果有什么影响由图2.3结果可知，平滑模板的大小影响处理后图片边缘的宽窄，阈值的大小影响处理后图像的清晰度。模版越大，处理后的图像越模糊，

12、阈值处理后得到细节越少，并且保留的线条越粗。而模版大小越小，处理后的图像模糊度越小，阈值处理后得到的细节越多，不过同时得到提取信息的效果就没有那么好。由图2.4结果，在控制Q不变时改变S，S越大，处理后的图像的边界越小，S越小，处理后的图像的边界越大。在控制S不变时改变Q，S越小，处理后的图像的边界越小，S越大，处理后的图像的边界越大。附录（实验代码）.pro程序：INCLUDEPATH+=C:Qtopencv2.2includeopencvC:Qtopencv2.2includeopencv2C:Qtopencv2.2includeLIBS+=C:Qtopencv2.2liblibcv.dl

13、l.aC:Qtopencv2.2liblibopencv_calib3d220.dll.aC:Qtopencv2.2liblibopencv_contrib220.dll.aC:Qtopencv2.2liblibopencv_core220.dll.aC:Qtopencv2.2liblibopencv_features2d220.dll.aC:Qtopencv2.2liblibopencv_flann220.dll.aC:Qtopencv2.2liblibopencv_gpu220.dll.aC:Qtopencv2.2liblibopencv_highgui220.dll.aC:Qtopenc

14、v2.2liblibopencv_imgproc220.dll.aC:Qtopencv2.2liblibopencv_legacy220.dll.aC:Qtopencv2.2liblibopencv_ml220.dll.aC:Qtopencv2.2liblibopencv_objdetect220.dll.aC:Qtopencv2.2liblibopencv_video220.dll.a2. 阈值处理：#include <stdio.h>#include <stdlib.h>#include <math.h>#include "cv.h"

15、#include "highgui.h"int main(void) IplImage *src_image = 0; IplImage *dst_image = 0; IplImage *dst_image_adaptive = 0; int c; int adaptive_method; int block_size; int offset; if( (src_image = cvLoadImage( "D:/3.tif", 0) = 0 ) return -1; dst_image = cvCreateImage(cvSize(src_image-

16、>width,src_image->height), IPL_DEPTH_8U, 1); cvNamedWindow( "Threshold", 1 ); dst_image_adaptive = cvCreateImage(cvSize(src_image->width,src_image->height), IPL_DEPTH_8U, 1); cvNamedWindow( "ThresholdAdaptive", 1); int thresh_val = 100; int thresh_valBef = 0; cvCreate

17、Trackbar( "thresh val", "Threshold", &thresh_val, 255, 0 ); adaptive_method = CV_ADAPTIVE_THRESH_MEAN_C; int adaptive_methodBef = 1; cvCreateTrackbar( "adaptive method", "ThresholdAdaptive", &adaptive_method, 1, 0 ); block_size = 1; int block_sizeBef =

18、 0; cvCreateTrackbar( "block size", "ThresholdAdaptive", &block_size, 60, 0 ); offset = 30; int offsetBef = 0; cvCreateTrackbar( "offset", "ThresholdAdaptive", &offset, 60, 0 ); for(;) if(thresh_valBef != thresh_val ) | (adaptive_methodBef != adaptive_

19、method ) | (block_sizeBef != block_size) | (offsetBef != offset) cvAdaptiveThreshold(src_image, dst_image_adaptive, (double)255, adaptive_method, CV_THRESH_BINARY, block_size*2+3, (double)(offset-30); cvShowImage("ThresholdAdaptive", dst_image_adaptive); cvThreshold(src_image, dst_image, (

20、double)thresh_val, (double)255, CV_THRESH_BINARY); cvShowImage("Threshold", dst_image); thresh_valBef = thresh_val; adaptive_methodBef = adaptive_method; block_sizeBef = block_size; offsetBef = offset; c = cvWaitKey(30); if( c = 'q' | c = 'Q' | (c & 255) = 27 | (c &

21、 255) = 32 ) break; cvDestroyWindow( "Threshold" );/销毁窗口 cvDestroyWindow( "ThresholdAdaptive" );/销毁窗口 cvReleaseImage( &src_image ); /释放图像 cvReleaseImage( &dst_image ); /释放图像 cvReleaseImage( &dst_image_adaptive ); /释放图像3. 全局阈值处理#include "cv.h"#include "c

22、xcore.h"#include "highgui.h"#include <iostream>using namespace cv;using namespace std;int hei;int wid;int BasicGlobalThreshold(IplImage* A) long N = hei * wid; int h256; for(int i = 0; i < 256; i+) hi = 0; uchar min=255,max=0; for(int i = 0; i < hei; i+) for(int j = 0; j &

23、lt; wid; j+) int k=(uchar*)(A->imageData + A->widthStep*i)j; if (k<min) min=k; if (k>max) max=k; hk+; int i,t,t1,t2,k1,k2; double u,u1,u2; t=0; u=0; for (i=min;i<max;i+) t+=hi; u+=i*hi; k2=(int) (u/t); / 计算此范围灰度的平均值 do k1=k2; t1=0; u1=0; for (i=min;i<=k1;i+) / 计算低灰度组的累加和 t1+=hi; u1

24、+=i*hi; t2=t-t1; u2=u-u1; if (t1) u1=u1/t1; / 计算低灰度组的平均值 else u1=0; if (t2) u2=u2/t2; / 计算高灰度组的平均值 else u2=0; k2=(int) (u1+u2)/2); / 得到新的阈值估计值 while(abs(k1-k2)>0.001); / 数据未稳定，继续 return k2;/直方图均衡变换函数void Hist_image(IplImage *target_image) float range = 0,255; /直方图的灰度级的范围（以8bit图为例）,range是一个数组 floa

25、t* ranges=range; /凡是变量前加了星号，都代表定义该变量是指针变量，就像ranges /不妨先定义我们将要画的直方图的尺寸，256×256 int hist_width=512; int hist_height=512; /原图转换为灰度图 IplImage* gray_image = cvCreateImage(cvGetSize(target_image),8,1); cvZero(gray_image); cvCvtColor(target_image,gray_image,CV_BGR2GRAY); /获取灰度图一维直方图的数据，统计图像在0 255像素值的数

26、据情况，也就是说灰度为0的像素有多少个 /，像素为1的像素点有多少个。如此类推，这些数据都将会存在gray_hist里面 CvHistogram* gray_hist = cvCreateHist(1,&hist_width,CV_HIST_ARRAY,ranges,1); cvCalcHist(&gray_image,gray_hist,0,0); /计算灰度图像的一维直方图,修改*gray_image可以获得原来图像的直方图 /归一化直方图 cvNormalizeHist(gray_hist,1.0); int scale = 1; /创建一张一维直方图的“图”，横坐标设为

27、灰度级别，纵坐标为像素个数（*scale） IplImage* hist_image = cvCreateImage(cvSize(hist_width*scale,hist_height),8,3); cvZero(hist_image); /统计直方图中的最大直方块 float max_value = 0; cvGetMinMaxHistValue(gray_hist, 0,&max_value,0,0); /分别将每个直方块的值绘制到hist_image中 float bin_val; for(int i=0;i<hist_width;i+) bin_val = cvQue

28、ryHistValue_1D(gray_hist,i); /这里可以把bin_val理解为灰度级为i的概率 int intensity = cvRound(bin_val*hist_height/max_value); /要绘制的高度 cvRectangle( hist_image, cvPoint(i*scale,hist_height-1), cvPoint(i+1)*scale - 1, hist_height-intensity), CV_RGB(255,255,0); cvNamedWindow("Histogram"); cvShowImage("Hi

29、stogram",hist_image);int main () IplImage* img = cvLoadImage("D:/2.tif",1); IplImage* img0 = cvCreateImage(cvGetSize(img),img->depth,img->nChannels); hei = img->height; wid = img->width; cvCopy(img,img0); CvScalar s; float k; k=BasicGlobalThreshold(img); cout<<k; int

30、 sum=0; for (int i=0;i<img->height;i+) for (int j=0;j<img->width;j+) s = cvGet2D(img,i,j); sum = (s.val0+s.val1+s.val2)/3; if (sum > k)/阈值 s.val0=s.val1=s.val2=255; cvSet2D(img,i,j,s); else s.val0=s.val1=s.val2=0; cvSet2D(img,i,j,s); Hist_image(img0); cvNamedWindow("Origin image&

31、quot;); cvShowImage("Origin image",img0); cvNamedWindow("Result"); cvShowImage("Result",img); cvWaitKey(0);4. 最佳阈值处理#include <iostream>#include <math.h>#include <cv.h>#include <highgui.h>#include <stdio.h>#include <windows.h>using nam

32、espace std;int hei;int wid;void otsu(IplImage* A, IplImage* B)long N = hei * wid;int h256;double p256,u256,w256;for(int i = 0; i < 256; i+) hi = 0; pi = 0; ui = 0; wi = 0;for(int i = 0; i < hei; i+) for(int j = 0; j < wid; j+) for(int k = 0; k < 256; k+) if(uchar*)(A->imageData + A-&g

33、t;widthStep*i)j = k) hk+; for( int i = 0; i < 256; i+) pi = hi / double(N);int T = 0;double uT,thegma2fang;double thegma2fang_max = -10000;for(int k = 0; k < 256; k+) uT = 0; for(int i = 0; i <= k; i+) uk += i*pi; wk += pi; for(int i = 0; i < 256; i+) uT += i*pi; thegma2fang = (uT*wk - u

34、k)*(uT*wk - uk) / (wk*(1-wk); if(thegma2fang > thegma2fang_max) thegma2fang_max = thegma2fang; T = k; cout<<T;/显示阈值for(int i = 0; i < hei; i+) for(int j = 0; j < wid; j+) if(uchar*)(A->imageData + A->widthStep*i)j > T) (uchar*)(B->imageData + B->widthStep*i)j = 255; els

35、e (uchar*)(B->imageData + B->widthStep*i)j = 0; int main(int argc, char* argv) IplImage* InPut = cvLoadImage( "D:/3.tif", 0 ); IplImage *OutPut = cvCreateImage(cvSize(InPut->width,InPut->height),IPL_DEPTH_8U,1); hei = InPut->height; wid = InPut->width; /计算阈值 otsu(InPut,Ou

36、tPut); cvNamedWindow( "Resource", 1 ); cvShowImage( "Resource", InPut ); cvNamedWindow( "Result", 1 ); cvShowImage( "Result", OutPut ); cvWaitKey(0);5. 梯度法检测边缘#include "cv.h"#include "cxcore.h"#include <iostream>#include <opencv2/c

37、ore/core.hpp>#include <opencv2/imgproc/imgproc.hpp>#include <opencv2/highgui/highgui.hpp>using namespace cv;using namespace std;int main(int argc, const char * argv) / 载入图像 Mat source_image = imread("D:/1.tif", CV_LOAD_IMAGE_GRAYSCALE), img,imgx,imgx2,imgy,imgy2,Mx; imshow(

38、"orgin img", source_image ); / 建立Sobel算子模板，gx,因为负值会变为0，因此还需要一个从下到上的 Mat gx = (Mat_<char>(3,3) << -1, -2, -1, 0, 0, 0, 1, 2, 1); Mat gx2 = (Mat_<char>(3,3) << 1, 2, 1, 0, 0, 0, -1, -2, -1); / 建立Sobel算子模板，gy Mat gy = (Mat_<char>(3,3) << -1, 0, 1, -2, 0, 2,

39、 -1, 0, 1); Mat gy2 = (Mat_<char>(3,3) << 1, 0, -1, 2, 0, -2, 1, 0, -1); /求偏导数 blur(source_image, img,Size(10,10);/用10乘10的模板模糊化 filter2D( img, imgx, -1, gx ); filter2D( img, imgy, img.depth(), gy ); filter2D( img, imgx2, -1, gx2 ); filter2D( img, imgy2, img.depth(), gy2 ); Mx=imgx+imgy+i

40、mgx2+imgy2;/ imshow("gx", imgx2);/ imshow("gy", imgy2); imshow("gradient", Mx); imwrite("D:/11.tif",Mx);/保存 /读取 IplImage* img0 = cvLoadImage("D:/11.tif",1); CvScalar s; int sum=0; for (int i=0;i<img0->height;i+) for (int j=0;j<img0->width

41、;j+) s = cvGet2D(img0,i,j); sum = (s.val0+s.val1+s.val2)/3; /这里可能考虑了三通道的缘故 if (sum > 128) s.val0=s.val1=s.val2=255; cvSet2D(img0,i,j,s); else s.val0=s.val1=s.val2=0; cvSet2D(img0,i,j,s); cvNamedWindow("Result"); cvShowImage("Result",img0); cvWaitKey(0);6. 区域生长#include <cxco

42、re.h>#include <cv.h>#include <stack>#include <stdio.h>#include <opencv2/core/core.hpp>#include <opencv2/imgproc/imgproc.hpp>#include <opencv2/highgui/highgui.hpp>typedef struct/保存种子像素 int x; int y;seedpoint;/区域生长，8连通区域，基于多个种子点，种子点可以自选，种子点在seed图像内设置为255void Grow

43、(IplImage* src,IplImage* seed, int t1)/ stack <seedpoint> seedd; seedpoint point; / 获取图像数据,保存种子区域 int height = seed->height; int width = seed->width; int step = seed->widthStep; uchar* seed_data = (uchar *)seed->imageData; uchar* src_data=(uchar *)src->imageData; int temp; /将种子点

44、压入堆栈 for(int i=0;i<height;i+) for(int j=0;j<width;j+) if(seed_datai*step+j=255) point.x=i; point.y=j; temp = src_datapoint.x*step+point.y; seedd.push(point); /取出种子点，分别作8邻域生长 seedpoint temppoint; while(!seedd.empty() point=seedd.top(); seedd.pop(); if(point.x>0)&&(point.x<(height-1)&&(point.y>0)&&(point.y<(width-1) if(seed_data(point.x-1)*step+point.y=0)&&(abs(src_data(point.x-1)*step+point.y-temp) <= t1) seed_data(point.x-1)*step+point.y=255; temppoint.x=point.x-1; temppoint.y=