From 2789acb5c59196abd3a94d0094cea9470667a67c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Th=C3=A9ophile=20Bastian?= <contact@tobast.fr>
Date: Tue, 13 Feb 2018 13:22:41 +0100
Subject: [PATCH] Implicit surfaces: add normal vector

---
 Implicit.cpp | 13 +++++++++++++
 Implicit.hpp |  3 +++
 2 files changed, 16 insertions(+)

diff --git a/Implicit.cpp b/Implicit.cpp
index 43462df..2ab78eb 100644
--- a/Implicit.cpp
+++ b/Implicit.cpp
@@ -3,3 +3,16 @@
 double ImplicitSurface::operator()(const Point& pt) const {
     return operator()(pt.x, pt.y, pt.z);
 }
+
+Point ImplicitSurface::normal_at(const Point& pt) const {
+    // this is simply the gradient.
+    static const double d_dist = 1e-8;
+
+    return Point(
+            (operator()(pt.x + d_dist, pt.y, pt.z) -
+                operator()(pt.x - d_dist, pt.y, pt.z)) / (2. * d_dist),
+            (operator()(pt.x, pt.y + d_dist, pt.z) -
+                operator()(pt.x, pt.y - d_dist, pt.z)) / (2. * d_dist),
+            (operator()(pt.x, pt.y, pt.z + d_dist) -
+                operator()(pt.x, pt.y, pt.z - d_dist)) / (2. * d_dist));
+}
diff --git a/Implicit.hpp b/Implicit.hpp
index 113c268..dc4a3d5 100644
--- a/Implicit.hpp
+++ b/Implicit.hpp
@@ -7,6 +7,9 @@ class ImplicitSurface {
         virtual double operator() (double x, double y, double z) const = 0;
         double operator()(const Point& pt) const;
         Point getCenter() const { return center;}
+
+        /// Compute the normal vector to the isosurface at `pt`
+        Point normal_at(const Point& pt) const;
     protected:
         Point center;
         ImplicitSurface(Point _center) : center(_center) {}