2019-10-25 23:21:22 qq_43294951 阅读数 24

由于OTSU比较简单,就是一个遍历来寻找最大的阈值分割的灰度值点,故本例直接给出代码,不做细评,代码如下,具体的写在代码注释

#include <iostream>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <math.h>
#include <stack>
#define _MATH_DEFINES_DEFINED
using namespace std;
using namespace cv;
using namespace Eigen;
float Otsu(const Mat& image,int T)
{
	int nr = image.rows;
	int nl = image.cols;
	Mat hist = Mat::zeros(1, 256, CV_32F);
	float* hi = hist.ptr<float>(0);
	float junzhi=0;//定义总均值,类一均值,类二均值
	float junzhi0 = 0;
	float junzhi1 = 0;
	float lei0 = 0;//定义类一,类二的概率
	float lei1 = 0;
	float fangcha0=0;//定义类一,类二的方差
	float fangcha1 = 0;
	float neifangcha = 0;//定义类内方差
	float jianfangcha = 0;//定义类间方差
	float zongfangcha = 0;//定义总方差,类内加类间
	float ppppp = 0;
	for (int i = 0; i < nr; i++)//这里是统计直方图
	{
		const float* im = image.ptr<float>(i);
		for (int j = 0; j < nl; j++)
		{
			hi[int(im[j])] = hi[int(im[j])] +1;
		}
	}
	//cout << hist << endl;
	for (int i = 0; i < 256; i++)//计算均值
	{
		junzhi = i * hi[i] / nr / nl;
	}
	for (int i = 0; i <= T; i++)//计算
	{
		lei0 = lei0 + hi[i] / nr / nl;
	}
	for (int i = T+1; i <256; i++)
	{
		lei1 = lei1 + hi[i] / nr / nl;
	}
	for (int i = 0; i <= T; i++)
	{
		junzhi0 = junzhi0 + i * hi[i] / nr / nl / lei0;
	}
	for (int i = T + 1; i < 256; i++)
	{
		junzhi1 = junzhi1 + i * hi[i] / nr / nl / lei1;
	}
	for (int i = 0; i <= T; i++)
	{
		fangcha0 = fangcha0 + (i -junzhi0)* (i - junzhi0) * hi[i] / nr / nl / lei0;
	}
	for (int i = T + 1; i < 256; i++)
	{
		fangcha1 = fangcha1 + (i - junzhi1) * (i - junzhi1) * hi[i] / nr / nl / lei1;
	}
	neifangcha = lei0 * fangcha0 + lei1 * fangcha1;
	jianfangcha = lei0 * lei1 * (junzhi1 - junzhi0) * (junzhi1 - junzhi0);
	zongfangcha = neifangcha + jianfangcha;
	ppppp= jianfangcha / zongfangcha;
	return ppppp;
}
int getT(const Mat& image, Mat& new_image)
{
	int temp = 0;
	float max = 0;
	max = Otsu(image, 20);
	for (int i =  0; i < 255; i++)
	{
		float gg = Otsu(image, i);
		if (max < gg)
		{
			temp = i;
			max =gg;
		}
	}//遍历得到最大的分割点temp
	for (int i = 0; i < image.rows; i++)
	{
		const float* im =image.ptr<float>(i);
		float* p = new_image.ptr<float>(i);
		for (int j = 0; j < image.cols; j++)
		{
			float ge = im[j];
			if (ge <= temp)
			{
				p[j] = 0;
			}
			else
			{
				p[j] = 255;
			}

		}
	}//图像进行二值化
	return temp;
}
int main()
{
	Mat image = imread("I:/C.jpg", IMREAD_GRAYSCALE);
	imshow("原图", image);
	Mat IM;
	image.convertTo(IM, CV_32F);
	//imshow("fsSSd", IM);
	Mat new_image = Mat::zeros(IM.rows, IM.cols, CV_32F);
	int T = getT(IM, new_image);
	cout << T << endl;
	imshow("分割以后", new_image);
	cv::waitKey(0);
}

ok,得到的阈值分割图如下,感觉海星
在这里插入图片描述

2019-04-22 16:24:12 Lemon_jay 阅读数 529

使用C++、opencv对图像进行简单的阈值分割

对图像进行颜色的阈值分割,往往不知道阈值设置为多少合适,需要不断测试,针对此问题,设置了阈值的滑动条,从而可以通过滑动滑动条快速找到阈值。

最常见的阈值分割即为R、G、B的阈值分割,即判断像素点的R、G、B值是否大于或小于某一阈值,满足条件的点就保留或删除。

本代码阈值分割的算法是直接按行按列依次遍历图像每个像素点,判断像素点的值是否满足分割条件,满足即将该点设置为黑色来达到分割的效果。代码为阈值、 R、G、B、和模式都设置了滑动条,其中模式有4种:

0----RGB        R、G、B值大于阈值的点设为黑色
1----|G-R|        G-R值的绝对值小于threshod_value的点设为黑色
2----|2G-R-B|        2G-R-B值的绝对值大于threshod_value的点设为黑色
3----G/R        G/R值小于threshod_value的点设为黑色


创建滑动条的API:

int createTrackbar(const string& trackbarname, const string& winname, int* value, int count, TrackbarCallback onChange=0,                                    void* userdata=0)

第一个参数,const string&类型的trackbarname,轨迹条的名字,用来代表我们创建的轨迹条。

第二个参数,const string&类型的winname,窗口的名字,表示这个轨迹条会依附到哪个窗口上,即对应named Window创建窗口时填的某一个窗口。

第三个参数,int*类型的value,一个指向整型的指针,表示滑块的位置。在创建时,滑块的初始位置就是该变量当前的值。

第四个参数,int类型的count,表示滑块可以达到的最大位置的值。滑块最小位置的值始终为0。

