Спектральные методы обработки изображений

Цифровая обработка изображений не ограничивается задачами трансформации, которые были рассмотрены в предыдущем разделе. Важной составляющей является кодирование информации, то есть представление изображения в виде, который позволяет эффективное хранение и передачу графической информации, выделить те или иные характеристики, которые нужны при анализе объектов на данном изображении, определить наиболее информативные характеристики и прочее. Для этой цели используются частотные методы анализа, прежде всего это анализ Фурье и анализ главных компонент. Анализ Фурье может использовать тригонометрический базис, или базис с компактным носителем, что приведет к вейвлет-методам.

Прямое дискретное преобразование Фурье (ДПФ) и обратное ему преобразования, имеют стандартную форму \[ F(\nu)=N^{-1}\sum_{\tau=0}^{N-1}f(\tau)\exp\left(-i\frac{2\pi\nu \tau}{N}\right), f(\tau)=\sum_{\nu=0}^{N-1}F(\nu)\exp\left(i\frac{2\pi\nu \tau}{N}\right). \] Часто сигнал (изображение в том числе) нужно преобразовать только в действительной области, без использования комплексных переменных. В этом случае традиционным методом, который разделяет частотные характеристики сигнала является дискретное косинус преобразования (ДКП), из которого и начнем.
Рассмотрим одномерное ДКП для N = 8. Выбор такого значения N не случаен, в формате Jpeg изображение разбивается на квадраты 8 × 8, что и привело к тому, что говоря об использовании частотных методов в обработке изображений используют как раз такое число составляющих дискретного преобразования.
Выпишем прямое дискретное косинус преобразование \[ F_\nu=\frac{1}{\sqrt{2}}C_\nu\sum_{i=0}^7p_i\cos\frac{(2i+1)\nu\pi}{16},\\ C_i=\left\{ \begin{array}{ll} \frac{1}{\sqrt{2}}, &\hbox{i=0;} \\ 1, &\hbox{ i>0}. \end{array} \right. \] Обратное дискретное косинус преобразование \[ p_i=\frac{1}{\sqrt{2}}C_i\sum_{\nu=0}^7F_\nu\cos\frac{(2i+1)\nu\pi}{16}. \] Понятное дело, что кроме ДКТ существуют другие тригонометрические преобразования, например ДСП - дискретное синус преобразования или ДПХ - преобразование Хартли.
Этот вид преобразования назван в честь Р.Хартли, который в 1942 г. опубликовал статью о паре интегральных преобразований - прямое и обратное, использующих введенную им функцию \(\mathrm {cas} \theta = \sin \theta + \cos \theta \). К началу 1980-х годов эти результаты оставались в забвении, пока к ним не привлек внимание исследователей Р.Брейсуэлл, разработавший основы теории как непрерывного, так и дискретного преобразования Хартли.
Дискретное преобразование Хартли (ДПХ) действительной функции f(τ) и соответствующее обратное преобразование определяются соотношениями \[ H(\nu)=\frac{1}{N}\sum_{\tau=0}^{N-1}f(\tau)\mathrm{cas } \frac{2\pi\nu\tau}{N}, \] \[ f(\tau)=\sum_{\nu=0}^{N-1}H(\nu)\mathrm{cas } \frac{2\pi\nu\tau}{N}, \] где используется обозначение cas θ= cos θ+sin θ, введенное Хартли.
Для доказательства данного факта воспользуемся свойством ортогональности \[ \sum_{\nu=0}^{N-1}\mathrm{cas } \frac{2\pi\nu\tau}{N}\mathrm{cas } \frac{2\pi\nu t}{N} =\left\{ \begin{array}{ll} N, &\tau=t;\\ 0, &\tau\ne t. \end{array} \right. \] Подставляя величину \[ \frac{1}{N}\sum_{\tau=0}^{N-1}f(\tau)\mathrm{cas } \frac{2\pi\nu\tau}{N}, \] которая определяет преобразование H (v) , в выражение \[ \sum_{\nu=0}^{N-1}H(\nu)\mathrm{cas } \frac{2\pi\nu\tau}{N}, \] получаем \[ \sum_{\nu=0}^{N-1}H(\nu)\mathrm{cas } \frac{2\pi\nu\tau}{N}= \sum_{\nu=0}^{N-1}\frac{1}{N}\sum_{t=0}^{N-1}f(t)\mathrm{cas } \frac{2\pi\nu t}{N}\mathrm{cas } \frac{2\pi\nu\tau}{N}= \frac{1}{N}\sum_{t=0}^{N-1}f(t)\sum_{\nu=0}^{N-1}\mathrm{cas } \frac{2\pi\nu t}{N}\mathrm{cas } \frac{2\pi\nu\tau}{N}= \frac{1}{N}\sum_{t=0}^{N-1}f(t)\times \left\{ \begin{array}{ll} N, &\tau=t;\\ 0, &\tau\ne t \end{array} \right. =f(\tau). \] что подтверждает справедливость обратного преобразования.
Коэффициент \(N^{-1} \) в ДПХ заимствуется из практики использования ДПФ, для которого величина F (0) равна постоянной составляющей функции \(f (\tau) \), то есть, ДПХ является симметричной процедурой. Кроме этого, ДПХ является действительным преобразованием, поскольку действительной есть функция \(f (\tau) \).
Проведем сравнение приведенных дискретных преобразований на примере множества точек p = {12,10,8,10,12,10,8,11}. Найдем ДКТ, ДСП и ДПХ этого множества. Проведем следующий эксперимент - отквантуем полученные результаты с точностью до десятков, а для синус преобразования, даже вдвое точнее.

dct={30,0,0,0,0,0,0,0}- черний цвет,
dst={0,25,0,10,0,10,0,10}- красный цвет,
dht={80,0,10,0,0,0,10,0}- зеленый цвет.

Применяя обратные преобразования, восстановим данные и найдем ошибку, которая визуализированна ниже Как видно даже из приведенного простого примера, худшее качество показывает синус-преобразование, а лучшее - преобразование Хартли. Но не все так просто, потому что количество ненулевых коэффициентов ДКП меньше чем у преобразования Хартли.
Полученные результаты говорят о перспективности преобразования Хартли. Для того, чтобы проверить и убедиться в эффективности дискретного преобразования Хартли перед косинус - и синус-преобразованием, построим нормы частотных доменов, полученных с использованием соответствующих преобразований.
Используя идеологию JPEG, разобьем изображение на квадраты со стороной в восемь пикселов и на каждом квадрате вычислим соответствующие частотные коэффициенты. Для ДКП это \[ c^k_{i,j}=\frac{1}{4}C_iC_j\sum_{x=0}^7\sum_{y=0}^7p_{x,y}\cos\frac{(2x+1)i\pi}{16}\cos\frac{(2y+1)j\pi}{16}, \] где \[ C_i=\left\{ \begin{array}{ll} \frac{1}{\sqrt{2}}, &i=0;\\ 1, &i\gt 0. \end{array} \right. \] для ДСП \[ s^k_{i,j}=\frac{1}{4}C_iC_j\sum_{x=0}^7\sum_{y=0}^7p_{x,y}\sin\frac{(2x+1)i\pi}{16}\sin\frac{(2y+1)j\pi}{16}, \] и для ДПХ \[ h^k_{i,j}=\sum_{x=0}^7\sum_{y=0}^7p_{x,y}\mathrm{cas }\frac{2xi\pi}{16}\mathrm{cas }\frac{2yj\pi}{16}, \] где k - номер квадрата.
Множество соответствующих частотных коэффициентов \(D (q_{i, j}) \left \{q^k_{i, j}| k = 0,1, ... \right \}\) будем называть частотным доменом. Норма частотного домена находится по правилу \(\| D (q_{i, j}) \| = \sqrt{\sum_k \left (q^k_{i, j} \right)^2} \). По сути, норма, это энергия, соответствующая коэффициентам соответствующего домена. Чем больше значение нормы, тем больше информации аккумулирует этот домен. Поэтому, об эффективности того или другого преобразования можно судить по тому, насколько быстро убывают нормы частотных доменов. Для тестового изображения Lena имеем следующие результаты

Наконец, преобразование Хартли иллюстрируют диаграммы

Чтобы сравнить эффективность кодирования ДКП и ДПХ (синус-преобразование сразу отбросим как неэффективное, этот факт сразу вытекает из анализа гистограмм энергии доменов), возьмем метод сжатия изображений Jpeg и заменим в нем ДКП на ДПХ.
Коротко остановимся на методе сжатия изображений JPEG. Схематично метод JPEG заключается в следующем:

1. Цветное изображение преобразуется из пространства смешивания цветов RGB в пространство YCrCb, где Y- люминесцентная составляющая (освещенность), Cr- компонента теплых оттенков и Cb- компонента холодных оттенков.
2. На втором этапе происходит прореживания пикселей (при высоких степенях сжатия из четырех значений Cr и Cb берутся не все, два или даже одно значение).
3. Группирование значений Y, Cr и Cb в блоки 8 × 8 с последующей фильтрацией, путем применения дискретного косинус-преобразования (ДКП).
4. Квантование полученных коэффициентов, путем деления каждого коэффициента косинус-преобразования Фурье на специальное число (QC коэффициент квантования), с последующим округлением до целого. Именно на этом этапе идет невосполнимая потеря информации, за счет чего, собственно, и реализуется существенное сжатие изображения.
5. Далее отквантованные значения упорядочиваются (согласно змейке JPEG) и сжимаются методом RLE (Run-lengthEncoding) - сжатие с использованием счетчиков нулей.
6. И на последнем этапе применяется префиксный или энтропийный метод сжатия данных без потерь.

В качестве критерия оценки восстановления исходного изображения \(I (i, j) \) и изображения \(K (i, j), (i = 0,1, ..., m-1, j = 0,1 ,. .., n-1) \), используем отношение сигнал / шум \[ PSNR(I,K)=20 \log_{10}\frac{255}{\sqrt{MSE(I,K)}} \] где \[ MSE(I,K)=\frac{1}{n\times m}\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}|I(i,j)-K(i,j)|^2.\] Для тестирования используем изображения тестовой базы TID2008 фирмы Kodak, к которым применим сжатие с двумя разными наборами QC - коэффициентами квантования. Результаты экспериментов объединим на одном графике.
Ось абсцисс соответствует размеру файла, ось ординат - значение PSNR. Точки синего цвета соответствуют ДКП, пурпурного - ДПХ. Видно, что они образуют по два кластеры, которые соответствуют режимам сжатия. Кроме того, приведены две линейные аппроксимации, соединяющие средние значения каждого кластера. Данная иллюстрация показывает существенное улучшение качества сжатия изображения при использовании ДПХ вместо ДКП. На основании изложенных результатов следует, что использование дискретного преобразования Хартли позволяет получить метод сжатия изображений, который превышает JPEG как по степени сжатия, так и по качеству восстановленных изображений.
Очевидно, что \[ \sum_{i=0}^{N-1}p_i\mathrm{cas } \frac{2\pi\nu i}{N}=\sqrt{2}\sum_{i=0}^{N-1}p_i\cos \left(\frac{2\pi\nu i}{N}-\frac{\pi}{4}\right) \] то есть, преобразование Хартли представляет собой преобразование Фурье по косинусу с фазовым сдвигом π/4 и это дает неплохие результаты. Таким образом возникает вопрос, а существует преобразования по косинусу с любым углом фазового сдвига.

Одномерное дискретное тригонометрическое преобразования со свободным фазовым сдвигом.

Теорема. Пусть \(\varphi\in\left(0,\frac{\pi}{2}\right)\), тогда для всех \(\{h_m\}_{m=0}^{N-1}\) таких, что \(-\infty\lt h_m\lt \infty,m=0,1,...,N-1\) и \[ H_k=\sum_{m=0}^{N-1}h_m\cos\left(\frac{2\pi mk}{N}-\varphi\right) \] имеет место равенство \[ h_n=\frac{2}{N\sin(2\varphi)}\sum_{k=0}^{N-1}H_k\sin\left(\frac{2\pi nk}{N}+\varphi\right). \]

Доказательство. Подставим во второе соотношение значение \(H_k \) из первого равенства, тогда получим \[ \frac{2}{N\sin(2\varphi)}\sum_{k=0}^{N-1}{\left(\sum_{m=0}^{N-1}{h_m\cos\left(\frac{2\pi m k}{N}-\varphi\right)}\right)\sin\left(\frac{2\pi n k}{N}+\varphi\right)}= \] \[ \frac{2}{N\sin(2\varphi)}\sum_{m=0}^{N-1}{h_m\sum_{k=0}^{N-1}{\cos\left(\frac{2\pi m k}{N}-\varphi\right)}\sin\left(\frac{2\pi n k}{N}+\varphi\right)}= \] \[ \frac{2}{N\sin(2\varphi)}\sum_{m=0}^{N-1}{h_m\sum_{k=0}^{N-1}{\left(\sin\left(\frac{2\pi k(n+m)}{N}\right)+\sin\left(\frac{2\pi k (n-m)}{N}+2\varphi\right)\right)}}. \] Заметим, что (см., например, [8] стр.82) \begin{equation} \label{1} \sum_{\nu=0}^{N-1}\sin(\alpha+\nu\beta)=\frac{1}{\sin\frac{\beta}{2}}\sin\left(\alpha+\frac{N-1}{2}\beta\right)\sin\frac{N\beta}{2}, \end{equation} Пусть, прежде всего, \(n + m \ne N \), тогда из (\ref{1}) получим, \[ \sum_{k=0}^{N-1}\sin\left(\frac{2\pi k(n+m)}{N}\right)= \frac{1}{\sin\left(\frac{2\pi (n+m)}{2N}\right)}\sin\left(\frac{2\pi(n+m)(N-1)}{2N}\right) \sin\left(\frac{2\pi(n+m)N}{2N}\right). \] Из условия \(n+m\ne N\) получаем \[ \sin\left(\frac{2\pi (n+m)}{N}\right)\ne 0, \] а, с другой стороны, \[ \sin\left(\frac{2\pi N}{2N}\right)=\sin(\pi(n+m))=0, \] то есть, для \(n+m\ne N\) \[ \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi k(n+m)}{N}\right)}=0. \] Если же \(n+m= N\), тогда \[ \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi k(n+m)}{N}\right)}= \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi kN}{N}\right)}= \sum_{k=0}^{N-1}{\sin\left(2\pi k\right)}=0. \] Таким образом для всех \(n,m=0,1,\ldots,N-1\) \[ \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi k(n+m)}{N}\right)}=0. \] Рассмотрим соотношение \[ \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi k(n-m)}{N}+2\varphi\right)}. \] Используя равенство (\ref{1}), получаем \[ \sum_{k=0}^{N-1}{\sin\left(\frac{2\pi k(n-m)}{N}+2\varphi\right)}= \frac{1}{\sin\left(\frac{2\pi (n-m)}{2N}\right)}\sin\left(\frac{2\pi (n-m)(N-1)}{2N}+2\varphi\right)\sin\left(\frac{2\pi (n-m)N}{2N}\right). \] Далее, замечая, что так как \(n,m=0,1,\ldots,N-1\), то \(n-m\ne N\). Пусть, \(n\ne m\) , тогда получаем, что \[ \sin\left(\frac{\pi (n-m)}{N}\right)\ne 0, \] но при этом \[ \sin\left(\frac{2\pi (n-m)N}{2N}\right)=\sin(\pi(n-m))=0 \] для всех \(n,m=0,1,\ldots,N-1\) и для \(n\ne m\), в том числе. И, наконец, пусть \(n=m\), тогда \[ \sum_{k=0}^{N-1}\sin\left(\frac{2\pi k(n-m)}{N}+2\varphi\right)=\sum_{k=0}^{N-1}\sin(0+2\varphi)=N\sin(2\varphi). \] Таким образом, получаем \[ \sum_{k=0}^{N-1}\left(\sin\left(\frac{2\pi k(n+m)}{N}\right)+\sin\left(\frac{2\pi k(n-m)}{N}+2\varphi\right)\right)= \left\{ \begin{array}{ll} N\sin(2\varphi), &n=m; \\ 0, &n,m=0,1,\ldots,N-1,n\ne m. \end{array} \right. \] Откуда получаем соотношение \[ \frac{1}{N\sin(2\varphi)}\sum_{m=0}^{N-1}h_m\sum_{k=0}^{N-1}\left(\sin\left(\frac{2\pi k(n+m)}{N}\right)+\sin\left(\frac{2\pi k(n-m)}{N}+2\varphi\right)\right)= \left\{ \begin{array}{ll} h_n, &n=m; \\ 0, &n,m=0,1,\ldots,N-1,n\ne m. \end{array} \right. \] Что и завершает доказательство теоремы.
Отметим, что если \(\varphi = \frac {\pi}{4} \), то мы получим дискретное преобразование Хартли (ДПХ), то есть ДПХ является частным случаем полученного ДТП.
Использование возможности использования фазового сдвига для улучшения восстановления данных позволяет, на основании полученного дискретного преобразования, строить адаптивные фильтры, подстраивая фильтр не только для входных данных, но и, например, для используемого метода квантования или природы шума, вносящего искажения в сигнал.

Двумерное дискретное тригонометрическое преобразование со свободным фазовым сдвигом.

Так как одним из популярных сфер использования дискретного косинус-преобразования является цифровая обработка двухмерных сигналов, то есть изображений, то для использования полученного дискретного тригонометрического преобразования в обработке изображений применим двухмерное преобразование:

прямой ход \[ h_{i,j}=\frac{2}{N\sin(2\varphi)}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}p_{n,m}\cos\left(\frac{2\pi ni}{N}-\varphi\right)\cos\left(\frac{2\pi mj}{N}-\psi\right), \] и обратный ход \[ p_{n,m}=\frac{2}{N\sin(2\psi)}\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}h_{i,j}\sin\left(\frac{2\pi ni}{N}+\psi\right)\sin\left(\frac{2\pi mj}{N}+\varphi\right). \]

