How can I find the lines that go through?

Question

I have a 80 x 80 x 50 dataset. From this dataset, I get some curves. These curves are randomly seeded, by which I mean that their starting points have been chosen randomly.

Below, you can find an example of points that form a curve :

points = {
  {81.7317, 64.6811, 14.7966}, {81.7316, 64.6588, 14.8941}, {81.7294, 64.6378, 14.9918},
  {81.7228, 64.6176, 15.0895}, {81.7087, 64.5979, 15.1865}, {81.6867, 64.5785, 15.2821},
  {81.6576, 64.5595, 15.3758}, {81.6223, 64.5413, 15.4676}, {81.5825, 64.5243, 15.5577},
  {81.5395, 64.5088, 15.6467}, {81.4945, 64.4953, 15.7349}, {81.4486, 64.4841, 15.8231},
  {81.4027, 64.4756, 15.9115}, {81.3576, 64.4704, 16.0006}, {81.3128, 64.4684, 16.09},
  {81.2674, 64.4693, 16.179}, {81.2217, 64.4729, 16.2679}, {81.1761, 64.4792, 16.3567},
  {81.1306, 64.488, 16.4453}, {81.0853, 64.4992, 16.5338}, {81.0403, 64.5127, 16.622},
  {80.9954, 64.5283, 16.71}, {80.9497, 64.5457, 16.7972}, {80.9025, 64.5649, 16.8833},
  {80.8539, 64.5861, 16.9681}, {80.8043, 64.6094, 17.0517}, {80.7555, 64.6349, 17.1352},
  {80.7075, 64.6628, 17.2184}, {80.6604, 64.693, 17.3012}, {80.6143, 64.7256, 17.3837},
  {80.569, 64.7605, 17.4658}, {80.5249, 64.798, 17.5473}, {80.4818, 64.838, 17.6282}, 
  {80.4398, 64.8809, 17.7082}, {80.399, 64.9267, 17.7872}, {80.3595, 64.9757, 17.8649},
  {80.3213, 65.0279, 17.9411}, {80.2848, 65.0819, 18.017}, {80.2493, 65.1363, 18.093},
  {80.2144, 65.1902, 18.1696}, {80.1795, 65.2433, 18.2469}, {80.1444, 65.295, 18.3249},
  {80.1088, 65.3453, 18.4037}, {80.0722, 65.3938, 18.4831}, {80.0344, 65.4404, 18.5631},
  {79.9951, 65.485, 18.6435}, {79.9539, 65.527, 18.7243}, {79.9104, 65.566, 18.8055}, 
  {79.8646, 65.6019, 18.8868}, {79.8163, 65.6347, 18.968}, {79.7657, 65.6645, 19.0489},
  {79.7136, 65.6922, 19.1296}, {79.6598, 65.7174, 19.2101}, {79.6044, 65.7397, 19.2903},
  {79.5475, 65.7588, 19.3703}, {79.489, 65.7743, 19.4499}, {79.4291, 65.7861, 19.5291},
  {79.3681, 65.7943, 19.6079}, {79.306, 65.7989, 19.6862}, {79.2432, 65.8002, 19.764},
  {79.1798, 65.7984, 19.8413}, {79.116, 65.794, 19.9181}, {79.0519, 65.7873, 19.9945}, 
  {78.9878, 65.7782, 20.0708}, {78.9239, 65.7658, 20.1467}, {78.86, 65.7502, 20.222},
  {78.7959, 65.7318, 20.2965}, {78.7317, 65.7111, 20.3703}, {78.6672, 65.6885, 20.4433},
  {78.6023, 65.6644, 20.5155}, {78.5371, 65.6392,20.587}, {78.4715, 65.6133, 20.6579},
  {78.4053, 65.587, 20.7281}, {78.3387, 65.5606, 20.7979}, {78.2716, 65.5344,20.8672},
  {78.204, 65.5085, 20.9361}, {78.1358, 65.4831, 21.0048}, {78.0668, 65.4575, 21.0725},
  {77.9969, 65.4312, 21.139}, {77.926, 65.4039, 21.204}, {77.8541, 65.3753, 21.2673}, 
  {77.7812, 65.3455, 21.3289}, {77.7071, 65.3148, 21.3887}, {77.6319, 65.2835, 21.4466},
  {77.5553, 65.2517, 21.5026}, {77.4776, 65.2198, 21.5568}, {77.3985, 65.188, 21.6091},
  {77.3183, 65.1564, 21.6597}, {77.2368, 65.1252, 21.7087}, {77.1543, 65.0947, 21.7562},
  {77.0708, 65.0648, 21.8023}, {76.9863, 65.0358, 21.8473}, {76.901, 65.0074, 21.8911},
  {76.8149, 64.9795, 21.9336}, {76.7278, 64.9523, 21.9745}, {76.6398, 64.9257, 22.0138},
  {76.5509, 64.8996, 22.0515}, {76.4613, 64.8738, 22.0877}, {76.371, 64.8486, 22.1225},
  {76.2801, 64.8238, 22.1558}, {76.1885, 64.7995, 22.1879}, {76.0964, 64.7757, 22.2187},
  {76.0038, 64.7524, 22.2483}, {75.9108, 64.7293, 22.2768}, {75.8174, 64.7061, 22.304},
  {75.7237, 64.6826, 22.3299}, {75.6297, 64.6588,22.3544}, {75.5354, 64.6346, 22.3774},
  {75.441, 64.6098, 22.3988}, {75.3463, 64.5844, 22.4184}, {75.2514, 64.5581, 22.436},
  {75.1565, 64.5307, 22.4515}, {75.0616, 64.502, 22.4645}, {74.9668, 64.4717, 22.4748},
  {74.8723, 64.4401, 22.4826}
  }

