2011-12-22

Filtr Lanczos

Filtr Lanczos jest zbudowany na bazie funkcji sinc, która jest odpowiedzią impulsową idealnego filtru dolnoprzepustowego. Dowolny sygnał możemy odfiltrować takim filtrem poprzez splot.

Idealny filtr dolnoprzepustowy jest nieciągły w przestrzeni częstotliwości i rozciąga się nieskończenie w czasie. W grafice komputerowej, w filtrowaniu próbek obrazu podczas jego skalowania lub resamplingu musimy się posłużyć przyciętą w przestrzeni wersją sinc. Dla likwidacji oscylacji w jego widmie biorących się z nieciągłości powstałych w wyniku jego przycięcia musi zostać wymnożony z oknem. Jako okna używa się tego stworzonego przez niejakiego pana Lanczos - stąd nazwa filtru. Dlaczego takim oknem, a nie innym nie wiem. Okazało się także w trakcie eksperymentów, że najlepsze wyniki daje okno o wielkości 2 i 3.

Jedno o czym musimy pamiętać, że tworzenie filtrów to w pewnej mierze połączenie sztuki z nauką i dodatkowo to, że każdy ma inne gusta i inną percepcję. Filtr Lanczos jest prawie idealnym filtrem dolnoprzepustowym. Jednak wcale nie znaczy to najlepszym. Mówiąc o filtrowaniu obrazu chodzi nam o odcięcie wysokich częstotliwości, ale należy pamiętać, że oko nie postrzega widma (w przestrzeni, ale także samej fali elektromagnetycznej), poza tym receptory w oku nie są rozłożone równomiernie by eliminować możliwość aliasingu, sam mózg dokłada swoje w usuwaniu z postrzeganego obrazu niechcianych zakłóceń. Każdy ma inaczej rozłożone czopki w oku i inaczej utkane neurony i inny bagaż doświadczeń. Czopki u każdego reagują na tą samą częstotliwość fali trochę inaczej. Każdy ma inną optykę oka, inne zmętnienie soczewki. I pewnie można by tak wymieniać.

Oba filtry o szerokościach 2 i 3 mają dla pewnym wartości ujemne wartości, stąd istnieje możliwość, że wynikowa wartość splotu będzie ujemna. Czyli należy przewidzieć możliwość radzenia sobie z ujemnymi wartościami koloru. Przy samplowaniu Jitter dla niektórych metod resamplowania możemy otrzymać wartość koloru większą niż jeden.

Niezależnie od tego jak wygląda widmo filtru, jeśli filtr ma ujemne wybrzuszenia oznacza to, że będzie on miał właściwości wyostrzające.

Funkcja okna ma postać:

$\displaystyle W(x) = \begin{cases} \frac{sin \displaystyle\frac{x \pi}{\tau}}{\displaystyle\frac{x \pi}{\tau}} & \text{dla} \;\;\; -\tau \le x \le \tau \\ \\\\\\\\\\\\\\ \qquad 0 & \text{dla} \;\qquad |x|>\tau \end{cases}$

, gdzie $\tau$ to szerokość okna.

Po złożeniu z funkcją sinc otrzymujemy filtr:

$F(x) = \displaystyle \begin{cases} \frac{sin x \cdot sin \displaystyle\frac{x \pi}{\tau}}{x \cdot \displaystyle\frac{x \pi}{\tau}} & \text{dla } \;\;\; -\tau \le x \le \tau \\ \\\\\\\\\\\\\\\\\\\\ \qquad\qquad 0 & \text{dla } \;\qquad |x|>\tau \end{cases}$

Dla zera $\lim_{x \rightarrow 0}{\frac{sinx}{x}}\rightarrow 1$.

Dla $\tau=3$ wykresy filtru, okna, samej funkcji sinc mają postać:


Pole powierzchni filtra to 1, czyli dokonując z nim splotu nie zmieniamy energii sygnału.

Porównanie widma funkcji przyciętej funkcji sinc bez i z oknem Lanczos:


Filtr dwuwymiarowy $L(x,y)=L(x)*L(y)$ ma wygląd:


Jak widać nie powstaje on przez obrót filtru 1D. Objętość pod powierzchnią filtra 2D to jeden.

Przykładowy kod w Matlabie użytego do wygenerowania wykresów:

function [] = plot_lanczos()

    x = linspace(-4, 4, 1000);

    [y, ysinc, ywin] = lanczos(x, 3);

    trapz(x, y)

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    plot(x, ywin, 'r o');
    plot(x, y, 'b.');
    plot(x, ysinc, 'g');
    title('lanczos - spatial');
    legend('window', 'lanczos', 'sinc');

    x = linspace(-3, 3, 100);
    [y, ysinc, ~] = lanczos(x, 3);

    Y = fft(y, 1000);
    YSINC = fft(ysinc, 1000);

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    title('lanczos - freq');
    plot(fftshift(abs(Y)), 'r');
    plot(fftshift(abs(YSINC)));

end

function [y, ysinc, ywin] = lanczos(tau)

    y = zeros(1, size(x, 2));
    ysinc = zeros(1, size(x, 2));
    ywin = zeros(1, size(x, 2));

    for i=1:size(y, 2)
        xx = abs(x(i));
        if (xx == 0)
            ysinc(i) = 1;
            ywin(i) = 1;
        else
            xx = xx * pi;
            xxx = xx / tau;
            ysinc(i) = sin(xx) / xx;
            if (xx/pi <= tau)
                ywin(i) = sin(xxx) / xxx;
            end  
        end

        y(i) = ysinc(i) * ywin(i);
    end
end

function [] = plot_lanczos_3d()

    N = 51;
    tau = 3;

    x = linspace(-tau, tau, N);
    y = linspace(-tau, tau, N);

    lx = linspace(-tau, tau, N);
    ly = linspace(-tau, tau, N);

    for i=1:N 
        lx(i) = lanczos(lx(i), tau);
        ly(i) = lanczos(ly(i), tau);
    end

    [lx, ly] = meshgrid(lx, ly);

    lz = lx .* ly;

    surface(x, y, lz);

    trapz(y, trapz(x, lz))

    function [y] = lanczos(x, tau)

        xx = abs(x);
        if (xx == 0)
            y = 1;
        else
            xx = xx * pi;
            xxx = xx / tau;
            ysinc = sin(xx) / xx;
            ywin = 0;
            if (xx/pi <= tau)
                ywin = sin(xxx) / xxx;
            end  
            y = ysinc * ywin;
        end
    end
end

Przykładowy kod klasy w C# implementujące filtr:

public class LanczosFilter : Filter
{
    public double Tau = 3;

    public override double Ray
    {
        get
        {
            return Tau;
        }
    }

    public override double Evaluate(double a_value)
    {
        double v = Math.Abs(a_value);

        if (v > Tau)
            return 0;

        if (v < Constants.DOUBLE_PRECISION)
            return 1;

        v = v * Constants.PI;
        double p3 = v / Tau;
        return Math.Sin(v) * Math.Sin(p3) / (v * p3);
    }
}

Resampling

To kombinacja upsamplingu i downsamplingu. Pozwala nam ona zmienić próbkowanie sygnału o niecałkowitą skale. Ponieważ jest to kombinacja upsamplingu i downsamplingu dla pewnych współczynników skalowania może się okazać, że współczynnik upsamplingu może być bardzo duży. Takie przypadki w moim przykładzie staram się odrzucać. No i dla współczynników skalowania całkowitych albo ich odwrotności musimy tylko skorzystać z upsamplingu lub downsamplingu.

Funkcja:
function [ tout, fout ] = resampling( fin, fsin, fsout )

m = lcm(fsin, fsout);
u = m / fsin;
d = m / fsout;

if (u > 1000)
    error('up scale factor too big');
end

if (u > 1)
    [tout, fout] = upsampling(fin, fsin, fsin*u);
    if (d > 1)
        [tout, fout] = downsampling(fout, fsin*u, fsout);
    end
elseif (d > 1)
    [tout, fout] = downsampling(fin, fsin, fsout);
else
    tout = (1/fsin) * linspace(0, size(fin, 2), size(fin, 2)); 
    fout = fin;
end

end

Ograniczenie na błąd powinniśmy sobie sami ustawić, w zależności jakiej długości sygnały przetwarzamy i ile mamy pamięci, ... czasu.

Upsampling

Upsampling czyli zwiększenie częstotliwości próbkowania sygnału o zadany całkowity współczynik przy zachowaniu jego widma.

Funkcja:
function [ tout, fout ] = upsampling( fin, fsin, fsout )

scale = fsout / fsin;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

fin = fin * scale;
fout = upsample(fin, scale);
tout =  linspace(0, (1 / fsout) * (size(fout, 2) - 1), size(fout, 2)); 

[~, k] = sinc_gen(fsout, fsin/2);
if (mod(size(k, 2), 2) == 0)
    [~, k] = sinc_gen(fsout, fsin/2, fix(size(k, 2) / 2) * 2 + 1);
end

fout = conv(fout, k, 'same');

fout = fout(1:1:end-scale+1);
tout = tout(1:1:end-scale+1);

end
Funkcja upsample uzupełnia fin zerami. Zera w ilości scale-1 są wstawiane między próbki i na koniec. My te końcowe zera usuwamy. Dzięki czemu ostatnia próbka sygnału wyjściowego odpowiada ostatniej próbce sygnału wejściowego.

Okno sinc musi mieć nieparzystą długość by uniknącć przesunięcia wyniku.

Sam splot to tzw. idealna rekonstrukcja, tutaj o tyle prostsza i szybsza, że współczynnik skalowania jest całkowity i sinc rekonstruujący może być wyliczony raz i używamy dla każdego punktu z osobna.

W stosunku do wbudowanej w Matlab funkcji interp błąd przez nas popełniany jest na tyle niewielki, że możemy powiedzieć, że wynika on z innej metody interpolacji.

Przykład:

%% input data

fsin = 6000;
fsout = 30000;
scale = fsout / fsin;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = upsampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

fout1 = interp(fin, scale);
fout1 = fout1(1:1:end-scale+1);

%% compare them

hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'interp');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'interp');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);



Downsampling

Downsampling czyli zmniejszenie częstotliwości próbkowania sygnału o zadany całkowity współczynnik przy zachowaniu jego widma (poniżej częstotliwości do której zmniejszamy). Ponieważ tutaj rekonstrukcja odbywa się z częstotliwością mniejszą o całkowitą wartość w przeciwieństwie do rekonstrukcji wystarczy nam tylko jeden impuls sinc. Cały kod jest znacznie szybszy.

Funkcja:

function [ fout, tout ] = downsampling( fin, fsin, fsout )

scale = fsin / fsout;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

% works only with odd window.
k = sinc_gen(fsin, fsout/2);
if (mod(size(k, 2), 2) == 0)
    k = sinc_gen(fsin, fsout/2, size(k, 2) + 1);
end

tout = linspace(0, (1/fsin) * (size(fin, 2)-1), ... 
                ceil(size(fin, 2)/scale));
fout = zeros(1, size(tout, 2));

for m=1:1:size(fout, 2)
    fout(m) = conv_point(fin, k, (m - 1) * scale + 1);
end

end

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 6000;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2.4; % czym dluzszy tym lepsze widmo

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% downsampling
[fout, tout] = downsampling(fin, fsin, fsout);

fout1 = decimate(fin, scale);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'decimate');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'decimate');

%% error;

error_calc(fout, forg);
error_calc(fout1, forg);




Błąd mojej wersji jest porównywalny do tej wbudowanej w Matlab. Ilość próbek wejściowych jest równa wielokrotnością skalowania. Jak widać obie wersje rozmijają się, Matlab na początku, ja na końcu. Dla innych liczb próbek obie funkcję potrafią zwrócić bardzo porównywalne wyniki.

Idealna rekonstrukcja sygnału

