Merge branch 'master' into stash_glyphcache

Conflicts:
	make/scripts/tests.sh
	src/jogl/classes/com/jogamp/graph/curve/OutlineShape.java
	src/jogl/classes/com/jogamp/graph/curve/Region.java
	src/jogl/classes/com/jogamp/graph/curve/opengl/GLRegion.java
	src/jogl/classes/com/jogamp/graph/curve/opengl/RegionRenderer.java
	src/jogl/classes/com/jogamp/graph/curve/opengl/Renderer.java
	src/jogl/classes/com/jogamp/graph/curve/opengl/TextRenderer.java
	src/jogl/classes/com/jogamp/graph/font/Font.java
	src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
	src/jogl/classes/jogamp/graph/curve/text/GlyphShape.java
	src/jogl/classes/jogamp/graph/curve/text/GlyphString.java
	src/jogl/classes/jogamp/graph/font/typecast/TypecastFont.java
	src/jogl/classes/jogamp/graph/font/typecast/TypecastGlyph.java
	src/jogl/classes/jogamp/graph/font/typecast/TypecastRenderer.java

11 files changed, 3233 insertions, 0 deletions

diff --git a/src/jogl/classes/com/jogamp/opengl/math/Binary16.java b/src/jogl/classes/com/jogamp/opengl/math/Binary16.java
new file mode 100644
index 000000000..33add46c2
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Binary16.java
@@ -0,0 +1,569 @@
+/**
+ * Copyright 2013 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+
+package com.jogamp.opengl.math;
+
+/**
+ * <p>
+ * Functions to convert values to/from the <code>binary16</code> format
+ * specified in <code>IEEE 754 2008</code>.
+ * </p>
+ */
+
+public final class Binary16
+{
+  /**
+   * The encoded form of negative infinity <code>-∞</code>.
+   */
+
+  public static final char NEGATIVE_INFINITY;
+
+  /**
+   * The encoded form of positive infinity <code>∞</code>.
+   */
+
+  public static final char POSITIVE_INFINITY;
+
+  /**
+   * The encoded form of positive zero <code>0</code>.
+   */
+
+  public static final char POSITIVE_ZERO;
+
+  /**
+   * The encoded form of negative zero <code>-0</code>.
+   */
+
+  public static final char NEGATIVE_ZERO;
+
+  /**
+   * The <i>bias</i> value used to offset the encoded exponent. A given
+   * exponent <code>e</code> is encoded as <code>{@link #BIAS} + e</code>.
+   */
+
+  public static final int  BIAS;
+
+  static {
+    NEGATIVE_INFINITY = 0xFC00;
+    POSITIVE_INFINITY = 0x7C00;
+    POSITIVE_ZERO = 0x0000;
+    NEGATIVE_ZERO = 0x8000;
+    BIAS = 15;
+  }
+
+  private static final int MASK_SIGN;
+  private static final int MASK_EXPONENT;
+  private static final int MASK_SIGNIFICAND;
+
+  static {
+    MASK_SIGN = 0x8000;
+    MASK_EXPONENT = 0x7C00;
+    MASK_SIGNIFICAND = 0x03FF;
+  }
+
+  /**
+   * One possible not-a-number value.
+   */
+
+  public static char exampleNaN()
+  {
+    final int n =
+      Binary16.packSetExponentUnbiasedUnchecked(16)
+        | Binary16.packSetSignificandUnchecked(1);
+    final char c = (char) n;
+    return c;
+  }
+
+  /**
+   * Return <code>true</code> if the given packed <code>binary16</code> value
+   * is infinite.
+   */
+
+  public static boolean isInfinite(
+    final char k)
+  {
+    if (Binary16.unpackGetExponentUnbiased(k) == 16) {
+      if (Binary16.unpackGetSignificand(k) == 0) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  /**
+   * Return <code>true</code> if the given packed <code>binary16</code> value
+   * is not a number (<code>NaN</code>).
+   */
+
+  public static boolean isNaN(
+    final char k)
+  {
+    final int e = Binary16.unpackGetExponentUnbiased(k);
+    final int s = Binary16.unpackGetSignificand(k);
+    return (e == 16) && (s > 0);
+  }
+
+  /**
+   * <p>
+   * Convert a double precision floating point value to a packed
+   * <code>binary16</code> value.
+   * </p>
+   * <p>
+   * For the following specific cases, the function returns:
+   * </p>
+   * <ul>
+   * <li><code>NaN</code> iff <code>isNaN(k)</code></li>
+   * <li>{@link #POSITIVE_INFINITY} iff
+   * <code>k == {@link Double#POSITIVE_INFINITY}</code></li>
+   * <li>{@link #NEGATIVE_INFINITY} iff
+   * <code>k == {@link Double#NEGATIVE_INFINITY}</code></li>
+   * <li>{@link #NEGATIVE_ZERO} iff <code>k == -0.0</code></li>
+   * <li>{@link #POSITIVE_ZERO} iff <code>k == 0.0</code></li>
+   * </ul>
+   * <p>
+   * Otherwise, the <code>binary16</code> value that most closely represents
+   * <code>k</code> is returned. This may obviously be an infinite value as
+   * the interval of double precision values is far larger than that of the
+   * <code>binary16</code> type.
+   * </p>
+   * 
+   * @see #unpackDouble(char)
+   */
+
+  public static char packDouble(
+    final double k)
+  {
+    if (Double.isNaN(k)) {
+      return Binary16.exampleNaN();
+    }
+    if (k == Double.POSITIVE_INFINITY) {
+      return Binary16.POSITIVE_INFINITY;
+    }
+    if (k == Double.NEGATIVE_INFINITY) {
+      return Binary16.NEGATIVE_INFINITY;
+    }
+    if (Double.doubleToLongBits(k) == Binary64.NEGATIVE_ZERO_BITS) {
+      return Binary16.NEGATIVE_ZERO;
+    }
+    if (k == 0.0) {
+      return Binary16.POSITIVE_ZERO;
+    }
+
+    final long de = Binary64.unpackGetExponentUnbiased(k);
+    final long ds = Binary64.unpackGetSign(k);
+    final long dn = Binary64.unpackGetSignificand(k);
+    final char rsr = Binary16.packSetSignUnchecked((int) ds);
+
+    /**
+     * Extract the 5 least-significant bits of the exponent.
+     */
+
+    final int rem = (int) (de & 0x001F);
+    final char rer = Binary16.packSetExponentUnbiasedUnchecked(rem);
+
+    /**
+     * Extract the 10 most-significant bits of the significand.
+     */
+
+    final long rnm = dn & 0xFFC0000000000L;
+    final long rns = rnm >> 42;
+    final char rnr = Binary16.packSetSignificandUnchecked((int) rns);
+
+    /**
+     * Combine the results.
+     */
+
+    return (char) (rsr | rer | rnr);
+  }
+
+  /**
+   * <p>
+   * Convert a single precision floating point value to a packed
+   * <code>binary16</code> value.
+   * </p>
+   * <p>
+   * For the following specific cases, the function returns:
+   * </p>
+   * <ul>
+   * <li><code>NaN</code> iff <code>isNaN(k)</code></li>
+   * <li>{@link #POSITIVE_INFINITY} iff
+   * <code>k == {@link Float#POSITIVE_INFINITY}</code></li>
+   * <li>{@link #NEGATIVE_INFINITY} iff
+   * <code>k == {@link Float#NEGATIVE_INFINITY}</code></li>
+   * <li>{@link #NEGATIVE_ZERO} iff <code>k == -0.0</code></li>
+   * <li>{@link #POSITIVE_ZERO} iff <code>k == 0.0</code></li>
+   * </ul>
+   * <p>
+   * Otherwise, the <code>binary16</code> value that most closely represents
+   * <code>k</code> is returned. This may obviously be an infinite value as
+   * the interval of single precision values is far larger than that of the
+   * <code>binary16</code> type.
+   * </p>
+   * 
+   * @see #unpackFloat(char)
+   */
+
+  public static char packFloat(
+    final float k)
+  {
+    if (Float.isNaN(k)) {
+      return Binary16.exampleNaN();
+    }
+    if (k == Float.POSITIVE_INFINITY) {
+      return Binary16.POSITIVE_INFINITY;
+    }
+    if (k == Float.NEGATIVE_INFINITY) {
+      return Binary16.NEGATIVE_INFINITY;
+    }
+    if (Float.floatToIntBits(k) == Binary32.NEGATIVE_ZERO_BITS) {
+      return Binary16.NEGATIVE_ZERO;
+    }
+    if (k == 0.0) {
+      return Binary16.POSITIVE_ZERO;
+    }
+
+    final long de = Binary32.unpackGetExponentUnbiased(k);
+    final long ds = Binary32.unpackGetSign(k);
+    final long dn = Binary32.unpackGetSignificand(k);
+    final char rsr = Binary16.packSetSignUnchecked((int) ds);
+
+    /**
+     * Extract the 5 least-significant bits of the exponent.
+     */
+
+    final int rem = (int) (de & 0x001F);
+    final char rer = Binary16.packSetExponentUnbiasedUnchecked(rem);
+
+    /**
+     * Extract the 10 most-significant bits of the significand.
+     */
+
+    final long rnm = dn & 0x7FE000L;
+    final long rns = rnm >> 13;
+    final char rnr = Binary16.packSetSignificandUnchecked((int) rns);
+
+    /**
+     * Combine the results.
+     */
+
+    return (char) (rsr | rer | rnr);
+  }
+
+  /**
+   * <p>
+   * Encode the unbiased exponent <code>e</code>. Values should be in the
+   * range <code>[-15, 16]</code> - values outside of this range will be
+   * truncated.
+   * </p>
+   * 
+   * @see #unpackGetExponentUnbiased(char)
+   */
+
+  public static char packSetExponentUnbiasedUnchecked(
+    final int e)
+  {
+    final int eb = e + Binary16.BIAS;
+    final int es = eb << 10;
+    final int em = es & Binary16.MASK_EXPONENT;
+    return (char) em;
+  }
+
+  /**
+   * <p>
+   * Encode the significand <code>s</code>. Values should be in the range
+   * <code>[0, 1023]</code>. Values outside of this range will be truncated.
+   * </p>
+   * 
+   * @see #unpackGetSignificand(char)
+   */
+
+  public static char packSetSignificandUnchecked(
+    final int s)
+  {
+    final int sm = s & Binary16.MASK_SIGNIFICAND;
+    return (char) sm;
+  }
+
+  /**
+   * <p>
+   * Encode the sign bit <code>s</code>. Values should be in the range
+   * <code>[0, 1]</code>, with <code>0</code> ironically denoting a positive
+   * value. Values outside of this range will be truncated.
+   * </p>
+   * 
+   * @see #unpackGetSign(char)
+   */
+
+  public static char packSetSignUnchecked(
+    final int s)
+  {
+    final int ss = s << 15;
+    final int sm = ss & Binary16.MASK_SIGN;
+    return (char) sm;
+  }
+
+  /**
+   * Show the given raw packed <code>binary16</code> value as a string of
+   * binary digits.
+   */
+
+  public static String toRawBinaryString(
+    final char k)
+  {
+    final StringBuilder b = new StringBuilder();
+    int z = k;
+    for (int i = 0; i < 16; ++i) {
+      if ((z & 1) == 1) {
+        b.insert(0, "1");
+      } else {
+        b.insert(0, "0");
+      }
+      z >>= 1;
+    }
+    return b.toString();
+  }
+
+  /**
+   * <p>
+   * Convert a packed <code>binary16</code> value <code>k</code> to a
+   * double-precision floating point value.
+   * </p>
+   * <p>
+   * The function returns:
+   * </p>
+   * <ul>
+   * <li><code>NaN</code> iff <code>isNaN(k)</code></li>
+   * <li>{@link Double#POSITIVE_INFINITY} iff
+   * <code>k == {@link #POSITIVE_INFINITY}</code></li>
+   * <li>{@link Double#NEGATIVE_INFINITY} iff
+   * <code>k == {@link #NEGATIVE_INFINITY}</code></li>
+   * <li><code>-0.0</code> iff <code>k == {@link #NEGATIVE_ZERO}</code></li>
+   * <li><code>0.0</code> iff <code>k == {@link #POSITIVE_ZERO}</code></li>
+   * <li><code>(-1.0 * n) * (2 ^ e) * 1.s</code>, for the decoded sign
+   * <code>n</code> of <code>k</code>, the decoded exponent <code>e</code> of
+   * <code>k</code>, and the decoded significand <code>s</code> of
+   * <code>k</code>.</li>
+   * </ul>
+   * 
+   * @see #packDouble(double)
+   */
+
+  public static double unpackDouble(
+    final char k)
+  {
+    if (Binary16.isNaN(k)) {
+      return Double.NaN;
+    }
+    if (k == Binary16.POSITIVE_INFINITY) {
+      return Double.POSITIVE_INFINITY;
+    }
+    if (k == Binary16.NEGATIVE_INFINITY) {
+      return Double.NEGATIVE_INFINITY;
+    }
+    if (k == Binary16.NEGATIVE_ZERO) {
+      return -0.0;
+    }
+    if (k == Binary16.POSITIVE_ZERO) {
+      return 0.0;
+    }
+
+    final long e = Binary16.unpackGetExponentUnbiased(k);
+    final long s = Binary16.unpackGetSign(k);
+    final long n = Binary16.unpackGetSignificand(k);
+
+    /**
+     * Shift the sign bit to the position at which it will appear in the
+     * resulting value.
+     */
+
+    final long rsr = s << 63;
+
+    /**
+     * 1. Bias the exponent.
+     * 
+     * 2. Shift the result left to the position at which it will appear in the
+     * resulting value.
+     */
+
+    final long reb = (e + Binary64.BIAS);
+    final long rer = reb << 52;
+
+    /**
+     * Shift the significand left to the position at which it will appear in
+     * the resulting value.
+     */
+
+    final long rnr = n << 42;
+    return Double.longBitsToDouble(rsr | rer | rnr);
+  }
+
+  /**
+   * <p>
+   * Convert a packed <code>binary16</code> value <code>k</code> to a
+   * single-precision floating point value.
+   * </p>
+   * <p>
+   * The function returns:
+   * </p>
+   * <ul>
+   * <li><code>NaN</code> iff <code>isNaN(k)</code></li>
+   * <li>{@link Float#POSITIVE_INFINITY} iff
+   * <code>k == {@link #POSITIVE_INFINITY}</code></li>
+   * <li>{@link Float#NEGATIVE_INFINITY} iff
+   * <code>k == {@link #NEGATIVE_INFINITY}</code></li>
+   * <li><code>-0.0</code> iff <code>k == {@link #NEGATIVE_ZERO}</code></li>
+   * <li><code>0.0</code> iff <code>k == {@link #POSITIVE_ZERO}</code></li>
+   * <li><code>(-1.0 * n) * (2 ^ e) * 1.s</code>, for the decoded sign
+   * <code>n</code> of <code>k</code>, the decoded exponent <code>e</code> of
+   * <code>k</code>, and the decoded significand <code>s</code> of
+   * <code>k</code>.</li>
+   * </ul>
+   * 
+   * @see #packFloat(float)
+   */
+
+  public static float unpackFloat(
+    final char k)
+  {
+    if (Binary16.isNaN(k)) {
+      return Float.NaN;
+    }
+    if (k == Binary16.POSITIVE_INFINITY) {
+      return Float.POSITIVE_INFINITY;
+    }
+    if (k == Binary16.NEGATIVE_INFINITY) {
+      return Float.NEGATIVE_INFINITY;
+    }
+    if (k == Binary16.NEGATIVE_ZERO) {
+      return -0.0f;
+    }
+    if (k == Binary16.POSITIVE_ZERO) {
+      return 0.0f;
+    }
+
+    final int e = Binary16.unpackGetExponentUnbiased(k);
+    final int s = Binary16.unpackGetSign(k);
+    final int n = Binary16.unpackGetSignificand(k);
+
+    /**
+     * Shift the sign bit to the position at which it will appear in the
+     * resulting value.
+     */
+
+    final int rsr = s << 31;
+
+    /**
+     * 1. Bias the exponent.
+     * 
+     * 2. Shift the result left to the position at which it will appear in the
+     * resulting value.
+     */
+
+    final int reb = (e + Binary32.BIAS);
+    final int rer = reb << 23;
+
+    /**
+     * Shift the significand left to the position at which it will appear in
+     * the resulting value.
+     */
+
+    final int rnr = n << 13;
+    return Float.intBitsToFloat(rsr | rer | rnr);
+  }
+
+  /**
+   * <p>
+   * Extract and unbias the exponent of the given packed <code>binary16</code>
+   * value.
+   * </p>
+   * <p>
+   * The exponent is encoded <i>biased</i> as a number in the range
+   * <code>[0, 31]</code>, with <code>0</code> indicating that the number is
+   * <i>subnormal</i> and <code>[1, 30]</code> denoting the actual exponent
+   * plus {@link #BIAS}. Infinite and <code>NaN</code> values always have an
+   * exponent of <code>31</code>.
+   * </p>
+   * <p>
+   * This function will therefore return:
+   * </p>
+   * <ul>
+   * <li>
+   * <code>0 - {@link #BIAS} = -15</code> iff the input is a <i>subnormal</i>
+   * number.</li>
+   * <li>An integer in the range
+   * <code>[1 - {@link #BIAS}, 30 - {@link #BIAS}] = [-14, 15]</code> iff the
+   * input is a <i>normal</i> number.</li>
+   * <li>
+   * <code>16</code> iff the input is {@link #POSITIVE_INFINITY},
+   * {@link #NEGATIVE_INFINITY}, or <code>NaN</code>.</li>
+   * </ul>
+   * 
+   * @see #packSetExponentUnbiasedUnchecked(int)
+   */
+
+  public static int unpackGetExponentUnbiased(
+    final char k)
+  {
+    final int em = k & Binary16.MASK_EXPONENT;
+    final int es = em >> 10;
+    return es - Binary16.BIAS;
+  }
+
+  /**
+   * Retrieve the sign bit of the given packed <code>binary16</code> value, as
+   * an integer in the range <code>[0, 1]</code>.
+   * 
+   * @see Binary16#packSetSignUnchecked(int)
+   */
+
+  public static int unpackGetSign(
+    final char k)
+  {
+    return (k & Binary16.MASK_SIGN) >> 15;
+  }
+
+  /**
+   * <p>
+   * Return the significand of the given packed <code>binary16</code> value as
+   * an integer in the range <code>[0, 1023]</code>.
+   * </p>
+   * 
+   * @see Binary16#packSetSignificandUnchecked(int)
+   */
+
+  public static int unpackGetSignificand(
+    final char k)
+  {
+    return k & Binary16.MASK_SIGNIFICAND;
+  }
+
+  private Binary16()
+  {
+    throw new AssertionError("Unreachable code, report this bug!");
+  }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Binary32.java b/src/jogl/classes/com/jogamp/opengl/math/Binary32.java
new file mode 100644
index 000000000..d98815d9f
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Binary32.java
@@ -0,0 +1,116 @@
+/**
+ * Copyright 2013 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+
+package com.jogamp.opengl.math;
+
+/**
+ * Functions for interrogating <code>binary32</code> (float) values.
+ */
+
+public final class Binary32
+{
+  static final int NEGATIVE_ZERO_BITS;
+  static final int MASK_SIGN;
+  static final int MASK_EXPONENT;
+  static final int MASK_SIGNIFICAND;
+  static final int BIAS;
+
+  static {
+    NEGATIVE_ZERO_BITS = 0x80000000;
+    MASK_SIGN = 0x80000000;
+    MASK_EXPONENT = 0x7ff00000;
+    MASK_SIGNIFICAND = 0x7fffff;
+    BIAS = 127;
+  }
+
+  /**
+   * <p>
+   * Extract and unbias the exponent of the given packed <code>float</code>
+   * value.
+   * </p>
+   * <p>
+   * The exponent is encoded <i>biased</i> as a number in the range
+   * <code>[0, 255]</code>, with <code>0</code> indicating that the number is
+   * <i>subnormal</i> and <code>[1, 254]</code> denoting the actual exponent
+   * plus {@link #BIAS}. Infinite and <code>NaN</code> values always have a
+   * biased exponent of <code>255</code>.
+   * </p>
+   * <p>
+   * This function will therefore return:
+   * </p>
+   * <ul>
+   * <li>
+   * <code>0 - {@link #BIAS} = -127</code> iff the input is a <i>subnormal</i>
+   * number.</li>
+   * <li>An integer in the range
+   * <code>[1 - {@link #BIAS}, 254 - {@link #BIAS}] = [-126, 127]</code> iff
+   * the input is a <i>normal</i> number.</li>
+   * <li>
+   * <code>255 - {@link #BIAS} = 128</code> iff the input is
+   * {@link #POSITIVE_INFINITY}, {@link #NEGATIVE_INFINITY}, or
+   * <code>NaN</code>.</li>
+   * </ul>
+   * 
+   * @see #packSetExponentUnbiasedUnchecked(int)
+   */
+
+  public static int unpackGetExponentUnbiased(
+    final float d)
+  {
+    final int b = Float.floatToRawIntBits(d);
+    final int em = b & Binary32.MASK_EXPONENT;
+    final int es = em >> 23;
+    return es - Binary32.BIAS;
+  }
+
+  /**
+   * <p>
+   * Return the sign of the given float value.
+   * </p>
+   */
+
+  public static int unpackGetSign(
+    final float d)
+  {
+    final int b = Float.floatToRawIntBits(d);
+    return ((b & Binary32.MASK_SIGN) >> 31) & 1;
+  }
+
+  /**
+   * <p>
+   * Return the significand of the given float value.
+   * </p>
+   */
+
+  public static int unpackGetSignificand(
+    final float d)
+  {
+    final int b = Float.floatToRawIntBits(d);
+    return b & Binary32.MASK_SIGNIFICAND;
+  }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Binary64.java b/src/jogl/classes/com/jogamp/opengl/math/Binary64.java
new file mode 100644
index 000000000..5efad433a
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Binary64.java
@@ -0,0 +1,116 @@
+/**
+ * Copyright 2013 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+
+package com.jogamp.opengl.math;
+
+/**
+ * Functions for interrogating <code>binary64</code> (double) values.
+ */
+
+public final class Binary64
+{
+  static final long NEGATIVE_ZERO_BITS;
+  static final long MASK_SIGN;
+  static final long MASK_EXPONENT;
+  static final long MASK_SIGNIFICAND;
+  static final long BIAS;
+
+  static {
+    NEGATIVE_ZERO_BITS = 0x8000000000000000L;
+    MASK_SIGN = 0x8000000000000000L;
+    MASK_EXPONENT = 0x7ff0000000000000L;
+    MASK_SIGNIFICAND = 0x000fffffffffffffL;
+    BIAS = 1023;
+  }
+
+  /**
+   * <p>
+   * Extract and unbias the exponent of the given packed <code>double</code>
+   * value.
+   * </p>
+   * <p>
+   * The exponent is encoded <i>biased</i> as a number in the range
+   * <code>[0, 2047]</code>, with <code>0</code> indicating that the number is
+   * <i>subnormal</i> and <code>[1, 2046]</code> denoting the actual exponent
+   * plus {@link #BIAS}. Infinite and <code>NaN</code> values always have a
+   * biased exponent of <code>2047</code>.
+   * </p>
+   * <p>
+   * This function will therefore return:
+   * </p>
+   * <ul>
+   * <li>
+   * <code>0 - {@link #BIAS} = -1023</code> iff the input is a
+   * <i>subnormal</i> number.</li>
+   * <li>An integer in the range
+   * <code>[1 - {@link #BIAS}, 2046 - {@link #BIAS}] = [-1022, 1023]</code>
+   * iff the input is a <i>normal</i> number.</li>
+   * <li>
+   * <code>2047 - {@link #BIAS} = 1024</code> iff the input is
+   * {@link #POSITIVE_INFINITY}, {@link #NEGATIVE_INFINITY}, or
+   * <code>NaN</code>.</li>
+   * </ul>
+   * 
+   * @see #packSetExponentUnbiasedUnchecked(int)
+   */
+
+  public static long unpackGetExponentUnbiased(
+    final double d)
+  {
+    final long b = Double.doubleToRawLongBits(d);
+    final long em = b & Binary64.MASK_EXPONENT;
+    final long es = em >> 52;
+    return es - Binary64.BIAS;
+  }
+
+  /**
+   * <p>
+   * Return the significand of the given double value.
+   * </p>
+   */
+  
+  public static long unpackGetSignificand(
+    final double d)
+  {
+    final long b = Double.doubleToRawLongBits(d);
+    return b & Binary64.MASK_SIGNIFICAND;
+  }
+
+  /**
+   * <p>
+   * Return the sign of the given double value.
+   * </p>
+   */
+  
+  public static long unpackGetSign(
+    final double d)
+  {
+    final long b = Double.doubleToRawLongBits(d);
+    return ((b & Binary64.MASK_SIGN) >> 63) & 1;
+  }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/FixedPoint.java b/src/jogl/classes/com/jogamp/opengl/math/FixedPoint.java
new file mode 100644
index 000000000..b7dbf183f
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/FixedPoint.java
@@ -0,0 +1,61 @@
+/*
+ * Copyright (c) 2009 Sun Microsystems, Inc. All Rights Reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * - Redistribution of source code must retain the above copyright
+ *   notice, this list of conditions and the following disclaimer.
+ *
+ * - Redistribution in binary form must reproduce the above copyright
+ *   notice, this list of conditions and the following disclaimer in the
+ *   documentation and/or other materials provided with the distribution.
+ *
+ * Neither the name of Sun Microsystems, Inc. or the names of
+ * contributors may be used to endorse or promote products derived from
+ * this software without specific prior written permission.
+ *
+ * This software is provided "AS IS," without a warranty of any kind. ALL
+ * EXPRESS OR IMPLIED CONDITIONS, REPRESENTATIONS AND WARRANTIES,
+ * INCLUDING ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS FOR A
+ * PARTICULAR PURPOSE OR NON-INFRINGEMENT, ARE HEREBY EXCLUDED. SUN
+ * MICROSYSTEMS, INC. ("SUN") AND ITS LICENSORS SHALL NOT BE LIABLE FOR
+ * ANY DAMAGES SUFFERED BY LICENSEE AS A RESULT OF USING, MODIFYING OR
+ * DISTRIBUTING THIS SOFTWARE OR ITS DERIVATIVES. IN NO EVENT WILL SUN OR
+ * ITS LICENSORS BE LIABLE FOR ANY LOST REVENUE, PROFIT OR DATA, OR FOR
+ * DIRECT, INDIRECT, SPECIAL, CONSEQUENTIAL, INCIDENTAL OR PUNITIVE
+ * DAMAGES, HOWEVER CAUSED AND REGARDLESS OF THE THEORY OF LIABILITY,
+ * ARISING OUT OF THE USE OF OR INABILITY TO USE THIS SOFTWARE, EVEN IF
+ * SUN HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
+ *
+ */
+
+package com.jogamp.opengl.math;
+
+public class FixedPoint {
+  public static final int toFixed(int value) {
+    if (value < -32768) value = -32768;
+    if (value > 32767) value = 32767;
+    return value * 65536;
+  }
+
+  public static final int toFixed(float value) {
+    if (value < -32768) value = -32768;
+    if (value > 32767) value = 32767;
+    return (int)(value * 65536.0f);
+  }
+
+  public static final float toFloat(int value) {
+    return (float)value/65536.0f;
+  }
+
+  public static final int mult(int x1, int x2) {
+    return (int) ( ((long)x1*(long)x2)/65536 );
+  }
+
+  public static final int div(int x1, int x2) {
+    return (int) ( (((long)x1)<<16)/x2 );
+  }
+}
+
diff --git a/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
new file mode 100644
index 000000000..191a83241
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
@@ -0,0 +1,573 @@
+/**
+ * Copyright 2010 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math;
+
+import java.nio.FloatBuffer;
+
+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.
+ * </p>
+ * <p>
+ * Derived from ProjectFloat.java - Created 11-jan-2004
+ * </p>
+ *
+ * @author Erik Duijs
+ * @author Kenneth Russell
+ * @author Sven Gothel
+ */
+public class FloatUtil {
+  private static final float[] IDENTITY_MATRIX =
+    new float[] {
+      1.0f, 0.0f, 0.0f, 0.0f,
+      0.0f, 1.0f, 0.0f, 0.0f,
+      0.0f, 0.0f, 1.0f, 0.0f,
+      0.0f, 0.0f, 0.0f, 1.0f };
+
+  private static final float[] ZERO_MATRIX =
+    new float[] {
+      0.0f, 0.0f, 0.0f, 0.0f,
+      0.0f, 0.0f, 0.0f, 0.0f,
+      0.0f, 0.0f, 0.0f, 0.0f,
+      0.0f, 0.0f, 0.0f, 0.0f };
+
+  /**
+   * Make matrix an identity matrix
+   */
+  public static final void makeIdentityf(float[] m, int offset) {
+    for (int i = 0; i < 16; i++) {
+      m[i+offset] = IDENTITY_MATRIX[i];
+    }
+  }
+
+  /**
+   * Make matrix an identity matrix
+   */
+  public static final void makeIdentityf(FloatBuffer m) {
+    final int oldPos = m.position();
+    m.put(IDENTITY_MATRIX);
+    m.position(oldPos);
+  }
+
+  /**
+   * Make matrix an zero matrix
+   */
+  public static final void makeZero(float[] m, int offset) {
+    for (int i = 0; i < 16; i++) {
+      m[i+offset] = 0;
+    }
+  }
+
+  /**
+   * Make matrix an zero matrix
+   */
+  public static final void makeZero(FloatBuffer m) {
+    final int oldPos = m.position();
+    m.put(ZERO_MATRIX);
+    m.position(oldPos);
+  }
+
+  /**
+   * @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
+   */
+  public static final void multMatrixf(final float[] a, int a_off, final float[] b, int b_off, float[] d, int d_off) {
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final float ai0=a[a_off+i+0*4],  ai1=a[a_off+i+1*4],  ai2=a[a_off+i+2*4],  ai3=a[a_off+i+3*4]; // row-i of a
+        d[d_off+i+0*4] = ai0 * b[b_off+0+0*4] + ai1 * b[b_off+1+0*4] + ai2 * b[b_off+2+0*4] + ai3 * b[b_off+3+0*4] ;
+        d[d_off+i+1*4] = ai0 * b[b_off+0+1*4] + ai1 * b[b_off+1+1*4] + ai2 * b[b_off+2+1*4] + ai3 * b[b_off+3+1*4] ;
+        d[d_off+i+2*4] = ai0 * b[b_off+0+2*4] + ai1 * b[b_off+1+2*4] + ai2 * b[b_off+2+2*4] + ai3 * b[b_off+3+2*4] ;
+        d[d_off+i+3*4] = ai0 * b[b_off+0+3*4] + ai1 * b[b_off+1+3*4] + ai2 * b[b_off+2+3*4] + ai3 * b[b_off+3+3*4] ;
+     }
+  }
+
+  /**
+   * @param a 4x4 matrix in column-major order (also result)
+   * @param b 4x4 matrix in column-major order
+   */
+  public static final void multMatrixf(final float[] a, int a_off, final float[] b, int b_off) {
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final int a_off_i = a_off+i;
+        final float ai0=a[a_off_i+0*4],  ai1=a[a_off_i+1*4],  ai2=a[a_off_i+2*4],  ai3=a[a_off_i+3*4]; // row-i of a
+        a[a_off_i+0*4] = ai0 * b[b_off+0+0*4] + ai1 * b[b_off+1+0*4] + ai2 * b[b_off+2+0*4] + ai3 * b[b_off+3+0*4] ;
+        a[a_off_i+1*4] = ai0 * b[b_off+0+1*4] + ai1 * b[b_off+1+1*4] + ai2 * b[b_off+2+1*4] + ai3 * b[b_off+3+1*4] ;
+        a[a_off_i+2*4] = ai0 * b[b_off+0+2*4] + ai1 * b[b_off+1+2*4] + ai2 * b[b_off+2+2*4] + ai3 * b[b_off+3+2*4] ;
+        a[a_off_i+3*4] = ai0 * b[b_off+0+3*4] + ai1 * b[b_off+1+3*4] + ai2 * b[b_off+2+3*4] + ai3 * b[b_off+3+3*4] ;
+     }
+  }
+
+  /**
+   * @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
+   */
+  public static final void multMatrixf(final float[] a, int a_off, final float[] b, int b_off, FloatBuffer d) {
+     final int dP = d.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final float ai0=a[a_off+i+0*4],  ai1=a[a_off+i+1*4],  ai2=a[a_off+i+2*4],  ai3=a[a_off+i+3*4]; // row-i of a
+        d.put(dP+i+0*4 , ai0 * b[b_off+0+0*4] + ai1 * b[b_off+1+0*4] + ai2 * b[b_off+2+0*4] + ai3 * b[b_off+3+0*4] );
+        d.put(dP+i+1*4 , ai0 * b[b_off+0+1*4] + ai1 * b[b_off+1+1*4] + ai2 * b[b_off+2+1*4] + ai3 * b[b_off+3+1*4] );
+        d.put(dP+i+2*4 , ai0 * b[b_off+0+2*4] + ai1 * b[b_off+1+2*4] + ai2 * b[b_off+2+2*4] + ai3 * b[b_off+3+2*4] );
+        d.put(dP+i+3*4 , ai0 * b[b_off+0+3*4] + ai1 * b[b_off+1+3*4] + ai2 * b[b_off+2+3*4] + ai3 * b[b_off+3+3*4] );
+     }
+  }
+
+  /**
+   * @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
+   */
+  public static final void multMatrixf(final FloatBuffer a, final float[] b, int b_off, FloatBuffer d) {
+     final int aP = a.position();
+     final int dP = d.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final float ai0=a.get(aP+i+0*4),  ai1=a.get(aP+i+1*4),  ai2=a.get(aP+i+2*4),  ai3=a.get(aP+i+3*4); // row-i of a
+        d.put(dP+i+0*4 , ai0 * b[b_off+0+0*4] + ai1 * b[b_off+1+0*4] + ai2 * b[b_off+2+0*4] + ai3 * b[b_off+3+0*4] );
+        d.put(dP+i+1*4 , ai0 * b[b_off+0+1*4] + ai1 * b[b_off+1+1*4] + ai2 * b[b_off+2+1*4] + ai3 * b[b_off+3+1*4] );
+        d.put(dP+i+2*4 , ai0 * b[b_off+0+2*4] + ai1 * b[b_off+1+2*4] + ai2 * b[b_off+2+2*4] + ai3 * b[b_off+3+2*4] );
+        d.put(dP+i+3*4 , ai0 * b[b_off+0+3*4] + ai1 * b[b_off+1+3*4] + ai2 * b[b_off+2+3*4] + ai3 * b[b_off+3+3*4] );
+     }
+  }
+
+  /**
+   * @param a 4x4 matrix in column-major order (also result)
+   * @param b 4x4 matrix in column-major order
+   */
+  public static final void multMatrixf(final FloatBuffer a, final float[] b, int b_off) {
+     final int aP = a.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final int aP_i = aP+i;
+        final float ai0=a.get(aP_i+0*4),  ai1=a.get(aP_i+1*4),  ai2=a.get(aP_i+2*4),  ai3=a.get(aP_i+3*4); // row-i of a
+        a.put(aP_i+0*4 , ai0 * b[b_off+0+0*4] + ai1 * b[b_off+1+0*4] + ai2 * b[b_off+2+0*4] + ai3 * b[b_off+3+0*4] );
+        a.put(aP_i+1*4 , ai0 * b[b_off+0+1*4] + ai1 * b[b_off+1+1*4] + ai2 * b[b_off+2+1*4] + ai3 * b[b_off+3+1*4] );
+        a.put(aP_i+2*4 , ai0 * b[b_off+0+2*4] + ai1 * b[b_off+1+2*4] + ai2 * b[b_off+2+2*4] + ai3 * b[b_off+3+2*4] );
+        a.put(aP_i+3*4 , ai0 * b[b_off+0+3*4] + ai1 * b[b_off+1+3*4] + ai2 * b[b_off+2+3*4] + ai3 * b[b_off+3+3*4] );
+     }
+  }
+
+  /**
+   * @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
+   */
+  public static final void multMatrixf(final FloatBuffer a, final FloatBuffer b, FloatBuffer d) {
+     final int aP = a.position();
+     final int bP = b.position();
+     final int dP = d.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final float ai0=a.get(aP+i+0*4),  ai1=a.get(aP+i+1*4),  ai2=a.get(aP+i+2*4),  ai3=a.get(aP+i+3*4); // row-i of a
+        d.put(dP+i+0*4 , ai0 * b.get(bP+0+0*4) + ai1 * b.get(bP+1+0*4) + ai2 * b.get(bP+2+0*4) + ai3 * b.get(bP+3+0*4) );
+        d.put(dP+i+1*4 , ai0 * b.get(bP+0+1*4) + ai1 * b.get(bP+1+1*4) + ai2 * b.get(bP+2+1*4) + ai3 * b.get(bP+3+1*4) );
+        d.put(dP+i+2*4 , ai0 * b.get(bP+0+2*4) + ai1 * b.get(bP+1+2*4) + ai2 * b.get(bP+2+2*4) + ai3 * b.get(bP+3+2*4) );
+        d.put(dP+i+3*4 , ai0 * b.get(bP+0+3*4) + ai1 * b.get(bP+1+3*4) + ai2 * b.get(bP+2+3*4) + ai3 * b.get(bP+3+3*4) );
+     }
+  }
+
+  /**
+   * @param a 4x4 matrix in column-major order (also result)
+   * @param b 4x4 matrix in column-major order
+   */
+  public static final void multMatrixf(final FloatBuffer a, final FloatBuffer b) {
+     final int aP = a.position();
+     final int bP = b.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final int aP_i = aP+i;
+        final float ai0=a.get(aP_i+0*4),  ai1=a.get(aP_i+1*4),  ai2=a.get(aP_i+2*4),  ai3=a.get(aP_i+3*4); // row-i of a
+        a.put(aP_i+0*4 , ai0 * b.get(bP+0+0*4) + ai1 * b.get(bP+1+0*4) + ai2 * b.get(bP+2+0*4) + ai3 * b.get(bP+3+0*4) );
+        a.put(aP_i+1*4 , ai0 * b.get(bP+0+1*4) + ai1 * b.get(bP+1+1*4) + ai2 * b.get(bP+2+1*4) + ai3 * b.get(bP+3+1*4) );
+        a.put(aP_i+2*4 , ai0 * b.get(bP+0+2*4) + ai1 * b.get(bP+1+2*4) + ai2 * b.get(bP+2+2*4) + ai3 * b.get(bP+3+2*4) );
+        a.put(aP_i+3*4 , ai0 * b.get(bP+0+3*4) + ai1 * b.get(bP+1+3*4) + ai2 * b.get(bP+2+3*4) + ai3 * b.get(bP+3+3*4) );
+     }
+  }
+
+  /**
+   * @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
+   */
+  public static final void multMatrixf(final FloatBuffer a, final FloatBuffer b, float[] d, int d_off) {
+     final int aP = a.position();
+     final int bP = b.position();
+     for (int i = 0; i < 4; i++) {
+        // one row in column-major order
+        final float ai0=a.get(aP+i+0*4),  ai1=a.get(aP+i+1*4),  ai2=a.get(aP+i+2*4),  ai3=a.get(aP+i+3*4); // row-i of a
+        d[d_off+i+0*4] = ai0 * b.get(bP+0+0*4) + ai1 * b.get(bP+1+0*4) + ai2 * b.get(bP+2+0*4) + ai3 * b.get(bP+3+0*4) ;
+        d[d_off+i+1*4] = ai0 * b.get(bP+0+1*4) + ai1 * b.get(bP+1+1*4) + ai2 * b.get(bP+2+1*4) + ai3 * b.get(bP+3+1*4) ;
+        d[d_off+i+2*4] = ai0 * b.get(bP+0+2*4) + ai1 * b.get(bP+1+2*4) + ai2 * b.get(bP+2+2*4) + ai3 * b.get(bP+3+2*4) ;
+        d[d_off+i+3*4] = ai0 * b.get(bP+0+3*4) + ai1 * b.get(bP+1+3*4) + ai2 * b.get(bP+2+3*4) + ai3 * b.get(bP+3+3*4) ;
+     }
+  }
+
+  /**
+   * 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
+   * @param v_out m_in * v_in
+   */
+  public static final void multMatrixVecf(float[] m_in, int m_in_off, float[] v_in, int v_in_off, float[] v_out, int v_out_off) {
+    for (int i = 0; i < 4; i++) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[i + v_out_off] =
+        v_in[0+v_in_off] * m_in[0*4+i+m_in_off] +
+        v_in[1+v_in_off] * m_in[1*4+i+m_in_off] +
+        v_in[2+v_in_off] * m_in[2*4+i+m_in_off] +
+        v_in[3+v_in_off] * m_in[3*4+i+m_in_off];
+    }
+  }
+
+  /**
+   * @param m_in 4x4 matrix in column-major order
+   * @param m_in_off
+   * @param v_in 4-component column-vector
+   * @param v_out m_in * v_in
+   */
+  public static final void multMatrixVecf(float[] m_in, float[] v_in, float[] v_out) {
+    for (int i = 0; i < 4; i++) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[i] =
+        v_in[0] * m_in[0*4+i] +
+        v_in[1] * m_in[1*4+i] +
+        v_in[2] * m_in[2*4+i] +
+        v_in[3] * m_in[3*4+i];
+    }
+  }
+
+  /**
+   * @param m_in 4x4 matrix in column-major order
+   * @param v_in 4-component column-vector
+   * @param v_out m_in * v_in
+   */
+  public static final void multMatrixVecf(FloatBuffer m_in, float[] v_in, int v_in_off, float[] v_out, int v_out_off) {
+    final int matrixPos = m_in.position();
+    for (int i = 0; i < 4; i++) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[i+v_out_off] =
+        v_in[0+v_in_off] * m_in.get(0*4+i+matrixPos) +
+        v_in[1+v_in_off] * m_in.get(1*4+i+matrixPos) +
+        v_in[2+v_in_off] * m_in.get(2*4+i+matrixPos) +
+        v_in[3+v_in_off] * m_in.get(3*4+i+matrixPos);
+    }
+  }
+
+  /**
+   * @param m_in 4x4 matrix in column-major order
+   * @param v_in 4-component column-vector
+   * @param v_out m_in * v_in
+   */
+  public static final void multMatrixVecf(FloatBuffer m_in, float[] v_in, float[] v_out) {
+    final int matrixPos = m_in.position();
+    for (int i = 0; i < 4; i++) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[i] =
+        v_in[0] * m_in.get(0*4+i+matrixPos) +
+        v_in[1] * m_in.get(1*4+i+matrixPos) +
+        v_in[2] * m_in.get(2*4+i+matrixPos) +
+        v_in[3] * m_in.get(3*4+i+matrixPos);
+    }
+  }
+
+  /**
+   * @param m_in 4x4 matrix in column-major order
+   * @param v_in 4-component column-vector
+   * @param v_out m_in * v_in
+   */
+  public static final void multMatrixVecf(FloatBuffer m_in, FloatBuffer v_in, FloatBuffer v_out) {
+    final int inPos = v_in.position();
+    final int outPos = v_out.position();
+    final int matrixPos = m_in.position();
+    for (int i = 0; i < 4; i++) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out.put(i + outPos,
+              v_in.get(0+inPos) * m_in.get(0*4+i+matrixPos) +
+              v_in.get(1+inPos) * m_in.get(1*4+i+matrixPos) +
+              v_in.get(2+inPos) * m_in.get(2*4+i+matrixPos) +
+              v_in.get(3+inPos) * m_in.get(3*4+i+matrixPos));
+    }
+  }
+
+  /**
+   * @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)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @param row row number to print
+   * @return matrix row string representation
+   */
+  public static StringBuilder matrixRowToString(StringBuilder sb, String f, FloatBuffer a, int aOffset, int rows, int columns, boolean rowMajorOrder, int row) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      final int a0 = aOffset + a.position();
+      if(rowMajorOrder) {
+          for(int c=0; c<columns; c++) {
+              sb.append( String.format( f+" ", a.get( a0 + row*columns + c ) ) );
+          }
+      } else {
+          for(int r=0; r<columns; r++) {
+              sb.append( String.format( f+" ", a.get( a0 + row + r*rows ) ) );
+          }
+      }
+      return sb;
+  }
+
+  /**
+   * @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)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @param row row number to print
+   * @return matrix row string representation
+   */
+  public static StringBuilder matrixRowToString(StringBuilder sb, String f, float[] a, int aOffset, int rows, int columns, boolean rowMajorOrder, int row) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      if(rowMajorOrder) {
+          for(int c=0; c<columns; c++) {
+              sb.append( String.format( f+" ", a[ aOffset + row*columns + c ] ) );
+          }
+      } else {
+          for(int r=0; r<columns; r++) {
+              sb.append( String.format( f+" ", a[ aOffset + row + r*rows ] ) );
+          }
+      }
+      return sb;
+  }
+
+  /**
+   * @param sb optional passed StringBuilder instance to be used
+   * @param rowPrefix optional prefix for each row
+   * @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)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @return matrix string representation
+   */
+  public static StringBuilder matrixToString(StringBuilder sb, String rowPrefix, String f, FloatBuffer a, int aOffset, int rows, int columns, boolean rowMajorOrder) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      final String prefix = ( null == rowPrefix ) ? "" : rowPrefix;
+      for(int i=0; i<rows; i++) {
+          sb.append(prefix).append("[ ");
+          matrixRowToString(sb, f, a, aOffset, rows, columns, rowMajorOrder, i);
+          sb.append("]").append(Platform.getNewline());
+      }
+      return sb;
+  }
+
+  /**
+   * @param sb optional passed StringBuilder instance to be used
+   * @param rowPrefix optional prefix for each row
+   * @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)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @return matrix string representation
+   */
+  public static StringBuilder matrixToString(StringBuilder sb, String rowPrefix, String f, float[] a, int aOffset, int rows, int columns, boolean rowMajorOrder) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      final String prefix = ( null == rowPrefix ) ? "" : rowPrefix;
+      for(int i=0; i<rows; i++) {
+          sb.append(prefix).append("[ ");
+          matrixRowToString(sb, f, a, aOffset, rows, columns, rowMajorOrder, i);
+          sb.append("]").append(Platform.getNewline());
+      }
+      return sb;
+  }
+
+  /**
+   * @param sb optional passed StringBuilder instance to be used
+   * @param rowPrefix optional prefix for each row
+   * @param f the format string of one floating point, i.e. "%10.5f", see {@link java.util.Formatter}
+   * @param a 4x4 matrix in column major order (OpenGL)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param b 4x4 matrix in column major order (OpenGL)
+   * @param bOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @return side by side representation
+   */
+  public static StringBuilder matrixToString(StringBuilder sb, String rowPrefix, String f, FloatBuffer a, int aOffset, FloatBuffer b, int bOffset, int rows, int columns, boolean rowMajorOrder) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      final String prefix = ( null == rowPrefix ) ? "" : rowPrefix;
+      for(int i=0; i<rows; i++) {
+          sb.append(prefix).append("[ ");
+          matrixRowToString(sb, f, a, aOffset, rows, columns, rowMajorOrder, i);
+          sb.append("=?= ");
+          matrixRowToString(sb, f, b, bOffset, rows, columns, rowMajorOrder, i);
+          sb.append("]").append(Platform.getNewline());
+      }
+      return sb;
+  }
+
+  /**
+   * @param sb optional passed StringBuilder instance to be used
+   * @param rowPrefix optional prefix for each row
+   * @param f the format string of one floating point, i.e. "%10.5f", see {@link java.util.Formatter}
+   * @param a 4x4 matrix in column major order (OpenGL)
+   * @param aOffset offset to <code>a</code>'s current position
+   * @param b 4x4 matrix in column major order (OpenGL)
+   * @param bOffset offset to <code>a</code>'s current position
+   * @param rows
+   * @param columns
+   * @param rowMajorOrder if true floats are layed out in row-major-order, otherwise column-major-order (OpenGL)
+   * @return side by side representation
+   */
+  public static StringBuilder matrixToString(StringBuilder sb, String rowPrefix, String f, float[] a, int aOffset, float[] b, int bOffset, int rows, int columns, boolean rowMajorOrder) {
+      if(null == sb) {
+          sb = new StringBuilder();
+      }
+      final String prefix = ( null == rowPrefix ) ? "" : rowPrefix;
+      for(int i=0; i<rows; i++) {
+          sb.append(prefix).append("[ ");
+          matrixRowToString(sb, f, a, aOffset, rows, columns, rowMajorOrder, i);
+          sb.append("=?= ");
+          matrixRowToString(sb, f, b, bOffset, rows, columns, rowMajorOrder, i);
+          sb.append("]").append(Platform.getNewline());
+      }
+      return sb;
+  }
+
+  public static final float E = 2.7182818284590452354f;
+
+  public static final float PI = 3.14159265358979323846f;
+
+  public static float abs(float a) { return (float) java.lang.Math.abs(a);  }
+
+  public static float pow(float a, float b) { return (float) java.lang.Math.pow(a, b);  }
+
+  public static float sin(float a) { return (float) java.lang.Math.sin(a);  }
+
+  public static float cos(float a) { return (float) java.lang.Math.cos(a);  }
+
+  public static float acos(float a) { return (float) java.lang.Math.acos(a);  }
+
+  public static float sqrt(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/Quaternion.java b/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java
new file mode 100644
index 000000000..52a59c599
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java
@@ -0,0 +1,414 @@
+/**
+ * Copyright 2010 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math;
+
+public class Quaternion {
+    protected float x, y, z, w;
+
+    public Quaternion() {
+        setIdentity();
+    }
+
+    public Quaternion(Quaternion q) {
+        x = q.x;
+        y = q.y;
+        z = q.z;
+        w = q.w;
+    }
+
+    public Quaternion(float x, float y, float z, float w) {
+        this.x = x;
+        this.y = y;
+        this.z = z;
+        this.w = w;
+    }
+
+    /**
+     * Constructor to create a rotation based quaternion from two vectors
+     *
+     * @param vector1
+     * @param vector2
+     */
+    public Quaternion(float[] vector1, float[] vector2) {
+        final float theta = FloatUtil.acos(VectorUtil.dot(vector1, vector2));
+        final float[] cross = VectorUtil.cross(vector1, vector2);
+        fromAxis(cross, theta);
+    }
+
+    /***
+     * Constructor to create a rotation based quaternion from axis vector and angle
+     * @param vector axis vector
+     * @param angle rotation angle (rads)
+     * @see #fromAxis(float[], float)
+     */
+    public Quaternion(float[] vector, float angle) {
+        fromAxis(vector, angle);
+    }
+
+    /***
+     * Initialize this quaternion with given axis vector and rotation angle
+     *
+     * @param vector axis vector
+     * @param angle rotation angle (rads)
+     */
+    public void fromAxis(float[] vector, float angle) {
+        final float halfangle = angle * 0.5f;
+        final float sin = FloatUtil.sin(halfangle);
+        final float[] nv = VectorUtil.normalize(vector);
+        x = (nv[0] * sin);
+        y = (nv[1] * sin);
+        z = (nv[2] * sin);
+        w = FloatUtil.cos(halfangle);
+    }
+
+    /**
+     * Transform the rotational quaternion to axis based rotation angles
+     *
+     * @return new float[4] with ,theta,Rx,Ry,Rz
+     */
+    public float[] toAxis() {
+        final float[] vec = new float[4];
+        final float scale = FloatUtil.sqrt(x * x + y * y + z * z);
+        vec[0] = FloatUtil.acos(w) * 2.0f;
+        vec[1] = x / scale;
+        vec[2] = y / scale;
+        vec[3] = z / scale;
+        return vec;
+    }
+
+    public float getW() {
+        return w;
+    }
+
+    public void setW(float w) {
+        this.w = w;
+    }
+
+    public float getX() {
+        return x;
+    }
+
+    public void setX(float x) {
+        this.x = x;
+    }
+
+    public float getY() {
+        return y;
+    }
+
+    public void setY(float y) {
+        this.y = y;
+    }
+
+    public float getZ() {
+        return z;
+    }
+
+    public void setZ(float z) {
+        this.z = z;
+    }
+
+    /**
+     * Add a quaternion
+     *
+     * @param q quaternion
+     */
+    public void add(Quaternion q) {
+        x += q.x;
+        y += q.y;
+        z += q.z;
+    }
+
+    /**
+     * Subtract a quaternion
+     *
+     * @param q quaternion
+     */
+    public void subtract(Quaternion q) {
+        x -= q.x;
+        y -= q.y;
+        z -= q.z;
+    }
+
+    /**
+     * Divide a quaternion by a constant
+     *
+     * @param n a float to divide by
+     */
+    public void divide(float n) {
+        x /= n;
+        y /= n;
+        z /= n;
+    }
+
+    /**
+     * Multiply this quaternion by the param quaternion
+     *
+     * @param q a quaternion to multiply with
+     */
+    public void mult(Quaternion q) {
+        final float w1 = w * q.w - x * q.x - y * q.y - z * q.z;
+
+        final float x1 = w * q.x + x * q.w + y * q.z - z * q.y;
+        final float y1 = w * q.y - x * q.z + y * q.w + z * q.x;
+        final float z1 = w * q.z + x * q.y - y * q.x + z * q.w;
+
+        w = w1;
+        x = x1;
+        y = y1;
+        z = z1;
+    }
+
+    /**
+     * Multiply a quaternion by a constant
+     *
+     * @param n a float constant
+     */
+    public void mult(float n) {
+        x *= n;
+        y *= n;
+        z *= n;
+    }
+
+    /***
+     * Rotate given vector by this quaternion
+     *
+     * @param vector input vector
+     * @return rotated vector
+     */
+    public float[] mult(float[] vector) {
+        // TODO : optimize
+        final float[] res = new float[3];
+        final Quaternion a = new Quaternion(vector[0], vector[1], vector[2], 0.0f);
+        final Quaternion b = new Quaternion(this);
+        final Quaternion c = new Quaternion(this);
+        b.inverse();
+        a.mult(b);
+        c.mult(a);
+        res[0] = c.x;
+        res[1] = c.y;
+        res[2] = c.z;
+        return res;
+    }
+
+    /**
+     * Normalize a quaternion required if to be used as a rotational quaternion
+     */
+    public void normalize() {
+        final float norme = (float) FloatUtil.sqrt(w * w + x * x + y * y + z * z);
+        if (norme == 0.0f) {
+            setIdentity();
+        } else {
+            final float recip = 1.0f / norme;
+
+            w *= recip;
+            x *= recip;
+            y *= recip;
+            z *= recip;
+        }
+    }
+
+    /**
+     * Invert the quaternion If rotational, will produce a the inverse rotation
+     */
+    public void inverse() {
+        final float norm = w * w + x * x + y * y + z * z;
+
+        final float recip = 1.0f / norm;
+
+        w *= recip;
+        x = -1 * x * recip;
+        y = -1 * y * recip;
+        z = -1 * z * recip;
+    }
+
+    /**
+     * Transform this quaternion to a 4x4 column matrix representing the
+     * rotation
+     *
+     * @return new float[16] column matrix 4x4
+     */
+    public float[] toMatrix() {
+        final float[] matrix = new float[16];
+        matrix[0] = 1.0f - 2 * y * y - 2 * z * z;
+        matrix[1] = 2 * x * y + 2 * w * z;
+        matrix[2] = 2 * x * z - 2 * w * y;
+        matrix[3] = 0;
+
+        matrix[4] = 2 * x * y - 2 * w * z;
+        matrix[5] = 1.0f - 2 * x * x - 2 * z * z;
+        matrix[6] = 2 * y * z + 2 * w * x;
+        matrix[7] = 0;
+
+        matrix[8] = 2 * x * z + 2 * w * y;
+        matrix[9] = 2 * y * z - 2 * w * x;
+        matrix[10] = 1.0f - 2 * x * x - 2 * y * y;
+        matrix[11] = 0;
+
+        matrix[12] = 0;
+        matrix[13] = 0;
+        matrix[14] = 0;
+        matrix[15] = 1;
+        return matrix;
+    }
+
+    /**
+     * Set this quaternion from a Sphereical interpolation of two param
+     * quaternion, used mostly for rotational animation.
+     * <p>
+     * Note: Method does not normalize this quaternion!
+     * </p>
+     * <p>
+     * See http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/
+     * quaternions/slerp/
+     * </p>
+     *
+     * @param a initial quaternion
+     * @param b target quaternion
+     * @param t float between 0 and 1 representing interp.
+     */
+    public void slerp(Quaternion a, Quaternion b, float t) {
+        final float cosom = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
+        final float t1 = 1.0f - t;
+
+        // if the two quaternions are close, just use linear interpolation
+        if (cosom >= 0.95f) {
+            x = a.x * t1 + b.x * t;
+            y = a.y * t1 + b.y * t;
+            z = a.z * t1 + b.z * t;
+            w = a.w * t1 + b.w * t;
+            return;
+        }
+
+        // the quaternions are nearly opposite, we can pick any axis normal to
+        // a,b
+        // to do the rotation
+        if (cosom <= -0.99f) {
+            x = 0.5f * (a.x + b.x);
+            y = 0.5f * (a.y + b.y);
+            z = 0.5f * (a.z + b.z);
+            w = 0.5f * (a.w + b.w);
+            return;
+        }
+
+        // cosom is now withion range of acos, do a SLERP
+        final float sinom = FloatUtil.sqrt(1.0f - cosom * cosom);
+        final float omega = FloatUtil.acos(cosom);
+
+        final float scla = FloatUtil.sin(t1 * omega) / sinom;
+        final float sclb = FloatUtil.sin(t * omega) / sinom;
+
+        x = a.x * scla + b.x * sclb;
+        y = a.y * scla + b.y * sclb;
+        z = a.z * scla + b.z * sclb;
+        w = a.w * scla + b.w * sclb;
+    }
+
+    /**
+     * Check if this quaternion represents an identity matrix for rotation,
+     * , ie (0,0,0,1).
+     *
+     * @return true if it is an identity rep., false otherwise
+     */
+    public boolean isIdentity() {
+        return w == 1 && x == 0 && y == 0 && z == 0;
+    }
+
+    /***
+     * Set this quaternion to identity (x=0,y=0,z=0,w=1)
+     */
+    public void setIdentity() {
+        x = y = z = 0;
+        w = 1;
+    }
+
+    /**
+     * compute the quaternion from a 3x3 column matrix
+     *
+     * @param m 3x3 column matrix
+     */
+    public void setFromMatrix(float[] m) {
+        final float T = m[0] + m[4] + m[8] + 1;
+        if (T > 0) {
+            final float S = 0.5f / (float) FloatUtil.sqrt(T);
+            w = 0.25f / S;
+            x = (m[5] - m[7]) * S;
+            y = (m[6] - m[2]) * S;
+            z = (m[1] - m[3]) * S;
+        } else {
+            if ((m[0] > m[4]) && (m[0] > m[8])) {
+                final float S = FloatUtil.sqrt(1.0f + m[0] - m[4] - m[8]) * 2f; // S=4*qx
+                w = (m[7] - m[5]) / S;
+                x = 0.25f * S;
+                y = (m[3] + m[1]) / S;
+                z = (m[6] + m[2]) / S;
+            } else if (m[4] > m[8]) {
+                final float S = FloatUtil.sqrt(1.0f + m[4] - m[0] - m[8]) * 2f; // S=4*qy
+                w = (m[6] - m[2]) / S;
+                x = (m[3] + m[1]) / S;
+                y = 0.25f * S;
+                z = (m[7] + m[5]) / S;
+            } else {
+                final float S = FloatUtil.sqrt(1.0f + m[8] - m[0] - m[4]) * 2f; // S=4*qz
+                w = (m[3] - m[1]) / S;
+                x = (m[6] + m[2]) / S;
+                y = (m[7] + m[5]) / S;
+                z = 0.25f * S;
+            }
+        }
+    }
+
+    /**
+     * Check if the the 3x3 matrix (param) is in fact an affine rotational
+     * matrix
+     *
+     * @param m 3x3 column matrix
+     * @return true if representing a rotational matrix, false otherwise
+     */
+    public boolean isRotationMatrix(float[] m) {
+        final float epsilon = 0.01f; // margin to allow for rounding errors
+        if (FloatUtil.abs(m[0] * m[3] + m[3] * m[4] + m[6] * m[7]) > epsilon)
+            return false;
+        if (FloatUtil.abs(m[0] * m[2] + m[3] * m[5] + m[6] * m[8]) > epsilon)
+            return false;
+        if (FloatUtil.abs(m[1] * m[2] + m[4] * m[5] + m[7] * m[8]) > epsilon)
+            return false;
+        if (FloatUtil.abs(m[0] * m[0] + m[3] * m[3] + m[6] * m[6] - 1) > epsilon)
+            return false;
+        if (FloatUtil.abs(m[1] * m[1] + m[4] * m[4] + m[7] * m[7] - 1) > epsilon)
+            return false;
+        if (FloatUtil.abs(m[2] * m[2] + m[5] * m[5] + m[8] * m[8] - 1) > epsilon)
+            return false;
+        return (FloatUtil.abs(determinant(m) - 1) < epsilon);
+    }
+
+    private float determinant(float[] m) {
+        return m[0] * m[4] * m[8] + m[3] * m[7] * m[2] + m[6] * m[1] * m[5]
+             - m[0] * m[7] * m[5] - m[3] * m[1] * m[8] - m[6] * m[4] * m[2];
+    }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
new file mode 100644
index 000000000..053876b56
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
@@ -0,0 +1,565 @@
+/**
+ * Copyright 2010 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math;
+
+import java.util.ArrayList;
+
+public class VectorUtil {
+
+    public enum Winding {
+        CW(-1), CCW(1);
+
+        public final int dir;
+
+        Winding(int dir) {
+            this.dir = dir;
+        }
+    }
+
+    public static final int COLLINEAR = 0;
+
+    /** 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)
+    {
+        return (vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2]);
+    }
+    /** Normalize a vector
+     * @param vector input vector
+     * @return normalized vector
+     */
+    public static float[] normalize(float[] vector)
+    {
+        final float[] newVector = new float[3];
+
+        final float d = FloatUtil.sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
+        if(d> 0.0f)
+        {
+            newVector[0] = vector[0]/d;
+            newVector[1] = vector[1]/d;
+            newVector[2] = vector[2]/d;
+        }
+        return newVector;
+    }
+
+    /** Scales a vector by param creating a new float[] for the result!
+     * @param vector input vector
+     * @param scale constant to scale by
+     * @return new scaled vector
+     * @deprecated Use {@link #scale(float[], float[], float)}
+     */
+    public static float[] scale(float[] vector, float scale)
+    {
+        final float[] newVector = new float[3];
+
+        newVector[0] = vector[0] * scale;
+        newVector[1] = vector[1] * scale;
+        newVector[2] = vector[2] * scale;
+        return newVector;
+    }
+
+    /** 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)
+    {
+        result[0] = vector[0] * scale;
+        result[1] = vector[1] * scale;
+        result[2] = vector[2] * scale;
+        return result;
+    }
+
+    /** Scales a vector by param using given result float[]
+     * @param result vector for the result
+     * @param vector input vector
+     * @param scale 3 component scale constant for each vector component
+     * @return given result vector
+     */
+    public static float[] scale(float[] result, float[] vector, float[] scale)
+    {
+        result[0] = vector[0] * scale[0];
+        result[1] = vector[1] * scale[1];
+        result[2] = vector[2] * scale[2];
+        return result;
+    }
+
+    /** Adds to vectors
+     * @param v1 vector 1
+     * @param v2 vector 2
+     * @return v1 + v2
+     */
+    public static float[] vectorAdd(float[] v1, float[] v2)
+    {
+        final float[] newVector = new float[3];
+
+        newVector[0] = v1[0] + v2[0];
+        newVector[1] = v1[1] + v2[1];
+        newVector[2] = v1[2] + v2[2];
+        return newVector;
+    }
+
+    /** cross product vec1 x vec2
+     * @param vec1 vector 1
+     * @param vec2 vecttor 2
+     * @return the resulting vector
+     */
+    public static float[] cross(float[] vec1, float[] vec2)
+    {
+        final float[] out = new float[3];
+
+        out[0] = vec2[2]*vec1[1] - vec2[1]*vec1[2];
+        out[1] = vec2[0]*vec1[2] - vec2[2]*vec1[0];
+        out[2] = vec2[1]*vec1[0] - vec2[0]*vec1[1];
+
+        return out;
+    }
+
+    /** Column Matrix Vector multiplication
+     * @param colMatrix column matrix (4x4)
+     * @param vec vector(x,y,z)
+     * @return result new float[3]
+     */
+    public static float[] colMatrixVectorMult(float[] colMatrix, float[] vec)
+    {
+        final float[] out = new float[3];
+
+        out[0] = vec[0]*colMatrix[0] + vec[1]*colMatrix[4] + vec[2]*colMatrix[8] + colMatrix[12];
+        out[1] = vec[0]*colMatrix[1] + vec[1]*colMatrix[5] + vec[2]*colMatrix[9] + colMatrix[13];
+        out[2] = vec[0]*colMatrix[2] + vec[1]*colMatrix[6] + vec[2]*colMatrix[10] + colMatrix[14];
+
+        return out;
+    }
+
+    /** Matrix Vector multiplication
+     * @param rawMatrix column matrix (4x4)
+     * @param vec vector(x,y,z)
+     * @return result new float[3]
+     */
+    public static float[] rowMatrixVectorMult(float[] rawMatrix, float[] vec)
+    {
+        final float[] out = new float[3];
+
+        out[0] = vec[0]*rawMatrix[0] + vec[1]*rawMatrix[1] + vec[2]*rawMatrix[2] + rawMatrix[3];
+        out[1] = vec[0]*rawMatrix[4] + vec[1]*rawMatrix[5] + vec[2]*rawMatrix[6] + rawMatrix[7];
+        out[2] = vec[0]*rawMatrix[8] + vec[1]*rawMatrix[9] + vec[2]*rawMatrix[10] + rawMatrix[11];
+
+        return out;
+    }
+
+    /** Calculate the midpoint of two values
+     * @param p1 first value
+     * @param p2 second vale
+     * @return midpoint
+     */
+    public static float mid(float p1, float p2)
+    {
+        return (p1+p2)/2.0f;
+    }
+
+    /** Calculate the midpoint of two points
+     * @param p1 first point
+     * @param p2 second point
+     * @return midpoint
+     */
+    public static float[] mid(float[] p1, float[] p2)
+    {
+        final float[] midPoint = new float[3];
+        midPoint[0] = (p1[0] + p2[0])*0.5f;
+        midPoint[1] = (p1[1] + p2[1])*0.5f;
+        midPoint[2] = (p1[2] + p2[2])*0.5f;
+
+        return midPoint;
+    }
+
+    /** 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)
+    {
+        final float w0 = point[0]-p0[0];
+        final float w1 = point[1]-p0[1];
+        final float w2 = point[2]-p0[2];
+
+        return FloatUtil.sqrt(w0*w0 + w1*w1 + w2*w2);
+    }
+
+    /**Check equality of 2 vec3 vectors
+     * @param v1 vertex 1
+     * @param v2 vertex 2
+     * @return
+     */
+    public static boolean checkEquality(float[] v1, float[] v2)
+    {
+        return Float.compare(v1[0], v2[0]) == 0 &&
+               Float.compare(v1[1], v2[1]) == 0 &&
+               Float.compare(v1[2], v2[2]) == 0 ;
+    }
+
+    /**Check equality of 2 vec2 vectors
+     * @param v1 vertex 1
+     * @param v2 vertex 2
+     * @return
+     */
+    public static boolean checkEqualityVec2(float[] v1, float[] v2)
+    {
+        return Float.compare(v1[0], v2[0]) == 0 &&
+               Float.compare(v1[1], v2[1]) == 0 ;
+    }
+
+    /** Compute the determinant of 3 vectors
+     * @param a vector 1
+     * @param b vector 2
+     * @param c vector 3
+     * @return the determinant value
+     */
+    public static float computeDeterminant(float[] a, float[] b, 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];
+    }
+
+    /** Check if three vertices are colliniear
+     * @param v1 vertex 1
+     * @param v2 vertex 2
+     * @param v3 vertex 3
+     * @return true if collinear, false otherwise
+     */
+    public static boolean checkCollinear(float[] v1, float[] v2, float[] v3)
+    {
+        return (computeDeterminant(v1, v2, v3) == VectorUtil.COLLINEAR);
+    }
+
+    /** Compute Vector
+     * @param vector storage for resulting Vector V1V2
+     * @param v1 vertex 1
+     * @param v2 vertex2 2
+     */
+    public static void computeVector(float[] vector, float[] v1, float[] v2) {
+        vector[0] = v2[0] - v1[0];
+        vector[1] = v2[1] - v1[1];
+        vector[2] = v2[2] - v1[2];
+    }
+
+    /** Check if vertices in triangle circumcircle
+     * @param a triangle vertex 1
+     * @param b triangle vertex 2
+     * @param c triangle vertex 3
+     * @param d vertex in question
+     * @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) {
+        final float[] A = a.getCoord();
+        final float[] B = b.getCoord();
+        final float[] C = c.getCoord();
+        final float[] D = d.getCoord();
+        return (A[0] * A[0] + A[1] * A[1]) * triArea(B, C, D) -
+               (B[0] * B[0] + B[1] * B[1]) * triArea(A, C, D) +
+               (C[0] * C[0] + C[1] * C[1]) * triArea(A, B, D) -
+               (D[0] * D[0] + D[1] * D[1]) * triArea(A, B, C) > 0;
+    }
+
+    /** Computes oriented area of a triangle
+     * @param a first vertex
+     * @param b second vertex
+     * @param c third vertex
+     * @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){
+        final float[] A = a.getCoord();
+        final float[] B = b.getCoord();
+        final float[] C = c.getCoord();
+        return (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
+    }
+
+    /** Computes oriented area of a triangle
+     * @param A first vertex
+     * @param B second vertex
+     * @param C third vertex
+     * @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){
+        return (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1])*(C[0] - A[0]);
+    }
+
+    /** Check if a vertex is in triangle using
+     * barycentric coordinates computation.
+     * @param a first triangle vertex
+     * @param b second triangle vertex
+     * @param c third triangle vertex
+     * @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){
+        // Compute vectors
+        computeVector(ac, a, c); //v0
+        computeVector(ab, a, b); //v1
+        computeVector(ap, a, p); //v2
+
+        // Compute dot products
+        final float dot00 = dot(ac, ac);
+        final float dot01 = dot(ac, ab);
+        final float dot02 = dot(ac, ap);
+        final float dot11 = dot(ab, ab);
+        final float dot12 = dot(ab, ap);
+
+        // Compute barycentric coordinates
+        final float invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
+        final float u = (dot11 * dot02 - dot01 * dot12) * invDenom;
+        final float v = (dot00 * dot12 - dot01 * dot02) * invDenom;
+
+        // Check if point is in triangle
+        return (u >= 0) && (v >= 0) && (u + v < 1);
+    }
+
+    /** Check if one of three vertices are in triangle using
+     * barycentric coordinates computation.
+     * @param a first triangle vertex
+     * @param b second triangle vertex
+     * @param c third triangle vertex
+     * @param p1 the vertex in question
+     * @param p2 the vertex in question
+     * @param p3 the vertex in question
+     * @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[] ac, float[] ab, float[] ap){
+        // Compute vectors
+        computeVector(ac, a, c); //v0
+        computeVector(ab, a, b); //v1
+
+        // Compute dot products
+        final float dotAC_AC = dot(ac, ac);
+        final float dotAC_AB = dot(ac, ab);
+        final float dotAB_AB = dot(ab, ab);
+
+        // Compute barycentric coordinates
+        final float invDenom = 1 / (dotAC_AC * dotAB_AB - dotAC_AB * dotAC_AB);
+        {
+            computeVector(ap, a, p1); //v2
+            final float dotAC_AP1 = dot(ac, ap);
+            final float dotAB_AP1 = dot(ab, ap);
+            final float u1 = (dotAB_AB * dotAC_AP1 - dotAC_AB * dotAB_AP1) * invDenom;
+            final float v1 = (dotAC_AC * dotAB_AP1 - dotAC_AB * dotAC_AP1) * invDenom;
+
+            // Check if point is in triangle
+            if ( (u1 >= 0) && (v1 >= 0) && (u1 + v1 < 1) ) {
+                return true;
+            }
+        }
+
+        {
+            computeVector(ap, a, p2); //v2
+            final float dotAC_AP2 = dot(ac, ap);
+            final float dotAB_AP2 = dot(ab, ap);
+            final float u = (dotAB_AB * dotAC_AP2 - dotAC_AB * dotAB_AP2) * invDenom;
+            final float v = (dotAC_AC * dotAB_AP2 - dotAC_AB * dotAC_AP2) * invDenom;
+
+            // Check if point is in triangle
+            if ( (u >= 0) && (v >= 0) && (u + v < 1) ) {
+                return true;
+            }
+        }
+
+        {
+            computeVector(ap, a, p3); //v2
+            final float dotAC_AP3 = dot(ac, ap);
+            final float dotAB_AP3 = dot(ab, ap);
+            final float u = (dotAB_AB * dotAC_AP3 - dotAC_AB * dotAB_AP3) * invDenom;
+            final float v = (dotAC_AC * dotAB_AP3 - dotAC_AB * dotAC_AP3) * invDenom;
+
+            // Check if point is in triangle
+            if ( (u >= 0) && (v >= 0) && (u + v < 1) ) {
+                return true;
+            }
+        }
+
+        return false;
+    }
+
+    /** Check if points are in ccw order
+     * @param a first vertex
+     * @param b second vertex
+     * @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){
+        return triArea(a,b,c) > 0;
+    }
+
+    /** Compute the winding of given points
+     * @param a first vertex
+     * @param b second vertex
+     * @param c third vertex
+     * @return Winding
+     */
+    public static Winding getWinding(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c) {
+        return triArea(a,b,c) > 0 ? Winding.CCW : Winding.CW ;
+    }
+
+    /** Computes the area of a list of vertices to check if ccw
+     * @param vertices
+     * @return positive area if ccw else negative area value
+     */
+    public static float area(ArrayList<? extends Vert2fImmutable> vertices) {
+        final int n = vertices.size();
+        float area = 0.0f;
+        for (int p = n - 1, q = 0; q < n; p = q++)
+        {
+            final float[] pCoord = vertices.get(p).getCoord();
+            final float[] qCoord = vertices.get(q).getCoord();
+            area += pCoord[0] * qCoord[1] - qCoord[0] * pCoord[1];
+        }
+        return area;
+    }
+
+    /** Compute the  general winding of the vertices
+     * @param vertices array of Vertices
+     * @return CCW or CW {@link Winding}
+     */
+    public static Winding getWinding(ArrayList<? extends Vert2fImmutable> vertices) {
+        return area(vertices) >= 0 ? Winding.CCW : Winding.CW ;
+    }
+
+
+    /** Compute intersection between two segments
+     * @param a vertex 1 of first segment
+     * @param b vertex 2 of first segment
+     * @param c vertex 1 of second segment
+     * @param d vertex 2 of second segment
+     * @return the intersection coordinates if the segments intersect, otherwise
+     * returns null
+     */
+    public static float[] seg2SegIntersection(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+        final float determinant = (a.getX()-b.getX())*(c.getY()-d.getY()) - (a.getY()-b.getY())*(c.getX()-d.getX());
+
+        if (determinant == 0)
+            return null;
+
+        final float alpha = (a.getX()*b.getY()-a.getY()*b.getX());
+        final float beta = (c.getX()*d.getY()-c.getY()*d.getY());
+        final float xi = ((c.getX()-d.getX())*alpha-(a.getX()-b.getX())*beta)/determinant;
+        final float yi = ((c.getY()-d.getY())*alpha-(a.getY()-b.getY())*beta)/determinant;
+
+        final float gamma = (xi - a.getX())/(b.getX() - a.getX());
+        final float gamma1 = (xi - c.getX())/(d.getX() - c.getX());
+        if(gamma <= 0 || gamma >= 1) return null;
+        if(gamma1 <= 0 || gamma1 >= 1) return null;
+
+        return new float[]{xi,yi,0};
+    }
+
+    /** Compute intersection between two segments
+     * @param a vertex 1 of first segment
+     * @param b vertex 2 of first segment
+     * @param c vertex 1 of second segment
+     * @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) {
+        final float[] A = a.getCoord();
+        final float[] B = b.getCoord();
+        final float[] C = c.getCoord();
+        final float[] D = d.getCoord();
+
+        final float determinant = (A[0]-B[0])*(C[1]-D[1]) - (A[1]-B[1])*(C[0]-D[0]);
+
+        if (determinant == 0) {
+            return false;
+        }
+
+        final float alpha = (A[0]*B[1]-A[1]*B[0]);
+        final float beta = (C[0]*D[1]-C[1]*D[1]);
+        final float xi = ((C[0]-D[0])*alpha-(A[0]-B[0])*beta)/determinant;
+
+        final float gamma = (xi - A[0])/(B[0] - A[0]);
+        final float gamma1 = (xi - C[0])/(D[0] - C[0]);
+        if(gamma <= 0 || gamma >= 1 || gamma1 <= 0 || gamma1 >= 1) {
+            return false;
+        }
+
+        return true;
+    }
+
+    /** Compute intersection between two lines
+     * @param a vertex 1 of first line
+     * @param b vertex 2 of first line
+     * @param c vertex 1 of second line
+     * @param d vertex 2 of second line
+     * @return the intersection coordinates if the lines intersect, otherwise
+     * returns null
+     */
+    public static float[] line2lineIntersection(Vert2fImmutable a, Vert2fImmutable b, Vert2fImmutable c, Vert2fImmutable d) {
+        final float determinant = (a.getX()-b.getX())*(c.getY()-d.getY()) - (a.getY()-b.getY())*(c.getX()-d.getX());
+
+        if (determinant == 0)
+            return null;
+
+        final float alpha = (a.getX()*b.getY()-a.getY()*b.getX());
+        final float beta = (c.getX()*d.getY()-c.getY()*d.getY());
+        final float xi = ((c.getX()-d.getX())*alpha-(a.getX()-b.getX())*beta)/determinant;
+        final float yi = ((c.getY()-d.getY())*alpha-(a.getY()-b.getY())*beta)/determinant;
+
+        return new float[]{xi,yi,0};
+    }
+
+    /** Check if a segment intersects with a triangle
+     * @param a vertex 1 of the triangle
+     * @param b vertex 2 of the triangle
+     * @param c vertex 3 of the triangle
+     * @param d vertex 1 of first segment
+     * @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){
+        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/Vert2fImmutable.java b/src/jogl/classes/com/jogamp/opengl/math/Vert2fImmutable.java
new file mode 100644
index 000000000..ec90b401f
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Vert2fImmutable.java
@@ -0,0 +1,39 @@
+/**
+ * Copyright 2012 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math;
+
+public interface Vert2fImmutable {
+    float getX();
+
+    float getY();
+
+    int getCoordCount();
+
+    float[] getCoord();
+
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Vert3fImmutable.java b/src/jogl/classes/com/jogamp/opengl/math/Vert3fImmutable.java
new file mode 100644
index 000000000..76bd02fbc
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Vert3fImmutable.java
@@ -0,0 +1,32 @@
+/**
+ * Copyright 2012 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math;
+
+public interface Vert3fImmutable extends Vert2fImmutable {
+    float getZ();
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java b/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java
new file mode 100644
index 000000000..5fbc28c60
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java
@@ -0,0 +1,360 @@
+/**
+ * Copyright 2010 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math.geom;
+
+import com.jogamp.opengl.math.VectorUtil;
+
+
+/**
+ * Axis Aligned Bounding Box. Defined by two 3D coordinates (low and high)
+ * The low being the the lower left corner of the box, and the high being the upper
+ * right corner of the box.
+ *
+ */
+public class AABBox implements Cloneable {
+    private float[] low = new float[3];
+    private float[] high = new float[3];
+    private float[] center = new float[3];
+
+    /** Create a Axis Aligned bounding box (AABBox)
+     * where the low and and high MAX float Values.
+     */
+    public AABBox() {
+        reset();
+    }
+
+    /** Create an AABBox specifying the coordinates
+     * of the low and high
+     * @param lx min x-coordinate
+     * @param ly min y-coordnate
+     * @param lz min z-coordinate
+     * @param hx max x-coordinate
+     * @param hy max y-coordinate
+     * @param hz max z-coordinate
+     */
+    public AABBox(float lx, float ly, float lz,
+                  float hx, float hy, float hz) {
+        setSize(lx, ly, lz, hx, hy, hz);
+    }
+
+    /** Create a AABBox defining the low and high
+     * @param low min xyz-coordinates
+     * @param high max xyz-coordinates
+     */
+    public AABBox(float[] low, float[] high) {
+        setSize(low[0],low[1],low[2], high[0],high[1],high[2]);
+    }
+
+    /** resets this box to the inverse low/high, allowing the next {@link #resize(float, float, float)} command to hit. */
+    public final void reset() {
+        setLow(Float.MAX_VALUE,Float.MAX_VALUE,Float.MAX_VALUE);
+        setHigh(-1*Float.MAX_VALUE,-1*Float.MAX_VALUE,-1*Float.MAX_VALUE);
+        center[0] = 0f;
+        center[1] = 0f;
+        center[2] = 0f;
+    }
+
+    /** Get the max xyz-coordinates
+     * @return a float array containing the max xyz coordinates
+     */
+    public final float[] getHigh() {
+        return high;
+    }
+
+    private final void setHigh(float hx, float hy, float hz) {
+        this.high[0] = hx;
+        this.high[1] = hy;
+        this.high[2] = hz;
+    }
+
+    /** Get the min xyz-coordinates
+     * @return a float array containing the min xyz coordinates
+     */
+    public final float[] getLow() {
+        return low;
+    }
+
+    private final void setLow(float lx, float ly, float lz) {
+        this.low[0] = lx;
+        this.low[1] = ly;
+        this.low[2] = lz;
+    }
+
+    private final void computeCenter() {
+        center[0] = (high[0] + low[0])/2;
+        center[1] = (high[1] + low[1])/2;
+        center[2] = (high[2] + low[2])/2;
+    }
+
+    /**
+     * Set size of the AABBox specifying the coordinates
+     * of the low and high.
+     *
+     * @param lx min x-coordinate
+     * @param ly min y-coordnate
+     * @param lz min z-coordinate
+     * @param hx max x-coordinate
+     * @param hy max y-coordinate
+     * @param hz max z-coordinate
+     */
+    public final void setSize(float lx, float ly, float lz,
+                              float hx, float hy, float hz) {
+        this.low[0] = lx;
+        this.low[1] = ly;
+        this.low[2] = lz;
+        this.high[0] = hx;
+        this.high[1] = hy;
+        this.high[2] = hz;
+        computeCenter();
+    }
+
+    /** Resize the AABBox to encapsulate another AABox
+     * @param newBox AABBox to be encapsulated in
+     */
+    public final void resize(AABBox newBox) {
+        float[] newLow = newBox.getLow();
+        float[] newHigh = newBox.getHigh();
+
+        /** test low */
+        if (newLow[0] < low[0])
+            low[0] = newLow[0];
+        if (newLow[1] < low[1])
+            low[1] = newLow[1];
+        if (newLow[2] < low[2])
+            low[2] = newLow[2];
+
+        /** test high */
+        if (newHigh[0] > high[0])
+            high[0] = newHigh[0];
+        if (newHigh[1] > high[1])
+            high[1] = newHigh[1];
+        if (newHigh[2] > high[2])
+            high[2] = newHigh[2];
+
+        computeCenter();
+    }
+
+    /** Resize the AABBox to encapsulate the passed
+     * xyz-coordinates.
+     * @param x x-axis coordinate value
+     * @param y y-axis coordinate value
+     * @param z z-axis coordinate value
+     */
+    public final void resize(float x, float y, float z) {
+        /** test low */
+        if (x < low[0])
+            low[0] = x;
+        if (y < low[1])
+            low[1] = y;
+        if (z < low[2])
+            low[2] = z;
+
+        /** test high */
+        if (x > high[0])
+            high[0] = x;
+        if (y > high[1])
+            high[1] = y;
+        if (z > high[2])
+            high[2] = z;
+
+        computeCenter();
+    }
+
+    /** Resize the AABBox to encapsulate the passed
+     * xyz-coordinates.
+     * @param xyz xyz-axis coordinate values
+     * @param offset of the array
+     */
+    public final void resize(float[] xyz, int offset) {
+        resize(xyz[0+offset], xyz[1+offset], xyz[2+offset]);
+    }
+
+    /** Check if the x & y coordinates are bounded/contained
+     *  by this AABBox
+     * @param x  x-axis coordinate value
+     * @param y  y-axis coordinate value
+     * @return true if  x belong to (low.x, high.x) and
+     * y belong to (low.y, high.y)
+     */
+    public final boolean contains(float x, float y) {
+        if(x<low[0] || x>high[0]){
+            return false;
+        }
+        if(y<low[1]|| y>high[1]){
+            return false;
+        }
+        return true;
+    }
+
+    /** Check if the xyz coordinates are bounded/contained
+     *  by this AABBox.
+     * @param x x-axis coordinate value
+     * @param y y-axis coordinate value
+     * @param z z-axis coordinate value
+     * @return true if  x belong to (low.x, high.x) and
+     * y belong to (low.y, high.y) and  z belong to (low.z, high.z)
+     */
+    public final boolean contains(float x, float y, float z) {
+        if(x<low[0] || x>high[0]){
+            return false;
+        }
+        if(y<low[1]|| y>high[1]){
+            return false;
+        }
+        if(z<low[2] || z>high[2]){
+            return false;
+        }
+        return true;
+    }
+
+    /** Check if there is a common region between this AABBox and the passed
+     *     2D region irrespective of z range
+     * @param x lower left x-coord
+     * @param y lower left y-coord
+     * @param w width
+     * @param h hight
+     * @return true if this AABBox might have a common region with this 2D region
+     */
+    public final boolean intersects(float x, float y, float w, float h) {
+        if (w <= 0 || h <= 0) {
+            return false;
+        }
+
+        final float _w = getWidth();
+        final float _h = getHeight();
+        if (_w <= 0 || _h <= 0) {
+            return false;
+        }
+
+        final float x0 = getMinX();
+        final float y0 = getMinY();
+        return (x + w > x0 &&
+                y + h > y0 &&
+                x < x0 + _w &&
+                y < y0 + _h);
+    }
+
+
+    /** Get the size of the Box where the size is represented by the
+     * length of the vector between low and high.
+     * @return a float representing the size of the AABBox
+     */
+    public final float getSize() {
+        return VectorUtil.computeLength(low, high);
+    }
+
+    /**Get the Center of the AABBox
+     * @return the xyz-coordinates of the center of the AABBox
+     */
+    public final float[] getCenter() {
+        return center;
+    }
+
+    /** Scale the AABBox by a constant
+     * @param size a constant float value
+     */
+    public final void scale(float size) {
+        float[] diffH = new float[3];
+        diffH[0] = high[0] - center[0];
+        diffH[1] = high[1] - center[1];
+        diffH[2] = high[2] - center[2];
+
+        diffH = VectorUtil.scale(diffH, size);
+
+        float[] diffL = new float[3];
+        diffL[0] = low[0] - center[0];
+        diffL[1] = low[1] - center[1];
+        diffL[2] = low[2] - center[2];
+
+        diffL = VectorUtil.scale(diffL, size);
+
+        high = VectorUtil.vectorAdd(center, diffH);
+        low = VectorUtil.vectorAdd(center, diffL);
+    }
+
+    public final float getMinX() {
+        return low[0];
+    }
+
+    public final float getMinY() {
+        return low[1];
+    }
+
+    public final float getMinZ() {
+        return low[2];
+    }
+
+    public final float getMaxX() {
+        return high[0];
+    }
+
+    public final float getMaxY() {
+        return high[1];
+    }
+
+    public final float getMaxZ() {
+        return high[2];
+    }
+
+    public final float getWidth(){
+        return high[0] - low[0];
+    }
+
+    public final float getHeight() {
+        return high[1] - low[1];
+    }
+
+    public final float getDepth() {
+        return high[2] - low[2];
+    }
+
+    @Override
+    public final AABBox clone() {
+        return new AABBox(this.low, this.high);
+    }
+
+    @Override
+    public final boolean equals(Object obj) {
+        if( obj == this ) {
+            return true;
+        }
+        if( null == obj || !(obj instanceof AABBox) ) {
+            return false;
+        }
+        final AABBox other = (AABBox) obj;
+        return VectorUtil.checkEquality(low, other.low) &&
+               VectorUtil.checkEquality(high, other.high) ;
+    }
+
+    @Override
+    public final String toString() {
+        return "[ dim "+getWidth()+" x "+getHeight()+" x "+getDepth()+
+               ", box "+low[0]+" / "+low[1]+" / "+low[2]+" .. "+high[0]+" / "+high[1]+" / "+high[2]+
+               ", ctr "+center[0]+" / "+center[1]+" / "+center[2]+" ]";
+    }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
new file mode 100644
index 000000000..fb311083f
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
@@ -0,0 +1,388 @@
+/**
+ * Copyright 2010 JogAmp Community. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without modification, are
+ * permitted provided that the following conditions are met:
+ *
+ *    1. Redistributions of source code must retain the above copyright notice, this list of
+ *       conditions and the following disclaimer.
+ *
+ *    2. Redistributions in binary form must reproduce the above copyright notice, this list
+ *       of conditions and the following disclaimer in the documentation and/or other materials
+ *       provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
+ * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+ * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * The views and conclusions contained in the software and documentation are those of the
+ * authors and should not be interpreted as representing official policies, either expressed
+ * or implied, of JogAmp Community.
+ */
+package com.jogamp.opengl.math.geom;
+
+import com.jogamp.common.os.Platform;
+
+/**
+ * Providing frustum {@link #getPlanes() planes} derived by different inputs
+ * ({@link #updateByPMV(float[], int) P*MV}, ..)
+ * used to {@link #classifySphere(float[], float) classify objects} and to test
+ * whether they are {@link #isOutside(AABBox) outside}.
+ *
+ * <p>
+ * Extracting the world-frustum planes from the P*Mv:
+ * <pre>
+ * Fast Extraction of Viewing Frustum Planes from the World-View-Projection Matrix
+ *   Gil Gribb <[email protected]>
+ *   Klaus Hartmann <[email protected]>
+ *   http://graphics.cs.ucf.edu/cap4720/fall2008/plane_extraction.pdf
+ * </pre>
+ * Classifying Point, Sphere and AABBox:
+ * <pre>
+ * Efficient View Frustum Culling
+ *   Daniel Sýkora <[email protected]>
+ *   Josef Jelínek <[email protected]>
+ *   http://www.cg.tuwien.ac.at/hostings/cescg/CESCG-2002/DSykoraJJelinek/index.html
+ * </pre>
+ * <pre>
+ * Lighthouse3d.com
+ * http://www.lighthouse3d.com/tutorials/view-frustum-culling/
+ * </pre>
+ *
+ * Fundamentals about Planes, Half-Spaces and Frustum-Culling:<br/>
+ * <pre>
+ * Planes and Half-Spaces,  Max Wagner <[email protected]>
+ * http://www.emeyex.com/site/tuts/PlanesHalfSpaces.pdf
+ * </pre>
+ * <pre>
+ * Frustum Culling,  Max Wagner <[email protected]>
+ * http://www.emeyex.com/site/tuts/FrustumCulling.pdf
+ * </pre>
+ * </p>
+ */
+public class Frustum {
+    /** Normalized planes[l, r, b, t, n, f] */
+	protected Plane[] planes = new Plane[6];
+
+	/**
+	 * Creates an undefined instance w/o calculating the frustum.
+	 * <p>
+	 * Use one of the <code>update(..)</code> methods to set the {@link #getPlanes() planes}.
+	 * </p>
+	 * @see #updateByPlanes(Plane[])
+	 * @see #updateByPMV(float[], int)
+	 */
+    public Frustum() {
+        for (int i = 0; i < 6; ++i) {
+            planes[i] = new Plane();
+        }
+    }
+
+	/**
+	 * Plane equation := dot(n, x - p) = 0 ->  ax + bc + cx + d == 0
+	 * <p>
+	 * In order to work w/ {@link Frustum#isOutside(AABBox) isOutside(..)} methods,
+	 * the normals have to point to the inside of the frustum.
+	 * </p>
+	 */
+    public static class Plane {
+        /** Normal of the plane */
+        public final float[] n = new float[3];
+
+        /** Distance to origin */
+        public float d;
+
+        /**
+         * Return signed distance of plane to given point.
+         * <ul>
+         *   <li>If dist &lt; 0 , then the point p lies in the negative halfspace.</li>
+         *   <li>If dist = 0 , then the point p lies in the plane.</li>
+         *   <li>If dist &gt; 0 , then the point p lies in the positive halfspace.</li>
+         * </ul>
+         * A plane cuts 3D space into 2 half spaces.
+         * <p>
+         * Positive halfspace is where the plane’s normals vector points into.
+         * </p>
+         * <p>
+         * Negative halfspace is the <i>other side</i> of the plane, i.e. *-1
+         * </p>
+         **/
+        public final float distanceTo(float x, float y, float z) {
+            return n[0] * x + n[1] * y + n[2] * z + d;
+        }
+
+        /** Return distance of plane to given point, see {@link #distanceTo(float, float, float)}. */
+        public final float distanceTo(float[] p) {
+            return n[0] * p[0] + n[1] * p[1] + n[2] * p[2] + d;
+        }
+
+        @Override
+        public String toString() {
+            return "Plane[ [ " + n[0] + ", " + n[1] + ", " + n[2] + " ], " + d + "]";
+        }
+    }
+
+    /** Index for left plane: {@value} */
+    public static final int LEFT   = 0;
+    /** Index for right plane: {@value} */
+    public static final int RIGHT  = 1;
+    /** Index for bottom plane: {@value} */
+    public static final int BOTTOM = 2;
+    /** Index for top plane: {@value} */
+    public static final int TOP    = 3;
+    /** Index for near plane: {@value} */
+    public static final int NEAR   = 4;
+    /** Index for far plane: {@value} */
+    public static final int FAR    = 5;
+
+    /**
+     * {@link Plane}s are ordered in the returned array as follows:
+     * <ul>
+     *   <li>{@link #LEFT}</li>
+     *   <li>{@link #RIGHT}</li>
+     *   <li>{@link #BOTTOM}</li>
+     *   <li>{@link #TOP}</li>
+     *   <li>{@link #NEAR}</li>
+     *   <li>{@link #FAR}</li>
+     * </ul>
+     * <p>
+     * {@link Plane}'s normals are pointing to the inside of the frustum
+     * in order to work w/ {@link #isOutside(AABBox) isOutside(..)} methods.
+     * </p>
+     *
+     * @return array of normalized {@link Plane}s, order see above.
+     */
+    public final Plane[] getPlanes() { return planes; }
+
+    /**
+     * Copy the given <code>src</code> planes into this this instance's planes.
+     * @param src the 6 source planes
+     */
+    public final void updateByPlanes(Plane[] src) {
+        for (int i = 0; i < 6; ++i) {
+            final Plane p0 = planes[i];
+            final float[] p0_n = p0.n;
+            final Plane p1 = src[i];
+            final float[] p1_n = p1.n;
+            p0_n[0] = p1_n[0];
+            p0_n[1] = p1_n[1];
+            p0_n[2] = p1_n[2];
+            p0.d = p1.d;
+        }
+    }
+
+    /**
+     * Calculate the frustum planes in world coordinates
+     * using the passed float[16] as premultiplied P*MV (column major order).
+     * <p>
+     * Frustum plane's normals will point to the inside of the viewing frustum,
+     * as required by this class.
+     * </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
+        {
+            final Plane p = planes[LEFT];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 0 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 0 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] + pmv[ pmv_off + 0 + 2 * 4 ];
+            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
+        {
+            final Plane p = planes[RIGHT];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 0 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 0 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] - pmv[ pmv_off + 0 + 2 * 4 ];
+            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
+        {
+            final Plane p = planes[BOTTOM];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 1 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 1 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] + pmv[ pmv_off + 1 + 2 * 4 ];
+            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
+        {
+            final Plane p = planes[TOP];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 1 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 1 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] - pmv[ pmv_off + 1 + 2 * 4 ];
+            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
+        {
+            final Plane p = planes[NEAR];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 2 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 2 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] + pmv[ pmv_off + 2 + 2 * 4 ];
+            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
+        {
+            final Plane p = planes[FAR];
+            final float[] p_n = p.n;
+            p_n[0] = pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 2 + 0 * 4 ];
+            p_n[1] = pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 2 + 1 * 4 ];
+            p_n[2] = pmv[ pmv_off + 3 + 2 * 4 ] - pmv[ pmv_off + 2 + 2 * 4 ];
+            p.d    = pmv[ pmv_off + 3 + 3 * 4 ] - pmv[ pmv_off + 2 + 3 * 4 ];
+        }
+
+        // Normalize all planes
+        for (int i = 0; i < 6; ++i) {
+            final Plane p = planes[i];
+            final float[] p_n = p.n;
+            final double invl = Math.sqrt(p_n[0] * p_n[0] + p_n[1] * p_n[1] + p_n[2] * p_n[2]);
+
+            p_n[0] /= invl;
+            p_n[1] /= invl;
+            p_n[2] /= invl;
+            p.d /= invl;
+        }
+    }
+
+	private static final boolean isOutsideImpl(Plane p, AABBox box) {
+	    final float[] low = box.getLow();
+	    final float[] high = box.getHigh();
+
+		if ( p.distanceTo(low[0],  low[1],  low[2])  > 0.0f ||
+		     p.distanceTo(high[0], low[1],  low[2])  > 0.0f ||
+		     p.distanceTo(low[0],  high[1], low[2])  > 0.0f ||
+		     p.distanceTo(high[0], high[1], low[2])  > 0.0f ||
+		     p.distanceTo(low[0],  low[1],  high[2]) > 0.0f ||
+		     p.distanceTo(high[0], low[1],  high[2]) > 0.0f ||
+		     p.distanceTo(low[0],  high[1], high[2]) > 0.0f ||
+		     p.distanceTo(high[0], high[1], high[2]) > 0.0f ) {
+			return false;
+		}
+		return true;
+	}
+
+	/**
+	 * Check to see if an axis aligned bounding box is completely outside of the frustum.
+	 * <p>
+	 * Note: If method returns false, the box may only be partially inside.
+	 * </p>
+	 */
+    public final boolean isAABBoxOutside(AABBox box) {
+        for (int i = 0; i < 6; ++i) {
+            if ( isOutsideImpl(planes[i], box) ) {
+                // fully outside
+                return true;
+            }
+        }
+        // We make no attempt to determine whether it's fully inside or not.
+        return false;
+    }
+
+
+    public static enum Location { OUTSIDE, INSIDE, INTERSECT };
+
+    /**
+     * Check to see if a point is outside, inside or on a plane of the frustum.
+     *
+     * @param p the point
+     * @return {@link Location} of point related to frustum planes
+     */
+    public final Location classifyPoint(float[] p) {
+        Location res = Location.INSIDE;
+
+        for (int i = 0; i < 6; ++i) {
+            final float d = planes[i].distanceTo(p);
+            if ( d < 0.0f ) {
+                return Location.OUTSIDE;
+            } else if ( d == 0.0f ) {
+                res = Location.INTERSECT;
+            }
+        }
+        return res;
+    }
+
+    /**
+     * Check to see if a point is outside of the frustum.
+     *
+     * @param p the point
+     * @return true if outside of the frustum, otherwise inside or on a plane
+     */
+    public final boolean isPointOutside(float[] p) {
+        return Location.OUTSIDE == classifyPoint(p);
+    }
+
+    /**
+     * Check to see if a sphere is outside, intersecting or inside of the frustum.
+     *
+     * @param p center of the sphere
+     * @param radius radius of the sphere
+     * @return {@link Location} of point related to frustum planes
+     */
+    public final Location classifySphere(float[] p, float radius) {
+        Location res = Location.INSIDE; // fully inside
+
+        for (int i = 0; i < 6; ++i) {
+            final float d = planes[i].distanceTo(p);
+            if ( d < -radius ) {
+                // fully outside
+                return Location.OUTSIDE;
+            } else if (d < radius ) {
+                // intersecting
+                res = Location.INTERSECT;
+            }
+        }
+        return res;
+    }
+
+    /**
+     * Check to see if a sphere is outside of the frustum.
+     *
+     * @param p center of the sphere
+     * @param radius radius of the sphere
+     * @return true if outside of the frustum, otherwise inside or intersecting
+     */
+    public final boolean isSphereOutside(float[] p, float radius) {
+        return Location.OUTSIDE == classifySphere(p, radius);
+    }
+
+    public StringBuilder toString(StringBuilder sb) {
+        if( null == sb ) {
+            sb = new StringBuilder();
+        }
+        sb.append("Frustum[ Planes[ ").append(Platform.NEWLINE)
+        .append(" L: ").append(planes[0]).append(", ").append(Platform.NEWLINE)
+        .append(" R: ").append(planes[1]).append(", ").append(Platform.NEWLINE)
+        .append(" B: ").append(planes[2]).append(", ").append(Platform.NEWLINE)
+        .append(" T: ").append(planes[3]).append(", ").append(Platform.NEWLINE)
+        .append(" N: ").append(planes[4]).append(", ").append(Platform.NEWLINE)
+        .append(" F: ").append(planes[5]).append("], ").append(Platform.NEWLINE)
+        .append("]");
+        return sb;
+    }
+
+	@Override
+	public String toString() {
+	    return toString(null).toString();
+	}
+}