Используя это преобразование для блока 8 × 8 пикселей, можно получить лучший результат чем дискретное косинус-преобразования, но для получения такого результата необходимо выбрать соответствующие фазовые сдвиги. Этим дискретным тригонометрическим преобразованием можно заменить дискретное косинус-преобразования в стандарте JPEG без изменения логики кодирования и квантования, только заменив преобразования.

Определение лучшего фазового сдвига.

Для определения наилучшего фазового сдвига, состоящего для двухмерного сигнала, из двух параметров \(\varphi \) и \(\psi \), предложен следующий подход:

1. Нужно вычислить значение Peak Signal to Noise Ration (PSNR) для восстановленного изображения, при почти крайних значениях сдвигов, например \(\varphi_1 = \pi / 100, \varphi_2 = \pi / 8 \) и \(\varphi_3 = \pi / 4 \), значение \(\psi \) лучше взять среднее, то есть \(\pi / 8 \).
2. Найти точку максимума параболы проведенной через точки A (PSNR ₁ ; φ ₁ ), B (PSNR ₂ ; φ ₂ ) и C (PSNR ₃ ; φ ₃ ).
3. Координата X найденной точки будет аппроксимацией искомого значения φ.
4. Перейти к шагу 1 и провести эту же операцию для ψ используя найденное значение φ, данную операцию можно повторять несколько раз, используя поочередно φ и ψ, но практика показала, что и одного раза достаточно.

