I have some problem with interpolation on scattered data. griddata(...) or interp2(...) functions do not work well. Delaunay function could triangulate data but maybe not well in all cases.

I have function to find nearest points. Interpolate between the three nearest point should if I think correct be equal to interpolation on delaunay triangulation but me a little bit stupid and can't immidiately figure out the equation. Anyone have it at hand?

I however found two other very interesting functions here http://fourier.eng.hmc.edu/e176/lectures/ch7/node7.html called "Radial Basis Function Method" and "Shepard method". I am however a little bit uncertain what happen then density vary for example then many points are tightly spaced as I expect them to pile up. Anyone used any of these?

On Thu, Mar 12, 2020, 10:49 AM Nicklas Karlsson <[hidden email]> wrote:

I have some problem with interpolation on scattered data. griddata(...) or interp2(...) functions do not work well. Delaunay function could triangulate data but maybe not well in all cases.

I have function to find nearest points. Interpolate between the three nearest point should if I think correct be equal to interpolation on delaunay triangulation but me a little bit stupid and can't immidiately figure out the equation. Anyone have it at hand?

I however found two other very interesting functions here http://fourier.eng.hmc.edu/e176/lectures/ch7/node7.html called "Radial Basis Function Method" and "Shepard method". I am however a little bit uncertain what happen then density vary for example then many points are tightly spaced as I expect them to pile up. Anyone used any of these?

I'm not familiar with those algorithms, but what difficulty are you having with a griddata method? What have you tried? Matlab has created a scatteredInterpolant function that works well for the whole process but but it is not yet implemented in octave.

I found a workaround for one of my data sets (a non uniform triangular mesh) following the example from:

If you do not have too many points (< 1000) you can use mesh-less
methods. they are known with many names, maongst th emost ppopular:
Gaussian processes (include the radial basis functions) or krigging.
See if regress_gp in the statistics package[1] solves your problem, if
not try the stk package[2] or gpml[3]

(careful, means the -th point in
your data, so x(i,:) or x(:,i) in your code)

If you're trying to use mesh techniques, I guess you will prefer
this solution : it's computationally less demanding, and it makes
the interpolant function smooth.

In this case, when the density is high, and points "pile up", then
the estimated will be
averaged over those points : it's up to you to decide whether this
corresponds to your model or not.

In the formula above, you can restrict to to the
3 nearest neighbours, or use the whole data set if you can afford
the computational load.

For the radial basis function, you can use for instance a gaussian
kernel :

if is low,
the interpolant will be "blurry", if is high it
will be "spiky".

On 12/03/2020 15:53, Nicklas Karlsson wrote:

I have some problem with interpolation on scattered
data. griddata(...) or interp2(...) functions do not work well.
Delaunay function could triangulate data but maybe not well in
all cases.

I have function to find nearest points. Interpolate between
the three nearest point should if I think correct be equal to
interpolation on delaunay triangulation but me a little bit
stupid and can't immidiately figure out the equation. Anyone
have it at hand?

I however found two other very interesting functions here http://fourier.eng.hmc.edu/e176/lectures/ch7/node7.html called
"Radial Basis Function Method" and "Shepard method". I am
however a little bit uncertain what happen then density vary
for example then many points are tightly spaced as I expect
them to pile up. Anyone used any of these?

Then thinking again. Sometimes point will end up outside the area in between the three nearest neighbours. Or maybe more correctly expressed in mathematical terms, point sometimes end up outside the convex hull spent up by the three nearest neighbours.

Maybe the weighted estimator will work, have to think. Whole data set will most probably work well, there are only about 100 points and guess plot figure use more exuction time. Averaging over points piled up is the correct solution.

For some reason Octave delaunay(...) function return a screwed up mesh, some angles are very small while others are very wide. This might be because most of the points are rather squared with a small scew while in one corner points in the outer corner is squeezed to a line because of limited signal. Then mesh is screwed up interpolation does not work well.

Thanks, Nicklas Karlsson

Den mån 16 mars 2020 kl 10:02 skrev Augustin Lefèvre <[hidden email]>:

You can use 3-nearest neighbour interpolation :

y=13∑j=13yjy = \frac{1}{3} \sum_{j=1}^3 y_j

where j is the index of the three nearest neighbours (based on xx values)

(careful, means the -th point in
your data, so x(i,:) or x(:,i) in your code)

If you're trying to use mesh techniques, I guess you will prefer
this solution : it's computationally less demanding, and it makes
the interpolant function smooth.

In this case, when the density is high, and points "pile up", then
the estimated will be
averaged over those points : it's up to you to decide whether this
corresponds to your model or not.

In the formula above, you can restrict to to the
3 nearest neighbours, or use the whole data set if you can afford
the computational load.

