The code was greatly simplified using only quaternion maths. Step 1) directly rotates Z around the axis defined by PxZ and no longer performs an intermediate rotation into the x,y plane and back . I have replaced the C++ code in the parent posting to use this approach.
For "fat bird" and the new bulbs posted in first post, P,Q,R should be {0,0,1} {0,1,0} {0,0,1}
For "plancton" P,Q,R can be {1,0,0},{0,1,0},{0,0,1}
I implemented three poles P(), Q() and R() as static functions. If you choose arbitrary axes, make sure P,Q,R have unit length (i.e. divide vectors by their norm).
You have 5 degrees of freedom for chosing P(), Q() and R() as long as you keep Q() orthogonal to R() (I think my implementation currently requires Q to be orthogonal to R or the calculations will be incorrect). However three degrees of freedom thereof effectively rotate the entire object in space. So the two relevant degrees of freedom are:
- how P relates to R (angle between these vectors)
- how Q is situated relative to P ( with the additional constraint Q _|_ R )
Now we can start to think about simplifying the entire computation.
I propose fixing R at {0,0,1}, Q at {1,0,0}, therefore Step 2) becomes the equivalent of the phi rotation in a spherical coordinate system of the "traditional" Mandelbulb. This simplifies Step 2) greatly (just rotate the x,y coordinates - ignore z).
Our two available degrees of freedom are therefore achieved by arbitrarily pointing the pole vector P on the unit sphere. The straightforward quaternion maths of Step 1) remain unchanged.
Here is an attempt at optimizing the code for fixed R,Q (as proposed) but with an arbitrary P vector. For integer powers 2-9 no trigonometry will be used. The most expensive computations are three sqrt() expressions per iteration (which could be trimmed to two, I suppose, by remembering the magnitude as member of the triplex object between iterations).
The more generic and simple code is in this posting
http://www.fractalforums.com/mandelbulb-renderings/a-new-class-of-bulb/msg9861/#msg9861 // euclidean vector utility class with the required operator overloads
struct triplex
{
double x,y,z;
// magnitude (length) of this object
inline double len() const { return sqrt(x*x+y*y+z*z); }
// vector addition between this object and another vector b
inline triplex operator + (const triplex &b) const
{
triplex result = {x + b.x, y + b.y, z + b.z};
return result;
}
// vector subtraction between this object and another vector b
inline triplex operator - (const triplex &b) const
{
triplex result = {x - b.x, y - b.y, z - b.z};
return result;
}
// vector division by a scalar b
inline triplex operator / (const double &b) const
{
triplex result = { x / b, y / b, z / b };
return result;
}
// vector multiplication with a scalar b
inline triplex operator * (const double &b) const
{
triplex result = { x * b, y * b, z * b };
return result;
}
// dot product between this object and another vector b
inline double operator * (const triplex &b) const
{
return x * b.x + y * b.y + z * b.z;
}
// cross product between this object and another vector b
inline triplex operator % (const triplex &b) const
{
triplex result = { y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x };
return result;
}
// return result of rotation around z axis
inline triplex rotz(double arg) const
{
triplex result = { cos(arg)*x - sin(arg)*y,
sin(arg)*x + cos(arg)*y,
z };
return result;
}
// quaternion helper class describing rotations
struct quaternion
{
// 1,i,j,k
double a,b,c,d;
// rotation around a given axis (given sine and cosine of HALF the rotation angle)
static inline quaternion describe_rotation(triplex v, double sina_2, double cosa_2)
{
quaternion result = { cosa_2, sina_2*v.x, sina_2*v.y, sina_2*v.z };
return result;
}
// rotation around a given axis (angle without range restriction)
static inline quaternion describe_rotation(triplex v, double angle)
{
double sina_2 = sin(angle/2);
double cosa_2 = cos(angle/2);
quaternion result = { cosa_2, sina_2*v.x, sina_2*v.y, sina_2*v.z };
return result;
}
// rotate a triplex in 3D space with the given quaternion
// see EN Wikipedia on Quaternions and spatial rotation
inline triplex rotate(triplex v) const
{
double t2 = a*b;
double t3 = a*c;
double t4 = a*d;
double t5 = -b*b;
double t6 = b*c;
double t7 = b*d;
double t8 = -c*c;
double t9 = c*d;
double t10 = -d*d;
triplex result = {
2*( (t8 + t10)*v.x + (t6 - t4)*v.y + (t3 + t7)*v.z ) + v.x,
2*( (t4 + t6)*v.x + (t5 + t10)*v.y + (t9 - t2)*v.z ) + v.y,
2*( (t7 - t3)*v.x + (t2 + t9)*v.y + (t5 + t8)*v.z ) + v.z
};
return result;
}
};
// quaternion with a fixed rotational axis of {0,0,1}
// simplified to only include the required variable members
// and to perform only a minimal amount of computations
struct zternion
{
// 1,k
double a,d;
// rotation around a given axis (given sine and cosine of HALF the rotation angle)
static inline zternion describe_rotation(double sina_2, double cosa_2)
{
zternion result = { cosa_2, sina_2 };
return result;
}
// rotation around a given axis (angle without range restriction)
static inline zternion describe_rotation(double angle)
{
double sina_2 = sin(angle/2);
double cosa_2 = cos(angle/2);
zternion result = { cosa_2, sina_2 };
return result;
}
// rotate a triplex in 3D space with the given zternion
// see EN Wikipedia on quaternions and spatial rotation
inline triplex rotate(triplex v) const
{
double t4 = a*d;
double t10 = -d*d;
triplex result = {
2*( t10 * v.x - t4 * v.y ) + v.x,
2*( t4 * v.x + t10 * v.y ) + v.y,
v.z
};
return result;
}
};
// raise this object to the given power relative to the pole vector
triplex operator ^ (const double &power) const
{
// storing the magnitude saves us multiple calls to len()
double magnitude = len();
// compute sine/cosine expressions of the angle theta between current vector
// and the pole vector - none of this requires actual sine/cosine calls.
// The most expensive function here is the sqrt()
double costheta = (*this * P()) / magnitude;
double costheta2 = costheta * costheta;
double sintheta2 = 1.0 - costheta2;
double sintheta = sqrt(sintheta2);
// compute normal vector that defines the plane in which the result shall live
triplex N = (P() % *this) / (magnitude * sintheta);
// start by creating a copy of this object
triplex Z_power = *this;
// Step 1)
// generate a result vector that is situated in the plane
// spanned by origin, pole and this object and has an resulting
// angle power*theta to the pole and is of the required magnitude.
{
int method = 0;
double scale;
quaternion rot;
// Handle accelerated cases for some integer powers (2-9)
// that do not require any extra trigonometry to work.
if (power == 2.0 || power == 3.0)
{
if (power == 2.0)
method = 1; // rotate P by angle arg=2*theta
else method = 2; // rotate Z by angle arg=2*theta
// sin(arg/2) = sin(theta)
// cos(arg/2) = cos(theta)
rot = quaternion::describe_rotation( N, sintheta, costheta );
scale = magnitude * magnitude; // magnitude^2
}
else if (power == 4.0 || power == 5.0)
{
if (power == 4.0)
method = 1; // rotate P by angle arg=4*theta
else method = 2; // rotate Z by angle arg=4*theta
// sin(arg/2) = sin(2*theta)=2*sin(theta)*cos(theta)
// cos(arg/2) = cos(2*theta)=cos^2(theta)-sin^2(theta)
rot = quaternion::describe_rotation( N, 2*sintheta*costheta, costheta2-sintheta2 );
double tmp = magnitude*magnitude; scale = tmp*tmp; // magnitude^4
}
else if (power == 6.0 || power == 7.0)
{
if (power == 6.0)
method = 1; // rotate P by angle arg=6*theta
else method = 2; // rotate Z by angle arg=6*theta
// sin(arg/2) = sin(3*theta)=3*sin(theta)-4*sin^3(theta)
// cos(arg/2) = cos(3*theta)=4*cos^3(theta)-3*cos(theta)
rot = quaternion::describe_rotation( N, 3*sintheta-4*sintheta2*sintheta,
4*costheta2*costheta-3*costheta );
double tmp = magnitude*magnitude; scale = tmp*tmp*tmp; // magnitude^6
}
else if (power == 8.0 || power == 9.0)
{
if (power == 8.0)
method = 1; // rotate P by angle arg=8*theta
else method = 2; // rotate Z by angle arg=8*theta
// sin(arg/2) = sin(4*theta)=8*cos^3(theta)*sin(theta)-4*cos(theta)*sin(theta)
// cos(arg/2) = cos(4*theta)=8*cos^4(theta)-8*cos^2(theta)+1
rot = quaternion::describe_rotation( N, 8*costheta2*costheta*sintheta-4*costheta*sintheta,
8*costheta2*costheta2-8*costheta2+1 );
double tmp = magnitude*magnitude; tmp = tmp*tmp; scale = tmp*tmp; // magnitude^8
}
switch(method)
{
case 1:
// scale and rotate P around axis N
Z_power = rot.rotate(P()) * scale;
break;
case 2:
// scale and rotate Z around axis N
Z_power = rot.rotate(Z_power) * scale;
break;
default:
// fall back to slower trigonometric method
// scale and rotate P around axis N
rot = quaternion::describe_rotation( N, power*acos(costheta) );
Z_power = rot.rotate(P()) * pow(magnitude, power);
break;
}
}
// Step2)
// Perform rotation around {0,0,1} z-axis by an angle of (power-1)*phi
{
double xy_len = sqrt(Z_power.x*Z_power.x+Z_power.y*Z_power.y);
double sinphi = Z_power.y / xy_len, sinphi2 = sinphi*sinphi;
double cosphi = Z_power.x / xy_len, cosphi2 = cosphi*cosphi;
int method = 0;
zternion zrot; // simplified quaternion with fixed rotational axis {0,0,1}
// Handle accelerated cases for some integer powers (2-9)
// that do not require any extra trigonometry to work.
if (power == 2.0 || power == 3.0)
{
if (power == 2.0)
method = 1; // rotate P by angle arg=2*phi
else method = 2; // rotate Z by angle arg=2*phi
// sin(arg/2) = sin(phi)
// cos(arg/2) = cos(phi)
zrot = zternion::describe_rotation( sinphi, cosphi );
}
else if (power == 4.0 || power == 5.0)
{
if (power == 4.0)
method = 1; // rotate P by angle arg=4*phi
else method = 2; // rotate Z by angle arg=4*phi
// sin(arg/2) = sin(2*phi)=2*sin(phi)*cos(phi)
// cos(arg/2) = cos(2*phi)=cos^2(phi)-sin^2(phi)
zrot = zternion::describe_rotation( 2*sinphi*cosphi, cosphi2-sinphi2 );
}
else if (power == 6.0 || power == 7.0)
{
if (power == 6.0)
method = 1; // rotate P by angle arg=6*phi
else method = 2; // rotate Z by angle arg=6*phi
// sin(arg/2) = sin(3*phi)=3*sin(phi)-4*sin^3(phi)
// cos(arg/2) = cos(3*phi)=4*cos^3(phi)-3*cos(phi)
zrot = zternion::describe_rotation( 3*sinphi-4*sinphi2*sinphi,
4*cosphi2*cosphi-3*cosphi );
}
else if (power == 8.0 || power == 9.0)
{
if (power == 8.0)
method = 1; // rotate P by angle arg=8*phi
else method = 2; // rotate Z by angle arg=8*phi
// sin(arg/2) = sin(4*phi)=8*cos^3(phi)*sin(phi)-4*cos(phi)*sin(phi)
// cos(arg/2) = cos(4*phi)=8*cos^4(phi)-8*cos^2(phi)+1
zrot = zternion::describe_rotation( 8*cosphi2*cosphi*sinphi-4*cosphi*sinphi,
8*cosphi2*cosphi2-8*cosphi2+1 );
}
switch(method)
{
case 1:
{ // rotate a vector facing direction {1,0,0} (direction of zero angle)
triplex Z_P = {xy_len,0,Z_power.z};
Z_power = zrot.rotate(Z_P); }
break;
case 2:
// rotate Z_power vector
Z_power = zrot.rotate(Z_power);
break;
default:
// generic rotation (use slow trigonometry)
Z_power = Z_power.rotz( (power-1)*atan2(Z_power.y, Z_power.x) );
break;
}
}
return Z_power;
}
// return the P pole vector which must have unit length
// results in a lathed Mandelbrot around axis P during Step 1
static inline triplex P()
{
triplex result = {0, 1, 0};
// return result / result.len(); // use this if the above is not unit lenth
return result;
}
};
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 ^ 3.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;
To get an actual *formula* out of it, pick a case (e.g. power=2.0) and see what is computed. Probably a software package like Mathematica would be very helpful to unravel this. You could simplify this even more by fixing also the vector P but you would lose all the degrees of freedom in your parameter space.
Christian