# 大津阈值分割（OSTU）

## 一、算法原理

p i = n i / N ， ∑ i = 0 L − 1 p i = 1 p_i = n_i /N， \sum_{i=0}^{L-1}p_i=1 p i ​=n i ​/N ，i =0 ∑L −1 ​p i ​=1

w 0 = P r ( C 0 ) = ∑ i = 0 k p i w_0=Pr(C_0)=\sum_{i=0}^{k}p_i w 0 ​=P r (C 0 ​)=i =0 ∑k ​p i ​

w 1 = P r ( C 1 ) = ∑ i = k + 1 L − 1 p i = 1 − w 0 w_1=Pr(C_1)=\sum_{i=k+1}^{L-1}p_i = 1-w_0 w 1 ​=P r (C 1 ​)=i =k +1 ∑L −1 ​p i ​=1 −w 0 ​

u 0 = ∑ i = 0 k i P r ( i ∣ C 0 ) = ∑ i = 0 k i p i / w 0 u_0=\sum_{i=0}^{k}iPr(i|C_0)=\sum_{i=0}^{k}ip_i/w_0 u 0 ​=i =0 ∑k ​i P r (i ∣C 0 ​)=i =0 ∑k ​i p i ​/w 0 ​

u 1 = ∑ i = k + 1 L − 1 i P r ( i ∣ C 1 ) = ∑ i = k + 1 L − 1 i p i / w 1 u_1=\sum_{i=k+1}^{L-1}iPr(i|C_1)=\sum_{i=k+1}^{L-1}ip_i/w_1 u 1 ​=i =k +1 ∑L −1 ​i P r (i ∣C 1 ​)=i =k +1 ∑L −1 ​i p i ​/w 1 ​

u T = ∑ i = 0 L − 1 i p i u_T=\sum_{i=0}^{L-1}ip_i u T ​=i =0 ∑L −1 ​i p i ​

w o u 0 + w 1 u 1 = u T , w 0 + w 1 = 1 w_ou_0+w_1u_1=u_T,w_0+w_1=1 w o ​u 0 ​+w 1 ​u 1 ​=u T ​,w 0 ​+w 1 ​=1
C_0和C_1类的方差可由下式求得：
σ 0 2 = ∑ i = 0 t ( i − u 0 ) 2 p i / w 0 \sigma_0^2=\sum_{i=0}^{t}(i-u_0)^2p_i/w_0 σ0 2 ​=i =0 ∑t ​(i −u 0 ​)2 p i ​/w 0 ​

σ 1 2 = ∑ i = t + 1 L − 1 ( i − u 1 ) 2 p i / w 1 \sigma_1^2=\sum_{i=t+1}^{L-1}(i-u_1)^2p_i/w_1 σ1 2 ​=i =t +1 ∑L −1 ​(i −u 1 ​)2 p i ​/w 1 ​

σ w 2 = w 0 σ 0 2 + w 1 σ 1 2 \sigma_w^2=w_0\sigma_0^2+w_1\sigma_1^2 σw 2 ​=w 0 ​σ0 2 ​+w 1 ​σ1 2 ​

σ B 2 = w 0 ( u 0 − u T ) 2 + w 1 ( u 1 − u T ) 2 = w 0 w 1 ( u 1 − u 0 ) 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 σB 2 ​=w 0 ​(u 0 ​−u T ​)2 +w 1 ​(u 1 ​−u T ​)2 =w 0 ​w 1 ​(u 1 ​−u 0 ​)2

σ T 2 = σ B 2 + σ w 2 \sigma_T^2=\sigma_B^2+\sigma_w^2 σT 2 ​=σB 2 ​+σw 2 ​

## 二、类间方差推导

σ B 2 = w 0 ( μ 0 − μ T ) 2 + w 1 ( μ 1 − μ T ) 2 = w 0 [ μ 0 − w 0 μ 0 − w 1 μ 1 ] 2 + w 1 [ μ 1 − w 0 μ 0 − w 1 μ 1 ] 2 = w 0 [ ( 1 − w 0 ) μ 0 − w 1 μ 1 ] 2 + w 1 [ ( 1 − w 1 ) μ 1 − w 0 μ 0 ] 2 = w 0 [ w 1 μ 0 − w 1 μ 1 ] 2 + w 1 [ w 0 μ 1 − w 0 μ 0 ] 2 = w 0 w 1 ( μ 0 − μ 1 ) 2 + w 0 w 1 ( μ 1 − μ 0 ) 2 = ( μ 1 − μ 0 ) 2 w 0 w 1 ( w 0 + w 1 ) = w 0 w 1 ( μ 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 σB ​2 =w 0 ​(μ0 ​−μT ​)2 +w 1 ​(μ1 ​−μT ​)2 =w 0 ​[μ0 ​−w 0 ​μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[μ1 ​−w 0 ​μ0 ​−w 1 ​μ1 ​]2 =w 0 ​[(1 −w 0 ​)μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[(1 −w 1 ​)μ1 ​−w 0 ​μ0 ​]2 =w 0 ​[w 1 ​μ0 ​−w 1 ​μ1 ​]2 +w 1 ​[w 0 ​μ1 ​−w 0 ​μ0 ​]2 =w 0 ​w 1 ​(μ0 ​−μ1 ​)2 +w 0 ​w 1 ​(μ1 ​−μ0 ​)2 =(μ1 ​−μ0 ​)2 w 0 ​w 1 ​(w 0 ​+w 1 ​)=w 0 ​w 1 ​(μ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++){

float class_0_prob = 0;

float class_1_prob = 0;
for(int i = 0; i  k; i++){
class_0_prob += pixelProb[i];
}
class_1_prob = 1 - class_0_prob;

float class_0_average = 0;

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));
}

float class_0_variance = 0;

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;

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;
else
tempImage.at<uchar>(i, j) = 255;
}
}


## 四、效果对比

Original: https://blog.csdn.net/qq_42780025/article/details/113000940
Author: nachr
Title: 大津阈值分割（OSTU）

(0)