Re: Unstructured Interpolation in V7?

*To*: mathgroup at smc.vnet.net*Subject*: [mg109698] Re: Unstructured Interpolation in V7?*From*: "Ingolf Dahl" <ingolf.dahl at telia.com>*Date*: Wed, 12 May 2010 07:34:38 -0400 (EDT)

Hi Luc and all you other, who have asked for scattered point interpolation. I am developing an interpolation package for Mathematica, named "Obtuse", which contains four methods for doing scattered point interpolation. Now I am looking for "test pilots" for my package. There are still some small details in the documentation that need to fixed, but I could send a preliminary version within a few days. Please send an email if you are interested (to ingolf.dahl at telia.com). If I want to interpolate using the testData set by Luc, I write Needs["Obtuse`"]; interpol = Interpolation[testData, Method -> "ObtuseAngle"]; interpol[{-122.000`, 38.000`}] Running this I obtain {0.0384003, 0.0490968} In the testData set from Luc, only x and y coordinates are specified. It would be possible to include the z or t coordinates also, or both z and t. The four interpolation methods are defined as extra options to the Interpolation function: The Delaunay interpolation method (option Method->"Delaunay") to the Interpolation command is a coordinate-based interpolation method, in my package allowing only 2D-points. The method is based on Delaunay triangulation (from the Computational Geometry Package) of the control points and triangular interpolation inside the Delaunay triangles. The Mathematica routines were initially written by Tom Wickham-Jones and modified by Daniel Huber (Switzerland) and Frank Iannarilli. I have modified them further, by introducing an extrapolation method to avoid error messages for extrapolated values, and by changing the calls to the routine to agree with the Interpolation command. The other three implemented methods are distance-based, and they interact with the control points via a distance matrix and a distance function. The distance-based interpolation methods could be applied to multidimensional points or to other types of objects, as long as the relevant distance function is specified. The Shepard interpolation (option Method->"Shepard") is included just for comparison, and should not be used for serious work without careful consideration. The Shepard interpolation is using the distance to a control point, raised to the inverse ShepardPower parameter, as a weight for the function value in the control point. The weights are then normalized to have 1 as their sum. If the value of ShepardPower is less than the dimension of the point set plus one, the distant points will dominate the interpolated value, and this makes it difficult to use the method if the dimensionality of the control point set is high. However, what is important here is the real dimensionality of the control point set, which might be lower than the dimensionality of the embedding coordinate system (if we have such a system). The method gives nice linear interpolation between two points if ShepardPower is one, but only for this value. Thus we might conclude that the Shepard interpolation is working best for interpolation in one-dimensional data sets. The Shepard interpolation method does not require any length scale parameter. The Radial Basis Function interpolation (option Method->"RBF") can be seen as an improvement of the Shepard interpolation. We let each control point be the center for one radial base function (a function where the value is determined by the distance to the center), and form a linear combination of all these functions. We choose the weights in this linear combination in such a way that it takes the function values in the control points. Since we have the same number of weight as we have function values, we should in general have a well-specified equation system to solve. The radial base functions could be chosen in various ways, but a fundamental property should be that they are finite but non-zero at origin and smoothly decreasing or increasing outwards. In the definition of the radial base function there is as a rule a length scale parameter included, so if the length scale vary strongly there might be a problem to choose one single appropriate value. For smooth functions the accuracy of the RBF interpolation can be amazingly good. The RBF interpolation method is translated from another system routines described by William Press in the course "Computational Statistics with Application to Bioinformatics". The ObtuseAngle interpolation (option Method->"ObtuseAngle") is based on the ObtuseAngle connection list. One starts with a weight function for each control point that is one everywhere. Then, if there is a connection from control point 1 to control point 2, the weight function for control point two is multiplied by a "shadow" function, varying only in a direction parallel to the connection. In the case that InterpolationOrder is one, this function is zero from point 1 and beyond, one from point 2 and beyond, and increases linearly from point 1 to point 2. In this way, the weight function will be zero in all points shadowed by point 1 from point 2, according to the Obtuse Angle shadowing rule. The same operation is applied to all connections. Then the weight for all control point are normalized to have 1 as their sum at the interpolation point. For InterpolationOrder equal to one, we always know that the interpolation coefficients are in the closed interval 0 to 1. For some applications this is essential. In the case that InterpolationOrder is two, the shadow function is chosen in such a way that it has a continuous derivative. At the same time it loses its ability to shadow points behind the from-point in the connections, and we have to introduce a second-level shadowing function for the second-level connections to zero the weight function for more distant points. If InterpolationOrder is three, the shadow function has second order continuous derivative, both the first-level, second-level and the third-level connections are used. As long as the options CutoffRadius and SmoothenDistance are not set, the ObtuseAngle interpolation method is scale-independent and does not need any scale parameter. The ObtuseAngle interpolation method is my invention, and the method has some qualities that should make it into an interesting alternative to RBF interpolation in some applications. The function values for the control points could be any kind of mathematical objects that are such that they can be linearly combined. The Delaunay and ObtuseAngle interpolation methods are typically local interpolation methods, where the function value in one control point often influence only the local neighborhood. The RBF interpolation is more of a global interpolation method, where the function value in one control point might influence the interpolated value almost everywhere. The Obtuse package also provides Mathematica functions to find and investigate networks created according to the obtuse-angle shadowing rule. Ingolf Dahl ingolf.dahl at telia.com > -----Original Message----- > From: Luc Barthelet [mailto:luc at tirnua.com] > Sent: den 11 maj 2010 12:25 > To: mathgroup at smc.vnet.net > Subject: [mg109658] Unstructured Interpolation in V7? > > I have a small data set (wind in the bay area, automatically extracted > by Mathematica) > > testData={{{-121.082`,38.955`},{- > 0.017250000000004206`,0.02987787643056= > 0946`}},{{- > 122.05`,37.98`},{0.03175000000000239`,0.054992613140314006`}},= > {{-122.12`,37.66`},{0.06795173495784468`,-0.011981724259015891`}},{{- > 121.= > 82`,37.7`},{0.06795173495784468`,0.011981724259015891`}},{{- > 121.3`,38.58`= > },{0.02599999999999625`,0.04503332099679369`}},{{- > 122.05`,37.43`},{0.0380= > 57551141832846`,0.01385181580469208`}},{{- > 122.28`,38.21`},{0.070476946558= > 94077`,0.02565151074942662`}},{{- > 122.556`,38.1436`},{0.05121000315664048`= > ,0.009029705238681629`}},{{- > 121.957`,38.378`},{0.03241939541710792`,0.011= > 79969494473454`}},{{- > 122.22`,37.7`},{0.09060231327711676`,0.0159756323453= > 5926`}},{{-122.12`,37.48`},{0.05662644579820153`,- > 0.009984770215851313`}}= > ,{{-121.6`,38.7`},{0.018640840680916426`,0.022215288850453874`}},{{- > 121.5= > `,38.52`},{0.02603289819229815`,0.031024799946315795`}},{{- > 122.25`,37.53`= > },{0.05121000315664048`,0.009029705238681629`}},{{- > 122.36`,37.62`},{0.051= > 21000315664048`,0.009029705238681629`}},{{- > 121.93`,37.37`},{0.06495190528= > 383432`,0.03750000000000142`}},{{- > 121.82`,37.34`},{0.06795173495784468`,0= > =2E011981724259015891`}},{{- > 122.81`,38.51`},{0.03983716857408126`,0.02300= > 000000000324`}},{{- > 121.25`,37.91`},{0.05121000315664048`,0.00902970523868= > 1629`}},{{- > 121.93`,38.28`},{0.09725818625133797`,0.03539908483420362`}},{= > {-122.465`,37.807`},{0.03464101615138304`,0.020000000000003126`}},{{- > 122.= > 21`,37.507`},{0.07386058147591257`,0.013023613325017891`}},{{- > 122.298`,37= > =2E772`},{0.05967048141990006`,-0.021718279101179405`}},{{- > 122.33`,37.8`}= > ,{0.1039999999999992`,0.`}},{{- > 122.4`,37.928`},{0.029877876430560946`,0.0= > 172499999999971`}},{{-122.975`,37.997`},{0.04530115663855838`,- > 0.00798781= > 617267963`}},{{- > 122.038`,38.057`},{0.048643822138060955`,0.04081701321509= > 712`}},{{-122.881`,37.361`},{0.06598211945181731`,- > 0.011634427903686628`}= > },{{-122.833`,37.759`},{0.055999999999997385`,0.`}},{{- > 122.833`,37.759`},= > {0.055999999999997385`,0.`}},{{- > 122.026`,38.223`},{0.07924132444627219`,0= > =2E04574999999999818`}},{{- > 122.447`,37.891`},{0.005035797152345367`,0.028= > 559424837354186`}}}; > > That interpolate, extrapolate nicely with > > ListVectorPlot[testData] > > Now, I would like to build the interpolated data outside of > ListVectorPlot, > > but Interpolation[testData] returns an error: > Interpolation::indim: The coordinates do not lie on a structured tensor > product grid > > Which is an old issue (for MathGroup at least). > > The solution are to use the IMS package imsUnstructuredInterpolation or > TWJ's TriangularInterpolation. (or Cases on ListVectorPlot.., which is > my current solution) > > But is there a new better way to do that in version 7? the > documentation > avoids the issue and nicely just uses Table to create dataset from a > function! (Which I guess is a hint that this is not resolved) > > Thanks, > Luc > > > > >