OTSU算法及其改进算法学习

这篇文章转自http://blog.csdn.net/jinshengtao/article/details/19506005点击打开链接

这篇文章还是来自斯坦福课后作业hw2_3，主要是结合一个例子介绍otsu算法【亦称为大律算法，小日本】及其改进算法。

   本文将先介绍老外的题目、解题思路及maltab解答，然后分析otsu算法步骤，末了给出opencv实现。

   

老外的题目：Binarization of Scanned Book Pages

题目大意：

网上图书服务，比如百度文库需要将大量藏书数字化。首先，书的每一页将被扫描。然后，这些扫描图片将被二值化，并通过字符识别引擎OCR处理，即图片转字符。对于传统书籍【由于装订原因，如果在不破坏书的情况下】，书的每一页被扫描时，由于纸张被弯曲导致扫描结果的光照不均匀。如下图所示：

现在要求：

1.     对每一幅图像使用otsu算法执行全局二值化处理，计算原始图像的直方图，并在该直方图上标注OTSU阈值

2.     对每一幅图像执行局部自适应阈值，根据局部变化区分对待均匀和非均匀区域。

 

解题思路：

第一题，直接使用matlab的graythresh函数，通过最大类间方差法【OTSU】找到图片的一个合适的阈值（threshold）。末了，用imhist求取直方图便是。

第二题，使用一个水平滑窗，大小为21列宽*图像原始高度，从左往右逐像素滑动。对于窗口内的像素，计算局部变化【方差或平均值，代码用的是方差】。若窗口内方差大于阈值，使用otsu算法计算窗口内局部阈值，并二值化该窗口内像素；若方差小于阈值，则是书页上的空白区域，将该窗口内所有像素设为白色。如下图所示：

【题外话，本题的意思就是在对图书使用OCR进行字符识别前，优化二值化结果，使得OCR结果更精确。】

matlab代码：

[plain] view plaincopy

1. clc; clear all;  
2. imageFiles = {'hw2_book_page_1.jpg', 'hw2_book_page_2.jpg'};  
3. for nImage = 1:length(imageFiles)  
4.     % Load image  
5.     img = im2double(imread(imageFiles{nImage}));  
6.     figure(1); clf;  
7.     imshow(img);  
8.     [height, width] = size(img);  
9.    % [pathStr, name, ext] = fileparts(imageFiles{nImage});  
10.     % Global thresholding  
11.     globalThresh = graythresh(img);  
12.     imgBinGlobal = im2bw(img, globalThresh);%Convert image to binary image, based on threshold  
13.     figure(2); clf;  
14.     imshow(imgBinGlobal);  
15.     figure(3); clf; set(gcf, 'Color', 'w');  
16.     imhist(img); hold on;  
17.     histCounts = imhist(img);  
18.     h = plot(globalThresh*ones(1,100), linspace(0,max(histCounts)), 'r-');  
19.     set(h, 'LineWidth', 2);  
20.     set(gca, 'FontSize', 26);  
21.     h = text(globalThresh+0.01, max(histCounts)/4, ...  
22.         sprintf('T = %.2f', globalThresh));  
23.     set(h, 'FontSize', 26);  
24.     ylabel('Frequency');  
25.    % imwrite(imgBinGlobal, ['Global_' name '.jpg']);  
26.     % Locally adaptive thresholding  
27.     imgBinLocal = imgBinGlobal;  
28.     winHalfWidth = 10;  
29.     localVarThresh = 0.002;  
30.     for col = 1:width  
31.         inCols = max(1,col-winHalfWidth) : min(width,col+winHalfWidth);  
32.         inRows = 1:height;  
33.         inTile = img(inRows, inCols);  
34.         localThresh = graythresh(inTile);  
35.         %localMean = mean2(inTile);  
36.         localVar = std(inTile(:))^2;    %方差  
37.         if localVar > localVarThresh  
38.             imgBinLocal(:,col) = im2bw(img(:,col), localThresh);  
39.         else  
40.             imgBinLocal(:,col) = 1;  
41.         end  
42.     end % col  
43.     figure(4); clf;  
44.     imshow(imgBinLocal);  
45.    % imwrite(imgBinLocal, ['Local_' name '.jpg']);  
46.     if nImage == 1  
47.         pause  
48.     end  
49. end % nImage  

实验结果：

hw2_book_page_1.jpg 原始图像：

