2012-02-02

How to center control in panel

This solution doesn't depend on any WinForms features like TableLayout, AutoScroll, Anchoring, etc. It requires panel with one control. It reacts on resizing of control and panel resizing. It doesn't change control size. When there is no sufficient place for control then scroll bar appears. Control left, top borders never go outside left, top borders of panel.

public static void CenterControlInPanel(Panel a_panel)
{
    if (a_panel.Controls.Count != 1)
        throw new InvalidOperationException();

    Control control = a_panel.Controls[0];
            
    a_panel.AutoScroll = false;
    a_panel.AutoSize = false;

    HScrollBar horz_bar = new HScrollBar();
    VScrollBar vert_bar = new VScrollBar();

    a_panel.Controls.Add(horz_bar);
    a_panel.Controls.Add(vert_bar);

    horz_bar.Dock = DockStyle.Bottom;
    vert_bar.Dock = DockStyle.Right;

    Panel panel = new Panel();
    a_panel.Controls.Add(panel);
    panel.Dock = DockStyle.Fill;

    a_panel.Controls.Remove(control);
    panel.Controls.Add(control);

    Action on_scroll = () =>
    {
        if (panel.ClientSize.Width < control.Width)
        {
            if (horz_bar.Value > 0)
                control.Left = -horz_bar.Value;
            else
                control.Left = 0;
        }
        else
            control.Left = (panel.ClientSize.Width - control.Width) / 2;


        if (panel.ClientSize.Height < control.Height)
        {
            if (vert_bar.Value > 0)
                control.Top = -vert_bar.Value;
            else
                control.Top = 0;
        }
        else
            control.Top = (panel.ClientSize.Height - control.Height) / 2;

    };

    Action on_resized = () =>
    {
        if (panel.ClientRectangle.Width < control.Width)
            horz_bar.Visible = true;
        else
            horz_bar.Visible = false;

        if (panel.ClientRectangle.Height < control.Height)
            vert_bar.Visible = true;
        else
            vert_bar.Visible = false;

        if (panel.ClientRectangle.Width < control.Width)
            horz_bar.Visible = true;
        else
            horz_bar.Visible = false;

        vert_bar.Minimum = 0;
        vert_bar.Maximum = control.Height;
        vert_bar.LargeChange = panel.ClientRectangle.Height;

        horz_bar.Minimum = 0;
        horz_bar.Maximum = control.Width;
        horz_bar.LargeChange = panel.ClientRectangle.Width;

        if (horz_bar.LargeChange + horz_bar.Value - 1 > control.Width)
            horz_bar.Value = control.Width - horz_bar.LargeChange + 1;
        if (vert_bar.LargeChange + vert_bar.Value - 1 > control.Height)
            vert_bar.Value = control.Height - vert_bar.LargeChange + 1;

        on_scroll();
    };

    panel.Resize += (s, e) => on_resized();
    horz_bar.Scroll += (s, e) => on_scroll();
    vert_bar.Scroll += (s, e) => on_scroll();
    control.Resize += (s, e) => on_resized();

    on_resized();
}

2012-02-01

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++;
        }
    }
}

Superellipsoid

Superellipsoida jest bryłą zamkniętą o powierzchni, która nie zawsze jest ciągła (zależnie od parametrów).

Równanie Superellipsoidy w postatci uwikłanej:

$ \left( \left|x\right|^{R} + \left|y\right|^{R} \right)^{\frac{T}{R}} + \left|z\right|^{T} = 1$

, gdzie $T \in \mathbb{R}$ i $R \in \mathbb{R}$

Bryła jest zamknięta w sześcianie (-1, 1). Bryła jest symetryczna względem płaszczyzn układu współrzędnych. Punkt znajduje się we wnętrzu bryły jeśli:

$ \left( \left|x\right|^{R} + \left|y\right|^{R} \right)^{\frac{T}{R}} + \left|z\right|^{T} < 1$

Gradient z funkcji uwikłanej definiującej powierzchnię:

$\displaystyle\frac{df(x,y,z)}{dx} = T\cdot\left(\left|x\right|^R+\left|y\right|^R\right)^{\frac{T}{R}-1}\cdot sign(x) \cdot \left|x\right|^{R-1}$

