Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""rbf - Radial basis functions for interpolation/smoothing scattered N-D data. 

2 

3Written by John Travers <jtravs@gmail.com>, February 2007 

4Based closely on Matlab code by Alex Chirokov 

5Additional, large, improvements by Robert Hetland 

6Some additional alterations by Travis Oliphant 

7Interpolation with multi-dimensional target domain by Josua Sassen 

8 

9Permission to use, modify, and distribute this software is given under the 

10terms of the SciPy (BSD style) license. See LICENSE.txt that came with 

11this distribution for specifics. 

12 

13NO WARRANTY IS EXPRESSED OR IMPLIED. USE AT YOUR OWN RISK. 

14 

15Copyright (c) 2006-2007, Robert Hetland <hetland@tamu.edu> 

16Copyright (c) 2007, John Travers <jtravs@gmail.com> 

17 

18Redistribution and use in source and binary forms, with or without 

19modification, are permitted provided that the following conditions are 

20met: 

21 

22 * Redistributions of source code must retain the above copyright 

23 notice, this list of conditions and the following disclaimer. 

24 

25 * Redistributions in binary form must reproduce the above 

26 copyright notice, this list of conditions and the following 

27 disclaimer in the documentation and/or other materials provided 

28 with the distribution. 

29 

30 * Neither the name of Robert Hetland nor the names of any 

31 contributors may be used to endorse or promote products derived 

32 from this software without specific prior written permission. 

33 

34THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 

35"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 

36LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 

37A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 

38OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 

39SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 

40LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 

41DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 

42THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 