第五个参数,Trackbar Callback类型的on Change,它有默认值0。这是一个指向回调函数的指针,每次滑块位置改变时,这个函数都会进行回调。并且这个函数的原型必须为void XXXX(int,void*),其中第一个参数是轨迹条的位置,第二个参数是用户数据(看下面的第六个参数)。如果回调是NULL指针,则表示没有回调函数的调用,仅第三个参数value有变化。

第六个参数,void*类型的userdata,也有默认值0。这个参数是用户传给回调函数的数据,用来处理轨迹条事件。如果使用的第三个参数value实参是全局变量的话,完全可以不去管这个userdata参数。


代码实现:

#include "stdafx.h"
#include <opencv2\opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
using namespace cv;

Mat src,dst;

//滑动条相关全局变量
int mode = 0;
int r = 255, g = 255, b = 255;
int threshod_value = 70;

//滑动条回调函数声明
void on_mode(int, void*);
void on_threshold(int, void*);

//阈值处理函数声明
void RGB_threshold();
void G_B_threshold();
void G_B_R_threshold();
void G_R_threshold();

//帮助文本
void show_help()
{
	std::cout << "---------help---------" << std::endl;
	std::cout << "mode:" << std::endl;
	std::cout << "0----RGB		R、G、B值大于阈值的点设为黑色" << std::endl;
	std::cout << "1----|G-R|		G-R值的绝对值小于threshod_value的点设为黑色" << std::endl;
	std::cout << "2----|2G-R-B|		2G-R-B值的绝对值大于threshod_value的点设为黑色" << std::endl;
	std::cout << "3----G/R		G/R值小于threshod_value的点设为黑色" << std::endl;
}
int main()
{
	system("color 02");

	show_help();

	//读取图像
	src = imread("D:\\哀.jpg");
	namedWindow("src", WINDOW_NORMAL);
	imshow("src", src);
	namedWindow("dst", WINDOW_NORMAL);

	//创建滑动条
	createTrackbar("mode", "src", &mode, 4, on_mode);
	on_mode(0, 0);

	createTrackbar("R", "src", &r, 255, on_threshold);
	on_threshold(0, 0);
	
	createTrackbar("G", "src", &g, 255, on_threshold);
	on_threshold(0, 0);
	
	createTrackbar("B", "src", &b, 255, on_threshold);
	on_threshold(0, 0);
	
	createTrackbar("thre_val", "src", &threshod_value, 255, on_threshold);
	on_threshold(0, 0);

	waitKey();
	return 0;
}

void on_mode(int, void*)
{
}
void on_threshold(int, void*)
{
	switch (mode)
	{
	case 0:
		RGB_threshold();
		break;
	case 1:
		G_B_threshold();
		break;
	case 2:
		G_B_R_threshold();
		break;
	case 3:
		G_R_threshold();
		break;
	}
}

void RGB_threshold()
{
	//复制源图到目标图像
	dst = src.clone();
	//遍历图像每个像素点
	for(int i=0;i<src.rows;i++)
		for (int j = 0; j < src.cols; j++)
		{
			//大于阈值的点设为黑色
			if (src.at<Vec3b>(i, j)[2] > r)
			{
				dst.at<Vec3b>(i, j)[0] = 0;
				dst.at<Vec3b>(i, j)[1] = 0;
				dst.at<Vec3b>(i, j)[2] = 0;
			}
			if (src.at<Vec3b>(i, j)[1] > g)
			{
				dst.at<Vec3b>(i, j)[0] = 0;
				dst.at<Vec3b>(i, j)[1] = 0;
				dst.at<Vec3b>(i, j)[2] = 0;
			}
			if (src.at<Vec3b>(i, j)[0] > b)
			{
				dst.at<Vec3b>(i, j)[0] = 0;
				dst.at<Vec3b>(i, j)[1] = 0;
				dst.at<Vec3b>(i, j)[2] = 0;
			}
		}
	imshow("dst", dst);
}
void G_B_threshold()
{
	dst = src.clone();
	for (int i = 0; i<src.rows; i++)
		for (int j = 0; j < src.cols; j++)
		{
			if ((    abs(src.at<Vec3b>(i, j)[1] - src.at<Vec3b>(i, j)[2])	<	threshod_value	))
			{
				dst.at<Vec3b>(i, j)[0] = 0;
				dst.at<Vec3b>(i, j)[1] = 0;
				dst.at<Vec3b>(i, j)[2] = 0;
			}
		}
	//开操作
	//Mat element = getStructuringElement(MORPH_RECT, Size(3, 3));
	//morphologyEx(dst, dst, MORPH_OPEN, element);
	imshow("dst", dst);
}
void G_B_R_threshold()
{
	dst = src.clone();
	for (int i = 0; i<src.rows; i++)
		for (int j = 0; j < src.cols; j++)
		{
			if ((abs(2*src.at<Vec3b>(i, j)[1] - src.at<Vec3b>(i, j)[2] - src.at<Vec3b>(i, j)[0])	>	threshod_value))
			{
				dst.at<Vec3b>(i, j)[0] = 0;
				dst.at<Vec3b>(i, j)[1] = 0;
				dst.at<Vec3b>(i, j)[2] = 0;
			}
		}
	//腐蚀操作
	//Mat element = getStructuringElement(MORPH_RECT, Size(3, 3));
	//morphologyEx(dst, dst, MORPH_ERODE, element);
	imshow("dst", dst);
}
void G_R_threshold()
{
	dst = src.clone();
	int height = dst.rows;
	int width = dst.cols;
	for (int row = 0; row < height; row++)
	{
		for (int col = 0; col < width; col++)
		{
			double a = (double)dst.at<Vec3b>(row, col)[1] / dst.at<Vec3b>(row, col)[2];
			a = a * 100;
			if (a < threshod_value)
			{
				dst.at<Vec3b>(row, col)[0] = 0;
				dst.at<Vec3b>(row, col)[1] = 0;
				dst.at<Vec3b>(row, col)[2] = 0;
			}
		}
	}
	imshow("dst", dst);
}

