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}