Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}
  }
share|improve this question
    
Better to create an array {x/h , y/h , z/h } initialized to all zero, then set to 1 each vozxel that a segment passes through. Unless of course its so big that you are memory limited. –  george2079 Jan 3 '13 at 15:23
    
Are you counting points or curves? –  belisarius Jan 3 '13 at 15:55
    
@belisarius I am counting curves. In the code written by me, First I count the points then I add them up and if they are greater than 0 then I say there is a curve here. –  cesm Jan 3 '13 at 15:59
1  
A spline curve could pass into a voxel that doesn't contain points .... –  belisarius Jan 3 '13 at 16:04
    
@belisarius You are right..Because of the distance between points that form spline curve is too small, I didn't give much possibility to this situation. But, I have to consider this possibility.. –  cesm Jan 3 '13 at 16:35
show 4 more comments

2 Answers

up vote 5 down vote accepted

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]}]]

2D voxels

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] }]

3D OP example

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]}]]
share|improve this answer
    
Could you explain the functions that you have written a bit more? Your code tries to find the voxels that the curve traverse, right? –  cesm Jan 4 '13 at 12:11
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.