diff --git a/CDT/include/CDTUtils.h b/CDT/include/CDTUtils.h
index dda19c4..f7a4cde 100644
--- a/CDT/include/CDTUtils.h
+++ b/CDT/include/CDTUtils.h
@@ -446,6 +446,23 @@ CDT_EXPORT bool isEncroachingOnEdge(
 template <typename T>
 CDT_EXPORT V2d<T> circumcenter(V2d<T> a, V2d<T> b, const V2d<T>& c);
 
+/// Doubled surface area of a triangle ABC
+template <typename T>
+CDT_EXPORT T doubledArea(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
+
+/// Surface area of a triangle ABC
+template <typename T>
+CDT_EXPORT T area(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
+
+/// Sine of smallest angle of triangle ABC
+template <typename T>
+CDT_EXPORT T
+sineOfSmallestAngle(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
+
+/// Smallest angle of triangle ABC in radians
+template <typename T>
+CDT_EXPORT T smallestAngle(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c);
+
 } // namespace CDT
 
 #ifndef CDT_USE_AS_COMPILED_LIBRARY
diff --git a/CDT/include/CDTUtils.hpp b/CDT/include/CDTUtils.hpp
index 62605ad..93ed4f9 100644
--- a/CDT/include/CDTUtils.hpp
+++ b/CDT/include/CDTUtils.hpp
@@ -331,4 +331,35 @@ V2d<T> circumcenter(V2d<T> a, V2d<T> b, const V2d<T>& c)
         c.y + (a.x * bLenSq - b.x * aLenSq) * denom);
 }
 
+template <typename T>
+T doubledArea(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
+{
+    return std::abs(orient2D(a, b, c));
+}
+
+template <typename T>
+T area(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
+{
+    return doubledArea(a, b, c) / T(2);
+}
+
+template <typename T>
+T sineOfSmallestAngle(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
+{
+    // find sides of the smallest angle using law of sines:
+    T sideA = distance(a, b), sideB = distance(b, c);
+    if(sideA > sideB)
+        std::swap(sideA, sideB);
+    sideA = std::max(sideA, distance(a, c));
+    return (doubledArea(a, b, c) / sideA) / sideB;
+}
+
+template <typename T>
+T smallestAngle(const V2d<T>& a, const V2d<T>& b, const V2d<T>& c)
+{
+    const T angleSine = sineOfSmallestAngle(a, b, c);
+    assert(angleSine >= -1 && angleSine <= 1);
+    return std::asin(angleSine);
+}
+
 } // namespace CDT
diff --git a/CDT/include/Triangulation.h b/CDT/include/Triangulation.h
index 119ae6c..473f08b 100644
--- a/CDT/include/Triangulation.h
+++ b/CDT/include/Triangulation.h
@@ -86,7 +86,7 @@ struct CDT_EXPORT IntersectingConstraintEdges
 /**
  * Enum of strategies for triangles refinement
  */
-struct CDT_EXPORT RefineTriangles
+struct CDT_EXPORT RefinementCriterion
 {
     /**
      * The Enum itself
@@ -94,8 +94,8 @@ struct CDT_EXPORT RefineTriangles
      */
     enum Enum
     {
-        ByAngle, ///< constraint minimum triangles angle
-        ByArea,  ///< constraint maximum triangles area
+        SmallestAngle, ///< constraint minimum triangles angle
+        LargestArea,   ///< constraint maximum triangles area
     };
 };
 
@@ -305,7 +305,8 @@ class CDT_EXPORT Triangulation
      * @param threshold threshold value for refinement
      */
     void refineTriangles(
-        RefineTriangles::Enum refinementConstrain = RefineTriangles::ByAngle,
+        RefinementCriterion::Enum refinementConstrain =
+            RefinementCriterion::SmallestAngle,
         T threshold = 20 / 180.0 * M_PI);
     /**
      * Erase triangles adjacent to super triangle
@@ -521,10 +522,10 @@ class CDT_EXPORT Triangulation
         VertInd iV3,
         VertInd iV4) const;
     TriInd edgeTriangle(Edge edge) const;
-    bool isBadTriangle(
+    bool isRefinementNeeded(
         const Triangle& tri,
-        RefineTriangles::Enum refinement,
-        T threshold) const;
+        RefinementCriterion::Enum refinementCriterion,
+        T refinementThreshold) const;
     /// Search in all fixed edges to find encroached edges, each fixed edge is
     /// checked against its opposite vertices
     /// Returns queue of encroached edges
@@ -540,7 +541,7 @@ class CDT_EXPORT Triangulation
         const V2d<T>& v = {},
         bool validV = false,
         bool fillBadTriangles = false,
-        RefineTriangles::Enum refinementConstrain = {},
+        RefinementCriterion::Enum refinementConstrain = {},
         T badTriangleThreshold = {});
     VertInd splitEncroachedEdge(Edge e, TriInd iT, TriInd iTopo);
     void changeNeighbor(TriInd iT, TriInd oldNeighbor, TriInd newNeighbor);
diff --git a/CDT/include/Triangulation.hpp b/CDT/include/Triangulation.hpp
index 9b05d46..19634bf 100644
--- a/CDT/include/Triangulation.hpp
+++ b/CDT/include/Triangulation.hpp
@@ -1286,45 +1286,23 @@ TriInd Triangulation<T, TNearPointLocator>::edgeTriangle(const Edge edge) const
 }
 
 template <typename T, typename TNearPointLocator>
-bool Triangulation<T, TNearPointLocator>::isBadTriangle(
+bool Triangulation<T, TNearPointLocator>::isRefinementNeeded(
     const Triangle& tri,
-    const RefineTriangles::Enum refinement,
-    const T threshold) const
-{
-    const V2d<T>& v1 = vertices[tri.vertices[0]];
-    const V2d<T>& v2 = vertices[tri.vertices[1]];
-    const V2d<T>& v3 = vertices[tri.vertices[2]];
-
-    const T twiceArea = std::abs(orient2D(v1, v2, v3));
-    bool ans = false;
-    switch(refinement)
-    {
-    case RefineTriangles::ByAngle: {
-        T opoLenV1 = distance(v2, v3);
-        T opoLenV2 = distance(v1, v3);
-        T opoLenV3 = distance(v1, v2);
-        // Let opoLenV1 is the smallest edge length
-        if((opoLenV2 < opoLenV1) && (opoLenV2 < opoLenV3))
-        {
-            std::swap(opoLenV1, opoLenV2);
-        }
-        else if(opoLenV3 < opoLenV1)
-        {
-            std::swap(opoLenV1, opoLenV3);
-        }
-        assert(opoLenV1 <= opoLenV3);
-        assert(opoLenV1 <= opoLenV2);
-        T samllestAngle = twiceArea / opoLenV3 / opoLenV2;
-        ans = samllestAngle < sin(threshold);
-        break;
-    }
-    case RefineTriangles::ByArea: {
-        T area = 0.5 * twiceArea;
-        ans = area > threshold;
-        break;
-    }
+    RefinementCriterion::Enum refinementCriterion,
+    const T refinementThreshold) const
+{
+    const V2d<T>& a = vertices[tri.vertices[0]];
+    const V2d<T>& b = vertices[tri.vertices[1]];
+    const V2d<T>& c = vertices[tri.vertices[2]];
+    switch(refinementCriterion)
+    {
+    case RefinementCriterion::SmallestAngle:
+        return smallestAngle(a, b, c) <= refinementThreshold;
+    case RefinementCriterion::LargestArea:
+        return area(a, b, c) >= refinementThreshold;
     }
-    return ans;
+    assert(false); // unreachable code
+    return false;
 }
 
 /// Search in all fixed edges to find encroached edges, each fixed edge is
@@ -1384,7 +1362,7 @@ Triangulation<T, TNearPointLocator>::resolveEncroachedEdges(
     const V2d<T>& v,
     const bool validV,
     const bool fillBadTriangles,
-    const RefineTriangles::Enum refinementConstrain,
+    const RefinementCriterion::Enum refinementConstrain,
     const T badTriangleThreshold)
 {
     IndexSizeType numOfSplits = 0;
@@ -1410,7 +1388,7 @@ Triangulation<T, TNearPointLocator>::resolveEncroachedEdges(
         {
             const Triangle& t = triangles[currTri];
             if(fillBadTriangles &&
-               isBadTriangle(t, refinementConstrain, badTriangleThreshold))
+               isRefinementNeeded(t, refinementConstrain, badTriangleThreshold))
             {
                 badTriangles.push_back(currTri);
             }
@@ -2297,7 +2275,7 @@ void Triangulation<T, TNearPointLocator>::tryInitNearestPointLocator()
 
 template <typename T, typename TNearPointLocator>
 void Triangulation<T, TNearPointLocator>::refineTriangles(
-    RefineTriangles::Enum refinementConstrain,
+    RefinementCriterion::Enum refinementConstrain,
     T threshold)
 {
     if(isFinalized())
@@ -2315,7 +2293,7 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
         if(t.vertices[0] < 3 || t.vertices[1] < 3 || t.vertices[2] < 3)
             continue;
 
-        if(isBadTriangle(t, refinementConstrain, threshold))
+        if(isRefinementNeeded(t, refinementConstrain, threshold))
         {
             const V2d<T> vert = circumcenter(
                 vertices[t.vertices[0]],
@@ -2335,7 +2313,7 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
         TriInd iT = badTriangles.front();
         const Triangle& t = triangles[iT];
         badTriangles.pop();
-        if(!isBadTriangle(t, refinementConstrain, threshold) ||
+        if(!isRefinementNeeded(t, refinementConstrain, threshold) ||
            numOfSteinerPoints >= maxSteinerPoints)
         {
             continue;
@@ -2373,7 +2351,7 @@ void Triangulation<T, TNearPointLocator>::refineTriangles(
             do
             {
                 const Triangle& t = triangles[currTri];
-                if(isBadTriangle(t, refinementConstrain, threshold))
+                if(isRefinementNeeded(t, refinementConstrain, threshold))
                 {
                     badTriangles.push(currTri);
                 }
diff --git a/visualizer/main.cpp b/visualizer/main.cpp
index fe74a0b..10fd21e 100644
--- a/visualizer/main.cpp
+++ b/visualizer/main.cpp
@@ -224,7 +224,7 @@ public slots:
             if(m_isRemoveOuterAndHoles)
             {
                 m_cdt.refineTriangles(
-                    CDT::RefineTriangles::ByAngle, 20 / 180.0 * M_PI);
+                    CDT::RefinementCriterion::SmallestAngle, 20 / 180.0 * M_PI);
                 m_cdt.eraseOuterTrianglesAndHoles();
             }
             else if(m_isRemoveOuter)
@@ -232,7 +232,7 @@ public slots:
             else if(m_isHideSuperTri)
             {
                 m_cdt.refineTriangles(
-                    CDT::RefineTriangles::ByAngle, 20 / 180.0 * M_PI);
+                    CDT::RefinementCriterion::SmallestAngle, 20 / 180.0 * M_PI);
                 m_cdt.eraseSuperTriangle();
             }
             const CDT::unordered_map<Edge, CDT::EdgeVec> tmp =