[vtkusers] odd vtkImageData interpolation behavior
David Gobbi
david.gobbi at gmail.com
Fri Jun 20 11:55:56 EDT 2008
Good work!
David
On Fri, Jun 20, 2008 at 11:53 AM, David Feng <dfeng at cs.unc.edu> wrote:
> It only took a little more digging, thankfully. The parametric coordinates
> needed to be shifted to ignore missing dimensions. I put this at the end of
> the ComputeStructuredCoordinates function, so both FindCell and
> FindAndGetCell can benefit. I've added the diff/patch to the notes on the
> bug.
>
> 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