biComplexMod5
preset values

Hypercomplex fractals

3D/4D fractals explored with stochastic IIM (Inverse Iterative Method)

Abuot Inverse Iteration Method:

Below I try to explain, in a very elementary way, the IIM (Inverse Iteration Method), and how I use it in glChAoS.p: this is important to know how to use the parameters

The classic Julia set is repulsive with respect to the polynomial fc(z) = z^2 + c, which means that the orbits of points near the set move farther away from it. With the more known DE (Distance Estimator) method, we explore, point by point a region of C space, and assign a color to the speed with which the point flee from orbit. Conversely, the Julia set is "attractive" with respect to the inverse function fc^-1 = sqrt(z-c), because the orbits of points far away from the Julia set come closer to it With the IIM, we choose a initial point on the C space, and we iterate until we reach the orbit.
There are 2 things to observe:
1) the fc^-1 = sqrt(z-c) formula produce two results fx^-1 = +sqrt(z-c) and fy^-1 = -sqrt(z-c) (sqrt property)
So if z0 is initial point { fx^-1(z0), fy^-1(z0) } are two points that will be closer to the Julia set.
From these, with next iteration, we have { fx^-1(fx^-1(z0)), fx^-1(fy^-1(z0)), fy^-1(fx^-1(z0)), fy^-1(fy^-1(z0)) }...
... and so on.
After k iteration, we have 2^(k-1) results: and a typical binary tree, of deeph = k-1, to explore.
glChAoS.P uses a stochastic approach obtained by choosing one of the two roots, random, at each stage of the recursion up to a determined deeph (selectable), and when we reach it, restart the iteration from starting point.
Statistically, sooner or later, we explore all the branches of the binary tree.
Formally, this can be considered a typical example for an iterated function system (IFS).
2) not all Julia sets generated by a single starting point can be visualized by the inverse iteration method, and a single point may not be enough to define the full object
In glChAoS.P, when the depth has been reached, the re-iteration restart from starting point plus a "possible" random displacement, if are defined "min" and "max" values.
For example, if starting point is p(0,0,0) and min = -1.0 and max = 1.0, the new starting point will be randomly generated in a sphere with center p(0,0,0) and radius = 1. In this way there are more possibilities to define the full set, but are also created some "garbage" starting points on the space: you can visually find a balance.


Quaternion Julia

    qJulia(const vec4 &p) {
        const int rnd = Random::get(int(0),int(INT_MAX));
        const float sign = (rnd&1) ? 1.f : -1.f)
        const float xQ = p.x * p.x, yQ = p.y * p.y, zQ = p.z * p.z, wQ = p.w * p.w;
        const float r = sqrtf(xQ + yQ + zQ + wQ);
        const float a = sqrtf((p.x+r)*.5);
        const float b = (r-p.x) * a / (yQ + zQ + wQ);
        pNew = sign * vec4(a, b*p.y, b*p.z , b*dim4D );
    }

 


juliaBulb

auto JuliaBulb(const vec3 &p, float sign1, float sign2) {
    const float xQ = p.x * p.x, yQ = p.y * p.y, zQ = p.z * p.z;
    float r = sqrtf(xQ+yQ+zQ);
    float a = sign1 * sqrtf(yQ + xQ*(2.f+zQ/(xQ+yQ)) - 2.f*p.x*r);
    float b = sign2 * sqrtf(r - p.x - a) * .5f;
    float c = yQ*yQ + xQ * (yQ - zQ);
    float d = a * (p.x * (r+p.x) + yQ);
    vp = vec3(b * ((p.x*yQ-d) * (xQ+yQ) - p.x*xQ*zQ) / (p.y*c), b, -p.z/sqrtf(2 * (r - d *(xQ+yQ)/c)));
}

const int rnd = Random::get(0,INT_MAX);
radiciEq(p - c, (rnd&1) ? 1.f : -1.f, (rnd&2) ? 1.f : -1.f);

 



JuliaBulb degree Nth

    JuliaBulbNth(const vec3 &p, int kTheta, int kPhi) {
        const float r = length(p);
        const int absOrder = abs(degreeN);
        const int k1 = (absOrder - (p.z<0 ? 0 : 1)) * .25, k2 = (3*absOrder + (p.z<0 ? 4 : 2)) * .25;

        const int dk = ((absOrder%2)==0 && kPhi>k1 && kPhi()) / float(degreeN);
        const float phi    = (asinf(p.z/r)    + (2 * kPhi   - dk) * glm::pi()) / float(degreeN);
        const float cosphi =cosf(phi);
        pNew = powf(r, 1.0f/float(degreeN)) * vec3(cosf(theta)*cosphi,sinf(theta)*cosphi,sinf(phi));
    }

    JuliaBulbNth(p - c, Random::get(0, INT_MAX) % degreeN, Random::get(0, INT_MAX) % degreeN);

 



Bicomplex

    Bicomplex(const vec4 &p, float sign1, float sign2)
    {
        const std::complex z1(p.x, p.y), z2(-p.w, p.z);
        const std::complex w1 = sign1 * sqrt(z1 - z2), w2 = sign2 *sqrt(z1 + z2);
        pNew = .5f * vec4(w1.real()+w2.real(), w1.imag()+w2.imag(), w2.imag()-w1.imag(), w1.real()-w2.real());
    }

    const int rnd = Random::get(int(0),int(INT_MAX));
    Bicomplex(p-c, (rnd&1) ? 1.f : -1.f, (rnd&2) ? 1.f : -1.f);

 

Previous formulas are a stochastic adaptation of Paul Nylander's Wolfram Mathematica code.

 

 

The following formulas were born by chance from BiComplex original formula:
I had not implemented yet 4D iteration, so I have set w parameter to start value, with no increment in iterated formula.
What appeared was pleasant, so I included them in the program


Bicomplex modified formula 0

    BicomplexMod0(const vec4 &p, float sign1, float sign2)
    {
        const std::complex z1(p.x, p.y), z2(-start.w, p.z);
        const std::complex w1 = sign1 * sqrt(z1 - z2), w2 = sign2 *sqrt(z1 + z2);
        pNew = .5f * vec4(w1.real()+w2.real(), w1.imag()+w2.imag(), w2.imag()-w1.imag(), w1.real()-w2.real());
    }

    const int rnd = Random::get(int(0),int(INT_MAX));
    BicomplexMod0(p - c, (rnd&1) ? 1.f : -1.f, (rnd&2) ? 1.f : -1.f);

 


Bicomplex modified formula 1

    BicomplexMod1(const vec4 &p, float sign1, float sign2)
    {
        const std::complex z1(p.x, p.y), z2(-p.w, start.z);
        const std::complex w1 = sign1 * sqrt(z1 - z2), w2 = sign2 *sqrt(z1 + z2);
        pNew = .5f * vec4(w1.real()+w2.real(), w1.imag()+w2.imag(), w2.imag()-w1.imag(), w1.real()-w2.real());
    }

    const int rnd = Random::get(int(0),int(INT_MAX));
    BicomplexMod1(p - c, (rnd&1) ? 1.f : -1.f, (rnd&2) ? 1.f : -1.f);