Flywheel/joml/sampling/BestCandidateSampling.java

/*
 * The MIT License
 *
 * Copyright (c) 2016-2021 JOML
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
package com.jozufozu.flywheel.repack.joml.sampling;

import java.nio.FloatBuffer;
import java.util.ArrayList;

import com.jozufozu.flywheel.repack.joml.Random;
import com.jozufozu.flywheel.repack.joml.Vector2f;
import com.jozufozu.flywheel.repack.joml.Vector3f;

/**
 * Creates samples using the "Best Candidate" algorithm.
 *
 * @author Kai Burjack
 */
public class BestCandidateSampling {

    private static final class IntHolder {
        int value;
    }

    /**
     * Generates Best Candidate samples on a unit sphere.
     * <p>
     * References:
     * <ul>
     * <li><a href="https://arxiv.org/ftp/cs/papers/0701/0701164.pdf">Indexing the Sphere with the Hierarchical Triangular Mesh</a>
     * <li><a href="http://math.stackexchange.com/questions/1244512/point-in-a-spherical-triangle-test">Point in a spherical triangle test</a>
     * </ul>
     *
     * @author Kai Burjack
     */
    public static class Sphere {
        /**
         * Implementation of a Hierarchical Triangular Mesh structure to index the sample points on the unit sphere for accelerating 1-nearest neighbor searches.
         *
         * @author Kai Burjack
         */
        private static final class Node {
            private static final int MAX_OBJECTS_PER_NODE = 32;

            private float v0x, v0y, v0z;
            private float v1x, v1y, v1z;
            private float v2x, v2y, v2z;
            private float cx, cy, cz;
            private float arc;

            private ArrayList objects;
            private Node[] children;

            Node() {
                this.children = new Node[8];
                float s = 1f;
                this.arc = 2.0f * (float) Math.PI;
                /*
                 * See: https://arxiv.org/ftp/cs/papers/0701/0701164.pdf
                 */
                this.children[0] = new Node(-s, 0, 0, 0, 0, s, 0, s, 0);
                this.children[1] = new Node(0, 0, s, s, 0, 0, 0, s, 0);
                this.children[2] = new Node(s, 0, 0, 0, 0, -s, 0, s, 0);
                this.children[3] = new Node(0, 0, -s, -s, 0, 0, 0, s, 0);
                this.children[4] = new Node(-s, 0, 0, 0, -s, 0, 0, 0, s);
                this.children[5] = new Node(0, 0, s, 0, -s, 0, s, 0, 0);
                this.children[6] = new Node(s, 0, 0, 0, -s, 0, 0, 0, -s);
                this.children[7] = new Node(0, 0, -s, 0, -s, 0, -s, 0, 0);
            }

            private Node(float x0, float y0, float z0, float x1, float y1, float z1, float x2, float y2, float z2) {
                this.v0x = x0;
                this.v0y = y0;
                this.v0z = z0;
                this.v1x = x1;
                this.v1y = y1;
                this.v1z = z1;
                this.v2x = x2;
                this.v2y = y2;
                this.v2z = z2;
                cx = (v0x + v1x + v2x) / 3.0f;
                cy = (v0y + v1y + v2y) / 3.0f;
                cz = (v0z + v1z + v2z) / 3.0f;
                float invCLen = Math.invsqrt(cx * cx + cy * cy + cz * cz);
                cx *= invCLen;
                cy *= invCLen;
                cz *= invCLen;
                // Compute maximum radius around triangle centroid
                float arc1 = greatCircleDist(cx, cy, cz, v0x, v0y, v0z);
                float arc2 = greatCircleDist(cx, cy, cz, v1x, v1y, v1z);
                float arc3 = greatCircleDist(cx, cy, cz, v2x, v2y, v2z);
                float dist = Math.max(Math.max(arc1, arc2), arc3);
                /*
                 * Convert radius into diameter, but also take into account the linear
                 * arccos approximation we use.
                 * This value was obtained empirically!
                 */
                dist *= 1.7f;
                this.arc = dist;
            }

