/*
* BioJava development code
*
* This code may be freely distributed and modified under the
* terms of the GNU Lesser General Public Licence. This should
* be distributed with the code. If you do not have a copy,
* see:
*
* http://www.gnu.org/copyleft/lesser.html
*
* Copyright for this code is held jointly by the individual
* authors. These should be listed in @author doc comments.
*
* For more information on the BioJava project and its aims,
* or to join the biojava-l mailing list, visit the home page
* at:
*
* http://www.biojava.org/
*
*/
package org.biojava.nbio.structure.symmetry.geometry;
import javax.vecmath.*;
import java.util.ArrayList;
import java.util.List;
/**
* Generate the permutations and sign changes for a Triple.
* For instance, (a,0,0) becomes (±a,0,0),(0,±a,0),(0,0,±a). In the general
* case 48 tuples are returned.
*/
class Permute {
private List<Point3i> triples = new ArrayList<Point3i>();
Permute(Point3i t) {
Point3i tmp = new Point3i();
tmp.x = t.x;
tmp.y = t.y;
tmp.z = t.z;
triples.add(tmp);
int n = 1;
// Unpermuted ±x,±y,±z
if (t.x != 0) {
for (int i = 0; i < n; ++i) {
Tuple3i m = triples.get(i);
triples.add(new Point3i(-m.x, m.y, m.z));
}
n *= 2;
}
if (t.y != 0) {
for (int i = 0; i < n; ++i) {
Point3i m = triples.get(i);
triples.add(new Point3i(m.x, -m.y, m.z));
}
n *= 2;
}
if (t.z != 0) {
for (int i = 0; i < n; ++i) {
Point3i m = triples.get(i);
triples.add(new Point3i(m.x, m.y, -m.z));
}
n *= 2;
}
if (t.x == t.y && t.y == t.z) {
return;
}
// At least one differing value
for (int i = 0; i < n; ++i) {
Point3i m = triples.get(i);
triples.add(new Point3i(m.y, m.z, m.x));
triples.add(new Point3i(m.z, m.x, m.y));
}
n *= 3;
if (t.x == t.y || t.y == t.z) {
return;
}
// At least two differing values
for (int i = 0; i < n; ++i) {
Point3i m = triples.get(i);
triples.add(new Point3i(m.y, m.x, m.z));
}
n *= 2;
}
public int size() {
return triples.size();
}
public Point3i get(int i) {
return triples.get(i);
}
};
/**
* Sample possible orientations. An orientation can be represented as a
* quaternion or as a rotation axis and angle. The orientations are somewhat
* evenly distributed over orientation space, but the current implementation
* is not suited for statistically uniform sampling.
* @author Peter
*/
public final class SphereSampler {
private static final List<Quat4d> orientations ;
// The rotational symmetries of the cube. (Not normalized, since
// PackSet.Add does this.)
private static final double[][] cubeSyms = {
{ 1, 0, 0, 0 },
// 180 deg rotations about 3 axes
{ 0, 1, 0, 0 },
{ 0, 0, 1, 0 },
{ 0, 0, 0, 1 },
// +/- 120 degree rotations about 4 leading diagonals
{ 1, 1, 1, 1 }, { 1, 1, 1, -1 }, { 1, 1, -1, 1 }, { 1, 1, -1, -1 },
{ 1, -1, 1, 1 }, { 1, -1, 1, -1 },
{ 1, -1, -1, 1 },
{ 1, -1, -1, -1 },
// +/- 90 degree rotations about 3 axes
{ 1, 1, 0, 0 }, { 1, -1, 0, 0 }, { 1, 0, 1, 0 }, { 1, 0, -1, 0 },
{ 1, 0, 0, 1 }, { 1, 0, 0, -1 },
// 180 degree rotations about 6 face diagonals
{ 0, 1, 1, 0 }, { 0, 1, -1, 0 }, { 0, 1, 0, 1 }, { 0, 1, 0, -1 },
{ 0, 0, 1, 1 }, { 0, 0, 1, -1 }, };
// Approximate spacing between grid points
private static final double delta = 0.15846;
// Amount of distortion towards grid corners to account for the projection back to a 4-sphere
private static final double sigma = 0.00;
private static final int ntot = 7416; // total number of grid points over the 24 cells
private static final int ncell = 309; // number of grid points per cell
private static final int nent = 18; // number of distinct grid points (with k>=l>=m)
// Why not (5,5,3) and (5,5,5)? Would they overlap with other cells? -SB
private static final int[] k = { 0, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5};
private static final int[] l = { 0, 1, 0, 2, 2, 1, 3, 3, 0, 2, 2, 4, 4, 4, 1, 3, 3, 5};
private static final int[] m = { 0, 1, 0, 0, 2, 1, 1, 3, 0, 0, 2, 0, 2, 4, 1, 1, 3, 1};
// number of permutations to expect for each (k,l,m) tuple
private static final int[] mult = { 1, 8, 6, 12, 8, 24, 24, 8, 6, 24, 24, 12, 24, 8, 24, 48, 24, 24};
static
{
List<Quat4d> myorientations = new ArrayList<Quat4d>();
// 60 equally spaced grid points covering all quaternions
// This could be removed, since the 48-cell bcc lattice has lower angle error -SB
for (int i = 0; i < IcosahedralSampler.getSphereCount(); i++) {
myorientations.add(IcosahedralSampler.getQuat4d(i));
}
// SB's interpretation of this code:
// Directly form a 4D grid over the 4-sphere of orientations. Start by
// making a 3D body-centered lattice with w=1. Then use the hypercube symmetries
// to extend this grid over the whole surface.
// The permuted (k,l,m) values together make a diagonal grid out to ±(5,5,5).
// The point spacing is distorted by the pind() function so that the
// projection of the points back to the 4-sphere will be more even.
// This is the c48u309 lattice from Karney 2006, with a max error of 10.07
// degrees.
List<Quat4d> grid = new ArrayList<Quat4d>();
int ncell1 = 0;
for (int n = 0; n < nent; ++n) { // for each tuple (k,l,m) above
Permute p = new Permute(new Point3i(k[n], l[n], m[n]));
assert (mult[n] == p.size());
for (int i = 0; i < mult[n]; ++i) { // for each permutation
Point3i t = p.get(i);
grid.add(new Quat4d(1.0, pind(0.5 * t.x, delta, sigma), pind(
0.5 * t.y, delta, sigma), pind(0.5 * t.z, delta, sigma)));
}
ncell1 += mult[n];
}
assert (ncell1 == ncell);
// Apply symmetry operations to distribute grid over all values of w
int nc = grid.size();
assert (nc == ncell);
for (int n = 1; n < 24; ++n) {
Quat4d q = new Quat4d(cubeSyms[n][0], cubeSyms[n][1],
cubeSyms[n][2], cubeSyms[n][3]);
for (int i = 0; i < nc; ++i) {
Quat4d qs = new Quat4d();
qs.mul(q, grid.get(i));
grid.add(qs);
// s.add(times(q, s.getOrientation(i)), s.getWeight(i)); // this data set has no weights
// Actually, it does have weights but we don't care since we don't need uniform sampling. -SB
}
}
assert (grid.size() == ntot);
myorientations.addAll(grid);
orientations = myorientations;
}
// this class cannot be instantiated
private SphereSampler() {
};
public static int getSphereCount() {
return orientations.size();
}
public static Quat4d getQuat4d(int index) {
return orientations.get(index);
}
public static void getAxisAngle(int index, AxisAngle4f axisAngle) {
axisAngle.set(orientations.get(index));
}
/**
* @param index Index between 0 and {@link #getSphereCount()}-1
* @param axisAngle (Output) modified to contain the index-th sampled orientation
*/
public static void getAxisAngle(int index, AxisAngle4d axisAngle) {
axisAngle.set(orientations.get(index));
}
// Convert from index to position. The sinh scaling tries to compensate
// for the bunching up that occurs when [1 x y z] is projected onto the
// unit sphere.
private static double pind(double ind, double delta, double sigma) {
return (sigma == 0) ? ind * delta : Math.sinh(sigma * ind * delta)
/ sigma;
}
}