Add class for Linearly Transformed Cosines evaluation

and refactor NelderMead3D

src/org/lwjgl/demo/util/LTC.java

@@ -0,0 +1,84 @@
+package org.lwjgl.demo.util;
+
+import org.joml.Vector3f;
+
+import static java.lang.Math.*;
+
+/**
+ * Linearly Transformed Cosines evaluation.
+ * <p>
+ * This code was adapted from: <a href=
+ * "https://github.com/selfshadow/ltc_code/blob/master/fit/LTC.h">https://github.com/selfshadow/ltc_code/blob/master/fit/LTC.h</a>
+ * 
+ * @author Kai Burjack
+ */
+public class LTC {
+    private static final float TWO_PI = (float) (2.0 * Math.PI);
+    private static final float ONE_OVER_PI = (float) (1.0 / Math.PI);
+
+    public float magnitude = 1;
+    public float fresnel = 1;
+    public float m11 = 1, m22 = 1, m13;
+    public float Xx = 1, Xz, Zx, Zz = 1;
+
+    private float absInvDet;
+    private float nm00, nm02, nm20, nm22;
+    private float im00, im02, im11, im20, im22;
+    {
+        update();
+    }
+
+    public void update() {
+        nm00 = Xx * m11;
+        nm02 = Xz * m11;
+        nm20 = Xx * m13 + Zx;
+        nm22 = Xz * m13 + Zz;
+        float s = 1.0f / (nm00 * (m22 * nm22) + nm20 * (-m22 * nm02));
+        im00 = (m22 * nm22) * s;
+        im02 = -(m22 * nm02) * s;
+        im11 = (nm00 * nm22 - nm20 * nm02) * s;
+        im20 = -(nm20 * m22) * s;
+        im22 = (nm00 * m22) * s;
+        absInvDet = abs(s);
+    }
+
+    public float eval(float x, float y, float z) {
+        float rx = im00 * x + im20 * z;
+        float ry = im11 * y;
+        float rz = im02 * x + im22 * z;
+        float invlen = 1.0f / (float) sqrt(rx * rx + ry * ry + rz * rz);
+        rx *= invlen;
+        ry *= invlen;
+        rz *= invlen;
+        float kx = nm00 * rx + nm20 * rz;
+        float ky = m22 * ry;
+        float kz = nm02 * rx + nm22 * rz;
+        float len = (float) sqrt(kx * kx + ky * ky + kz * kz);
+        return ONE_OVER_PI * magnitude * max(0.0f, rz) * (len * len * len) * absInvDet;
+    }
+
+    public Vector3f sampleUniform(float cosTheta, float phi01, Vector3f dest) {
+        return sample(cosTheta, phi01, dest);
+    }
+
+    public Vector3f sampleCosineWeighted(float cosTheta, float phi01, Vector3f dest) {
+        return sample((float) sqrt(cosTheta), phi01, dest);
+    }
+
+    private Vector3f sample(float c, float phi, Vector3f dest) {
+        float s = (float) sqrt(1.0 - c * c);
+        float p = TWO_PI * phi;
+        float sinp = (float) sin(p);
+        float cosp = (float) cos(p);
+        float x = s * cosp, y = s * sinp, z = c;
+        float kx = nm00 * x + nm20 * z;
+        float ky = m22 * y;
+        float kz = nm02 * x + nm22 * z;
+        float invlen = 1.0f / (float) sqrt(kx * kx + ky * ky + kz * kz);
+        dest.x = kx * invlen;
+        dest.y = ky * invlen;
+        dest.z = kz * invlen;
+        return dest;
+    }
+
+}

src/org/lwjgl/demo/util/NelderMead3D.java

@@ -1,5 +1,7 @@
 package org.lwjgl.demo.util;
 