            private void split() {
                float w0x = v1x + v2x;
                float w0y = v1y + v2y;
                float w0z = v1z + v2z;
                float len0 = Math.invsqrt(w0x * w0x + w0y * w0y + w0z * w0z);
                w0x *= len0;
                w0y *= len0;
                w0z *= len0;
                float w1x = v0x + v2x;
                float w1y = v0y + v2y;
                float w1z = v0z + v2z;
                float len1 = Math.invsqrt(w1x * w1x + w1y * w1y + w1z * w1z);
                w1x *= len1;
                w1y *= len1;
                w1z *= len1;
                float w2x = v0x + v1x;
                float w2y = v0y + v1y;
                float w2z = v0z + v1z;
                float len2 = Math.invsqrt(w2x * w2x + w2y * w2y + w2z * w2z);
                w2x *= len2;
                w2y *= len2;
                w2z *= len2;
                children = new Node[4];
                /*
                 * See: https://arxiv.org/ftp/cs/papers/0701/0701164.pdf
                 */
                children[0] = new Node(v0x, v0y, v0z, w2x, w2y, w2z, w1x, w1y, w1z);
                children[1] = new Node(v1x, v1y, v1z, w0x, w0y, w0z, w2x, w2y, w2z);
                children[2] = new Node(v2x, v2y, v2z, w1x, w1y, w1z, w0x, w0y, w0z);
                children[3] = new Node(w0x, w0y, w0z, w1x, w1y, w1z, w2x, w2y, w2z);
            }

            private void insertIntoChild(Vector3f o) {
                for (int i = 0; i < children.length; i++) {
                    Node c = children[i];
                    /*
                     * Idea: Test whether ray from origin towards point cuts triangle
                     *
                     * See: http://math.stackexchange.com/questions/1244512/point-in-a-spherical-triangle-test
                     */
                    if (isPointOnSphericalTriangle(o.x, o.y, o.z, c.v0x, c.v0y, c.v0z, c.v1x, c.v1y, c.v1z, c.v2x, c.v2y, c.v2z, 1E-6f)) {
                        c.insert(o);
                        return;
                    }
                }
            }

            void insert(Vector3f object) {
                if (children != null) {
                    insertIntoChild(object);
                    return;
                }
                if (objects != null && objects.size() == MAX_OBJECTS_PER_NODE) {
                    split();
                    for (int i = 0; i < MAX_OBJECTS_PER_NODE; i++)
                        insertIntoChild((Vector3f) objects.get(i));
                    objects = null;
                    insertIntoChild(object);
                } else {
                    if (objects == null)
                        objects = new ArrayList(MAX_OBJECTS_PER_NODE);
                    objects.add(object);
                }
            }

            /**
             * This is essentially a ray cast from the origin of the sphere to the point to test and then checking whether that ray goes through the triangle.
             * <p>
             * Reference: <a href="http://www.graphics.cornell.edu/pubs/1997/MT97.pdf">Fast,
             * Minimum Storage Ray/Triangle Intersection</a>
             */
            private static boolean isPointOnSphericalTriangle(float x, float y, float z, float v0X, float v0Y, float v0Z, float v1X, float v1Y, float v1Z, float v2X, float v2Y,
                    float v2Z, float epsilon) {
                float edge1X = v1X - v0X;
                float edge1Y = v1Y - v0Y;
                float edge1Z = v1Z - v0Z;
                float edge2X = v2X - v0X;
                float edge2Y = v2Y - v0Y;
                float edge2Z = v2Z - v0Z;
                float pvecX = y * edge2Z - z * edge2Y;
                float pvecY = z * edge2X - x * edge2Z;
                float pvecZ = x * edge2Y - y * edge2X;
                float det = edge1X * pvecX + edge1Y * pvecY + edge1Z * pvecZ;
                if (det > -epsilon && det < epsilon)
                    return false;
                float tvecX = -v0X;
                float tvecY = -v0Y;
                float tvecZ = -v0Z;
                float invDet = 1.0f / det;
                float u = (tvecX * pvecX + tvecY * pvecY + tvecZ * pvecZ) * invDet;
                if (u < 0.0f || u > 1.0f)
                    return false;
                float qvecX = tvecY * edge1Z - tvecZ * edge1Y;
                float qvecY = tvecZ * edge1X - tvecX * edge1Z;
                float qvecZ = tvecX * edge1Y - tvecY * edge1X;
                float v = (x * qvecX + y * qvecY + z * qvecZ) * invDet;
                if (v < 0.0f || u + v > 1.0f)
                    return false;
                float t = (edge2X * qvecX + edge2Y * qvecY + edge2Z * qvecZ) * invDet;
                return t >= epsilon;
            }

