大津阈值分割（OSTU）

文章首发于我的个人博客

大津法是一种灰度图像自适应阈值分割算法，是日本学者Ostu于1979年提出，又称类间方差阈值分割法。大津法根据图像的灰度分布将图像分为前景和背景两部分，前景是我们分割出来的部分。前景和背景的分割值就是我们要通过类间方差法求出的阈值。

一、算法原理

设图像灰度级为L，灰度级为i的像素点数为ni，那么直方图分布为：
pi=ni/N，∑i=0L−1pi=1p_i = n_i /N， \sum_{i=0}^{L-1}p_i=1 pi=ni/N，i=0∑L−1pi=1
按灰度级用阈值t划分为两类：C_0=(0,1,…,t)和C_1=(t+1,t+2…,L-1)。因此C_0类和C_1类的概率级均值分别由一下各式给出：
w0=Pr(C0)=∑i=0kpiw_0=Pr(C_0)=\sum_{i=0}^{k}p_i w0=Pr(C0)=i=0∑kpi

w1=Pr(C1)=∑i=k+1L−1pi=1−w0w_1=Pr(C_1)=\sum_{i=k+1}^{L-1}p_i = 1-w_0 w1=Pr(C1)=i=k+1∑L−1pi=1−w0

u0=∑i=0kiPr(i∣C0)=∑i=0kipi/w0u_0=\sum_{i=0}^{k}iPr(i|C_0)=\sum_{i=0}^{k}ip_i/w_0 u0=i=0∑kiPr(i∣C0)=i=0∑kipi/w0

u1=∑i=k+1L−1iPr(i∣C1)=∑i=k+1L−1ipi/w1u_1=\sum_{i=k+1}^{L-1}iPr(i|C_1)=\sum_{i=k+1}^{L-1}ip_i/w_1 u1=i=k+1∑L−1iPr(i∣C1)=i=k+1∑L−1ipi/w1

uT=∑i=0L−1ipiu_T=\sum_{i=0}^{L-1}ip_i uT=i=0∑L−1ipi

可以看出，对任何t值，下式都能成立：
wou0+w1u1=uT,w0+w1=1w_ou_0+w_1u_1=u_T,w_0+w_1=1 wou0+w1u1=uT,w0+w1=1
C_0和C_1类的方差可由下式求得：
σ02=∑i=0t(i−u0)2pi/w0\sigma_0^2=\sum_{i=0}^{t}(i-u_0)^2p_i/w_0 σ02=i=0∑t(i−u0)2pi/w0

σ12=∑i=t+1L−1(i−u1)2pi/w1\sigma_1^2=\sum_{i=t+1}^{L-1}(i-u_1)^2p_i/w_1 σ12=i=t+1∑L−1(i−u1)2pi/w1

由此便可定义，类内方差(within)：
σw2=w0σ02+w1σ12\sigma_w^2=w_0\sigma_0^2+w_1\sigma_1^2 σw2=w0σ02+w1σ12
类间方差(Between):
σB2=w0(u0−uT)2+w1(u1−uT)2=w0w1(u1−u0)2\sigma_B^2=w_0(u_0-u_T)^2+w_1(u_1-u_T)^2=w_0w_1(u_1-u_0)^2 σB2=w0(u0−uT)2+w1(u1−uT)2=w0w1(u1−u0)2
总体方差：
σT2=σB2+σw2\sigma_T^2=\sigma_B^2+\sigma_w^2 σT2=σB2+σw2
何为类间方差？我们将灰度级L按t非为两类，反应在直方图上就是两类灰度。我们要做的就是找到能使两边灰度差的最大的阈值，根据这个阈值对灰度图像进行分割。因此我们要做的就是找到能使类间方差最大的t值。

二、类间方差推导

σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0[μ0−w0μ0−w1μ1]2+w1[μ1−w0μ0−w1μ1]2=w0[(1−w0)μ0−w1μ1]2+w1[(1−w1)μ1−w0μ0]2=w0[w1μ0−w1μ1]2+w1[w0μ1−w0μ0]2=w0w1(μ0−μ1)2+w0w1(μ1−μ0)2=(μ1−μ0)2w0w1(w0+w1)=w0w1(μ1−μ0)2{\sigma_B}^2=w_0(\mu_0-\mu_T)^2+w_1(\mu_1-\mu_T)^2 \\ = w_0[\mu_0-w_0\mu_0-w_1\mu_1]^2+w_1[\mu_1-w_0\mu_0-w_1\mu_1]^2 \\ = w_0[(1-w_0)\mu_0-w_1\mu_1]^2+w_1[(1-w_1)\mu_1-w_0\mu_0]^2 \\ =w_0[w_1\mu_0-w_1\mu_1]^2+w_1[w_0\mu_1-w_0\mu_0]^2 \\ =w_0w_1(\mu_0-\mu_1)^2+w_0w_1(\mu_1-\mu_0)^2 \\ =(\mu_1-\mu_0)^2w_0w_1(w_0+w_1)\\ =w_0w_1(\mu_1-\mu_0)^2 σB2=w0(μ0−μT)2+w1(μ1−μT)2=w0[μ0−w0μ0−w1μ1]2+w1[μ1−w0μ0−w1μ1]2=w0[(1−w0)μ0−w1μ1]2+w1[(1−w1)μ1−w0μ0]2=w0[w1μ0−w1μ1]2+w1[w0μ1−w0μ0]2=w0w1(μ0−μ1)2+w0w1(μ1−μ0)2=(μ1−μ0)2w0w1(w0+w1)=w0w1(μ1−μ0)2

三、c++实现（核心代码）

 auto tempImage = originImage.clone();std::vector<int> countPixel(256, 0);//直方图灰度分布概率std::vector<float> pixelProb;std::vector<int> pixel(256);for(int i = 0; i < 256; i++){pixel[i] = i;}for(int i = 0; i < tempImage.rows; i++){for(int j = 0; j < tempImage.cols; j++){countPixel[tempImage.at<uchar>(i, j)]++;}}//计算直方图分布float sum_pixel = tempImage.rows * tempImage.cols;for(auto iter = countPixel.begin(); iter != countPixel.end(); iter++){pixelProb.push_back(*iter / sum_pixel);}//类内方差std::vector<float> class_within_variance_vec;//类间方差std::vector<float> class_between_variance_vec;for(int k = 0; k < 256; k++){//类0的概率float class_0_prob = 0;//类1的概率float class_1_prob = 0;for(int i = 0; i <= k; i++){class_0_prob += pixelProb[i];}class_1_prob = 1 - class_0_prob;//std::cout<<"class_0_prob:"<<class_0_prob<<" class_1_prob:"<<class_1_prob<<std::endl;//类0（前景）灰度均值float class_0_average = 0;//类1（背景）灰度均值float class_1_average = 0;for(int i = 0; i <= k; i++){class_0_average += i * (pixelProb[i] / (class_0_prob+0.000001));}for(int i = k+1; i <= 255; i++){class_1_average += i * (pixelProb[i] / (class_1_prob+0.000001));}//类0（前景）方差float class_0_variance = 0;//类1（背景）方差float class_1_variance = 0;for(int i = 0; i <= k; i++){class_0_variance += std::pow((i-class_0_average),2)*(pixelProb[i]/(class_0_prob+0.000001));}for(int i = k+1; i <= 255; i++){class_1_variance += std::pow((i-class_1_average),2)*(pixelProb[i]/(class_1_prob+0.000001));}//类内方差float class_within_variance = class_0_prob*std::pow(class_0_variance,2) + class_1_prob*std::pow(class_1_variance,2);class_within_variance_vec.push_back(class_within_variance);//类间方差float class_between_variance = class_0_prob * class_1_prob * std::pow((class_1_average-class_0_average),2);class_between_variance_vec.push_back(class_between_variance);}auto max_value_iter = std::max_element(class_between_variance_vec.begin(), class_between_variance_vec.end());//得到类间方差最大值int k_value = std::distance(class_between_variance_vec.begin(), max_value_iter);std::cout<<"k_value:"<< k_value<<"  max_var:"<<*max_value_iter<<std::endl;/* auto min_value_iter = std::min_element(class_within_variance_vec.begin(), class_within_variance_vec.end());int kk_value = std::distance(class_within_variance_vec.begin(), min_value_iter);std::cout<<"kk_value:"<< kk_value<<"  max_var:"<<*min_value_iter<<std::endl;*/for(int i = 0; i < tempImage.rows; i++){for(int j = 0; j < tempImage.cols; j++){if(tempImage.at<uchar>(i, j) <= k_value)tempImage.at<uchar>(i, j) = 0;elsetempImage.at<uchar>(i, j) = 255;}}

四、效果对比

下图从左到右依次为，原灰度图、OpenCV提供OSTU算法分割结果、自己实现的OSTU分割效果。

下图为调用matlab函数实现的大津法阈值分割效果

必不可少的Lenna。