geom/ellipse: add sample_surface function

geom/ellipse.cpp

@@ -11,7 +11,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2015-2017 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2015-2018 Danny Robson <danny@nerdcruft.net>
  */
 
 #include "ellipse.hpp"
@@ -46,9 +46,6 @@ util::geom::intersects (A<S,T> a, B<S,T> b)
 }
 
 
-//-----------------------------------------------------------------------------
-
-
 ///////////////////////////////////////////////////////////////////////////////
 template <size_t DimensionV, typename ValueT>
 util::point<DimensionV,ValueT>
@@ -103,7 +100,6 @@ util::geom::bounds (K<S,T> k)
 }
 
 
-
 ///////////////////////////////////////////////////////////////////////////////
 #define INSTANTIATE_S_T(S,T) \
 template bool util::geom::intersects (ellipse<S,T>, util::point<S,T>); \

geom/ellipse.hpp

@@ -11,19 +11,21 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  *
- * Copyright 2015 Danny Robson <danny@nerdcruft.net>
+ * Copyright 2015-2018 Danny Robson <danny@nerdcruft.net>
  */
 
 #ifndef __UTIL_GEOM_ELLIPSE_HPP
 #define __UTIL_GEOM_ELLIPSE_HPP
 
-#include <cstdlib>
-
 #include "fwd.hpp"
 
 #include "../point.hpp"
 #include "../vector.hpp"
 
+#include <cstdlib>
+#include <random>
+
+
 namespace util::geom {
     ///////////////////////////////////////////////////////////////////////////
     template <size_t S, typename ValueT>
@@ -38,6 +40,40 @@ namespace util::geom {
         util::vector<S,ValueT> up;
     };
 
+    using ellipse3f = ellipse<3,float>;
+
+
+    /// returns the approximate surface area of the ellipse
+    ///
+    /// the general form involves _substantially_ more expensive and
+    /// complicated maths which is prohibitive right now.
+    ///
+    /// the relative error should be at most 1.061%
+    inline float
+    area (ellipse3f self)
+    {
+        auto const semiprod = self.radius * self.radius.indices<1,2,0> ();
+        auto const semipow  = pow (semiprod, 1.6f);
+
+        return 4 * pi<float> * std::pow (sum (semipow) / 3, 1/1.6f);
+    }
+
+
+    template <typename T>
+    T
+    area (ellipse<2,T> self)
+    {
+        return pi<T> * product (self.radius);
+    }
+
+
+    template <size_t S, typename T>
+    T
+    volume (ellipse<S,T> self)
+    {
+        return 4 / T{3} * pi<T> * product (self.radius);
+    }
+
 
     template <size_t DimensionV, typename ValueT>
     point<DimensionV,ValueT>
@@ -57,7 +93,35 @@ namespace util::geom {
         ellipse<DimensionV,ValueT>
     );
 
-    using ellipse3f = ellipse<3,float>;
+
+    /// returns a point that is uniformly distributed about the ellipse
+    /// surface.
+    ///
+    /// NOTE: I don't have a strong proof that the below is in fact properly
+    /// uniformly distributed, so if you need a strong guarantee for your work
+    /// then it might not be the best option. But visual inspection appears to
+    /// confirm there aren't obvious patterns.
+    ///
+    /// the concept was taken from: https://math.stackexchange.com/a/2514182
+    template <typename RandomT>
+    util::point3f
+    sample_surface (ellipse3f self, RandomT &generator)
+    {
+        const auto [a, b, c] = self.radius;
+        const auto a2 = a * a;
+        const auto b2 = b * b;
+        const auto c2 = c * c;
+
+        // generate a direction vector from a normally distributed random variable
+        auto const x = std::normal_distribution<float> (0, a2) (generator);
+        auto const y = std::normal_distribution<float> (0, b2) (generator);
+        auto const z = std::normal_distribution<float> (0, c2) (generator);
+
+        // find the distance to the surface along the direction vector
+        auto const d = std::sqrt (x * x / a2 + y * y / b2 + z * z / c2);
+
+        return self.origin + util::vector3f {x,y,z} / d;
+    }
 }