Logo by mrob - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Visit us on facebook
 
*
Welcome, Guest. Please login or register. March 29, 2024, 03:48:18 PM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1] 2   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: A new class of bulb?  (Read 5696 times)
Description: see for yourself
0 Members and 1 Guest are viewing this topic.
cbuchner1
Fractal Phenom
******
Posts: 443


« on: December 12, 2009, 06:20:39 PM »

Hi,

this method is comprised of two steps to raise a point Z=(x,y,z) to a given power in 3D space.

Step 1)

First we perform a 2D complex operation inside the plane spanned up by the coordinate origin, the end point of the pole vector P (an arbitrary pointed unit vector) and the point Z. The easiest way to do this is to determine the angle theta between P and Z and to simply raise this angle to the required power i.e. rotating Z by (power-1)*theta around the axis N=PxZ. We also scale the resulting vector with the scalar |Z|^(power-1).

When you perform only this step for power 2 you get a boring rotated Mandelbrot ("lathed" mandelbrot) - because each plane that contains the pole vector also contains the 2D Mandelbrot.

Step 2)

Choose another unit vector Q perpendicular to another pole vector R. Now I rotate the previously resulting point Z around the pole axis by an amount of ((power-1) * arg), where arg is defined as the angle 0...2*PI between the vector Q and the projection of the vector Z onto the plane that is normal to R. Q and the origin are by definition contained in this plane.

EDIT: P can be identical to R, or one can also choose a different rotational axis R in Step 2. This creates much larger variety of objects depending on how these axes are situated to each other.



To generate the bulb, iterate over the z^power+c Mandelbrot-like formula until the norm of Z exceeds radius 2. Repeat for every point in space within a certain distance (e.g. 2) from the origin.

Do you think this creates a new class of Mandelbulb or has this been covered in any of Bugman's generalisations? Step 1) is an equivalent of scaling theta and the magnitude (norm) of Z in spherical coordinates, whereas step 2) performs an operation on the angle phi (but not in spherical coordinates, but rather cylindrical ones)

What you see is a power 6 bulb where P=R








Christian
« Last Edit: December 15, 2009, 05:51:00 PM by cbuchner1 » Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #1 on: December 12, 2009, 06:43:47 PM »


And a power 8





Logged
bugman
Conqueror
*******
Posts: 122



WWW
« Reply #2 on: December 12, 2009, 09:03:12 PM »

Can you express this power formula in {x, y, z} notation?
Logged
twinbee
Fractal Fertilizer
*****
Posts: 383



WWW
« Reply #3 on: December 12, 2009, 11:13:12 PM »

Fascinating - wish I could see the object more clearly. Have you tried rendering the quadratic version?
Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #4 on: December 12, 2009, 11:24:57 PM »

Fascinating - wish I could see the object more clearly. Have you tried rendering the quadratic version?

The quadratic version looks a bit like a fat bird with etremely short wings wink
Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #5 on: December 13, 2009, 02:01:46 PM »


The quadratic version looks a bit like a fat bird with etremely short wings wink


After careful consideration and comparison with bugmans renders I conclude that this is indeed a new kind of bulb. It consists of two chained operations (step 1 and step 2) with cylindrical coordinates or symmetries - and does not operate on spherical coordinates.

For reference I am including a piece of C++ code (minus the outer loop across voxel space x,y,z) to reproduce the "fat bird". Get higher power bulbs by adjusting the exponent inside the iteration loop. Both the implementation of step1 and step2 are based on quaternions because these greatly simplify rotations around arbitrary axes.

I am sure this could be simplified into a single formula like {x,y,z} = f(x,y,z), in particular when choosing pole vectors along the axes. But right now I prefer to try new combinations and operations instead of simplifying the existing ones.

By the way I am hoping for the zoomed fractal detail of this object to be less stretched in some places because the transformations do not have the same kind of distortions near the poles as created by the spherical coordinate system.

