[MOAB-dev] r6027 - MOAB/trunk/tools/mbcslam

[iulian at mcs.anl.gov]

Author: iulian
Date: 2013-03-08 18:38:13 -0600 (Fri, 08 Mar 2013)
New Revision: 6027

Modified:
   MOAB/trunk/tools/mbcslam/CslamUtils.cpp
   MOAB/trunk/tools/mbcslam/case1_test.cpp
Log:
increase the robustness of l'Huiller itself
we compute angles between 2 vectors, and our usual 
angle(CartVect & u, CartVect & v) is not robust enough
inline double angle( const CartVect& u, const CartVect& v )
  { return std::acos( (u % v) / std::sqrt((u % u) * (v % v)) ); }

when u is close to v, (u % v) / std::sqrt((u % u) * (v % v)) is 
slightly larger than 1, and std::acos will give nan
which is a shame :(

this should have been taken case by stdlib (or math?) library.
maybe there are some better implementations of std::acos



Modified: MOAB/trunk/tools/mbcslam/CslamUtils.cpp
===================================================================
--- MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-03-07 19:28:15 UTC (rev 6026)
+++ MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-03-09 00:38:13 UTC (rev 6027)
@@ -697,13 +697,23 @@
  *
  *  E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
  */
+
+//! Interior angle between two vectors
+inline double angleSafe( const CartVect& u, const CartVect& v )
+  {
+    double a = ( (u % v) / std::sqrt((u % u) * (v % v)) );
+    if (a> 1.) a= 1.;
+    if (a<-1.) a =-1.;
+    return std::acos(a);
+  }
+
 double area_spherical_triangle_lHuiller(double * ptA, double * ptB, double * ptC, double Radius)
 {
   // now, a is angle BOC, O is origin
   CartVect vA(ptA), vB(ptB), vC(ptC);
-  double a = angle(vB, vC);
-  double b = angle(vC, vA);
-  double c = angle(vA, vB);
+  double a = angleSafe(vB, vC);
+  double b = angleSafe(vC, vA);
+  double c = angleSafe(vA, vB);
   int sign =1;
   if ((vA*vB)%vC < 0)
     sign = -1;
@@ -712,6 +722,8 @@
   if (tmp<0.)
     tmp = 0.;
   double E = 4*atan(sqrt(tmp));
+  if (E!=E)
+    std::cout << "bad nan here \n";
   return sign * E * Radius*Radius;
 }
 

Modified: MOAB/trunk/tools/mbcslam/case1_test.cpp
===================================================================
--- MOAB/trunk/tools/mbcslam/case1_test.cpp	2013-03-07 19:28:15 UTC (rev 6026)
+++ MOAB/trunk/tools/mbcslam/case1_test.cpp	2013-03-09 00:38:13 UTC (rev 6027)
@@ -118,6 +118,7 @@
 {
 
 
+  const char *filename_mesh1 = STRINGIFY(SRCDIR) "/eulerHomme.vtk";
   if (argc > 1)