结果图:

2013-01-24 14:14:54 ily6418031hwm 阅读数 13028

       由于阈值处理直观、实现简单且计算速度快,因此图像阈值处理在图像分割中处于核心地位,下面我会为大家介绍阈值处理的方法,并用OpenCV给出实现的代码。

第一种:图像阈值分割---基本的全局阈值

1.  处理流程:

             1.为全局阈值选择一个初始估计值T(图像的平均灰度)。
             2.用T分割图像。产生两组像素:G1有灰度值大于T的像素组成,G2有小于等于T像素组成。
             3.计算G1和G2像素的平均灰度值m1和m2;
             4.计算一个新的阈值:T = (m1 + m2) / 2;
             5.重复步骤2和4,直到连续迭代中的T值间的差为零。

             此种方法主要适用于图像直方图有明显波谷。

2.源代码:      

/******************************************************************************* 
函数名称  : BasicGlobalThreshold 
函数描述  : 图像阈值分割---基本的全局阈值(适合图像直方图有明显波谷)
输入参数  : N/A
输出参数  : N/A
返 回 值    : int k1---------------------------------迭代得到的阈值
作     者    : N/A    2013-01-22
处理流程:1.为全局阈值选择一个初始估计值T(图像的平均灰度)。
          2.用T分割图像。产生两组像素:G1有灰度值大于T的像素组成,G2有小于等于T像素组成。
          3.计算G1和G2像素的平均灰度值m1和m2;
          4.计算一个新的阈值:T = (m1 + m2) / 2;
          5.重复步骤2和4,直到连续迭代中的T值间的差为零。
*******************************************************************************/ 
 int CThreasholdProcess::BasicGlobalThreshold(int*pg,int start,int end)
 {
	 int  i,t,t1,t2,k1,k2;
	 double u,u1,u2;    
	 t=0;     
	 u=0;

	 for (i=start;i<end;i++) 
	 {
		 t+=pg[i];        
		 u+=i*pg[i];
	 }

	 k2=(int) (u/t);                          //  计算此范围灰度的平均值    
	 do 
	 {
		 k1=k2;
		 t1=0;    
		 u1=0;
		 for (i=start;i<=k1;i++) 
		 {             //  计算低灰度组的累加和
			 t1+=pg[i];    
			 u1+=i*pg[i];
		 }

		 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(k1!=k2);                           //  数据未稳定,继续

	 return k1;                              //  返回阈值
 }
3.演示结果

a).被处理的源图像                                                                                                                         b).基本的全局阈值处理后图像                                                                      c)源图直方图

  

           

2012-02-11 16:18:32 renshengrumenglibing 阅读数 48377

在图像处理时,受外界光线的干扰一般比较大,假如在阈值分割时采用固


定阈值,那么在环境改变时分割效果受影响极大,那么为了避免此影响就


必须采用动态阈值,自动求出合适的阈值进行分割。
本文的介绍几种主要的图像分割方法,并给出自动阈值分割的源代码




图像分割是图像处理与计算机视觉领域低层次视觉中最为基础和重要的领域之一,它是对图像进行视觉分析和模式识别的基本前提.阈值法是一种传统的图像分割方法,因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术.已被应用于很多的领域。本文是在阅读大量国内外相关文献的基础上,对阈值分割技术稍做总结,分三个大类综述阈值选取方法,然后对阈值化算法的评估做简要介绍。


1.引言
所谓图像分割是指根据灰度、彩色、空间纹理、几何形状等特征把图像划分成若干个互不相交的区域,使得这些特征在同一区域内,表现出一致性或相似性,而在不同区域间表现出明显的不同[37].简单的讲,就是在一幅图像中,把目标从背景中分离出来,以便于进一步处理。图像分割是图像处理与计算机视觉领域低层次视觉中最为基础和重要的领域之一,它是对图像进行视觉分析和模式识别的基本前提.同时它也是一个经典难题,到目前为止既不存在一种通用的图像分割方法,也不存在一种判断是否分割成功的客观标准。
阈值法是一种传统的图像分割方法,因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术.已被应用于很多的领域,例如,在红外技术应用中,红外无损检测中红外热图像的分割,红外成像跟踪系统中目标的分割;在遥感应用中,合成孔径雷达图像中目标的分割等;在医学应用中,血液细胞图像的分割,磁共振图像的分割;在农业工程应用中,水果品质无损检测过程中水果图像与背景的分割。在工业生产中,机器视觉运用于产品质量检测等等。在这些应用中,分割是对图像进一步分析、识别的前提,分割的准确性将直接影响后续任务的有效性,其中阈值的选取是图像阈值分割方法中的关键技术。


2.阈值分割的基本概念
图像阈值化分割是一种最常用,同时也是最简单的图像分割方法,它特别适用于目标和背景占据不同灰度级范围的图像[1]。它不仅可以极大的压缩数据量,而且也大大简化了分析和处理步骤,因此在很多情况下,是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。图像阈值化的目的是要按照灰度级,对像素集合进行一个划分,得到的每个子集形成一个与现实景物相对应的区域,各个区域内部具有一致的属性,而相邻区域布局有这种一致属性。这样的划分可以通过从灰度级出发选取一个或多个阈值来实现。
阈值分割法是一种基于区域的图像分割技术,其基本原理是:通过设定不同的特征阈值,把图像像素点分为若干类.常用的特征包括:直接来自原始图像的灰度或彩色特征;由原始灰度或彩色值变换得到的特征.设原始图像为f(x,y),按照一定的准则在f(x,y)中找到特征值T,将图像分割为两个部分,分割后的图像为


 
若取 :b0=0(黑),b1=1(白),即为我们通常所说的图像二值化。


   
          (原始图像)    (阈值分割后的二值化图像)


