<p>You may be able to eliminate the triangular artifacts by using Natural Neighbor Interpolation.  You can find a description of the technique and an associated description of the algorithm at <a href="https://gwlucastrig.github.io/TinfourDocs/NaturalNeighborIntro/index.html" rel="nofollow noreferrer">An Introduction to Natural Neighbor Interpolation</a>.   I've got an open-source implementation in Java.  And while it won't solve your problem (since you're working in C/C++), you could run your data through the Tinfour Viewer demonstration application to see how it looks under Natural Neighbor Interpolation. That would, at least, tell you if it was worth further consideration.  If you're interested, you can find more information at the <a href="https://github.com/gwlucastrig/Tinfour" rel="nofollow noreferrer">Tinfour Project</a>.</p>
<p>It looks like your input data is from a Lidar source.  So here's a picture of some Lidar data taken over Bear Mountain in Salisbury, CT, USA. Interpolation was performed using NNI. NNI works fine when your data represents a coherent surface.  For Lidar, that means sticking to all ground points or all first returns.  If you start mixing in things like tree branches and other vegetation, you'd be better served by something like an Inverse Distance Weighting solution.</p>
<p><a href="https://i.sstatic.net/sghlz.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/sghlz.jpg" alt="Lidar sample over Bear Mountain" /></a></p>

CC BY-SA 4.0
2022-04-27T12:48:27.407Z
If you are familiar with the code in the linked article could you link the specific part of the code that computes the weights based on the Delaunay triangulation? As I understood it it is supposed to compute the overlapping areas if a point were to be introduced, though I am not certain on the details of that. They mention the Bowyer-Watson tessellation algorithm, and as I understand it, one would have to recompute the Delaunay triangulation change (in terms of 1-ring neighbours) for each pixel.
lightxbulb
CC BY-SA 4.0
2022-04-27T14:40:12.303Z
Sure.  If you click on the link marked "Implementation Notes" at the top of the article I cited, you'll find an article describing the algorithm in detail.  It also explains how it is possible to perform a Natural Neighbor Interpolation without modifying the Delaunay triangulation (an important part of making the thing run fast).  You can find it at https://gwlucastrig.github.io/TinfourDocs/NaturalNeighborTinfourAlgorithm/index.html .  Also, if you download the code from Github, there is a Java class called NaturalNeighborInterpolator that implements the algorithm.
Gary Lucas
CC BY-SA 4.0
2022-04-27T14:49:35.363Z
Also, you may find it useful to look at Liang and Hale's algorithm.  Although the algorithm I use is largely my own, their earlier work used a lot of the same ideas I applied. You may find their approach more to your liking (in either case, I will be interested in hearing more about what you do). See https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.174.2191&rep=rep1&type=pdf
Gary Lucas
CC BY-SA 4.0
2022-04-27T20:24:14.583Z
Thank you. I have found this too: http://alexbeutel.com/papers/tsas-terranni.pdf but it seems fairly involved. What I do is sparse data interpolation, and I have been using some PDEs for that, but I realize that generalized barycentric coordinates based approaches can be faster, and have other properties. I was interested in case I ever get to implementing it, but it's currently not high up on my priority queue.
lightxbulb
CC BY-SA 4.0
2022-04-27T12:05:40.493Z
2022-04-27T12:05:40.493Z
Gary Lucas
<p>I am looking for an algorithm which would smoothly interpolate triangles of a mesh (computed by Delaunay triangulation) where each vertex has some value (elevation in my case). I need it for <a href="https://pdal.io" rel="nofollow noreferrer">PDAL</a> where I want to try to implement it. Its current implementation uses linear interpolation which makes the result not so pleasant:</p>
<p><a href="https://i.sstatic.net/m9T2P.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/m9T2P.png" alt="enter image description here" /></a></p>
<p>(image is already post-processed to display the result as hill-shading)</p>
<p>Is there some public algorithm which would interpolate it smoothly, without sharp edges? Result should be a raster (rectangular grid of points).</p>

CC BY-SA 4.0
2022-04-11T18:55:55.597Z
https://en.m.wikipedia.org/wiki/Phong_shading
lightxbulb
CC BY-SA 4.0
2022-04-12T05:57:40.840Z
@lightxbulb I actually need just to interpolate elevations, not to directly produce visual shading and playing with light direction
Martin Ždila
CC BY-SA 4.0
2022-04-12T08:36:26.167Z
You can try this: https://en.m.wikipedia.org/wiki/Point-normal_triangle
Also if the result is supposed to be on a raster grid and you have sparse data, is there a reason you compute the Delaunay triangulation instead of applying interpolation directly on the sparse data?
lightxbulb
CC BY-SA 4.0
2022-04-12T08:58:26.607Z
@lightxbulb thanks! I will take a look on PN triangle. And actually you are right - I have sparse points and want them to produce a smooth surface. Can you recommend me a method which could do it?
Martin Ždila
CC BY-SA 4.0
2022-04-12T09:03:37.970Z
I am assuming those are irregularly distributed? If that is the case you can solve the harmonic or biharmonic equation with Dirichlet boundary conditions defined by the point values.  I admit that this may be harder though if you lack the mathematical background. You can also look into  literature on radial basis functions for interpolation.
lightxbulb
CC BY-SA 4.0
2022-04-12T09:12:03.937Z
@lightxbulb Yes, irregulary distributed. Also there are billions of points (LiDAR data) so I think the computation would require supercomputers ;-). I have no math background and planned to give it not much effort. Meanwhile I found https://github.com/wangbolun300/Scattered_Points_Interpolation which could be a good solution.
Martin Ždila
CC BY-SA 4.0
2022-04-12T10:48:02.177Z
How do you even store the image if the sparse data is billions of points? Or do you generate the interpolation only for a specific region? Also Delaunay on billion of points should also be very very slow.
lightxbulb
<mesh><interpolation><triangulation>
Smooth triangular mesh interpolation
Martin Ždila
# Smooth triangular mesh interpolation

