Welcome to Fractal Forums

i am wondering how to calculate a DistanceEstimator value for that Amazin Fractal
http://www.fractalforums.com/3d-fractal-generation/amazing-fractal/

i am a bit stuck

when thinking about it, a distance estimator is not an option here ?
i have to binary search the ray for intersection ?
please correct me if i am wrong  O0

Well, DE for the Mset is based on the Koebe 1/4 Theorem (yey 4 wikipedia!), which does not necessarily hold in the complex plane. That doesn't mean DE wouldn't work in vector space. If we have the mapping Hⁿ_c(x), we could attempt distance estimation with ln(abs(d(Hⁿ_c(x))/dc)) * abs(Hⁿ_c(x)) / abs(d(Hⁿ_c(x))/dc). The derivative of the function H with respect to c... isn't easy. d(H¹_c(x))/dc is one, that's simple enough, but then we use the chain rule to extend the derivative with each iteration... if we have DE_n, z_n, z_n+1 as our variables, we could calculate DE_n+1 using some pretty weird conditional hooey.

And you'd need to define multiplication of these 3D numbers.

Hehe, I wrote this post in the hope I could think of something, and I've convinced myself in the meantime that it would utterly fail  ;D

It would have to be a conditional derivative, as the function is non-continuous... but even then for:
point.x > 1
point.x = 2-point.x    which is basically   f(x) = 2 - x  which makes f ' (x) = -1  which may or may not be useful for DE. 

Hmm... yeah, and for the spherical inversion you'd have a derivative of 0 for radii of less than 1/2, a derivative of 1 for radii greater than 1, and -1/x^2 for radii in between. Scaling just gives a derivative of 2. So I guess if we say H(x)=M(I(C(x))) + pixel (M means magnify, I is invert, C is cube), the derivative of H(x) = C'(x)*I'(C(x))*M'(I(C(x))), which could be calculated with the following code:

loop:
 derivOld=deriv
 deriv=point.x
 if point.x > 1 || point.x < 1 then
  deriv=-1
 else
  deriv=1
 endif
%Code for folding goes here...
%Code for calculating radius goes here...
 if length<0.5
  deriv=0
 elseif length<1
  deriv=-deriv/(point.x*point.x)
 endif
%Code for spherical inversion goes here...
 deriv=deriv*2
%Code for scaling goes here...
%Code for constant addition goes here...
 deriv=deriv*derivOld+1

I think this would work. Because DE isn't really legit for 3D space, you might do something like computing separate DE's using .x, .y, and .z values and averaging them or something...  :hmh:

With respect to the analytical DE, it's probably better for this type of fractal to consider the analytical DE for IFS rather than standard escape-time fractals.
For affine IFS (assuming uniform scaling and no skewing) then the DE is simply:
 cabs(final magnitude)/total scale
Which actually makes sense since here "total scale" is the derivative, though quite why multiplying by log(cabs(final magnitude)) is definitely not correct in this case I'm not sure :)
However for IFS with non-uniform scaling and/or skewing etc. then you have to get the magnitude of the determinant of the product of the transformation matrices used and of course use the correct matrix from each iteration - in the case of rendering escape-time IFS where the entire tree is traversed per point then you have to store the path of the closest route and compute the product of the matrices for that route after you've finished traversal - but in the case of the Mandelbox/Juliabox you only have one route.

One of the great things about UltraFractal 5 is that the architecture of the programming language allows plugins. My 3D Fractal Raytrace (UF5) was designed to handle virtually any 3D or 4D non-IFS fractal. The first plugins were for quaternion and FractInt style hypercomplex fractals, but it works with every other type I have tried. I have already prsesented some Mandelbulb examples, and now here is an Amazing Fractal example. 3D Fractal Raytrace (UF5) needs the derivative for the potential function for approaching the surface. For the Amazing Fractal I used the Delta DE potential method. DE in this context stands for DErivative, not Distance Estimate. In the Delta DE method the derivative is estimated at each iteration by adding a small delta to each of the x, y and z components and then subtracting the resultant value to get dx, dy and dz. This method works very well for the Amazing Fractal and for all other types I have tried. My other potential methods, which are often faster, require an explicit derivative calculation. I found that Timeroot's approach works reasonably well for my other potential methods.

Here is Cathedral, which uses scale = 3

Unfortunately the Ultrafractal database is down right now, and I don't know when it will be back up. I will upload the plugin for Amazing Fractal when the database is back up.

ok, i have to ask again, about the Delta DE Method


my understanding of the algorithm for Delta DE:

z=initialzvalue
while length(z)<bailout and itercount<maxIter
   z=fractalfunc(z)
   zderive = fractalfunc(z+delta)-z ?!??!?!???????????
end

You missed a rather important factor (1/delta) and of course you don't add "delta" itself to z, but a vector of length delta in the appropriate direction (or directions if using Buddhi's method) i.e.  zderive = (fractalfunc(z+delta*vector)-z)/delta

Edit: correction that should be:  zderive = (fractalfunc(initialz+delta*vector)-z)/delta

is that vector the delta from the last z value and the current one ?
or how is that vector achieved ?!?

The vector (or vectors) is a unit vector in an appropriate direction, exactly what depending on whether you impliment my version or Buddhi's version of the method.

In my version the vector is the unit direction vector of the viewing ray for the current pixel i.e. you iterate at position "step" on the ray and at position "step+delta". This method suffers from the same problem as my older delta DE method that uses the difference in the smooth iteration values rather than the final z value in that rays travelling parallel to the surface will produce very small zderiv values and so you have to program around this (for example using the maxdist array method that I've outlined in a number of places).

In Buddhi's version you actually use 3 unit vectors (1,0,0), (0,1,0) and (0,0,1) i.e. your four locations to iterate are "step" position (px,py,pz) plus (px+delta,py,pz), (px,py+delta,pz) and (px,py,pz+delta). From these you calculate 3 values of zderiv and combine them using zderivfinal = sqrt(zderiv1^2+zderiv2^2+zderiv3^2) (at least that's what I believe Buddhi is doing). This method avoids the directionality problems associated with using the viewing ray vector. I haven't actually implimented this myself yet as obviously this requires 4 iterations per ray step rather than just 2. I should add that an alternative to using the pythagoras is (possibly) to simply use the largest of the 3 (absolute) values of zderiv or to use the quick and dirty methods (using max,med and min) for 3D pythagoras.

Hm...can't say I've encountered this method before. I understand what you're doing to get zderive, but what is it used for? In Trifox's pseudocode, it's not actually applied.

It's simply applied to get a distance estimate in place of the analytical derivative i.e.

DE = 0.5*log(cabs(z))*cabs(z)/zderiv

Note that for Buddhi's method you can avoid the triple divide by delta by calculating zderiv without using the divide and then use:

DE = 0.5*delta*log(cabs(z))*cabs(z)/zderiv

With respect to the analytical derivative for the Mandelbox can someone more maths savvy than myself please answer the question as to what is the derivative of f(x,y,z) where f(x,y,z)=(-x,y,z) or f(x,y,z)=(x,-y,z) or f(x,y,z)=(x,y,-z) as we need the derivative for the rectangular folding ?
Can we just use the obvious e.g. if the formula is say folding the x as (v-x) then we simply negate dx ?

I should point out that that's how the delta z method is applied, my older delta DE method using the smooth iteration values is slightly different.

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."

For a noticeable speed-up, change this


to this:


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".

Oh, wow, thanks folks. I understand what you're doing now. I'll have to give that a shot in my own code. :D

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.  O0

one more question:

does the DeltaDE gets negative values when inside ?
for isosurfacing of the formula

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;

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.

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.

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 :D
ok, i will try that ... thank you very much for helping again! :angel1:

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.

I think the answer is in your question - "depending on the method i think ( julia/mandelbrot )" :)

how does your binary search works then ?!
i used to have -1 for inside, and 1 for outside :D

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.

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 ;)

It just struck me that a more GPU friendly manner to accomplish the same thing as the binary search is to:

1. Step until DE<closeness or "inside" - set "count" to zero, set teststep to laststepdistance/n
2. Step backwards by teststep and increment "count" and get DE for new position
3. if DE<closeness or "inside" and count<n goto 2
4. point found

Use larger values of n for more accuracy, I would guess around 10 would be fine in most cases.
Technically you shouldn't need the "count" at all but....