Klasa liczby zespolonej

Moja implementacja klasy liczby zespolonej, nie jest może wszystko-mająca, ale do FFT wystarczy. Do tego testy. Publikuje to dlatego, że przechodzę na implementację klasy Complex z System.Numerics. Jedyną różnicę jaką zauważyłem to zachowanie się w przypadku liczenia przesunięcia fazowego z (0,0). Zgodnie z definicją z wikipedii i nazwijmy to zdrowym rozsądkiem wartość ta powinna być nieokreślona, tymczasem w System.Numerics jest to zero.

public struct Complex
{
    public static readonly Complex ZERO = new Complex(0, 0);
    public static readonly Complex ONE = new Complex(1, 0);
    public static readonly Complex i = new Complex(0, 1);
    public static readonly Complex UNDEFINED = new Complex(Double.NaN, Double.NaN);

    public readonly double Real;
    public readonly double Image;

    public Complex(double a_real, double a_image = 0)
    {
        Real = a_real;
        Image = a_image;
    }

    public static bool operator ==(Complex a_a, Complex a_b)
    {
        return (a_a.Real == a_b.Real) & (a_a.Image == a_b.Image);
    }

    public static bool operator !=(Complex a_a, Complex a_b)
    {
        return (a_a.Real != a_b.Real) | (a_a.Image != a_b.Image);
    }

    public static Complex operator +(Complex a_a)
    {
        return a_a;
    }

    public static Complex FromPolar(double a_magnitude, double a_phase)
    {
        return new Complex(a_magnitude * Math.Cos(a_phase), 
            a_magnitude * Math.Sin(a_phase));
    }

    public static Complex FromPolar(double a_phase)
    {
        return new Complex(Math.Cos(a_phase), Math.Sin(a_phase));
    }

    public static Complex operator -(Complex a_a)
    {
        return new Complex(-a_a.Real, -a_a.Image);
    }

    public static Complex operator +(Complex a_a, Complex a_b)
    {
        return new Complex(
            a_a.Real + a_b.Real, 
            a_a.Image + a_b.Image);
    }

    public static Complex operator -(Complex a_a, Complex a_b)
    {
        return new Complex(
            a_a.Real - a_b.Real, 
            a_a.Image - a_b.Image);
    }

    public static Complex operator *(Complex a_a, Complex a_b)
    {
        return new Complex(
            a_a.Real * a_b.Real - a_a.Image * a_b.Image, 
            a_a.Real * a_b.Image + a_a.Image * a_b.Real);
    }

    public static Complex operator *(Complex a_a, double a_b)
    {
        return new Complex(
            a_a.Real * a_b,
            a_a.Image * a_b);
    }

    public static Complex operator *(double a_a, Complex a_b)
    {
        return new Complex(
            a_b.Real * a_a,
            a_b.Image * a_a);
    }

    public static Complex operator /(Complex a_a, Complex a_b)
    {
        double den = a_b.Real * a_b.Real + a_b.Image * a_b.Image;

        return new Complex(
            (a_a.Real * a_b.Real + a_a.Image * a_b.Image) / den,
            (a_a.Image * a_b.Real - a_a.Real * a_b.Image) / den);
    }

    public Complex Scale(double a_scale)
    {
        return new Complex(
            Real * a_scale,
            Image * a_scale);
    }

    public bool IsAlmostEquals(Complex a_a, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        return Real.IsAlmostEquals(a_a.Real, a_precision) &&
            Image.IsAlmostEquals(a_a.Image, a_precision);
    }

    public bool IsAlmostRelativeEquals(Complex a_a, 
        double a_precision = Constants.DOUBLE_PRECISION)
    {
        return Real.IsAlmostRelativeEquals(a_a.Real, a_precision) &&
            Image.IsAlmostRelativeEquals(a_a.Image, a_precision);
    }

    public override int GetHashCode()
    {
        return Real.GetHashCode() ^ Image.GetHashCode();
    }

    public override bool Equals(object a_obj)
    {
        if (a_obj == null)
            return false;

        if (!typeof(Complex).Equals(a_obj.GetType()))
            return false;
        Complex complex = (Complex)a_obj;

        return (complex.Real == Real) && (complex.Image == Image);
    }

    public bool Equals(Complex a_complex)
    {
        return (a_complex.Real == Real) && (a_complex.Image == Image);
    }

    [DebuggerHidden]
    public Complex Conjugate
    {
        get
        {
            return new Complex(Real, -Image);
        }
    }

    public double Abs
    {
        get
        {
            return MathExtensions.Hypot(Real, Image);
        }
    }

    public double Phase
    {
        get
        {
            if (IsAlmostRelativeEquals(ZERO))
                return Double.NaN;

            return Math.Atan2(Image, Real);
        }
    }

    public override string ToString()
    {
        if ((Real != 0) && (Image != 0))
            return String.Format("({0} + {1}i)", Real, Image);
        else if (Real != 0)
            return Real.ToString();
        else if (Image != 0)
            return String.Format("{1}i", Image);
        else
            return "0";
    }
}

Testy:

[TestMethod]
public void Complex_Test()
{
    // Constructor
    {
        Assert.AreEqual(4, new Complex(4).Real);
        Assert.AreEqual(0, new Complex(4).Image);
        Assert.AreEqual(5, new Complex(5, 6).Real);
        Assert.AreEqual(6, new Complex(5, 6).Image);
    }

    // ==, !=
    {
        Assert.IsTrue(new Complex(5, 4) == new Complex(5, 4));
        Assert.IsTrue(new Complex(5.3, 4.7) == new Complex(5.3, 4.7));
        Assert.IsFalse(new Complex(5, 4) == new Complex(5, 5));
        Assert.IsFalse(new Complex(5, 4) == new Complex(4, 4));
        Assert.IsFalse(new Complex(5, 5) == new Complex(5, 4));
        Assert.IsFalse(new Complex(4, 4) == new Complex(5, 4));

        Assert.IsFalse(new Complex(5, 4) != new Complex(5, 4));
        Assert.IsTrue(new Complex(5, 4) != new Complex(5, 5));
        Assert.IsTrue(new Complex(5, 4) != new Complex(4, 4));
        Assert.IsTrue(new Complex(5, 5) != new Complex(5, 4));
        Assert.IsTrue(new Complex(4, 4) != new Complex(5, 4));
    }

    // unary -
    {
        Assert.IsTrue(new Complex(-1, 1) == -new Complex(1, -1));
        Assert.IsTrue(new Complex(-1.6, 1.5) == -new Complex(1.6, -1.5));
        Assert.IsTrue(new Complex(-1, -1) == -new Complex(1, 1));
        Assert.IsTrue(new Complex(1, -1) == -new Complex(-1, 1));
        Assert.IsTrue(new Complex(1, 0) == -new Complex(-1, 0));
        Assert.IsTrue(new Complex(0, 0) == -new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 1) == -new Complex(0, -1));

        Assert.IsTrue(-new Complex(-1, 1) == - -new Complex(1, -1));
        Assert.IsTrue(-new Complex(-1, -1) == - -new Complex(1, 1));
        Assert.IsTrue(-new Complex(1, -1) == - -new Complex(-1, 1));
        Assert.IsTrue(-new Complex(1, 0) == - -new Complex(-1, 0));
        Assert.IsTrue(-new Complex(0, 0) == - -new Complex(0, 0));
        Assert.IsTrue(-new Complex(0, 1) == - -new Complex(0, -1));
    }

    // unary +
    {
        Assert.IsTrue(-new Complex(-1, 1) == +new Complex(1, -1));
        Assert.IsTrue(-new Complex(-1, -1) == +new Complex(1, 1));
        Assert.IsTrue(-new Complex(-1.2, -1.5) == +new Complex(1.2, 1.5));
        Assert.IsTrue(-new Complex(1, -1) == +new Complex(-1, 1));
        Assert.IsTrue(-new Complex(1, 0) == +new Complex(-1, 0));
        Assert.IsTrue(-new Complex(0, 0) == +new Complex(0, 0));
        Assert.IsTrue(-new Complex(0, 1) == +new Complex(0, -1));
    }

    // +
    {
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) + new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) + new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, 0) + new Complex(0, 1));
        Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0) + new Complex(0, 0));
        Assert.IsTrue(new Complex(1, 0) == new Complex(0, 0) + new Complex(1, 0));
        Assert.IsTrue(new Complex(1, 1) == new Complex(1, 0) + new Complex(0, 1));
        Assert.IsTrue(new Complex(1, 1) == new Complex(0, 1) + new Complex(1, 0));
        Assert.IsTrue(new Complex(3, 3) == new Complex(1, 2) + new Complex(2, 1));
        Assert.IsTrue(new Complex(4.1, 3.7) == new Complex(1.4, 2.6) + 
            new Complex(2.7, 1.1));
    }

    // -
    {
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) - -new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) - -new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, 0) - -new Complex(0, 1));
        Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0) - -new Complex(0, 0));
        Assert.IsTrue(new Complex(1, 0) == new Complex(0, 0) - -new Complex(1, 0));
        Assert.IsTrue(new Complex(1, 1) == new Complex(1, 0) - -new Complex(0, 1));
        Assert.IsTrue(new Complex(1, 1) == new Complex(0, 1) - -new Complex(1, 0));
        Assert.IsTrue(new Complex(3, 3) == new Complex(1, 2) - -new Complex(2, 1));
        Assert.IsTrue(new Complex(3.8, 3.8) == new Complex(1.7, 2.6) - 
            -new Complex(2.1, 1.2));
    }

    // * Complex
    {
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 1) * new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(0, 1));
        Assert.IsTrue(new Complex(0, 0) == new Complex(1, 0) * new Complex(0, 0));
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0) * new Complex(1, 0));
        Assert.IsTrue(new Complex(0, 1) == new Complex(1, 0) * new Complex(0, 1));
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, 1) * new Complex(1, 0));
        Assert.IsTrue(new Complex(0, 5) == new Complex(1, 2) * new Complex(2, 1));
        Assert.IsTrue(new Complex(-14, 31) == new Complex(2, 3) * 
            new Complex(5, 8));
        AreAlmostRelativeEquals(new Complex(-14.32, 37.47), 
            new Complex(2.3, 3.2) * new Complex(5.6, 8.5));
    }

    // Complex * double, double * Complex, Scale
    {
        Action<Complex, Complex, double> scale_test = (c1, c2, s) =>
        {
            Assert.IsTrue(c1 == c2 * s);
            Assert.IsTrue(c1 == s * c2);
            Assert.IsTrue(c1 == c2.Scale(s));
        };

        scale_test(new Complex(0, 0), new Complex(0, 0), 0);
        scale_test(new Complex(0, 0), new Complex(4, 4), 0);
        scale_test(new Complex(0, 0), new Complex(4, 0), 0);
        scale_test(new Complex(0, 0), new Complex(0, 4), 0);
        scale_test(new Complex(0, 0), new Complex(0, 0), 4);
        scale_test(new Complex(16, 16), new Complex(4, 4), 4);
        scale_test(new Complex(16, 0), new Complex(4, 0), 4);
        scale_test(new Complex(0, 16), new Complex(0, 4), 4);
        scale_test(new Complex(0, 0), new Complex(0, 0), -4);
        scale_test(new Complex(-16, -16), new Complex(4, 4), -4);
        scale_test(new Complex(-16, 0), new Complex(4, 0), -4);
        scale_test(new Complex(0, -16), new Complex(0, 4), -4);
        scale_test(new Complex(-16, -12), new Complex(4, 3), -4);
        scale_test(new Complex(-16, 12), new Complex(4, -3), -4);
        scale_test(new Complex(16, -12), new Complex(-4, 3), -4);
        scale_test(new Complex(17.63, -13.12), new Complex(-4.3, 3.2), -4.1);
    }

    // /
    {
        AreAlmostRelativeEquals(new Complex(23.0 / 41.0, 2.0 / 41.0), 
            new Complex(2, 3) / new Complex(4, 5));
        AreAlmostRelativeEquals(new Complex(23.0 / 13.0, -2.0 / 13.0), 
            new Complex(4, 5) / new Complex(2, 3));
        AreAlmostRelativeEquals(
            new Complex(1.4896226415094338, -0.13867924528301892), 
            new Complex(4.4, 5.3) / new Complex(2.6, 3.8));
    }

    // IsAlmostEquals
    {
        Complex c = new Complex(0, 0);

        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, -0.5) * 
            Constants.DOUBLE_PRECISION));

        c = new Complex(40, -40);

        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(2, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-2, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(0.5, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, 2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsFalse(c.IsAlmostEquals(c + new Complex(-0.5, -2) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(0.5, -0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, 0.5) * 
            Constants.DOUBLE_PRECISION));
        Assert.IsTrue(c.IsAlmostEquals(c + new Complex(-0.5, -0.5) * 
            Constants.DOUBLE_PRECISION));
    }

    // IsAlmostRelativeEquals
    {
        Complex c = new Complex(1 / Constants.DOUBLE_PRECISION, 
            1 / Constants.DOUBLE_PRECISION);

        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, 2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, -2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, 2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, -2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, 0.5)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, 0.5)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(2, -0.5)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-2, -0.5)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(0.5, 2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(0.5, -2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-0.5, 2)));
        Assert.IsFalse(c.IsAlmostRelativeEquals(c + new Complex(-0.5, -2)));
        Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(0.5, 0.5)));
        Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(0.5, -0.5)));
        Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(-0.5, 0.5)));
        Assert.IsTrue(c.IsAlmostRelativeEquals(c + new Complex(-0.5, -0.5)));
    }

    // Equals(obj), Equals(complex)
    {
        Assert.AreEqual(new Complex(5, 4), new Complex(5, 4));
        Assert.AreEqual(new Complex(5.4, 4.4), new Complex(5.4, 4.4));
        Assert.AreNotEqual(new Complex(5, 4), new Complex(5, 5));
        Assert.AreNotEqual(new Complex(5, 4), new Complex(4, 4));
        Assert.AreNotEqual(new Complex(5, 5), new Complex(5, 4));
        Assert.AreNotEqual(new Complex(4, 4), new Complex(5, 4));

        Assert.IsTrue(new Complex(5, 4).Equals(new Complex(5, 4)));
        Assert.IsTrue(new Complex(5.3, 4.3).Equals(new Complex(5.3, 4.3)));
        Assert.IsFalse(new Complex(5, 4).Equals(new Complex(5, 5)));
        Assert.IsFalse(new Complex(5, 4).Equals(new Complex(4, 4)));
        Assert.IsFalse(new Complex(5, 5).Equals(new Complex(5, 4)));
        Assert.IsFalse(new Complex(4, 4).Equals(new Complex(5, 4)));
    }

    // Conjugate
    {
        Assert.IsTrue(new Complex(0, 0) == new Complex(0, 0).Conjugate);
        Assert.IsTrue(new Complex(1, 0) == new Complex(1, 0).Conjugate);
        Assert.IsTrue(new Complex(0, -1) == new Complex(0, 1).Conjugate);
        Assert.IsTrue(new Complex(0, 1) == new Complex(0, -1).Conjugate);
        Assert.IsTrue(new Complex(2, -3) == new Complex(2, 3).Conjugate);
        Assert.IsTrue(new Complex(2, 3) == new Complex(2, -3).Conjugate);
        Assert.IsTrue(new Complex(-2, -3) == new Complex(-2, 3).Conjugate);
        Assert.IsTrue(new Complex(-2, 3) == new Complex(-2, -3).Conjugate);
        Assert.IsTrue(new Complex(-2.2, 3.2) == new Complex(-2.2, -3.2).Conjugate);
    }

    // Abs
    {
        Assert.IsTrue(0 == new Complex(0, 0).Abs);
        Assert.IsTrue(1 == new Complex(0, 1).Abs);
        Assert.IsTrue(1 == new Complex(0, -1).Abs);
        Assert.IsTrue(1 == new Complex(1, 0).Abs);
        Assert.IsTrue(1 == new Complex(-1, 0).Abs);
        Assert.IsTrue(5 == new Complex(3, 4).Abs);
        Assert.IsTrue(5 == new Complex(-3, 4).Abs);
        Assert.IsTrue(5 == new Complex(3, -4).Abs);
        Assert.IsTrue(5 == new Complex(-3, -4).Abs);
        Assert.IsTrue(2*MathExtensions.SQRT2 == new Complex(2, 2).Abs);
    }

    // Phase
    {
        Assert.AreEqual(Double.NaN,
            MathExtensions.ToDeg(new Complex(0, 0).Phase));
        Assert.AreEqual(Double.NaN,
            MathExtensions.ToDeg(new Complex(Constants.DOUBLE_PRECISION / 2, 
            Constants.DOUBLE_PRECISION / 2).Phase));
        Assert.AreEqual(0,
            MathExtensions.ToDeg(new Complex(1, 0).Phase));
        Assert.AreEqual(45,
            MathExtensions.ToDeg(new Complex(1, 1).Phase));
        Assert.AreEqual(90,
            MathExtensions.ToDeg(new Complex(0, 1).Phase));
        Assert.AreEqual(135,
            MathExtensions.ToDeg(new Complex(-1, 1).Phase));
        Assert.AreEqual(180,
            MathExtensions.ToDeg(new Complex(-1, 0).Phase));
        Assert.AreEqual(225 - 360,
            MathExtensions.ToDeg(new Complex(-1, -1).Phase));
        Assert.AreEqual(270 - 360,
            MathExtensions.ToDeg(new Complex(0, -1).Phase));
        Assert.AreEqual(315 - 360,
            MathExtensions.ToDeg(new Complex(1, -1).Phase));
    }
}