$\displaystyle\frac{df(x,y,z)}{dx} = T\cdot\left(\left|x\right|^R+\left|y\right|^R\right)^{\frac{T}{R}-1}\cdot sign(x) \cdot \left|x\right|^{R-1}$

$\displaystyle\frac{df(x,y,z)}{dx} = T\cdot\left(\left|x\right|^R+\left|y\right|^R\right)^{\frac{T}{R}-1}\cdot sign(y) \cdot \left|y\right|^{R-1}$


$\displaystyle\frac{df(x,y,z)}{dx} = T \cdot sign(z) \cdot \left|z\right|^{T-1}$

Gradient dla punktu powierzchni bryły jest równy normalnej w tym punkcie.

Równanie parametryczne ma postać:

$\begin{align} x(u,v) &{}= A c\left(v,\frac{2}{t}\right) c\left(u,\frac{2}{r}\right) \\ y(u,v) &{}= B c\left(v,\frac{2}{t}\right) s\left(u,\frac{2}{r}\right) \\ z(u,v) &{}= C s\left(v,\frac{2}{t}\right) \\ & -\pi/2 \le v \le \pi/2, \quad -\pi \le u < \pi ,\end{align}$

,gdzie:

$\begin{align} c(\omega,m) &{}= \text{sign}(\cos \omega) |\cos \omega|^m \\ s(\omega,m) &{}= \text{sign}(\sin \omega) |\sin \omega|^m\end{align}$

Kod generujący z tego równania siatkę w Matlabie:

function [ ] = plot_se( )

R = 2.3;
T = 0.5;
scale = [1 1 1];

n = 200;
etamax = pi / 2;
etamin = -pi / 2;
wmax = pi;
wmin = -pi;
delta = (etamax - etamin) / n;
dw = (wmax - wmin) / n;
[i,j] = meshgrid(1:n+1, 1:n+1);
eta = etamin + (i-1) * delta;
w = wmin + (j-1) * dw;
x = scale(1) .* sign(cos(eta)) .* abs(cos(eta)).^(2/T) .* ...
sign(cos(w)) .* abs(cos(w)).^(2/R);
y = scale(2) .* sign(cos(eta)) .* abs(cos(eta)).^(2/T) .* ...
sign(sin(w)) .* abs(sin(w)).^(2/R);
z = scale(3) .* sign(sin(eta)) .* abs(sin(eta)).^(2/T);

z = z - 2;

meshc(x,y,z);
title(sprintf('R: %.1f, T: %.1f', R, T));

end


Przykłady:



Cały kod jest podobny do kodu dla SuperQuadric. Zmieniamy tylko zawartość funkcji Norm i Gradient i zmieniay kod tak by zamiast Exponent wykorzystywane był dwie zmienne T i R.

Kod:

public class SuperellipsoidObject : RenderableObject
{
    private RootFinder m_root_finder;
    private double m_R = 1;
    private double m_T = 1;
    private double m_TR = 1;

    private double Norm(Vector3 a_point)
    {
        return Math.Pow(Math.Pow(Math.Abs(a_point.X), R) + 
            Math.Pow(Math.Abs(a_point.Y), R), TR) + 
            Math.Pow(Math.Abs(a_point.Z), T);
    }

    [YAXNode]
    public double R
    {
        get
        {
            return m_R;
        }
        set
        {
            m_R = value;
            m_TR = T / R;
        }
    }

    [YAXNode]
    public double T
    {
        get
        {
            return m_T;
        }
        set
        {
            m_T = value;
            m_TR = T / R;
        }
    }

    private double TR
    {
        get
        {
            return m_TR;
        }
    }

    public SuperellipsoidObject(Vector3 a_right, Vector3 a_up) :
        base(a_right, a_up)
    {
        Name = "Superellipsoid";
        Closed = true;
        m_uv_mapper = new SphericalUVMapper();
        Update(UpdateFlags.All);
    }

    private Vector3 Gradient(Vector3 a_point)
    {
        double ax = Math.Abs(a_point.X);
        double ay = Math.Abs(a_point.Y);
        double xR1 = Math.Pow(ax, R - 1);
        double yR1 = Math.Pow(ay, R - 1);
        double xryr = Math.Pow(xR1 * ax + yR1 * ay, TR - 1);

        return new Vector3(
            xryr * Math.Sign(a_point.X) * T * xR1,
            xryr * Math.Sign(a_point.Y) * T * yR1, 
            Math.Sign(a_point.Z) * T * Math.Pow(Math.Abs(a_point.Z), T - 1));   
    }

