Метод Оцу

Граничний приклад зображення з використанням алгоритму Оцу.

Оригінальне зображення.

У комп'ютерному баченні та обробці зображень метод Оцу, названий на честь Шаблон:Нп (大津展之 Ōtsu Nobuyuki), є методом, заснованим на кластеризації, для автоматичного обчислення порогового зображення, або зведення сірого зображення до бінарного зображення. Алгоритм передбачає, що зображення містить два класи пікселів, наступної бі-модальної гістограми: пікселі переднього плану і пікселі тла, потім обчислюється оптимальний поріг, що розділяє два класи, так, що їх комбінований діапазон (дисперсія кластера) є мінімальною або рівноцінною (тому що сума попарних квадратичних відстаней постійна), так, що їх міжкластерна дисперсія є максимальною. Отже, метод Оцу є приблизно одновимірним дискретним аналогом дискримінантного аналізу Фішера. Метод Оцу також безпосередньо пов'язаний з Шаблон:Нп.

Розширення вихідного методу до багаторівневого порогового значення називається багатоточковим методом Оцу.

Метод

У методі Оцу ми вичерпно шукаємо поріг, що мінімізує дисперсію всередині класу, яка визначається як зважена сума дисперсій двох класів:

${\sigma }_w^2(t)={\omega }_0(t){\sigma }_0^2(t)+{\omega }_1(t){\sigma }_1^2(t)$

Вага ${\omega }_0$ та ${\omega }_1$ ймовірності двох класів, розділені порогом $t$ і ${\sigma }_0^2$ та ${\sigma }_1^2$ — відхилення цих двох класів.

Клас ймовірностей ${\omega }_{0,1}(t)$ обчислюється з $L$ гістограм:

$\begin{array}{cc}{\omega }_0(t)& =\sum _{i=0}^{t-1}p(i)\\ {\omega }_1(t)& =\sum _{i=t}^{L-1}p(i)\end{array}$

Оцу показує, що мінімізація дисперсії всередині класу збігається з максимізацією міжкласової дисперсії:

$\begin{array}{cc}{\sigma }_b^2(t)& ={\sigma }^2-{\sigma }_w^2(t)={\omega }_0({\mu }_0-{\mu }_T{)}^2+{\omega }_1({\mu }_1-{\mu }_T{)}^2\\ & ={\omega }_0(t){\omega }_1(t){[{\mu }_0(t)-{\mu }_1(t)]}^2\end{array}$

що виражається в термінах ймовірностей класу $\omega$ та класу значень $\mu$.

коли клас ${\mu }_{0,1,T}(t)$ означає:

$\begin{array}{cc}{\mu }_0(t)& =\sum _{i=0}^{t-1}i\frac{p(i)}{{\omega }_0}\\ {\mu }_1(t)& =\sum _{i=t}^{L-1}i\frac{p(i)}{{\omega }_1}\\ {\mu }_T& =\sum _{i=0}^{L-1}ip(i)\end{array}$

Наступні співвідношення можна легко перевірити:

$\begin{array}{cc}{\omega }_0{\mu }_0+{\omega }_1{\mu }_1& ={\mu }_T\\ {\omega }_0+{\omega }_1& =1\end{array}$

Ймовірності класу і засоби класу можуть бути обчислені ітеративно. Ця ідея дає ефективний алгоритм.

Алгоритм

1. Обчисліть гістограму та ймовірність кожного рівня інтенсивності
2. Установіть початкові ${\omega }_i(0)$ та ${\mu }_i(0)$
3. Пройдіть усі можливі порогові значення $t=1,\dots$ максимальна інтенсивність
   1. Оновіть ${\omega }_i$ та ${\mu }_i$
   2. Обчисліть ${\sigma }_b^2(t)$
4. Бажаний поріг відповідає максимуму ${\sigma }_b^2(t)$

MATLAB реалізація

histogramCounts — це 256-елементна гістограма зображення в градаціях сірого, різних рівнях сірого (типові для 8-бітних зображень). Рівень — це поріг зображення.

function level = otsu(histogramCounts)
total = sum(histogramCounts); % '''total''' is the number of pixels in the given image. 
%% OTSU automatic thresholding method
sumB = 0;
wB = 0;
maximum = 0.0;
sum1 = dot( (0:255), histogramCounts); 
for ii=1:256
    wB = wB + histogramCounts(ii);
    wF = total - wB;
    if (wB == 0 || wF == 0)
        continue;
    end
    sumB = sumB +  (ii-1) * histogramCounts(ii);
    mF = (sum1 - sumB) / wF;
    between = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF);
    if ( between >= maximum )
        level = ii;
        maximum = between;
    end
