.nr PS 12
.nr VS 12
.\"fspecial TR CTR
.\"fspecial TI CTI
.\"fspecial TB CTB
.\"kern 1 
.special SS S
.fspecial SS
.ds ABSTRACT РЕФЕРАТ
.ce 1
.TL
.PP
.AI
Министерство образования и науки Российской Федерации
.PP
.AI
Федеральное государственное бюджетное образовательное учреждение
.PP
.AI
высшего профессионального образования
.PP
.AI
СЫКТЫВКАРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
.PP
.AI
Институт точных наук и информационных технологий
.PP
.AI
Кафедра прикладной математики
.PP
.AI
Моделирование фона видеоизображений с помощью гауссовых распределений.

.AI
Сыктывкарский Государственный Университет
.EQ
delim ##
define vdots '\[u22EE]'
define ddots '\[u22F1]'
define rdtos '\[u22F0]'
.EN
.bp \" page 1
Содержание
.PP
1) Введение
.PP
2) Реализация
.RS
.RS
.IP 2.1)
Построение функции правдоподобия. EM алгоритм.
.IP 2.2)
Формула log-sum-exp
.IP 2.3)
Эффективное вычисление функции плотности
.IP 2.4)
Практическая реализация EM алгоритма.
.RS
.IP 2.4.1)
Условия задачи
.IP 2.4.2)
Формулы максимизации
.IP 2.4.3)
Моделирование фона
.IP 2.4.4)
Выбор начальных средних с помощью алгоритма k-means (k средних).
.bp
.IP 2.4.5)
Оптимальный выбор средних для случая k=2. 
.RE
.RE
.RE
3) Список литературы.
.PP
.bp \" page 1
.NH
Введение.
.PP
В задачах распознавания образов есть класс задач, в которых требуется определить
движущиеся объекты на видео последовательности для дальнейшего определения к какому классу принадлежит
объект, подсчёта объектов или расчёта скорости и траектории движения объекта. Можно вести поиск интересующего нас
объекта в каждом кадре видео по всем его пикселям, но для облегчения задачи можно
определять сразу относится ли пиксель к фону, и тогда его можно не рассматривать в расчётах.
Для решения задачи разделения каждого пикселя видео на фон и не фон, можно применить несколько способов:
.PP
.IP 1)
Вычитать из текущего кадра предыдущий и оценивать по некоторому порогу относится пиксель к фону, либо нет.
Данный способ будет давать большую погрешность, когда фон не является статичным 
(например: падающий снег, колышущиеся деревья и трава, пузырьки в воде рис. 1) либо
камера сильно "шумит".
.PP
.IP 2)
Строить некоторую статистическую модель для фона, в которой мы можем учитывать изменения фона по времени
и шум камеры.
.PP
В данной работе рассматривается моделирование фона видео изображения смесью гауссиан при помощи EM и k-means
алгоритмов, построение модели фона,
ускорение стандартных алгоритмов вычисления фона.
.PP
Данная работа актуальна для автоматизации учёта, робототехники и систем безопасности.
.PP
.PP
.PSPIC -C fish.eps 5i 5i
.RS
.ce 1
Рис. 1 Пример нестатичного фона
.RE
.bp \" page 1
.NH 1
Реализация
.LP
.NH 2
Построение функции правдоподобия. EM алгоритм.
.PP
Для набора данных (пикселей) будем искать функцию правдоподобия в виде суммы гауссиан с помощью EM - алгоритма
(expectation - maximization - algorithm). В общем случае идея алгоритма заключается в следующем.
.PP
Пусть плотность распределения на множестве #bold X# имеет вид смеси #bold k# распределений:
.PP
.EQ
p( bold x) = sum from {j=1} to k
  omega sub j ^ p sub j ( bold x), sum from {j=1} to k omega sub j ^ = 1, omega sub j ^ >= 0.