+import static java.lang.Math.*;
+
 /**
  * Nelder-Mead algorithm for non-linear optimization.
  * <p>
@@ -16,9 +18,23 @@ public interface FitnessFunction {
 
     private final float[] s = new float[4 * 3];
     private final float[] f = new float[4];
+    private final float reflect;
+    private final float expand;
+    private final float contract;
+    private final float shrink;
+
+    public NelderMead3D() {
+        this(1.0f, 2.0f, 0.5f, 0.5f);
+    }
+
+    public NelderMead3D(float reflect, float expand, float contract, float shrink) {
+        this.reflect = reflect;
+        this.expand = expand;
+        this.contract = contract;
+        this.shrink = shrink;
+    }
 
     public float minimize(float[] init, float d, float minDf, int N, FitnessFunction fn, float[] outMin) {
-        float reflect = 1.0f, expand = 2.0f, contract = 0.5f, shrink = 0.5f;
         s[0] = init[0];
         s[1] = init[1];
         s[2] = init[2];
@@ -33,107 +49,86 @@ public float minimize(float[] init, float d, float minDf, int N, FitnessFunction
         int lo = 0, hi, nh, i;
         for (i = 0; i < N; i++) {
             lo = hi = nh = 0;
-            if (f[1] < f[lo])
-                lo = 1;
-            if (f[1] > f[hi]) {
-                nh = hi;
-                hi = 1;
-            } else if (f[1] > f[nh])
-                nh = 1;
-            if (f[2] < f[lo])
-                lo = 2;
-            if (f[2] > f[hi]) {
-                nh = hi;
-                hi = 2;
-            } else if (f[2] > f[nh])
-                nh = 2;
-            if (f[3] < f[lo])
-                lo = 3;
-            if (f[3] > f[hi]) {
-                nh = hi;
-                hi = 3;
-            } else if (f[3] > f[nh])
-                nh = 3;
-            float a = Math.abs(f[lo]), b = Math.abs(f[hi]);
-            if (2.0f * Math.abs(a - b) <= (a + b) * minDf)
+            for (int j = 0; j < 3; j++) {
+                if (f[j + 1] < f[lo])
+                    lo = j + 1;
+                if (f[j + 1] > f[hi])
+                    hi = j + 1;
+                else if (f[j + 1] > f[nh])
+                    nh = j + 1;
+            }
+            float a = abs(f[lo]), b = abs(f[hi]);
+            if (2.0f * abs(a - b) <= (a + b) * minDf)
                 break;
             float ox = 0.0f, oy = 0.0f, oz = 0.0f;
-            if (0 != hi) {
-                ox += s[0];
-                oy += s[1];
-                oz += s[2];
-            }
-            if (1 != hi) {
-                ox += s[3];
-                oy += s[4];
-                oz += s[5];
-            }
-            if (2 != hi) {
-                ox += s[6];
-                oy += s[7];
-                oz += s[8];
-            }
-            if (3 != hi) {
-                ox += s[9];
-                oy += s[10];
-                oz += s[11];
+            for (int j = 0; j < 4; j++) {
+                if (j == hi)
+                    continue;
+                ox += s[3 * j];
+                oy += s[3 * j + 1];
+                oz += s[3 * j + 2];
             }
             ox *= 0.33333333333f;
             oy *= 0.33333333333f;
             oz *= 0.33333333333f;
-            float rx, ry, rz;
-            rx = ox + reflect * (ox - s[3 * hi]);
-            ry = oy + reflect * (oy - s[3 * hi + 1]);
-            rz = oz + reflect * (oz - s[3 * hi + 2]);
+            float rx = ox + reflect * (ox - s[3 * hi]);
+            float ry = oy + reflect * (oy - s[3 * hi + 1]);
+            float rz = oz + reflect * (oz - s[3 * hi + 2]);
             float fr = fn.f(rx, ry, rz);
             if (fr < f[nh]) {
-                minimize(fr, rx, ry, rz, ox, oy, oz, expand, lo, hi, fn);
-                continue;
+                if (!expand(fr, ox, oy, oz, lo, hi, fn))
+                    reflect(fr, rx, ry, rz, hi);
+            } else {
+                if (!contract(ox, oy, oz, hi, fn))
+                    shrink(lo, fn);
             }
-            minimize(ox, oy, oz, contract, lo, hi, shrink, fn);
         }
         outMin[0] = s[3 * lo];
         outMin[1] = s[3 * lo + 1];
         outMin[2] = s[3 * lo + 2];
         return f[lo];
     }
 
-    private void minimize(float fr, float rx, float ry, float rz, float ox, float oy, float oz, float expand, int lo,
-            int hi, FitnessFunction fn) {
+    private boolean expand(float fr, float ox, float oy, float oz, int lo, int hi, FitnessFunction fn) {
         if (fr < f[lo]) {
-            float ex, ey, ez;
-            ex = ox + expand * (ox - s[3 * hi]);
-            ey = oy + expand * (oy - s[3 * hi + 1]);
-            ez = oz + expand * (oz - s[3 * hi + 2]);
+            float ex = ox + expand * (ox - s[3 * hi]);
+            float ey = oy + expand * (oy - s[3 * hi + 1]);
+            float ez = oz + expand * (oz - s[3 * hi + 2]);
             float fe = fn.f(ex, ey, ez);
             if (fe < fr) {
                 s[3 * hi] = ex;
                 s[3 * hi + 1] = ey;
                 s[3 * hi + 2] = ez;
                 f[hi] = fe;
-                return;
+                return true;
             }
         }
+        return false;
+    }
+
+    private void reflect(float fr, float rx, float ry, float rz, int hi) {
         s[3 * hi] = rx;
         s[3 * hi + 1] = ry;
         s[3 * hi + 2] = rz;
         f[hi] = fr;
     }
 
-    private void minimize(float ox, float oy, float oz, float contract, int lo, int hi, float shrink,
-            FitnessFunction fn) {
-        float cx, cy, cz;
-        cx = ox - contract * (ox - s[3 * hi]);
-        cy = oy - contract * (oy - s[3 * hi + 1]);
-        cz = oz - contract * (oz - s[3 * hi + 2]);
+    private boolean contract(float ox, float oy, float oz, int hi, FitnessFunction fn) {
+        float cx = ox - contract * (ox - s[3 * hi]);
+        float cy = oy - contract * (oy - s[3 * hi + 1]);
+        float cz = oz - contract * (oz - s[3 * hi + 2]);
         float fc = fn.f(cx, cy, cz);
         if (fc < f[hi]) {
             s[3 * hi] = cx;
             s[3 * hi + 1] = cy;
             s[3 * hi + 2] = cz;
             f[hi] = fc;
-            return;
+            return true;
         }
+        return false;
+    }
+
+    private void shrink(int lo, FitnessFunction fn) {
         for (int k = 0; k < 4; k++) {
             if (k == lo)
                 continue;