MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

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]


*)




  • Prev by Date: Re: Question about subscripts and polynomial
  • Next by Date: Bug in Mathematica: regular expression applied to very long string
  • Previous by thread: Re: Can Mathematica interpolate non-uniform scatter data?
  • Next by thread: Re: Using variable name as string AND value ?