Metatorus and other iso

I've seen a lot of articles about iso surface generation and triangulation, but all of them explain how to make blobs, or how to do marching cubes, or about normals generation. Imho it's trivial nowadays, to write such articles. This looks like: "mom, I know how to make metaballs!".

In this article I want to explain how to make a metatorus, metacube and how to make a 6DOF rotation of meta primitives. This doc isn't for beginners!

Metatorus.

Metatorus formula is simple (free to say more...):

(x^2 + y^2 + z^2 + R^2 - r^2)^2 = 4*(R^2)*(x^2 + y^2)

(I've found it in POV-ray doc's =))

Here R - radius of revolution, and r - disk radius.

Metacube.

This is for "brain dead" ones! ;) In fact this is a dirty hack, just to make "looks good". To make an "original" meta cube - just fill "cube" in meta field with 1 (or another iso-level value).

(x^n + y^n + z^n) = 1

n - curve factor (8-12 for good result). The greater n is - the sharper the cube edge is.

Rotating meta primitives.

It's just as simple as metatorus. =) Normally meta primitives are defined as F(x,y,z), sampled with step (grid size) in cube with dimensions (W,H,D) => to rotate a meta primitive, just give rotated coordinates to its equation!

Instead of:

for(z=0;z<D;z++) for(y=0;y<H;y++) for(x=0;x<W;x++) field[x + y*W + z*(W*H)] = isoSurface(x - W/2,y - H/2,z - D/2);

write:

point_t t; for(z=0;z<D;z++) for(y=0;y<H;y++) for(x=0;x<W;x++){ t.x = x - W/2; //or something.. t.y = y - H/2; t.z = z - D/2; Rotate(&t,&Matrix); field[x + y*W + z*(W*H)] = isoSurface(t.x,t.y,t.z); }

An experienced eye says: "Cool, but rotating one point is about 9 muls and 6 adds! And if the dimensions of the field cube are 32x32x32 (most common) you have 294912 addition muls, and 196608 addition adds for one meta primitive!! Isn't it too much?!" Yea, this is tooooo much, but who cares.... on a P4 - 1.9GHz you'll, maybe, get >70fps with 32x32x32 field cube and 8 metatoruses. I didn't check this.. =)

So.. and what to do...? The answer on this question is in cube dimensions and grid size: dimensions of the field cube are fixed on constant values, and so is the grid size! => We can replace the rotation with only one vector addition!

Look: The rotation matrix is:

|Xx Yx Zx 0| |Xy Yy Zy 0| |Xz Yz Zz 0| |0 0 0 1|

Here

{Xx,Xy,Xz} - is X axis of a new coordinate system,

{Yx,Yy,Yz} - is Y axis of a new coordinate system,

{Zx,Zy,Zz} - is Z axis of a new coordinate system.

In particular, "axis" is a vector that defines the direction, so to move in space in X axis direction we can just multiply the length we want to move on its direction ({Xx,Yx,Zx}). For moving on a constant length (our grid size) we can just add premultiplied direction to the current value. The same rule is applied for Y and Z.

We know the corner value of the field cube - (-W/2,H/2,D/2), that becomes our "key point", the point from that we "moving".

Now let's rewrite the previous loop:

//some old style C code =) point_t Pz,Pyz,Pxyz; Pz.x = -W/2; Pz.y = -H/2; Pz.z = -D/2; Rotate(&Pz,&Matrix); for(z=0;z<D;z++){ Pyz = Pz; for(y=0;y<H;y++){ Pxyz = Pyz; for(x=0;x<W;x++){ field[x + y*W + z*(W*H)] = isoSurface(Pxyz.x,Pxyz.y,Pxyz.z); Pxyz.x += Matrix[0]; //or something.. =) Pxyz.y += Matrix[4]; Pxyz.z += Matrix[8]; } Pyz.x += Matrix[1]; Pyz.y += Matrix[5]; Pyz.z += Matrix[9]; } Pz.x += Matrix[2]; Pz.y += Matrix[6]; Pz.z += Matrix[10]; }

Pz - point that moves along Z axis, and gives rotated Z to the result, Pyz - point that moves along Y axis and adds rotated Y to the result, Pxyz - point that moves along X axis, and it's(with Pz and Pyz) our result value with rotated X,Y,Z.

And don't forget to calculate vertex normal vectors using meta field data instead of averaging of face normals!!

PS. In this doc "grid size" == 1.

Happy coding.