how can i "gridify" a cube and find distance between the cubes that have points inside them? thank you!

I have a cube. in the cube there are thousands of points and a randomly oriented line. I broke down the cube into smaller cubes using an incremented distance in the x y and z directions then found the intersection points of the line and the planes of the cube. eventually I want to use the smaller cubes as a generalization for the points so the line has something to hit. when the line does hit a smaller cube with points in it i need to find the length of the line inside the smaller cube and sum these up. I have the coordinates of the intersection points and the cubes that have points inside them i am just having a hard time relating the two. any ideas would be appreciated. Thank you!

    cut = 10;
  cubes = cut^3;

k1=0;
for x1 = 1:cut
    Xline(x1) = xlo + (xlen/cut) * (x1);  % breaking the cube up in the x direction 

      for y1 = 1:cut
          Yline(y1) = ylo + (ylen/cut) * (y1); % y direction 

          for z1 = 1:cut
              Zline(z1) = zlo + (zlen/cut) * (z1); % x direction 

              %finding if there is a point in the cube
              [rowxyz] = find(Xline(x1)-xlen/cut < x & x < Xline(x1)& Yline(y1)-ylen/cut < y & y < Yline(y1)& Zline(z1)-zlen/cut < z & z < Zline(z1)); 

              if (rowxyz ~= 0)  
              k1 = k1 +1;
              row1xyz{k1} = {rowxyz,x1,y1,z1};  
              xlinmin(k1) = Xline(x1)-xlen/cut;  % saved the xtremes of each smaller box with points inside to compare to the intersection points
              xlinmax(k1) = Xline(x1);
              ylinmin(k1) = Yline(y1)-ylen/cut;
              ylinmax(k1) = Yline(y1);
              zlinmin(k1) = Zline(z1)-zlen/cut;
              zlinmax(k1) = Zline(z1);
              box{x1,y1,z1} = 1;     % saves a 3d matrix of 1's and 0's (1 means there is a point in the box)
              else box{x1,y1,z1} = 0; % 0 means no point inside box

              end

          end
      end
  end

for k = 1:cut % finds intersections along each x plane
nx = [1,0,0];
V0x = [Xline(k),0,0];
[xint,checkx ]= plane_line_intersect(nx,V0x,plane1pt,plane2pt);
if checkx == 1 
    xint1 = xint;
else xint1 = [0,0,0];

end
ny = [0,1,0];        %finds intersection along each y plane
V0y = [0,Yline(k),0];
[yint,checky] = plane_line_intersect(ny,V0y,plane1pt,plane2pt);
if checky == 1 
    yint1 = yint;
else yint1 = [0,0,0];

end
nz = [0,0,1];    % finds intersection along each z plane
V0z = [0,0,Zline(k)];
[zint,checkz] = plane_line_intersect(nz,V0z,plane1pt,plane2pt);
if checkz == 1 
    zint1 = zint;
else zint1 = [0,0,0];

end
allint{k} = [xint1;yint1;zint1];
end

allint = cell2mat(allint');
allint( ~any(allint,2), : ) = []; % matrix with all intersection points