For the radial basis function, you can use for instance a gaussian
kernel :

if is low,
the interpolant will be "blurry", if is high it
will be "spiky".

On 12/03/2020 15:53, Nicklas Karlsson wrote:

I have some problem with interpolation on scattered
data. griddata(...) or interp2(...) functions do not work well.
Delaunay function could triangulate data but maybe not well in
all cases.

I have function to find nearest points. Interpolate between
the three nearest point should if I think correct be equal to
interpolation on delaunay triangulation but me a little bit
stupid and can't immidiately figure out the equation. Anyone
have it at hand?

I however found two other very interesting functions here http://fourier.eng.hmc.edu/e176/lectures/ch7/node7.html called
"Radial Basis Function Method" and "Shepard method". I am
however a little bit uncertain what happen then density vary
for example then many points are tightly spaced as I expect
them to pile up. Anyone used any of these?

Regards Nicklas Karlsson

Den mån 16 mars 2020 kl 10:02 skrev Augustin Lefèvre <[hidden email]>:

You can use 3-nearest neighbour interpolation :

y=13∑j=13yjy = \frac{1}{3} \sum_{j=1}^3 y_j

where j is the index of the three nearest neighbours (based on xx values)

(careful, means the -th point in
your data, so x(i,:) or x(:,i) in your code)

If you're trying to use mesh techniques, I guess you will prefer
this solution : it's computationally less demanding, and it makes
the interpolant function smooth.

In this case, when the density is high, and points "pile up", then
the estimated will be
averaged over those points : it's up to you to decide whether this
corresponds to your model or not.

In the formula above, you can restrict to to the
3 nearest neighbours, or use the whole data set if you can afford
the computational load.

For the radial basis function, you can use for instance a gaussian
kernel :

if is low,
the interpolant will be "blurry", if is high it
will be "spiky".

On 12/03/2020 15:53, Nicklas Karlsson wrote:

I have some problem with interpolation on scattered
data. griddata(...) or interp2(...) functions do not work well.
Delaunay function could triangulate data but maybe not well in
all cases.

I have function to find nearest points. Interpolate between
the three nearest point should if I think correct be equal to
interpolation on delaunay triangulation but me a little bit
stupid and can't immidiately figure out the equation. Anyone
have it at hand?

I however found two other very interesting functions here http://fourier.eng.hmc.edu/e176/lectures/ch7/node7.html called
"Radial Basis Function Method" and "Shepard method". I am
however a little bit uncertain what happen then density vary
for example then many points are tightly spaced as I expect
them to pile up. Anyone used any of these?

On Mon, Mar 23, 2020 at 6:08 AM Nicklas Karlsson <[hidden email]> wrote:

Then thinking again. Sometimes point will end up outside the area in between the three nearest neighbours. Or maybe more correctly expressed in mathematical terms, point sometimes end up outside the convex hull spent up by the three nearest neighbours.

Maybe the weighted estimator will work, have to think. Whole data set will most probably work well, there are only about 100 points and guess plot figure use more exuction time. Averaging over points piled up is the correct solution.

For some reason Octave delaunay(...) function return a screwed up mesh, some angles are very small while others are very wide. This might be because most of the points are rather squared with a small scew while in one corner points in the outer corner is squeezed to a line because of limited signal. Then mesh is screwed up interpolation does not work well.

Matlab has a 'scatteredInterpolant' function that is still in the wish list for Octave. But Octave does have griddata. you give the function your scattered data points as three vectors x, y, f(x,y), and request output points (or grid) xi, yi, and it will apply one of the methods it has implemented. Currently only
'nearest neighbor' interpolation and linear interpolation are implemented in the octave releases. I just implemented the Biharmonic spline interpolation (smooth interpolant, the 'v4' option in matlab) in bug #33539 but haven't worked up a full patch yet. Also looking for proper matlab compatible cubic and 'natural' method interpolations.

linear, cubic, and natural methods are all triangulation based. they take the scattered points, use octave's delauney triangulation, then calculate local interpolants within the triangles. linear (actually bilinear i guess) is simply a weighted some of triangle areas (the interpolation point creates three sub-triangles within the mesh triangle, and value is a weighted projection inversely proportional to sub-triangle size, where a big sub-triangle has you more strongly weighted toward the point not associated with that sub-triangle). the v4 method does not triangulate, and does a weighted interpolation over the entire mesh. that can get slow and memory intensive, but it's currently the only 'smooth' interpolation function Octave has implemented (sort of).

it sounds like griddata should do what you want, unless I've misread the thread.

