Math and PMVMatrix: Cleanup and Refine

- Added final qualifier where possible

- Refined API doc

- FloatUtil:
  - Add machine EPSILON
    - fixed value and runtime computed (real machEps)
    - incl. isZero(..), isEqual(..)

  - Add makeRotationAxis(..)
    - Moved from PMVMatrix for reusage

  - Add makeRotationEuler(..)
    - New, not recommended due to Gimbal-Lock

  - Add copyMatrix[Column|Row](..)

  - Add more PI variations and trigo-func float mappings

  - Removed cross and normalize, use VectorUtil!

VectorUtil:
  - Add copyVec*
  - Add equals and isZero w/ and w/o EPSILON
  - Add distance[Square]
  - Add length[Square]

PMVMatrix:
  - Removed 'destroy' method in favor of making most fields 'final'.

AffineTransform:
  - Added AABBox transform
  - Public multiply

3 files changed, 460 insertions, 163 deletions

diff --git a/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
index d2976357d..6177a6b2d 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
@@ -37,23 +37,23 @@ import com.jogamp.common.os.Platform;
  * Basic Float math utility functions.
  * <p>
  * Implementation assumes linear matrix layout in column-major order
- * matching OpenGL's implementation, translation matrix example:
+ * matching OpenGL's implementation, illustration:
  * <pre>
-   Row-Major Order:
-     1 0 0 x
-     0 1 0 y
-     0 0 1 z
-     0 0 0 1
- * </pre>
- * <pre>
-   Column-Major Order:
-     1 0 0 0
-     0 1 0 0
-     0 0 1 0
-     x y z 1
+  Row-Major                    Column-Major (OpenGL):
+
+        | 0  1  2  3  |            | 0  4  8  12 |
+        |             |            |             |
+        | 4  5  6  7  |            | 1  5  9  13 |
+    M = |             |        M = |             |
+        | 8  9  10 11 |            | 2  6  10 14 |
+        |             |            |             |
+        | 12 13 14 15 |            | 3  7  11 15 |
  * </pre>
  * </p>
  * <p>
+ * See <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html">Matrix-FAQ</a>
+ * </p>
+ * <p>
  * Derived from ProjectFloat.java - Created 11-jan-2004
  * </p>
  *
