This is a diff of the original fractal.cpp from mandelbulber 1.08 and the generalized box fold changes. The code essentially replaces the Tglad mandelbox box fold with the generalized fold. I had trouble implementing the "folding limit" function, so I have used this parameter to select between a number of polyhedron.
Folding limit polyhedron
default tet
2 cube
3 oct
4 dodeca
5 turncated oct.
Disclaimer. The cube selection does not reproduce the original box fold. This points to a bug in either the theory or the implementation. However it produces some interesting fractals as is. I expect to find the bug when I try to implement the melting poly version suggested by Hobold.
I would be honored if Budda wants to integrate this into the mandelbulber, but he might want to wait a couple of days to see if I can fix the bug and melt some polys.
I did not realize that the code option compressed the file to a nice little window. Here is the entire fractower.cpp with the generalized box fold changes.
/********************************************************
/ MANDELBULBER *
/ *
/ author: Krzysztof Marczak *
/ contact: buddhi1980@gmail.com *
/ licence: GNU GPL *
********************************************************/
/*
* fractal.cpp
*
* Created on: 2010-01-23
* Author: krzysztof
*/
#include "Render3D.h"
#include "interface.h"
#include "primitives.h"
inline CVector3 Cross(CVector3 A , CVector3 B){
CVector3 C;
C.x = A.y*B.z - A.z*B.y;
C.y = -A.x*B.z + A.z*B.x;
C.z = A.x*B.y - A.y*B.x;
return (C);
}
// The box used by the generalized box fold is specified by a set
// of normalized unit vectors. Each vector is a normal that defines a plane.
// Xp dot N = foldingvalue.
#define sqrt_i3 .57735
CVector3 Nv_tet[] =
{
{ sqrt_i3, sqrt_i3,-sqrt_i3},
{ sqrt_i3,-sqrt_i3, sqrt_i3},
{-sqrt_i3, sqrt_i3, sqrt_i3},
{-sqrt_i3,-sqrt_i3,-sqrt_i3}
};
int sides_tet = 4;
CVector3 Nv_cube[] =
{
{ 1, 0, 0},
{-1, 0, 0},
{ 0, 1, 0},
{ 0,-1, 0},
{ 0, 0, 1},
{ 0, 0,-1}
};
int sides_cube = 6;
CVector3 Nv_oct[] =
{
{ sqrt_i3, sqrt_i3,-sqrt_i3},
{ sqrt_i3,-sqrt_i3, sqrt_i3},
{-sqrt_i3, sqrt_i3, sqrt_i3},
{-sqrt_i3,-sqrt_i3,-sqrt_i3},
{ sqrt_i3, sqrt_i3, sqrt_i3},
{-sqrt_i3,-sqrt_i3, sqrt_i3},
{-sqrt_i3, sqrt_i3,-sqrt_i3},
{ sqrt_i3,-sqrt_i3,-sqrt_i3}
};
int sides_oct = 8;
CVector3 Nv_oct_cube[] =
{
{ sqrt_i3, sqrt_i3,-sqrt_i3},
{ sqrt_i3,-sqrt_i3, sqrt_i3},
{-sqrt_i3, sqrt_i3, sqrt_i3},
{-sqrt_i3,-sqrt_i3,-sqrt_i3},
{ sqrt_i3, sqrt_i3, sqrt_i3},
{-sqrt_i3,-sqrt_i3, sqrt_i3},
{-sqrt_i3, sqrt_i3,-sqrt_i3},
{ sqrt_i3,-sqrt_i3,-sqrt_i3},
{ 1, 0, 0},
{-1, 0, 0},
{ 0, 1, 0},
{ 0,-1, 0},
{ 0, 0, 1},
{ 0, 0,-1}
};
int sides_oct_cube = 14;
#define aa ((1.0+sqrt(5.0))/2.0)
#define bb (1.0/sqrt(aa*aa+1))
CVector3 Nv_dodeca[] =
{
{ 0, bb, aa*bb},
{ 0, bb,-aa*bb},
{ 0,-bb, aa*bb},
{ 0,-bb,-aa*bb},
{ bb, aa*bb, 0},
{ bb,-aa*bb, 0},
{-bb, aa*bb, 0},
{-bb,-aa*bb, 0},
{ aa*bb, 0, bb},
{-aa*bb, 0, bb},
{ aa*bb, 0,-bb},
{-aa*bb, 0,-bb}
};
int sides_dodeca = 12;
/**
* Compute the fractal at the point, in one of the various modes
*
* Mode: normal: Returns distance
* fake_ao: Returns minimum radius
* colouring: Returns colour index
* delta_DE1, delta_DE2: Returns radius
*/
template<int Mode>
double Compute(CVector3 z, const sFractal &par, int *iter_count)
{
int L;
double distance = 0;
double w = 0;
double constantw = 0;
CVector3 dz(1.0, 0.0, 0.0);
CVector3 one(1.0, 0.0, 0.0);
double r_dz = 1;
double ph_dz = 0;
double th_dz = 0;
double p = par.doubles.power; //mandelbulb power
CVector3 constant;
double fixedRadius = par.mandelbox.doubles.foldingSphericalFixed;
double fR2 = fixedRadius * fixedRadius;
double minRadius = par.mandelbox.doubles.foldingSphericalMin;
double mR2 = minRadius * minRadius;
double tglad_factor1 = fR2 / mR2;
double tgladDE = 1.0;
double scale = par.mandelbox.doubles.scale;
enumFractalFormula actualFormula = par.formula;
if (actualFormula == kaleidoscopic || actualFormula == menger_sponge)
tgladDE = 1.0;
double tgladColor = 1.0;
if (par.juliaMode)
{
constant = par.doubles.julia;
}
else
{
constant = z * par.doubles.constantFactor;
}
bool hybridEnabled = false;
if (actualFormula == hybrid) hybridEnabled = true;
double r = z.Length();
double min = 1e200;
for (L = 0; L < par.N; L++)
{
if (hybridEnabled)
{
int tempL = L;
if(tempL > (int)par.formulaSequence.size()-1) tempL = (int)par.formulaSequence.size()-1;
actualFormula = par.formulaSequence[tempL];
p = par.hybridPowerSequence[tempL];
scale = p;
}
if (par.IFS.foldingMode)
{
if (par.IFS.absX) z.x = fabs(z.x);
if (par.IFS.absY) z.y = fabs(z.y);
if (par.IFS.absZ) z.z = fabs(z.z);
for (int i = 0; i < IFS_VECTOR_COUNT; i++)
{
if (par.IFS.enabled[i])
{
z = par.IFS.rot[i].RotateVector(z);
double length = z.Dot(par.IFS.doubles.direction[i]);
if (length < par.IFS.doubles.distance[i])
{
z -= par.IFS.doubles.direction[i] * 2.0 * (length - par.IFS.doubles.distance[i]);
}
}
}
z = par.IFS.mainRot.RotateVector(z - par.IFS.doubles.offset) + par.IFS.doubles.offset;
z *= par.IFS.doubles.scale;
z -= par.IFS.doubles.offset * (par.IFS.doubles.scale - 1.0);
r = z.Length();
}
if (par.tgladFoldingMode)
{
if (z.x > par.doubles.foldingLimit)
{
z.x = par.doubles.foldingValue - z.x;
tgladColor *= 0.9;
}
else if (z.x < -par.doubles.foldingLimit)
{
z.x = -par.doubles.foldingValue - z.x;
tgladColor *= 0.9;
}
if (z.y > par.doubles.foldingLimit)
{
z.y = par.doubles.foldingValue - z.y;
tgladColor *= 0.9;
}
else if (z.y < -par.doubles.foldingLimit)
{
z.y = -par.doubles.foldingValue - z.y;
tgladColor *= 0.9;
}
if (z.z > par.doubles.foldingLimit)
{
z.z = par.doubles.foldingValue - z.z;
tgladColor *= 0.9;
}
else if (z.z < -par.doubles.foldingLimit)
{
z.z = -par.doubles.foldingValue - z.z;
tgladColor *= 0.9;
}
r = z.Length();
}
if (par.sphericalFoldingMode)
{
double fR2_2 = par.doubles.foldingSphericalFixed * par.doubles.foldingSphericalFixed;
double mR2_2 = par.doubles.foldingSphericalMin * par.doubles.foldingSphericalMin;
double r2_2 = r * r;
double tglad_factor1_2 = fR2_2 / mR2_2;
if (r2_2 < mR2_2)
{
z = z * tglad_factor1_2;
tgladDE *= tglad_factor1_2;
tgladColor += 1;
}
else if (r2_2 < fR2_2)
{
double tglad_factor2_2 = fR2_2 / r2_2;
z = z * tglad_factor2_2;
tgladDE *= tglad_factor2_2;
tgladColor += 10;
}
r = z.Length();
}
switch (actualFormula)
{
case trig_DE:
{
double r1 = pow(r, p - 1);
double r2 = r1 * r;
double th = z.GetAlfa();
double ph = -z.GetBeta();
if (Mode == 0)
{
double p_r1_rdz = p * r1 * r_dz;
double ph_phdz = (p - 1.0) * ph + ph_dz;
double th_thdz = (p - 1.0) * th + th_dz;
CVector3 rot(th_thdz, ph_phdz);
dz = rot * p_r1_rdz + one;
r_dz = dz.Length();
th_dz = dz.GetAlfa();
ph_dz = -dz.GetBeta();
}
CVector3 rot(p * th, p * ph);
z = rot * r2 + constant;
r = z.Length();
break;
}
case trig_optim:
{
//optimisation based on: http://www.fractalforums.com/mandelbulb-implementation/realtime-renderingoptimisations/
double th0 = asin(z.z / r);
double ph0 = atan2(z.y, z.x);
double rp = pow(r, p - 1.0);
double th = th0 * p;
double ph = ph0 * p;
double cth = cos(th);
r_dz = rp * r_dz * p + 1.0;
rp *= r;
z = CVector3(cth * cos(ph), cth * sin(ph), sin(th)) * rp + constant;
r = z.Length();
break;
}
case mandelbulb2:
{
double temp, tempR;
tempR = sqrt(z.x * z.x + z.y * z.y);
z *= (1.0 / tempR);
temp = z.x * z.x - z.y * z.y;
z.y = 2.0 * z.x * z.y;
z.x = temp;
z *= tempR;
tempR = sqrt(z.y * z.y + z.z * z.z);
z *= (1.0 / tempR);
temp = z.y * z.y - z.z * z.z;
z.z = 2.0 * z.y * z.z;
z.y = temp;
z *= tempR;
tempR = sqrt(z.x * z.x + z.z * z.z);
z *= (1.0 / tempR);
temp = z.x * z.x - z.z * z.z;
z.z = 2.0 * z.x * z.z;
z.x = temp;
z *= tempR;
z = z * r;
z += constant;
r = z.Length();
break;
}
case mandelbulb3:
{
double temp, tempR;
double sign = 1.0;
double sign2 = 1.0;
if (z.x < 0) sign2 = -1.0;
tempR = sqrt(z.x * z.x + z.y * z.y);
z *= (1.0 / tempR);
temp = z.x * z.x - z.y * z.y;
z.y = 2.0 * z.x * z.y;
z.x = temp;
z *= tempR;
if (z.x < 0) sign = -1.0;
tempR = sqrt(z.x * z.x + z.z * z.z);
z *= (1.0 / tempR);
temp = z.x * z.x - z.z * z.z;
z.z = 2.0 * z.x * z.z * sign2;
z.x = temp * sign;
z *= tempR;
z = z * r;
z += constant;
r = z.Length();
break;
}
case mandelbulb4:
{
double rp = pow(r, p - 1);
double angZ = atan2(z.y, z.x);
double angY = atan2(z.z, z.x);
double angX = atan2(z.z, z.y);
CRotationMatrix rotM;
rotM.RotateX(angX * (p - 1));
rotM.RotateY(angY * (p - 1));
rotM.RotateZ(angZ * (p - 1));
z = rotM.RotateVector(z) * rp + constant;
r = z.Length();
break;
}
case xenodreambuie:
{
double rp = pow(r, p);
double th = atan2(z.y, z.x);
double ph = acos(z.z / r);
if (ph > 0.5 * M_PI)
{
ph = M_PI - ph;
}
else if (ph < -0.5 * M_PI)
{
ph = -M_PI - ph;
}
z.x = rp * cos(th * p) * sin(ph * p);
z.y = rp * sin(th * p) * sin(ph * p);
z.z = rp * cos(ph * p);
z = z + constant;
r = z.Length();
break;
}
case fast_trig:
{
double x2 = z.x * z.x;
double y2 = z.y * z.y;
double z2 = z.z * z.z;
double temp = 1.0 - z2 / (x2 + y2);
double newx = (x2 - y2) * temp;
double newy = 2.0 * z.x * z.y * temp;
double newz = -2.0 * z.z * sqrt(x2 + y2);
z.x = newx + constant.x;
z.y = newy + constant.y;
z.z = newz + constant.z;
r = z.Length();
break;
}
case minus_fast_trig:
{
double x2 = z.x * z.x;
double y2 = z.y * z.y;
double z2 = z.z * z.z;
double temp = 1.0 - z2 / (x2 + y2);
double newx = (x2 - y2) * temp;
double newy = 2.0 * z.x * z.y * temp;
double newz = 2.0 * z.z * sqrt(x2 + y2);
z.x = newx + constant.x;
z.y = newy + constant.y;
z.z = newz + constant.z;
r = z.Length();
break;
}
case hypercomplex:
{
CVector3 newz(z.x * z.x - z.y * z.y - z.z * z.z - w * w, 2.0 * z.x * z.y - 2.0 * w * z.z, 2.0 * z.x * z.z - 2.0 * z.y * w);
double neww = 2.0 * z.x * w - 2.0 * z.y * z.z;
z = newz + constant;
w = neww;
r = sqrt(z.x * z.x + z.y * z.y + z.z * z.z + w * w);
break;
}
case quaternion:
{
CVector3 newz(z.x * z.x - z.y * z.y - z.z * z.z - w * w, 2.0 * z.x * z.y, 2.0 * z.x * z.z);
double neww = 2.0 * z.x * w;
z = newz + constant;
w = neww;
r = sqrt(z.x * z.x + z.y * z.y + z.z * z.z + w * w);
break;
}
case menger_sponge:
{
double temp;
z.x = fabs(z.x);
z.y = fabs(z.y);
z.z = fabs(z.z);
if (z.x - z.y < 0)
{
temp = z.y;
z.y = z.x;
z.x = temp;
}
if (z.x - z.z < 0)
{
temp = z.z;
z.z = z.x;
z.x = temp;
}
if (z.y - z.z < 0)
{
temp = z.z;
z.z = z.y;
z.y = temp;
}
if (Mode == colouring)
{
double length2 = z.Length();
if (length2 < min) min = length2;
}
z *= 3.0;
z.x -= 2.0;
z.y -= 2.0;
if (z.z > 1.0) z.z -= 2.0;
r = z.Length();
tgladDE *= 3.0;
break;
}
case tglad:
{
if (par.mandelbox.rotationsEnabled)
{
bool lockout = false;
z = par.mandelbox.rot[0][0].RotateVector(z);
if (z.x > par.mandelbox.doubles.foldingLimit)
{
z.x = par.mandelbox.doubles.foldingValue - z.x;
tgladColor += par.mandelbox.doubles.colorFactorX;
lockout = true;
}
z = par.mandelbox.rotinv[0][0].RotateVector(z);
z = par.mandelbox.rot[1][0].RotateVector(z);
if (!lockout && z.x < -par.mandelbox.doubles.foldingLimit)
{
z.x = -par.mandelbox.doubles.foldingValue - z.x;
tgladColor += par.mandelbox.doubles.colorFactorX;
}
z = par.mandelbox.rotinv[1][0].RotateVector(z);
lockout = false;
z = par.mandelbox.rot[0][1].RotateVector(z);
if (z.y > par.mandelbox.doubles.foldingLimit)
{
z.y = par.mandelbox.doubles.foldingValue - z.y;
tgladColor += par.mandelbox.doubles.colorFactorY;
lockout = true;
}
z = par.mandelbox.rotinv[0][1].RotateVector(z);
z = par.mandelbox.rot[1][1].RotateVector(z);
if (!lockout && z.y < -par.mandelbox.doubles.foldingLimit)
{
z.y = -par.mandelbox.doubles.foldingValue - z.y;
tgladColor += par.mandelbox.doubles.colorFactorY;
}
z = par.mandelbox.rotinv[1][1].RotateVector(z);
lockout = false;
z = par.mandelbox.rot[0][2].RotateVector(z);
if (z.z > par.mandelbox.doubles.foldingLimit)
{
z.z = par.mandelbox.doubles.foldingValue - z.z;
tgladColor += par.mandelbox.doubles.colorFactorZ;
lockout = true;
}
z = par.mandelbox.rotinv[0][2].RotateVector(z);
z = par.mandelbox.rot[1][2].RotateVector(z);
if (!lockout && z.z < -par.mandelbox.doubles.foldingLimit)
{
z.z = -par.mandelbox.doubles.foldingValue - z.z;
tgladColor += par.mandelbox.doubles.colorFactorZ;
}
z = par.mandelbox.rotinv[1][2].RotateVector(z);
}
else
{
int i;
CVector3 *Nv;
int sides;
// HACK HACK HACK. I high jacked the foldingLimit parameter to chose the poly type.
// Do not try this at home.
Nv = Nv_tet;
sides = sides_tet;
if(int(par.mandelbox.doubles.foldingLimit) == 2){
Nv = Nv_cube;
sides = sides_cube;
}
if(int(par.mandelbox.doubles.foldingLimit) == 3){
Nv = Nv_oct;
sides = sides_oct;
}
if(int(par.mandelbox.doubles.foldingLimit) == 4){
Nv = Nv_dodeca;
sides = sides_dodeca;
}
if(int(par.mandelbox.doubles.foldingLimit) == 5){
Nv = Nv_oct_cube;
sides = sides_oct_cube;
}
int sort[3];
int tmp_sort;
double Z_Dot_Nv[3];
double tmp_Z_Dot_Nv;
// Find the three closest normal vectors to z as defined by max dot product.
sort[0] = 0;
Z_Dot_Nv[0] = z.Dot(Nv[0]);
sort[1]= 1;
Z_Dot_Nv[1] = z.Dot(Nv[1]);
sort[2]= 2;
Z_Dot_Nv[2] = z.Dot(Nv[2]);
if(Z_Dot_Nv[1]>Z_Dot_Nv[0]){
tmp_sort = sort[0];
tmp_Z_Dot_Nv = Z_Dot_Nv[0];
sort[0] = sort[1];
Z_Dot_Nv[0]=Z_Dot_Nv[1];
sort[1]=tmp_sort;
Z_Dot_Nv[1]= tmp_Z_Dot_Nv;
}
if(Z_Dot_Nv[2]>Z_Dot_Nv[1]){
tmp_sort = sort[1];
tmp_Z_Dot_Nv = Z_Dot_Nv[1];
sort[1] = sort[2];
Z_Dot_Nv[1]=Z_Dot_Nv[2];
sort[2]=tmp_sort;
Z_Dot_Nv[2]=tmp_Z_Dot_Nv;
}
if(Z_Dot_Nv[1]>Z_Dot_Nv[0]){
tmp_sort = sort[0];
tmp_Z_Dot_Nv = Z_Dot_Nv[0];
sort[0] = sort[1];
Z_Dot_Nv[0]=Z_Dot_Nv[1];
sort[1]=tmp_sort;
Z_Dot_Nv[1]= tmp_Z_Dot_Nv;
}
for(i=3;i<sides;i++){
tmp_Z_Dot_Nv = z.Dot(Nv[i]);
tmp_sort = i;
if(tmp_Z_Dot_Nv>Z_Dot_Nv[2]){
sort[2] = tmp_sort;
Z_Dot_Nv[2]= tmp_Z_Dot_Nv;
if(tmp_Z_Dot_Nv>Z_Dot_Nv[1]){
sort[2] = sort[1];
Z_Dot_Nv[2]=Z_Dot_Nv[1];
sort[1]=tmp_sort;
Z_Dot_Nv[1]= tmp_Z_Dot_Nv;
if(tmp_Z_Dot_Nv>Z_Dot_Nv[0]){
sort[1] = sort[0];
Z_Dot_Nv[1]=Z_Dot_Nv[0];
sort[0]=tmp_sort;
Z_Dot_Nv[0]= tmp_Z_Dot_Nv;
}
}
}
}
CVector3 Nv0=Nv[sort[0]];
CVector3 Nv1=Nv[sort[1]];
CVector3 Nv2=Nv[sort[2]];
CVector3 new_z;
double new_z_sqr;
CVector3 Zm;
// Assume z inside the poly and we are wasting our time.
new_z = z;
new_z_sqr = new_z.Dot(new_z);
// Find reflection point to closest plain.
Zm = z - (Nv0 + Nv0)*(z.Dot(Nv0) - par.mandelbox.doubles.foldingValue);
if(new_z_sqr > Zm.Dot(Zm)){
new_z = Zm;
new_z_sqr = new_z.Dot(new_z);
tgladColor += par.mandelbox.doubles.colorFactorX;
}
// Find rotation point to closest line.
CVector3 T01, L01;
L01 = Cross(Nv0,Nv1);
L01 = L01 * (1.0/L01.Length());
T01 = (Nv0+Nv1) * ( (par.mandelbox.doubles.foldingValue)/(1 + Nv0.Dot(Nv1)));
CVector3 Zr;
Zr = (T01 + L01 * z.Dot(L01)) * 2 - z;
if(new_z_sqr > Zr.Dot(Zr)){
new_z = Zr;
new_z_sqr = new_z.Dot(new_z);
tgladColor += par.mandelbox.doubles.colorFactorY;
}
// Find inversion point to closest vert.
CVector3 Zi;
double a;
a = ((par.mandelbox.doubles.foldingValue)-T01.Dot(Nv2))/(L01.Dot(Nv2));
Zi = (L01*a + T01) * 2 - z;
if(new_z_sqr > Zi.Dot(Zi)){
new_z = Zi;
new_z_sqr = new_z.Dot(new_z);
tgladColor += par.mandelbox.doubles.colorFactorZ;
}
z = new_z;
}
r = z.Length();
double r2 = r * r;
z += par.mandelbox.doubles.offset;
if (r2 < mR2)
{
z *= tglad_factor1;
tgladDE *= tglad_factor1;
tgladColor += par.mandelbox.doubles.colorFactorSp1;
}
else if (r2 < fR2)
{
double tglad_factor2 = fR2 / r2;
z *= tglad_factor2;
tgladDE *= tglad_factor2;
tgladColor += par.mandelbox.doubles.colorFactorSp2;
}
z -= par.mandelbox.doubles.offset;
z = par.mandelbox.mainRot.RotateVector(z);
z = z * scale + constant;
tgladDE = tgladDE*fabs(scale)+1.0;
r = z.Length();
break;
}
case smoothMandelbox:
{
double sm = par.mandelbox.doubles.sharpness;
double zk1 = SmoothConditionAGreaterB(z.x, par.mandelbox.doubles.foldingLimit,sm);
double zk2 = SmoothConditionALessB(z.x, -par.mandelbox.doubles.foldingLimit,sm);
z.x = z.x * (1.0 - zk1) + (par.mandelbox.doubles.foldingValue - z.x) * zk1;
z.x = z.x * (1.0 - zk2) + (-par.mandelbox.doubles.foldingValue - z.x) * zk2;
tgladColor += (zk1 + zk2) * par.mandelbox.doubles.colorFactorX;
double zk3 = SmoothConditionAGreaterB(z.y, par.mandelbox.doubles.foldingLimit,sm);
double zk4 = SmoothConditionALessB(z.y, -par.mandelbox.doubles.foldingLimit,sm);
z.y = z.y * (1.0 - zk3) + (par.mandelbox.doubles.foldingValue - z.y) * zk3;
z.y = z.y * (1.0 - zk4) + (-par.mandelbox.doubles.foldingValue - z.y) * zk4;
tgladColor += (zk3 + zk4) * par.mandelbox.doubles.colorFactorY;
double zk5 = SmoothConditionAGreaterB(z.z, par.mandelbox.doubles.foldingLimit,sm);
double zk6 = SmoothConditionALessB(z.z, -par.mandelbox.doubles.foldingLimit,sm);
z.z = z.z * (1.0 - zk5) + (par.mandelbox.doubles.foldingValue - z.z) * zk5;
z.z = z.z * (1.0 - zk6) + (-par.mandelbox.doubles.foldingValue - z.z) * zk6;
tgladColor += (zk5 + zk6) * par.mandelbox.doubles.colorFactorZ;
r = z.Length();
double r2 = r * r;
double tglad_factor2 = fR2 / r2;
double rk1 = SmoothConditionALessB(r2, mR2, sm);
double rk2 = SmoothConditionALessB(r2, fR2, sm);
double rk21 = (1.0 - rk1) * rk2;
z = z * (1.0 - rk1) + z * (tglad_factor1 * rk1);
z = z * (1.0 - rk21) + z * (tglad_factor2 * rk21);
tgladDE = tgladDE * (1.0 - rk1) + tgladDE * (tglad_factor1 * rk1);
tgladDE = tgladDE * (1.0 - rk21) + tgladDE * (tglad_factor2 * rk21);
tgladColor += rk1 * par.mandelbox.doubles.colorFactorSp1;
tgladColor += rk21 * par.mandelbox.doubles.colorFactorSp2;
z = par.mandelbox.mainRot.RotateVector(z);
z = z * scale + constant;
tgladDE = tgladDE * fabs(scale) + 1.0;
r = z.Length();
break;
}
case foldingIntPow2:
{
if (z.x > par.doubles.FoldingIntPowFoldFactor) z.x = par.doubles.FoldingIntPowFoldFactor * 2.0 - z.x;
else if (z.x < -par.doubles.FoldingIntPowFoldFactor) z.x = -par.doubles.FoldingIntPowFoldFactor * 2.0 - z.x;
if (z.y > par.doubles.FoldingIntPowFoldFactor) z.y = par.doubles.FoldingIntPowFoldFactor * 2.0 - z.y;
else if (z.y < -par.doubles.FoldingIntPowFoldFactor) z.y = -par.doubles.FoldingIntPowFoldFactor * 2.0 - z.y;
if (z.z > par.doubles.FoldingIntPowFoldFactor) z.z = par.doubles.FoldingIntPowFoldFactor * 2.0 - z.z;
else if (z.z < -par.doubles.FoldingIntPowFoldFactor) z.z = -par.doubles.FoldingIntPowFoldFactor * 2.0 - z.z;
r = z.Length();
double fR2_2 = 1.0;
double mR2_2 = 0.25;
double r2_2 = r * r;
double tglad_factor1_2 = fR2_2 / mR2_2;
if (r2_2 < mR2_2)
{
z = z * tglad_factor1_2;
}
else if (r2_2 < fR2_2)
{
double tglad_factor2_2 = fR2_2 / r2_2;
z = z * tglad_factor2_2;
}
z = z * 2.0;
double x2 = z.x * z.x;
double y2 = z.y * z.y;
double z2 = z.z * z.z;
double temp = 1.0 - z2 / (x2 + y2);
double newx = (x2 - y2) * temp;
double newy = 2.0 * z.x * z.y * temp;
double newz = -2.0 * z.z * sqrt(x2 + y2);
z.x = newx + constant.x;
z.y = newy + constant.y;
z.z = newz + constant.z;
z.z *= par.doubles.FoldingIntPowZfactor;
r = z.Length();
break;
}
case kaleidoscopic:
{
if (par.IFS.absX) z.x = fabs(z.x);
if (par.IFS.absY) z.y = fabs(z.y);
if (par.IFS.absZ) z.z = fabs(z.z);
for (int i = 0; i < IFS_VECTOR_COUNT; i++)
{
if (par.IFS.enabled[i])
{
z = par.IFS.rot[i].RotateVector(z);
double length = z.Dot(par.IFS.doubles.direction[i]);
if (length < par.IFS.doubles.distance[i])
{
z -= par.IFS.doubles.direction[i] * (2.0 * (length - par.IFS.doubles.distance[i]) * par.IFS.doubles.intensity[i]);
}
}
}
z = par.IFS.mainRot.RotateVector(z - par.IFS.doubles.offset) + par.IFS.doubles.offset;
if(par.IFS.doubles.edge.x > 0) z.x = par.IFS.doubles.edge.x - fabs(par.IFS.doubles.edge.x - z.x);
if(par.IFS.doubles.edge.y > 0) z.y = par.IFS.doubles.edge.y - fabs(par.IFS.doubles.edge.y - z.y);
if(par.IFS.doubles.edge.z > 0) z.z = par.IFS.doubles.edge.z - fabs(par.IFS.doubles.edge.z - z.z);
if (Mode == colouring)
{
double length2 = z.Length();
if (length2 < min) min = length2;
}
z *= par.IFS.doubles.scale;
if(par.IFS.mengerSpongeMode)
{
z.x -= par.IFS.doubles.offset.x * (par.IFS.doubles.scale - 1.0);
z.y -= par.IFS.doubles.offset.y * (par.IFS.doubles.scale - 1.0);
if (z.z > 0.5 * par.IFS.doubles.offset.z * (par.IFS.doubles.scale - 1.0)) z.z -= par.IFS.doubles.offset.z * (par.IFS.doubles.scale - 1.0);
}
else
{
z -= par.IFS.doubles.offset * (par.IFS.doubles.scale - 1.0);
}
tgladDE *= par.IFS.doubles.scale;
r = z.Length();
break;
}
case mandelboxVaryScale4D:
{
scale = scale + par.mandelbox.doubles.vary4D.scaleVary * (fabs(scale) - 1.0);
CVector3 oldz = z;
z.x = fabs(z.x + par.mandelbox.doubles.vary4D.fold) - fabs(z.x - par.mandelbox.doubles.vary4D.fold) - z.x;
z.y = fabs(z.y + par.mandelbox.doubles.vary4D.fold) - fabs(z.y - par.mandelbox.doubles.vary4D.fold) - z.y;
z.z = fabs(z.z + par.mandelbox.doubles.vary4D.fold) - fabs(z.z - par.mandelbox.doubles.vary4D.fold) - z.z;
w = fabs(w + par.mandelbox.doubles.vary4D.fold) - fabs(w - par.mandelbox.doubles.vary4D.fold) - w;
if(z.x != oldz.x) tgladColor += par.mandelbox.doubles.colorFactorX;
if(z.y != oldz.y) tgladColor += par.mandelbox.doubles.colorFactorY;
if(z.z != oldz.z) tgladColor += par.mandelbox.doubles.colorFactorZ;
double rr = pow(z.x * z.x + z.y * z.y + z.z * z.z + w * w, par.mandelbox.doubles.vary4D.rPower);
double m = scale;
if (rr < par.mandelbox.doubles.vary4D.minR * par.mandelbox.doubles.vary4D.minR)
{
m = scale / (par.mandelbox.doubles.vary4D.minR * par.mandelbox.doubles.vary4D.minR);
tgladColor += par.mandelbox.doubles.colorFactorSp1;
}
else if (rr < 1.0)
{
m = scale / rr;
tgladColor += par.mandelbox.doubles.colorFactorSp2;
}
z = z * m + constant;
w = w * m + par.mandelbox.doubles.vary4D.wadd;
tgladDE = tgladDE * fabs(m) + 1.0;
r = sqrt(z.x * z.x + z.y * z.y + z.z * z.z + w * w);
break;
}
case aexion:
{
if(L == 0)
{
double cx = fabs(constant.x + constant.y + constant.z) + par.doubles.cadd;
double cy = fabs(-constant.x - constant.y + constant.z) + par.doubles.cadd;
double cz = fabs(-constant.x + constant.y - constant.z) + par.doubles.cadd;
double cw = fabs(constant.x - constant.y - constant.z) + par.doubles.cadd;
constant.x = cx;
constant.y = cy;
constant.z = cz;
constantw = cw;
double tempx = fabs(z.x + z.y + z.z) + par.doubles.cadd;
double tempy = fabs(-z.x - z.y + z.z) + par.doubles.cadd;
double tempz = fabs(-z.x + z.y - z.z) + par.doubles.cadd;
double tempw = fabs(z.x - z.y - z.z) + par.doubles.cadd;
z.x = tempx;
z.y = tempy;
z.z = tempz;
w = tempw;
}
double tempx = z.x * z.x - z.y * z.y + 2.0 * w * z.z + constant.x;
double tempy = z.y * z.y - z.x * z.x + 2.0 * w * z.z + constant.y;
double tempz = z.z * z.z - w * w + 2.0 * z.x * z.y + constant.z;
double tempw = w * w - z.z * z.z + 2.0 * z.x * z.y + constantw;
z.x = tempx;
z.y = tempy;
z.z = tempz;
w = tempw;
r = sqrt(z.x * z.x + z.y * z.y + z.z * z.z + w * w);
break;
}
case benesi:
{
double r1 = z.y*z.y + z.z*z.z;
double newx = 0;
if(constant.x < 0 || z.x < sqrt(r1))
{
newx = z.x*z.x - r1;
}
else
{
newx = -z.x*z.x + r1;
}
r1 = - 1.0/sqrt(r1) * 2.0 * fabs(z.x);
double newy = r1 * (z.y*z.y - z.z*z.z);
double newz = r1 * 2.0 * z.y * z.z;
z.x = newx + constant.x;
z.y = newy + constant.y;
z.z = newz + constant.z;
r = z.Length();
break;
}
case bristorbrot:
{
double newx = z.x*z.x - z.y*z.y - z.z*z.z;
double newy = z.y * (2.0 * z.x - z.z);
double newz = z.z * (2.0 * z.x + z.y);
z.x = newx + constant.x;
z.y = newy + constant.y;
z.z = newz + constant.z;
r = z.Length();
break;
}
case invertX:
{
z.x = z.x >= 0.0 ? z.x*z.x/(fabs(z.x) + p) : -z.x*z.x/(fabs(z.x) + p);
r = z.Length();
break;
}
case invertY:
{
z.y = z.y >= 0.0 ? z.y*z.y/(fabs(z.y) + p) : -z.y*z.y/(fabs(z.y) + p);
r = z.Length();
break;
}
case invertZ:
{
z.z = z.z >= 0.0 ? z.z*z.z/(fabs(z.z) + p) : -z.z*z.z/(fabs(z.z) + p);
r = z.Length();
break;
}
case invertR:
{
double rInv = r*r/(r + p);
z.x = z.x / r * rInv;
z.y = z.y / r * rInv;
z.z = z.z / r * rInv;
r = z.Length();
break;
}
case sphericalFold:
{
double rr = r*r;
double pp = p*p;
if (rr < pp)
{
z.x = 1.0 / pp;
z.y = 1.0 / pp;
z.z = 1.0 / pp;
}
else if (rr < pp*4.0)
{
z.x = 1.0 / rr;
z.y = 1.0 / rr;
z.z = 1.0 / rr;
}
r = z.Length();
break;
}
case powXYZ:
{
z.x = z.x >= 0 ? pow(z.x,p) : -pow(-z.x,p);
z.y = z.y >= 0 ? pow(z.y,p) : -pow(-z.y,p);
z.z = z.z >= 0 ? pow(z.z,p) : -pow(-z.z,p);
r = z.Length();
break;
}
case scaleX:
{
z.x = z.x * p;
r = z.Length();
break;
}
case scaleY:
{
z.y = z.y * p;
r = z.Length();
break;
}
case scaleZ:
{
z.z = z.z * p;
r = z.Length();
break;
}
case offsetX:
{
z.x = z.x + p;
r = z.Length();
break;
}
case offsetY:
{
z.y = z.y + p;
r = z.Length();
break;
}
case offsetZ:
{
z.z = z.z + p;
r = z.Length();
break;
}
case angleMultiplyX:
{
double angle = atan2(z.z,z.y)*p;
double tempR = sqrt(z.z*z.z + z.y*z.y);
z.y = tempR * cos(angle);
z.z = tempR * sin(angle);
r = z.Length();
break;
}
case angleMultiplyY:
{
double angle = atan2(z.z,z.x)*p;
double tempR = sqrt(z.z*z.z + z.x*z.x);
z.x = tempR * cos(angle);
z.z = tempR * sin(angle);
r = z.Length();
break;
}
case angleMultiplyZ:
{
double angle = atan2(z.y,z.x)*p;
double tempR = sqrt(z.x*z.x + z.y*z.y);
z.x = tempR * cos(angle);
z.y = tempR * sin(angle);
r = z.Length();
break;
}
case hybrid:
break;
case none:
break;
}
//************************** iteration terminate conditions *****************
if (Mode == deltaDE1)
{
if (r > 1e10)
break;
}
else if (Mode == deltaDE2)
{
if (L == *iter_count)
break;
}
if (actualFormula == menger_sponge || actualFormula == kaleidoscopic)
{
if (r > 1000)
{
distance = (r - 2.0) / tgladDE;
break;
}
}
else if (actualFormula == tglad || actualFormula == smoothMandelbox || actualFormula == mandelboxVaryScale4D)
{
if (r > 1024)
{
distance = r / fabs(tgladDE);
break;
}
}
else
{
if (Mode == normal) //condition for all other trigonometric and hypercomplex fractals
{
if (r > 1e2)
{
distance = 0.5 * r * log(r) / r_dz;
break;
}
}
else if (Mode == fake_AO) //mode 2
{
if (r < min) min = r;
if (r > 1e15)
{
distance = min;
break;
}
}
else if (Mode == colouring) //mode 1
{
distance = z.Length();
if (distance < min) min = distance;
if (distance > 1e15)
{
distance = min;
break;
}
}
}
}
//************ return values *****************
N_counter += L + 1;
Loop_counter++;
if (L < 64)
histogram[L]++;
else
histogram[63]++;
if (iter_count != NULL)
*iter_count = L;
if (Mode == normal)
{
if (L == par.N)
distance = 0;
return distance;
}
if (Mode == deltaDE1 || Mode == deltaDE2)
return r;
if (Mode == fake_AO)
return distance;
if (Mode == colouring)
{
if (par.formula == hybrid)
{
if (min > 100)
min = 100;
if (distance > 20)
distance = 20;
if (tgladColor > 1000)
tgladColor = 1000;
return distance * 5000.0 + tgladColor * 100.0 + min * 1000.0;
}
else if (actualFormula == tglad || actualFormula == smoothMandelbox || actualFormula == mandelboxVaryScale4D)
return tgladColor * 100.0 + z.Length()*par.mandelbox.doubles.colorFactorR;
else if (actualFormula == kaleidoscopic || actualFormula == menger_sponge)
return min * 1000.0;
else
return distance * 5000.0;
}
}
//******************* Calculate distance *******************8
double CalculateDistance(CVector3 point, sFractal ¶ms, bool *max_iter)
{
int L;
double distance;
params.specialColour = 0;
if (params.limits_enabled)
{
bool limit = false;
double distance_a = 0;
double distance_b = 0;
double distance_c = 0;
if (point.x < params.doubles.amin)
{
distance_a = fabs(params.doubles.amin - point.x) + params.doubles.detailSize;
limit = true;
}
if (point.x > params.doubles.amax)
{
distance_a = fabs(params.doubles.amax - point.x) + params.doubles.detailSize;
limit = true;
}
if (point.y < params.doubles.bmin)
{
distance_a = fabs(params.doubles.bmin - point.y) + params.doubles.detailSize;
limit = true;
}
if (point.y > params.doubles.bmax)
{
distance_b = fabs(params.doubles.bmax - point.y) + params.doubles.detailSize;
limit = true;
}
if (point.z < params.doubles.cmin)
{
distance_c = fabs(params.doubles.cmin - point.z) + params.doubles.detailSize;
limit = true;
}
if (point.z > params.doubles.cmax)
{
distance_c = fabs(params.doubles.cmax - point.z) + params.doubles.detailSize;
limit = true;
}
if (limit)
{
if (max_iter != NULL)
*max_iter = false;
distance = dMax(distance_a, distance_b, distance_c);
return distance;
}
}
if (params.analitycDE)
{
distance = Compute<normal>(point, params, &L);
if (max_iter != NULL)
{
if (L == params.N)
*max_iter = true;
else
*max_iter = false;
}
if (L < params.minN && distance < params.doubles.detailSize)
distance = params.doubles.detailSize;
if (params.interiorMode)
{
if (distance < 0.5 * params.doubles.detailSize || L == params.N)
{
distance = params.doubles.detailSize;
if (max_iter != NULL)
*max_iter = false;
}
}
}
else
{
double deltaDE = 1e-10;
double r = Compute<deltaDE1>(point, params, &L);
int retval = L;
point.x += deltaDE;
point.y += 0;
point.z += 0;
double r2 = Compute<deltaDE2>(point, params, &L);
double dr1 = fabs(r2 - r) / deltaDE;
point.x -= deltaDE;
point.y += deltaDE;
point.z += 0;
r2 = Compute<deltaDE2>(point, params, &L);
double dr2 = fabs(r2 - r) / deltaDE;
point.x += 0;
point.y -= deltaDE;
point.z += deltaDE;
r2 = Compute<deltaDE2>(point, params, &L);
double dr3 = fabs(r2 - r) / deltaDE;
double dr = sqrt(dr1 * dr1 + dr2 * dr2 + dr3 * dr3);
if(params.linearDEmode)
{
distance = 0.5 * r / dr;
}
else
{
distance = 0.5 * r * log(r) / dr;
}
if (retval == params.N)
{
if (max_iter != NULL)
*max_iter = true;
distance = 0;
}
else if (max_iter != NULL)
*max_iter = false;
if (L < params.minN && distance < params.doubles.detailSize)
distance = params.doubles.detailSize;
if (params.interiorMode)
{
if (distance < 0.5 * params.doubles.detailSize || retval == 256)
{
distance = params.doubles.detailSize;
if (max_iter != NULL)
*max_iter = false;
}
}
}
//plane
if (params.primitives.planeEnable)
{
double planeDistance = PrimitivePlane(point, params.doubles.primitives.planeCentre, params.doubles.primitives.planeNormal);
if(planeDistance < distance) params.specialColour = 253;
distance = (planeDistance < distance) ? planeDistance : distance;
}
//box
if (params.primitives.boxEnable)
{
double boxDistance = PrimitiveBox(point, params.doubles.primitives.boxCentre, params.doubles.primitives.boxSize);
if(boxDistance < distance) params.specialColour = 252;
distance = (boxDistance < distance) ? boxDistance : distance;
}
//inverted box
if (params.primitives.invertedBoxEnable)
{
double boxDistance = PrimitiveInvertedBox(point, params.doubles.primitives.invertedBoxCentre, params.doubles.primitives.invertedBoxSize);
if(boxDistance < distance) params.specialColour = 251;
distance = (boxDistance < distance) ? boxDistance : distance;
}
//sphere
if (params.primitives.sphereEnable)
{
double sphereDistance = PrimitiveSphere(point, params.doubles.primitives.sphereCentre, params.doubles.primitives.sphereRadius);
if(sphereDistance < distance) params.specialColour = 250;
distance = (sphereDistance < distance) ? sphereDistance : distance;
}
//invertedSphere
if (params.primitives.invertedSphereEnable)
{
double sphereDistance = PrimitiveInvertedSphere(point, params.doubles.primitives.invertedSphereCentre, params.doubles.primitives.invertedSphereRadius);
if(sphereDistance < distance) params.specialColour = 249;
distance = (sphereDistance < distance) ? sphereDistance : distance;
}
//water
if (params.primitives.waterEnable)
{
double waterDistance = PrimitiveWater(point, params.doubles.primitives.waterHeight, params.doubles.primitives.waterAmplitude,
params.doubles.primitives.waterLength, params.doubles.primitives.waterRotation, params.primitives.waterIterations);
if(waterDistance < distance) params.specialColour = 248;
distance = (waterDistance < distance) ? waterDistance : distance;
}
if (distance < 0) distance = 0;
return distance;
}
// force template instantiation
template double Compute<normal>(CVector3, const sFractal&, int*);
template double Compute<colouring>(CVector3, const sFractal&, int*);
template double Compute<fake_AO>(CVector3, const sFractal&, int*);
template double Compute<deltaDE1>(CVector3, const sFractal&, int*);
template double Compute<deltaDE2>(CVector3, const sFractal&, int*);