Code:
                   // 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;
                        }

                        // quaternion helper class describing rotations
                        struct quaternion
                        {
                            //     1,i,j,k
                            double a,b,c,d;

                            // 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;
                            }
                        };

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

                            {
                                // NOTE: for integer powers 2,3,4... we will probably find formulas easily
                                // that do not require trigonometry at all.
                                // However the quaternion expects us to put in sine and cosine
                                // of HALF the rotation angle, which is no fun. The formulas for
                                // sin,cos(angle/2) as a function of cos(angle) have square roots in
                                // them and the sign needs to be inverted according to which quadrant
                                // the half angle is in.
                                double theta = acos(costheta);

                                // rotation axis is N, amount is (power-1)*theta
                                quaternion rot = quaternion::describe_rotation( N, (power-1)*theta );

                                // rotate and scale Z_power
                                Z_power = rot.rotate(Z_power) * pow(magnitude, power-1);

                                // NOTE: Equivalent transformation (maybe much easier to simplify)
                                // rotate P by power*theta and scale with magnitude^power
                                // quaternion rot = quaternion::describe_rotation( N, power*theta );
                                // Z_power = rot.rotate(P()) * pow(magnitude, power);
                            }

                            // Step 2)
                            // Perform rotation around pole axis by an angle of (power-1)*phi
                            {
                                // orthogonal projection of Z_power onto the plane normal to R
                                triplex Z_p = R() * (Z_power * R());
                                triplex Z_normal = Z_power - Z_p;

                                // determine angle between Q and Z_normal (Q is assumed to be in the plane already)
                                // and this into a range of phi = 0...2PI with respect to Q, R
                                double cosphi = (Z_normal * Q()) / Z_normal.len();
                                double phi = acos(cosphi);
                                if ((Q() % Z_normal) * R() < 0)
                                    phi = M_PI*2-phi;

                                // rotate by an angle of (power-1) * phi around the pole axis R
                                quaternion rot2 = quaternion::describe_rotation( R(), power*phi-phi );
                                Z_power = rot2.rotate(Z_power);
                            }

                            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, 0, 1};  // use {1,0,0} to get the Plancton bulbs
//                            return result / result.len();  // use this if the above is not unit lenth
                            return result;
                        }

                        // return the Q pole vector which must have unit length and be perpendicular to R
                        // defines the angular basis for rotation around vector R
                        static inline triplex Q()
                        {
                            triplex result = {0, 1, 0};
//                            return result / result.len();  // use this if the above is not unit lenth
                            return result;
                        }

                        // return the R pole vector
                        // rotational axis during Step 2. Can be identical to P but does not have to be.
                        static inline triplex R()
                        {
                            triplex result = {0, 0, 1};
//                            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 ^ 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;
« Last Edit: December 15, 2009, 01:30:51 PM by cbuchner1 » Logged
TedWalther
Guest
« Reply #6 on: December 13, 2009, 11:16:29 PM »

Wow, I'd love to see some of twinbees high-res renders on this; as soon as I saw your power-6, it looked exactly like an almond/apple blossom.  What does pow-3 look like?

Here is another idea I've been thinking about, but haven't figured out how to do it.  What about volume?  With the 2-D version of mandelbrot, you could say that the area swept out by the point Z is doubled, then the magnitude of Z is doubled.  How to extend that to 3-D so that the volume between Z, the X axis, and the Z axis, is doubled, then the resulting Z has its magnitude doubled?

Ted
Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #7 on: December 14, 2009, 03:33:15 AM »

Wow, I'd love to see some of twinbees high-res renders on this; as soon as I saw your power-6, it looked exactly like an almond/apple blossom.  What does pow-3 look like?

Unfortunately High res and voxels are a bit contradictory (unless you think of the great GigaVoxels project on the GPU, which allows one to dynamically zoom into the voxel space which is re-computed on demand).

But after simplifying the formulas they will most likely allow for distance estimation and raymarching approaches. Right now my formulas are a bit too generic and it would be a major pain to find the correct derivative term needed for raymarching.

Oh wow. The objects get even more interesting when I use two separate rotational axes.

Step1) rotates around P(), e.g. [1,0,0]
Step2) rotates around R(), e.g. [0,0,1] where the vector Q() e.g [0,1,0] defines the angular reference for the rotation.

I am getting the most amazing voxel objects out of that. Ready for this? And that's just one viewing angle.

This is a power 3 with above axes



And here is a power 4



A power 5... Is this thing alive?



There you have it. It's voxel plancton.

I am beginning to see that the choice of spherical coordinates for the "original" Mandelbulb was seriously limiting the variety of bulb "species" that you can get out of it - the spherical coords are using a fixed set of rotational axes. Making these axes freely definable is what will create a much larger variety of objects.

I really need to shift my voxel computation to the GPU. Right now I am waiting 5 minutes for a 512x512x512 object to be generated on the CPU.

Christian
« Last Edit: December 14, 2009, 04:26:12 AM by cbuchner1 » Logged
TedWalther
Guest
« Reply #8 on: December 14, 2009, 07:46:49 AM »

Now if only there was some way to modify Q at every iteration...
Logged
kram1032
Fractal Senior
******
Posts: 1863


« Reply #9 on: December 14, 2009, 03:11:38 PM »

really nice ULOs you have there cheesy
Unknown Living Objects smiley

Actually, they could just as well be UFOs... smiley

If you get a realtime-ish approach out of this, it would be cool to see animated axes smiley
Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #10 on: December 14, 2009, 03:26:27 PM »

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

Code:
                   // 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
« Last Edit: December 15, 2009, 01:35:38 PM by cbuchner1 » Logged
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #11 on: December 15, 2009, 02:48:43 AM »

A CUDA enabled plancton app was published here: http://forums.nvidia.com/index.php?showtopic=153257

