Title: Generalized Box Fold.
Post by: fractower on August 25, 2010, 07:09:02 AM
I have taken up Tglad's Mandelbox advanced folding challenge and think I have an algorithm for a generalized conformal box fold. Consider any convex box (cube, tet, potato etc.). Identity map any point inside the box. For any point outside the box find the closest point on the box and perform a point inversion. Because the box is convex the closest point is unique.
Now apply this to a cube. Any localized group of points where the closest point is a plane will be reflected. Any localized group of points which are closest to a line will be rotated by Pi around the line. Any localized group of points which are closest to a vertex will be inverted through the vertex. This performs the same operation as the box fold for a greatly increased computational cost.
Now apply this algorithm to any convex box consisting of planes, edges and vortexes. The result is a kind of conformal convex box inversion which preserves the symmetry of the box. I suspect that any resulting images will be less interesting than those obtained using the kaleidoscope method since symmetries are often boring, but Who CAres WheN WE ARE TAlKInG ABOUT MATH PURITY... OOPS, I tend to over capitalize when I get excited.
I am not sure that applying this to a box with curves will not be conformal, but it may result in interesting pictures.
Title: Re: Generalized Box Fold.
Post by: hobold on August 25, 2010, 03:13:26 PM
I agree, this algorithm really does lead to a contiguous conformal mapping. I am not sure if it is really equivalent to the Mandelbox, though.
There is one kind of "curvy boxes" to which this construction can be easily generalized. As an example, take the set of points that are exactly at some given distance d from a given cube. This gets you a rounded cube: faces are still flat squares, but edges are now cylindrical and corners are spherical surface patches (both of radius d). The original reflection with respect to the the closest point can still be done with this new surface, and still leads to a contiguous conformal mapping.
Title: Re: Generalized Box Fold.
Post by: fractower on August 25, 2010, 08:29:30 PM
The box fold equivalence probably needs a bit more explanation. For the box fold the points outside the box will be reflected 1, 2, or 3 times by orthogonal planes. A single reflection is just a reflection. Two reflections is equivalent to a rotation by Pi around the line of intersection. Three reflections is equivalent to a point inversion through the point of intersection of the three planes. The points which are reflected only once are in a prism projection orthogonal to a cube face which coincidentally are the points whose closest point to the cube is that cube face. A similar wedge and pyramid exist for the 2 and 3 reflection points.
The melted cube will be relatively easy to implement. The result of the 6 comparisons used to determine the box folds can also be used to determine the whether a reflection, rotation or inversion is to be performed. It will also be nice to compare with the original Mandelbox.
Title: Re: Generalized Box Fold.
Post by: Tglad on August 26, 2010, 03:50:28 AM
Nice way to think of the box fold and it sounds like it does indeed generalise to any convex polyhedron fold, and it isn't a sequence of folds so is symmetric. Its quite similar to msltoe's fractals in the way it reflects around nearest points on a polyhedron http://www.fractalforums.com/3d-fractal-generation/sierpinski-like-fractals-using-an-iterative-function/ (http://www.fractalforums.com/3d-fractal-generation/sierpinski-like-fractals-using-an-iterative-function/).
Can't do curvy folds without losing conformality (or whatever the word is), but might make nice, stretchy results none-the-less. Also, it you couldn't use this directly on a spherical fold as the distance to edge isn't maintained in conformal sphere inversions, so that's one reason why it couldn't apply to inflated polyhedrons.
Title: Re: Generalized Box Fold.
Post by: hobold on August 26, 2010, 03:46:18 PM
Dang, you are right. :-( Inversions aren't just reflections on spheres. And there seems to be no way to properly blend a true inversion into such a piecewise map. Still, it might be enlightening to continuously control the distortion of angles and get a feel for the amount and type of whipped cream that results from that particular flaw.
Title: Re: Generalized Box Fold.
Post by: fractower on September 21, 2011, 08:52:47 AM
I finally had time to program the generalized box fold and produce some renders. I hope to get the the melted cube version that Hobold suggested soon.
(http://nocache-nocookies.digitalgott.com/gallery/8/2792_21_09_11_8_47_33.jpeg)
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 21, 2011, 09:04:53 AM
Formulaaaaaaaa! ;D
Title: Re: Generalized Box Fold.
Post by: fractower on September 21, 2011, 03:17:25 PM
It is more of an algorithm than a formula. If z is inside the polyhedron, do nothing. If z is outside the box, find the closest point on the surface of the box and do a point inversion.
If I can figure out how to post code on the thread, I will. I only touch fractal.cpp, so anyone who can build mandelbulber 1.08 can play with it.
Title: Re: Generalized Box Fold.
Post by: cKleinhuis on September 21, 2011, 04:06:29 PM
just use
the
code [ code ] tag
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 21, 2011, 04:18:11 PM
Do you mean it's not an escape-time formula? :'(
Title: Re: Generalized Box Fold.
Post by: taurus on September 21, 2011, 04:24:13 PM
@ fractower
ever thought of cooperating with buddhi? afaik he is in search of people helping him coding new formulas in mandelbulber.
Title: Re: Generalized Box Fold.
Post by: fractower on September 21, 2011, 06:22:55 PM
Do you mean it's not an escape-time formula? It is essentially Tglad's Mandelbox with the cubic box fold replaced with a generalized box fold. Is there a way to add a source file to a thread?
Title: Re: Generalized Box Fold.
Post by: Buddhi on September 21, 2011, 06:24:21 PM
@ fractower
ever thought of cooperating with buddhi? afaik he is in search of people helping him coding new formulas in mandelbulber.
Of course I'm interested to implement some new formulas
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 21, 2011, 06:49:25 PM
Do you mean it's not an escape-time formula? It is essentially Tglad's Mandelbox with the cubic box fold replaced with a generalized box fold. Is there a way to add a source file to a thread? Yes, so paste directly the code please! :dink:
Title: Re: Generalized Box Fold.
Post by: fractower on September 22, 2011, 04:50:51 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. 19a20,111
>
>
> 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;
>
>
456,484c548,626
< if (z.x > par.mandelbox.doubles.foldingLimit)
< {
< z.x = par.mandelbox.doubles.foldingValue - z.x;
< tgladColor += par.mandelbox.doubles.colorFactorX;
< }
< else if (z.x < -par.mandelbox.doubles.foldingLimit)
< {
< z.x = -par.mandelbox.doubles.foldingValue - z.x;
< tgladColor += par.mandelbox.doubles.colorFactorX;
< }
< if (z.y > par.mandelbox.doubles.foldingLimit)
< {
< z.y = par.mandelbox.doubles.foldingValue - z.y;
< tgladColor += par.mandelbox.doubles.colorFactorY;
< }
< else if (z.y < -par.mandelbox.doubles.foldingLimit)
< {
< z.y = -par.mandelbox.doubles.foldingValue - z.y;
< tgladColor += par.mandelbox.doubles.colorFactorY;
< }
< if (z.z > par.mandelbox.doubles.foldingLimit)
< {
< z.z = par.mandelbox.doubles.foldingValue - z.z;
< tgladColor += par.mandelbox.doubles.colorFactorZ;
< }
< else if (z.z < -par.mandelbox.doubles.foldingLimit)
< {
< z.z = -par.mandelbox.doubles.foldingValue - z.z;
< tgladColor += par.mandelbox.doubles.colorFactorZ;
---
> 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;
485a628,674
> }
> }
> }
> 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;
487c676
<
---
>
Title: Re: Generalized Box Fold.
Post by: fractower 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*);
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 22, 2011, 04:58:52 PM
It looks far too complicated to translate into raw assembly code for me :-\ :'( I pass the ball :tongue1:
Title: Re: Generalized Box Fold.
Post by: marius on September 22, 2011, 07:23:39 PM
It looks far too complicated to translate into raw assembly code for me :-\ :'( I pass the ball :tongue1:
With the right incantations about position independence and calling conventions, and discipline regarding globals, you can write the m3b formulae in C, no?
Title: Re: Generalized Box Fold.
Post by: Aexion on September 22, 2011, 08:42:26 PM
None was interested (on implementing it), but still here are the multidimensional box fold: :)
(ok,ok, they where overly complex.. but it was a nice try!
also, there are other shapes, such as donuts and diamonds)
http://www.fractalforums.com/3d-fractal-generation/platonic-dimensions/ (http://www.fractalforums.com/3d-fractal-generation/platonic-dimensions/)
(http://www.rfractals.net/share/Icosahedron.png)
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 22, 2011, 09:15:34 PM
Your formula has a different look from this one ;)
And for the question to me; no! It is too time consuming passing to C code to assembly. Anyway if you think it is simple... I will pass u the addresses and you cqn try yourself ;)
Title: Re: Generalized Box Fold.
Post by: Aexion on September 22, 2011, 11:44:31 PM
Your formula has a different look from this one ;)
Ok.ok.. I know.. I just was mentioning them (my formulas are always a mess :) ).. Anyways, if you have nothing to do, you can just play with this GyroidFold formula, made of nested Gyroids, complete with my best wishes to the user: bool BGyroidfold(double x, double y, double z,int miter){
const double ctx =x;
const double cty =y;
const double ctz =z;
x=0;
y=0;
z=0;
for(unsigned char Counter=0;Counter<miter;Counter++){
if(x>1)x=2-x;
else
if(x<-1)x=-2-x;
if(y>1)y=2-y;
else
if(y<-1)y=-2-y;
if(z>1)z=2-z;
else
if(z<-1)z=-2-z;
const float rw=(cos(x) * sin(y) + cos(y) * sin(z) + cos(z) * sin(x))*3.141592653;//Gyroid
const float vm=fabs(rw);
if(vm<0.5){
x=x*4;
y=y*4;
z=z*4;
} else
if(vm<1){
const float vsq=vm*vm;
x/=vsq;
y/=vsq;
z/=vsq;
}
x=x*2+ctx;
y=y*2+cty;
z=z*2+ctz;
if ((x*x+y*y+z*z)>3600) break;
}
if((x*x+y*y+z*z)<3600)
return true;
else
return false;
}
:devil:
Title: Re: Generalized Box Fold.
Post by: marius on September 23, 2011, 01:05:22 AM
And for the question to me; no! It is too time consuming passing to C code to assembly. Anyway if you think it is simple... I will pass u the addresses and you cqn try yourself ;)
darkbeam, how about something like: // From fractalforums.com/mandelbulb-3d/asking-permission/msg27618/#msg27618
__attribute__((packed)) struct TIteration3Dext {
double J4,Rold,RStopD,x,y,z,w;//: Double; // use with neg indizes before C1
double C1, C2, C3;//: Double; //
double J1, J2, J3;//: Double; //+24 these are the constants for adding. (x0, y0, z0) or (Cx, Cy, Cz). w0 is J4 (-56)
void* PVar;//: Pointer; //+48 the actual formulas adress of constants and vars
float SmoothItD;//: Single; //+52
double Rout;//: Double; //+56 the square of the vector length, can be used as input
int ItResultI;//: Integer; //+64
int maxIt;//: Integer; //+68
float RStop;//: Single; //+72
int nHybrid[5];//: array[0..5] of Integer; //+76 Hybrid counts or single weights in interpolhybrid
void* fHPVar[5];//: array[0..5] of Pointer; //+100 pointer to constants+vars, PVars-8=0.5, dOption values below -8
void* fHybrid[5];//: array[0..5] of ThybridIteration2; //+124 adresses of formulas
int CalcSIT;//: LongBool; //+148
int DoJulia;//: LongBool; //+152
float LNRStop;//: Single; //+156
int DEoption;//: Integer; //+160
float fHln[5];//: array[0..5] of Single; //+164 for SmoothIts
int iRepeatFrom;//: Integer; //+188
double OTrap;//: Double; //+192
double VaryScale;//: Double; //+200 to use in vary by its
int bFirstIt;//: Integer; //+208 used also as iteration count, is set to 0 on it-start
int bTmp;//: Integer; //+212 tmpBuf, free of use.
double Dfree1;//: Double; //+216
double Dfree2;//: Double; //+224
double Deriv1;//: Double; //+232 for 4D first deriv or as full derivs
double Deriv2;//: Double; //+240
double Deriv3;//: Double; //+248
double SMatrix4[4][4];//: TSMatrix4; //+256
};
// fastcall is not quite delphi fastcall.
// first two args are ok, third is in ecx in delphi, on stack here.
void __attribute__((fastcall)) formula(
double* x, // [eax]
double* y, // [edx]
double* arg, // [ebp+8], points to TIteration3Dext.C1
void* dummy // so we end w/ ret 8 as delphi expects
) {
// Compute ptr to proper start of TIteration3Dext struct.
struct TIteration3Dext* cfg = (struct TIteration3Dext*)(arg-7);
// Read / write some fields as demonstration.. not a real formula ;-)
double r = cfg->x + cfg->y - cfg->z * cfg->w;
cfg->Deriv1 = r;
}
compile with 'gcc -c -Os -m32'. The objdump -D of the generated code shows: Disassembly of section .text:
00000000 <formula>:
0: 55 push %ebp
1: 89 e5 mov %esp,%ebp
3: 8b 45 08 mov 0x8(%ebp),%eax
6: 83 e8 38 sub $0x38,%eax
9: dd 40 18 fldl 0x18(%eax)
c: dc 40 20 faddl 0x20(%eax)
f: dd 40 28 fldl 0x28(%eax)
12: dc 48 30 fmull 0x30(%eax)
15: de e9 fsubrp %st,%st(1)
17: dd 98 10 01 00 00 fstpl 0x110(%eax)
1d: 5d pop %ebp
1e: c2 08 00 ret $0x8
That should be fairly compatible with m3b's expectations, eyeballing suggests. Some definitions need polishing and further explanations (perhaps Jesse can chime in, but you probably know the details already). But after that this should make it a whole lot easier to code m3f files than in asm. And you get compiler optimizations for free ;D
Title: Re: Generalized Box Fold.
Post by: DarkBeam on September 23, 2011, 09:04:46 AM
Never really tried that way :o ... It can be a good idea ty!
Luca
Title: Re: Generalized Box Fold.
Post by: fractower on September 23, 2011, 09:41:59 AM
I was able to get the melting option to work.
(http://nocache-nocookies.digitalgott.com/gallery/8/2792_23_09_11_9_36_55.jpeg)
Title: Re: Generalized Box Fold.
Post by: Jesse on September 23, 2011, 09:40:43 PM
That should be fairly compatible with m3b's expectations, eyeballing suggests.
Some definitions need polishing and further explanations (perhaps Jesse can chime in, but you probably know the details already).
But after that this should make it a whole lot easier to code m3f files than in asm. And you get compiler optimizations for free ;D
Thx, i added a C compatible struct in the message you mentioned: www.fractalforums.com/mandelbulb-3d/asking-permission/msg27618/#msg27618 (http://www.fractalforums.com/mandelbulb-3d/asking-permission/msg27618/#msg27618) But you always have to check whether calls or absolute addresses are used by the compiler. For constants and user variables you could also do the trick like in "struct TIteration3Dext* cfg = .." with the PVar pointer. Luca, hijacking is stearing the thread in a not intended direction like m3d intern struct code in a generalized box fold topic thread :dink: (just because you wanted to know this in another thread, you are not guilty or so!)
Title: Re: Generalized Box Fold.
Post by: knighty on September 23, 2011, 11:31:57 PM
Nice sequence fractower. A video?
|