一般意义下,阈值运算可以看作是对图像中某点的灰度、该点的某种局部特性以及该点在图像中的位置的一种函数,这种阈值函数可记作
    T(x,y,N(x,y),f(x,y))
式中,f(x,y)是点(x,y)的灰度值;N(x,y)是点(x,y)的局部邻域特性.根据对T的不同约束,可以得到3种不同类型的阈值[37],即
    点相关的全局阈值T=T(f(x,y))  
(只与点的灰度值有关)
区域相关的全局阈值T=T(N(x,y),f(x,y))  
(与点的灰度值和该点的局部邻域特征有关)
    局部阈值或动态阈值T=T(x,y,N(x,y),f(x,y))
(与点的位置、该点的灰度值和该点邻域特征有关)


图像阈值化这个看似简单的问题,在过去的四十年里受到国内外学者的广泛关注,产生了数以百计的阈值选取方法[2-9],但是遗憾的是,如同其他图像分割算法一样,没有一个现有方法对各种各样的图像都能得到令人满意的结果,甚至也没有一个理论指导我们选择特定方法处理特定图像。
所有这些阈值化方法,根据使用的是图像的局部信息还是整体信息,可以分为上下文无关(non-contextual)方法(也叫做基于点(point-dependent)的方法)和上下文相关(contextual)方法(也叫做基于区域(region-dependent)的方法);根据对全图使用统一阈值还是对不同区域使用不同阈值,可以分为全局阈值方法(global thresholding)和局部阈值方法(local thresholding,也叫做自适应阈值方法adaptive thresholding);另外,还可以分为双阈值方法(bilever thresholding)和多阈值方法(multithresholding)
本文分三大类对阈值选取技术进行综述:
1) 基于点的全局阈值方法;
2) 基于区域的全局阈值方法
3) 局部阈值方法和多阈值方法




3.基于点的全局阈值选取方法
3.1  p-分位数法
1962年Doyle[10]提出的p-分位数法(也称p-tile法)可以说是最古老的一种阈值选取方法。该方法使目标或背景的像素比例等于其先验概率来设定阈值,简单高效,但是对于先验概率难于估计的图像却无能为力。
例如,根据先验知识,知道图像目标与背景象素的比例为PO/PB,则可根据此条件直接在图像直方图上找到合适的阈值T,使得f(x,y)>=T的象素为目标,f(x,y)<T的象素为背景。


3.2  迭代方法选取阈值[11]
初始阈值选取为图像的平均灰度T0,然后用T0将图像的象素点分作两部分,计算两部分各自的平均灰度,小于T0的部分为TA,大于T0的部分为TB
计算   ,将T1 作为新的全局阈值代替T0,重复以上过程,如此迭代,直至TK 收敛,即TK+1 =TK 
经试验比较,对于直方图双峰明显,谷底较深的图像,迭代方法可以较快地获得满意结果。但是对于直方图双峰不明显,或图像目标和背景比例差异悬殊,迭代法所选取的阈值不如最大类间方差法。


3.3  直方图凹面分析法
从直观上说,图像直方图双峰之间的谷底,应该是比较合理的图像分割阈值,但是实际的直方图是离散的,往往十分粗糙、参差不齐,特别是当有噪声干扰时,有可能形成多个谷底。从而难以用既定的算法,实现对不同类型图像直方图谷底的搜索。
Rosenfeld和Torre[12]提出可以构造一个包含直方图 的最小凸多边形 ,由集差 确定 的凹面。若 和 分别表示 与 在灰度级之处的高度,则 取局部极大值时所对应的灰度级可以作为阈值。也有人使用低通滤波的方法平滑直方图,但是滤波尺度的选择并不容易[13]。
但此方法仍然容易受到噪声干扰,对不同类型的图像,表现出不同的分割效果。往往容易得到假的谷底。但此方法对某些只有单峰直方图的图像,也可以作出分割。如:
 


3.4 最大类间方差法
由Otsu[14]于1978年提出的最大类间方差法以其计算简单、稳定有效,一直广为使用。从模式识别的角度看,最佳阈值应当产生最佳的目标类与北京类的分离性能,此性能我们用类别方差来表征,为此引入类内方差 、类间方差 和总体方差 ,并定义三个等效的准则测量:
 ,  ,  .                 (3)
鉴于计算量的考量,人们一般通过优化第三个准则获取阈值。此方法也有其缺陷,kittler和Illingworth[15]的实验揭示:当图像中目标与背景的大小之比很小时方法失效。
 在实际运用中,往往使用以下简化计算公式:
   (T) = WA(μa-μ)2  + Wb(μb-μ)2 
其中, 为两类间最大方差,WA 为A类概率,μa为A类平均灰度,Wb 为B类概率,μb为B类平均灰度,μ为图像总体平均灰度。
即阈值T将图像分成A,B两部分,使得两类总方差 (T)取最大值的T,即为最佳分割阈值。


3.5 熵方法
八十年代以来,许多学者将Shannon信息熵的概念应用于图像阈值化,其基本思想都是利用图像的灰度分布密度函数定义图像的信息熵,根据假设的不同或视角的不同提出不同的熵准则,最后通过优化该准则得到阈值。Pun[16]通过使后验熵的上限最大来确定阈值。Kapur等人[17]的方法假定目标和背景服从两个不同的概率分布 和 定义
                 (4)
使得熵
                           (5)
达到最大求得最佳阈值。
此方法又称为KSW熵方法。


3.6 最小误差阈值
此方法来源于Bayes最小误差分类方法。
 
Eb(T)是目标类错分到背景类的概率,Eo(T)是背景类错分到目标类的概率
总的误差概率 E(T) = Eb(T) + Eo(T)
使E(T)取最小值,即为最优分类方法。


在Kittler和Illingworth[18]于1986年提出的最小误差法中,直方图被视为目标与背景混合集概率密度函数 的估计
                 (9)