    private double GetDistToNearestAxisOrOneOneBoxPlaneAndInsideOneOneBox(
        Vector3 a_ray_start, Vector3 a_ray_dir)
    {
        double dist_min = Double.PositiveInfinity;

        double dist = -(a_ray_start.Z / a_ray_dir.Z);
        if (dist &gt; 0)
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() &lt;= Vector3.ONE)
                dist_min = dist;
        }

        dist = -(a_ray_start.X / a_ray_dir.X);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() &lt;= Vector3.ONE)
                dist_min = dist;
        }

        dist = -(a_ray_start.Y / a_ray_dir.Y);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() &lt;= Vector3.ONE)
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.XAXIS, Vector3.XAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.XAXIS, -Vector3.XAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.YAXIS, Vector3.YAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.YAXIS, -Vector3.YAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.ZAXIS, Vector3.ZAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.ZAXIS, -Vector3.ZAXIS);
        if ((dist &gt; 0) &amp;&amp; (dist &lt; dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        return dist_min;
    }

    private double GetMinimalDistanceFromRayToSurface(Vector3 a_ray_start, 
        Vector3 a_ray_dir, double a_a, double a_b)
    {
        bool start_inside = Inside(a_ray_start);
        bool inside = start_inside;
        double min_t = Double.PositiveInfinity;

        for (int i = 0; i &lt; Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
        {
            if (Math.Abs(a_a - a_b) &lt; Constants.ROOT_FINDER_ABSOLUTE_ERROR)
                break;

            min_t = (a_a + a_b) / 2;
            Vector3 point = a_ray_start + a_ray_dir * min_t;
            Vector3 mid_grad = Gradient(point);

            double norm = Norm(point);

            if (start_inside)
            {
                if (norm &gt; 1)
                {
                    inside = false;
                    break;
                }
            }
            else
            {
                if (norm &lt; 1)
                {
                    inside = true;
                    break;
                }
            }

            double mid_dot = mid_grad * a_ray_dir;

            if (mid_dot &lt; 0)
                a_a = min_t;
            else
                a_b = min_t;
        }

        if (inside == start_inside)
            return Double.PositiveInfinity;

        return min_t;
    }

    private bool Inside(Vector3 a_local_point)
    {
        return !Norm(a_local_point).IsAlmostGreaterThen(1);
    }

    internal override Intersection GetIntersection(
        Intersection a_source_ray_intersection, Ray a_ray)
    {
        Vector3 local_dir;
        Vector3 local_start;
        double near_t;
        double far_t;

        TransformToLocal(a_ray, out local_dir, out local_start);

        if (!GetLocalBoundBox().GetIntersections(
            local_start, local_dir, out near_t, out far_t))
        {
            return Scene.NoIntersection;
        }

        Func<double, double=""> f = (t) =&gt; Norm(local_start + local_dir * t) - 1;

        bool back_hit = false;

        if (near_t &lt; 0)
        {
            if (Inside(local_start))
                back_hit = true;
            if (Norm(local_start).IsAlmostEquals(1))
            {
                Vector3 normal = -Gradient(local_start).Normalized;
                if (normal * local_dir &gt; 0)
                    back_hit = true;

                local_start = local_start + local_dir * Constants.MINIMAL_DISTANT;
                far_t -= Constants.MINIMAL_DISTANT;
                near_t -= Constants.MINIMAL_DISTANT;
            }
        }

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

        if (near_t &gt; 0)
        {
            local_start = local_start + local_dir * 
                (near_t - Constants.MINIMAL_DISTANT * 2);
            far_t = far_t - (near_t - Constants.MINIMAL_DISTANT * 2);
            near_t = Constants.MINIMAL_DISTANT * 2;
        }

        double local_dist = Double.PositiveInfinity;

        for (; ; )
        {
            double plane_t = GetDistToNearestAxisOrOneOneBoxPlaneAndInsideOneOneBox(
                local_start, local_dir);

            if (plane_t == Double.PositiveInfinity)
                break;

            double min_t = GetMinimalDistanceFromRayToSurface(
                local_start, local_dir, 0, plane_t);

            if (min_t != Double.PositiveInfinity)
            {
                local_dist = m_root_finder.FindRoot(0, min_t, f);

                if (local_dist == Double.PositiveInfinity)
                    local_dist = m_root_finder.FindRoot(min_t, plane_t, f);
            }
            else
                local_dist = m_root_finder.FindRoot(0, plane_t, f);

            if (local_dist != Double.PositiveInfinity)
                break;

            local_start = local_start + local_dir * plane_t;
            far_t -= plane_t;
        }

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

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

        Vector3 local_norm = -Gradient(local_pos).Normalized;

        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 &gt;= 0)
            return Scene.NoIntersection;

        return intersection;
    }

    protected override AABB GetLocalBoundBox()
    {
        return new AABB(-Vector3.ONE, Vector3.ONE);
    }

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

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

    public override string ToString()
    {
        return String.Format("Superellipsoid: {0}, R: {1}, T: {2}", Name, R, T);
    }

    internal override void GetTangents(Intersection a_intersection, 
        out Vector3 a_tangent_x, out Vector3 a_tangent_y)
    {
        a_tangent_x = Vector3.CrossProduct(Up, a_intersection.Pos).Normalized;

        if (a_intersection.BackHit)
        {
            a_tangent_y = Vector3.CrossProduct(a_intersection.Normal, 
                a_tangent_x).Normalized;
        }
        else
        {
            a_tangent_y = Vector3.CrossProduct(a_tangent_x, 
                a_intersection.Normal).Normalized;
        }
    }

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

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            m_root_finder = RootFinder.Create(
                Scene.RenderOptions.RootFindingMethod);
        }
    }

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

        base.ScaleAbsolute(a_scale);
    }
}

