Mercurial > hg > octave-nkf
annotate scripts/geometry/griddata.m @ 10793:be55736a0783
Grammarcheck the documentation from m-files.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sun, 18 Jul 2010 20:35:16 -0700 |
parents | 95c3e38098bf |
children | aa40bdbfa478 |
rev | line source |
---|---|
8920 | 1 ## Copyright (C) 1999, 2000, 2007, 2008, 2009 Kai Habel |
6823 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
6823 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
6823 | 18 |
19 ## -*- texinfo -*- | |
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
10549
diff
changeset
|
20 ## @deftypefn {Function File} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) |
6823 | 21 ## @deftypefnx {Function File} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) |
22 ## | |
23 ## Generate a regular mesh from irregular data using interpolation. | |
24 ## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}. | |
25 ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If | |
26 ## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh. | |
27 ## | |
6826 | 28 ## The interpolation method can be @code{"nearest"}, @code{"cubic"} or |
29 ## @code{"linear"}. If method is omitted it defaults to @code{"linear"}. | |
6823 | 30 ## @seealso{delaunay} |
31 ## @end deftypefn | |
32 | |
10549 | 33 ## Author: Kai Habel <kai.habel@gmx.de> |
6823 | 34 ## Adapted-by: Alexander Barth <barth.alexander@gmail.com> |
35 ## xi and yi are not "meshgridded" if both are vectors | |
36 ## of the same size (for compatibility) | |
37 | |
6826 | 38 function [rx, ry, rz] = griddata (x, y, z, xi, yi, method) |
10549 | 39 |
6823 | 40 if (nargin == 5) |
6826 | 41 method = "linear"; |
6823 | 42 endif |
43 if (nargin < 5 || nargin > 7) | |
6826 | 44 print_usage (); |
6823 | 45 endif |
46 | |
6826 | 47 if (ischar (method)) |
48 method = tolower (method); | |
6823 | 49 endif |
6826 | 50 if (! all (size (x) == size (y) & size (x) == size (z))) |
51 error ("griddata: x, y, and z must be vectors of same length"); | |
6823 | 52 endif |
53 | |
9677
8cf522ce9c4d
fix griddata with vectors
Jaroslav Hajek <highegg@gmail.com>
parents:
9501
diff
changeset
|
54 ## Meshgrid xi and yi if they are a row and column vector. |
8cf522ce9c4d
fix griddata with vectors
Jaroslav Hajek <highegg@gmail.com>
parents:
9501
diff
changeset
|
55 if (rows (xi) == 1 && columns (yi) == 1) |
6826 | 56 [xi, yi] = meshgrid (xi, yi); |
6823 | 57 endif |
58 | |
9677
8cf522ce9c4d
fix griddata with vectors
Jaroslav Hajek <highegg@gmail.com>
parents:
9501
diff
changeset
|
59 if (! size_equal (xi, yi)) |
6826 | 60 error ("griddata: xi and yi must be vectors or matrices of same size"); |
6823 | 61 endif |
62 | |
6826 | 63 [nr, nc] = size (xi); |
6823 | 64 |
9501 | 65 x = x(:); |
66 y = y(:); | |
67 z = z(:); | |
9495 | 68 |
8507 | 69 ## Triangulate data. |
6826 | 70 tri = delaunay (x, y); |
10548
479536c5bb10
Replace lowercase nan with NaN for visual cue in scripts
Rik <code@nomad.inbox5.com>
parents:
9677
diff
changeset
|
71 zi = NaN (size (xi)); |
6823 | 72 |
6826 | 73 if (strcmp (method, "cubic")) |
8664 | 74 error ("griddata: cubic interpolation not yet implemented"); |
6823 | 75 |
6826 | 76 elseif (strcmp (method, "nearest")) |
8507 | 77 ## Search index of nearest point. |
6826 | 78 idx = dsearch (x, y, tri, xi, yi); |
79 valid = !isnan (idx); | |
6823 | 80 zi(valid) = z(idx(valid)); |
81 | |
6826 | 82 elseif (strcmp (method, "linear")) |
8507 | 83 ## Search for every point the enclosing triangle. |
84 tri_list = tsearch (x, y, tri, xi(:), yi(:)); | |
6823 | 85 |
8507 | 86 ## Only keep the points within triangles. |
9496 | 87 valid = !isnan (tri_list); |
88 tri_list = tri_list(valid); | |
6826 | 89 nr_t = rows (tri_list); |
6823 | 90 |
9496 | 91 tri = tri(tri_list,:); |
92 | |
8507 | 93 ## Assign x,y,z for each point of triangle. |
9496 | 94 x1 = x(tri(:,1)); |
95 x2 = x(tri(:,2)); | |
96 x3 = x(tri(:,3)); | |
6826 | 97 |
9496 | 98 y1 = y(tri(:,1)); |
99 y2 = y(tri(:,2)); | |
100 y3 = y(tri(:,3)); | |
6826 | 101 |
9496 | 102 z1 = z(tri(:,1)); |
103 z2 = z(tri(:,2)); | |
104 z3 = z(tri(:,3)); | |
6823 | 105 |
8507 | 106 ## Calculate norm vector. |
6826 | 107 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); |
9496 | 108 ## Normalize. |
109 N = diag (norm (N, "rows")) \ N; | |
6823 | 110 |
8507 | 111 ## Calculate D of plane equation |
112 ## Ax+By+Cz+D = 0; | |
6823 | 113 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); |
114 | |
8507 | 115 ## Calculate zi by solving plane equation for xi, yi. |
9677
8cf522ce9c4d
fix griddata with vectors
Jaroslav Hajek <highegg@gmail.com>
parents:
9501
diff
changeset
|
116 zi(valid) = -(N(:,1).*xi(:)(valid) + N(:,2).*yi(:)(valid) + D) ./ N(:,3); |
6823 | 117 |
118 else | |
6826 | 119 error ("griddata: unknown interpolation method"); |
6823 | 120 endif |
121 | |
6826 | 122 if (nargout == 3) |
6823 | 123 rx = xi; |
124 ry = yi; | |
125 rz = zi; | |
6826 | 126 elseif (nargout == 1) |
6823 | 127 rx = zi; |
6826 | 128 elseif (nargout == 0) |
129 mesh (xi, yi, zi); | |
6823 | 130 endif |
131 endfunction | |
132 | |
8153
ec0a13863eb7
Only run tests that depend on HDF5 and QHull if Octave was actually
Soren Hauberg <hauberg@gmail.com>
parents:
7017
diff
changeset
|
133 %!testif HAVE_QHULL |
6823 | 134 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
135 %! x = xx(:); | |
136 %! x = x + 10 * (2 * round(rand(size(x))) - 1) * eps; | |
137 %! y = yy(:); | |
138 %! y = y + 10 * (2 * round(rand(size(y))) - 1) * eps; | |
139 %! z = sin(2*(x.^2+y.^2)); | |
140 %! zz = griddata(x,y,z,xx,yy,'linear'); | |
141 %! zz2 = sin(2*(xx.^2+yy.^2)); | |
142 %! zz2(isnan(zz)) = NaN; | |
143 %! assert (zz, zz2, 100 * eps) | |
144 | |
145 %!demo | |
146 %! x=2*rand(100,1)-1; | |
147 %! y=2*rand(size(x))-1; | |
148 %! z=sin(2*(x.^2+y.^2)); | |
149 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
150 %! griddata(x,y,z,xx,yy); | |
151 %! title('nonuniform grid sampled at 100 points'); | |
152 | |
153 %!demo | |
154 %! x=2*rand(1000,1)-1; | |
155 %! y=2*rand(size(x))-1; | |
156 %! z=sin(2*(x.^2+y.^2)); | |
157 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
158 %! griddata(x,y,z,xx,yy); | |
159 %! title('nonuniform grid sampled at 1000 points'); | |
160 | |
161 %!demo | |
162 %! x=2*rand(1000,1)-1; | |
163 %! y=2*rand(size(x))-1; | |
164 %! z=sin(2*(x.^2+y.^2)); | |
165 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
166 %! griddata(x,y,z,xx,yy,'nearest'); | |
167 %! title('nonuniform grid sampled at 1000 points with nearest neighbor'); |