其中, 为先验概率, ,求解下列方程可得到Bayes最小误差阈值
                        (10)
遗憾的是上式中 , 和 通常是未知的,Nakagawa和Rosenfeld[19]提倡用拟合方法从直方图中估计这些参数,但是算法相当复杂,不易实现。




3.7 矩量保持法
矩量保持(moment-preserving)法[20] ,即矩守恒阈值法,是1985年提出的,其基本思想是最佳的阈值应该使分割前后图像的矩量保持不变,由此可以得到一组矩量保持方程,求解该方程组就可以得到最佳阈值。


3.8 模糊集方法
模糊集理论较好的描述了人类视觉中的模糊性和随机性,因此在图像阈值化领域受到了广泛的关注。模糊集阈值化方法的基本思想是,选择一种S状的隶属度函数定义模糊集,隶属度为0.5的灰度级对应了阈值,当然在上述隶属度函数的表达式中阈值是一个未知的参数;然后在此模糊集上定义某种准则函数(例如整个图像的总体模糊度),通过优化准则函数来确定最佳阈值。
Pal等[21]首先,他们把一幅具有 个灰度级的 图像看作一个模糊集 ,其中隶属函数 定义如下:
               (11)
参数 称之为交叉点(即 )。由此从图像 的空间 平面得到模糊特性 平面。然后,基于此模糊集定义了图像的线性模糊度 、二次模糊度 和模糊熵 ,使这三个量取最小值时的交叉点 即为最佳阈值。
文献[21]指出模糊隶属度函数在该算法中的作用仅在于将图像由灰度数据空间转换为模糊空间 ,其函数的形式对增强结果几乎没有影响。这就使我们有理由使用一些形式简单的函数形式。例如国内学者发表的一种模糊阈值方法[22]:
 
隶属度μ(x)表示灰度x具有明亮特性的程度,c为隶属函数窗宽,q对应隶属度为0.5的灰度级。设灰度级 的模糊率为:
  = min{μ(l),1-μ(l)}
则得到整幅图像的模糊率[44] 
 
其中,MN为图像尺寸,L为图像总灰度级, 图像中灰度为 的象素个数。
对应于不同的q值,就可以计算出相应的图像模糊率,选取使得 最小的q值,作为图像分割的最佳阈值即可。


3.9 小结
对于基于点的全局阈值选取方法,除上述主要几种之外还许多,但大多都是以上述基本方法为基础,做出的改进方法或者对算法的优化,如使用递推方法以降低算法复杂性。
例如在文献[42]中,提出一种使目标和背景差距最大的阈值求取方法,类似于最大类间方差阈值法。是它的一种简化算法。
又如1984年Dunn等人[23]提出了均匀化误差阈值选取方法,这种方法实质上是要使将背景点误分为目标点的概率等于将目标点误分为背景点的概率。类似于最小误差阈值法。
近年来有一些新的研究手段被引入到阈值选取中。比如人工智能,在文献[24] 中,描述了如何用人工智能的方法,寻找直方图的谷底点,作为全局阈值分割。其它如神经网络,数学形态学[39][46],小波分析与变换[40]等等。
总的来说,基于点的全局阈值算法,与其它几大类方法相比,算法时间复杂度较低,易于实现,适合应用于在线实时图像处理系统。由于我的研究方向为机器视觉,所作的项目要求算法具有良好的实时性,因此针对基于点的全局阈值方法,阅读了较多的文献,在综述里叙述也相对比较详细。


4 基于区域的全局阈值选取方法
对一幅图像而言,不同的区域,比如说目标区域或背景区域,同一区域内的象素,在位置和灰度级上同时具有较强的一致性和相关性。
而在上述基于点的全局阈值选取方法中,有一个共同的弊病,那就是它们实际上只考虑了直方图提供的灰度级信息,而忽略了图像的空间位置细节,其结果就是它们对于最佳阈值并不是反映在直方图的谷点的情况会束手无策,不幸我们通常遇到的很多图像恰恰是这种情况。另一方面,完全不同的两幅图片却可以有相同的直方图,所以即使对于峰谷明显的情况,这些方法也不能保证你得到合理的阈值。于是,人们又提出了很多基于空间信息的阈值化方法。
可以说,局域区域的全局阈值选取方法,是基于点的方法,再加上考虑点领域内象素相关性质组合而成,所以某些方法常称为“二维xxx方法”。由于考虑了象素领域的相关性质,因此对噪声有一定抑止作用[41]。
4.1 二维熵阈值分割方法[25]
使用灰度级-局域平均灰度级形成的二维灰度直方图[43]进行阈值选取,这样就得到二维熵阈值化方法。
 
(二维灰度直方图: 灰度-领域平均灰度)
如图,在0区和1区,象素的灰度值与领域平均灰度值接近,说明一致性和相关性较强,应该大致属于目标或背景区域;2区和3区一致性和相关性较弱,可以理解为噪声或边界部分。二维熵阈值分割,就是选择(S,T)对,使得目标类和背景类的后验熵最大。(具体方法是一维熵阈值分割的推广,可参见上一节)
Abutaleb[26],和Pal]结合Kapur]和Kirby的方法,分别提出了各自的二维熵阈值化方法,其准则函数都是使目标熵和背景熵之和最大化。Brink[27]的方法则是使这两者中的较小者最大化,该方法的计算复杂度为 ,后来有人改进为递推快速算法将时间复杂度降为 (其中 为最大灰度级数)。 


4.2  简单统计法
Kittler等人[28],[29]提出一种基于简单的图像统计的阈值选取方法。使用这种方法,阈值可以直接计算得到,从而避免了分析灰度直方图,也不涉及准则函数的优化。该方法的计算公式为
                         (19)
其中, 
    
    
因为e(x,y)表征了点(x,y)领域的性质,因此本方法也属于基于区域的全局阈值法。