Mamy sygnał fin, który został spróbkowany z częstotliwością fsin. Chcemy wygenerować sygnał, który został spróbkowany z częstotliwością fsout. Jeśli fsout > fsin w sygnale wyjściowym nie powinno być żadnej częstotliwości powyżej fsin. Jeśli fsout < fsin w sygnale wyjściowym powinny się znaleźć tylko częstotliwości mniejsze od fsout z sygnału wejściowego.

Wszelkie interpolacje (spline, bicubic, itp) w dziedzinie czasu odpadają, gdyż wprowadzimy tym częstotliwości do sygnału, których w nim nie było oryginalnie. Jedne sprawują się gorzej, inne lepiej, ale generalnie interpolacje odtworzą sygnał wizualnie niby dobrze, częstotliwościowo nie będzie to to samo. Ma to szczególne znaczenie kiedy pracujemy z sygnałem dźwiękowym. No dobra za wyjątkiem podanej niżej interpolacji. Rekonstrukcję możemy też przeprowadzić z dziedzinie częstotliwości.

Aby odtworzyć oryginalny sygnał potrzeba nam minimum dwie próbki na na okres, w praktyce powinno być ich więcej. Możemy próbkować dokładnie w miejscach gdzie sinusoida przechodzi przez zero. A nawet jeśli jesteśmy minimalnie minięci podczas próbkowania z zerem, pobierane przez nas wartości są małe i błąd przez to jest duży.

Teoria mówi, że sygnał oryginalny może zostać dokładnie zrekonstruowany ze swych próbek, jeśli próbki były pobierane w odstępach czasu T, zaś sygnał próbkowany nie zawiera częstotliwości większych niż połowa częstotliwości próbkowania. Czyli jeśli próbkujemy z częstotliwością 44100Hz, to sygnał nie powinien zawierać częstotliwości większych niż 22050Hz. W przeciwnym razie częstotliwości wyższe od limitu pojawią się w sygnale zrekonstruowanym zniekształcając go (aliasign). Twierdzenie Whittaker–Shannon'a mówi:

