221 views

Skip to first unread message

Jan 26, 2010, 6:26:08 AM1/26/10

to

Hi all,

Can Mathematica interpolate non-uniform scatter data?

Like ListInterpolation or Interpolation functions...

z = {{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, ,

0.12467473338522769`, 0.18640329676226988`,

0.24740395925452294`, , 0.36627252908604757`,

0.42367625720393803`, 0.479425538604203`}, {0.`,

0.12467473338522769`, 0.24740395925452294`, 0.36627252908604757`,

0.479425538604203`, 0.5850972729404622`, 0.6816387600233341`, ,

0.8414709848078965`}, {0.`, 0.18640329676226988`,

0.36627252908604757`, 0.5333026735360201`, 0.6816387600233341`,

0.806081108260693`, 0.9022675940990952`, 0.9668265566961802`,

0.9974949866040544`}, {0.`, 0.24740395925452294`,

0.479425538604203`, 0.6816387600233341`, 0.8414709848078965`,

0.9489846193555862`, 0.9974949866040544`, 0.9839859468739369`,

0.9092974268256817`}};

g = ListInterpolation[z, {{0, 1}, {0, 2}}];

g[0.5, 0.5]

However, the previous interpolation will give us a lot of NULL in the

expression...

Jan 27, 2010, 1:42:45 AM1/27/10

to

to do it. Anyway, here is it:

now we can do the interpolation:

Jan 27, 2010, 6:23:51 AM1/27/10

to

Hi DH,

Did you miss out anything in your reply?

How to do the interpolation with non-uniform scattered data (with

holes)?

On Jan 27, 2:42 pm, dh <d...@metrohm.com> wrote:

> to do it. Anyway, here is it:

>

> now we can do the interpolation:

>

> Lawrence Teo wrote:

> > Hi all,

>

> > Can Mathematica interpolate non-uniform scatter data?

> > Like ListInterpolation or Interpolation functions...

>

