Aternative 3D Mandelbrot ?

[David Makin]

Apologies if anyone's tried this before, I would do so myself but I'm sick of the noise of my Windows PC and don't have any development code for 3D fractals on my Mac yet.
Anyway for a 3D point (x,y,z) with x horizontal left to right, y vertical bottom to top and z horizontal near to far:
To square (x,y,z) rotate (x,y,z) around the x axis by theta=-atan2(y/z)? such that the point ends up in the x/z plane (with z positive) (i.e. newx,0,newz) then rotate by the angle phi=atan2(newz/newx)? then rotate by -theta around the x axis - and obviously scale by the magnitude
Haven't tried it but I think it is either an alternative algorithm for quaternions or it may be altogether more interesting

Edit: Apologies for the typo - I should have said the angle phi around the y axis, though I think the atan is correct !!!

[cKleinhuis]

you have a typo: what is your second axis around, surely not x-axis a second time

i think it is covered in a generalized approach, i would like to define the mandelbulb formula by simple 2 base vectors, which serve as additional 3 parameters, they do not need to be orthogonal ( although it seems to be the most used approach )

with the 2 base vectors defining the coordinate system you can experiment with every direction or axis combination you like, this would get rid of always define new formulas for just changing one axis ( from forward to backward, or the order )

spherical coordinate systems can be easily extended to work with 2 arbitrary functions, just the theta and phi value calculation would have to be adjusted, but this would be a straight forward task, i am right now having access to maple, and can support you with formulas for arbitrary axis calculation

standard mandelbulb formula ( nylanders method) would have base vectors alpha/beta of (0,0,1) <- the complex plane, rotation around z-axis, and (0,-1,0) <- the y axis reversed ( because he is rotating clockwise )

and your above formulation would be for axes alpha and beta (1,0,0) and (0,1,0) for y axis beta, or (0,0,1) for z axis beta

clearly the base vectors have to be normalized!

[cKleinhuis]

oh, i see you have a different thought, by taking the transformed values for rotating, this is new, please try it out

[David Makin]

you have a typo: what is your second axis around, surely not x-axis a second time

Corrected in the edit - I of course meant phi as the angle around the y axis from the positive x axis i.e. the angle of (xnew,0,znew) to the positive x axis

[David Makin]

Just had a rethink about that method - an easier way to impliment an equivalent method is as follows:

With (x,y,z) such that x is on the horizontal left to right, y is the vertical bottom to top and z is the horizontal near to far:

First rotate around the x axis by angle theta such that we get a new point (x,ynew,0) where ynew is positive then simply treat (x,ynew) as a complex number and perform any function on it giving a new complex value (xnew,ynew2) then rotate (xnew,ynew2,0)  by -theta around the x axis to give the new 3D value (xnew,ynew3,znew).

So in actual fact that's:

(x,y,z)
1. ynew = magnitude of (y,z)
2. (xnew,ynew2) = complex fn(x,ynew)
3. (ynew3,znew) = (y,z)*ynew2/ynew
New point is (xnew,ynew3,znew)

OK, so for (x,y,z)^2 then that's

(x,y,z) -> (x,sqrt(y*y+z*z),0) -> (x*x-y*y-z*z, 2*x*sqrt(y*y+z*z),0) -> (x*x-y*y-z*z, 2*x*y, 2*x*z)

So it is just another way of doing 3D quaternions *however* the interesting part is that *any* function (not involving translation) of this 3D quaternion form can be accomplished by performing the 3 stages above such that the 3D function is applied as the equivalent complex function on stage 2.

Unfortunately I already knew this - in fact exactly the same thing applies to full 4D quaternions
i.e.
Quaternionic fn(x,y,z,w)
(x,ynew) = (x,mag(y,z,w))
(xnew,ynew2) = complex fn(x,ynew)
(ynew3, znew, wnew) = (y,z,w)*ynew2/ynew
=(xnew, ynew3, znew,wnew)
 
Of course that doesn't cover z^p where both z *and* p are fully quaternionic but is fine when p is complex, float or int.

[David Makin]

I was thinking about this some more and went through the foillowing reasoning:

The method for quaternions gives the known issue with multiplication in that multiplication is not commutative - the issue here being that in the form I described then to perform a multiplication there are *two* separate theta angles which are used to put the two 3D values into two complex values which are then operated on but we now don't know which of the two theta angles to use after the complex operation.

However if the final operation (rotating back by the original theta) is changed for multiplication to rotating back by the sum of both thetas then we have a new number form that is in fact commutative, not only that but we can also modify the method such that this new final rotation is consistent with other functions such as z^p as follows:

Instead of:

(x,y,z)
1. ynew = magnitude of (y,z)
2. (xnew,ynew2) = complex fn(x,ynew)
3. (ynew3,znew) = (y,z)*ynew2/ynew
New point is (xnew,ynew3,znew)

