Math: Complete Matrix4f w/ Vec[234]f and adopt it throughout Quaternion, Ray, AABBox, Frustum, Stereo*, ... adding hook to PMVMatrix

Motivation was to simplify matrix + vector math usage, ease review and avoid usage bugs.

Matrix4f implementation uses dedicated float fields instead of an array.

Performance didn't increase much,
as JVM >= 11(?) has some optimizations to drop the array bounds check.

AMD64 + OpenJDK17
- Matrix4f.mul(a, b) got a roughly ~10% enhancement over FloatUtil.multMatrix(a, b, dest)
- Matrix4f.mul(b) roughly ~3% slower than FloatUtil.multMatrix(a, b, dest)
- FloatUtil.multMatrix(a, a_off, b, b_off, dest) is considerable slower than all
- Matrix4f.invert(..) roughly ~3% slower than FloatUtil.invertMatrix(..)

RaspberryPi 4b aarch64 + OpenJDK17
- Matrix4f.mul(a, b) got a roughly ~10% enhancement over FloatUtil.multMatrix(a, b, dest)
- Matrix4f.mul(b) roughly ~20% slower than FloatUtil.multMatrix(a, b)
- FloatUtil.multMatrix(a, a_off, b, b_off, dest) is considerable slower than all
- Matrix4f.invert(..) roughly ~4% slower than FloatUtil.invertMatrix(..)

Conclusion
- Matrix4f.mul(b) needs to be revised (esp for aarch64)
- Matrix4f.invert(..) should also not be slower ..

20 files changed, 3238 insertions, 871 deletions

diff --git a/src/jogl/classes/com/jogamp/graph/font/Font.java b/src/jogl/classes/com/jogamp/graph/font/Font.java
index 4399bbad7..3d7e9f6a1 100644
--- a/src/jogl/classes/com/jogamp/graph/font/Font.java
+++ b/src/jogl/classes/com/jogamp/graph/font/Font.java
@@ -29,6 +29,7 @@ package com.jogamp.graph.font;
 
 import com.jogamp.graph.curve.OutlineShape;
 import com.jogamp.graph.geom.plane.AffineTransform;
+import com.jogamp.opengl.math.Vec3f;
 import com.jogamp.opengl.math.geom.AABBox;
 
 /**
@@ -151,7 +152,7 @@ public interface Font {
          * @param dest AABBox instance set to this metrics boundary in font em-size [0..1]
          * @return the given and set AABBox 'dest' in font units
          */
-        AABBox getBounds(final AABBox dest, final float[] tmpV3);
+        AABBox getBounds(final AABBox dest);
     }
 
     /**
@@ -193,10 +194,9 @@ public interface Font {
         /**
          * Return the AABBox in font em-size [0..1], copying into given dest.
          * @param dest AABBox instance set to this metrics boundary in font em-size [0..1]
-         * @param tmpV3 caller provided temporary 3-component vector
          * @return the given and set AABBox 'dest' in font em-size [0..1]
          */
-        AABBox getBounds(final AABBox dest, float[] tmpV3);
+        AABBox getBounds(final AABBox dest);
 
         /**
          * Return the AABBox in font em-size [0..1], creating a new copy.
diff --git a/src/jogl/classes/com/jogamp/graph/geom/plane/AffineTransform.java b/src/jogl/classes/com/jogamp/graph/geom/plane/AffineTransform.java
index 62cda0322..74036f97d 100644
--- a/src/jogl/classes/com/jogamp/graph/geom/plane/AffineTransform.java
+++ b/src/jogl/classes/com/jogamp/graph/geom/plane/AffineTransform.java
@@ -16,6 +16,7 @@
  */
 /**
  * @author Denis M. Kishenko
+ * @author Sven Gothel, (c) 2010-2023
  */
 package com.jogamp.graph.geom.plane;
 
@@ -23,6 +24,8 @@ package com.jogamp.graph.geom.plane;
 
 import com.jogamp.graph.geom.Vertex;
 import com.jogamp.opengl.math.FloatUtil;