Этот метод также подходит и для одномерного сигнала. Данная операция занимает больше времени, чем применение Jpeg, так как изображение полностью кодируется и декодируется минимум 6 раз, но она позволяет получить почти наилучшие значения φ и ψ. Также можно рассчитать значение PSNR для всех значений φ и ψ с определенным шагом и выбрать наилучший результат, но приведенный метод позволяет значительно сократить количество вычислений преобразования и предоставляет более точный результат.
Для тестирования были взяты тестовые изображения из базы данных изображений TID2008 фирмы Kodak. Для каждого изображения были определены лучшие значения сдвигов, а также построены графики каждого значения PSNR в зависимости от фазовых сдвигов

Рассмотрим следующую задачу: как выбрать базисные векторы, чтобы минимальным числом базисных векторов можно было наилучшим образом приблизить данные из некоторого набора.
Решению этой задачи отвечает метод главных компонент - ортогональное линейное преобразование базиса, при котором первый вектор нового базиса соответствует направлению максимальной дисперсии данных, второй вектор - в следующем направлении максимальной дисперсии и так далее.

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

Метод главных компонент.

Метод главных компонент (МГК) - один из основных способов сокращения размерности данных с уменьшением количества информации, разработанный Карлом Пирсоном (Karl Pearson) в 1901 году. Используется во многих областях, таких как распознавание образов, компьютерное зрение, сжатие данных и др. (Более подробно с МГК можно ознакомиться, если перейти по ссылке). Нахождение главных компонент сводится к вычислению собственных векторов и значений ковариационной матрицы входных данных.
Перейдем к рассмотрению метода МГК (PCA - Principal Component Analysis). Сначала проведем центрирование данных, для чего найдем константу μ, которая наилучшим образом описывает входные данные \[ \varepsilon \left( \mu \right) = \sum\limits_{i = 1}^n {\left( {x_i - \mu } \right)^2 \to \mathop {\min }\limits_\mu } . \] Для нахождения минимума найдем производную и приравняем ее к нулю, после чего найдем значение μ, которое дает минимум \[ \frac{d}{d\mu }\varepsilon \left( \mu \right) = - 2\sum\limits_{i = 1}^n {\left( {x_i - \mu } \right) = 0 \Rightarrow } \sum\limits_{i = 1}^n {\mu = } \sum\limits_{i = 1}^n {x_i \Rightarrow } \mu n = \sum\limits_{i = 1}^n {x_i \Rightarrow } \mu = \frac{1}{n}\sum\limits_{i = 1}^n {x_i } . \] Далее проведем центрирование данных, то есть, переопределим входные значения \(x_{new} = x_{old} - \mu \), вычитая из каждого полученное значение μ. Ясно, что новые данные имеют математическое ожидание, равное нулю \[ E\left( {X - E\left( X \right)} \right) = E\left( X \right) - E\left( X\right) = 0. \] Таким образом, мы сделали параллельный сдвиг в существующий системе координат.
В дальнейшем будем считать, что наши данные центрированы и мы хотим найти наиболее точное представление данных \(D = \left\{ {x_1 ,...,x_n } \right\}\) в некотором подпространстве W с размерностью \(k \lt n \).
Пусть \(\left \{ {e_1 ..., e_k} \right \} \) ортогональный базис W . Произвольный вектор из W может быть записан в виде линейной комбинации векторов базиса, то есть, \(x_1 \) можно поставить в соответствие некоторый вектор \(\sum \limits_{i = 1}^k {\alpha_{1, i}} e_i \) из W . Погрешность между ними находится следующим образом \[ \varepsilon_1 = \left\| {x_1 - \sum\limits_{i = 1}^k {\alpha_{1,i} } e_i } \right\|_2^2 = \left\langle {x_1 - \sum\limits_{i = 1}^k {\alpha_{1,i} } e_i ,x_1 - \sum\limits_{i = 1}^k {\alpha_{1,i} } e_i } \right\rangle . \] Иллюстрация погрешности восстановления вектора. Чтобы найти полную погрешность, найдем сумму величин погрешностей по всем \(x_j \) \begin{equation} \label{2.1} \varepsilon \underbrace {\left( {e_1 ,...,e_k ,\alpha_{1,1} ,...,\alpha_{n,k} } \right)}_{unknowns} = \sum\limits_{j = 1}^n {\varepsilon_j } = \sum\limits_{j = 1}^n {\left\| {x_j - \sum\limits_{i = 1}^k {\alpha_{j,i} } e_i } \right\|_2^2 } . \end{equation} Чтобы минимизировать погрешность, нужно взять частные производные и учесть ограничения на ортогональность \(\left\{ {e_1 ,...,e_k } \right\}\).
Вначале упростим соотношение (\ref{2.1}) \[ \varepsilon \left( {e_1 ,...,e_k ,\alpha_{1,1} ,...,\alpha_{n,k} } \right) = \sum\limits_{j = 1}^n {\left\| {x_j - \sum\limits_{i = 1}^k {\alpha_{j,i} } e_i } \right\|_2^2 = }\] \[= \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - 2\sum\limits_{j = 1}^n {x_j^T \sum\limits_{i = 1}^k {\alpha_{j,i} } e_i + \sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\alpha_{j,i}^2 } = } } } \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - 2\sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\alpha_{j,i} x_j^T } e_i + \sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\alpha_{j,i}^2 } .} } } \] Тогда \[ \frac{\partial }{\partial \alpha_{m,\ell} }\varepsilon \left( {e_1 ,...,e_k,\alpha_{1,1} ,...,\alpha_{n,k} } \right) = - 2x_m^T e_\ell + 2\alpha_{m,\ell}. \] Необходимое и достаточное условие экстремума будет иметь вид \[ - 2x_m^T e_\ell + 2\alpha_{m,\ell} = 0 \Rightarrow \alpha_{m,\ell} = x_m^T e_\ell . \] Таким образом, погрешность (\ref{2.1}) будет иметь вид \[ \varepsilon \left( {e_1 ,...,e_k } \right) = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - 2\sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\left( {x_j^T e_i } \right)x_j^T } e_i + \sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\left( {x_j^T e_i } \right)^2} .} } } \] Проводя упрощения, имеем \begin{equation} \label{2.2} \varepsilon \left( {e_1 ,...,e_k } \right) = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{j = 1}^n {\sum\limits_{i = 1}^k {\left( {x_j^T e_i } \right)^2} .} } \end{equation} Учитывая, что \(\left\langle {a,b} \right\rangle = a^Tb\) и \(\left\langle {b,a} \right\rangle = \left\langle {a,b} \right\rangle \), получаем \[ \left( {a^Tb} \right)^2 = \left( {a^Tb} \right)\left( {a^Tb} \right) = \left( {b^Ta} \right)\left( {a^Tb} \right) = b^T\left( {aa^T} \right)b, \] то есть, \[ \varepsilon \left( {e_1 ,...,e_k } \right) = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{i = 1}^k {e_i^T \left( {\sum\limits_{j = 1}^n {\left( {x_j x_j^T } \right)} } \right)e_i } } = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{i = 1}^k {e_i^T Se_i } } , \] где \(S = \sum\limits_{j = 1}^n {\left( {x_j x_j^T } \right)}\) ковариационная матрица.
Следующим шагом будет минимизация \[\varepsilon \left( {e_1 ,...,e_k } \right) = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{i = 1}^k {e_i^T Se_i } }\] при условии \(e_i^T e_i = 1 \) для всех i . Для решения этой задачи будем использовать метод неопределенных множителей Лагранжа (Lagrange). Введем множители \(\lambda_1 ,...,\lambda_k\), отмечая, что \(\sum\limits_{j =1}^n {\left\| {x_j } \right\|_2^2 } \equiv Const\) и выпишем функцию цели \[ L\left( {e_1 ,...,e_k } \right) = \sum\limits_{i = 1}^k {e_i^T Se_i } - \sum\limits_{i = 1}^k {\lambda_i \left( {e_i^T e_i - 1} \right)} . \] Замечая, что \(\frac{d}{dX}\left( {X^TX} \right) = \frac{d}{dX}\left\langle {X,X} \right\rangle = 2X\) и если A симметричная матрица, то \(\frac{d}{dX}\left( {X^TAX} \right) = 2AX\), имеем, \[ \frac{\partial }{\partial e_m }L \left( {e_1 ,...,e_k } \right) = 2Se_m - 2\lambda_m e_m = 0, \] то есть, \(Se_m = \lambda_m e_m\). Таким образом, необходимо найти решение уравнения \(\left( {S - \lambda I} \right)e = 0\) (где I единичная матрица), что эквивалентно тому, что \(\lambda_m \) и \(e_m\) есть собственные значения и собственные векторы ковариационной матрицы S.
При этом погрешность имеет вид \begin{equation}\label{2.3} \varepsilon \left( {e_1 ,...,e_k } \right) = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{i = 1}^k {\lambda_i\left\| {e_i } \right\|_2^2 } } = \sum\limits_{j = 1}^n {\left\| {x_j } \right\|_2^2 - \sum\limits_{i = 1}^k {\lambda_i} } . \end{equation} Минимизация (\ref{2.3}) состоит в выборе базиса W из k собственных векторов матрицы S , которые соответствуют k крупнейшим собственным значениям. Большее собственное значение S дает большую вариацию в направлении соответствующего собственного вектора. Этот результат можно переформулировать следующим образом - проекция X на подпространство размера k , которая обеспечивает наибольшую вариацию. Таким образом МГК может трактоваться следующим образом - берем ортогональный базис и вращаем его пока на одном из направлений не получим максимальную вариацию. Фиксируем это направление и вращаем другие, пока не найдем второе направление и так далее.
Пусть \(\left \{ {e_1 ..., e_n} \right \} \) все собственные векторы матрицы S , отсортированные в порядке уменьшения соответствующего собственного значения, тогда для произвольного \[ x_i = \sum\limits_{j = 1}^n {\alpha_{i,j} } e_j = \underbrace {\alpha _{i,1} e_1 + ... + \alpha_{i,k} e_k }_{approximation} + \overbrace {\alpha _{i,k + 1} e_{k + 1} + ... + \alpha_{i,n} e_n }^{error} \]
коэффициенты \(\alpha_{m, \ell} = x_m^T e_\ell \) являются координатами главных компонент, и чем больше значение k , тем лучшую аппроксимацию получаем. Приведем алгоритмизацию МГК.
Пусть \(D = \left \{ {x_1^0, ..., x_n^0} \right \} \) входные (оригинальные) данные, каждый из этих векторов \(x_i^0 \) имеет размер N

