Re: Can Mathematica interpolate non-uniform scatter data?
- To: mathgroup at smc.vnet.net
- Subject: [mg107466] Re: Can Mathematica interpolate non-uniform scatter data?
- From: Frank Iannarilli <frankeye at cox.net>
- Date: Sat, 13 Feb 2010 05:24:13 -0500 (EST)
- References: <hjmjgg$ee$1@smc.vnet.net> <hjon95$4ga$1@smc.vnet.net>
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]
*)