全局阈值化处理结果，即全局otsu结果：

 

全局OTSU结果在灰度直方图中的位置【注意这里所有的灰度都被缩放到0-1之间，包括阈值才0.65，我后面自己实现的要167】：

局部OTSU效果：

可以明显发现大片黑色木有了，得，另一副图片的结果，大家自己去斯坦福下载学习吧。

下面着重介绍OTSU算法原理及实现：

内容参考原文《A Threshold Selection Method from Gray-Level Histograms》

最大类间方差是由日本学者大津(Nobuyuki Otsu)于1979年提出，是一种自适应的阈值确定方法。算法假设图像像素能够根据阈值，被分成背景[background]和目标[objects]两部分。然后，计算该最佳阈值来区分这两类像素，使得两类像素区分度最大【用方差表达，具体公式见后】。OTSU的扩展算法，可进行多级阈值处理，称为“Multi Otsu method”【题外话】

设原始灰度级为M，灰度级为i的像素点个数为n_i，对灰度直方图进行归一化：

 

opencv实现代码：

[cpp] view plaincopy

1. // m_otsu.cpp : 定义控制台应用程序的入口点。  
2. //  
3.   
4. #include "stdafx.h"  
5. #include "cv.h"  
6. #include "highgui.h"  
7.   
8. int _tmain(int argc, _TCHAR* argv[])  
9. {  
10.     int i,j,nThresh;  
11.     int nHistogram[256] = {0};  
12.     double fStdHistogram[256] = {0.0};  
13.     double fGrayAccu[256] = {0.0};  
14.     double fGrayAve[256] = {0.0};  
15.     double fAverage = 0;  
16.     double fTemp = 0;  
17.     double fMax = 0;  
18.     IplImage *src,*dst;  
19.       
20.       
21.     src = cvLoadImage("test.jpg",CV_LOAD_IMAGE_GRAYSCALE);  
22.     dst = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);  
23.   
24.     //统计直方图  
25.     // 每行  
26.     for(i = 0; i < src->height; i++)  
27.     {  
28.         // 每列  
29.         for(j = 0; j < src->width; j++)  
30.         {  
31.           
32.             nHistogram[(unsigned char)src->imageData[i*src->width+j]] ++;  
33.         }  
34.     }  
35.       
36.     //归一化直方图  
37.     for(i = 0; i <= 255;i++)  
38.     {  
39.         fStdHistogram[i] = nHistogram[i]/(double)(src->width * src->height);  //Pi  
40.         //printf("%f\n",fStdHistogram[i]);  
41.     }  
42.       
43.     for(i=0;i<=255;i++)  
44.     {  
45.         for(j=0;j<=i;j++)  
46.         {  
47.             fGrayAccu[i] += fStdHistogram[j];                                   //所有灰度级，关于w0的数组                           
48.             fGrayAve[i] += j*fStdHistogram[j];                                  //所有灰度级，关于u(t)的数组  
49.         }  
50.   
51.         fAverage += i*fStdHistogram[i];                                         //uT  
52.         //printf("%f\n",fAverage);  
53.     }  
54.   
55.   
56.     //计算OSTU  
57.     for(i=0;i<=255;i++)  
58.     {  
59.         fTemp=(fAverage*fGrayAccu[i]-fGrayAve[i])*(fAverage*fGrayAccu[i]-fGrayAve[i])/(fGrayAccu[i]*(1-fGrayAccu[i]));  
60.   
61.         if(fTemp>fMax)  
62.         {  
63.             fMax=fTemp;  
64.             nThresh=i;  
65.         }  
66.     }  
67.   
68.     //计算二值图像  
69.     for (i=0;i<src->height;i++)  
70.     {  
71.         for (j=0;j<src->width;j++)  
72.         {  
73.             if ((unsigned char)src->imageData[i*src->width+j]<nThresh)  
74.             {  
75.                 dst->imageData[i*src->width+j] = 0;  
76.             }else{  
77.                 dst->imageData[i*src->width+j] = 255;  
78.             }  
79.         }  
80.     }  
81.     printf("%d",nThresh);  
82.     cvNamedWindow("otsu",0);  
83.     cvShowImage("otsu",dst);  
84.     cvSaveImage("otsu_result.jpg",dst);  
85.     cvWaitKey(0);  
86.     return 0;  
87. }  

实验结果：【全局OTSU】

阈值为167.