Mercurial > octave

comparison scripts/geometry/griddatan.m @ 28211:289882040316

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
griddatan.m: Rewrite docstring and redo input validation. * griddatan.m: Rewrite docstring. Redo input validation to make it more strict and error messages more specific. Add BIST tests for new input validation.
author Rik <rik@octave.org>
date Sun, 12 Apr 2020 18:14:05 -0700
parents bb929d5a34cb
children bb50407f9c60
comparison
equal deleted inserted replaced
28210:bb929d5a34cb 28211:289882040316
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}) 27 ## @deftypefn {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi})
28 ## @deftypefnx {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}) 28 ## @deftypefnx {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method})
29 ## @deftypefnx {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}, @var{options}) 29 ## @deftypefnx {} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}, @var{options})
30 ## 30 ##
31 ## Generate a regular mesh from irregular data using interpolation. 31 ## Interpolate irregular source data @var{x}, @var{y} at points specified by
32 ## @var{xi}.
32 ## 33 ##
33 ## The function is defined by @code{@var{y} = f (@var{x})}. 34 ## The input @var{x} is an MxN matrix representing M points in an N-dimensional
34 ## The interpolation points are all @var{xi}. 35 ## space. The input @var{y} is a single-valued column vector (Mx1)
36 ## representing a function evaluated at the points @var{x}, i.e.,
37 ## @code{@var{y} = fcn (@var{x})}. The input @var{xi} is a list of points
38 ## for which the function output @var{yi} should be approximated through
39 ## interpolation. @var{xi} must have the same number of columns (@var{N})
40 ## as @var{x} so that the dimensionality matches.
35 ## 41 ##
36 ## The interpolation method can be @qcode{"nearest"} or @qcode{"linear"}. 42 ## The optional input interpolation @var{method} can be @qcode{"nearest"} or
37 ## If method is omitted it defaults to @qcode{"linear"}. 43 ## @qcode{"linear"}. When the method is @qcode{"nearest"}, the output @var{yi}
44 ## will be the closest point in the original data @var{x} to the query point
45 ## @var{xi}. When the method is @qcode{"linear"}, the output @var{yi} will
46 ## be a linear interpolation between the two closest points in the original
47 ## source data. If @var{method} is omitted or empty, it defaults to
48 ## @qcode{"linear"}.
38 ## 49 ##
39 ## The optional argument @var{options} is passed directly to Qhull when 50 ## The optional argument @var{options} is passed directly to Qhull when
40 ## computing the Delaunay triangulation used for interpolation. See 51 ## computing the Delaunay triangulation used for interpolation. See
41 ## @code{delaunayn} for information on the defaults and how to pass different 52 ## @code{delaunayn} for information on the defaults and how to pass different
42 ## values. 53 ## values.
54 ##
55 ## Example
56 ##
57 ## @example
58 ## @group
59 ## ## Evaluate sombrero() function at irregular data points
60 ## x = 16*gallery ("uniformdata", [200,1], 1) - 8;
61 ## y = 16*gallery ("uniformdata", [200,1], 11) - 8;
62 ## z = sin (sqrt (x.^2 + y.^2)) ./ sqrt (x.^2 + y.^2);
63 ## ## Create a regular grid and interpolate data
64 ## [xi, yi] = ndgrid (linspace (-8, 8, 50));
65 ## zi = griddatan ([x, y], z, [xi(:), yi(:)]);
66 ## zi = reshape (zi, size (xi));
67 ## ## Plot results
68 ## clf ();
69 ## plot3 (x, y, z, "or");
70 ## hold on
71 ## surf (xi, yi, zi);
72 ## legend ("Original Data", "Interpolated Data");
73 ## @end group
74 ## @end example
75 ##
76 ## Programming Notes: If the input is complex the real and imaginary parts
77 ## are interpolated separately. For 2-D and 3-D data additional interpolation
78 ## methods are available by using the @code{griddata} function.
43 ## @seealso{griddata, griddata3, delaunayn} 79 ## @seealso{griddata, griddata3, delaunayn}
44 ## @end deftypefn 80 ## @end deftypefn
45 81
46 function yi = griddatan (x, y, xi, method = "linear", varargin) 82 function yi = griddatan (x, y, xi, method = "linear", varargin)
47 83
48 if (nargin < 3) 84 if (nargin < 3)
49 print_usage (); 85 print_usage ();
50 endif 86 endif
51 87
52 if (ischar (method))
53 method = tolower (method);
54 endif
55
56 [m, n] = size (x); 88 [m, n] = size (x);
57 [mi, ni] = size (xi); 89 [mi, ni] = size (xi);
58 90
59 if (n != ni || rows (y) != m || columns (y) != 1) 91 if (m < n + 1)
60 error ("griddatan: dimensional mismatch"); 92 error ("griddatan: number of points in X (rows of X) must be greater than dimensionality of data + 1 (columns of X + 1)");
93 endif
94 if (! iscolumn (y) || rows (y) != m)
95 error ("griddatan: Y must be a column vector with the same number of points (rows) as X");
96 endif
97 if (n != ni)
98 error ("griddatan: dimension of query data XI (columns) must match X");
99 endif
100
101 if (nargin > 3)
102 if (isempty (method))
103 method = "linear";
104 elseif (! ischar (method))
105 error ("griddatan: METHOD must be a string");
106 endif
107
108 method = tolower (method);
109 if (strcmp (method, "linear") || strcmp (method, "nearest"))
110 ## Do nothing, these are implemented methods
111 elseif (any (strcmp (method, {"cubic", "v4"})))
112 error ('griddatan: "%s" METHOD is available for 2-D inputs by using "griddata"', method);
113
114 elseif (strcmp (method, "natural"))
115 ## FIXME: Remove when griddata.m supports "natural" method.
116 error ('griddatan: "natural" interpolation METHOD not yet implemented');
117
118 else
119 error ('griddatan: unknown interpolation METHOD: "%s"', method);
120 endif
121
61 endif 122 endif
62 123
63 ## triangulate data 124 ## triangulate data
64 tri = delaunayn (x, varargin{:}); 125 tri = delaunayn (x, varargin{:});
65 126
66 yi = NaN (mi, 1); 127 yi = NaN (mi, 1);
67 128
68 if (strcmp (method, "nearest")) 129 if (strcmp (method, "linear"))
69 ## search index of nearest point
70 idx = dsearchn (x, tri, xi);
71 valid = ! isnan (idx);
72 yi(valid) = y(idx(valid));
73
74 elseif (strcmp (method, "linear"))
75 ## search for every point the enclosing triangle 130 ## search for every point the enclosing triangle
76 [tri_list, bary_list] = tsearchn (x, tri, xi); 131 [tri_list, bary_list] = tsearchn (x, tri, xi);
77 132
78 ## only keep the points within triangles. 133 ## only keep the points within triangles.
79 valid = ! isnan (tri_list); 134 valid = ! isnan (tri_list);
87 yi(valid) = sum (y(tri(tri_list,:)).' .* bary_list, 2); 142 yi(valid) = sum (y(tri(tri_list,:)).' .* bary_list, 2);
88 else 143 else
89 yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2); 144 yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2);
90 endif 145 endif
91 146
92 elseif (any (strcmp (method, {"cubic", "v4"}))) 147 else
93 error ('griddata: "%s" METHOD only valid for 2-D inputs using "griddata"', method); 148 ## search index of nearest point
149 idx = dsearchn (x, tri, xi);
150 valid = ! isnan (idx);
151 yi(valid) = y(idx(valid));
94 152
95 elseif (strcmp (method, "natural"))
96 ## FIXME: implement missing interpolation method 'natural' for 3-D inputs.
97 error ('griddatan: "natural" interpolation METHOD not yet implemented');
98
99 else
100 error ('griddatan: unknown interpolation METHOD: "%s"', method);
101 endif 153 endif
102 154
103 endfunction 155 endfunction
104 156
105 157
107 %! [xx,yy] = meshgrid (linspace (-1,1,32)); 159 %! [xx,yy] = meshgrid (linspace (-1,1,32));
108 %! xi = [xx(:), yy(:)]; 160 %! xi = [xx(:), yy(:)];
109 %! x = 2*rand (100,2) - 1; 161 %! x = 2*rand (100,2) - 1;
110 %! x = [x;1,1;1,-1;-1,-1;-1,1]; 162 %! x = [x;1,1;1,-1;-1,-1;-1,1];
111 %! y = sin (2 * sum (x.^2,2)); 163 %! y = sin (2 * sum (x.^2,2));
112 %! zz = griddatan (x,y,xi,"linear"); 164 %! zz = griddatan (x,y,xi, "linear");
113 %! zz2 = griddata (x(:,1),x(:,2),y,xi(:,1),xi(:,2),"linear"); 165 %! zz2 = griddata (x(:,1),x(:,2),y,xi(:,1),xi(:,2), "linear");
114 %! assert (zz, zz2, 1e-10); 166 %! assert (zz, zz2, 1e-10);
115 167
116 %!testif HAVE_QHULL 168 %!testif HAVE_QHULL
117 %! [xx,yy] = meshgrid (linspace (-1,1,32)); 169 %! [xx,yy] = meshgrid (linspace (-1,1,32));
118 %! xi = [xx(:), yy(:)]; 170 %! xi = [xx(:), yy(:)];
119 %! x = 2*rand (100,2) - 1; 171 %! x = 2*rand (100,2) - 1;
120 %! x = [x;1,1;1,-1;-1,-1;-1,1]; 172 %! x = [x;1,1;1,-1;-1,-1;-1,1];
121 %! y = sin (2*sum (x.^2,2)); 173 %! y = sin (2*sum (x.^2,2));
122 %! zz = griddatan (x,y,xi,"nearest"); 174 %! zz = griddatan (x,y,xi, "nearest");
123 %! zz2 = griddata (x(:,1),x(:,2),y,xi(:,1),xi(:,2),"nearest"); 175 %! zz2 = griddata (x(:,1),x(:,2),y,xi(:,1),xi(:,2), "nearest");
124 %! assert (zz, zz2, 1e-10); 176 %! assert (zz, zz2, 1e-10);
125 177
126 %!testif HAVE_QHULL <*56515> 178 %!testif HAVE_QHULL <*56515>
127 %! x = [ 0, 0; 1, 1; 0, 1; 1, 0 ]; 179 %! x = [ 0, 0; 1, 1; 0, 1; 1, 0 ];
128 %! y = [ 1; 2; 3; 4 ]; 180 %! y = [ 1; 2; 3; 4 ];
131 183
132 ## Test input validation 184 ## Test input validation
133 %!error griddatan () 185 %!error griddatan ()
134 %!error griddatan (1) 186 %!error griddatan (1)
135 %!error griddatan (1,2) 187 %!error griddatan (1,2)
136 %!error griddatan (1,2,3) 188 %!error <number of points in X> griddatan (1,2,3)
137 %!error <OPTIONS argument must be a string> griddatan (1,2,3,4,5) 189 %!error <Y must be a column vector> griddatan ([1;2],[3,4], 1)
138 %!error <unknown interpolation METHOD> griddatan (1,2,3,4) 190 %!error <Y must .* same number of points .* as X> griddatan ([1;2],[3;4;5], 1)
139 %!#error <"v4" METHOD only valid for 2-D inputs> 191 %!error <dimension of .* XI .* must match X> griddatan ([1;2],[3;4], [1, 2])
140 %! griddatan (ones(2,2,2), 2*ones(2,2,2), 3, "v4") 192 %!error <METHOD must be a string> griddatan ([1;2],[3;4], 1, 5)
141 %!error <"cubic" METHOD only valid for> griddatan (1, 2, 3, "cubic") 193 %!error <"v4" METHOD is available for 2-D> griddatan ([1;2],[3;4], 1, "v4")
142 %!error <"natural" .* not yet implemented> griddatan (1, 2, 3, "natural") 194 %!error <"cubic" METHOD is available> griddatan ([1;2],[3;4], 1, "cubic")
143 %!error <unknown interpolation METHOD: "foobar"> griddatan (1, 2, 3, "foobar") 195 %!error <"natural" .* not yet implemented> griddatan ([1;2],[3;4], 1, "natural")
196 %!error <unknown .* METHOD: "foobar"> griddatan ([1;2],[3;4], 1, "foobar")