Przykłady:

R=2.3, T=0.5


R=0.5, T=0.5

2012-01-30

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.

Przybliżenie wartości pochodnej

Dla niektórych metod poszukiwania miejsc zerowych potrzebujemy wartości pochodnych, z reguły 1, 2 i 3 rzędu. Samą badaną funkcję nie zawsze mamy daną wzorem analitycznym by te pochodne policzyć, czasami jest to zbyt kłopotliwe. Wtedy nie pozostaje nam nic innego jak przybliżyć wartość pochodnej korzystając z wzoru:

$f'(x_{0}) = \displaystyle \lim _{h\to 0}\frac{f(x_{0}+h) - f(x_{0})}{h}$

Wartość drugiej otrzymamy poprzez podstawienie do podanego wzoru wartości pierwszej pochodnej.

Możemy wyróżnić trzy metody doboru przedziału. Centralnie wokół punktu x, na prawo od punktu x i na lewo od punktu x:

$\delta_h[f](x) = f(x+\tfrac12h)-f(x-\tfrac12h)$

$\nabla_h[f](x) = f(x) - f(x-h)$

$\delta_h[f](x) = f(x+\tfrac12h)-f(x-\tfrac12h)$

Ja skorzystałem z metody pierwszej. Kod na pierwszą i drugą pochodną jest następujący:

public class CentralDifferenceMethod
{
    private double m_h;
    private Func<double, double> m_func;

    public CentralDifferenceMethod()
        : base()
    {
    }
        
    public CentralDifferenceMethod(Func<double, double> a_func, 
        double a_h = Constants.DOUBLE_PRECISION * 2)
    {
        m_h = a_h;
        m_func = a_func;
    }

    public double Derive1(double a_x)
    {
        return (m_func(a_x + m_h) - m_func(a_x - m_h)) / (2 * m_h);
    }

    public double Derive2(double a_x)
    {
        return (m_func(a_x + m_h) - 2 * m_func(a_x) + 
            m_func(a_x - m_h)) / (m_h * m_h);
    }

