News: Visit the official fractalforums.com Youtube Channel  ## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages:    Go Down       Author Topic: 3D Mandelbrot, pole vector method (see Theory forum)  (Read 3096 times) Description: 0 Members and 1 Guest are viewing this topic.
cbuchner1
Fractal Phenom      Posts: 443 « on: December 10, 2009, 12:53:16 AM »

Rendered in voxel technique with some transparency. The pole vector is pointing (0,1,0), i.e. in y direction. For a brief description of my technique check out the Triplex Algebra thread in the Theory forum and/or look at the source code in the third post.

From the first angle this looks almost like your usual mandelbrot. But watch what happens when I rotate.     Oh wow, this was unexpected. I guess this qualifies as a bulb, huh? Now I really need an analytical method for this (and the derivative) to do some realtime raytracing.
Voxels are nice, but it's hard to make out surface details.

Christian
 « Last Edit: December 10, 2009, 02:41:11 AM by cbuchner1 » Logged
cbuchner1
Fractal Phenom      Posts: 443 « Reply #1 on: December 10, 2009, 01:43:01 AM »

Pole vector towards [1,1,1] and normalized. Bizarre creature.  « Last Edit: December 10, 2009, 01:44:47 AM by cbuchner1 » Logged
cbuchner1
Fractal Phenom      Posts: 443 « Reply #2 on: December 10, 2009, 02:19:54 AM »

partial C++ source code included. To use the code you should include outer loops to cover all voxels in x,y,z coordinates (double variables) and adjust the method by which the volume texture is stored to suit your needs

The best rendition of a mandelbrot is obtained by using a pole vector of {1.0, 0.0, 0.0}

Or maybe you're just interested in the maths.

Code:
// euclidean vector utility class with the required operator overloads
struct triplex
{
double x,y,z;

// magnitude (length) of this object
double len() { return sqrt(x*x+y*y+z*z); }

// return result of rotation around x axis
triplex rotx(double arg)
{
triplex result = { x,
cos(arg)*y - sin(arg)*z,
sin(arg)*y + cos(arg)*z };
return result;
}

// return result of rotation around y axis
triplex roty(double arg)
{
triplex result = { cos(arg)*x - sin(arg)*z,
y,
sin(arg)*x + cos(arg)*z };
return result;
}

// return result of rotation around z axis
triplex rotz(double arg)
{
triplex result = { cos(arg)*x - sin(arg)*y,
sin(arg)*x + cos(arg)*y,
z                                };
return result;
}

// vector addition between this object and another vector b
triplex operator + (const triplex &b)
{
triplex result = {x + b.x, y + b.y, z + b.z};
return result;
}

// vector division by a scalar b
triplex operator / (const double &b)
{
triplex result = { x / b, y / b, z / b };
return result;
}

// vector multiplication with a scalar b
triplex operator * (const double &b)
{
triplex result = { x * b, y * b, z * b };
return result;
}

// dot product between this object and another vector b
double operator * (const triplex &b)
{
return x * b.x + y * b.y + z * b.z;
}

// cross product between this object and another vector b
triplex operator % (const triplex &b)
{
triplex result =  { y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x };
return result;
}

// raise this object to the power b relative to the pole vector
triplex operator ^ (const double &b)
{
triplex Z_test = *this;

// storing the magnitude saves us multiple calls to len()
double magnitude = len();

// compute normal vector that defines the plane in which the result shall live
triplex N = (pole() % *this) / magnitude;

// compute angle theta between current vector and the pole vector
// which is of range 0 to PI
double theta = acos( (*this * pole()) / magnitude );

// generate a result vector that is situated in the plane
// spanned by origin, pole and this object and has an angle theta
// to the pole and is of the required magnitude.

// rotate Z space such that normal vector points "up" (in z direction)
double N_phi   = atan2(N.y, N.x);
double N_theta = atan2(-N.x, N.z);
Z_test = Z_test.rotz(-N_phi);
Z_test = Z_test.roty(-N_theta);

// z coord of Z_test is now expected to be 0, i.e.
// all following operations take place in the x,y plane

// rotate Z_test by theta*(b-1.0)
Z_test = Z_test.rotz(theta * (b-1.0));

// scale magnitude of Z_test with magnitude raised to the power b-1
// so the result is the original length raised to the power b
Z_test = Z_test * pow(magnitude, (b-1.0));

// reverse the previous rotations
Z_test = Z_test.roty(N_theta);
Z_test = Z_test.rotz(N_phi);

return Z_test;
}

// return the pole vector (which must have unit length)
static triplex pole()
{
triplex result = {1.0, 1.0, 1.0};
return result / result.len();
}
};