(Asked by: username_0 on 2022-04-11)

I am looking for an algorithm which would smoothly interpolate triangles of a mesh (computed by Delaunay triangulation) where each vertex has some value (elevation in my case). I need it for [PDAL](https://pdal.io) where I want to try to implement it. Its current implementation uses linear interpolation which makes the result not so pleasant:

[![enter image description here](https://i.sstatic.net/m9T2P.png)](https://i.sstatic.net/m9T2P.png)

(image is already post-processed to display the result as hill-shading)

Is there some public algorithm which would interpolate it smoothly, without sharp edges? Result should be a raster (rectangular grid of points).

**Tags:** mesh, interpolation, triangulation

> **Comments:**
> **username_1**: https://en.m.wikipedia.org/wiki/Phong_shading  
> **username_0**: @username_1 I actually need just to interpolate elevations, not to directly produce visual shading and playing with light direction  
> **username_1**: You can try this: https://en.m.wikipedia.org/wiki/Point-normal_triangle Also if the result is supposed to be on a raster grid and you have sparse data, is there a reason you compute the Delaunay triangulation instead of applying interpolation directly on the sparse data?  
> **username_0**: @username_1 thanks! I will take a look on PN triangle. And actually you are right - I have sparse points and want them to produce a smooth surface. Can you recommend me a method which could do it?  
> **username_1**: I am assuming those are irregularly distributed? If that is the case you can solve the harmonic or biharmonic equation with Dirichlet boundary conditions defined by the point values.  I admit that this may be harder though if you lack the mathematical background. You can also look into  literature on radial basis functions for interpolation.  
> **username_0**: @username_1 Yes, irregulary distributed. Also there are billions of points (LiDAR data) so I think the computation would require supercomputers ;-). I have no math background and planned to give it not much effort. Meanwhile I found https://github.com/wangbolun300/Scattered_Points_Interpolation which could be a good solution.  
> **username_1**: How do you even store the image if the sparse data is billions of points? Or do you generate the interpolation only for a specific region? Also Delaunay on billion of points should also be very very slow.  

---

## 1 answers

### Answer by username_2 on 2022-04-27. Score: 3 | Accepted answer

You may be able to eliminate the triangular artifacts by using Natural Neighbor Interpolation. You can find a description of the technique and an associated description of the algorithm at [An Introduction to Natural Neighbor Interpolation](https://gwlucastrig.github.io/TinfourDocs/NaturalNeighborIntro/index.html). I've got an open-source implementation in Java. And while it won't solve your problem (since you're working in C/C++), you could run your data through the Tinfour Viewer demonstration application to see how it looks under Natural Neighbor Interpolation. That would, at least, tell you if it was worth further consideration. If you're interested, you can find more information at the [Tinfour Project](https://github.com/gwlucastrig/Tinfour).

It looks like your input data is from a Lidar source. So here's a picture of some Lidar data taken over Bear Mountain in Salisbury, CT, USA. Interpolation was performed using NNI. NNI works fine when your data represents a coherent surface. For Lidar, that means sticking to all ground points or all first returns. If you start mixing in things like tree branches and other vegetation, you'd be better served by something like an Inverse Distance Weighting solution.

[![Lidar sample over Bear Mountain](https://i.sstatic.net/sghlz.jpg)](https://i.sstatic.net/sghlz.jpg)

> **Comments:**
> **username_1**: If you are familiar with the code in the linked article could you link the specific part of the code that computes the weights based on the Delaunay triangulation? As I understood it it is supposed to compute the overlapping areas if a point were to be introduced, though I am not certain on the details of that. They mention the Bowyer-Watson tessellation algorithm, and as I understand it, one would have to recompute the Delaunay triangulation change (in terms of 1-ring neighbours) for each pixel.  
> **username_2**: Sure.  If you click on the link marked "Implementation Notes" at the top of the article I cited, you'll find an article describing the algorithm in detail.  It also explains how it is possible to perform a Natural Neighbor Interpolation without modifying the Delaunay triangulation (an important part of making the thing run fast).  You can find it at https://gwlucastrig.github.io/TinfourDocs/NaturalNeighborTinfourAlgorithm/index.html .  Also, if you download the code from Github, there is a Java class called NaturalNeighborInterpolator that implements the algorithm.  
> **username_2**: Also, you may find it useful to look at Liang and Hale's algorithm.  Although the algorithm I use is largely my own, their earlier work used a lot of the same ideas I applied. You may find their approach more to your liking (in either case, I will be interested in hearing more about what you do). See https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.174.2191&rep=rep1&type=pdf  
> **username_1**: Thank you. I have found this too: http://alexbeutel.com/papers/tsas-terranni.pdf but it seems fairly involved. What I do is sparse data interpolation, and I have been using some PDEs for that, but I realize that generalized barycentric coordinates based approaches can be faster, and have other properties. I was interested in case I ever get to implementing it, but it's currently not high up on my priority queue.  

---

