[vtkusers] odd vtkImageData interpolation behavior

[David Gobbi]

Hi David,

If you add this to the bugtracker, http://www.vtk.org/Bug/view_all_bug_page.php,
then it will be assigned to someone, but it's best if you keep digging yourself
if you need the fix soon.   The guts of vtkImageData are very old, so it isn't
likely that anyone is fully up-to-speed on that code.

   David


On Fri, Jun 20, 2008 at 10:31 AM, David Feng <dfeng at cs.unc.edu> wrote:
> 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
>>>
>>>
>
>
>