From 0cfeda19c1b4de0c0e84f24eef4883687aa471dd Mon Sep 17 00:00:00 2001
From: "T.Alderighi" <thomasalderighi@gmail.com>
Date: Thu, 21 Nov 2019 17:14:34 +0100
Subject: [PATCH] fix bug on rotationmatrix computation corner case

---
 vcg/math/matrix33.h | 44 ++++++++++++++++++++++++++++++++------------
 1 file changed, 32 insertions(+), 12 deletions(-)

diff --git a/vcg/math/matrix33.h b/vcg/math/matrix33.h
index f74551ff..efa797cf 100644
--- a/vcg/math/matrix33.h
+++ b/vcg/math/matrix33.h
@@ -518,19 +518,39 @@ Matrix33<S> RotationMatrix(vcg::Point3<S> v0,vcg::Point3<S> v1,bool normalized=t
             rotM.SetIdentity();
             return rotM;
         }
-        if (dot<(-(S)1+epsilon))
-        {
-            rotM.SetZero();
-            rotM[0][0]=-1;
-            rotM[1][1]=-1;
-            rotM[2][2]=-1;
-            return (rotM);
-        }
-        ///find the axis of rotation
-        CoordType axis;
-        axis=v0^v1;
-        axis.Normalize();
 
+        //find the axis of rotation
+        CoordType axis;
+
+        //if dot = -1 rotating to opposite vertex
+        //the problem is underdefined, so choose axis such that division is more stable
+        //alternative solution at http://cs.brown.edu/research/pubs/pdfs/1999/Moller-1999-EBA.pdf
+        if (dot < (S)-1 + epsilon)
+        {
+            S max = std::numeric_limits<S>::min();
+            int maxInd = 0;
+            for (int i = 0; i < 3; ++i)
+            {
+                if (v0[i] > max)
+                {
+                    max    = v0[i];
+                    maxInd = i;
+                }
+            }
+
+            axis[maxInd] = - (v0[(maxInd+2) % 3] / v0[maxInd]);
+            axis[(maxInd+1) % 3] = 0;
+            axis[(maxInd+2) % 3] = 1;
+
+            dot = (S)-1;
+        } 
+        else 
+        {
+            axis=v0^v1;
+        }
+
+        axis.Normalize();
+    
         ///construct rotation matrix
         S u=axis.X();
         S v=axis.Y();