[vtkusers] odd vtkImageData interpolation behavior
David Feng
dfeng at cs.unc.edu
Fri Jun 20 11:53:41 EDT 2008
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