1. Найдем среднее \(\mu = \frac{1}{n} \sum \limits_{i = 1}^n {x_i^0} \).
2. Вычтем среднее по каждому вектору \(x_i = x_i^0 - \mu \).
3. Найдем ковариационную матрицу \(S = \sum \limits_{j = 1}^n {x_jx_j^T} \).
4. Вычислим собственные векторы \(\left \{ {e_1 ..., e_k} \right \} \), которые соответствуют k самым большим собственным значениям S.
5. Пусть \(\left \{ {e_1 ..., e_k} \right \} \) составляют матрицу \(E = \left [{e_1 ... e_k} \right] \).
6. Тогда ближе всего к x лежит \(z^TE + \mu \).

Итерационный алгоритм вычисления главных компонент.

Описанный метод получения главных компонент достаточно ресурсоемок и нестабильный, особенно в случае, если собственные значения матрицы близки к нулю.
Поэтому более эффективным является использование итерационного метода нахождения главных компонент. Для этого рассмотрим задачу (\ref{2.1}) с другой точки зрения.
В случае i = 1 задача (\ref{2.1}) сводится к нахождению одной компоненты \(e_1 \), которая наилучшим образом восстанавливает все исходные данные \(\left \{ {x_1, ... , x_n} \right \} \). \begin{equation} \label{2.4} \varepsilon \left( {e_1 ,\alpha_{1,1} ,...,\alpha_{n,1} } \right) = \sum\limits_{j = 1}^n {\left\| {x_j - \alpha_{j,1} e_1 } \right\|_2^2 } \to \min \end{equation} по всем \(e_1 \) i \(\left\{ {\alpha_{i,1} } \right\}_{i = 1}^n \) при условии \(\sum\limits_{i = 1}^n {\alpha_{i,1}^2 = 1.} \)
Если \(\left\{ {\tilde{\alpha }_{i,1} } \right\}_{i = 1}^n \) и \(\tilde{e}_1\) есть решение этой задачи, то \(\Delta x_j = x_j - \tilde {\alpha }_{j,1}\tilde {e}_1 \) - погрешность восстановления данных одной первой главной компонентой, тогда, решая задачу \[ \sum\limits_{j = 1}^n {\left\| {\Delta x_j - \alpha_{j,2} e_2 } \right\|_2^2 } \to \min \] по всем \(e_2\) та \(\left\{ {\alpha_{i,2} } \right\}_{i = 1}^n \) при условии \(\sum\limits_{i = 1}^n {\alpha_{i,1}^2 = 1} \), получаем вторую главную компоненту \(\tilde {e}_2 \) и соответствующий вектор \(\left\{ {\tilde {\alpha }_{i,2} } \right\}_{i = 1}^n \) и т.д.
При фиксованных \(\left\{ {\alpha_{i,1} } \right\}_{i = 1}^n \) задача (\ref{2.4}) решается методом наименьших квадратов. В силу того, что функцией цели является квадратичный функционал, необходимое и достаточное условие экстремума совпадают. Таким образом, решение задачи сводится к поиску решения уравнения \[ \frac{\partial }{\partial e_1 }\varepsilon \left( {e_1 ,\alpha_{1,1} ,...,\alpha_{n,1} } \right) = - 2\sum\limits_{j = 1}^n {\left( {x_j - \alpha_{j,1} e_1 } \right)\alpha_{j,1} } = - 2\left( {\sum\limits_{j =1}^n {x_j \alpha_{j,1} } - \sum\limits_{j = 1}^n {\alpha_{j,1}^2 e_1 } } \right). \] Отсюда имеем \[ \frac{\partial }{\partial e_1 }\varepsilon \left( {e_1 ,\alpha_{1,1} ,...,\alpha_{n,1} } \right) = - 2\sum\limits_{j = 1}^n {\left( {x_j - \alpha_{j,1} e_1 } \right)\alpha_{j,1} } = - 2\left( {\sum\limits_{j =1}^n {x_j \alpha_{j,1} } - \sum\limits_{j = 1}^n {\alpha_{j,1}^2 e_1 } } \right). \] учитывая условие нормирования единицей, то есть \(\sum\limits_{i = 1}^n {\alpha_{i,1}^2 = 1}\), имеем \[ e_1 = \sum\limits_{j = 1}^n {x_j \alpha_{j,1} } . \] Следующий шаг будем делать из условия, что в задаче (\ref{2.4}) нам известна компонента \(e_1 \) и требуется найти экстремум по \(\left\{ {\alpha_{i,1} } \right\}_{i = 1}^n\) \[ \frac{\partial }{\partial \alpha_{\nu ,1} }\varepsilon \left( {e_1 ,\alpha_{1,1} ,...,\alpha_{n,1} } \right) = - 2\left( {x_\nu - \alpha_{\nu ,1} e_1 } \right)e_1 = - 2\left( {\left\langle {x_\nu ,e_1 } \right\rangle - \alpha_{\nu ,1} \left\langle {e_1 ,e_1 } \right\rangle } \right) = 0, \] то есть \[ \alpha_{\nu ,1} = \frac{\left\langle {x_\nu ,e_1 } \right\rangle }{\left\langle {e_1 ,e_1 } \right\rangle }, \] где, как обычно, \(\left\langle {x,y} \right\rangle\)- скалярное произведение векторов x и y.
Далее, считая найденные \(\left\{ {\alpha_{i,1} } \right\}_{i = 1}^n \) уже известными, повторяем весь процесс, пока не достигнем стабилизации погрешности.
Полученные \(e_1 \) будем считать первой главной компонентой \(\tilde {e}_1 \). Тогда \(\Delta x_j = x_j - \tilde {\alpha }_{j,1} \tilde {e}_1\) - погрешность восстановления данных одной первой главной компонентой.
Применяя этот метод к погрешности восстановления \(\Delta x_j \), находим вторую главную компоненту \(e_2 \) вместе с коэффициентами \(\alpha_{j, 2} \), и так далее.
Приведем алгоритмизацию этого алгоритма.
Во-первых центрируем данные, вычитая из входных данных среднее значение и в дальнейшем считаем, что данные в среднем равны нулю.

1. Пусть номер итерации ν = 1.
2. Выбираем стартовые значения \(\left \{ {\alpha_{i, 1}^\nu} \right\}_{i = 1}^n \), например, пусть все они между собой равны, то есть \(\alpha_{i, 1}^\nu = \frac{1}{\sqrt n}, i = 1,2, ..., n \).
3. Вычислим \(e_1^\nu = \sum \limits_{j = 1}^n {x_j \alpha_{j, 1}^\nu} \).
4. Далее найдем \(\beta_i = \frac {\left \langle {x_i, e_1^\nu} \right \rangle} {\left \langle {e_1^\nu, e_1^\nu} \right \rangle} \) и, после нормирования единицей, имеем \[ \alpha_{i,1}^{\nu + 1} = \frac{\beta_i }{\sqrt {\sum\limits_{j = 1}^n {\beta_j^2 } } }. \]
5. Пусть ν = ν+ 1.
6. Проведем проверку критерия остановки, в качестве которого может быть либо стабилизация коэффициентов \(\left\{ {\alpha_{i,1}^\nu } \right\}_{i = 1}^n\), либо стабилизация главной компоненты \(e_1^\nu \) или проверка на заданное фиксированное число итераций. Если условие остановки итерационного процесса не выполнено, то переходим к пункту 3.

Переход от модели RGB к оптимальной трехкомпонентной модели.

