Хвильки (вейвлети) Хаара в 1D

Хвильки (вейвлети) - це математичне приладдя для ієрархіного розкладення зображень. У цій статті я розповім про найпростішу форму хвильок - базис Хаара в одновимірному просторі.

Щоб отримати базове розуміння як працюють хвильки розглянемо простий приклад. Припустимо, що ми маємо одновимірне зображення з такими значеннями:

$$[9\ 7\ 3\ 5]$$

Ми можемо зберегти це зображення залишивши лише середні значення для кожної двійки пікселів:

$$[8\ 4]$$

Очевидно, що частина інформації загубилась. Щоб отримати початкові чотири пікселі нам також треба зберегати коефіцієнти деталізації. Тут такими коефіцієнтами є $1$ для першої пари і $-1$ для другої пари значень.

Роздільність	Середнє значення	Коеф. деталізації
$4$	$[\ 9\ 7\ 3\ 5\ ]$
$2$	$[\ 8\ 4\ ]$	$[\ 1\ {-1}\ ]$
$1$	$[\ 6\ ]$	$[\ 2\ ]$

Отже, для одновимірного базиса Хаара хвильковие перетворення нашого зображення виглядає так: $$[\ 6\ 2\ 1\ -1\ ]$$ Зауважте, що інформація не була ані втрачена, ані набута. Початкове зображення мало 4 коефіцієнти і перетворене також. Початкове зображення можна легко відновити з його перетворення.

Базис хвилькового перетворення

Розглянемо дві ортогональні функції:

1. Маштабувальна - $\psi(x)$.
2. Хвилькова (материнська) - $\phi(x)$.

Які задаються так: $$\begin{equation} \phi(s) := \begin{cases} 1 & 0 \le x < 1\\ 0 & \mbox{інакше.} \end{cases} \end{equation}$$ $$\begin{equation} \psi(s) := \begin{cases} 1 & 0 \le x < \frac{1}{2}\\ -1 & \frac{1}{2} \le x < 1\\ 0 & \mbox{інакше.} \end{cases} \end{equation}$$ Фактично, $[\ 8\ 4\ 1\ -1\ ]$ кодує $[\ 9\ 7\ 3\ 5\ ]$ у такому вигляді: $$\begin{equation} 8\phi(x) + \psi(x) + 4\phi(x-1) - \psi(x-1).\label{dec9735pre} \end{equation}$$ УВАГА! Зараз ми зробимо невеличкий виверт. Ми розглядатимемо зображення як кусково-неперервні функції на проміжку $[0, 1)$. Наприклад, одновимірне зображення $[\ 9\ 7\ 3\ 5\ ]$ виглядатиме так:

Отже, тепер ми вже не можемо сказати, що $\eqref{dec9735pre}$ предсталяє $[\ 9\ 7\ 3\ 5\ ]$. Тепер нам потрібно відмаштабувати наші функції. У загальному випдаку для зображення розміру $2^n$ на першому кроці треба використовувати $\phi_n^i(x) = \phi(2^nx-i)$ і $\psi_n^i(x) = \psi(2^nx-i)$.

Базис Хаара має чудову властивість ортогональності, тобто всі базисні функції - $\phi_0^0, \psi_0^0, \phi_1^1, \psi_1^1, \dots$ ортогональні. Нам також хотілось би нормалізувати базис.

Означення

Базисна функція $u(x)$ нормалізована, якщо $$\langle u|u\rangle = \int_0^1 u(x)u(x)dx = 1.$$

Базис Хаара можна нормалізувати підправивши наявне визначення базисних функцій так: $$\begin{align} \phi_i^j(x):=2^{j/2}\phi(2^jx-i),\\ \psi_i^j(x):=2^{j/2}\psi(2^jx-i). \end{align}$$ В такому базисі ненормалізовані коефіцієнти $[\ 6\ 2\ 1\ -1\ ]$ стають нормалізованими $$[\quad 6\quad 2\quad \frac{1}{\sqrt{2}}\quad -\frac{1}{\sqrt{2}}\quad ].$$ Ніж очислювати ненормалізовані коейіцієнти і потім нормалізовувати їх, можна включити нормалізацію в алгоритм.