    public static  Func<double, double> Derive1(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive2(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, 
        Func<double, double> a_d1, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_d1, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }
}

Zauważmy, że pierwsza pochodna wymaga dwóch punktów, druga trzech, itd.

Superquadric

Superquadric to rodzaj bryły, której powierzchnia nie jest gładka, zamkniętej w przedziale (-1,1). Równanie jej powierzchni:

$\left|x\right|^{\displaystyle r }+ \left|y\right|^{\displaystyle s} + \left|z\right|^{\displaystyle t}=1$

Wektor (r,s,t) definiuje nam kształt bryły. Generalnie eksponenta < 1 to wklęsłość, e = 1 to linia prosta. e = 2 to sfera, dla dużych e bryła przybliża się do sześcianu. Jeśli $ \left|x\right|^r + \left|y\right|^s + \left|z\right|^t < 1$ to punkt znajduje się wewnątrz bryły. Bryła jest zawsze symetryczna wobec płaszczyzn układu współrzędnych. Równanie parametryczne powierzchni:

$ \begin{align} x(u,v) &{}= A c\left(v,\frac{2}{r}\right) c\left(u,\frac{2}{r}\right) \\ y(u,v) &{}= B c\left(v,\frac{2}{s}\right) s\left(u,\frac{2}{s}\right) \\ z(u,v) &{}= C s\left(v,\frac{2}{t}\right) \\ & -\frac{\pi}{2} \le v \le \frac{\pi}{2}, \quad -\pi \le u < \pi ,\end{align} $

, gdzie:

$\begin{align} c(\omega,m) &{}= \text{sign}(\cos \omega) |\cos \omega|^m \\ s(\omega,m) &{}= \text{sign}(\sin \omega) |\sin \omega|^m\end{align}$

Kod Matlaba generujący siatkę z równania parametrycznego:
function [ ] = plot_sq(  )

  exponent = [ 0.5, 4, 1]; 
  scale = [1 1 1]; 

  n = 50;
  etamax = pi / 2;
  etamin = -pi / 2;
  wmax = pi;
  wmin = -pi;
  delta = (etamax - etamin) / n;
  dw = (wmax - wmin) / n;
  [i,j] = meshgrid(1:n+1, 1:n+1);
  eta = etamin + (i-1) * delta;
  w = wmin + (j-1) * dw;
  x = scale(1) .* sign(cos(eta)) .* abs(cos(eta)).^(2/exponent(1)) .* ... 
      sign(cos(w)) .* abs(cos(w)).^(2/exponent(1));
  y = scale(2) .* sign(cos(eta)) .* abs(cos(eta)).^(2/exponent(2)) .* ... 
      sign(sin(w)) .* abs(sin(w)).^(2/exponent(2));
  z = scale(3) .* sign(sin(eta)) .* abs(sin(eta)).^(2/exponent(3));
  
  z = z - 2;
  
  meshc(x,y,z);
 
end

Przykłady:

(0.5, 4, 1)


(0.5, 400, 0.5)


(1, 1, 0.5)


Najwięcej czasu zajęło mi napisanie prawidłowego kodu na intersekcję. Po pierwsze korzystamy z faktu, że nasza płaszczyzny układu współrzędnych oraz płaszczyzny BoundBox-a dzielą bryłę na 8 ćwiartek. Dla każdej ćwiartki znajdujemy miejsce zerowe osobno. W razie jego braku przemieszczamy się do następnej ćwiartki. Robimy tak długo aż opuścimy BoundBox (3 razy maksymalnie). Tym sposobem ograniczamy możliwość znalezienia nie tego miejsca zerowego co trzeba. Ale ciągle jeśli np. ćwiartka jest wypukła lub wkłęsła promień może ją przeciąć w dwóch miejscach i RootFinder, który znajduje jedno z miejsc zerowych i nie jest powiedziane, które to będzie, znajdzie to drugie, dalsze. Musimy znowu podzielić naszą przestrzeń poszukiwań. Służy do tego funkcja GetMinimalDistanceFromRayToSurface. Jest to trochę zmodyfikowany BisectionRootFinder. Poszukujemy punktu w którym promień jest najbliżej powierzchni bryły. Poszukujemy w niej minimum iloczynu skalarnego normalnej Superquadric (liczonej poza powierzchnią bryły wzdłuż promienia) i kierunku promienia. W miejscu gdzie te dwa wektory są prostopadłe, czyli czubku wypukłości iloczyn ten wynosi zero. Równocześnie na bieżąco w każdej iteracji sprawdzamy, czy czasem nie przeszliśmy przez powierzchnię bryły, co może się zdarzyć jako że zmierzamy do punktu wypukłości. Jeśli tak przerywamy iterację, gdyż mamy zakres do poszukiwań w którym promień przechodzi przez bryłę. No i ostatnią modyfikacją by uniknąć detekcji uderzenia w punkt wyjścia było przesunięcie minimalnego dystansu z zero na Constants.MINIMAL_DISTANT.

Gradient dla punktu poza powierzchnią bryły to taki wektor, że prosta będąca jego przedłużeniem przetnie powierzchnię pod kątem prostym. Wektor gradientu jest w punkcie przecięcia prostej normalną powierzchni.

Jak na mój gust cały kod chyba nie robi błędów. Jak promień nie przecinałby ćwiartki odległość promienia od powierzchni będzie funkcją z jednym ekstremum. Tak mi się wydaje....

Kod klasy:

public class SuperquadricObject : RenderableObject
{
    private RootFinder m_root_finder;