end
end

Matlab має вбудовані функції graythresh() та multithresh() в Image Processing Toolbox, які реалізуються за допомогою методу Оцу та методу Мульти-Оцу, відповідно.

Інший підхід з векторизованним методом (може бути легко перетворений в версію матриці масиву python для обробки GPU)

function  [threshold_otsu] = Thresholding_Otsu( Image)
%Intuition:
%(1)pixels are divided into two groups
%(2)pixels within each group are very similar to each other 
%   Parameters:
%   t : threshold 
%   r : pixel value ranging from 1 to 255
%   q_L, q_H : the number of lower and higher group respectively
%   sigma : group variance
%   miu : group mean
%   Author: Lei Wang
%   Date  : 22/09/2013
%   References : Wikepedia, 
%   for multi children Otsu method, please visit : https://drive.google.com/file/d/0BxbR2jt9XyxteF9fZ0NDQ0dKQkU/view?usp=sharing
%   This is my original work




nbins = 256;
counts = imhist(Image,nbins);
p = counts / sum(counts);




for t = 1 : nbins
   q_L = sum(p(1 : t));
   q_H = sum(p(t + 1 : end));
   miu_L = sum(p(1 : t) .* (1 : t)') / q_L;
   miu_H = sum(p(t + 1 : end) .* (t + 1 : nbins)') / q_H;
   sigma_b(t) = q_L * q_H * (miu_L - miu_H)^2;
end




[~,threshold_otsu] = max(sigma_b(:));
end

Реалізація має невелику надмірність обчислень. Але оскільки метод Оцу виконується швидко, реалізація прийнята і зрозуміла. Хоча в деяких середовищах, оскільки ми використовуємо форму векторизації, обчислення циклів може бути швидше. Цей метод може бути легко перетворений в многопороговий метод з мінімальними мітками для heap—children labels.

Python реалізація

import numpy as np
from collections import Counter

def Thresholding_Otsu(img):
    nbins = 256
    pixel_counts  = Counter(img.ravel())
    counts = np.array([0 for x in range(256)])
    for c in sorted(pixel_counts):
        counts[c] = pixel_counts[c]
    p = counts/sum(counts)
    sigma_b = np.zeros((256,1))
    for t in range(nbins):
        q_L = sum(p[:t]) 
        q_H = sum(p[t:]) 
        if q_L ==0 or q_H == 0:
            continue
            
        miu_L = sum(np.dot(p[:t],np.transpose(np.matrix([i for i in range(t)]) )))/q_L
        miu_H = sum(np.dot(p[t:],np.transpose(np.matrix([i for i in range(t,256)]))))/q_H
        sigma_b[t] = q_L*q_H*(miu_L-miu_H)**2
        
    return np.argmax(sigma_b)

JavaScript реалізація

Варіант 1

Вхідний аргумент total — це кількість пікселів в даному зображенні. Вхідний аргумент histogram являє собою 256-елементну гістограму зображення в градаціях сірого, різних рівнів сірого (типові для 8-бітних зображень). Ця функція виводить поріг для зображення.

function otsu(histogram, total) {
    var sum = 0;
    for (var i = 1; i < histogram.length + 1; ++i) //normally it will be 255 but sometimes we want to change step
        sum += i * histogram[i];
    var sumB = 0;
    var wB = 0;
    var wF = 0;
    var mB;
    var mF;
    var max = 0.0;
    var between = 0.0;
    var threshold1 = 0.0;
    var threshold2 = 0.0;
    for (var i = 0; i < histogram.length; ++i) {
        wB += histogram[i];
        if (wB == 0)
            continue;
        wF = total - wB;
        if (wF == 0)
            break;
        sumB += i * histogram[i];
        mB = sumB / wB;
        mF = (sum - sumB) / wF;
        between = wB * wF * (mB - mF) * (mB - mF);
        if ( between >= max ) {
            threshold1 = i;
            if ( between > max ) {
                threshold2 = i;
            }
            max = between;            
        }
    }
    return ( threshold1 + threshold2 ) / 2.0;
}

Варіант 2

function otsu(histogram, pixelsNumber) {
	var sum = 0
	  , sumB = 0
	  , wB = 0
	  , wF = 0
      , mB
	  , mF
	  , max = 0
	  , between
	  , threshold = 0;
	for (var i = 0; i < 256; ++i) {
      wB += histogram[i];
      if (wB == 0)
        continue;
      wF = pixelsNumber - wB;
      if (wF == 0)
        break;
      sumB += i * histogram[i];
      mB = sumB / wB;
      mF = (sum - sumB) / wF;
      between = wB * wF * Math.pow(mB - mF, 2);
      if (between > max) {
        max = between;
        threshold = i;
      }
    }
    return threshold;
}




// To test: open any image in browser and run code in console
var im = document.getElementsByTagName('img')[0]
  , cnv = document.createElement('canvas')
  , ctx = cnv.getContext('2d');
cnv.width = im.width;
cnv.height = im.height;
ctx.drawImage(im, 0, 0);
var imData = ctx.getImageData(0, 0, cnv.width, cnv.height)
  , histogram = Array(256)
  , i
  , red
  , green
  , blue
  , gray;
for (i = 0; i < 256; ++i)
	histogram[i] = 0;
for (i = 0; i < imData.data.length; i += 4) {
  red = imData.data[i];
  blue = imData.data[i + 1];
  green = imData.data[i + 2];
  // alpha = imData.data[i + 3];
  // https://en.wikipedia.org/wiki/Grayscale
  gray = red * .2126 + green * .07152 + blue * .0722;
  histogram[Math.round(gray)] += 1;
}
var threshold = otsu(histogram, imData.data.length / 4);
console.log("threshold = %s", threshold);
for (i = 0; i < imData.data.length; i += 4) {
  imData.data[i] = imData.data[i + 1] = imData.data[i + 2] =
    imData.data[i] >= threshold ? 255 : 0;
  // opacity 255 = 100%
  imData.data[i + 3] = 255;
}
ctx.putImageData(imData, 0, 0);
document.body.appendChild(cnv);
console.log("finished");

Result.

Обмеження

Метод Оцу демонструє відносно хорошу продуктивність, якщо припустити, що гістограма має бімодальний розподіл і передбачає наявність глибокої і гострої долини між двома піками. Але якщо площа об'єкта мала в порівнянні з фоновою областю, гістограма більше не має бімодального розподілу. І якщо дисперсія об'єкта та інтенсивність фону великі в порівнянні з середньою різницею, або зображення сильно пошкоджено адитивним шумом, гостра долина гістограми рівня сірого погіршується. Тоді можливий неправильний поріг, визначений методом Оцу, призводить до помилки сегментації. (Тут ми визначаємо розмір об'єкта як відношення площі об'єкта до всієї області зображення, а середня різниця — різниця середніх інтенсивностей об'єкта і фону)

Показано, що з результатів експериментів ефективність глобальних методів порога, включаючи метод Оцу, обмежена розміром невеликого об'єкта, малої середньої різниці, великими дисперсіями об'єкта та інтенсивністю фону, великою кількістю доданого шуму і т. д.

Поліпшення

Існує безліч поліпшень, спрямованих на різні обмеження методу Оцу. Один відомий і ефективний спосіб відомий як двовимірний метод Оцу. У цьому підході вивчається значення рівня сірого для кожного пікселя, а також середнє значення його безпосередній околиці, так що результати бінаризації значно поліпшуються, особливо для зображень, спотворених шумом.

На кожному пікселі обчислюється середнє значення сірого рівня околиці. Нехай сірий рівень даного зображення буде розділений на значення $L$ та середній рівень сірого також буде розділений на ті ж самі значення $L$. Потім формується пара: рівень сірого пікселя і середнє значення околиці. Кожна пара належить двовимірному скриньки. Загальна кількість ящиків, очевидно, $L\times L$. Загальна кількість входжень (частота), $f_{ij}$, пари $(i,j)$ поділеній на загальну кількість пікселів на зображенні $N$, визначає спільну функцію ймовірності ймовірності в двовимірної гістограмі:

$P_{ij}=\frac{f_{ij}}{N},\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}{P}_{ij}=1$

Двомірний метод Оцу буде розроблений на основі двовимірної гістограми наступним чином.

Ймовірності двох класів можна позначити як:

$\begin{array}{cc}{\omega }_0& =\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}{P}_{ij}\\ {\omega }_1& =\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}{P}_{ij}\end{array}$