Use:

(x,y,z)
1. ynew = magnitude of (y,z)
2. (xnew,ynew2) = complex fn(x,ynew)
3. (ynew3,znew) = complex fn(y,z)
4. (ynew4,znew2) = (ynew3,znew)*ynew2/mag(ynew3,znew)
New point is (xnew,ynew4,znew2)

This not only makes the multiplications commutative but also gives us a way of applying any (non-translational) function to this number form.
Without checking carefully I think this is simply an alternative to the Mandelbulb that will give similar results, but I believe the above method gives us a (reasonably optimum) way of obtaining higher functions such as transcendentals etc.

Multiply: (x,y,z)*(a,b,c)

1. ynew=mag(y,z), bnew=mag(b,c)
2. (xnew,ynew2) = (x,ynew)*(a,bnew)
3. (ynew3,znew) = (y,z)*(b,c)
4. (ynew4,znew2) = (ynew3,znew)*ynew2/mag(ynew3,znew)
New point is (xnew,ynew4,znew2)

Edit: just to clarify that of course in stage 2 and 3 of the multiplication the calculations are complex ones (so assume (x+ynew*i)*(a+bnew*i) on step 2 for example).

Ediit: Note that the above in this post has been modified so that the rescale in step 4 has *ynew2 instead of *ynew as I wrote originally !!

[David Makin]

OK - I confess I'm rediscovering the wheel as usual - see:

http://www.fractalforums.com/theory/triplex-algebra/msg21795/#msg21795

However I decided to play wiith log and exp and found that for this form they're straightforward (using C's atan2):

log(x,y,z) = (0.5*log(x*x+y*y+z*z), atan2(sqrt(x*x+y*y),x), atan2(z,y))

exp(x,y,z) = (exp(x)*cos(y), exp(x)*sin(y)*cos(z), exp(x)*sin(y)*sin(z))

Here's some "C" code:

Code:

void tExp(double &x, double &y, double &z)
{
	double m;
	x = exp(x);
	if (z==0.0f)
	{
		if (y!=0.0f)
		{
			m = x;
			x = x*cos(y);
			y = m*sin(y);
		}
	}
	else if (y==0.0)
	{
		y = z = 0.0f;
	}
	else
	{
		m = x*sin(y);
		x = x*cos(y);
		y = m*cos(z);
		z = m*sin(z);
	}
}

void tLog(double &x, double &y, double &z)
{
	double m,m1;
	if (z==0.0f)
	{
		if (y==0.0f)
		{
			if (x<0.0f)
			{
				y = -3.1415926f;
			}
			x = log(fabs(x));
		}
		else
		{
			m = x*x + y*y;
			y = atan2(y,x);
			x = 0.5*log(m);
		}
	}
	else if (y==0.0f)
	{
		m = x*x + z*z;
		y = atan2(fabs(z),x);
		if (z>=0.0f)
		{
			z = 0.5f*3.1415926f;
		}
		else
		{
			z = -0.5f*3.1415926f;
		}
		x = 0.5f*log(m);
	}
	else
	{
		m1 = sqrt(m = y*y + z*z);
		m += x*x;
		z = atan2(z,y);
		y = atan2(m1,x);
		x = 0.5f*log(m);
	}
}

Test calls:

	double x,y,z,x1,y1,z1;
	int count=0;
	srandom(round([[NSProcessInfo processInfo] systemUptime]));
	do
	{
		x = (float)random()/(65536.0f*256.0f)-64.0f;
		y = (float)random()/(65536.0f*256.0f)-64.0f;
		z = (float)random()/(65536.0f*256.0f)-64.0f;
		x1 = x;
		y1 = y;
		z1 = z;
		tLog(x, y, z);
		tExp(x, y, z);
		if ((fabs(x-x1)/(fabs(x)+1e-20)>0.001f)
			||(fabs(y-y1)/(fabs(y)+1e-20)>0.001f)
			||(fabs(z-z1)/(fabs(z)+1e-20)>0.001f))
		{
			x1 = fabs(x-x1)/(fabs(x)+1e-20);
			y1 = fabs(y-y1)/(fabs(y)+1e-20);
			z1 = fabs(z-z1)/(fabs(z)+1e-20);
		}
		x1 = x;
		y1 = y;
		z1 = z;
		tLog(x, y, z);
		tExp(x, y, z);
		if ((fabs(x-x1)/(fabs(x)+1e-20)>0.001f)
			||(fabs(y-y1)/(fabs(y)+1e-20)>0.001f)
			||(fabs(z-z1)/(fabs(z)+1e-20)>0.001f))
		{
			x1 = fabs(x-x1)/(fabs(x)+1e-20);
			y1 = fabs(y-y1)/(fabs(y)+1e-20);
			z1 = fabs(z-z1)/(fabs(z)+1e-20);
		}
	} while (++count<5000);