[vtk-developers] Add GetPoint function to vtkStructuredGrid

[David Gobbi]

On Fri, Feb 5, 2010 at 3:23 PM, David Doria <daviddoria+vtk at gmail.com> wrote:
> On Fri, Feb 5, 2010 at 5:15 PM, David Gobbi <david.gobbi at gmail.com> wrote:
>>
>> The extents can become negative as a result of vtkImagePad, for
>> instance.  The only way to use it to pad on the bottom or the left is
>> by setting the OutputWholeExtent so that it starts at a negative
>> number.  Yikes!  Also any filters that do "extent translation" can
>> lead to negative extents.
>>
>> The tests don't cover negative extents, so when they occur, it's
>> guaranteed that some downstream filter is gonna choke on them because
>> it "assumes" that extents are always positive.
>>
>>   David
>
> Ah yes yes, extents can be negative - how about this then?
>
> void vtkStructuredGrid::GetPoint(int x, int y, int z, double p[3])
> {
>   int extent[6];
>   this->GetExtent(extent);
>   if(x < extent[0] || x > extent[1] ||
>      y < extent[2] || y > extent[3] ||
>      z < extent[4] || z > extent[5])
>     {
>     return; // out of bounds!
>     }
>
>   //get relative coordinates to the "lower left" grid corner
>   int xrel = x - extent[0];
>   int yrel = y - extent[2];
>   int zrel = z - extent[4];
>   int ind = zrel*dims[1]*dims[0] + yrel*dims[0] + xrel;
>
>   grid->GetPoint(ind, p);
> }
> Thanks,
>
> David

Close but still no cigar.  Don't use "int", use vtkIdType instead.
Otherwise your code will not be 64-bit safe.  In fact, to make it
64-bit safe you need to do this:

   vtkIdType xrel = x - extent[0];
   vtkIdType yrel = y - extent[2];
   vtkIdType zrel = z - extent[4];
   vtkIdType ind = zrel*dims[1]*dims[0] + yrel*dims[0] + xrel;

Unless xrel, yrel, and zrel are also vtkIdType, the mathematical
computation of "zrel*dims[1]*dims[0] + yrel*dims[0] + xrel" will be
done in 32-bit regardless the type of the variable it gets stored in.

This is all kind of moot, though, because you are kinda reinventing
the wheel here.  After you do the bounds check you should just call
the following function that will do the rest for you:

vtkStructuredData::ComputePointIdForExtent(int extent[6], int ijk[3]);

It already does everything but the bounds check.

    David