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.
   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

CROSS-REFERENCE INFORMATION ^

This function calls: 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.
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

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