Here's an example with simple symmetric parameters, which could easily extend to 3d.
<Quoted Image Removed>
shifts is 0, bends are (1,0) and (0,1), and flips are (-1,1) and (1,-1). Giving simpler code:
getPixelCount(int &count, Vector2 &point, int depth)
for (i = 0; i<2; i++)
pos = point - bends[i];
pos /= dot(pos, pos);
pos += bends[i];
pos -= 2.0 * flips[i] * pos.dot(flips[i]) / dot(flips[i], flips[i]);
pos *= scale;
pos += C;
if (dot(pos, pos) > 4)
continue;
if (depth == 0)
count++;
else
getPixelCount(count, pos, depth - 1);
In this simplified code you reduce the matrix outcomes to simpler vector sums using normalised vectors derived from the pixel position. However you do not define scale anywhere I can see.
Anyway, taking point to mean a structure holding at least the 2 screen pixel coordinates you subtract the unit vectors and normalise to give you a unit vector. There is no zero test here shown prior to the normalisation.
You then start to manipulate this unit vector, but I do not follow the code at this point as the for instruction is spent as I understand it. Bends[] are rotations and at this point we would double the rotation for a Mandelbrot outcome. However, you have no machinery in this code to facilitate that, so I suggest using (-1,0 ) and (0,-1).
To generalise the rotation here I would use bends[] as
(cos ( arctan( y/x)), sin( arctan(y/x)))
where y and x are derived from point, possibly point.x and point.y. And also instance a matrix rotation procedure. The addition of the vectors here is only a partial rotation and introduces some Shift.
The general doubling could then be achieved by doubling the arctan, or arctan2.
You then do some normalised flips[], which in general are reflections in the origin. 2 such reflections are equivalent to a rotation through 2pi, but in this case reflection in the origin is always equivalent to a rotation by pi. It is reflection in the axes that has to be repeated to obtain a rotation. Your flips are rotations of pi, but again, if you add them you introduce shift into the vector resultants. To avoid this requires a matrix product procedure. Doubling it only increases this shift drift!
Then you scale. Here I would suggest a scale that is derived from point.x*point.y*2.
These suggestions are clearly based on the Mandelbrot formulae, and it is the last one that is problematic with regard to being conformal. It may be that this scale factor is what ultimately ruins the project, but rather than giving in maybe the modulo arithmetic can be used to to constrain it within bounds in a pseudo normalisation.
Thinking this through makes me suspicious about this last constraint being physically relevant in modeling dynamic systems. But what do I know?