>linear, cubic, and natural methods are all triangulation based. they take
>the scattered points, use octave's delauney triangulation, then calculate
>local interpolants within the triangles. linear (actually bilinear i guess)
>is simply a weighted some of triangle areas (the interpolation point
>creates three sub-triangles within the mesh triangle, and value is a
>weighted projection inversely proportional to sub-triangle size, where a
>big sub-triangle has you more strongly weighted toward the point not
>associated with that sub-triangle). the v4 method does not triangulate,
>and does a weighted interpolation over the entire mesh. that can get slow
>and memory intensive, but it's currently the only 'smooth' interpolation
>function Octave has implemented (sort of).

Have you considered Kriging? I have used the stk package from
Octave-forge with profit.

--
Francesco Potortì (ricercatore) Voice: +39.050.621.3058
ISTI - Area della ricerca CNR Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa Skype: wnlabisti
(gate 20, 1st floor, room C71) Web: http://fly.isti.cnr.it

As there are several values for each point I where able to get a good triangulation of area by supplying other values for the same indexes. Maybe this is a good time implementing the missing function, it should not be to hard for linear interpolation.

Interpolation using delaunay triangulation is useful then working with data both from measured data and from solved partial differential equations, for example using the Finite Element Method.

Nicklas Karlsson

Den mån 23 mars 2020 kl 14:57 skrev Francesco Potortì <[hidden email]>:

>linear, cubic, and natural methods are all triangulation based. they take
>the scattered points, use octave's delauney triangulation, then calculate
>local interpolants within the triangles. linear (actually bilinear i guess)
>is simply a weighted some of triangle areas (the interpolation point
>creates three sub-triangles within the mesh triangle, and value is a
>weighted projection inversely proportional to sub-triangle size, where a
>big sub-triangle has you more strongly weighted toward the point not
>associated with that sub-triangle). the v4 method does not triangulate,
>and does a weighted interpolation over the entire mesh. that can get slow
>and memory intensive, but it's currently the only 'smooth' interpolation
>function Octave has implemented (sort of).

Have you considered Kriging? I have used the stk package from
Octave-forge with profit.

--
Francesco Potortì (ricercatore) Voice: +39.050.621.3058
ISTI - Area della ricerca CNR Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa Skype: wnlabisti
(gate 20, 1st floor, room C71) Web: http://fly.isti.cnr.it

On Mon, Mar 23, 2020 at 10:45 AM Nicklas Karlsson
<[hidden email]> wrote:
>
> As there are several values for each point I where able to get a good triangulation of area by supplying other values for the same indexes. Maybe this is a good time implementing the missing function, it should not be to hard for linear interpolation.

If you are referring to the ScatteredInterpolant function, I'm not
sure if Octave can yet support the classdef features required to make
it fully compatible. Someone did create a temporary 'wrapper' that
never got evaluated/implemented here:
https://savannah.gnu.org/bugs/?35821

It would take some work to make sure it still applies to the current
codebase and actually works.

Regarding griddata, the linear interpolation based on delaunay
triangulation already exists. it's the 'cubic' and 'natural' methods
still missing.

Le 23/03/2020 à 14:57, Francesco Potortì a écrit :

>> linear, cubic, and natural methods are all triangulation based. they take
>> the scattered points, use octave's delauney triangulation, then calculate
>> local interpolants within the triangles. linear (actually bilinear i guess)
>> is simply a weighted some of triangle areas (the interpolation point
>> creates three sub-triangles within the mesh triangle, and value is a
>> weighted projection inversely proportional to sub-triangle size, where a
>> big sub-triangle has you more strongly weighted toward the point not
>> associated with that sub-triangle). the v4 method does not triangulate,
>> and does a weighted interpolation over the entire mesh. that can get slow
>> and memory intensive, but it's currently the only 'smooth' interpolation
>> function Octave has implemented (sort of).
> Have you considered Kriging? I have used the stk package from
> Octave-forge with profit.

Happy to read that ;-)

I'm sorry I haven't been following this thread carefully, but if anybody
needs help to perform scattered-data interpolation with stk, feel free
to ask here or on [hidden email].

It's not possible to send your own triangulation to griddata(...) and
I need to do because griddata(...) function does not work, it actually
crash Octave. Data come from measurement on many different point in
system and hopefully same triangulation work for interpolation for all
values.
>
> Think the best option is to modify griddata(...) function to my need. It only took a minute or so to find the griddata function on filesystem, it is a rather short function. Need to supply my own triangulation and it is also useful to show triangulation used so to check it is reasonable.
>
> Only drawback modifying griddata(...) function is I have to make it available to others to fulfill license requirements. ;)
>

