geom/tri: add triangle area formula

geom/tri.cpp

@@ -11,7 +11,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2015 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2016-2018 Danny Robson <danny@nerdcruft.net>
  */
 
-
+#include "tri.hpp"

geom/tri.hpp

@@ -11,7 +11,68 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2015 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2016-2018 Danny Robson <danny@nerdcruft.net>
  */
 
+#ifndef CRUFT_GEOM_TRI_HPP
+#define CRUFT_GEOM_TRI_HPP
 
+#include <cstddef>
+
+#include "../point.hpp"
+
+namespace util::geom {
+    ///////////////////////////////////////////////////////////////////////////
+    // n-dimensional triangle area
+    template <std::size_t DimensionV, typename ValueT>
+    ValueT
+    area (util::point<DimensionV,ValueT> a,
+          util::point<DimensionV,ValueT> b,
+          util::point<DimensionV,ValueT> c)
+    {
+        // heron's formula
+        const auto ab = util::distance (a, b);
+        const auto bc = util::distance (b, c);
+        const auto ca = util::distance (c, a);
+
+        const auto s = (ab + bc + ca) / 2;
+
+        return std::sqrt (s * (s - ab) * (s - bc) * (s - ca));
+    }
+
+
+    //-------------------------------------------------------------------------
+    // 2-dimension triangle area
+    template <typename T>
+    T
+    area (util::point<2,T> a, util::point<2,T> b, util::point<2,T> c)
+    {
+        //                | x1 y1 1 |
+        // area = 0.5 det | x2 y2 1 |
+        //                | x3 y3 1 |
+
+        return std::abs (
+            -b.x * a.y
+            +c.x * a.y
+            +a.x * b.y
+            -c.x * b.y
+            -a.x * c.y
+            +b.x * c.y
+        ) / 2;
+    }
+
+
+    //-------------------------------------------------------------------------
+    // 3-dimension triangle area
+    template <typename T>
+    T
+    area (util::point<3,T> a, util::point<3,T> b, util::point<3,T> c)
+    {
+        const auto ab = a - b;
+        const auto ac = a - c;
+
+        return norm (cross (ab, ac)) / 2;
+    }
+};
+
+#endif