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 */
021
022package org.biojava.nbio.structure.geometry;
023
024import org.biojava.nbio.structure.jama.EigenvalueDecomposition;
025import org.biojava.nbio.structure.jama.Matrix;
026
027import javax.vecmath.Matrix3d;
028import javax.vecmath.Point3d;
029import javax.vecmath.Vector3d;
030
031import java.util.ArrayList;
032import java.util.List;
033
034/**
035 * The moment of inertia, otherwise known as the angular mass or rotational
036 * inertia, of a rigid body determines the torque needed for a desired angular
037 * acceleration about a rotational axis. It depends on the body's mass
038 * distribution and the axis chosen, with larger moments requiring more torque
039 * to change the body's rotation.
040 * <p>
041 * More in https://en.wikipedia.org/wiki/Moment_of_inertia.
042 *
043 * @author Peter Rose
044 * @author Aleix Lafita
045 *
046 */
047public class MomentsOfInertia {
048
049        private List<Point3d> points = new ArrayList<Point3d>();
050        private List<Double> masses = new ArrayList<Double>();
051
052        private boolean modified = true;
053
054        private double[] principalMomentsOfInertia = new double[3];
055        private Vector3d[] principalAxes = new Vector3d[3];
056
057        private Matrix3d orientation = new Matrix3d();
058
059        public enum SymmetryClass {
060                LINEAR, PROLATE, OBLATE, SYMMETRIC, ASYMMETRIC
061        };
062
063        /** Creates a new empty instance of MomentsOfInertia */
064        public MomentsOfInertia() {
065        }
066
067        public void addPoint(Point3d point, double mass) {
068                points.add(point);
069                masses.add(mass);
070                modified = true;
071        }
072
073        public Point3d getCenterOfMass() {
074
075                if (points.size() == 0) {
076                        throw new IllegalStateException(
077                                        "MomentsOfInertia: no points defined");
078                }
079
080                Point3d center = new Point3d();
081                double totalMass = 0.0;
082
083                for (int i = 0, n = points.size(); i < n; i++) {
084                        double mass = masses.get(i);
085                        totalMass += mass;
086                        center.scaleAdd(mass, points.get(i), center);
087                }
088                center.scale(1.0 / totalMass);
089                return center;
090        }
091
092        public double[] getPrincipalMomentsOfInertia() {
093
094                if (modified) {
095                        diagonalizeTensor();
096                        modified = false;
097                }
098                return principalMomentsOfInertia;
099        }
100
101        /**
102         * The principal axes of intertia
103         *
104         * @return
105         */
106        public Vector3d[] getPrincipalAxes() {
107
108                if (modified) {
109                        diagonalizeTensor();
110                        modified = false;
111                }
112                return principalAxes;
113        }
114
115        /**
116         * The orientation Matrix is a 3x3 Matrix with a column for each principal
117         * axis. It represents the orientation (rotation) of the principal axes with
118         * respect to the axes of the coordinate system (unit vectors [1,0,0],
119         * [0,1,0] and [0,0,1]).
120         * <p>
121         * The orientation matrix indicates the rotation to bring the coordinate
122         * axes to the principal axes, in this direction.
123         *
124         * @return the orientation Matrix as a Matrix3d object
125         */
126        public Matrix3d getOrientationMatrix() {
127
128                if (modified) {
129                        diagonalizeTensor();
130                        modified = false;
131                }
132                return orientation;
133        }
134
135        /**
136         * The effective value of this distance for a certain body is known as its
137         * radius of / gyration with respect to the given axis. The radius of
138         * gyration corresponding to Ijj / is defined as /
139         * http://www.eng.auburn.edu/~marghitu/MECH2110/C_4.pdf / radius of gyration
140         * k(j) = sqrt(I(j)/m)
141         */
142        public double[] getElipsisRadii() {
143                if (modified) {
144                        diagonalizeTensor();
145                        modified = false;
146                }
147                double m = 0;
148                for (int i = 0, n = points.size(); i < n; i++) {
149                        m += masses.get(i);
150                }
151                double[] r = new double[3];
152                for (int i = 0; i < 3; i++) {
153                        r[i] = Math.sqrt(principalMomentsOfInertia[i] / m);
154                }
155                return r;
156        }
157
158        public double getRadiusOfGyration() {
159                Point3d c = getCenterOfMass();
160                Point3d t = new Point3d();
161                double sum = 0;
162                for (int i = 0, n = points.size(); i < n; i++) {
163                        t.set(points.get(i));
164                        sum += t.distanceSquared(c);
165                }
166                sum /= points.size();
167                return Math.sqrt(sum);
168        }
169
170        public SymmetryClass getSymmetryClass(double threshold) {
171                if (modified) {
172                        diagonalizeTensor();
173                        modified = false;
174                }
175                double ia = principalMomentsOfInertia[0];
176                double ib = principalMomentsOfInertia[1];
177                double ic = principalMomentsOfInertia[2];
178                boolean c1 = (ib - ia) / (ib + ia) < threshold;
179                boolean c2 = (ic - ib) / (ic + ib) < threshold;
180
181                if (c1 && c2) {
182                        return SymmetryClass.SYMMETRIC;
183                }
184                if (c1) {
185                        return SymmetryClass.OBLATE;
186                }
187                if (c2) {
188                        return SymmetryClass.PROLATE;
189                }
190                return SymmetryClass.ASYMMETRIC;
191        }
192
193        public double symmetryCoefficient() {
194                if (modified) {
195                        diagonalizeTensor();
196                        modified = false;
197                }
198                double ia = principalMomentsOfInertia[0];
199                double ib = principalMomentsOfInertia[1];
200                double ic = principalMomentsOfInertia[2];
201                double c1 = 1.0f - (ib - ia) / (ib + ia);
202                double c2 = 1.0f - (ic - ib) / (ic + ib);
203                return Math.max(c1, c2);
204        }
205
206        public double getAsymmetryParameter(double threshold) {
207                if (modified) {
208                        diagonalizeTensor();
209                        modified = false;
210                }
211                if (getSymmetryClass(threshold).equals(SymmetryClass.SYMMETRIC)) {
212                        return 0.0;
213                }
214                double a = 1.0 / principalMomentsOfInertia[0];
215                double b = 1.0 / principalMomentsOfInertia[1];
216                double c = 1.0 / principalMomentsOfInertia[2];
217                return (2 * b - a - c) / (a - c);
218        }
219
220        public double[][] getInertiaTensor() {
221
222                Point3d p = new Point3d();
223                double[][] tensor = new double[3][3];
224
225                // calculate the inertia tensor at center of mass
226                Point3d com = getCenterOfMass();
227
228                for (int i = 0, n = points.size(); i < n; i++) {
229                        double mass = masses.get(i);
230                        p.sub(points.get(i), com);
231
232                        tensor[0][0] += mass * (p.y * p.y + p.z * p.z);
233                        tensor[1][1] += mass * (p.x * p.x + p.z * p.z);
234                        tensor[2][2] += mass * (p.x * p.x + p.y * p.y);
235
236                        tensor[0][1] -= mass * p.x * p.y;
237                        tensor[0][2] -= mass * p.x * p.z;
238                        tensor[1][2] -= mass * p.y * p.z;
239                }
240
241                tensor[1][0] = tensor[0][1];
242                tensor[2][0] = tensor[0][2];
243                tensor[2][1] = tensor[1][2];
244
245                return tensor;
246        }
247
248        private void diagonalizeTensor() {
249
250                Matrix m = new Matrix(getInertiaTensor());
251                EigenvalueDecomposition eig = m.eig();
252
253                principalMomentsOfInertia = eig.getRealEigenvalues();
254                double[][] eigenVectors = eig.getV().getArray();
255
256                // Get the principal axes from the eigenVectors
257                principalAxes[0] = new Vector3d(eigenVectors[0][0], eigenVectors[1][0],
258                                eigenVectors[2][0]);
259                principalAxes[1] = new Vector3d(eigenVectors[0][1], eigenVectors[1][1],
260                                eigenVectors[2][1]);
261                principalAxes[2] = new Vector3d(eigenVectors[0][2], eigenVectors[1][2],
262                                eigenVectors[2][2]);
263
264                // Convert the principal axes into a rotation matrix
265                for (int i = 0; i < 3; i++)
266                        orientation.setColumn(i, principalAxes[i]);
267                orientation.negate();
268
269        }
270}