@@ -113,6 +113,92 @@ public class FloatUtil {
   }
 
   /**
+   * Make a rotation matrix from the given axis and angle in radians.
+   * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
+   */
+  public static final void makeRotationAxis(final float angrad, float x, float y, float z, final float[] mat, final int mat_offset, final float[] tmpVec3f) {
+        final float c = cos(angrad);
+        final float ic= 1.0f - c;
+        final float s = sin(angrad);
+
+        tmpVec3f[0]=x; tmpVec3f[1]=y; tmpVec3f[2]=z;
+        VectorUtil.normalize(tmpVec3f);
+        x = tmpVec3f[0]; y = tmpVec3f[1]; z = tmpVec3f[2];
+
+        // Rotation matrix (Row Order):
+        //      xx(1-c)+c  xy(1-c)+zs xz(1-c)-ys 0
+        //      xy(1-c)-zs yy(1-c)+c  yz(1-c)+xs 0
+        //      xz(1-c)+ys yz(1-c)-xs zz(1-c)+c  0
+        //      0          0          0          1
+        final float xy = x*y;
+        final float xz = x*z;
+        final float xs = x*s;
+        final float ys = y*s;
+        final float yz = y*z;
+        final float zs = z*s;
+        mat[0+0*4+mat_offset] = x*x*ic+c;
+        mat[1+0*4+mat_offset] = xy*ic+zs;
+        mat[2+0*4+mat_offset] = xz*ic-ys;
+
+        mat[0+1*4+mat_offset] = xy*ic-zs;
+        mat[1+1*4+mat_offset] = y*y*ic+c;
+        mat[2+1*4+mat_offset] = yz*ic+xs;
+
+        mat[0+2*4+mat_offset] = xz*ic+ys;
+        mat[1+2*4+mat_offset] = yz*ic-xs;
+        mat[2+2*4+mat_offset] = z*z*ic+c;
+  }
+
+  /**
+   * Make a concatenated rotation matrix in column-major order from the given Euler rotation angles in radians.
+   * <p>
+   * The rotations are applied in the given order:
+   * <ul>
+   *  <li>y - heading</li>
+   *  <li>z - attitude</li>
+   *  <li>x - bank</li>
+   * </ul>
+   * </p>
+   * @param bankX the Euler pitch angle in radians. (rotation about the X axis)
+   * @param headingY the Euler yaw angle in radians. (rotation about the Y axis)
+   * @param attitudeZ the Euler roll angle in radians. (rotation about the Z axis)
+   * <p>
+   * Implementation does not use Quaternion and hence is exposed to
+   * <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q34">Gimbal-Lock</a>
+   * </p>
+   * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q36">Matrix-FAQ Q36</a>
+   * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToMatrix/index.htm">euclideanspace.com-eulerToMatrix</a>
+   */
+  public static final void makeRotationEuler(final float bankX, final float headingY, final float attitudeZ, final float[] mat, final int mat_offset) {
+      // Assuming the angles are in radians.
+      final float ch = cos(headingY);
+      final float sh = sin(headingY);
+      final float ca = cos(attitudeZ);
+      final float sa = sin(attitudeZ);
+      final float cb = cos(bankX);
+      final float sb = sin(bankX);
+
+      mat[0+0*4+mat_offset] =  ch*ca;
+      mat[0+1*4+mat_offset] =  sh*sb    - ch*sa*cb;
+      mat[0+2*4+mat_offset] =  ch*sa*sb + sh*cb;
+      mat[1+0*4+mat_offset] =  sa;
+      mat[1+1*4+mat_offset] =  ca*cb;
+      mat[1+2*4+mat_offset] = -ca*sb;
+      mat[2+0*4+mat_offset] = -sh*ca;
+      mat[2+1*4+mat_offset] =  sh*sa*cb + ch*sb;
+      mat[2+2*4+mat_offset] = -sh*sa*sb + ch*cb;
+
+      mat[3+0*4+mat_offset] =  0;
+      mat[3+1*4+mat_offset] =  0;
+      mat[3+2*4+mat_offset] =  0;
+
+      mat[0+3*4+mat_offset] =  0;
+      mat[1+3*4+mat_offset] =  0;
+      mat[2+3*4+mat_offset] =  0;
+      mat[3+3*4+mat_offset] =  1;
+  }
+
+  /**
    * @param a 4x4 matrix in column-major order
    * @param b 4x4 matrix in column-major order
    * @param d result a*b in column-major order
@@ -252,79 +338,6 @@ public class FloatUtil {
   }
 
   /**
-   * Normalize vector
-   *
-   * @param v makes len(v)==1
-   */
-  public static final void normalize(float[] v) {
-    float r = (float) Math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
-
-    if ( r == 0.0 || r == 1.0) {
-      return;
-    }
-
-    r = 1.0f / r;
-
-    v[0] *= r;
-    v[1] *= r;
-    v[2] *= r;
-  }
-
-  /**
-   * Normalize vector
-   *
-   * @param v makes len(v)==1
-   */
-  public static final void normalize(FloatBuffer v) {
-    final int vPos = v.position();
-
-    float r = (float) Math.sqrt(v.get(0+vPos) * v.get(0+vPos) +
-                                v.get(1+vPos) * v.get(1+vPos) +
-                                v.get(2+vPos) * v.get(2+vPos));
-
-    if ( r == 0.0 || r == 1.0) {
-      return;
-    }
-
-    r = 1.0f / r;
-
-    v.put(0+vPos, v.get(0+vPos) * r);
-    v.put(1+vPos, v.get(1+vPos) * r);
-    v.put(2+vPos, v.get(2+vPos) * r);
-  }
-
-
-  /**
-   * Calculate cross-product of 2 vector
-   *
-   * @param v1 3-component vector
-   * @param v2 3-component vector
-   * @param result v1 X v2
-   */
-  public static final void cross(float[] v1, float[] v2, float[] result) {
-    result[0] = v1[1] * v2[2] - v1[2] * v2[1];
-    result[1] = v1[2] * v2[0] - v1[0] * v2[2];
-    result[2] = v1[0] * v2[1] - v1[1] * v2[0];
-  }
-
-  /**
-   * Calculate cross-product of 2 vector
-   *
-   * @param v1 3-component vector
-   * @param v2 3-component vector
-   * @param result v1 X v2
-   */
-  public static final void cross(FloatBuffer v1, FloatBuffer v2, FloatBuffer result) {
-    final int v1Pos = v1.position();
-    final int v2Pos = v2.position();
-    final int rPos  = result.position();
-
-    result.put(0+rPos, v1.get(1+v1Pos) * v2.get(2+v2Pos) - v1.get(2+v1Pos) * v2.get(1+v2Pos));
-    result.put(1+rPos, v1.get(2+v1Pos) * v2.get(0+v2Pos) - v1.get(0+v1Pos) * v2.get(2+v2Pos));
-    result.put(2+rPos, v1.get(0+v1Pos) * v2.get(1+v2Pos) - v1.get(1+v1Pos) * v2.get(0+v2Pos));
-  }
-
-  /**
    * @param m_in 4x4 matrix in column-major order
    * @param m_in_off
    * @param v_in 4-component column-vector
@@ -412,6 +425,46 @@ public class FloatUtil {
   }
 
   /**
+   * Copy the named column of the given column-major matrix to v_out.
+   * <p>
+   * v_out may be 3 or 4 components long, hence the 4th row may not be stored.
+   * </p>
+   * @param m_in input column-major matrix
+   * @param m_in_off offset to input matrix
+   * @param column named column to copy
+   * @param v_out the column-vector storage, at least 3 components long
+   * @param v_out_off offset to storage
+   */
+  public static final void copyMatrixColumn(final float[] m_in, final int m_in_off, final int column, final float[] v_out, final int v_out_off) {
+      v_out[0+v_out_off]=m_in[0+column*4+m_in_off];
+      v_out[1+v_out_off]=m_in[1+column*4+m_in_off];
+      v_out[2+v_out_off]=m_in[2+column*4+m_in_off];
+      if( v_out.length > 3+v_out_off ) {
+          v_out[3+v_out_off]=m_in[3+column*4+m_in_off];
+      }
+  }
+
+  /**
+   * Copy the named row of the given column-major matrix to v_out.
+   * <p>
+   * v_out may be 3 or 4 components long, hence the 4th column may not be stored.
+   * </p>
+   * @param m_in input column-major matrix
+   * @param m_in_off offset to input matrix
+   * @param row named row to copy
+   * @param v_out the row-vector storage, at least 3 components long
+   * @param v_out_off offset to storage
+   */
+  public static final void copyMatrixRow(final float[] m_in, final int m_in_off, final int row, final float[] v_out, final int v_out_off) {
+      v_out[0+v_out_off]=m_in[row+0*4+m_in_off];
+      v_out[1+v_out_off]=m_in[row+1*4+m_in_off];
+      v_out[2+v_out_off]=m_in[row+2*4+m_in_off];
+      if( v_out.length > 3+v_out_off ) {
+          v_out[3+v_out_off]=m_in[row+3*4+m_in_off];
+      }
+  }
+
+  /**
    * @param sb optional passed StringBuilder instance to be used
    * @param f the format string of one floating point, i.e. "%10.5f", see {@link java.util.Formatter}
    * @param a mxn matrix (rows x columns)
@@ -570,20 +623,129 @@ public class FloatUtil {
       return sb;
   }
 
+  @SuppressWarnings("unused")
+  private static void calculateMachineEpsilonFloat() {
+      final long t0;
+      if( DEBUG_EPSILON ) {
+          t0 = Platform.currentTimeMillis();
+      }
+      float machEps = 1.0f;
+      int i=0;
+      do {
+          machEps /= 2.0f;
+          i++;
+      } while (1.0f + (machEps / 2.0f) != 1.0f);
+      machEpsilon = machEps;
+      if( DEBUG_EPSILON ) {
+          final long t1 = Platform.currentTimeMillis();
+          System.err.println("MachineEpsilon: "+machEpsilon+", in "+i+" iterations within "+(t1-t0)+" ms");
+      }
+  }
+  private static volatile boolean machEpsilonAvail = false;
+  private static float machEpsilon = 0f;
+  private static final boolean DEBUG_EPSILON = false;
+
+  /**
+   * Return computed machine Epsilon value.
+   * <p>
+   * The machine Epsilon value is computed once.
+   * </p>
+   * <p>
+   * On a reference machine the result was {@link #EPSILON} in 23 iterations.
+   * </p>
+   * @see #EPSILON
+   */
+  public static float getMachineEpsilon() {
+      if( !machEpsilonAvail ) {
+          synchronized(FloatUtil.class) {
+              if( !machEpsilonAvail ) {
+                  machEpsilonAvail = true;
+                  calculateMachineEpsilonFloat();
+              }
+          }
+      }
+      return machEpsilon;
+  }
+
   public static final float E = 2.7182818284590452354f;
 
