Opencv源码解析------二值化
二值化的4种方法：
（1） 对于RGB图像进行灰度化，然后扫描图像的每个像素。大于阈值的设为255（白色）。小于阈值的设为0（黑色）。
（2） 计算图像像素的均值K，然后扫描图像的每个像素，像素值大于K的设为255（白色），小于K的设为黑色（0）。
（3） 大津法：利用直方图的形式统计像素的分布情况然后，求得图像直方图的两个最高峰然后将两峰之间的峰谷最低处的值作为阈值进行分割。
（4） 方法四：http://en.wikipedia.org/wiki/Thresholding_(image_processing)使用近似一维Means方法寻找二值化阈值，该方法的大致步骤如下：
1. 一个初始化阈值T，可以自己设置或者根据随机方法生成。
2. 根据阈值图每个像素数据P(n,m)分为对象像素数据G1与背景像素数据G2。(n为行，m为列)
3. G1的平均值是m1, G2的平均值是m2
4. 一个新的阈值T’ = (m1 +m2)/2
5. 回到第二步，用新的阈值继续分像素数据为对象与北京像素数据，继续2～4步，直到计算出来的新阈值等于上一次阈值。
源码实现：
//Haartest.cpp : 定义控制台应用程序的入口点。
#include "stdafx.h"
#include <iostream>
#include <opencv2\opencv.hpp>
using namespace std;
using namespace cv;
//voidcolor(Mat &image)
//{
// for(int i=0;i<image.rows;i++)
// {
// for(int j=0;j<image.cols;j++)
// {
// if(image.at<Vec3b>(i,j)[0]==255&&image.at<Vec3b>(i,j)[1]==255&&image.at<Vec3b>(i,j)[0]==255)
// {
// image.at<Vec3b>(i,j)[0]=0;
// image.at<Vec3b>(i,j)[1]=255;
// image.at<Vec3b>(i,j)[2]=0;
// }
//
// }
// }
//}
static double
getThreshVal_Otsu_8u(const Mat& _src ) //大津法求值
{
Size size = _src.size();
if(_src.isContinuous() ) //判断是否继续
{
size.width *= size.height;
size.height = 1;
}
const int N = 256; //256个值统计像素的分布情况
int i, j,h[N] = {0};
for( i = 0;i < size.height; i++ )
{
const uchar* src =_src.data + _src.step*i; //利用指针的方式遍历图像
j = 0;
#ifCV_ENABLE_UNROLLED
for( ;j <= size.width - 4; j += 4 )
{
intv0 = src[j], v1 = src[j+1];
h[v0]++; h[v1]++; //统计像素值出现的次数。。。
v0 = src[j+2]; v1 = src[j+3];
h[v0]++; h[v1]++;
}
#endif
for( ;j < size.width; j++ )//统计像素出现的次数
h[src[j]]++;
}
double mu =0, scale = 1./(size.width*size.height); //计算尺度的大小
for( i = 0;i < N; i++ )
mu += i*(double)h[i]; //计算所有像素的和
mu *= scale; //计算像素均值，图像的总平均灰度为mu
double mu1= 0, q1 = 0;
doublemax_sigma = 0, max_val = 0;
for( i = 0;i < N; i++ )
{
doublep_i, q2, mu2, sigma;
p_i = h[i]*scale; //像素i出现的概率 前景像素占图像的比例
mu1 *= q1;
q1 += p_i; //像素i出现的概率的和
q2 = 1 - q1; //背景像素占图像的比例
if(std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )
continue;
mu1 = (mu1 + i*p_i)/q1; //前景像素的均值
mu2 = (mu - q1*mu1)/q2; //背景像素的均值
//************类间方差sigma的计算公式为sigma =q1*q2*(mu1-mu2)^2*****************//
sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
if(sigma > max_sigma )
{
max_sigma = sigma;
max_val = i; //求得使类间方差最大的阈值
}
}
returnmax_val;
}
static void
thresh_8u( const Mat& _src, Mat& _dst, uchar thresh, ucharmaxval, int type )
{
int i, j,j_scalar = 0;
uchartab[256];
Size roi = _src.size();
roi.width *= _src.channels();
if(_src.isContinuous() && _dst.isContinuous() )
{
roi.width *= roi.height;
roi.height = 1;
}
#ifdef HAVE_TEGRA_OPTIMIZATION
if(tegra::thresh_8u(_src, _dst, roi.width, roi.height, thresh, maxval, type))
return;
#endif
switch(type )
{
caseTHRESH_BINARY: //边缘性dst(x,y) = max_value, if src(x,y)>threshold 0, otherwise.
for( i= 0; i <= thresh; i++ )
tab[i] = 0;
for( ;i < 256; i++ )
tab[i] = maxval;
break;
caseTHRESH_BINARY_INV://dst(x,y) = 0, ifsrc(x,y)>threshold; dst(x,y) = max_value, otherwise.
for( i= 0; i <= thresh; i++ )
tab[i] = maxval;
for( ;i < 256; i++ )
tab[i] = 0;
break;
caseTHRESH_TRUNC: //dst(x,y)= threshold, if src(x,y)>threshold; dst(x,y) = src(x,y), otherwise.
for( i= 0; i <= thresh; i++ )
tab[i] = (uchar)i;
for( ;i < 256; i++ )
tab[i] = thresh;
break;
caseTHRESH_TOZERO://dst(x,y) = src(x,y), if(x,y)>threshold ; dst(x,y) = 0,otherwise.
for( i= 0; i <= thresh; i++ )
tab[i] = 0;
for( ;i < 256; i++ )
tab[i] = (uchar)i;
break;
caseTHRESH_TOZERO_INV://dst(x,y) = 0, ifsrc(x,y)>threshold ; dst(x,y) =src(x,y), otherwise.
for( i= 0; i <= thresh; i++ )
tab[i] = (uchar)i;
for( ;i < 256; i++ )
tab[i] = 0;
break;
default:
CV_Error( CV_StsBadArg, "Unknown threshold type" );
}
#if CV_SSE2//如果条件编译符号CV_SSE被定义过，则调用checkHardwareSupport(CV_CPU_SSE).获取
//你的电脑上的CPU是否支持CPU_SSE指令集优化特性，如果为真支持，为假不支持
if(checkHardwareSupport(CV_CPU_SSE2) )
{
__m128i_x80 = _mm_set1_epi8('\x80');
__m128ithresh_u = _mm_set1_epi8(thresh);
__m128ithresh_s = _mm_set1_epi8(thresh ^ 0x80);
__m128imaxval_ = _mm_set1_epi8(maxval);
j_scalar = roi.width & -8;
for( i= 0; i < roi.height; i++ )
{
constuchar* src = (constuchar*)(_src.data + _src.step*i);
uchar*dst = (uchar*)(_dst.data + _dst.step*i);
switch(type )
{
caseTHRESH_BINARY:
for(j = 0; j <= roi.width - 32; j += 32 )
{
__m128iv0, v1;
v0 = _mm_loadu_si128( (const __m128i*)(src +j) );
v1 = _mm_loadu_si128( (const __m128i*)(src +j + 16) );
v0 = _mm_cmpgt_epi8( _mm_xor_si128(v0,_x80), thresh_s );
v1 = _mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s );
v0 = _mm_and_si128( v0,maxval_ );
v1 = _mm_and_si128( v1,maxval_ );
_mm_storeu_si128( (__m128i*)(dst + j), v0 );
_mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
}
for(; j <= roi.width - 8; j += 8 )
{
__m128iv0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
v0 = _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s );
v0 = _mm_and_si128( v0,maxval_ );
_mm_storel_epi64( (__m128i*)(dst + j), v0 );
}
break;
caseTHRESH_BINARY_INV:
for(j = 0; j <= roi.width - 32; j += 32 )
{
__m128iv0, v1;
v0 = _mm_loadu_si128( (const __m128i*)(src +j) );
v1 = _mm_loadu_si128( (const __m128i*)(src +j + 16) );
v0 = _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s );
v1 = _mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s );
v0 = _mm_andnot_si128( v0,maxval_ );
v1 = _mm_andnot_si128( v1,maxval_ );
_mm_storeu_si128( (__m128i*)(dst + j), v0 );
_mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
}
for(; j <= roi.width - 8; j += 8 )
{
__m128iv0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
v0 = _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s );
v0 = _mm_andnot_si128( v0, maxval_ );
_mm_storel_epi64( (__m128i*)(dst + j), v0 );
}
break;
caseTHRESH_TRUNC:
for(j = 0; j <= roi.width - 32; j += 32 )
{
__m128iv0, v1;
v0 = _mm_loadu_si128( (const __m128i*)(src +j) );
v1 = _mm_loadu_si128( (const __m128i*)(src +j + 16) );
v0 = _mm_subs_epu8( v0,_mm_subs_epu8( v0, thresh_u ));
v1 = _mm_subs_epu8( v1,_mm_subs_epu8( v1, thresh_u ));
_mm_storeu_si128( (__m128i*)(dst + j), v0 );
_mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
}
for(; j <= roi.width - 8; j += 8 )
{
__m128iv0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
v0 = _mm_subs_epu8( v0,_mm_subs_epu8( v0, thresh_u ));
_mm_storel_epi64( (__m128i*)(dst + j), v0 );
}
break;
caseTHRESH_TOZERO:
for(j = 0; j <= roi.width - 32; j += 32 )
{
__m128iv0, v1;
v0 = _mm_loadu_si128( (const __m128i*)(src +j) );
v1 = _mm_loadu_si128( (const __m128i*)(src +j + 16) );
v0 = _mm_and_si128( v0,_mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ));
v1 = _mm_and_si128( v1,_mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s ));
_mm_storeu_si128( (__m128i*)(dst + j), v0 );
_mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
}
for(; j <= roi.width - 8; j += 8 )
{
__m128iv0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
v0 = _mm_and_si128( v0,_mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ));
_mm_storel_epi64( (__m128i*)(dst + j), v0 );
}
break;
caseTHRESH_TOZERO_INV:
for(j = 0; j <= roi.width - 32; j += 32 )
{
__m128iv0, v1;
v0 = _mm_loadu_si128( (const __m128i*)(src +j) );
v1 = _mm_loadu_si128( (const __m128i*)(src +j + 16) );
v0 = _mm_andnot_si128(_mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ), v0 );
v1 = _mm_andnot_si128(_mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s ), v1 );
_mm_storeu_si128( (__m128i*)(dst + j), v0 );
_mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
}
for(; j <= roi.width - 8; j += 8 )
{
__m128iv0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
v0 = _mm_andnot_si128(_mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ), v0 );
_mm_storel_epi64( (__m128i*)(dst + j), v0 );
}
break;
}
}
}
#endif
if(j_scalar < roi.width )
{
for( i= 0; i < roi.height; i++ )
{
constuchar* src = (constuchar*)(_src.data + _src.step*i);
uchar*dst = (uchar*)(_dst.data + _dst.step*i);
j = j_scalar;
#if CV_ENABLE_UNROLLED //这段代码没看懂什么意思但是发现没多大影响
for(; j <= roi.width - 4; j += 4 )
{
uchart0 = tab[src[j]];
uchart1 = tab[src[j+1]];
dst[j] = t0;
dst[j+1] = t1;
t0 = tab[src[j+2]];
t1 = tab[src[j+3]];
dst[j+2] = t0;
dst[j+3] = t1;
}
#endif
for(; j < roi.width; j++ )
dst[j] = tab[src[j]];
}
}
}
class ThresholdRunner : publicParallelLoopBody //进行并行计算提高算法的速度
{
public:
ThresholdRunner(Mat _src, Mat _dst, double _thresh, double_maxval, int _thresholdType)
{
src = _src;
dst = _dst;
thresh = _thresh;
maxval = _maxval;
thresholdType = _thresholdType;
}
void operator () ( constRange& range ) const //赋值函数的形式
{
introw0 = range.start;
introw1 = range.end;
Mat srcStripe = src.rowRange(row0, row1);
Mat dstStripe = dst.rowRange(row0,row1);
if(srcStripe.depth() == CV_8U)
{
thresh_8u( srcStripe, dstStripe, (uchar)thresh, (uchar)maxval,thresholdType );
}
}
private:
Mat src;
Mat dst;
intnStripes;
doublethresh;
doublemaxval;
intthresholdType;
};
double threshold1( Mat src, Mat dst, double thresh, doublemaxval, int type )
{
booluse_otsu = (type & THRESH_OTSU) != 0;
type &= THRESH_MASK;
if(use_otsu ) //是否应用大津法求阈值
{
CV_Assert( src.type() == CV_8UC1 ); //若括号中的表达式为false则返回一个错误信息
thresh =getThreshVal_Otsu_8u(src); //大津法得到阈值
}
dst.create( src.size(), src.type() );
if( src.depth()== CV_8U )
{
intithresh = cvFloor(thresh);
//**********cvRound 返回和参数最接近的整数值。 cvFloor 返回不大于参数的最大整数值
//cvCeil 返回不小于参数的最小整数值。//
thresh = ithresh;
intimaxval = cvRound(maxval);
if(type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<uchar>(imaxval);
//saturate_cast防止数据溢出大致原理如下：
/*if(data<0)
data=0;
else if(data>255)
data=255;*/
if(ithresh < 0 || ithresh >= 255 )
{
if(type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type== THRESH_TOZERO_INV) && ithresh < 0) ||
(type == THRESH_TOZERO&& ithresh >= 255) )
{
intv = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
type ==THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
returnthresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_16S )
{
intithresh = cvFloor(thresh);
thresh = ithresh;
intimaxval = cvRound(maxval);
if(type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<short>(imaxval);
if(ithresh < SHRT_MIN || ithresh >= SHRT_MAX )
{
if(type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type== THRESH_TOZERO_INV) && ithresh < SHRT_MIN)||
(type == THRESH_TOZERO&& ithresh >= SHRT_MAX) )
{
intv = type == THRESH_BINARY ? (ithresh >= SHRT_MAX? 0 : imaxval) :
type == THRESH_BINARY_INV ?(ithresh >= SHRT_MAX ? imaxval : 0) :
/*type== THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
returnthresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_32F )
;
else
CV_Error( CV_StsUnsupportedFormat, "" );
parallel_for_(Range(0, dst.rows),
ThresholdRunner(src, dst,thresh, maxval, type),
dst.total()/(double)(1<<16));
returnthresh;
}
int _tmain(int argc,_TCHAR* argv[])
{
Mat img = imread("lena.bmp");
Mat dst;
doublevalue=0.0;
cvtColor(img,dst,CV_BGR2GRAY);
value=threshold1(dst,dst,124,255,THRESH_BINARY);
imshow("src",img); //二值化的原图
imshow("dst",dst);
waitKey(0);
return 0;
}
实验结果：
4051
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
