[vtkusers] odd vtkImageData interpolation behavior

[David Feng]

David,

I'm not terribly familiar with the inner workings of much of VTK, but 
I'm guessing that somewhere along the line vtkVoxel interpolation is 
used because the z dimension for these data sets is non-zero.  I haven't 
delved any deeper into the code, thought, so I can't be sure.  If I were 
to guess, I'd say that the following lines are related:

  // compute the pixel coordinates and parametric coordinates of point x
  if ( this->ComputeStructuredCoordinates(x, loc, pcoords) == 0 )
    {
    return -1;
    }
  // compute interpolation weights
  this->Voxel->InterpolateFunctions(pcoords,weights);

Looking quickly at the DataDescription variable, a switch statement 
there that goes over the different types (VTK_XZ_PLANE, for example) and 
interpolates with the correct cell type might fix the problem.  
Something like this:

switch (this->DataDescription)
    {
      case VTK_EMPTY:
      default:
        return -1;
        break;
       
      case VTK_SINGLE_POINT: // cellId can only be = 0
        vtkVertex::InterpolationFunctions(pcoords,weights);     
        break;

      case VTK_X_LINE:
      case VTK_Y_LINE:   
      case VTK_Z_LINE:
        vtkLine::InterpolationFunctions(pcoords,weights);
        break;

    
    case VTK_XY_PLANE:     
    case VTK_YZ_PLANE:
    case VTK_XZ_PLANE:
      vtkPixel::InterpolationFunctions(pcoords,weights);
      break;

    case VTK_XYZ_GRID:
     vtkVoxel::InterpolationFunctions(pcoords,weights);
     break;
}

I threw this in there instead of the vtkVoxel interpolation call, and 
this improved the situation a bit, but there were still some apparent 
off-by-one issues on the borders of the interpolated images.  Is there 
someone that usually maintains this file?  If so, s/he can probably fix 
this faster than I can.  If not, I can keep digging.

Thanks,
David

David Gobbi wrote:
> Hi David,
>
> I haven't had time to look at the pertinent sections of vtkImageData.cxx in
> detail, but it certainly sounds like there is a problem.
>
> From my glances at the code, it seems that the cell type should be
> determined from the DataDescription member of vtkImageData, and
> for 2D data, interpolation _should_ be done with vtkPixel cells instead
> of with vtkVoxel cells, i.e. if the code is working the way that it should,
> then your if statements should not be required.  So the question is,
> why is it interpolating with vtkVoxel if the data is 2D?
>
> About the confusion with vtkVoxel having data values at the corners
> instead of at the center (the way that people usually think of voxels),
> that is a VTK-ism and perhaps the name "vtkVoxel" for that cell type
> is misleading.  The "data" in vtkImageData is associated with the
> points in the data, not with the cells, and in most imaging filters
> the cells are ignored completely.
>
>    David
>
>
>
>
> On Thu, Jun 19, 2008 at 3:06 PM, David Feng <dfeng at cs.unc.edu> wrote:
>   
>> I tracked down the problem in vtkImageData's FindCell.  There are actually
>> two:
>>
>> 1 - Even when a dimension is of size 1 (the Y dimension in my case), there
>> is still a voxel cell in that dimension.  The indexing code automatically
>> sets all slice indexing to 0 when it multiplies the z voxel index by
>> (dims[1]-1).  This is fixed with a simple check on the first two dimensions.
>>  If they are of size one, set them to size two.  Obviously it doesn't matter
>> if the dimensions are used in order.
>>
>> 2 - The interpolation weights that a computed before returning from FromCell
>> assume that the dimensions will be used in order.  The values need to be
>> shifted if dimensions are to be ignored.  The shifting is a bit ugly.
>> Both of these problems are fixed by inserting the following code after the
>> InterpolateFunctions function call inside of FindCell, somewhere around line
>> 769 of vtkimagedata.cxx:
>>
>> //
>>  // dimensions are used for cell indexing later.  if a dimension has only
>> one
>>  // voxel, we need to make sure it still has one cell.  Also,
>>  // the coordinate weights need to be shifted to ignore that dimension.
>>  // This is ugly.  The order of these ifs makes a difference.
>>  //
>>
>>  //
>>  // is y trivial?  if so, keep all even blocks of two weights
>>  //
>>  if (dims[1] == 1)
>>   {
>>   dims[1] = 2;
>>   weights[2] = weights[4];
>>   weights[3] = weights[5];
>>   }
>>
>>  //
>>  // is x trivial? if so, keep all even indexed weights
>>  //
>>  if (dims[0] == 1)
>>   {
>>   dims[0] = 2;
>>   weights[1] = weights[2];
>>   weights[2] = weights[4];
>>   weights[3] = weights[6];
>>   }
>>
>> I'll try to make a patch and submit it.
>>
>> David
>>
>> David Feng wrote:
>>     
>>> I was indeed missing something.
>>> I thought that VTK was using cell-centered data, which means that each
>>> voxel contains a single data point at its center.  In fact VTK puts a data
>>> point at every vertex of the voxel.  This means that there is one fewer
>>> voxel in each dimension than there are data points, so the cell indexing was
>>> in fact correct.
>>> I still don't know why my XZ image is not interpolating correctly, though.
>>>
>>> David
>>>
>>> David Feng wrote:
>>>       
>>>> I noticed some odd behavior when interpolating a vtkImageData object.
>>>>  Rather than make a 2D image and rotate it to the XZ plane, I made it a 3D
>>>> image with dimensions (x,1,y).  When I call FindCell simply using the
>>>> world-space pixel coordinates, the interpolated values were not what I
>>>> expected.  Taking a closer look at the function, I noticed the following:
>>>>
>>>> return loc[2] * (dims[0]-1)*(dims[1]-1) +
>>>>         loc[1] * (dims[0]-1) + loc[0];
>>>>
>>>> ...where loc is the interpolated image-space pixel coordinate and dims
>>>> are the dimensions of the image.  The return value is supposed to the 1D
>>>> array index within the image.  With this formula, a pixel location of
>>>> (0,0,5) and (0,0,0) have the same 1D index for an image such as mine (XZ
>>>> instead of XY, with dimensions (x,1,y)).   I haven't looked more closely,
>>>> but indexing into a 3D should actually look like this, if I'm not mistaken:
>>>>
>>>> return loc[2] * dims[0] * dims[1] +
>>>>         loc[1] * dims[0] + loc[0];
>>>>
>>>> I'm not sure why the dimensions are being decremented in the original
>>>> code.  It would seem to make the indexing off by several pixels for each row
>>>> and slice.  When I made the change, I started getting the values I expected.
>>>>  A simple way to test this is to use a vtkProbeFilter to probe an XZ image
>>>> using itself as the input sample points -- without the change, my input is
>>>> different from my output, which should not be the case.   Am I missing
>>>> something?
>>>>
>>>> Thanks,
>>>> David
>>>>
>>>>
>>>>
>>>>         
>>>       
>> _______________________________________________
>> This is the private VTK discussion list.
>> Please keep messages on-topic. Check the FAQ at:
>> http://www.vtk.org/Wiki/VTK_FAQ
>> Follow this link to subscribe/unsubscribe:
>> http://www.vtk.org/mailman/listinfo/vtkusers
>>
>>