Any good SDF distance estimators for a paraboloid
category: code [glöplog]
For no specific reason at all (ahem), I would like to have a true paraboloid SDF in an upcoming raytrace scene. However, I'm having a bit of trouble building a good distance estimator, since the naive one (distance from the point to the paraboloid surface along paraboloid axis) tends to overshoot for obvious reasons.
I made an attempt to calculate an analytical distance. If we assume the paraboloid is in point POS, and its axis points to a direction DIR, then define D_ORTHO as the distance orthogonal to the paraboloid axis and D_AX as the distance along the paraboloid axis, we can walk the paraboloid surface along the orthogonal direction and form the actual distance function quite easily. I'm getting something like DIST = sqrt((X^2-D_AX)^2 + (-X + D_ORTHO)^2), where X is the distance we walk into the axial direction from the paraboloid axis.
Then, derivating and solving for the roots of X to find the minimum distance, I'm getting 4X^3 + (-4D_AX+2)X - 2D_ORTHO = 0, which is irritatingly a third degree polynomial. The solution seems to be correct, however, at least on cursory inspection and using a random cubic polynomial solver from Wolfram.
Example with solution: https://i.imgur.com/Xd23YjI.png
Random cubic solver from Wolfram: http://www.wolframalpha.com/widgets/view.jsp?id=3f4366aeb9c157cf9a30c90693eafc55
In the above image, I simply have an (x^2) paraboloid with a (2, -1) point that I'm solving for for demonstration purposes.
Now, a cubic solver is general knowledge, and one can of course be implemented, but it contains a ton of cubic roots so it isn't going to be particularly fast, and there are other irritating side effects like having to check which of the solutions will be real in the first place. Before I delve into this (which I will do anyway tomorrow after work), I'd like to ask whether anyone here has encountered this previously and developed a simpler, yet reasonably efficient distance estimator for paraboloids.
I made an attempt to calculate an analytical distance. If we assume the paraboloid is in point POS, and its axis points to a direction DIR, then define D_ORTHO as the distance orthogonal to the paraboloid axis and D_AX as the distance along the paraboloid axis, we can walk the paraboloid surface along the orthogonal direction and form the actual distance function quite easily. I'm getting something like DIST = sqrt((X^2-D_AX)^2 + (-X + D_ORTHO)^2), where X is the distance we walk into the axial direction from the paraboloid axis.
Then, derivating and solving for the roots of X to find the minimum distance, I'm getting 4X^3 + (-4D_AX+2)X - 2D_ORTHO = 0, which is irritatingly a third degree polynomial. The solution seems to be correct, however, at least on cursory inspection and using a random cubic polynomial solver from Wolfram.
Example with solution: https://i.imgur.com/Xd23YjI.png
Random cubic solver from Wolfram: http://www.wolframalpha.com/widgets/view.jsp?id=3f4366aeb9c157cf9a30c90693eafc55
In the above image, I simply have an (x^2) paraboloid with a (2, -1) point that I'm solving for for demonstration purposes.
Now, a cubic solver is general knowledge, and one can of course be implemented, but it contains a ton of cubic roots so it isn't going to be particularly fast, and there are other irritating side effects like having to check which of the solutions will be real in the first place. Before I delve into this (which I will do anyway tomorrow after work), I'd like to ask whether anyone here has encountered this previously and developed a simpler, yet reasonably efficient distance estimator for paraboloids.
Actual image instead of link:
You should read up on the so called "lipschitz"-constant. As long as your object-calculation doesn´t deliver a lipschitz-correct distance you will overshoot as you said!
You can overcome this problem in some cheaty ways if you can´t come up with such a correct distance also:
1) always only return a fraction of the calculated distance...this way you may still overshoot, but not that visibly (in lower amounts), just let your DE for this object return it´s distance like "0.2*distance" or even lower maybe! This obviously slows your whole shader down a lot! ;)
2) make use of bounding-spheres/boxes...march the distance of these, and only if you are inside the bounding volume return the distance of your object...this way you won´t miss the object that easily, or also move inside of it by overshooting. This could even speed up your shader a bit! ;)
3) A combination of 2)+3), which should give you some graphics with very few artefacts, as soon as you have parametrized it in a nice way.
Maybe this helps, but you should of course try to find a lipshitz-conform DE first, if possible, while the above tricks are used quite commonly in demoscene-raymarching, as we just want to have nice graphics and dont need them to be correct in math, which would only slow our effects down in most cases! :D
You can overcome this problem in some cheaty ways if you can´t come up with such a correct distance also:
1) always only return a fraction of the calculated distance...this way you may still overshoot, but not that visibly (in lower amounts), just let your DE for this object return it´s distance like "0.2*distance" or even lower maybe! This obviously slows your whole shader down a lot! ;)
2) make use of bounding-spheres/boxes...march the distance of these, and only if you are inside the bounding volume return the distance of your object...this way you won´t miss the object that easily, or also move inside of it by overshooting. This could even speed up your shader a bit! ;)
3) A combination of 2)+3), which should give you some graphics with very few artefacts, as soon as you have parametrized it in a nice way.
Maybe this helps, but you should of course try to find a lipshitz-conform DE first, if possible, while the above tricks are used quite commonly in demoscene-raymarching, as we just want to have nice graphics and dont need them to be correct in math, which would only slow our effects down in most cases! :D
I'd cheat by calculating the distance to a cone (yay stackoverflow), then taking the square root of its axis (eg.
Code:
if the cone is along the x-axis, where p is the current point) before calculating the distance. It's an ugly hack, I haven't tried it, but it might work.p.x=sqrt(p.x);
Ok the code there got messed up a little.
And what hardy said, it doesn't have to be exact.
And what hardy said, it doesn't have to be exact.
just so everyone can follow:
DE = Distance Estimation
DE = Distance Estimation
An even worse hack: just use a very large ellipsoid (i.e. a scaled sphere)
I know a good anagram of paraboloid. ;)
Joking aside, if you only need the distance for marching, you can just as well use the true distance to the paraboloid along the ray, which has a much simpler formula (only second degree). In other words, just intersect your ray with the paraboloid directly.
If you also need the distance for other things (such as fake AO or some such), that doesn't help you much of course, though you could choose to use the complicated nearest distance formula just for that (if the structure/constaints of your code allow it).
Joking aside, if you only need the distance for marching, you can just as well use the true distance to the paraboloid along the ray, which has a much simpler formula (only second degree). In other words, just intersect your ray with the paraboloid directly.
If you also need the distance for other things (such as fake AO or some such), that doesn't help you much of course, though you could choose to use the complicated nearest distance formula just for that (if the structure/constaints of your code allow it).
@ɧคɾɗվ: Already doing 2) Also compartmentalizing the shape itself with multiple recursive bounding boxes to avoid checking all shapes at all positions (at the moment this object has ~30 primitives and I'm probably going up to 40-50 eventually). Also, already tried doing 1) but the slowdown was a bit unacceptable, considering how many objects I have.
@porocyon: A scaled sphere isn't the right shape, but the DE is really easy and it might not matter here. I'm at least checking this possibility out before doing the cubic solver.
@Blueberry: Firstly, I missed this possibility entirely, and didn't even consider that the ray-specific formula might be easier. However, there is a slight problem with this in the sense, that our (i.e. Trilkk's) framework does not pass the ray direction to the SDF functions, but only the current position of the ray. The aim is to then return the largest safe distance we can traverse before we hit anything and march that + overshoot, which produces a reasonably good estimate in most cases.
@porocyon: A scaled sphere isn't the right shape, but the DE is really easy and it might not matter here. I'm at least checking this possibility out before doing the cubic solver.
@Blueberry: Firstly, I missed this possibility entirely, and didn't even consider that the ray-specific formula might be easier. However, there is a slight problem with this in the sense, that our (i.e. Trilkk's) framework does not pass the ray direction to the SDF functions, but only the current position of the ray. The aim is to then return the largest safe distance we can traverse before we hit anything and march that + overshoot, which produces a reasonably good estimate in most cases.
...I can make it use the ray direction too if you want. I'd prefer not to unless it's absolutely necessary though, since that's some amount of rewriting.
Out of context, weren't Gamma2 blobs faked by tracing a paraboloid?
Just make the rayDirection a global variable and there´s about no amount of refactoring any code! ;)
I was too tired yesterday but today I investigated this further and implemented a part of the cubic solver. However, there is a problem with the implementation of the pow function in GLSL, as I have to take the 3rd root of a negative number. While this is possible in my specific corner case and produces a real result, the implementation of pow inside glsl is incapable of handling it due to not handling complex numbers by design (obviously).
In other words, unless we change the framework and include the ray direction, this can't be done. I think I'm falling back to an ellipsoid.
Thanks everyone.
The initial implementation was like this for anyone interested:
And yes, I am aware that I'm only investigating the always real root of the polynomial (positive discriminant) and the above function could thus fail inside the cup of the paraboloid, where there are additional roots, even if I had a pow function which could handle this. However, it's enough for you to test and see that it can't be done this way.
In other words, unless we change the framework and include the ray direction, this can't be done. I think I'm falling back to an ellipsoid.
Thanks everyone.
The initial implementation was like this for anyone interested:
Code:
float cubicsolver(float a, float b, float c, float d)
{
float x1 = -b/(3*a);
float x2 = -b*b*b/(27*a*a*a) + b*c/(6*a*a) - d/(2*a);
float x3 = sqrt( pow(x2, 2) + pow(c/(3*a) - b*b/(9*a*a), 3) );
return x1 + pow(x2 + x3, 1/3) + pow(x2 - x3, 1/3);
}
/// SDF of a paraboloid.
/// \param pos Position of ray.
/// \param p1 Origin of paraboloid.
/// \param p2 Axis vector of paraboloid.
/// \param radius Cutoff of paraboloid in axis-ortho direction.
/// \param steepness Paraboloid steepness (i.e. steepness*x^2).
/// \param wall_thickness Paraboloid wall thickness.
float sdf_paraboloid(vec3 pos, vec3 p1, vec3 p2, float radius, float steepness, float wall_thickness)
{
vec3 dir_pos = pos - p1;
vec3 proj_p2 = p2*(dot(dir_pos, p2)/dot(p2,p2));
vec3 proj_p2_ortho = dir_pos - proj_p2;
float len_proj_p2_ortho = length(proj_p2_ortho);
float len_proj_p2 = length(proj_p2);
float x = cubicsolver(steepness*steepness*4, 0, -steepness*4*len_proj_p2+2, -2*len_proj_p2_ortho);
float dist = sqrt(pow(x*x-len_proj_p2, 2) + pow(len_proj_p2_ortho-x, 2));
if(len_proj_p2_ortho > radius)
{
dist = max(dist, len_proj_p2_ortho - radius);
}
return dist;
}
And yes, I am aware that I'm only investigating the always real root of the polynomial (positive discriminant) and the above function could thus fail inside the cup of the paraboloid, where there are additional roots, even if I had a pow function which could handle this. However, it's enough for you to test and see that it can't be done this way.
Quote:
While this is possible in my specific corner case and produces a real result, the implementation of pow inside glsl is incapable of handling it due to not handling complex numbers by design (obviously).
In other words, unless we change the framework and include the ray direction, this can't be done.
It is, of course, impossible to compute complex cubic roots without native support for it in the language's standard library. :-)
Interestingly I implemented something like this not too long ago and got it to work with a cubic root solver. Nevertheless, the solution was not really numerically stable and had issues, especially when I wanted to compute the gradients/normals.
I ended up using a really basic, hacked together approach (just divide by the implicit function by the length of its gradient):
If you already have something similar... You might want to improve your sphere tracing code - the above function just worked surprisingly well for my use case. It's maybe not optimal in terms of the returned distances, but most probably much faster to evaluate than any cubic root solver based variant.
I ended up using a really basic, hacked together approach (just divide by the implicit function by the length of its gradient):
Code:
float fParaboloid(float3 p) {
return max(p.z-1, (p.x*p.x + p.y*p.y - p.z) / sqrt(4*p.x*p.x + 4*p.y*p.y + 1));
}
If you already have something similar... You might want to improve your sphere tracing code - the above function just worked surprisingly well for my use case. It's maybe not optimal in terms of the returned distances, but most probably much faster to evaluate than any cubic root solver based variant.
This is close enough for the sphere tracing algorithms I use:
damn, this bigger screenshot with its debug-colors is mesmerizing!
Look at its middle for fifteen seconds, then watch away from your screen!! ;)
(outcome may vary depending on applied resolution and used page-zoom!)
Look at its middle for fifteen seconds, then watch away from your screen!! ;)
(outcome may vary depending on applied resolution and used page-zoom!)
distanceToParabolaSegment is identical to distanceToQuadtaticBezier, without the boundaries of the bezier.
shadertoy.com has functions, that let you input a quadratic Bezier, and it converts that to a parabola.
quadraticBedier to point shortestDistance is overdone on shadertoy.
there are also many fast (non itterative) rough estimate/bound solutions.
The PRECISE complexity of this is low, it computes the 2 local minima of a quadratic function (the quadratic here always has 1 or 2 minima and 0 or 1 maxima, and it just is never the local maxima)
the 2 end points are some dot() bound? i forgot the details, its VERY common, and just copy pasted a lot.
bonus points, there are even a lot of parametrizations of quadraticBeziers on shadertoy (== with uniform distance of points along barabola/bezier)
shadertoy.com has functions, that let you input a quadratic Bezier, and it converts that to a parabola.
quadraticBedier to point shortestDistance is overdone on shadertoy.
there are also many fast (non itterative) rough estimate/bound solutions.
The PRECISE complexity of this is low, it computes the 2 local minima of a quadratic function (the quadratic here always has 1 or 2 minima and 0 or 1 maxima, and it just is never the local maxima)
the 2 end points are some dot() bound? i forgot the details, its VERY common, and just copy pasted a lot.
bonus points, there are even a lot of parametrizations of quadraticBeziers on shadertoy (== with uniform distance of points along barabola/bezier)
This seems to work well: https://www.shadertoy.com/view/3lSczz. Just geed it with vec2( length(p.xz), p.y ) to make it revolve into a paraboloid.
*feed