Mercurial > hg > octave-nkf

comparison scripts/geometry/griddata.m @ 6826:8618f29520c6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-08-24 16:02:07 by jwe]
author jwe
date Fri, 24 Aug 2007 16:02:07 +0000
parents 9fddcc586065
children 93c65f2a5668
comparison
equal deleted inserted replaced
6825:59e22e30aff8 6826:8618f29520c6
24 ## Generate a regular mesh from irregular data using interpolation. 24 ## Generate a regular mesh from irregular data using interpolation.
25 ## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}. 25 ## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}.
26 ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If 26 ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If
27 ## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh. 27 ## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh.
28 ## 28 ##
29 ## The interpolation method can be 'nearest', 'cubic' or 'linear'. 29 ## The interpolation method can be @code{"nearest"}, @code{"cubic"} or
30 ## If method is omitted it defaults to 'linear'. 30 ## @code{"linear"}. If method is omitted it defaults to @code{"linear"}.
31 ## @seealso{delaunay} 31 ## @seealso{delaunay}
32 ## @end deftypefn 32 ## @end deftypefn
33 33
34 ## Author: Kai Habel <kai.habel@gmx.de> 34 ## Author: Kai Habel <kai.habel@gmx.de>
35 ## Adapted-by: Alexander Barth <barth.alexander@gmail.com> 35 ## Adapted-by: Alexander Barth <barth.alexander@gmail.com>
36 ## xi and yi are not "meshgridded" if both are vectors 36 ## xi and yi are not "meshgridded" if both are vectors
37 ## of the same size (for compatibility) 37 ## of the same size (for compatibility)
38 38
39 function [rx, ry, rz] = griddata (x,y,z,xi,yi,method) 39 function [rx, ry, rz] = griddata (x, y, z, xi, yi, method)
40 40
41 if (nargin == 5) 41 if (nargin == 5)
42 method="linear"; 42 method = "linear";
43 endif 43 endif
44 if (nargin < 5 || nargin > 7) 44 if (nargin < 5 || nargin > 7)
45 print_usage(); 45 print_usage ();
46 endif 46 endif
47 47
48 if (ischar(method)) 48 if (ischar (method))
49 method=tolower(method); 49 method = tolower (method);
50 endif 50 endif
51 if (!all( (size(x)==size(y)) & (size(x)==size(z)))) 51 if (! all (size (x) == size (y) & size (x) == size (z)))
52 error("griddata: x,y,z must be vectors of same length"); 52 error ("griddata: x, y, and z must be vectors of same length");
53 endif 53 endif
54 54
55 ## meshgrid xi and yi if they are vectors unless they 55 ## meshgrid xi and yi if they are vectors unless they
56 ## are vectors of the same length 56 ## are vectors of the same length
57 if (isvector(xi) && isvector(yi) && numel(xi) ~= numel(yi)) 57 if (isvector (xi) && isvector (yi) && numel (xi) != numel (yi))
58 [xi,yi]=meshgrid(xi,yi); 58 [xi, yi] = meshgrid (xi, yi);
59 endif 59 endif
60 60
61 if (any(size(xi) != size(yi))) 61 if (any (size (xi) != size (yi)))
62 error("griddata: xi and yi must be vectors or matrices of same size"); 62 error ("griddata: xi and yi must be vectors or matrices of same size");
63 endif 63 endif
64 64
65 [nr,nc] = size(xi); 65 [nr, nc] = size (xi);
66 66
67 ## triangulate data 67 ## triangulate data
68 tri = delaunay(x,y); 68 tri = delaunay (x, y);
69 zi = nan(size(xi)); 69 zi = nan (size (xi));
70 70
71 if strcmp(method,"cubic") 71 if (strcmp (method, "cubic"))
72 error("griddata(...,'cubic') cubic interpolation not yet implemented\n") 72 error ("griddata: cubic interpolation not yet implemented")
73 73
74 elseif strcmp(method,'nearest') 74 elseif (strcmp (method, "nearest"))
75 ## search index of nearest point 75 ## search index of nearest point
76 idx = dsearch(x,y,tri,xi,yi); 76 idx = dsearch (x, y, tri, xi, yi);
77 valid = !isnan(idx); 77 valid = !isnan (idx);
78 zi(valid) = z(idx(valid)); 78 zi(valid) = z(idx(valid));
79 79
80 elseif strcmp(method,'linear') 80 elseif (strcmp (method, "linear"))
81 ## search for every point the enclosing triangle 81 ## search for every point the enclosing triangle
82 tri_list= tsearch(x,y,tri,xi(:),yi(:)); 82 tri_list= tsearch (x, y, tri, xi(:), yi(:));
83 83
84 ## only keep the points within triangles. 84 ## only keep the points within triangles.
85 valid = !isnan(reshape(tri_list,size(xi))); 85 valid = !isnan (reshape (tri_list, size (xi)));
86 tri_list = tri_list(!isnan(tri_list)); 86 tri_list = tri_list(!isnan (tri_list));
87 nr_t = rows(tri_list); 87 nr_t = rows (tri_list);
88 88
89 ## assign x,y,z for each point of triangle 89 ## assign x,y,z for each point of triangle
90 x1 = x(tri(tri_list,1));y1=y(tri(tri_list,1));z1=z(tri(tri_list,1)); 90 x1 = x(tri(tri_list,1));
91 x2 = x(tri(tri_list,2));y2=y(tri(tri_list,2));z2=z(tri(tri_list,2)); 91 x2 = x(tri(tri_list,2));
92 x3 = x(tri(tri_list,3));y3=y(tri(tri_list,3));z3=z(tri(tri_list,3)); 92 x3 = x(tri(tri_list,3));
93
94 y1 = y(tri(tri_list,1));
95 y2 = y(tri(tri_list,2));
96 y3 = y(tri(tri_list,3));
97
98 z1 = z(tri(tri_list,1));
99 z2 = z(tri(tri_list,2));
100 z3 = z(tri(tri_list,3));
93 101
94 ## calculate norm vector 102 ## calculate norm vector
95 N = cross([x2-x1, y2-y1, z2-z1],[x3-x1, y3-y1, z3-z1]); 103 N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]);
96 N_norm = sqrt(sumsq(N,2)); 104 N_norm = sqrt (sumsq (N, 2));
97 N = N ./ N_norm(:,[1,1,1]); 105 N = N ./ N_norm(:,[1,1,1]);
98 106
99 ## calculate D of plane equation 107 ## calculate D of plane equation
100 ## Ax+By+Cz+D=0; 108 ## Ax+By+Cz+D=0;
101 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); 109 D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1);
102 110
103 ## calculate zi by solving plane equation for xi,yi 111 ## calculate zi by solving plane equation for xi,yi
104 zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3); 112 zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3);
105 113
106 else 114 else
107 error("griddata: unknown interpolation method"); 115 error ("griddata: unknown interpolation method");
108 endif 116 endif
109 117
110 if nargout == 3 118 if (nargout == 3)
111 rx = xi; 119 rx = xi;
112 ry = yi; 120 ry = yi;
113 rz = zi; 121 rz = zi;
114 elseif nargout == 1 122 elseif (nargout == 1)
115 rx = zi; 123 rx = zi;
116 elseif nargout == 0 124 elseif (nargout == 0)
117 mesh(xi, yi, zi); 125 mesh (xi, yi, zi);
118 endif 126 endif
119 endfunction 127 endfunction
120 128
121 %!test 129 %!test
122 %! [xx,yy]=meshgrid(linspace(-1,1,32)); 130 %! [xx,yy]=meshgrid(linspace(-1,1,32));