from qgriddata by Tim Hattrell
Quick data gridding for unevenly distributed data sets

qgriddata(xmin,xmax,dx,ymin,ymax,dy,zmin,zmax,dz,x,y,z,v);
function [xg,yg,zg,vg] = qgriddata(xmin,xmax,dx,ymin,ymax,dy,zmin,zmax,dz,x,y,z,v);
% [XG,YG,ZG,VG] = QGRIDDATA(XMIN,XMAX,DX,YMIN,YMAX,DY,ZMIN,ZMAX,DZ,X,Y,Z,V)
%
%  Quick data gridding.
%
%  This function takes unevenly distributed data in vectors X,Y,Z,V and
%  returns data in arrays XG,YG,ZG,VG which are suitable for use in Matlab
%  functions such as SLICE.
%
%  For any valid I, X(I),Y(I),Z(I),V(I) defines the x,y and z coordinates
%  of data with a scalar value V(I) associated with that coordinate. The
%  arguments XMIN,XMAX,DX... define the minimum, maximum and stepsize to
%  resample data on for a particular axis.
%
%  This function is much faster than using GRIDDATA for large data sets as
%  it does not attempt any interpolation. Instead, the function searches
%  for data points which fall in a particular cell of dimensions DX,DY,DZ.
%  The mean of all the points in the cell is returned in VG. If there are
%  no data points inside that volume then NaN is returned, i.e. THERE IS NO
%  INTERPOLATION. For this reason, if DX,DY or DZ are set too small the 
%  function will not be effective.
%
%  Tim Hattrell March 2006

[xg,yg,zg] = meshgrid(xmin:dx:xmax,ymin:dy:ymax,zmin:dz:zmax); % Make grid for data
vg = NaN(size(xg)); % Allocate space for vg

% Initialise counters
i = 1;
j = 1;
k = 1;

h = waitbar(0,'Processing...'); % Progress indicator

for i = 1:size(xg,1)
    for j = 1:size(xg,2)
        for k = 1:size(xg,3)
            xc = xg(i,j,k); % Current x,y and z coordinates
            yc = yg(i,j,k);
            zc = zg(i,j,k);
            mask = find( (abs(x-xc) <= dx) & (abs(y-yc) <= dy) & (abs(z-zc) <= dz) ); % Find data points in current volume
            if isempty(mask)
                vg(i,j,k) = NaN; % If no data points then NaN
            else
                vg(i,j,k) = mean(v(mask)); % If data points present then set to mean
            end
        end
    end
    waitbar(i/size(xg,1),h) % Update progress indicator
end

close(h) % Close progress indicator

Contact us at files@mathworks.com