$x(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot {\rm sinc}\left(\frac{t - nT}{T}\right)\,$

Inaczej sygnał oryginalny to splot impulsów z próbek z funkcją sinc o odpowiednim charakterze.

Teraz mając odtworzony sygnał ciągły, próbkujemy go ponownie z nową częstotliwością. Jeśli z zmniejszą to najpierw pozbywamy się z niego częstotliwości większych od połowy nowej częstotliwości próbkowania.

Całość robi poniższy kod. Jeśli mamy zrekonstruować sygnał z niższą częstotliwością to filtrowane są próbki wejściowe. Taki sygnał podlega rekonstrukcji zgodnie z wyżej podaną formułą tylko w interesujących nas czasach.

Kod funkcji:

function [ fout, tout ] = reconstruction( fin, fsin, fsout )

Tout = 1 / fsout;
Tin = 1 / fsin;

if (fsin == fsout)
    fout = fin;
    tout = Tout* linspace(0, size(fout, 2), size(fout, 2)); 
    return;
end

scale = fsout / fsin;

if (scale < 1)
    % for compatibility with downsampling and decimate
    fout = zeros(1, ceil(size(fin, 2)* scale));
else
    fout = zeros(1, fix(size(fin, 2) * scale));
end

tout = linspace(0, (1/fsin) * (size(fin, 2) - 1), size(fout, 2));            

% when downsampling cut frequencies below fsout/2 from fin 
% function.
if (fsin > fsout)
    k = sinc_gen(fsin, fsout/2);
    fin = conv(fin, k, 'same');
end

W = fix(160 / (fsout*2 / fsin)); % experimental
W = min(W, size(fin, 2));

tin = (0:1:size(fin, 2)-1) * Tin;

if (W > size(fin, 2))
    W = size(fin, 2);
end

for m=1:size(fout, 2)
    [~, min_index] = min(abs(tout(m) - tin));
    min_index = fix(min_index - W/2);
    if (min_index < 1)
        min_index = 1;
    end
    max_index = min_index + W - 1;
    if (max_index > size(fin, 2))
        max_index = size(fin, 2);
    end
    sc = sinc((tout(m) - (min_index-1:1:max_index-1)*Tin) / Tin);
    fout(m) = sum(fin(min_index:1:max_index) .* sc);
end

% for compatibility with upsampling
if (scale > 1)
    fout = fout(1:1:end-scale+1);
    tout = tout(1:1:end-scale+1);
end

Dla idealnej rekonstrukcji okno W powinno być równe długości funkcji, w praktyce może być znacznie krótsze. Jeśli mamy do czynienia z downsamplingiem to najpierw sygnał trzeba przefiltrować filtrem dolnoprzepustowym. Inaczej częstotliwości powyżej fsout, które były w fin odbiją się i pojawią się jako zakłócenia w sygnale. Z tak zakłóconym sygnałem już niż nie możemy zrobić - nie da się stwierdzić co jest oryginałem, a co zakłóceniem.

Rekonstrukcja taka jest dokładna, ale wyjątkowo czasochłonna, dla każdej próbki wyjściowej trzeba policzyć osobny impuls sinc.

Normalnie okno sinc powinno być równe wektorowi wejściowemu, tutaj jest krótsze, o długości powiedzmy eksperymentalnej dla której nie widać różnicy.

Jest też naniesionych parę poprawek na kompatybilność z innymi funkcjami dotyczącymi długości wektora wyjściowego.

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 7300;
Tin = 1/fsin;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

%build input func.
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = reconstruction(fin, fsin, fsout);

%resample
[fout1, tout1] = resampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(0:Tin/20:tend, func(0:Tin/20:tend), 'g');
plot(tout, fout, 'b');
plot(tout1, fout1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% compare fq

FORG = fftshift(abs(fft(forg)));
FIN = fftshift(abs(fft(func(0:Tin:tend))));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqout1 = fsout/2*linspace(-1, 1, size(FOUT1, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

s = size(FIN, 2) / size(FOUT, 2);

hold on;
plot(fqout, s*FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, s*FOUT, 'b');
plot(fqout1, s*FOUT1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);



Filtr FIR o minimalnym przesunięciu fazowym

Weźmy wszystkie filtry o takiej samej charakterystyce częstotliwościowej G, ale różniących się charakterystyką przesunięcia fazowe. Filtr o minimalnym przesunięciu fazowym to taki filtr, że w dowolnym okresie czasu $<0, t>$ energia sygnału wyjściowego jest większa lub równa dowolnemu innemu filtrowi z G. Czyli sygnał w takim filtrze najszybciej narasta i najszybciej maleje. Sygnał jest maksymalnie przesunięty ku zeru.

Dla filtrów FIR Matlab posiada funkcję do znajdowania MPF (minimum-phase filter) rceps. Niestety nie działa ona za dobrze. Rezultat nie wygląda tak jak powinien, a na dodatek dla niektórych długości wektora policzenie wyniku jest niemożliwe.

Moja wersja:


function [ minblip ] = fir_min_phase( impulse )

impulse1 = [impulse, zeros(1, size(impulse, 2)*7)];

y = ifft(log_safe(abs(fft(impulse1))));

%y(1) = y(1) * 1;
y(2:fix(end/2)) = y(2:fix(end/2)) * 2;
%y(fix(end/2+1) = y(fix(end/2+1) * 1;
y(fix(end/2)+2:end) = 0;

minblip = real(ifft(exp(fft(y))));

minblip = minblip(1:1:size(impulse, 2));

minblip = minblip / trapz(minblip);

end

function [ r ] = log_safe( v )

r = log(v);

for k=1:1:size(r, 2)
if (r(k) == -inf)
r(k) = -1e-100;
end
end

end


Kompletna magia dla mnie - matematycznie nie wiem jak to działa. Zmniejszając początkową rozdzielczość okna, albo źle dobierając okno +/-1 możemy otrzymać wersję Matlabową.

Porównanie impulsu sinc i jego wersji MPF.


fs = 7093790/2;
d1 = firwin(fs, 21000, 2047);
d2 = fir_min_phase(d1);
[~, d3] = rceps(d1);

hold on;
plot(d1, 'r');
plot(d2, 'g');
plot(d3, 'b');
legend('org', 'my', 'matlab');

D1 = abs(fftshift(fft(d1)));
D2 = abs(fftshift(fft(d2)));
D3 = abs(fftshift(fft(d3)));
fx = linspace(-fs/2, fs/2, 2047);

hold on;
plot(fx, D1, '*');
plot(fx, D2, 'g');
plot(fx, D3, 'o');
legend('org', 'my', 'matlab');




Widzimy, że we wszystkich 3 przypadkach charakterystyka częstotliwościowa filtrów jest taka sama.

2011-12-19

Idealny filtr dolnoprzepustowy

Idealny filtr dolnoprzepustowy powinien po stronie częstotliwości być linią prostą o wartości dodatniej od f=0 do częstotliwości odcięcia i powyżej nie przyjmować wartość zero. Ponieważ w transformacie Fouriera odwrotnej chcąc po stronie czasu uzyskać wartości rzeczywiste prostokąt ten powinien być symetryczny względem f=0.

Najpierw zdefiniujmy sobie tzw. funkcję prostokątną:

$\text{rect}(t) = \begin{cases}0 & \text{if } \; |t| \;> \;^1/_2 \\ ^1/_2 & \text{if }\; |t| \;= \;^1/_2 \\1 & \text{if }\; |t| \;<\; ^1/_2 \\\end{cases}$

W naszym wypadku przyjmując za częstotliwość odcięcia B filtr po stronie częstotliwości powinien mieć postać $\displaystyle rect(\frac{f}{2B})$.

Jego transformatą po stronie czasu jest:

$\displaystyle f(t)=\mathcal{F}^{-1}\{ rect(\frac{f}{2B})\} = 2Bsinc(2Bt)$

Funkcja generująca taki filtr:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

  T = 1 / sampling_rate; 

  if (nargin == 2)
    samples = 160 / (cutfq*2 / sampling_rate); % experimental
  end

  tout = T * (1:samples); 
  tout = tout - T*(samples+1)/2;
  B = 2*cutfq;
  k = sinc(B*tout);
  k = k / trapz(k);

  if (nargout == 2)
      t = tout;
  end

end
Jeśli nie podamy ilości próbek jest ona ustalana na najlepszą (według mnie) wartość, która przy filtrowaniu nie powoduje błędów.

Zmienna sampling_rate to częstotliwość próbkowania sygnału. sampling_rate/2 to największa częstotliwość w sygnale. Czyli np: jeśli sampling_rate=44100, to $f_{max}=22050$. Na samym końcu impuls jest normalizowany, tak by jego pole powierzchni było równe jeden, dzięki czemu splot z sygnałem zachowuje jego energię.

Zauważmy, że dla takiej samej wartości samples jeśli stosunek sampling_rate do cutoff_fq FIR będzie taki sam to odpowiedź impulsowa będzie taka sama (w sensie wartości w tablicy).

Przykładowy filtr sinc:

[~, w3] = sinc_gen(44100, 4000, 200);
tx = linspace(0, 200/44100, 200);
plot(tx, w3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');



Filtr jest symetryczny i zmierza do zera jak czaszmierza do nieskończoności. Dla parzystej liczby sampling_rate dwie środkowe próbki mają takie same wartości. Dla nieparzystej największą wartość ma środkowa próbka (ten jest bardziej pożądany, gdyż nie przesuwa sygnału podczas splotu).

Wraz ze zwiększaniem sampling rate i utrzymywaniem wysokiej cutoff_fq kształt filtru będzie coraz bardziej przypominać impuls.

W dziedzinie częstotliwości filtr taki wygląda tak:

f3 = abs(fftshift(fft(w3)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');



Oscylacje na brzegach są spowodowane tym, że nas filtr jest gwałtownie się urywa na brzegach. Aby tego uniknąć na filtry nakłada się okno, jest ich wiele rodzajów, i nie zamierzam wnikać, który się do czego się nadaje.

w4 = w3 .* blackman(200)';
f3 = abs(fftshift(fft(w4)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');



Moja wyspecjalizowana funkcja generująca sinca:
function [ k, t ] = firwin( fs, fcut, len, prevent_shift )

k = sinc_gen(fs, fcut);

if (sum(size(len) - [1 1]) == 0)
    L = len;
else
    L = size(len, 2);
end

if (size(k, 2) > L)
    k = sinc_gen(fs, fcut, L);
end;

p = true;
if (nargin == 4)
    p = prevent_shift;
end

% prevent time shifting
if (p)
    if (mod(size(k, 2), 2) == 0)
        k = sinc_gen(fs, fcut, size(k, 2) - 1);
    end
end

% apply window
w = kaiser(size(k, 2), 9);
k = k .* w';

% normalize
k = k / trapz(k);

if (nargout == 2)
    t = linspace(0, (size(k, 2)-1)/fs, size(k, 2));
end

end

Korzysta ona z ustawionego na stałe okna (Kaiser, 9). Można także zapobiegać generacji sinca o parzystej liczbie próbek, splot z takim sincem w czasie przesuwa nam rezultat o jeden.

Porównanie do funkcji wbudowanej w Matlaba:

spectrum_cut = 50000;
Fs = 350000;
len = 100;

k1 = fir1(len-1, 2 * spectrum_cut / Fs, kaiser(len, 9));
k1 = k1 / trapz(k1);

k2 = firwin(Fs, spectrum_cut, len, false);

hold on;
plot(k1, '*');
plot(k2, 'r');

err = error_calc(k1, k2);



Funkcja pomocnicza:

function [ r ] = error_calc( v1, v2 )

res = sum(abs((v1 - v2) / v2) / size(v1, 2));

if (nargout == 0)
    display(sprintf('error: %f', res));
else
    r = res;
end

end

Resampler LocalGrid

Budujemy dwuwymiarową (może też być jednowymiarowa) siatkę wag filtra. Środek filtru jest pozycji (0,0). Stąd dla filtru o promieniu 3 siatka musi mieć wymiary 7x7. Odległości z których pobieramy wartości filtra to -3, -2, -1, 0, 1, 2, 3. Zakładając, że wartości brzegowe w naszej siatce to zera możemy wymiary siatki zmniejszyć o 2.

Taką siatkę centrujemy w punkcie padania promienia i dla każdego piksela, na który nakłada się siatka dodajemy wartość koloru promienia ważoną wartością z siatki.

W przypadku supersamlingu wymiary siatki mnożymy przez pierwiastek z ilości promieni przewidywanych na piksel.

Bardzo ważne jest by uwzględnić wartości zwracane przez Sampler, tzn. zamiast ... 13.5, 14.5, 15.5, ... możemy dostać ... 13.49, 14.5, 15.51, ... (oczywiście z dużo mniejszym błędem). Tego rodzaju błędów nie da się uniknąć - zaokrąglanie liczb binarnych do postaci dziesiętnej.

W samym Resamplerze dodałem właściwość Subreslution, która dodatkowo zwiększa wymiary siatki. Jej zadaniem jest zwiększenie dokładności resamplowania, tak, by metoda ta w miarę zwiększania właściwości coraz bardziej przypominała resampler Exact, a coraz mniej Resizer. By kod mógł współdziałać zarówno z filtrem Box o promieniu 0.5 jak i z innymi filtrami o promieniach całkowitych Subreslution musi być nieparzyste. Inaczej wagi albo jednego albo drugiego filtru będą oscylować na granicy promienia w zależności od błędów numerycznych Samplera.

W analizie kodu najważniejsze jest dobre zrozumienie wszystkich $\pm\frac{1}{2}$, $\pm 1$, Floor, Ceiling, Round.

Kod:

public struct PixelInterpolator
{
    private ColorFloat m_numerator;
    private double m_denumerator;

    /// <summary>
    /// Note: a_ray_color in SRGB space.
    /// </summary>
    /// <param name="a_ray_color"></param>
    /// <param name="a_filter_weight"></param>
    public void Collect(ColorFloat a_ray_color, double a_filter_weight)
    {
        m_numerator += a_ray_color * a_filter_weight;
        m_denumerator += a_filter_weight;
    }

    /// <summary>
    /// Note: color in SRGB space.
    /// </summary>
    public ColorFloat Color
    {
        get
        {
            if (m_denumerator == 0)
                return ColorFloat.Black;

            return m_numerator / m_denumerator;
        }
    }

    public override string ToString()
    {
        return String.Format("num: {0}, denum: {1}", m_numerator, m_denumerator);
    }
}

public abstract class InterpolateResampler : FilterResampler
{
    protected PixelInterpolator[,] m_pixels;
    protected OverlayCorrector m_overlay_corrector;

    protected override ColorFloat ResamplePixel(int a_x, int a_y)
    {
        return Gamma.SRGBToLinear(m_pixels[a_x, a_y].Color);
    }

    internal override void RenderStart(RenderStartPhase a_phase)
    {
        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            m_pixels = new PixelInterpolator[Film.Width, Film.Height];
            m_overlay_corrector = OverlayCorrector.Create(
                OverlayMethod.Mirror, Film.Width, Film.Height, 
                Filter.Ray.Ceiling());
        }
    }

}

public class LocalGridResampler : InterpolateResampler
{
    private double[,] m_weights;

    [YAXNode]
    public int Subresolution = 3;

    private double m_scale;

    internal override void RenderStart(RenderStartPhase a_phase)
    {
        base.RenderStart(a_phase);

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            var sampler = Film.Sampler as UniformSampler;

            if (sampler == null)
                throw new InvalidOperationException();
                    
            if (Subresolution % 2 == 0)
                throw new InvalidOperationException();

            m_scale = Subresolution * sampler.Subresolution;
            PrepareWeights();
        }
    }

    private void PrepareWeights()
    {
        List<List<PixelWeight>> horz_weights, vert_weights;
        ResizerResampler.PrecalculateWeights((int)m_scale, (int)m_scale, 1, 1, 
            out horz_weights, out vert_weights, Filter, false);
        double[] weights = horz_weights[0].Select(pw => pw.Weight).ToArray();

        Debug.Assert(weights.Length % 2 == 1);

        m_weights = new double[weights.Length, weights.Length];

        for (int x = 0; x < weights.Length; x++)
        {
            double wx = weights[x];

            for (int y = 0; y < weights.Length; y++)
            {
                double wy = weights[y];
                m_weights[x, y] = wx * wy;
            }
        }
    }

    private double GetWeight(double a_dx, double a_dy)
    {
        int x = (a_dx * m_scale).Round();
        x += (m_weights.GetLength(0) - 1) / 2;

        if (x < 0)
            return 0;
        else if (x >= m_weights.GetLength(0))
            return 0;

        int y = (a_dy * m_scale).Round();
        y += (m_weights.GetLength(1) - 1) / 2;

        if (y < 0)
            return 0;
        else if (y >= m_weights.GetLength(1))
            return 0;

        return m_weights[x, y];
    }

    internal override void SetRayColor(Vector2 a_pixel, ColorFloat a_color)
    {
        a_color = Gamma.LinearToSRGB(a_color);

        int start_x = (int)Math.Floor(a_pixel.X - Filter.Ray);
        int end_x = (int)Math.Ceiling(a_pixel.X + Filter.Ray) - 1;

        int start_y = (int)Math.Floor(a_pixel.Y - Filter.Ray);
        int end_y = (int)Math.Ceiling(a_pixel.Y + Filter.Ray) - 1;

        for (int y = start_y; y <= end_y; y++)
        {
            int cy = m_overlay_corrector.CorrectY(y);

            for (int x = start_x; x <= end_x; x++)
            {
                int cx = m_overlay_corrector.CorrectX(x);
                double w = GetWeight((x + 0.5) - a_pixel.X, (y + 0.5) - a_pixel.Y);
                m_pixels[cx, cy].Collect(a_color, w);
            }
        }

        #if DEBUG
        for (int y = start_y - 1; y <= end_y + 1; y++)
        {
            Debug.Assert(GetWeight((start_x + 0.5 - 1) - a_pixel.X, (y + 0.5) - 
                a_pixel.Y).IsAlmostEquals(0));
            Debug.Assert(GetWeight((end_x + 0.5 + 1) - a_pixel.X, (y + 0.5) - 
                a_pixel.Y).IsAlmostEquals(0));
        }
        for (int x = start_x; x <= end_x; x++)
        {
            Debug.Assert(GetWeight((x + 0.5) - a_pixel.X, (start_y + 0.5 - 1) - 
                a_pixel.Y).IsAlmostEquals(0));
            Debug.Assert(GetWeight((x + 0.5) - a_pixel.X, (end_y + 0.5 + 1) - 
                a_pixel.Y).IsAlmostEquals(0));
        }
        #endif
    }

    public override ResamplerType ResamplerType
    {
        get
        {
            return ResamplerType.LocalGrid;
        }
    }
}

Resampler Exact

Dla każdego promienia skanujemy piksele wokół (uwzględniamy promień filtra). Dla każdego piksela wagę liczymy dokładnie z filtra. Nie korzystamy z żadnych prekalkulowanych tabelek. Metoda ta powinna zapewniać największą dokładność.

W analizie kodu najważniejsze jest dobre zrozumienie wszystkich $\pm\frac{1}{2}$, $\pm 1$, Floor, Ceiling, Round.

Kod:
public struct ColorInterpolator
{
    private ColorFloat m_numerator;
    private double m_denumerator;

    public void Collect(ColorFloat a_ray_color, double a_filter_weight)
    {
        m_numerator += a_ray_color * a_filter_weight;
        m_denumerator += a_filter_weight;
    }

    public ColorFloat Color
    {
        get
        {
            if (m_denumerator == 0)
                return ColorFloat.Black;

            return m_numerator / m_denumerator;
        }
    }
}

public abstract class Resampler
{
    public Filter Filter = new BoxFilter();
    internal Film Film;

    public static Resampler Create(ResamplerType a_type)
    {
        if (a_type == ResamplerType.Exact)
            return new ExactResampler();
        else if (a_type == ResamplerType.LocalGrid)
            return new LocalGridResampler();
        else if (a_type == ResamplerType.Resizer)
            return new ResizerResampler();
        else
            throw new NotImplementedException();
    }

    public Rectangle GetAffectedRect(Rectangle a_src_rect)
    {
        int border = (int)Math.Ceiling(Filter.Ray) + 2;
        a_src_rect.Inflate(border, border);
        a_src_rect.Intersect(new Rectangle(0, 0, Film.Width, Film.Height));
        return a_src_rect;
    }

    public void Resample(Rectangle a_rect, ColorArrayBase a_dest)
    {
        for (int y = a_rect.Top; y < a_rect.Bottom; y++)
        {
            for (int x = a_rect.Left; x < a_rect.Right; x++)
            {
                a_dest.SetColor(x, y, ResamplePixel(x, y));
            }
        }
    }

    public abstract void RenderStart(RenderStartPhase a_phase);
    protected abstract ColorFloat ResamplePixel(int a_x, int a_y);
    public abstract void SetRayColor(Vector2 a_pixel, ColorFloat a_color);
}

public class ExactResampler : InterpolateResampler
{
    public override void SetRayColor(Vector2 a_pixel, ColorFloat a_color)
    {
        int start_x = (int)Math.Floor(a_pixel.X - Filter.Ray);
        int end_x = (int)Math.Ceiling(a_pixel.X + Filter.Ray) - 1;

        int start_y = (int)Math.Floor(a_pixel.Y - Filter.Ray);
        int end_y = (int)Math.Ceiling(a_pixel.Y + Filter.Ray) - 1;

        for (int y = start_y; y <= end_y; y++)
        {
            int cy = m_overlay_corrector.CorrectY(y);
            double wy = Filter.Evaluate((y + 0.5) - a_pixel.Y);

            for (int x = start_x; x <= end_x; x++)
            {
                int cx = m_overlay_corrector.CorrectX(x);
                double wx = Filter.Evaluate((x + 0.5) - a_pixel.X);

                m_rays[cx, cy].Collect(a_color, wy * wx);
            }
        }

        Debug.Assert(Filter.Evaluate(
            a_pixel.X - (start_x + 0.5 - 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.X - (end_x + 0.5 + 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.Y - (start_y + 0.5 - 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.Y - (end_y + 0.5 + 1)).IsAlmostEquals(0));
    }
}

Resamplowanie Exact i LocalGrid - porównanie

Resampling LocalGrid jest nieznacznie szybszy. Dla filtru Lanczos który wymaga policzenia dwóch funkcji sinus i jest chyba najkosztowniejszym obliczeniowo filtrem jest to 10%.

Jeśli chodzi o jakość to już Resampler LocalGrid z minimalną Subresolution równą 2 jest prawie identyczny jak Exact. Dalsze zwiększanie Subresolution nie poprawia jakości generowanego obrazu.

2011-12-18

Splot

Poniżej podane są funkcje w Matlabie realizujące splot, działające jak ich odpowiedniki w Matlabie:
function [ fout ] = conv_full( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  fout = zeros(1,s1+s2-1);
  
  for m=1:size(fout, 2)
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

function [ fout ] = conv_same( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  fout = zeros(1, H-L+1);
  
  for m=L:H
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m-L+1) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

Funkcja conv_same przydaje się w przetwarzaniu sygnałów, wtedy z reguły u traktujemy jako sygnał zaś v jako odwrotną transformatę Fouriera jakiegoś filtru, którego chcemy użyć. Splotowi u i v w czasie odpowiada mnożenie po stronie częstotliwości. Wynik conv_same będzie takiego samego rozmiaru jak u.

Funkcja v może być także odpowiedzią impulsową funkcji przejścia H(s) dowolnego układu. Znając H(s) liczymy jaka jest odpowiedź układu na impuls (1 po stronie s), czyli jest to transformata odwrotna 1*H(s). Mówimy wtedy o odpowiedzi impulsowej filtru. Po stronie czasu mamy układ o funkcji przejścia H(s), podajemy na jego wejście impuls (po stronie częstotliwości widmo impulsu jest ciągłe od minus nieskończoności do plus nieskończoności i o stałej amplitudzie). To co otrzymujemy na wyjściu to odpowiedź w czasie przefiltrowanych wszystkich możliwych częstotliwości. Transformatą Fouriera tej odpowiedzi wyznacza nam charakterystykę częstotliwością filtru.

I jeszcze przykład funkcji realizującej splot dla danego punktu:
function [ r ] = conv_point( fin, k, point )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  if (point < 1)
    error('bad point index');
  end
  
  if (point > size(fin, 2))
    error('bad point index');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  m = point + L - 1;
  
  fin1 = max(m-s2+1, 1);
  fin2 = min(m, s1);
  k1 = max(m-s1+1, 1);
  k2 = min(m, s2);
  r = sum(fin(fin1:fin2) .* k(k2:-1:k1));

end

Porównanie wersji same i full:

fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
fs = 24000; % czestotliwosc probkowania

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

fcut = 1000; % czestotliwosc odciecia
tend = 3/fsin1; % dlugosc funkcji
sinc_samples = 200;

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

if (mod(size(tx, 2), 2) == 1)
tx = tx(1:1:end-1);
end

% wejscie
fin = func_in(tx)/15;

% filtr
w1 = sinc_gen(fs, fcut, sinc_samples);

fout1 = conv(fin, w1, 'same');
fout2 = conv(fin, w1, 'full');

figure(1);
hold on;
plot(fin, 'r');
plot(w1, 'g');
plot(fout1, 'b');
plot(fout2, 'c');
legend('wejscie', 'filtr', 'conv full', 'conv same');



Przykład:

fs = 30000; % czestotliwosc probkowania
fcut = 1000; % czestotliwosc odciecia
fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
tend = 1/fsin1*3.35; % dlugosc funkcji
sinc_samples = 1000;

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

% funkcja wyjsciowa - idealne odciecie.
func_out = @(ti) sin(2*pi*fsin1*ti+fsin1_phase);

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

% oblicz wartosci obu funkcji w czasie.
fin = func_in(tx);
fout = func_out(tx);

%[~, fin] = sinc_gen(fs, (fsin1+fsin2)*1.12, size(tx, 2));
%[~, fout] = sinc_gen(fs, fsin1, size(tx, 2));

% przygotuj sinca
[~, w1] = sinc_gen(fs, fcut, sinc_samples);

% oba okna w czestotliwosci.
fw1 = fft(w1);
fwx = linspace(-fs/2, fs/2, size(fw1, 2));

% funkcja wejsciowa w czestotliwosci.
ffin = fftshift(abs(fft(fin)));
ffinx = linspace(-fs/2, fs/2, size(ffin, 2));

% przyciecie poprzez splot
fout2 = conv(fin, w1, 'same');

figure(1);
hold on;
plot(fwx, fftshift(abs(fw1)), 'r');
legend('sinc');
title('filtr dolnoprzepustowy w przestrzeni czestotliwosci');

figure(2);
hold on;
plot(fin, 'c');
plot(fout);
plot(fout2, 'r');
title('porowanie w czasie');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

figure(3);
hold on;
plot(ffinx, fftshift(abs(fft(fin))), 'c');
plot(ffinx, fftshift(abs(fft(fout))));
plot(ffinx, fftshift(abs(fft(fout2))), 'r');
title('porowanie w czestotliwosci');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');





Po odkomentowaniu linii zastępujących nasze funkcje sincami, wyniki przedstawiają się tak:



Test naszych funkcji:

% przetestujmy poprawnosc dzialania naszej funkcji conv_cut.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = conv_same(a, a);
y12 = conv_same(a, b);
y13 = conv_same(a, c);
y14 = conv_same(a, d);

y21 = conv_same(b, a);
y22 = conv_same(b, b);
y23 = conv_same(b, c);
y24 = conv_same(b, d);

y31 = conv_same(c, a);
y32 = conv_same(c, b);
y33 = conv_same(c, c);
y34 = conv_same(c, d);

y41 = conv_same(d, a);
y42 = conv_same(d, b);
y43 = conv_same(d, c);
y44 = conv_same(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

% przetestujmy poprawnosc dzialania naszej funkcji conv_full.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'full');
x12 = conv(a, b, 'full');
x13 = conv(a, c, 'full');
x14 = conv(a, d, 'full');

x21 = conv(b, a, 'full');
x22 = conv(b, b, 'full');
x23 = conv(b, c, 'full');
x24 = conv(b, d, 'full');

x31 = conv(c, a, 'full');
x32 = conv(c, b, 'full');
x33 = conv(c, c, 'full');
x34 = conv(c, d, 'full');

x41 = conv(d, a, 'full');
x42 = conv(d, b, 'full');
x43 = conv(d, c, 'full');
x44 = conv(d, d, 'full');

y11 = conv_full(a, a);
y12 = conv_full(a, b);
y13 = conv_full(a, c);
y14 = conv_full(a, d);

y21 = conv_full(b, a);
y22 = conv_full(b, b);
y23 = conv_full(b, c);
y24 = conv_full(b, d);

y31 = conv_full(c, a);
y32 = conv_full(c, b);
y33 = conv_full(c, c);
y34 = conv_full(c, d);

y41 = conv_full(d, a);
y42 = conv_full(d, b);
y43 = conv_full(d, c);
y44 = conv_full(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

%% przetestujmy poprawnosc dzialania naszej funkcji conv_point.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = zeros(1, size(a, 2));
y12 = zeros(1, size(a, 2));
y13 = zeros(1, size(a, 2));
y14 = zeros(1, size(a, 2));
for i=1:1:size(a, 2)
y11(i) = conv_point(a, a, i);
y12(i) = conv_point(a, b, i);
y13(i) = conv_point(a, c, i);
y14(i) = conv_point(a, d, i);
end

y21 = zeros(1, size(b, 2));
y22 = zeros(1, size(b, 2));
y23 = zeros(1, size(b, 2));
y24 = zeros(1, size(b, 2));
for i=1:1:size(b, 2)
y21(i) = conv_point(b, a, i);
y22(i) = conv_point(b, b, i);
y23(i) = conv_point(b, c, i);
y24(i) = conv_point(b, d, i);
end

y31 = zeros(1, size(c, 2));
y32 = zeros(1, size(c, 2));
y33 = zeros(1, size(c, 2));
y34 = zeros(1, size(c, 2));
for i=1:1:size(c, 2)
y31(i) = conv_point(c, a, i);
y32(i) = conv_point(c, b, i);
y33(i) = conv_point(c, c, i);
y34(i) = conv_point(c, d, i);
end

y41 = zeros(1, size(d, 2));
y42 = zeros(1, size(d, 2));
y43 = zeros(1, size(d, 2));
y44 = zeros(1, size(d, 2));
for i=1:1:size(d, 2)
y41(i) = conv_point(d, a, i);
y42(i) = conv_point(d, b, i);
y43(i) = conv_point(d, c, i);
y44(i) = conv_point(d, d, i);
end

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);


Funkcja pomocnicza:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

T = 1 / sampling_rate;

if (nargin == 2)
samples = 160 / (cutfq*2 / sampling_rate); % experimental
end

tout = T * (1:samples);
tout = tout - T*(samples+1)/2;
B = 2*cutfq;
k = sinc(B*tout);
k = k / trapz(k);

if (nargout == 2)
t = tout;
end

end

Impuls, Splot, rodzaje filtrów, odpowiedź impulsowa

W poprzednich postach podałem przykład funkcji sinc i jej transformacji. W miarę jak będziemy brać coraz większy zakres czasu prostokąt po stronie częstotliwości będzie mieć coraz mniejsze oscylacje w punkcie skoku. W miarę jak nasza funkcja sinc będzie coraz węższa nasz prostokąt będzie coraz szerszy. Nasza funkcja sinc będzie coraz bardziej przypominać impuls. Innymi słowy transformatą impulsu jest nieskończone widmo częstotliwości o takiej samej amplitudzie. Impuls jest paczką sinusoid o wszystkich częstotliwościach i takiej samej amplitudzie.

Splot funkcji definiujemy jako:

$f(t)*h(t) = \int_{-\infty}^{\infty} f(\tau)h(t-\tau)d\tau$

Istnieje także jego dyskretny odpowiednik. W transformacie Fouriera splot ma bardzo ważną cechę. Mianowicie mnożeniu funkcji po stronie częstotliwości odpowiada splot po stronie czasu i na odwrót.

Co się stanie jeśli na jakiś filtr, czyli obiekt który odpowiednio filtruje częstotliwości sygnały wejściowego podamy impuls. Impuls to nieskończone widmo częstotliwości. Odpowiedz którą otrzymamy to suma odfiltrowanych sinusoid. Taka odpowiedź daje nam jednoznaczną informację o charakterystyce częstotliwościowej filtru.

W zależności od cech odpowiedzi możemy różnie klasyfikować filtry. Ze względu na energię sygnału wyjściowego: aktywne (dodają energii) i bierne (nie dodają energii, inaczej pasywne). Ze względu na czas po jakim następuje odpowiedź na sygnał wejściowy, czy jest on zawsze taki sam, czy od czegoś zależy np. od sygnału wcześniejszego. W szczególności mówimy, że filtr jest niezmienny w czasie (time-invariant) jeśli na każdą pojedyńczą sinusoidę odpowiada on w taki sam sposób i po takim samym czasie. Inną cechą filtrów jest ich linearność, bądź i nie. Filtr jest linearny jeśli dowolne rozłożenie sygnału wejściowego, w szczególności na pojedyncze sinusoidy i ich przefiltrowanie, a następnie złożenie sygnału wyjściowego równe jest odpowiedzi filtru na złożony sygnał wejściowy. Inaczej filtr filtruje każdą częstotliwość niezależnie od innych. Innym podziałem jest FIR (finite impulse response) i IIR (infinite impulse response). Odpowiedź filtru FIR zanika w skończonym czasie, zaś odpowiedź filtru IIR nie - asymptotycznie zmierza do zera z reguły oscylując. Filtr IIR możemy przyciąć do FIR zwłaszcza jeśli posługujemy się jego dyskretną postacią.

Wszystkie filtry o których będziemy mówić są pasywne, liniowe, time-invariant, FIR.

Robiąc z filtru IIR filtr FIR wprowadzamy błąd. Dobierając odpowiednio długość FIR możemy te zniekształcenia przenieść poniżej rozdzielczości sygnału. Lub powyżej częstotliwości poniżej których nie powinno ich być (jak w sygnale audio). Część błędu, czyli oscylacje na brzegach możemy zmniejszyć splatając odpowiedzieć impulsową z oknem. Funkcji okien jest bardzo wiele - po stronie częstotliwości mają one kształt okna, najprostszym jest RECT.

Jak wcześniej powiedzieliśmy mnożeniu po stronie częstotliwości odpowiada splot po stronie czasu. W tym przypadku splot funkcji i odpowiedzi impulsowej. Jest to bardzo ważna cecha pozwalająca nam bez użycia transformat filtrować sygnał.

Przesunięcie w czasie, a FFT

Przesunięcie w czasie nie wpływa na absolutne wartości amplitudy w FFT. Zmienia się tylko faza składowych sinusoid. Dla symetrycznego wykresu względem jego środka po stronie czasu część urojona amplitudy jest zerowa. Połówka symetrycznego wykresu będzie miała takie same absolutne amplitudy po stronie częstotliwości jak cały symetryczny wykres, różnica będzie w przesunięciach fazowych. Niesymetryczność po stronie FFT amplitud lub faz, powoduje, że po stronie czasu nasz wykres nie jest
rzeczywisty. Z uwagi na błędy zaokrąglenia rozsądnie jest jeśli spodziewamy się symetryczności dla bezpieczeństwa brać wartości bezwzględne
przy konwersji z czasu do częstotliwości, a przy konwersji odwrotnej część rzeczywistą wyniku.

Wpływ przesunięcia w czasie na FFT:
function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

  T = 1 / sampling_rate; 

  if (nargin == 2)
    samples = 160 / (cutfq*2 / sampling_rate); % experimental
  end

  tout = T * (1:samples); 
  tout = tout - T*(samples+1)/2;
  B = 2*cutfq;
  k = sinc(B*tout);
  k = k / trapz(k);

  if (nargout == 2)
    t = tout;
  end

end

spectrum_cut = 20000;
fs = 3.5e6;

[k, t] = firwin(fs, spectrum_cut, 2047);

k1 = [zeros(1, 700), k];

K = fftshift(fft(k, 2047));
K1 = fftshift(fft(k1, 2747));
fx = linspace(-fs/2, fs/2, 2047);
fx1 = linspace(-fs/2, fs/2, 2747);

figure(1);
hold on;
plot(fx, abs(K), '.');
plot(fx1, abs(K1), 'r');
legend('org', 'shifted');

figure(2);
hold on;
plot(fx, unwrap(angle(K)), '.');
plot(fx1, unwrap(angle(K1)), 'r');
legend('org', 'shifted');



Niesymetryczność w czasie to niezerowa część urojona po stronie FFT:

figure(3);
hold on;
plot(fx, imag(K), '.');
plot(fx1, imag(K1), 'r');
legend('org', 'shifted');


2011-12-12

Transformata Fouriera

Istnieje wiele transformat, które generalnie można powiedzieć przekształcają funkcję z jednej przestrzeni w inną. Z reguły istnieje także transformata odwrotna. Nie wszystkie funkcje da się także poddać takiej transformacji. Transformata Fouriera to takie przekształcenie, które możemy zinterpretować jako przekształcenie funkcji w dziedzinie czasie w funkcję w dziedzinie częstotliwości.

$f(\omega) = \int_{-\infty}^{\infty} f(x)\ e^{- 2\pi i x t}dt$

Transformata odwrotna:

$f(t) = \int_{-\infty}^{\infty} f(\omega)\ e^{2 \pi i x \omega}d\omega$

W wyniku transformacji funkcji o wartościach rzeczywistych zmiennej rzeczywistej transformata Fouriera jest symetryczna względem zera. Tylko taka symetria po stronie częstotliwości gwarantuje nam otrzymanie funkcji o wartościach rzeczywistych podczas transformacji odwrotnej. Ma to szczególne znaczenie podczas dyskretnej transformaty Fouriera kiedy niedokładności w symetrii są częste.

Po stronie częstotliwości widmo sygnału interpretujemy jako sumę sinusoid o danych częstotliwościach, przesunięciu fazowym i amplitudzie. Ponieważ we wzorach mamy całkę, a nie sumę, musiała by to być suma nieskończona. Przechodząc z całek na sumy przechodzimy na dyskretną transformatę Fouriera. Wynikiem transformaty jest funkcja liczby zespolonej - jej moduł to amplituda dla danej częstotliwości, jej przesunięcie fazowe to przesunięcie fazowe sinusoidy o danej częstotliwości.

Ponieważ z reguły jesteśmy zainteresowaniu amplitudą widma jej wykresy najczęściej podaje się po stronie częstotliwości. Przesunięcia fazowe z reguły pomija się. Zauważmy, że dla różnie poprzesuwanych składowych sinusoid widmo amplitud jest takie same. I drugiej strony przesuniecie funkcji w czasie nie zmienia jej widma. Właściwość ta jest wykorzystywana podczas kompresji obrazów i filmów. Szczególnym przypadkiem dobrania przesunięć fazowych jest funkcja o minimalnym przesunięciu fazowym (jest to taka funkcja w której w dowolnym czasie $<0,t>$ energia sygnału jest maksymalna).

W wyniku tego, że dyskretna transformata Fouriera operuje na dyskretnych wartościach wprowadza ona błędy. Wspomnianą wcześniej asymetrię funkcji rzeczywistej i nierzeczywistą postać jej transformacji odwrotnej. Innym rodzajem błędu jest to, że pojedyncza sinusoida może być reprezentowana poprzez rozmyty pik po stronie częstotliwości. Dlatego po retransformatą warto wziąć jej część rzeczystą (nie moduł).

W dalszej części będziemy mówili tylko o transformacji dyskretnej.

Oto przykład transformaty funkcji sinc:


Teoretycznie jej transformatą powinien być prostokąt. W rzeczywistości ponieważ wartości są dyskretne, a funkcja nie jest nieskończona, transformata nie jest idealnym prostokątem i wykazuje oscylacje na brzegach skoku. Są one tym mniejsze im nasz sinc jest dłuższy.

Dla transformaty odwrotnej: sinc po stronie częstotliwości będzie prostokątem po stronie czasu z takimi samymi oscylacjami na brzegach. Ponieważ nasz sinc jest skończony w czasie brakuje nam widma dla wyższych częstotliwości i one uniemożliwiają nam zreprodukowanie skoku i zlikwidowanie oscylacji - można powiedzieć ich wygaszenie.

Przykład jak może wyglądać transformacja złożenia dwóch sinusoid o f=650 i f=1500. Oczywiście w zależności jak dobierzemy krok próbkowania naszych ciągłych funkcji, ilości próbek czasu i ilości próbek częstotliwości możemy otrzymać wykres bardziej lub mniej przypominający piki.

Warto podkreślić, że dyskretna transformata jest bezwymiarowa podobnie jak wartości wejściowych próbek. Chcąc otrzymać wartości częstotliwości na osi musimy je policzyć.

Skończony sygnał w czasie nieokresowy wymaga do opisu nieskończonego widma częstotliwości. Podobnie jest z nieciągłościami w czasie. I drugą stronę: skończone widmo nie powoli nam na opisanie skończonego w czasie sygnału. Te dwa ograniczenia silnie wiążą się z DFFT, której liczba próbek jest skończona.

Metody resamplowania

Nie chodzi mi tyle o filtry, ale o sposób samego resamplowania.

Sposób pierwszy. Robimy podobnie jak podczas zmiany rozdzielczości obrazka (tak działał resampler po staremu) Tworzymy tablice subpixeli, każdy subpixel przyciągamy do odpowiedniego subpixela. Taką tablicę resamplujemy. Wady to duża ilość pamięci, niedokładność związana z nieuwzględnianiem prawdziwych współrzędnych promienia.

Do metody drugiej i trzeciej musimy użyć innego sposobu pamiętania wartości kolorów promieni.

Ponieważ rezygnujemy z bitmapy subpixeli, która także nawiasem mówiąc byłaby kłopotem podczas resamplownania adaptywnego, które zamierzam dodać, musimy jakoś pamiętać wartości kolorów promieni, które będą potrzebne do zresamplowania próbek które jeszcze nie zostały zrenderowane - z innych prostokątów. Przez prostokąt rozumiem tutaj obszar znajdujący się obok, który nie został jeszcze zrenderowany. W obecnej chwili scena jest dzielona na w miarę równe prostokąty o określonej licznie subpikseli i wszystkie są one renderowane równolegle. Ponieważ w procesie resamplingu wartość subpikseli granicznych może być potrzeba do wyliczenia pikseli w sąsiednim prostokącie musimy je jakoś zapamiętać. Dodatkowo aby uniknąć problemów w wyścigiem zabronione jest renderowanie sąsiednich, wzajemnie na siebie wpływających prostokątów. Możemy na przykład pamiętać je w tablicy hashującej, liście o wielkości równej liczbie pikseli na scenie - każdy element listy byłby null-em, ewentualnie bardzo rzadko tablicą wartości i współrzędnych promieni. Z tym, że podczas zmiany sposobu renderowanie np. z prostokątów na linia po linii mogło by się okazać, że kod wyliczający czy dany piksel będzie jeszcze wykorzystywany mógłby okazać się nie uniwersalny. No i jeśli renderujemy prostokątami to subpiksele brzegowe trzeba pamiętać maksymalnie dla 3 kwadratów - potrzebny byłby jakiś licznik referencji - brzegowe 2, rogowe - 3, albo 2 jeśli są na brzegu sceny. Dosyć kłopotliwe... Dlatego zdecydowałem się na rozwiązanie wymagające więcej pamięci, ale za to bardziej uniwersalne. Resamplowanie działa w 2D. Zmiana rozdzielczości obrazka za pomocą filtru 2D jest wolniejsza niż dwukrotne użycie filtru 1D z pionie i w poziome. Wynik jest ten sam. Kształt filtru powstaje poprzez obrót filtra 1D wokół Y, jeśli tak możemy nieprecyzyjnie powiedzieć. Ogólny wzór na wartość zrekonstruowanego piksela ma postać:

$\displaystyle C(x,y) = \frac{\sum f(x-x_i,y-y_i)c(x_i,y_i) }{\sum f(x-x_i,y-y_i)}$

Wartości wag filtra wyliczone z tego wzory nie wymagają normalizacji, nawet jeśli pole pod filtrem jest różne od 1, a ni tylko jeśli wyliczone wagi obarczone są błędem numerycznym. Normalizacja jest wzór - normalizacja to podzielenie wszystkich wag przez sumę wag, tak by ta suma wynosiła 1.

$C(x,y)$ to wynikowy kolor. $c(x_i,y_i)$ to kolor promienia. $f(x-x_i,y-y_i)$ to wartość filtru. Zgodnie z tym wzorem korzystnie jest pamiętać kolor jako iloraz dwóch liczb, dzięki czemu możemy dodawać do piksela wartości nowych promieni. Znika nam także problem pamiętania wartości wyliczonych promieni na potrzeby sąsiednich prostokątów.

Pomimo normalizacji może dojść do przestrzelenia koloru. Weźmy dwie wagi: 2 i -0.5. Niech oba kolory mają wartość 1, wtedy wartość zrekonstruowana to 1, ale jeśli drugi kolor ma w tym punkcie wartość 0, to otrzymamy 2/1.5. Widać to zwłaszcza podczas Samplowania Jitter 1x1 na Lanczos, który ma przedziały ujemne.

No i teraz dwie pozostałe metody.

Druga to liczenie za każdym razem wagi filtru. W metodzie trzeciej budujemy subtablicę wag o określonej wielkości. Centrujemy ją na pikselu wynikowym. Patrzymy w którą komórkę subtablicy trafił promień i bierzemy z niej wagę. Przy czym subtablica może mieć większą rozdzielczość niż wartość ilości próbek na piksel, dzięki czemu może być dokładniejsza niż metoda pierwsza. Przy odpowiednio dużej wielkości powinna być tak samo efektywna jak metoda druga, a przy tym szybsza.

Metody 2 i 3 w stosunku do 1 generują obraz nieznacznie przesunięty. Przesunięcie jest znacznie mniejsze od 0.5 i ... nie wiem skąd się bierze.

2011-12-11

Zmiany w kodzie

Cały kod obecnie odpowiadający za resamplowanie generuje obrazy nie tak dobrej jakości jakby mógł, gdyż zamiast opierać się na prawdziwej odległości pozycji promienia na płaszczyźnie projekcji od pozycji piksela który chcemy zrekonstruować, obecnie najpierw generowana jest bitmapa subpixeli i ona jest reskalowana. Ta metoda ma też inną wadę, wymaga strasznie dużo pamięci. Tak więc wraz z przepisaniem samplowania (zmianie musi ulec sposób generowania próbek) i resamplowania (sposób rekonstrukcji) zamiast bitmapy subpixeli będziemy tylko pamiętali tyle subpixeli ile potrzebne nam jest do przyszłej rekonstrukcji sąsiednich pixeli.

Niejaki przy okazji pojedyńczy obiekt kamery zostanie zastąpiony kolekcją kamer z jedną aktywną. Każda kamera zostanie połączona z obiektem klasy Film. Klasy Film utworzą swoją własną kolekcję, tak więc wiele kamer będzie mogło korzystać z tego samego filmy. Do klasy Film zostanie zaś przeniesione większość opcji z RenderOptions. W tym wszelkie informacje pozwalające zbudować łańcuch postprocesorów obrazu.

InterestRects zostały przeniesione do Filmu. ... I w trakcie zrezygnowałem z tego. Prawie tego nie używałem. Poza tym dla o ile na początku się starałem by rzeczywiście dobrze definiować obszary zainteresowania to później nie miałem już na to sił - liczba testowych scen to już przeszło tysiąc. Poza tym tak naprawdę to w wyniku zmian w kodzie popsuć się może wszystko i wszędzie nie tylko w obszarach zainteresowań.

SceneElement został zastąpiony przez SceneActor. SceneActor to baza dla wszystkich obiektów mogących występować na scenie - posiadających pozycję. ScemeElement jest od teraz klasą bazową dla SceneActor - posiada nazwę i właściwość Scene. Jest też dolnym ograniczeniem parametru generycznego dla ObjectList.

Klasa Material dziedziczy z SceneElement. Klasa MaterialsList dziedziczy z ObjectsList.

Dodałem do kamery dwie nowe metody UVToPixel i PixelToUV. Przestrzeń UV to zawierający się od 0 do 1 zbiór punktów z których może zostać wypuszczony promień. Dla kamer które wysyłają promienie z prostokątnego obszaru płaszczyzny wartość UV nie różnią niczym w sensie znaczenia od koordynatów tekstury. Co innego kiedy dojdą kamery które wysyłają promienie z wycinka walca (panoramiczne) albo z wycinka sfery (tzw. fish-eye). Procedura Camera.CreateRay wymaga do wygenerowania promienia współrzędnych UV.

Samplery generują próbki zawsze w zakresie 0 do 1. Parametry do generacji do liczba wierszy i kolumn. Ich iloczyn to ilość punktów do wygenerowania. Dla samplerów jednorodnych na każdy prostokąt będzie przypadać jedna próbka, dla samplerów niejednorodnych niekoniecznie.

Postprocesory tworzą łańcuch ustalany w runtime. Parametry dla nich (gamma correction, tone mapping) zniknęły z klasy Film. Łańcuch postprocesorów podłączony jest do klasy Film. Na razie zrezygnowałem z utworzenia osobnej kolekcji łańcuchów postprocesorów. Taki łańcuch postprocesorów jest serializowany do xml bezpośrednio.

Parametrami dla postprocesingu jest jest bitmapa wejściowa i wyjściowa. Z reguły mogą być one takie same. Jeśli tak nie jest to od postprocesora wymaga się by podał jaki ma być jej rozmiar. O ile ma się on zmienić. Wymaga się też od niego by prawidłowo przeliczył na nowy obszar prostokątu. Sama klasa odpowiedzialna za wywołania postprocesorów będzie się starała utrzymać taką bitmapę wyjścia do następnego razu aby nie tworzyć ich za każdym razem. Taka bitmapa wyjściowa będzie mu później wielokrotnie w trakcie generowania poszczególnych obszarów przesyłana mu do uzupełnienia tych obszarów i będzie wykorzystana w dalszym łańcuchu. Jeśli postprocesor nie skaluje obrazu bitmapa wyjściowa będzie równa wejściowej. Posprocesor winien w niej zmodyfikować tylko zadany obszar.

Ponieważ kod resamplujący ulegnie zmianie postanowiłem obecny kod, który jest tak naprawdę kodem na zmniejszanie wielkości obrazka o SampleSize przekształcić w uniwersalny postprocesor pozwalający na skalowanie wyrenderowanej bitmapy.

Problem z mapowaniem

Przy okazji mapowania na torus powierzchni danych wielomianem pojawił się problem, który na razie odłożyłem w przyszłość, ale prędzej czy później trzeba go będzie jakoś rozwiązać. Mianowicie torus dany bezpośrednio mapowany jest w sześcian UVW od 0 do 1. W przypadku powierzchni wyższego rzędu mapowanej na torus nie wiemy nic o geometrii bryły i jest ona mapowana w UVW proporcjonalnie do $r/(R+r)$. Na razie poradziłem sobie z tym rozciągając powierzchnię wielomianową w Y do sześcianu, dzięki czemu dobrze współpracuje z mapperem. A następnie ręcznie ją skalując by przypominała torus. Niestety trzeba to robić ręcznie.

Wyobraźmy sobie teraz siatkę przypominającą spłaszczoną sferę. Chcemy ją zamapować sferycznie. Teraz ponieważ sfera wymaga UVW w pełnym sześcianie musielibyśmy najpierw przeskalować siatkę by zawierała się w sześcianie. A następnie spłaszczyć ją skalowaniem.

Rozwiązanie pierwsze to skalowanie UVW zanim zostanie przekształcone w UV. Z tym, że na razie (nie wiem czy słusznie) rozdzialiłem świat UVW i UV. Nie ma przejścia od Volume do Layer. I raczej nie chciałbym tego robić póki nie będzie to absolutnie potrzebne.

Rozwiązanie drugie to dodanie do mapperów skalowania i translacji, być może nawet macierzy, którą UVW jest transformowane przez mapowaniem.

Problem z nieskończonościami

Użycie nieskończoności (np. jako granice na AABB) powoduje wiele problemów. Niekończoności nie porównują się dobrze ze sobą, wiele standardowych funkcji (Math.Sign) zwraca wyjątek po ich napotkaniu. Zastąpienie ich Double.MaxValue i Double.MinValue jest dużo lepszym pomysłem. Z jednym wyjątkiem liczby takie konwertują się na tekst, ale z powrotem już nie, gdyż otrzymujemy Overflow Exception. Tak więc musimy wprowadzić nasze własne minima i maksima i się ich trzymać: Double.MaxValue/10 i Double.MinValue/10

Powierzchnie dane wielomianem

A dokładniej wielomiany w funkcji trzech zmiennych x, y, z, przeważnie w postaci uwikłanej. Powierzchnie, zwłaszcza te zamknięte, to z reguły wielomiany stopnia 4 lub 6. Przy obliczaniu intersekcji z nimi bardzo ważne jest przybliżenie punktu startu promienia do AABB, by zmniejszyć błędy obliczeniowe, które mogą się pojawić. Dla niektórych powierzchni AABB nie da się określić wprost, gdyż zależy ona od parametrów takiego wielomianu (np. promień dla sfery).

W przypadku niektórych kształtów likwidowałem parametry funkcji, gdyż odpowiadały one tak naprawdę za skalowanie. Czasami też przesuwałem powierzchnię (manipulując x, y, z) by znalazła ona się w środku układu współrzędnych oraz by zawierała się w przedziale od -1 do 1.

Na ten moment sprawdziłem renderowanie prawie wszystkich powierzchni zamkniętych podanych na MathWorld. I wystarczy tego. Cel jakim było w ogóle zrenderowanie takich powierzchni został osiągnięty. Na jego potrzeby musiałem napisać całkiem spory kawałek kodu, który bierze wielomian w postaci tekstu, parsuje do tablicy wyrażeń, buduje z niej RPN, z RPN drzewo. Drzewo wyrażeń można zoptymalizować - wymnażanie przez stała, dodawanie, mnożenie wyrazów, uproszczenia związane z 0 i 1. Oprócz tego wyrażenie można rozszerzyć, czyli wymnożyć, wypotęgować co się da - zbudować jedną wielką sumę wyrazów. Dzielenie przez inne wielomiany, czyli funkcje wymierne nie jest obsługiwane, podobnie jak wykładniki potęgowania nie będące liczbami. Dla zmiennych wykładniki te powinny być dodatnimi liczbami całkowitymi. Brakuje kodu na redukcję wyrażeń, czyli czegoś odwrotnego do rozszerzania. Gradient jest liczony wprost z funkcji bez jej rozszerzania. Dla wyliczenia intersekcji podstawiamy do wielomianu równanie kierunkowe linii, rozszerzamy je i grupujemy po stopniu t. Każda taka grupa to współczynnik dla wielomianu intersekcji. Następnie kompilujemy kod który zwraca nam taki wielomian z już wyliczonymi współczynnikami po podaniu mu jako parametrów kierunku i punktu startu promienia. Kompilacja następuje poprzez zbudowanie LINQ Expression. Kod wynikowy jest bardzo szybki.

Szczegółów kodu nie zamierzam podawać. Ale kilka rad po napisaniu czegoś takiego mam. Drzewo wyrażenia powinno być niezmienne. Totalnie niezmienne drzewo będzie wymagało skopiowania całego drzewa przy dowolnej zmianie, a będziemy to robić bardzo często. Najlepiej kod drzewo umieścić w osobnej przestrzeni nazw i poukrywać wszystkie do manipulowania jego strukturą bezpośrednio (rodzic, dzieci) i na zewnątrz wystawić tylko zbiór ustalonych metod. Dlatego najlepiej metody zmieniające drzewo jakoś zaizolować, a na zewnątrz wystawić tylko metody nie modyfikujące drzewa, ale zwracające nowe zmienione drzewo. Co do elementów drzewa to powinno się ono składać tylko z symboli, dodawania, mnożenia, dzielenia, potęgowania. Unarne minusy i odejmowanie powinniśmy zastąpić dodawaniem i mnożeniem przez -1. Dzielenie także jak najszybciej konwertować na mnożenie (obsługujemy dzielenie tylko przez stałe). Przez symbol rozumiemy wyrażenie A*N^B, gdzie A to stała, N to nazwa symbolu, B to jego potęga. A=0 to zero, B=0 to liczby. Dodatkowo warto od początku założyć, że dodawanie i mnożenie może mieć więcej dzieci, dzięki czemu znajdowanie grup typu ABCD+EFGH będzie szybkie i proste. Dzięki temu, że potęgowanie mamy od razu w symbolu nie musimy przeszukiwać drzewa nadmiernie. A stała w symbolu przez którą mnożymy sprawi, że często nie będziemy musieli tworzyć nowych symboli liczb. Na samym początku redukujemy liczby i przenosimy potęgowanie do symboli tak długo jak pozbędziemy się potęgowania i dzielenia. Pozbywamy się także minusów i odejmowania. Spłaszczamy drzewo tak by dodawania i mnożenie nie posiadało innych dzieci - odpowiednio dodawania i mnożenia i w takiej postaci staramy się utrzymywać drzewo.

Sam kod intersekcji nie różni się zbytnio od kodu intersekcji dla torusa. Z tym, że inaczej jest budowany wielomian intersekcji.

Można oczywiście policzyć gradient i współczynniki wielomianu intersekcji na piechotę. Jednakże jest to bardzo pracochłonne. Poza tym mało uniwersalne - mając cały kod przetwarzający możemy implementować dowolne powierzchnie dane wielomianem.

W moim przypadku największe problemy sprawił mi kształt serca: $(x^2+\frac{9}{4} y^2+z^2-1)^3-x^2z^3-\frac{9}{80}y^2z^3$

Serce bez AABB:


Serce z AABB - start promienia jest przybliżany do AABB i dla takiego punktu startu budowany jest wielomian intersekcji.


Serce z AABB - likwidacja miejsc zerowych.


Sprawdziłem jak wyglądają współczynniki i same miejsca zerowego dla promienia odbitego od błędnych punktów. Okazało się, że następuje ponowna detekcja uderzenia w powierzchnie od której się odbijamy. Możemy powiększyć Constants.MINIMIAL_DISTANS = 1e-8 - miejsce od którego szukamy pierwiastków. Ale w tym przypadku ponieważ pierwszy współczynnik wielomianu intersekcji okazał się mniejszy od dokładności Constants.DOUBLE_PRECISION = 1e-12 postanowiłem napisać kod eliminujący takie współczynniki, przesuwający pozostałe współczynniki i zmiejszający stopień wielomianu. Tak mały współczynnik może odpowiadać za miejsce zerowe w granicach 1e-5. Taki zabieg nigdzie nie wyprodukował dodatkowych błędów. Błędy nie znikły całkowicie, ale są już prawie niewidoczne, a na pewno rozmyją przy supersamplowaniu.

Serce z przesuwaniem punktu startu wzdłuż normalnej.


Serce ciągle ma dwa piksele za ciemne. Nie pozostało mi nic innego jak odsunąć punkt startu od powierzchni wzdłuż normalnej. Wartość przy której znikły ostatnie błędy to Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA = 1e-5. Jest to bardzo dużo.

Przykład powierzchni jabłka z refrakcją:


Ponieważ dysponuję kodem na obiekt torusa liczonego bezpośrednio i teraz mam do dyspozycji torus dany przez wielomian mogę porównać ich jakość i szybkość. Implementacja oparta na wielomianie jest ponad dwa razy wolniejsza. Przede wszystkim największy współczynnik wielomianu, który po uproszczeniu jeden jedynką tutaj jest wyliczany. Prawdopodobnie pozostałe współczynniki też nie są w optymalnej postaci. Coś za coś. Np. dla serca liczba współczynników to 217 po pełnym rozwinięciu. To daje obraz zabawy podczas wyliczania współczynników wielomianu.

Kod obiektu sceny:

public class PolynomialSurface : RenderableObject
{
    protected PolynomialFunction m_polynomial;
    private Func<Vector3, Vector3> m_gradient; 
    private Func<Vector3, Vector3, Polynomial> m_intersection;
    private PolynomialSolver m_solver;

    [YAXNode]
    public string Equation;

    [YAXNode]
    public string Parameter1 = "";

    [YAXNode]
    public string Parameter2 = "";

    [YAXNode]
    public string Parameter3 = "";

    [YAXNode]
    public AABB LocalBoundBox;

    private Vector3 m_aabb_minimum;
    private Vector3 m_aabb_maximum;

    public PolynomialSurface(Vector3 a_right, Vector3 a_up) :
        base(a_right, a_up)
    {
        Name = "Polynomial surface";
        m_aabb_minimum = Vector3.MAXIMUM;
        m_aabb_maximum = Vector3.MINIMUM;
        Closed = false;
        LocalBoundBox = AABB.INFINITY;
        m_uv_mapper = new SphericalUVMapper();

        Update(UpdateFlags.All);
    }

    public override Intersection GetIntersection(
        Intersection a_source_ray_intersection, Ray a_ray)
    {
        Vector3 local_dir;
        Vector3 local_start;
        bool back_hit = false;

        if (!TransformToLocalToBoundBox(a_ray, out local_dir, out local_start))
            return Scene.NoIntersection;

        if ((a_source_ray_intersection != null) && 
            (a_source_ray_intersection.SceneObject == this))
        {
            if (a_ray.RaySurfaceSide == RaySurfaceSide.RefractedSide)
            {
                local_start += ((WorldToLocalNormal * 
                    a_source_ray_intersection.Normal).Normalized) *
                    -Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
            }
            else
            {
                local_start += ((WorldToLocalNormal * 
                    a_source_ray_intersection.Normal).Normalized) *
                    Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
            }
        }

        Polynomial p = m_intersection(local_start, local_dir);

        double dist = m_solver.Solve(p);

        if (dist == Double.PositiveInfinity)
            return Scene.NoIntersection;

        Vector3 local_pos = local_start + local_dir * dist;

        Vector3 normal = m_gradient(local_pos);

        back_hit = normal * local_dir > 0;

        if (back_hit && OneSide)
            return Scene.NoIntersection;

        Vector3 world_pos = LocalToWorld * local_pos;
        double world_dist = (world_pos - a_ray.Start).Length;

        Intersection intersection = new Intersection()
        {
            PrevIntersection = a_source_ray_intersection, 
            SceneObject = this,
            SourceRay = a_ray,
            LocalPos = local_pos,
            Dist = world_dist,
            Scene = Scene,
            BackHit = back_hit,
            Pos = world_pos
        };

        // Numerical errors.
        if (intersection.Normal * intersection.SourceRay.Dir >= 0)
            return Scene.NoIntersection;

            m_aabb_minimum = Vector3.Minimize(local_pos, m_aabb_minimum);
            m_aabb_maximum = Vector3.Maximize(local_pos, m_aabb_maximum);

        return intersection;
    }

    public override Vector3 GetNormal(Intersection a_intersection)
    {
        if (a_intersection.BackHit)
        {
            return (LocalToWorldNormal * 
                -m_gradient(a_intersection.LocalPos)).Normalized;
        }
        else
        {
            return (LocalToWorldNormal * 
                m_gradient(a_intersection.LocalPos)).Normalized;
        }
    }

    public override Vector2 GetUV(Intersection a_intersection)
    {
        return UVMapper.Map(a_intersection.UVW);  
    }

    public override Vector3 GetUVW(Intersection a_intersection)
    {
        return base.GetUVW(a_intersection) * 0.5 + Vector3.HALF;
    }

    public override string ToString()
    {
        return String.Format("PolynomialSurface: {0}", Name);
    }

    public override void GetTangents(Intersection a_intersection, 
        out Vector3 a_tangent_x, out Vector3 a_tangent_y)
    {
        if (a_intersection.BackHit)
        {
            a_tangent_x = Vector3.CrossProduct(
                Up, a_intersection.Normal).Normalized;
            a_tangent_y = Vector3.CrossProduct(
                a_tangent_x, a_intersection.Normal).Normalized;
        }
        else
        {
            a_tangent_x = Vector3.CrossProduct(
                a_intersection.Normal, Up).Normalized;
            a_tangent_y = Vector3.CrossProduct(
                a_intersection.Normal, a_tangent_x).Normalized;
        }
    }

    public override Scene Scene
    {
        get
        {
            return base.Scene;
        }
        set
        {
            if (Scene != null)
                Scene.RenderEnd -= RenderEnd;

            base.Scene = value;

            if (Scene != null)
                Scene.RenderEnd += RenderEnd;
        }
    }

    private void RenderEnd(bool a_all)
    {
        if (!a_all)
            return;

        if (LocalBoundBox != AABB.INFINITY)
        {
            if ((!LocalBoundBox.Minimum.IsAlmostRelativeEquals(m_aabb_minimum)) &&
                (!LocalBoundBox.Maximum.IsAlmostRelativeEquals(m_aabb_maximum)))
            {
                LogManager.GetCurrentClassLogger().Error("Bound box for {0}", Name);
                LogManager.GetCurrentClassLogger().Error("Actually: {0}", 
                    LocalBoundBox);
                LogManager.GetCurrentClassLogger().Error("Should be: {0}", 
                    new AABB(m_aabb_minimum, m_aabb_maximum));
            }
        }
    }

    protected override AABB GetLocalBoundBox()
    {
        return LocalBoundBox;
    }

    protected override void RenderStart(RenderStartPhase a_phase)
    {
        base.RenderStart(a_phase);

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            string equation = ApplyParameters(Equation);
            m_polynomial = new PolynomialFunction(equation);
            m_gradient = m_polynomial.Gradient;
            m_intersection = m_polynomial.Intersection;

            m_solver = new PolynomialSolver(PolynomialSolveMethod.Sturm,
                Scene.RenderOptions.RootFindingMethod);

            m_aabb_minimum = Vector3.MAXIMUM;
            m_aabb_maximum = Vector3.MINIMUM;
        }
    }

    private string ApplyParameters(string a_equation)
    {
        a_equation = ApplyParameter(a_equation, Parameter1);
        a_equation = ApplyParameter(a_equation, Parameter2);
        a_equation = ApplyParameter(a_equation, Parameter3);

        return a_equation;
    }

    private string ApplyParameter(string a_equation, string a_parameter)
    {
        if (String.IsNullOrWhiteSpace(a_parameter))
            return a_equation;

        string[] ar = a_parameter.Split('=');

        if (ar.Length < 2)
            return a_equation;

        a_equation = a_equation.Replace(ar[0], "(" + ar[1].Replace(",", ".") + ")");

        return a_equation;

    }

    public override void ScaleAbsolute(double a_scale)
    {
        Scale *= a_scale;

        base.ScaleAbsolute(a_scale);
    }
}

Twierdzenie Sturma

Twierdzenie Sturma pozwala nam stwierdzić ile rzeczywistych pierwiastków wielomianu $f(x)$ znajduje się w określonym przedziale $\left (a, b \right )$, gdzie $a<b$, $f(a) \neq 0$, $f(b) \neq 0$. Liczba pierwiastków rzeczywistych w takim przedziale to $\sigma(a)-\sigma(b)$, gdzie $\sigma(x)$ to liczba zmian znaków w tzw. sekwencji Sturma. Sekwencje tą budujemy następująco:

$p_0(x)=f(x)$
$p_1(x)=f'(x)$
$p_2(x)=-rem(p_0, p_1)=p_1 q_0 - p_0$
$p_3(x)=-rem(p_1, p_2)=p_2 q_1 - p_1$
...
$p_m(x)=-rem(p_{m-1}, p_{m-2}) = p_{m-1} q_{m-2} - p_{m-2}$

$rem$ oznacza resztę z dzielenia wielomianów.

Przykładowo jeśli obliczeniach otrzymamy: ++--, to $\sigma=1$, jeśli +-+- to $\sigma=2$.

Twierdzenie Sturma działa tylko dla pierwiastków pojedynczych. Jeśli wielomian posiada podwójne pierwiastki musimy się ich najpierw pozbyć budując wielomian bez nich: $\displaystyle ^P/_{GCD(P, P')}$. Przez dzielenie rozumiemy dzielenie wielomianów. GCD to największy wspólny dzielnik. W moim przypadku sekwencje tą wykorzystujemy do obliczeń związanych z raytracingiem. Tak więc póki co nie musiałem się przejmować przypadkami wielokrotnych pierwiastków. Jedynym jest torus zredukowany do sfery.

Tak więc zakładamy, że nasz wielomian posiada tylko pojedyncze pierwiastki.

Dla wielokrotnych twierdzenie nie będzie działało dobrze. Zmiana w znaku zostanie wykryta zawsze dla pierwiastka nieparzystego. Ale może zostać także stwierdzona dla parzystego. Przynajmniej tak to działa do wielomianów stopnia czwartego - sprawdzone doświadczalnie. Innymi słowy nigdy nie dostaniemy fałszywej liczby pierwiastków ale możemy nie dostać prawdziwej. W przypadku torusa zredukowanego do sfery problemem są raczej błędy numeryczne niż ten błąd.

Sekwencja Sturma posiada następujące właściwości:

$1^\circ \quad p_m(x) \neq 0$ \: dla $x \in (a, b)$
$2^\circ \quad \text{jeśli} \:\: 0<i<m-1$ \text{to} $p_{i-1}(x) = p_{i+1}(x) = 0$ \: \text{dla} \: $x \in (a, b)$
$3^\circ \quad \text{jeśli} \:\: 0<i<m$ to $p_{i-1}(x) \neq 0$ i $p_{i}(x) \neq 0$ \: \text{dla} \: $x \in (a, b)$
$4^\circ \quad$ istnieje takie $\alpha$, że dla $p_0(x)=0$ zachodzi $sign(p_0(x-\alpha)) = sign(p_0(x+\alpha))$

$1^\circ$ Wynika z definicji reszty. Musi ona być co najmniej o stopień mniejsza od dzielnej. Stopień $p_1$ też jest o stopień niższy od $p_0$. Stopień wielomianu maleje wraz z wyrazami sekwencji. Jeśli wielomian jest stopnia zerowego to nasza sekwencja liczy tylko jeden wyraz. Dla wielomianu pierwszego stopnia jego pochodna zawsze będzie różna od zera. Dla wielomianów wyższych stopni, ponieważ pierwiastki są pojedyncze to $p_0(x)$ i $p_1(x)$ nie mają wspólnych miejsc zerowych, reszta $p_2$ musi być różna od zera. I nie ma miejsc zerowych wspólnych z $p_1$ i $p_0$. W konsekwencji kolejne reszty nie mogą być zerowe, ich stopień z braku wspólnych miejsc zerowych będzie się zmniejszał zawsze o jeden do $p_{m-1}=ax-b$, $p_m=a$.

$2^\circ$ Sekwencja musi mieć co najmniej 3 wyrazy. Jeśli $p_1=0$ to z definicji wynika, że $p_2=p_1 q_0 - p_0 = -p_0$. Jeśli $p_i=0$, gdzie $i>1$ to z definicji na ogólny wyraz ciągu: $p_i(x)= p_{i-1} q_{i-2} - p_{i-2}$, wynika, że $p_{i+1}(x)= p_{i} q_{i-1} - p_{i-1}$, dla $p_i=0$ mamy: $p_{i+1}=-p_{i-1}$.

$3^\circ$ Wynika z rozważań z punktu $1^\circ$. Stopień wielomianu w sekwencji zawsze maleje od jeden do stałej i nie mają one wspólnych miejsc zerowych.

$4^\circ$ Wynika z definicji pochodnej. Pochodna wielomianu osiąga zera tylko w ekstremach wielomianu. Pomiędzy nimi, czyli w obszarze występowania miejsc zerowych wielomianu jest monotoniczna. O ile miejsca zerowe nie są wielokrotne, wtedy miejsce zerowe może być jednocześnie ekstremum. Dla pojedynczych miejsc zerowych wielomian zawsze przechodzi przez miejsce zerowe i zmienia znak.

Wszystkie te cztery właściwości pozwalają nam udowodnić poprawność twierdzenia. Udowodnimy, że funkcja $\sigma$ wyznaczająca ilość zmian znaków jest malejąca i maleje o jeden tylko przy przejściu przez miejsce zerowe wielomianu.

Sekwencja składa się z wielomianów, których miejsca zerowe są pojedyncze. Jeśli przechodzimy przez zero w wyrazie sekwencji wyraz ten zmienia znak. Z $3^\circ$ wynika, że taki wyraz otoczony jest przez wyrazy niezerowe. Ale nie jest nigdzie powiedziane, że inne wyrazy nie mogą być zerami.

Przypuśćmy, że miejsce zerowe nastąpiło w jednym lub więcej wyrazach różnych od $p_0$. Zgodnie z $4^\circ$ i $2^\circ$ mamy cztery kombinacje znaków przed miejscem zerowym: -++, --+, +--, ++- i po:
--+, -++, ++-, +--. tak więc liczba znaków w sekwencji nie ulega zmianie.

Przypuśćmy, że miejsce zerowe wystąpiło w wyrazie $p_0$ (może też wystąpić w innych wyrazach sekwencji). Zgodnie z $4^\circ$ $p_1$ nie zmienia znaku w otoczeniu miejsca zerowego. Same wielomian $p_0$ zmienia znak z uwagi na to, że posiada on tylko pojedyncze miejsca zerowe. Jeśli wielomian $p_0$ zmienia znak z + na - to znak pochodnej to -, tak więc możliwa kombinacja znaków przed to: +-, -+, zaś po: --, ++.

Widzimy więc, że tylko podczas przejścia przez miejsca zerowe wielomianu $p_0$ liczba zmian znaków w sekwencji maleje o jeden. Możemy więc powiedzieć, że liczba zmian znaków w przedziale $(a, b)$ równa się $\sigma(a) - \sigma(b)$. Ponieważ zmiana znaków następuje przy przejściu przez miejsce zerowe wielomianu $p_0$ jest to tym samym liczba miejsc zerowych w tym przedziale.

Twierdzenie Sturma wykorzystuje się do znajdowania miejsc zerowych wielomianów. Dla wielomianów stopnia 5 i wyższego nie istnieje rozwiązanie dane wzorami wprost. Poza tym wzory na pierwiastki wielomianu stopnia 3 i 4 są niestabilne numerycznie (dokładność double jest niewystarczająca by policzyć je dokładnie, lub by policzyć je dobrze kiedy są blisko siebie - wtedy często wyliczamy jako pierwiastek ekstremum między nimi). Najpierw metodą np. binarnego poszukiwania wyodrębniamy obszary w których zmiana znaku jest pojedyncza, a następnie jedną z metod aproksymacji (Newton, ModifedRegulaFasli, Bisection) przybliżamy się do miejsca zerowego, korzystając z faktu, że jest pojedyncze, i znak z jednej jego strony jest przeciwny do znaku z drugiej strony.

W przypadku obliczeń do raytracingu, interesuje nas tylko pierwsze rozwiązanie większe od zera, a ściśle od Constants.MINIMUM_DISTANCE.

Kod klasy implementującej sekwencje Sturma:

public struct SturmSequence
{
    private readonly Polynomial[] m_seq;
    public readonly int Length;

    public SturmSequence(Polynomial a_poly)
    {
        a_poly = a_poly.RemoveZeroes();

        m_seq = new Polynomial[a_poly.Order + 1];

        Length = 1;
        m_seq[0] = a_poly;

        if (a_poly.Order == 0)
            return;

        Length++;
        m_seq[1] = a_poly.Differentiate();

        while (m_seq[Length - 1].Order != 0)
        {
            Length++;
            m_seq[Length - 1] = -(m_seq[Length - 3] % m_seq[Length - 2]);
            if (m_seq[Length - 1].Order == 0)
                return;
        }
    }

    public static Polynomial SingleRoots(Polynomial a_poly, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        return a_poly / Polynomial.GCD(a_poly, a_poly.Differentiate(), a_precision);
    }

    public Polynomial this[int a_index]
    {
        get
        {
            if (a_index >= Length)
                throw new IndexOutOfRangeException();

            return m_seq[a_index];
        }
    }

    public int NumberOfSignChangesAt(double a_x)
    {
        int sc = 0;

        double c1 = m_seq[0].Evaluate(a_x);

        for (int i = 1; i < Length; i++)
        {
            double c2 = m_seq[i].Evaluate(a_x);

            if ((c1 * c2 < 0) || (c1 == 0.0))
                sc++;

            c1 = c2;
        }

        return sc;
    }

    public int NumberOfSignChangesAtPositiveInfinity()
    {
        int sc = 0;

        double c1 = m_seq[0].LeadCoefficient;

        for (int i = 1; i < Length; i++)
        {
            double c2 = m_seq[i].LeadCoefficient;

            if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                sc++;

            c1 = c2;
        }

        return sc;
    }

    public int NumberOfSignChangesAtNegativeInfinity()
    {
        int sc = 0;

        double c1 = -m_seq[0].LeadCoefficient;
        if (m_seq[0].Order % 2 == 0)
            c1 = -c1;

        for (int i = 1; i < Length; i++)
        {
            double c2 = -m_seq[i].LeadCoefficient;
            if (m_seq[i].Order % 2 == 0)
                c2 = -c2;

            if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
                sc++;

            c1 = c2;
        }

        return sc;
    }
}
Kod klasy solvera wykorzystującego powyższą sekwencje znajdujący tylko jedno miejsce zerowe większe od Constants.MINIMUM_DISTANCE większe od zera:
public static class Constants
{
    public const int STURM_ANY_ROOTS_MAXIMUM_ITERATIONS = 32;
    public const int STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS = 800;

    public const double MINIMAL_DISTANT = 1e-8;
}

public struct PolynomialSolverAprox
{
    private readonly SturmSequence m_sturm_sequence;
    private RootFinder m_root_finder;

    public PolynomialSolverAprox(Polynomial a_polynomial, RootFinder a_root_finder)
    {
        m_sturm_sequence = new SturmSequence(a_polynomial);
        m_root_finder = a_root_finder;
    }

    public double FindOneRootInRange(double a_a, double a_b, int a_num_roots_at_a)
    {
        for (int i = 0; 
             i < Constants.STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS; i++)
        {
            double mid = (a_a + a_b) / 2;
            int atmid = m_sturm_sequence.NumberOfSignChangesAt(mid);
            int nroots = a_num_roots_at_a - atmid;

            if (nroots == 1)
            {
                Polynomial p = m_sturm_sequence[0];
                return m_root_finder.FindRoot(a_a, mid, (x) => p.Evaluate(x));
            }
            else if (nroots == 0)
                a_a = mid;
            else
                a_b = mid;
        }

        return a_a;
    }

    public double SolveOne()
    {
        int nroots_at_a = m_sturm_sequence.NumberOfSignChangesAt(
           Constants.MINIMAL_DISTANT);
        int nroots_at_pos_inf = 
            m_sturm_sequence.NumberOfSignChangesAtPositiveInfinity();

        int nroots = nroots_at_a - nroots_at_pos_inf;

        if (nroots == 0)
            return Double.PositiveInfinity;

        double a = Constants.MINIMAL_DISTANT;
        double b = a + 1.5;
        int nroots_at_b = nroots_at_a;

        for (int its = 0; its < Constants.STURM_ANY_ROOTS_MAXIMUM_ITERATIONS; its++)
        {
            nroots_at_b = m_sturm_sequence.NumberOfSignChangesAt(b);

            if (nroots_at_a - nroots_at_b >= 1)
                break;

            b *= 2.0;
        }

        if (nroots_at_a == nroots_at_b)
            return Double.PositiveInfinity;

        return FindOneRootInRange(a, b, nroots_at_a);
    }
}