Samplowanie Poissona

Samplowanie Poissona powinno nam dać równomiernie losowo rozłożone próbki. Próbki są rozłożone losowo. W przeciwieństwie do samplowania losowego (Jitter) próbki są rozłożone równomiernie (są w miarę równych odstępach od siebie). Sampler ten nie generuje ściśle określonej liczby próbek na piksel - może być ich więcej, albo mniej. Czym odstępstwo od średniej większe tym rzadsze.

Generowanie próbek tą metodą jest wyjątkowo wolne. Złożoność obliczeniowa to $N log N$ W kodzie używam i tak już zmodyfikowanej metody, której złożoność obliczeniowa to N.

Rozkład punktów wygenerowanych przez sampler jest porównywalny do rozkładu fotoreceptorów siatkówki.

Aby wygenerować zbiór punktów o rozkładzie Poisson-Disc bierzemy losowy punk i kładziemy go w przestrzeni którą chcemy zapełnić (może być wielowymiarowa). Następnie losujemy pewien zbiór punktów (liczba prób zakończona w kodzie jako Repeats) i wybieramy ten, który jest najbardziej oddalony od już postawionych (minimalny dystans do już postawionych jest maksymalny) oraz nie jest bliżej niż pewien promień $r_{MIN}$, który w praktyce wyznacza nam gęstość rozkładu. Algorytm ten jest bardzo wolny. Po skończonym czasie może się zdarzyć (powinno się zdarzyć jeśli mamy na myśli skończony bardzo długi czasu, że nie możemy już dodać żadnych nowych punktów. W praktyce algorytm trzeba przerwać dużo wcześniej, co z drugiej strony może dać nam rozkład z dziurami. Możemy sobie wyobrazić metodę w której losowanie odbywa się tylko w obszarach, w których mamy pewność, że dziury są.

W szybszej wersji tworzymy siatkę o rozmiarze tak dobranym by średnio jeden punkt wpadał w jedno oczko siatki i w przypadku wylosowania punktu sprawdzamy tylko obszar 3x3 w takiej siatce. Rozmiar oczka siatki wyznaczany jest przez $r_{MIN}$. Dwa punkty w sąsiednich oczkach siatki (także po przekątnych) maksymalnie od siebie oddalone nie zostawią w środku miejsca na kolejny punkt. Ponieważ odległość taka jest największa pod przekątnej rozmiar siatki definiuje się jako: $\frac{r_{MIN}}{\sqrt{2}}$.

Na początku losujemy punkt i dodajemy go do siatki i do listy aktywnej. Następnie tak długo jak lista aktywna nie jest pusta losujemy z niej jeden punkt A. Dla tego punktu losujemy Repeats punktów C wokół niego oddalonych od $0.5 \cdot r_{MIN}$ do $r_{MIN}$. Następnie dla każdego punktu B z C sprawdzamy siatkę 3x3 wokół punktu B. Jeśli punkt B jest oddalony od wszystkich punktów w siatce 3x3 o więcej niż $r_{MIN}$ wstawiamy go siatki i do listy aktywnej. W przypadku gdy nie udało nam się dodać żadnego nowego punktu wokół A usuwamy go z listy aktywnej. Tym samym czym większa liczba prób tym mniej potencjalnych dziur, które da się jeszcze coś wepchnąć.

Przykłady rozkładu próbek:

Kod samplera Poisson:

public class PoissonSampler : NonUniformSampler
{
    public int Repeats = 20;
    public double Density = 1;

    private double m_cell_size;
    private double m_radius;
    private int m_grid_width;
    private int m_grid_height;
    private Random m_random;
    private int m_width;
    private int m_height;
    private int[,] m_grid;
    List<int> m_actives = new List<int>();
    List<Vector2> m_samples = new List<Vector2>();

    private void PrepareSamples()
    {
        m_width = Film.Width;
        m_height = Film.Height;

        if (Film.Scene.RenderOptions.Renderer is SamplesVisualizerRenderer)
        {
            double scale = (Film.Scene.RenderOptions.Renderer as 
                SamplesVisualizerRenderer).Scale;
            m_width = (m_width / scale).Round();
            m_height = (m_height / scale).Round();
        }

        m_radius = 1.0 / Density;
        m_cell_size = m_radius / Math.Sqrt(2);
        m_grid_width = (int)(m_width / m_cell_size) + 1;
        m_grid_height = (int)(m_height / m_cell_size) + 1;

        m_grid = new int[m_grid_width, m_grid_height];
        m_grid.Fill(-1);

        AddFirstPoint();

        while (m_actives.Count != 0)
        {
            int index = m_random.Next(m_actives.Count);

            bool found = false;

            for (int i = 0; i < Repeats; i++)
                found |= AddNextPoint(m_actives[index]);

            if (!found)
                m_actives.RemoveAt(index);
        }
    }

    private Vector2 GenerateRandomAround(int a_index)
    {
        double radius = (m_radius + m_radius * m_random.NextDouble());
        double angle = 2 * MathExtensions.PI * m_random.NextDouble();
        Vector2 center = m_samples[a_index];
        return new Vector2(center.X + radius * Math.Sin(angle), 
            center.Y + radius * Math.Cos(angle));
    }

    private bool AddNextPoint(int a_index)
    {
        Vector2 q = GenerateRandomAround(a_index);

        if ((q.X < 0) || (q.X >= m_width) || (q.Y < 0) || (q.Y >= m_height))
            return false;

        int index_x = GetXIndex(q);
        int index_y = GetYIndex(q);

        for (int x = Math.Max(0, index_x - 2); 
             x < Math.Min(m_grid_width, index_x + 3); x++)
        {
            for (int y = Math.Max(0, index_y - 2); 
                 y < Math.Min(m_grid_height, index_y + 3); y++)
            {
                int index = m_grid[x, y];

                if (index != -1)
                {
                    if ((q - m_samples[index]).Length < m_radius)
                        return false;
                }
            }
        }

        m_grid[index_x, index_y] = m_samples.Count;
        m_actives.Add(m_samples.Count);
        m_samples.Add(q);

        return true;
    }

    private int GetXIndex(Vector2 a_point)
    {
        return (int)(a_point.X / m_cell_size);
    }

    private int GetYIndex(Vector2 a_point)
    {
        return (int)(a_point.Y / m_cell_size);
    }

    private void AddFirstPoint()
    {
        Vector2 p = new Vector2(m_width * m_random.NextDouble(), 
            m_height * m_random.NextDouble());

        m_grid[GetXIndex(p), GetYIndex(p)] = m_samples.Count;
        m_actives.Add(m_samples.Count);
        m_samples.Add(p);
    }

    public override IEnumerable<Vector2> GetSamples(Rectangle a_rect)
    {
        int top = (int)(a_rect.Top / m_cell_size) + 1;
        int bottom = (int)(a_rect.Bottom / m_cell_size) + 1;
        int left = (int)(a_rect.Left / m_cell_size) + 1;
        int right = (int)(a_rect.Right / m_cell_size) + 1;

        for (int x = Math.Max(0, left - 2); 
            x <= Math.Min(m_grid_width - 1, right + 2); x++)
        {
            for (int y = Math.Max(0, top - 2); 
                y <= Math.Min(m_grid_height - 1, bottom + 2); y++)
            {
                if (m_grid[x, y] != -1)
                {
                    Vector2 p = m_samples[m_grid[x, y]];

                    if (p.X < a_rect.Left)
                        continue;
                    if (p.X > a_rect.Right)
                        continue;
                    if (p.Y < a_rect.Top)
                        continue;
                    if (p.Y > a_rect.Bottom)
                        continue;

                    yield return p;
                }
            }
        }
    }
}

Parametr skalowania można pominąć, jest on używany tylko podczas wizualizacji rozkładu próbek i w periodogramie.

Periodogram

Periodogram to zsumowana transformata DFFT serii obrazów. DFFT z pojedynczego obrazu może nam nie pozwolić wychwycić jakiegoś trendu w obrazie, rzeczy które powtarzają się z pewną częstotliwością. Co innego suma takich DFFT.

Periodogram jest przydatny do oceny jakości samplerów. Idealnie wymagamy od samplera by próbki były rozmieszczone w miarę równo od siebie i losowo. W miarę równo po stronie FFT oznacza w miarę równo od siebie odległe piki, tym węższe (mniej rozmyte) czym równiej od siebie rozmieszone. Ponieważ mamy do czynienia z próbkami w 2D piki te powinny tworzyć okręgi, co oznacza, że jeśli weźmiemy próbkę A i kąty wszystkich próbek wokół niej rozmieszczonych to są one równomierny rozmieszczone.

Punkty do generowania perodiogramów musimy przeskalować. Normalnie mamy próbkę na piksel albo i więcej jeśli mamy subsampling. Obraz musi też być w miarę duży by wychwycić wzorzec częstotliwości w pojedynczym obrazie. Odstępami pomiędzy próbkami regulujemy odstępy pomiędzy okręgami na periodogramie. Trochę zabawy i można dobrać parametry jak należy.

Interesuje nas tylko kilka pierwszych okręgów, czym dalej tym bardziej znikają one w szumie i potrzeba więcej obrazów z większą ilością próbek.

Do generacji periodogramów służy osobna klasa renderera:

public class SamplesPeriodogramRenderer : Renderer
{
    [YAXNode]
    public int Repeats = 100;

    [YAXNode]
    public int ScaleSamples = 25;

    [YAXNode]
    public double ScaleValues = 0.01;

    ColorArrayFloat m_periodogram;

    public override void Render()
    {
        m_periodogram = new ColorArrayFloat(Film.Width, Film.Height);

        for (int i = 0; i < Repeats; i++)
        {
            m_rects_to_update.Enqueue(Film.GetFilmRect());
            DoUpdateCallback();

            Scene scene = Scene.DeepClone();
            scene.RenderOptions.Renderer = 
                Renderer.Create(RendererType.SamplesVisualizer);
            (scene.RenderOptions.Renderer as SamplesVisualizerRenderer).Scale = 
                ScaleSamples;
            scene.ActiveCamera.Film.Postprocessors.Add(new FFT2Transform());

            if (i != 0)
                (scene.ActiveCamera.Film.Resampler as ExactResampler).SuppressLog = 
                    true;

            if ((i % 5 == 0) && (i > 0))
                Loggers.Raytracer.Info("Periodogram progress: {0}", i * Repeats / 
                    100);

            Renderer renderer = scene.CreateRenderer();
            renderer.ShouldStop += () => DoShouldStop();
            renderer.Render();

            var part = renderer.Canvas.PostprocessedArray;

            foreach (var p in part.Pixels())
                m_periodogram[p.X, p.Y] += part[p.X, p.Y] * ScaleValues;

            m_rects_to_update.Enqueue(Film.GetFilmRect());
            DoUpdateCallback();

            if (DoShouldStop())
                break;
        }
    }

    public override void UpdateBitmap(Bitmap a_bmp)
    {
        using (Graphics g = Graphics.FromImage(a_bmp))
        {
            Rectangle rect;
            while (m_rects_to_update.TryDequeue(out rect))
            {
                Bitmap bmp = m_periodogram.GetBitmap(rect);
                g.DrawImage(bmp, rect);
            }
        }
    }

    public override RendererType RendererType
    {
        get
        {
            return RendererType.SamplerPeriodogram;
        }
    }
}

Jej zadaniem jest generować tylko punkty z powierzchni projekcji kamery, same promienie nie są śledzone. Same próbki są skalowane, by było widać odstępy między nimi. Skalowaniu podlegają też końcowe wartości kolorów na periodogramie, chodzi o to by nie był on zbyt jasny ani zbyt ciemny.

Za samą generację DFFT2 odpowiada postprocesor.

Dla mnie najważniejsze było potwierdzenie tego, że mój kod generuje rozkład Poissona, a nie coś innego. Moje periodogramy wyglądają tak jak inne, które można znaleźć w sieci, tak więc w kodzie raczej nie ma błędów. Same periodogramy wyglądają jak powinny, czyli w kodzie samplera Poisonna nie ma błędów.

W przyszłości zamierzam dodać samplowanie adaptatywne. Być może będę generować próbki tylko dla małego kwadratu i wykorzystywać je do całego renderowane obrazu (z jakimś mniej lub bardziej losowym mieszaniem). Obecnie renderowanie odbywa się kafelkami, zamierzam dla każdego kafelka osobno generować punkty do samplowania i zwalniać je jak nie są już potrzebne. Pozwoli to zaoszczędzić trochę pamięci, zwłaszcza przy samplowaniu adaptatywnym. Dla samplowania adaptatywnego będzie to chyba naturalna opcja dla zsamplowania jakiegoś obszaru z wyższą częstotliwością poprzez utworzenie nowej pary sampler/resampler. Z tym zwalnianiem to też nie taka prosta sprawa, dla Poisonna potrzebuje brzegowe próbki z sąsiednich kafelków, tak więc tak od razu nie możemy się pozbyć próbek z kafelków. Periodogram na pewno pozwoli mi stwierdzić, że wszystko działa.

Filtr Nearest Neighbour

Filtr ten nie w zasadzie filtrem. Dla wartości zero zwraca jeden, dla pozostałych zero. Ma charakter impulsu. Promień może być dowolnie mały ale większy od zera, ustawiłem go na 0.001. Ze względu na charakter filtra zablokowałem jego używanie dla resamplowania - nie ma to sensu. Możemy go tylko wykorzystać w postporcesingu podczas zmiany rozdzielczości.

Dla jasności przy pomniejszaniu z źródła branych jest dokładnie tyle pikseli ile ma obraz docelowy. Podczas powiększania zachowuje się on tak jak box.

public class NearestNeighbourFilter : Filter
{
    public override double Ray
    {
        get
        {
            return 0.001;
        }
    }

    public override double Evaluate(double a_value)
    {
        if (a_value.IsAlmostEquals(0))
            return 1; // TODO: mniejsze.
        else
            return 0;
    }
}

Ponieważ nie jest to filtr w pełnym tego sensie, do procedury obliczania wag musiałem dodać małą poprawkę, tak by dla tego filtra zawsze generował się impuls.

Zmiana systemu cachowania

Zadaniem cachowania jest przetrzymywanie w pamięci załadowanych zasobów tak długo jak są używane i przez pewien czas jak już nie są używane. Tak, że podczas ponownego żądania ich użycia nie musimy tracić czasu na ich załadowanie. Kiedy np. dwa materiały proszą o taką samą bitmapę z pliku to dostają zawsze taką samą z cachu. Dla bitmap ładowanych z pliku zakładamy, że nie można ich modyfikować.

Rozmiar pamięci cache jest określony sztywno z góry, jest to też rozmiar pamięci jaki możemy przeznaczyć na tekstury, mapy normalnych, mapowanie wypukłości itp. Kiedy kolejne ładowanie z pliku doprowadziło do przekroczenia rozmiaru cachu poszukujemy w nim nieużywanych wpisów i zwalniamy je. Jeśli mimo to przekroczyliśmy rozmiar cache raportowany jest błąd. Po zakończeniu renderowania obiekty zwalniają zasoby pobrane z cache i raportują CacheManager o zakończeniu ich używania. Jeśli od zasobu w cache odłączyli się wszyscy jego użytkownicy jest on wolny. Taki zasobów zwalniamy natychmiast kiedy brakuje nam pamięci, albo po określonym czasie nieużywania.

Próbowałem także zaimplementować rozwiązanie, które w takiej sytuacji wyładowało by zasób obecnie używany ale najdawniej użyty. Sprawa niby atrakcyjna okazała się kłopotliwa. Taki zasób trzeba synchronizować przy każdym użyciu. Każde pobranie koloru dla piksela. Nie da się tego przeskoczyć. Idealnie chcemy takiego efektu: wiele wątków renderujących uzyskuje dostęp równoległy do danych. Ale jeśli CacheManager oznaczył obiekt do zwolnienia to zakłada swoją blokadę i czeka aż wątki robocze go opuszczą. Wtedy go zwalnia. Jeśli w tym momencie wątek roboczy próbuje użyć zasobu musi on czekać na zwolnienie go przez CacheManager, zajmuje go swoją blokadą (dla wątków roboczych) i sprawdza (nie ma go, więc ponownie ładuje). Potrzebujemy dwóch obiektów do synchronizacji. Dałem sobie z tym spokój.

Wszystkie zasoby są ładowane w momencie startu renderowania.

Limit pamięci później będzie można rozbudować o bardziej mądre, agresywne zajmowanie wolnej pamięci.

Od razu odrzucam korzystanie z weak references, one się do tego nie nadają. Po przejściu obiekt w stan niewykorzystywania (tylko słaba referencja na niego wskazuje) jest on bardzo szybko zwalniany przez GC.

Nasz system cache musi mieć w osobnym wątku zaimplementowaną logikę zwalniającą pliki nieużywane przez dłuższy czas - niezależnie od tego czy pamięć jest, czy jej nie ma.

System cache ma mieć charakter globalny. Powinien być całkowicie oderwany od konkretnej sceny.

W przyszłości dodane zostaną mip-mapy. Warto to mieć na uwadze.

Obecny system cache zostanie przerobiony tak by spełniał powyższe wymagania.

Jeden obiekt może być związany tylko z jednym wpisem w cache tego samego typu. To ograniczenie wyszło podczas pisania i nie wydaje mi się, by była potrzeba jego rewidowania.

System cache w stosunku do poprzedniego stał się bardziej uniwersalny. Obiekt który chce być cachowany musi zaimplementować interfejs ICacheable.

Pewna dodatkowa logika była potrzebna ba zabezpieczenie się przed podwójnym ładowaniem tego samego zasobu. W takim wypadku jeden wątek ładuje zasób, drugi (reszta) czeka na zakończenie operacji. Sytuacja taka może się zdarzyć jeśli więcej wątków próbuje pobrać zasób jednocześnie.

Kod interfejsu ICacheable:

public interface ICacheable
{
    int MemoryUsage { get; }
    object Data { get; }
    string FileName { get; }
}

Kod menadżera pamięci tymczasowej CacheManager:

public class CacheManager
{
    #region CacheEntry
    private class CacheEntry
    {
        public object Data;
        public int MemoryUsage;
        public string FileName;
        public Type DataType;
        public DateTime ReleasedTime;
        public List<ICacheable> Users = new List<ICacheable>();

        public bool ContainsUser(ICacheable a_user)
        {
            for (int i = 0; i < Users.Count; i++)
            {
                if (Object.ReferenceEquals(Users[i], a_user))
                    return true;
            }

            return false;
        }
    }
    #endregion

    public static long CACHE_SIZE = 1024 * 1024 * 1024;
    public static TimeSpan NOT_USED_PERIOD = new TimeSpan(0, 0, seconds: 10);
    public static int SLEEP_PERIOD_MS = 1000;

    private static long s_memory_usage;

    private static List<CacheEntry> s_cache = new List<CacheEntry>();
    private static List<CacheEntry> s_loading = new List<CacheEntry>();

    private static Logger Logger = LogManager.GetLogger("Cache");

    static CacheManager()
    {
        Task.Factory.StartNew(() => 
        {
            for (; ; )
            {
                Thread.Sleep(SLEEP_PERIOD_MS);

                lock (s_cache)
                {
                    for (int i = s_cache.Count - 1; i >= 0; i--)
                    {
                        CacheEntry ce = s_cache[i];

                        if (ce.Users.Count != 0)
                            continue;

                        if (DateTime.Now - ce.ReleasedTime > NOT_USED_PERIOD)
                        {
                            s_cache.RemoveAt(i);
                            Interlocked.Add(ref s_memory_usage, -ce.MemoryUsage);

                            Logger.Info("Release not used data for: {0}", 
                                ce.FileName);

                            GC.Collect();
                        }
                    }
                }
            }
        });
    }

    public static long MemoryUsage
    {
        get
        {
            return s_memory_usage;
        }
    }

    private static CacheEntry FindEntry(ICacheable a_user)
    {
        lock (s_cache)
        {
            foreach (var ce in s_cache)
            {
                if (Object.ReferenceEquals(ce.Data, a_user.Data))
                {
                    Debug.Assert(ce.Data != null);
                    Debug.Assert(ce.ContainsUser(a_user));

                    return ce;
                }
            }

            return null;
        }
    }

    public static T GetOrLoad<T>(ICacheable a_user, Action a_loader) where T : class
    {
        CacheEntry nce = null;

        lock (s_cache)
        {
            foreach (var ce in s_cache)
            {
                if (ce.FileName != a_user.FileName)
                    continue;

                if (!(ce.Data is T))
                    continue;

                nce = ce;
                break;
            }

            if (nce == null)
            {
                lock (s_loading)
                {
                    foreach (var ce in s_loading)
                    {
                        if (ce.FileName != a_user.FileName)
                            continue;

                        if (!(ce.Data is T))
                            continue;

                        nce = ce;
                        break;
                    }


                    if (nce == null)
                    {
                        nce = new CacheEntry();

                        nce.DataType = typeof(T);
                        nce.FileName = a_user.FileName;

                        s_loading.Add(nce);
                    }
                }
            }
        }

        lock (nce)
        {
            if (nce.Data != null)
            {
                Logger.Info("Cached data finded: {0}", nce.FileName);
                Debug.Assert(!nce.ContainsUser(a_user));
                nce.Users.Add(a_user);

                return (T)nce.Data;
            }
            else
            {
                Logger.Info("Adding to cache: {0}", a_user.FileName);

                a_loader();

                nce.Data = a_user.Data;
                nce.MemoryUsage = a_user.MemoryUsage;

                Debug.Assert(!nce.ContainsUser(a_user));
                nce.Users.Add(a_user);

                Debug.Assert(nce.Data != null);
                Debug.Assert(!String.IsNullOrWhiteSpace(nce.FileName));
                Debug.Assert(nce.DataType == nce.Data.GetType());
                Debug.Assert(nce.MemoryUsage > 0);

                while (nce.MemoryUsage + MemoryUsage > CACHE_SIZE)
                {
                    Logger.Info("Memory usage exceeded: {0} > {1}", 
                        nce.MemoryUsage + MemoryUsage, CACHE_SIZE);

                    if (!RemoveEntry())
                        throw new OutOfMemoryException("Increase cache size");

                    if (MemoryUsage == 0)
                        break;
                }

                Interlocked.Add(ref s_memory_usage, nce.MemoryUsage);

                lock (s_cache)
                {
                    s_cache.Add(nce);
                }

                lock (s_loading)
                {
                    s_loading.Remove(nce);
                }

                return (T)nce.Data;
            }
        }
    }

    private static bool RemoveEntry()
    {
        CacheEntry ce = null;

        lock (s_cache)
        {
            DateTime last_usage_time = new DateTime(0);

            foreach (var c in s_cache)
            {
                if (c.Users.Count != 0)
                    continue;

                if (c.ReleasedTime >= last_usage_time)
                {
                    last_usage_time = c.ReleasedTime;
                    ce = c;
                }
            }
            
            if (ce == null)
                return false;

            Logger.Info("Releasing not used cache data: {0}", ce.FileName);

            s_cache.Remove(ce);
            Interlocked.Add(ref s_memory_usage, -ce.MemoryUsage);
            return true;
        }
    }

    public static void Unbind(ICacheable a_user)
    {
        Logger.Info("Unbinding: {0}", a_user.FileName);

        lock (s_cache)
        {
            CacheEntry ce = FindEntry(a_user);

            Debug.Assert(ce != null);
            Debug.Assert(ce.ContainsUser(a_user));

            for (int i = 0; i < ce.Users.Count; i++)
            {
                if (Object.ReferenceEquals(ce.Users[i], a_user))
                {
                    ce.Users.RemoveAt(i);
                    break;
                }
            }

            if (ce.Users.Count == 0)
                ce.ReleasedTime = DateTime.Now;
        }
    }
}

Wybrane fragmenty z obiektu implementującego interfejs ICacheable:

Filtr trójkątny

Filtr ten to nic innego jak interpolacja bilinearna - wartość pośrednia zmienia się liniowo. Pole pod powierzchnią filtra to 1. Zaś sam filtr jak nazwa wskazuje ma postać trójkąta o wysokości 1 i szerokości podstawy 2 (promień 1). Przykładowy kod w C#:

public class TriangleFilter : Filter
{
    public override double Ray
    {
        get
        {
            return 1;
        }
    }
 
    public override double Evaluate(double a_value)
    {
        if (a_value >= 1)
            return 0;
        if (a_value < -1)
            return 0;

        return 1 - Math.Abs(a_value);
    }
}

Filtr Gauss'a

Filtr Gaussa ma właściwości rozmywające.

Kształt filtru wywodzi się z funkcji na rozkład normalny:

$f(x)=e^{\displaystyle -\alpha x^2}$

Dla x dążącego do nieskończoności wartość filtru dąży do zera. W naszym przypadku wartość filtra poza promieniem powinna wynosić zera i najlepiej dla eliminacji oscylacji jak x dąży do granicy promienia wartość filtra powinna dążyć do zera. Najprościej możemy to zrobić przesuwając odpowiednio wartości filtru:

$f(x)=e^{\displaystyle -\alpha x^2} - e^{\displaystyle -\alpha r^2}$

Ważną jego cechą jest to, że jest on transformatą Fouriera samego siebie - kształt jest taki sam w czasie jak i w częstotliwości.

Kod w Matlabie potrzebny do generacji wykresów:

function [] = plot_gaussian()

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

    ym1 = gaussian(x, 0.5, 2);
    ym2 = gaussian(x, 1/sqrt(2), 2);
    ym3 = gaussian(x, 1, 2);
    ym4 = gaussian(x, 0.5, 4);

    trapz(x, ym1)
    trapz(x, ym2)
    trapz(x, ym3)
    trapz(x, ym4)

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    title('gaussian - time');
    plot(x, ym1, 'r');
    plot(x, ym2, 'g');
    plot(x, ym3, 'b');
    plot(x, ym4, 'c');
    legend('alpha = 0.5, ray = 2', 'alpha = 0.7, ray = 2', ... 
        'alpha = 1, ray = 2', 'alpha = 0.5, ray = 4');

    YM1 = fft(ym1, 100000);
    YM1 = [YM1(1:1:1000),YM1(end-1000:1:end)];

    YM2 = fft(ym2, 100000);
    YM2 = [YM2(1:1:1000),YM2(end-1000:1:end)];

    YM3 = fft(ym3, 100000);
    YM3 = [YM3(1:1:1000),YM3(end-1000:1:end)];

    YM4 = fft(ym4, 100000);
    YM4 = [YM4(1:1:1000),YM4(end-1000:1:end)];

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    title('gaussian - freq');
    plot(fftshift(abs(YM1)), 'r');
    plot(fftshift(abs(YM2)), 'g');
    plot(fftshift(abs(YM3)), 'b');
    plot(fftshift(abs(YM4)), 'c');
    legend('alpha = 0.5, ray = 2', 'alpha = 0.7, ray = 2', ... 
        'alpha = 1, ray = 2', 'alpha = 0.5, ray = 4');

end

function [y] = gaussian(x, alpha, ray)

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

    for i=1:size(y, 2)

        if (abs(x(i)) > ray)
            y(i) = 0;
        else
            y(i) = exp(-alpha * x(i) * x(i)) - exp(-alpha * ray * ray);
        end
    end
end

Przykładowy kod filtra w C#:

public class GaussianFilter : Filter
{
    private double m_exp;
    private double m_alpha = 0.5;

    [YAXNode]
    public double Alpha
    {
        get
        {
            return m_alpha;
        }
        set
        {
            m_alpha = value;
            m_exp = Math.Exp(-m_alpha * Ray * Ray);
        }
    }

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

    public override double Evaluate(double a_value)
    {
        if (a_value >= 2)
            return 0;
        if (a_value < -2)
            return 0;

        return Math.Exp(-m_alpha * a_value * a_value) - m_exp;
    }
}

Dla naszych zastosowań używanie innego większego promienia nie ma sensu. Filtr taki tylko coraz bardziej rozmywa obraz.

Filtr prostokątny

Chyba najprostszy filtr o ile go tak można nazwać. Uśrednia wartości próbek w obszarze splotu. Jego promień to 0.5. Co oznacza, że przy przekształcaniu obrazu 1:1 pod uwagę zostanie wzięta wartość tylko jednego piksela. W dziedzinie częstotliwości odpowiada mu funkcja sinc o widmie rozciągającym się do nieskończoności. Filtr taki generuje silny aliasing. Wykresem jest linia prosta o wartości 1 od -0.5 do 0.5. Pole powierzchni to 1. Filtr ten w żaden sposób nie wygładza pomniejszonego lub powiększonego obrazu. Przy pomniejszaniu pod uwagę brany jest pojedynczy piksel obrazu wynikowego. Przy powiększaniu brany jest jeden piksel obrazu źródłowego.

Przykładowy kod filtru w C#:

public class BoxFilter : Filter
{
    public override double Ray
    {
        get
        {
            return 0.5;
        }
    }
 
    public override double Evaluate(double a_value)
    {
        a_value = Math.Abs(a_value);

        if (a_value >= 0.5)
            return 0;
        if (a_value < -0.5) 
            return 0;

        return 1;
    }
}

Box w przeciwieństwie do innych filtrów nie ma na brzegach wartości zerowych. Stąd by uniknąć by promień przy samplowaniu 1x1 (szczególnie Jitter) nie był uwzględniany w dwóch sąsiednich pikselach punkt z jednego brzegu uwzględniamy, z drugiego nie. Lewy brzeg jest domknięty, prawy nie. Taki sposób domknięcia bardziej naturalnie współgra z zaokrąglaniem liczb. Dla samplowania Grid domknięcie przedziałów nie ma znaczenia bo zawsze pobieramy wartość ze środków.

Filtr Mitchell

Autorami pracy, która opisuje filtr są Mitchell i Netravali. Filtr zwyczajowo jednak określa się mianem filtru Mitchell-a. Równania filtru:

$ F(x) = \frac{1}{6} \begin{cases} \ (12-9B-6C)|x|^3 + (12B+6C-18)x^2 + (6-2B) & | x | < 1 \\ (-B-6C)|x|^3 + (6B+30C)x^2 + (-12B-48C)|x| + (8B+24C) & 1 \le |x|<2 \\ 0 & |x|>2 \end{cases}$

Filtr jest parametryzowny poprzez B i C. Zaleca się by B i C spełniało równanie: $B+2C=1$. Typowe wartości powszechnie używane to $B=^1/_3$ i $C=^1/_3$. Stanowią one dobry kompromis pomiędzy rozmyciem, a oscylacjami (ringing). Poniżej rysunek z oryginalnej pracy na której zaznaczono eksperymentalne wyniki subiektywnych odczuć grupy badanej.

Punktem wyjścia do genezy filtru było wyjście od generalnej postaci filtru bicubicznego, czyli złożenia dwóch wielomianów trzeciego stopnia dla zakresów wartości jak w definicji filtru. Każdy taki wielomian ma cztery współczynniki. Następnie zakładając, że w punkcie x=1 filtr powinien być ciągły, pochodna powinna być ciągła, pole powierzchni pod filtrem powinno być równe jeden $\int_{-2}^2F(x)d x=1$. Po nałożeniu tych warunków liczba parametrów redukuje się do dwóch: B i C.

Filtr ten obszar o ujemnych wartościach. Czyli posiada właściwości wyostrzające. W specyficznych warunkach wartość odfiltrowanego punktu może być ujemna.

Wygląd filtru dla kilku typowych wartości B i C:

Pole pod wykresem filtrów to 1. Dzięki czemu filtr ten nie zmienia energii sygnału.

Kod użyty do generacji tego wykresu:

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

ym1 = mitchell(1/3, 1/3);
ym2 = mitchell(1, 0);
ym3 = mitchell(0, 0.5);
ym4 = mitchell(3/2, 1/3);

trapz(x, ym1)
trapz(x, ym2)
trapz(x, ym3)
trapz(x, ym4)

figure(3);
hold on;
title('mitchell - spatial');
plot(x, ym1, 'r');
plot(x, ym2, 'g');
plot(x, ym3, 'b');
plot(x, ym4, 'c');
legend('B=1/3, C = 1/3', 'B=1, C=0', 'B=0, C=0.5', 'A=3/2, B=1/3');

function [y] = mitchell(B, C)

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

    for i=1:size(y, 2)

        xx = abs(x(i));
        x2 = xx * xx;

        if (xx < 1)

            y(i)= (12 - 9 * B - 6 * C) * x2 * xx + ...
                (12 * B + 6 * C - 18) * x2 + ...
                (6 - 2 * B);

        elseif (xx <= 2)

            y(i)= (-B - 6 * C) * x2 * xx + ...
                (6 * B + 30 * C) * x2 + ...
                (-12 * B - 48 * C) * xx + ...
                (8 * B + 24 * C);
        end

        y(i) = y(i) / 6;
    end
end

Przykładowy kod w C#:

public class MitchellFilter : Filter
{
    public double B = 1.0 / 3.0;
    public double C = 1.0 / 3.0;

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

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

        if (a_value >= 2)
            return 0;
        if (a_value < -2)
            return 0;

        double v2 = v * v;

        if (v > 1)
        {
            return (1.0 / 6.0) * (
                (-B - 6 * C) * v2 * v +
                (6 * B + 30 * C) * v2 +
                (-12 * B - 48 * C) * v +
                8 * B + 24 * C);
        }
        else
        {
            return (1.0 / 6.0) * (
                (12 - 9 * B - 6 * C) * v2 * v +
                (12 * B + 6 * C - 18) * v2 +
                6 - 2 * B);
        }
    }
}