Дискретизация. Антиалиасинг. Геометрические преобразования растровых изображений

Теорема Найквиста-Котельникова

Теорема Найквиста-Котельникова дает ответ на вопрос, какой частоты дискретизации fs достаточно для того, чтобы не произошло потери информации, т.е. чтобы по дискретизованному сигналу можно было восстановить исходный. Применительно к изображениям это грубо (поскольку еще не ясно, как происходит восстановление, которое зависит от устройства отображения) можно понимать так: "Какая разрешающая способность должна быть у растра, чтобы он сохранил все детали исходного аналогового изображения". Хотя потеря информации даже в случае соблюдения условий теоремы Котельникова произойдет из-за того, что значения дискретизованной функции (растрового изображения) в компьютере сами хранятся с ограниченной точностью. Передача цветов и оттенков лучшим образом при ограниченном диапазоне значений является задачей квантования, которая рассмотрена в лекции 12.

(рис 7.3) Срез изображения как сигнал и его частотный спектр.(рис 7.2) Гребенчатый фильтр и его преобразование Фурье.

В доказательстве теоремы и далее будет использоваться операция свертки функций I(x), J(x), определяемая так:

$$(I * J)(x) \stackrel{def}{ = }\int_{-\infty}^{\infty} I(u) \cdot J(x - u)du.$$

Теорема 7.1.1 (Найквиста-Котельникова). Для того чтобы сигнал I(x) можно было восстановить по его дискретному образу, его спектр должен быть ограничен максимальной частотой fH и частота дискретизации fs должна быть более 2fH.

Доказательство использует факты из математического и функционального анализа (см. например [3]). Пусть Is(x) - дискретный образ исходного сигнала I(x), как обычно, T = 1/fs - период дискретизации, тогда

$$I_s(x) = I(x)\cdot Comb(x) = I(x) \sum\limits_{n=-\infty}^{+\infty}\sigma (x-nT) = \\ \sum\limits_{n=-\infty}^{+\infty}I(nT)\cdot \sigma (x-nT).$$

Образом функции Comb в частотной области является функция

$$FComb(f) = f_s \sum\limits_{k=-\infty}^{+\infty}\sigma (f - kf_s),$$

а Фурье-образ I(x) по-прежнему будем обозначать F(f). Умножение функций в пространственной области соответствует их свертке (будем обозначать ее $$*$$ ) в частотной и наоборот. Соответственно, рассмотрим свертку F и FComb, являющуюся Фурье-образом Is(x) (обозначим его Fs(f) ):

$$F_s(f) = F * FComb(f) = F(f) * f_s \sum\limits_{k=-\infty}^{+\infty}\sigma (f - kf_s) \stackrel{(1)}{ = } \\ f_s \sum\limits_{k=-\infty}^{+\infty}F(f-kf_s),$$

где переход (1) произошел благодаря сдвигающему свойству дельта-функции при свертке. Как видно из последнего выражения, Fs(f) представляет собой бесконечную сумму функций F(f), умноженных на fs и сдвинутых на fs относительно друг друга, поэтому при условии fs > 2fH носители соседних сдвинутых версий не пересекаются, и отдельно, взяв центральную копию F(f) (k = 0) и применив к ней обратное преобразование Фурье, можно получить исходный сигнал I(x). Центральная копия берется путем умножения Fs(f) на прямоугольную функцию $$T \cdot Rect(T \cdot f)$$, где

