[vtkusers] Uniformly subdividing a 3D space with points

[David Doria]

On Mon, Jan 18, 2010 at 2:22 AM, Bryn Lloyd <blloyd at vision.ee.ethz.ch> wrote:
> David Doria wrote:
>>
>> On Sun, Jan 17, 2010 at 10:08 AM, Bryn Lloyd <blloyd at vision.ee.ethz.ch>
>> wrote:
>>>
>>> Hi David,
>>>
>>>
>>> You could try the following:
>>>
>>> 1. put your regular grid (e.g. vtkImageData) into a vtkCellLocator.
>>>
>>> 2. query for each of your points in the point cloud:
>>> vtkCellLocator::FindCell(...)
>>>
>>> 3. store this information in a way suitable for you, e.g. vtkDataArray
>>> for
>>> the point cloud.
>>>
>>>
>>> BTW. vtkDataSet (and all sub-classes) has its own "FindCell" method, so
>>> you
>>> could skip the vtkCellLocator.
>>>
>>>
>>> Cheers
>>> Bryn
>>>
>>>
>>>
>>>
>>>
>>> David Doria wrote:
>>>>
>>>> I have a point cloud. I would like to break it down into a uniform
>>>> grid of cubic voxels (axis aligned is fine). That is, I would like to
>>>> be able to ask "which points are in voxel (i,j,k)?". The only way I
>>>> know how to do this is to actually create a cube to represent each
>>>> voxel and then, for each cube, use a vtkSelectEnclosedPoints filter.
>>>> This seems very "manual" - I'd need to come up with a data structure
>>>> to store the cubes in this (i,j,k) type format, etc. It seems like
>>>> there should be a way to divide the space into a grid of cubes all in
>>>> one shot.
>>>>
>>>> Does anyone know of a filter that does this efficiently?
>>>>
>>>> Thanks,
>>>>
>>>> David
>>>> _______________________________________________
>>
>> Ok, seems like we're headed in the right direction. I wanted to make
>> sure I was using vtkCellLocator correctly, so I tried to use the
>> GenerateRepresentation function to visualize it. However, GetLevel()
>> seems to return -1?
>>
>> http://www.cmake.org/Wiki/VTK/Examples/CellLocator
>>
>> Assuming it was working, I started building the code for my original
>> question:
>> http://www.cmake.org/Wiki/VTK/Examples/GridPointCloud
>>
>> It seems to kind of work... except it always says point 0 is in cell
>> -1. There shouldn't be negative cell indices, right? So you're
>> recommending asking which cell every point is in, then storing the
>> cellId as an array in the PointData? Then I would somehow (how?) query
>> this PointData array to get all of the point indices that are in a
>> particular cell?
>>
>> Also, then I could traverse the cells in linear order, but is there a
>> way to access them in (i,j,k) notation without simply writing my own
>> (linear index)->(i,j,k) function?
>>
>> Thanks for your help,
>>
>> David
>
>
> If FindCell returns -1, then it could not find a cell which contains the
> points. The point might be outside the regular grid (image)...
>
>
> BTW. I just checked vtkImageData again. You might be interested in the
> function "ComputeStructuredCoordinates". It will compute the (i,j,k)
> coordinates of the cell from your 3D point (x,y,z).

Ah, awesome. I didn't realize ImageData created cells (exactly the
voxels I was looking for) automatically.

So now my plan was to do this:

1) for each point, compute which cell (i,j,k) it is in and store it in
an array in the PointData
2) the main loop over the voxels would then, for each voxel V,
traverse the new PointData array and collect indices of points which
belong to voxel V and operate on them

I don't have much experience with algorithm efficiency stuff, but it
seems to me like this could get slow as the number of voxels gets
large. Or maybe this array traversal+compare is not a problem?

Thanks again,

David