Mercurial > octave-nkf
annotate scripts/geometry/griddata.m @ 9501:3c40d81c197f
ChangeLog and style fixes
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 06 Aug 2009 14:41:24 -0400 |
| parents | 079c06f37f17 |
| children | 8cf522ce9c4d |
| 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 -*- | |
| 20 ## @deftypefn {Function File} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) | |
| 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 | |
| 33 ## Author: Kai Habel <kai.habel@gmx.de> | |
| 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) |
| 6823 | 39 |
| 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 | |
| 8507 | 54 ## Meshgrid xi and yi if they are vectors unless they |
| 55 ## are vectors of the same length. | |
| 6826 | 56 if (isvector (xi) && isvector (yi) && numel (xi) != numel (yi)) |
| 57 [xi, yi] = meshgrid (xi, yi); | |
| 6823 | 58 endif |
| 59 | |
| 6826 | 60 if (any (size (xi) != size (yi))) |
| 61 error ("griddata: xi and yi must be vectors or matrices of same size"); | |
| 6823 | 62 endif |
| 63 | |
| 6826 | 64 [nr, nc] = size (xi); |
| 6823 | 65 |
| 9501 | 66 x = x(:); |
| 67 y = y(:); | |
| 68 z = z(:); | |
| 9495 | 69 |
| 8507 | 70 ## Triangulate data. |
| 6826 | 71 tri = delaunay (x, y); |
| 72 zi = nan (size (xi)); | |
| 6823 | 73 |
| 6826 | 74 if (strcmp (method, "cubic")) |
| 8664 | 75 error ("griddata: cubic interpolation not yet implemented"); |
| 6823 | 76 |
| 6826 | 77 elseif (strcmp (method, "nearest")) |
| 8507 | 78 ## Search index of nearest point. |
| 6826 | 79 idx = dsearch (x, y, tri, xi, yi); |
| 80 valid = !isnan (idx); | |
| 6823 | 81 zi(valid) = z(idx(valid)); |
| 82 | |
| 6826 | 83 elseif (strcmp (method, "linear")) |
| 8507 | 84 ## Search for every point the enclosing triangle. |
| 85 tri_list = tsearch (x, y, tri, xi(:), yi(:)); | |
| 6823 | 86 |
| 8507 | 87 ## Only keep the points within triangles. |
| 9496 | 88 valid = !isnan (tri_list); |
| 89 tri_list = tri_list(valid); | |
| 6826 | 90 nr_t = rows (tri_list); |
| 6823 | 91 |
| 9496 | 92 tri = tri(tri_list,:); |
| 93 | |
| 8507 | 94 ## Assign x,y,z for each point of triangle. |
| 9496 | 95 x1 = x(tri(:,1)); |
| 96 x2 = x(tri(:,2)); | |
| 97 x3 = x(tri(:,3)); | |
| 6826 | 98 |
| 9496 | 99 y1 = y(tri(:,1)); |
| 100 y2 = y(tri(:,2)); | |
| 101 y3 = y(tri(:,3)); | |
| 6826 | 102 |
| 9496 | 103 z1 = z(tri(:,1)); |
| 104 z2 = z(tri(:,2)); | |
| 105 z3 = z(tri(:,3)); | |
| 6823 | 106 |
| 8507 | 107 ## Calculate norm vector. |
| 6826 | 108 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]); |
| 9496 | 109 ## Normalize. |
| 110 N = diag (norm (N, "rows")) \ N; | |
| 6823 | 111 |
| 8507 | 112 ## Calculate D of plane equation |
| 113 ## Ax+By+Cz+D = 0; | |
| 6823 | 114 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); |
| 115 | |
| 8507 | 116 ## Calculate zi by solving plane equation for xi, yi. |
| 117 zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D) ./ N(:,3); | |
| 6823 | 118 |
| 119 else | |
| 6826 | 120 error ("griddata: unknown interpolation method"); |
| 6823 | 121 endif |
| 122 | |
| 6826 | 123 if (nargout == 3) |
| 6823 | 124 rx = xi; |
| 125 ry = yi; | |
| 126 rz = zi; | |
| 6826 | 127 elseif (nargout == 1) |
| 6823 | 128 rx = zi; |
| 6826 | 129 elseif (nargout == 0) |
| 130 mesh (xi, yi, zi); | |
| 6823 | 131 endif |
| 132 endfunction | |
| 133 | |
|
8153
ec0a13863eb7
Only run tests that depend on HDF5 and QHull if Octave was actually
Soren Hauberg <hauberg@gmail.com>
parents:
7017
diff
changeset
|
134 %!testif HAVE_QHULL |
| 6823 | 135 %! [xx,yy]=meshgrid(linspace(-1,1,32)); |
| 136 %! x = xx(:); | |
| 137 %! x = x + 10 * (2 * round(rand(size(x))) - 1) * eps; | |
| 138 %! y = yy(:); | |
| 139 %! y = y + 10 * (2 * round(rand(size(y))) - 1) * eps; | |
| 140 %! z = sin(2*(x.^2+y.^2)); | |
| 141 %! zz = griddata(x,y,z,xx,yy,'linear'); | |
| 142 %! zz2 = sin(2*(xx.^2+yy.^2)); | |
| 143 %! zz2(isnan(zz)) = NaN; | |
| 144 %! assert (zz, zz2, 100 * eps) | |
| 145 | |
| 146 %!demo | |
| 147 %! x=2*rand(100,1)-1; | |
| 148 %! y=2*rand(size(x))-1; | |
| 149 %! z=sin(2*(x.^2+y.^2)); | |
| 150 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
| 151 %! griddata(x,y,z,xx,yy); | |
| 152 %! title('nonuniform grid sampled at 100 points'); | |
| 153 | |
| 154 %!demo | |
| 155 %! x=2*rand(1000,1)-1; | |
| 156 %! y=2*rand(size(x))-1; | |
| 157 %! z=sin(2*(x.^2+y.^2)); | |
| 158 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
| 159 %! griddata(x,y,z,xx,yy); | |
| 160 %! title('nonuniform grid sampled at 1000 points'); | |
| 161 | |
| 162 %!demo | |
| 163 %! x=2*rand(1000,1)-1; | |
| 164 %! y=2*rand(size(x))-1; | |
| 165 %! z=sin(2*(x.^2+y.^2)); | |
| 166 %! [xx,yy]=meshgrid(linspace(-1,1,32)); | |
| 167 %! griddata(x,y,z,xx,yy,'nearest'); | |
| 168 %! title('nonuniform grid sampled at 1000 points with nearest neighbor'); |