Метод Оцу

Метод Оцу (англ. Otsu's method) — это алгоритм вычисления порога бинаризации для полутонового изображения, используемый в области компьютерного распознавания образов и обработки изображений для получения чёрно-белых изображений. Назван в честь доктора инженерии Токийского университета Нобуюки Отсу^[англ.].

Алгоритм позволяет разделить пиксели двух классов («полезные» и «фоновые»), рассчитывая такой порог, чтобы внутриклассовая дисперсия была минимальной[1]. Метод Оцу также имеет улучшенную версию для поддержки нескольких уровней изображения[2], который получил название мульти-Оцу метод.

В различных русскоязычных источниках можно встретить различные способы написания фамилии автора, например, можно встретить Метод Отса или Метод Отсу.

Математический метод

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

${\displaystyle \sigma _{w}^{2}(t)=\omega _{1}(t)\sigma _{1}^{2}(t)+\omega _{2}(t)\sigma _{2}^{2}(t)}$,

где ${\displaystyle i}$ — число, указывающее на класс пикселя, веса ${\displaystyle \omega _{i}}$ — это вероятности двух классов, разделённых порогом ${\displaystyle t}$, ${\displaystyle \sigma _{i}^{2}}$ — дисперсия этих классов.

Оцу показал, что минимизация дисперсии внутри класса равносильна максимизации дисперсии между классами:[1]

${\displaystyle \sigma _{b}^{2}(t)=\sigma ^{2}-\sigma _{w}^{2}(t)=\omega _{1}(t)\omega _{2}(t)\left[\mu _{1}(t)-\mu _{2}(t)\right]^{2}}$

которая выражается в терминах вероятности ${\displaystyle \omega _{i}}$ и среднего арифметического класса ${\displaystyle \mu _{i}}$, которое, в свою очередь, может обновляться итеративно. Эта идея привела к эффективному алгоритму.

Алгоритм

Пусть дано полутоновое изображение ${\displaystyle G(i,j),i={\overline {1,Height}},j={\overline {1,Width}}.}$ Счётчик повторений ${\displaystyle k=0.}$

1. Вычислить гистограмму ${\displaystyle p(l)}$ изображения и частоту ${\displaystyle N(l)}$ для каждого уровня интенсивности изображения ${\displaystyle G}$.
2. Вычислить начальные значения для ${\displaystyle \omega _{1}(0),\omega _{2}(0)}$ и ${\displaystyle \mu _{1}(0),\mu _{2}(0).}$
3. Для каждого значения ${\displaystyle t={\overline {1,max(G)}}}$ — полутона — горизонтальная ось гистограммы:
   • Обновляем ${\displaystyle \omega _{1},\omega _{2}}$ и ${\displaystyle \mu _{1},\mu _{2}}$;
   • Вычисляем минимизацию дисперсии внутри класса: ${\displaystyle \sigma _{b}^{2}(t)=\omega _{1}(t)\omega _{2}(t)\left[\mu _{1}(t)-\mu _{2}(t)\right]^{2}}$;
   • Если ${\displaystyle \sigma _{b}^{2}(t)}$ больше, чем имеющееся, то запоминаем ${\displaystyle \sigma _{b}^{2}}$ и значение порога ${\displaystyle t}$, ведь искомый порог соответствует максимуму , иначе повторяем действие 3.
4. Искомый порог соответствует максимуму ${\displaystyle \sigma _{b}^{2}(t).}$
5. Вычисляем итеративно обновляемое среднее арифметическое класса ${\displaystyle \mu _{i}}$:
   • ${\displaystyle N_{T}=\sum _{i=0}^{max(G)}{p(i)}}$,
   • ${\displaystyle \omega _{1}(t)={{\sum _{i=0}^{t-1}{p(i)}} \over {N_{T}}}={\sum _{i=0}^{t-1}{N(i)}},\quad \omega _{2}(t)=1-\omega _{1}(t)}$,
   • ${\displaystyle \mu _{T}={{\sum _{i=0}^{max(G)}{i\cdot p(i)}} \over {N_{T}}}={\sum _{i=0}^{max(G)}{i\cdot N(i)}},}$
   • ${\displaystyle \mu _{1}(t)={{\sum _{i=0}^{t-1}{i\cdot p(i)}} \over {N_{T}\cdot \omega _{1}(t)}}={{\sum _{i=0}^{t-1}{i\cdot N(i)}} \over {\omega _{1}(t)}},\quad \mu _{2}(t)={{\mu _{T}-\mu _{1}(t)\cdot \omega _{1}(t)} \over {\omega _{2}(t)}}.}$

Исходное изображение

Бинаризация с порогом по Оцу

Программная реализация

JavaScript

В данной функции аргумент pixelsNumber является общим числом пикселей изображения, а аргумент histogram — гистограмма 8-битового полутонового изображения, представленная в виде одномерного массива, в котором номер элемента кодирует номер полутона, а значение поля — количество пикселей с этим полутоном.

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 * .7152 + 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");

При выполнении этого кода в браузере рядом с оригиналом появится чёрно-белое изображение, которое было получено при следовании этого алгоритма. Примерный результат можно найти здесь.

Реализация на языке C

// Количество уровней интенсивности у изображения.
// Стандартно для серого изображения -- это 256. От 0 до 255.
const int INTENSITY_LAYER_NUMBER = 256;

// Возвращает гистограмму по интенсивности изображения от 0 до 255 включительно
void calculateHist(const IMAGE *image, const int size, int* hist) {
    // Иницилизируем гистограмму нулями
    memset(hist, 0, INTENSITY_LAYER_NUMBER * sizeof(*hist));

    // Вычисляем гистограмму
    for (int i = 0; i < size; ++i) {
        ++hist[image[i]];
    }
}

// Посчитать сумму всех интенсивностей
int calculateIntensitySum(const IMAGE *image, const int size) {
    int sum = 0;
    for (int i = 0; i < size; ++i) {
        sum += image[i];
    }

    return sum;
}

// Функция возвращает порог бинаризации для полутонового изображения image с общим числом пикселей size.
// const IMAGE *image содержит интенсивность изображения от 0 до 255 включительно.
// size -- количество пикселей в image.
int otsuThreshold(const IMAGE *image, const int size) {
    int hist[INTENSITY_LAYER_NUMBER];
    calculateHist(image, size, hist);

    // Необходимы для быстрого пересчёта разности дисперсий между классами
    int all_pixel_count = size;
    int all_intensity_sum = calculateIntensitySum(image, size);

    int best_thresh = 0;
    double best_sigma = 0.0;

    int first_class_pixel_count = 0;
    int first_class_intensity_sum = 0;

    // Перебираем границу между классами
    // thresh < INTENSITY_LAYER_NUMBER - 1, т.к. при 255 в ноль уходит знаменатель внутри for
    for (int thresh = 0; thresh < INTENSITY_LAYER_NUMBER - 1; ++thresh) {
        first_class_pixel_count += hist[thresh];
        first_class_intensity_sum += thresh * hist[thresh];

        double first_class_prob = first_class_pixel_count / (double) all_pixel_count;
        double second_class_prob = 1.0 - first_class_prob;

        double first_class_mean = first_class_intensity_sum / (double) first_class_pixel_count;
        double second_class_mean = (all_intensity_sum - first_class_intensity_sum) 
            / (double) (all_pixel_count - first_class_pixel_count);

        double mean_delta = first_class_mean - second_class_mean;

        double sigma = first_class_prob * second_class_prob * mean_delta * mean_delta;

        if (sigma > best_sigma) {
            best_sigma = sigma;
            best_thresh = thresh;
        }
    }

   return best_thresh;
}

Примечания

1. N. Otsu. A threshold selection method from gray-level histograms (англ.) // IEEE Trans. Sys., Man., Cyber. : journal. — 1979. — Vol. 9. — P. 62—66.
2. Ping-Sung Liao and Tse-Sheng Chen and Pau-Choo Chung. A Fast Algorithm for Multilevel Thresholding (неопр.) // J. Inf. Sci. Eng.. — 2001. — Т. 17. — С. 713—727.

Ссылки

• Lecture notes on thresholding — covers the Otsu method.