diff --git a/images/interpolation_nn_c1.img b/images/interpolation_nn_c1.img
index 0c16345..dcc8e10 100644
--- a/images/interpolation_nn_c1.img
+++ b/images/interpolation_nn_c1.img
@@ -1 +1 @@
-
\ No newline at end of file
+
\ No newline at end of file
diff --git a/images/interpolation_nn_c1.jpg b/images/interpolation_nn_c1.jpg
index bb196da..8a4591b 100644
Binary files a/images/interpolation_nn_c1.jpg and b/images/interpolation_nn_c1.jpg differ
diff --git a/src/delaunay_core/interpolation.rs b/src/delaunay_core/interpolation.rs
index 4a9829b..4c699c9 100644
--- a/src/delaunay_core/interpolation.rs
+++ b/src/delaunay_core/interpolation.rs
@@ -314,9 +314,11 @@ where
/// The value to interpolate is given by the `i` parameter. The gradient that defines the derivative at
/// each input vertex is given by the `g` parameter.
///
- /// The `flatness` parameter blends between an interpolation that ignores the given gradients (value 0.0)
- /// or adheres to it strongly (values larger than ~2.0) in the vicinity of any vertex. When in doubt, using
- /// a value of 1.0 should result in a good interpolation and is also the fastest.
+ /// The `flatness` parameter (should be >= 0.0) blends between an interpolation that adheres less to the gradient (small values)
+ /// or adheres to it strongly (values larger than ~2.0) in the vicinity of any vertex. A value of 0.5 reproduces
+ /// the original Sibson C1 interpolant. A value of 1.0 is the fastest and works well in many situations.
+ /// A value of 0.0 reproduces the I1 interpolant introduced in Flötotto's thesis chapter 4.4.
+ /// If both the flatness and the gradients are 0.0 Sibsons C0 interpolant is recreated and [Self::interpolate] should be used.
///
/// Returns `None` for any point outside the triangulation's convex hull.
///
@@ -358,6 +360,10 @@ where
///
/// This method uses the C1 extension proposed by Sibson in
/// "A brief description of natural neighbor interpolation, R. Sibson, 1981"
+ ///
+ /// It is also worth looking at Julia Flötotto's thesis
+ /// A Coordinate System associated to a Point Cloud issued from a Manifold: Definition, Properties and Applications
+ /// chapter 4.1 and 4.4
pub fn interpolate_gradient(
&self,
i: I,
@@ -390,14 +396,15 @@ where
let handle = self.triangulation.vertex(*handle);
let pos_i = handle.position();
let h_i = i(handle);
- let diff = pos_i.sub(position);
+ let diff = position.sub(pos_i);
let r_i2 = diff.length2();
let r_i = r_i2.powf(flatness);
let c1_weight_i = *weight / r_i;
let grad_i = g(handle);
let zeta_i = h_i + diff.dot(grad_i.into());
- alpha = alpha + c1_weight_i * r_i;
- beta = beta + c1_weight_i * r_i2;
+
+ alpha = alpha + c1_weight_i * r_i2;
+ beta = beta + *weight * r_i2;
sum_c1_weights = sum_c1_weights + c1_weight_i;
sum_c1 = sum_c1 + zeta_i * c1_weight_i;
sum_c0 = sum_c0 + h_i * *weight;