Перейдем к использованию метода главных компонент в компьютерной графике. Здесь важную роль играет конвертация изображения из пространства равноправных цветных характеристик в пространство неравноправных (см. раздел "Свет и цвет"). Произвольное изображение визуализируется с использованием смешивания равноправных цветных компонент - красной, зеленой и синей составляющих - модель RGB. В качестве неравноправных цветных компонент, как правило, используют, соответственно, также три компонента - значение освещенности (люминесцентную составляющую), характеристику теплых тонов и характеристику холодных тонов. Использование неравноправных цветных компонент используют для сжатия изображений и видео-потоков, применяя к каждой из неравноправных компонент свой метод сжатия.
Можно подходить к проблеме построения неравноправного цветового пространства с другой точки зрения, исходя из максимальной информативности каждой компоненты. Применим метод главных компонент для получения неравноправной трехкомпонентной модели оптимальной с точки зрения минимизации среднеквадратичной погрешности восстановления данного изображения. Таким образом, первая из полученных цветных компонент будет нести наибольшую информацию относительно изображения между всех полученных цветных компонент, а вторая будет иметь наибольшую информацию среди тех, кто остались.
Таким образом, в нашей терминологии, задача (\ref{2.1}) будет иметь вид \[ \left\| {R - \sum\limits_{i = 1}^3 {\alpha_{r,i} } e_i } \right\|_2^2 + \left\| {G - \sum\limits_{i = 1}^3 {\alpha_{g,i} } e_i } \right\|_2^2 + \left\| {B - \sum\limits_{i = 1}^3 {\alpha_{b,i} } e_i } \right\|_2^2 \to \min , \] где минимум берется по всем \(\alpha_{r, i} \alpha_{g, i} \alpha_{b, i} \) и \(e_i \quad i = 1,2,3. \) В качестве данных рассмотрим тестовое изображение Lena. Применяя метод главных компонент, имеем \[ \left( { {\begin{array}{*{20}c} R \hfill \\G \hfill \\B \hfill \\ \end{array} }} \right) = \left( {\begin{array}{l} \mbox{ 0.767785 0.45439 0.4517034} \\ \mbox{ - 0.6164395 0.716085 0.3274513} \\ \mbox{ - 0.174667 - 0.5298553 0.8299064} \\ \end{array}} \right)\left( { {\begin{array}{*{20}c} {e_1 } \hfill \\ {e_2 } \hfill \\ {e_3 } \hfill \\ \end{array} }} \right). \] Здесь в первом столбце матрицы стоят коэффициенты \(\alpha_{r, 1}, \alpha_{r, 2}, \alpha_{r, 3} \), во втором - \(\alpha_{g, 1} , \alpha_{g, 2}, \alpha_{g, 3} \) и в третьем - \(\alpha_{b, 1}, \alpha_{b, 2}, \alpha_{b, 3} \). Тогда восстановления тестового изображения по одной компоненте можно записать в виде \[ R_{i,j}=\mbox{0.767785 }Y_{i,j} , \\ G_{i,j}=\mbox{ 0.45439} Y_{i,j}, \\ B_{i,j}=\mbox{0.4517034} Y_{i,j}, \] где \(Y_{i,j}\) – значение первой главной компоненты \(e_1 \), которые соответствуют пикселю с координатами (i, j) . Сравним восстановления изображения компонентами различных пространств цветов Синим цветом обозначена погрешность восстановления изображения по одной, по двум или полная реконструкция тремя компонентами пространства YUV, желтым цветом -YСrCb, а пурпурным - по главным компонентам.

Оптимальный метод кодирования изображения.

В рассмотренном выше примере изображение было разложено на три оптимальных (с точки зрения энергии слоев) информационных домена (слоя). Обобщим данный подход используя для разбиения изображения на \(3 \times n^2 (n \ge 1) \) информационных доменов с возможностью динамической подгрузки (по аналогии с прогрессивным Jpeg). То есть, реализовать оптимальный метод прогрессивной развертки изображения на \(3 \times n^2 (n \ge 1) \) слоев.
Суть метода состоит в следующем. Разобьем наше изображение размером \(H \times W \) на квадраты \(P_{x, y} \) размером \(n \times n (n \ge 1) \). Для простоты будем считать, что \(H = N \times n, W = M \times n \). Каждому пикселю \(p_{nx + i, ny + j} (i, j = 0,1, ..., n-1) \) данного квадрата соответствуют три цветные компоненты \(r{nx + i, ny + j}, g_{nx + i, ny + j}, b_{nx + i, ny + j} \) (красная, зеленая и синяя составляющие). То есть, общая информация на квадрате \(P_{x, y} \) описывается данным \[ (r_{nx+i,ny+j},g_{nx+i,ny+j},b_{nx+i,ny+j}),(r_{nx+i,ny+j},g_{nx+i,ny+j},b_{nx+i,ny+j}),...,(r_{nx+i,ny+j},g_{nx+i,ny+j},b_{nx+i,ny+j}). \] Сформируем домены входящей информации, например, следующим образом \[ X_0=\{r_{0,0},r_{n,0},...,r_{nx,ny},...,r_{H-n-1,W-n-1}\}, \] \[ X_1=\{g_{0,0},g_{n,0},...,g_{nx,ny},...,g_{H-n-1,W-n-1}\}, \] \[ X_2=\{b_{0,0},b_{n,0},...,b_{nx,ny},...,b_{H-n-1,W-n-1}\}, \] \[ X_3=\{r_{1,0},r_{n+1,0},...,r_{nx+1,ny},...,r_{H-n,W-n-1}\}, \] \[ X_4=\{g_{1,0},g_{n+1,0},...,g_{nx+1,ny},...,g_{H-n,W-n-1}\}, \] \[ X_5=\{b_{1,0},b_{n+1,0},...,b_{nx+1,ny},...,b_{H-n,W-n-1}\}, \] \[\ldots\] \[ X_{3\times n^2-2}=\{r_{n-1,n-1},r_{2n-1,2n-1},...,r_{(n+1)x-1,(n+1)y-1},...,r_{H,W}\}, \] \[ X_{3\times n^2-1}=\{g_{n-1,n-1},g_{2n-1,2n-1},...,g_{(n+1)x-1,(n+1)y-1},...,g_{H,W}\}, \] \[ X_{3\times n^2}=\{b_{n-1,n-1},b_{2n-1,2n-1},...,b_{(n+1)x-1,(n+1)y-1},...,b_{H,W}\}. \] Рассмотрим задачу \[ \sum_{j=0}^{3\times n^2}\left\|X_j-\sum_{i=0}^k\alpha_{j,i}Y_j\right\|^2_2\to\min, \] где минимум берется по всем \(\left\{Y_i\right\}_{i=0}^k, \alpha_{j,i}(i=0,1,\ldots,k,j=0,1,...,3\times n^2).\)
Как видно, для решения этой задачи можно использовать итерационный метод главных компонент, алгоритм которого приведен выше.
Таким образом, пусть мы получили все \(\left\{Y_i\right\}_{i=0}^k, \alpha_{j,i}(i=0,1,\ldots,k,j=0,1,...,3\times n^2).\) Опишем алгоритм прогрессивного восстановления исходных данных.
Понятное дело, что множество \(\left \{X_j \right \}_{j = 0}^{3 \times n^2} \) описывает все пиксели изображения, то есть, если \[X^0_j=\alpha_{j,0}Y_0,\] то \(\left\{X_j^0\right\}_{j=0}^{3\times n^2}\) является первой итерацией восстановления нашего изображения. На втором шаге \( X^1_j=\alpha_{j,0}Y_0+\alpha_{j,1}Y_1 \), то есть \(X^1_j=X^0_j+\alpha_{j,1}Y_1\) и \(\left\{X_j^1\right\}_{j=0}^{3\times n^2}\) будет второй итерацией. Таким образом имеем \(X^k_j=\sum_{i=0}^k\alpha_{j,i}Y_i=X^{k-1}_j+\alpha_{j,k}Y_k\), а \(\left\{X_j^k\right\}_{j=0}^{3\times n^2}\) будет k -й итерацией восстановления нашего изображения. Таким образом имеем \(X^k_j=\sum_{i=0}^k\alpha_{j,i}Y_i=X^{k-1}_j+\alpha_{j,k}Y_k\), а \(\left\{X_j^k\right\}_{j=0}^{3\times n^2}\) будет k -й итерацией восстановления нашего изображения. Отметим, что в отличие от фильтрации на основе преобразования Фурье, форма примитива, по которому делается разложение методом главных компонент, может быть любой формы при условии плотной упаковки плоскости, очевидно, что это может быть прямоугольник, шестиугольник, который наиболее плотно упаковывает плоскость, или даже такая экзотическая фигура, как фрактал-дракон.

Приведем некоторые примеры.

Пусть задан любой двумерный паркет, покрывающий всю плоскость. Каждый элемент паркета состоит из n пикселей. Перенумеруем все пиксели элемента паркета произвольным, но фиксированным способом. Выбирая из каждого элемента паркета пиксель с одним и тем же номером, получим n компонент. Заметим, что результат оптимальной фильтрации не зависит от способа нумерации точек на элементе паркета.
Для простейшего случая - квадратного паркета имеем Например, для паркета с элементами в виде шестиугольных сот с 76 точек можно использовать следующую нумерацию Применим МГК к изображению Lena в градациях серого Выпишем нумерацию для паркета «Дракон» с 64 пикселей Сделаем сравнения качества фильтрации по разному паркета размером 64 пикселя. Фиолетовым цветом обозначена погрешность восстановления с использованием фильтрации по паркету «дракон», цветом «морской волны» - по шестиугольным сотам, темно-синим - по квадратам 8 × 8, желтым - по прямоугольной решетке 16 × 4 и красным - по решетке 4 × 16 . Приведем фрагмент этих графиков в другом масштабе.

Вейвлеты в компьютерной графике.

