W przypadku niektórych kształtów likwidowałem parametry funkcji, gdyż odpowiadały one tak naprawdę za skalowanie. Czasami też przesuwałem powierzchnię (manipulując x, y, z) by znalazła ona się w środku układu współrzędnych oraz by zawierała się w przedziale od -1 do 1.
Na ten moment sprawdziłem renderowanie prawie wszystkich powierzchni zamkniętych podanych na MathWorld. I wystarczy tego. Cel jakim było w ogóle zrenderowanie takich powierzchni został osiągnięty. Na jego potrzeby musiałem napisać całkiem spory kawałek kodu, który bierze wielomian w postaci tekstu, parsuje do tablicy wyrażeń, buduje z niej RPN, z RPN drzewo. Drzewo wyrażeń można zoptymalizować - wymnażanie przez stała, dodawanie, mnożenie wyrazów, uproszczenia związane z 0 i 1. Oprócz tego wyrażenie można rozszerzyć, czyli wymnożyć, wypotęgować co się da - zbudować jedną wielką sumę wyrazów. Dzielenie przez inne wielomiany, czyli funkcje wymierne nie jest obsługiwane, podobnie jak wykładniki potęgowania nie będące liczbami. Dla zmiennych wykładniki te powinny być dodatnimi liczbami całkowitymi. Brakuje kodu na redukcję wyrażeń, czyli czegoś odwrotnego do rozszerzania. Gradient jest liczony wprost z funkcji bez jej rozszerzania. Dla wyliczenia intersekcji podstawiamy do wielomianu równanie kierunkowe linii, rozszerzamy je i grupujemy po stopniu t. Każda taka grupa to współczynnik dla wielomianu intersekcji. Następnie kompilujemy kod który zwraca nam taki wielomian z już wyliczonymi współczynnikami po podaniu mu jako parametrów kierunku i punktu startu promienia. Kompilacja następuje poprzez zbudowanie LINQ Expression. Kod wynikowy jest bardzo szybki.
Szczegółów kodu nie zamierzam podawać. Ale kilka rad po napisaniu czegoś takiego mam. Drzewo wyrażenia powinno być niezmienne. Totalnie niezmienne drzewo będzie wymagało skopiowania całego drzewa przy dowolnej zmianie, a będziemy to robić bardzo często. Najlepiej kod drzewo umieścić w osobnej przestrzeni nazw i poukrywać wszystkie do manipulowania jego strukturą bezpośrednio (rodzic, dzieci) i na zewnątrz wystawić tylko zbiór ustalonych metod. Dlatego najlepiej metody zmieniające drzewo jakoś zaizolować, a na zewnątrz wystawić tylko metody nie modyfikujące drzewa, ale zwracające nowe zmienione drzewo. Co do elementów drzewa to powinno się ono składać tylko z symboli, dodawania, mnożenia, dzielenia, potęgowania. Unarne minusy i odejmowanie powinniśmy zastąpić dodawaniem i mnożeniem przez -1. Dzielenie także jak najszybciej konwertować na mnożenie (obsługujemy dzielenie tylko przez stałe). Przez symbol rozumiemy wyrażenie A*N^B, gdzie A to stała, N to nazwa symbolu, B to jego potęga. A=0 to zero, B=0 to liczby. Dodatkowo warto od początku założyć, że dodawanie i mnożenie może mieć więcej dzieci, dzięki czemu znajdowanie grup typu ABCD+EFGH będzie szybkie i proste. Dzięki temu, że potęgowanie mamy od razu w symbolu nie musimy przeszukiwać drzewa nadmiernie. A stała w symbolu przez którą mnożymy sprawi, że często nie będziemy musieli tworzyć nowych symboli liczb. Na samym początku redukujemy liczby i przenosimy potęgowanie do symboli tak długo jak pozbędziemy się potęgowania i dzielenia. Pozbywamy się także minusów i odejmowania. Spłaszczamy drzewo tak by dodawania i mnożenie nie posiadało innych dzieci - odpowiednio dodawania i mnożenia i w takiej postaci staramy się utrzymywać drzewo.
Sam kod intersekcji nie różni się zbytnio od kodu intersekcji dla torusa. Z tym, że inaczej jest budowany wielomian intersekcji.
Można oczywiście policzyć gradient i współczynniki wielomianu intersekcji na piechotę. Jednakże jest to bardzo pracochłonne. Poza tym mało uniwersalne - mając cały kod przetwarzający możemy implementować dowolne powierzchnie dane wielomianem.
W moim przypadku największe problemy sprawił mi kształt serca: $(x^2+\frac{9}{4} y^2+z^2-1)^3-x^2z^3-\frac{9}{80}y^2z^3$
Serce bez AABB:
Serce z AABB - start promienia jest przybliżany do AABB i dla takiego punktu startu budowany jest wielomian intersekcji.
Serce z AABB - likwidacja miejsc zerowych.
Sprawdziłem jak wyglądają współczynniki i same miejsca zerowego dla promienia odbitego od błędnych punktów. Okazało się, że następuje ponowna detekcja uderzenia w powierzchnie od której się odbijamy. Możemy powiększyć
Constants.MINIMIAL_DISTANS = 1e-8 - miejsce od którego szukamy pierwiastków. Ale w tym przypadku ponieważ pierwszy współczynnik wielomianu intersekcji okazał się mniejszy od dokładności Constants.DOUBLE_PRECISION = 1e-12 postanowiłem napisać kod eliminujący takie współczynniki, przesuwający pozostałe współczynniki i zmiejszający stopień wielomianu. Tak mały współczynnik może odpowiadać za miejsce zerowe w granicach 1e-5. Taki zabieg nigdzie nie wyprodukował dodatkowych błędów. Błędy nie znikły całkowicie, ale są już prawie niewidoczne, a na pewno rozmyją przy supersamplowaniu.Serce z przesuwaniem punktu startu wzdłuż normalnej.
Serce ciągle ma dwa piksele za ciemne. Nie pozostało mi nic innego jak odsunąć punkt startu od powierzchni wzdłuż normalnej. Wartość przy której znikły ostatnie błędy to
Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA = 1e-5. Jest to bardzo dużo.Przykład powierzchni jabłka z refrakcją:
Ponieważ dysponuję kodem na obiekt torusa liczonego bezpośrednio i teraz mam do dyspozycji torus dany przez wielomian mogę porównać ich jakość i szybkość. Implementacja oparta na wielomianie jest ponad dwa razy wolniejsza. Przede wszystkim największy współczynnik wielomianu, który po uproszczeniu jeden jedynką tutaj jest wyliczany. Prawdopodobnie pozostałe współczynniki też nie są w optymalnej postaci. Coś za coś. Np. dla serca liczba współczynników to 217 po pełnym rozwinięciu. To daje obraz zabawy podczas wyliczania współczynników wielomianu.
Kod obiektu sceny:
public class PolynomialSurface : RenderableObject
{
protected PolynomialFunction m_polynomial;
private Func<Vector3, Vector3> m_gradient;
private Func<Vector3, Vector3, Polynomial> m_intersection;
private PolynomialSolver m_solver;
[YAXNode]
public string Equation;
[YAXNode]
public string Parameter1 = "";
[YAXNode]
public string Parameter2 = "";
[YAXNode]
public string Parameter3 = "";
[YAXNode]
public AABB LocalBoundBox;
private Vector3 m_aabb_minimum;
private Vector3 m_aabb_maximum;
public PolynomialSurface(Vector3 a_right, Vector3 a_up) :
base(a_right, a_up)
{
Name = "Polynomial surface";
m_aabb_minimum = Vector3.MAXIMUM;
m_aabb_maximum = Vector3.MINIMUM;
Closed = false;
LocalBoundBox = AABB.INFINITY;
m_uv_mapper = new SphericalUVMapper();
Update(UpdateFlags.All);
}
public override Intersection GetIntersection(
Intersection a_source_ray_intersection, Ray a_ray)
{
Vector3 local_dir;
Vector3 local_start;
bool back_hit = false;
if (!TransformToLocalToBoundBox(a_ray, out local_dir, out local_start))
return Scene.NoIntersection;
if ((a_source_ray_intersection != null) &&
(a_source_ray_intersection.SceneObject == this))
{
if (a_ray.RaySurfaceSide == RaySurfaceSide.RefractedSide)
{
local_start += ((WorldToLocalNormal *
a_source_ray_intersection.Normal).Normalized) *
-Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
}
else
{
local_start += ((WorldToLocalNormal *
a_source_ray_intersection.Normal).Normalized) *
Constants.POLYNOMIAL_SURFACE_FROM_MOVE_DELTA;
}
}
Polynomial p = m_intersection(local_start, local_dir);
double dist = m_solver.Solve(p);
if (dist == Double.PositiveInfinity)
return Scene.NoIntersection;
Vector3 local_pos = local_start + local_dir * dist;
Vector3 normal = m_gradient(local_pos);
back_hit = normal * local_dir > 0;
if (back_hit && OneSide)
return Scene.NoIntersection;
Vector3 world_pos = LocalToWorld * local_pos;
double world_dist = (world_pos - a_ray.Start).Length;
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;
m_aabb_minimum = Vector3.Minimize(local_pos, m_aabb_minimum);
m_aabb_maximum = Vector3.Maximize(local_pos, m_aabb_maximum);
return intersection;
}
public override Vector3 GetNormal(Intersection a_intersection)
{
if (a_intersection.BackHit)
{
return (LocalToWorldNormal *
-m_gradient(a_intersection.LocalPos)).Normalized;
}
else
{
return (LocalToWorldNormal *
m_gradient(a_intersection.LocalPos)).Normalized;
}
}
public override Vector2 GetUV(Intersection a_intersection)
{
return UVMapper.Map(a_intersection.UVW);
}
public override Vector3 GetUVW(Intersection a_intersection)
{
return base.GetUVW(a_intersection) * 0.5 + Vector3.HALF;
}
public override string ToString()
{
return String.Format("PolynomialSurface: {0}", Name);
}
public override void GetTangents(Intersection a_intersection,
out Vector3 a_tangent_x, out Vector3 a_tangent_y)
{
if (a_intersection.BackHit)
{
a_tangent_x = Vector3.CrossProduct(
Up, a_intersection.Normal).Normalized;
a_tangent_y = Vector3.CrossProduct(
a_tangent_x, a_intersection.Normal).Normalized;
}
else
{
a_tangent_x = Vector3.CrossProduct(
a_intersection.Normal, Up).Normalized;
a_tangent_y = Vector3.CrossProduct(
a_intersection.Normal, a_tangent_x).Normalized;
}
}
public override Scene Scene
{
get
{
return base.Scene;
}
set
{
if (Scene != null)
Scene.RenderEnd -= RenderEnd;
base.Scene = value;
if (Scene != null)
Scene.RenderEnd += RenderEnd;
}
}
private void RenderEnd(bool a_all)
{
if (!a_all)
return;
if (LocalBoundBox != AABB.INFINITY)
{
if ((!LocalBoundBox.Minimum.IsAlmostRelativeEquals(m_aabb_minimum)) &&
(!LocalBoundBox.Maximum.IsAlmostRelativeEquals(m_aabb_maximum)))
{
LogManager.GetCurrentClassLogger().Error("Bound box for {0}", Name);
LogManager.GetCurrentClassLogger().Error("Actually: {0}",
LocalBoundBox);
LogManager.GetCurrentClassLogger().Error("Should be: {0}",
new AABB(m_aabb_minimum, m_aabb_maximum));
}
}
}
protected override AABB GetLocalBoundBox()
{
return LocalBoundBox;
}
protected override void RenderStart(RenderStartPhase a_phase)
{
base.RenderStart(a_phase);
if (a_phase == RenderStartPhase.PrepareObjectToRender)
{
string equation = ApplyParameters(Equation);
m_polynomial = new PolynomialFunction(equation);
m_gradient = m_polynomial.Gradient;
m_intersection = m_polynomial.Intersection;
m_solver = new PolynomialSolver(PolynomialSolveMethod.Sturm,
Scene.RenderOptions.RootFindingMethod);
m_aabb_minimum = Vector3.MAXIMUM;
m_aabb_maximum = Vector3.MINIMUM;
}
}
private string ApplyParameters(string a_equation)
{
a_equation = ApplyParameter(a_equation, Parameter1);
a_equation = ApplyParameter(a_equation, Parameter2);
a_equation = ApplyParameter(a_equation, Parameter3);
return a_equation;
}
private string ApplyParameter(string a_equation, string a_parameter)
{
if (String.IsNullOrWhiteSpace(a_parameter))
return a_equation;
string[] ar = a_parameter.Split('=');
if (ar.Length < 2)
return a_equation;
a_equation = a_equation.Replace(ar[0], "(" + ar[1].Replace(",", ".") + ")");
return a_equation;
}
public override void ScaleAbsolute(double a_scale)
{
Scale *= a_scale;
base.ScaleAbsolute(a_scale);
}
}
Brak komentarzy:
Prześlij komentarz