eljungsk / ScatteredInterpolation.jl

Interpolation of scattered data
Other
21 stars 10 forks source link

Scattered linear interpolation? #22

Open marianoarnaiz opened 3 years ago

marianoarnaiz commented 3 years ago

Hi everyone. Not to complaint about the package it is actually my favourite interpolation in Julia. I was wondering it there is a way to using linear interpolation with it? I have some data that cannot overshoot and linear would be fantastic! Kind regards. M.

eljungsk commented 3 years ago

Hi, There is unfortunately no linear interpolation implemented in the package at the moment. It's in the works though, but I can't make any promises for when it will be done. /Emil

pranavkrishna30 commented 2 years ago

Have there been any updates on the addition of linear interpolation? -Pranav

marianoarnaiz commented 2 years ago

I am still waiting too hahaha. Maybe we can help!

marianoarnaiz commented 2 years ago

If any one has the equations aI will give it a try

pranavkrishna30 commented 2 years ago

I believe that you would need to use the Delaunay triangulation and find simplex around the queried point. You can then use barycentric interpolation like done in https://github.com/johwing/matlab_scatteredInterpn/blob/master/scatteredInterpn.m. Here is a pdf (https://infoscience.epfl.ch/record/99883?ln=en) that also outlines more of the theory behind the triangulation and has the matrix formulations for the calculations to calculate interpolated values.

I think the issue is that there is not a native Julia Delaunay package for higher dimensions. https://github.com/JuliaGeometry/VoronoiDelaunay.jl only does 2D to my understanding and https://github.com/eschnett/Delaunay.jl seems to just use pycall (at which point you can just pycall linear scattered interpolants from python).

I am still pretty new to a lot of the theory behind scattered interpolation so pardon any errors.

eljungsk commented 2 years ago

Hi guys,

Thanks for volunteering to help out here. There is actually already an implementation in the linear branch (see PR #23), but I haven't gotten around to writing tests and documentation yet.

Can you try that branch and see if it works for you?

marianoarnaiz commented 2 years ago

Hi @eljungsk can you let me know how to call it? I have lots of data waiting for this, I can help with the documentation and examples hahaha.

eljungsk commented 2 years ago

Hi @marianoarnaiz It should be the same API as the other methods, but with Linear() as the method:

itp = interpolate(Linear(), points, sample's)

Help with tests and docs would be very appreciated!

marianoarnaiz commented 2 years ago

Great! Give a week or so! I will message you here.

It is on the new version

pranavkrishna30 commented 2 years ago

When I was testing, I am unable to recover some of the original data points that I put into the interpolant. For example, I made up some engine thrust data that you might typically get. I get a domain error when querying some of the original datapoints at the boundary. Here is a link to the very simple script with the data: https://github.com/pranavkrishna30/testing_repo/blob/main/julia/ScatterdInterp/linear_scattered_interp.jl

eljungsk commented 2 years ago

Thanks for testing, @pranavkrishna30! I've updated the code to allow for some rounding error at the domain boundaries. Your testcase should work now. Maybe you can give it another try?

pranavkrishna30 commented 2 years ago

I gave it another try it seems to error on another point (0.9, 20000) when trying to recover the original points. If you change the tolerance to 8*eps() in find_simplex(), it resolves the error but there is definitely some small rounding error on the output. For my understanding, what is causing the rounding error/need for tolerance in the first place?

carlocastoldi commented 1 year ago

I am pretty new to the problem, but I tried to use ScatteredInterpolation.jl with the linear branch. I register there is a considerable difference in computation time with Linear() method rather then Multiquadratic(). Moreover, distributing a package with a dependency on python+scipy is quite hard.

As I have stated earlier I know nothing about the need of Delaunay algorithm for computing a linear interpolation, but will the final version skip python+scipy all together? Perhaps interfacing with QHull directly. Could QHull.jl help in any way?

Lastly, I know that ScatteredInterpolation.jl is about interpolating in any dimension, but would it be possible to create a custom (fast-er) 2D scattered linear interpolator using the above-mentioned VoronoiDelaunay.jl?

carlocastoldi commented 1 year ago

After a couple of benchmarks in my use case, I see that Linear() method with #23 is about 50 to 60 times slower than Multiquadratic()

carlocastoldi commented 1 year ago

Excuse me, it may sound stupid but isn't the linear interpolation just a Polyharmonic with k=1?

marianoarnaiz commented 12 months ago

Hi everyone. Sorry I have been away. I cannot make the Linear() option work.

iPDF =  interpolate(Linear(), d[:,1:3]', d[:,4]);
ERROR: UndefVarError: `Linear` not defined
Stacktrace:
 [1] top-level scope

Any idea? Also, is Polyharmonic(1) really linear like @carlocastoldi says hehehe. Let me know

carlocastoldi commented 12 months ago

From what i see here, i would say so: https://en.wikipedia.org/wiki/Polyharmonic_spline#Examples

marianoarnaiz commented 11 months ago

From what i see here, i would say so: https://en.wikipedia.org/wiki/Polyharmonic_spline#Examples

Great! Thanks!