Всплески (wavelets) появились в начале 1980-х годов на стыке теории функций, функционального анализа, обработки сигналов и изображений, квантовой теории поля. Развитие теории всплесков связано с именами А. Гроссмана, Ж. Морлета, Ж. Стремберга, Я. Мейера, С. Малла, И. Добеши и многих др. (Более подробно с вейвлетами можно ознакомиться, если перейти по ссылке)
Подобно анализу Фурье теория всплесков находит применение в исследовании частотных характеристик функции f с помощью всплеск-преобразования. При описании функций с конечным носителем, сигналов, с переменными во времени частотами с помощью анализа Фурье, возникают трудности, которые объясняются свойствами тригонометрических функций на оси \(\mathbb{R} \). Например, преобразование Фурье функции f в случае δ-функции Дирака, которая имеет в качестве носителя единую точку, вырождается в функцию \(\hat{\delta} (\omega) = 1 \), распространенную на всю числовую ось. При обработке сигналов с переменными частотами традиционный анализ Фурье не позволяет выделить частоты, составляющие сигнал в окрестности произвольного момента времени. Аппарат всплесков позволяет автоматически осуществлять локализацию как по времени, так и в частотной области.
Приведем один алгоритм фильтрации векторного поля, то есть набора векторов одинаковой размерности.
Эти результаты являются основополагающими для всех алгоритмов, связанных с оптимальной фильтрацией фильтрами, подобными фильтрам Хаара и Уолша-Адамара одной и многих переменных.
В численном и функциональном анализе дискретные вейвлет-преобразования (ДВП) относятся к преобразованиям, в которых вейвлеты представлены дискретными сигналами (выборками).
Впервые ДВП было придумано венгерским математиком Альфредом Хааром. Для входного сигнала, представленного массивом 2 ⁿ чисел, вейвлет-преобразование Хаара просто группирует элементы по 2 и образует из них суммы и разности. Группировка сумм проводится рекурсивно для следующего уровня разложения. В итоге получается 2 ⁿ -1 разность и одна общая сумма.
Это простое ДВП иллюстрирует общие полезные свойства вейвлетов. Во-первых, преобразование можно выполнить используя \(n \log_2n \) операций. Во-вторых, оно не только раскладывает сигнал на некоторое подобие частотных полос (путем анализа его в разных масштабах), но и представляет временную область, то есть моменты возникновения тех или иных частот в сигнале.
Приведем определение и обозначение, необходимые нам в дальнейшем.
Множество векторов \(X^k=\left\{x_i^k\right\}(k = 0, .., n-1, i = 0,1, .., N-1)\) одинаковой розмерности N будем называть компонентами векторного поля или просто компонентами. Под фильтром \(H^\nu=\{h_i^\nu \} (\nu = 0, .., n-1)\) будем понимать вектор с n координатами.
Систему векторов \(H^\nu (\nu = 0, .., n-1)\) будем называть ортонормованной, если \[ H^\nu\bot H^\mu \Leftrightarrow \left\lt H^\nu,H^\mu \right\gt=0 \hbox{ и } \left\lt H^\nu,H^\nu \right\gt=1. \] Здесь <a, b> традиционное скалярное произведение векторов a и b.
Под фильтрацией компонентов \(X^k (k = 0, .., n-1)\) системой фильтров \(H^k (k = 0, .., n-1)\) будем понимать вычисления множества (частотного домена) по правилу \begin{equation}\label{w:1} Y^\nu=\sum_{k=0}^{n-1}h_\nu^k X^k. \end{equation} В координатной форме это соотношение будет иметь вид \[ y_i^\nu=\sum_{k=0}^{n-1}h_\nu^k x_i^k. \] Величину \begin{equation}\label{w:2} X_m^\nu=\sum_{k=0}^{m-1}h_\nu^k Y^k \end{equation} будем называть восстановлением компоненты \(X^\nu\) по m частотным доменам. Если m = n, то получаем восстановление компоненты по полному набору частотных доменов.
Первое соотношение называется декомпозицией, а последнее - реконструкцией данных.
В случае, если для данного массива \(F = \{f_i \} \) и для n = 4 компоненты определены правилом \[ X^0=\{x_i^0 \}=\{f_{4i} \},X^1=\{x_i^1 \}=\{f_{4i+1} \},\\ X^2=\{x_i^2 \}=\{f_{4i+2} \},X^3=\{x_i^3\}=\{f_{4i+3} \}, \] для фильтров Уолша-Адамара \[ H^0=0.5 \{1,1,1,1\},H^1=0.5 \{1,1,-1,-1\},\\ H^2=0.5\{1,-1,1,-1\},H^3=0.5\{1,-1,-1,1\} \] процедура декомпозиции (\ref{w:1}) будет иметь вид \[ y_i^0=0.5 (f_{4i}+f_{4i+1}+f_{4i+2}+f_{4i+3} ),\\ y_i^1=0.5(f_{4i}+f_{4i+1}-f_{4i+2}-f_{4i+3} ),\\ y_i^2=0.5 (f_{4i}-f_{4i+1}+f_{4i+2}-f_{4i+3} ),\\ y_i^3=0.5 (f_{4i}-f_{4i+1}-f_{4i+2}+f_{4i+3} ), \] а формулы реконструкции (\ref{w:2}), например, для m = 2 \[ f_{4i}=0.5 (y_i^0+y_i^1 ),\\ f_{4i+1}=0.5 (y_i^0+y_i^1 ),\\ f_{4i+2}= 0.5 (y_i^0-y_i^1 ),\\ f_{4i+3}=0.5(y_i^0-y_i^1 ). \] Полная реконструкция (то есть m = n ) определяется соотношениями \[ f_{4i}=\frac{1}{2} (y_i^0+y_i^1+y_i^2+y_i^3 ),\\ f_{4i+1}=\frac{1}{2}(y_i^0+y_i^1-y_i^2-y_i^3 ),\\ f_{4i+2}=\frac{1}{2}(y_i^0-y_i^1+y_i^2-y_i^3 ),\\ f_{4i+3}=\frac{1}{2}(y_i^0-y_i^1-y_i^2+y_i^3 ). \] Таким образом, задача использования фильтров разбивается на две. Первая из которых - сведение сигнала (в общем случае многокомпонентного и многоразмерного) к системе компонент одинаковой размерности. В частности, для фильтрации изображений это 3-компонентный двумерный сигнал, для видео - трехкомпонентный трехмерный и так далее. Вторая задача состоит в получении системы лучших (в том или ином смысле) фильтров для n компонент.
Фундаментальным понятием теории вейвлет есть понятие кратномасштабного анализа (КМА).
Кратномасштабный анализ в \(L^2 (\mathbb{R}^n) \) - это последовательность замкнутых подпространств \[ \ldots \subset V^{-1}\subset V^0 \subset V^1 \subset \] для которых выполняются условия

1. \(\overline{\bigcup_{j\in Z}V^j}=L^2 (\mathbb{R}^n );\)
2. \(\bigcap_{j\in Z}V^j=\{0\};\)
3. \(f\in V^j \Leftrightarrow f\left(2^{-j\circ}\right)\in V^0;\)
4. Найдется такая функция \(\varphi\in V^0\)(масштабирующая функция), множество ее сдвигов φ (x-n) создает ортогональный базис пространства \(V^0\).

Данная последовательность пространств гарантирует, что для произвольной интегрированной функции f и заранее заданного числа \(\varepsilon \gt 0 \) существует пространство \(V^k \) и элемент \(F_k \in V^k \) такие, что \(\| f-f_k \| \lt \varepsilon \). Отметим, что в евклидовом пространстве лучшим приближением функции является ортогональная проекция функции на это пространство, поэтому в качестве \(f_k \) будем брать ортогональную проекцию f на это пространство.
Кроме того, из свойств КМА следует, что для любого \(j \in Z \) множество функций \(\varphi^j_k(x)=2^{j/2}\varphi(2^jx-k)\) \(k=0,\pm 1, \ldots\) создают ортогональный базис. Множитель \(2^{j / 2} \) является результатом того, что замена x на 2x уменьшает норму в \(2^{1/2} \) раз .
Так как \(V^0 \subset V^1 \), то функция φ является линейной комбинацией функций \(\{\varphi^1_n \} \). Таким образом, найдется множество коэффициентов \(\{h_n \} \) таких, что \begin{equation}\label{w.1} \varphi(x)=\sqrt{2}\sum_{n\in Z}h_n\varphi(2x-n), \end{equation} Из равенства Парсеваля имеем \(\sum_nh^2_n = 1. \)
Функция φ (x) называется масштабирующей (скелинг) функцией, а функция \begin{equation}\label{w.2} \psi(x)=\sqrt{2}\sum_{n\in Z}(-1)^{1-n}h_{1-n}\varphi(2x-n), \end{equation} называется вейвлетом или всплеском.
Учитывая, что \(V^{- 1} \subset V^0 \), можно записать \(V^{- 1} \oplus W^{- 1} = V^0. \) В качестве всплеска, который определяет сдвиги базиса в пространствах \(W^j \) возьмем функцию ψ и положим \[ \psi_n^j (x)=2^{j/2} \psi (2^j x-n). \] В таком случае \[ \sum_{n\in Z}c_n^{-1}\varphi_n^{-1}(x)+ \sum_{n\in Z}d_n^{-1}\psi_n^{-1}(x)= \sum_{n\in Z}c_n^{0}\varphi_n^{0}(x). \] Отсюда, применяя операцию скалярного произведения, получаем \[ c^{-1}_k=\sum_{n\in Z}c_n^0\left\lt \varphi^0_n(x),\varphi^{-1}_k(x)\right>= \sum_{n\in Z}c_n^0\left\lt \varphi^0_n(x), \sum_{m\in Z}h_{m-2k}\varphi^{-1}_m(x)\right>= \sum_{n\in Z}c_n^0{h_{n-2k}}. \] и \[ d^{-1}_k= \sum_{n\in Z}c_n^0{g_{n-2k}}. \] где \(g_n=(-1)^{1-n} h_{1-n}\).
Ясно, что в общем случае \(V^j \oplus W^j=V^{j+1}\) и \[\sum_{n\in Z}c_n^{j}\varphi_n^{j}(x)+ \sum_{n\in Z}d_n^{j}\psi_n^{j}(x)= \sum_{n\in Z}c_n^{j+1}\varphi_n^{j+1}(x)\] и общие формулы декомпозиции \[ c^{j}_k= \sum_{n\in Z}c_n^{j+1}{h_{n-2k}}, \] и \[ d^{j}_k= \sum_{n\in Z}c_n^{j+1}{g_{n-2k}}. \] Обобщая выше сказанное, имеем разложение \(V^0 \) в ортогональную сумму \[ V^0=W^{-1}\oplus W^{-2}\oplus \ldots W^{-k}\oplus V^{-k} \] и формула реконструкции на каждом шаге будет иметь вид \[ c^{j+1}_k= \sum_{n\in Z}c_n^{j}h_{k-2n}+\sum_{n\in Z}d_n^{j}g_{k-2n}.\] Приведенный алгоритм называется каскадным.