.EN
.PP
где # p sub j# - функция правдоподобия # j #-й компоненты смеси, # omega sub j# - ее априорная вероятность.
Пусть функции правдоподобия принадлежат параметрическому семейству распределений
.EQ
phi ( bold x, ETA )
.EN
и отличаются только значениями параметра
.EQ
p sub j = phi ( bold x ;  ETA sub j )
.EN
Задача разделения смеси заключается в том, чтобы, имея выборку #bold X sup m# случайных и независимых наблюдений из смеси #p( bold x )#, зная число #k# и функцию #phi#, оценить вектор параметров #THETA = ( omega sub 1,..., omega sup k; ETA sub 1,..., ETA sub k )#
.PP
Искусственно вводится вспомогательный вектор скрытых переменных #G#, обладающий двумя замечательными свойствами.
.PP
.IP а)
С одной стороны, он может быть вычислен, если известен вектор вектор параметров #THETA#
.IP б)
С другой стороны, поиск максимума правдоподобия сильно упрощается, если известны значения скрытых переменных.
.PP
EM - алгоритм состоит из итерационного повторения двух шагов:
.IP
На # bold E-шаге # вычисляется ожидаемое значение (expectation) вектора скрытых переменных #G# по текущему приближению вектора параметров #THETA#.
Обозначим через #p ( bold x , ETA sub j )# плотность вероятности того, что объект # bold x # получен из #j#-й компоненты смеси. По формуле условной вероятности
.EQ
p ( bold x, ETA sub j ) = p( bold x ) ^ P( ETA sub j | bold x ) = omega sub j ^ p sub j ( bold x )
.EN
Введём обозначение
.EQ
g sub ij = P( ETA sub j | bold x sub i )
.EN
Это неизвестная апостериорная вероятность того, что обучающий объект # bold x sub i # получен из #j#-й компоненты смеси. Возьмём эти величины в качестве скрытых переменных.
.EQ
sum from {j = 1} to k g sub ij = 1,
.EN
для любого #i=1..m#, так как имеет смысл полной вероятности принадлежать объекту #bold x sub i # одной из #bold k# компонент смеси. Из формулы Байеса:
.EQ
g sub ij = { omega sub j ^ p sub j ( x sub i ) } over { sum from {s = 1} to k omega sub s ^ p sub s (x sub i ) }
.EN
.IP
На # bold M-шаге# решается задача максимизации правдоподобия (maximization) и находится следующее приближение вектора #THETA# по текущим значениям векторов #G# и #THETA#.
Будем максимизировать логарифм полного правдоподобия:
.EQ
Q ( THETA ) = ln prod from {i = 1} to m p ( x sub i ) = sum from {i = 1} to m sum from {j = 1} to k omega sub j ^ p sub j ( x sub i )
.EN
Решая оптимизационную задачу Лагранжа с ограничением на сумму #omega sub j#, находим:
.EQ
omega sub j = 1 over m sum from {i = 1} to m g sub ij, j = 1..k 
.EN
.EQ
Eta sub j = argmax sub { ETA } sum from {i = 1} from m g sub ij ln phi ( x, ETA ) , j = 1..k
.EN
.PP
Рассмотрим данный алгоритм более подробно для смеси гауссовых распределений (mixture of gaussian):
.PP
Формула функции плотности многомерного гауссовского
распределения имеет вид:
.EQ (1)
eta ( bold x ~|~ mu ,~SIGMA ) = 
 1 over {(2 pi ) sup {M/2} left | SIGMA right | sup {1/2} }
 italic e sup { -{{1 over 2} ( bold x ~-~ bold mu ) sup T ~
  SIGMA ~ sup -1 ( bold x ~-~ bold mu )}}.
.EN
.ne 3
Здесь:
.RS
.IP # bold x #
вектор #(x sub 1 , x sub 2 , x sub 3 , ..., x sub M )#,
размерности \fIM\fP;
.IP #SIGMA#
ковариационная матрица, а #left | SIGMA right |# \(em её
детерминант;
.IP #bold mu #
вектор средних значений # bold x #.
.RE
.LP
Обозначение ``# eta ( bold x ~|~ bold mu ,~SIGMA )#''
указывает на то, что функция получает
вектор #bold x#, как аргумент, и параметризуется значениями
#SIGMA# и #bold mu#.
.PP
Плотность вероятности для смеси из #K# гауссовских
распределений, вес каждой из которых задан #omega sub k#,
задаётся так:
.EQ (2)
f sub bold x ( bold x ) = sum from {k=1} to K
  omega sub k ^ eta ( bold x ~|~ bold mu sub k,~SIGMA sub k ).