            private int child(float x, float y, float z) {
                for (int i = 0; i < children.length; i++) {
                    Node c = children[i];
                    if (isPointOnSphericalTriangle(x, y, z, c.v0x, c.v0y, c.v0z, c.v1x, c.v1y, c.v1z, c.v2x, c.v2y, c.v2z, 1E-5f))
                        return i;
                }
                // No child found. This can happen in 'nearest()' when querying possible nearby nodes
                return 0;
            }

            /**
             * Reference: <a href="https://en.wikipedia.org/wiki/Great-circle_distance#Vector_version">https://en.wikipedia.org/</a>
             */
            private float greatCircleDist(float x1, float y1, float z1, float x2, float y2, float z2) {
                float dot = x1 * x2 + y1 * y2 + z1 * z2;
                /*
                 * Just use a linear function, because we (mostly) do less-than comparisons on the result.
                 * We just need a linear function which:
                 * f(-1) = PI
                 * f(0)  = PI/2
                 * f(1)  = 0
                 */
                return (float) (-Math.PIHalf * dot + Math.PIHalf);
                //return (float) Math.acos(dot);
            }

            float nearest(float x, float y, float z) {
                return nearest(x, y, z, Float.POSITIVE_INFINITY);
            }
            float nearest(float x, float y, float z, float n) {
                float gcd = greatCircleDist(x, y, z, cx, cy, cz);
                /*
                 * If great-circle-distance between query point and centroid is larger than the current smallest distance 'n' plus the great circle diameter 'arc', we abort here,
                 * because then it is not possible for any point in the triangle patch to be closer to the query point than 'n'.
                 */
                /*
                 * Yes, we are subtracting two great-circle distances from one another here, which we did not even compute correctly
                 * using our overly linear arccos approximation. But the 1.7 factor above will take care that we still stay conservative
                 * enough here and not rejecting triangle patches which would contain samples nearer than 'n'.
                 */
                if (gcd - arc > n)
                    return n;
                float nr = n;
                if (children != null) {
                    int num = children.length, mod = num-1;
                    for (int i = child(x, y, z), c = 0; c < num; i = (i + 1) & mod, c++) {
                        float n1 = children[i].nearest(x, y, z, nr);
                        nr = Math.min(n1, nr);
                    }
                    return nr;
                }
                for (int i = 0; objects != null && i < objects.size(); i++) {
                    Vector3f o = (Vector3f) objects.get(i);
                    float d = greatCircleDist(o.x, o.y, o.z, x, y, z);
                    if (d < nr)
                        nr = d;
                }
                return nr;
            }
        }

        private boolean onHemisphere;
        private int numSamples;
        private int numCandidates = 60; // <- use a reasonable default
        private long seed;

