大津演算法（Otsu）是一種自動影像二值化方法，通過分析影像的灰度分佈，自動找到最適合的閾值，突顯出影像中的目標特徵。這個演算法在處理影像時非常有用，因為使用這個演算法可以幫助我們快速且準確的將影像二值化。
期望值是概率論和統計學中的一個重要概念，用來衡量隨機變數的平均值u或中心趨勢。它是在概率分佈下所有可能的值乘以其對應的概率後的加權平均。
在數學上，對於離散型隨機變數 X，期望值E[X] 的計算公式如下：
當我們用一個六面骰子（六面的正方體，每個面上標有 1 到 6 的數字）投擲時，我們可以計算這個隨機事件的期望值。
假設我們想要計算骰子的點數的期望值。每個點數有等可能性（1/6）出現，因為骰子是均勻的。
因此，期望值計算如下：
在這個例子中，骰子的點數的期望值是 3.5。這意味著，如果我們丟了無限次的骰子，平均來說，我們可以期望得到接近 3.5 的點數。
標準差是一種統計量，用來衡量一組數據的散佈程度或變異程度。當數據中的值之間差異較大時，標準差會較大；而當數據的值彼此間的差異較小時，標準差會較小。
標準差的數學公式如下：
我們可以通過觀察下圖來理解這一點，這是一張平均值為0的常態分佈圖，當常態分佈的標準差（用符號σ表示）越小，這代表數據相對集中。換句話說，大部分的數據點都會集中在平均值附近，較少的數據點會偏離平均值。反之，當常態分佈的標準差越大，這代表數據相對分散，數據點會廣泛地分佈在平均值周圍。
變異數（Variance）是用來衡量一組數據的散佈程度或變異程度。它代表著數據點與其平均值之間的差異程度的平均值。變異數可以幫助我們瞭解數據的分佈情況，以及數據點與平均值之間的偏離程度。
變異數可以透過將標準差平方後得到。如果我們已經知道一組數據的標準差，我們可以通過對其進行平方運算，得到數據的變異數。
這種關係可以用以下的數學公式表示：
在大津演算法中，我們的目標是找到一個最適合的閾值，以便將影像的像素分成前景和背景兩個類別。大津演算法告訴我們，我們希望讓同一類的數據更加相似，而不同類之間的數據更加不同。
在介紹大津演算法是怎麼實現之前，我們需要搞懂一些名詞：
以下是大津演算法中用來計算類間變異數的公式：
可以看到，因為σ 平方是整張圖片的變異數是常數，當類內變異數最小時，類間變異數會最大。只要將兩個類之間的差異(類間變異數)最大化，這兩個類分別的內部差異(類內變異數)就自動的被最小化。換句話說，我們不需要同時追求最大化類之間的差異並最小化內部差異。
大津演算法的原理是搜尋所有閾值t，找出0~255中出現最大類間變異數所在的閾值為多少，實現使同一類的數據更加相似，而不同類之間的數據更加不同的效果。
為了找到類間變異數最大值，我們需要透過上方公式求出類間變異數。
我們需要先算出這張圖片灰階值的平均，假設我們今天已經將圖片的直方圖求出機率質量函數H(X=xi)，H(X=xi)代表整張圖片灰階值為xi的像素個數除以圖片總像素(或是灰階值xi在這張圖片出現的機率)，最大為1。
這個直方圖的期望值(或平均值)可以用數學公式表示：
背景的權重 w1(t)表示的是在二值化閾值 t情況下，像素出現灰階值小於等於 t 的總和機率。反之，前景的權重w2(t) 則是像素出現灰階值大於二值化閾值 t的總和機率。因為兩個權重的總和為1，我們可以藉由w1(t)推出w2(t)的機率。這兩個函數的性質適合使用CDF函數來表示累加的範圍。
接下來，我們要求出背景內的平均值和前景內的平均值，我們可以藉由μ1(t)期望值來推出μ2(t)期望值，以減少不必要的迴圈運算
大津演算法的整體步驟如下：
這個程式碼提供了兩種運算方式，一個是使用OpenCV內建的函式，另一個是使用自己寫的大津演算法，兩者效果一樣。你可以透過設定USE_OPENCV 成1或0來決定你要使用前者還是後者的實現方式。
#define USE_OPENCV 1
使用OpenCV 函式庫中的 cv::threshold 函數，用於對灰階影像進行二值化。
double threshold = cv::threshold(grayImage, binaryImage, 0,255, cv::THRESH_OTSU);
printf("threshold:%d",(int)threshold);
grayImage：這是輸入的灰階影像，表示你要對這張影像進行二值化。binaryImage：這是輸出的二值化影像，函數會將二值化後的影像存儲在這個變數中。0：這是二值化的閾值，因為我們使用大津演算法，不需要提供閾值，所以輸入0。255：這是閾值的上界，也就是將所有高於閾值的像素設為前景（白色）。cv::THRESH_OTSU：這是指定的二值化方法，使用了大津演算法來自動找到適合的閾值，以最大化類間變異數。cv::threshold 函數回傳otsu演算法的結果。calcHist 函數用來計算一張影像的灰階直方圖。它遍歷影像的每個像素，根據像素的灰階值將對應的直方圖的計數值增加。
接下來使用otsu 函式實現了大津演算法，尋找閾值。
cv::Mat hist=calcHist(grayImage);
uchar threshold=otsu(hist);
printf("threshold:%d",threshold);
cv::threshold(grayImage, binaryImage, threshold, 255, cv::THRESH_BINARY);
otsu 函式實現了大津演算法，具體實現的步驟如下：
首先將直方圖轉換成機率質量函數 hist_pmf，並除以圖片總像素進行正規化，以便進行計算。
接著計算整張影像的平均灰階值 u。
在迴圈中，從閾值 0 到 255 逐個測試。對每個閾值 t，計算出在閾值左側的背景權重 w1，並計算出背景的平均灰階值 u1。同時計算在閾值右側的前景權重 w2，以及前景的平均灰階值 u2。然後計算出類間變異數 sigma_b_2。
每次計算出的 sigma_b_2 與先前的最大值 max_sigma_b_2 進行比較，如果當前的 sigma_b_2 更大，則更新最大值並記錄對應的閾值 t。
最終，函式返回找到的最適合的閾值。
// 使用OTSU方法計算閾值
uchar otsu(cv::Mat hist)
{
cv::Mat hist_pmf;
hist.convertTo(hist_pmf, CV_32F);
hist_pmf /= cv::sum(hist)[0];
float u = 0.0f;
for (int i = 0; i < 256; i++)
{
u += i * hist_pmf.at<float>(i);
}
float w1 = 0.0f;
float sum1 = 0.0f;
float max_sigma_b_2 = -1.0f;
uchar threshold = 0;
for (int t = 0; t < 256; t++)
{
//累加背景權重
w1 += hist_pmf.at<float>(t);
if (w1 == 0.0f)
continue;
//累加前景權重
float w2 = 1.0f - w1;
if (w2 == 0.0f)
continue;
sum1 += t * hist_pmf.at<float>(t);
//計算背景平均值
float u1 = sum1 / w1;
float sum2 = u - sum1;
//計算前景平均值
float u2 = sum2 / w2;
//計算類間變異數
float sigma_b_2 = w1 * w2 * (u1 - u2) * (u1 - u2);
if (sigma_b_2 > max_sigma_b_2)
{
max_sigma_b_2 = sigma_b_2;
threshold = t;
}
}
return threshold;
}
#include <iostream>
#include <inttypes.h>
#include <opencv2/opencv.hpp>
#include <opencv2/core/utils/logger.hpp>
#define USE_OPENCV 1
using namespace std;
// 函數原型
uchar otsu(cv::Mat hist);
cv::Mat calcHist(cv::Mat image);
int main()
{
cv::utils::logging::setLogLevel(cv::utils::logging::LOG_LEVEL_ERROR);
// 讀取灰度影像
cv::Mat grayImage = cv::imread("C:\\Users\\vince\\Downloads\\Lenna.png", cv::IMREAD_GRAYSCALE);
cv::Mat binaryImage;
#if USE_OPENCV
// 使用OpenCV的OTSU閾值化
double threshold = cv::threshold(grayImage, binaryImage, 0, 255, cv::THRESH_OTSU);
printf("OTSU threshold:%d", (int)threshold);
#else
// 計算影像的直方圖
cv::Mat hist = calcHist(grayImage);
// 使用OTSU方法計算閾值
uchar threshold = otsu(hist);
printf("OTSU threshold:%d", threshold);
// 根據計算得到的閾值進行二值化
cv::threshold(grayImage, binaryImage, threshold, 255, cv::THRESH_BINARY);
#endif
// 顯示二值化影像
cv::imshow("Binary Image", binaryImage);
cv::waitKey(0);
return 0;
}
// 計算灰度影像的直方圖
cv::Mat calcHist(cv::Mat image)
{
cv::Mat hist = cv::Mat::zeros(cv::Size(1, 256), CV_32SC1);
for (int y = 0; y < image.rows; y++)
{
for (int x = 0; x < image.cols; x++)
{
uchar value = image.at<uchar>(y, x);
hist.at<int>(value)++;
}
}
return hist;
}
// 使用OTSU方法計算閾值
uchar otsu(cv::Mat hist)
{
cv::Mat hist_pmf;
hist.convertTo(hist_pmf, CV_32F);
hist_pmf /= cv::sum(hist)[0];
float u = 0.0f;
for (int i = 0; i < 256; i++)
{
u += i * hist_pmf.at<float>(i);
}
float w1 = 0.0f;
float sum1 = 0.0f;
float max_sigma_b_2 = -1.0f;
uchar threshold = 0;
for (int t = 0; t < 256; t++)
{
//累加背景權重
w1 += hist_pmf.at<float>(t);
if (w1 == 0.0f)
continue;
//累加前景權重
float w2 = 1.0f - w1;
if (w2 == 0.0f)
continue;
sum1 += t * hist_pmf.at<float>(t);
//計算背景平均值
float u1 = sum1 / w1;
float sum2 = u - sum1;
//計算前景平均值
float u2 = sum2 / w2;
//計算類間變異數
float sigma_b_2 = w1 * w2 * (u1 - u2) * (u1 - u2);
if (sigma_b_2 > max_sigma_b_2)
{
max_sigma_b_2 = sigma_b_2;
threshold = t;
}
}
return threshold;
}
可以看到，即使我們沒有提供閾值，大津演算法幫我們找出了最適合這張圖片的閾值。大津演算法的這個特性使得影像處理變得更加自動化，無需事先了解影像的特性，就能夠獲得最佳的結果。這對於處理大量影像或需要快速處理影像的場景非常有用。
iThome鐵人賽
看影片追技術