Avoid Numerical Error with Strict Floating Point Types

Code is much better off when it takes the discreteness of numerical computation seriously, here is a robust and flexible "Real" class as an example (The following is done using C# style code... but it can be done in any other language just as easily).

It should also be noted that the epsilon value chosen in the example below is somewhat arbitrary. The value at a absolute maximum should be 1E-6, but it is acceptable for values all the way to 1E-12 (and possibly smaller). Just remember that the whole point of it is to provide a work around for the numerical error the results from storing values in a binary mantissa, with a binary exponent. I don't doubt that there is/are thes(is)/(es) on this subject with a statistical study of what a good epsilon is, but for the purposes of a non-critical 3D application I would say you can just fudge it.

struct Real{
   
   public const float EPSILON = 0.0000001f;
   private float number;
   
   public float Value{
      get{ return this.number; }
      set{ this.number = value; }
   }

   public Real(float num){
      this.number = num;
   }

   public static bool operator==(Real lhs, Real rhs){
      if (Abs(lhs.number - rhs.number) < EPSILON){
         return true;
      }
      return false;
   }
   
   // Implement "Equals" if Java or C# using == operator
   // ...

   public static Real operator/(Real lhs, Real rhs){
      ASSERT(Abs(rhs.number) >= EPSILON);
      Real result = new Real(lhs.number/rhs.number);
      return result;
   }
   public static Real operator*(Real lhs, Real rhs){
      Real result = new Real(lhs*rhs);
      return result;
   }
   // ... all other operators   

}

Calculating Normals

It happens time and again where I have to read in a certain file format or I'm creating the initial code for a 3D engine and I completely forget how to go about calculating normals given a set of vertices and faces. So I've decided to write down the pseudo code for this to avoid future aggravation.

The following code assumes that a vertex array exists, a face array exists (where faces are simply a set of indices into the vertex array) and that suitable data structures for holding 3D vectors, faces etc. exist.

// Step 1: Initialize all the vertices to (0,0,0)
foreach vertex v{
   v = Vertex3D(0,0,0);
}

// Step 2: Assign normals...
foreach face f1{
   Vector faceNormal = CrossProduct(f1.v1-f1.v0, f1.v2-f1.v0);
   foreach vertex v1 in f1{
      if (v1.normal is not initialized){
         v1.normal = faceNormal;
      }
      else{
         if (condition for not smoothing){
               /* duplicate v1 */
               /* set the index in f1 to the index 
                   of the duplicated vertex*/
               /* set the normal of the duplicated vertex to 
                   be the normal of f1*/
         }
         else{
            v1.normal += faceNormal;
         }
      }
   }
}

// Step 3: Average normals in vertices that were stupidly 
// duplicated (optional - but necessary for .3DS file loading)
foreach vertex v1{
   foreach vertex v2{
      if (v1 == v2){
         continue;
      }
      if (v1.position == v2.position && condition for smoothing 
          e.g., v1.smoothgrp == v2.smoothgrp){
         v1.position = v2.position;
         v1.normal += v2.normal;
         v2.normal += v1.normal;
      }
   }
}

// Step 4: Normalize every normal
foreach vertex v{
   v.normalize();
}

Explanation:
In step 1 we make sure all normals are initialized to a null vector so that we can add onto them without any consequences (as we do in the next part of the algorithm). In step 2 we iterate through each face and calculate its face normal (this would be the normal of all its vertices in the case of a purely faceted mesh). With this normal we can now pick and choose what vertices will average that with other face normals and which will actually be faceted/unsmoothed.

I have placed the phrase "condition for not smoothing", typically this can be derived from smoothing groups (e.g., in the case of .3DS file loading we could check to see if the smoothing group associated with v1 isn't the same as f1's) or it could be derived from the angle between v1's normal and the current faceNormal (e.g., if the angle between the two normals is greater than 75 degrees then we don't average the normal). In the case where the smoothing group isn't the same or the angle is greater then we have to make another vertex and ensure that f1 uses that vertex instead (this is because the normal must be different to allow for a faceted/unsmoothed look).

Step 3 is a bit of a caveat, in certain cases we have file formats (*cough .3DS), that decided to duplicate certain vertices in the same smoothing group (or within a degree/radian range that should be considered smooth). As a result this will give the appearence of seams on certain rounded objects - most notably spheres. The seam is a result of the point where multiple normals are present at a single vertex; these normals are pointing in the directions of their faces, which leads to multiple normals at a single location in space which should be smooth. The effect of this is faceted faces around the affected vertices. By averaging out the normals of these vertices we gain a normal which is appropriate and eliminates the hideous seams in the mesh. Another important point to mention is that when doing the comparsion between the positions of the two vertices (i.e., v1.position == v2.position, we should be using an "EPSILON" value) - numerical compuation is a slippery slope and we need to ensure that we take the floating point accuracy of the two positions very seriously (floats can very often not be completely equal even when they are supposed to be). Thus the "equals" (==) computation should be carried out as follows:

bool operator==(vector v1, vector v2){
   if (v1 - EPSILON  <= v2 <= v1 + EPSILON){
      return true;
   }
   return false;
}

Where an operation like v1 + EPSILON is interpretted:

(v1.x + EPSILON, v1.y + EPSILON, v1.z + EPSILON)

(Note that there are better ways to do this... epsilon should be a well established concept in your type sturctures/classes/libraries). With that said I would like to point out that we are very lucky that this is a loading procedure and not a real-time, per frame procedure... because O(n^2) algorithms suck.

The final step (step 4) in the algorithm is to simply normalize the normals for every vertex. In doing this we average the normals for any vertices where multiple normals were added on and make the normals acceptable by any typical graphics API (e.g., OpenGL, D3D).

Further Points:

• It might be a good idea to simply remove duplicated vertices within the same smoothing group altogether - this is especially true if you plan on deforming surfaces or manipulating the normals in the future (... and in real-time especially).


• The above procedure could be accomplished as follows:
  

  foreach smoothing group (or some other condition) s{
     foreach vertex v1 in s{
        foreach vertex v2 in s{
           if (v1 == v2){  continue; }
           // As before, use EPSILON in the following
           if (v1.position == v2.position){
              // We must now eliminate either v1 or v2 
              // and point (i.e., change the reference on) 
              // all faces that referenced that vertex towards 
              // the one we don't eliminate
              
              // 1. Eliminate v2
              // 2. Foreach face with a previous 
              // reference to v2, reference it to v1
           }
        }
     }
  }



• Don't fall into the floating point error trap, becareful of reading in very small vertex coordinates as well as calculating very small normals (usually < 1E-06 can be a bad thing, especially when you're dividing)