CC BY
247
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
Трехмерные и двумерные изображения: модели, алгоритмы и области анализа
П.А. Чочиа
Аннотация — Рассматриваются вопросы модификации двухмасштабной модели и алгоритмов обработки при переходе от двумерных к трехмерным изображениям. Показаны изменения области анализа, алгоритмов фильтрации, декомпозиции изображений и обнаружения объектов. Предложены быстрые алгоритмы вычисления локального среднего и порядковых статистик по скользящему окну для 3Б-изображений.
Ключевые слова — Обработка изображений, трехмерное изображение, модель изображений, алгоритм обработки, область анализа, быстрые алгоритмы.
I. Введение
Когда употребляют термин «трехмерное изображение» зачастую понимают совершенно различные виды данных [1]-[3], основные из которых следующие.
1. Данные, задаваемые функцией трех координат и являющиеся гомеоморфным отображением некоторого объемного участка трехмерного пространства, включая все содержащиеся в нем объекты.
2. Стереоскопическое изображение, состоящее из пары двумерных изображений, за счет диспаратности дающих наблюдателю представление о расположении объектов.
3. Изображение, являющееся проекцией трехмерной сцены (например, аксонометрической), позволяющее оценить форму и расположение объектов, но при этом остающееся двумерным.
4. Двумерное изображение, каждая точка которого соответствует некоторым координатам в трехмерном пространстве, например дальности или рельефу.
5. Особым способом сформированные изображения, создающие образы объектов, например голограммы.
6. Видеопоследовательности, содержащие набор кадров объектов. Такие данные могут представляться в виде трехмерного массива, но одна из координат при этом является не пространственной, а координатой времени.
Далее под трехмерным или ЗБ-изображением (непрерывным или дискретным) будут пониматься изображения исключительно первого типа. По существу, 3D-изображение представляет собой расширение обычного 2D-изображения путем добавления еще одного пространственного измерения.
Способы формирования трехмерного изображения могут быть различными; наиболее известным является
Статья получена 20 ноября 2014.
П.А. Чочиа, к.т.н., с.н.с., Институт проблем передачи информации РАН, г. Москва. (e-mail: chochia@iitp.ru)).
томографическое сканирование — рентгеновская компьютерная томография или магнитно-резонансная томография. Возможно получение 3D-изображения в результате сейсморазведки при геологических исследованиях, в микроскопии при использовании объектива с переменным фокусным расстоянием, при компьютерном моделировании трехмерных объектов и сцен, или каким-то иным образом.
Обработка и анализ трехмерных изображений играют в настоящее время существенную роль во многих областях исследований, особенно в медицине и геологии. В настоящей работе рассмотрены вопросы расширения модели двумерного изображения [4] и применения ее к трехмерным изображениям, вопросы модификации операций частотной и пространственной фильтраций при переходе в 3D, вопросы сглаживания и декомпозиции изображений [5],[11], фильтрации помех, обнаружения контуров и объектов, а также вычислительные аспекты реализации некоторых алгоритмов для трехмерных изображений [16].
Большинство алгоритмов фильтрации при переходе от двумерных к трехмерным сигналам модифицируются сравнительно просто. Это будет показано на примерах наиболее распространенных алгоритмов, основанных на частотной и пространственной фильтрациях.
II. Особенности трехмерного изображения В дискретном виде 3D-изображение представляется массивом X = [xmnk] размерами MxNxK. Как и в 2D, значение каждого элемента xmnk есть квантованное на (xmax+1) градаций значение логарифма яркости (энергии) 0 < xmnk < xmax, которое для краткости будем называть просто яркостью. Дискретный элемент 3D-изображения принято называть вексель. Трехмерное изображение, отображающее некоторую сцену, можно рассматривать состоящим из плотно упакованных связных трехмерных областей (объектов), соответствующих деталям сцены. Областью или объектом будем называть максимальное по размеру связное множество элементов изображения, имеющих близкие, возможно плавно меняющиеся значения яркости. Области могут соприкасаться произвольным образом, в том числе одна область может быть полностью окружена другой. На границах соседних областей значения яркости должны заметно различаться. Не соприкасающиеся области могут иметь произвольные, в том числе и совпадающие яркости. Пространственные границы между соседними областями, как объектами различающейся яркости, называют контурами. 3D-изображения характеризуются
1
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
следующими свойствами:
— 3D-изображение — совокупность объектов, плотно заполняющих пространство изображения;
— контуры в 3D-изображении — пространственные границы между объектами;
— сечение 3D-изображения плоскостью и проекция его на плоскость любого направления дают двумерный сигнал со всеми свойствами обычного 2D-изображения.
III. Области анализа
Область анализа — подмножество исходных данных, используемое в оценке параметров. Методы, в которых для каждой точки (или малых фрагментов) изображения используются свои параметры обработки, определяемые по ограниченной и, как правило, центрированной в данной точке области анализа, называют локальными.
Рассмотрим связное множество элементов xijle Vd(xmnk) таких, которые отстоят от центрального элемента xmnk на расстояние не далее, чем d и вместе составляют фигуру некоторой задаваемой формы. При d < 2a3 множество Vd(xmnk), окружающее центральный элемент (воксель) xmnk, будем называть окрестностью и обозначать Vmnk, а при d >>1 — фрагментом и обозначать Wmnk. Отметим, что в зависимости от выполняемых операций сам центральный элемент xmnk может как принадлежать, так и не принадлежать Vd(x). Соответственно, операции вида
ymnk = f{xijl 1 xijl e Vd(xmnk) }, (1)
в которых результат в каждой точке (m,n,k) зависит лишь от значений элементов xijl, входящих в Vd(xmnk), называются локальными операциями.
При переходе из 2D в 3D варианты симметричных окрестностей и соседства элементов претерпевают следующие изменения. Окрестность из 2x2 элементов (4 пикселя) становится окрестностью из 2x2x2 элементов (8 вокселей), в которой каждый воксель соседствует с каждым. В двумерной окрестности из 3x3 элементов (9 пикселей), как известно, можно рассматривать два варианта соседства элементов: 4-соседство (только по сторонам пикселей) и 8-соседство (по сторонам и вершинам пикселей) [3]. Аналогом первого из них в 3D будет окрестность с 6-соседством вокселей (Рис. 1,а). Аналогом второго — окрестность с 26-соседством вокселей (Рис. 1,в). Возможен промежуточный вариант с 18-соседством вокселей (Рис. 1,б). Выбор варианта окрестности обычно определяется контекстом задачи и используемым алгоритмом.
б
в
Рис. 1. Окрестности и соседства вокселей в 3D: а) 6-соседство; б) 18-соседство; в) 26-соседство.
В некоторых случаях нас интересует не весь набор точек, попадающих в область анализа, а лишь некоторое его подмножество, включающее центральный элемент, которое будем называть областью принадлежности. Способ выбора области принадлежности зависит от
задачи; некоторые варианты рассмотрены в разделе VI.
IV. Двухмасштабная многокомпонентная модель
ТРЕХМЕРНОГО ИЗОБРАЖЕНИЯ
Для формулировки сведений об основных свойствах изображений — топологических (форм, размеров областей и контурных перепадов между ними) и статистических (взаимосвязи значений элементов) необходима соответствующая модель изображения. Чтобы быть полезной, она должна описывать свойства изображений на расстояниях, обусловленных особенностями задач, а также давать возможности для построения эффективных алгоритмов обработки и анализа изображений.
Статистические взаимосвязи элементов изображения, находящихся на больших расстояниях, существенно отличаются от аналогичных свойств близлежащих элементов и их не удается описать одними и теми же соотношениями. Для обычного двумерного изображения была разработана двухмасштабная многокомпонентная модель, достаточно хорошо описывающая взаимосвязи элементов и на малых расстояниях в несколько шагов дискретизации, и на больших — соразмерных размерам объектов изображения [4]. Она с успехом может быть перенесена на 3D-изображения.
Значения элементов 3D-изображения X = [xmnk] при этом будут представляться в виде суммы статистически независимых компонент:
xmnk Smnk ^ tmnk + ^mnk (2)
Первый член суммы — кусочно-гладкая компонента Smnk, определяющая уровни яркости протяженных областей изображения; tmnk — текстурно-детальная компонента, несущая информацию о текстуре и мелких деталях; £,mnk — шумовая компонента, определяемая шумами регистратора, аналого-цифрового преобразователя и др. Все компоненты предполагаются
независимыми и аддитивными, а tmnk и £,mnk — нормально распределенными и несмещенными.
A. Масштаб малого размера
На масштабе малого размера (масштабе элементов окрестности) рассматривается сравнительно небольшое связное множество элементов, расположенных на расстоянии нескольких шагов дискретизации. Как и в двумерной модели [4], элементы трехмерного изображения разделяются на два непересекающихся множества: попадающие на граничные участки (контурные) и не попадающие (внутренние), составляющие вместе полное изображение. Окрестность Vmnk элемента xmnk рассматривается как группа из R элементов xrmnk e Vmnk, r = 1,...,R, ближайших к xmnk, и попадающих в то же множество (контурное или внутреннее), что и элемент xmnk (Рис. 2).
Методом наименьших квадратов проводится гиперплоскость, наиболее близкая значениям элементов из Vmnk, составляющая с гиперплоскостью, ориентированной вдоль осей координат MNK, некоторый угол, величина и направление которого в точке (m,n,k) характеризуется вектором gmnk. В точке r окрестности проведенная гиперплоскость отличается от значения
2
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
xrmnk на случайную величину yrmnk . Такое представление
позволяет связать значения элементов окрестности
XLk ^ Vmnk фоРМУлой:
Xmnk ^mnk + Р gmnk + У mnk , (3)
где pmnk — значение проведенной гиперплоскости в центральной точке окрестности (m,n,k), pr — расстояние между центральным элементом Xmnk и xrmnk, grmnk — величина проекции gmnk на вектор из Xmnk в xrmnk, а yrmnk —
случайная величина.
Вводится понятие контурной маски E = [emnk]: emnk = 1 для контурных и emnk = 0 для внутренних элементов. Обозначая для контурных и внутренних элементов grmnk
через Vmnk и Vmnk , а Ymnk через Cmnk и hmnk
соответственно, представим grmnk и jrmnk в виде сумм
grmnk = emnk <nk +1-eLk )Ymn
и Ymnk emnk Zmnk +(1 emnk )nmnk ' В
результате получим формулу модели трехмерной окрестности:
Xmnk ^ mnk ~+emnk (ф mnk pr +C'mnkMKnWmnkpr Wmk) • (4)
Здесь Zmnk — стохастическое возбуждение в точке r
окрестности для контурных, а nrmnk — для внутренних
элементов. Случайные величины ф,^, Wmnk, Zmnk, и Pmnk считаются некоррелированными и несмещенными, а шумовые составляющие Zmnk, и pmnk — нормально распределенными. Проведенные эксперименты [4],[5] показывают, что значения дисперсий компонент ф и ф различаются в 10-100 раз.
B. Масштаб большого размера
На масштабе большого размера (масштабе объектов фрагмента) полагается, что гладкие составляющие Sv тех частей областей uv (v = 1,...,V), которые попадают во фрагмент Wmnk, задаются полиномом степени не выше, чем ю. Тогда Sii внутри фрагмента Wmnk может быть описана формулой:
V ю ю-p ю-p-q
Sj (Wmnk ) = Ё §Bv Ё Ё Ё аР^Р ^ (5)
v=1 p=0 q=0 r=0
здесь (i,j,l) — точка области u во фрагменте Wmnk; 5 v = 1, если точка (i,j,l) e u и 5 v = 0 в остальных
U J u
случаях. Добавляя в (5) текстурную tj и шумовую ^ составляющие, получим выражение для значения элементов внутри фрагмента:
^ ю-p ю-p-q
^ = Ё 5uv ЁЁ Ё ^РЛГ + j + ^
(6)
V
v=1 у p=0 q=0 r =0 J
Это общая формула модели, описывающей значения элементов областей внутри фрагмента. Компоненты tj
и считаются нормально распределенными, но
значения дисперсий tv вообще говоря различаются от области к области. Двумерный вариант изображения, содержащего нескольких областей, а также окрестность и фрагмент вокруг элемента Xmnk, показаны на Рис. 2. Для 3Б-изображения как окрестность с фрагментом, так и сами области будут трехмерными фигурами.
На многих реальных изображениях области, в
пределах типичного фрагмента анализа, имеют приблизительно постоянные яркости, не меняющиеся заметно. Поэтому во многих случаях допустимо выбрать минимальную степень полинома: ю = 0. Т огда
Sj (Wmnk) = Svmnk и (6) упрощается до следующего вида:
X
ijl
V
ЁЬ ( S'mnk + t, +Ц ) .
(7)
Это формула кусочно-постоянной модели фрагмента для представления участков областей изображения, попадающих во фрагмент Wmnk.
На масштабе большого размера обычно предполагается, что фрагмент Wmnk является прямоугольным параллелепипедом; это вызвано особенностями алгоритмической реализации методов обработки трехмерных изображений.
Рис. 2. Участок изображения; обрабатываемый элемент окружен границами окрестности и фрагмента.
V. Изменение частотных алгоритмов в 3D
Одним из эффективных подходов к обработке сигналов является применение методов фильтрации, использующих ортогональные преобразования: Фурье, Уолша-Адамара, косинусное, и др. [6]. Наиболее распространен класс преобразований, обеспечивающих разложение сигналов по гармоническим функциям, из которых важнейшим является преобразование Фурье; на его примере и покажем модификацию преобразования.
При переходе от 2D- к 3D-изображениям трехмерными становятся как пространственная, так и частотная области представления данных, т.е. вместо двумерного получаем трехмерное Фурье-преобразование; основные свойства преобразования Фурье при этом сохраняются.
Пусть f(x,y,z) — непрерывная функция трех
переменных x, y и z. Пара трехмерных непрерывных преобразований Фурье (прямое и обратное) задается следующими выражениями:
F(u, v, w) = J J J f (x, y, z)e-i2n(xu+yv+zw)dxdydz (8)
и f (x,y,z) = J J J F(u,v,w)ei2n(xu+yv+zw)dudvdw, (9)
где u, v и w — непрерывные частотные, а x, y и z — непрерывные пространственные переменные.
Трехмерное дискретное преобразование Фурье (ДПФ)
3
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
может быть записано в виде:
M -1 N-1 K-1
F(u,v, w) = XXS f (m,n,k)e-nM3m,M+vn/N+wk/K> , (10)
m=0 n=0 k=0
где f(m,n,k) — трехмерный массив размерами MxNxK. Обратное трехмерное дискретное преобразование Фурье будет иметь вид:
1 M -1 N-1 K-1
f (m,n,k) = —L_Y У У F(u,v,w) e2n(um/M+vn/N+wk/K) (11)
J MNKU v=0w=0
для 0 < m < M, 0 < n < N и 0 < k < K, а также 0 < u < M, 0 < v < N и 0 < w < K.
Как и в двумерном случае, частотная фильтрация трехмерного сигнала осуществляется выполнением прямого преобразования Фурье (10), модификацией полученного спектра F(u,v,w), и последующим обратным преобразованием (11). По аналогии с преобразованием Фурье, другие ортогональные преобразования при переходе к трехмерному сигналу также несложно видоизменяются добавлением размерности.
VI. Изменение пространственных алгоритмов в 3D Пространственная фильтрация обычно осуществляется локальными операторами согласно формуле (1) с определенными ограничениями на размеры локальной области анализа Vd(x), поскольку в случае функции f(v) в (1) общего вида, объем необходимых вычислений возрастает как минимум пропорционально числу элементов, попадающих в Vd(x). Окрестности, которые в 2D-случае содержат порядка (2d)2 пикселей, в 3D будут содержать уже порядка (2d)3 вокселей. Тем не менее для некоторых преобразований в трехмерном случае, как и в двумерном, удается построить быстрые алгоритмы.
A. Алгоритмы, использующие оценку среднего по фрагменту
Ряд алгоритмов обработки двумерных изображений, связанных со сглаживанием, выделением низко- или высокочастотных составляющих, улучшением изображения, обобщаются следующей формулой [5],[7]:
ymn f(xmn Smn,Dmn) + bSmn + c,
(12)
где xmn — значение центрального элемента, Smn — оценка сглаженного значения для точки (m,n) (например, значение арифметического среднего или медианы по фрагменту), Dmn — оценка дисперсии по фрагменту, f(u,v) — зависимость усиления контраста, b и c — параметры преобразования. Так, при f(u,v) = 0 и c = 0 получаем выделение низкочастотной составляющей (сглаживание), при f(u,v) = u, b = 0 и c = 0,5 — выделение высокочастотной составляющей, а при f(u,v) > u, b < 1 и c ~ (1-b)/2 — улучшение изображения [5],[7].
Преобразования вида (12) легко переносятся на 3D-изображения добавлением размерности, но с учетом того, что Smnk, а также u и v в f(u,v) должны определяться уже не по двумерному, а по трехмерному фрагменту (прямоугольному параллелепипеду) с центром в точке (m,n,k):
ymnk f(xmnk Smnk,Dmnk) + bSmnk + c. (13)
Модификации алгоритмов быстрого вычисления арифметического среднего или медианы Smnk, а также
дисперсии Dmnk рассмотрены в конце статьи.
B. Операторы контурных перепадов В основе большинства алгоритмов обнаружения контуров на изображении лежат операторы контурных перепадов, использующие оценки значений либо первой, либо второй производной [3],[8]. Для 2D-изображений контуры объектов суть линии, разделяющие плоские области; в случае, когда линии не имеют разрывов, получаем контурную карту изображаемой сцены.
В трехмерном пространстве контурное изображение представляет собой множество разделяющих объекты поверхностей произвольной формы, и здесь возможно возникновение различных конфигураций контуров, например, становятся допустимыми самопересечения и узлы. Сечение 3D-контурного изображения плоскостью представляет собой двумерную карту контуров. Однако такая карта не обязана совпадать с картой, получаемой проведением контуров по 2D-изображению, являющемуся сечением 3D-изображения той же плоскостью.
Рассмотрим модификации базовых операторов контурных перепадов на основе первой и второй производных при переходе от 2D- к 3D-изображениям в рамках рассмотренных выше трехмерных окрестностей.
1) На основе первой производной: оператор Робертса Оператор Робертса для двумерного изображения формулируется как сумма модулей разностей значений диагональных элементов по квадрату 2x2 пикселей: y(m,n)={lx(m,n)-x(m+1,n+1)l+lx(m+1,n)-x(m,n+1)l}/2. (14)
Для трехмерного изображения оператор Робертса [9] можно представить в виде суммы модулей разностей элементов в диагоналях куба из 2x2x2 пикселей: y(m,n,k) = (15)
={lx(m,n,k)-x(m+ 1,n+1,k+1)l+ lx(m+ 1,n,k)-x(m,n+ 1,k+1)l+ +lx(m,n+1,k)-x(m+1,n,k+1)l+lx(m,n,k+1)-x(m+1,n+1,k)l}/4.
2) На основе первой производной: оператор Собела Для двумерного изображения оператор Собела задается выражением:
y(m,n) = {lx(m-1,n-1)+2x(m-1,n)+x(m-1,n+1) -
- x(m+1,n-1)-2x(m+1,n)-x(m+1,n+1)l +
+ lx(m-1,n-1)+2x(m,n-1)+x(m+1,n-1) - (16)
- x(m+ 1,n-1 )-2x(m+1 ,n)-x(m+ 1,n+1) l}/8.
Чтобы построить оператор Собела [3] (и аналогичные
ему) для 3D-изображения, первоначально определим частные отклики как модули разностей значений элементов для каждого из направлений m, n и k, что соответствует выражению под одним из знаков модуля в (16). Для направления m такую зависимость ym(m,n,k) с точностью до коэффициентов a, b и c можно выразить формулой:
ym(m,n,k) =
=lax(m- 1,n- 1,k-1)+bx(m- 1,n- 1,k)+ax(m-1 ,n-1 ,k+1)+
+ bx(m-1,n,k-1) +cx(m-1,n,k) +bx(m-1,n,k+1)+
+ ax(m- 1,n+ 1,k-1)+bx(m- 1,n+ 1,k)+ax(m-1 ,n+1 ,k+1)-
- ax(m+ 1,n-1,k-1 )-bx(m+ 1,n-1 ,k)-ax(m+1,n- 1,k+1)-
- bx(m+1,n,k-1) -cx(m+1,n,k) -bx(m+1,n,k+1)-
- ax(m+ 1,n+1,k-1 )-bx(m+ 1,n+1 ,k)-ax(m+1,n+ 1,k+1)l/
/(4a+4b+c). (17)
Аналогичным образом записываются модули разностей yn(m,n,k) и yk(m,n,k) для составляющих n и k.
4
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
Трехмерный оператор, аналогичный оператору Собела, выразим через корень из суммы квадратов откликов по трем направлениям:
y(m,n,k)= {[ym(m,n,k)]2+ [y„(m,n,k)]2+[yk(m,n,k)]2}1/2. (18) Коэффициенты a, b и c рекомендуется выбирать как 1, 2 и 3 соответственно.
3) На основе второй производной: оператор Лапласа Оператор Лапласа (лапласиан) для двумерного
изображения задается следующей формулой [3]: y(m,n) =
{x(m-1 ,n-1)+x(m- 1,n)+x(m- 1 ,n+1)+x(m,n-1)+x(m,n+1) + +x(m+1,n-1)+x(m+1,n)+x(m+1,n+1)}/8 - x(m,n). (19)
Для 3D-изображения модификация лапласиана может быть выражена через набор элементов окрестности:
У mnk
Е
у QV xijleVmnk
- x„
(20)
где Vmnk — окрестность точки xmnk, не включающая саму центральную точку xmnk, QV — число точек в такой окрестности, {x^} — набор элементов окрестности Vmnk. В качестве Vmnk может выступать одна из окрестностей, показанных на Рис. 1; для окрестности (а) QV = 6, для (б) Qv = 18, для (в) Qv = 26.
4) Другие контурные операторы
Как и для двумерного случая, возможно построение и других детекторов, позволяющих обнаруживать контурные перепады, например, разности максимума и минимума значений элементов окрестности Vmnk:
x
ijl
ymnk = max {xijl I xijl e Vmnk} - min{ xtji I xj e Vmnk}.
C. Фильтрация импульсных помех Импульсными помехами называют сравнительно редкие искажения отдельных элементов изображения, когда значения помехи значительно отличаются от истинных значений сигнала и не коррелированны с ними. Модель искажения изображения импульсными помехами проста. Значение каждого из элементов xmnk с вероятностью p заменяется на случайное значение 4mnk независимо от значений остальных элементов. Обозначим через X’ = [x'mnk] исходное неискаженное изображение, через X = [xmnk] — искаженное импульсной помехой, а через Y = [ymnk] — результат фильтрации. Процесс искажения представится в виде:
Г xm nk с вероятностью (1 - p);
xmnk=
4 m
с вероятностью p.
(21)
Как правило полагается, что значения импульсных помех 4mnk распределены равномерно в диапазоне яркостей [0,xmax]. Фильтрация импульсных помех сводиться к обнаружению помех и коррекции искаженных отсчетов яркости.
Наиболее распространенные алгоритмы фильтрации основываются на предсказании значения элемента xmnk по окружающей его окрестности Vmnk. При этом используются локальные корреляционные связи близлежащих элементов изображения и учитывается, что шум пространственно декоррелирован. При обнаружении сравниваются наблюдаемое xmnk и предсказываемое xmnk значения. Если они отличаются более, чем на величину некоторого порога обнаружения
5, считается, что xmnk — помеха и осуществляется ее исправление на значение xmnk [10].
Общая формула вычисления xmnk в целом повторяет формулу (1):
xmnk = /{xijl 1 xijl e Vmnk}. (22)
Значение функции /{•} будем находить следующим образом. Выберем из набора элементов окрестности xijle Vmnk усеченное множество {xj} путем отбрасывания 2n точек с наименьшими и наибольшими значениями (n ~ QV/8). Множество точек {x'j} будет соответствовать так называемой области принадлежности. Проведем через точки {x'j} ближайшую гиперплоскость и возьмем ее значение в точке (m,n,k). Полученная величина и будет искомым предсказываемым значением xmnk.
Обобщенная формула фильтрации для большинства алгоритмов будет следующей:
Гxmnk , есЛИ Knk - xmnk 1 < 5; ^
ymnk = 1 ~ , - , . с (23)
если Ixmnk - xmnk 1 ^ 5.
mnk
При больших p, когда вероятность искажения пары соседних вокселей достаточно высока, рекомендуется итеративный процесс фильтрации: первоначально с большим значением порога 5, затем с постепенным его уменьшением. Рекомендуется для первой итерации применять окрестность на Рис. 1,а, а для последующих — окрестности на Рис. 1,б и (в).
D. Декомпозиция изображения
Для двумерного изображения вопросы декомпозиции, означающей в терминах модели (2) разделение изображения на сглаженную S и текстурно-детальную компоненты (t + 4), были рассмотрены в [5],[7],[11]. Суть изложенного в них алгоритма заключается в том, что для каждой точки изображения (m,n) на основе последовательного анализа окрестности Vmn и фрагмента Wmn, окружающих элемент xmn, выбирается часть точек Wmn в качестве области принадлежности, и значение Smn оценивается по выбранному множеству.
Алгоритм декомпозиции использует анализ распределений вероятностей значений элементов по окрестности Vmn и фрагменту Wmn. Пример фрагмента, охватывающего три области, а также распределение значений элементов по фрагменту (гистограмма) показаны на Рис. 3. Областью принадлежности в данном случае является множество точек объекта и.
Рис. 3. а) Область анализа (фрагмент); б) распределение вероятностей значений элементов по фрагменту.
Поскольку используемый подход основан на анализе распределений вероятностей, то он позволяет легко
5
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
перейти от 2D- к 3D-изображению — достаточно лишь вместо двумерных окрестности Vmn и фрагмента Wmn подставить в формулы трехмерные окрестность Vmnk и фрагмент Wmnk. В нотации [12], алгоритм декомпозиции трехмерного изображения формулируется следующим образом.
При заданных размерах txlxl окрестности Vmnk и LxLxL фрагмента Wmnk (l < L), центрированных в точке (m,n,k), ширине яркостных интервалов анализа А и А , а также ранговых параметрах nV < l3/2 и nW < L3/2 соответственно, значение сглаженной компоненты Smnk в (2) находится выполнением следующих операций для каждой точки (m,n,k) изображения:
1. Подсчитываются вероятности распределения (гистограммы) по трехмерным окрестности Н^ = { hVnk (i) } и фрагменту Н^к = { hWnk (i) } с центром в (m,n,k).
2. По гистограмме окрестности HmVnk и параметру nV находятся значения R = RV(nV/l3) и R = RV(1 - nV/l3); здесь ранговые значения R(x) определяются как решения
уравнения £ *=0 ^k (i) = -Г где hL (i) — гистограмма значений элементов в окрестности Vmnk. Промежуточное значение x находится сравнением -mnk с R1 и R2 :
XV = -mnk, если R1 < -mnk < R1 '; XV = R , если -mnk < R1V ;
V V V
и X = R2 , если Xmnk > R1 .
3. Из элементов Vmnk выбираются z таких xrmnk e V,
mnk
(r = 1,..., z), что попадают в интервал (XV - AV, XV + AV), где AV — его полуширина. По значениям Xrmnk из данного интервала, подсчитывается среднее [13]:
n л- , = 1 —
Xmnk A( V mnk , Xmnk ,
X„
= A(Vmnk , Xmnk , n , А ) = - Ё X
Z r = 1
r
mnk
(24)
где XV - AV < Xmnk < XV + AV.
4. Аналогично п. 2, по гистограмме фрагмента H'Wnk и
W r»W r>W/ W/т 3\
заданному n находятся значения R1 = R (n /L ) и rw = rw(1 - nW/L3). Значение XW находится сравнением
X„„k с RW и RW :
~W — X = X,
mnk, если RW < Xmrk < RV ; XW = RW ,
--RW , если X„„k > RW
mnk 1
если Xmnk < R1W ; и X'
5. Сглаженное значение Smnk находится как медиана по усеченной гистограмме фрагмента HmWnk — той ее части,
которая расположена в интервале (X - A , X + A ):
Smnk = med(Wm„k, Xmnk , nW, AW).
(25)
Рис. 4. Результат декомпозиции: исходное изображение (слева), сглаженная компонента S„n (справа) и графики отмеченной строки.
Полученное значение Smnk считается искомой сглаженной компонентой. Упрощенный вариант данного алгоритма был позже опубликован в [14] под названием билатеральная фильтрация. Пример декомпозиции двумерного изображения размерами 256x256 элементов с фрагментом сглаживания Wmnk размерами 15x15 элементов показан на Рис. 4.
E. Обнаружение объектов заданного объема В [12] показано, что изложенный выше алгоритм декомпозиции можно использовать для обнаружения объектов на изображении. Аналогично 2D-изображе-ниям, для которых решается задача обнаружения объектов по их площади, в трехмерной модификации ставится задача обнаружения объектов по их объему.
По аналогии с двумерной, трехмерная задача также допускает формулировку в трех вариантах: обнаружение объектов с объемом (т.е. числом элементов) N j меньше заданного T1, больше заданного T2, и обнаружение объектов, имеющих объем в интервале T1 < N] < T2.
1) Обнаружение объектов с N1 > T Предполагается, что изображение состоит из достаточно ровного фона (большая область U0), на которой имеется набор небольших областей U1,...,UJ, отстоящих друг от друга достаточно далеко, и можно выбрать некоторый размер фрагмента L (L3/2 > T) такой, что в любой фрагмент Wmnk попадает не более одной области c N] > T, либо несколько меньших, но при условии EN1 < T (U1 _ W). В п. 5 алгоритма декомпозиции (25) выберем nW= T, a R'W = R(T/L3) и R2W = R(1-T/L3). Обработкой изображения с данными значениями R1W и R2W получим:
ymnk Smnk, (26)
т.е. сглаженную компоненту исходного изображения, на которой остались области с N] > T, обнаруживаемые детектором со значением порога S(U0) ± 5, где S(U0) — средняя яркость фона, а 5 < mini{IS(U1) - S(U 0)l}; (S(U 1) — яркости соответствующих областей).
2) Обнаружение объектов с N1 < T
Сглаженная компонента Smnk, в (25) содержит лишь области с N1 > T, а области с N1 < T содержатся в разностной компоненте tmn = (Xmn - ^mn) - Smn. Детекция объектов с N1 < T достигается пороговым обнаружением в точках, где ltmnl > 5 (5 — порог обнаружения).
3) Обнаружение объектов с Tj < N] < T2 Возможны два варианта решения.
В первом случае сначала выберем nW = T1. Тогда сглаженная компонента Smnk в (25) будет содержать объекты с N] > T1. Осуществим повторную ее обработку алгоритмом (25) с nW = T2 (T2 > T1). Очевидно, во вновь полученной сглаженной компоненте S’mnk будут содержаться лишь объекты с N1 > T1. Взяв разность ymnk = ISmnk - S’mnkl получим сигнал, содержащий объекты в диапазоне T1 < N] < T2. Недостаток данного решения — алгоритм получается двухпроходовым.
Второй вариант. Обратим внимание, что при анализе гистограмм по окрестности и фрагменту используются два разных порога (nV и nW). Выберем размеры окрестности l и фрагмента L больше обычного — такими, чтобы l3 > 2T1 и L3 > 2T2. Задав R1V и RV, как
6
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
RV = RV(T1/l2) и RV = RV(1 - ТУ/2), после операции (24) будем иметь xmnk, которое уже не содержит области с N] < Т1. Далее в п. 5 алгоритма декомпозиции, при анализе HW, зададим RW и R2W как RW = RW(T2/L3) и rw = rw(1 - t2/l3). Получив значение Smnk в (25), возьмем разность ymnk = I xmnk - Smnk\, на которой объекты выделяются пороговым обнаружением.
Отметим, что как понятие площади в двумерном случае, так и понятие объема в трехмерном варианте алгоритма используются в несколько необычном смысле — как «локальный» объем, т.е. объем той части объекта, которая попадает внутрь фрагмента Wmnk.
VII. Некоторые вычислительные алгоритмы
A. Сумма по прямоугольному параллелепипеду Введем обозначение для суммы по фрагменту 2Б-изображения:
S
(ij)(mn)
m-1 n-1
EE x»,
т.е. S(j)(mn) — сумма значений элементов xuv, попадающих в прямоугольный фрагмент, диагональные точки которого имеют координаты (i,j) и (m-1,n-1). Обратим внимание, что фрагмент при этом не включает точку с координатами (m,n) и соответствующие ей строку и столбец. Аналогично для трехмерного изображения, S(ijl)(mnk) — сумма значений элементов xijl в прямоугольном параллелепипеде с координатами углов (i,j,l) и (m-1,n-1,k-1):
S
{ijl){mnk)
m-1 n-1 k-1
ZEE;
uvw
u=i v=j w=l
Для двумерного изображения классический способ вычисления суммы S(mn)(m+Hn+L) по скользящему прямоугольному фрагменту размерами HxL элементов при переходе от элемента (m,n) к элементу (m,n+1) сводится к формуле
S(m,n+1)(m+H,n+L+1) S (m,n)(m+H,n+L) S(m,n)(m+H,n+1) +S(m,n+L)(m+H,n+L+1), (26)
где последние два члена — суммы элементов по левому (удаляемому) и правому (добавляемому) столбцам фрагмента. Алгоритм требует 4 операции независимо от размеров фрагмента — 2 операции в выражении (26) и две операции на пересчет каждой из сумм по столбцу S(m,n)(m+H,n+V) при переходе от строки m к строке m+1. Дополнительно требуется N ячеек для хранения сумм по столбцам.
При переходе к 3Б-изображению, формула (26) будет модифицирована для скользящего прямоугольного параллелепипеда размерами HxLxJ:
S(m,n,k+1)(m+H,n+L,k+J+1) S(m,n,k)(m+H,n+L,k+J) S(m,n,k)(m+H,n+L,k+1)+ + S(m,n,k+J)(m+H,n+L,k+J+1), (27)
где последние два члена — суммы элементов по левой (удаляемой) и правой (добавляемой) граням
параллелепипеда. Данный алгоритм требует уже 6 арифметических операции независимо от размеров фрагмента — 2 операции в выражении (27), две операции на пересчет каждой из сумм по граням и две на пересчет сумм по столбцам. Кроме того, требуется K
ячеек для хранения массива сумм по граням и NxK ячеек для хранения сумм по столбцам.
Для двумерного изображения известен и другой алгоритм вычисления суммы по прямоугольнику произвольного размера. Пусть для каждой точки (m,n) подсчитаны суммы Smn = S(0,0)(m,n) по прямоугольнику с диагональными элементами x00 и xm-1,n-1. Тогда сумма S(ij)(mn) значений элементов внутри прямоугольника с угловыми координатами (i, j) и (m-1,n-1) будет:
S(ij)(mn) Smn Sm,j Si,
+ Sij
(28)
что в среднем для каждого элемента изображения требует 2 операции для вычисления суммы Smn и 3 операции для вычислений по формуле (28). Однако для хранения сумм Smn требуется уже MxN ячеек, равное размеру изображения. Преимуществом является то, что за те же 3 операции можно вычислять S(ij)(mn) для любых значений координат, а не только по скользящему фрагменту.
Алгоритм (28) также может быть модифицирован для трехмерного изображения. Пусть для каждой точки (m,n,k) подсчитаны суммы Smnk по прямоугольному параллелепипеду с диагональными элементами x000 и xm-1,n-1,k-1, т.е. Smnk = S(ooo)(mnk). Нетрудно показать, что в таком случае сумма S(ijl)(mnk) значений элементов внутри параллелепипеда с угловыми координатами (i, j, l) и (m-1,n-1,k-1) вычисляется при помощи следующей операции:
S(ijl)(mnk) = Smnk -Smjk -Smnl -Sink +Smjl +Sijk +Sinl -Sijl. (29)
Таким образом с учетом того, что для вычисления каждого из значений Smnk требуется 3 арифметических операции, для вычисления суммы S(ijl)(mnk) для каждого элемента трехмерного изображения потребуется в среднем 10 арифметических операций. Объем дополнительно требуемой памяти составит MxNxK ячеек с разрядностью, достаточной для значений сумм Smnk.
Аналогично можно находить дисперсии по фрагменту D(ijl)(mnk), вычисляя значения сумм квадратов S(mnk)(x2) для каждой из точек (m,n,k) и для прямоугольного параллелепипеда S(ijl)(mnk)(x2), и затем пользуясь формулой
D(ijl)(mnk) = {S(ijl)(mnk)(x ) - (S(ijl)(mnk)) }/N(ijl)(mnk), (30)
где S(ijl)(mnk)(x2) — сумма квадратов значений элементов, попадающих в параллелепипед, а общее число точек в параллелепипеде равно N(ijl)(mnk) = (m-i)x(n-j)x(k-l).
B. Порядковые статистики по прямоугольному параллелепипеду
Основой для вычисления порядковых статистик по фрагменту W как 2D- так и 3Б-изображения служит гистограмма распределения значений яркости по этому фрагменту hW(x), а также ее интегральная характеристика FW(x):
x
FW (x) = Z hW (i) ; FW(xmax) = NW, (3 1)
i=0
где xmax — максимально возможное значение яркости, а NW — число точек во фрагменте W. Порядковые статистики вида RW(n), где 0 < n < NW, представляют собой зависимость:
RW(n) = z, если FW(z-1) < n < FW(z). (32)
Алгоритм скользящего вычисления гистограммы по
7
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
фрагменту строится аналогично формулам (26) и (27), т.е. при смещении фрагмента к следующей точке производится удаление точек на одной грани фрагмента и добавление точек на противоположной грани [15]. В двумерном случае алгоритм скользящего вычисления гистограммы по фрагменту при переходе от точки (m,n) к соседней точке (m,n+1) требует в среднем 2H числа операций (H — число строк во фрагменте). В трехмерном случае при переходе от точки (m,n,k) к точке (m,n,k+1) потребуется уже 2HxL операций, где HxL — число точек в грани параллелепипеда, перпендикулярной направлению смещения K.
В трехмерном случае число требуемых операций растет пропорционально произведению HxL. Однако, если (HxL) > (xmax+1), то вместо операций со значениями отдельных точек выгодно заранее сформировать гистограммы для граней параллелепипеда hj)(mnk+1)(x), а
затем осуществлять операции вычитания и прибавления таких гистограмм:
hW (x) = hW (x) - hF (x) +
h(i,j,k+1Xm,n,k+J+1)(x) h(ijk )(m,n,k+J)(X) h(ijk )(m,n,k+1)(x)
+ h(i,j,k+J)(m,n,k+J+1)(X) ,
(33)
где J — размер фрагмента в направлении смещения K. Действия по формуле (33) требуют в среднем 2(xmax+1) арифметических операций на одну точку изображения независимо от размера параллелепипеда. Для пересчета гистограмм по граням параллелепипеда hfjk)(mnk+1)(x)
нужно дополнительно в среднем 2L операций на точку и (xmax+1)xK ячеек памяти для хранения K гистограмм.
Число операций (xmax+1), требуемое для прибавле-ния/вычитания каждой из гистограмм по граням параллелепипеда hFjk)(mnk +1)(x), можно уменьшить, если
воспользоваться пространственной корреляцией обрабатываемых данных. На участках с медленным изменением яркости, которых на реальных изображениях обычно большинство, размах значений элементов сравнительно невелик — в несколько раз меньше полного диапазона в (xmax+1) градаций. Добавив к каждой из гистограмм hFjk)(mnk +1)(x) по 2 ячейки для
запоминания минимального и максимального значений распределения, и, соответственно, обрабатывая лишь указываемый диапазон градаций, удается дополнительно в несколько раз сократить общее число операций.
VIII. Заключение
Рассмотрены вопросы модификации двумерных модели и алгоритмов обработки изображений в применении к трехмерным изображениям. Показаны варианты изменения области анализа, алгоритмов фильтрации, декомпозиции изображений и обнаружения объектов, которые сравнительно несложно модифицируются при переходе от 2D- к 3Б-изображениям. Предложены вычислительные алгоритмы, сокращающие количество операций при определении значений среднего и порядковых статистик по скользящему фрагменту для 3D- изображения.
Библиография
[1] Toriwaki J., Yoshida H. Fundamentals of Three-Dimensional Digital Image Processing. New York: Springer, 2009.
[2] Красильников Н. Цифровая обработка 2D и 3D изображений. СПб.: BHV-Петербург, 2011.
[3] Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2012.
[4] Чочиа П.А. “Двухмасштабная модель изображения”. В кн.
Кодирование и обработка изображений. М.: Наука, 1988, С. 69-87.
[5] Чочиа П.А. Обработка и анализ изображений на основе двухмасштабной модели: Препринт ИППИ АН СССР. М.: ВИНИТИ, 1986.
[6] Ахмед Н., Рао К. Ортогональные преобразования при обработке цифровых сигналов. М.: Связь, 1980.
[7] Chochia P.A. “Image Enhancement Using Sliding Histograms”. Computer Vision, Graphics, Image Processing, 1988, vol. 44, no. 2, pp. 211-229.
[8] Прэтт У. Цифровая обработка изображений. М.: Мир, 1982. Т. 1, 2.
[9] Робертс Л. “Автоматическое восприятие трехмерных объектов”. В кн. Интегральные роботы. М.: Мир, 1973, С. 162-208.
[10] Чочиа П.А. “Цифровая фильтрация импульсных помех на телевизионных изображениях”. Техника средств связи: сер. Техника телевидения, 1984, вып.1, C. 26-36.
[11] Чочиа П.А. “Сглаживание изображения при сохранении контуров”. В кн. Кодирование и обработка изображений. М.: Наука, 1988, С. 87-98.
[12] Чочиа П.А. “Некоторые алгоритмы обнаружения объектов на основе двухмасштабной модели изображения”. Информационные процессы, 2014, Т. 14, № 2, С. 117-136.
[13] Lee J.S. “Digital Image Smoothing and the Sigma Filter”. Computer Vision, Graphics, Image Processing, 1983, vol. 24, no. 2. pp. 255-269.
[14] Tomasi C., Manduchi R. “Bilateral filtering for gray and color images”. Proc. IEEE 6th Int. Conf. on Computer Vision, Bombay, India, 1998, pp. 839-846
[15] Чочиа П.А. “Параллельный алгоритм вычисления скользящей гистограммы”. Автометрия, 1990, № 2, С. 40-44.
[16] Чочиа П.А. “Модификация модели и алгоритмов обработки при переходе от двумерных к трехмерным изображениям.” В кн. IX
Международная научно-практическая конференция “Современные информационные технологии и ИТ-образование." Сборник избранных трудов. М.: МГУ, 2014, С. 820-833.
8
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 2, no. 11, 2014
Three-dimensional and two-dimensional images: modification of image model, analysis area processing algorithms
P.A. Chochia
Abstract - The specificities of three dimensional images are formulated. The adaptation of analysis area and two-scale image model to 3D-images is studied. The modifications of various image processing algorithms to 3D-images are demonstrated. Fast algorithms for calculating local average and order statistics in the moving window for 3D-images are proposed.
Keywords — Image processing, three-dimensional image, image model, image processing algorithm, analysis area, fast algorithms.
9
