/*
* (C) Copyright 2016-2017 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 org.joml.sampling;
import java.util.ArrayList;
import org.joml.Vector2f;
import org.joml.Vector3f;
/**
* Creates samples using the "Best Candidate" algorithm.
*
* @author Kai Burjack
*/
public class BestCandidateSampling {
/**
* Generates Best Candidate samples on a unit sphere.
* <p>
* References:
* <ul>
* <li><a href="https://www.microsoft.com/en-us/research/wp-content/uploads/2005/09/tr-2005-123.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://www.microsoft.com/en-us/research/wp-content/uploads/2005/09/tr-2005-123.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 = 1.0f / (float) Math.sqrt(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 = 1.0f / (float) Math.sqrt(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 = 1.0f / (float) Math.sqrt(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 = 1.0f / (float) Math.sqrt(w2x * w2x + w2y * w2y + w2z * w2z);
w2x *= len2;
w2y *= len2;
w2z *= len2;
children = new Node[4];
/*
* See: https://www.microsoft.com/en-us/research/wp-content/uploads/2005/09/tr-2005-123.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.cs.virginia.edu/~gfx/Courses/2003/ImageSynthesis/papers/Acceleration/Fast%20MinimumStorage%20RayTriangle%20Intersection.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, 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 final Random rnd;
private final Node otree;
/**
* Create a new instance of {@link Sphere}, initialize the random number generator with the given <code>seed</code> and generate <code>numSamples</code>
* number of 'best candidate' sample positions on the unit sphere with each sample tried <code>numCandidates</code> number of times, and call the given
* <code>callback</code> for each sample generate.
*
* @param seed
* the seed to initialize the random number generator with
* @param numSamples
* the number of samples to generate
* @param numCandidates
* the number of candidates to test for each sample
* @param callback
* will be called for each sample generated
*/
public Sphere(long seed, int numSamples, int numCandidates, Callback3d callback) {
this.rnd = new Random(seed);
this.otree = new Node();
compute(numSamples, numCandidates, callback);
}
private void compute(int numSamples, int numCandidates, Callback3d callback) {
for (int i = 0; i < numSamples; i++) {
float bestX = 0, bestY = 0, bestZ = 0, 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);
float minDist = otree.nearest(x, y, z, Float.POSITIVE_INFINITY);
if (minDist > bestDist) {
bestDist = minDist;
bestX = x;
bestY = y;
bestZ = z;
}
}
callback.onNewSample(bestX, bestY, bestZ);
otree.insert(new Vector3f(bestX, bestY, bestZ));
}
}
}
/**
* 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 final Random rnd;
private final QuadTree qtree;
/**
* Create a new instance of {@link Disk}, initialize the random number generator with the given <code>seed</code> and generate <code>numSamples</code>
* number of 'best candidate' sample positions on the unit disk with each sample tried <code>numCandidates</code> number of times, and call the given <code>callback</code>
* for each sample generate.
*
* @param seed
* the seed to initialize the random number generator with
* @param numSamples
* the number of samples to generate
* @param numCandidates
* the number of candidates to test for each sample
* @param callback
* will be called for each sample generated
*/
public Disk(long seed, int numSamples, int numCandidates, Callback2d callback) {
this.rnd = new Random(seed);
this.qtree = new QuadTree(-1, -1, 2);
generate(numSamples, numCandidates, callback);
}
private void generate(int numSamples, int numCandidates, Callback2d callback) {
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));
}
}
}
/**
* Generates Best Candidate samples on a unit quad.
*
* @author Kai Burjack
*/
public static class Quad {
private final Random rnd;
private final QuadTree qtree;
/**
* Create a new instance of {@link Quad}, initialize the random number generator with the given <code>seed</code> and generate <code>numSamples</code>
* number of 'best candidate' sample positions on the unit quad with each sample tried <code>numCandidates</code> number of times, and call the given <code>callback</code>
* for each sample generate.
*
* @param seed
* the seed to initialize the random number generator with
* @param numSamples
* the number of samples to generate
* @param numCandidates
* the number of candidates to test for each sample
* @param callback
* will be called for each sample generated
*/
public Quad(long seed, int numSamples, int numCandidates, Callback2d callback) {
this.rnd = new Random(seed);
this.qtree = new QuadTree(-1, -1, 2);
generate(numSamples, numCandidates, callback);
}
private void generate(int numSamples, int numCandidates, Callback2d callback) {
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));
}
}
}
/**
* 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 final Random rnd;
private final Octree octree;
/**
* Create a new instance of {@link Cube}, initialize the random number generator with the given <code>seed</code> and generate <code>numSamples</code>
* number of '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.
*
* @param seed
* the seed to initialize the random number generator with
* @param numSamples
* the number of samples to generate
* @param numCandidates
* the number of candidates to test for each sample
* @param callback
* will be called for each sample generated
*/
public Cube(long seed, int numSamples, int numCandidates, Callback3d callback) {
this.rnd = new Random(seed);
this.octree = new Octree(-1, -1, -1, 2);
generate(numSamples, numCandidates, callback);
}
private void generate(int numSamples, int numCandidates, Callback3d callback) {
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));
}
}
}
}