Binning and averaging huge pointcloud in Matlab

Question

I have a huge 3D pointcloud, [3x40e6]. It's a pointcloud of a flat surface, so it's more like 2.5D in that sense. I would like to bin the pointcloud into a fixed mesh range so that i can just put it into a 2D matrix and view it with imagesc(mymap).

I did solve this problem, but it takes too long. This is what i got so far. It works fine but takes 10 minutes. It interpolates all kind of jive and i don't need that.

xlin=linspace(min(xx),max(xx),meshsz*dxxyy);
ylin=linspace(min(yy),max(yy),meshsz);
[X,Y]=meshgrid(xlin,ylin);
disp('+ Flattening (X,Y,Z) information into 2D (X,Y)(Z) mesh..')
%Fit to 2D grid (takes a long time)
Z=griddata(xx,yy,zz,X,Y);%,'cubic');

What i would really want to do is bin all the data from my [3xN] vector into my specific 2D range map. I guess i would have to bin all the values into that map (there are some excellent and fast binning algorithms) but i would also need the specific Z-coordinate scalars in there because that's what i want to average in each bin.

Thanks!

Answer 1

Solved. Takes 200ms per million points.

yidx=[min(yy):dy:max(yy)];
xidx=[min(xx):dx:max(xx)];
ZmapSum=zeros(length(yidx),length(xidx));
ZmapIdx=zeros(size(ZmapSum));

[nx,binx] = histc(xx,xidx);
[ny,biny] = histc(yy,yidx);
%bin==0 means the value is out of range
binx=binx+1; biny=biny+1;
%binzero=( (binx==0) | (biny==0) );
%binx(binzero) = [];
%biny(binzero) = [];
%xx(binzero) = [];
%yy(binzero) = [];
%zz(binzero) = [];

%binx and biny give their respective bin locations
for i=1:1:length(xx)
    ZmapSum(biny(i),binx(i))=ZmapSum(biny(i),binx(i))+zz(i);
    ZmapIdx(biny(i),binx(i))=ZmapIdx(biny(i),binx(i))+1;
end

Zmap=ZmapSum./ZmapIdx;

Answer 2

Here's how you can do it in one step. accumarray allows you to easily swap the function you use for combining the data, so that you ca use e.g. std to see the local variability, or numel to see the counts in each bin.

%# transform your x,y coordinates to pixel (=bin) values
minX = min(xx);
maxX = max(xx);
minY = min(yy);
maxY = max(yy);

targetSize = [512 512];

xxBin = round( (xx-minX)/(maxX-minX)*(targetSize(1)-1) ) +1;
yyBin = round( (yy-minY)/(maxY-minY)*(targetSize(2)-1) ) +1;

%# map by using accumarray, take the mean of each bin 
map = accumarray([xxBin(:),yyBin(:)],zz,targetSize,@mean,0);