.EN
Так как функция #f sub bold x ( bold x )# есть плотность
вероятности, её интеграл должен быть равен единице,
очевидно, что сумма всех весов #sum omega# тоже должна
равняться единице, и все\ #omega sub k > 0#.
.PP
Функция правдоподобия рассчитывается всегда для конечного
набора значений #bold x sub n#, как произведение
их плотностей вероятности:
.EQ
L =  f sub bold x ( bold x sub 1 ) 
 f sub bold x ( bold x sub 2 ) ...
 f sub bold x ( bold x sub N ) =
 prod from {n=1} to N f sub bold x ( bold x sub n ).
.EN
.PP
Функция правдоподобия для смеси из #K# гауссовских
распределений для всего набора значений
#( bold x sub 1 , bold x sub 2 , ..., bold x sub N )# имеет вид:
.EQ (3)
L = prod from {n=1} to N sum from {k=1} to K
  omega sub k eta ( bold x sub n ~|~ bold mu sub k,~SIGMA sub k ).
.EN
В практическом смысле, в ходе выполнения какого-то алгоритма,
нас чаще всего интересует как изменялась функция
правдоподобности, а не её абсолютное значение. Поэтому
вместо #L# можно взять её логарифм, при этом
#log L# превратит произведение в сумму:
.DS L
.EQ
  log L =
   log
    left [
     prod from {n=1} to N
      sum from {k=1} to K
       omega sub k eta ( bold x sub n ~|~ bold mu sub k,~SIGMA sub k )
    right ]
.EN
.EQ
  lineup = 
  log sum from {k=1} to K
  omega sub k 
  eta ( bold x sub 1 ~|~ bold mu sub k,~SIGMA sub k ) ~+~
  ... ~+~
 log sum from {k=1} to K
  omega sub k
  eta ( bold x sub N ~|~ bold mu sub k,~SIGMA sub k )
.EN
.EQ
 lineup =
 sum from {n=1} to N log sum from {k=1} to K
  omega sub k 
  eta ( bold x sub n ~|~ bold mu sub k,~SIGMA sub k ).
.EN
.DE
.PP
Существует один важный практический аспект. При подсчёте
сумм плотностей,
часто значения плотностей оказываются настолько малыми,
что в операциях с плавающей точкой происходит исчезновение
порядка (floating point underflow).
В арифметике IEEE-754 исчезновение порядка приводит к тому
что число денормализуется и плавно, а не резко, приближается
к нулю.
Большинство сопроцессоров операции с денормализованными
числами выполняют через микропрограмму, а не аппаратно,
что обычно медленнее на порядок.
Потому целесообразно и при вычислении\ # eta 
( bold x ~|~ bold mu ,~SIGMA )#
переходить к логарифмированию:
.EQ
log eta ( bold x ~|~ bold mu ,~SIGMA ) mark~= 
 ~-~ M over 2 log {(2 pi )}
 ~-~ {1 over 2} log left | SIGMA right |
 ~-~ {1 over 2} ( bold x ~-~ bold mu ) sup T ~
  SIGMA ~ sup -1 ( bold x ~-~ bold mu )
.EN
.EQ
lineup~= -{1 over 2} left ( M log (2 pi )
 ~+~ log left | SIGMA right |
 ~+~ ( bold x ~-~ bold mu ) sup T ~ 
    SIGMA ~ sup -1 ( bold x ~-~ bold mu )
 right )^.
.EN
После логарифмирования, значения должны быть восстановлены
экспоненцированием.
.PP
.LP
.NH 2
# bold Формула# # bold log-sum-exp #
.PP
В формуле для смеси гауссовских распределений каждая
компонента плотности имеет вес #omega sub k#. Будем
обозначать логарифм плотности распределения с весом
#omega# как\ #gamma #:
.EQ
gamma = log eta ( bold x ~|~ bold mu ,~SIGMA ) + log omega
.EN
Формула для расчёта #log L# тогда примет вид:
.EQ
 log L = sum from {n=1} to N log
 left (^
  sum from {k=1} to K
   exp^( gamma sub k )
 ^ right )^.
.EN
.PP
Значение #exp( gamma sub k )# может опять таки оказаться
настолько малым, что это приведёт к исчезновению порядка.
Для того чтобы избежать этого,
при вычислении логарифма суммы сначала находится самое
большое #gamma sub max#, а затем применяется формула
.EQ
log left (^
  sum from {k=1} to K
   exp^( gamma sub k )
 ^right )^ =
gamma sub max + 
 log left (^
  sum from {k=1} to K
   exp^( gamma sub k - gamma sub max )
 ^right )^,