        /**
         * Create a new instance of {@link Sphere} to configure and generate 'best candidate' sample positions on the unit sphere.
         */
        public Sphere() {}

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xyzs</code> float array.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xyzs
         *            will hold the x, y and z coordinates of all samples in the order <code>XYZXYZXYZ...</code>.
         *            This array must have a length of at least <code>numSamples</code>
         * @return this
         */
        public Sphere generate(final float[] xyzs) {
            final IntHolder i = new IntHolder();
            return generate(new Callback3d() {
                public void onNewSample(float x, float y, float z) {
                    xyzs[3 * i.value + 0] = x;
                    xyzs[3 * i.value + 1] = y;
                    xyzs[3 * i.value + 2] = z;
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xyzs</code> FloatBuffer.
         * <p>
         * The samples will be written starting at the current position of the FloatBuffer. The position of the FloatBuffer will not be modified.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xyzs
         *            will hold the x, y and z coordinates of all samples in the order <code>XYZXYZXYZ...</code>.
         *            This FloatBuffer must have at least <code>numSamples</code> remaining elements.
         *            The position of the buffer will not be modified by this method
         * @return this
         */
        public Sphere generate(final FloatBuffer xyzs) {
            final IntHolder i = new IntHolder();
            final int pos = xyzs.position();
            return generate(new Callback3d() {
                public void onNewSample(float x, float y, float z) {
                    xyzs.put(pos + 3 * i.value + 0, x);
                    xyzs.put(pos + 3 * i.value + 1, y);
                    xyzs.put(pos + 3 * i.value + 2, z);
                    i.value++;
                }
            });
        }

        /**
         * Set the seed to initialize the pseudo-random number generator with.
         *
         * @param seed
         *          the seed value
         * @return this
         */
        public Sphere seed(long seed) {
            this.seed = seed;
            return this;
        }

        /**
         * Set the number of samples to generate.
         *
         * @param numSamples
         *          the number of samples
         * @return this
         */
        public Sphere numSamples(int numSamples) {
            this.numSamples = numSamples;
            return this;
        }

        /**
         * Set the number of candidates to try for each generated sample.
         *
         * @param numCandidates
         *          the number of candidates to try
         * @return this
         */
        public Sphere numCandidates(int numCandidates) {
            this.numCandidates = numCandidates;
            return this;
        }

        /**
         * Set whether to generate samples on a hemisphere around the <code>+Z</code> axis.
         * <p>
         * The default is <code>false</code>, which will generate samples on the whole unit sphere.
         *
         * @param onHemisphere
         *          whether to generate samples on the hemisphere
         * @return this
         */
        public Sphere onHemisphere(boolean onHemisphere) {
            this.onHemisphere = onHemisphere;
            return this;
        }

        /**
         * Generate 'best candidate' sample call the given <code>callback</code> for each generated sample.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param callback
         *            will be called with the coordinates of each generated sample position
         * @return this
         */
        public Sphere generate(Callback3d callback) {
            Random rnd = new Random(seed);
            Node otree = new Node();
            for (int i = 0; i < numSamples; i++) {
                float bestX = Float.NaN, bestY = Float.NaN, bestZ = Float.NaN, bestDist = 0.0f;
                for (int c = 0; c < numCandidates; c++) {
                    /*
                     * Random point on sphere
                     *
                     * Reference: <a href="http://mathworld.wolfram.com/SpherePointPicking.html">http://mathworld.wolfram.com/</a>
                     */
                    float x1, x2;
                    do {
                        x1 = rnd.nextFloat() * 2.0f - 1.0f;
                        x2 = rnd.nextFloat() * 2.0f - 1.0f;
                    } while (x1 * x1 + x2 * x2 > 1.0f);
                    float sqrt = (float) Math.sqrt(1.0 - x1 * x1 - x2 * x2);
                    float x = 2 * x1 * sqrt;
                    float y = 2 * x2 * sqrt;
                    float z = 1.0f - 2.0f * (x1 * x1 + x2 * x2);
                    if (onHemisphere) {
                        z = Math.abs(z);
                    }
                    float minDist = otree.nearest(x, y, z);
                    if (minDist > bestDist) {
                        bestDist = minDist;
                        bestX = x;
                        bestY = y;
                        bestZ = z;
                    }
                }
                callback.onNewSample(bestX, bestY, bestZ);
                otree.insert(new Vector3f(bestX, bestY, bestZ));
            }
            return this;
        }
    }

    /**
     * Simple quatree that can handle points and 1-nearest neighbor search.
     *
     * @author Kai Burjack
     */
    private static class QuadTree {
        private static final int MAX_OBJECTS_PER_NODE = 32;

        // Constants for the quadrants of the quadtree
        private static final int PXNY = 0;
        private static final int NXNY = 1;
        private static final int NXPY = 2;
        private static final int PXPY = 3;

        private float minX, minY, hs;
        private ArrayList objects;
        private QuadTree[] children;

        QuadTree(float minX, float minY, float size) {
            this.minX = minX;
            this.minY = minY;
            this.hs = size * 0.5f;
        }

        private void split() {
            children = new QuadTree[4];
            children[NXNY] = new QuadTree(minX, minY, hs);
            children[PXNY] = new QuadTree(minX + hs, minY, hs);
            children[NXPY] = new QuadTree(minX, minY + hs, hs);
            children[PXPY] = new QuadTree(minX + hs, minY + hs, hs);
        }

        private void insertIntoChild(Vector2f o) {
            children[quadrant(o.x, o.y)].insert(o);
        }

        void insert(Vector2f object) {
            if (children != null) {
                insertIntoChild(object);
                return;
            }
            if (objects != null && objects.size() == MAX_OBJECTS_PER_NODE) {
                split();
                for (int i = 0; i < objects.size(); i++)
                    insertIntoChild((Vector2f) objects.get(i));
                objects = null;
                insertIntoChild(object);
            } else {
                if (objects == null)
                    objects = new ArrayList(MAX_OBJECTS_PER_NODE);
                objects.add(object);
            }
        }

        private int quadrant(float x, float y) {
            if (x < minX + hs) {
                if (y < minY + hs)
                    return NXNY;
                return NXPY;
            }
            if (y < minY + hs)
                return PXNY;
            return PXPY;
        }

        float nearest(float x, float y, float lowerBound, float upperBound) {
            float ub = upperBound;
            if (x < minX - upperBound || x > minX + hs * 2 + upperBound || y < minY - upperBound
                    || y > minY + hs * 2 + upperBound)
                return ub;
            if (children != null) {
                for (int i = quadrant(x, y), c = 0; c < 4; i = (i + 1) & 3, c++) {
                    float n1 = children[i].nearest(x, y, lowerBound, ub);
                    ub = Math.min(n1, ub);
                    if (ub <= lowerBound)
                        return lowerBound;
                }
                return ub;
            }
            float ub2 = ub * ub;
            float lb2 = lowerBound * lowerBound;
            for (int i = 0; objects != null && i < objects.size(); i++) {
                Vector2f o = (Vector2f) objects.get(i);
                float d = o.distanceSquared(x, y);
                if (d <= lb2)
                    return lowerBound;
                if (d < ub2)
                    ub2 = d;
            }
            return (float) Math.sqrt(ub2);
        }
    }

    /**
     * Generates Best Candidate samples on a unit disk.
     *
     * @author Kai Burjack
     */
    public static class Disk {
        private int numSamples;
        private int numCandidates = 60; // <- use a reasonable default
        private long seed;

        /**
         * Create a new instance of {@link Disk} to configure and generate 'best candidate' sample positions on the unit disk.
         */
        public Disk() {}

        /**
         * Set the seed to initialize the pseudo-random number generator with.
         *
         * @param seed
         *          the seed value
         * @return this
         */
        public Disk seed(long seed) {
            this.seed = seed;
            return this;
        }

        /**
         * Set the number of samples to generate.
         *
         * @param numSamples
         *          the number of samples
         * @return this
         */
        public Disk numSamples(int numSamples) {
            this.numSamples = numSamples;
            return this;
        }

        /**
         * Set the number of candidates to try for each generated sample.
         *
         * @param numCandidates
         *          the number of candidates to try
         * @return this
         */
        public Disk numCandidates(int numCandidates) {
            this.numCandidates = numCandidates;
            return this;
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xys</code> float array.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xys
         *            will hold the x and y coordinates of all samples in the order <code>XYXYXY...</code>.
         *            This array must have a length of at least <code>numSamples</code>
         * @return this
         */
        public Disk generate(final float[] xys) {
            final IntHolder i = new IntHolder();
            return generate(new Callback2d() {
                public void onNewSample(float x, float y) {
                    xys[2 * i.value + 0] = x;
                    xys[2 * i.value + 1] = y;
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xys</code> FloatBuffer.
         * <p>
         * The samples will be written starting at the current position of the FloatBuffer. The position of the FloatBuffer will not be modified.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xys
         *            will hold the x and y coordinates of all samples in the order <code>XYXYXY...</code>. This FloatBuffer must have at least <code>numSamples</code> remaining elements. The
         *            position of the buffer will not be modified by this method
         * @return this
         */
        public Disk generate(final FloatBuffer xys) {
            final IntHolder i = new IntHolder();
            final int pos = xys.position();
            return generate(new Callback2d() {
                public void onNewSample(float x, float y) {
                    xys.put(pos + 3 * i.value + 0, x);
                    xys.put(pos + 3 * i.value + 1, y);
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and call the given <code>callback</code> for each generated sample.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param callback
         *            will be called with the coordinates of each generated sample position
         * @return this
         */
        public Disk generate(Callback2d callback) {
            QuadTree qtree = new QuadTree(-1, -1, 2);
            Random rnd = new Random(seed);
            for (int i = 0; i < numSamples; i++) {
                float bestX = 0, bestY = 0, bestDist = 0.0f;
                for (int c = 0; c < numCandidates; c++) {
                    float x, y;
                    do {
                        x = rnd.nextFloat() * 2.0f - 1.0f;
                        y = rnd.nextFloat() * 2.0f - 1.0f;
                    } while (x * x + y * y > 1.0f);
                    float minDist = qtree.nearest(x, y, bestDist, Float.POSITIVE_INFINITY);
                    if (minDist > bestDist) {
                        bestDist = minDist;
                        bestX = x;
                        bestY = y;
                    }
                }
                callback.onNewSample(bestX, bestY);
                qtree.insert(new Vector2f(bestX, bestY));
            }
            return this;
        }
    }

    /**
     * Generates Best Candidate samples on a unit quad.
     *
     * @author Kai Burjack
     */
    public static class Quad {
        private int numSamples;
        private int numCandidates = 60; // <- use a reasonable default
        private long seed;

        /**
         * Create a new instance of {@link Quad} to configure and generate 'best candidate' sample positions on the unit quad.
         */
        public Quad() {}

        /**
         * Set the seed to initialize the pseudo-random number generator with.
         *
         * @param seed
         *          the seed value
         * @return this
         */
        public Quad seed(long seed) {
            this.seed = seed;
            return this;
        }

        /**
         * Set the number of samples to generate.
         *
         * @param numSamples
         *          the number of samples
         * @return this
         */
        public Quad numSamples(int numSamples) {
            this.numSamples = numSamples;
            return this;
        }

        /**
         * Set the number of candidates to try for each generated sample.
         *
         * @param numCandidates
         *          the number of candidates to try
         * @return this
         */
        public Quad numCandidates(int numCandidates) {
            this.numCandidates = numCandidates;
            return this;
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xyzs</code> float array.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xyzs
         *            will hold the x, y and z coordinates of all samples in the order <code>XYZXYZXYZ...</code>.
         *            This array must have a length of at least <code>numSamples</code>
         * @return this
         */
        public Quad generate(final float[] xyzs) {
            final IntHolder i = new IntHolder();
            return generate(new Callback2d() {
                public void onNewSample(float x, float y) {
                    xyzs[2 * i.value + 0] = x;
                    xyzs[2 * i.value + 1] = y;
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xys</code> FloatBuffer.
         * <p>
         * The samples will be written starting at the current position of the FloatBuffer. The position of the FloatBuffer will not be modified.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xys
         *            will hold the x and y coordinates of all samples in the order <code>XYXYXY...</code>. This FloatBuffer must have at least <code>numSamples</code> remaining elements. The position of
         *            the buffer will not be modified by this method
         * @return this
         */
        public Quad generate(final FloatBuffer xys) {
            final IntHolder i = new IntHolder();
            final int pos = xys.position();
            return generate(new Callback2d() {
                public void onNewSample(float x, float y) {
                    xys.put(pos + 3 * i.value + 0, x);
                    xys.put(pos + 3 * i.value + 1, y);
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and call the given <code>callback</code> for each generated sample.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param callback
         *            will be called with the coordinates of each generated sample position
         * @return this
         */
        public Quad generate(Callback2d callback) {
            QuadTree qtree = new QuadTree(-1, -1, 2);
            Random rnd = new Random(seed);
            for (int i = 0; i < numSamples; i++) {
                float bestX = 0, bestY = 0, bestDist = 0.0f;
                for (int c = 0; c < numCandidates; c++) {
                    float x = rnd.nextFloat() * 2.0f - 1.0f;
                    float y = rnd.nextFloat() * 2.0f - 1.0f;
                    float minDist = qtree.nearest(x, y, bestDist, Float.POSITIVE_INFINITY);
                    if (minDist > bestDist) {
                        bestDist = minDist;
                        bestX = x;
                        bestY = y;
                    }
                }
                callback.onNewSample(bestX, bestY);
                qtree.insert(new Vector2f(bestX, bestY));
            }
            return this;
        }
    }

    /**
     * Simple octree for points and 1-nearest neighbor distance query.
     *
     * @author Kai Burjack
     */
    private static class Octree {
        private static final int MAX_OBJECTS_PER_NODE = 32;

        // Constants for the octants of the octree
        private static final int PXNYNZ = 0;
        private static final int NXNYNZ = 1;
        private static final int NXPYNZ = 2;
        private static final int PXPYNZ = 3;
        private static final int PXNYPZ = 4;
        private static final int NXNYPZ = 5;
        private static final int NXPYPZ = 6;
        private static final int PXPYPZ = 7;

        private float minX, minY, minZ, hs;
        private ArrayList objects;
        private Octree[] children;

        Octree(float minX, float minY, float minZ, float size) {
            this.minX = minX;
            this.minY = minY;
            this.minZ = minZ;
            this.hs = size * 0.5f;
        }

        private void split() {
            children = new Octree[8];
            children[NXNYNZ] = new Octree(minX, minY, minZ, hs);
            children[PXNYNZ] = new Octree(minX + hs, minY, minZ, hs);
            children[NXPYNZ] = new Octree(minX, minY + hs, minZ, hs);
            children[PXPYNZ] = new Octree(minX + hs, minY + hs, minZ, hs);
            children[NXNYPZ] = new Octree(minX, minY, minZ + hs, hs);
            children[PXNYPZ] = new Octree(minX + hs, minY, minZ + hs, hs);
            children[NXPYPZ] = new Octree(minX, minY + hs, minZ + hs, hs);
            children[PXPYPZ] = new Octree(minX + hs, minY + hs, minZ + hs, hs);
        }

        private void insertIntoChild(Vector3f o) {
            children[octant(o.x, o.y, o.z)].insert(o);
        }

        void insert(Vector3f object) {
            if (children != null) {
                insertIntoChild(object);
                return;
            }
            if (objects != null && objects.size() == MAX_OBJECTS_PER_NODE) {
                split();
                for (int i = 0; i < objects.size(); i++)
                    insertIntoChild((Vector3f) objects.get(i));
                objects = null;
                insertIntoChild(object);
            } else {
                if (objects == null)
                    objects = new ArrayList(MAX_OBJECTS_PER_NODE);
                objects.add(object);
            }
        }

        private int octant(float x, float y, float z) {
            if (x < minX + hs)
                if (y < minY + hs) {
                    if (z < minZ + hs)
                        return NXNYNZ;
                    return NXNYPZ;
                } else if (z < minZ + hs)
                    return NXPYNZ;
                else
                    return NXPYPZ;
            else if (y < minY + hs) {
                if (z < minZ + hs)
                    return PXNYNZ;
                return PXNYPZ;
            } else if (z < minZ + hs)
                return PXPYNZ;
            else
                return PXPYPZ;
        }

        float nearest(float x, float y, float z, float lowerBound, float upperBound) {
            float up = upperBound;
            if (x < minX - upperBound || x > minX + hs * 2 + upperBound || y < minY - upperBound || y > minY + hs * 2 + upperBound ||
                    z < minZ - upperBound || z > minZ + hs * 2 + upperBound)
                return up;
            if (children != null) {
                for (int i = octant(x, y, z), c = 0; c < 8; i = (i + 1) & 7, c++) {
                    float n1 = children[i].nearest(x, y, z, lowerBound, up);
                    up = Math.min(n1, up);
                    if (up <= lowerBound)
                        return lowerBound;
                }
                return up;
            }
            float up2 = up * up;
            float lb2 = lowerBound * lowerBound;
            for (int i = 0; objects != null && i < objects.size(); i++) {
                Vector3f o = (Vector3f) objects.get(i);
                float d = o.distanceSquared(x, y, z);
                if (d <= lb2)
                    return lowerBound;
                if (d < up2)
                    up2 = d;
            }
            return (float) Math.sqrt(up2);
        }
    }

    /**
     * Generates Best Candidate samples inside a unit cube.
     *
     * @author Kai Burjack
     */
    public static class Cube {
        private int numSamples;
        private int numCandidates = 60; // <- use a reasonable default
        private long seed;

        /**
         * Create a new instance of {@link Cube} to configure and generate 'best candidate' sample positions
         * on the unit cube with each sample tried <code>numCandidates</code> number of times, and call the
         * given <code>callback</code> for each sample generate.
         */
        public Cube() {}

        /**
         * Set the seed to initialize the pseudo-random number generator with.
         *
         * @param seed
         *          the seed value
         * @return this
         */
        public Cube seed(long seed) {
            this.seed = seed;
            return this;
        }

        /**
         * Set the number of samples to generate.
         *
         * @param numSamples
         *          the number of samples
         * @return this
         */
        public Cube numSamples(int numSamples) {
            this.numSamples = numSamples;
            return this;
        }

        /**
         * Set the number of candidates to try for each generated sample.
         *
         * @param numCandidates
         *          the number of candidates to try
         * @return this
         */
        public Cube numCandidates(int numCandidates) {
            this.numCandidates = numCandidates;
            return this;
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xyzs</code> float array.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xyzs
         *            will hold the x, y and z coordinates of all samples in the order <code>XYZXYZXYZ...</code>.
         *            This array must have a length of at least <code>numSamples</code>
         * @return this
         */
        public Cube generate(final float[] xyzs) {
            final IntHolder i = new IntHolder();
            return generate(new Callback3d() {
                public void onNewSample(float x, float y, float z) {
                    xyzs[3 * i.value + 0] = x;
                    xyzs[3 * i.value + 1] = y;
                    xyzs[3 * i.value + 2] = z;
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and store the coordinates of all generated samples into the given <code>xyzs</code> FloatBuffer.
         * <p>
         * The samples will be written starting at the current position of the FloatBuffer. The position of the FloatBuffer will not be modified.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param xyzs
         *            will hold the x, y and z coordinates of all samples in the order <code>XYZXYZXYZ...</code>.
         *            This FloatBuffer must have at least <code>numSamples</code> remaining elements.
         *            The position of the buffer will not be modified by this method
         * @return this
         */
        public Cube generate(final FloatBuffer xyzs) {
            final IntHolder i = new IntHolder();
            final int pos = xyzs.position();
            return generate(new Callback3d() {
                public void onNewSample(float x, float y, float z) {
                    xyzs.put(pos + 3 * i.value + 0, x);
                    xyzs.put(pos + 3 * i.value + 1, y);
                    xyzs.put(pos + 3 * i.value + 2, z);
                    i.value++;
                }
            });
        }

        /**
         * Generate 'best candidate' sample positions and call the given <code>callback</code> for each generated sample.
         * <p>
         * <em>This method performs heap allocations, so should be used sparingly.</em>
         *
         * @param callback
         *            will be called with the coordinates of each generated sample position
         * @return this
         */
        public Cube generate(Callback3d callback) {
            Octree octree = new Octree(-1, -1, -1, 2);
            Random rnd = new Random(seed);
            for (int i = 0; i < numSamples; i++) {
                float bestX = 0, bestY = 0, bestZ = 0, bestDist = 0.0f;
                for (int c = 0; c < numCandidates; c++) {
                    float x = rnd.nextFloat() * 2.0f - 1.0f;
                    float y = rnd.nextFloat() * 2.0f - 1.0f;
                    float z = rnd.nextFloat() * 2.0f - 1.0f;
                    float minDist = octree.nearest(x, y, z, bestDist, Float.POSITIVE_INFINITY);
                    if (minDist > bestDist) {
                        bestDist = minDist;
                        bestX = x;
                        bestY = y;
                        bestZ = z;
                    }
                }
                callback.onNewSample(bestX, bestY, bestZ);
                octree.insert(new Vector3f(bestX, bestY, bestZ));
            }
            return this;
        }
    }

}