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.symmetry.geometry;
023
024import org.biojava.nbio.structure.jama.EigenvalueDecomposition;
025import org.biojava.nbio.structure.jama.Matrix;
026
027import javax.vecmath.Point3d;
028import javax.vecmath.Vector3d;
029import java.util.ArrayList;
030import java.util.List;
031
032/**
033 *
034 * @author Peter
035 */
036public class MomentsOfInertia {
037        private List<Point3d> points = new ArrayList<Point3d>();
038        private List<Double> masses = new ArrayList<Double>();
039
040        private boolean modified = true;
041
042        private double[] principalMomentsOfInertia = new double[3];
043        private Vector3d[] principalAxes = new Vector3d[3];
044
045        public enum SymmetryClass {LINEAR, PROLATE, OBLATE, SYMMETRIC, ASYMMETRIC};
046
047        /** Creates a new instance of MomentsOfInertia */
048        public MomentsOfInertia() {
049        }
050
051        public void addPoint(Point3d point, double mass) {
052                points.add(point);
053                masses.add(mass);
054                modified = true;
055        }
056
057        public Point3d centerOfMass() {
058                if (points.size() == 0) {
059                        throw new IllegalStateException("MomentsOfInertia: no points defined");
060                }
061
062                Point3d center = new Point3d();
063                double totalMass = 0.0;
064
065                for (int i = 0, n = points.size(); i < n; i++) {
066                        double mass = masses.get(i);
067                        totalMass += mass;
068                        center.scaleAdd(mass, points.get(i), center);
069                }
070                center.scale(1.0/totalMass);
071                return center;
072        }
073
074        public double[] getPrincipalMomentsOfInertia() {
075                if (modified) {
076                        diagonalizeTensor();
077                        modified = false;
078                }
079                return principalMomentsOfInertia;
080        }
081
082        public Vector3d[] getPrincipalAxes() {
083                if (modified) {
084                        diagonalizeTensor();
085                        modified = false;
086                }
087                return principalAxes;
088        }
089        // The effective value of this distance for a certain body is known as its radius of
090        // gyration with respect to the given axis. The radius of gyration corresponding to Ijj
091        // is defined as
092        // http://www.eng.auburn.edu/~marghitu/MECH2110/C_4.pdf
093        // radius of gyration k(j) = sqrt(I(j)/m)
094        public double[] getElipsisRadii() {
095                if (modified) {
096                        diagonalizeTensor();
097                        modified = false;
098                }
099                double m = 0;
100                for (int i = 0, n = points.size(); i < n; i++) {
101                         m += masses.get(i);
102                }
103                double[] r = new double[3];
104                for (int i = 0; i < 3; i++) {
105                        r[i] = Math.sqrt(principalMomentsOfInertia[i]/m);
106                }
107                return r;
108        }
109
110        public double getRadiusOfGyration() {
111                Point3d c = centerOfMass();
112                Point3d t = new Point3d();
113                double sum = 0;
114                for (int i = 0, n = points.size(); i < n; i++) {
115                        t.set(points.get(i));
116                        sum += t.distanceSquared(c);
117                }
118                sum /= points.size();
119                return Math.sqrt(sum);
120        }
121
122        public SymmetryClass getSymmetryClass(double threshold) {
123                if (modified) {
124                        diagonalizeTensor();
125                        modified = false;
126                }
127                double ia = principalMomentsOfInertia[0];
128                double ib = principalMomentsOfInertia[1];
129                double ic = principalMomentsOfInertia[2];
130                boolean c1 = (ib - ia) / (ib + ia) < threshold;
131                boolean c2 = (ic - ib) / (ic + ib) < threshold;
132
133                if (c1 && c2) {
134                        return SymmetryClass.SYMMETRIC;
135                }
136                if (c1) {
137                        return SymmetryClass.OBLATE;
138                }
139                if (c2) {
140                        return SymmetryClass.PROLATE;
141                }
142                return SymmetryClass.ASYMMETRIC;
143        }
144
145        public double symmetryCoefficient() {
146                if (modified) {
147                        diagonalizeTensor();
148                        modified = false;
149                }
150                double ia = principalMomentsOfInertia[0];
151                double ib = principalMomentsOfInertia[1];
152                double ic = principalMomentsOfInertia[2];
153                double c1 = 1.0f - (ib - ia) / (ib + ia);
154                double c2 = 1.0f - (ic - ib) / (ic + ib);
155                return Math.max(c1, c2);
156        }
157
158        public double getAsymmetryParameter(double threshold) {
159           if (modified) {
160                        diagonalizeTensor();
161                        modified = false;
162                }
163           if (getSymmetryClass(threshold).equals(SymmetryClass.SYMMETRIC)) {
164                   return 0.0;
165           }
166           double a = 1.0/principalMomentsOfInertia[0];
167           double b = 1.0/principalMomentsOfInertia[1];
168           double c = 1.0/principalMomentsOfInertia[2];
169           return (2 * b - a - c) / (a - c);
170        }
171
172        public double[][] getInertiaTensor() {
173                Point3d p = new Point3d();
174                double[][] tensor = new double[3][3];
175
176                // calculate the inertia tensor at center of mass
177                Point3d com = centerOfMass();
178
179                for (int i = 0, n = points.size(); i < n; i++) {
180                        double mass = masses.get(i);
181                        p.sub(points.get(i), com);
182
183                        tensor[0][0] += mass * (p.y * p.y + p.z * p.z);
184                        tensor[1][1] += mass * (p.x * p.x + p.z * p.z);
185                        tensor[2][2] += mass * (p.x * p.x + p.y * p.y);
186
187                        tensor[0][1] -= mass * p.x * p.y;
188                        tensor[0][2] -= mass * p.x * p.z;
189                        tensor[1][2] -= mass * p.y * p.z;
190                }
191
192                tensor[1][0] = tensor[0][1];
193                tensor[2][0] = tensor[0][2];
194                tensor[2][1] = tensor[1][2];
195
196                return tensor;
197        }
198
199        private void diagonalizeTensor() {
200                Matrix m = new Matrix(getInertiaTensor());
201                EigenvalueDecomposition eig = m.eig();
202
203                principalMomentsOfInertia = eig.getRealEigenvalues();
204                double[][] eigenVectors = eig.getV().getArray();
205                principalAxes[0] = new Vector3d(eigenVectors[0][0], eigenVectors[1][0],eigenVectors[2][0]);
206                principalAxes[1] = new Vector3d(eigenVectors[0][1], eigenVectors[1][1],eigenVectors[2][1]);
207                principalAxes[2] = new Vector3d(eigenVectors[0][2], eigenVectors[1][2],eigenVectors[2][2]);
208        }
209}