Answer 1

One interpretation of you question is how to find voxels intersected by line segments, where we may assume the segments may be larger than a voxel so simply looking for points in voxles wont cut it. You may find this useful:

lineseg = Table[RandomReal[{-5, 5}, 3], {2}];
dx = {.25, .3, 1};
splitline[line_, index_, del_] := Module[{x0, x1, lenlinex},
   {x0, x1} = Sort[#[[index]] & /@ line ];
   lenlinex = (line[[2, index]] - line[[1, index]]);
   If[lenlinex != 0,
    Table[ (i del - line[[1, index]])/lenlinex   , {i, 
      Ceiling[x0/del], Floor[x1/del]}], {Floor[x0/del]}]];
voxelhit[line_, dx_] := Module[{c, ci, p},
         c = Join[{0}, 
                 Sort[Flatten[splitline[line, #, dx[[#]]] & /@ {1, 2, 3}]],
                  {1}];
        Table[
         ci = (c[[i]] + c[[i + 1]])/2  ; 
         p = line[[1]] + ci (line[[2]] - line[[1]]) ; 
         Table[ Floor[p[[i]]/dx[[i]]] , {i, 3}], {i, Length[c] - 1}]];
res = voxelhit[lineseg, dx];

res is a list of voxel indices

This is a (2D) plot for verification:

  xrange = Union[#[[1]] & /@ res][[{1, -1}]] + {0, 1};
  yrange = Union[#[[2]] & /@ res][[{1, -1}]] + {0, 1};
  grid = Graphics[{
  Table[ 
    Line[{ { i dx[[1]] , yrange[[1]]  dx[[2]] } , { i dx[[1]] , 
       yrange[[2]] dx[[2]] } }], {i, xrange[[1]], xrange[[2]]}], 
  Table[ 
    Line[{ {xrange[[1]]  dx[[1]], i dx[[2]]  } , { 
        xrange[[2]] dx[[1]] , i dx[[2]] } }], {i, yrange[[1]], 
      yrange[[2]]}]}];
Show[grid, 
  Graphics[{Rectangle[  dx[[1 ;; 2]]  #[[1 ;; 2]], 
      dx[[1 ;; 2]] ({1, 1} + #[[1 ;; 2]])] & /@ res, Red, 
   Line[#[[1 ;; 2]] & /@ lineseg]}]]

finish up by calling voxelhit for each segment of each curve. Note we never loop over the huge dx/dy/dz space so it should be much faster than your original

Here is the result for your example (points as defined in the question)

Graphics3D[{Cuboid[ {.2, .2, .2} # , {.2, .2, .2} (# + {1,1,1}) ] & /@ 
     Flatten[ 
        Table[ voxelhit[points[[i ;; i + 1]], {.2, .2, .2}] , {i, 
               Length[points] - 1}], 1] , Red, Line[points] }]

Edit, a somewhat tighter version with some explination:

lineseg = Table[RandomReal[{-5, 5}, 3], {2}];
dx = {.25, 1, 1};

splitline takes a single line segment and considers a single family of planes based on index(=1,2,or 3) with spacing del. Use a parametric description of the line segment, ie x=p1 + c (p2-p1) find all c ( 0 lt c lt 1 ) such that x[[index]] == i*del for any integr i.

Note this is generalzed for non-cubic voxels, ie dx is differenc in each direction. For completenss you might want the include an offset to the voxel grid origin.

splitline[line_, index_, del_] := Module[{x0, x1, lenlinex},
   {x0, x1} = Sort[#[[index]] & /@ line ];
   lenlinex = (line[[2, index]] - line[[1, index]]);
   If[lenlinex != 0,
     (Range[Ceiling[x0/del], Floor[x1/del]] del - line[[1, index]])/lenlinex,
     {Floor[x0/del]}]];

voxelhit assembles the results of splitline called for each direction, ie we slice up the line everywhere it breaks any of the three sets of planes. Now you have the line sliced up into segments each entirely within a voxel. Return the list of voxel indices.

voxelhit[line_, dx_] := Module[{lastc},
   lastc = 0;
   Reap[( 
     Sow[Floor[(line[[1]] + ((lastc + #)/2) (line[[2]] - 
            line[[1]]))/dx]]; lastc = #) & /@ Append[
            Union[Flatten[splitline[line, #, dx[[#]]] & /@ {1, 2, 3}]],1]][[2, 1]]];
res = voxelhit[lineseg, dx];

tighter version of the 2d plot

xrange = Union[#[[1]] & /@ res][[{1, -1}]] + {-1, 2};
yrange = Union[#[[2]] & /@ res][[{1, -1}]] + {-1, 2};
grid = Graphics[{
    Line[dx[[1 ;; 2]] #  & /@ Function[y, {#, y}] /@ yrange & /@ Range[xrange /. List -> Sequence]],
    Line[dx[[1 ;; 2]] #  & /@ Function[x, {x, #}] /@ xrange & /@ Range[yrange /. List -> Sequence]]}];
rectangle[lst_] := Rectangle[lst[[1]], lst[[2]]];
Show[grid, 
     Graphics[{rectangle[ Function[v, dx[[1 ;; 2]] (#[[1 ;; 2]] + v)] /@  {0, 1}]  & /@ 
    res, Red, Line[#[[1 ;; 2]] & /@ lineseg]}]]

Answer 2

Perhaps this is similar and faster than your code:

numcurves = 10;
curves = RandomReal[{1, 100}, {numcurves, 115, 3}];
resolution = 1/.2;
k = curves resolution;
Tally[Flatten[((DeleteDuplicates@Round@#) & /@ k)/resolution, 1]]

But beware! As I wrote in a comment, you may have a spline curve passing through a voxel that isn't occupied by any point.