Home > wafo > kdetools > ffndgrid.m

# ffndgrid

## PURPOSE

Fast 'n' Furious N-D data gridding.

## SYNOPSIS

[zzgrid, xvec] = ffndgrid(x, f, delta,limits,aver);

## DESCRIPTION

``` 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.
dx1 to dxD and Nx1 to NxD defines the stepsize of grids and
number of bins, respectively, in the D dimesional space. Empty
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

## CROSS-REFERENCE INFORMATION

This function calls:
 cell Create cell array. colorbar Display color bar (color scale) colormap Color look-up table. error Display message and abort function. full Convert sparse matrix to full matrix. gca Get handle to current axis. hot Black-red-yellow-white color map. imagesc Scale data and display as image. inputname Input argument name. pcolor Pseudocolor (checkerboard) plot. permute Permute array dimensions. set Set object properties. shading Color shading mode. sparse Create sparse matrix. sprintf Write formatted data to string. title Graph title. xlabel X-axis label. ylabel Y-axis label.
This function is called by:

## SOURCE CODE

```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.
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
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 %
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);
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));
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
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)
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```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003