.EN
которая гарантирует,
что хотя бы одно значение не вызовет исчезновения порядка,
а остальные значения в таком случае
всё равно не окажут влияния на сумму в силу своей малости.
Такой трюк известен под названием
\(Fo\fIформулa\ log-sum-exp\fP\(Fc.
.PP
.LP
.NH 2
# bold Эффективное # # bold вычисление # # bold функции # # bold плотности #
.PP
Функция плотности гауссовского распределения (1)
содержит скаляр
.EQ
( bold x ~-~ bold mu ) sup T ~ 
  SIGMA ~ sup -1 ( bold x ~-~ bold mu ) =
 theta sup T ~ SIGMA ~ sup -1 theta .
.EN
Конечно, можно вычислить это значение,
как говорится, \(Foв лоб\(Fc,
найдя обратную ковариационную матрицу и выполнив две операции
умножения \(Foматрица-вектор\(Fc.
Тем не менее, существует способ сделать это проще и
компактнее, что не маловажно для систем,
обладающих скромными ресурсами. Свойства ковариационной
матрицы позволяют осуществить упрощение, потому что эта
матрица \fIсимметричная\fP и,
в большинстве практических случаев,
\fIположительно-определённая\fP,
а значит её можно факторизовать,
пользуясь \fIразложением Холецкого\fP,
на две треугольные матрицы
.EQ
bold SIGMA = bold {CC} sup T ,
.EN
где:
.RS
.IP #bold C#
нижняя треугольная матрица со строго положительными
диагональными элементами;
.IP #bold C sup T#
верхняя треугольная матрица.
.LP
.RE
Любую симметричную положительно-определённую матрицу можно
представить в таком виде.
Взяв нижнюю треугольную матрицу #bold C#,
мы можем очень быстро,
через прямую подстановку,
найти вектор #u# такой,
что будет выполняться равенство
.EQ
bold C u = theta .
.EN
Если теперь заменить вектор #theta# на #bold C u#,
получим:
.EQ
 theta sup T ~ SIGMA ~ sup -1 theta ~mark = 
  ( bold C u) sup T ^( bold C bold C sup T ) sup -1 ^bold C u
.EN
.EQ
 ~lineup =
   u sup T bold C sup T
   bold C sup -T bold C sup -1
   bold C u
.EN
.EQ
~lineup = u sup T u = u cdot u.
.EN
Алгоритм разложения Холецкого достаточно прост и эффективен
в реализации на ЭВМ,
и,
в том случае если размерность вектора
#bold x# невелика и известна заранее,
все циклы алгоритма,
вычисляющие коэффициенты матриц разложения,
можно развёрнуть.
.PP

.PP
.LP
.NH 2
Практическая реализация EM алгоритма.
.PP
.LP
.NH 3
# bold Условия # # bold задачи#
.PP
Дано #N# точек
#( bold x sub 1 , bold x sub 2 , ... , bold x sub N )#
размерности #M#.
Известно,
что все точки принадлежат #K#
многомерным гауссовским распределениям вида
#omega sub k^ eta ( bold x ~|~ bold mu sub k,~SIGMA sub k )#
с неизвестными параметрами
#bold mu sub k#,
#SIGMA sub k#,
#omega sub k#:
.RS
.IP
#bold mu sub k#
вектор средних значений для #k#-го распределения;
.IP
#SIGMA sub k#
ковариационная матрица #M times M# для #k#-го распределения;
.IP
#omega sub k#
вес #k#-го распределения, #sum from k omega sub k = 1#.
.RE
.PP
Каждая точка #bold x sub n# имеет некоторую вероятность
присутствия в распределении #k#.
Эту вероятность мы будем обозначать как #P sub {nk}#.
Все параметры #P sub {nk}# для всех точек удобно свести
в матрицу #P sub {N times K}#:
.EQ
P sub {N times K} ~=~
 left [
  pile {
   P sub {11} above
   P sub {21} above
   vdots above 
   P sub {N1}
  } ~
  pile {
   P sub {12} above
   P sub {22} above
   vdots above
   P sub {N2}
  } ~
  pile {
   P sub {13} above
   P sub {23} above
   vdots above
   P sub {N3}
  } ~ 
  pile {
   ... above
   ... above
   ddots above
   ...
  } ~
  pile {
   P sub {1K} above
   P sub {2K} above
   vdots above
   P sub {NK}
  }
 right ] .
.EN
Правдоподобность и вероятности #P sub {nk} #
.PP
Вероятностную плотность для точки с координатами 
#bold x# будем обозначать как\ #P( bold x )#.
Это вероятность попасть в точку по этим координатам,
если случайным образом выбирать её из всей выборки.
Правдоподобность параметров распределения будем обозначать
как L.
Наша задача максимизировать L,
подобрав параметры #bold mu sub k#, #SIGMA sub k#,
#omega sub k#.
Так как мы полагаем точки с данными независимыми,
их L будет произведением вероятностей присутствия точки из
выборки в #bold x sub n#:
.EQ
 L~=~ prod from n=1 to N P( bold x sub n ) =
  P( bold x sub 1 ) ^times^
  P( bold x sub 2 ) ^times^ ... ^times^
  P( bold x sub N ).
.EN
Для того,
чтобы посчитать L,
нам нужно найти плотности распределения в точках
#bold x sub n# и перемножить их.
Для каждой точки #bold x sub n# плотность определяется
как сумма плотностей из #K# распределений,
пропорционально весу #k#-го распределения:
.EQ
P( bold x sub n ) = sum from {k=1} to K
  omega sub k ^ p sub j ^ ( bold x).
.EN
Так как значения #P sub {nk}#,
то есть значения строки #n# матрицы #P sub {N times K}#,
есть вероятности вхождения #bold x sub n# в
распределение #k#,
их сумма должна равняться единице:
.EQ
sum from k=1 to k P sub {nk} = 1.
.EN
Значит #P( bold x sub n )# это сумма
.EQ
P( bold x sub n ) =
 P sub {n1} cdot P( bold x sub n ) ~+~
 P sub {n2} cdot P( bold x sub n ) ~+~ ... ~+~
 P sub {NK} cdot P( bold x sub n ),
.EN
но нам также известно,
что
.EQ
P( bold x sub n ) =
 omega sub 1 ^ eta ( bold x sub n ~|~ bold mu sub 1 ,~SIGMA sub 1 )
 ~+~
 omega sub 2 ^ eta ( bold x sub n ~|~ bold mu sub 2 ,~SIGMA sub 2 )
 ~+~ ... ~+~
 omega sub K ^ eta ( bold x sub n ~|~ bold mu sub K ,~SIGMA sub K ).
.EN
Последнее означает,
что каждое значение #P sub {nk}# можно вычислить по формуле:
.EQ
P sub {nk} =
 {omega sub k ^ eta ( bold x sub n ~|~ bold mu sub k ,~SIGMA sub k )} over 
 {P( bold x sub n )}.
.EN
При заполнении матрицы #P sub {N times K}#
значениями #P sub {nk}#,
вычисления рационально построить по следующей схеме:
.RS
Для каждой строки #n# матрицы #P sub {N times K}#:
.IP а)
Заполнить строку значениями
#P sub {nk} cdot P( bold x sub n ) =
omega sub k ^ eta ( bold x sub n ~|~ bold mu sub k ,~SIGMA sub k ).#
.IP б)
Подсчитать сумму всех элементов в строке #sum P sub {nk}#.
.IP в)
Поделить все элементы строки на сумму,
полученную на предыдущем шаге.
.RE
.LP
.bp \"page 6
.NH 3
#bold Формулы# #bold максимизации#
.PP
Значения #bold mu bar sub k ,~SIGMA bar sub k ,~omega bar sub k#
на шаге максимизации вычисляются по следующим формулам:
.DS
.EQ (1)
omega bar sub k ~=~ 1 over N sum from n P sub {nk}
.EN
.EQ (2)
bold mu bar sub k ~=~
 {sum from n P sub {nk} cdot bold x sub n } over {sum from n to nothing P sub {nk} }
.EN
.EQ (3)
SIGMA bar sub k ~=~
 {
  sum from n
   P sub {nk} cdot ( bold x sub n - italic {bold mu} bar sub k )
    ( bold x sub n - bold mu bar sub k ) sup T 
 } over
 {
  sum from n to nothing
   P sub {nk}
 }
.EN
.DE
Пошаговое выполнение алгоритма
.PP
.nr step 1 1
Пользователь предоставляет данные #bold x sub n#,
размерности #N#, #M# и #K#,
значение #epsilon# и,
возможно,
вектор начальных значений для #bold mu sub k#.
Реализация алгоритма может состоять из следующих шагов:
.IP \n[step]. 3
На первом шаге:
.RS
.IP а)
Установить веса #omega sub k#, как #omega sub k = 1 over k#.
.IP б)
Установить #SIGMA sub k# в
#I = size -2 
{
  left [ matrix { 
  ccol { 1 above 0 above size -2 {vdots} above 0 }
  ccol { 0 above 1 above size -2 {vdots} above 0 }
  ccol { size -1 {...} above size -1 {...} above size -2 {ddots} above size -1 {...} }
  ccol { 0 above 0 above size -2 {vdots} above 1 }
 } right ]
}.#
.IP в)
Выбрать #bold mu sub k#,
как #K# случайно взятых значения из набора #bold x sub n#,
или же взять начальные значения,
предоставленные пользователем.
.RE
.IP \n+[step].
На втором шаге:
.RS
.IP а)
Выполнить шаг \fIE\fP:
.RS
.IP а.1)
вычислить #P sub {nk}# и #P( bold x sub n )#;
.IP а.2)
вычислить сумму всех логарифмов #P( bold x sub n )#;
.IP а.3)
сохранить #size -2 {sum from n} log P( bold x sub n )#
в переменную #loglike_prev#.
.RE
.RE
.RS
.IP б)
Выполнить шаг \fIM\fP:
.RS
.IP б.1)
пользуясь формулой (1) и весами #P sub {N times K}#,
вычислить #omega bar sub k#;
.IP б.2)
пользуясь формулой (2) и весами #P sub {N times K}#,
вычислить #bold mu bar sub k#;
.IP б.3)
пользуясь формулой (3) и весами #P sub {N times K}#,
вычислить #SIGMA bar sub k#.
.RE
.RE
.IP \n+[step].
На третьем шаге:
.RS
.IP а)
Выполнить шаг \fIE\fP,
как в 2а,
но сохраняя результат 2а.3 в переменную \fIloglike\fP.
.IP б)
Вычислить #abs(loglike - loglike_prev)#,
и сравнить с #epsilon#,
заданным пользователем.
Если #abs(loglike - loglike_prev) < epsilon#,
тогда перейти на шаг 4.
.IP в)
Выполнить шаг \fIM\fP, как в 2б.
.IP г)
Перейти на шаг 3.
.RE
.IP \n+[step].
Завершить выполнение алгоритма с успешным статусом,
а значения\ #bold mu bar sub k,~SIGMA bar sub k,~omega bar sub k#
вернуть пользователю.
.PSPIC -C Hist.eps 3i 3i
.RS
.ce 1
Рис. 2 Сумма гауссиан построенных по выборке
.RE
.PP
 Примеры с кластеризацией данных с помощью EM алгоритма.