Інтенсивність означає вектори значень двох класів і загальної середньої вектора, які можуть бути виражені таким чином:

$\begin{array}{cc}{\mu }_0& =[{\mu }_{0i},{\mu }_{0j}{]}^T={[\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}i\frac{P_{ij}}{{\omega }_0},\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}j\frac{P_{ij}}{{\omega }_0}]}^T\\ {\mu }_1& =[{\mu }_{1i},{\mu }_{1j}{]}^T={[\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}i\frac{P_{ij}}{{\omega }_1},\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}j\frac{P_{ij}}{{\omega }_1}]}^T\\ {\mu }_T& =[{\mu }_{Ti},{\mu }_{Tj}{]}^T={[\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}i{P}_{ij},\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}j{P}_{ij}]}^T\end{array}$

У більшості випадків ймовірність недіагональної оцінки буде незначною, тому легко перевірити:

${\omega }_0+{\omega }_1\cong 1$

${\omega }_0{\mu }_0+{\omega }_1{\mu }_1\cong {\mu }_T$

Міжкласова дискретна матриця визначається як

$S_b=\sum _{k=0}^1{\omega }_k[({\mu }_k-{\mu }_T)({\mu }_k-{\mu }_T{)}^T]$

Слід дискретної матриці може бути виражений як