    [YAXNode]
    public Vector3 Exponents = new Vector3(1, 1, 1);

    public SuperquadricObject(Vector3 a_right, Vector3 a_up) :
        base(a_right, a_up)
    {
        Name = "Superquadric";
        Closed = true;
        m_uv_mapper = new SphericalUVMapper();
        Update(UpdateFlags.All);
    }

    private double Norm(Vector3 a_point)
    {
        return Math.Pow(Math.Abs(a_point.X), Exponents.X) +
                Math.Pow(Math.Abs(a_point.Y), Exponents.Y) +
                Math.Pow(Math.Abs(a_point.Z), Exponents.Z);
    }

    private Vector3 Gradient(Vector3 a_point)
    {
        return new Vector3(
            Exponents.X * Math.Sign(a_point.X) * 
                Math.Pow(Math.Abs(a_point.X), Exponents.X - 1),
            Exponents.Y * Math.Sign(a_point.Y) * 
                Math.Pow(Math.Abs(a_point.Y), Exponents.Y - 1),
            Exponents.Z * Math.Sign(a_point.Z) * 
                Math.Pow(Math.Abs(a_point.Z), Exponents.Z - 1));
    }

    private double GetDistToNearestAxisOrOneOneBoxPlaneAndInsideOneOneBox(
        Vector3 a_ray_start, Vector3 a_ray_dir)
    {
        double dist_min = Double.PositiveInfinity;

        double dist = -(a_ray_start.Z / a_ray_dir.Z);
        if (dist > 0)
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() <= Vector3.ONE)
                dist_min = dist;
        }

