Sekwencja Hammersley'a

Jest to pewna modyfikacja sekwencji Halton'a. Dla jednego z wymiarów sekwencja ma postać: $\displaystyle\frac{n-0.5}{N}$, gdzie n to numer kolejnej próbki. N to liczba próbek. Sekwencja ta ma mniejszą dyspersję niż sekwencja Haltona, za cenę określenia z góry ilości generowanych próbek.

-0.5 w wzorze zapewnia nam centrowanie punktów na środkach pikseli.

Sekwencja ta degraduje się dla wyższych wymiarów. W przypadku 2D bardzo szybko.

Przykłady:

2

3

4

5

Widzimy, że być może i jest to sekwencja o niskiej dyspersji, ale punkty zdecydowanie układają się w wzór, który jak każda powtarzalność w próbkowaniu podatny jest na aliasing. Widać to też na poniższym periodogramie.

Przykład pierodogramu dla bazy 3:

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 class HammersleySampler : LowDiscrepancySequenceSampler
{
    private int m_n = 1;
    private List<Vector2>[,] m_samples;
    private int m_count;

    public int BaseY = 3;

    private double NextX()
    {
        return (m_n - 0.5) / m_count;
    }

    public override IEnumerable<Vector2> GetSamples(Rectangle a_rect)
    {
        foreach (var p in a_rect.EnumPixels())
        {
            if (m_samples[p.X, p.Y] == null)
                continue;

            foreach (var s in m_samples[p.X, p.Y])
                yield return s;
        }
    }

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

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

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
            PrepareSamples();
    }

    private void PrepareSamples()
    {
        m_samples = new List<Vector2>[Film.Width, Film.Height];
        m_count = Film.Width * Film.Height * Subresolution * Subresolution;

        for (int i=0; i<m_count; i++)
        {
            Vector2 s = new Vector2(
                NextX() * Film.Width, 
                LowDiscrepancyMath.RadicalInverse(m_n, BaseY) * Film.Height);
            Point p = new Point((int)s.X, (int)s.Y);

            if (m_samples[p.X, p.Y] == null)
                m_samples[p.X, p.Y] = new List<Vector2>();
            m_samples[p.X, p.Y].Add(s);

            m_n++;
        }
    }
}