CC BY
29
Чочиа П.А.
Институт проблем передачи информации им. А.А. Харкевича Российской академии наук (ИППИ РАН), г. Москва, к.т.н., старший научный сотрудник, chochia @iitp . ru
Модификация модели и алгоритмов обработки при переходе от двумерных к трехмерным изображениям
КЛЮЧЕВЫЕ СЛОВА
Обработка изображений, трехмерное изображение, модель изображения, большие данные.
АННОТАЦИЯ
Рассматриваются вопросы модификации двухмасштабной модели и алгоритмов обработки при переходе от двумерных к трехмерным изображениям. Показываются способы сокращения объема вычислений при обработке больших данных.
1. Особенности трехмерного изображения
Под термином трехмерное или 3D-изображение зачастую понимаются совершенно различные виды данных [1,2]. Ниже под трехмерным изображением (непрерывным или дискретным) будут пониматься исключительно изображения, задаваемые функцией трех координат, которая гомеоморфно отображает некоторый объемный участок трехмерного пространства, включая все содержащиеся в нем объекты. 3D-изображение можно рассматривать как расширение 2D-изображения путем добавления еще одного пространственного измерения. Предполагается, что значения точек исходного пространства с координатами (m,n,k) описываются функцией f(m,n,k), имеющей свойства плотности некоторой выбранной физической характеристики.
В дискретном виде трехмерное изображение представляется массивом X = [xmnk] размерами MxNxK. Как и в 2D-изображении, значение каждого элемента xmnk есть квантованное на (xmax+1) градаций значение логарифма яркости (энергии) 0 0 xmnk 0 xmax, которое для краткости будем называть просто яркостью. Дискретный элемент 3D-изображения принято называть воксель (англ. voxel — volumetric pixel). Трехмерное изображение, отображающее некоторую сцену, можно рассматривать состоящим из плотно упакованных связных трехмерных областей (объектов), соответствующих деталям сцены. Областью или объектом будем называть максимальное по размеру связное множество элементов изображения, имеющих близкие, возможно плавно меняющиеся значения яркости. Области могут соприкасаться произвольным образом, в том числе одна область может быть полностью окружена другой. На границах соседних областей значения яркости должны заметно различаться. Не
соприкасающиеся области могут иметь произвольные, в том числе и совпадающие яркости. Пространственные границы между соседними областями, как объектами различающейся яркости, называют контурами. Трехмерные изображения характеризуются следующими свойствами:
— в трехмерном изображении объект (^-объект) — пространственная фигура, а само 3D-изображение — пространственно ограниченная совокупность плотно упакованных 3D-объектов, заполняющих отображаемый участок пространства;
— контуры в 3D-изображении — пространственные границы между объектами;
— сечение 3D-изображения плоскостью любого направления, а также проекция 3D-изображения на плоскость любого направления дают двумерный сигнал со всеми свойствами обычного 2D-изображения.
1.1. Локальные операции, окрестности и соседство элементов
Рассматриваемая ниже двухмасштабная модель изображения является по сути локальной, т.е. такой, которая для каждого произвольно взятого элемента изображения Хтпк описывает его взаимосвязи лишь с некоторым пространственно ограниченным множеством элементов Vd(Хmnk) или просто Vd(x), которые его окружают и находятся от центрального элемента не далее расстояния d. Множество элементов Vd(x), окружающее центральным элемент (воксель) Хтпк, при dЦ 23 будем называть окрестностью и обозначать далее как Утпк, а при d Ш 1 — фрагментом и обозначать Штпк- Отметим также, что в зависимости от выполняемых операций сам центральный элемент Хтпк может как принадлежать, так и не принадлежать Vd(x). Соответственно, операции вида
Утпк = ЛХд1 | Х^ □ У^Хтпк)}, (1)
в которых результат в каждой точке (т,п,к) зависит лишь от значений элементов Хщ, входящих в Vd(Хmnk), называются локальными операциями.
При переходе из 2D в 3D, варианты симметричных окрестностей и соседства элементов претерпевают следующие изменения. Окрестность из 22 элементов (4 пикселя) становится окрестностью из 222 элементов (8 вокселей), в которой каждый воксель соседствует с каждым. В двумерной окрестности из 33 элементов (9 пикселей), как известно, можно рассматривать два варианта соседства элементов: 4-соседство (только по граням пикселей) и 8-соседство (по граням и вершинам пикселей) [3]. Аналогом первого из них в 3D будет окрестность с 6-соседством вокселей (Рис. 1,а). Аналогом второго — окрестность с 26-соседством вокселей (Рис. 1,в). Возможен промежуточный вариант с 18-соседством вокселей (Рис. 1,б). Выбор варианта окрестности обычно определяется контекстом задачи и используемым алгоритмом.
Рис. 1. Варианты окрестностей и соседства элементов (вокселей) в 3D: а) 6-соседство вокселей; б) 18-соседство вокселей; в) 26-соседство вокселей
2. Двухмасштабная модель трехмерного изображения
Аналогично модели обычного 2D-изображения [4], 3D-изображение также удается описать двухмасштабной моделью, представляющей статистические взаимосвязи элементов как на масштабе малого размера в несколько шагов дискретизации, так и на масштабе большого размера, соразмерного объектам изображения. Значения элементов 3D-изображения X = [Хтпк] при этом представляются в виде суммы статистически независимых компонент:
Хтпк = Smnk + tmnk + ^тпк. (2)
Первый член суммы — кусочно-гладкая компонента Smnk, определяющая средние уровни яркости протяженных областей изображения; тк — текстурно-детальная компонента, несущая информацию о текстуре и мелких деталях; ^тпк — шумовая компонента, определяемая шумами регистратора, аналого-цифрового преобразователя и др. Все компоненты предполагаются независимыми и аддитивными, а тк и ^тпк — нормально распределенными и несмещенными.
2.1. Масштаб малого размера
На масштабе малого размера (масштабе элементов окрестности) рассматривается сравнительно небольшое связное множество элементов, расположенных на расстоянии нескольких шагов дискретизации. Как и в двумерной модели [4], элементы трехмерного изображения разделяются на два непересекающихся множества: попадающие на контурные участки (контурные) и не попадающие (фоновые), составляющие вместе полное изображение. Вводится понятие окрестности Утпк элемента Хтпк как группы из R элементов е Vmnk, г = 1,...Д, ближайших к Хтпк, и попадающих в то же множество (контурное или фоновое), что и центральный элемент Хтпк.
Методом наименьших квадратов проводится гиперплоскость, наиболее близкая значениям элементов из Утпк, составляющая с гиперплоскостью, ориентированной вдоль осей координат MNK, некоторый угол, величина и направление которого в точке (т,п,к) характеризуется вектором gmnk. В точке г окрестности проведенная гиперплоскость отличается от значения хгттк на случайную величину чгтпк. Такое
представление позволяет связать значения элементов окрестности X mnk еVmnk формулой:
Xmnk ^ mnk + Р gmnk + У mnk , (3)
где ^mnk — значение проведенной плоскости в центральной точке окрестности (m,n,k), pr — расстояние между центральным элементом Xmnk и хГтл, grmnk — величина проекции gmnk на вектор из Xmnk в хгтЛ, а ym* — случайная величина.
Вводится понятие контурной маски E = [emnk]: emnk = 1 для контурных и emnk = 0 для фоновых элементов. Обозначая для контурных и фоновых
r r r r у r r
элементов gmnk через jmnk и у mnk, а у mnk через Z mnk и 4mnk соответственно, представим grmnk и у1^ в виде сумм grmnk = егтЛ ф!, +(1 - ^л Wmnk и fmnk = ernnkz rmnk +(1 - Xmnk ■ В результате получим формулу модели трехмерной окрестности:
Xmnk = ^mnk + emnk (ф mnkР + Zmnk ) + (1 — emnk )(у mnkР + Пmnk ) ■ (4)
Здесь Zrm„k — стохастическое возбуждение в точке r окрестности для контурных, а nrmnk — для фоновых элементов. Случайные величины фmnk, ^mnk, Zmnk, и nmnk считаются некоррелированными и несмещенными, а шумовые составляющие Zmnk, и ^mnk — нормально распределенными. Эксперименты [4,5], показывают, что значения дисперсий компонент ф и ф различаются в 10-100 раз.
2.2. Масштаб большого размера
На масштабе большого размера (масштабе объектов фрагмента) полагается, что гладкие составляющие S тех частей областей uv (v = 1,...,V), которые попадают во фрагмент Wmnk, задаются полиномом степени не выше, чем Тогда Sj внутри фрагмента Wmnk может быть описана
формулой:
V 00 OO-JiOO-i-i
WJ=Z5,ZZ Z
v-l p-fi i-0 r-0 (5)
здесь (i,j,l) — точка области uv во фрагменте Wmnk; = 1, если точка (Щ) О uv и = 0 в остальных случаях. Добавляя в (5) текстурную tiji и шумовую составляющие, получим выражение для значения элементов внутри фрагмента:
оо OO-JiOO-Ji-tf
(6)
т-1 Г-О
Это общая формула модели, описывающей значения элементов областей внутри фрагмента. Значения ^ и ^ считаются нормально распределенными, но значения дисперсий V вообще говоря различаются от области к области.
На многих реальных изображениях области, в пределах типичного фрагмента анализа, имеют приблизительно постоянные яркости, не
меняющиеся заметно. Поэтому во многих случаях допустимо выбрать минимальную степень полинома: ы = 0. Тогда S¡]1 (Жтпк) = Svmnk и (6) упрощается до следующего вида:
V
4=2 \ (^ +4 ■
(7)
Это формула кусочно-постоянной модели фрагмента для представления участков областей изображения, попадающих во фрагмент
Ж„пк-
На масштабе большого размера обычно предполагается, что фрагмент Штпк является прямоугольным параллелепипедом; это вызвано особенностями алгоритмической реализации методов обработки трехмерных изображений.
3. Изменение алгоритмов фильтрации при переходе в 3D
Большинство алгоритмов фильтрации при переходе от двумерных к трехмерным сигналам модифицируются сравнительно просто. Покажем это на примере наиболее распространенных алгоритмов, основанных на частотной и пространственной фильтрациях.
3.1. Алгоритмы частотной фильтрации
Одним из эффективных подходов в обработке сигналов является применение глобальных методов фильтрации, использующих ортогональные преобразования: Фурье, Уолша-Адамара, косинусное, и пр. [6]. Наиболее распространен класс преобразований, обеспечивающих разложение сигналов по гармоническим функциям, из которых важнейшим является преобразование Фурье; на его примере и покажем модификацию преобразования.
При переходе от 2D- к 3D-изображениям трехмерными становятся как пространственная, так и частотная области представления данных, т.е. вместо двумерного получаем трехмерное Фурье-преобразование; все основные свойства преобразования Фурье при этом сохраняются.
Пусть f(x,y,z) — непрерывная функция трех переменных х, у и z. Пара трехмерных непрерывных преобразований Фурье (прямое и обратное)
задается следующими выражениями:
00 00 00
-00 -00 -00
00 00 00
-00 -00 -00 РЗ
где и, V и w — непрерывные частотные, ах, у и z — непрерывные пространственные переменные. Трехмерное дискретное преобразование Фурье (ДПФ) может быть записано в виде:
т=0 н=0 (10)
где f[m,n,k) — трехмерный цифровой массив размерами MxNxK. Обратное
трехмерное дискретное преобразование Фурье будет иметь вид:
1 М-1 N-1 К-1
f{m,n,k) = —— ^ ^ ^ F(u,v,w)en<umlM+vr!lN+wkSK)
для 0 < m < M, 0 < n < N и 0 < k < K, а также 0 < u < M, 0 < v < N и 0 < w < K.
Как и в двумерном случае, частотная фильтрация трехмерного сигнала производится выполнением прямого преобразования Фурье (10), модификацией полученного спектра F(u,v,w), и последующим обратным преобразованием (11). По аналогии с преобразованием Фурье, другие ортогональные преобразования при переходе к трехмерному сигналу также несложно видоизменяются добавлением размерности.
3.2. Алгоритмы пространственной фильтрации
Пространственная фильтрация обычно осуществляется локальными операторами согласно формуле (1) с определенными ограничениями на размеры локальной области анализа Vd(x), поскольку в случае функции Д-/) в (1) общего вида объем необходимых вычислений возрастает как минимум пропорционально числу элементов, попадающих в Vd(x). Окрестности, которые в 2D случае содержат порядка (2d)2 пикселей, в 3D будут содержать уже порядка (2d)3 вокселей. Тем не менее для некоторых преобразований в трехмерном случае, как и в двумерном, удается построить быстрые алгоритмы.
3.2.1. Алгоритмы, использующие оценку среднего по фрагменту
Ряд алгоритмов обработки двумерных изображений, связанных со сглаживанием, выделением низко- или высокочастотных составляющих, улучшением изображения, обобщаются следующей формулой [5,7]:
ymn = f(xmn — Smn,Dmn
) + bSmn + С, (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 Q 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 рассмотрены в конце статьи.
3.2.2. Модификация операторов контурных перепадов
В основе большинства алгоритмов обнаружения контуров на изображении лежат операторы контурных перепадов, использующие оценки значений либо первой, либо второй производной [3,8]. Для 2D-изображений контуры объектов суть линии, разделяющие плоские области; в случае, когда линии не имеют разрывов, получаем контурную карту изображаемой сцены.
В 3D-пространстве контурное изображение представляет собой множество разделяющих объекты поверхностей произвольной формы, и здесь возможно возникновение самых различных конфигураций контуров, например, становятся допустимыми самопересечения и узлы. Сечение 3D-контурного изображения плоскостью представляет собой двумерную карту контуров. Однако такая карта не обязана совпадать с картой, получаемой проведением контуров по 2D-изображению, являющемуся сечением 3D-изображения той же плоскостью.
Рассмотрим модификации базовых операторов контурных перепадов на основе первой и второй производных при переходе от 2D- к 3D-изображениям в рамках трехмерных окрестностей, рассмотренных выше.
На основе первой производной: оператор Робертса
Оператор Робертса для двумерного изображения описывается формулой:
у(т,п) = {|х(т,п) - х(т+1,п+1)| + |х(т+1,п) - х(т,п+1)|}/2, (14)
т.е. как сумма модулей разностей значений диагональных элементов по квадрату 2x2 пикселей. Для 3D-изображения оператор Робертса [9] будет представлять сумму модулей разностей элементов в диагоналях куба из 2x2x2 пикселей:
у(т,п,к) = {|х(т,п,к) - х(т+1,п+1,к+1)| + |х(т+1,п,к) - х(т,п+1,к+1)| +
+ |х(т,п+1,к) - х(т+1,п,к+1)| + |х(т,п,к+1) - х(т+1,п+1,к)|}/4. (15)
На основе первой производной: оператор Собела
Оператор Собела для двумерного изображения описывается формулой:
у(т,п) = {|х(т-1,п-1) + 2х(т-1,п) + х(т-1,п+1) -
- х(т+1,п-1) - 2х(т+1,п) - х(т+1,п+1)| +
+ |х(т-1,п-1) + 2х(т,п-1) + х(т+1,п-1) - (16)
- х(т+1,п-1) - 2х(т+1,п) - х(т+1,п+1)|}/8.
Чтобы задать оператор Собела [3] (и аналогичные ему) для 3D-изображения, первоначально определим частные отклики как модули разностей значений элементов для каждого из направлений т, п и к, что соответствует выражению под одним из знаков модуля в (16). Для направления т такую зависимостьут(т,п,к) с точностью до коэффициентов а, Ь и с можно выразить формулой:
ут(т,п,к) = |ах(т-1,п-1,к-1) + Ьх(т-1,п-1,к) + ах(т-1,п-1,к+1) +
+ Ьх(т-1,п,к-1) + сх(т-1,п,к) + Ьх(т-1,п,к+1) +
+ ах(т-1,п+1,к-1) + Ьх(т-1,п+1,к) + ах(т-1,п+1,к+1) -
- ах(ш+1,п-1,к-1') - Ьх{ш+1,п-1,к) - ax(m+1,n-1,k+1) - (17)
- Ьх(т+1,п,к-1) - сх(т+1,п,к) - Ьх(т+1,п,к+1) -
- ax(m+1,n+1,k-1) - Ьх(т+1,п+1,к) - ax(m+1,n+1,k+1)|/(4a+4b+c). Аналогичным образом записываются модули разностей уп(т,п,к) и
Ук(т,п,к) для составляющих пик. Трехмерный оператор, аналогичный оператору Собела, выразим через корень из суммы квадратов откликов по трем направлениям:
у(т,п,к) = {[ут(т,п,к)]2 + [у„(т,п,к)]2 + [ук(т,п,к)]2}1/2. (18)
Коэффициенты а, Ь и с рекомендуется выбирать как 1, 2 и 3 соответственно.
На основе второй производной: оператор Лапласа Оператор Лапласа для двумерного изображения описывается формулой:
у(т,п) = {х(т-1,п-1) + х(т-1,п) + х(т-1,п+1) + х(т,п-1) + х(т,п+1) +
+ х(т+1,п-1) + х(т+1,п) + х(т+1,п+1)}/8 - х(т,п). (19)
Модификацию оператора Лапласа [3] для 3D-изображения можно выразить как
УмнК-
1 ^ бг * 'У
\
■ ТГ
(20)
где Vmnk — окрестность точки хтпк, не включающая саму центральную точку хтпк, Qv — число точек в такой окрестности, {х^} — набор элементов окрестности Vmnk■ В качестве Vmnk может выступать одна из окрестностей, показанных на Рис. 1; для окрестности (а) Qv = 6, для (б) Qv = 18, для (в) Qv = 26.
Другие контурные операторы
Как и для двумерного случая, возможно построение и других детекторов, позволяющих обнаруживать контурные перепады, например, разности максимума и минимума значений элементов окрестности Vmnk:
у(т,п,к) = тах{хУ1 | х^ □ Vmnk} - тт{ х^ | хП1 □ Vmnk}■
3.2.3. Фильтрация импульсных помех
Импульсными помехами называют сравнительно редкие искажения отдельных элементов изображения, когда значения помехи значительно отличаются от истинных значений сигнала и не коррелированны с ними. Модель искажения изображения импульсными помехами несложна. Значение каждого из элементов хт^ изображения с вероятностью р заменяется на случайное значение ^к независимо от значений остальных элементов. Обозначим через X = [х т^] исходное неискаженное изображение, через X = [хт^] — искаженное импульсной помехой, а через У = \утпк] — результат фильтрации. Процесс искажения представится в виде: Г^ с вероятностью для неискаженного элемента;
Л ТН ТН — |
с вероятностью^ для элемента, искаженного помехой. (21)
Как правило полагается, что значения импульсных помех ^к распределены равномерно в диапазоне яркостей [0, хтах]. Фильтрация
импульсных помех сводиться к обнаружению помех и коррекции искаженных отсчетов яркости.
Наиболее распространенные алгоритмы фильтрации основываются на предсказании значения элемента хтпк по окружающей его окрестности Утпк. При этом используются локальные корреляционные связи близлежащих элементов изображения и учитывается, что шум пространственно декоррелирован. При обнаружении сравниваются наблюдаемое Хтпк и предсказываемое хтпк значения. Если они отличаются более, чем на величину некоторого порога обнаружения б, считается, что Хтпк — помеха и осуществляется ее исправление на значение хтпк. Общая формула вычисления хтпк при этом в целом повторяет формулу (1):
ХтА = АХф | Хф □ Утпк}- (22)
В качестве функции Д-} в (22) могут выступать различные линейные или нелинейные операторы — как каузальные, так и некаузальные [10]. Обобщенная формула фильтрации для большинства алгоритмов будет следующей:
? I им Л- I ?
Хин*.> еСЛН ГииЛ _ — (23)
Как было показано в [10], наилучшие результаты фильтрации импульсных помех достигаются алгоритмами, основанными на использовании в качестве Д^} в (22) значения медианы по симметричной окрестности обрабатываемого элемента: ДУтпк) = med{xг:/7 □ Утпк}. При больших р, когда вероятность искажения пары соседних вокселей достаточно высока, рекомендуется итеративный процесс фильтрации: первоначально с большим значением порога б, затем с постепенным его уменьшением. Рекомендуется для первой итерации применять окрестность на Рис. 1,а, а для последующих — окрестности на Рис. 1,б и (в).
3.2.4. Декомпозиция изображения
Вопросы декомпозиции, означающей в терминах модели (2) разделение изображения на сглаженную 51 и текстурно-детальную компоненты (Ь + были рассмотрены в [5,7,11,12]. Согласно изложенному в них алгоритму, для каждой точки изображения (т,п) производится последовательный анализ изображения сначала по окрестности Утп, а затем по фрагменту Жтп, окружающим центральный элемент Хтп. При этом по сути анализируются лишь распределения значений элементов по окрестности Утп и фрагменту Жтп.
Такой подход позволяет перейти от 2D- к 3D-изображению — достаточно вместо двумерных окрестности Утп и фрагмента Ж™ подставить в формулы трехмерные окрестность Утпк и фрагмент Жтпк. В нотации [12], алгоритм декомпозиции трехмерного изображения формулируется следующим образом.
При заданных размерах 1x1x1 окрестности Утпк и LxLxL фрагмента Жтпк
(/ < V), центрированных в точке (т,п,к), ширине яркостных интервалов анализа ДУ и ДЖ, а также ранговых параметрах пУ < /3/2 и пЖ < ¿3/2 соответственно, значение сглаженной компоненты Smnk в (2) находится выполнением следующих операций для каждой точки (т,п,к) изображения:
1. Подсчитываются гистограммы (вероятности распределения) по трехмерным окрестности Н^ = { ^ (г)} и фрагменту Н1Л = { ЪшпЛ (г)} с центром в (т,п,к).
2. По гистограмме окрестности Н^пк и заданному значению пУ находятся ранговые параметры ^ = /?У(>7У//3) и ^ = /?у(1 - пу//3); здесь /?(х)
определяется как решение уравнения —, где —
гистограмма значений элементов в окрестности Утпк. Сравнением значения центрального элемента Хтпк с и определяется промежуточное усеченное значение хГ: хГ = Хтпк, если < Хтпк < ^; хГ = ^, если Хтпк < ^Г; и хГ = ^, если Хтпк > К.
3. Из элементов окрестности Утпк выбираются z значений хгтпк □ Утпк (г = 1,..., z), попадающих в интервал (хГ - ДУ, хГ + ДУ), где ДУ — его полуширина. По значениям хгтпк, попадающим в данный интервал, подсчитывается среднее [13]:
2 ^ (24)
4. Аналогично п. 2, по гистограмме фрагмента Н1"пк и заданному пЖ находятся ранговые значения В" = ИЖ(пш/Ь3) и К" = ЯЖ(1 - nW/L3). Сравнением значения хтпк с и Щ определяется значение х": х" = хтпк, если < хтпк < ; ХЖ = , если хтпк < ; и х" = ^ , если хтпк > ^ .
5. Сглаженное значение Smnk находится по гистограмме фрагмента Н1"пк как медиана значений Нт1пк, попадающих в интервал (х" - ДЖ, х" + ДЖ):
Smnk = те4(Жтпк, хтпк,л", А"). (25)
Полученное значение Smnk считается искомой сглаженной компонентой. Упрощенный вариант данного алгоритма позже был опубликован в [14] под названием билатеральная фильтрация.
3.2.5. Обнаружение объектов заданного объема
В [12] показано, что изложенный выше алгоритм декомпозиции можно использовать для обнаружения объектов на изображении. Аналогично 2D-изображениям, для которых решается задача обнаружения объектов по их площади, в трехмерной модификации ставится задача обнаружения объектов по их объему.
Данная задача также формулируется в трех вариантах: обнаружение объектов с объемом (т.е. числом элементов) Q, больше заданного Т меньше заданного Т, и обнаружение объектов, имеющих объем в интервале
Т < Qi < Т2. Алгоритм обнаружения для 3D-изображений принципиально не отличается от алгоритма для 2D, который изложен в [12]; достаточно лишь в обозначения координат добавить третье измерение, а плоский прямоугольный фрагмент Штп заменить на прямоугольный параллелепипед Нт^. Как понятие площади в двумерном случае, так и понятие объема в трехмерном варианте алгоритма используется несколько в необычном смысле — как «локальный» объем, т.е. объем той части объекта, которая попадает внутрь фрагмента Нтпк-4. Некоторые вычислительные алгоритмы 4.1. Сумма по прямоугольному параллелепипеду Введем обозначение для суммы по фрагменту 20-изображения:
гп-1 п-1 ^СУХти) = >
т.е. S(iJ■)(mn) — сумма значений элементов х^, попадающих в прямоугольный фрагмент, диагональные точки которого имеют координаты (у и (т-1,п-1). Обратим внимание, что фрагмент при этом не включает точку с координатами (т,п) и соответствующие ей строку и столбец. Аналогично для трехмерного изображения, Smmnk) — сумма значений элементов Ху в прямоугольном параллелепипеде с координатами диагональных углов (iJ,I) и (т-1,п-1,к-1)\
т-1 и-1 к-1
о _
(!]1)(тпк) . , . , . , "^ЦУ^ 11=1 У?=1
Для 2D-изображения классический способ вычисления суммы S(mn) [т+нп+1) по скользящему прямоугольному фрагменту размерами Ях£ элементов при переходе от элемента (т,п) к элементу (т,п+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+н,n+l) при переходе от строки т к строке т+1. Дополнительно требуется N ячеек для хранения сумм по столбцам.
При переходе к 3D-изображению, формула (26) будет модифицирована для скользящего прямоугольного параллелепипеда размерами Нх£х/:
S (m,n,к+1)(m+H,n+L,к+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), две операции на пересчет каждой из сумм по граням и две на пересчет сумм по столбцам. Кроме того, требуется N ячеек
для хранения массива сумм по граням и N*K ячеек для хранения сумм по столбцам.
Для двумерного изображения известен и другой алгоритм вычисления суммы по прямоугольнику произвольного размера. Пусть для каждой точки (m,n) подсчитаны суммы Smn по прямоугольнику с диагональными элементами xoo и Xm-i,n-i, т.е. Smn = S(o,o)(m,n). Тогда сумма S^mn) значений элементов внутри прямоугольника с угловыми координатами (i,j) и (m-1,n-1) будет:
S(ij)(mn) = Smn — Sm,j — Si,n + Sij, (28)
что в среднем для каждого элемента изображения требует 2 операции для вычисления суммы Smn и 3 операции для вычислений по формуле (28). Однако для хранения сумм Smn требуется уже MxN ячеек, равное размеру изображения.
Алгоритм (28) также может быть модифицирован для трехмерного изображения. Пусть для каждой точки (m,n,k) подсчитаны суммы Smnk по прямоугольному параллелепипеду с диагональными элементами X000 и Xm-I,n-I,k-1, т.е. Smnk = S(ooo)(mnk). Нетрудно показать, что в таком случае сумма Sm (mnk) значений элементов внутри параллелепипеда с угловыми координатами (i,j, 7) и (m-1,n-1,k-1) вычисляется при помощи следующей операции:
S(ij7)(mnk) = Smnk — Smjk — Smnl — Sink + Smjl + Sijk + Sinl — Sijl. (29)
Таким образом с учетом того, что для вычисления каждого из значений Smnk требуется 3 арифметических операции, для вычисления суммы S(iji)(mnk) для каждого элемента трехмерного изображения потребуется в среднем 10 арифметических операций. Объем дополнительной требуемой памяти составит MxNxK ячеек с разрядностью достаточной для хранения значений сумм Smnk.
Аналогично можно находить дисперсии по фрагменту D^mnk), вычисляя значения сумм квадратов S(mnk)(x2) для каждой из точек (m,n,k) и для прямоугольного параллелепипеда S(j/)(mnk)(x2), и затем пользуясь формулой
D(jl)(mnk) = {S(j7)(mnk)(X2) - (S(j7)(mnk))2}/N(j7)(mnk), (30)
где S(iji)(mnk)(x2) — сумма квадратов значений элементов, попадающих в параллелепипед, а N^mnk) = (m-i)x(n-j)x(k-7) — число точек в параллелепипеде.
4.2. Порядковые статистики по прямоугольному параллелепипеду
Основой для вычисления порядковых статистик по фрагменту W как 2D- так и 3D-изображения служит гистограмма распределения значений яркости по этому фрагменту hW(x), а также ее интегральная характеристика Fw{x):
1=0 (31)
где Xmax — максимально возможное значение яркости, а NW — число точек во фрагменте W Порядковые статистики вида RW(n), где 0 < n < NW, представляют собой зависимость:
RW(n) = г, если FW(z-1) < n < FW(z). (32)
Алгоритм скользящего вычисления гистограммы по фрагменту строится аналогично формулам (26) и (27), т.е. при смещении фрагмента к следующей точке производится удаление точек на одной грани фрагмента и добавление точек на противоположной грани [15]. В двумерном случае алгоритм скользящего вычисления гистограммы по фрагменту при переходе от точки (m,n) к соседней точке (m,n+1) требует в среднем 2H числа операций (H — число строк во фрагменте). В трехмерном случае при переходе от точки (m,n,k) к точке (m,n,k+1) потребуется уже 2HxL операций, где HxL — число точек в грани параллелепипеда, перпендикулярной направлению смещения K.
В трехмерном случае число требуемых операций растет пропорционально произведению HxL. Однако, если (HxL) > (Xmax+1), то вместо операций со значениями отдельных точек выгодно заранее сформировать гистограммы для граней параллелепипеда hFjkXm,nM1)(x) и осуществлять уже операции вычитания и прибавления таких гистограмм:
h(i,j, k+1)( m,n,k+J+l)(X) = \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) арифметических операций на одну точку изображения независимо от размера параллелепипеда. Для пересчета гистограмм по граням параллелепипеда hFjkXmnM1)(x) нужно дополнительно в среднем 2L операций на точку и (Xmax+1)xK ячеек памяти для хранения K гистограмм.
Число операций (Xmax+1), требуемое для прибавления/вычитания каждой из гистограмм по граням параллелепипеда hj)(m,n,k+1)(x), можно значительно уменьшить, если воспользоваться свойством пространственной корреляции обрабатываемых данных. На участках с медленным изменением яркости, которых на реальных изображениях обычно большинство, размах значений элементов сравнительно невелик — в несколько раз меньше полного диапазона в (Xmax+1) градаций. Добавив к каждой из гистограмм h(jks(mnM1)(x) по 2 ячейки для запоминания минимального и максимального значений распределения, и, соответственно, обрабатывая лишь указываемый диапазон градаций, удается дополнительно в несколько раз сократить общее число операций.
Литература
1. Toriwaki J., Yoshida H. Fundamentals of Three-Dimensional Digital Image Processing. NY.: 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.