Ортогональные вейвлеты.

Пусть \(\varphi \) является ортогональной функцией на решетке с шагом 2, то есть \[ \left\lt \varphi,\varphi(\cdot-2k)\right>=\delta_{0,k} \] Или, что то же самое, \[ \sum_ih_ih_{i-2k}=\delta_{0,k}.\] Известен факт, что если ортогональная функция имеет компактный носитель, то этот носитель четный. Для простейшего случая он равен двум \[ h_0=\frac{1}{\sqrt{2}},h_1=\frac{1}{\sqrt{2}}. \] Соответствующая базовая функция называется функцией Хаара Найдем тензорное произведение различных вариантов этих функций\(\varphi(x)\varphi(y),\varphi(x)\psi(y),\psi(x)\varphi(y)\) и \(\psi(x)\psi(y)\), в результате получим четверку фильтров (фильтров Хаара) \[ H^{++}=\frac{1}{2}\left( \begin{array}{cc} 1 &1\\ 1 &1 \\ \end{array} \right), H^{\pm}=\frac{1}{2}\left( \begin{array}{rr} 1 &1\\ -1 &-1 \\ \end{array} \right), H^{+-}=\frac{1}{2}\left( \begin{array}{cc} 1 &-1\\ 1 &-1 \\ \end{array} \right), H^{--}=\frac{1}{2}\left( \begin{array}{rr} 1 &-1\\ -1 &1 \\ \end{array} \right). \] Результатом фильтраци четверки пикселов \(\left\{p_{n,m} \right\}_{n=2i,m=2j}^{2i+1,2j+1}\) будет \[ c_{i,j}^{++}=\frac{1}{2} (p_{2i,2j}+p_{2i+1,2j}+p_{2i,2j+1}+p_{i+1,2j+1} ),\\ c_{i,j}^{\pm}=\frac{1}{2} (p_{2i,2j}+p_{2i+1,2j}-p_{2i,2j+1}-p_{i+1,2j+1} ),\\ c_{i,j}^{+-}=\frac{1}{2} (p_{2i,2j}-p_{2i+1,2j}+p_{2i,2j+1}-p_{i+1,2j+1} ),\\ c_{i,j}^{--}=\frac{1}{2} (p_{2i,2j}-p_{2i+1,2j}-p_{2i,2j+1}+p_{i+1,2j+1} ). \] Полученные множества называются частотными доменами,\(C^{++}=\left\{c_{i,j}^{++} \right\}\) - низко-частотный домен, \(C^{\pm}=\left\{c_{i,j}^{\pm} \right\}\) та \(C^{+-}=\left\{c_{i,j}^{+-} \right\}\) средне-частотные домены и \(C^{--}=\left\{c_{i,j}^{--} \right\}\)–высоко-частотный домен.
Из них наиболее информативным является низко-частотный домен, а наименее - высоко-частотный.
Применим данное преобразование (проведение декомпозиции) к тестовому изображению Lena и упорядочим следующим образом \[ \left[ \begin{array}{ll} C^{++} &C^{+-}\\ C^{\pm} &C^{--} \\ \end{array} \right]. \] После нормализации результат фильтрации будет выглядеть следующим образом Учитывая тот факт, что фильтры Хаара создают полную систему, можно выписать реконструкцию данных \[ p_{2i,2j}=\frac{1}{2}(c_{i,j}^{++}+c_{i,j}^\pm+c_{i,j}^{+ -}+c_{i,j}^{--}),\\ p_{2i,2j+1}=\frac{1}{2}(c_{i,j}^{++}-c_{i,j}^\pm+c_{i,j}^{+ -}-c_{i,j}^{--}),\\ p_{2i+1,2j}=\frac{1}{2}(c_{i,j}^{++}+c_{i,j}^\pm-c_{i,j}^{+ -}-c_{i,j}^{--}),\\ p_{2i+1,2j+1}=\frac{1}{2}(c_{i,j}^{++}-c_{i,j}^\pm-c_{i,j}^{+ -}+c_{i,j}^{--}). \] Последовательное применение декомпозиции к низко-частотному домену называется каскадной схемой А применение декомпозиции к всему множеству доменов - схемой полного пакета Например, схему полного пакета для двух шагов декомпозиции можно проиллюстрировать следующим образом Для трех шагов Качество фильтрации изображения оценивается, впереди всего, тем, насколько более информативный низко-частотный домен, сравнивая его со средне и высоко-частотными доменами.
Среди всех ортогональных вейвлетов базис Хаара имеет худшие показатели.
Можно сравнить разложение константы методом ДКП (дискретное косинус-преобразования) и по Хаару Для получения вейвлетов более сложной конструкции условия ортогональности недостаточно. Ингрид Добеши предложила добавить условия равенства нулю моментов \[ \int x^k\psi(x)=0, k=0,1,\ldots,n-1, \] что эквивалентно условию \[ \sum_i i^k(-1)^{i}h_{i}=0, k=0,1,\ldots,n-1. \] В случае, когда длина носителя равна 4, условия ортогональности будут иметь вид \[ \begin{array}{l} h_0^2+h_1^2+h_2^2+h_3^2=1, \\ h_0h_2+h_1h_3=0. \end{array} \] Добавляем два условия равенства моментов нулю \[ \begin{array}{l} h_0-h_1+h_2-h_3=0, \\ -h_1+2h_2-3h_3=0. \end{array} \] Решая полученную систему уравнений, получим коэффициенты ортогонального фильтра Добеши порядка 4: \[ h_0=\frac{1+\sqrt{3}}{4\sqrt{2}}, h_1=\frac{3+\sqrt{3}}{4\sqrt{2}}, h_2=\frac{3-\sqrt{3}}{4\sqrt{2}}, h_3=\frac{1-\sqrt{3}}{4\sqrt{2}}.\] В матричном виде процесс декомпозиции можно записать так \[ \left[ \begin{array}{c} c^{-1}_{0}\\c^{-1}_{1}\\c^{-1}_{2}\\c^{-1}_{3}\\ d^{-1}_{0}\\d^{-1}_{1}\\d^{-1}_{2}\\d^{-1}_{3} \end{array} \right] = \left[ \begin{array}{cccccccc} h_0&h_1&h_2&h_3\\ &&h_0&h_1&h_2&h_3\\ &&&&h_0&h_1&h_2&h_3\\ h_2&h_3&&&&&h_0&h_1\\ h_3&-h_2&h_1&-h_0\\ &&h_3&-h_2&h_1&-h_0\\ &&&&h_3&-h_2&h_1&-h_0\\ h_1&-h_0&&&&&h_3&-h_2 \end{array} \right] \left[ \begin{array}{c} c^{0}_{0}\\c^{0}_{1}\\c^{0}_{2}\\c^{0}_{3}\\ c^{0}_{4}\\c^{0}_{5}\\c^{0}_{6}\\c^{0}_{7} \end{array} \right].\] А реконструкция будет выглядеть следующим образом \[ \left[ \begin{array}{c} c^{0}_{0}\\c^{0}_{1}\\c^{0}_{2}\\c^{0}_{3}\\ c^{0}_{4}\\c^{0}_{5}\\c^{0}_{6}\\c^{0}_{7} \end{array} \right] = \left[ \begin{array}{cccccccc} h_0&&&h_2&h_3&&&h_1\\ h_1&&&h_3&-h_2&&&-h_0\\ h_2&h_0&&&h_1&h_3&&\\ h_3&h_1&&&-h_0&-h_2&&\\ &h_2&h_0&&&h_1&h_3&\\ &h_3&h_1&&&-h_0&-h_2&\\ &&h_2&h_0&&&h_1&h_3\\ &&h_3&h_1&&&-h_0&-h_2 \end{array} \right] \left[ \begin{array}{c} c^{-1}_{0}\\c^{-1}_{1}\\c^{-1}_{2}\\c^{-1}_{3}\\ d^{-1}_{0}\\d^{-1}_{1}\\d^{-1}_{2}\\d^{-1}_{3} \end{array} \right] .\] Для n = 6 соответствующий фильтр имеет вид \[ h_0=\frac{1}{16\sqrt{2}}\left(1+\sqrt{10}+\sqrt{5+2\sqrt{10}}\right), \] \[ h_1=\frac{1}{16\sqrt{2}}\left(5+\sqrt{10}+3\sqrt{5+2\sqrt{10}}\right), \] \[ h_2=\frac{1}{16\sqrt{2}}\left(10-2\sqrt{10}+2\sqrt{5+2\sqrt{10}}\right), \] \[ h_3=\frac{1}{16\sqrt{2}}\left(10-2\sqrt{10}-\sqrt{5+2\sqrt{10}}\right), \] \[ h_4=\frac{1}{16\sqrt{2}}\left(5+\sqrt{10}-3\sqrt{5+2\sqrt{10}}\right), \] \[ h_5=\frac{1}{16\sqrt{2}}\left(1+\sqrt{10}-\sqrt{5+2\sqrt{10}}\right), \] Использование вейвлет Добеши к фильтрации изображения сводится к последовательному применению фильтров сначала к каждой строке изображения, тем самым проводя деление на низко и высокочастотный домены, а потом применим те же фильтры к каждому из полученных доменов уже по столбцам. Такая последовательность обработки называется схемой Малла (Mallat). Ортогональные фильтры Добеши порядка 4 используются в методе сжатия DJVU.
Но самыми распространенными вейвлет-фильтрами являются биортогональные, которые используются в широком спектре разных задач, так или иначе связанных с частотным анализом сигналов. Относительно нашей темы, то надо отметить, что метод JPEG2000 основан на биортогональных вейвлетах.
Биортогональные преобразования являются обратимыми, но не обязательно ортогональным, что дает возможность получить дополнительные преимущества, так сделать вейвлет-фильтры симметричными, в отличие от ортогональных, которые сделать симметричными невозможно.
Пара масштабируемых функций \[ \varphi^{(0)}(x)=\sqrt{2}\sum_{n\in Z}h^{(0)}_n\varphi^{(0)}(2x-n), \] \[ \varphi^{(1)}(x)=\sqrt{2}\sum_{n\in Z}h^{(1)}_n\varphi^{(1)}(2x-n), \] вместе с соответствующей парой вейвлетов \[ \psi^{(0)}(x)=\sqrt{2}\sum_{n\in Z}g^{(0)}_n\varphi^{(0)}(2x-n), \] \[ \psi^{(1)}(x)=\sqrt{2}\sum_{n\in Z}g^{(1)}_n\varphi^{(1)}(2x-n), \] где \(g_n^{(m)}=(-1)^{1-n} h_{1-n}^{(m)},m=0,1,\) называется биортогональной, если \[ \left\lt \varphi^{(m)}\varphi^{(m)}(\circ -2k)\right\gt=\delta_{0,k},(m,n=0,1;n\ne m;k=0,\pm 1,…), \] тогда имеет место реконструкция сигнала \[ c^{j+1}_k=\sum_{n\in Z}c^j_nh^{(1)}_{k-2n}+\sum_{n\in Z}d^j_ng^{(0)}_{k-2n} \] где формулы декомпозиции такие \[ c^j_k=\sum_{n\in Z}c^{j+1}_nh^{(0)}_{k-2n}, d^j_k=\sum_{n\in Z}c^{j+1}_ng^{(1)}_{k-2n}. \] Как в ортогональном случае наиболее простыми являются фильтры Хаара, так для биортогональных, самый простой вариант дают фильтры 1-3, которые имеют следующий вид \[ h_0^{(0)}=1, h_0^{(1)}=1, h_{\pm 1}^{(1)}=\frac{1}{2}. \] В таком случае декомпозиция сигнала выглядит следующим образом \[ с_k^j=с_{2k}^{j+1}, d_k^j=-\frac{1}{2}c_{2k}^{j+1}+с_{2k+1}^{j+1}-\frac{1}{2} с_{2k+2}^{j+1}, \] а реконструкция \[ с_{2k}^{j+1}=с_k^j, с_{2k+1}^{j+1}=d_k^j+\frac{1}{2} с_k^j+\frac{1}{2} с_{k+1}^j. \] Приведем расчет классической биортогональной пары 3-5 Коэна-Добеши- Фово.
Пусть n = 1 и выпишем условие точности на константе, тогда \[ h_0^{(0)}=\frac{1}{\sqrt{2}}, h_{\pm 1}^{(0)}=\frac{1}{2\sqrt{2}}, \] а для n=2 \[ h_0^{(1)}+2h_{\pm 2}^{(1)}=\frac{1}{\sqrt{2}}, h_{\pm 1}^{(1)}=\frac{1}{2\sqrt{2}}, \] и условие биортогональности имеет вид \[ \frac{1}{8}+\frac{\sqrt{2}}{2}h^{(1)}_{\pm 2}=0, \frac{3}{4}-\sqrt{2}h^{(1)}_{\pm 2}=1. \] Отсюда сразу получаем биортогональную пару \[ h_0^{(0)}=\frac{\sqrt{2}}{2}, h^{(0)}_{\pm 1}=\frac{\sqrt{2}}{4},\\ h_0^{(1)}=\frac{3\sqrt{2}}{4}, h^{(1)}_{\pm 1}=\frac{\sqrt{2}}{4}, h^{(1)}_{\pm 2}=-\frac{\sqrt{2}}{8}. \] Прямой ход вейвлет-преобразования массива \(F=\{f_i\}\) в точке \(2i\) выглядит следующим образом \[ \hat{f}_{2i}=\frac{\sqrt{2}}{4}(f_{2i-1}+2f_{2i}+f_{2i+1}),\\ \hat{f}_{2i+1}=-\frac{\sqrt{2}}{8}(f_{2i-1}+2f_{2i}-6f_{2i+1}+2f_{2i+2}+f_{2i+3}), \] а обратный \[ f_{2i+1}=\frac{\sqrt{2}}{4}(\hat{f}_{2i}+2\hat{f}_{2i+1}+\hat{f}_{2i+2}),\\ f_{2i}=-\frac{\sqrt{2}}{8}(\hat{f}_{2i-2}+2\hat{f}_{2i-1}-6\hat{f}_{2i}+2\hat{f}_{2i+1}+\hat{f}_{2i+2}). \] Приведем самые популярные биортогональные фильтры 5-7 \[ h_0^{(0)}=\frac{3\sqrt{2}}{2}, h^{(0)}_{\pm 1}=\frac{\sqrt{2}}{4},h^{(0)}_{\pm 2}=-\frac{\sqrt{2}}{20},\\ h_0^{(1)}=\frac{17\sqrt{2}}{28}, h^{(1)}_{\pm 1}=\frac{73\sqrt{2}}{280}, h^{(1)}_{\pm 2}=-\frac{3\sqrt{2}}{56}, h^{(1)}_{\pm 3}=-\frac{3\sqrt{2}}{280}. \] Отметим, что в случае известной природы шума, вейвлет анализ является эффективным инструментом его удаления. В этом случае домены частот, соответствующих шуму, можно просто обратить в ноль, в результате чего изображение будет немного сглажено, но без соответствующих помех, для решения этой задачи вовсе не обязательно делать разложение в полный пакет.