$$Rect(f) =\left\{ \begin{array}{cc} 1, |f| \le 1 \\ 0, |f| > 1 \\ \end{array} \right.$$

Т.е. $$F(f) = F_s(f) \cdot T \cdot Rect(T \cdot f)$$, - образ исходной функции получен. Заметим, этому умножению в частотной области соответствует свертка в пространственной области. Применив обратное преобразование Фурье к $$T \cdot Rect(T \cdot f)$$, получим функцию

$$sinc(\pi x/T)\stackrel{def}{ = }\frac{\sin \pi x/T}{\pi x/T}.$$

Применив свертку с Is(x), получаем

$$I(x) = I_s(x) * sinc(\pi x/T) \\ = \sum\limits_{n=-\infty}^{+\infty}\left( I(nT) \cdot \sigma (x - nT)\right) * sinc((x - nT)/T) \\ \stackrel{(2)}{ = }\sum\limits_{n=-\infty}^{+\infty}I(nT) \cdot sinc((x - nT)/T),$$

где переход (2) также произошел благодаря сдвигающему свойству дельта-функции при свертке. Последняя формула называется Интерполяционной формулой Найквиста-Шеннона.

Для завершения доказательства осталось показать, что невозможно однозначно восстановить сигнал при $$f_s \le 2f_H$$.

Приведем соответствующий пример. Зафиксируем две частоты - fs и fH, $$f_s \le 2f_H$$ ; для упрощения рассуждений предположим, что fs > fH (в общем случае может быть более двух наложений сдвинутых образов, что усложнит построение контрпримера). Из формулы (7.6) следует, что Фурье-образ Is(x) является периодической функцией с периодом fs, поэтому вся информация для восстановления содержится в одном периоде (например $$f \in [-\frac{f_s}{2}, \frac{f_s}{2}]$$ ). Рассмотрим две функции, Фурье-образы которых равны

$$F_1(f) =\left\{ \begin{array}{ccc} f_H - |f|, |f| \in [0, f_s - f_H) \\ \frac{2f_H-f_s}{2}, |f| \in [f_s - f_H, f_H] \\ 0, |f| > f_H \\ \end{array} \right.$$ $$F_2(f) = \left\{ \begin{array}{cc} f_H - |f|, |f| \in [0, f_H] \\ 0, |f| > f_H \\ \end{array} \right.$$

(функция однозначно задается своим Фурье-образом). При дискретизации с частотой fs в соответствии с формулой (7.6) Фурье-образ в интервале $$[-\frac{f_s}{2}, \frac{f_s}{2}]$$ для обеих функций будет равен

$$F_s(f) = \left\{ \begin{array}{cc} f_s(f_H - |f|), |f| \in [0, f_s - f_H) \\ f_s(2f_H - f_s), |f| \in [f_s - f_H, f_s/2] \\ \end{array} \right.$$ (см. рис. 7.4).

Таким образом, в этом случае однозначная реконструкция невозможна.

(рис 7.5) Пример двух функций, дискретизированный образ которых совпадает(рис 7.4) Адекватная и неадекватная частоты дискретизации

Искажение сигнала и борьба с этим эффектом

Как было показано при доказательстве теоремы Найквиста-Котельникова, при недостаточной частоте дискретизации восстановленный сигнал будет искажен. Это наглядно можно наблюдать в частотной области (см. рис. 7.6), т.к. при этом копии частотного спектра исходного сигнала будут суммироваться в пересекающихся областях, что даст мнимое увеличение веса компонент с этими частотами в спектре, некую подмену высокочастотных компонент низкочастотными, эффект, известный как алиасинг (англ. aliasing). Особенно наглядно этот эффект можно наблюдать на резких контрастных изменениях яркости (см., например, изображение на рис. 7.7). Этот эффект легко объясним, если рассмотреть прямоугольный сигнал и его Фурье-образ ( sinc ), который имеет спектр бесконечной ширины (см. рис. 7.9, прямоугольный фильтр).

Для борьбы с подобными явлениями применяют префильтрацию - свертку с некой функцией фильтра (о фильтрации в общем случае см. соответствующую лекцию 8), что эквивалентно умножению на Фурье-образ функции-фильтра в частотной области, перед тем как производится дискретизация. Цель префильтрации - заранее отсечь высокочастотные компоненты, которые могут привести к алиасингу. В идеале, для этого следовало бы применять функцию sinc (см. формулу (7.7), рис. 7.8), действие которой как раз и состоит в отсечении высоких частот, но она имеет бесконечный носительНоситель функции - максимальная область, в которой она принимает значения отличные от 0., что затрудняет ее применение в пространственной области. На практике используются различные аппроксимации с ограниченным носителем (см. таблицу 7.1 и рис. 7.9 и 7.10). Гауссовский фильтр обычно применяется с $$\sigma = 0, 5$$ (так он лучше всего приближает sinc в частотной области), R берется порядка 2-3. У кубического фильтра присутствует параметр $$a \in [-3, 0]$$, с помощью которого можно регулировать степень размытия (увеличение a ) или, наоборот, более четкой передачи краев несмотря на некоторый алиасинг (уменьшение a ), стандарным компромиссным значением является a = -1. Для Lanzcos, который представляет собой обрезанный и сглаженный sinc, радиус R обычно равен 2 (так называемый Lanzcos2 ), реже - до 4. Большие значения не используются ввиду того, что носитель, а значит, и область интегрирования при свертке, будет слишком большой.

Понятно, что для дискретизации достаточно знать значения префильтрованной функции только в точках дискретизации. Это замечание служит основой для рассуждений в следующих разделах.

(рис 7.6) Эффект алиасинга при недостаточной частоте дискретизации.

В двумерном случае фильтрация описывается следующим образом: интенсивность пикселя I(x, y) с центром в точке (x, y) будет определяться формулой двумерной свертки

$$I(x, y) = \int\limits_{supp(F)}^{}F(s, t)Img(x - s, y - t) ds dt,$$

где F(s, t) - функция фильтра с центром в 0, supp(F) - ее носитель (область, где она не равна 0 ), Img(x, y) - непрерывное аналоговое изображение.

(рис 7.8) Видимый эффект алиасинга на изображении: слева - исходное изображение без антиалиасинга, справа - с антиалиасингом. (Изображения получены с помощью POV-Ray.)(рис 7.7) Функция-фильтр Sinc.

Одномерные функции-фильтры для антиалиасинга

Название	Функция фильтра F(x)
Импульсный (pulse)	$$F_p(x) = \left\{ \begin{array}{cc} 1, |x| \le 1/2 \\ 0, |x| > 1/2 \\ \end{array} \right.$$
Треугольный (triangle)	$$F_t(x) = \max \{1 - |x|, 0\}$$
Гауссовский (Gaussian)	$$F_{G(\sigma ,R)}(x) = \left\{ \begin{array}{cc} \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}}, |x| \le R \\ 0, |x| > R \\ \end{array} \right.$$
Кубический (cubic)	$$F_{c(a)}(x) = \left\{ \begin{array}{ccc} (a + 2)|x|^3 - (a + 3)|x|^2 + 1, 0 \le |x| < 1 \\ a|x| ^3 - 5a|x|^2 + 8a|x| - 4a, 1 \le |x| < 2 \\ 0, |x| > 2 \\ \end{array} \right.$$
Ланцоша (Lanzcos)	$$F_{L(R)}(x) = \left\{ \begin{array}{cc} sinc(x/R) \cdot sinc(x), 0 \le |x| \le R \\ 0, |x| > R \\ \end{array} \right.$$

Двумерные аналоги одномерного фильтра f11(x) строятся двумя путями:

как функция от радиуса: $$f_2(x, y) = f_1(r) = f_1(\sqrt{x^2 + y^2})$$ ;

как произведение: $$f_2(x, y) = f_1(x) \cdot f_1(y)$$.

Первый вариант - более корректный, но второй обладает свойством сепарабельности, т.е. в выражении (7.8) можно двумерное интегрирование разбить на два одномерных:

$$I(x, y) = \int\limits_{supp(F_1)\times supp(F_1)}^{}F_1(s)F_1(t)Img(x - s, y - t) ds dt \\ = \int\limits_{supp(F_1)}^{} ds F_1(s)\left( \int\limits_{supp(F_1)}^{}F_1(t)Img(x - s, y - t) dt \right)$$

Как будет показано в разделе 7.3, в случае передискретизации интегрирование заменяется суммированием и в этом случае быстрее производить вычисления в два этапа, каждый из которых состоит в дискретной свертке по одному измерению и может быть вычислен достаточно эффективно (в том числе и с применением параллельных алгоритмов). В алгоритмах растеризации с префильтрацией более популярны фильтры первого типа, а для геометрических преобразований, где происходит передискретизация - второго.

(рис 7.10) Фильтры для антиалиасинга.(рис 7.9) Фильтры для антиалиасинга (продолжение).

Двумерные аналоги фильтров, представленных в таблице 7.1 и на рис. 7.9 и 7.10, приведены в таблицах 7.2, 7.3 (для простоты частота дискретизации равна 1, произвольный случай получается масштабированием по осям). Прямоугольный фильтр соответствует либо цилиндру (первый тип), либо параллелепипеду с квадратным основанием (англ. Box filter) (второй тип); треугольный соответствует либо конусу, либо пирамиде соответственно. По смыслу прямоугольные и треугольные аналоги соответствуют так называемым невзвешенной и взвешенной площадным выборкам. Если рассмотреть Фурье-образы их одномерных аналогов на рис. 7.9, можно заметить, что взвешенная площадная выборка лучше отфильтровывает высокие частоты, поэтому рекомендуется применять именно ее.

Радиально-симметричные фильтры для антиалиасинга

Название	Функция фильтра $$F(s, t) = F_1(r)$$ ( $$r = \sqrt{s^2 + t^2}$$ )	Изображение
Цилиндрический	$$F_1(r) = F_p(r)$$
Конусообразный	$$F_1(r) = F_t(r)$$
Кубический	$$F_1(r) = F_{c(a)}(r)$$
Ланцоша	$$F_1(r) = F_{L(R)}(r)$$

Проблема растеризации с учетом антиалиасинга рассмотрена в разделе 7.2. В разделе 7.3 рассмотрены проблемы получения дискретных изображений по другим дискретным изображениям.

Одномерные функции-фильтры для антиалиасинга

Название	Функция фильтра $$F(s, t) = F_1(s) \cdot F_1(t)$$	Изображение
Параллелепипед	$$F_1(x) = F_p(x)$$
Пирамидальный	$$F_1(x) = F_t(x)$$
Кубический	$$F_1(x) = F_{c(a)}(x)$$
Ланцоша	$$F_1(x) = F_{L(R)}(x)$$

7.2. Антиалиасинг

Антиалиасинг или Фильтрация-сглаживание (англ. antialiasing) - устранение видимых артефактов дискретизации, особенно для изображений линий или краев объектов. Возможны два основных подхода.

Учет эффектов дискретизации в алгоритмах построения изображений. То есть, фактически, встраивание в них префильтрации.

Адаптивная дискретизация с постфильтрацией. Получение нового дискретизированного изображения по другому дискретизированному с помощью дискретной фильтрации, т.н. постфильтрации. Данный подход может быть эффективен при переходе от промежуточной большей частоты дискретизации к меньшей. Процесс растеризации с увеличенной по отношению к требуемой частотой первичной дискретизации с последующей постфильтрацией называется супердискретизацией (англ. supersampling). Первичная частота дискретизации может варьироваться в разных частях изображения в зависимости от его свойств.

Растеризация с антиалиасингом

При растеризации абстрактно определенных примитивов (точек, линий, многоугольников) можно рассматривать их как двумерный сигнал и применить к нему префильтрацию как описано в предыдущем разделе. В лекции 1 в качестве одной из моделей понятия растра рассматривалось замощение плоскости, где каждый пиксель представляет собой квадрат. При растеризации этому квадрату нужно приписать цвет так, чтобы при взгляде на изображение растра на каком-либо устройстве вывода воспринимаемая картинка была наиболее похожа на исходный непрерывный двумерный сигнал.

Заметим, что выражение (7.8) линейно по Img и, следовательно, обладает свойством суперпозиции по отношению к Img. Т.е. если Img состоит из нескольких непересекающихся объектов, то в (7.8) достаточно будет просуммировать вклад от этих объектов. Конечно, надо оговориться, что это возможно только в том случае, если атрибут пикселя принимает достаточное количество значений (для монохромных изображений антиалиасинг неприменим). Для цветных изображений проводимые рассуждения справедливы для каждого канала.

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

В данном разделе мы ограничимся лишь растеризацией базовых объектов - отрезков и общих кривых. Впрочем, подобные методы с незначительными модификациями применимы для антиалиасинга краев (которые также являются отрезками или общими кривыми) более сложных объектов (см., например, [29]).

(рис 7.11) Фильтрация для прямой.

Алгоритм Гупты-Спрулла

Данный алгоритм для растеризации отрезков с целочисленными координатами концов был предложен Гуптой и Спруллом в статье [34]. В нем предполагается использование радиально-симметричного фильтра (такого как конический). Пусть его радиус равен R. Для того чтобы выяснить, каким цветом следует закрасить пиксель, необходимо посчитать значение свертки фильтра с функцией прямой ( I = 1 в вытянутом прямоугольнике толщиной t и 0 вне него), см. рис. 7.11. Обозначим результат этой операции FR,t(r), где r - расстояние от центра закрашиваемого пикселя до центральной оси прямой. Заметим, что FR,t(r) = 0 для r > R + t/2.

Идея алгоритма состоит в модификации алгоритма Брезенхема (см. раздел 3.2) для одновременной закраски нескольких пикселей разной интенсивности в столбце (см. раздел 3.2).

(рис 7.12) Алгоритм Гупты-Спрулла.

Количество пикселей, которые могут быть закрашены в одном столбце зависит от t и R, а также от наклона прямой. Канонически будем полагать радиус фильтра R = 1 и толщину отрезка t = 1. Тогда максимально может быть закрашено 4 пикселя при наклоне, близком к $$45^{\circ}$$. Если рассматривать конический фильтр, то в случае 4 пикселей возможный вклад самого верхнего или нижнего будет достаточно мал, поэтому для упрощения алгоритма будем закрашивать только ближайшие 3 к пересечению прямой с вертикалью пикселя.

Пусть, как и в разделе 3.2, строится отрезок $$(0,0) \to (a,b)$$.

Обозначим нижний текущий (как в алгоритме Брезенхема) пиксель P2, нижний - P1, а верхний - P3 (см. рис. 7.12). Пусть v - расстояние от центральной точки между рассматриваемыми пикселями до выбранного пикселя P2. Оно может быть как положительным (если выбран шаг s ), так и отрицательным (если выбран шаг d ). Тогда соответствующие расстояния до оси прямой равны

$$r_1 = (1 + v) \cos \theta, r_2 = |v| \cos \theta, r_3 = (1 - v) \cos \theta,$$

где $$\theta$$ равен углу наклона прямой (см. рис. 7.12), т.е. $$\cos \theta = a/ \sqrt{a^2 + b^2}$$.

Если обратиться к исходному алгоритму Брезенхема (см. раздел 3.2), несложно заметить, что $$$e = (v - {\frac{1}{2}}) \cdot 2a$$$ или $$2av = e+a$$. Тогда $$v \cos \theta = \frac{e+a}{2\sqrt{a^2+b^2}}$$.

Для ускорения вычислений построим дискретную аппроксимацию FR,t(r), разбив интервал [0,R + t/2] на равные части и аппроксимируя значение F в каждом интервале значением в центре этого интервала. Получим таблицу значений интенсивности ITable.

Пусть функция round находит номер интервала по действительному значению в [0,R + t/2]. Тогда, модифицируя алгоритм Брезенхема, получаем следующий алгоритм:

// plot(x,y,I) закрашивает пиксель (x,y) с интенсивностью I

e = 2b - a;
Delta eS = 2b;
Delta eD = 2b - 2a;

denom = 1/(2*sqrt(a*a+b*b));
fixPart = 2a*denom; // фикс. часть r_1 и r_3

// (x,y) - Координаты текущей точки
x= 0; y = 0;

while( x < a )
{
      av2 = (e+a)*denom;

      plot(x, y-1, ITable[round(fixPart+av2)] );
      plot(x, y, ITable[round(abs(av2))] );
      plot(x, y+1, ITable[round(fixPart-av2)] );

      if( e > 0 )
      {
            // d : диагональное смещение
            x++; y++;
            e += Delta eD;
      }
      else
      {
            // s : горизонтальное смещение
            x++;
            e += Delta eS;
      }
}

Конечно, еще желательно специально обрабатывать начало и конец отрезка, но подробное обсуждение этого выходит за рамки данной книги. Подробности можно узнать в [34].

Алгоритм Ву

Ву разработал альтернативный алгоритм для растеризации кривых со встроенным антиалиасингом [54]. В нем он предложил несколько другой критерий для выбора значений интенсивности пикселей.

Пусть мы хотим растеризовать некоторую кривую. Без ограничения общности будем считать, что она локально наклонена горизонтально, т.е. ее можно представить как функцию y = f(x) и f'(x) < 1. При классической растеризации без антиалиасинга для каждого столбца пикселей с абсциссой i просто закрашивается ближайший к f(i) пиксель с ординатой $$Y_i = \lfloor f(i) + 1/2 \rfloor $$, т.е. минимизируется ошибка |Yi - f(i)|. В то же время Ву замечает, что при такой растеризации возможны неприятные визуальные эффекты (см. рис. 7.13), и предлагает минимизировать динамическую ошибку, возникающую при растеризации Eij = |(f(i) - f(j)) - (Yi - Yj)| и показывающую насколько искажается наклон кривой. Если часть кривой соответствует абсциссам в отрезке [0, a], то предлагается минимизировать норму матрицы $$E_{ij}, i, j \in \overline{0, a}$$, что в общем случае является вычислительно сложной задачей комбинаторной оптимизации.

Если же рассматривать возможность закрашивать пиксели с несколькими уровнями интенсивности, то можно в каждом столбце распределять исходную интенсивность в точке (i, f(i)) по двум соседним пикселям так, чтобы центр тяжести (если интенсивность рассматривать как массу) приходился на точку (i, f(i)) ; таким образом с учетом усреднения получаем E = 0. Пусть исходная интенсивность кривой равна I0, тогда

$$I(i, \lfloor f(i) \rfloor ) = I_{-}, I(i, \lceil f(i) \rceil ) = I_{+},$$

так что

$$I_0 = I_{-} + I_{+} \\ I_0f(i) = I_{-}\lfloor f(i) \rfloor + I_{+}\lceil f(i) \rceil .$$ (рис 7.13) Погрешность при аппроксимации кривой.

Из этой системы легко получается решение:

$$I_{-} = I_0(\lceil f(i) \rceil - f_i),$$ $$I_{+} = I_0(f_i - \lfloor f(i) \rfloor ).$$

Если рассмотреть данный подход с точки зрения фильтрации, то сигнал здесь представлен объединением квадратов

$$|x - i| \le \frac{1}{2}, |y - f(i)| \le \frac{1}{2}, i \in \overline{0, a},$$

и используется фильтр-параллелепипед радиуса 1/2 (см. рис. 7.14).

Рассмотрим применение данного подхода к растеризации отрезка с целочисленными координатами концов. Пусть строится отрезок $$(0,0) \to (a,b)$$. Тогда прямая определяется формулой y = kx, где k = b/a. Для максимального быстродействия будем проводить все вычисления приближенно, используя лишь целочисленную арифметику. Пусть D отвечает за дробную часть y. $$D = \{f(i)\} \cdot 2^n$$, где n - число разрядов в машинном представлении целых чисел. На каждом шаге будем прибавлять к D аппроксимацию k (обозначим ее d ):

$$D = D + d, d = \lfloor k2^n + 1/2\rfloor .$$ (рис 7.14) Аппроксимация кривой по Ву.

Ошибка аппроксимации e = k - d2-n не превосходит по модулю 2-n. D естественным образом будет содержать дробную часть f(i), умноженную на 2n, т.к. машинная арифметика автоматически осуществляет сложение по модулю 2n. Единственное, что нужно в тех случаях, когда происходит переполнениеОбычно процессор содержит специальный флаг, который указывает на то, произошло ли переполнение при предыдущей арифметической операции., увеличивать текущее значение y (которое отвечает за целую часть kx ) на 1 (т.е. диагональный сдвиг).

Пусть значения интенсивности пикселя определяются m -разрядным числом, тогда максимальное значение I0 = 2m-1. m, как правило, много меньше n (обычно, m = 8, n = 32 ). В соответствии с уравнениями 7.9 для столбца с абсциссой x,

$$I_{+} = I(x, y + 1) = I_0(kx - \lfloor kx\rfloor ) \\ = (2^m - 1)(D2^{-n} + ex) \\ = D2^{m-n} + (2^m - 1)ex - D2^{-n}.$$

Учитывая, что I+ должна быть целочисленной, пренебрежем последними двумя членами (в статье [54] подробно обосновывается почему ошибка будет пренебрежимо малой): I+ = D2n-m. Фактически это целое, полученное из m старших двоичных разрядов D. Отсюда довольно легко получить и $$I_{-} = I_0 - I_{+} = \overline{I_{+}}$$ ( $$\overline{\phantom{I_{+}}}$$ здесь означает двоичное дополнение). Это следует из того, что битовое представление 2m - 1 - D2m-n является двоичным дополнением D2m-n.

Также используя оптимизацию, заключающуюся в одновременной закраске с двух концов (см. раздел 3.2), получаем следующий алгоритм.

// Координаты концов отрезка - (0,0) и (a,b)
// plot(x,y,I) закрашивает пиксель (x,y) с интенсивностью I
// I0 - максимальная интенсивность (2^m-1)

x0 = 0; x1 = a; y0 = 0; y1 = b;
plot(x0,y0,I0);
plot(x1,y1,I0);

D = 0;
d = floor( (b/a)*2^n + 0.5 );

while( x0 < x1 )
{
      D = D + d;
      if( произошло переполнение D )
      {
            y0++; y1--;
      }

      I1 = D / 2^(n-m); // битовый сдвиг вправо на n-m
      I2 = двоичное_доп( I1 );
      plot(x0,y0,I1);
      plot(x1,y1,I1);
      plot(x0,y0+1,I2);
      plot(x1,y1-1,I2);
}

Хотя алгоритм Ву теоретически не очень качественный, т.к. использует грубую аппроксимацию кривой и грубый фильтр-параллелепипед, но на практике он успешно справляется с задачей антиалиасинга, обладая при этом большей простотой реализации (в частности, аппаратно) и скоростью, чем алгоритм Гупты-Спрулла.

7.3. Геометрические преобразования растровых изображений

Построение дискретного изображения, которое было подвергнуто геометрическому преобразованию, при условии наличия дискретизации исходного изображения является фактически проблемой передискретизации.

Пусть исходное (непрерывное) изображение задано функцией I(x, y), тогда его дискретизация Is(x, y) (считаем частоту равной fs ) будет получена умножением на функцию двумерной "решетки":

$$I_s(x, y) = I(x, y) \cdot Comb_{f_s}(x, y) \\ = I(x, y) \cdot \sum\limits_{i,j}^{}\sigma \left( x - \frac{i}{f_s}, y - \frac{j}{f_s}\right), (i, j) \in \mathbb{Z}^2 \\ = \sum\limits_{i,j}^{}I \left( \frac{i}{f_s}, \frac{j}{f_s}\right) \cdot \sigma \left( x - \frac{i}{f_s}, y - \frac{j}{f_s}\right) .$$

Следовательно, мы получили растр, с атрибутами

$$I_s(i, j)\stackrel{def}{ = }I \left( \frac{i}{f_s}, \frac{j}{f_s}\right) .$$

Пусть преобразование задано функцией $$T \colon \mathbb{R}^2 \to \mathbb{R}^2$$:

$$\left( \begin{array}{c} x' \\ y' \end{array} \right) = T \left( \begin{array}{c} x \\ y \end{array} \right) .$$

Тогда преобразованное изображение равно I'(x', y') = I(T-1(x', y')), дискретизация (c частотой f's ) преобразованного изображения должна быть получена как

$$I' (x', y') = I'(x', y') \cdot Comb_{f_s'}(x', y') \\ = I' (x', y') \cdot \sum\limits_{i',j'}^{}\sigma \left( x' - \frac{i'}{f_s'}, y' - \frac{j'}{f_s'}\right), (i', j') \in \mathbb{Z}^2 \\ = \sum\limits_{i',j'}^{}I' \left( \frac{i'}{f_s'}, \frac{j'}{f_s'}\right) \cdot \sigma \left( x' - \frac{i'}{f_s'}, y' - \frac{j'}{f_s'}\right) \\ = \sum\limits_{i',j'}^{}I \left( T^{-1}\left( \frac{i'}{f_s'}, \frac{j'}{f_s'}\right)\right) \cdot \sigma \left( x' - \frac{i'}{f_s'}, y' - \frac{j'}{f_s'}\right).$$

Мы рассматриваем задачу, когда нам неизвестно исходное изображение, а известна только его дискретизация Is. В таком случае в качестве I в (7.12) следует подставить реконструированное изображение Ir. Реконструкция производится с помощью реконструирующего фильтра RFilter:

$$I_r(x, y) = \int_{}^{} I_s(s, t) \cdot RFilter(x - s, y - t) ds dt$$ $$\stackrel{7.11}{ = } \sum_{i,j}^{} I_s(i, j) \cdot RFilter \left( x - \frac{i}{f_s}, y - \frac{j}{f_s}\right) .$$

В идеальном случае, когда начальная частота дискретизации была больше частоты Найквиста и в качестве RFilter выступает функция sinc изображение будет реконструировано точно, т.е. Ir = I. На практике же в качестве RFilter выступают различные аппроксимации sinc с локальным носителем (см. таблицу 7.1), что дает некоторые искажения (в частности, размытие).

Если в выражение (7.12) подставить Ir, вычисленное по формуле (7.14) вместо I и воспользоваться определением (7.11) для атрибутов нового растра I's(i', j') то получаем, что

$$I' (x', y') = \sum\limits_{i,j}^{} I_s(i, j) \cdot RFilter \left( T ^{-1} \left( \left( \begin{array}{c} \frac{i'}{f_s'}, \\ \frac{j'}{f_s'} \end{array} \right) \right) - \left( \begin{array}{c} \frac{i}{f_s}, \\ \frac{j}{f_s} \end{array} \right) \right) .$$

Таким образом, задача передискретизации сводится к применению дискретной свертки (фактически суммированию) исходного дискретизованного изображения Is с функцией фильтра RFilter, центрированной в прообразе нового пикселя при преобразовании T (см. рис. 7.15).

(рис 7.15) Передискретизация.

Линейные фильтры, используемые в качестве RFilter, фактически осуществляют локальную интерполяцию или аппроксимацию поверхности Ir по точкам дискретизации. Если мы просто будем считать Ir некоторой поверхностью, то в этом случае

$$I_s'(i', j') = I_r \left( T^{-1}\left( \frac{i'}{f_s'}, \frac{j'}{f_s'}\right) \right) .$$

Если в качестве фильтра выступает функция-параллелепипед, то при передискретизации новое значение будет просто средним по пикселям, попавшим в область носителя (для малого радиуса ( 1/2 ) - просто значение ближайшего пикселя). Пирамидальному фильтру соответствует билинейная интерполяция. Более качественные результаты можно получить при использовании бикубической интерполяции, которая соответствует сепарабельному кубическому фильтру (см. таблицу 7.3). Она является стандартной в популярном растровом редакторе Adobe Photoshop.

(рис 7.16) Аффинные преобразования.

Предметом нашего рассмотрения будут в основном аффинные преобразования:

$$\left( \begin{array}{c} x' \\ y' \end{array} \right) = \left( \begin{array}{cc} a_{11} a_{12} \\ a_{21} a_{22} \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) + \left( \begin{array}{c} b_1 \\ b_2 \end{array} \right) ,$$

частным случаем которых являются сдвиги, растяжения, скосы и повороты (см. рис. 7.16).

В трехмерной графике подобные проблемы возникают при текстурировании, т.е. наложении искаженных растровых изображений на поверхности объектов, которые затем проецируются на экран. Там, как правило, речь идет о перспективных преобразованиях (задаются невырожденной трехмерной матрицей $$= \{a_{ij}\}_{i=1,j=1}^{3,3}$$ ):

$$x' = \frac{a_{11}x + a_{12}y + a_{13}}{a_{31}x + a_{32}y + a_{33}},$$ $$y' = \frac{a_{21}x + a_{22}y + a_{23}}{a_{31}x + a_{32}y + a_{33}}.$$ (рис 7.17) Супердискретизация

При перспективных преобразованиях возможны значительные искажения, при которых эффективная частота дискретизации, обратно пропорциональная расстоянию между переведенными точками $$T_{-1}\left( \frac{i'}{f_s'}, \frac{j'}{f_s'}\right)$$, может быть значительно больше исходной, которая полагается равной 1 (в пикселях исходного изображения). В этом случае прибегают к супердискретизации. Локально для одного пикселя в I' строится промежуточная дискретизация с частотой, большей чем f' (обычно в целое число раз, которое зависит от степени искаженияВопрос о том, как приблизительно вычислить это число, выходит за рамки данной книги. Довольно хорошо разработан данный вопрос для задач текстурирования. ), а затем по ней с помощью второго этапа дискретной фильтрации (зачастую здесь применяется простейший усредняющий фильтр-параллелепипед) получаем значение интенсивности результирующего пикселя (см. рис. 7.17). Более подробно о супердискретизации можно узнать из литературы, посвященной текстурированию.

Как всегда, на практике, пожертвовав точностью, можно получить более быстродействующие алгоритмы.

Подход Веймана

Вейманом был предложен подход [50], основанный на кодах Ротштейна (см. разделе 3.2, это эквивалентно кодированию смещений при растеризации отрезка. Данный подход позволяет осуществлять такие преобразования, как скос и масштабирование с коэффициентами, заданными рациональными числами. С помощью комбинации этих преобразований можно представить произвольное аффинное преобразование (см. ниже).

(рис 7.18) Код Ротштейна.

Рассмотрим задачу масштабирования по горизонтали из изображения шириной q пикселей в изображение шириной p пикселей, p > q (т.е. коэффициент масштабирования равен p/q ). Для начала строится код Ротштейна (например, алгоритмом Брезенхема (см. раздел 3.2)). Наиболее просто скопировать q исходных колонок в те колонки, которые помечены кодом 1. Но это приведет к пробелам в тех столбцах, где значение кода равно 0. Поэтому данное распределение по колонкам применяется к циклическим перестановкам кода Ротштейна; при этом происходит суммирование в каждой колонке, а затем результат усредняется (делится на q, т.к. всего было p циклических перестановок и каждая колонка была закрашена q раз).

Для случая сжатия (когда p < q ) строится код длины q, и теперь уже, наоборот, 1 -цы соответствуют выбору тех исходных столбцов, которые копируются в результирующие p. Аналогично, с целью устранения алиасинга, результат усредняется по всем циклическим перестановкам кода (деление снова на q, т.к. прибавка к каждой колонке происходит при каждой перестановке).

// q - исходная ширина, p - ширина результата (p > q)
// I0 - исходное изображение (q x n)
// I1 - результирующее изображение (p x n)

I1 = 0; // установим все пиксели в I1 равными 0
code = ComputeCode( p, q ); // код длины p

// по всем циклическим перестановкам
foreach( perm in 1...p )
{
      code = CyclicShift( code );
      srcI = 0;

      foreach( dstI in 1...p )
            if( code[dstI] == 1 )
            {
                  foreach( J in 1...p )
                        I1(dstI,J) += I0(srcI,J);
                  srcI++;
            }
}
// усредняем
foreach( pixel in I1 )
      I1(pixel) = I1(pixel) / q;

Для вертикального скоса с коэффициентом p/q < 1 применяется модификация, где 1 в коде указывает на то, что текущий столбец следует сдвинуть вверх на 1 пиксель (см. рис. 7.19). Так же, как и для растяжения, с целью фильтрации результат усредняется по циклическим перестановкам кода.

(рис 7.19) Алгоритм Веймана для вертикального скоса.

Разложение преобразований в композицию более простых

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

Самый простой и полезный пример - представление общего преобразования T в виде композиции преобразований по столбцам и по строкам. Преобразования, сохраняющие столбцы или строки, можно эффективно распараллелить (по столбцам и строкам соответственно).

Пусть $$T = R \circ C$$, где C сохраняет столбцы, а R - строки. Пусть $$T \colon (x, y) \mapsto (x'', y'')$$, а промежуточный результат - $$(x', y') (C \colon (x, y) \mapsto (x', y'), R \colon (x', y') \mapsto (x'', y''))$$. Тогда

$$T = \left( \begin{array}{c} t_1(x, y) \\ t_2(x, y) \end{array} \right); C = \left( \begin{array}{c} x \\ c_2(x, y) \end{array} \right); R = \left( \begin{array}{c} r_1(x', y') \\ y' \end{array} \right) .$$

Отсюда получаем, что

$$T = \left( \begin{array}{c} t_1(x, y) \\ t_2(x, y) \end{array} \right) = \left( \begin{array}{c} r_1(x, c_2(x, y)) \\ c_2(x, y) \end{array} \right) .$$

Таким образом, c2(x, y) = t2(x, y) и r1(x, t2(x, y)) = t1(x, y). Для того, чтобы найти r1, надо в выражении t1 вычленить все подвыражения, содержащие y, и привести их к виду t2(x, y). После этого, заменив эти подвыражения на y, получим выражение для r1.

Рассмотрим эту процедуру на примере поворота на угол $$\theta$$:

$$T(x, y) = \left( \begin{array}{cc} \cos \theta -\sin \theta \\ \sin \theta \cos \theta \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) .$$

Тогда

$$t_1(x, y) = \cos \theta \cdot x - \sin \theta \cdot y, \\ c_2(x, y) = t_2(x, y) = \sin \theta \cdot x + \cos \theta \cdot y,$$

тогда $$y = (t_2(x, y) - \sin \theta \cdot x)/ \cos \theta$$ (это, конечно, возможно, когда $$\theta \ne \pm 90^{\circ})$$. Отсюда

$$t_1(x, y) = \cos \theta \cdot x - \sin \theta \cdot \left(\frac{t_2(x, y) - \sin \theta \cdot x)}{\cos \theta} \right) \\ = \sec \theta \cdot x - \tan \theta \cdot t_2(x, y).$$

Получаем, что $$r_1(x, y) = \sec \theta \cdot x - \tan \theta \cdot y$$. Итак, получено следующее разложение (наглядно на рис. 7.20):

$$\stackrel{T}{\left( \begin{array}{cc} \cos \theta -\sin \theta \\ \sin \theta \cos \theta \end{array} \right)} = \stackrel{R}{ \left( \begin{array}{cc} \sec \theta - \tan \theta \\ 0 1 \end{array} \right)} \cdot \stackrel{C}{ \left( \begin{array}{cc} 1 0 \\ \sin \theta \cos \theta \end{array} \right)}$$

В вырожденном случае ( $$\theta = \pm 90^{\circ}$$ ) такое разложение не получится, но в этом случае всё T сведется к замене x на y и возможному отражению (когда $$\theta = -90^{\circ}$$ ). Следует также отметить, что вращений на углы $$\theta$$, близких к $$\pm 90^{\circ}$$, следует избегать, т.к. при их применении произойдут слишком большие искажения по вертикали (слишком сильное сжатие). В этом случае корректнее произвести поворот сначала на $$\pm 90^{\circ}$$, а затем уже на $$\theta - (\pm 90^{\circ})$$ (знак соответствует близкому к $$\theta$$ углу). Преобразования C и R представляют собой композицию сжатия и скоса по вертикали и горизонтали соответственно. Поэтому данная декомпозиция позволяет реализовать поворот с помощью нескольких последовательных применений алгоритма Веймана для скосов и масштабирования.

(рис 7.20) Разложение вращения на скос и смещение по вертикали (C), а затем по горизонтали (R).

Существует альтернативное разложение на 3 скоса для матрицы поворота (наглядно на рис. 7.21), которое позволяет избавиться от существенных искажений при $$\theta \approx \pm 90^{\circ}$$, присущих вышеизложенному алгоритму:

$$\left( \begin{array}{cc} \cos \theta -\sin \theta \\ \sin \theta \cos \theta \end{array} \right) = \left( \begin{array}{cc} 1 - \tan {\frac{\theta}{2}} \\ 0 1 \end{array} \right) \left( \begin{array}{cc} 1 0 \\ \sin \theta 1 \end{array} \right) \left( \begin{array}{cc} 1 - \tan {\frac{\theta}{2}} \\ 0 1 \end{array} \right)$$ (рис 7.21) Разложение вращения на 3 скоса.

Алгоритм Веймана рекомендуется применять именно вместе с таким разложением.

7.4. Заключение

Тема дискретизации, представления и обработки изображений (англ. image processing) является самостоятельной научной дисциплиной. Более подробно с ней можно ознакомиться по книгам [44], [32], [9]. Тема геометрических преобразований подробно раскрыта в книге [52].