fractower
Iterator
Posts: 173
|
|
« Reply #15 on: September 22, 2011, 04:56:03 PM » |
|
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*);
|