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)