Can Mathematica interpolate non-uniform scatter data?
Lawrence Teo
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...
dh
to do it. Anyway, here is it:
now we can do the interpolation:
Lawrence Teo
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=
dh
lately my post get mutilated, we are still seraching the reason. try again:
JH
> Hi DH,
>
> Did you miss out anything in your reply?
> How to do the interpolation with non-uniform scattered data (with
> holes)?
>
>
> > 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
Lawrence Teo
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`=
JH
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
Frank Iannarilli
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]
*)