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}