+import com.jogamp.opengl.math.Vec2f;
+import com.jogamp.opengl.math.Vec3f;
 import com.jogamp.opengl.math.geom.AABBox;
 
 public class AffineTransform implements Cloneable {
@@ -53,17 +56,17 @@ public class AffineTransform implements Cloneable {
     /**
      * The values of transformation matrix
      */
-    float m00;
-    float m10;
-    float m01;
-    float m11;
-    float m02;
-    float m12;
+    private float m00;
+    private float m10;
+    private float m01;
+    private float m11;
+    private float m02;
+    private float m12;
 
     /**
      * The transformation <code>type</code>
      */
-    transient int type;
+    private transient int type;
 
     public AffineTransform() {
         setToIdentity();
@@ -404,10 +407,10 @@ public class AffineTransform implements Cloneable {
      * @return dst for chaining
      */
     public final AABBox transform(final AABBox src, final AABBox dst) {
-        final float[] srcLo = src.getLow();
-        final float[] srcHi = src.getHigh();
-        dst.setSize(srcLo[0] * m00 + srcLo[1] * m01 + m02, srcLo[0] * m10 + srcLo[1] * m11 + m12, srcLo[2],
-                    srcHi[0] * m00 + srcHi[1] * m01 + m02, srcHi[0] * m10 + srcHi[1] * m11 + m12, srcHi[2]);
+        final Vec3f lo = src.getLow();
+        final Vec3f hi = src.getHigh();
+        dst.setSize(lo.x() * m00 + lo.y() * m01 + m02, lo.x() * m10 + lo.y() * m11 + m12, lo.z(),
+                    hi.x() * m00 + hi.y() * m01 + m02, hi.x() * m10 + hi.y() * m11 + m12, hi.z());
         return dst;
     }
 
@@ -475,6 +478,33 @@ public class AffineTransform implements Cloneable {
     }
 
     /**
+     * @param src source of transformation
+     * @param dst destination of transformation, maybe be equal to <code>src</code>
+     * @return dst for chaining
+     */
+    public final Vec2f transform(final Vec2f src, final Vec2f dst) {
+        final float x = src.x();
+        final float y = src.y();
+        dst.setX( x * m00 + y * m01 + m02 );
+        dst.setY( x * m10 + y * m11 + m12 );
+        return dst;
+    }
+
+    /**
+     * @param src source of transformation
+     * @param dst destination of transformation, maybe be equal to <code>src</code>
+     * @return dst for chaining
+     */
+    public final Vec3f transform(final Vec3f src, final Vec3f dst) {
+        final float x = src.x();
+        final float y = src.y();
+        dst.setX( x * m00 + y * m01 + m02 );
+        dst.setY( x * m10 + y * m11 + m12 );
+        dst.setZ( src.z() ); // just copy z
+        return dst;
+    }
+
+    /**
      *
      * @param src
      * @param dst
diff --git a/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
index 7aba7fa73..f793629d6 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/FloatUtil.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2010 JogAmp Community. All rights reserved.
+ * Copyright 2010-2023 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:
@@ -44,20 +44,26 @@ import com.jogamp.opengl.math.geom.Frustum;
  * Implementation assumes linear matrix layout in column-major order
  * matching OpenGL's implementation, illustration:
  * <pre>
-    Row-Major                  Column-Major (OpenGL):
+    Row-Major                       Column-Major (OpenGL):
 
-        |  0   1   2   3 |         |  0   4   8  12 |
-        |                |         |                |
-        |  4   5   6   7 |         |  1   5   9  13 |
-    M = |                |     M = |                |
-        |  8   9  10  11 |         |  2   6  10  14 |
-        |                |         |                |
-        | 12  13  14  15 |         |  3   7  11  15 |
+        |  0   1   2  tx |
+        |                |
+        |  4   5   6  ty |
+    M = |                |
+        |  8   9  10  tz |
+        |                |
+        | 12  13  14  15 |
 
-           C   R                      C   R
+           R   C                      R   C
          m[0*4+3] = tx;             m[0+4*3] = tx;
          m[1*4+3] = ty;             m[1+4*3] = ty;
          m[2*4+3] = tz;             m[2+4*3] = tz;
+
+          RC (std subscript order)   RC (std subscript order)
+         m03 = tx;                  m03 = tx;
+         m13 = ty;                  m13 = ty;
+         m23 = tz;                  m23 = tz;
+
  * </pre>
  * </p>
  * <p>
@@ -71,7 +77,7 @@ import com.jogamp.opengl.math.geom.Frustum;
  * Implementation utilizes unrolling of small vertices and matrices wherever possible
  * while trying to access memory in a linear fashion for performance reasons, see:
  * <ul>
- *   <li><a href="https://code.google.com/p/java-matrix-benchmark/">java-matrix-benchmark</a></li>
+ *   <li><a href="https://lessthanoptimal.github.io/Java-Matrix-Benchmark/">java-matrix-benchmark</a></li>
  *   <li><a href="https://github.com/lessthanoptimal/ejml">EJML Efficient Java Matrix Library</a></li>
  * </ul>
  * </p>
@@ -580,35 +586,6 @@ public final class FloatUtil {
   }
 
   /**
-   * Make given matrix the perspective {@link #makeFrustum(float[], int, boolean, float, float, float, float, float, float) frustum}
-   * matrix based on given parameters.
-   * <p>
-   * All matrix fields are only set if <code>initM</code> is <code>true</code>.
-   * </p>
-   *
-   * @param m 4x4 matrix in column-major order (also result)
-   * @param m_offset offset in given array <i>m</i>, i.e. start of the 4x4 matrix
-   * @param initM if true, given matrix will be initialized w/ identity matrix,
-   *              otherwise only the frustum fields are set.
-   * @param fovhv {@link FovHVHalves} field of view in both directions, may not be centered, either in radians or tangent
-   * @param zNear
-   * @param zFar
-   * @return given matrix for chaining
-   * @throws GLException if {@code zNear <= 0} or {@code zFar <= zNear}
-   * @see #makeFrustum(float[], int, boolean, float, float, float, float, float, float)
-   * @see Frustum#updateByFovDesc(float[], int, boolean, Frustum.FovDesc)
-   */
-  public static float[] makePerspective(final float[] m, final int m_offset, final boolean initM,
-                                        final FovHVHalves fovhv, final float zNear, final float zFar) throws GLException {
-      final FovHVHalves fovhvTan = fovhv.toTangents();  // use tangent of half-fov !
-      final float top    =         fovhvTan.top    * zNear;
-      final float bottom = -1.0f * fovhvTan.bottom * zNear;
-      final float left   = -1.0f * fovhvTan.left   * zNear;
-      final float right  =         fovhvTan.right  * zNear;
-      return makeFrustum(m, m_offset, initM, left, right, bottom, top, zNear, zFar);
-  }
-
-  /**
    * Make given matrix the <i>look-at</i> matrix based on given parameters.
    * <p>
    * Consist out of two matrix multiplications:
@@ -1072,12 +1049,12 @@ public final class FloatUtil {
    * @param vec4Tmp2 4 component vector for temp storage
    * @return true if successful, otherwise false (z is 1)
    */
-  public static boolean mapObjToWinCoords(final float objx, final float objy, final float objz,
-                                          final float[] modelMatrix, final int modelMatrix_offset,
-                                          final float[] projMatrix, final int projMatrix_offset,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float[] win_pos, final int win_pos_offset,
-                                          final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
+  public static boolean mapObjToWin(final float objx, final float objy, final float objz,
+                                    final float[] modelMatrix, final int modelMatrix_offset,
+                                    final float[] projMatrix, final int projMatrix_offset,
+                                    final int[] viewport, final int viewport_offset,
+                                    final float[] win_pos, final int win_pos_offset,
+                                    final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
       vec4Tmp1[0] = objx;
       vec4Tmp1[1] = objy;
       vec4Tmp1[2] = objz;
@@ -1120,18 +1097,15 @@ public final class FloatUtil {
    * @param objz
    * @param mat4PMv [projection] x [modelview] matrix, i.e. P x Mv
    * @param viewport 4 component viewport vector
-   * @param viewport_offset
    * @param win_pos 3 component window coordinate, the result
-   * @param win_pos_offset
    * @param vec4Tmp1 4 component vector for temp storage
    * @param vec4Tmp2 4 component vector for temp storage
    * @return true if successful, otherwise false (z is 1)
    */
-  public static boolean mapObjToWinCoords(final float objx, final float objy, final float objz,
-                                          final float[/*16*/] mat4PMv,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float[] win_pos, final int win_pos_offset,
-                                          final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
+  public static boolean mapObjToWin(final float objx, final float objy, final float objz,
+                                    final float[/*16*/] mat4PMv,
+                                    final int[] viewport, final float[] win_pos,
+                                    final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
       vec4Tmp2[0] = objx;
       vec4Tmp2[1] = objy;
       vec4Tmp2[2] = objz;
@@ -1152,9 +1126,9 @@ public final class FloatUtil {
       vec4Tmp1[2] = vec4Tmp1[2] * vec4Tmp1[3] + 0.5f;
 
       // Map x,y to viewport
-      win_pos[0+win_pos_offset] = vec4Tmp1[0] * viewport[2+viewport_offset] + viewport[0+viewport_offset];
-      win_pos[1+win_pos_offset] = vec4Tmp1[1] * viewport[3+viewport_offset] + viewport[1+viewport_offset];
-      win_pos[2+win_pos_offset] = vec4Tmp1[2];
+      win_pos[0] = vec4Tmp1[0] * viewport[2] + viewport[0];
+      win_pos[1] = vec4Tmp1[1] * viewport[3] + viewport[1];
+      win_pos[2] = vec4Tmp1[2];
 
       return true;
   }
@@ -1180,12 +1154,12 @@ public final class FloatUtil {
    * @param mat4Tmp2 16 component matrix for temp storage
    * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
    */
-  public static boolean mapWinToObjCoords(final float winx, final float winy, final float winz,
-                                          final float[] modelMatrix, final int modelMatrix_offset,
-                                          final float[] projMatrix, final int projMatrix_offset,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float[] obj_pos, final int obj_pos_offset,
-                                          final float[/*16*/] mat4Tmp1, final float[/*16*/] mat4Tmp2) {
+  public static boolean mapWinToObj(final float winx, final float winy, final float winz,
+                                    final float[] modelMatrix, final int modelMatrix_offset,
+                                    final float[] projMatrix, final int projMatrix_offset,
+                                    final int[] viewport, final int viewport_offset,
+                                    final float[] obj_pos, final int obj_pos_offset,
+                                    final float[/*16*/] mat4Tmp1, final float[/*16*/] mat4Tmp2) {
     // mat4Tmp1 = P x Mv
     multMatrix(projMatrix, projMatrix_offset, modelMatrix, modelMatrix_offset, mat4Tmp1, 0);
 
@@ -1242,11 +1216,11 @@ public final class FloatUtil {
    * @param vec4Tmp2 4 component vector for temp storage
    * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
    */
-  public static boolean mapWinToObjCoords(final float winx, final float winy, final float winz,
-                                          final float[/*16*/] mat4PMvI,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float[] obj_pos, final int obj_pos_offset,
-                                          final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
+  public static boolean mapWinToObj(final float winx, final float winy, final float winz,
+                                    final float[/*16*/] mat4PMvI,
+                                    final int[] viewport, final int viewport_offset,
+                                    final float[] obj_pos, final int obj_pos_offset,
+                                    final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
     vec4Tmp1[0] = winx;
     vec4Tmp1[1] = winy;
     vec4Tmp1[2] = winz;
@@ -1296,19 +1270,17 @@ public final class FloatUtil {
    * @param vec4Tmp2 4 component vector for temp storage
    * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
    */
-  public static boolean mapWinToObjCoords(final float winx, final float winy, final float winz1, final float winz2,
-                                          final float[/*16*/] mat4PMvI,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float[] obj1_pos, final int obj1_pos_offset,
-                                          final float[] obj2_pos, final int obj2_pos_offset,
-                                          final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
+  public static boolean mapWinToObj(final float winx, final float winy, final float winz1, final float winz2,
+                                    final float[/*16*/] mat4PMvI, final int[] viewport,
+                                    final Vec3f objPos1, final Vec3f objPos2,
+                                    final float[/*4*/] vec4Tmp1, final float[/*4*/] vec4Tmp2) {
     vec4Tmp1[0] = winx;
     vec4Tmp1[1] = winy;
     vec4Tmp1[3] = 1.0f;
 
     // Map x and y from window coordinates
-    vec4Tmp1[0] = (vec4Tmp1[0] - viewport[0+viewport_offset]) / viewport[2+viewport_offset];
-    vec4Tmp1[1] = (vec4Tmp1[1] - viewport[1+viewport_offset]) / viewport[3+viewport_offset];
+    vec4Tmp1[0] = (vec4Tmp1[0] - viewport[0]) / viewport[2];
+    vec4Tmp1[1] = (vec4Tmp1[1] - viewport[1]) / viewport[3];
 
     // Map to range -1 to 1
     vec4Tmp1[0] = vec4Tmp1[0] * 2 - 1;
@@ -1329,15 +1301,14 @@ public final class FloatUtil {
 
     vec4Tmp2[3] = 1.0f / vec4Tmp2[3];
 
-    obj1_pos[0+obj1_pos_offset] = vec4Tmp2[0] * vec4Tmp2[3];
-    obj1_pos[1+obj1_pos_offset] = vec4Tmp2[1] * vec4Tmp2[3];
-    obj1_pos[2+obj1_pos_offset] = vec4Tmp2[2] * vec4Tmp2[3];
+    objPos1.set( vec4Tmp2[0] * vec4Tmp2[3],
+                 vec4Tmp2[1] * vec4Tmp2[3],
+                 vec4Tmp2[2] * vec4Tmp2[3] );
 
     //
     // winz2
     //
-    vec4Tmp1[2] = winz2;
-    vec4Tmp1[2] = vec4Tmp1[2] * 2 - 1;
+    vec4Tmp1[2] = winz2 * 2 - 1;
 
     // object raw coords = Inv(P x Mv) *  winPos  -> mat4Tmp2
     multMatrixVec(mat4PMvI, vec4Tmp1, vec4Tmp2);
@@ -1348,9 +1319,9 @@ public final class FloatUtil {
 
     vec4Tmp2[3] = 1.0f / vec4Tmp2[3];
 
-    obj2_pos[0+obj2_pos_offset] = vec4Tmp2[0] * vec4Tmp2[3];
-    obj2_pos[1+obj2_pos_offset] = vec4Tmp2[1] * vec4Tmp2[3];
-    obj2_pos[2+obj2_pos_offset] = vec4Tmp2[2] * vec4Tmp2[3];
+    objPos2.set( vec4Tmp2[0] * vec4Tmp2[3],
+                 vec4Tmp2[1] * vec4Tmp2[3],
+                 vec4Tmp2[2] * vec4Tmp2[3] );
 
     return true;
   }
@@ -1379,13 +1350,13 @@ public final class FloatUtil {
    * @param mat4Tmp2 16 component matrix for temp storage
    * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
    */
-  public static boolean mapWinToObjCoords(final float winx, final float winy, final float winz, final float clipw,
-                                          final float[] modelMatrix, final int modelMatrix_offset,
-                                          final float[] projMatrix, final int projMatrix_offset,
-                                          final int[] viewport, final int viewport_offset,
-                                          final float near, final float far,
-                                          final float[] obj_pos, final int obj_pos_offset,
-                                          final float[/*16*/] mat4Tmp1, final float[/*16*/] mat4Tmp2) {
+  public static boolean mapWinToObj4(final float winx, final float winy, final float winz, final float clipw,
+                                     final float[] modelMatrix, final int modelMatrix_offset,
+                                     final float[] projMatrix, final int projMatrix_offset,
+                                     final int[] viewport, final int viewport_offset,
+                                     final float near, final float far,
+                                     final float[] obj_pos, final int obj_pos_offset,
+                                     final float[/*16*/] mat4Tmp1, final float[/*16*/] mat4Tmp2) {
     // mat4Tmp1 = P x Mv
     multMatrix(projMatrix, projMatrix_offset, modelMatrix, modelMatrix_offset, mat4Tmp1, 0);
 
@@ -1397,7 +1368,7 @@ public final class FloatUtil {
     mat4Tmp2[0] = winx;
     mat4Tmp2[1] = winy;
     mat4Tmp2[2] = winz;
-    mat4Tmp2[3] = 1.0f;
+    mat4Tmp2[3] = clipw;
 
     // Map x and y from window coordinates
     mat4Tmp2[0] = (mat4Tmp2[0] - viewport[0+viewport_offset]) / viewport[2+viewport_offset];
@@ -1417,8 +1388,6 @@ public final class FloatUtil {
       return false;
     }
 
-    mat4Tmp2[3+raw_off] = 1.0f / mat4Tmp2[3+raw_off];
-
     obj_pos[0+obj_pos_offset] = mat4Tmp2[0+raw_off];
     obj_pos[1+obj_pos_offset] = mat4Tmp2[1+raw_off];
     obj_pos[2+obj_pos_offset] = mat4Tmp2[2+raw_off];
@@ -1459,7 +1428,7 @@ public final class FloatUtil {
   public static boolean mapWinToRay(final float winx, final float winy, final float winz0, final float winz1,
                                     final float[] modelMatrix, final int modelMatrix_offset,
                                     final float[] projMatrix, final int projMatrix_offset,
-                                    final int[] viewport, final int viewport_offset,
+                                    final int[] viewport,
                                     final Ray ray,
                                     final float[/*16*/] mat4Tmp1, final float[/*16*/] mat4Tmp2, final float[/*4*/] vec4Tmp2) {
       // mat4Tmp1 = P x Mv
@@ -1469,11 +1438,9 @@ public final class FloatUtil {
       if ( null == invertMatrix(mat4Tmp1, mat4Tmp1) ) {
           return false;
       }
-      if( mapWinToObjCoords(winx, winy, winz0, winz1, mat4Tmp1,
-                            viewport, viewport_offset,
-                            ray.orig, 0, ray.dir, 0,
-                            mat4Tmp2, vec4Tmp2) ) {
-          VectorUtil.normalizeVec3( VectorUtil.subVec3(ray.dir, ray.dir, ray.orig) );
+      if( mapWinToObj(winx, winy, winz0, winz1, mat4Tmp1, viewport,
+                      ray.orig, ray.dir, mat4Tmp2, vec4Tmp2) ) {
+          ray.dir.sub(ray.orig).normalize();
           return true;
       } else {
           return false;
@@ -1485,9 +1452,8 @@ public final class FloatUtil {
    * @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
-   * @return given result matrix <i>d</i> for chaining
    */
-  public static float[] multMatrix(final float[] a, final int a_off, final float[] b, final int b_off, final float[] d, final int d_off) {
+  public static void multMatrix(final float[] a, final int a_off, final float[] b, final int b_off, final float[] d, final int d_off) {
       final float b00 = b[b_off+0+0*4];
       final float b10 = b[b_off+1+0*4];
       final float b20 = b[b_off+2+0*4];
@@ -1540,8 +1506,6 @@ public final class FloatUtil {
       d[d_off+3+1*4] = ai0 * b01  +  ai1 * b11  +  ai2 * b21  +  ai3 * b31 ;
       d[d_off+3+2*4] = ai0 * b02  +  ai1 * b12  +  ai2 * b22  +  ai3 * b32 ;
       d[d_off+3+3*4] = ai0 * b03  +  ai1 * b13  +  ai2 * b23  +  ai3 * b33 ;
-
-      return d;
   }
 
   /**
@@ -1612,9 +1576,8 @@ public final class FloatUtil {
    * Multiply matrix: [a] = [a] x [b]
    * @param a 4x4 matrix in column-major order (also result)
    * @param b 4x4 matrix in column-major order
-   * @return given result matrix <i>a</i> for chaining
    */
-  public static float[] multMatrix(final float[] a, final int a_off, final float[] b, final int b_off) {
+  public static void multMatrix(final float[] a, final int a_off, final float[] b, final int b_off) {
       final float b00 = b[b_off+0+0*4];
       final float b10 = b[b_off+1+0*4];
       final float b20 = b[b_off+2+0*4];
@@ -1667,8 +1630,6 @@ public final class FloatUtil {
       a[a_off+3+1*4] = ai0 * b01  +  ai1 * b11  +  ai2 * b21  +  ai3 * b31 ;
       a[a_off+3+2*4] = ai0 * b02  +  ai1 * b12  +  ai2 * b22  +  ai3 * b32 ;
       a[a_off+3+3*4] = ai0 * b03  +  ai1 * b13  +  ai2 * b23  +  ai3 * b33 ;
-
-      return a;
   }
 
   /**
@@ -1778,11 +1739,10 @@ public final class FloatUtil {
    * @param m_in_off
    * @param v_in 4-component column-vector
    * @param v_out m_in * v_in
-   * @return given result vector <i>v_out</i> for chaining
    */
-  public static float[] multMatrixVec(final float[] m_in, final int m_in_off,
-                                      final float[] v_in, final int v_in_off,
-                                      final float[] v_out, final int v_out_off) {
+  public static void multMatrixVec(final float[] m_in, final int m_in_off,
+                                   final float[] v_in, final int v_in_off,
+                                   final float[] v_out, final int v_out_off) {
       // (one matrix row in column-major order) X (column vector)
       v_out[0 + v_out_off] = v_in[0+v_in_off] * m_in[0*4+m_in_off  ]  +  v_in[1+v_in_off] * m_in[1*4+m_in_off  ] +
                              v_in[2+v_in_off] * m_in[2*4+m_in_off  ]  +  v_in[3+v_in_off] * m_in[3*4+m_in_off  ];
@@ -1798,8 +1758,31 @@ public final class FloatUtil {
       final int m_in_off_3 = 3+m_in_off;
       v_out[3 + v_out_off] = v_in[0+v_in_off] * m_in[0*4+m_in_off_3]  +  v_in[1+v_in_off] * m_in[1*4+m_in_off_3] +
                              v_in[2+v_in_off] * m_in[2*4+m_in_off_3]  +  v_in[3+v_in_off] * m_in[3*4+m_in_off_3];
+  }
 
-      return v_out;
+  /**
+   * @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 void multMatrixVec(final float[] m_in, final int m_in_off,
+                                   final float[] v_in, final float[] v_out) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[0] = v_in[0] * m_in[0*4+m_in_off  ]  +  v_in[1] * m_in[1*4+m_in_off  ] +
+                 v_in[2] * m_in[2*4+m_in_off  ]  +  v_in[3] * m_in[3*4+m_in_off  ];
+
+      final int m_in_off_1 = 1+m_in_off;
+      v_out[1] = v_in[0] * m_in[0*4+m_in_off_1]  +  v_in[1] * m_in[1*4+m_in_off_1] +
+                 v_in[2] * m_in[2*4+m_in_off_1]  +  v_in[3] * m_in[3*4+m_in_off_1];
+
+      final int m_in_off_2 = 2+m_in_off;
+      v_out[2] = v_in[0] * m_in[0*4+m_in_off_2]  +  v_in[1] * m_in[1*4+m_in_off_2] +
+                 v_in[2] * m_in[2*4+m_in_off_2]  +  v_in[3] * m_in[3*4+m_in_off_2];
+
+      final int m_in_off_3 = 3+m_in_off;
+      v_out[3] = v_in[0] * m_in[0*4+m_in_off_3]  +  v_in[1] * m_in[1*4+m_in_off_3] +
+                 v_in[2] * m_in[2*4+m_in_off_3]  +  v_in[3] * m_in[3*4+m_in_off_3];
   }
 
   /**
@@ -1845,46 +1828,59 @@ public final class FloatUtil {
   }
 
   /**
-   * Copy the named column of the given column-major matrix to v_out.
-   * <p>
-   * v_out may be 3 or 4 components long, hence the 4th row may not be stored.
-   * </p>
-   * @param m_in input column-major matrix
-   * @param m_in_off offset to input matrix
-   * @param column named column to copy
-   * @param v_out the column-vector storage, at least 3 components long
-   * @param v_out_off offset to storage
+   * Affine 3f-vector transformation by 4x4 matrix
+   *
+   * 4x4 matrix multiplication with 3-component vector,
+   * using {@code 1} for for {@code v_in[3]} and dropping {@code v_out[3]},
+   * which shall be {@code 1}.
+   *
+   * @param m_in 4x4 matrix in column-major order
+   * @param m_in_off
+   * @param v_in 3-component column-vector
+   * @param v_out m_in * v_in, 3-component column-vector
    * @return given result vector <i>v_out</i> for chaining
    */
-  public static float[] copyMatrixColumn(final float[] m_in, final int m_in_off, final int column, final float[] v_out, final int v_out_off) {
-      v_out[0+v_out_off]=m_in[0+column*4+m_in_off];
-      v_out[1+v_out_off]=m_in[1+column*4+m_in_off];
-      v_out[2+v_out_off]=m_in[2+column*4+m_in_off];
-      if( v_out.length > 3+v_out_off ) {
-          v_out[3+v_out_off]=m_in[3+column*4+m_in_off];
-      }
+  public static float[] multMatrixVec3(final float[] m_in, final int m_in_off,
+                                       final float[] v_in, final float[] v_out) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[0] = v_in[0] * m_in[0*4+m_in_off  ]  +  v_in[1] * m_in[1*4+m_in_off  ] +
+                 v_in[2] * m_in[2*4+m_in_off  ]  +       1f * m_in[3*4+m_in_off  ];
+
+      final int m_in_off_1 = 1+m_in_off;
+      v_out[1] = v_in[0] * m_in[0*4+m_in_off_1]  +  v_in[1] * m_in[1*4+m_in_off_1] +
+                 v_in[2] * m_in[2*4+m_in_off_1]  +       1f * m_in[3*4+m_in_off_1];
+
+      final int m_in_off_2 = 2+m_in_off;
+      v_out[2] = v_in[0] * m_in[0*4+m_in_off_2]  +  v_in[1] * m_in[1*4+m_in_off_2] +
+                 v_in[2] * m_in[2*4+m_in_off_2]  +       1f * m_in[3*4+m_in_off_2];
+
       return v_out;
   }
 
   /**
-   * Copy the named row of the given column-major matrix to v_out.
-   * <p>
-   * v_out may be 3 or 4 components long, hence the 4th column may not be stored.
-   * </p>
-   * @param m_in input column-major matrix
-   * @param m_in_off offset to input matrix
-   * @param row named row to copy
-   * @param v_out the row-vector storage, at least 3 components long
-   * @param v_out_off offset to storage
+   * Affine 3f-vector transformation by 4x4 matrix
+   *
+   * 4x4 matrix multiplication with 3-component vector,
+   * using {@code 1} for for {@code v_in[3]} and dropping {@code v_out[3]},
+   * which shall be {@code 1}.
+   *
+   * @param m_in 4x4 matrix in column-major order
+   * @param m_in_off
+   * @param v_in 3-component column-vector
+   * @param v_out m_in * v_in, 3-component column-vector
    * @return given result vector <i>v_out</i> for chaining
    */
-  public static float[] copyMatrixRow(final float[] m_in, final int m_in_off, final int row, final float[] v_out, final int v_out_off) {
-      v_out[0+v_out_off]=m_in[row+0*4+m_in_off];
-      v_out[1+v_out_off]=m_in[row+1*4+m_in_off];
-      v_out[2+v_out_off]=m_in[row+2*4+m_in_off];
-      if( v_out.length > 3+v_out_off ) {
-          v_out[3+v_out_off]=m_in[row+3*4+m_in_off];
-      }
+  public static float[] multMatrixVec3(final float[] m_in, final float[] v_in, final float[] v_out) {
+      // (one matrix row in column-major order) X (column vector)
+      v_out[0] = v_in[0] * m_in[0*4  ]  +  v_in[1] * m_in[1*4  ] +
+                 v_in[2] * m_in[2*4  ]  +       1f * m_in[3*4  ];
+
+      v_out[1] = v_in[0] * m_in[0*4+1]  +  v_in[1] * m_in[1*4+1] +
+                 v_in[2] * m_in[2*4+1]  +       1f * m_in[3*4+1];
+
+      v_out[2] = v_in[0] * m_in[0*4+2]  +  v_in[1] * m_in[1*4+2] +
+                 v_in[2] * m_in[2*4+2]  +       1f * m_in[3*4+2];
+
       return v_out;
   }
 
@@ -2276,7 +2272,7 @@ public final class FloatUtil {
   }
 
   /**
-   * Return true if value is zero, i.e. it's absolute value < <code>epsilon</code>.
+   * Return true if value is zero, i.e. it's absolute value < {@link #EPSILON}.
    * @see #EPSILON
    */
   public static boolean isZero(final float a) {
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Matrix4.java b/src/jogl/classes/com/jogamp/opengl/math/Matrix4.java
deleted file mode 100644
index a080d4442..000000000
--- a/src/jogl/classes/com/jogamp/opengl/math/Matrix4.java
+++ /dev/null
@@ -1,172 +0,0 @@
-/**
- * Copyright 2014 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 com.jogamp.opengl.GLException;
-import com.jogamp.opengl.fixedfunc.GLMatrixFunc;
-
-import com.jogamp.opengl.util.PMVMatrix;
-
-/**
- * Simple float array-backed float 4x4 matrix
- * exposing {@link FloatUtil} matrix functionality in an object oriented manner.
- * <p>
- * Unlike {@link PMVMatrix}, this class only represents one single matrix
- * without a complete {@link GLMatrixFunc} implementation,
- * allowing this class to be more lightweight.
- * </p>
- * <p>
- * Implementation is not mature - WIP and subject to change.
- * </p>
- */
-public class Matrix4 {
-
-    public Matrix4() {
-          matrix     = new float[16];
-          matrixTxSx = new float[16];
-          mat4Tmp1   = new float[16];
-          vec4Tmp1   = new float[4];
-          FloatUtil.makeIdentity(matrixTxSx);
-          loadIdentity();
-    }
-
-    public final float[] getMatrix() {
-        return matrix;
-    }
-
-    public final void loadIdentity() {
-        FloatUtil.makeIdentity(matrix);
-    }
-
-    /**
-     * Multiply matrix: [this] = [this] x [m]
-     * @param m 4x4 matrix in column-major order
-     */
-    public final void multMatrix(final float[] m, final int m_offset) {
-        FloatUtil.multMatrix(matrix, 0, m, m_offset);
-    }
-
-    /**
-     * Multiply matrix: [this] = [this] x [m]
-     * @param m 4x4 matrix in column-major order
-     */
-    public final void multMatrix(final float[] m) {
-        FloatUtil.multMatrix(matrix, m);
-    }
-
-    /**
-     * Multiply matrix: [this] = [this] x [m]
-     * @param m 4x4 matrix in column-major order
-     */
-    public final void multMatrix(final Matrix4 m) {
-        FloatUtil.multMatrix(matrix, m.getMatrix());
-    }
-
-    /**
-     * @param v_in 4-component column-vector
-     * @param v_out this * v_in
-     */
-    public final void multVec(final float[] v_in, final float[] v_out) {
-        FloatUtil.multMatrixVec(matrix, v_in, v_out);
-    }
-
-    /**
-     * @param v_in 4-component column-vector
-     * @param v_out this * v_in
-     */
-    public final void multVec(final float[] v_in, final int v_in_offset, final float[] v_out, final int v_out_offset) {
-        FloatUtil.multMatrixVec(matrix, 0, v_in, v_in_offset, v_out, v_out_offset);
-    }
-
-    public final void translate(final float x, final float y, final float z) {
-        multMatrix(FloatUtil.makeTranslation(matrixTxSx, false, x, y, z));
-    }
-
-    public final void scale(final float x, final float y, final float z) {
-        multMatrix(FloatUtil.makeScale(matrixTxSx, false, x, y, z));
-    }
-
-    public final void rotate(final float angrad, final float x, final float y, final float z) {
-        multMatrix(FloatUtil.makeRotationAxis(mat4Tmp1, 0, angrad, x, y, z, vec4Tmp1));
-    }
-
-    /**
-     * Rotate the current matrix with the given {@link Quaternion}'s rotation {@link Quaternion#toMatrix(float[], int) matrix representation}.
-     */
-    public final void rotate(final Quaternion quat) {
-        multMatrix(quat.toMatrix(mat4Tmp1, 0));
-    }
-
-    public final void transpose() {
-        System.arraycopy(matrix, 0, mat4Tmp1, 0, 16);
-        FloatUtil.transposeMatrix(mat4Tmp1, matrix);
-    }
-
-    public final float determinant() {
-        return FloatUtil.matrixDeterminant(matrix);
-    }
-
-    public final boolean invert() {
-        return null != FloatUtil.invertMatrix(matrix, matrix);
-    }
-
-    public final void makeOrtho(final float left, final float right, final float bottom, final float top, final float zNear, final float zFar) {
-        multMatrix( FloatUtil.makeOrtho(mat4Tmp1, 0, true, left, right, bottom, top, zNear, zFar) );
-    }
-
-    /**
-     * @param left
-     * @param right
-     * @param bottom
-     * @param top
-     * @param zNear
-     * @param zFar
-     * @throws GLException if {@code zNear <= 0} or {@code zFar <= zNear}
-     *                     or {@code left == right}, or {@code bottom == top}.
-     * @see FloatUtil#makeFrustum(float[], int, boolean, float, float, float, float, float, float)
-     */
-    public final void makeFrustum(final float left, final float right, final float bottom, final float top, final float zNear, final float zFar) throws GLException {
-        multMatrix( FloatUtil.makeFrustum(mat4Tmp1, 0, true, left, right, bottom, top, zNear, zFar) );
-    }
-
-    /**
-     * @param fovy_rad
-     * @param aspect
-     * @param zNear
-     * @param zFar
-     * @throws GLException if {@code zNear <= 0} or {@code zFar <= zNear}
-     * @see FloatUtil#makePerspective(float[], int, boolean, float, float, float, float)
-     */
-    public final void makePerspective(final float fovy_rad, final float aspect, final float zNear, final float zFar) throws GLException {
-        multMatrix( FloatUtil.makePerspective(mat4Tmp1, 0, true, fovy_rad, aspect, zNear, zFar) );
-    }
-
-    private final float[] matrix, matrixTxSx;
-    private final float[] mat4Tmp1, vec4Tmp1;
-}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Matrix4f.java b/src/jogl/classes/com/jogamp/opengl/math/Matrix4f.java
new file mode 100644
index 000000000..5951c7d98
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Matrix4f.java
@@ -0,0 +1,1878 @@
+/**
+ * Copyright 2014-2023 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.opengl.math.geom.AABBox;
+import com.jogamp.opengl.math.geom.Frustum;
+import com.jogamp.opengl.math.geom.Frustum.Plane;
+
+/**
+ * Basic 4x4 float matrix implementation using fields for intensive use-cases (host operations).
+ * <p>
+ * Implementation covers {@link FloatUtil} matrix functionality, exposed in an object oriented manner.
+ * </p>
+ * <p>
+ * Unlike {@link com.jogamp.opengl.util.PMVMatrix PMVMatrix}, this class only represents one single matrix
+ * without a complete {@link com.jogamp.opengl.fixedfunc.GLMatrixFunc GLMatrixFunc} implementation.
+ * </p>
+ * <p>
+ * For array operations the layout is expected in column-major order
+ * matching OpenGL's implementation, illustration:
+ * <pre>
+    Row-Major                       Column-Major (OpenGL):
+
+        |  0   1   2  tx |
+        |                |
+        |  4   5   6  ty |
+    M = |                |
+        |  8   9  10  tz |
+        |                |
+        | 12  13  14  15 |
+
+           R   C                      R   C
+         m[0*4+3] = tx;             m[0+4*3] = tx;
+         m[1*4+3] = ty;             m[1+4*3] = ty;
+         m[2*4+3] = tz;             m[2+4*3] = tz;
+
+          RC (std subscript order)   RC (std subscript order)
+         m03 = tx;                  m03 = tx;
+         m13 = ty;                  m13 = ty;
+         m23 = tz;                  m23 = tz;
+
+ * </pre>
+ * </p>
+ * <p>
+ * <ul>
+ *   <li><a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html">Matrix-FAQ</a></li>
+ *   <li><a href="https://en.wikipedia.org/wiki/Matrix_%28mathematics%29">Wikipedia-Matrix</a></li>
+ *   <li><a href="http://www.euclideanspace.com/maths/algebra/matrix/index.htm">euclideanspace.com-Matrix</a></li>
+ * </ul>
+ * </p>
+ * <p>
+ * Implementation utilizes unrolling of small vertices and matrices wherever possible
+ * while trying to access memory in a linear fashion for performance reasons, see:
+ * <ul>
+ *   <li><a href="https://lessthanoptimal.github.io/Java-Matrix-Benchmark/">java-matrix-benchmark</a></li>
+ *   <li><a href="https://github.com/lessthanoptimal/ejml">EJML Efficient Java Matrix Library</a></li>
+ * </ul>
+ * </p>
+ * @see com.jogamp.opengl.util.PMVMatrix
+ * @see FloatUtil
+ */
+public class Matrix4f {
+
+    /**
+     * Creates a new identity matrix.
+     */
+    public Matrix4f() {
+        m00 = m11 = m22 = m33 = 1.0f;
+        // remaining fields have default init to zero
+    }
+
+    /**
+     * Creates a new matrix copying the values of the given {@code src} matrix.
+     */
+    public Matrix4f(final Matrix4f src) {
+        load(src);
+    }
+
+    /**
+     * Creates a new matrix based on given float[4*4] column major order.
+     * @param m 4x4 matrix in column-major order
+     */
+    public Matrix4f(final float[] m) {
+        load(m);
+    }
+
+    /**
+     * Creates a new matrix based on given float[4*4] column major order.
+     * @param m 4x4 matrix in column-major order
+     * @param m_off offset for matrix {@code m}
+     */
+    public Matrix4f(final float[] m, final int m_off) {
+        load(m, m_off);
+    }
+
+    //
+    // Write to Matrix via load(..)
+    //
+
+    /**
+     * Set this matrix to identity.
+     * <pre>
+      Translation matrix (Column Order):
+      1 0 0 0
+      0 1 0 0
+      0 0 1 0
+      0 0 0 1
+     * </pre>
+     * @return this matrix for chaining
+     */
+    public final Matrix4f loadIdentity() {
+       m00 = m11 = m22 = m33 = 1.0f;
+       m01 = m02 = m03 =
+       m10 = m12 = m13 =
+       m20 = m21 = m23 =
+       m30 = m31 = m32 = 0.0f;
+       return this;
+    }
+
+    /**
+     * Load the values of the given matrix {@code b} to this matrix.
+     * @param src the source values
+     * @return this matrix for chaining
+     */
+    public Matrix4f load(final Matrix4f src) {
+        m00 = src.m00; m10 = src.m10; m20 = src.m20; m30 = src.m30;
+        m01 = src.m01; m11 = src.m11; m21 = src.m21; m31 = src.m31;
+        m02 = src.m02; m12 = src.m12; m22 = src.m22; m32 = src.m32;
+        m03 = src.m03; m13 = src.m13; m23 = src.m23; m33 = src.m33;
+        return this;
+    }
+
+    /**
+     * Load the values of the given matrix {@code src} to this matrix.
+     * @param src 4x4 matrix float[16] in column-major order
+     * @return this matrix for chaining
+     */
+    public Matrix4f load(final float[] src) {
+        m00 = src[0+0*4]; // column 0
+        m10 = src[1+0*4];
+        m20 = src[2+0*4];
+        m30 = src[3+0*4];
+        m01 = src[0+1*4]; // column 1
+        m11 = src[1+1*4];
+        m21 = src[2+1*4];
+        m31 = src[3+1*4];
+        m02 = src[0+2*4]; // column 2
+        m12 = src[1+2*4];
+        m22 = src[2+2*4];
+        m32 = src[3+2*4];
+        m03 = src[0+3*4]; // column 3
+        m13 = src[1+3*4];
+        m23 = src[2+3*4];
+        m33 = src[3+3*4];
+        return this;
+    }
+
+    /**
+     * Load the values of the given matrix {@code src} to this matrix.
+     * @param src 4x4 matrix float[16] in column-major order
+     * @param src_off offset for matrix {@code src}
+     * @return this matrix for chaining
+     */
+    public Matrix4f load(final float[] src, final int src_off) {
+        m00 = src[src_off+0+0*4];
+        m10 = src[src_off+1+0*4];
+        m20 = src[src_off+2+0*4];
+        m30 = src[src_off+3+0*4];
+        m01 = src[src_off+0+1*4];
+        m11 = src[src_off+1+1*4];
+        m21 = src[src_off+2+1*4];
+        m31 = src[src_off+3+1*4];
+        m02 = src[src_off+0+2*4];
+        m12 = src[src_off+1+2*4];
+        m22 = src[src_off+2+2*4];
+        m32 = src[src_off+3+2*4];
+        m03 = src[src_off+0+3*4];
+        m13 = src[src_off+1+3*4];
+        m23 = src[src_off+2+3*4];
+        m33 = src[src_off+3+3*4];
+        return this;
+    }
+
+    /**
+     * Load the values of the given matrix {@code src} to this matrix.
+     * <p>
+     * Implementation uses relative {@link FloatBuffer#get()},
+     * hence caller may want to issue {@link FloatBuffer#reset()} thereafter.
+     * </p>
+     * @param src 4x4 matrix {@link FloatBuffer} in column-major order
+     * @return this matrix for chaining
+     */
+    public Matrix4f load(final FloatBuffer src) {
+        m00 = src.get();
+        m10 = src.get();
+        m20 = src.get();
+        m30 = src.get();
+        m01 = src.get();
+        m11 = src.get();
+        m21 = src.get();
+        m31 = src.get();
+        m02 = src.get();
+        m12 = src.get();
+        m22 = src.get();
+        m32 = src.get();
+        m03 = src.get();
+        m13 = src.get();
+        m23 = src.get();
+        m33 = src.get();
+        return this;
+    }
+
+    //
+    // Read out Matrix via get(..)
+    //
+
+    /** Gets the ith component, 0 <= i < 16 */
+    public float get(final int i) {
+        switch (i) {
+            case 0+4*0: return m00;
+            case 1+4*0: return m10;
+            case 2+4*0: return m20;
+            case 3+4*0: return m30;
+
+            case 0+4*1: return m01;
+            case 1+4*1: return m11;
+            case 2+4*1: return m21;
+            case 3+4*1: return m31;
+
+            case 0+4*2: return m02;
+            case 1+4*2: return m12;
+            case 2+4*2: return m22;
+            case 3+4*2: return m32;
+
+            case 0+4*3: return m03;
+            case 1+4*3: return m13;
+            case 2+4*3: return m23;
+            case 3+4*3: return m33;
+
+            default: throw new IndexOutOfBoundsException();
+        }
+    }
+
+    /**
+     * Get the named column of the given column-major matrix to v_out.
+     * @param column named column to copy
+     * @param v_out the column-vector storage
+     * @return given result vector <i>v_out</i> for chaining
+     */
+    public Vec4f getColumn(final int column, final Vec4f v_out) {
+        v_out.set( get(0+column*4),
+                   get(1+column*4),
+                   get(2+column*4),
+                   get(3+column*4) );
+        return v_out;
+    }
+
+    /**
+     * Get the named column of the given column-major matrix to v_out.
+     * @param column named column to copy
+     * @param v_out the column-vector storage
+     * @return given result vector <i>v_out</i> for chaining
+     */
+    public Vec3f getColumn(final int column, final Vec3f v_out) {
+        v_out.set( get(0+column*4),
+                   get(1+column*4),
+                   get(2+column*4) );
+        return v_out;
+    }
+
+    /**
+     * Get the named row of the given column-major matrix to v_out.
+     * @param row named row to copy
+     * @param v_out the row-vector storage
+     * @return given result vector <i>v_out</i> for chaining
+     */
+    public Vec4f getRow(final int row, final Vec4f v_out) {
+        v_out.set( get(row+0*4),
+                   get(row+1*4),
+                   get(row+2*4),
+                   get(row+3*4) );
+        return v_out;
+    }
+
+    /**
+     * Get the named row of the given column-major matrix to v_out.
+     * @param row named row to copy
+     * @param v_out the row-vector storage
+     * @return given result vector <i>v_out</i> for chaining
+     */
+    public Vec3f getRow(final int row, final Vec3f v_out) {
+        v_out.set( get(row+0*4),
+                   get(row+1*4),
+                   get(row+2*4) );
+        return v_out;
+    }
+
+    /**
+     * Get this matrix into the given float[16] array at {@code dst_off} in column major order.
+     *
+     * @param dst float[16] array storage in column major order
+     * @param dst_off offset
+     */
+    public void get(final float[] dst, final int dst_off) {
+        dst[dst_off+0+0*4] = m00;
+        dst[dst_off+1+0*4] = m10;
+        dst[dst_off+2+0*4] = m20;
+        dst[dst_off+3+0*4] = m30;
+        dst[dst_off+0+1*4] = m01;
+        dst[dst_off+1+1*4] = m11;
+        dst[dst_off+2+1*4] = m21;
+        dst[dst_off+3+1*4] = m31;
+        dst[dst_off+0+2*4] = m02;
+        dst[dst_off+1+2*4] = m12;
+        dst[dst_off+2+2*4] = m22;
+        dst[dst_off+3+2*4] = m32;
+        dst[dst_off+0+3*4] = m03;
+        dst[dst_off+1+3*4] = m13;
+        dst[dst_off+2+3*4] = m23;
+        dst[dst_off+3+3*4] = m33;
+    }
+
+    /**
+     * Get this matrix into the given float[16] array in column major order.
+     *
+     * @param dst float[16] array storage in column major order
+     */
+    public void get(final float[] dst) {
+        dst[0+0*4] = m00;
+        dst[1+0*4] = m10;
+        dst[2+0*4] = m20;
+        dst[3+0*4] = m30;
+        dst[0+1*4] = m01;
+        dst[1+1*4] = m11;
+        dst[2+1*4] = m21;
+        dst[3+1*4] = m31;
+        dst[0+2*4] = m02;
+        dst[1+2*4] = m12;
+        dst[2+2*4] = m22;
+        dst[3+2*4] = m32;
+        dst[0+3*4] = m03;
+        dst[1+3*4] = m13;
+        dst[2+3*4] = m23;
+        dst[3+3*4] = m33;
+    }
+
+    /**
+     * Get this matrix into the given {@link FloatBuffer} in column major order.
+     * <p>
+     * Implementation uses relative {@link FloatBuffer#put(float)},
+     * hence caller may want to issue {@link FloatBuffer#reset()} thereafter.
+     * </p>
+     *
+     * @param dst {@link FloatBuffer} array storage in column major order
+     */
+    public void get(final FloatBuffer dst) {
+        dst.put( m00 );
+        dst.put( m10 );
+        dst.put( m20 );
+        dst.put( m30 );
+        dst.put( m01 );
+        dst.put( m11 );
+        dst.put( m21 );
+        dst.put( m31 );
+        dst.put( m02 );
+        dst.put( m12 );
+        dst.put( m22 );
+        dst.put( m32 );
+        dst.put( m03 );
+        dst.put( m13 );
+        dst.put( m23 );
+        dst.put( m33 );
+    }
+
+    //
+    // Basic matrix operations
+    //
+
+    /**
+     * Returns the determinant of this matrix
+     * @return the matrix determinant
+     */
+    public float determinant() {
+        float ret = 0;
+        ret += m00 * ( + m11*(m22*m33 - m23*m32) - m12*(m21*m33 - m23*m31) + m13*(m21*m32 - m22*m31));
+        ret -= m01 * ( + m10*(m22*m33 - m23*m32) - m12*(m20*m33 - m23*m30) + m13*(m20*m32 - m22*m30));
+        ret += m02 * ( + m10*(m21*m33 - m23*m31) - m11*(m20*m33 - m23*m30) + m13*(m20*m31 - m21*m30));
+        ret -= m03 * ( + m10*(m21*m32 - m22*m31) - m11*(m20*m32 - m22*m30) + m12*(m20*m31 - m21*m30));
+        return ret;
+    }
+
+    /**
+     * Transpose this matrix.
+     *
+     * @return this matrix for chaining
+     */
+    public final Matrix4f transpose() {
+        float tmp;
+
+        tmp = m10;
+        m10 = m01;
+        m01 = tmp;
+
+        tmp = m20;
+        m20 = m02;
+        m02 = tmp;
+
+        tmp = m30;
+        m30 = m03;
+        m03 = tmp;
+
+        tmp = m21;
+        m21 = m12;
+        m12 = tmp;
+
+        tmp = m31;
+        m31 = m13;
+        m13 = tmp;
+
+        tmp = m32;
+        m32 = m23;
+        m23 = tmp;
+
+        return this;
+    }
+
+    /**
+     * Transpose the given {@code src} matrix into this matrix.
+     *
+     * @param src source 4x4 matrix
+     * @return this matrix (result) for chaining
+     */
+    public final Matrix4f transpose(final Matrix4f src) {
+        if( src == this ) {
+            return transpose();
+        }
+        m00 = src.m00;
+        m10 = src.m01;
+        m20 = src.m02;
+        m30 = src.m03;
+
+        m01 = src.m10;
+        m11 = src.m11;
+        m21 = src.m12;
+        m31 = src.m13;
+
+        m02 = src.m20;
+        m12 = src.m21;
+        m22 = src.m22;
+        m32 = src.m23;
+
+        m03 = src.m30;
+        m13 = src.m31;
+        m23 = src.m32;
+        m33 = src.m33;
+        return this;
+    }
+
+    /**
+     * Invert this matrix.
+     * @return false if this matrix is singular and inversion not possible, otherwise true
+     */
+    public boolean invert() {
+        final float scale;
+        {
+            float a = Math.abs(m00);
+            float max = a;
+
+            a = Math.abs(m01); if( a > max ) max = a;
+            a = Math.abs(m02); if( a > max ) max = a;
+            a = Math.abs(m03); if( a > max ) max = a;
+
+            a = Math.abs(m10); if( a > max ) max = a;
+            a = Math.abs(m11); if( a > max ) max = a;
+            a = Math.abs(m12); if( a > max ) max = a;
+            a = Math.abs(m13); if( a > max ) max = a;
+
+            a = Math.abs(m20); if( a > max ) max = a;
+            a = Math.abs(m21); if( a > max ) max = a;
+            a = Math.abs(m22); if( a > max ) max = a;
+            a = Math.abs(m23); if( a > max ) max = a;
+
+            a = Math.abs(m30); if( a > max ) max = a;
+            a = Math.abs(m31); if( a > max ) max = a;
+            a = Math.abs(m32); if( a > max ) max = a;
+            a = Math.abs(m33); if( a > max ) max = a;
+
+            if( 0 == max ) {
+                return false;
+            }
+            scale = 1.0f/max;
+        }
+
+        final float a00 = m00*scale;
+        final float a10 = m10*scale;
+        final float a20 = m20*scale;
+        final float a30 = m30*scale;
+
+        final float a01 = m01*scale;
+        final float a11 = m11*scale;
+        final float a21 = m21*scale;
+        final float a31 = m31*scale;
+
+        final float a02 = m02*scale;
+        final float a12 = m12*scale;
+        final float a22 = m22*scale;
+        final float a32 = m32*scale;
+
+        final float a03 = m03*scale;
+        final float a13 = m13*scale;
+        final float a23 = m23*scale;
+        final float a33 = m33*scale;
+
+        final float b00 = + a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+        final float b01 = -( + a10*(a22*a33 - a23*a32) - a12*(a20*a33 - a23*a30) + a13*(a20*a32 - a22*a30));
+        final float b02 = + a10*(a21*a33 - a23*a31) - a11*(a20*a33 - a23*a30) + a13*(a20*a31 - a21*a30);
+        final float b03 = -( + a10*(a21*a32 - a22*a31) - a11*(a20*a32 - a22*a30) + a12*(a20*a31 - a21*a30));
+
+        final float b10 = -( + a01*(a22*a33 - a23*a32) - a02*(a21*a33 - a23*a31) + a03*(a21*a32 - a22*a31));
+        final float b11 = + a00*(a22*a33 - a23*a32) - a02*(a20*a33 - a23*a30) + a03*(a20*a32 - a22*a30);
+        final float b12 = -( + a00*(a21*a33 - a23*a31) - a01*(a20*a33 - a23*a30) + a03*(a20*a31 - a21*a30));
+        final float b13 = + a00*(a21*a32 - a22*a31) - a01*(a20*a32 - a22*a30) + a02*(a20*a31 - a21*a30);
+
+        final float b20 = + a01*(a12*a33 - a13*a32) - a02*(a11*a33 - a13*a31) + a03*(a11*a32 - a12*a31);
+        final float b21 = -( + a00*(a12*a33 - a13*a32) - a02*(a10*a33 - a13*a30) + a03*(a10*a32 - a12*a30));
+        final float b22 = + a00*(a11*a33 - a13*a31) - a01*(a10*a33 - a13*a30) + a03*(a10*a31 - a11*a30);
+        final float b23 = -( + a00*(a11*a32 - a12*a31) - a01*(a10*a32 - a12*a30) + a02*(a10*a31 - a11*a30));
+
+        final float b30 = -( + a01*(a12*a23 - a13*a22) - a02*(a11*a23 - a13*a21) + a03*(a11*a22 - a12*a21));
+        final float b31 = + a00*(a12*a23 - a13*a22) - a02*(a10*a23 - a13*a20) + a03*(a10*a22 - a12*a20);
+        final float b32 = -( + a00*(a11*a23 - a13*a21) - a01*(a10*a23 - a13*a20) + a03*(a10*a21 - a11*a20));
+        final float b33 = + a00*(a11*a22 - a12*a21) - a01*(a10*a22 - a12*a20) + a02*(a10*a21 - a11*a20);
+
+        final float det = (a00*b00 + a01*b01 + a02*b02 + a03*b03) / scale;
+
+        if( 0 == det ) {
+            return false;
+        }
+
+        m00 = b00 / det;
+        m10 = b01 / det;
+        m20 = b02 / det;
+        m30 = b03 / det;
+
+        m01 = b10 / det;
+        m11 = b11 / det;
+        m21 = b12 / det;
+        m31 = b13 / det;
+
+        m02 = b20 / det;
+        m12 = b21 / det;
+        m22 = b22 / det;
+        m32 = b23 / det;
+
+        m03 = b30 / det;
+        m13 = b31 / det;
+        m23 = b32 / det;
+        m33 = b33 / det;
+        return true;
+    }
+
+    /**
+     * Invert the {@code src} matrix values into this matrix
+     * @param src the source matrix, which values are to be inverted
+     * @return false if {@code src} matrix is singular and inversion not possible, otherwise true
+     */
+    public boolean invert(final Matrix4f src) {
+        return load(src).invert();
+    }
+
+    /**
+     * Multiply matrix: [this] = [this] x [b]
+     * <p>
+     * Roughly 15% slower than {@link #mul(Matrix4f, Matrix4f)}
+     * Roughly  3% slower than {@link FloatUtil#multMatrix(float[], float[])}
+     * </p>
+     * @param b 4x4 matrix
+     * @return this matrix for chaining
+     * @see #mul(Matrix4f, Matrix4f)
+     */
+    public final Matrix4f mul(final Matrix4f b) {
+        // return mul(new Matrix4f(this), b); // <- roughly half speed
+        float ai0=m00; // row-0, m[0+0*4]
+        float ai1=m01;
+        float ai2=m02;
+        float ai3=m03;
+        m00 = ai0 * b.m00  +  ai1 * b.m10  +  ai2 * b.m20  +  ai3 * b.m30 ;
+        m01 = ai0 * b.m01  +  ai1 * b.m11  +  ai2 * b.m21  +  ai3 * b.m31 ;
+        m02 = ai0 * b.m02  +  ai1 * b.m12  +  ai2 * b.m22  +  ai3 * b.m32 ;
+        m03 = ai0 * b.m03  +  ai1 * b.m13  +  ai2 * b.m23  +  ai3 * b.m33 ;
+
+        ai0=m10; //row-1, m[1+0*4]
+        ai1=m11;
+        ai2=m12;
+        ai3=m13;
+        m10 = ai0 * b.m00  +  ai1 * b.m10  +  ai2 * b.m20  +  ai3 * b.m30 ;
+        m11 = ai0 * b.m01  +  ai1 * b.m11  +  ai2 * b.m21  +  ai3 * b.m31 ;
+        m12 = ai0 * b.m02  +  ai1 * b.m12  +  ai2 * b.m22  +  ai3 * b.m32 ;
+        m13 = ai0 * b.m03  +  ai1 * b.m13  +  ai2 * b.m23  +  ai3 * b.m33 ;
+
+        ai0=m20; // row-2, m[2+0*4]
+        ai1=m21;
+        ai2=m22;
+        ai3=m23;
+        m20 = ai0 * b.m00  +  ai1 * b.m10  +  ai2 * b.m20  +  ai3 * b.m30 ;
+        m21 = ai0 * b.m01  +  ai1 * b.m11  +  ai2 * b.m21  +  ai3 * b.m31 ;
+        m22 = ai0 * b.m02  +  ai1 * b.m12  +  ai2 * b.m22  +  ai3 * b.m32 ;
+        m23 = ai0 * b.m03  +  ai1 * b.m13  +  ai2 * b.m23  +  ai3 * b.m33 ;
+
+        ai0=m30; // row-3, m[3+0*4]
+        ai1=m31;
+        ai2=m32;
+        ai3=m33;
+        m30 = ai0 * b.m00  +  ai1 * b.m10  +  ai2 * b.m20  +  ai3 * b.m30 ;
+        m31 = ai0 * b.m01  +  ai1 * b.m11  +  ai2 * b.m21  +  ai3 * b.m31 ;
+        m32 = ai0 * b.m02  +  ai1 * b.m12  +  ai2 * b.m22  +  ai3 * b.m32 ;
+        m33 = ai0 * b.m03  +  ai1 * b.m13  +  ai2 * b.m23  +  ai3 * b.m33 ;
+        return this;
+    }
+
+    /**
+     * Multiply matrix: [this] = [a] x [b]
+     * <p>
+     * Roughly 13% faster than {@link #mul(Matrix4f)}
+     * Roughly 11% faster than {@link FloatUtil#multMatrix(float[], float[])}
+     * </p>
+     * @param a 4x4 matrix, can't be this matrix
+     * @param b 4x4 matrix, can't be this matrix
+     * @return this matrix for chaining
+     * @see #mul(Matrix4f)
+     */
+    public final Matrix4f mul(final Matrix4f a, final Matrix4f b) {
+        // row-0, m[0+0*4]
+        m00 = a.m00 * b.m00  +  a.m01 * b.m10  +  a.m02 * b.m20  +  a.m03 * b.m30 ;
+        m01 = a.m00 * b.m01  +  a.m01 * b.m11  +  a.m02 * b.m21  +  a.m03 * b.m31 ;
+        m02 = a.m00 * b.m02  +  a.m01 * b.m12  +  a.m02 * b.m22  +  a.m03 * b.m32 ;
+        m03 = a.m00 * b.m03  +  a.m01 * b.m13  +  a.m02 * b.m23  +  a.m03 * b.m33 ;
+
+        //row-1, m[1+0*4]
+        m10 = a.m10 * b.m00  +  a.m11 * b.m10  +  a.m12 * b.m20  +  a.m13 * b.m30 ;
+        m11 = a.m10 * b.m01  +  a.m11 * b.m11  +  a.m12 * b.m21  +  a.m13 * b.m31 ;
+        m12 = a.m10 * b.m02  +  a.m11 * b.m12  +  a.m12 * b.m22  +  a.m13 * b.m32 ;
+        m13 = a.m10 * b.m03  +  a.m11 * b.m13  +  a.m12 * b.m23  +  a.m13 * b.m33 ;
+
+        // row-2, m[2+0*4]
+        m20 = a.m20 * b.m00  +  a.m21 * b.m10  +  a.m22 * b.m20  +  a.m23 * b.m30 ;
+        m21 = a.m20 * b.m01  +  a.m21 * b.m11  +  a.m22 * b.m21  +  a.m23 * b.m31 ;
+        m22 = a.m20 * b.m02  +  a.m21 * b.m12  +  a.m22 * b.m22  +  a.m23 * b.m32 ;
+        m23 = a.m20 * b.m03  +  a.m21 * b.m13  +  a.m22 * b.m23  +  a.m23 * b.m33 ;
+
+        // row-3, m[3+0*4]
+        m30 = a.m30 * b.m00  +  a.m31 * b.m10  +  a.m32 * b.m20  +  a.m33 * b.m30 ;
+        m31 = a.m30 * b.m01  +  a.m31 * b.m11  +  a.m32 * b.m21  +  a.m33 * b.m31 ;
+        m32 = a.m30 * b.m02  +  a.m31 * b.m12  +  a.m32 * b.m22  +  a.m33 * b.m32 ;
+        m33 = a.m30 * b.m03  +  a.m31 * b.m13  +  a.m32 * b.m23  +  a.m33 * b.m33 ;
+
+        return this;
+    }
+
+    /**
+     * @param v_in 4-component column-vector
+     * @param v_out this * v_in
+     * @returns v_out for chaining
+     */
+    public final float[] mulVec4f(final float[/*4*/] v_in, final float[/*4*/] v_out) {
+        // (one matrix row in column-major order) X (column vector)
+        final float x = v_in[0], y = v_in[1], z = v_in[2], w = v_in[3];
+        v_out[0] = x * m00 + y * m01 + z * m02 + w * m03;
+        v_out[1] = x * m10 + y * m11 + z * m12 + w * m13;
+        v_out[2] = x * m20 + y * m21 + z * m22 + w * m23;
+        v_out[3] = x * m30 + y * m31 + z * m32 + w * m33;
+        return v_out;
+    }
+
+    /**
+     * @param v_in 4-component column-vector
+     * @param v_out this * v_in
+     * @returns v_out for chaining
+     */
+    public final Vec4f mulVec4f(final Vec4f v_in, final Vec4f v_out) {
+        // (one matrix row in column-major order) X (column vector)
+        final float x = v_in.x(), y = v_in.y(), z = v_in.z(), w = v_in.w();
+        v_out.set( x * m00 + y * m01 + z * m02 + w * m03,
+                   x * m10 + y * m11 + z * m12 + w * m13,
+                   x * m20 + y * m21 + z * m22 + w * m23,
+                   x * m30 + y * m31 + z * m32 + w * m33 );
+        return v_out;
+    }
+
+    /**
+     * Affine 3f-vector transformation by 4x4 matrix
+     *
+     * 4x4 matrix multiplication with 3-component vector,
+     * using {@code 1} for for {@code v_in[3]} and dropping {@code v_out[3]},
+     * which shall be {@code 1}.
+     *
+     * @param v_in 3-component column-vector
+     * @param v_out m_in * v_in, 3-component column-vector
+     * @returns v_out for chaining
+     */
+    public final float[] mulVec3f(final float[/*3*/] v_in, final float[/*3*/] v_out) {
+        // (one matrix row in column-major order) X (column vector)
+        final float x = v_in[0], y = v_in[1], z = v_in[2];
+        v_out[0] = x * m00 + y * m01 + z * m02 + 1f * m03;
+        v_out[1] = x * m10 + y * m11 + z * m12 + 1f * m13;
+        v_out[2] = x * m20 + y * m21 + z * m22 + 1f * m23;
+        return v_out;
+    }
+
+    /**
+     * Affine 3f-vector transformation by 4x4 matrix
+     *
+     * 4x4 matrix multiplication with 3-component vector,
+     * using {@code 1} for for {@code v_in.w()} and dropping {@code v_out.w()},
+     * which shall be {@code 1}.
+     *
+     * @param v_in 3-component column-vector {@link Vec3f}
+     * @param v_out m_in * v_in, 3-component column-vector {@link Vec3f}
+     * @returns v_out for chaining
+     */
+    public final Vec3f mulVec3f(final Vec3f v_in, final Vec3f v_out) {
+        // (one matrix row in column-major order) X (column vector)
+        final float x = v_in.x(), y = v_in.y(), z = v_in.z();
+        v_out.set( x * m00 + y * m01 + z * m02 + 1f * m03,
+                   x * m10 + y * m11 + z * m12 + 1f * m13,
+                   x * m20 + y * m21 + z * m22 + 1f * m23 );
+        return v_out;
+    }
+
+    //
+    // Matrix setTo...(), affine + basic
+    //
+
+    /**
+     * Set this matrix to translation.
+     * <pre>
+      Translation matrix (Column Order):
+      1 0 0 0
+      0 1 0 0
+      0 0 1 0
+      x y z 1
+     * </pre>
+     * @param x x-axis translate
+     * @param y y-axis translate
+     * @param z z-axis translate
+     * @return this matrix for chaining
+     */
+    public final Matrix4f setToTranslation(final float x, final float y, final float z) {
+        m00 = m11 = m22 = m33 = 1.0f;
+        m03 = x;
+        m13 = y;
+        m23 = z;
+        m01 = m02 =
+        m10 = m12 =
+        m20 = m21 =
+        m30 = m31 = m32 = 0.0f;
+        return this;
+    }
+
+    /**
+     * Set this matrix to translation.
+     * <pre>
+      Translation matrix (Column Order):
+      1 0 0 0
+      0 1 0 0
+      0 0 1 0
+      x y z 1
+     * </pre>
+     * @param t translate Vec3f
+     * @return this matrix for chaining
+     */
+    public final Matrix4f setToTranslation(final Vec3f t) {
+        return setToTranslation(t.x(), t.y(), t.z());
+    }
+
+    /**
+     * Set this matrix to scale.
+     * <pre>
+      Scale matrix (Any Order):
+      x 0 0 0
+      0 y 0 0
+      0 0 z 0
+      0 0 0 1
+     * </pre>
+     * @param x x-axis scale
+     * @param y y-axis scale
+     * @param z z-axis scale
+     * @return this matrix for chaining
+     */
+    public final Matrix4f setToScale(final float x, final float y, final float z) {
+        m33 = 1.0f;
+        m00 = x;
+        m11 = y;
+        m22 = z;
+        m01 = m02 = m03 =
+        m10 = m12 = m13 =
+        m20 = m21 = m23 =
+        m30 = m31 = m32 = 0.0f;
+        return this;
+    }
+
+    /**
+     * Set this matrix to rotation from the given axis and angle in radians.
+     * <pre>
+        Rotation matrix (Column Order):
+        xx(1-c)+c  xy(1-c)+zs xz(1-c)-ys 0
+        xy(1-c)-zs yy(1-c)+c  yz(1-c)+xs 0
+        xz(1-c)+ys yz(1-c)-xs zz(1-c)+c  0
+        0          0          0          1
+     * </pre>
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
+     * @param ang_rad angle in radians
+     * @param x x of rotation axis
+     * @param y y of rotation axis
+     * @param z z of rotation axis
+     * @return this matrix for chaining
+     */
+    public final Matrix4f setToRotationAxis(final float ang_rad, float x, float y, float z) {
+        final float c = FloatUtil.cos(ang_rad);
+        final float ic= 1.0f - c;
+        final float s = FloatUtil.sin(ang_rad);
+
+        final float[] tmpVec3f = { x, y, z };
+        VectorUtil.normalizeVec3(tmpVec3f);
+        x = tmpVec3f[0]; y = tmpVec3f[1]; z = tmpVec3f[2];
+
+        final float xy = x*y;
+        final float xz = x*z;
+        final float xs = x*s;
+        final float ys = y*s;
+        final float yz = y*z;
+        final float zs = z*s;
+        m00 = x*x*ic+c;
+        m10 = xy*ic+zs;
+        m20 = xz*ic-ys;
+        m30 = 0;
+
+        m01 = xy*ic-zs;
+        m11 = y*y*ic+c;
+        m21 = yz*ic+xs;
+        m31 = 0;
+
+        m02 = xz*ic+ys;
+        m12 = yz*ic-xs;
+        m22 = z*z*ic+c;
+        m32 = 0;
+
+        m03 = 0f;
+        m13 = 0f;
+        m23 = 0f;
+        m33 = 1f;
+
+        return this;
+    }
+
+    /**
+     * Set this matrix to rotation from the given axis and angle in radians.
+     * <pre>
+        Rotation matrix (Column Order):
+        xx(1-c)+c  xy(1-c)+zs xz(1-c)-ys 0
+        xy(1-c)-zs yy(1-c)+c  yz(1-c)+xs 0
+        xz(1-c)+ys yz(1-c)-xs zz(1-c)+c  0
+        0          0          0          1
+     * </pre>
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
+     * @param ang_rad angle in radians
+     * @param axis rotation axis
+     * @return this matrix for chaining
+     */
+    public final Matrix4f setToRotationAxis(final float ang_rad, final Vec3f axis) {
+        return setToRotationAxis(ang_rad, axis.x(), axis.y(), axis.z());
+    }
+
+    /**
+     * Set this matrix to rotation from the given Euler rotation angles in radians.
+     * <p>
+     * The rotations are applied in the given order:
+     * <ul>
+     *  <li>y - heading</li>
+     *  <li>z - attitude</li>
+     *  <li>x - bank</li>
+     * </ul>
+     * </p>
+     * @param bankX the Euler pitch angle in radians. (rotation about the X axis)
+     * @param headingY the Euler yaw angle in radians. (rotation about the Y axis)
+     * @param attitudeZ the Euler roll angle in radians. (rotation about the Z axis)
+     * @return this matrix for chaining
+     * <p>
+     * Implementation does not use Quaternion and hence is exposed to
+     * <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q34">Gimbal-Lock</a>,
+     * consider using {@link #setToRotation(Quaternion)}.
+     * </p>
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q36">Matrix-FAQ Q36</a>
+     * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToMatrix/index.htm">euclideanspace.com-eulerToMatrix</a>
+     * @see #setToRotation(Quaternion)
+     */
+    public Matrix4f setToRotationEuler(final float bankX, final float headingY, final float attitudeZ) {
+        // Assuming the angles are in radians.
+        final float ch = FloatUtil.cos(headingY);
+        final float sh = FloatUtil.sin(headingY);
+        final float ca = FloatUtil.cos(attitudeZ);
+        final float sa = FloatUtil.sin(attitudeZ);
+        final float cb = FloatUtil.cos(bankX);
+        final float sb = FloatUtil.sin(bankX);
+
+        m00 =  ch*ca;
+        m10 =  sa;
+        m20 = -sh*ca;
+        m30 =  0;
+
+        m01 =  sh*sb    - ch*sa*cb;
+        m11 =  ca*cb;
+        m21 =  sh*sa*cb + ch*sb;
+        m31 =  0;
+
+        m02 =  ch*sa*sb + sh*cb;
+        m12 = -ca*sb;
+        m22 = -sh*sa*sb + ch*cb;
+        m32 =  0;
+
+        m03 =  0;
+        m13 =  0;
+        m23 =  0;
+        m33 =  1;
+
+        return this;
+    }
+
+    /**
+     * Set this matrix to rotation using the given Quaternion.
+     * <p>
+     * Implementation Details:
+     * <ul>
+     *   <li> makes identity matrix if {@link #magnitudeSquared()} is {@link FloatUtil#isZero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     *   <li> The fields [m00 .. m22] define the rotation</li>
+     * </ul>
+     * </p>
+     *
+     * @param q the Quaternion representing the rotation
+     * @return this matrix for chaining
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q54">Matrix-FAQ Q54</a>
+     * @see Quaternion#toMatrix(float[], int)
+     * @see #getRotation()
+     */
+    public final Matrix4f setToRotation(final Quaternion q) {
+        // pre-multiply scaled-reciprocal-magnitude to reduce multiplications
+        final float norm = q.magnitudeSquared();
+        if ( FloatUtil.isZero(norm, FloatUtil.EPSILON) ) {
+            // identity matrix -> srecip = 0f
+            loadIdentity();
+            return this;
+        }
+        final float srecip;
+        if ( FloatUtil.isEqual(1f, norm, FloatUtil.EPSILON) ) {
+            srecip = 2f;
+        } else {
+            srecip = 2.0f / norm;
+        }
+
+        final float x = q.x();
+        final float y = q.y();
+        final float z = q.z();
+        final float w = q.w();
+
+        final float xs = srecip * x;
+        final float ys = srecip * y;
+        final float zs = srecip * z;
+
+        final float xx = x  * xs;
+        final float xy = x  * ys;
+        final float xz = x  * zs;
+        final float xw = xs * w;
+        final float yy = y  * ys;
+        final float yz = y  * zs;
+        final float yw = ys * w;
+        final float zz = z  * zs;
+        final float zw = zs * w;
+
+        m00 = 1f - ( yy + zz );
+        m01 =      ( xy - zw );
+        m02 =      ( xz + yw );
+        m03 = 0f;
+
+        m10 =      ( xy + zw );
+        m11 = 1f - ( xx + zz );
+        m12 =      ( yz - xw );
+        m13 = 0f;
+
+        m20 =      ( xz - yw );
+        m21 =      ( yz + xw );
+        m22 = 1f - ( xx + yy );
+        m23 = 0f;
+
+        m30 = m31 = m32 = 0f;
+        m33 = 1f;
+        return this;
+    }
+
+    /**
+     * Returns the rotation [m00 .. m22] fields converted to a Quaternion.
+     * @param res resulting Quaternion
+     * @return the resulting Quaternion for chaining.
+     * @see Quaternion#setFromMatrix(float, float, float, float, float, float, float, float, float)
+     * @see #setToRotation(Quaternion)
+     */
+    public final Quaternion getRotation(final Quaternion res) {
+        res.setFromMatrix(m00, m01, m02, m10, m11, m12, m20, m21, m22);
+        return res;
+    }
+
+    /**
+     * Set this matrix to orthogonal projection.
+     * <pre>
+      Ortho matrix (Column Order):
+      2/dx  0     0    0
+      0     2/dy  0    0
+      0     0     2/dz 0
+      tx    ty    tz   1
+     * </pre>
+     * @param left
+     * @param right
+     * @param bottom
+     * @param top
+     * @param zNear
+     * @param zFar
+     * @return this matrix for chaining
+     */
+    public Matrix4f setToOrtho(final float left, final float right,
+                               final float bottom, final float top,
+                               final float zNear, final float zFar) {
+        {
+            // m00 = m11 = m22 = m33 = 1f;
+            m10 = m20 = m30 = 0f;
+            m01 = m21 = m31 = 0f;
+            m02 = m12 = m32 = 0f;
+            // m03 = m13 = m23 = 0f;
+        }
+        final float dx=right-left;
+        final float dy=top-bottom;
+        final float dz=zFar-zNear;
+        final float tx=-1.0f*(right+left)/dx;
+        final float ty=-1.0f*(top+bottom)/dy;
+        final float tz=-1.0f*(zFar+zNear)/dz;
+
+        m00 =  2.0f/dx;
+        m11 =  2.0f/dy;
+        m22 = -2.0f/dz;
+
+        m03 = tx;
+        m13 = ty;
+        m23 = tz;
+        m33 = 1f;
+
+        return this;
+    }
+
+    /**
+     * Set this matrix to frustum.
+     * <pre>
+      Frustum matrix (Column Order):
+      2*zNear/dx   0          0   0
+      0            2*zNear/dy 0   0
+      A            B          C  -1
+      0            0          D   0
+     * </pre>
+     * @param left
+     * @param right
+     * @param bottom
+     * @param top
+     * @param zNear
+     * @param zFar
+     * @return this matrix for chaining
+     * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
+     *                                  or {@code left == right}, or {@code bottom == top}.
+     */
+    public Matrix4f setToFrustum(final float left, final float right,
+                                 final float bottom, final float top,
+                                 final float zNear, final float zFar) throws IllegalArgumentException {
+        if( zNear <= 0.0f || zFar <= zNear ) {
+            throw new IllegalArgumentException("Requirements zNear > 0 and zFar > zNear, but zNear "+zNear+", zFar "+zFar);
+        }
+        if( left == right || top == bottom) {
+            throw new IllegalArgumentException("GL_INVALID_VALUE: top,bottom and left,right must not be equal");
+        }
+        {
+            // m00 = m11 = m22 = m33 = 1f;
+            m10 = m20 = m30 = 0f;
+            m01 = m21 = m31 = 0f;
+            m03 = m13 = 0f;
+        }
+        final float zNear2 = 2.0f*zNear;
+        final float dx=right-left;
+        final float dy=top-bottom;
+        final float dz=zFar-zNear;
+        final float A=(right+left)/dx;
+        final float B=(top+bottom)/dy;
+        final float C=-1.0f*(zFar+zNear)/dz;
+        final float D=-2.0f*(zFar*zNear)/dz;
+
+        m00 = zNear2/dx;
+        m11 = zNear2/dy;
+
+        m02 = A;
+        m12 = B;
+        m22 = C;
+        m32 = -1.0f;
+
+        m23 = D;
+        m33 = 0f;
+
+        return this;
+    }
+
+    /**
+     * Set this matrix to perspective {@link #setToFrustum(float, float, float, float, float, float) frustum} projection.
+     *
+     * @param fovy_rad angle in radians
+     * @param aspect aspect ratio width / height
+     * @param zNear
+     * @param zFar
+     * @return this matrix for chaining
+     * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
+     * @see #setToFrustum(float, float, float, float, float, float)
+     */
+    public Matrix4f setToPerspective(final float fovy_rad, final float aspect, final float zNear, final float zFar) throws IllegalArgumentException {
+        final float top    =  FloatUtil.tan(fovy_rad/2f) * zNear; // use tangent of half-fov !
+        final float bottom =  -1.0f * top;    //          -1f * fovhvTan.top * zNear
+        final float left   = aspect * bottom; // aspect * -1f * fovhvTan.top * zNear
+        final float right  = aspect * top;    // aspect * fovhvTan.top * zNear
+        return setToFrustum(left, right, bottom, top, zNear, zFar);
+    }
+
+    /**
+     * Set this matrix to perspective {@link #setToFrustum(float, float, float, float, float, float) frustum} projection.
+     *
+     * @param fovhv {@link FovHVHalves} field of view in both directions, may not be centered, either in radians or tangent
+     * @param zNear
+     * @param zFar
+     * @return this matrix for chaining
+     * @throws IllegalArgumentException if {@code zNear <= 0} or {@code zFar <= zNear}
+     * @see #setToFrustum(float, float, float, float, float, float)
+     * @see Frustum#updateByFovDesc(float[], int, boolean, Frustum.FovDesc)
+     */
+    public Matrix4f setToPerspective(final FovHVHalves fovhv, final float zNear, final float zFar) throws IllegalArgumentException {
+        final FovHVHalves fovhvTan = fovhv.toTangents();  // use tangent of half-fov !
+        final float top    =         fovhvTan.top    * zNear;
+        final float bottom = -1.0f * fovhvTan.bottom * zNear;
+        final float left   = -1.0f * fovhvTan.left   * zNear;
+        final float right  =         fovhvTan.right  * zNear;
+        return setToFrustum(left, right, bottom, top, zNear, zFar);
+    }
+
+    /**
+     * 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 updateFrustumPlanes(final Frustum frustum) {
+        // Left:   a = m41 + m11, b = m42 + m12, c = m43 + m13, d = m44 + m14  - [1..4] column-major
+        // Left:   a = m30 + m00, b = m31 + m01, c = m32 + m02, d = m33 + m03  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.LEFT];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 + m00,
+                     m31 + m01,
+                     m32 + m02 );
+            p.d    = m33 + m03;
+        }
+
+        // Right:  a = m41 - m11, b = m42 - m12, c = m43 - m13, d = m44 - m14  - [1..4] column-major
+        // Right:  a = m30 - m00, b = m31 - m01, c = m32 - m02, d = m33 - m03  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.RIGHT];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 - m00,
+                     m31 - m01,
+                     m32 - m02 );
+            p.d    = m33 - m03;
+        }
+
+        // Bottom: a = m41m21, b = m42m22, c = m43m23, d = m44m24  - [1..4] column-major
+        // Bottom: a = m30m10, b = m31m11, c = m32m12, d = m33m13  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.BOTTOM];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 + m10,
+                     m31 + m11,
+                     m32 + m12 );
+            p.d    = m33 + m13;
+        }
+
+        // Top:   a = m41 - m21, b = m42 - m22, c = m43 - m23, d = m44 - m24  - [1..4] column-major
+        // Top:   a = m30 - m10, b = m31 - m11, c = m32 - m12, d = m33 - m13  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.TOP];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 - m10,
+                     m31 - m11,
+                     m32 - m12 );
+            p.d    = m33 - m13;
+        }
+
+        // Near:  a = m41m31, b = m42m32, c = m43m33, d = m44m34  - [1..4] column-major
+        // Near:  a = m30m20, b = m31m21, c = m32m22, d = m33m23  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.NEAR];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 + m20,
+                     m31 + m21,
+                     m32 + m22 );
+            p.d    = m33 + m23;
+        }
+
+        // Far:   a = m41 - m31, b = m42 - m32, c = m43 - m33, d = m44 - m34  - [1..4] column-major
+        // Far:   a = m30 - m20, b = m31 - m21, c = m32m22, d = m33m23  - [0..3] column-major
+        {
+            final Frustum.Plane p = frustum.getPlanes()[Frustum.FAR];
+            final Vec3f p_n = p.n;
+            p_n.set( m30 - m20,
+                     m31 - m21,
+                     m32 - m22 );
+            p.d    = m33 - m23;
+        }
+
+        // Normalize all planes
+        for (int i = 0; i < 6; ++i) {
+            final Plane p = frustum.getPlanes()[i];
+            final Vec3f p_n = p.n;
+            final float invLen = 1f / p_n.length();
+            p_n.scale(invLen);
+            p.d *= invLen;
+        }
+    }
+
+    /**
+     * Make given matrix the <i>look-at</i> matrix based on given parameters.
+     * <p>
+     * Consist out of two matrix multiplications:
+     * <pre>
+     *   <b>R</b> = <b>L</b> x <b>T</b>,
+     *   with <b>L</b> for <i>look-at</i> matrix and
+     *        <b>T</b> for eye translation.
+     *
+     *   Result <b>R</b> can be utilized for <i>projection or modelview</i> multiplication, i.e.
+     *          <b>M</b> = <b>M</b> x <b>R</b>,
+     *          with <b>M</b> being the <i>projection or modelview</i> matrix.
+     * </pre>
+     * </p>
+     * @param eye 3 component eye vector
+     * @param center 3 component center vector
+     * @param up 3 component up vector
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public Matrix4f setToLookAt(final Vec3f eye, final Vec3f center, final Vec3f up, final Matrix4f tmp) {
+        // normalized forward!
+        final Vec3f fwd = new Vec3f( center.x() - eye.x(),
+                                     center.y() - eye.y(),
+                                     center.z() - eye.z() ).normalize();
+
+        /* Side = forward x up, normalized */
+        final Vec3f side = fwd.cross(up).normalize();
+
+        /* Recompute up as: up = side x forward */
+        final Vec3f up2 = side.cross(fwd);
+
+        m00 = side.x();
+        m10 = up2.x();
+        m20 = -fwd.x();
+        m30 = 0;
+
+        m01 = side.y();
+        m11 = up2.y();
+        m21 = -fwd.y();
+        m31 = 0;
+
+        m02 = side.z();
+        m12 = up2.z();
+        m22 = -fwd.z();
+        m32 = 0;
+
+        m03 = 0;
+        m13 = 0;
+        m23 = 0;
+        m33 = 1;
+
+        return mul( tmp.setToTranslation( -eye.x(), -eye.y(), -eye.z() ) );
+    }
+
+    //
+    // Matrix affine operations using setTo..()
+    //
+
+    /**
+     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(float, float, float, float) axis-rotation matrix}.
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
+     * @param angrad angle in radians
+     * @param x x of rotation axis
+     * @param y y of rotation axis
+     * @param z z of rotation axis
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f rotate(final float ang_rad, final float x, final float y, final float z, final Matrix4f tmp) {
+        return mul( tmp.setToRotationAxis(ang_rad, x, y, z) );
+    }
+
+    /**
+     * Rotate this matrix about give axis and angle in radians, i.e. multiply by {@link #setToRotationAxis(float, Vec3f) axis-rotation matrix}.
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38">Matrix-FAQ Q38</a>
+     * @param angrad angle in radians
+     * @param axis rotation axis
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f rotate(final float ang_rad, final Vec3f axis, final Matrix4f tmp) {
+        return mul( tmp.setToRotationAxis(ang_rad, axis) );
+    }
+
+    /**
+     * Rotate this matrix with the given {@link Quaternion}, i.e. multiply by {@link #setToRotation(Quaternion) Quaternion's rotation matrix}.
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f rotate(final Quaternion quat, final Matrix4f tmp) {
+        return mul( tmp.setToRotation(quat) );
+    }
+
+    /**
+     * Translate this matrix, i.e. multiply by {@link #setToTranslation(float, float, float) translation matrix}.
+     * @param x x translation
+     * @param y y translation
+     * @param z z translation
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f translate(final float x, final float y, final float z, final Matrix4f tmp) {
+        return mul( tmp.setToTranslation(x, y, z) );
+    }
+
+    /**
+     * Translate this matrix, i.e. multiply by {@link #setToTranslation(Vec3f) translation matrix}.
+     * @param t translation Vec3f
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f translate(final Vec3f t, final Matrix4f tmp) {
+        return mul( tmp.setToTranslation(t) );
+    }
+
+    /**
+     * Scale this matrix, i.e. multiply by {@link #setToScale(float, float, float) scale matrix}.
+     * @param x x scale
+     * @param y y scale
+     * @param z z scale
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f scale(final float x, final float y, final float z, final Matrix4f tmp) {
+        return mul( tmp.setToScale(x, y, z) );
+    }
+
+    /**
+     * Scale this matrix, i.e. multiply by {@link #setToScale(float, float, float) scale matrix}.
+     * @param s scale for x-, y- and z-axis
+     * @param tmp temporary Matrix4f used for multiplication
+     * @return this matrix for chaining
+     */
+    public final Matrix4f scale(final float s, final Matrix4f tmp) {
+        return mul( tmp.setToScale(s, s, s) );
+    }
+
+    //
+    // Matrix Stack
+    //
+
+    /**
+     * Push the matrix to it's stack, while preserving this matrix values.
+     * @see #pop()
+     */
+    public final void push() {
+        stack.push(this);
+    }
+
+    /**
+     * Pop the current matrix from it's stack, replacing this matrix values.
+     * @see #push()
+     */
+    public final void pop() {
+        stack.pop(this);
+    }
+
+    //
+    // equals
+    //
+
+    /**
+     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @param epsilon consider using {@link FloatUtil#EPSILON}
+     * @return true if all components differ less than {@code epsilon}, otherwise false.
+     */
+    public boolean isEqual(final Matrix4f o, final float epsilon) {
+        if( this == o ) {
+            return true;
+        } else {
+            return FloatUtil.isEqual(m00, o.m00, epsilon) &&
+                   FloatUtil.isEqual(m01, o.m01, epsilon) &&
+                   FloatUtil.isEqual(m02, o.m02, epsilon) &&
+                   FloatUtil.isEqual(m03, o.m03, epsilon) &&
+                   FloatUtil.isEqual(m10, o.m10, epsilon) &&
+                   FloatUtil.isEqual(m11, o.m11, epsilon) &&
+                   FloatUtil.isEqual(m12, o.m12, epsilon) &&
+                   FloatUtil.isEqual(m13, o.m13, epsilon) &&
+                   FloatUtil.isEqual(m20, o.m20, epsilon) &&
+                   FloatUtil.isEqual(m21, o.m21, epsilon) &&
+                   FloatUtil.isEqual(m22, o.m22, epsilon) &&
+                   FloatUtil.isEqual(m23, o.m23, epsilon) &&
+                   FloatUtil.isEqual(m30, o.m30, epsilon) &&
+                   FloatUtil.isEqual(m31, o.m31, epsilon) &&
+                   FloatUtil.isEqual(m32, o.m32, epsilon) &&
+                   FloatUtil.isEqual(m33, o.m33, epsilon);
+        }
+    }
+
+    /**
+     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @return true if all components differ less than {@link FloatUtil#EPSILON}, otherwise false.
+     */
+    public boolean isEqual(final Matrix4f o) {
+        return isEqual(o, FloatUtil.EPSILON);
+    }
+
+    @Override
+    public boolean equals(final Object o) {
+        if( o instanceof Matrix4f ) {
+            return isEqual((Matrix4f)o, FloatUtil.EPSILON);
+        } else {
+            return false;
+        }
+    }
+
+    //
+    // Static multi Matrix ops
+    //
+
+    /**
+     * Map object coordinates to window coordinates.
+     * <p>
+     * Traditional <code>gluProject</code> implementation.
+     * </p>
+     *
+     * @param obj object position, 3 component vector
+     * @param mMv modelview matrix
+     * @param mP projection matrix
+     * @param viewport 4 component viewport vector
+     * @param winPos 3 component window coordinate, the result
+     * @return true if successful, otherwise false (z is 1)
+     */
+    public static boolean mapObjToWin(final Vec3f obj, final Matrix4f mMv, final Matrix4f mP,
+                                      final int[] viewport, final float[] winPos)
+    {
+        final Vec4f vec4Tmp1 = new Vec4f(obj, 1f);
+
+        // vec4Tmp2 = Mv * o
+        // rawWinPos = P  * vec4Tmp2
+        // rawWinPos = P * ( Mv * o )
+        // rawWinPos = P * Mv * o
+        final Vec4f vec4Tmp2 = mMv.mulVec4f(vec4Tmp1, new Vec4f());
+        final Vec4f rawWinPos = mP.mulVec4f(vec4Tmp2, vec4Tmp1);
+
+        if (rawWinPos.w() == 0.0f) {
+            return false;
+        }
+
+        final float s = ( 1.0f / rawWinPos.w() ) * 0.5f;
+
+        // Map x, y and z to range 0-1 (w is ignored)
+        rawWinPos.scale(s).add(0.5f, 0.5f, 0.5f, 0f);
+
+        // Map x,y to viewport
+        winPos[0] = rawWinPos.x() * viewport[2] + viewport[0];
+        winPos[1] = rawWinPos.y() * viewport[3] + viewport[1];
+        winPos[2] = rawWinPos.z();
+
+        return true;
+    }
+
+    /**
+     * Map object coordinates to window coordinates.
+     * <p>
+     * Traditional <code>gluProject</code> implementation.
+     * </p>
+     *
+     * @param obj object position, 3 component vector
+     * @param mPMv [projection] x [modelview] matrix, i.e. P x Mv
+     * @param viewport 4 component viewport vector
+     * @param winPos 3 component window coordinate, the result
+     * @return true if successful, otherwise false (z is 1)
+     */
+    public static boolean mapObjToWin(final Vec3f obj, final Matrix4f mPMv,
+                                      final int[] viewport, final float[] winPos)
+    {
+        final Vec4f vec4Tmp2 = new Vec4f(obj, 1f);
+
+        // rawWinPos = P * Mv * o
+        final Vec4f rawWinPos = mPMv.mulVec4f(vec4Tmp2, new Vec4f());
+
+        if (rawWinPos.w() == 0.0f) {
+            return false;
+        }
+
+        final float s = ( 1.0f / rawWinPos.w() ) * 0.5f;
+
+        // Map x, y and z to range 0-1 (w is ignored)
+        rawWinPos.scale(s).add(0.5f, 0.5f, 0.5f, 0f);
+
+        // Map x,y to viewport
+        winPos[0] = rawWinPos.x() * viewport[2] + viewport[0];
+        winPos[1] = rawWinPos.y() * viewport[3] + viewport[1];
+        winPos[2] = rawWinPos.z();
+
+        return true;
+    }
+
+    /**
+     * Map window coordinates to object coordinates.
+     * <p>
+     * Traditional <code>gluUnProject</code> implementation.
+     * </p>
+     *
+     * @param winx
+     * @param winy
+     * @param winz
+     * @param mMv 4x4 modelview matrix
+     * @param mP 4x4 projection matrix
+     * @param viewport 4 component viewport vector
+     * @param objPos 3 component object coordinate, the result
+     * @param mat4Tmp 16 component matrix for temp storage
+     * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
+     */
+    public static boolean mapWinToObj(final float winx, final float winy, final float winz,
+                                      final Matrix4f mMv, final Matrix4f mP,
+                                      final int[] viewport,
+                                      final Vec3f objPos,
+                                      final Matrix4f mat4Tmp)
+    {
+        // invPMv = Inv(P x Mv)
+        final Matrix4f invPMv = mat4Tmp.mul(mP, mMv);
+        if( !invPMv.invert() ) {
+            return false;
+        }
+
+        final Vec4f winPos = new Vec4f(winx, winy, winz, 1f);
+
+        // Map x and y from window coordinates
+        winPos.add(-viewport[0], -viewport[1], 0f, 0f).scale(1f/viewport[2], 1f/viewport[3], 1f, 1f);
+
+        // Map to range -1 to 1
+        winPos.scale(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
+
+        // rawObjPos = Inv(P x Mv) *  winPos
+        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+
+        if ( rawObjPos.w() == 0.0f ) {
+            return false;
+        }
+        objPos.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
+
+        return true;
+    }
+
+    /**
+     * Map window coordinates to object coordinates.
+     * <p>
+     * Traditional <code>gluUnProject</code> implementation.
+     * </p>
+     *
+     * @param winx
+     * @param winy
+     * @param winz
+     * @param invPMv inverse [projection] x [modelview] matrix, i.e. Inv(P x Mv)
+     * @param viewport 4 component viewport vector
+     * @param objPos 3 component object coordinate, the result
+     * @param mat4Tmp 16 component matrix for temp storage
+     * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
+     */
+    public static boolean mapWinToObj(final float winx, final float winy, final float winz,
+                                      final Matrix4f invPMv,
+                                      final int[] viewport,
+                                      final Vec3f objPos,
+                                      final Matrix4f mat4Tmp)
+    {
+        final Vec4f winPos = new Vec4f(winx, winy, winz, 1f);
+
+        // Map x and y from window coordinates
+        winPos.add(-viewport[0], -viewport[1], 0f, 0f).scale(1f/viewport[2], 1f/viewport[3], 1f, 1f);
+
+        // Map to range -1 to 1
+        winPos.scale(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
+
+        // rawObjPos = Inv(P x Mv) *  winPos
+        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+
+        if ( rawObjPos.w() == 0.0f ) {
+            return false;
+        }
+        objPos.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
+
+        return true;
+    }
+
+    /**
+     * Map two window coordinates to two object coordinates,
+     * distinguished by their z component.
+     * <p>
+     * Traditional <code>gluUnProject</code> implementation.
+     * </p>
+     *
+     * @param winx
+     * @param winy
+     * @param winz1
+     * @param winz2
+     * @param invPMv inverse [projection] x [modelview] matrix, i.e. Inv(P x Mv)
+     * @param viewport 4 component viewport vector
+     * @param objPos1 3 component object coordinate, the result
+     * @param mat4Tmp 16 component matrix for temp storage
+     * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
+     */
+    public static boolean mapWinToObj(final float winx, final float winy, final float winz1, final float winz2,
+                                      final Matrix4f invPMv,
+                                      final int[] viewport,
+                                      final Vec3f objPos1, final Vec3f objPos2,
+                                      final Matrix4f mat4Tmp)
+    {
+        final Vec4f winPos = new Vec4f(winx, winy, winz1, 1f);
+
+        // Map x and y from window coordinates
+        winPos.add(-viewport[0], -viewport[1], 0f, 0f).scale(1f/viewport[2], 1f/viewport[3], 1f, 1f);
+
+        // Map to range -1 to 1
+        winPos.scale(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
+
+        // rawObjPos = Inv(P x Mv) *  winPos1
+        final Vec4f rawObjPos = invPMv.mulVec4f(winPos, new Vec4f());
+
+        if ( rawObjPos.w() == 0.0f ) {
+            return false;
+        }
+        objPos1.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
+
+        //
+        // winz2
+        //
+        // Map Z to range -1 to 1
+        winPos.setZ( winz2 * 2f - 1f );
+
+        // rawObjPos = Inv(P x Mv) *  winPos2
+        invPMv.mulVec4f(winPos, rawObjPos);
+
+        if ( rawObjPos.w() == 0.0f ) {
+            return false;
+        }
+        objPos2.set( rawObjPos.scale( 1f / rawObjPos.w() ) );
+
+        return true;
+    }
+
+    /**
+     * Map window coordinates to object coordinates.
+     * <p>
+     * Traditional <code>gluUnProject4</code> implementation.
+     * </p>
+     *
+     * @param winx
+     * @param winy
+     * @param winz
+     * @param clipw
+     * @param mMv 4x4 modelview matrix
+     * @param mP 4x4 projection matrix
+     * @param viewport 4 component viewport vector
+     * @param near
+     * @param far
+     * @param obj_pos 4 component object coordinate, the result
+     * @param mat4Tmp 16 component matrix for temp storage
+     * @return true if successful, otherwise false (failed to invert matrix, or becomes infinity due to zero z)
+     */
+    public static boolean mapWinToObj4(final float winx, final float winy, final float winz, final float clipw,
+                                       final Matrix4f mMv, final Matrix4f mP,
+                                       final int[] viewport,
+                                       final float near, final float far,
+                                       final Vec4f objPos,
+                                       final Matrix4f mat4Tmp)
+    {
+        // invPMv = Inv(P x Mv)
+        final Matrix4f invPMv = mat4Tmp.mul(mP, mMv);
+        if( !invPMv.invert() ) {
+            return false;
+        }
+
+        final Vec4f winPos = new Vec4f(winx, winy, winz, clipw);
+
+        // Map x and y from window coordinates
+        winPos.add(-viewport[0], -viewport[1], -near, 0f).scale(1f/viewport[2], 1f/viewport[3], 1f/(far-near), 1f);
+
+        // Map to range -1 to 1
+        winPos.scale(2f, 2f, 2f, 1f).add(-1f, -1f, -1f, 0f);
+
+        // objPos = Inv(P x Mv) *  winPos
+        invPMv.mulVec4f(winPos, objPos);
+
+        if ( objPos.w() == 0.0f ) {
+            return false;
+        }
+        return true;
+    }
+
+    /**
+     * Map two window coordinates w/ shared X/Y and distinctive Z
+     * to a {@link Ray}. The resulting {@link Ray} maybe used for <i>picking</i>
+     * using a {@link AABBox#getRayIntersection(Ray, float[]) bounding box}.
+     * <p>
+     * Notes for picking <i>winz0</i> and <i>winz1</i>:
+     * <ul>
+     *   <li>see {@link FloatUtil#getZBufferEpsilon(int, float, float)}</li>
+     *   <li>see {@link FloatUtil#getZBufferValue(int, float, float, float)}</li>
+     *   <li>see {@link FloatUtil#getOrthoWinZ(float, float, float)}</li>
+     * </ul>
+     * </p>
+     * @param winx
+     * @param winy
+     * @param winz0
+     * @param winz1
+     * @param mMv 4x4 modelview matrix
+     * @param mP 4x4 projection matrix
+     * @param viewport 4 component viewport vector
+     * @param ray storage for the resulting {@link Ray}
+     * @param mat4Tmp1 16 component matrix for temp storage
+     * @param mat4Tmp2 16 component matrix for temp storage
+     * @return true if successful, otherwise false (failed to invert matrix, or becomes z is infinity)
+     */
+    public static boolean mapWinToRay(final float winx, final float winy, final float winz0, final float winz1,
+                                      final Matrix4f mMv,
+                                      final Matrix4f mP,
+                                      final int[] viewport,
+                                      final Ray ray,
+                                      final Matrix4f mat4Tmp1, final Matrix4f mat4Tmp2) {
+        // invPMv = Inv(P x Mv)
+        final Matrix4f invPMv = mat4Tmp1.mul(mP, mMv);
+        if( !invPMv.invert() ) {
+            return false;
+        }
+
+        if( mapWinToObj(winx, winy, winz0, winz1, invPMv, viewport,
+                        ray.orig, ray.dir, mat4Tmp2) ) {
+            ray.dir.sub(ray.orig).normalize();
+            return true;
+        } else {
+            return false;
+        }
+    }
+
+    //
+    // String and internals
+    //
+
+    /**
+     * @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}
+     * @return matrix string representation
+     */
+    public StringBuilder toString(final StringBuilder sb, final String rowPrefix, final String f) {
+        final float[] tmp = new float[16];
+        this.get(tmp);
+        return FloatUtil.matrixToString(sb, rowPrefix, f,tmp, 0, 4, 4, false /* rowMajorOrder */);
+    }
+
+    @Override
+    public String toString() {
+        return toString(null, null, "%10.5f").toString();
+    }
+
+    private float m00, m10, m20, m30;
+    private float m01, m11, m21, m31;
+    private float m02, m12, m22, m32;
+    private float m03, m13, m23, m33;
+
+    final Stack stack = new Stack(0, 16*16); // start w/ zero size, growSize is half GL-min size (32)
+
+    private static class Stack {
+        private int position;
+        private float[] buffer;
+        private final int growSize;
+
+        /**
+         * @param initialSize initial size
+         * @param growSize grow size if {@link #position()} is reached, maybe <code>0</code>
+         *        in which case an {@link IndexOutOfBoundsException} is thrown.
+         */
+        public Stack(final int initialSize, final int growSize) {
+            this.position = 0;
+            this.growSize = growSize;
+            this.buffer = new float[initialSize];
+        }
+
+        private final void growIfNecessary(final int length) throws IndexOutOfBoundsException {
+            if( position + length > buffer.length ) {
+                if( 0 >= growSize ) {
+                    throw new IndexOutOfBoundsException("Out of fixed stack size: "+this);
+                }
+                final float[] newBuffer =
+                        new float[buffer.length + growSize];
+                System.arraycopy(buffer, 0, newBuffer, 0, position);
+                buffer = newBuffer;
+            }
+        }
+
+        public final Matrix4f push(final Matrix4f src) throws IndexOutOfBoundsException {
+            growIfNecessary(16);
+            src.get(buffer, position);
+            position += 16;
+            return src;
+        }
+
+        public final Matrix4f pop(final Matrix4f dest) throws IndexOutOfBoundsException {
+            position -= 16;
+            dest.load(buffer, position);
+            return dest;
+        }
+    }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java b/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java
index 849477f54..2bb0f96c6 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/Quaternion.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2010 JogAmp Community. All rights reserved.
+ * Copyright 2010-2023 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:
@@ -99,7 +99,7 @@ public class Quaternion {
         return FloatUtil.sqrt(magnitudeSQ);
     }
 
-    public final float getW() {
+    public final float w() {
         return w;
     }
 
@@ -107,7 +107,7 @@ public class Quaternion {
         this.w = w;
     }
 
-    public final float getX() {
+    public final float x() {
         return x;
     }
 
@@ -115,7 +115,7 @@ public class Quaternion {
         this.x = x;
     }
 
-    public final float getY() {
+    public final float y() {
         return y;
     }
 
@@ -123,7 +123,7 @@ public class Quaternion {
         this.y = y;
     }
 
-    public final float getZ() {
+    public final float z() {
         return z;
     }
 
@@ -142,15 +142,15 @@ public class Quaternion {
      * Returns the dot product of this quaternion with the given quaternion
      */
     public final float dot(final Quaternion quat) {
-        return dot(quat.getX(), quat.getY(), quat.getZ(), quat.getW());
+        return dot(quat.x(), quat.y(), quat.z(), quat.w());
     }
 
     /**
      * Returns <code>true</code> if this quaternion has identity.
      * <p>
      * Implementation uses {@link FloatUtil#EPSILON epsilon} to compare
-     * {@link #getW() W} {@link FloatUtil#isEqual(float, float, float) against 1f} and
-     * {@link #getX() X}, {@link #getY() Y} and {@link #getZ() Z}
+     * {@link #w() W} {@link FloatUtil#isEqual(float, float, float) against 1f} and
+     * {@link #x() X}, {@link #y() Y} and {@link #z() Z}
      * {@link FloatUtil#isZero(float, float) against zero}.
      * </p>
      */
@@ -404,12 +404,12 @@ public class Quaternion {
      * </ul>
      * </p>
      * For details see {@link #rotateByEuler(float, float, float)}.
-     * @param angradXYZ euler angel array in radians
+     * @param angradXYZ euler angle array in radians
      * @return this quaternion for chaining.
      * @see #rotateByEuler(float, float, float)
      */
-    public final Quaternion rotateByEuler(final float[] angradXYZ) {
-        return rotateByEuler(angradXYZ[0], angradXYZ[1], angradXYZ[2]);
+    public final Quaternion rotateByEuler(final Vec3f angradXYZ) {
+        return rotateByEuler(angradXYZ.x(), angradXYZ.y(), angradXYZ.z());
     }
 
     /**
@@ -450,48 +450,42 @@ public class Quaternion {
 
     /***
      * Rotate the given vector by this quaternion
+     * @param vecIn vector to be rotated
+     * @param vecOut result storage for rotated vector, maybe equal to vecIn for in-place rotation
      *
-     * @param vecOut result float[3] storage for rotated vector, maybe equal to vecIn for in-place rotation
-     * @param vecOutOffset offset in result storage
-     * @param vecIn float[3] vector to be rotated
-     * @param vecInOffset offset in vecIn
      * @return the given vecOut store for chaining
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q63">Matrix-FAQ Q63</a>
      */
-    public final float[] rotateVector(final float[] vecOut, final int vecOutOffset, final float[] vecIn, final int vecInOffset) {
-        if ( VectorUtil.isVec3Zero(vecIn, vecInOffset, FloatUtil.EPSILON) ) {
-            vecOut[0+vecOutOffset] = 0f;
-            vecOut[1+vecOutOffset] = 0f;
-            vecOut[2+vecOutOffset] = 0f;
+    public final Vec3f rotateVector(final Vec3f vecIn, final Vec3f vecOut) {
+        if( vecIn.isZero() ) {
+            vecOut.set(0, 0, 0);
         } else {
-            final float vecX = vecIn[0+vecInOffset];
-            final float vecY = vecIn[1+vecInOffset];
-            final float vecZ = vecIn[2+vecInOffset];
+            final float vecX = vecIn.x();
+            final float vecY = vecIn.y();
+            final float vecZ = vecIn.z();
             final float x_x = x*x;
             final float y_y = y*y;
             final float z_z = z*z;
             final float w_w = w*w;
 
-            vecOut[0+vecOutOffset] =   w_w * vecX
-                                     + x_x * vecX
-                                     - z_z * vecX
-                                     - y_y * vecX
-                                     + 2f * ( y*w*vecZ - z*w*vecY + y*x*vecY + z*x*vecZ );
+            vecOut.setX(   w_w * vecX
+                         + x_x * vecX
+                         - z_z * vecX
+                         - y_y * vecX
+                         + 2f * ( y*w*vecZ - z*w*vecY + y*x*vecY + z*x*vecZ ) );
                                      ;
 
-            vecOut[1+vecOutOffset] =   y_y * vecY
-                                     - z_z * vecY
-                                     + w_w * vecY
-                                     - x_x * vecY
-                                     + 2f * ( x*y*vecX + z*y*vecZ + w*z*vecX - x*w*vecZ );
-                                     ;
-
-            vecOut[2+vecOutOffset] =   z_z * vecZ
-                                     - y_y * vecZ
-                                     - x_x * vecZ
-                                     + w_w * vecZ
-                                     + 2f * ( x*z*vecX + y*z*vecY - w*y*vecX + w*x*vecY );
-                                     ;
+            vecOut.setY(   y_y * vecY
+                         - z_z * vecY
+                         + w_w * vecY
+                         - x_x * vecY
+                         + 2f * ( x*y*vecX + z*y*vecZ + w*z*vecX - x*w*vecZ ) );;
+
+            vecOut.setZ(   z_z * vecZ
+                         - y_y * vecZ
+                         - x_x * vecZ
+                         + w_w * vecZ
+                         + 2f * ( x*z*vecX + y*z*vecY - w*y*vecX + w*x*vecY ) );
         }
         return vecOut;
     }
@@ -593,21 +587,19 @@ public class Quaternion {
      * @return this quaternion for chaining.
      * @see <a href="http://www.euclideanspace.com/maths/algebra/vectors/lookat/index.htm">euclideanspace.com-LookUp</a>
      */
-    public Quaternion setLookAt(final float[] directionIn, final float[] upIn,
-                                final float[] xAxisOut, final float[] yAxisOut, final float[] zAxisOut) {
+    public Quaternion setLookAt(final Vec3f directionIn, final Vec3f upIn,
+                                final Vec3f xAxisOut, final Vec3f yAxisOut, final Vec3f zAxisOut) {
         // Z = norm(dir)
-        VectorUtil.normalizeVec3(zAxisOut, directionIn);
+        zAxisOut.set(directionIn).normalize();
 
         // X = upIn x Z
         //     (borrow yAxisOut for upNorm)
-        VectorUtil.normalizeVec3(yAxisOut, upIn);
-        VectorUtil.crossVec3(xAxisOut, yAxisOut, zAxisOut);
-        VectorUtil.normalizeVec3(xAxisOut);
+        yAxisOut.set(upIn).normalize();
+        xAxisOut.cross(yAxisOut, zAxisOut).normalize();
 
         // Y = Z x X
         //
-        VectorUtil.crossVec3(yAxisOut, zAxisOut, xAxisOut);
-        VectorUtil.normalizeVec3(yAxisOut);
+        yAxisOut.cross(zAxisOut, xAxisOut).normalize();
 
         /**
             final float m00 = xAxisOut[0];
@@ -642,42 +634,42 @@ public class Quaternion {
      * </p>
      * @param v1 not normalized
      * @param v2 not normalized
-     * @param tmpPivotVec float[3] temp storage for cross product
-     * @param tmpNormalVec float[3] temp storage to normalize vector
+     * @param tmpPivotVec temp storage for cross product
+     * @param tmpNormalVec temp storage to normalize vector
      * @return this quaternion for chaining.
      */
-    public final Quaternion setFromVectors(final float[] v1, final float[] v2, final float[] tmpPivotVec, final float[] tmpNormalVec) {
-        final float factor = VectorUtil.normVec3(v1) * VectorUtil.normVec3(v2);
+    public final Quaternion setFromVectors(final Vec3f v1, final Vec3f v2, final Vec3f tmpPivotVec, final Vec3f tmpNormalVec) {
+        final float factor = v1.length() * v2.length();
         if ( FloatUtil.isZero(factor, FloatUtil.EPSILON ) ) {
             return setIdentity();
         } else {
-            final float dot = VectorUtil.dotVec3(v1, v2) / factor; // normalize
+            final float dot = v1.dot(v2) / factor; // normalize
             final float theta = FloatUtil.acos(Math.max(-1.0f, Math.min(dot, 1.0f))); // clipping [-1..1]
 
-            VectorUtil.crossVec3(tmpPivotVec, v1, v2);
+            tmpPivotVec.cross(v1, v2);
 
-            if ( dot < 0.0f && FloatUtil.isZero( VectorUtil.normVec3(tmpPivotVec), FloatUtil.EPSILON ) ) {
+            if ( dot < 0.0f && FloatUtil.isZero( tmpPivotVec.length(), FloatUtil.EPSILON ) ) {
                 // Vectors parallel and opposite direction, therefore a rotation of 180 degrees about any vector
                 // perpendicular to this vector will rotate vector a onto vector b.
                 //
                 // The following guarantees the dot-product will be 0.0.
                 int dominantIndex;
-                if (Math.abs(v1[0]) > Math.abs(v1[1])) {
-                    if (Math.abs(v1[0]) > Math.abs(v1[2])) {
+                if (Math.abs(v1.x()) > Math.abs(v1.y())) {
+                    if (Math.abs(v1.x()) > Math.abs(v1.z())) {
                         dominantIndex = 0;
                     } else {
                         dominantIndex = 2;
                     }
                 } else {
-                    if (Math.abs(v1[1]) > Math.abs(v1[2])) {
+                    if (Math.abs(v1.y()) > Math.abs(v1.z())) {
                         dominantIndex = 1;
                     } else {
                         dominantIndex = 2;
                     }
                 }
-                tmpPivotVec[dominantIndex]           = -v1[(dominantIndex + 1) % 3];
-                tmpPivotVec[(dominantIndex + 1) % 3] =  v1[dominantIndex];
-                tmpPivotVec[(dominantIndex + 2) % 3] =  0f;
+                tmpPivotVec.set( dominantIndex,           -v1.get( (dominantIndex + 1) % 3 ) );
+                tmpPivotVec.set( (dominantIndex + 1) % 3,  v1.get( dominantIndex ) );
+                tmpPivotVec.set( (dominantIndex + 2) % 3,  0f );
             }
             return setFromAngleAxis(theta, tmpPivotVec, tmpNormalVec);
         }
@@ -698,41 +690,41 @@ public class Quaternion {
      * </p>
      * @param v1 normalized
      * @param v2 normalized
-     * @param tmpPivotVec float[3] temp storage for cross product
+     * @param tmpPivotVec temp storage for cross product
      * @return this quaternion for chaining.
      */
-    public final Quaternion setFromNormalVectors(final float[] v1, final float[] v2, final float[] tmpPivotVec) {
-        final float factor = VectorUtil.normVec3(v1) * VectorUtil.normVec3(v2);
+    public final Quaternion setFromNormalVectors(final Vec3f v1, final Vec3f v2, final Vec3f tmpPivotVec) {
+        final float factor = v1.length() * v2.length();
         if ( FloatUtil.isZero(factor, FloatUtil.EPSILON ) ) {
             return setIdentity();
         } else {
-            final float dot = VectorUtil.dotVec3(v1, v2) / factor; // normalize
+            final float dot = v1.dot(v2) / factor; // normalize
             final float theta = FloatUtil.acos(Math.max(-1.0f, Math.min(dot, 1.0f))); // clipping [-1..1]
 
-            VectorUtil.crossVec3(tmpPivotVec, v1, v2);
+            tmpPivotVec.cross(v1, v2);
 
-            if ( dot < 0.0f && FloatUtil.isZero( VectorUtil.normVec3(tmpPivotVec), FloatUtil.EPSILON ) ) {
+            if ( dot < 0.0f && FloatUtil.isZero( tmpPivotVec.length(), FloatUtil.EPSILON ) ) {
                 // Vectors parallel and opposite direction, therefore a rotation of 180 degrees about any vector
                 // perpendicular to this vector will rotate vector a onto vector b.
                 //
                 // The following guarantees the dot-product will be 0.0.
                 int dominantIndex;
-                if (Math.abs(v1[0]) > Math.abs(v1[1])) {
-                    if (Math.abs(v1[0]) > Math.abs(v1[2])) {
+                if (Math.abs(v1.x()) > Math.abs(v1.y())) {
+                    if (Math.abs(v1.x()) > Math.abs(v1.z())) {
                         dominantIndex = 0;
                     } else {
                         dominantIndex = 2;
                     }
                 } else {
-                    if (Math.abs(v1[1]) > Math.abs(v1[2])) {
+                    if (Math.abs(v1.y()) > Math.abs(v1.z())) {
                         dominantIndex = 1;
                     } else {
                         dominantIndex = 2;
                     }
                 }
-                tmpPivotVec[dominantIndex]           = -v1[(dominantIndex + 1) % 3];
-                tmpPivotVec[(dominantIndex + 1) % 3] =  v1[dominantIndex];
-                tmpPivotVec[(dominantIndex + 2) % 3] =  0f;
+                tmpPivotVec.set( dominantIndex,           -v1.get( (dominantIndex + 1) % 3 ) );
+                tmpPivotVec.set( (dominantIndex + 1) % 3,  v1.get( dominantIndex ) );
+                tmpPivotVec.set( (dominantIndex + 2) % 3,  0f );
             }
             return setFromAngleNormalAxis(theta, tmpPivotVec);
         }
@@ -748,14 +740,14 @@ public class Quaternion {
      * </p>
      * @param angle rotation angle (rads)
      * @param vector axis vector not normalized
-     * @param tmpV3f float[3] temp storage to normalize vector
+     * @param tmpV3f temp storage to normalize vector
      * @return this quaternion for chaining.
      *
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q56">Matrix-FAQ Q56</a>
-     * @see #toAngleAxis(float[])
+     * @see #toAngleAxis(Vec3f)
      */
-    public final Quaternion setFromAngleAxis(final float angle, final float[] vector, final float[] tmpV3f) {
-        VectorUtil.normalizeVec3(tmpV3f, vector);
+    public final Quaternion setFromAngleAxis(final float angle, final Vec3f vector, final Vec3f tmpV3f) {
+        tmpV3f.set(vector).normalize();
         return setFromAngleNormalAxis(angle, tmpV3f);
     }
 
@@ -772,17 +764,17 @@ public class Quaternion {
      * @return this quaternion for chaining.
      *
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q56">Matrix-FAQ Q56</a>
-     * @see #toAngleAxis(float[])
+     * @see #toAngleAxis(Vec3f)
      */
-    public final Quaternion setFromAngleNormalAxis(final float angle, final float[] vector) {
-        if ( VectorUtil.isVec3Zero(vector, 0, FloatUtil.EPSILON) ) {
+    public final Quaternion setFromAngleNormalAxis(final float angle, final Vec3f vector) {
+        if( vector.isZero() ) {
             setIdentity();
         } else {
             final float halfangle = angle * 0.5f;
             final float sin = FloatUtil.sin(halfangle);
-            x = vector[0] * sin;
-            y = vector[1] * sin;
-            z = vector[2] * sin;
+            x = vector.x() * sin;
+            y = vector.y() * sin;
+            z = vector.z() * sin;
             w = FloatUtil.cos(halfangle);
         }
         return this;
@@ -791,24 +783,22 @@ public class Quaternion {
     /**
      * Transform the rotational quaternion to axis based rotation angles
      *
-     * @param axis float[3] storage for computed axis
+     * @param axis storage for computed axis
      * @return the rotation angle in radians
-     * @see #setFromAngleAxis(float, float[], float[])
+     * @see #setFromAngleAxis(float, Vec3f, Vec3f)
      */
-    public final float toAngleAxis(final float[] axis) {
+    public final float toAngleAxis(final Vec3f axis) {
         final float sqrLength = x*x + y*y + z*z;
         float angle;
         if ( FloatUtil.isZero(sqrLength, FloatUtil.EPSILON) ) { // length is ~0
             angle = 0.0f;
-            axis[0] = 1.0f;
-            axis[1] = 0.0f;
-            axis[2] = 0.0f;
+            axis.set( 1.0f, 0.0f, 0.0f );
         } else {
             angle = FloatUtil.acos(w) * 2.0f;
             final float invLength = 1.0f / FloatUtil.sqrt(sqrLength);
-            axis[0] = x * invLength;
-            axis[1] = y * invLength;
-            axis[2] = z * invLength;
+            axis.set( x * invLength,
+                      y * invLength,
+                      z * invLength );
         }
         return angle;
     }
@@ -816,7 +806,7 @@ public class Quaternion {
     /**
      * Initializes this quaternion from the given Euler rotation array <code>angradXYZ</code> in radians.
      * <p>
-     * The <code>angradXYZ</code> array is laid out in natural order:
+     * The <code>angradXYZ</code> vector is laid out in natural order:
      * <ul>
      *  <li>x - bank</li>
      *  <li>y - heading</li>
@@ -824,12 +814,12 @@ public class Quaternion {
      * </ul>
      * </p>
      * For details see {@link #setFromEuler(float, float, float)}.
-     * @param angradXYZ euler angel array in radians
+     * @param angradXYZ euler angle vector in radians holding x-bank, y-heading and z-attitude
      * @return this quaternion for chaining.
      * @see #setFromEuler(float, float, float)
      */
-    public final Quaternion setFromEuler(final float[] angradXYZ) {
-        return setFromEuler(angradXYZ[0], angradXYZ[1], angradXYZ[2]);
+    public final Quaternion setFromEuler(final Vec3f angradXYZ) {
+        return setFromEuler(angradXYZ.x(), angradXYZ.y(), angradXYZ.z());
     }
 
     /**
@@ -857,7 +847,7 @@ public class Quaternion {
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q60">Matrix-FAQ Q60</a>
      * @see <a href="http://vered.rose.utoronto.ca/people/david_dir/GEMS/GEMS.html">Gems</a>
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToQuaternion/index.htm">euclideanspace.com-eulerToQuaternion</a>
-     * @see #toEuler(float[])
+     * @see #toEuler(Vec3f)
      */
     public final Quaternion setFromEuler(final float bankX, final float headingY, final float attitudeZ) {
         if ( VectorUtil.isZero(bankX, headingY, attitudeZ, FloatUtil.EPSILON) ) {
@@ -889,58 +879,45 @@ public class Quaternion {
 
     /**
      * Transform this quaternion to Euler rotation angles in radians (pitchX, yawY and rollZ).
+     * <p>
+     * The <code>result</code> array is laid out in natural order:
+     * <ul>
+     *  <li>x - bank</li>
+     *  <li>y - heading</li>
+     *  <li>z - attitude</li>
+     * </ul>
+     * </p>
      *
-     * @param result the float[] array storing the computed angles.
-     * @return the double[] array, filled with heading, attitude and bank in that order..
+     * @param result euler angle result vector for radians x-bank, y-heading and z-attitude
+     * @return the Vec3f `result` filled with x-bank, y-heading and z-attitude
      * @see <a href="http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/index.htm">euclideanspace.com-quaternionToEuler</a>
      * @see #setFromEuler(float, float, float)
      */
-    public float[] toEuler(final float[] result) {
+    public Vec3f toEuler(final Vec3f result) {
         final float sqw = w*w;
         final float sqx = x*x;
         final float sqy = y*y;
         final float sqz = z*z;
-        final float unit = sqx + sqy + sqz + sqw; // if normalized is one, otherwise
-        // is correction factor
+        final float unit = sqx + sqy + sqz + sqw; // if normalized is one, otherwise, is correction factor
         final float test = x*y + z*w;
 
         if (test > 0.499f * unit) { // singularity at north pole
-            result[0] =  0f;
-            result[1] =  2f * FloatUtil.atan2(x, w);
-            result[2] =  FloatUtil.HALF_PI;
+            result.set( 0f,                          // x-bank
+                        2f * FloatUtil.atan2(x, w),  // y-heading
+                        FloatUtil.HALF_PI );         // z-attitude
         } else if (test < -0.499f * unit) { // singularity at south pole
-            result[0] =  0f;
-            result[1] = -2 * FloatUtil.atan2(x, w);
-            result[2] = -FloatUtil.HALF_PI;
+            result.set( 0f,                          // x-bank
+                       -2 * FloatUtil.atan2(x, w),   // y-heading
+                       -FloatUtil.HALF_PI );         // z-attitude
         } else {
-            result[0] = FloatUtil.atan2(2f * x * w - 2 * y * z, -sqx + sqy - sqz + sqw);
-            result[1] = FloatUtil.atan2(2f * y * w - 2 * x * z,  sqx - sqy - sqz + sqw);
-            result[2] = FloatUtil.asin( 2f * test / unit);
+            result.set( FloatUtil.atan2(2f * x * w - 2 * y * z, -sqx + sqy - sqz + sqw), // x-bank
+                        FloatUtil.atan2(2f * y * w - 2 * x * z,  sqx - sqy - sqz + sqw), // y-heading
+                        FloatUtil.asin( 2f * test / unit) );                             // z-attitude
         }
         return result;
     }
 
     /**
-     * Initializes this quaternion from a 4x4 column rotation matrix
-     * <p>
-     * See <a href="ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/quatut.ps.Z">Graphics Gems Code</a>,<br/>
-     * <a href="http://mathworld.wolfram.com/MatrixTrace.html">MatrixTrace</a>.
-     * </p>
-     * <p>
-     * Buggy <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q55">Matrix-FAQ Q55</a>
-     * </p>
-     *
-     * @param m 4x4 column matrix
-     * @return this quaternion for chaining.
-     * @see #toMatrix(float[], int)
-     */
-    public final Quaternion setFromMatrix(final float[] m, final int m_off) {
-        return setFromMatrix(m[0+0*4+m_off], m[0+1*4+m_off], m[0+2*4+m_off],
-                             m[1+0*4+m_off], m[1+1*4+m_off], m[1+2*4+m_off],
-                             m[2+0*4+m_off], m[2+1*4+m_off], m[2+2*4+m_off]);
-    }
-
-    /**
      * Compute the quaternion from a 3x3 column rotation matrix
      * <p>
      * See <a href="ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/quatut.ps.Z">Graphics Gems Code</a>,<br/>
@@ -951,7 +928,7 @@ public class Quaternion {
      * </p>
      *
      * @return this quaternion for chaining.
-     * @see #toMatrix(float[], int)
+     * @see #setFromMatrix(Matrix4f)
      */
     public Quaternion setFromMatrix(final float m00, final float m01, final float m02,
                                     final float m10, final float m11, final float m12,
@@ -996,6 +973,24 @@ public class Quaternion {
     }
 
     /**
+     * Compute the quaternion from a 3x3 column rotation matrix
+     * <p>
+     * See <a href="ftp://ftp.cis.upenn.edu/pub/graphics/shoemake/quatut.ps.Z">Graphics Gems Code</a>,<br/>
+     * <a href="http://mathworld.wolfram.com/MatrixTrace.html">MatrixTrace</a>.
+     * </p>
+     * <p>
+     * Buggy <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q55">Matrix-FAQ Q55</a>
+     * </p>
+     *
+     * @return this quaternion for chaining.
+     * @see Matrix4f#getRotation(Quaternion)
+     * @see #setFromMatrix(float, float, float, float, float, float, float, float, float)
+     */
+    public Quaternion setFromMatrix(final Matrix4f m) {
+        return m.getRotation(this);
+    }
+
+    /**
      * Transform this quaternion to a normalized 4x4 column matrix representing the rotation.
      * <p>
      * Implementation Details:
@@ -1008,7 +1003,8 @@ public class Quaternion {
      * @param mat_offset
      * @return the given matrix store
      * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q54">Matrix-FAQ Q54</a>
-     * @see #setFromMatrix(float[], int)
+     * @see #setFromMatrix(Matrix4f)
+     * @see #setFromMatrix(float, float, float, float, float, float, float, float, float)
      */
     public final float[] toMatrix(final float[] matrix, final int mat_offset) {
         // pre-multiply scaled-reciprocal-magnitude to reduce multiplications
@@ -1061,58 +1057,22 @@ public class Quaternion {
     }
 
     /**
-     * @param index the 3x3 rotation matrix column to retrieve from this quaternion (normalized). Must be between 0 and 2.
-     * @param result the vector object to store the result in.
-     * @return the result column-vector for chaining.
+     * Transform this quaternion to a normalized 4x4 column matrix representing the rotation.
+     * <p>
+     * Implementation Details:
+     * <ul>
+     *   <li> makes identity matrix if {@link #magnitudeSquared()} is {@link FloatUtil#isZero(float, float) is zero} using {@link FloatUtil#EPSILON epsilon}</li>
+     * </ul>
+     * </p>
+     *
+     * @param matrix store for the resulting normalized column matrix 4x4
+     * @return the given matrix store
+     * @see <a href="http://web.archive.org/web/20041029003853/http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q54">Matrix-FAQ Q54</a>
+     * @see #setFromMatrix(float, float, float, float, float, float, float, float, float)
+     * @see Matrix4f#setToRotation(Quaternion)
      */
-    public float[] copyMatrixColumn(final int index, final float[] result, final int resultOffset) {
-        // pre-multipliy scaled-reciprocal-magnitude to reduce multiplications
-        final float norm = magnitudeSquared();
-        final float srecip;
-        if ( FloatUtil.isZero(norm, FloatUtil.EPSILON) ) {
-            srecip= 0f;
-        } else if ( FloatUtil.isEqual(1f, norm, FloatUtil.EPSILON) ) {
-            srecip= 2f;
-        } else {
-            srecip= 2.0f / norm;
-        }
-
-        // compute xs/ys/zs first to save 6 multiplications, since xs/ys/zs
-        // will be used 2-4 times each.
-        final float xs = x * srecip;
-        final float ys = y * srecip;
-        final float zs = z * srecip;
-        final float xx = x * xs;
-        final float xy = x * ys;
-        final float xz = x * zs;
-        final float xw = w * xs;
-        final float yy = y * ys;
-        final float yz = y * zs;
-        final float yw = w * ys;
-        final float zz = z * zs;
-        final float zw = w * zs;
-
-        // using s=2/norm (instead of 1/norm) saves 3 multiplications by 2 here
-        switch (index) {
-        case 0:
-            result[0+resultOffset] = 1.0f - (yy + zz);
-            result[1+resultOffset] = xy + zw;
-            result[2+resultOffset] = xz - yw;
-            break;
-        case 1:
-            result[0+resultOffset] = xy - zw;
-            result[1+resultOffset] = 1.0f - (xx + zz);
-            result[2+resultOffset] = yz + xw;
-            break;
-        case 2:
-            result[0+resultOffset] = xz + yw;
-            result[1+resultOffset] = yz - xw;
-            result[2+resultOffset] = 1.0f - (xx + yy);
-            break;
-        default:
-            throw new IllegalArgumentException("Invalid column index. " + index);
-        }
-        return result;
+    public final Matrix4f toMatrix(final Matrix4f matrix) {
+        return matrix.setToRotation(this);
     }
 
     /**
@@ -1126,10 +1086,10 @@ public class Quaternion {
      * @param zAxis vector representing the <i>orthogonal</i> z-axis of the coordinate system.
      * @return this quaternion for chaining.
      */
-    public final Quaternion setFromAxes(final float[] xAxis, final float[] yAxis, final float[] zAxis) {
-        return setFromMatrix(xAxis[0], yAxis[0], zAxis[0],
-                             xAxis[1], yAxis[1], zAxis[1],
-                             xAxis[2], yAxis[2], zAxis[2]);
+    public final Quaternion setFromAxes(final Vec3f xAxis, final Vec3f yAxis, final Vec3f zAxis) {
+        return setFromMatrix(xAxis.x(), yAxis.x(), zAxis.x(),
+                             xAxis.y(), yAxis.y(), zAxis.y(),
+                             xAxis.z(), yAxis.z(), zAxis.z());
     }
 
     /**
@@ -1140,11 +1100,11 @@ public class Quaternion {
      * @param zAxis vector representing the <i>orthogonal</i> z-axis of the coordinate system.
      * @param tmpMat4 temporary float[4] matrix, used to transform this quaternion to a matrix.
      */
-    public void toAxes(final float[] xAxis, final float[] yAxis, final float[] zAxis, final float[] tmpMat4) {
-        toMatrix(tmpMat4, 0);
-        FloatUtil.copyMatrixColumn(tmpMat4, 0, 2, zAxis, 0);
-        FloatUtil.copyMatrixColumn(tmpMat4, 0, 1, yAxis, 0);
-        FloatUtil.copyMatrixColumn(tmpMat4, 0, 0, xAxis, 0);
+    public void toAxes(final Vec3f xAxis, final Vec3f yAxis, final Vec3f zAxis, final Matrix4f tmpMat4) {
+        tmpMat4.setToRotation(this);
+        tmpMat4.getColumn(2, zAxis);
+        tmpMat4.getColumn(1, yAxis);
+        tmpMat4.getColumn(0, xAxis);
     }
 
     /**
@@ -1154,6 +1114,7 @@ public class Quaternion {
      * @param m 3x3 column matrix
      * @return true if representing a rotational matrix, false otherwise
      */
+    @Deprecated
     public final boolean isRotationMatrix3f(final 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)
@@ -1171,6 +1132,7 @@ public class Quaternion {
         return (FloatUtil.abs(determinant3f(m) - 1) < epsilon);
     }
 
+    @Deprecated
     private final float determinant3f(final 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];
@@ -1193,10 +1155,10 @@ public class Quaternion {
             return false;
         }
         final Quaternion comp = (Quaternion) o;
-        return Math.abs(x - comp.getX()) <= ALLOWED_DEVIANCE &&
-               Math.abs(y - comp.getY()) <= ALLOWED_DEVIANCE &&
-               Math.abs(z - comp.getZ()) <= ALLOWED_DEVIANCE &&
-               Math.abs(w - comp.getW()) <= ALLOWED_DEVIANCE;
+        return Math.abs(x - comp.x()) <= ALLOWED_DEVIANCE &&
+               Math.abs(y - comp.y()) <= ALLOWED_DEVIANCE &&
+               Math.abs(z - comp.z()) <= ALLOWED_DEVIANCE &&
+               Math.abs(w - comp.w()) <= ALLOWED_DEVIANCE;
     }
     @Override
     public final int hashCode() {
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Ray.java b/src/jogl/classes/com/jogamp/opengl/math/Ray.java
index a528f0763..25a7d9a70 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/Ray.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/Ray.java
@@ -47,14 +47,14 @@ import com.jogamp.opengl.math.geom.AABBox;
  * </p>
  */
 public class Ray {
-    /** Origin of Ray, float[3]. */
-    public final float[] orig = new float[3];
+    /** Origin of Ray. */
+    public final Vec3f orig = new Vec3f();
 
-    /** Normalized direction vector of ray, float[3]. */
-    public final float[] dir = new float[3];
+    /** Normalized direction vector of ray. */
+    public final Vec3f dir = new Vec3f();
 
     @Override
     public String toString() {
-        return "Ray[orig["+orig[0]+", "+orig[1]+", "+orig[2]+"], dir["+dir[0]+", "+dir[1]+", "+dir[2]+"]]";
+        return "Ray[orig["+orig+"], dir["+dir+"]]";
     }
 }
\ No newline at end of file
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Vec2f.java b/src/jogl/classes/com/jogamp/opengl/math/Vec2f.java
index 0b67102fb..0c7854216 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/Vec2f.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/Vec2f.java
@@ -48,6 +48,11 @@ public final class Vec2f {
         set(o);
     }
 
+    /** Creating new Vec2f using Vec3f, dropping z. */
+    public Vec2f(final Vec3f o) {
+        set(o);
+    }
+
     public Vec2f copy() {
         return new Vec2f(this);
     }
@@ -66,6 +71,12 @@ public final class Vec2f {
         this.y = o.y;
     }
 
+    /** this = o while dropping z, returns this. */
+    public void set(final Vec3f o) {
+        this.x = o.x();
+        this.y = o.y();
+    }
+
     /** this = { x, y }, returns this. */
     public void set(final float x, final float y) {
         this.x = x;
@@ -88,6 +99,13 @@ public final class Vec2f {
         }
     }
 
+    /** xy = this, returns xy. */
+    public float[] get(final float[/*2*/] xy) {
+        xy[0] = this.x;
+        xy[1] = this.y;
+        return xy;
+    }
+
     /** Gets the ith component, 0 <= i < 2 */
     public float get(final int i) {
         switch (i) {
@@ -108,10 +126,17 @@ public final class Vec2f {
         return new Vec2f(this).scale(val);
     }
 
-    /** this = this * val, returns this. */
-    public Vec2f scale(final float val) {
-        x *= val;
-        y *= val;
+    /** this = this * s, returns this. */
+    public Vec2f scale(final float s) {
+        x *= s;
+        y *= s;
+        return this;
+    }
+
+    /** this = this * { sx, sy }, returns this. */
+    public Vec2f scale(final float sx, final float sy) {
+        x *= sx;
+        y *= sy;
         return this;
     }
 
@@ -158,6 +183,7 @@ public final class Vec2f {
         return this;
     }
 
+    /** Return true if all components are zero, i.e. it's absolute value < {@link #EPSILON}. */
     public boolean isZero() {
         return FloatUtil.isZero(x) && FloatUtil.isZero(y);
     }
@@ -282,6 +308,53 @@ public final class Vec2f {
         return true;
     }
 
+    /**
+     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @param epsilon consider using {@link FloatUtil#EPSILON}
+     * @return true if all components differ less than {@code epsilon}, otherwise false.
+     */
+    public boolean isEqual(final Vec2f o, final float epsilon) {
+        if( this == o ) {
+            return true;
+        } else {
+            return FloatUtil.isEqual(x, o.x, epsilon) &&
+                   FloatUtil.isEqual(y, o.y, epsilon);
+        }
+    }
+
+    /**
+     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @return true if all components differ less than {@link FloatUtil#EPSILON}, otherwise false.
+     */
+    public boolean isEqual(final Vec2f o) {
+        return isEqual(o, FloatUtil.EPSILON);
+    }
+
+    @Override
+    public boolean equals(final Object o) {
+        if( o instanceof Vec2f ) {
+            return isEqual((Vec2f)o, FloatUtil.EPSILON);
+        } else {
+            return false;
+        }
+    }
+
     @Override
     public String toString() {
         return x + " / " + y;
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Vec3f.java b/src/jogl/classes/com/jogamp/opengl/math/Vec3f.java
index 9cead54c2..d5c725ad7 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/Vec3f.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/Vec3f.java
@@ -35,6 +35,12 @@ package com.jogamp.opengl.math;
  * and its data layout from JOAL's Vec3f.
  */
 public final class Vec3f {
+    public static final Vec3f ONE = new Vec3f(VectorUtil.VEC3_ONE);
+    public static final Vec3f UNIT_Y = new Vec3f(VectorUtil.VEC3_UNIT_Y);
+    public static final Vec3f UNIT_Y_NEG = new Vec3f(VectorUtil.VEC3_UNIT_Y_NEG);
+    public static final Vec3f UNIT_Z = new Vec3f(VectorUtil.VEC3_UNIT_Z);
+    public static final Vec3f UNIT_Z_NEG = new Vec3f(VectorUtil.VEC3_UNIT_Z_NEG);
+
     private float x;
     private float y;
     private float z;
@@ -45,6 +51,16 @@ public final class Vec3f {
         set(o);
     }
 
+    /** Creating new Vec3f using Vec4f, dropping w. */
+    public Vec3f(final Vec4f o) {
+        set(o);
+    }
+
+    /** Creating new Vec3f using { Vec2f, z}. */
+    public Vec3f(final Vec2f o, final float z) {
+        set(o, z);
+    }
+
     public Vec3f copy() {
         return new Vec3f(this);
     }
@@ -65,6 +81,22 @@ public final class Vec3f {
         return this;
     }
 
+    /** this = { o, z }, returns this. */
+    public Vec3f set(final Vec2f o, final float z) {
+        this.x = o.x();
+        this.y = o.y();
+        this.z = z;
+        return this;
+    }
+
+    /** this = o while dropping w, returns this. */
+    public Vec3f set(final Vec4f o) {
+        this.x = o.x();
+        this.y = o.y();
+        this.z = o.z();
+        return this;
+    }
+
     /** this = { x, y, z }, returns this. */
     public Vec3f set(final float x, final float y, final float z) {
         this.x = x;
@@ -91,6 +123,14 @@ public final class Vec3f {
         }
     }
 
+    /** xyz = this, returns xyz. */
+    public float[] get(final float[/*3*/] xyz) {
+        xyz[0] = this.x;
+        xyz[1] = this.y;
+        xyz[2] = this.z;
+        return xyz;
+    }
+
     /** Gets the ith component, 0 <= i < 3 */
     public float get(final int i) {
         switch (i) {
@@ -114,11 +154,19 @@ public final class Vec3f {
         return new Vec3f(this).scale(val);
     }
 
-    /** this = this * val, returns this. */
-    public Vec3f scale(final float val) {
-        x *= val;
-        y *= val;
-        z *= val;
+    /** this = this * s, returns this. */
+    public Vec3f scale(final float s) {
+        x *= s;
+        y *= s;
+        z *= s;
+        return this;
+    }
+
+    /** this = this * { sx, sy, sz }, returns this. */
+    public Vec3f scale(final float sx, final float sy, final float sz) {
+        x *= sx;
+        y *= sy;
+        z *= sz;
         return this;
     }
 
@@ -169,7 +217,8 @@ public final class Vec3f {
         return this;
     }
 
-    public boolean is_zero() {
+    /** Return true if all components are zero, i.e. it's absolute value < {@link #EPSILON}. */
+    public boolean isZero() {
         return FloatUtil.isZero(x) && FloatUtil.isZero(y) && FloatUtil.isZero(z);
     }
 
@@ -266,8 +315,57 @@ public final class Vec3f {
     public boolean intersects(final Vec3f o) {
         if( Math.abs(x-o.x) >= FloatUtil.EPSILON || Math.abs(y-o.y) >= FloatUtil.EPSILON || Math.abs(z-o.z) >= FloatUtil.EPSILON ) {
             return false;
+        } else {
+            return true;
+        }
+    }
+
+    /**
+     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @param epsilon consider using {@link FloatUtil#EPSILON}
+     * @return true if all components differ less than {@code epsilon}, otherwise false.
+     */
+    public boolean isEqual(final Vec3f o, final float epsilon) {
+        if( this == o ) {
+            return true;
+        } else {
+            return FloatUtil.isEqual(x, o.x, epsilon) &&
+                   FloatUtil.isEqual(y, o.y, epsilon) &&
+                   FloatUtil.isEqual(z, o.z, epsilon);
+        }
+    }
+
+    /**
+     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @return true if all components differ less than {@link FloatUtil#EPSILON}, otherwise false.
+     */
+    public boolean isEqual(final Vec3f o) {
+        return isEqual(o, FloatUtil.EPSILON);
+    }
+
+    @Override
+    public boolean equals(final Object o) {
+        if( o instanceof Vec3f ) {
+            return isEqual((Vec3f)o, FloatUtil.EPSILON);
+        } else {
+            return false;
         }
-        return true;
     }
 
     @Override
diff --git a/src/jogl/classes/com/jogamp/opengl/math/Vec4f.java b/src/jogl/classes/com/jogamp/opengl/math/Vec4f.java
new file mode 100644
index 000000000..1a20015a9
--- /dev/null
+++ b/src/jogl/classes/com/jogamp/opengl/math/Vec4f.java
@@ -0,0 +1,348 @@
+/**
+ * Copyright 2022-2023 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;
+
+/**
+ * 4D Vector based upon four float components.
+ *
+ * Implementation borrowed from [gfxbox2](https://jausoft.com/cgit/cs_class/gfxbox2.git/tree/include/pixel/pixel3f.hpp#n29)
+ * and its data layout from JOAL's Vec3f.
+ */
+public final class Vec4f {
+    private float x;
+    private float y;
+    private float z;
+    private float w;
+
+    public Vec4f() {}
+
+    public Vec4f(final Vec4f o) {
+        set(o);
+    }
+
+    /** Creating new Vec4f using { o, w }. */
+    public Vec4f(final Vec3f o, final float w) {
+        set(o, w);
+    }
+
+    public Vec4f copy() {
+        return new Vec4f(this);
+    }
+
+    public Vec4f(final float[/*4*/] xyzw) {
+        set(xyzw);
+    }
+
+    public Vec4f(final float x, final float y, final float z, final float w) {
+        set(x, y, z, w);
+    }
+
+    /** this = o, returns this. */
+    public Vec4f set(final Vec4f o) {
+        this.x = o.x;
+        this.y = o.y;
+        this.z = o.z;
+        this.w = o.w;
+        return this;
+    }
+
+    /** this = { o, w }, returns this. */
+    public Vec4f set(final Vec3f o, final float w) {
+        this.x = o.x();
+        this.y = o.y();
+        this.z = o.z();
+        this.w = w;
+        return this;
+    }
+
+    /** this = { x, y, z, w }, returns this. */
+    public Vec4f set(final float x, final float y, final float z, final float w) {
+        this.x = x;
+        this.y = y;
+        this.z = z;
+        this.w = w;
+        return this;
+    }
+
+    /** this = xyzw, returns this. */
+    public Vec4f set(final float[/*4*/] xyzw) {
+        this.x = xyzw[0];
+        this.y = xyzw[1];
+        this.z = xyzw[2];
+        this.w = xyzw[3];
+        return this;
+    }
+
+    /** Sets the ith component, 0 <= i < 4 */
+    public void set(final int i, final float val) {
+        switch (i) {
+            case 0: x = val; break;
+            case 1: y = val; break;
+            case 2: z = val; break;
+            case 3: w = val; break;
+            default: throw new IndexOutOfBoundsException();
+        }
+    }
+
+    /** xyzw = this, returns xyzw. */
+    public float[] get(final float[/*4*/] xyzw) {
+        xyzw[0] = this.x;
+        xyzw[1] = this.y;
+        xyzw[2] = this.z;
+        xyzw[3] = this.w;
+        return xyzw;
+    }
+
+    /** Gets the ith component, 0 <= i < 4 */
+    public float get(final int i) {
+        switch (i) {
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            case 3: return w;
+            default: throw new IndexOutOfBoundsException();
+        }
+    }
+
+    public float x() { return x; }
+    public float y() { return y; }
+    public float z() { return z; }
+    public float w() { return w; }
+
+    public void setX(final float x) { this.x = x; }
+    public void setY(final float y) { this.y = y; }
+    public void setZ(final float z) { this.z = z; }
+    public void setW(final float w) { this.w = w; }
+
+    /** Returns this * val; creates new vector */
+    public Vec4f mul(final float val) {
+        return new Vec4f(this).scale(val);
+    }
+
+    /** this = this * s, returns this. */
+    public Vec4f scale(final float s) {
+        x *= s;
+        y *= s;
+        z *= s;
+        w *= s;
+        return this;
+    }
+
+    /** this = this * { sx, sy, sz, sw }, returns this. */
+    public Vec4f scale(final float sx, final float sy, final float sz, final float sw) {
+        x *= sx;
+        y *= sy;
+        z *= sz;
+        w *= sw;
+        return this;
+    }
+
+    /** Returns this + arg; creates new vector */
+    public Vec4f plus(final Vec4f arg) {
+        return new Vec4f(this).add(arg);
+    }
+
+    /** this = this + { dx, dy, dz, dw }, returns this. */
+    public Vec4f add(final float dx, final float dy, final float dz, final float dw) {
+        x += dx;
+        y += dy;
+        z += dz;
+        w += dw;
+        return this;
+    }
+
+    /** this = this + b, returns this. */
+    public Vec4f add(final Vec4f b) {
+        x += b.x;
+        y += b.y;
+        z += b.z;
+        w += b.w;
+        return this;
+    }
+
+    /** Returns this - arg; creates new vector */
+    public Vec4f minus(final Vec4f arg) {
+        return new Vec4f(this).sub(arg);
+    }
+
+    /** this = this - b, returns this. */
+    public Vec4f sub(final Vec4f b) {
+        x -= b.x;
+        y -= b.y;
+        z -= b.z;
+        w -= b.w;
+        return this;
+    }
+
+    /** Return true if all components are zero, i.e. it's absolute value < {@link #EPSILON}. */
+    public boolean isZero() {
+        return FloatUtil.isZero(x) && FloatUtil.isZero(y) && FloatUtil.isZero(z) && FloatUtil.isZero(w);
+    }
+
+    /**
+     * Return the length of this vector, a.k.a the <i>norm</i> or <i>magnitude</i>
+     */
+    public float length() {
+        return (float) Math.sqrt(lengthSq());
+    }
+
+    /**
+     * Return the squared length of this vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
+     */
+    public float lengthSq() {
+        return x*x + y*y + z*z + w*w;
+    }
+
+    /**
+     * Normalize this vector in place
+     */
+    public Vec4f normalize() {
+        final float lengthSq = lengthSq();
+        if ( FloatUtil.isZero( lengthSq ) ) {
+            x = 0.0f;
+            y = 0.0f;
+            z = 0.0f;
+            w = 0.0f;
+        } else {
+            final float invSqr = 1.0f / (float)Math.sqrt(lengthSq);
+            x *= invSqr;
+            y *= invSqr;
+            z *= invSqr;
+            w *= invSqr;
+        }
+        return this;
+    }
+
+    /**
+     * Return the squared distance between this vector and the given one.
+     * <p>
+     * When comparing the relative distance between two points it is usually sufficient to compare the squared
+     * distances, thus avoiding an expensive square root operation.
+     * </p>
+     */
+    public float distSq(final Vec4f o) {
+        final float dx = x - o.x;
+        final float dy = y - o.y;
+        final float dz = z - o.z;
+        final float dw = w - o.w;
+        return dx*dx + dy*dy + dz*dz + dw*dw;
+    }
+
+    /**
+     * Return the distance between this vector and the given one.
+     */
+    public float dist(final Vec4f o) {
+        return (float)Math.sqrt(distSq(o));
+    }
+
+
+    /**
+     * Return the dot product of this vector and the given one
+     * @return the dot product as float
+     */
+    public float dot(final Vec4f o) {
+        return x*o.x + y*o.y + z*o.z + w*o.w;
+    }
+
+    /**
+     * Return the cosines of the angle between two vectors
+     */
+    public float cosAngle(final Vec4f o) {
+        return dot(o) / ( length() * o.length() ) ;
+    }
+
+    /**
+     * Return the angle between two vectors in radians
+     */
+    public float angle(final Vec4f o) {
+        return (float) Math.acos( cosAngle(o) );
+    }
+
+    public boolean intersects(final Vec4f o) {
+        if( Math.abs(x-o.x) >= FloatUtil.EPSILON || Math.abs(y-o.y) >= FloatUtil.EPSILON || Math.abs(z-o.z) >= FloatUtil.EPSILON ||
+            Math.abs(w-o.w) >= FloatUtil.EPSILON) {
+            return false;
+        }
+        return true;
+    }
+
+    /**
+     * Equals check using a given {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @param epsilon consider using {@link FloatUtil#EPSILON}
+     * @return true if all components differ less than {@code epsilon}, otherwise false.
+     */
+    public boolean isEqual(final Vec4f o, final float epsilon) {
+        if( this == o ) {
+            return true;
+        } else {
+            return FloatUtil.isEqual(x, o.x, epsilon) &&
+                   FloatUtil.isEqual(y, o.y, epsilon) &&
+                   FloatUtil.isEqual(z, o.z, epsilon) &&
+                   FloatUtil.isEqual(w, o.w, epsilon);
+        }
+    }
+
+    /**
+     * Equals check using {@link FloatUtil#EPSILON} value and {@link FloatUtil#isEqual(float, float, float)}.
+     * <p>
+     * Implementation considers following corner cases:
+     * <ul>
+     *    <li>NaN == NaN</li>
+     *    <li>+Inf == +Inf</li>
+     *    <li>-Inf == -Inf</li>
+     * </ul>
+     * @param o comparison value
+     * @return true if all components differ less than {@link FloatUtil#EPSILON}, otherwise false.
+     */
+    public boolean isEqual(final Vec4f o) {
+        return isEqual(o, FloatUtil.EPSILON);
+    }
+
+    @Override
+    public boolean equals(final Object o) {
+        if( o instanceof Vec4f ) {
+            return isEqual((Vec4f)o, FloatUtil.EPSILON);
+        } else {
+            return false;
+        }
+    }
+
+    @Override
+    public String toString() {
+        return x + " / " + y + " / " + z + " / " + w;
+    }
+}
diff --git a/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
index a07153155..8edbd0cd7 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/VectorUtil.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2010 JogAmp Community. All rights reserved.
+ * Copyright 2010-2023 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:
@@ -34,7 +34,6 @@ import com.jogamp.graph.geom.plane.Winding;
 public final class VectorUtil {
 
     public static final float[] VEC3_ONE = { 1f, 1f, 1f };
-    public static final float[] VEC3_ZERO = { 0f, 0f, 0f };
     public static final float[] VEC3_UNIT_Y = { 0f, 1f, 0f };
     public static final float[] VEC3_UNIT_Y_NEG = { 0f, -1f, 0f };
     public static final float[] VEC3_UNIT_Z = { 0f, 0f, 1f };
@@ -1025,14 +1024,15 @@ public final class VectorUtil {
      * @param epsilon
      * @return resulting intersecting if exists, otherwise null
      */
-    public static float[] line2PlaneIntersection(final float[] result, final Ray ray, final float[/*4*/] plane, final float epsilon) {
-        final float tmp = dotVec3(ray.dir, plane) ;
+    public static Vec3f line2PlaneIntersection(final Vec3f result, final Ray ray, final Vec4f plane, final float epsilon) {
+        final Vec3f plane3 = new Vec3f(plane);
+        final float tmp = ray.dir.dot(plane3);
 
         if ( Math.abs(tmp) < epsilon ) {
             return null; // ray is parallel to plane
         }
-        scaleVec3 ( result, ray.dir, -( dotVec3(ray.orig, plane) + plane[3] ) / tmp ) ;
-        return addVec3(result, result, ray.orig);
+        result.set( ray.dir );
+        return result.scale( -( ray.orig.dot(plane3) + plane.w() ) / tmp ).add(ray.orig);
     }
 
     /** Compute intersection between two segments
diff --git a/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java b/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java
index f858b1c0d..9e8edfbf0 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/geom/AABBox.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2010 JogAmp Community. All rights reserved.
+ * Copyright 2010-2023 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:
@@ -29,9 +29,11 @@ package com.jogamp.opengl.math.geom;
 
 import com.jogamp.graph.geom.plane.AffineTransform;
 import com.jogamp.opengl.math.FloatUtil;
+import com.jogamp.opengl.math.Matrix4f;
 import com.jogamp.opengl.math.Quaternion;
 import com.jogamp.opengl.math.Ray;
-import com.jogamp.opengl.math.VectorUtil;
+import com.jogamp.opengl.math.Vec3f;
+import com.jogamp.opengl.util.PMVMatrix;
 
 
 /**
@@ -51,9 +53,9 @@ import com.jogamp.opengl.math.VectorUtil;
  */
 public class AABBox {
     private static final boolean DEBUG = FloatUtil.DEBUG;
-    private final float[] low = new float[3];
-    private final float[] high = new float[3];
-    private final float[] center = new float[3];
+    private final Vec3f low = new Vec3f();
+    private final Vec3f high = new Vec3f();
+    private final Vec3f center = new Vec3f();
 
     /**
      * Create an Axis Aligned bounding box (AABBox) with the
@@ -109,42 +111,34 @@ public class AABBox {
     public final AABBox 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;
+        center.set( 0f, 0f, 0f);
         return this;
     }
 
     /** Get the max xyz-coordinates
-     * @return a float array containing the max xyz coordinates
+     * @return max xyz coordinates
      */
-    public final float[] getHigh() {
+    public final Vec3f getHigh() {
         return high;
     }
 
     private final void setHigh(final float hx, final float hy, final float hz) {
-        this.high[0] = hx;
-        this.high[1] = hy;
-        this.high[2] = hz;
+        this.high.set(hx, hy, hz);
     }
 
     /** Get the min xyz-coordinates
-     * @return a float array containing the min xyz coordinates
+     * @return min xyz coordinates
      */
-    public final float[] getLow() {
+    public final Vec3f getLow() {
         return low;
     }
 
     private final void setLow(final float lx, final float ly, final float lz) {
-        this.low[0] = lx;
-        this.low[1] = ly;
-        this.low[2] = lz;
+        this.low.set(lx, ly, lz);
     }
 
     private final void computeCenter() {
-        center[0] = (high[0] + low[0])/2f;
-        center[1] = (high[1] + low[1])/2f;
-        center[2] = (high[2] + low[2])/2f;
+        center.set(high).add(low).scale(1f/2f);
     }
 
     /**
@@ -154,9 +148,9 @@ public class AABBox {
      * @return this AABBox for chaining
      */
     public final AABBox copy(final AABBox src) {
-        System.arraycopy(src.low, 0, low, 0, 3);
-        System.arraycopy(src.high, 0, high, 0, 3);
-        System.arraycopy(src.center, 0, center, 0, 3);
+        low.set(src.low);
+        high.set(src.high);
+        center.set(src.center);
         return this;
     }
 
@@ -186,12 +180,23 @@ public class AABBox {
      */
     public final AABBox setSize(final float lx, final float ly, final float lz,
                                 final float hx, final float hy, final 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;
+        this.low.set(lx, ly, lz);
+        this.high.set(hx, hy, hz);
+        computeCenter();
+        return this;
+    }
+
+    /**
+     * Set size of the AABBox specifying the coordinates
+     * of the low and high.
+     *
+     * @param low min xyz-coordinates
+     * @param high max xyz-coordinates
+     * @return this AABBox for chaining
+     */
+    public final AABBox setSize(final Vec3f low, final Vec3f high) {
+        this.low.set(low);
+        this.high.set(high);
         computeCenter();
         return this;
     }
@@ -202,25 +207,30 @@ public class AABBox {
      * @return this AABBox for chaining
      */
     public final AABBox resize(final AABBox newBox) {
-        final float[] newLow = newBox.getLow();
-        final float[] newHigh = newBox.getHigh();
+        final Vec3f newLow = newBox.getLow();
+        final Vec3f 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];
+        if (newLow.x() < low.x()) {
+            low.setX( newLow.x() );
+        }
+        if (newLow.y() < low.y()) {
+            low.setY( newLow.y() );
+        }
+        if (newLow.z() < low.z()) {
+            low.setZ( newLow.z() );
+        }
 
         /** 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];
-
+        if (newHigh.x() > high.x()) {
+            high.setX( newHigh.x() );
+        }
+        if (newHigh.y() > high.y()) {
+            high.setY( newHigh.y() );
+        }
+        if (newHigh.z() > high.z()) {
+            high.setZ( newHigh.z() );
+        }
         computeCenter();
         return this;
     }
@@ -229,34 +239,32 @@ public class AABBox {
      * Resize the AABBox to encapsulate another AABox, which will be <i>transformed</i> on the fly first.
      * @param newBox AABBox to be encapsulated in
      * @param t the {@link AffineTransform} applied on <i>newBox</i> on the fly
-     * @param tmpV3 temp float[3] storage
+     * @param tmpV3 temporary storage
      * @return this AABBox for chaining
      */
-    public final AABBox resize(final AABBox newBox, final AffineTransform t, final float[] tmpV3) {
+    public final AABBox resize(final AABBox newBox, final AffineTransform t, final Vec3f tmpV3) {
         /** test low */
         {
-            final float[] newBoxLow = newBox.getLow();
+            final Vec3f newBoxLow = newBox.getLow();
             t.transform(newBoxLow, tmpV3);
-            tmpV3[2] = newBoxLow[2];
-            if (tmpV3[0] < low[0])
-                low[0] = tmpV3[0];
-            if (tmpV3[1] < low[1])
-                low[1] = tmpV3[1];
-            if (tmpV3[2] < low[2])
-                low[2] = tmpV3[2];
+            if (tmpV3.x() < low.x())
+                low.setX( tmpV3.x() );
+            if (tmpV3.y() < low.y())
+                low.setY( tmpV3.y() );
+            if (tmpV3.z() < low.z())
+                low.setZ( tmpV3.z() );
         }
 
         /** test high */
         {
-            final float[] newBoxHigh = newBox.getHigh();
+            final Vec3f newBoxHigh = newBox.getHigh();
             t.transform(newBoxHigh, tmpV3);
-            tmpV3[2] = newBoxHigh[2];
-            if (tmpV3[0] > high[0])
-                high[0] = tmpV3[0];
-            if (tmpV3[1] > high[1])
-                high[1] = tmpV3[1];
-            if (tmpV3[2] > high[2])
-                high[2] = tmpV3[2];
+            if (tmpV3.x() > high.x())
+                high.setX( tmpV3.x() );
+            if (tmpV3.y() > high.y())
+                high.setY( tmpV3.y() );
+            if (tmpV3.z() > high.z())
+                high.setZ( tmpV3.z() );
         }
 
         computeCenter();
@@ -273,25 +281,25 @@ public class AABBox {
      */
     public final AABBox resize(final float x, final float y, final float z) {
         /** test low */
-        if (x < low[0]) {
-            low[0] = x;
+        if (x < low.x()) {
+            low.setX( x );
         }
-        if (y < low[1]) {
-            low[1] = y;
+        if (y < low.y()) {
+            low.setY( y );
         }
-        if (z < low[2]) {
-            low[2] = z;
+        if (z < low.z()) {
+            low.setZ( z );
         }
 
         /** test high */
-        if (x > high[0]) {
-            high[0] = x;
+        if (x > high.x()) {
+            high.setX( x );
         }
-        if (y > high[1]) {
-            high[1] = y;
+        if (y > high.y()) {
+            high.setY( y );
         }
-        if (z > high[2]) {
-            high[2] = z;
+        if (z > high.z()) {
+            high.setZ( z );
         }
 
         computeCenter();
@@ -320,6 +328,16 @@ public class AABBox {
     }
 
     /**
+     * Resize the AABBox to encapsulate the passed
+     * xyz-coordinates.
+     * @param xyz xyz-axis coordinate values
+     * @return this AABBox for chaining
+     */
+    public final AABBox resize(final Vec3f xyz) {
+        return resize(xyz.x(), xyz.y(), xyz.z());
+    }
+
+    /**
      * Check if the x & y coordinates are bounded/contained
      * by this AABBox
      * @param x  x-axis coordinate value
@@ -328,10 +346,10 @@ public class AABBox {
      * y belong to (low.y, high.y)
      */
     public final boolean contains(final float x, final float y) {
-        if(x<low[0] || x>high[0]){
+        if(x<low.x() || x>high.x()){
             return false;
         }
-        if(y<low[1]|| y>high[1]){
+        if(y<low.y()|| y>high.y()){
             return false;
         }
         return true;
@@ -347,13 +365,13 @@ public class AABBox {
      * y belong to (low.y, high.y) and  z belong to (low.z, high.z)
      */
     public final boolean contains(final float x, final float y, final float z) {
-        if(x<low[0] || x>high[0]){
+        if(x<low.x() || x>high.x()){
             return false;
         }
-        if(y<low[1]|| y>high[1]){
+        if(y<low.y()|| y>high.y()){
             return false;
         }
-        if(z<low[2] || z>high[2]){
+        if(z<low.z() || z>high.z()){
             return false;
         }
         return true;
@@ -390,13 +408,13 @@ public class AABBox {
     /**
      * Check if {@link Ray} intersects this bounding box.
      * <p>
-     * Versions uses the SAT[1], testing 6 axes.
+     * Versions uses the SAT.y(), testing 6 axes.
      * Original code for OBBs from MAGIC.
-     * Rewritten for AABBs and reorganized for early exits[2].
+     * Rewritten for AABBs and reorganized for early exits.z().
      * </p>
      * <pre>
-     * [1] SAT = Separating Axis Theorem
-     * [2] http://www.codercorner.com/RayAABB.cpp
+     * .y() SAT = Separating Axis Theorem
+     * .z() http://www.codercorner.com/RayAABB.cpp
      * </pre>
      * @param ray
      * @return
@@ -405,19 +423,19 @@ public class AABBox {
         // diff[XYZ] -> VectorUtil.subVec3(diff, ray.orig, center);
         //  ext[XYZ] -> extend VectorUtil.subVec3(ext, high, center);
 
-        final float dirX  = ray.dir[0];
-        final float diffX = ray.orig[0] - center[0];
-        final float extX  = high[0] - center[0];
+        final float dirX  = ray.dir.x();
+        final float diffX = ray.orig.x() - center.x();
+        final float extX  = high.x() - center.x();
         if( Math.abs(diffX) > extX && diffX*dirX >= 0f ) return false;
 
-        final float dirY  = ray.dir[1];
-        final float diffY = ray.orig[1] - center[1];
-        final float extY  = high[1] - center[1];
+        final float dirY  = ray.dir.y();
+        final float diffY = ray.orig.y() - center.y();
+        final float extY  = high.y() - center.y();
         if( Math.abs(diffY) > extY && diffY*dirY >= 0f ) return false;
 
-        final float dirZ  = ray.dir[2];
-        final float diffZ = ray.orig[2] - center[2];
-        final float extZ  = high[2] - center[2];
+        final float dirZ  = ray.dir.z();
+        final float diffZ = ray.orig.z() - center.z();
+        final float extZ  = high.z() - center.z();
         if( Math.abs(diffZ) > extZ && diffZ*dirZ >= 0f ) return false;
 
         final float absDirY = Math.abs(dirY);
@@ -442,7 +460,7 @@ public class AABBox {
      * or null if none exist.
      * <p>
      * <ul>
-     *  <li>Original code by Andrew Woo, from "Graphics Gems", Academic Press, 1990 [2]</li>
+     *  <li>Original code by Andrew Woo, from "Graphics Gems", Academic Press, 1990 .z()</li>
      *  <li>Optimized code by Pierre Terdiman, 2000 (~20-30% faster on my Celeron 500)</li>
      *  <li>Epsilon value added by Klaus Hartmann.</li>
      * </ul>
@@ -458,8 +476,8 @@ public class AABBox {
      * Report bugs: [email protected] (original author)
      * </p>
      * <pre>
-     * [1] http://www.codercorner.com/RayAABB.cpp
-     * [2] http://tog.acm.org/resources/GraphicsGems/gems/RayBox.c
+     * .y() http://www.codercorner.com/RayAABB.cpp
+     * .z() http://tog.acm.org/resources/GraphicsGems/gems/RayBox.c
      * </pre>
      * @param result vec3
      * @param ray
@@ -467,45 +485,45 @@ public class AABBox {
      * @param assumeIntersection if true, method assumes an intersection, i.e. by pre-checking via {@link #intersectsRay(Ray)}.
      *                           In this case method will not validate a possible non-intersection and just computes
      *                           coordinates.
-     * @param tmp1V3 temp vec3
-     * @param tmp2V3 temp vec3
-     * @param tmp3V3 temp vec3
      * @return float[3] result of intersection coordinates, or null if none exists
      */
-    public final float[] getRayIntersection(final float[] result, final Ray ray, final float epsilon,
-                                            final boolean assumeIntersection,
-                                            final float[] tmp1V3, final float[] tmp2V3, final float[] tmp3V3) {
+    public final Vec3f getRayIntersection(final Vec3f result, final Ray ray, final float epsilon,
+                                          final boolean assumeIntersection) {
         final float[] maxT = { -1f, -1f, -1f };
 
-        final float[] origin = ray.orig;
-        final float[] dir = ray.dir;
+        final Vec3f origin = ray.orig;
+        final Vec3f dir = ray.dir;
 
         boolean inside = true;
 
         // Find candidate planes.
         for(int i=0; i<3; i++) {
-            if(origin[i] < low[i]) {
-                result[i] = low[i];
+            final float origin_i = origin.get(i);
+            final float dir_i = dir.get(i);
+            final float low_i = low.get(i);
+            final float high_i = high.get(i);
+            if(origin_i < low_i) {
+                result.set(i, low_i);
                 inside    = false;
 
                 // Calculate T distances to candidate planes
-                if( 0 != Float.floatToIntBits(dir[i]) ) {
-                    maxT[i] = (low[i] - origin[i]) / dir[i];
+                if( 0 != Float.floatToIntBits(dir_i) ) {
+                    maxT[i] = (low_i - origin_i) / dir_i;
                 }
-            } else if(origin[i] > high[i]) {
-                result[i] = high[i];
+            } else if(origin_i > high_i) {
+                result.set(i, high_i);
                 inside    = false;
 
                 // Calculate T distances to candidate planes
-                if( 0 != Float.floatToIntBits(dir[i]) ) {
-                    maxT[i] = (high[i] - origin[i]) / dir[i];
+                if( 0 != Float.floatToIntBits(dir_i) ) {
+                    maxT[i] = (high_i - origin_i) / dir_i;
                 }
             }
         }
 
         // Ray origin inside bounding box
         if(inside) {
-            System.arraycopy(origin, 0, result, 0, 3);
+            result.set(origin);
             return result;
         }
 
@@ -530,22 +548,22 @@ public class AABBox {
             } */
             switch( whichPlane ) {
                 case 0:
-                    result[1] = origin[1] + maxT[whichPlane] * dir[1];
-                    if(result[1] < low[1] - epsilon || result[1] > high[1] + epsilon) { return null; }
-                    result[2] = origin[2] + maxT[whichPlane] * dir[2];
-                    if(result[2] < low[2] - epsilon || result[2] > high[2] + epsilon) { return null; }
+                    result.setY( origin.y() + maxT[whichPlane] * dir.y() );
+                    if(result.y() < low.y() - epsilon || result.y() > high.y() + epsilon) { return null; }
+                    result.setZ( origin.z() + maxT[whichPlane] * dir.z() );
+                    if(result.z() < low.z() - epsilon || result.z() > high.z() + epsilon) { return null; }
                     break;
                 case 1:
-                    result[0] = origin[0] + maxT[whichPlane] * dir[0];
-                    if(result[0] < low[0] - epsilon || result[0] > high[0] + epsilon) { return null; }
-                    result[2] = origin[2] + maxT[whichPlane] * dir[2];
-                    if(result[2] < low[2] - epsilon || result[2] > high[2] + epsilon) { return null; }
+                    result.setX( origin.x() + maxT[whichPlane] * dir.x() );
+                    if(result.x() < low.x() - epsilon || result.x() > high.x() + epsilon) { return null; }
+                    result.setZ( origin.z() + maxT[whichPlane] * dir.z() );
+                    if(result.z() < low.z() - epsilon || result.z() > high.z() + epsilon) { return null; }
                     break;
                 case 2:
-                    result[0] = origin[0] + maxT[whichPlane] * dir[0];
-                    if(result[0] < low[0] - epsilon || result[0] > high[0] + epsilon) { return null; }
-                    result[1] = origin[1] + maxT[whichPlane] * dir[1];
-                    if(result[1] < low[1] - epsilon || result[1] > high[1] + epsilon) { return null; }
+                    result.setX( origin.x() + maxT[whichPlane] * dir.x() );
+                    if(result.x() < low.x() - epsilon || result.x() > high.x() + epsilon) { return null; }
+                    result.setY( origin.y() + maxT[whichPlane] * dir.y() );
+                    if(result.y() < low.y() - epsilon || result.y() > high.y() + epsilon) { return null; }
                     break;
                 default:
                     throw new InternalError("XXX");
@@ -553,16 +571,16 @@ public class AABBox {
         } else {
             switch( whichPlane ) {
                 case 0:
-                    result[1] = origin[1] + maxT[whichPlane] * dir[1];
-                    result[2] = origin[2] + maxT[whichPlane] * dir[2];
+                    result.setY( origin.y() + maxT[whichPlane] * dir.y() );
+                    result.setZ( origin.z() + maxT[whichPlane] * dir.z() );
                     break;
                 case 1:
-                    result[0] = origin[0] + maxT[whichPlane] * dir[0];
-                    result[2] = origin[2] + maxT[whichPlane] * dir[2];
+                    result.setX( origin.x() + maxT[whichPlane] * dir.x() );
+                    result.setZ( origin.z() + maxT[whichPlane] * dir.z() );
                     break;
                 case 2:
-                    result[0] = origin[0] + maxT[whichPlane] * dir[0];
-                    result[1] = origin[1] + maxT[whichPlane] * dir[1];
+                    result.setX( origin.x() + maxT[whichPlane] * dir.x() );
+                    result.setY( origin.y() + maxT[whichPlane] * dir.y() );
                     break;
                 default:
                     throw new InternalError("XXX");
@@ -577,14 +595,14 @@ public class AABBox {
      * @return a float representing the size of the AABBox
      */
     public final float getSize() {
-        return VectorUtil.distVec3(low, high);
+        return low.dist(high);
     }
 
     /**
      * Get the Center of this AABBox
      * @return the xyz-coordinates of the center of the AABBox
      */
-    public final float[] getCenter() {
+    public final Vec3f getCenter() {
         return center;
     }
 
@@ -594,24 +612,17 @@ public class AABBox {
      * high and low is recomputed by scaling its distance to fixed center.
      * </p>
      * @param size a constant float value
-     * @param tmpV3 caller provided temporary 3-component vector
      * @return this AABBox for chaining
      * @see #scale2(float, float[])
      */
-    public final AABBox scale(final float size, final float[] tmpV3) {
-        tmpV3[0] = high[0] - center[0];
-        tmpV3[1] = high[1] - center[1];
-        tmpV3[2] = high[2] - center[2];
-
-        VectorUtil.scaleVec3(tmpV3, tmpV3, size); // in-place scale
-        VectorUtil.addVec3(high, center, tmpV3);
+    public final AABBox scale(final float size) {
+        final Vec3f tmp = new Vec3f();
+        tmp.set(high).sub(center).scale(size);
+        high.set(center).add(tmp);
 
-        tmpV3[0] = low[0] - center[0];
-        tmpV3[1] = low[1] - center[1];
-        tmpV3[2] = low[2] - center[2];
+        tmp.set(low).sub(center).scale(size);
+        low.set(center).add(tmp);
 
-        VectorUtil.scaleVec3(tmpV3, tmpV3, size); // in-place scale
-        VectorUtil.addVec3(low, center, tmpV3);
         return this;
     }
 
@@ -621,13 +632,12 @@ public class AABBox {
      * high and low is scaled and center recomputed.
      * </p>
      * @param size a constant float value
-     * @param tmpV3 caller provided temporary 3-component vector
      * @return this AABBox for chaining
      * @see #scale(float, float[])
      */
-    public final AABBox scale2(final float size, final float[] tmpV3) {
-        VectorUtil.scaleVec3(high, high, size); // in-place scale
-        VectorUtil.scaleVec3(low, low, size); // in-place scale
+    public final AABBox scale2(final float size) {
+        high.scale(size);
+        low.scale(size);
         computeCenter();
         return this;
     }
@@ -637,9 +647,9 @@ public class AABBox {
      * @param t the float[3] translation vector
      * @return this AABBox for chaining
      */
-    public final AABBox translate(final float[] t) {
-        VectorUtil.addVec3(low, low, t); // in-place translate
-        VectorUtil.addVec3(high, high, t); // in-place translate
+    public final AABBox translate(final Vec3f t) {
+        low.add(t);
+        high.add(t);
         computeCenter();
         return this;
     }
@@ -650,46 +660,46 @@ public class AABBox {
      * @return this AABBox for chaining
      */
     public final AABBox rotate(final Quaternion quat) {
-        quat.rotateVector(low, 0, low, 0);
-        quat.rotateVector(high, 0, high, 0);
+        quat.rotateVector(low, low);
+        quat.rotateVector(high, high);
         computeCenter();
         return this;
     }
 
     public final float getMinX() {
-        return low[0];
+        return low.x();
     }
 
     public final float getMinY() {
-        return low[1];
+        return low.y();
     }
 
     public final float getMinZ() {
-        return low[2];
+        return low.z();
     }
 
     public final float getMaxX() {
-        return high[0];
+        return high.x();
     }
 
     public final float getMaxY() {
-        return high[1];
+        return high.y();
     }
 
     public final float getMaxZ() {
-        return high[2];
+        return high.z();
     }
 
     public final float getWidth(){
-        return high[0] - low[0];
+        return high.x() - low.x();
     }
 
     public final float getHeight() {
-        return high[1] - low[1];
+        return high.y() - low.y();
     }
 
     public final float getDepth() {
-        return high[2] - low[2];
+        return high.z() - low.z();
     }
 
     @Override
@@ -701,29 +711,65 @@ public class AABBox {
             return false;
         }
         final AABBox other = (AABBox) obj;
-        return VectorUtil.isVec2Equal(low, 0, other.low, 0, FloatUtil.EPSILON) &&
-               VectorUtil.isVec3Equal(high, 0, other.high, 0, FloatUtil.EPSILON) ;
+        return low.isEqual(other.low) && high.isEqual(other.high);
     }
     @Override
     public final int hashCode() {
         throw new InternalError("hashCode not designed");
     }
 
+    public AABBox transform(final AABBox result, final float[/*16*/] mat4, final int mat4_off,
+                            final float[] vec3Tmp0, final float[] vec3Tmp1) {
+        result.reset();
+        FloatUtil.multMatrixVec3(mat4, mat4_off, low.get(vec3Tmp0), vec3Tmp1);
+        result.resize(vec3Tmp1);
+
+        FloatUtil.multMatrixVec3(mat4, mat4_off, high.get(vec3Tmp0), vec3Tmp1);
+        result.resize(vec3Tmp1);
+
+        result.computeCenter();
+        return result;
+    }
+
+    public AABBox transformMv(final AABBox result, final PMVMatrix pmv,
+                              final float[] vec3Tmp0, final float[] vec3Tmp1) {
+        result.reset();
+        pmv.multMvMatVec3f(low.get(vec3Tmp0), vec3Tmp1);
+        result.resize(vec3Tmp1);
+
+        pmv.multMvMatVec3f(high.get(vec3Tmp0), vec3Tmp1);
+        result.resize(vec3Tmp1);
+
+        result.computeCenter();
+        return result;
+    }
+
+    public AABBox transform(final AABBox result, final Matrix4f mat,
+                            final Vec3f vec3Tmp) {
+        result.reset();
+        result.resize( mat.mulVec3f(low, vec3Tmp) );
+
+        result.resize( mat.mulVec3f(high, vec3Tmp) );
+
+        result.computeCenter();
+        return result;
+    }
+
     /**
      * Assume this bounding box as being in object space and
      * compute the window bounding box.
      * <p>
      * If <code>useCenterZ</code> is <code>true</code>,
-     * only 4 {@link FloatUtil#mapObjToWinCoords(float, float, float, float[], int[], int, float[], int, float[], float[]) mapObjToWinCoords}
+     * only 4 {@link FloatUtil#mapObjToWin(float, float, float, float[], int[], float[], float[], float[]) mapObjToWinCoords}
      * operations are made on points [1..4] using {@link #getCenter()}'s z-value.
-     * Otherwise 8 {@link FloatUtil#mapObjToWinCoords(float, float, float, float[], int[], int, float[], int, float[], float[]) mapObjToWinCoords}
+     * Otherwise 8 {@link FloatUtil#mapObjToWin(float, float, float, float[], int[], float[], float[], float[]) mapObjToWinCoords}
      * operation on all 8 points are performed.
      * </p>
      * <pre>
-     *  [2] ------ [4]
+     *  .z() ------ [4]
      *   |          |
      *   |          |
-     *  [1] ------ [3]
+     *  .y() ------ [3]
      * </pre>
      * @param mat4PMv P x Mv matrix
      * @param view
@@ -736,42 +782,42 @@ public class AABBox {
     public AABBox mapToWindow(final AABBox result, final float[/*16*/] mat4PMv, final int[] view, final boolean useCenterZ,
                               final float[] vec3Tmp0, final float[] vec4Tmp1, final float[] vec4Tmp2) {
         {
-            // System.err.printf("AABBox.mapToWindow.0: view[%d, %d, %d, %d], this %s%n", view[0], view[1], view[2], view[3], toString());
-            final float objZ = useCenterZ ? center[2] : getMinZ();
-            FloatUtil.mapObjToWinCoords(getMinX(), getMinY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            // System.err.printf("AABBox.mapToWindow.p1: %f, %f, %f -> %f, %f, %f%n", getMinX(), getMinY(), objZ, vec3Tmp0[0], vec3Tmp0[1], vec3Tmp0[2]);
+            // System.err.printf("AABBox.mapToWindow.0: view[%d, %d, %d, %d], this %s%n", view.x(), view.y(), view.z(), view[3], toString());
+            final float objZ = useCenterZ ? center.z() : getMinZ();
+            FloatUtil.mapObjToWin(getMinX(), getMinY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            // System.err.printf("AABBox.mapToWindow.p1: %f, %f, %f -> %f, %f, %f%n", getMinX(), getMinY(), objZ, vec3Tmp0.x(), vec3Tmp0.y(), vec3Tmp0.z());
             // System.err.println("AABBox.mapToWindow.p1:");
             // System.err.println(FloatUtil.matrixToString(null, "  mat4PMv", "%10.5f", mat4PMv, 0, 4, 4, false /* rowMajorOrder */));
 
             result.reset();
-            result.resize(vec3Tmp0, 0);
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMinX(), getMaxY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            // System.err.printf("AABBox.mapToWindow.p2: %f, %f, %f -> %f, %f, %f%n", getMinX(), getMaxY(), objZ, vec3Tmp0[0], vec3Tmp0[1], vec3Tmp0[2]);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMinX(), getMaxY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            // System.err.printf("AABBox.mapToWindow.p2: %f, %f, %f -> %f, %f, %f%n", getMinX(), getMaxY(), objZ, vec3Tmp0.x(), vec3Tmp0.y(), vec3Tmp0.z());
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMaxX(), getMinY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            // System.err.printf("AABBox.mapToWindow.p3: %f, %f, %f -> %f, %f, %f%n", getMaxX(), getMinY(), objZ, vec3Tmp0[0], vec3Tmp0[1], vec3Tmp0[2]);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMaxX(), getMinY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            // System.err.printf("AABBox.mapToWindow.p3: %f, %f, %f -> %f, %f, %f%n", getMaxX(), getMinY(), objZ, vec3Tmp0.x(), vec3Tmp0.y(), vec3Tmp0.z());
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMaxX(), getMaxY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            // System.err.printf("AABBox.mapToWindow.p4: %f, %f, %f -> %f, %f, %f%n", getMaxX(), getMaxY(), objZ, vec3Tmp0[0], vec3Tmp0[1], vec3Tmp0[2]);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMaxX(), getMaxY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            // System.err.printf("AABBox.mapToWindow.p4: %f, %f, %f -> %f, %f, %f%n", getMaxX(), getMaxY(), objZ, vec3Tmp0.x(), vec3Tmp0.y(), vec3Tmp0.z());
+            result.resize(vec3Tmp0);
         }
 
         if( !useCenterZ ) {
             final float objZ = getMaxZ();
-            FloatUtil.mapObjToWinCoords(getMinX(), getMinY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMinX(), getMinY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMinX(), getMaxY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMinX(), getMaxY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMaxX(), getMinY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMaxX(), getMinY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            result.resize(vec3Tmp0);
 
-            FloatUtil.mapObjToWinCoords(getMaxX(), getMaxY(), objZ, mat4PMv, view, 0, vec3Tmp0, 0, vec4Tmp1, vec4Tmp2);
-            result.resize(vec3Tmp0, 0);
+            FloatUtil.mapObjToWin(getMaxX(), getMaxY(), objZ, mat4PMv, view, vec3Tmp0, vec4Tmp1, vec4Tmp2);
+            result.resize(vec3Tmp0);
         }
         if( DEBUG ) {
             System.err.printf("AABBox.mapToWindow: view[%d, %d], this %s -> %s%n", view[0], view[1], toString(), result.toString());
@@ -782,7 +828,6 @@ public class AABBox {
     @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]+" ]";
+               ", box "+low+" .. "+high+", ctr "+center+" ]";
     }
 }
diff --git a/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
index 8b0fa559e..4d098cb72 100644
--- a/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
+++ b/src/jogl/classes/com/jogamp/opengl/math/geom/Frustum.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2010 JogAmp Community. All rights reserved.
+ * Copyright 2010-2023 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:
@@ -29,9 +29,11 @@ package com.jogamp.opengl.math.geom;
 
 import jogamp.common.os.PlatformPropsImpl;
 
-import com.jogamp.common.os.Platform;
 import com.jogamp.opengl.math.FloatUtil;
 import com.jogamp.opengl.math.FovHVHalves;
+import com.jogamp.opengl.math.Matrix4f;
+import com.jogamp.opengl.math.Vec3f;
+import com.jogamp.opengl.math.geom.Frustum.FovDesc;
 
 /**
  * Providing frustum {@link #getPlanes() planes} derived by different inputs
@@ -103,6 +105,7 @@ public class Frustum {
             this.zNear = zNear;
             this.zFar = zFar;
         }
+        @Override
         public final String toString() {
             return "FrustumFovDesc["+fovhv.toStringInDegrees()+", Z["+zNear+" - "+zFar+"]]";
         }
@@ -134,7 +137,7 @@ public class Frustum {
 	 */
     public static class Plane {
         /** Normal of the plane */
-        public final float[] n = new float[3];
+        public final Vec3f n = new Vec3f();
 
         /** Distance to origin */
         public float d;
@@ -155,17 +158,22 @@ public class Frustum {
          * </p>
          **/
         public final float distanceTo(final float x, final float y, final float z) {
-            return n[0] * x + n[1] * y + n[2] * z + d;
+            return n.x() * x + n.y() * y + n.z() * z + d;
         }
 
         /** Return distance of plane to given point, see {@link #distanceTo(float, float, float)}. */
         public final float distanceTo(final float[] p) {
-            return n[0] * p[0] + n[1] * p[1] + n[2] * p[2] + d;
+            return n.x() * p[0] + n.y() * p[1] + n.z() * p[2] + d;
+        }
+
+        /** Return distance of plane to given point, see {@link #distanceTo(float, float, float)}. */
+        public final float distanceTo(final Vec3f p) {
+            return n.x() * p.x() + n.y() * p.y() + n.z() * p.z() + d;
         }
 
         @Override
         public String toString() {
-            return "Plane[ [ " + n[0] + ", " + n[1] + ", " + n[2] + " ], " + d + "]";
+            return "Plane[ [ " + n + " ], " + d + "]";
         }
     }
 
@@ -221,9 +229,9 @@ public class Frustum {
      * Operation Details:
      * <ul>
      *   <li>The given {@link FovDesc} will be transformed
-     *       into the given float[16] as a perspective matrix (column major order) first,
-     *       see {@link FloatUtil#makePerspective(float[], int, boolean, FovHVHalves, float, float)}.</li>
-     *   <li>Then the float[16] perspective matrix is used to {@link #updateByPMV(float[], int)} this instance.</li>
+     *       into the given perspective matrix (column major order) first,
+     *       see {@link Matrix4f#setToPerspective(FovHVHalves, float, float)}.</li>
+     *   <li>Then the perspective matrix is used to {@link Matrix4f#updateFrustumPlanes(Frustum)} this instance.</li>
      * </ul>
      * </p>
      * <p>
@@ -232,17 +240,15 @@ public class Frustum {
      * </p>
      *
      * @param m 4x4 matrix in column-major order (also result)
-     * @param m_offset offset in given array <i>m</i>, i.e. start of the 4x4 matrix
-     * @param initM if true, given matrix will be initialized w/ identity matrix,
-     *              otherwise only the frustum fields are set.
      * @param fovDesc {@link Frustum} {@link FovDesc}
      * @return given matrix for chaining
-     * @see FloatUtil#makePerspective(float[], int, boolean, FovHVHalves, float, float)
+     * @see Matrix4f#setToPerspective(FovHVHalves, float, float)
+     * @see Matrix4f#updateFrustumPlanes(Frustum)
+     * @see Matrix4f#getFrustum(Frustum, FovDesc)
      */
-    public float[] updateByFovDesc(final float[] m, final int m_offset, final boolean initM,
-                                   final FovDesc fovDesc) {
-        FloatUtil.makePerspective(m, m_offset, initM, fovDesc.fovhv, fovDesc.zNear, fovDesc.zFar);
-        updateByPMV(m, 0);
+    public Matrix4f updateByFovDesc(final Matrix4f m, final FovDesc fovDesc) {
+        m.setToPerspective(fovDesc.fovhv, fovDesc.zNear, fovDesc.zFar);
+        m.updateFrustumPlanes(this);
         return m;
     }
 
@@ -259,10 +265,10 @@ public class Frustum {
         // Left:   a = m30 + m00, b = m31 + m01, c = m32 + m02, d = m33 + m03  - [0..3] column-major
         {
             final Plane p = planes[LEFT];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 0 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 0 + 1 * 4 ],
+                     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 ];
         }
 
@@ -270,10 +276,10 @@ public class Frustum {
         // Right:  a = m30 - m00, b = m31 - m01, c = m32 - m02, d = m33 - m03  - [0..3] column-major
         {
             final Plane p = planes[RIGHT];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 0 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 0 + 1 * 4 ],
+                     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 ];
         }
 
@@ -281,10 +287,10 @@ public class Frustum {
         // Bottom: a = m30 + m10, b = m31 + m11, c = m32 + m12, d = m33 + m13  - [0..3] column-major
         {
             final Plane p = planes[BOTTOM];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 1 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 1 + 1 * 4 ],
+                     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 ];
         }
 
@@ -292,10 +298,10 @@ public class Frustum {
         // Top:   a = m30 - m10, b = m31 - m11, c = m32 - m12, d = m33 - m13  - [0..3] column-major
         {
             final Plane p = planes[TOP];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 1 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 1 + 1 * 4 ],
+                     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 ];
         }
 
@@ -303,10 +309,10 @@ public class Frustum {
         // Near:  a = m30 + m20, b = m31 + m21, c = m32 + m22, d = m33 + m23  - [0..3] column-major
         {
             final Plane p = planes[NEAR];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] + pmv[ pmv_off + 2 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] + pmv[ pmv_off + 2 + 1 * 4 ],
+                     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 ];
         }
 
@@ -314,38 +320,35 @@ public class Frustum {
         // Far:   a = m30 - m20, b = m31 - m21, c = m32 + m22, d = m33 + m23  - [0..3] column-major
         {
             final Plane p = planes[FAR];
-            final float[] p_n = p.n;
-            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 ];
+            final Vec3f p_n = p.n;
+            p_n.set( pmv[ pmv_off + 3 + 0 * 4 ] - pmv[ pmv_off + 2 + 0 * 4 ],
+                     pmv[ pmv_off + 3 + 1 * 4 ] - pmv[ pmv_off + 2 + 1 * 4 ],
+                     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;
+            final Vec3f p_n = p.n;
+            final float invLen = 1f / p_n.length();
+            p_n.scale(invLen);
+            p.d *= invLen;
         }
     }
 
 	private static final boolean isOutsideImpl(final Plane p, final 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 ) {
+	    final Vec3f lo = box.getLow();
+	    final Vec3f hi = box.getHigh();
+
+		if ( p.distanceTo(lo.x(), lo.y(), lo.z()) > 0.0f ||
+		     p.distanceTo(hi.x(), lo.y(), lo.z()) > 0.0f ||
+		     p.distanceTo(lo.x(), hi.y(), lo.z()) > 0.0f ||
+		     p.distanceTo(hi.x(), hi.y(), lo.z()) > 0.0f ||
+		     p.distanceTo(lo.x(), lo.y(), hi.z()) > 0.0f ||
+		     p.distanceTo(hi.x(), lo.y(), hi.z()) > 0.0f ||
+		     p.distanceTo(lo.x(), hi.y(), hi.z()) > 0.0f ||
+		     p.distanceTo(hi.x(), hi.y(), hi.z()) > 0.0f ) {
 			return false;
 		}
 		return true;
diff --git a/src/jogl/classes/com/jogamp/opengl/util/PMVMatrix.java b/src/jogl/classes/com/jogamp/opengl/util/PMVMatrix.java
index f8f1d3930..1aa305c2e 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/PMVMatrix.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/PMVMatrix.java
@@ -46,6 +46,7 @@ import jogamp.common.os.PlatformPropsImpl;
 import com.jogamp.common.nio.Buffers;
 import com.jogamp.common.util.FloatStack;
 import com.jogamp.opengl.math.FloatUtil;
+import com.jogamp.opengl.math.Matrix4f;
 import com.jogamp.opengl.math.Quaternion;
 import com.jogamp.opengl.math.Ray;
 import com.jogamp.opengl.math.geom.AABBox;
@@ -208,9 +209,9 @@ public final class PMVMatrix implements GLMatrixFunc {
           // Mvit Modelview-Inverse-Transpose
           matrixArray = new float[5*16];
 
-          mP_offset   = 0*16;
-          mMv_offset  = 1*16;
-          mTex_offset = 4*16;
+          // mP_offset   = 0*16;
+          // mMv_offset  = 1*16;
+          // mTex_offset = 4*16;
 
           matrixPMvMvit = Buffers.slice2Float(matrixArray,  0*16, 4*16);  // P + Mv + Mvi + Mvit
           matrixPMvMvi  = Buffers.slice2Float(matrixArray,  0*16, 3*16);  // P + Mv + Mvi
@@ -241,6 +242,8 @@ public final class PMVMatrix implements GLMatrixFunc {
      * Issues {@link #glLoadIdentity()} on all matrices,
      * i.e. {@link GLMatrixFunc#GL_MODELVIEW GL_MODELVIEW}, {@link GLMatrixFunc#GL_PROJECTION GL_PROJECTION} or {@link GL#GL_TEXTURE GL_TEXTURE}
      * and resets all internal states.
+     *
+     * Leaves {@link GLMatrixFunc#GL_MODELVIEW GL_MODELVIEW} the active matrix mode.
      */
     public final void reset() {
         FloatUtil.makeIdentity(matrixArray, mMv_offset);
@@ -403,7 +406,6 @@ public final class PMVMatrix implements GLMatrixFunc {
         }
     }
 
-
     /**
      * Multiplies the {@link #glGetPMatrixf() P} and {@link #glGetMvMatrixf() Mv} matrix, i.e.
      * <pre>
@@ -432,6 +434,71 @@ public final class PMVMatrix implements GLMatrixFunc {
         return mat4MvP;
     }
 
+    /**
+     * v_out = Mv * v_in
+     * @param v_in float[4] input vector
+     * @param v_out float[4] output vector
+     */
+    public final void multMvMatVec4f(final float[/*4*/] v_in, final float[/*4*/] v_out) {
+        FloatUtil.multMatrixVec(matrixArray, mMv_offset, v_in, v_out);
+    }
+
+    /**
+     * v_out = Mv * v_in
+     *
+     * Affine 3f-vector transformation by 4x4 matrix, see {@link FloatUtil#multMatrixVec3(float[], int, float[], float[])}.
+     *
+     * @param v_in float[3] input vector
+     * @param v_out float[3] output vector
+     */
+    public final void multMvMatVec3f(final float[/*3*/] v_in, final float[/*3*/] v_out) {
+        FloatUtil.multMatrixVec3(matrixArray, mMv_offset, v_in, v_out);
+    }
+
+    /**
+     * v_out = P * v_in
+     * @param v_in float[4] input vector
+     * @param v_out float[4] output vector
+     */
+    public final void multPMatVec4f(final float[/*4*/] v_in, final float[/*4*/] v_out) {
+        FloatUtil.multMatrixVec(matrixArray, v_in, v_out); // mP_offset := 0
+    }
+
+    /**
+     * v_out = P * v_in
+     *
+     * Affine 3f-vector transformation by 4x4 matrix, see {@link FloatUtil#multMatrixVec3(float[], int, float[], float[])}.
+     *
+     * @param v_in float[3] input vector
+     * @param v_out float[3] output vector
+     */
+    public final void multPMatVec3f(final float[/*3*/] v_in, final float[/*3*/] v_out) {
+        FloatUtil.multMatrixVec3(matrixArray, v_in, v_out); // mP_offset := 0
+    }
+
+    /**
+     * v_out = P * Mv * v_in
+     * @param v_in float[4] input vector
+     * @param v_out float[4] output vector
+     */
+    public final void multPMvMatVec4f(final float[/*4*/] v_in, final float[/*4*/] v_out) {
+        FloatUtil.multMatrixVec(matrixArray, mMv_offset, v_in, mat4Tmp1);
+        FloatUtil.multMatrixVec(matrixArray, mat4Tmp1, v_out); // mP_offset := 0
+    }
+
+    /**
+     * v_out = P * Mv * v_in
+     *
+     * Affine 3f-vector transformation by 4x4 matrix, see {@link FloatUtil#multMatrixVec3(float[], int, float[], float[])}.
+     *
+     * @param v_in float[3] input vector
+     * @param v_out float[3] output vector
+     */
+    public final void multPMvMatVec3f(final float[/*3*/] v_in, final float[/*3*/] v_out) {
+        FloatUtil.multMatrixVec3(matrixArray, mMv_offset, v_in, mat4Tmp1);
+        FloatUtil.multMatrixVec3(matrixArray, mat4Tmp1, v_out); // mP_offset := 0
+    }
+
     //
     // GLMatrixFunc implementation
     //
@@ -534,6 +601,24 @@ public final class PMVMatrix implements GLMatrixFunc {
     }
 
     /**
+     * Load the current matrix with the values of the given {@link Matrix4f}.
+     */
+    public final void glLoadMatrixf(final Matrix4f m) {
+        if(matrixMode==GL_MODELVIEW) {
+            m.get(matrixArray, mMv_offset);
+            dirtyBits |= DIRTY_INVERSE_MODELVIEW | DIRTY_INVERSE_TRANSPOSED_MODELVIEW | DIRTY_FRUSTUM ;
+            modifiedBits |= MODIFIED_MODELVIEW;
+        } else if(matrixMode==GL_PROJECTION) {
+            m.get(matrixArray, mP_offset);
+            dirtyBits |= DIRTY_FRUSTUM ;
+            modifiedBits |= MODIFIED_PROJECTION;
+        } else if(matrixMode==GL.GL_TEXTURE) {
+            m.get(matrixArray, mTex_offset);
+            modifiedBits |= MODIFIED_TEXTURE;
+        }
+    }
+
+    /**
      * Load the current matrix with the values of the given {@link Quaternion}'s rotation {@link Quaternion#toMatrix(float[], int) matrix representation}.
      */
     public final void glLoadMatrix(final Quaternion quat) {
@@ -633,6 +718,21 @@ public final class PMVMatrix implements GLMatrixFunc {
         }
     }
 
+    public final void glMultMatrixf(final Matrix4f m) {
+        if(matrixMode==GL_MODELVIEW) {
+            new Matrix4f(matrixArray, mMv_offset).mul(m).get(matrixArray, mMv_offset);
+            dirtyBits |= DIRTY_INVERSE_MODELVIEW | DIRTY_INVERSE_TRANSPOSED_MODELVIEW | DIRTY_FRUSTUM ;
+            modifiedBits |= MODIFIED_MODELVIEW;
+        } else if(matrixMode==GL_PROJECTION) {
+            new Matrix4f(matrixArray, mP_offset).mul(m).get(matrixArray, mP_offset);
+            dirtyBits |= DIRTY_FRUSTUM ;
+            modifiedBits |= MODIFIED_PROJECTION;
+        } else if(matrixMode==GL.GL_TEXTURE) {
+            new Matrix4f(matrixArray, mTex_offset).mul(m).get(matrixArray, mTex_offset);
+            modifiedBits |= MODIFIED_TEXTURE;
+        }
+    }
+
     @Override
     public final void glTranslatef(final float x, final float y, final float z) {
         glMultMatrixf(FloatUtil.makeTranslation(matrixTxSx, false, x, y, z), 0);
@@ -645,7 +745,7 @@ public final class PMVMatrix implements GLMatrixFunc {
 
     @Override
     public final void glRotatef(final float ang_deg, final float x, final float y, final float z) {
-        glMultMatrixf(FloatUtil.makeRotationAxis(mat4Tmp1, 0, ang_deg * FloatUtil.PI / 180.0f, x, y, z, mat4Tmp2), 0);
+        glMultMatrixf(FloatUtil.makeRotationAxis(mat4Tmp1, 0, FloatUtil.adegToRad(ang_deg), x, y, z, mat4Tmp2), 0);
     }
 
     /**
@@ -728,7 +828,7 @@ public final class PMVMatrix implements GLMatrixFunc {
     public final boolean gluProject(final float objx, final float objy, final float objz,
                                     final int[] viewport, final int viewport_offset,
                                     final float[] win_pos, final int win_pos_offset ) {
-        return FloatUtil.mapObjToWinCoords(objx, objy, objz,
+        return FloatUtil.mapObjToWin(objx, objy, objz,
                           matrixArray, mMv_offset,
                           matrixArray, mP_offset,
                           viewport, viewport_offset,
@@ -754,7 +854,7 @@ public final class PMVMatrix implements GLMatrixFunc {
     public final boolean gluUnProject(final float winx, final float winy, final float winz,
                                       final int[] viewport, final int viewport_offset,
                                       final float[] obj_pos, final int obj_pos_offset) {
-        return FloatUtil.mapWinToObjCoords(winx, winy, winz,
+        return FloatUtil.mapWinToObj(winx, winy, winz,
                                            matrixArray, mMv_offset,
                                            matrixArray, mP_offset,
                                            viewport, viewport_offset,
@@ -788,7 +888,7 @@ public final class PMVMatrix implements GLMatrixFunc {
                                  final int[] viewport, final int viewport_offset,
                                  final float near, final float far,
                                  final float[] obj_pos, final int obj_pos_offset ) {
-        return FloatUtil.mapWinToObjCoords(winx, winy, winz, clipw,
+        return FloatUtil.mapWinToObj4(winx, winy, winz, clipw,
                 matrixArray, mMv_offset,
                 matrixArray, mP_offset,
                 viewport, viewport_offset,
@@ -842,14 +942,12 @@ public final class PMVMatrix implements GLMatrixFunc {
      * @return true if successful, otherwise false (failed to invert matrix, or becomes z is infinity)
      */
     public final boolean gluUnProjectRay(final float winx, final float winy, final float winz0, final float winz1,
-                                         final int[] viewport, final int viewport_offset,
+                                         final int[] viewport,
                                          final Ray ray) {
         return FloatUtil.mapWinToRay(winx, winy, winz0, winz1,
                 matrixArray, mMv_offset,
-                matrixArray, mP_offset,
-                viewport, viewport_offset,
-                ray,
-                mat4Tmp1, mat4Tmp2, mat4Tmp3);
+                matrixArray, mP_offset, viewport,
+                ray, mat4Tmp1, mat4Tmp2, mat4Tmp3);
     }
 
     public StringBuilder toString(StringBuilder sb, final String f) {
@@ -1066,8 +1164,10 @@ public final class PMVMatrix implements GLMatrixFunc {
         return res;
     }
 
+    private static final int mP_offset   = 0*16;
+    private static final int mMv_offset  = 1*16;
+    private static final int mTex_offset = 4*16;
     private final float[] matrixArray;
-    private final int mP_offset, mMv_offset, mTex_offset;
     private final FloatBuffer matrixPMvMvit, matrixPMvMvi, matrixPMv, matrixP, matrixTex, matrixMv, matrixMvi, matrixMvit;
     private final float[] matrixTxSx;
     private final float[] mat4Tmp1, mat4Tmp2, mat4Tmp3;
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/EyeParameter.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/EyeParameter.java
index 43a6cfc58..e0f465da7 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/EyeParameter.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/EyeParameter.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2014 JogAmp Community. All rights reserved.
+ * Copyright 2014-2023 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:
@@ -28,6 +28,7 @@
 package com.jogamp.opengl.util.stereo;
 
 import com.jogamp.opengl.math.FovHVHalves;
+import com.jogamp.opengl.math.Vec3f;
 
 /**
  * Constant single eye parameter of the viewer, relative to its {@link ViewerPose}.
@@ -36,8 +37,8 @@ public final class EyeParameter {
     /** Eye number, <code>0</code> for the left eye and <code>1</code> for the right eye. */
     public final int number;
 
-    /** float[3] eye position vector used to define eye height in meter relative to <i>actor</i>. */
-    public final float[] positionOffset;
+    /** eye position vector used to define eye height in meter relative to <i>actor</i>. */
+    public final Vec3f positionOffset;
 
     /** Field of view in both directions, may not be centered, either {@link FovHVHalves#inTangents} or radians. */
     public final FovHVHalves fovhv;
@@ -51,18 +52,18 @@ public final class EyeParameter {
     /** Z-axis eye relief in meter. */
     public final float eyeReliefZ;
 
-    public EyeParameter(final int number, final float[] positionOffset, final FovHVHalves fovhv,
+    public EyeParameter(final int number, final Vec3f positionOffset, final FovHVHalves fovhv,
                         final float distNoseToPupil, final float verticalDelta, final float eyeRelief) {
         this.number = number;
-        this.positionOffset = new float[3];
-        System.arraycopy(positionOffset, 0, this.positionOffset, 0, 3);
+        this.positionOffset = new Vec3f(positionOffset);
         this.fovhv = fovhv;
         this.distNoseToPupilX = distNoseToPupil;
         this.distMiddleToPupilY = verticalDelta;
         this.eyeReliefZ = eyeRelief;
     }
+    @Override
     public final String toString() {
-        return "EyeParam[num "+number+", posOff["+positionOffset[0]+", "+positionOffset[1]+", "+positionOffset[2]+"], "+fovhv+
-                      ", distPupil[noseX "+distNoseToPupilX+", middleY "+distMiddleToPupilY+", reliefZ "+eyeReliefZ+"]]";
+        return "EyeParam[num "+number+", posOff["+positionOffset+"], "+fovhv+
+               ", distPupil[noseX "+distNoseToPupilX+", middleY "+distMiddleToPupilY+", reliefZ "+eyeReliefZ+"]]";
     }
 }
\ No newline at end of file
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/LocationSensorParameter.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/LocationSensorParameter.java
index b795927cd..6294adee1 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/LocationSensorParameter.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/LocationSensorParameter.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2015 JogAmp Community. All rights reserved.
+ * Copyright 2015-2023 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:
@@ -27,6 +27,7 @@
  */
 package com.jogamp.opengl.util.stereo;
 
+import com.jogamp.opengl.math.Matrix4f;
 import com.jogamp.opengl.math.geom.Frustum;
 
 /**
@@ -38,13 +39,15 @@ public final class LocationSensorParameter {
     /** The {@link Frustum}'s {@link Frustum.FovDesc} description of the location sensor. */
     public final Frustum.FovDesc frustumDesc;
     /** The {@link Frustum}'s float[16] projection matrix of the location sensor. */
-    public final float[] frustumProjMat;
+    public final Matrix4f frustumProjMat;
 
     public LocationSensorParameter(final Frustum.FovDesc fovDesc) {
         this.frustumDesc = fovDesc;
         this.frustum = new Frustum();
-        this.frustumProjMat = frustum.updateByFovDesc(new float[16], 0, true, fovDesc);
+        this.frustumProjMat = frustum.updateByFovDesc(new Matrix4f(), fovDesc);
     }
+
+    @Override
     public final String toString() {
         return "LocationSensor["+frustumDesc+"]";
     }
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoDevice.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoDevice.java
index b6112650a..85e752302 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoDevice.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoDevice.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2014 JogAmp Community. All rights reserved.
+ * Copyright 2014-2023 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:
@@ -33,6 +33,7 @@ import com.jogamp.nativewindow.util.PointImmutable;
 import jogamp.opengl.Debug;
 
 import com.jogamp.opengl.math.FovHVHalves;
+import com.jogamp.opengl.math.Vec3f;
 
 /**
  * Interface describing a native stereoscopic device
@@ -94,7 +95,7 @@ public interface StereoDevice {
     public int getRequiredRotation();
 
     /**
-     * Return the device default eye position offset for {@link #createRenderer(int, int, float[], FovHVHalves[], float)}.
+     * Return the device default eye position offset for {@link #createRenderer(int, int, Vec3f, FovHVHalves[], float)}.
      * <p>
      * Result is an array of float values for
      * <ul>
@@ -105,7 +106,7 @@ public interface StereoDevice {
      * </p>
      * @return
      */
-    public float[] getDefaultEyePositionOffset();
+    public Vec3f getDefaultEyePositionOffset();
 
     /**
      * Returns the device default {@link FovHVHalves} for all supported eyes
@@ -198,7 +199,7 @@ public interface StereoDevice {
      * Returns the supported distortion compensation of the {@link StereoDeviceRenderer},
      * e.g. {@link StereoDeviceRenderer#DISTORTION_BARREL}, {@link StereoDeviceRenderer#DISTORTION_CHROMATIC}, etc.
      * @see StereoDeviceRenderer#getDistortionBits()
-     * @see #createRenderer(int, int, float[], FovHVHalves[], float, int)
+     * @see #createRenderer(int, int, Vec3f, FovHVHalves[], float, int)
      * @see #getRecommendedDistortionBits()
      * @see #getMinimumDistortionBits()
      */
@@ -212,7 +213,7 @@ public interface StereoDevice {
      * User shall use the recommended distortion compensation to achieve a distortion free view.
      * </p>
      * @see StereoDeviceRenderer#getDistortionBits()
-     * @see #createRenderer(int, int, float[], FovHVHalves[], float, int)
+     * @see #createRenderer(int, int, Vec3f, FovHVHalves[], float, int)
      * @see #getSupportedDistortionBits()
      * @see #getMinimumDistortionBits()
      */
@@ -227,7 +228,7 @@ public interface StereoDevice {
      * @see #getSupportedDistortionBits()
      * @see #getRecommendedDistortionBits()
      * @see StereoDeviceRenderer#getDistortionBits()
-     * @see #createRenderer(int, int, float[], FovHVHalves[], float, int)
+     * @see #createRenderer(int, int, Vec3f, FovHVHalves[], float, int)
      */
     public int getMinimumDistortionBits();
 
@@ -245,6 +246,6 @@ public interface StereoDevice {
      * @return
      */
     public StereoDeviceRenderer createRenderer(final int distortionBits,
-                                               final int textureCount, final float[] eyePositionOffset,
+                                               final int textureCount, final Vec3f eyePositionOffset,
                                                final FovHVHalves[] eyeFov, final float pixelsPerDisplayPixel, final int textureUnit);
 }
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoUtil.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoUtil.java
index b6f76a343..63cb3e1e7 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoUtil.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/StereoUtil.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2014 JogAmp Community. All rights reserved.
+ * Copyright 2014-2023 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:
@@ -27,9 +27,9 @@
  */
 package com.jogamp.opengl.util.stereo;
 
-import com.jogamp.opengl.math.FloatUtil;
+import com.jogamp.opengl.math.Matrix4f;
 import com.jogamp.opengl.math.Quaternion;
-import com.jogamp.opengl.math.VectorUtil;
+import com.jogamp.opengl.math.Vec3f;
 import com.jogamp.opengl.util.CustomGLEventListener;
 import com.jogamp.opengl.util.stereo.StereoDeviceRenderer.Eye;
 
@@ -137,24 +137,24 @@ public class StereoUtil {
      * @param eye
      * @param zNear frustum near value
      * @param zFar frustum far value
-     * @param mat4Projection float[16] projection matrix result
-     * @param mat4Modelview float[16] modelview matrix result
+     * @param mat4Projection projection matrix result
+     * @param mat4Modelview modelview matrix result
      */
     public static void getSBSUpstreamPMV(final ViewerPose viewerPose, final Eye eye,
                                          final float zNear, final float zFar,
-                                         final float[] mat4Projection, final float[] mat4Modelview) {
-        final float[] mat4Tmp1 = new float[16];
-        final float[] mat4Tmp2 = new float[16];
-        final float[] vec3Tmp1 = new float[3];
-        final float[] vec3Tmp2 = new float[3];
-        final float[] vec3Tmp3 = new float[3];
+                                         final Matrix4f mat4Projection, final Matrix4f mat4Modelview) {
+        final Matrix4f mat4Tmp1 = new Matrix4f();
+        final Matrix4f mat4Tmp2 = new Matrix4f();
+        final Vec3f vec3Tmp1 = new Vec3f();
+        final Vec3f vec3Tmp2 = new Vec3f();
+        final Vec3f vec3Tmp3 = new Vec3f();
 
         final EyeParameter eyeParam = eye.getEyeParameter();
 
         //
         // Projection
         //
-        FloatUtil.makePerspective(mat4Projection, 0, true, eyeParam.fovhv, zNear, zFar);
+        mat4Projection.setToPerspective(eyeParam.fovhv, zNear, zFar);
 
         //
         // Modelview
@@ -162,21 +162,17 @@ public class StereoUtil {
         final Quaternion rollPitchYaw = new Quaternion();
         // private final float eyeYaw = FloatUtil.PI; // 180 degrees in radians
         // rollPitchYaw.rotateByAngleY(eyeYaw);
-        final float[] shiftedEyePos = rollPitchYaw.rotateVector(vec3Tmp1, 0, viewerPose.position, 0);
-        VectorUtil.addVec3(shiftedEyePos, shiftedEyePos, eyeParam.positionOffset);
+        final Vec3f shiftedEyePos = rollPitchYaw.rotateVector(viewerPose.position, vec3Tmp1).add(eyeParam.positionOffset);
 
         rollPitchYaw.mult(viewerPose.orientation);
-        final float[] up = rollPitchYaw.rotateVector(vec3Tmp2, 0, VectorUtil.VEC3_UNIT_Y, 0);
-        final float[] forward = rollPitchYaw.rotateVector(vec3Tmp3, 0, VectorUtil.VEC3_UNIT_Z_NEG, 0);
-        final float[] center = VectorUtil.addVec3(forward, shiftedEyePos, forward);
+        final Vec3f up = rollPitchYaw.rotateVector(Vec3f.UNIT_Y, vec3Tmp2);
+        final Vec3f forward = rollPitchYaw.rotateVector(Vec3f.UNIT_Z_NEG, vec3Tmp3); // -> center
+        final Vec3f center = forward.add(shiftedEyePos);
 
-        final float[] mLookAt = FloatUtil.makeLookAt(mat4Tmp2, 0, shiftedEyePos, 0, center, 0, up, 0, mat4Tmp1);
-        final float[] mViewAdjust = FloatUtil.makeTranslation(mat4Modelview, true,
-                                                              eyeParam.distNoseToPupilX,
-                                                              eyeParam.distMiddleToPupilY,
-                                                              eyeParam.eyeReliefZ);
-
-        /* mat4Modelview = */ FloatUtil.multMatrix(mViewAdjust, mLookAt);
+        final Matrix4f mLookAt = mat4Tmp2.setToLookAt(shiftedEyePos, center, up, mat4Tmp1);
+        mat4Modelview.mul( mat4Tmp1.setToTranslation( eyeParam.distNoseToPupilX,
+                                                      eyeParam.distMiddleToPupilY,
+                                                      eyeParam.eyeReliefZ ), mLookAt);
     }
 
 }
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/ViewerPose.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/ViewerPose.java
index 10ee4c994..5d2cf925c 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/ViewerPose.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/ViewerPose.java
@@ -1,5 +1,5 @@
 /**
- * Copyright 2014 JogAmp Community. All rights reserved.
+ * Copyright 2014-2023 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:
@@ -28,13 +28,14 @@
 package com.jogamp.opengl.util.stereo;
 
 import com.jogamp.opengl.math.Quaternion;
+import com.jogamp.opengl.math.Vec3f;
 
 /**
  * {@link #position} and {@link #orientation} of viewer.
  */
 public final class ViewerPose {
     /**
-     * float[3] position of viewer in meter.
+     * position of viewer in meter.
      * <p>
      * Apply the following to resolve the actual eye position:
      * <ul>
@@ -43,13 +44,13 @@ public final class ViewerPose {
      * </ul>
      * </p>
      */
-    public final float[] position;
+    public final Vec3f position;
 
     /** Orientation of viewer. */
     public final Quaternion orientation;
 
     public ViewerPose() {
-        this.position = new float[3];
+        this.position = new Vec3f();
         this.orientation = new Quaternion();
     }
     public ViewerPose(final float[] position, final Quaternion orientation) {
@@ -64,11 +65,14 @@ public final class ViewerPose {
     }
     /** Set position and orientation of this instance. */
     public final void setPosition(final float posX, final float posY, final float posZ) {
-        position[0] = posX;
-        position[1] = posY;
-        position[2] = posZ;
+        position.set( posX, posY, posZ );
     }
+    /** Set position and orientation of this instance. */
+    public final void setPosition(final Vec3f pos) {
+        position.set( pos );
+    }
+    @Override
     public final String toString() {
-        return "ViewerPose[pos["+position[0]+", "+position[1]+", "+position[2]+"], "+orientation+"]";
+        return "ViewerPose[pos["+position+"], "+orientation+"]";
     }
 }
\ No newline at end of file
diff --git a/src/jogl/classes/com/jogamp/opengl/util/stereo/generic/GenericStereoDeviceFactory.java b/src/jogl/classes/com/jogamp/opengl/util/stereo/generic/GenericStereoDeviceFactory.java
index 957758e78..0cdef8770 100644
--- a/src/jogl/classes/com/jogamp/opengl/util/stereo/generic/GenericStereoDeviceFactory.java
+++ b/src/jogl/classes/com/jogamp/opengl/util/stereo/generic/GenericStereoDeviceFactory.java
@@ -35,6 +35,7 @@ import com.jogamp.nativewindow.util.Dimension;
 import com.jogamp.nativewindow.util.DimensionImmutable;
 import com.jogamp.opengl.math.FloatUtil;
 import com.jogamp.opengl.math.FovHVHalves;
+import com.jogamp.opengl.math.Vec3f;
 import com.jogamp.opengl.util.stereo.EyeParameter;
 import com.jogamp.opengl.util.stereo.StereoDevice;
 import com.jogamp.opengl.util.stereo.StereoDeviceConfig;
@@ -53,7 +54,7 @@ public class GenericStereoDeviceFactory extends StereoDeviceFactory {
      */
     public static GenericStereoDeviceConfig createMono(final String name,
                              final DimensionImmutable surfaceSizeInPixel, final float[] screenSizeInMeters,
-                             final float[] defaultEyePositionOffset) {
+                             final Vec3f defaultEyePositionOffset) {
         final float pupilCenterFromScreenTopInMeters = screenSizeInMeters[1] / 2f;
         final float d2r = FloatUtil.PI / 180.0f;
         return new GenericStereoDeviceConfig(
@@ -90,7 +91,7 @@ public class GenericStereoDeviceFactory extends StereoDeviceFactory {
     public static GenericStereoDeviceConfig createStereoSBS(final String name,
                              final DimensionImmutable surfaceSizeInPixel, final float[] screenSizeInMeters,
                              final float interpupillaryDistanceInMeters, final float fovy,
-                             final float[] defaultEyePositionOffset) {
+                             final Vec3f defaultEyePositionOffset) {
         final float pupilCenterFromScreenTopInMeters = screenSizeInMeters[1] / 2f;
         final float d2r = FloatUtil.PI / 180.0f;
 
@@ -138,7 +139,7 @@ public class GenericStereoDeviceFactory extends StereoDeviceFactory {
                              final DimensionImmutable surfaceSizeInPixel, final float[] screenSizeInMeters,
                              final float interpupillaryDistanceInMeters, final float fovy,
                              final DimensionImmutable eyeTextureSize,
-                             final float[] defaultEyePositionOffset) {
+                             final Vec3f defaultEyePositionOffset) {
         DistortionMesh.Producer lenseDistMeshProduce = null;
         try {
             lenseDistMeshProduce =