How to find number of points in 3d gridded data

I have a large dataset of 3D points (x,y,z) and need to divide them into a regular grid of voxel and to count the number of points that fall into each voxel. Because the dataset is very large I need to use some efficient built-in function and avoid to iterate along each cell comparirg the coordinates.

How can I create a 3D grid of voxels? How can I compute the number of points for each cell?

Any help it is appriciate. Thanks in advance

Answer 1

[saumya jayasundara]

prompt='What is the threshould value for cell size?'; T=input(prompt);

 load('Point2.mat');
 [m,n] = size(Point);
 count=1;
 Xmax=Point(count,1);
 Xmin=Point(count,1);
 Ymax=Point(count,2);
 Ymin=Point(count,2);
 while count< m
    count=count+1;
    if Xmax<Point(count,1)
        Xmax=Point(count,1);
    end
    if Xmin>Point(count,1)
        Xmin=Point(count,1);
    end
    if Ymax<Point(count,2)
        Ymax=Point(count,2);
    end
    if Ymin>Point(count,2)
        Ymin=Point(count,2);
    end
 end

 Xmax;
 Xmin;
 Ymax;
 Ymin;

 corner1=[Xmin,Ymin];
 corner2=[Xmax,Ymin];
 corner3=[Xmax,Ymax];
 corner4=[Xmin,Ymax];

 cellx = floor(((Xmax-Xmin)/T)+T);
 celly = floor(((Ymax-Ymin)/T)+T);
 cell= zeros(cellx,celly);

 cellcount=1;
 for ycount = 1:celly
     for pcount = 1:cellx
         cell(pcount,ycount)= cellcount;

         cellcount=cellcount+1;
     end

 end
 attribute1=zeros(m,2);
 Point=[Point(:,1:3) attribute1];

 'Devide Grid'

 %for cellcount2= 1:cellcount
 for xcount= 1:m
  if Point(xcount,4)==0
     for cellxc =1:cellx
         l=cellxc*T;
         if(Xmin+l-T)<=Point(xcount,1)
             if Point(xcount,1)<(Xmin+l)
                 for cellyc= 1:celly
                     k=cellyc*T;
                     if (Ymin+k-T)<=Point(xcount,2)
                         if Point(xcount,2)<(Ymin+k)
                                Point(xcount,4)=cell(cellxc,cellyc);
                                break
                         else
                             continue
                         end
                     else
                         continue
                     end

                 end

             else

                continue
             end
         else

             continue
         end
     end
  else
     continue
  end
 end
 'Devided'
 %Point(1:100,:)
 %file=fopen('point.txt','w');
 %fprintf(file,Point);

 %sum=0;
 %sumdiv=0;
 %for cellcnt= 1:cellcount
 %   for Pcount2=1:m
  %      if Point(Pcount2,4)==cellcnt
   %       sum=sum+Point(Pcount2,3);
    %      sumdiv=sumdiv+1;
     %   else
      %      continue
       % end
    %end
    %avg=sum/sumdiv;
    %for Pcount2=1:m
     %   if Point(Pcount2,4)==cellcnt
      %    Point(Pcount2,5)=avg;
       % else
        %    continue
        %end
    %end
 %end