It takes half a second on a nVidia GTX 260 to render a fractal with new parameters. Rotating, zooming and panning however is interactive frame rates. Fractal parameters can be controlled with keypresses. The Mandelbulbs from the top posting in this thread appear when theta is set to 0 or 180 degrees, i.e. the P vector is pointing upwards or downwards (I use spherical coordinates to control where the P vector is pointing). Various plancton-like creatures lurk in between.

I even found a Donald Duck somewhere wink



Christian


« Last Edit: December 15, 2009, 01:28:57 PM by cbuchner1 » Logged
David Makin
Global Moderator
Fractal Senior
******
Posts: 2286



Makin' Magic Fractals
WWW
« Reply #12 on: December 15, 2009, 02:05:24 PM »

Hi Christian, would I be correct in assuming that you are simply applying two rotations in the Mandelbulb manner but the rotations concerned are around user-definable axes ?
In which case if you want the axes concerned to be perpendicular than you can simply transform x,y,z by a standard 3D rotation then apply the Mandelbulb calculation then rotate back by the inverse standard 3D rotation - is that correct ?
Obviously if you want the two axes concerned to be independant then it's slightly more complicated but can be done essentially the same way.

Have you tried rotation around independant axes that do not go through the origin ?
« Last Edit: December 15, 2009, 02:07:19 PM by David Makin » Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
cbuchner1
Fractal Phenom
******
Posts: 443


« Reply #13 on: December 15, 2009, 03:44:58 PM »

Hi Christian, would I be correct in assuming that you are simply applying two rotations in the Mandelbulb manner but the rotations concerned are around user-definable axes ?
In which case if you want the axes concerned to be perpendicular than you can simply transform x,y,z by a standard 3D rotation then apply the Mandelbulb calculation then rotate back by the inverse standard 3D rotation - is that correct ?

Not quite. The first step is in essence the Mandelbulb formula with a theta rotation only - however I believe my definition of theta differs somewhat from the definition in a spherical coordinate system (the rotation is only identical when phi is 0).  The second step is only rotating phi (cylindrical rotation). The rotational axes are not necessarily perpendicular. Note that in the spherical coordinate system there is a convergence towards the poles (due to x=cos(theta)*cos(phi), y=cos(theta)*sin(phi), z=sin(theta) where cos(theta) becomes 0 for theta converging towards +/-PI. In my two-step approach I avoid this problem. So my resulting Mandelbulbs look somewhat different (see first post)

A simplification is achieved by making one of the rotation axes align with a coordinate axis. In my simplification approach I've taken I've aligned the phi rotation with the z axis, but the first step then requires quaternion maths to do the theta rotation for an arbitrary axis (which defines the angle theta against Z)

Have you tried rotation around independant axes that do not go through the origin ?

Not yet, need to brush up my math skills first and learn if it can be done with some help by quaternions, too. wink

Christian
« Last Edit: December 15, 2009, 08:35:52 PM by cbuchner1 » Logged
BradC
Safarist
******
Posts: 85



« Reply #14 on: December 15, 2009, 05:18:42 PM »

Have you tried rotation around independant axes that do not go through the origin ?
Not yet...

I believe you can do the following: To rotate a point (x,y,z) by an angle \theta about an axis aligned with a vector \overset{\rightharpoonup }{p} but passing through point (a,b,c), first subtract (a,b,c) from your point (x,y,z) to move the axis so that it passes through the origin, then do the quaternion q \overset{\rightharpoonup }{v} q^{-1} step with q=\cos \left(\frac{\theta }{2}\right)+\overset{\rightharpoonup }{p} \sin \left(\frac{\theta }{2}\right) as usual, and then finally add (a,b,c) back in again to move everything back where it should be.

BTW, I really like this class of fractals you've come up with. I'm looking forward to exploring them some. smiley
« Last Edit: December 15, 2009, 05:26:41 PM by BradC, Reason: Oops, I had two different variables both called p. » Logged
Pages: [1] 2   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
P 8 bulb inside a P 2 bulb. probably not... Mandelbulb Renderings KRAFTWERK 12 4301 Last post March 08, 2010, 09:45:21 AM
by kram1032
A Klingon D7-class battlecruiser Images Showcase (Rate My Fractal) thom 0 1233 Last post April 27, 2012, 11:32:09 PM
by thom
Epirial Pizza Delivery, Dreadnought Class Images Showcase (Rate My Fractal) John Smith 0 1089 Last post August 04, 2012, 02:42:58 AM
by John Smith
class at itztapalapa - mexiko Fractal Forums News cKleinhuis 1 1112 Last post October 24, 2012, 03:42:25 AM
by LMarkoya
Sudarsky Class I Images Showcase (Rate My Fractal) Pauldelbrot 4 1354 Last post March 10, 2014, 04:29:17 PM
by Pauldelbrot

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.328 seconds with 30 queries. (Pretty URLs adds 0.011s, 2q)