Sekwencje Van Der Corput'a i Halton'a

W kolejnych postach zostanie omówionych kilka sekwencji, użytecznych dla samplowania. Wszystkie są quasi-losowe, czyli nie losowe w sensie pseudo-losowe, ale powtarzalne, wykazujące pewien wzór, ale wyglądające na losowe. Wszystkie sekwencje generują liczby z zakresu $<0, 1)$. Punkty sekwencji są w miarę równomiernie od siebie oddalone i w miarę równomiernie wypełniają przestrzeń (niska dyspersja). Warunek ten jest spełniony po wygenerowaniu dowolnych n wyrazów sekwenji. Większość sekwencji liczona jest wprost, jedynym parametrem jest numer wyrazu sekwencji.

Pierwszą z sekwencji jest sekwencja Van Der Corput'a. Parametrem do generacji serii punktów jest Base - podstawa systemu liczbowego dla której będzie generowana sekwencja, czyli powinna być większa od 1 (2 dla binarnego, 10 dla dziesiętnego).

Sekwencja wielowymiarowe, których bazy są liczbami relatywnie pierwszymi dla siebie to sekwencja Haltona.

Parametrem do generacji poszczególnych liczb w sekwencji jest $n \in {1,2,3,4,...}$. Dla sekwencji wielowymiarowych dla każdego wymiaru powinniśmy użyć inną bazę by uniknąć degradacji. Poza tym by uniknąć degradacji bazy sekwencji wielowymiarowych nie powinny mieć wspólnych podzielników - powinny być liczbami relatywnie pierwszymi (coprime).

Generowanie sekwencji opisowo wygląda tak. Weźmy kolejną liczbę dodatnią N. Potraktujmy ją jako liczbę w systemie o zadanej bazie $abcd_{BASE}$. Odpowiednik w sekwencji ma postać: $\frac{d}{BASE}+\frac{c}{BASE^2}+\frac{b}{BASE^3}+\frac{a}{BASE^4}$. Np. dla systemu dziesiętnego: $214 \rightarrow 0.412$.

Kod:

public static class LowDiscrepancyMath
{
    public static double RadicalInverse(int a_n, int a_base)
    {
        double y = 0;
        double b = a_base;
        while (a_n > 0)
        {
            int d_i = a_n % a_base;
            y += d_i / b;
            a_n /= a_base;
            b *= a_base;
        }
        return y;
    }  
}

public abstract class LowDiscrepancySequenceSampler : NonUniformSampler
{
    public int Subresolution = 3;

    public override IEnumerable<Vector2> GetSamples(Rectangle a_rect)
    {
        int samples_x = a_rect.Width * Subresolution;
        int samples_y = a_rect.Height * Subresolution;

        for (int y = 0; y < samples_y; y++)
        {
            for (int x = 0; x < samples_x; x++)
            {
                Vector2 sample = NextSample();
                yield return new Vector2(
                    a_rect.Left + sample.X * a_rect.Width,
                    a_rect.Top + sample.Y * a_rect.Height);
            }
        }
    }

    protected abstract Vector2 NextSample();
}

public class HaltonSampler : LowDiscrepancySequenceSampler
{
    protected int m_n = 1;

    public int BaseX = 2;
    public int BaseY = 3;

    protected override Vector2 NextSample()
    {
        int n = m_n;
        m_n++;

        return new Vector2(
            LowDiscrepancyMath.RadicalInverse(n, BaseX),
            LowDiscrepancyMath.RadicalInverse(n, BaseY));
    }

    public override SamplerType SamplerType
    {
        get
        {
            return SamplerType.Halton;
        }
    }
}

Przykłady sekwencji:

(2,3)

(2,4)

(2,5)

Poniższy przykład zakłada sekwencję, w której na każdy piksel przypada jedna próbka. W przypadku samplowania typu Grid byłby to biały prostokąt. Tutaj możemy zauważyć jak dla tej sekwencji co pewną ilość punktów następuje duża zmiana w rozkładzie próbek utrzymująca się przez pewien czas.

Periodogram dla (2,3):

Widać na nim pewną powtarzalność i tym samym podatność na aliasing.