Can Mathematica interpolate non-uniform scatter data?
in [Mathematica]
|
Prev: Sum and NSum problem
Next: Access to Open Excel Cell
From: Lawrence Teo on 1 Feb 2010 06:13 Hi DH, This is so upsetting. I need the answer so urgently but just can't see your reply. :) On 29 Jan, 20:43, dh <d...(a)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...(a)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`= , > >>> 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...
From: JH on 5 Feb 2010 07:13 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. JH
From: Frank Iannarilli on 13 Feb 2010 05:24
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] *) |