// perform a loop over x,y,z here and fix how the result is stored according to your needs.

triplex Z = {x,y,z};
triplex C = Z;

double mag;
int iter = 0;

// iteration loop. This couldn't be any easier, right?
while((mag = Z.len()) < 2.0)
{
Z = (Z ^ 2.0)  + C;

if (++iter == 100) break;
}

// store the result in a volume texture. This has been tweaked to result
// in the following graphically displayed right handed coordinate system
// x points right
// y points "away"
// z points up
mem[volumeSize.width*volumeSize.height*(volumeSize.depth-1-iy) + volumeSize.width*iz + ix] = iter;

I would be in particular interested in a non-trigonometric approach to this, i.e. a solution to the equations

Code:
1) x^2 + y^2 + z^2 = (power * |Z|)^2
2) x*px + y*py + z*pz = power * |Z| * cos(power * theta)
3) x*nx + y*ny + z*nz = 0

theta is the angle between Z and the pole vector P=[px,py,pz]
N=[nx,ny,nz] is the normal vector  (normalized cross product of pole vector and Z)

[x,y,z] are then the resulting coordinates of Z^power

The first equation describes the required length of the resulting vector Z^power
The second equation describes the required angle between the pole vector and Z^power
The third equation describes the plane in which the solution needs to fall into.

Note that x,y,z are the components of the Z^power vector which we try to determine. Of the original vector Z we only
need its length |Z| which is known.

I would expect this to return two results (roots of a quadratic equation), one of which is the point I need.

Can someone help me out with Mathematica or Maple and solve this equation system for x,y,z please? I tried to do this on a sheet of paper but it got too messy and so I went for implementation of the trigonometric version first.
 « Last Edit: December 12, 2009, 03:15:54 PM by cbuchner1 » Logged
cbuchner1
Fractal Phenom      Posts: 443 « Reply #3 on: December 11, 2009, 12:33:09 AM »

Code:
1) x^2 + y^2 + z^2 = (power * |Z|)^2
2) x*px + y*py + z*pz = power * |Z| * cos(power * theta)
3) x*nx + y*ny + z*nz = 0

No more need for help here. I got myself a free 15 day trial of Mathematica and it helped me to solve equation 2) and 3) as a function of the z coordinate. This in turn I can insert into equation 1) and get a solution. The harder part will be to reject the "wrong" solution of the quadratic equation. hmm....

Christian
 « Last Edit: December 11, 2009, 12:37:34 AM by cbuchner1 » Logged
cbuchner1
Fractal Phenom      Posts: 443 « Reply #4 on: December 11, 2009, 08:25:51 PM »

The posted rotation code is buggy. It should not compute N_theta based on the original N, but on z-axis rotated N around by am amount of -N_phi.

Code:
// rotate Z space such that normal vector points "up" (in z direction)
double N_phi   = atan2(N.y, N.x);
double N_theta = atan2(-N.x, N.z);
Z_test = Z_test.rotz(-N_phi);
Z_test = Z_test.roty(-N_theta);

Correct code looks like this

Code:
// rotate Z space such that normal vector points "up" (in z direction)
double N_phi   = atan2(N.y, N.x);
triplex Nxz = N.rotz(-N_phi);                // insert this line for correct rotation
double N_theta = atan2(-Nxz.x, Nxz.z); // and base N_theta calculation on above result
Z_test = Z_test.rotz(-N_phi);
Z_test = Z_test.roty(-N_theta);

As an effect of this bug you get nicely asymmetrical 3D Mandelbrots.

And the result of the bugfix will be a boring symmetrical fractal, lathed simply around its axis defined by the pole vector.
 « Last Edit: December 12, 2009, 03:17:27 PM by cbuchner1 » Logged
 Pages:    Go Down
 Related Topics Subject Started by Replies Views Last post  A very strange idea; vector to the Mandelbrot Let's collaborate on something! Timeroot 8 3222 January 23, 2010, 07:44:58 AM by Timeroot  The South Pole Mandelbulber Gallery mclarekin 0 478 January 10, 2014, 01:19:30 PM by mclarekin  The Mandelbrot Set. Perturbation Theory. UltraFractal SeryZone 0 3192 October 14, 2014, 07:54:14 PM by SeryZone  Fastest mandelbrot software? Perturbation method? Fractal Programs « 1 2 » recursiveidentity 25 6564 August 02, 2017, 04:19:11 AM by Dinkydau  North Pole Images Showcase (Rate My Fractal) Lois 1 1315 November 24, 2017, 03:11:26 AM by Lois