+  /** The value PI, i.e. 180 degrees in radians. */
   public static final float PI = 3.14159265358979323846f;
 
-  public static float abs(float a) { return java.lang.Math.abs(a);  }
+  /** The value 2PI, i.e. 360 degrees in radians. */
+  public static final float TWO_PI = 2f * PI;
+
+  /** The value PI/2, i.e. 90 degrees in radians. */
+  public static final float HALF_PI = PI / 2f;
+
+  /** The value PI/4, i.e. 45 degrees in radians. */
+  public static final float QUARTER_PI = PI / 4f;
+
+  /** The value PI^2. */
+  public final static float SQUARED_PI = PI * PI;
+
+  /**
+   * Epsilon for floating point {@value}, as once computed via {@link #getMachineEpsilon()} on an AMD-64 CPU.
+   * <p>
+   * Definition of machine epsilon guarantees that:
+   * <pre>
+   *        1.0f + EPSILON != 1.0f
+   * </pre>
+   * In other words: <i>machEps</i> is the maximum relative error of the chosen rounding procedure.
+   * </p>
+   * <p>
+   * A number can be considered zero if it is in the range (or in the set):
+   * <pre>
+   *    <b>MaybeZeroSet</b> e ]-<i>machEps</i> .. <i>machEps</i>[  <i>(exclusive)</i>
+   * </pre>
+   * While comparing floating point values, <i>machEps</i> allows to clip the relative error:
+   * <pre>
+   *    boolean isZero    = afloat < EPSILON;
+   *    boolean isNotZero = afloat >= EPSILON;
+   *
+   *    boolean isEqual    = abs(bfloat - afloat) < EPSILON;
+   *    boolean isNotEqual = abs(bfloat - afloat) >= EPSILON;
+   * </pre>
+   * </p>
+   * @see #equals(float, float, float)
+   * @see #isZero(float, float)
+   */
+  public static final float EPSILON = 1.1920929E-7f; // Float.MIN_VALUE == 1.4e-45f ; double EPSILON 2.220446049250313E-16d
+
+  /**
+   * Return true if both values are equal, i.e. their absolute delta < <code>epsilon</code>.
+   * @see #EPSILON
+   */
+  public static boolean equals(final float a, final float b, final float epsilon) {
+      return Math.abs(a - b) < epsilon;
+  }
+
+  /**
+   * Return true if value is zero, i.e. it's absolute value < <code>epsilon</code>.
+   * @see #EPSILON
+   */
+  public static boolean isZero(final float a, final float epsilon) {
+      return Math.abs(a) < epsilon;
+  }
+
+  public static float abs(final float a) { return java.lang.Math.abs(a);  }
+
+  public static float pow(final float a, final float b) { return (float) java.lang.Math.pow(a, b);  }
+
+  public static float sin(final float a) { return (float) java.lang.Math.sin(a);  }
+
+  public static float asin(final float a) { return (float) java.lang.Math.asin(a);  }
+
+  public static float cos(final float a) { return (float) java.lang.Math.cos(a);  }
 