$\begin{array}{cc}& \mathrm{tr}(S_b)\\ =& {\omega }_0[({\mu }_{0i}-{\mu }_{Ti}{)}^2+({\mu }_{0j}-{\mu }_{Tj}{)}^2]+{\omega }_1[({\mu }_{1i}-{\mu }_{Ti}{)}^2+({\mu }_{1j}-{\mu }_{Tj}{)}^2]\\ =& \frac{({\mu }_{Ti}{\omega }_0-{\mu }_i{)}^2+({\mu }_{Tj}{\omega }_0-{\mu }_j{)}^2}{{\omega }_0(1-{\omega }_0)}\end{array}$

де

${\mu }_i=\sum _{i=0}^s\sum _{j=0}^{t}i{P}_{ij}$

${\mu }_j=\sum _{i=0}^s\sum _{j=0}^{t}j{P}_{ij}$

Подібно одновимірному методу Оцу, оптимальний поріг $(s,t)$ виходить шляхом максимізації $\mathrm{tr}(S_b)$.

Алгоритм

$s$ і $t$ отримуються ітеративно, що дуже схоже на одновимірний Метод Оцу. Значення $s$ і $t$ змінюються доки ми не досягнемо максимуму $\mathrm{tr}(S_b)$, що є

max,s,t = 0;
for ss: 0 to L-1 do
    for tt: 0 to L-1 do
        evaluate tr(S_b);
        if tr(S_b) > max
            max = tr(S,b);
            s = ss;
            t = tt;
        end if
    end for
end for
return s,t;

Matlab implementation

функції входу та виходу:

function threshold = 2D_otsu(hists, total)
maximum = 0.0;
threshold = 0;
helperVec = 0:255;
mu_t0 = sum(sum(repmat(helperVec',1,256).*hists));
mu_t1 = sum(sum(repmat(helperVec,256,1).*hists));
p_0 = zeros(256);
mu_i = p_0;
mu_j = p_0;
for ii = 1:256
    for jj = 1:256
        if jj == 1
            if ii == 1
                p_0(1,1) = hists(1,1);
            else
                p_0(ii,1) = p_0(ii-1,1) + hists(ii,1);
                mu_i(ii,1) = mu_i(ii-1,1)+(ii-1)*hists(ii,1);
                mu_j(ii,1) = mu_j(ii-1,1);
            end
        else
            p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj);
            mu_i(ii,jj) = mu_i(ii,jj-1)+mu_i(ii-1,jj)-mu_i(ii-1,jj-1)+(ii-1)*hists(ii,jj);
            mu_j(ii,jj) = mu_j(ii,jj-1)+mu_j(ii-1,jj)-mu_j(ii-1,jj-1)+(jj-1)*hists(ii,jj);
        end

        if (p_0(ii,jj) == 0)
            continue;
        end
        if (p_0(ii,jj) == total)
            break;
        end
        tr = ((mu_i(ii,jj)-p_0(ii,jj)*mu_t0)^2 + (mu_j(ii,jj)-p_0(ii,jj)*mu_t0)^2)/(p_0(ii,jj)*(1-p_0(ii,jj)));

        if ( tr >= maximum )
            threshold = ii;
            maximum = tr;
        end
    end
end
end

Див. також

• Кластеризація зображення за фрагментами інтенсивності
• Implementation of Otsu's thresholding method as GIMP-plugin using Script-Fu (a Scheme-based language).
• Lecture notes on thresholding — covers the Otsu method.
• A plugin for ImageJ using Otsu's method to do the threshold.
• A full explanation of Otsu's method with a working example and Java implementation.
• Implementation of Otsu's method in ITK
• Otsu Thresholding in C# A straightforward C# implementation with explanation.
• Otsu's method using MATLAB