Java Code Examples for javax.vecmath.AxisAngle4d#set()

/**
 * Given a rotation matrix calculates the rotation axis and angle for it.
 * The angle is calculated from the trace, the axis from the eigenvalue
 * decomposition.
 * If given matrix is improper rotation or identity matrix then
 * axis (0,0,0) and angle 0 are returned.
 * @param m
 * @return
 * @throws IllegalArgumentException if given matrix is not a rotation matrix (determinant not 1 or -1)
 */
public static AxisAngle4d getRotAxisAndAngle(Matrix3d m) {
	double determinant = m.determinant();

	if (!(Math.abs(determinant)-1.0<DELTA)) throw new IllegalArgumentException("Given matrix is not a rotation matrix");

	AxisAngle4d axisAndAngle = new AxisAngle4d(new Vector3d(0,0,0),0);

	double[] d = {m.m00,m.m10,m.m20,
			m.m01,m.m11,m.m21,
			m.m02,m.m12,m.m22};

	Matrix r = new Matrix(d,3);

	if (!deltaComp(r.det(), 1.0, DELTA)) {
		// improper rotation: we return axis 0,0,0 and angle 0
		return axisAndAngle;
	}

	EigenvalueDecomposition evd = new EigenvalueDecomposition(r);

	Matrix eval = evd.getD();
	if (deltaComp(eval.get(0, 0),1.0,DELTA) && deltaComp(eval.get(1, 1),1.0,DELTA) && deltaComp(eval.get(2, 2),1.0,DELTA)) {
		// the rotation is an identity: we return axis 0,0,0 and angle 0
		return axisAndAngle;
	}
	int indexOfEv1;
	for (indexOfEv1=0;indexOfEv1<3;indexOfEv1++) {
		if (deltaComp(eval.get(indexOfEv1, indexOfEv1),1,DELTA)) break;
	}
	Matrix evec = evd.getV();
	axisAndAngle.set(new Vector3d(evec.get(0,indexOfEv1), evec.get(1, indexOfEv1), evec.get(2, indexOfEv1)),
			Math.acos((eval.trace()-1.0)/2.0));

	return axisAndAngle;
}
 

public static void getAxisAngle(int index, AxisAngle4d axisAngle) {
	quat.set(orientations[index]);
	axisAngle.set(quat);
}
 

private boolean evaluatePermutation(List<Integer> permutation) {
	// permutate subunits
	for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) {
		transformedCoords[j].set(originalCoords[permutation.get(j)]);
	}

	int fold = PermutationGroup.getOrder(permutation);

	// TODO implement this piece of code using at origin superposition
	Quat4d quat = UnitQuaternions.relativeOrientation(
			originalCoords, transformedCoords);
	AxisAngle4d axisAngle = new AxisAngle4d();
	Matrix4d transformation = new Matrix4d();

	transformation.set(quat);
	axisAngle.set(quat);

	Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
	if (axis.lengthSquared() < 1.0E-6) {
		axisAngle.x = 0;
		axisAngle.y = 0;
		axisAngle.z = 1;
		axisAngle.angle = 0;
	} else {
		axis.normalize();
		axisAngle.x = axis.x;
		axisAngle.y = axis.y;
		axisAngle.z = axis.z;
	}

	CalcPoint.transform(transformation, transformedCoords);

	double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);

	if (subunitRmsd <parameters.getRmsdThreshold()) {
		// transform to original coordinate system
		combineWithTranslation(transformation);
		QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
		if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) {
			return false;
		}

		scores.setRmsdCenters(subunitRmsd);
		Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
		rotations.addRotation(symmetryOperation);
		return true;
	}
	return false;
}
 

private void solve() {
		initialize();
		Vector3d trans = new Vector3d(subunits.getCentroid());
		trans.negate();
		List<Point3d[]> traces = subunits.getTraces();

//		Point3d[] x = SuperPosition.clonePoint3dArray(traces.get(0));
//		SuperPosition.center(x);
//		Point3d[] y = SuperPosition.clonePoint3dArray(traces.get(1));
//		SuperPosition.center(y);

		Point3d[] x = CalcPoint.clonePoint3dArray(traces.get(0));
		CalcPoint.translate(trans, x);
		Point3d[] y = CalcPoint.clonePoint3dArray(traces.get(1));
		CalcPoint.translate(trans, y);

		// TODO implement this piece of code using at origin superposition
		Quat4d quat = UnitQuaternions.relativeOrientation(
				x, y);
		AxisAngle4d axisAngle = new AxisAngle4d();
		Matrix4d transformation = new Matrix4d();

		transformation.set(quat);
		axisAngle.set(quat);

		Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
		if (axis.lengthSquared() < 1.0E-6) {
			axisAngle.x = 0;
			axisAngle.y = 0;
			axisAngle.z = 1;
			axisAngle.angle = 0;
		} else {
			axis.normalize();
			axisAngle.x = axis.x;
			axisAngle.y = axis.y;
			axisAngle.z = axis.z;
		}

		CalcPoint.transform(transformation, y);

		// if rmsd or angle deviation is above threshold, stop
		double angleThresholdRadians = Math.toRadians(parameters.getAngleThreshold());
		double deltaAngle = Math.abs(Math.PI-axisAngle.angle);

		if (deltaAngle > angleThresholdRadians) {
			rotations.setC1(subunits.getSubunitCount());
			return;
		}

		// add unit operation
		addEOperation();

		// add C2 operation
		int fold = 2;
		combineWithTranslation(transformation);
		List<Integer> permutation = Arrays.asList(1,0);
		QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
		scores.setRmsdCenters(0.0); // rmsd for superposition of two subunits centers is zero by definition

		if (scores.getRmsd() > parameters.getRmsdThreshold() || deltaAngle > angleThresholdRadians) {
			rotations.setC1(subunits.getSubunitCount());
			return;
		}

		Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
		rotations.addRotation(symmetryOperation);
	}
 