(FYI please keep the mailing list in the conversation, and it is the
convention of this mailing list is to bottom post. see
http://www.idallen.com/topposting.html )

can you provide a sample of the data that is crashing Octave? I'd be
curious to see what the specific problem is.

Also, I'm no expert on licensing, but my understanding is you can make
and keep private any changes to the code you want unless you
distribute your modified program in some fashion, in which case you
need to make the source code available along with it.

That said, if it improves the function while maintaining compatibility
I'm sure we'd love to consider any improvement to the code.

> Il giorno 23 mar 2020, alle ore 18:09, Nicholas Jankowski <[hidden email]> ha scritto:
>
> It's not possible to send your own triangulation to griddata(...) and
> I need to do because griddata(...) function does not work, it actually
> crash Octave. Data come from measurement on many different point in
> system and hopefully same triangulation work for interpolation for all
> values.
>>
>> Think the best option is to modify griddata(...) function to my need. It only took a minute or so to find the griddata function on filesystem, it is a rather short function. Need to supply my own triangulation and it is also useful to show triangulation used so to check it is reasonable.
>>
>> Only drawback modifying griddata(...) function is I have to make it available to others to fulfill license requirements. ;)
>>
>
> (FYI please keep the mailing list in the conversation, and it is the
> convention of this mailing list is to bottom post. see
> http://www.idallen.com/topposting.html )
>
> can you provide a sample of the data that is crashing Octave? I'd be
> curious to see what the specific problem is.
>
> Also, I'm no expert on licensing, but my understanding is you can make
> and keep private any changes to the code you want unless you
> distribute your modified program in some fashion, in which case you
> need to make the source code available along with it.
>
> That said, if it improves the function while maintaining compatibility
> I'm sure we'd love to consider any improvement to the code.
>

I actually already have a modified version which I used locally.
I attach it in case it is of use to anyone else.
c.

> I actually already have a modified version which I used locally.
> I attach it in case it is of use to anyone else.

interesting and fairly simple extension. I think a bit about how that
could be added in a compatible fashion, it may be worth including as
an extension of the function. (input handling gets just a touch more
complicated since a recent patch makes griddata accept 3D coordinates
for matlab compatibility, but should be manageable)

Den mån 23 mars 2020 kl 18:30 skrev c. <[hidden email]>:

> Il giorno 23 mar 2020, alle ore 18:09, Nicholas Jankowski <[hidden email]> ha scritto:
>
> It's not possible to send your own triangulation to griddata(...) and
> I need to do because griddata(...) function does not work, it actually
> crash Octave. Data come from measurement on many different point in
> system and hopefully same triangulation work for interpolation for all
> values.
>>
>> Think the best option is to modify griddata(...) function to my need. It only took a minute or so to find the griddata function on filesystem, it is a rather short function. Need to supply my own triangulation and it is also useful to show triangulation used so to check it is reasonable.
>>
>> Only drawback modifying griddata(...) function is I have to make it available to others to fulfill license requirements. ;)
>>
>
> (FYI please keep the mailing list in the conversation, and it is the
> convention of this mailing list is to bottom post. see
> http://www.idallen.com/topposting.html )
>
> can you provide a sample of the data that is crashing Octave? I'd be
> curious to see what the specific problem is.
>
> Also, I'm no expert on licensing, but my understanding is you can make
> and keep private any changes to the code you want unless you
> distribute your modified program in some fashion, in which case you
> need to make the source code available along with it.
>
> That said, if it improves the function while maintaining compatibility
> I'm sure we'd love to consider any improvement to the code.
>

I actually already have a modified version which I used locally.
I attach it in case it is of use to anyone else.
c.

Tried it and it works!

griddata(...) function work on gridded data while this function work on triangulation. Will write a few a words about it tonight and suggest it is added to Octave among the "tri" functions as it belong there.

There already are triplot(...) trimesh(...) trisurf(...) and suggest tridata(...) would be a good name for this function. It is probably a good idea to use similar parameters as on the other tri functions tridata(tri, x, y, xi, yi) but are a little bit unsure about method if higher degree was ever to be added, know lagrange polynoms are commonly used for partial differential equations.

Nicklas Karlsson:
>> Only drawback modifying griddata(...) function is I have to make it
>> available to others to fulfill license requirements. ;)

Nicklas, no free software licence requires that. You can modify
whatever you want for in-house purposes.

Nicholas Jankowski:
>Also, I'm no expert on licensing, but my understanding is you can make
>and keep private any changes to the code you want unless you
>distribute your modified program in some fashion, in which case you
>need to make the source code available along with it.

Yes

--
Francesco Potortì (ricercatore) Voice: +39.050.621.3058
ISTI - Area della ricerca CNR Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa Skype: wnlabisti
(gate 20, 1st floor, room C71) Web: http://fly.isti.cnr.it