diff --git a/src/interpolation/mod.rs b/src/interpolation/mod.rs
index b748cdd..b5a4c41 100755
--- a/src/interpolation/mod.rs
+++ b/src/interpolation/mod.rs
@@ -97,7 +97,7 @@ impl Interpolant for Laplace {
             let loc = dt.locate(p[0], p[1]);
             match loc {
                 Ok(_tr) => {
-                    match dt.insert_one_pt(p[0], p[1], 0.) {
+                    match dt.insert_one_pt_interpol(p[0], p[1]) {
                         Ok(pi) => {
                             //-- no extrapolation
                             if dt.is_vertex_convex_hull(pi) {
@@ -134,7 +134,9 @@ impl Interpolant for Laplace {
                                 re.push(Ok(z / sumweights));
                             }
                         }
-                        Err(e) => re.push(Ok(dt.stars[e.0].pt[2])),
+                        Err((pi, _updated)) => {
+                            re.push(Ok(dt.stars[pi].pt[2]));
+                        }
                     }
                 }
                 Err(_e) => re.push(Err(StartinError::OutsideConvexHull)),
@@ -243,7 +245,7 @@ impl Interpolant for NNI {
             let loc = dt.locate(p[0], p[1]);
             match loc {
                 Ok(_tr) => {
-                    match dt.insert_one_pt(p[0], p[1], 0.) {
+                    match dt.insert_one_pt_interpol(p[0], p[1]) {
                         Ok(pi) => {
                             //-- no extrapolation
                             if dt.is_vertex_convex_hull(pi) {
diff --git a/src/lib.rs b/src/lib.rs
index 4c99c88..3433217 100755
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -582,6 +582,48 @@ impl Triangulation {
         Ok(pi)
     }
 
+    pub fn insert_one_pt_interpol(&mut self, px: f64, py: f64) -> Result<usize, (usize, bool)> {
+        let pz = 0.0;
+        if !self.is_init {
+            return self.insert_one_pt_init_phase(px, py, pz);
+        }
+        //-- walk
+        let p: [f64; 3] = [px, py, pz];
+        let tr = self.walk(&p);
+        if geom::distance2d_squared(&self.stars[tr.v[0]].pt, &p) <= (self.snaptol * self.snaptol) {
+            return Err((tr.v[0], false));
+        }
+        if geom::distance2d_squared(&self.stars[tr.v[1]].pt, &p) <= (self.snaptol * self.snaptol) {
+            return Err((tr.v[1], false));
+        }
+        if geom::distance2d_squared(&self.stars[tr.v[2]].pt, &p) <= (self.snaptol * self.snaptol) {
+            return Err((tr.v[2], false));
+        }
+        //-- ok we now insert the point in the data structure
+        let pi: usize;
+        if self.removed_indices.is_empty() {
+            self.stars.push(Star::new(px, py, pz));
+            pi = self.stars.len() - 1;
+        } else {
+            // self.stars.push(Star::new(px, py, pz));
+            pi = self.removed_indices.pop().unwrap();
+            self.stars[pi].pt[0] = px;
+            self.stars[pi].pt[1] = py;
+            self.stars[pi].pt[2] = pz;
+        }
+        //-- flip13()
+        self.flip13(pi, &tr);
+        //-- update_dt()
+        self.update_dt(pi);
+        self.cur = pi;
+        //-- extra attributes
+        match &mut self.attributes {
+            Some(x) => x.push(json!({})),
+            _ => (),
+        }
+        Ok(pi)
+    }
+
     fn update_z_value_duplicate(&mut self, vi: usize, newz: f64) -> bool {
         let mut re = false;
         match self.duplicates_handling {
diff --git a/tests/interpolate.rs b/tests/interpolate.rs
index 1d5c1cf..649ff55 100755
--- a/tests/interpolate.rs
+++ b/tests/interpolate.rs
@@ -128,6 +128,22 @@ fn existing_point() {
     assert_eq!(Ok(11.1), interpolate(&i_idw, &mut dt, &vec![[5.0, 5.0]])[0]);
 }
 
+#[test]
+fn existing_point_highest() {
+    let mut pts: Vec<[f64; 3]> = Vec::new();
+    pts.push([0.0, 0.0, 1.0]);
+    pts.push([10.0, 0.0, 2.0]);
+    pts.push([10.0, 10.0, 3.0]);
+    pts.push([0.0, 10.0, 4.0]);
+    let mut dt = startin::Triangulation::new();
+    dt.set_duplicates_handling(startin::DuplicateHandling::Lowest);
+    dt.insert(&pts, startin::InsertionStrategy::AsIs);
+    let _re = dt.insert_one_pt(5.0, 5.0, 11.1);
+
+    let i_lap = startin::interpolation::Laplace {};
+    assert_eq!(Ok(11.1), interpolate(&i_lap, &mut dt, &vec![[5.0, 5.0]])[0]);
+}
+
 #[test]
 fn middle() {
     let mut dt = four_points();