/**
 * Returns the AxisAngle of the helix transformation
 * @param transformation helix transformation
 * @return
 */
public AxisAngle4d getAxisAngle() {
	AxisAngle4d axis = new AxisAngle4d();
	axis.set(this.transformation);
	return axis;
}
 

/**
 * Superimpose subunits based on the given permutation. Then check whether
 * the superposition passes RMSD thresholds and create a Rotation to
 * represent it if so.
 * @param permutation A list specifying which subunits should be aligned by the current transformation
 * @return A Rotation representing the permutation, or null if the superposition did not meet thresholds.
 */
private Rotation superimposePermutation(List<Integer> permutation) {
	// permutate subunits
	for (int j = 0, n = subunits.getSubunitCount(); j < n; j++) {
		transformedCoords[j].set(originalCoords[permutation.get(j)]);
	}

	int fold = PermutationGroup.getOrder(permutation);

	// get optimal transformation and axisangle by subunit superposition
	// TODO implement this piece of code using at origin superposition
	Quat4d quat = UnitQuaternions.relativeOrientation(
			originalCoords, transformedCoords);
	AxisAngle4d axisAngle = new AxisAngle4d();
	Matrix4d transformation = new Matrix4d();

	transformation.set(quat);
	axisAngle.set(quat);

	Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
	if (axis.lengthSquared() < 1.0E-6) {
		axisAngle.x = 0;
		axisAngle.y = 0;
		axisAngle.z = 1;
		axisAngle.angle = 0;
	} else {
		axis.normalize();
		axisAngle.x = axis.x;
		axisAngle.y = axis.y;
		axisAngle.z = axis.z;
	}

	CalcPoint.transform(transformation, transformedCoords);
	double subunitRmsd = CalcPoint.rmsd(transformedCoords, originalCoords);

	if (subunitRmsd < parameters.getRmsdThreshold()) {
		combineWithTranslation(transformation);

		// evaluate superposition of CA traces
		QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
		if (scores.getRmsd() < 0.0 || scores.getRmsd() > parameters.getRmsdThreshold()) {
			return null;
		}

		scores.setRmsdCenters(subunitRmsd);
		Rotation symmetryOperation = createSymmetryOperation(permutation, transformation, axisAngle, fold, scores);
		return symmetryOperation;
	}
	return null;
}
 

/**
 * Test {@link UnitQuaternions#orientation(javax.vecmath.Point3d[])}.
 * <p>
 * Tests the identity orientation, orientation around one coordinate axis
 * and orientation around a non-coordinate axis.
 *
 * @throws StructureException
 * @throws IOException
 */
@Test
public void testOrientation() throws IOException, StructureException {

	// Get points from a structure. It is difficult to generate points
	// with no bias in their distribution (too uniform, ie).
	Structure pdb = StructureIO.getStructure("4hhb.A");
	Point3d[] cloud = Calc.atomsToPoints(StructureTools
			.getRepresentativeAtomArray(pdb));

	// Center the cloud at the origin
	CalcPoint.center(cloud);

	// Orient its principal axes to the coordinate axis
	Quat4d orientation = UnitQuaternions.orientation(cloud);
	Matrix4d transform = new Matrix4d();
	transform.set(orientation);
	transform.invert();
	CalcPoint.transform(transform, cloud);

	// The orientation found now should be 0 (it has been re-oriented)
	orientation = UnitQuaternions.orientation(cloud);
	AxisAngle4d axis = new AxisAngle4d();
	axis.set(orientation);

	// No significant rotation
	assertEquals(orientation.x, 0.0, 0.01);
	assertEquals(orientation.y, 0.0, 0.01);
	assertEquals(orientation.z, 0.0, 0.01);
	assertEquals(axis.angle, 0.0, 0.01);

	// Now try to recover an orientation
	Quat4d quat = new Quat4d(0.418, 0.606, 0.303, 0.606);

	Matrix4d mat = new Matrix4d();
	mat.set(quat);

	CalcPoint.transform(mat, cloud);

	orientation = UnitQuaternions.orientation(cloud);

	// Test recovering the quaternion (q and -q same rotation)
	assertEquals(Math.abs(orientation.x), quat.x, 0.01);
	assertEquals(Math.abs(orientation.y), quat.y, 0.01);
	assertEquals(Math.abs(orientation.z), quat.z, 0.01);
	assertEquals(Math.abs(orientation.w), quat.w, 0.01);
}
 

/**
 * Calculate the rotation angle component of the input unit quaternion.
 *
 * @param q
 *            unit quaternion Quat4d
 * @return the angle in radians of the input quaternion
 */
public static double angle(Quat4d q) {
	AxisAngle4d axis = new AxisAngle4d();
	axis.set(q);
	return axis.angle;
}
 

/**
 * Returns the AxisAngle of the helix transformation
 *
 * @param transformation
 *            helix transformation
 * @return
 */
private static AxisAngle4d getAxisAngle(Matrix4d transformation) {
	AxisAngle4d axis = new AxisAngle4d();
	axis.set(transformation);
	return axis;
}