FFNDGRID Fast 'n' Furious N-D data gridding. CALL: [fgrid, xvec] = ffndgrid(x,f, delta, limits, aver ); fgrid = Matrix of gridded data. xvec = cellarray of gridvectors xl1:dx1:xu1 or linspace(xl1,xu1,Nx1) depending on delta. x = [x1 x2,...xD] coordinate matrix for unevenly spaced data, f. size NxD. f = f(x), vector of function values length N. delta = [dx1+i*pad, dx2 ,...,dxD] or [-Nx1+i*pad -Nx2,...,NxD], where dx1 to dxD and Nx1 to NxD defines the stepsize of grids and number of bins, respectively, in the D dimesional space. Empty gridpoints are padded with IMAG(delta(1))=pad. (default -75) limits = [xl1 xu1 xl2 ...xuN fl fu n0], contain the limits in the D-dimensional x1-x2...xN-f-plane of where to grid. The last parameter, n0, removes outliers from the data set by ignoring grid points with n0 or less observations. When n0 is negative it is treated as the percentage of the total number of data points. (default [min(x1),max(x1),min(x2),.... max(xN),min(f),max(f),0]) aver = 0 If each value of fgrid is the sum of all points falling within each cell. 1 If each value of fgrid is the average of all points falling within each cell. (default) FFNDGRID grids unevenly spaced data in vector f into a matrix fgrid. NOTE: - The vector limits can be padded with NaNs if only certain limits are desired, e g if xl1 and fl are wanted: ffndgrid(x, f, [],[.5 nan nan nan 45]) - If no output arguments are given, FFGRID will plot the gridded function with the prescribed axes using PCOLOR. Examples: N = 500;D=2; sz = [N ,D ]; x = wraylrnd(1,sz); z = ones(sz(1),1); [nc, xv] = ffndgrid(x,z,-15,[],0); % Histogram pcolor(xv{:},nc) % [XV,YV]=meshgrid(xv{:}); text(XV(:),YV(:),int2str(nc(:))) dx = [diff(xv{1}(1:2)) diff(xv{2}(1:2))]; contourf(xv{:}, nc/(N*prod(dx))) % 2-D probability density plot. colorbar colormap jet See also griddata
Create cell array. | |
Display color bar (color scale) | |
Color look-up table. | |
Display message and abort function. | |
Convert sparse matrix to full matrix. | |
Get handle to current axis. | |
Black-red-yellow-white color map. | |
Scale data and display as image. | |
Input argument name. | |
Pseudocolor (checkerboard) plot. | |
Permute array dimensions. | |
Set object properties. | |
Color shading mode. | |
Create sparse matrix. | |
Write formatted data to string. | |
Graph title. | |
X-axis label. | |
Y-axis label. |
001 function [zzgrid, xvec] = ffndgrid(x, f, delta,limits,aver); 002 %FFNDGRID Fast 'n' Furious N-D data gridding. 003 % 004 % CALL: [fgrid, xvec] = ffndgrid(x,f, delta, limits, aver ); 005 % 006 % fgrid = Matrix of gridded data. 007 % xvec = cellarray of gridvectors 008 % xl1:dx1:xu1 or linspace(xl1,xu1,Nx1) depending on delta. 009 % x = [x1 x2,...xD] coordinate matrix for unevenly spaced data, f. 010 % size NxD. 011 % f = f(x), vector of function values length N. 012 % delta = [dx1+i*pad, dx2 ,...,dxD] or [-Nx1+i*pad -Nx2,...,NxD], where 013 % dx1 to dxD and Nx1 to NxD defines the stepsize of grids and 014 % number of bins, respectively, in the D dimesional space. Empty 015 % gridpoints are padded with IMAG(delta(1))=pad. (default -75) 016 % limits = [xl1 xu1 xl2 ...xuN fl fu n0], contain the limits in the 017 % D-dimensional x1-x2...xN-f-plane of where to grid. The last 018 % parameter, n0, removes outliers from the data set by ignoring 019 % grid points with n0 or less observations. When n0 is negative 020 % it is treated as the percentage of the total number of data points. 021 % (default [min(x1),max(x1),min(x2),.... max(xN),min(f),max(f),0]) 022 % aver = 0 If each value of fgrid is the sum of all points falling 023 % within each cell. 024 % 1 If each value of fgrid is the average of all points falling 025 % within each cell. (default) 026 % 027 % FFNDGRID grids unevenly spaced data in vector f into a matrix fgrid. 028 % 029 % NOTE: - The vector limits can be padded with NaNs if only 030 % certain limits are desired, e g if xl1 and fl are wanted: 031 % 032 % ffndgrid(x, f, [],[.5 nan nan nan 45]) 033 % 034 % - If no output arguments are given, FFGRID will plot the gridded 035 % function with the prescribed axes using PCOLOR. 036 % 037 % Examples: 038 % N = 500;D=2; sz = [N ,D ]; 039 % x = wraylrnd(1,sz); z = ones(sz(1),1); 040 % [nc, xv] = ffndgrid(x,z,-15,[],0); % Histogram 041 % pcolor(xv{:},nc) % 042 % [XV,YV]=meshgrid(xv{:}); 043 % text(XV(:),YV(:),int2str(nc(:))) 044 % dx = [diff(xv{1}(1:2)) diff(xv{2}(1:2))]; 045 % contourf(xv{:}, nc/(N*prod(dx))) % 2-D probability density plot. 046 % colorbar 047 % colormap jet 048 % 049 % See also griddata 050 051 % Tested on: MatLab 4.2, 5.0, 5.1, 5.2 and 5.3. 052 % History: 053 % revised pab 02.08.2001 054 % - made it general for D dimensions + changed name from ffgrid to ffndgrid 055 % -added nargchk + examples. 056 % -updated help header to wafo-style 057 % - moved dx and dy into delta =[dx,dy] 058 % -removed call to bin 059 % modified by Per A. Brodtkorb 060 % 05.10.98 secret option: aver 061 % optionally do not take average of values for each point 062 % 12.06-98 063 % by 064 % 28.7.97, Oyvind.Breivik@gfi.uib.no. 065 % 066 % Oyvind Breivik 067 % Department of Geophysics 068 % University of Bergen 069 % NORWAY 070 071 072 073 074 error(nargchk(2,5,nargin)) 075 076 [r, c] = size(x); 077 if r==1,% Make sure x is a column vector. 078 x = x(:); 079 end 080 081 [N,D]=size(x); 082 f = f(:); 083 if length(f)==1, f = f(ones(N,1),:) ; 084 elseif length(f)~=N,error('The length of f must equal size(x,1)!'),end 085 086 % Default values 087 dx = repmat(-75,1,D); 088 pad = 0; 089 xyz = [zeros(1,2*D) min(f), max(f), 0]; 090 xyz(1:2:2*D) = min(x,[],1); 091 xyz(2:2:2*D) = max(x,[],1); 092 093 094 % take average of values for each point (default) 095 if (nargin < 5)|isempty(aver), aver = 1; end 096 if (nargin >= 4) & ~isempty(limits), 097 nlim = length(limits); 098 ind = find(~isnan(limits(1:min(7,nlim)))); 099 if any(ind), xyz(ind) = limits(ind);end 100 end 101 if nargin>=3&~isempty(delta), 102 Nd =length(delta); 103 delta = delta(1:min(Nd,D)); 104 if Nd ==1, delta = repmat(delta,1,D);end 105 ind = find(~(isnan(delta)|delta==0)); 106 if any(ind), 107 dx(ind) = real(delta(ind)); 108 pad = imag(delta(1)); 109 end 110 end 111 112 xL = xyz(1:2:2*D); 113 xU = xyz(2:2:2*D); 114 115 fL = xyz(end-2); 116 fU = xyz(end-1); 117 n0 = xyz(end); 118 119 ind = find(dx<0); 120 if any(ind), 121 if any(dx(ind)~=round(dx(ind))), 122 error('Some of Nx1,...NxD in delta are not an integer!'), 123 end 124 dx(ind) = (xU(ind)-xL(ind))./(abs(dx(ind))-1); 125 end 126 127 128 % bin data in D-dimensional-space 129 binx = round((x - xL(ones(N,1),:))./dx(ones(N,1),:)) +1; 130 131 fgsiz = ones(1,max(D,2)); 132 xvec = cell(1,D); 133 for iy=1:D, 134 xvec{iy} = xL(iy):dx(iy):xU(iy); 135 fgsiz(iy) = length(xvec{iy}); 136 end 137 if D>1 138 in = all((binx >= 1) & (binx <= fgsiz(ones(N,1),:)),2) & (fL <= f) & (f <= fU); 139 else 140 in = (binx >= 1) & (binx <= fgsiz(1)) & (fL <= f) & (f <= fU); 141 end 142 binx = binx(in,:); 143 f = f(in); 144 N = length(binx); % how many datapoints are left now? 145 146 Nf = prod(fgsiz); 147 148 if D>1, 149 fact = [1 cumprod(fgsiz(1:D-1))]; 150 binx = sum((binx-1).*fact(ones(N,1),:),2)+1; % linear index to fgrid 151 end 152 % Fast gridding 153 fgrid = sparse(binx,1,f,Nf,1);% z-coordinate 154 155 if n0~=0 | aver, 156 ngrid = sparse(binx,1,ones(N,1),Nf, 1); % no. of obs per grid cell 157 if(n0 < 0), n0 = -n0*N; end % n0 is given as percentage 158 159 if n0~=0, 160 % Remove outliers 161 fgrid(ngrid <= n0) = 0; 162 ngrid(ngrid <= n0) = 0; 163 N = full(sum(ngrid(:))); % how many datapoints are left now? 164 end 165 end 166 167 ind = full(find(fgrid)); % find nonzero values 168 169 if aver, 170 fgrid(ind) = fgrid(ind)./ngrid(ind); % Make average of values for each point 171 end 172 173 if pad~=0, 174 Nil=find(fgrid==0); 175 fgrid(Nil) = full(fgrid(Nil))+pad; % Empty grid points are set to pad 176 end 177 178 rho = length(ind)/(Nf); % density of nonzero values in the grid 179 if rho>0.4, 180 fgrid = full(fgrid); 181 end 182 if r==1, 183 fgrid = fgrid.'; 184 else 185 fgrid = reshape(fgrid,fgsiz); 186 switch D % make sure fgrid is stored in the same way as meshgrid 187 case 2, fgrid=fgrid.'; 188 case 3, fgrid=permute(fgrid,[2 1 3]); 189 end 190 end 191 192 193 194 195 if (nargout > 0) 196 zzgrid = fgrid; 197 elseif D==2,% no output, then plot 198 colormap(flipud(hot)) %colormap jet 199 if 1, 200 %figure('Position', [100 100 size(fgrid)]) 201 imagesc(xvec{:}, fgrid) 202 set(gca,'YDir','normal') 203 else 204 pcolor(xvec{:}, fgrid) 205 shading flat %interp 206 end 207 colorbar 208 xlabel(inputname(1)) 209 ylabel(inputname(1)) 210 zstr=inputname(2); 211 dum = full(size(fgrid')); 212 if (~isempty(zstr)) % all this vital information ... 213 str = sprintf('Color scale: %s, %d data points, grid: %dx%d, density: %4.2f', ... 214 inputname(3), N, dum(1), dum(2), rho); 215 else 216 str = sprintf('%d data points, grid: %dx%d, density: %4.2f', ... 217 N, dum(1), dum(2), rho); 218 end 219 title(str); 220 end 221 222 223 224
Comments or corrections to the WAFO group