cKleinhuis
|
|
« Reply #15 on: March 25, 2010, 12:58:26 AM » |
|
i got this pseudocode from fractalrebel via private mail, and i think it is worth to make public, to have a decent description about what is going on
fractalrebel: "It went faster than I thought. I hope this will be satisfactory:
This is a metacode approach to the main calculation sequence in REB_3DFractalRaytrace. There are four nested loops in the main calculation. Before the loops begin calculations are carried out to determine the minimum distance (mindist) and the maximum distance (maxdist). Mindist is a distance between the camera and the object which acts as the distance starting point. Maxdist lies on the other side of the object and determines when distance calculations stop if the ray misses the object. The mindist and maxdist values come from two planes whose position is defined by the user.
The loops fill two complex arrays, one for the x-y final values for the fractal surface and the other array for the z values for the fractal surface. The imaginary portion of the second array holds information that aids in coloring the fractal.
REPEAT ; loops over the y screen/pixel coordinates
REPEAT ; loops over the x screen/pixel coordinates In this section initial x, y and z values are calculated based upon the target screen pixel, the camera position, the mindist value and the camera vector.
REPEAT ; iterates until the surface or maxdist is reached In this section x, y and z values updated to reflect the new distance calculated in the inner most loop. The innermost loop distance calculation (called d) is added to the current distance value (called dist). This loop continues until d < epsilon or dist > maxdist. REPEAT ; iterates x, y and z until maxiter reached or bailout In this section the fractal is iterated, as is the derivative of the fractal. Only Delta DE will be discussed in detail here. Consider the vector R = (x,y,z,w). We can effectively ignore w for most purposes, so consider the vector to be R = (x,y,z). With each iteration the length of R is calculated, with bailout occuring when length(R) = L > bailout. Since we are ignoring the 4th dimension and differences between Julia type and Mandlel type fractals, we only need to calculate three additional lengths:
Lx = length(R(x+delta,y,z)) Ly = length(R(x,y+delta,z)) Lz = length(R(x,y,z+delta))
From L, Lx, Ly and Lz we can estimate the derivative at the point (x,y,z):
dx = (L-Lx)/delta dy = (L-Ly)/delta dz = (L-Lz)/delta
This is a vector, and we need the magnitude of the vector:
dr = length(dx,dy,dz)
We now have everything we need to estimate the distance from the current postion to the surface. Each time through the loop the estimate will be smaller, as we will be closer to the surface. The stepsize parameter is to insure that any given step will not be too large and go through the surface. Typically a stepsize of 1.0 works very well.
d = L*log(L)/dr*(stepsize/2)
d is added to dist to get the current distance estimate from the camera to the surface. When d < epsilon the loop ends. Epsilon typically has dimensions on the order of pixel size. UNTIL ; iter = maxiter or sqrt(x^2+y^2+z^2) > bailout UNTIL ; distance < epsilon UNTIL ; all x pixel values UNTIL ; all x pixel values
Downstream there are calculations relating to other objects such a planes, shadows on the planes and shadows on the object, raytracing vs distance color calculations, etc.
I recently plugged in Mstloe's new algorithm involving applications of symmetry. It took about 10 minutes worth of coding, and it works just fine. I am now working on solving the problems with y -z values which he noted.
The plugins are all 4D vector format, and I use a couple of helper functions to convert between a 4D vector and two complex numbers."
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
David Makin
|
|
« Reply #16 on: March 25, 2010, 01:04:54 AM » |
|
For a noticeable speed-up, change this Lx = length(R(x+delta,y,z)) Ly = length(R(x,y+delta,z)) Lz = length(R(x,y,z+delta))
From L, Lx, Ly and Lz we can estimate the derivative at the point (x,y,z):
dx = (L-Lx)/delta dy = (L-Ly)/delta dz = (L-Lz)/delta
This is a vector, and we need the magnitude of the vector:
dr = length(dx,dy,dz)
We now have everything we need to estimate the distance from the current postion to the surface. Each time through the loop the estimate will be smaller, as we will be closer to the surface. The stepsize parameter is to insure that any given step will not be too large and go through the surface. Typically a stepsize of 1.0 works very well.
d = L*log(L)/dr*(stepsize/2)
to this: Lx = length(R(x+delta,y,z)) Ly = length(R(x,y+delta,z)) Lz = length(R(x,y,z+delta))
From L, Lx, Ly and Lz we can estimate the derivative at the point (x,y,z):
dx = (L-Lx) dy = (L-Ly) dz = (L-Lz)
This is a vector, and we need the magnitude of the vector:
dr = length(dx,dy,dz)
We now have everything we need to estimate the distance from the current postion to the surface. Each time through the loop the estimate will be smaller, as we will be closer to the surface. The stepsize parameter is to insure that any given step will not be too large and go through the surface. Typically a stepsize of 1.0 works very well.
d = delta*L*log(L)/dr*(stepsize/2)
Also I find it better to separate the distance estimate from the distance to step so that one can modify the closeness separately from the accuracy used, this can be performed simply by changing the last bit to: de = 0.5*delta*L*log(L)/dr Where de is the distance to test against the closeness threshold, and then (if not close enough) use: step = de/@accuracy Where an accuracy of 1.0 will be OK but larger values can be used for problematical renders - increasing the accuracy at a speed cost *but* independant of the closeness to "inside".
|
|
« Last Edit: March 25, 2010, 01:11:51 AM by David Makin »
|
Logged
|
|
|
|
reesej2
Guest
|
|
« Reply #17 on: March 25, 2010, 01:49:04 AM » |
|
Oh, wow, thanks folks. I understand what you're doing now. I'll have to give that a shot in my own code.
|
|
|
Logged
|
|
|
|
Timeroot
|
|
« Reply #18 on: March 25, 2010, 01:52:59 AM » |
|
Ahh, I was just about to post about this when David posted the question about derivative of (-x,y,z). I learned something interesting last night while watching Alfred Hitchcock. Since a holomorphic function must satisfy the Riemann-Cauchy equations, and anything that satisfies them is conformal, and Liouville's theorem limits the 3D conformals, the only things you can take the derivative of in 3D space - any function, any coordinate system - are scaling, translation, inversion (that is, reciprocal of the radius, along with any one reflection), and reflection. And reflection doesn't work either, as can easily be shown that the complex conjugate function doesn't satisfy the Riemann-Cauchy equations either. So, no, you can't take the derivative of (-x,y,z), or anything else, really.... this might explain why the method for DE I suggested earlier worked! With the exception of the reflection, which used a bit of fudging with a conditional derivative. One sad result of this is that it will probably never be possible to construct a triplex equivalent of the Newton fractal. Of course, you could fudge to find your "derivative" in other ways; viewing it as a vector-valued function, you get the derivative of (-1,1,1). How you will use that, though... I don't know. Another possibility is to take a directional derivative. But then you need to decide on a direction. One possibility would be taking the average derivative over all possible directions which would be... 2/3. Really useful there, yeah. In my method for computing the DE, I calculated each variable's derivative separately, it worked fine.
|
|
|
Logged
|
Someday, man will understand primary theory; how every aspect of our universe has come about. Then we will describe all of physics, build a complete understanding of genetic engineering, catalog all planets, and find intelligent life. And then we'll just puzzle over fractals for eternity.
|
|
|
cKleinhuis
|
|
« Reply #19 on: March 25, 2010, 05:04:55 PM » |
|
one more question:
does the DeltaDE gets negative values when inside ? for isosurfacing of the formula
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
cKleinhuis
|
|
« Reply #20 on: March 25, 2010, 08:59:00 PM » |
|
the Delta De function as implemented by me on gpu does always return 0 distance fracfunc is actually the mandelbox .... input is the 3d coordinate we want to measure the distance for ... seed is the global defined julia/pertubation value for that fractal float delta=0.10; float stepsize=1.0; float iterateDeltaDE(vec3 pixel){
vec3 z=seed ; vec3 c=pixel; int i=0; vec3 L; vec3 d1; float Llength; float d; float dr;
while((i++<=maxIterationDepth)&& ( length(z)<16.0) ) { z=fracFunc(z,c); Llength=length(z); L.x = length(vec3(z.x+delta,z.y,z.z)); L.y = length(vec3(z.x,z.y+delta,z.z)); L.z = length(vec3(z.x,z.y,z.z+delta)); d1.x = (Llength-L.x)/delta; d1.y = (Llength-L.y)/delta; d1.z = (Llength-L.z)/delta; dr= length(d1); d = Llength*log(Llength)/dr*(stepsize/2.0); } return d;
}
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
David Makin
|
|
« Reply #21 on: March 25, 2010, 09:45:41 PM » |
|
the Delta De function as implemented by me on gpu does always return 0 distance You're misunderstanding the code I think - you have to do the full iteration for *4* points. Edit: corrected to use the lengths instead of final z values ! Get z and c for the base point (x,y,z) and iterate to get r=length(final z). Now in turn get z and c for the points (x+delta,y,z), (x,y+delta,z) and (x,y,z+delta) and iterate each giving final z lengths rx,ry and rz. Then your derivative estimate is length(r-rx,r-ry,r-rz)/delta.
|
|
« Last Edit: March 25, 2010, 11:00:17 PM by David Makin »
|
Logged
|
|
|
|
David Makin
|
|
« Reply #22 on: March 25, 2010, 09:48:24 PM » |
|
one more question:
does the DeltaDE gets negative values when inside ? for isosurfacing of the formula
No, if you detect "inside" you treat it the same way as if the DE value is less than the closeness threshold i.e. you either just take that as the solid point or you instigate a binary search between that point and the previous step position to get the exact solid boundary.
|
|
|
Logged
|
|
|
|
cKleinhuis
|
|
« Reply #23 on: March 25, 2010, 09:51:48 PM » |
|
omg, thank you for clarifying ... what do you mean by get c for each point ? c is the same for all ? depending on the method i think ( julia/mandelbrot ) ok, 1. Iterate Point for ( x,y,z ), and save resulting z value, rebel called it L 2. Iterate Point for ( x+delta,y,z ), and save resulting zx value, rebel called it Lx 3. Iterate Point for ( x,y+delta,z ), and save resulting zy value, rebel called it Ly 4. Iterate Point for ( x,y,z+delta ), and save resulting zz value, rebel called it Lz omg, i hoped it would be something faster ok, i will try that ... thank you very much for helping again!
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
cKleinhuis
|
|
« Reply #24 on: March 25, 2010, 09:53:45 PM » |
|
one more question:
does the DeltaDE gets negative values when inside ? for isosurfacing of the formula
No, if you detect "inside" you treat it the same way as if the DE value is less than the closeness threshold i.e. you either just take that as the solid point or you instigate a binary search between that point and the previous step position to get the exact solid boundary. ok, thank you so, what is your opinion about binary searching when this happens ? because normally when marching the ray, and decreasing the step size in respect of the distance this should never happen ?!
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
David Makin
|
|
« Reply #25 on: March 25, 2010, 10:15:20 PM » |
|
one more question:
does the DeltaDE gets negative values when inside ? for isosurfacing of the formula
No, if you detect "inside" you treat it the same way as if the DE value is less than the closeness threshold i.e. you either just take that as the solid point or you instigate a binary search between that point and the previous step position to get the exact solid boundary. ok, thank you so, what is your opinion about binary searching when this happens ? because normally when marching the ray, and decreasing the step size in respect of the distance this should never happen ?! The binary search is only really needed when the "closeness" is set quite large and/or the "accuracy" setting in my implimentation is set very low but it allows the accuracy setting to be set lower than otherwise and thus can be used for optimising render times.
|
|
|
Logged
|
|
|
|
David Makin
|
|
« Reply #26 on: March 25, 2010, 10:19:32 PM » |
|
what do you mean by get c for each point ? c is the same for all ? depending on the method i think ( julia/mandelbrot )
I think the answer is in your question - "depending on the method i think ( julia/mandelbrot )"
|
|
|
Logged
|
|
|
|
cKleinhuis
|
|
« Reply #27 on: March 25, 2010, 10:32:43 PM » |
|
how does your binary search works then ?! i used to have -1 for inside, and 1 for outside
|
|
|
Logged
|
---
divide and conquer - iterate and rule - chaos is No random!
|
|
|
David Makin
|
|
« Reply #28 on: March 25, 2010, 10:33:37 PM » |
|
I'll just add that using a "closeness" value that is accurate enough for rendering the shape to a given resolution will often be too large to be accurate enough for correct lighting normals - this is another reason for using a binary search since having to reduce the closeness value further (or accuracy) just to get the correct normals is not optimum.
|
|
|
Logged
|
|
|
|
David Makin
|
|
« Reply #29 on: March 25, 2010, 10:41:18 PM » |
|
how does your binary search works then ?! i used to have -1 for inside, and 1 for outside Step along the ray using the step values you calculate from the DE, if at a new step position either the DE is less than the closeness threshold *or* the point is inside then set "binary step" to "last step distance/2", step backwards by that distance and continue through your rendering loop as normal but with a "binary search" flag enabled and in your code if this is enabled after calculating the new DE value divide "binary step" by 2 then decide whether to go forwards (if DE>closeness) or backwards (DE<closeness or "inside") and step by "binary step" in the appropriate direction. Obviously there are a number of different ways of terminating the binary search, e.g. have a binary search count and stop when a certain number of binary search steps has been performed or stop when the "binary step" value is less than a given threshold or stop when the new DE value is very close to the closeness value. I should add that implimenting it in a GPU friendly manner is not easy
|
|
« Last Edit: March 25, 2010, 10:44:55 PM by David Makin »
|
Logged
|
|
|
|
|