Целочисленное вейвлет-преобразование.

Существует целый ряд задач, когда частотное преобразование сигнала (в нашем случае, изображение) не должно вносить ошибку округления. И это принципиально. В таком случае, все, что было изложено ранее, не подходит.
Использование дискретного тригонометрического преобразования Фурье, так же, как и ортогонального (или биортогональных) вейвлет преобразования, приводит к появлению ошибок округления. Поэтому для целочисленного частотного анализа изображения было предложено использовать целочисленное вейвлет-преобразования. Традиционно для этой цели используется сепарабельном схема или схема Малла. В частности, в стандарте JPEG2000 применяется следующая конструкция.
Пусть \(c^0_i \) входные данные, тогда прямой ход выглядит следующим образом \[ c_{2n}^{1}=c^{0}_{2n}-\left[\frac{c^{0}_{2n-1}+c^{0}_{2n+1}+2}{4}\right], c_{2n+1}^{1}=c^{0}_{2n+1}+\left[\frac{c^{1}_{2n}+c^{1}_{2n+2}}{2}\right], \] а обратный – \[ c_{2n+1}^{0}=c^{1}_{2n+1}-\left[\frac{c^{1}_{2n}+c^{1}_{2n+2}}{2}\right], c_{2n}^{0}=c^{1}_{2n}+\left[\frac{c^{0}_{2n-1}+c^{0}_{2n+1}+2}{4}\right]. \] Множество коэффициентов с четными индексами образует низкочастотный домен (L), а с нечетными - высокочастотный (H). Последовательное применение того же алгоритма позволяет разделить исходные данные на четыре частотных домена, из которых домен LL наиболее чувствителен к ошибке, а домен HH наименее чувствителен. Как правило, частотные преобразования используются не для компонентов равноправного пространства цветов (R, G, B), а для неравноправных, как YCrCb, YIQ и пр.
Сам процесс преобразования цветовых компонентов из равноправных (тип данных byte) к неравноправным (тип данных double или float), уже вносит ошибку округления. В случае, если точность преобразования является принципиальной, все следует сделать в целых числах. Соответствующее преобразование можно провести следующим образом

Y = [(R + 2 G + B) / 4]; U = (R - G); V = (B - G);

и обратный ход

G = (Y - [(U + V) / 4]); R = U + G; B = V + G;

где [.] целая часть числа.
Полученные целочисленные цветовые компоненты YUV не столь эффективны, как YCrCb или YIQ, зато лишены такой проблемы, как ошибки округления.
Заметим, что сепарабельних схема изначально ориентирована на одномерные данные, поэтому использование двумерной вейвлет-схемы для обработки изображений является лучшим. Для фильтрации данных \(p_{i, j} \) была предложена несепарабельна (что непредставимо в виде тензорного произведения одномерных) вейвлет-схема с той же сложностью, то есть с равным числом арифметических операций.
На первом этапе проводится децимация данных с четными индексами \[ s^0_{i,j}=p_{2i,2j}. \] На втором шаге по децимованим значениям строим линейные восстановления в точках с одним нечетным индексом и вычисляем соответствующие частотные коэффициенты \[ s^1_{i,j}=p_{2i+1,2j}-\left[\frac{p_{2i,2j}+p_{2i+2,2j}}{2}\right], s^2_{i,j}=p_{2i,2j+1}-\left[\frac{p_{2i,2j}+p_{2i,2j+2}}{2}\right], \] и, наконец, на третьем этапе получаем частотные коэффициенты с нечетными индексами \[ s^3_{i,j}=p_{2i,2j+1}-\left[\frac{s^0_{i,j}+s^0_{i,j+1}+s^0_{i+1,j}+s^0_{i+1,j+1}}{4}\right]-\left[\frac{s^1_{i,j}+s^1_{i,j+1}+s^2_{i,j}+s^2_{i+1,j}}{4}\right]. \] По аналогии со схемой Малла можно записать \[ HH=\left\{s^3_{i,j}\right\}, HL=\left\{s^2_{i,j}\right\}, LH=\left\{s^1_{i,j}\right\}, LL=\left\{s^0_{i,j}\right\}. \]

---