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));
endPrzykł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 > 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);
double norm = Norm(point);
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("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
Brak komentarzy:
Prześlij komentarz