43(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 

44OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 

45""" 

46import numpy as np 

47 

48from scipy import linalg 

49from scipy.special import xlogy 

50from scipy.spatial.distance import cdist, pdist, squareform 

51 

52__all__ = ['Rbf'] 

53 

54 

55class Rbf(object): 

56 """ 

57 Rbf(*args) 

58 

59 A class for radial basis function interpolation of functions from 

60 N-D scattered data to an M-D domain. 

61 

62 Parameters 

63 ---------- 

64 *args : arrays 

65 x, y, z, ..., d, where x, y, z, ... are the coordinates of the nodes 

66 and d is the array of values at the nodes 

67 function : str or callable, optional 

68 The radial basis function, based on the radius, r, given by the norm 

69 (default is Euclidean distance); the default is 'multiquadric':: 

70 

71 'multiquadric': sqrt((r/self.epsilon)**2 + 1) 

72 'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1) 

73 'gaussian': exp(-(r/self.epsilon)**2) 

74 'linear': r 

75 'cubic': r**3 

76 'quintic': r**5 

77 'thin_plate': r**2 * log(r) 

78 

79 If callable, then it must take 2 arguments (self, r). The epsilon 

80 parameter will be available as self.epsilon. Other keyword 

81 arguments passed in will be available as well. 

82 

83 epsilon : float, optional 

84 Adjustable constant for gaussian or multiquadrics functions 

85 - defaults to approximate average distance between nodes (which is 

86 a good start). 

87 smooth : float, optional 

88 Values greater than zero increase the smoothness of the 

89 approximation. 0 is for interpolation (default), the function will 

90 always go through the nodal points in this case. 

91 norm : str, callable, optional 

92 A function that returns the 'distance' between two points, with 

93 inputs as arrays of positions (x, y, z, ...), and an output as an 

94 array of distance. E.g., the default: 'euclidean', such that the result 

95 is a matrix of the distances from each point in ``x1`` to each point in 

96 ``x2``. For more options, see documentation of 

97 `scipy.spatial.distances.cdist`. 

98 mode : str, optional 

99 Mode of the interpolation, can be '1-D' (default) or 'N-D'. When it is 

100 '1-D' the data `d` will be considered as 1-D and flattened 

101 internally. When it is 'N-D' the data `d` is assumed to be an array of 

102 shape (n_samples, m), where m is the dimension of the target domain. 

103 

104 

105 Attributes 

106 ---------- 

107 N : int 

108 The number of data points (as determined by the input arrays). 

109 di : ndarray 

110 The 1-D array of data values at each of the data coordinates `xi`. 

111 xi : ndarray 

112 The 2-D array of data coordinates. 

113 function : str or callable 

114 The radial basis function. See description under Parameters. 

115 epsilon : float 

116 Parameter used by gaussian or multiquadrics functions. See Parameters. 

117 smooth : float 

118 Smoothing parameter. See description under Parameters. 

119 norm : str or callable 

120 The distance function. See description under Parameters. 

121 mode : str 

122 Mode of the interpolation. See description under Parameters. 

123 nodes : ndarray 

124 A 1-D array of node values for the interpolation. 

125 A : internal property, do not use 

126 

127 Examples 

128 -------- 

129 >>> from scipy.interpolate import Rbf 

130 >>> x, y, z, d = np.random.rand(4, 50) 

131 >>> rbfi = Rbf(x, y, z, d) # radial basis function interpolator instance 

132 >>> xi = yi = zi = np.linspace(0, 1, 20) 

133 >>> di = rbfi(xi, yi, zi) # interpolated values 

134 >>> di.shape 

135 (20,) 

136 

137 """ 

138 # Available radial basis functions that can be selected as strings; 

139 # they all start with _h_ (self._init_function relies on that) 

140 def _h_multiquadric(self, r): 

141 return np.sqrt((1.0/self.epsilon*r)**2 + 1) 

142 

143 def _h_inverse_multiquadric(self, r): 

144 return 1.0/np.sqrt((1.0/self.epsilon*r)**2 + 1) 

145 

146 def _h_gaussian(self, r): 

147 return np.exp(-(1.0/self.epsilon*r)**2) 

148 

149 def _h_linear(self, r): 

150 return r 

151 

152 def _h_cubic(self, r): 

153 return r**3 

154 

155 def _h_quintic(self, r): 

156 return r**5 

157 

158 def _h_thin_plate(self, r): 

159 return xlogy(r**2, r) 

160 

161 # Setup self._function and do smoke test on initial r 

162 def _init_function(self, r): 

163 if isinstance(self.function, str): 

164 self.function = self.function.lower() 

165 _mapped = {'inverse': 'inverse_multiquadric', 

166 'inverse multiquadric': 'inverse_multiquadric', 

167 'thin-plate': 'thin_plate'} 

168 if self.function in _mapped: 

169 self.function = _mapped[self.function] 

170 

171 func_name = "_h_" + self.function 

172 if hasattr(self, func_name): 

173 self._function = getattr(self, func_name) 

174 else: 

175 functionlist = [x[3:] for x in dir(self) 

176 if x.startswith('_h_')] 

177 raise ValueError("function must be a callable or one of " + 

178 ", ".join(functionlist)) 

179 self._function = getattr(self, "_h_"+self.function) 

180 elif callable(self.function): 

181 allow_one = False 

182 if hasattr(self.function, 'func_code') or \ 

183 hasattr(self.function, '__code__'): 

184 val = self.function 

185 allow_one = True 

186 elif hasattr(self.function, "__call__"): 

187 val = self.function.__call__.__func__ 

188 else: 

189 raise ValueError("Cannot determine number of arguments to " 

190 "function") 

191 

192 argcount = val.__code__.co_argcount 

193 if allow_one and argcount == 1: 

194 self._function = self.function 

195 elif argcount == 2: 

196 self._function = self.function.__get__(self, Rbf) 

197 else: 

198 raise ValueError("Function argument must take 1 or 2 " 

199 "arguments.") 

200 

201 a0 = self._function(r) 

202 if a0.shape != r.shape: 

203 raise ValueError("Callable must take array and return array of " 

204 "the same shape") 

205 return a0 

206 

207 def __init__(self, *args, **kwargs): 

208 # `args` can be a variable number of arrays; we flatten them and store 

209 # them as a single 2-D array `xi` of shape (n_args-1, array_size), 

210 # plus a 1-D array `di` for the values. 

211 # All arrays must have the same number of elements 

212 self.xi = np.asarray([np.asarray(a, dtype=np.float_).flatten() 

213 for a in args[:-1]]) 

214 self.N = self.xi.shape[-1] 

215 

216 self.mode = kwargs.pop('mode', '1-D') 

217 

218 if self.mode == '1-D': 

219 self.di = np.asarray(args[-1]).flatten() 

220 self._target_dim = 1 

221 elif self.mode == 'N-D': 

222 self.di = np.asarray(args[-1]) 

223 self._target_dim = self.di.shape[-1] 

224 else: 

225 raise ValueError("Mode has to be 1-D or N-D.") 

226 

227 if not all([x.size == self.di.shape[0] for x in self.xi]): 

228 raise ValueError("All arrays must be equal length.") 

229 

230 self.norm = kwargs.pop('norm', 'euclidean') 

231 self.epsilon = kwargs.pop('epsilon', None) 

232 if self.epsilon is None: 

233 # default epsilon is the "the average distance between nodes" based 

234 # on a bounding hypercube 

235 ximax = np.amax(self.xi, axis=1) 

236 ximin = np.amin(self.xi, axis=1) 

237 edges = ximax - ximin 

238 edges = edges[np.nonzero(edges)] 

239 self.epsilon = np.power(np.prod(edges)/self.N, 1.0/edges.size) 

240 

241 self.smooth = kwargs.pop('smooth', 0.0) 

242 self.function = kwargs.pop('function', 'multiquadric') 

243 

244 # attach anything left in kwargs to self for use by any user-callable 

245 # function or to save on the object returned. 

246 for item, value in kwargs.items(): 

247 setattr(self, item, value) 

248 

249 # Compute weights 

250 if self._target_dim > 1: # If we have more than one target dimension, 

251 # we first factorize the matrix 

252 self.nodes = np.zeros((self.N, self._target_dim), dtype=self.di.dtype) 

253 lu, piv = linalg.lu_factor(self.A) 

254 for i in range(self._target_dim): 

255 self.nodes[:, i] = linalg.lu_solve((lu, piv), self.di[:, i]) 

256 else: 

257 self.nodes = linalg.solve(self.A, self.di) 

258 

259 @property 

260 def A(self): 

261 # this only exists for backwards compatibility: self.A was available 

262 # and, at least technically, public. 

263 r = squareform(pdist(self.xi.T, self.norm)) # Pairwise norm 

264 return self._init_function(r) - np.eye(self.N)*self.smooth 

265 

266 def _call_norm(self, x1, x2): 

267 return cdist(x1.T, x2.T, self.norm) 

268 

269 def __call__(self, *args): 

270 args = [np.asarray(x) for x in args] 

271 if not all([x.shape == y.shape for x in args for y in args]): 

272 raise ValueError("Array lengths must be equal") 

273 if self._target_dim > 1: 

274 shp = args[0].shape + (self._target_dim,) 

275 else: 

276 shp = args[0].shape 

277 xa = np.asarray([a.flatten() for a in args], dtype=np.float_) 

278 r = self._call_norm(xa, self.xi) 

279 return np.dot(self._function(r), self.nodes).reshape(shp)