> > z = {{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, ,

> > 0.12467473338522769`, 0.18640329676226988`,

> > 0.24740395925452294`, , 0.36627252908604757`,

> > 0.42367625720393803`, 0.479425538604203`}, {0.`,

> > 0.12467473338522769`, 0.24740395925452294`, 0.36627252908604757=

Jan 29, 2010, 7:43:53 AM1/29/10

to

lately my post get mutilated, we are still seraching the reason. try again:

Jan 29, 2010, 7:50:20 AM1/29/10

to

On 27 ene, 12:23, Lawrence Teo <lawrence...@yahoo.com> wrote:

> Hi DH,

>

> Did you miss out anything in your reply?

> How to do the interpolation with non-uniform scattered data (with

> holes)?

>

> Hi DH,

>

> Did you miss out anything in your reply?

> How to do the interpolation with non-uniform scattered data (with

> holes)?

>

> > to do it. Anyway, here is it:

>

> > now we can do the interpolation:

>

> > Lawrence Teo wrote:

> > > Hi all,

>

> > > Can Mathematica interpolate non-uniform scatter data?

> > > Like ListInterpolation or Interpolation functions...

>

> > > z = {{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, ,

> > > 0.12467473338522769`, 0.18640329676226988`,

> > > 0.24740395925452294`, , 0.36627252908604757`,

> > > 0.42367625720393803`, 0.479425538604203`}, {0.`,

> > > 0.12467473338522769`, 0.24740395925452294`, 0.366272529086047=>

> > now we can do the interpolation:

>

> > Lawrence Teo wrote:

> > > Hi all,

>

> > > Can Mathematica interpolate non-uniform scatter data?

> > > Like ListInterpolation or Interpolation functions...

>

> > > z = {{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, ,

> > > 0.12467473338522769`, 0.18640329676226988`,

> > > 0.24740395925452294`, , 0.36627252908604757`,

> > > 0.42367625720393803`, 0.479425538604203`}, {0.`,

57=

> `,

> > > 0.479425538604203`, 0.5850972729404622`, 0.6816387600233341`,=

,

> > > 0.8414709848078965`}, {0.`, 0.18640329676226988`,

`,

> > > 0.806081108260693`, 0.9022675940990952`, 0.9668265566961802`,

> > > 0.9974949866040544`}, {0.`, 0.24740395925452294`,

> > > 0.479425538604203`, 0.6816387600233341`, 0.8414709848078965`,

,

> > > 0.9092974268256817`}};

> > > g = ListInterpolation[z, {{0, 1}, {0, 2}}];

> > > g[0.5, 0.5]

>

> > > However, the previous interpolation will give us a lot of NULL in the

> > > expression...

Hello, Lawrence:

As far as I know, the answer is NOT YET (for more than one independent

variable). The subsequent question is obvious: WHY? (I don't know,

it's strange to me).

When dealing with data you know we have several options:

1. "Total" interpolation (look for a unique function --generally a

polynomial-- that passes thru all the points, as the Lagrange or

Hermite polynomials (1D) or generalizations for more dimensions, as

Shepard's. In Mathematica, for 1 indep. vble. we have

InterpolatingPolynomial[]).

2. "Local" or piecewise interpolation (as the Mathematica:

Interpolation[])

3. Fit (approximation of data). In Mathematica: Fit, FindFit, ...

For the first type, I have implemented some algorithms due to Shepard

(in 1968, not so modern!). Here I copy my implementation of one:

points= Table[xi = \[Pi] (2 Random[] - 1);yi = \[Pi] (2 Random[] - 1);

{xi, yi, Sin[xi] Cos[yi]},{9}]; note: few points, please

graphPoints = ListPointPlot3D[points, PlotStyle -> PointSize[0.04]];

Clear[n, vi, v];

n := Dimensions[points][[1]];

Array[vi, n];

Do[

vi[i] = (\!\(

\*UnderoverscriptBox[\(\[Product]\), \(j = 1\), \(i - 1\)]

\*SuperscriptBox[\((

\*SuperscriptBox[\((x - \(points[[j]]\)[[1]])\), \(2\)] +

\*SuperscriptBox[\((y - \(points[[j]]\)[[

2]])\), \(2\)])\), \(2\)]\)) (\!\(

\*UnderoverscriptBox[\(\[Product]\), \(j = i + 1\), \(n\)]

\*SuperscriptBox[\((

\*SuperscriptBox[\((x - \(points[[j]]\)[[1]])\), \(2\)] +

\*SuperscriptBox[\((y - \(points[[j]]\)[[

2]])\), \(2\)])\), \(2\)]\)), {i, 1, n}];

v = \!\(

\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(n\)]\(vi[i]\)\);

funcionShepard2[xx_, yy_] := \!\(

\*UnderoverscriptBox[\(\[Sum]\), \(i =

1\), \(n\)]\((\(points[[i]]\)[[3]]\

\*FractionBox[\(vi[i]\), \(v\)])\)\) /. {x -> xx, y -> yy}

graphSup =

Plot3D[funcionShepard2[xx,

yy], {xx, -\[Pi], \[Pi]}, {yy, -\[Pi], \[Pi]},

PlotRange -> {{-\[Pi], \[Pi]}, {-\[Pi], \[Pi]}, {-1, 1}},

ClippingStyle -> None]

Show[graphPoints, graphSup]

_____________________

For nD (n>1) Piecewise Interpolation of scattered (irregular) data,

one way is triangulation (in Mathematica we have DelaunayTriangulation

[], see Computational Geometry Package). But the interpolation

algorithm is not included in Mathem (I think).

_____________________

In a way, we can say Mathem. DOES HAVE (local) interpolation of

scattered data (2 indep. vbles), in List3DPlot[...], we can choose the

interpolation order, and Mathem. interpolates for a pretty drawing:

BUT THEY DON'T LET US USING THAT INTERPOLATION (Why???)

_____________________

I hope this has been useful for you.

JH

Feb 1, 2010, 6:13:17 AM2/1/10

to

Hi DH,

This is so upsetting. I need the answer so urgently but just can't see

your reply. :)

This is so upsetting. I need the answer so urgently but just can't see

your reply. :)

On 29 Jan, 20:43, dh <d...@metrohm.com> wrote:

> lately my post get mutilated, we are still seraching the reason. try a=

gain:

>

> to do it. Anyway, here is it:

>

> now we can do the interpolation:

>

> Lawrence Teo wrote:

> > On Jan 27, 2:42 pm, dh <d...@metrohm.com> wrote:

> >> to do it. Anyway, here is it:

> >> now we can do the interpolation:

> >> Lawrence Teo wrote:

> >>> Hi all,

> >>> Can Mathematica interpolate non-uniform scatter data?

> >>> Like ListInterpolation or Interpolation functions...

> >>> z = {{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, ,

> >>> 0.12467473338522769`, 0.18640329676226988`,

> >>> 0.24740395925452294`, , 0.36627252908604757`,

> >>> 0.42367625720393803`, 0.479425538604203`}, {0.`,

> >>> 0.12467473338522769`, 0.24740395925452294`, 0.366272529086047=

57=

> > `,

> >>> 0.479425538604203`, 0.5850972729404622`, 0.6816387600233341`,=

,

> >>> 0.8414709848078965`}, {0.`, 0.18640329676226988`,

> >>> 0.36627252908604757`, 0.5333026735360201`, 0.6816387600233341=

`,

> >>> 0.806081108260693`, 0.9022675940990952`, 0.9668265566961802`,

> >>> 0.9974949866040544`}, {0.`, 0.24740395925452294`,

> >>> 0.479425538604203`, 0.6816387600233341`, 0.8414709848078965`,

> >>> 0.9489846193555862`, 0.9974949866040544`, 0.9839859468739369`=

Feb 5, 2010, 7:13:46 AM2/5/10

to

I was a bit mistaken: Mathematica does 'global' interpolation for more

than 1 indep. vble.: with the same command as for 1 i.v.:

InterpolatingPolynomial. But this command does not work for your data,

and my implementation of Shepard gets stuck. As your grid is regular,

and the "only" problem is that you have a few holes inside, perhaps

you can repare those values fixing one of the independent variables.

than 1 indep. vble.: with the same command as for 1 i.v.:

InterpolatingPolynomial. But this command does not work for your data,

and my implementation of Shepard gets stuck. As your grid is regular,

and the "only" problem is that you have a few holes inside, perhaps

you can repare those values fixing one of the independent variables.

JH

Feb 13, 2010, 5:24:03 AM2/13/10

to

Hi,

Here is what you want, at least for scattered {x,y,z} data:

I grabbed Tom Wickham-Jones' ExtendGraphics package, specifically the

package TriangularInterpolate. That used to internally invoke qhull

for Delaunay via MathLink. But now that Mathematica (version 6+)

includes Delauney, I simply modified Tom's function into a standalone

which calls the "native" Delaunay. Hence, my modification of his

package called TriangularInterpolateNative.

I think I tried submitting this to MathSource a while ago, with no

effect. So forgive the ugliness: I'll paste the text of the

TriangularInterpolateNative.m file later below, following an example

"notebook" of its use:

NOTEBOOK:

-------------------

Needs["ExtendGraphics`TriangularInterpolateNative`"]

?TriangularInterpolateNative

TriangularInterpolateNative[ pts] constructs a

TriangularInterpolation

function which represents an approximate function that interpolates

the data. The data must have the form {x,y,z} and do not have to

be regularly spaced.

To be able to raster-scan later, define bounding region coordinates/

values:

xyzTab = Join[Table[10 RandomReal[], {10}, {3}], {{0, 0, 0}, {0, 10,

0}, {10, 0, 0}, {10, 10, 0}}];

ListPlot3D[xyzTab, PlotRange -> All]

xyInterp = TriangularInterpolateNative[xyzTab];

xyInterp[6., 3.]

Now raster-scan this triangular interpolation into a regular grid that

Interpolation[] can handle.

xyInterpRegTab = Flatten[Table[{i, j, xyInterp[i, j]}, {i, 0, 10,

0.5}, {j, 0, 10, 0.5}], 1];

ListPlot3D[xyInterpRegTab, PlotRange -> All]

xyInterpReg = Interpolation[Map[({{#[[1]], #[[2]]}, #[[3]]}) &,

xyInterpRegTab]]

ListPlot3D[Flatten[Table[{i, j, xyInterpReg[i, j]}, {i, 0, 10, 0.5},

{j, 0, 10, 0.5}], 1], PlotRange -> All]

Timing Efficiency Comparison

Timing[Flatten[Table[{i, j, xyInterp[i, j]}, {i, 0, 10, 0.1}, {j, 0,

10, 0.1}], 1];]

{4., Null}

Timing[Flatten[Table[{i, j, xyInterpReg[i, j]}, {i, 0, 10, 0.1}, {j,

0, 10, 0.1}], 1];]

{0.125, Null}

PACKAGE: Save below into file "TriangularInterpolateNative.m" in

folder ExtendGraphics in Applications folder

==========

(* ::Package:: *)

(* :Name: TriangularInterpolateNative *)

(* :Title: TriangularInterpolateNative *)

(* :Author: Tom Wickham-Jones - made Native by F. Iannarilli *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 2.2 - Native: Version 6.X *)

(* :Summary:

This package provides functions for triangular interpolation.

*)

(* :History:

Created summer 1993 by Tom Wickham-Jones.

This package is described in the book

Mathematica Graphics: Techniques and Applications.

Tom Wickham-Jones, TELOS/Springer-Verlag 1994.

*)

(*:Warnings:

This package uses Mathematica-Native Computational Geometry package

for DelaunayTriangulation.

*)

BeginPackage[ "ExtendGraphics`TriangularInterpolateNative`",

{"ComputationalGeometry`"}]

TriangularInterpolateNative::usage =

"TriangularInterpolateNative[ pts] constructs a

TriangularInterpolation

function which represents an approximate function that interpolates

the data. The data must have the form {x,y,z} and do not have to

be regularly spaced."

TriangularInterpolating::usage =

"TriangularInterpolating[ data] represents an approximate function

whose values are found by interpolation."

Begin[ "`Private`"]

TriangularInterpolateNative[ pnts_List /;

MatrixQ[ N[ pnts], NumberQ] &&

Length[ First[ pnts]] === 3, opts___Rule]:=

TriangularInterpolateNative[ N[pnts],

DelaunayTriangulation[ Map[Take[#, 2]&, N[ pnts]], Hull->True],

opts]

TriangularInterpolateNative[ pnts_List /;

MatrixQ[ N[pnts], NumberQ] &&

Length[ First[ pnts]] === 3,

{adj_, hull_}, opts___Rule] :=

Block[{t, tri, num, work, minx, maxx, miny, maxy, tris, ext},

(* adjacency list -> triangles, courtesy dh in MathGroup:

"If we partition the neighbours into all possible pairs and add

the vertex,we obtain alltriangels that contain this vertex.

However,this gives the same triangle several times.

Therefore,we sort the triangle points and use

Unique to get rid of the multiplicity."

*)

t=Function[x,Prepend[#,x[[1]]]&/@Partition[x[[2]],2,1]]/@ adj;

tri=Sort/@Flatten[t,1] //Union;

num = Range[ Length[ tri]] ;

tris = Map[ Part[ pnts, #]&, tri] ;

work = Map[ Min[ Part[ Transpose[ #], 1]]&, tris] ;

minx = SortFun[ work] ;

work = Map[ Max[ Part[ Transpose[ #], 1]]&, tris] ;

maxx = SortFun[ work] ;

work = Map[ Min[ Part[ Transpose[ #], 2]]&, tris] ;

miny = SortFun[ work] ;

work = Map[ Max[ Part[ Transpose[ #], 2]]&, tris] ;

maxy = SortFun[ work] ;

TriangularInterpolating[

{hull, tri, pnts, minx, maxx, miny, maxy}]

]

(*

The bounding box info is stored in minx maxx miny and maxy

The format is (eg minx) { tri-number, xval}

the values always increase as you down the list.

*)

(*

Colinear data point problem

*)

SortFun[ data_] :=

Sort[

Transpose[ { Range[ Length[ data]], data}],

Less[ Part[ #1, 2], Part[#2, 2]]&]

Format[ t_TriangularInterpolating] := "TriangularInterpolating[ <> ]"

(*

FindFirstMin[ test, imin, imax, list]

list {{t1, v1}, {t2, v2}, {t3, v3}, ...}

The vi are the minimum x(y) values of the triangles.

Return the index of the first for which the test is greater

or equal than.

Return Take[ list, pos] where pos

test >= Part[ list, pos, 2] &&

test < Part[ list, pos+1, 2]

FindFirstMax[ test, imin, imax, list]

Return the index pos

test > Part[ list, pos, 2] &&

test <= Part[ list, pos+1, 2]

*)

FindFirstMin[ x_, imin_, imax_, list_] :=

Block[ {pos},

Which[

x < Part[ list, 1, 2], {},

x >= Part[ list, -1, 2], Map[ First, list],

True,

pos = FindFirstMin1[ x, imin, imax, list] ;

If[ x < Part[ list, pos, 2] ||

(pos < Length[ list] && x >= Part[ list, pos+1,2]),

Print[ "Error in Bounding Box calc"]] ;

Map[ First, Take[ list, pos]]]

]

FindFirstMin1[ x_, imin_, imax_, list_] :=

Block[{pos},

If[ imin === imax,

imin,

pos = Floor[ (imin+imax)/2] ;

If[ x < Part[ list, pos, 2],

FindFirstMin1[ x, imin, pos, list],

If[ x < Part[ list, pos+1, 2],

pos,

FindFirstMin1[ x, pos+1, imax, list]

]

]

]

]

FindFirstMax[ x_, imin_, imax_, list_] :=

Block[ {pos},

Which[

x <= Part[ list, 1, 2], Map[ First, list],

x > Part[ list, -1, 2], {},

True,

pos = FindFirstMax1[ x, imin, imax, list] ;

If[ x <= Part[ list, pos, 2] ||

(pos < Length[ list] && x > Part[ list, pos+1,2]),

Print[ "Error in Bounding Box calc"]] ;

Map[ First, Drop[ list, pos]]]

]

FindFirstMax1[ x_, imin_, imax_, list_] :=

Block[{pos},

If[ imin === imax,

imin,

pos = Floor[ (imin+imax)/2] ;

If[ x <= Part[ list, pos, 2],

FindFirstMax1[ x, imin, pos, list],

If[ x <= Part[ list, pos+1, 2],

pos,

FindFirstMax1[ x, pos+1, imax, list]

]

]

]

]

PointInTri[

{x_, y_},

{{x1_,y1_,z1_},

{x2_,y2_,z2_},

{x3_,y3_,z3_}}] :=

Block[{t1,t2,t3,eps = 5 $MachineEpsilon},

t1 = -x1 y + x2 y + x y1 - x2 y1 - x y2 + x1 y2 ;

t2 = -x2 y + x3 y + x y2 - x3 y2 - x y3 + x2 y3 ;

t3 = x1 y - x3 y - x y1 + x3 y1 + x y3 - x1 y3 ;

If[ t1 > -eps && t2 > -eps && t3 > -eps ||

t1 < eps && t2 < eps && t3 < eps,

True,

False]

]

FindTri[ pt_, tris_, pts_, rng_] :=

Block[{tst},

tst = Part[ pts, Part[ tris, First[ rng]]] ;

If[ PointInTri[ pt, tst],

First[ rng],

If[ Length[ rng] === 1,

{},

FindTri[ pt, tris, pts, Rest[ rng]]]]

]

InterpolatePointInTri[

{x_, y_},

{{x1_,y1_,z1_},

{x2_,y2_,z2_},

{x3_,y3_,z3_}}] :=

Block[{den, a, b, c},

den = x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 - x2 y3 ;

a = -(y2 z1) + y3 z1 + y1 z2 - y3 z2 - y1 z3 + y2 z3 ;

b = x2 z1 - x3 z1 - x1 z2 + x3 z2 + x1 z3 - x2 z3 ;

c = x3 y2 z1 - x2 y3 z1 - x3 y1 z2 + x1 y3 z2 + x2 y1 z3 - x1 y2

z3 ;

a/den x + b/den y + c/den

]

TriangularInterpolating::dmval =

"Input value `1` lies outside domain of the interpolating function."

TriangularInterpolating[

{hull_, tris_, pts_, minx_, maxx_, miny_, maxy_}][ x_?NumberQ, y_?

NumberQ] :=

Block[{x0, y0, x1, y1, tri, len},

len = Length[ minx] ;

x0 = FindFirstMin[ x, 1, len, minx] ; (* In x0 and above *)

x1 = FindFirstMax[ x, 1, len, maxx] ; (* In x1 and below *)

y0 = FindFirstMin[ y, 1, len, miny] ;

y1 = FindFirstMax[ y, 1, len, maxy] ;

rng = Intersection[ x0, x1, y0, y1] ;

If[ rng =!= {},

tri = FindTri[ {x, y}, tris, pts, rng],

tri = {}] ;

If[ tri === {},

Message[ TriangularInterpolating::dmval, {x,y}];

Indeterminate

(* else *) ,

tri = Part[ pts, Part[ tris, tri]] ;

InterpolatePointInTri[ {x, y}, tri]]

]

FixInd[ pts_] :=

Block[{num},

res = Select[ pts, FreeQ[ #, Indeterminate]&] ;

If[ Length[ res] === 3,

res,

{}]

]

MinMax[ data_] :=

Block[{ x},

x = Map[ First, data] ;

{Min[x], Max[x]}

]

MinMax[ data_, pos_] :=

Block[{ d},

d = Map[ Part[ #, pos]&, data] ;

{Min[d], Max[d]}

]

CheckNum[ x_ /; (Head[x] === Integer && x > 1), len_] := x

CheckNum[ x_, len_] := Ceiling[ N[ Sqrt[ len] +2]]

End[]

EndPackage[]

(*:Examples:

<<ExtendGraphics`TriangularInterpolateNative`

pnts = Table[ { x = Random[], y = Random[], x y}, {50}];

fun = TriangularInterpolate[ pnts]

fun[ .5,.7]

*)

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu