001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.symmetry.geometry;
022
023import javax.vecmath.*;
024import java.util.ArrayList;
025import java.util.List;
026
027
028/**
029 * Generate the permutations and sign changes for a Triple.
030 * For instance, (a,0,0) becomes (±a,0,0),(0,±a,0),(0,0,±a). In the general
031 * case 48 tuples are returned.
032 */
033class Permute {
034
035        private List<Point3i> triples = new ArrayList<Point3i>();
036
037        Permute(Point3i t) {
038                Point3i tmp = new Point3i();
039                tmp.x = t.x;
040                tmp.y = t.y;
041                tmp.z = t.z;
042                triples.add(tmp);
043                int n = 1;
044
045                // Unpermuted ±x,±y,±z
046                if (t.x != 0) {
047                        for (int i = 0; i < n; ++i) {
048                                Tuple3i m = triples.get(i);
049
050                                triples.add(new Point3i(-m.x, m.y, m.z));
051                        }
052                        n *= 2;
053                }
054
055                if (t.y != 0) {
056                        for (int i = 0; i < n; ++i) {
057                                Point3i m = triples.get(i);
058                                triples.add(new Point3i(m.x, -m.y, m.z));
059                        }
060                        n *= 2;
061                }
062
063                if (t.z != 0) {
064                        for (int i = 0; i < n; ++i) {
065                                Point3i m = triples.get(i);
066                                triples.add(new Point3i(m.x, m.y, -m.z));
067                        }
068                        n *= 2;
069                }
070                if (t.x == t.y && t.y == t.z) {
071                        return;
072                }
073
074                // At least one differing value
075                for (int i = 0; i < n; ++i) {
076                        Point3i m = triples.get(i);
077                        triples.add(new Point3i(m.y, m.z, m.x));
078                        triples.add(new Point3i(m.z, m.x, m.y));
079                }
080                n *= 3;
081
082                if (t.x == t.y || t.y == t.z) {
083                        return;
084                }
085
086                // At least two differing values
087                for (int i = 0; i < n; ++i) {
088                        Point3i m = triples.get(i);
089                        triples.add(new Point3i(m.y, m.x, m.z));
090                }
091                n *= 2;
092        }
093
094        public int size() {
095                return triples.size();
096        }
097
098        public Point3i get(int i) {
099                return triples.get(i);
100        }
101};
102
103/**
104 * Sample possible orientations. An orientation can be represented as a
105 * quaternion or as a rotation axis and angle. The orientations are somewhat
106 * evenly distributed over orientation space, but the current implementation
107 * is not suited for statistically uniform sampling.
108 * @author Peter
109 */
110public final class SphereSampler {
111
112        private static final List<Quat4d> orientations ;
113
114
115
116        // The rotational symmetries of the cube. (Not normalized, since
117        // PackSet.Add does this.)
118        private static final double[][] cubeSyms = {
119                { 1, 0, 0, 0 },
120                // 180 deg rotations about 3 axes
121                { 0, 1, 0, 0 },
122                { 0, 0, 1, 0 },
123                { 0, 0, 0, 1 },
124                // +/- 120 degree rotations about 4 leading diagonals
125                { 1, 1, 1, 1 }, { 1, 1, 1, -1 }, { 1, 1, -1, 1 }, { 1, 1, -1, -1 },
126                { 1, -1, 1, 1 }, { 1, -1, 1, -1 },
127                { 1, -1, -1, 1 },
128                { 1, -1, -1, -1 },
129                // +/- 90 degree rotations about 3 axes
130                { 1, 1, 0, 0 }, { 1, -1, 0, 0 }, { 1, 0, 1, 0 }, { 1, 0, -1, 0 },
131                { 1, 0, 0, 1 }, { 1, 0, 0, -1 },
132                // 180 degree rotations about 6 face diagonals
133                { 0, 1, 1, 0 }, { 0, 1, -1, 0 }, { 0, 1, 0, 1 }, { 0, 1, 0, -1 },
134                { 0, 0, 1, 1 }, { 0, 0, 1, -1 }, };
135
136
137
138        // Approximate spacing between grid points
139        private static final double delta = 0.15846;
140        // Amount of distortion towards grid corners to account for the projection back to a 4-sphere
141        private static final double sigma = 0.00;
142        private static final int ntot = 7416; // total number of grid points over the 24 cells
143        private static final int ncell = 309; // number of grid points per cell
144        private static final int nent = 18; // number of distinct grid points (with k>=l>=m)
145
146        // Why not (5,5,3) and (5,5,5)? Would they overlap with other cells? -SB
147        private static final int[] k = { 0, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5};
148        private static final int[] l = { 0, 1, 0, 2, 2, 1, 3, 3, 0, 2, 2, 4, 4, 4, 1, 3, 3, 5};
149        private static final int[] m = { 0, 1, 0, 0, 2, 1, 1, 3, 0, 0, 2, 0, 2, 4, 1, 1, 3, 1};
150
151        // number of permutations to expect for each (k,l,m) tuple
152        private static final int[] mult = { 1, 8, 6, 12, 8, 24, 24, 8, 6, 24, 24, 12, 24, 8, 24, 48, 24, 24};
153
154        static
155        {
156
157                List<Quat4d> myorientations = new ArrayList<Quat4d>();
158
159                // 60 equally spaced grid points covering all quaternions
160                // This could be removed, since the 48-cell bcc lattice has lower angle error -SB
161                for (int i = 0; i < IcosahedralSampler.getSphereCount(); i++) {
162                        myorientations.add(IcosahedralSampler.getQuat4d(i));
163                }
164
165                // SB's interpretation of this code:
166                // Directly form a 4D grid over the 4-sphere of orientations. Start by
167                // making a 3D body-centered lattice with w=1. Then use the hypercube symmetries
168                // to extend this grid over the whole surface.
169                // The permuted (k,l,m) values together make a diagonal grid out to ±(5,5,5).
170                // The point spacing is distorted by the pind() function so that the
171                // projection of the points back to the 4-sphere will be more even.
172                
173                // This is the c48u309 lattice from Karney 2006, with a max error of 10.07
174                // degrees.
175                
176                List<Quat4d> grid = new ArrayList<Quat4d>();
177                int ncell1 = 0;
178                for (int n = 0; n < nent; ++n) { // for each tuple (k,l,m) above
179                        Permute p = new Permute(new Point3i(k[n], l[n], m[n]));
180                        assert (mult[n] == p.size());
181                        for (int i = 0; i < mult[n]; ++i) { // for each permutation
182                                Point3i t = p.get(i);
183                                grid.add(new Quat4d(1.0, pind(0.5 * t.x, delta, sigma), pind(
184                                                0.5 * t.y, delta, sigma), pind(0.5 * t.z, delta, sigma)));
185                        }
186                        ncell1 += mult[n];
187                }
188                assert (ncell1 == ncell);
189                // Apply symmetry operations to distribute grid over all values of w
190                int nc = grid.size();
191                assert (nc == ncell);
192                for (int n = 1; n < 24; ++n) {
193                        Quat4d q = new Quat4d(cubeSyms[n][0], cubeSyms[n][1],
194                                        cubeSyms[n][2], cubeSyms[n][3]);
195                        for (int i = 0; i < nc; ++i) {
196                                Quat4d qs = new Quat4d();
197                                qs.mul(q, grid.get(i));
198                                grid.add(qs);
199                                //      s.add(times(q, s.getOrientation(i)), s.getWeight(i)); // this data set has no weights
200                                // Actually, it does have weights but we don't care since we don't need uniform sampling. -SB
201                        }
202                }
203                assert (grid.size() == ntot);
204                myorientations.addAll(grid);
205
206                orientations = myorientations;
207
208        }
209
210        // this class cannot be instantiated
211        private SphereSampler() {
212        };
213
214        public static int getSphereCount() {
215
216                return orientations.size();
217        }
218
219        public static Quat4d getQuat4d(int index) {
220
221                return orientations.get(index);
222        }
223
224        public static void getAxisAngle(int index, AxisAngle4f axisAngle) {
225
226                axisAngle.set(orientations.get(index));
227        }
228
229        /**
230         * @param index Index between 0 and {@link #getSphereCount()}-1
231         * @param axisAngle (Output) modified to contain the index-th sampled orientation
232         */
233        public static void getAxisAngle(int index, AxisAngle4d axisAngle) {
234
235                axisAngle.set(orientations.get(index));
236        }
237
238        // Convert from index to position. The sinh scaling tries to compensate
239        // for the bunching up that occurs when [1 x y z] is projected onto the
240        // unit sphere.
241        private static double pind(double ind, double delta, double sigma) {
242                return (sigma == 0) ? ind * delta : Math.sinh(sigma * ind * delta)
243                                / sigma;
244        }
245
246}