.PP
.PSPIC -C clusters1.eps 3i 3i
.RS
.ce 1
Рис. 3 Кластеризация
.RE
.PP
.PSPIC -C clusters2.eps 3i 3i
.RS
.ce 1
Рис. 4 Кластеризация
.RE
.PP
.bp \"page 8
.LP
.NH 3
#bold Моделирование# # bold фона .#
.PP
Рассмотрим теперь моделирование фона видеоизображения: пусть #V sub {n times s} - # область видеоизображения размером #n# на #s#. Распишем пошагово наши действия:
.IP a)
По нескольким первым кадрам с некоторым прореживанием собрать иcторию пикселей (набрать выборку) по области #V#.
.PSPIC -C Bufer.eps
.RS
.ce 1
Рис. 5 Схема сбора данных
.RE
.PP
.IP b)
Вычислить модель фона с помощью #EM-алгоритма# для выборки.
.IP c)
Для каждого пикселя области #V# последующих кадров:
.RS
.IP c.1)
Вычислить по некоторому порогу, принадлежит ли данный пиксель модели фона, если нет, то  отнести этот пиксель к движению.
.IP c.2)
Обновить модель фона по данному пикселю с помощью IIR фильтра.
.IP c.3)
Добавить данный пиксель в историю (выборку).
.IP c.4)
При необходимости пересчитать модель фона по области #V#.
.RE