4.3  直方图变化法
从理论上说,直方图的谷底是非常理想的分割阈值,然后在实际应用中,图像常常受到噪声等的影响而使其直方图上原本分离的峰之间的谷底被填充,或者目标和背景的峰相距很近或者大小差不多,要检测他们的谷底就很难了。
在上一节基于点的全局阈值方法中,我们知道直方图凹面分析法的弊病是容易受到噪声干扰,对不同类型的图像,表现出不同的分割效果。往往容易得到假的谷底。这是由于原始的直方图是离散的,而且含噪声,没有考虑利用象素领域性质。
而直方图变化法,就是利用一些象素领域的局部性质变换原始的直方图为一个新的直方图。这个新的直方图与原始直方图相比,或者峰之间的谷底更深,或者谷转变成峰从而更易于检测。这里的象素领域局部性质,在很多方法中经常用的是象素的梯度值。
 例如,由于目标区的象素具有一定的一致性和相关性,因此梯度值应该较小,背景区也类似。而边界区域或者噪声,就具有较大的梯度值。最简单的直方图变换方法,就是根据梯度值加权,梯度值小的象素权加大,梯度值大的象素权减小。这样,就可以使直方图的双峰更加突起,谷底更加凹陷。


4.4 其它基于区域的全局阈值法
松弛法利用邻域约束条件迭代改进线性方程系统的收敛特性,当用于图像阈值化时其思想是:首先根据灰度级按概率将像素分为“亮”和“暗”两类,然后按照领域像素的概率调整每个像素的概率,调整过程迭代进行,使得属于亮(暗)区域的像素“亮(暗)”的概率变得更大。
其它还有许多方法利用灰度值和梯度值散射图,或者利用灰度值和平均灰度值散射图。


5 局部阈值法和多阈值法


5.1 局部阈值(动态阈值)
当图像中有如下一些情况:有阴影,照度不均匀,各处的对比度不同,突发噪声,背景灰度变化等,如果只用一个固定的全局阈值对整幅图像进行分割,则由于不能兼顾图像各处的情况而使分割效果受到影响。有一种解决办法就是用与象素位置相关的一组阈值(即阈值使坐标的函数)来对图像各部分分别进行分割。这种与坐标相关的阈值也叫动态阈值,此方法也叫变化阈值法,或自适应阈值法。这类算法的时间复杂性可空间复杂性比较大,但是抗噪能力强,对一些用全局阈值不易分割的图像有较好的效果。
例如,一幅照度不均(左边亮右边暗)的原始图像为:

如果只选择一个全局阈值进行分割,那么将出现下面两种情况,都不能得到满意的效果。


                
(阈值低,对亮区效果好,则暗区差)          (阈值高,对暗区效果好,则亮区差)


若使用局部阈值,则可分别在亮区和暗区选择不同的阈值,使得整体分割效果较为理性。


 
(按两个区域取局部阈值的分割结果)
进一步,若每个数字都用不同的局部阈值,则可达到更理想的分割效果。

/************************************************************************/
/* 全局阈值分割  自动求取阈值        */
/************************************************************************/
//自动求取阈值,增加对场景的适应性
//只需求取一次,之后就可以一直使用
#include<cv.h>
#include <highgui.h>
#include <iostream>
#include <math.h>
using namespace std;
int main(){
	IplImage * image,* image2;
	image = cvLoadImage("E:\\image\\dowels.tif",0);
	cvNamedWindow("image",1);
	cvShowImage("image",image);
	image2 = cvCreateImage(cvSize(image->width,image->height),image->depth,1);
	double T = 0;
	double dT0 = 1.0;//阈值求取结束标志
	double dT = 255.0;

	//求取平均灰度,作为阈值T的初始值T0
	  int i, j;
     double T0 = 0,T1 = 0,T2 = 0;//初始阈值
	 int count1,count2;
	 unsigned char * ptr,*dst;
	 for (i = 0 ; i< image->height ; i++)
	 {
		 for (j =0 ; j < image->width;j++)
		 {
			 ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
			 T0 += ((double)(*ptr))/image->width/image->height;
		 }
	 }
	 cout<<"T0:     "<<T0<<endl;
	 T = (int)(T0 + 0.5); 
	 //计算T两侧的灰度平均值,然后把二者的均值赋值给T
	 while (dT > dT0)
	 {

		 T1 = 0;
		 T2 = 0;
		 count1 = 0;
		 count2 = 0;
		 for (i = 0 ; i< image->height ; i++)
		 {
			 for (j =0 ; j < image->width;j++)
			 {
				 ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
				if (*ptr > T)
				{
					T1 += ((double)(*ptr))/image->width/image->height;
					count1++;
				} 
				else if(*ptr < T)
				{
					T2 +=  ((double)(*ptr))/image->width/image->height;
					count2++;
				}
			 }
		 }

		 T1 = T1*image->width*image->height/count1;
		 T2 = T2*image->width*image->height/count2;
		 dT = fabs(T - (T1 + T2)/2);
		  
		 cout<<"T1"<<T1<<endl;
		 cout<<"T2"<<T2<<endl;
		 cout<<"dT  " << dT<<endl;
		 T = (T1 + T2)/2;
		 cout<<"T:     "<<T<<endl;
	 }
	 


	 //根据求取的阈值进行分割
	 for (i = 0 ; i< image2->height ; i++)
	 {
		 for (j =0 ; j < image2->width;j++)
		 {
			 ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
			  dst = (unsigned char *)image2->imageData+i*image2->widthStep+j;
			if (*ptr > T)
			{
				*dst = 255;
			} 
			else
			{
				*dst =0;
			}
		 }
	 }

	 cvNamedWindow("image2",1);
	 cvShowImage("image2",image2);
	 cvSaveImage("E:\\image\\dowels2.tif",image2);
	 cvWaitKey(0);
	 return 0;
}




