Convert to MHL's Euler definitions (-x, -z, -x).

git-svn-id: svn+ssh://git-open/scm/biosvn/bug2drug@6109 1ff99547-d4dc-0310-9227-d4293942b958

ExpandSet.cpp

@@ -331,47 +331,37 @@ void Quaternion::PrintEuler(ostream& s) const {
   // Convert to rotation matrix (assume quaternion is already
   // normalized)
   double
-	//	m00 = 1 - 2*y*y - 2*z*z,
-	m01 = 	  2*x*y - 2*z*w,
-	m02 = 	  2*x*z + 2*y*w,
-	//	m10 = 	  2*x*y + 2*z*w,
-	m11 = 1 - 2*x*x - 2*z*z,
-	m12 = 	  2*y*z - 2*x*w,
-	m20 = 	  2*x*z - 2*y*w,
-	m21 = 	  2*y*z + 2*x*w,
-	m22 = 1 - 2*x*x - 2*y*y;
+    m00 = 1 - 2*y*y - 2*z*z,
+    m01 =     2*x*y - 2*z*w,
+    m02 =     2*x*z + 2*y*w,
+    m10 =     2*x*y + 2*z*w,
+    // m11 = 1 - 2*x*x - 2*z*z,
+    m12 =     2*y*z - 2*x*w,
+    m20 =     2*x*z - 2*y*w,
+    // m21 =     2*y*z + 2*x*w,
+    m22 = 1 - 2*x*x - 2*y*y;
   // Taken from Ken Shoemake, "Euler Angle Conversion", Graphics Gems
   // IV, Academic 1994.
   //
   //    http://vered.rose.utoronto.ca/people/david_dir/GEMS/GEMS.html
-  double sy = sqrt(m02*m02 + m12*m12);
+  double sy = sqrt(m10*m10 + m20*m20);
   double a,  b, c;
-  b = atan2(sy, m22);
+  b = atan2(sy, m00);
   if (sy > 16 * numeric_limits<double>::epsilon()) {
-	a = atan2(m21, m20);
-	c = atan2(m12, -m02);
+	a = atan2(m02, m01);
+	c = atan2(m20, -m10);
   } else {
 	a = 0;
-	c = atan2(m01, m11);
+	c = atan2(m12, m22);
   }
-  // b is already in [0, pi]. Fold a, c to [0, 2*pi].
-  if (a < 0)
-	a += 2*M_PI;
-  if (c < 0)
-	c += 2*M_PI;
-  // Convert -0 to 0.
-  if (a < numeric_limits<double>::min())
-	a = 0;
-  if (c < numeric_limits<double>::min())
-	c = 0;
-  s << fixed << setprecision(9) << setw(11) << a << " "
-    << setw(11) << b << " " << setw(11) << c;
+  s << fixed << setprecision(9) << setw(12) << a << " "
+    << setw(12) << b << " " << setw(12) << c;
 
 #if !defined(NDEBUG)
   // Sanity check.  Convert from Euler angles back to a quaternion, q
-  Quaternion q = Quaternion(cos(c/2), 0, 0, -sin(c/2)). // -c about z
-    Times(Quaternion(cos(b/2), 0, -sin(b/2), 0). // -b about y
-	  Times(Quaternion(cos(a/2), 0, 0, -sin(a/2)))); // -a about z
+  Quaternion q = Quaternion(cos(c/2), -sin(c/2), 0, 0). // -c about x
+    Times(Quaternion(cos(b/2), 0, 0, -sin(b/2)). // -b about z
+	  Times(Quaternion(cos(a/2), -sin(a/2), 0, 0))); // -a about x
   // and check that q is parallel to *this.
   double t = abs(q.w * w + q.x * x + q.y * y + q.z * z);
   assert(t > 1 - 16 * numeric_limits<double>::epsilon());

makefile → Makefile.distrib

@@ -1,6 +1,8 @@
 # $Id$
 distrib:
 	cd data && make
+
+publish: distrib
 	rsync -av index.html ExpandSet.cpp /home/ckarney/web/orientation/
-	rsync -av data/orient.tgz data/*.quat data/*.grid data/*.euler /home/ckarney/web/orientation/data/
+	rsync -av data/orient*.tgz data/*.quat data/*.grid data/*.euler /home/ckarney/web/orientation/data/
 	rsync -av --delete /home/ckarney/web/orientation/ petrel.petrel.org:web/orientation/

PackSet.cpp

@@ -404,7 +404,7 @@ double PackSet::MinMaxRadius(size_t k, double eps) {
 }
 
 double PackSet::Volume(size_t ind, double radius, size_t n,
-		       double& maxrad, const Quaternion& tx) {
+		       double& maxrad, const Quaternion& tx, bool docomp) {
   PackSet temp;
   temp.Add(Quaternion(1,0,0,0));
   Quaternion q = Member(ind);
@@ -415,6 +415,28 @@ double PackSet::Volume(size_t ind, double radius, size_t n,
     if (Dist(Closeness(q, Member(i))) <= 2 * radius)
       temp.Add(Member(i) * qc, 0, true);
   }
+  if (!docomp) {
+    cout << "x=[" << setprecision(13);
+    for (size_t i = 0;;) {
+      bool turnp = false;
+      Vector3D<double> t =
+	tx.transformPoint(temp.Member(i).RotateVector(turnp));
+      cout << t.x << " " << t.y << " " << t.z;
+      if (++i == temp.Number()) {
+	cout << "];" << endl;
+	break;
+      } else
+	cout << ";";
+    }
+    cout << "[v, c] = voronoin(x);" << endl;
+    cout << "y = v(c{1},:);" << endl;
+    cout << "[k, volume] = convhulln(y);" << endl;
+    cout << "rad = max(sqrt(sum(v(c{1},:).^2,2)));" << endl;
+    cout << "i = i+1;" << endl;
+    cout << "vol(i) = volume;" << endl;
+    cout << "maxrad(i) = rad;" << endl;
+    return 0;
+  }
   // Integrate over sphere of radius radt by constructing cubical grid
   // and count cubes.  V = int WithinCell(x) dV.  Better to integrate
   // over the surface of the unit sphere and compute V = int l(Omega)^3

PackSet.h

@@ -36,7 +36,9 @@ class PackSet {
   size_t MonteCarlo(size_t num, double delta, double beta);
   size_t MonteCarloA(size_t num, double delta, double beta);
   double Volume(size_t ind, double radius, size_t n,
-		double& maxrad, const Quaternion& tx = Quaternion(1,0,0,0));
+		double& maxrad,
+		const Quaternion& tx = Quaternion(1,0,0,0),
+		bool docomp = true);
 private:
   vector<Quaternion> m_set;
   vector<unsigned char> m_cat;