-  public static float pow(float a, float b) { return (float) java.lang.Math.pow(a, b);  }
+  public static float acos(final float a) { return (float) java.lang.Math.acos(a);  }
 
-  public static float sin(float a) { return (float) java.lang.Math.sin(a);  }
+  public static float tan(final float a) { return (float) java.lang.Math.tan(a); }
 
-  public static float cos(float a) { return (float) java.lang.Math.cos(a);  }
+  public static float atan(final float a) { return (float) java.lang.Math.atan(a); }
 
-  public static float acos(float a) { return (float) java.lang.Math.acos(a);  }
+  public static float atan2(final float y, final float x) { return (float) java.lang.Math.atan2(y, x); }
 
-  public static float sqrt(float a) { return (float) java.lang.Math.sqrt(a);  }
+  public static float sqrt(final float a) { return (float) java.lang.Math.sqrt(a);  }
 
 }
\ No newline at end of file
diff --git a/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
index 734b7459b..ff16fba79 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
@@ -31,6 +31,9 @@ import java.util.ArrayList;
 
 public class VectorUtil {
 
+    /** Zero vector */
+    public static final float[] ZERO       = new float[] {  0f,  0f,  0f };
+
     public enum Winding {
         CW(-1), CCW(1);
 
@@ -43,41 +46,171 @@ public class VectorUtil {
 
     public static final int COLLINEAR = 0;
 
+    /**
+     * Copies a vector of length 3
+     * @param dst output vector
+     * @param dstOffset offset of dst in array
+     * @param src input vector
+     * @param srcOffset offset of src in array
+     * @return copied output vector for chaining
+     */
+    public static float[] copyVec3(final float[] dst, int dstOffset, final float[] src, int srcOffset)
+    {
+        System.arraycopy(src, srcOffset, dst, dstOffset, 3);
+        return dst;
+    }
+
+    /**
+     * Copies a vector of length 4
+     * @param dst output vector
+     * @param dstOffset offset of dst in array
+     * @param src input vector
+     * @param srcOffset offset of src in array
+     * @return copied output vector for chaining
+     */
+    public static float[] copyVec4(final float[] dst, int dstOffset, final float[] src, int srcOffset)
+    {
+        System.arraycopy(src, srcOffset, dst, dstOffset, 4);
+        return dst;
+    }
+
+    /**
+     * Return true if both vectors are equal, no {@link FloatUtil#EPSILON} is taken into consideration.
+     */
+    public static boolean equals(final float[] vec1, int vec1Offset, final float[] vec2, int vec2Offset) {
+        return vec1[0+vec1Offset] == vec2[0+vec2Offset] &&
+               vec1[1+vec1Offset] == vec2[1+vec2Offset] &&
+               vec1[2+vec1Offset] == vec2[2+vec2Offset];
+    }
+
+    /**
+     * Return true if both vectors are equal, i.e. their absolute delta < <code>epsilon</code>.
+     * @see FloatUtil#EPSILON
+     */
+    public static boolean equals(final float[] vec1, int vec1Offset, final float[] vec2, int vec2Offset, final float epsilon) {
+        return Math.abs(vec1[0+vec1Offset] - vec2[0+vec2Offset]) < epsilon &&
+               Math.abs(vec1[1+vec1Offset] - vec2[1+vec2Offset]) < epsilon &&
+               Math.abs(vec1[2+vec1Offset] - vec2[2+vec2Offset]) < epsilon ;
+    }
+
+    /**
+     * Return true if vector is zero, no {@link FloatUtil#EPSILON} is taken into consideration.
+     */
+    public static boolean isZero(final float[] vec, final int vecOffset) {
+        return 0f == vec[0+vecOffset] && 0f == vec[1+vecOffset] && 0f == vec[2+vecOffset];
+    }
+
+    /**
+     * Return true if vector is zero, i.e. it's absolute components < <code>epsilon</code>.
+     * @see FloatUtil#EPSILON
+     */
+    public static boolean isZero(final float[] vec, final int vecOffset, final float epsilon) {
+        return isZero(vec[0+vecOffset], vec[1+vecOffset], vec[2+vecOffset], epsilon);
+    }
+
+    /**
+     * Return true if all three vector components are zero, i.e. it's their absolute value < <code>epsilon</code>.
+     * @see FloatUtil#EPSILON
+     */
+    public static boolean isZero(final float x, final float y, final float z, final float epsilon) {
+        return Math.abs(x) < epsilon &&
+               Math.abs(y) < epsilon &&
+               Math.abs(z) < epsilon ;
+    }
+
+    /**
+     * Return the squared distance between the given two points described vector v1 and v2.
+     * <p>
+     * When comparing the relative distance between two points it is usually sufficient to compare the squared
+     * distances, thus avoiding an expensive square root operation.
+     * </p>
+     */
+    public static float distanceSquared(final float[] v1, final float[] v2) {
+        final float dx = v1[0] - v2[0];
+        final float dy = v1[1] - v2[1];
+        final float dz = v1[2] - v2[2];
+        return dx * dx + dy * dy + dz * dz;
+    }
+
+    /**
+     * Return the distance between the given two points described vector v1 and v2.
+     */
+    public static float distance(final float[] v1, final float[] v2) {
+        return FloatUtil.sqrt(distanceSquared(v1, v2));
+    }
+
     /** compute the dot product of two points
      * @param vec1 vector 1
      * @param vec2 vector 2
      * @return the dot product as float
      */
-    public static float dot(float[] vec1, float[] vec2)
+    public static float dot(final float[] vec1, final float[] vec2)
     {
         return (vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2]);
     }
 
     /**
+     * Compute the squared length of a vector, a.k.a the squared <i>norm</i>
+     */
+    public static float lengthSquared(final float[] vec) {
+        return vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2];
+    }
+    /**
+     * Compute the length of a vector, a.k.a the <i>norm</i>
+     */
+    public static float length(final float[] vec) {
+        return FloatUtil.sqrt(lengthSquared(vec));
+    }
+
+    /**
      * Normalize a vector
+     * @param result output vector
      * @param vector input vector
-     * @return normalized vector
+     * @return normalized output vector
      */
-    public static float[] normalize(final float[] result, float[] vector)
-    {
-        final float d = FloatUtil.sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
-        if(d> 0.0f)
-        {
-            result[0] = vector[0]/d;
-            result[1] = vector[1]/d;
-            result[2] = vector[2]/d;
+    public static float[] normalize(final float[] result, final float[] vector) {
+        final float lengthSq = lengthSquared(vector);
+        if ( FloatUtil.isZero(lengthSq, FloatUtil.EPSILON) ) {
+            result[0] = 0f;
+            result[1] = 0f;
+            result[2] = 0f;
+        } else {
+            final float invSqr = 1f / FloatUtil.sqrt(lengthSq);
+            result[0] = vector[0] * invSqr;
+            result[1] = vector[1] * invSqr;
+            result[2] = vector[2] * invSqr;
         }
         return result;
     }
 
     /**
+     * Normalize a vector in place
+     * @param result output vector
+     * @param vector input vector
+     * @return normalized output vector
+     */
+    public static float[] normalize(final float[] vector) {
+        final float lengthSq = lengthSquared(vector);
+        if ( FloatUtil.isZero(lengthSq, FloatUtil.EPSILON) ) {
+            vector[0] = 0f;
+            vector[1] = 0f;
+            vector[2] = 0f;
+        } else {
+            final float invSqr = 1f / FloatUtil.sqrt(lengthSq);
+            vector[0] *= invSqr;
+            vector[1] *= invSqr;
+            vector[2] *= invSqr;
+        }
+        return vector;
+    }
+
+    /**
      * Scales a vector by param using given result float[]
      * @param result vector for the result
      * @param vector input vector
      * @param scale single scale constant for all vector components
      */
-    public static float[] scale(float[] result, float[] vector, float scale)
-    {
+    public static float[] scale(final float[] result, final float[] vector, final float scale) {
         result[0] = vector[0] * scale;
         result[1] = vector[1] * scale;
         result[2] = vector[2] * scale;
@@ -90,7 +223,7 @@ public class VectorUtil {
      * @param scale 3 component scale constant for each vector component
      * @return given result vector
      */
-    public static float[] scale(float[] result, float[] vector, float[] scale)
+    public static float[] scale(final float[] result, final float[] vector, final float[] scale)
     {
         result[0] = vector[0] * scale[0];
         result[1] = vector[1] * scale[1];
@@ -99,13 +232,12 @@ public class VectorUtil {
     }
 
     /**
-     * Adds to vectors
+     * Adds two vectors
      * @param v1 vector 1
      * @param v2 vector 2
      * @return v1 + v2
      */
-    public static float[] vectorAdd(float[] result, float[] v1, float[] v2)
-    {
+    public static float[] vectorAdd(final float[] result, final float[] v1, final float[] v2) {
         result[0] = v1[0] + v2[0];
         result[1] = v1[1] + v2[1];
         result[2] = v1[2] + v2[2];
@@ -113,17 +245,29 @@ public class VectorUtil {
     }
 
     /**
+     * Subtracts two vectors
+     * @param v1 vector 1
+     * @param v2 vector 2
+     * @return v1 - v2
+     */
+    public static float[] vectorSub(final float[] result, final float[] v1, final float[] v2) {
+        result[0] = v1[0] - v2[0];
+        result[1] = v1[1] - v2[1];
+        result[2] = v1[2] - v2[2];
+        return result;
+    }
+
+    /**
      * cross product vec1 x vec2
-     * @param vec1 vector 1
-     * @param vec2 vecttor 2
+     * @param v1 vector 1
+     * @param v2 vector 2
      * @return the resulting vector
      */
-    public static float[] cross(final float[] result, float[] vec1, float[] vec2)
+    public static float[] cross(final float[] result, final float[] v1, final float[] v2)
     {
-        result[0] = vec2[2]*vec1[1] - vec2[1]*vec1[2];
-        result[1] = vec2[0]*vec1[2] - vec2[2]*vec1[0];
-        result[2] = vec2[1]*vec1[0] - vec2[0]*vec1[1];
-
+        result[0] = v1[1] * v2[2] - v1[2] * v2[1];
+        result[1] = v1[2] * v2[0] - v1[0] * v2[2];
+        result[2] = v1[0] * v2[1] - v1[1] * v2[0];
         return result;
     }
 
@@ -132,7 +276,7 @@ public class VectorUtil {
      * @param vec vector(x,y,z)
      * @return result
      */
-    public static float[] colMatrixVectorMult(final float[] result, float[] colMatrix, float[] vec)
+    public static float[] colMatrixVectorMult(final float[] result, final float[] colMatrix, final float[] vec)
     {
         result[0] = vec[0]*colMatrix[0] + vec[1]*colMatrix[4] + vec[2]*colMatrix[8] + colMatrix[12];
         result[1] = vec[0]*colMatrix[1] + vec[1]*colMatrix[5] + vec[2]*colMatrix[9] + colMatrix[13];
@@ -146,7 +290,7 @@ public class VectorUtil {
      * @param vec vector(x,y,z)
      * @return result
      */
-    public static float[] rowMatrixVectorMult(final float[] result, float[] rawMatrix, float[] vec)
+    public static float[] rowMatrixVectorMult(final float[] result, final float[] rawMatrix, final float[] vec)
     {
         result[0] = vec[0]*rawMatrix[0] + vec[1]*rawMatrix[1] + vec[2]*rawMatrix[2] + rawMatrix[3];
         result[1] = vec[0]*rawMatrix[4] + vec[1]*rawMatrix[5] + vec[2]*rawMatrix[6] + rawMatrix[7];
@@ -160,7 +304,7 @@ public class VectorUtil {
      * @param p2 second vale
      * @return midpoint
      */
-    public static float mid(float p1, float p2)
+    public static float mid(final float p1, final float p2)
     {
         return (p1+p2)/2.0f;
     }
@@ -171,7 +315,7 @@ public class VectorUtil {
      * @param p2 second point
      * @return midpoint
      */
-    public static float[] mid(final float[] result, float[] p1, float[] p2)
+    public static float[] mid(final float[] result, final float[] p1, final float[] p2)
     {
         result[0] = (p1[0] + p2[0])*0.5f;
         result[1] = (p1[1] + p2[1])*0.5f;
@@ -180,22 +324,13 @@ public class VectorUtil {
         return result;
     }
 
-    /** Compute the norm of a vector
-     * @param vec vector
-     * @return vorm
-     */
-    public static float norm(float[] vec)
-    {
-        return FloatUtil.sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
-    }
-
     /** Compute distance between 2 points
      * @param p0 a ref point on the line
      * @param vec vector representing the direction of the line
      * @param point the point to compute the relative distance of
      * @return distance float
      */
-    public static float computeLength(float[] p0, float[] point)
+    public static float computeLength(final float[] p0, final float[] point)
     {
         final float w0 = point[0]-p0[0];
         final float w1 = point[1]-p0[1];
@@ -209,7 +344,7 @@ public class VectorUtil {
      * @param v2 vertex 2
      * @return
      */
-    public static boolean checkEquality(float[] v1, float[] v2)
+    public static boolean checkEquality(final float[] v1, final float[] v2)
     {
         return Float.compare(v1[0], v2[0]) == 0 &&
                Float.compare(v1[1], v2[1]) == 0 &&
@@ -221,7 +356,7 @@ public class VectorUtil {
      * @param v2 vertex 2
      * @return
      */
-    public static boolean checkEqualityVec2(float[] v1, float[] v2)
+    public static boolean checkEqualityVec2(final float[] v1, final float[] v2)
     {
         return Float.compare(v1[0], v2[0]) == 0 &&
                Float.compare(v1[1], v2[1]) == 0 ;
@@ -233,7 +368,7 @@ public class VectorUtil {
      * @param c vector 3
      * @return the determinant value
      */
-    public static float computeDeterminant(float[] a, float[] b, float[] c)
+    public static float computeDeterminant(final float[] a, final float[] b, final float[] c)
     {
         return a[0]*b[1]*c[2] + a[1]*b[2]*c[0] + a[2]*b[0]*c[1] - a[0]*b[2]*c[1] - a[1]*b[0]*c[2] - a[2]*b[1]*c[0];
     }
@@ -244,7 +379,7 @@ public class VectorUtil {
      * @param v3 vertex 3
      * @return true if collinear, false otherwise
      */
-    public static boolean checkCollinear(float[] v1, float[] v2, float[] v3)
+    public static boolean checkCollinear(final float[] v1, final float[] v2, final float[] v3)
     {
         return (computeDeterminant(v1, v2, v3) == VectorUtil.COLLINEAR);
     }
@@ -254,7 +389,7 @@ public class VectorUtil {
      * @param v1 vertex 1
      * @param v2 vertex2 2
      */
-    public static void computeVector(float[] vector, float[] v1, float[] v2) {
+    public static void computeVector(final float[] vector, final float[] v1, final float[] v2) {
         vector[0] = v2[0] - v1[0];
         vector[1] = v2[1] - v1[1];
         vector[2] = v2[2] - v1[2];
@@ -268,7 +403,7 @@ public class VectorUtil {
      * @return true if the vertex d is inside the circle defined by the
      * vertices a, b, c. from paper by Guibas and Stolfi (1985).
      */
-    public static boolean inCircle(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+    public static boolean inCircle(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
         final float[] A = a.getCoord();
         final float[] B = b.getCoord();
         final float[] C = c.getCoord();
@@ -286,7 +421,7 @@ public class VectorUtil {
      * @return compute twice the area of the oriented triangle (a,b,c), the area
      * is positive if the triangle is oriented counterclockwise.
      */
-    public static float triArea(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c){
+    public static float triArea(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c){
         final float[] A = a.getCoord();
         final float[] B = b.getCoord();
         final float[] C = c.getCoord();
@@ -300,7 +435,7 @@ public class VectorUtil {
      * @return compute twice the area of the oriented triangle (a,b,c), the area
      * is positive if the triangle is oriented counterclockwise.
      */
-    public static float triArea(float[] A, float[] B, float[] C){
+    public static float triArea(final float[] A, final float[] B, final float[] C){
         return (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1])*(C[0] - A[0]);
     }
 
@@ -312,9 +447,9 @@ public class VectorUtil {
      * @param p the vertex in question
      * @return true if p is in triangle (a, b, c), false otherwise.
      */
-    public static boolean vertexInTriangle(float[] a, float[]  b, float[]  c,
-                                           float[] p,
-                                           float[] ac, float[] ab, float[] ap){
+    public static boolean vertexInTriangle(final float[] a, final float[]  b, final float[]  c,
+                                           final float[] p,
+                                           final float[] ac, final float[] ab, final float[] ap){
         // Compute vectors
         computeVector(ac, a, c); //v0
         computeVector(ab, a, b); //v1
@@ -350,9 +485,9 @@ public class VectorUtil {
      * @param tmpAP
      * @return true if p1 or p2 or p3 is in triangle (a, b, c), false otherwise.
      */
-    public static boolean vertexInTriangle3(float[] a, float[]  b, float[]  c,
-                                            float[] p1, float[] p2, float[] p3,
-                                            float[] tmpAC, float[] tmpAB, float[] tmpAP){
+    public static boolean vertexInTriangle3(final float[] a, final float[]  b, final float[]  c,
+                                            final float[] p1, final float[] p2, final float[] p3,
+                                            final float[] tmpAC, final float[] tmpAB, final float[] tmpAP){
         // Compute vectors
         computeVector(tmpAC, a, c); //v0
         computeVector(tmpAB, a, b); //v1
@@ -412,7 +547,7 @@ public class VectorUtil {
      * @param c third vertex
      * @return true if the points a,b,c are in a ccw order
      */
-    public static boolean ccw(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c){
+    public static boolean ccw(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c){
         return triArea(a,b,c) > 0;
     }
 
@@ -422,7 +557,7 @@ public class VectorUtil {
      * @param c third vertex
      * @return Winding
      */
-    public static Winding getWinding(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c) {
+    public static Winding getWinding(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c) {
         return triArea(a,b,c) > 0 ? Winding.CCW : Winding.CW ;
     }
 
@@ -430,7 +565,7 @@ public class VectorUtil {
      * @param vertices
      * @return positive area if ccw else negative area value
      */
-    public static float area(ArrayList<? extends Vert2fImmutable> vertices) {
+    public static float area(final ArrayList<? extends Vert2fImmutable> vertices) {
         final int n = vertices.size();
         float area = 0.0f;
         for (int p = n - 1, q = 0; q < n; p = q++)
@@ -446,7 +581,7 @@ public class VectorUtil {
      * @param vertices array of Vertices
      * @return CCW or CW {@link Winding}
      */
-    public static Winding getWinding(ArrayList<? extends Vert2fImmutable> vertices) {
+    public static Winding getWinding(final ArrayList<? extends Vert2fImmutable> vertices) {
         return area(vertices) >= 0 ? Winding.CCW : Winding.CW ;
     }
 
@@ -458,7 +593,7 @@ public class VectorUtil {
      * @param d vertex 2 of second segment
      * @return the intersection coordinates if the segments intersect, otherwise returns null
      */
-    public static float[] seg2SegIntersection(final float[] result, Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+    public static float[] seg2SegIntersection(final float[] result, final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
         final float determinant = (a.getX()-b.getX())*(c.getY()-d.getY()) - (a.getY()-b.getY())*(c.getX()-d.getX());
 
         if (determinant == 0)
@@ -487,7 +622,7 @@ public class VectorUtil {
      * @param d vertex 2 of second segment
      * @return true if the segments intersect, otherwise returns false
      */
-    public static boolean testSeg2SegIntersection(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+    public static boolean testSeg2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
         final float[] A = a.getCoord();
         final float[] B = b.getCoord();
         final float[] C = c.getCoord();
@@ -520,7 +655,7 @@ public class VectorUtil {
      * @return the intersection coordinates if the lines intersect, otherwise
      * returns null
      */
-    public static float[] line2lineIntersection(final float[] result, Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+    public static float[] line2lineIntersection(final float[] result, final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
         final float determinant = (a.getX()-b.getX())*(c.getY()-d.getY()) - (a.getY()-b.getY())*(c.getX()-d.getX());
 
         if (determinant == 0)
@@ -545,7 +680,7 @@ public class VectorUtil {
      * @param e vertex 2 of first segment
      * @return true if the segment intersects at least one segment of the triangle, false otherwise
      */
-    public static boolean testTri2SegIntersection(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d, Vert2fImmutable e){
+    public static boolean testTri2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d, final Vert2fImmutable e){
         return testSeg2SegIntersection(a, b, d, e) ||
                testSeg2SegIntersection(b, c, d, e) ||
                testSeg2SegIntersection(a, c, d, e) ;
diff --git a/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
index fb311083f..8b4eba222 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
@@ -186,8 +186,8 @@ public class Frustum {
      * </p>
      */
     public void updateByPMV(float[] pmv, int pmv_off) {
-        // Left:   a = m41 + m11, b = m42 + m12, c = m43 + m13, d = m44 + m14  - [1..4] row-major
-        // Left:   a = m30 + m00, b = m31 + m01, c = m32 + m02, d = m33 + m03  - [0..3] row-major
+        // Left:   a = m41 + m11, b = m42 + m12, c = m43 + m13, d = m44 + m14  - [1..4] column-major
+        // Left:   a = m30 + m00, b = m31 + m01, c = m32 + m02, d = m33 + m03  - [0..3] column-major
         {
             final Plane p = planes[LEFT];
             final float[] p_n = p.n;
@@ -197,8 +197,8 @@ public class Frustum {
             p.d    = pmv[ pmv_off + 3 + 3 * 4 ] + pmv[ pmv_off + 0 + 3 * 4 ];
         }
 
-        // Right:  a = m41 - m11, b = m42 - m12, c = m43 - m13, d = m44 - m14  - [1..4] row-major
-        // Right:  a = m30 - m00, b = m31 - m01, c = m32 - m02, d = m33 - m03  - [0..3] row-major
+        // Right:  a = m41 - m11, b = m42 - m12, c = m43 - m13, d = m44 - m14  - [1..4] column-major
+        // Right:  a = m30 - m00, b = m31 - m01, c = m32 - m02, d = m33 - m03  - [0..3] column-major
         {
             final Plane p = planes[RIGHT];
             final float[] p_n = p.n;
@@ -208,8 +208,8 @@ public class Frustum {
             p.d    = pmv[ pmv_off + 3 + 3 * 4 ] - pmv[ pmv_off + 0 + 3 * 4 ];
         }
 
-        // Bottom: a = m41 + m21, b = m42 + m22, c = m43 + m23, d = m44 + m24  - [1..4] row-major
-        // Bottom: a = m30 + m10, b = m31 + m11, c = m32 + m12, d = m33 + m13  - [0..3] row-major
+        // Bottom: a = m41 + m21, b = m42 + m22, c = m43 + m23, d = m44 + m24  - [1..4] column-major
+        // Bottom: a = m30 + m10, b = m31 + m11, c = m32 + m12, d = m33 + m13  - [0..3] column-major
         {
             final Plane p = planes[BOTTOM];
             final float[] p_n = p.n;
@@ -219,8 +219,8 @@ public class Frustum {
             p.d    = pmv[ pmv_off + 3 + 3 * 4 ] + pmv[ pmv_off + 1 + 3 * 4 ];
         }
 
-        // Top:   a = m41 - m21, b = m42 - m22, c = m43 - m23, d = m44 - m24  - [1..4] row-major
-        // Top:   a = m30 - m10, b = m31 - m11, c = m32 - m12, d = m33 - m13  - [0..3] row-major
+        // Top:   a = m41 - m21, b = m42 - m22, c = m43 - m23, d = m44 - m24  - [1..4] column-major
+        // Top:   a = m30 - m10, b = m31 - m11, c = m32 - m12, d = m33 - m13  - [0..3] column-major
         {
             final Plane p = planes[TOP];
             final float[] p_n = p.n;
@@ -230,8 +230,8 @@ public class Frustum {
             p.d    = pmv[ pmv_off + 3 + 3 * 4 ] - pmv[ pmv_off + 1 + 3 * 4 ];
         }
 
-        // Near:  a = m41 + m31, b = m42 + m32, c = m43 + m33, d = m44 + m34  - [1..4] row-major
-        // Near:  a = m30 + m20, b = m31 + m21, c = m32 + m22, d = m33 + m23  - [0..3] row-major
+        // Near:  a = m41 + m31, b = m42 + m32, c = m43 + m33, d = m44 + m34  - [1..4] column-major
+        // Near:  a = m30 + m20, b = m31 + m21, c = m32 + m22, d = m33 + m23  - [0..3] column-major
         {
             final Plane p = planes[NEAR];
             final float[] p_n = p.n;
@@ -241,8 +241,8 @@ public class Frustum {
             p.d    = pmv[ pmv_off + 3 + 3 * 4 ] + pmv[ pmv_off + 2 + 3 * 4 ];
         }
 
-        // Far:   a = m41 - m31, b = m42 - m32, c = m43 - m33, d = m44 - m34  - [1..4] row-major
-        // Far:   a = m30 - m20, b = m31 - m21, c = m32 + m22, d = m33 + m23  - [0..3] row-major
+        // Far:   a = m41 - m31, b = m42 - m32, c = m43 - m33, d = m44 - m34  - [1..4] column-major
+        // Far:   a = m30 - m20, b = m31 - m21, c = m32 + m22, d = m33 + m23  - [0..3] column-major
         {
             final Plane p = planes[FAR];
             final float[] p_n = p.n;