.PP
.LP
.\".bp \"page 12
.NH 3
Выбор начальных средних с помощью алгоритма k-means (k средних).
.PP
Проблемой EM алгоритма является, то что он может на некоторых наборах данных
долго сходиться и давать плохие результаты. Поэтому важным является выбор начальных
значений для вектора средних.
.PP
Рассмотрим применение алгоритма k-means к начальному выбору средних.
Действие алгоритма таково, что он стремится минимизировать суммарное квадратичное отклонение точек кластеров от центров этих кластеров:
.EQ
V={sum from i=1 to k}{sum from {x sub j \[mo] S sub j}} (x sub j - mu sub j ) sup 2
.EN
где #k - # число кластеров, #S sub j - # полученные кластеры, и #mu sub j -# центры масс векторов 
.PP
Алгоритм нахождения центров масс представляет собой версию EM-алгоритма, применяемого также для разделения смеси гауссиан. Он разбивает множество элементов векторного пространства на заранее известное число кластеров k.
Основная идея заключается в том, что на каждой итерации перевычисляется центр масс для каждого кластера, полученного на предыдущем шаге, затем векторы разбиваются на кластеры вновь в соответствии с тем, какой из новых центров оказался ближе по выбранной метрике.
Алгоритм завершается, когда на какой-то итерации не происходит изменения кластеров. Это происходит за конечное число итераций, так как количество возможных разбиений конечного множество конечно, а на каждом шаге суммарное квадратичное уклонение V уменьшается, поэтому зацикливание невозможно.
.LP
.\".bp \"page 12
.NH 3
Оптимальный выбор средних для случая k = 2.
.PP
Для более быстрой и точной сходимости EM алгоритма для случая ограниченного пространства RGB можно выбрать начальное приближение для среднего более точно. Для этого рассмотрим одномерное пространство Y, где 
.EQ
Y = round(0.299 cdot R + 0.587 cdot G + 0.114 cdot B), Y \[mo] [0:255]
.EN
Будем считать, что Y принимает целые значения.
.PP
Если мы рассмотрим гистограмму нашей выборки пикселей в Y, нам интуитивно будет понятно где брать значения средних. Но как объяснить это компьютеру? Для начала найдём точку на гистограмме, где будет разделение гауссовых распределений. Гистограмму, которая приведена ниже, можно было бы разделить где-то на промежутке от 6 до 15, а потом выбрать за средние максимумы в точках 5 и 17, а потом обратно перейти в формат RGB, запомнив заранее какой нибудь пиксель попавший в выбранные начения Y. Для выбора точки разделения рассмотрим функцию
.EQ
T(y) = F(x) cdot y + (1 - F(y + 1)) cdot (255 - y), y \[mo] [0:255]
.EN
где F(y) - функция распределения дискретной величины y.
.PP
Будем брать точкой разделения точку минимума функции #T(y)#. Далее будем брать за первое среднее точку Y с наибольшим значением гистограммы слева от этой точки, а за второе - с наибольшим значением гистограммы справа от этой точки. Данный алгоритм можно обобщить на случай k больших 2-х, либо использовать данный алгоритм для пространства RGB. 
.PP
Данный выбор средних даёт ускорение сходимости в несколько раз большее, чем при случайном выборе и выборе с помощью алгоритма k-means.
.PSPIC -C hist.eps 5i 5i
.RS
.ce 1
Рис. 6 Гистограмма данных
.RE

.PSPIC -C value.eps 5i 5i
.RS
.ce 1
Рис. 7 Функция #T(y)#
.RE

.PP 
Заметим, что данный алгоритм можно улучшать и модифицировать для достижения большей точности и высокой скорости работы.
.PSPIC -C fish_recognition.eps 5i 5i
.RS
.ce 1
Рис. 8 Разделение на фон и "не фон".
.RE
.PSPIC -C auto_recognition.eps 5i 5i
.RS
.ce 1
Рис. 9 Разделение на фон и "не фон".
.RE
.bp \"page 14
.NH 1
Список литературы:
.PP
.IP 1)
Background Modeling using Mixture of Gaussians for Foreground
Detection - A Survey, T. Bouwmans, F. El Baf, B. Vachon
.IP 2)
The Expectation Maximization Algorithm. A short tutorial. Sean Borman
.IP 3)
http://ru.wikipedia.org/wiki/EM-алгоритм
.IP 4)
http://www.machinelearning.ru/wiki/index.php?title=EM_алгоритм_(пример)
.IP 5)
http://ru.wikipedia.org/wiki/K-means


