[vtkusers] odd vtkImageData interpolation behavior

David Feng dfeng at cs.unc.edu
Fri Jun 20 10:46:28 EDT 2008


I've added it:
http://www.vtk.org/Bug/view.php?id=7213

I'll keep updating it as time permits.  Thanks for your help.

David

David Gobbi wrote:
> 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
>>>>
>>>>
>>>>         
>>
>>     





More information about the vtkusers mailing list