5.1.1 阈值插值法
 首先将图像分解成系列子图,由于子图相对原图很小,因此受阴影或对比度空间变化等带来的问题的影响会比较小。然后对每个子图计算一个局部阈值(此时的阈值可用任何一种固定阈值选取方法)。通过对这些子图所得到的阈值进行插值,就可以得到对原图中每个象素进行分割所需要的合理阈值。这里对应每个象素的阈值合起来构成的一个曲面,叫做阈值曲面。


5.1.2 水线阈值算法
水线(也称分水岭或流域,watershed)阈值算法可以看成是一种特殊的自适应迭代阈值方法,它的基本思想是:初始时,使用一个较大的阈值将两个目标分开,但目标间的间隙很大;在减小阈值的过程中,两个目标的边界会相向扩张,它们接触前所保留的最后像素集合就给出了目标间的最终边界,此时也就得到了阈值。


5.1.3 其它的局部阈值法
文献[30]提出了一种基于阈值曲面的二维遗传算法。遗传算法是基于进化论中自然选择机理的、并行的、统计的随机化搜索方法,所以在图像处理中常用来确定分割阈值。
 文献[31] [32]中提出一种基于局部梯度最大值的插值方法。首先平滑图像,并求得具有局部梯度最大值的像素点,然后利用这些像素点的位置和灰度在图像上内插,得到灰度级阈值表面。
除此之外,典型的局部阈值方法还有White和Rohrer[33]的加权移动平均阈值方法,Perez和Gonzalez[34]的适用于非均匀照射下图像的局部阈值方法以及Shio[35]的与照射无关的对比度度量阈值方法等。总的来说,这类算法的时间和空间复杂度都较大,但是抗噪能力强,对一些使用全局阈值法不宜分割的图像具有较好的效果。


5.2 多阈值法
很显然,如果图像中含有占据不同灰度级区域的几个目标,则需要使用多个阈值才能将它们分开。其实多域值分割,可以看作单阈值分割的推广,前面讨论的大部分阈值化技术,诸如Otsu的最大类间方差法, Kapur的最大熵方法、矩量保持法和最小误差法等等都可以推广到多阈值的情形。以下介绍另外几种多阈值方法。
5.2.1 基于小波的多域值方法。
小波变换的多分辨率分析能力也可以用于直方图分析[36],一种基于直方图分析的多阈值选取方法思路如下:首先在粗分辨率下,根据直方图中独立峰的个数确定分割区域的类数,这里要求独立峰应该满足三个条件:(1)具有一定的灰度范围;(2)具有一定的峰下面积;(3)具有一定的峰谷差。然后,在相邻峰之间确定最佳阈值,这一步可以利用多分辨的层次结构进行。首先在最低分辨率一层进行,然后逐渐向高层推进,直到最高分辨率。可以基于最小距离判据对在最低层选取的所有阈值逐层跟踪,最后以最高分辨率层的阈值为最佳阈值。
5.2.2 基于边界点的递归多域值方法。
这是一种递归的多阈值方法。首先,将象素点分为边界点和非边界点两类,边界点再根据它们的邻域的亮度分为较亮的边界点和较暗的边界点两类,然后用这两类边界点分别作直方图,取两个直方图中的最高峰多对应的灰度级作为阈值。接下去,再分别对灰度级高于和低于此阈值的像素点递归的使用这一方法,直至得到预定的阈值数。


5.2.3 均衡对比度递归多域值方法。
首先,对每一个可能阈值计算它对应于它的平均对比度
                                         
其中, 是阈值为 时图像总的对比度, 是阈值 检测到的边界点的数目。然后,选择 的直方图上的峰值所对应的灰度级为最佳阈值。对于多阈值情形,首先用这种方法确定一个初始阈值,接着,去掉初始阈值检测到的边界点的贡献再做一次 的直方图,并依据新的直方图选择下一个阈值。这一过程可以这样一直进行下去,直到任何阈值的最大平均对比度小于某个给定的限制为止。


6 阈值化算法评价简介
尽管人们在图像分割方面做了许多研究工作,但由于尚无通用的分割理论,现已提出的分割算法大都是针对具体问题的,并没有一种适合于所有图像的通用的分割算法。另一方面,给定一个实际图像分割问题要选择合用的分割算法也还没有标准的方法。为解决这些问题需要研究对图像分割的评价问题。分割评价是改进和提高现有算法性能、改善分割质量和指导新算法研究的重要手段。
 然而,如同所有的图像分割方法一样,阈值化结果的评价是一个比较困难的问题。事实上对图像分割本身还缺乏比较系统的精确的研究,因此对其评价则更差一些。人们先后已经提出了几十个评价准则。这些准则中又有定性的,也有定量的;有分析算法的,也有检测实验结果的,文献[37]将它们大致分为13类。
文献[4] 中选择摄影师、建筑物和模特三幅图像作为标准图像,并采用趋于一致性度量和形状参数对几种常用的全局阈值方法的分割结果进行了评价。结果表明对于这三幅图像,如果希望得到的二值图像比较均匀且目标的形状较好,推荐使用最大熵方法、矩量保持方法和最大类间方差法。
文献[38] 中以磁盘及鹤模型作标准图像,在噪声条件下用错分概率、形状和均匀性度量作为标准评估了五种常见的整体阈值选取方法的性能。这五种方法是四元树方法、矩量保持法、最大类间方差法、最大熵方法和简单统计法。结果表明各种方法的性能不仅与所处理的图像有关,而且也和所选用的准则有关。该文献也指出,对于一般实时应用来说,可以选择最大类间方差方法和简单统计法。
最后,评价的目的是为了能指导、改进和提高分割,如何把评价和分割应用联系起来尚有许多工作要做。一个可能的方法是结合人工智能技术,建立分割专家系统[45],以有效的利用评价结果进行归纳推理,从而把对图像的分割由目前比较盲目的试验阶段推进到系统地实现的阶段。
2019-07-06 21:20:18 weixin_43124455 阅读数 348