Реалізація

Наступний код створює масив 512 символів завдовжки і заповнює його значеннями функції Веєрштраса на проміжку $[-1, 1]$:

#include <algorithm>
#include <cassert>
#include <vector>
#include <fstream>

double Weierstrass(double x) {
    const int summandsCount = 100;
    const double a = 0.7;
    const double b = 11;
    double aRaised = 1;
    double bRaised = 1;
    double res = 0;
    for (int i = 0; i < summandsCount; ++i) {
        aRaised *= a;
        bRaised *= b;
        res += aRaised * cos(bRaised * atan(1) * 4 * x);
    }
    return res;
}

const auto kLevelsCount = 9;
const auto kN = 1 << 9;
const auto kSqrt2 = sqrt(2);

void save(const std::string& name, const std::vector<double>& data) {
    std::ofstream out(name);
    for (int i = 0; i < (int)data.size(); ++i) out << data[i] << std::endl;
}

void DecompositionStep(std::vector<double>& data, std::vector<double> storage, int size) {
    for (int i = 0; i < size; i = i + 2) {
        storage[i / 2] = (data[i] + data[i + 1]) / kSqrt2;
        storage[size / 2 + i / 2] = (data[i] - data[i + 1]) / kSqrt2;
    }
    std::copy_n(storage.cbegin(), size, data.begin());
}

void Decomposition(std::vector<double>& data, std::vector<double>& storage) {
    assert(data.size() == storage.size() && "At the end lod levels take 1 position and differences (data.size() - 1) positions");
    std::transform(data.cbegin(), data.cend(), data.begin(), [s = data.size()](double y) { return y / sqrt(s); });
    auto currSize = (int)data.size();
    while (currSize > 1) {
        DecompositionStep(data, storage, currSize);
        currSize /= 2;
    }
}

int main()
{
    std::vector<double> data(kN);
    const auto interval = std::make_pair(-0.5, 0.5);
    const auto length = interval.second - interval.first;
    for (int i = 0; i < kN; ++i) {
        double x = interval.first + length * i / (kN - 1);
        data[i] = Weierstrass(x);
    }
    
    std::vector<double> storage(data.size());
    Decomposition(data, storage);
    save("compressed.dat", data);
    return 0;
}

А цей код MATLAB'у зчитує утворений файл і відновлює початкові дані починаючи з найгрубішого рівня. По мірі відновлення він записує зображення кожного рівня в файл анімації.

h = figure;
axis tight manual
filename = 'haarWeierstrass.gif';
hold on

xStart = -1;
xEnd = 1;

load 'C:\Yola\VC\Test17\Test17\compressed.dat' compressed -ascii
L = length(compressed);
lodLevelsCount = log(L)/log(2) + 1;
storage = zeros(size(compressed));
for i = 1:lodLevelsCount
    lodStart = 1;
    detailStart = 2^(i-1) + 1;
    localY = compressed(1:detailStart - 1);
    x = linspace(xStart, xEnd, length(localY));
    plot(x, localY);
    title(['Рівень деталізації ' num2str(i-1) ', ' num2str(2^(i-1)) ' точок'])
    drawnow
    frame = getframe(h);
    im = frame2im(frame); 
    [imind,cm] = rgb2ind(im,256); 
    if i == 1 
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf); 
    else 
        imwrite(imind,cm,filename,'gif','WriteMode','append'); 
    end 
    
    if (i < lodLevelsCount)
        % відновлення наступного рівня
        m = 2^((i-1)/2);
        for j = 1:detailStart - 1
            storage(2*j - 1) = compressed(j) + m * compressed(detailStart - 1 + j);
            storage(2*j) = compressed(j) - m * compressed(detailStart - 1 + j);
        end
        compressed(1:2*(detailStart-1)) = storage(1:2*(detailStart-1));
        pause(1);
        clf();
    end
end