* @author John L. Moreland * @see edu.sdsc.mbt.Bond * @see edu.sdsc.mbt.util.BondFactory */ public class Algebra { public static final int X_AXIS = 0; public static final int Y_AXIS = 1; public static final int Z_AXIS = 2; /** * Compute the distance between two virtex coordinates. *
*
* v1--v2 ** @param v1 First virtex coordinate as double[3]. * @param v2 Second virtex coordinate as double[3]. * @return the computed distence. */ public static double distance( double v1[], double v2[] ) { return Math.sqrt( (v2[0] - v1[0]) * (v2[0] - v1[0]) + (v2[1] - v1[1]) * (v2[1] - v1[1]) + (v2[2] - v1[2]) * (v2[2] - v1[2]) ); } /** * Calculate the angle formed by three virtex coordinates. *
*
* v1 v3 * \ / * v2 **
* @param v1 An end virtex coordinate as double[3]. * @param v2 The middle virtex coordinate as double[3]. * @param v3 An end virtex coordinate as double[3]. * @return angle in degrees. */ public static double angle( double v1[], double v2[], double v3[] ) { // See: http://mathworld.wolfram.com/Line-LineAngle.html // Compute the numerator. double a = (v1[0] - v2[0]) * (v3[0] - v2[0]) + (v1[1] - v2[1]) * (v3[1] - v2[1]) + (v1[2] - v2[2]) * (v3[2] - v2[2]); // Compute the distance between v1 and v2. double b = Math.sqrt( (v1[0] - v2[0]) * (v1[0] - v2[0]) + (v1[1] - v2[1]) * (v1[1] - v2[1]) + (v1[2] - v2[2]) * (v1[2] - v2[2]) ); // Compute the distance between v3 and v2. double c = Math.sqrt( (v3[0] - v2[0]) * (v3[0] - v2[0]) + (v3[1] - v2[1]) * (v3[1] - v2[1]) + (v3[2] - v2[2]) * (v3[2] - v2[2]) ); // Calculate the angle between the two vectors. double cos = a / (b * c); double rad = Math.acos( cos ); double degrees = Math.toDegrees( rad ); return degrees; } /** * Calculate the dihedral angle defined by four virtex coordinates. *
*
* v1 v4 * \ / * v2--v3 **
* @param v1 An end virtex coordinate as double[3]. * @param v2 An interior virtex coordinate as double[3]. * @param v3 An interior virtex coordinate as double[3]. * @param v4 An end virtex coordinate as double[3]. * @return dihedral angle in degrees. */ public static double dihedralAngle( double v1[], double v2[], double v3[], double v4[] ) { // See: http://mathworld.wolfram.com/DihedralAngle.html // Calculate the normal components for plane 1 (defined by v1-v2-v3). double a1 = (v3[1] - v2[1]) * (v1[2] - v2[2]) - (v1[1] - v2[1]) * (v3[2] - v2[2]); double b1 = (v1[0] - v2[0]) * (v3[2] - v2[2]) - (v3[0] - v2[0]) * (v1[2] - v2[2]); double c1 = (v3[0] - v2[0]) * (v1[1] - v2[1]) - (v1[0] - v2[0]) * (v3[1] - v2[1]); // Calculate the normal components for plane 2 (defined by v2-v3-v4). double a2 = (v3[1] - v2[1]) * (v4[2] - v2[2]) - (v4[1] - v2[1]) * (v3[2] - v2[2]); double b2 = (v4[0] - v2[0]) * (v3[2] - v2[2]) - (v3[0] - v2[0]) * (v4[2] - v2[2]); double c2 = (v3[0] - v2[0]) * (v4[1] - v2[1]) - (v4[0] - v2[0]) * (v3[1] - v2[1]); // Calculate the dihedral angle between the two plane's normals. double a = a1*a2 + b1*b2 + c1*c2; double b = Math.sqrt( a1*a1 + b1*b1 + c1*c1 ); double c = Math.sqrt( a2*a2 + b2*b2 + c2*c2 ); double cos = a / (b * c); double rad = Math.acos( cos ); double degrees = Math.toDegrees( rad ); //get the angle //first, compute the back bond normal double[] center = new double[]{ v2[0] - v3[0], v2[1] - v3[1], v2[2] - v3[2] }; double[] back = new double[]{ v4[0] - v3[0], v4[1] - v3[1], v4[2] - v3[2] }; double[] normal = Algebra.crossProduct(center, back); //compute the projection of the front bond to get the correct angle between the projection and //the normal: if the angle < 90, then the angle is positive double[] bond = new double[]{ v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2] }; double[] cross = Algebra.crossProduct(center, bond); double[] projection = Algebra.crossProduct(cross, center); //now get the angle if (Algebra.getAngleBetweenVectors(projection, normal) > 90){ degrees = -degrees; } return degrees; } /** * Compute the midpoint between two points. */ public static Vec3d midpoint(Vec3d a, Vec3d b){ double x = a.value[0] + (b.value[0] - a.value[0])/2; double y = a.value[1] + (b.value[1] - a.value[1])/2; double z = a.value[2] + (b.value[2] - a.value[2])/2; return new Vec3d(x, y, z); } /** * Compute the distance between coordinates using double values */ public static double distance( Vec3d v1, Vec3d v2 ) { double dist = 0.0; dist += (v2.value[0] - v1.value[0])*(v2.value[0] - v1.value[0]); dist += (v2.value[1] - v1.value[1])*(v2.value[1] - v1.value[1]); dist += (v2.value[2] - v1.value[2])*(v2.value[2] - v1.value[2]); dist = Math.sqrt( dist ); return dist; } /** * Compute the distance between coordinates using float values */ public static float distance( Vec3f v1, Vec3f v2 ) { float dist = 0.0f; dist += (v2.value[0] - v1.value[0])*(v2.value[0] - v1.value[0]); dist += (v2.value[1] - v1.value[1])*(v2.value[1] - v1.value[1]); dist += (v2.value[2] - v1.value[2])*(v2.value[2] - v1.value[2]); dist = (float)Math.sqrt( dist ); return dist; } /** * Returns coordinates of the given Atom in 3d space. * @param atom * @return */ public static Vec3d getAtomCoordinates(Atom atom){ double x = atom.coordinate[0]; double y = atom.coordinate[1]; double z = atom.coordinate[2]; return new Vec3d(x, y, z); } /** * Calculates the angle formed by the three points. * @param v2 coordinates of the middle point * @param v1 coordinates of an end point * @param v3 coordinates of an end point * @return angle in degrees */ public static double getAngle(Vec3d v2, Vec3d v1, Vec3d v3){ //get the angle based on the coordinates of three points expressed in whole degrees //v1 is the central atom at the vertex double a = (v2.value[0] - v1.value[0])*(v3.value[0] - v1.value[0]) + (v2.value[1] - v1.value[1])*(v3.value[1] - v1.value[1]) + (v2.value[2] - v1.value[2])*(v3.value[2] - v1.value[2]); double b = Math.sqrt((v2.value[0] - v1.value[0])*(v2.value[0] - v1.value[0]) + (v2.value[1] - v1.value[1])*(v2.value[1] - v1.value[1]) + (v2.value[2] - v1.value[2])*(v2.value[2] - v1.value[2])); double c = Math.sqrt((v3.value[0] - v1.value[0])*(v3.value[0] - v1.value[0]) + (v3.value[1] - v1.value[1])*(v3.value[1] - v1.value[1]) + (v3.value[2] - v1.value[2])*(v3.value[2] - v1.value[2])); double cos = a/(b*c); double rad = Math.acos(cos); double degrees = Math.toDegrees(rad); return degrees; } /** * Calculates the dihedral angle defined by four points * @param v3 coordinates of an end point * @param v1 coordinates of the second point in the sequence * @param v2 coordinates of the third point * @param v4 coordinates of the other end point * @return */ public static double getDihedralAngle(Vec3d v3, Vec3d v1, Vec3d v2, Vec3d v4){ //calculate the components of normals to the planes double a1 = (v2.value[1] - v1.value[1])*(v3.value[2] - v1.value[2]) - (v3.value[1] - v1.value[1])*(v2.value[2] - v1.value[2]); double b1 = (v3.value[0] - v1.value[0])*(v2.value[2] - v1.value[2]) - (v2.value[0] - v1.value[0])*(v3.value[2] - v1.value[2]); double c1 = (v2.value[0] - v1.value[0])*(v3.value[1] - v1.value[1]) - (v3.value[0] - v1.value[0])*(v2.value[1] - v1.value[1]); double a2 = (v2.value[1] - v1.value[1])*(v4.value[2] - v1.value[2]) - (v4.value[1] - v1.value[1])*(v2.value[2] - v1.value[2]); double b2 = (v4.value[0] - v1.value[0])*(v2.value[2] - v1.value[2]) - (v2.value[0] - v1.value[0])*(v4.value[2] - v1.value[2]); double c2 = (v2.value[0] - v1.value[0])*(v4.value[1] - v1.value[1]) - (v4.value[0] - v1.value[0])*(v2.value[1] - v1.value[1]); double a = a1*a2 + b1*b2 + c1*c2; double b = Math.sqrt(a1*a1 + b1*b1 + c1*c1); double c = Math.sqrt(a2*a2 + b2*b2 + c2*c2); double cos = a/(b*c); double rad = Math.acos(cos); double degrees = Math.toDegrees(rad); return degrees; } /** * Compute the vector cross product between v1 and v2 storing the * output in result. *
* @param v1 First vector argument. * @param v2 Second vector argument. * @param result The vector in which the result will be stored. */ public static void crossProduct( double[] v1, double[] v2, double[] result ) { if ( v1.length != 3 ) throw new IllegalArgumentException( "v1.length != 3" ); if ( v2.length != 3 ) throw new IllegalArgumentException( "v2.length != 3" ); if ( result.length != 3 ) throw new IllegalArgumentException( "result.length != 3" ); if ( v1 == v2 ) throw new IllegalArgumentException( "v1 == v2" ); if ( result == v1 ) throw new IllegalArgumentException( "result == v1" ); if ( result == v2 ) throw new IllegalArgumentException( "result == v2" ); result[0] = + ( ( v1[1] * v2[2] ) - ( v2[1] * v1[2] ) ); result[1] = - ( ( v1[0] * v2[2] ) - ( v2[0] * v1[2] ) ); result[2] = + ( ( v1[0] * v2[1] ) - ( v2[0] * v1[1] ) ); } public static void crossProduct( float[] v1, float[] v2, float[] result ) { if ( v1.length != 3 ) throw new IllegalArgumentException( "v1.length != 3" ); if ( v2.length != 3 ) throw new IllegalArgumentException( "v2.length != 3" ); if ( result.length != 3 ) throw new IllegalArgumentException( "result.length != 3" ); if ( v1 == v2 ) throw new IllegalArgumentException( "v1 == v2" ); if ( result == v1 ) throw new IllegalArgumentException( "result == v1" ); if ( result == v2 ) throw new IllegalArgumentException( "result == v2" ); result[0] = + ( ( v1[1] * v2[2] ) - ( v2[1] * v1[2] ) ); result[1] = - ( ( v1[0] * v2[2] ) - ( v2[0] * v1[2] ) ); result[2] = + ( ( v1[0] * v2[1] ) - ( v2[0] * v1[1] ) ); } public static double[] crossProduct( double[] v1, double[] v2 ) { double[] result = new double[3]; if ( v1.length != 3 ) throw new IllegalArgumentException( "v1.length != 3" ); if ( v2.length != 3 ) throw new IllegalArgumentException( "v2.length != 3" ); result[0] = + ( ( v1[1] * v2[2] ) - ( v2[1] * v1[2] ) ); result[1] = - ( ( v1[0] * v2[2] ) - ( v2[0] * v1[2] ) ); result[2] = + ( ( v1[0] * v2[1] ) - ( v2[0] * v1[1] ) ); return result; } /** * Return the vector dot product between v1 and v2. *
* @param v1 First vector argument.
* @param v2 Second vector argument.
* @return The dot product value.
*/
public static double dotProduct( double[] v1, double[] v2 )
{
if ( v1.length != v2.length )
throw new IllegalArgumentException( "v1.length != v2.length" );
double sumOfSquares = 0.0;
for ( int i=0; i
* @param matrix The 16-element rotation matrix.
* @param point The point to be rotated in-place.
*/
public static void matrixRotate( double[] matrix, double[] point )
{
if ( matrix == null )
throw new NullPointerException( "null matrix" );
if ( matrix.length < 16 )
throw new IllegalArgumentException( "matrix.length < 16" );
if ( point == null )
throw new NullPointerException( "null point" );
if ( point.length < 3 )
throw new IllegalArgumentException( "point.length < 3" );
double[] m = matrix;
double x = point[0];
double y = point[1];
double z = point[2];
double w = 1.0;
point[0] = m[0]*x + m[1]*y + m[2]*z + m[3]*w;
point[1] = m[4]*x + m[5]*y + m[6]*z + m[7]*w;
point[2] = m[8]*x + m[9]*y + m[10]*z + m[11]*w;
}
/**
* Convert the angle-axis {angle,x,y,z} rotation vector to a 16-element
* rotation matrix.
*
* @param rotation The angle-axis rotation vector.
* @param matrix The resulting 16-element rotation matrix.
*/
public static void rotationToMatrix( double[] rotation, double[] matrix )
{
if ( rotation == null )
throw new NullPointerException( "null rotation" );
if ( rotation.length < 4 )
throw new IllegalArgumentException( "rotation.length < 4" );
if ( matrix == null )
throw new NullPointerException( "null matrix" );
if ( matrix.length < 16 )
throw new IllegalArgumentException( "matrix.length < 16" );
// Make sure that the rotation vector is unit length.
double len = Math.sqrt(
rotation[1] * rotation[1] +
rotation[2] * rotation[2] +
rotation[3] * rotation[3] );
double norm = 1.0 / len;
if ( norm == 0.0 ) norm = 1.0;
double a = rotation[0];
double x = rotation[1] * norm;
double y = rotation[2] * norm;
double z = rotation[3] * norm;
double c = Math.cos( a );
double s = Math.sin( a );
double t = 1.0 - c;
matrix[0] = t*x*x + c;
matrix[1] = t*x*y - s*z;
matrix[2] = t*x*z + s*y;
matrix[3] = 0.0;
matrix[4] = t*x*y + s*z;
matrix[5] = t*y*y + c;
matrix[6] = t*y*z - s*x;
matrix[7] = 0.0;
matrix[8] = t*x*z - s*y;
matrix[9] = t*y*z + s*x;
matrix[10] = t*z*z + c;
matrix[11] = 0.0;
matrix[12] = 0.0;
matrix[13] = 0.0;
matrix[14] = 0.0;
matrix[15] = 1.0;
}
/**
* Calculate the position of the view point on Z axis, given the center of the
* viewed structure and the required distance
* @param target
* @param distance
* @return
*/
public static double[] getViewPoint(double[] target, double distance){
double result = target[2] + Math.sqrt(Math.pow(distance, 2) - Math.pow(target[0], 2) - Math.pow(target[1], 2) );
return new double[] { 0, 0, result };
}
public static double[] multiplyMatrices(double[][] a, double[][] b){
double[] A = new double[16];
double[] B = new double[16];
A[0] = a[0][0];
A[1] = a[0][1];
A[2] = a[0][2];
A[3] = 0;
A[4] = a[1][0];
A[5] = a[1][1];
A[6] = a[1][2];
A[7] = 0;
A[8] = a[2][0];
A[9] = a[2][1];
A[10] = a[2][2];
A[11] = 0;
A[12] = 0;
A[13] = 0;
A[14] = 0;
A[15] = 1;
B[0] = b[0][0];
B[1] = b[0][1];
B[2] = b[0][2];
B[3] = 0;
B[4] = b[1][0];
B[5] = b[1][1];
B[6] = b[1][2];
B[7] = 0;
B[8] = b[2][0];
B[9] = b[2][1];
B[10] = b[2][2];
B[11] = 0;
B[12] = 0;
B[13] = 0;
B[14] = 0;
B[15] = 1;
return Algebra.multiplySquareMatrices(A, B);
}
/**
* Calculate a product of two matrices and return the result
* @param first multiplier
* @param second multiplier matrix
* @return result
*/
public static double[] multiplySquareMatrices(double[] a, double[] b){
if (a.length != b.length) throw new IllegalArgumentException("The matrices are not of equal size.");
if (a.length != 16) throw new IllegalArgumentException("Matrix length must be 16.");
double[] out = new double[16];
out[0] = a[0]*b[0] + a[4]*b[1] + a[8]*b[2] + a[12]*b[3];
out[1] = a[1]*b[0] + a[5]*b[1] + a[9]*b[2] + a[13]*b[3];
out[2] = a[2]*b[0] + a[6]*b[1] + a[10]*b[2] + a[14]*b[3];
out[3] = a[3]*b[0] + a[7]*b[1] + a[11]*b[2] + a[15]*b[3];
out[4] = a[0]*b[4] + a[4]*b[5] + a[8]*b[6] + a[12]*b[7];
out[5] = a[1]*b[4] + a[5]*b[5] + a[9]*b[6] + a[13]*b[7];
out[6] = a[2]*b[4] + a[6]*b[5] + a[10]*b[6] + a[14]*b[7];
out[7] = a[3]*b[4] + a[7]*b[5] + a[11]*b[6] + a[15]*b[7];
out[8] = a[0]*b[8] + a[4]*b[9] + a[8]*b[10] + a[12]*b[11];
out[9] = a[1]*b[8] + a[5]*b[9] + a[9]*b[10] + a[13]*b[11];
out[10] = a[2]*b[8] + a[6]*b[9] + a[10]*b[10] + a[14]*b[11];
out[11] = a[3]*b[8] + a[7]*b[9] + a[11]*b[10] + a[15]*b[11];
out[12] = a[0]*b[12] + a[4]*b[13] + a[8]*b[14] + a[12]*b[15];
out[13] = a[1]*b[12] + a[5]*b[13] + a[9]*b[14] + a[13]*b[15];
out[14] = a[2]*b[12] + a[6]*b[13] + a[10]*b[14] + a[14]*b[15];
out[15] = a[3]*b[12] + a[7]*b[13] + a[11]*b[14] + a[15]*b[15];
return out;
}
/**
* This method solves a system of three equations with three unknowns and returns
* an array of x,y,z values
* @param system linear array of all 12 constituents: a1, b1, c1, d1, a2, b2.. etc
* @return
*/
public static double[] solveSystemOfEquations(double[] system){
double[] out = new double[3];
double a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3;
a1 = system[0];
b1 = system[1];
c1 = system[2];
d1 = system[3];
a2 = system[4];
b2 = system[5];
c2 = system[6];
d2 = system[7];
a3 = system[8];
b3 = system[9];
c3 = system[10];
d3 = system[11];
/* System.out.println("a1 = " + a1);
System.out.println("b1 = " + b1);
System.out.println("c1 = " + c1);
System.out.println("d1 = " + d1);
System.out.println("a2 = " + a2);
System.out.println("b2 = " + b2);
System.out.println("c2 = " + c2);
System.out.println("d2 = " + d2);
System.out.println("a3 = " + a3);
System.out.println("b3 = " + b3);
System.out.println("c3 = " + c3);
System.out.println("d3 = " + d3);
*/
//solve the system using the matrix method
double determinant = a1*b2*c3 - a1*b3*c2 + a2*b3*c1 - a2*b1*c3 + a3*b1*c2 - a3*b2*c1;
if (determinant == 0) return null;
//find x
out[0] = (d1*b2*c3 - d1*b3*c2 + d2*b3*c1 - d2*b1*c3 + d3*b1*c2 - d3*b2*c1)/determinant;
out[1] = (a1*d2*c3 - a1*d3*c2 + a2*d3*c1 - a2*d1*c3 + a3*d1*c2 - a3*d2*c1)/determinant;
out[2] = (a1*b2*d3 - a1*b3*d2 + a2*b3*d1 - a2*b1*d3 + a3*b1*d2 - a3*b2*d1)/determinant;
return out;
}
public static double getAngleBetweenVectors(double[] a, double[] b){
double[] ta = getNormalizedVector(a);
double[] tb = getNormalizedVector(b);
double cos = dotProduct(ta, tb)/(vectorLength(ta) * vectorLength(tb));
double angle = Math.acos(cos);
return Math.toDegrees(angle);
}
public static double[] getVector(double[] start, double[] end){
double[] out = new double[3];
out[0] = end[0] - start[0];
out[1] = end[1] - start[1];
out[2] = end[2] - start[2];
return out;
}
public static double[] getNormalizedVector(double[] start, double[] end){
double[] out = new double[3];
out[0] = end[0] - start[0];
out[1] = end[1] - start[1];
out[2] = end[2] - start[2];
return Algebra.getNormalizedVector(out);
}
public static void main(String[] args){
/*
double[] a = new double[]{ 1.5470, -0.0260 , -0.0000};
double[] b = new double[]{ 1.9350 , 1.4340, -0.0000};
double[] out = Algebra.getNormalizedVector(a, b);
System.out.println(out[0] + "\t" + out[1] + "\t" + out[2]);
*/
double[][] a = new double[3][3];
a[0][0] = 1;
a[1][0] = 1;
a[2][0] = 0;
a[0][1] = 1;
a[1][1] = 2;
a[2][1] = 0;
a[0][2] = 1;
a[1][2] = 2;
a[2][2] = 1;
double[][] res = Algebra.getInverseMatrix(a);
for (int i = 0; i < 3; i++){
for (int j = 0; j < 3; j++){
System.out.println("res[" + i + "][" + j + "] = " + res[i][j]);
}
}
}
public static double[][] getInverseMatrix(double[][] a){
if (a[0].length != 3 || a[1].length != 3 || a[2].length != 3 || a.length != 3) return null;
//compute determinant
double det = a[0][0]*(a[1][1]*a[2][2] - a[2][1]*a[1][2]) - a[0][1]*(a[1][0]*a[2][2] - a[2][0]*a[1][2]) +
a[0][2]*(a[1][0]*a[2][1] - a[2][0]*a[1][1]);
System.out.println("det = " + det);
double[][] out = new double[3][3];
//calculate the elements of the new matrix
out[0][0] = (a[1][1]*a[2][2] - a[2][1]*a[1][2])/det;
out[0][1] = -(a[1][0]*a[2][2] - a[2][0]*a[1][2])/det;
out[0][2] = (a[1][0]*a[2][1] - a[2][0]*a[1][1])/det;
out[1][0] = -(a[0][1]*a[2][2] - a[2][1]*a[0][2])/det;
out[1][1] = (a[0][0]*a[2][2] - a[2][0]*a[0][2])/det;
out[1][2] = -(a[0][0]*a[2][1] - a[2][0]*a[0][1])/det;
out[2][0] = (a[0][1]*a[1][2] - a[1][1]*a[0][2])/det;
out[2][1] = -(a[0][0]*a[1][2] - a[1][0]*a[0][2])/det;
out[2][2] = (a[0][0]*a[1][1] - a[1][0]*a[0][1])/det;
return out;
}
/**
* Multiply a vector times a scalar.
*
* @param vector The vector to be scaled in-place.
*/
public static void scalarMultiply( float[] vector, float scalar )
{
for ( int i=0; i
* @param vector1 The vector to be modified.
* @param vector2 The vector to subtract.
*/
public static void subtractVectors( double[] vector1, double[] vector2 )
{
if ( vector1 == null )
throw new NullPointerException( "null vector1" );
if ( vector2 == null )
throw new NullPointerException( "null vector2" );
if ( vector1.length != vector2.length )
throw new IllegalArgumentException( "vector length mismatch" );
for ( int i=0; i