otsu 双阈值算法 图像分割

理论

网上有很多博客都写了理论,但是很少有代码 就算有也是opencv c++的代码
所以贴一个自己写的matlab版本
otsu理论的话, 看冈萨雷斯的《数字图像处理》 第10.3.6节 多阈值处理
看完书,按照书上的理论来写的

MATLAB代码

otsu 双阈值函数:



function [t1,t2]=DoubleOtsuThresh(img)
%
%  Otsu 双阈值求解
%  输入 图像img,输出 最优阈值t1和t2(归一化,范围在[0,1])
%
%


BinsNum = 256;  
hist = imhist(img,BinsNum);   
p = hist / sum(hist);          % 直方图的概率密度函数
mG= sum(p .* (1:BinsNum)');     % 全局均值

P1 = cumsum(p);                 % 概率分布
m1 = cumsum(p .* (1:BinsNum)')./P1; % 256*1  每个阈值的前景平均灰度

% 根据算法理论,从k2+1累加到L-1,可以先倒着累加再翻转回来
P3= cumsum(flip(p));          
m3 = cumsum(flip(p) .* flip(1:BinsNum)')./P3;
P3=flip(P3);
P3=[P3(2:end) ;0];             %P3的索引用k2,则P3(1)实际上是从p(2)加到p(256),所以要移动一个
                               %移动后也用不到P3(1)这个值,因为k2>k1>1
m3=flip(m3);
m3=[m3(2:end) ;0];

% m2 P2为 k1*k2的索引矩阵,为 256*256 的上三角矩阵 因为k1<k2
m2=zeros(BinsNum,BinsNum);    
P2=zeros(BinsNum,BinsNum);


for k1=1:BinsNum-2
    for k2=k1+1:BinsNum-1
        P2(k1,k2)=1-P1(k1)-P3(k2);
        m2(k1,k2) = sum( (k1+1:k2)' .* p(k1+1:k2) )./ P2(k1,k2) ;        
    end
end


% 遍历k1,k2各种组合  求方差
variance=zeros(BinsNum,BinsNum);
for k1=1:BinsNum-2
    for k2=k1+1:BinsNum-1
        
    % variance 为256*256 的上三角矩阵 因为l1<k2
    variance(k1,k2) = P1(k1)*(m1(k1)-mG)^2 + ...
                      P2(k1,k2)*(m2(k1,k2)-mG)^2 + ...
                      P3(k2)*(m3(k2)-mG)^2;
   
    end
end

% 最大方差点即为最优阈值
[~,index] = max(variance(:));
[index_row,index_col] = ind2sub(size(variance),index);

%灰度值为[0,255] 因此需要-1
t1= (index_row-1)./(BinsNum-1);
t2= (index_col-1)./(BinsNum-1);

end

根据上面求得的两个阈值,把图分成0 127 255三个灰度值显示的函数:


function out=img2gray(img,t1,t2)
 [M,N]=size(img);
 out=zeros(M,N);
  
 for m=1:M
     for n=1:N
         if img(m,n)<t1
             out(m,n)=0;
         else
             if( img(m,n)>t1 && img(m,n)<t2)
                 out(m,n)=127;  %灰度级可调
             else
                out(m,n)=255;	%灰度级可调
             end
         end
     end
 end
end

测试函数

%%B 彩色图
 imgB =imread('window.jpg');  
 imgB=rgb2gray(imgB);
 imgB=im2double(imgB);


%% 先高斯滤波平滑图B,去高频,再用阈值处理 
 gaussH=fspecial('gaussian',[5 5],10);
 smoothB=imfilter(imgB,gaussH);

 [t1B,t2B] =DoubleOtsuThresh(smoothB);    

% 如果imgB是unit8类型,则灰度范围在[0,255] 要把算法得到的两个阈值乘255
% 如果imgB是[0,1]范围的double类型则不用乘
%  T1B=t1B*255;   
%  T2B=t2B*255;

 T1B=t1B;  
 T2B=t2B;

J9=img2gray(smoothB,T1B,T2B); %转为灰度图

 figure 
 imshow(J9,[]);title('平滑+双阈值分割'); % imshow里面J9后面的[]相当于直方图均衡效果,让灰度值均匀分布,可以去掉对比一下

图B平滑滤波后,右边是灰度直方图:
在这里插入图片描述
双阈值处理:
在这里插入图片描述


%%
% —— 图D —— %


  imgD=imread('ladder.jpg');
%   imgD=rgb2gray(imgD);
 [t1D,t2D] =DoubleOtsuThresh(imgD);    

 T1D=t1D*255;
 T2D=t2D*255;
  % 分为 0 127 255 显示
  
J8=img2gray(imgD,T1D,T2D);

 
 figure 
 subplot 121
 imshow(imgD);title('原图D');

subplot 122
 imshow(J8,[]);title('双阈值分割');
 

处理结果:
在这里插入图片描述


%% 用书上的图测试算法,输出的两个阈值和书上结果一样 80177

imgC=imread('iceberg.tif');  % 用书上例子的图片测试,验证算法正确与否
 

 [t1,t2] =DoubleOtsuThresh(imgC); 
 T1=t1*255;
 T2=t2*255;
J7=img2gray(imgC,T1,T2);
 figure
 subplot 121
 imshow(imgC);title('原图C');
 subplot 122
 imshow(J7,[]);title('双阈值');
 disp(['     otsu双阈值为:   ' ,'t1=',num2str(T1),'     t2=',num2str(T2)])

书上的例子: 在这里插入图片描述

The end .
( markdown如何支持matlab的代码片啊 ,为什么matlab这个key没有效果,只能沦落到用js这个key来高亮代码片)

没有更多推荐了,返回首页