        dist = -(a_ray_start.X / a_ray_dir.X);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() <= Vector3.ONE)
                dist_min = dist;
        }

        dist = -(a_ray_start.Y / a_ray_dir.Y);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs() <= Vector3.ONE)
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.XAXIS, Vector3.XAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.XAXIS, -Vector3.XAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.YAXIS, Vector3.YAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.YAXIS, -Vector3.YAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.ZAXIS, Vector3.ZAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        dist = Plane.GetDist(a_ray_start, a_ray_dir, Vector3.ZAXIS, -Vector3.ZAXIS);
        if ((dist > 0) && (dist < dist_min))
        {
            Vector3 p = (a_ray_start + a_ray_dir * dist);
            if (p.Abs().IsAlmostLessThen(Vector3.ONE))
                dist_min = dist;
        }

        return dist_min;
    }

    private double GetMinimalDistanceFromRayToSurface(
        Vector3 a_ray_start, Vector3 a_ray_dir,
        double a_a, double a_b)
    {
        bool start_inside = Inside(a_ray_start);
        bool inside = start_inside;
        double min_t = Double.PositiveInfinity;

        for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
        {
            if (Math.Abs(a_a - a_b) < Constants.ROOT_FINDER_ABSOLUTE_ERROR)
                break;

            min_t = (a_a + a_b) / 2;
            Vector3 point = a_ray_start + a_ray_dir * min_t;
            Vector3 mid_grad = Gradient(point);

            // Norm(point)
            double norm = (mid_grad * (point / Exponents)); 

            if (start_inside)
            {
                if (norm > 1)
                {
                    inside = false;
                    break;
                }
            }
            else
            {
                if (norm < 1)
                {
                    inside = true;
                    break;
                }
            }

            double mid_dot = mid_grad * a_ray_dir;

            if (mid_dot < 0)
                a_a = min_t;
            else
                a_b = min_t;
        }

        if (inside == start_inside)
            return Double.PositiveInfinity;

        return min_t;
    }

    private bool Inside(Vector3 a_local_point)
    {
        return !Norm(a_local_point).IsAlmostGreaterThen(1);
    }

    internal override Intersection GetIntersection(
        Intersection a_source_ray_intersection, Ray a_ray)
    {
        Vector3 local_dir;
        Vector3 local_start;
        double near_t;
        double far_t;

        TransformToLocal(a_ray, out local_dir, out local_start);

        if (!GetLocalBoundBox().GetIntersections(local_start, 
            local_dir, out near_t, out far_t))
        {
            return Scene.NoIntersection;
        }

        Func<double, double> f = (t) => Norm(local_start + local_dir * t) - 1;

        bool back_hit = false;

        if (near_t < 0)
        {
            if (Inside(local_start))
                back_hit = true;
            if (Norm(local_start).IsAlmostEquals(1))
            {
                Vector3 normal = -Gradient(local_start).Normalized;
                if (normal * local_dir > 0)
                    back_hit = true;

                local_start = local_start + local_dir * Constants.MINIMAL_DISTANT;
                far_t -= Constants.MINIMAL_DISTANT;
                near_t -= Constants.MINIMAL_DISTANT;
            }
        }

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

        if (near_t > 0)
        {
            local_start = local_start + local_dir * 
                (near_t - Constants.MINIMAL_DISTANT * 2);
            far_t = far_t - (near_t - Constants.MINIMAL_DISTANT * 2);
            near_t = Constants.MINIMAL_DISTANT * 2;
        }

        double local_dist = Double.PositiveInfinity;

        for (; ; )
        {
            double plane_t = 
                GetDistToNearestAxisOrOneOneBoxPlaneAndInsideOneOneBox(
                    local_start, local_dir);

            if (plane_t == Double.PositiveInfinity)
                break;

            double min_t = 
                GetMinimalDistanceFromRayToSurface(
                    local_start, local_dir, 0, plane_t);

            if (min_t != Double.PositiveInfinity)
            {
                local_dist = m_root_finder.FindRoot(0, min_t, f);

                if (local_dist == Double.PositiveInfinity)
                    local_dist = m_root_finder.FindRoot(min_t, plane_t, f);
            }
            else
                local_dist = m_root_finder.FindRoot(0, plane_t, f);

            if (local_dist != Double.PositiveInfinity)
                break;

            local_start = local_start + local_dir * plane_t;
            far_t -= plane_t;
        }

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

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

        Vector3 local_norm = -Gradient(local_pos).Normalized;

        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;

        return intersection;
    }

    protected override AABB GetLocalBoundBox()
    {
        return new AABB(-Vector3.ONE, Vector3.ONE);
    }

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

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

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

    internal override void GetTangents(Intersection a_intersection, 
        out Vector3 a_tangent_x, out Vector3 a_tangent_y)
    {
        a_tangent_x = Vector3.CrossProduct(Up, a_intersection.Pos).Normalized;

        if (a_intersection.BackHit)
        {
            a_tangent_y = Vector3.CrossProduct(a_intersection.Normal, 
                a_tangent_x).Normalized;
        }
        else
        {
            a_tangent_y = Vector3.CrossProduct(a_tangent_x, 
                a_intersection.Normal).Normalized;
        }
    }

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

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            m_root_finder = RootFinder.Create(
                Scene.RenderOptions.RootFindingMethod);
        }
    }

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

        base.ScaleAbsolute(a_scale);
    }
}

Dla przyspieszenia obliczeń, zrezygnowałem z detekcji intersekcji bryły z AABB w współrzędnych świata, zamiast tego robię to z AABB w układzie lokalnym, który dużo lepiej opina bryłę. Cały kod przyspiesza o jakieś 25% zależnie od położenia Superquadric.

Normalna liczona jest z definicji gradientu.

Mam nadzieję, że reszta kodu jest jasna. Także dla mnie. Co okaże się za jakieś pół roku, gdy trzeba będzie do tego wrócić.

Przykłady wygenerowane przy użyciu powyższego kodu:

(4.4, 4.2, 0.5)


(0.5, 0.5, 0.5)


(2.3, 0.5, 1)