[vtkusers] vtkContourFilter rounding issue with vtkIntArray

Cory Quammen cory.quammen at kitware.com
Thu Jul 23 14:45:10 EDT 2015


Awesome! Thanks in advance for contributing!

Cory

On Thu, Jul 23, 2015 at 2:34 PM, Bengt Rosenberger <bengt at ctech.com> wrote:

>  Cory,
>
> I can confirm a call to vtkDataArrayRoundIfNecessary fixes the issue. I
> will prepare and submit a patch.
>
> Thank you,
> Bengt
>
>
>   Bengt,
>
>  That's what it looks like to me. Could you try to add a call to vtkDataArrayRoundIfNecessary there
> and report back if it solves your problem? If so, would you be interested
> in contributing a patch to VTK? Developer instructions are available here:
>
>  <https://gitlab.kitware.com/vtk/vtk/blob/master/CONTRIBUTING.md>
> https://gitlab.kitware.com/vtk/vtk/blob/master/CONTRIBUTING.md
>
>  Feel free to add me as a reviewer: @cory-quammen is my username.
>
>  Thanks,
> Cory
>
> On Thu, Jul 23, 2015 at 1:01 PM, Bengt Rosenberger <bengt at ctech.com>
> wrote:
>
>>  Small correction: vtkDataArrayInterpolateTuple is either called from
>> vtkDataArray, ln 496 or ln 515. I can't exactly tell because I'm unable to
>> step through that macro. However, I could see that
>> vtkDataArrayInterpolateTuple, ln 103 was called.
>>
>> Is there just a vtkDataArrayRoundIfNecessary missing?
>>
>> Bye,
>> Bengt
>>
>> Am 23.07.2015 um 18:50 schrieb Bengt Rosenberger:
>>
>> Hi Cory,
>>
>> thanks for the reply.
>>
>> When you run the provided code with the provided data, the interpolation
>> for the int data ultimately happens in vtkDataArray::InterpolateEdge, ln
>> 520, which in turn calls vtkDataArrayInterpolateTuple, ln 103. This one
>> performs a multiplication with double and then casts back to the type
>> Scalar, it does not round the final interpolant. The overload above it does
>> tho, but it's not called.
>>
>> Hope that helps,
>> Bengt
>>
>> Am 23.07.2015 um 15:30 schrieb Cory Quammen:
>>
>> Bengt,
>>
>>  From what I can tell looking at the interpolating code in
>> vtkDataArray.cxx, it looks like integer types should be rounded, and it
>> looks like this support has been in VTK for a long time. Can you point out
>> where in VTK you are seeing the casting back to integer type instead of
>> rounding? You may have discovered a bug.
>>
>>  Thanks,
>> Cory
>>
>> On Wed, Jul 22, 2015 at 6:05 PM, Bengt Rosenberger < <bengt at ctech.com>
>> bengt at ctech.com> wrote:
>>
>>> Hello everyone!
>>>
>>> Sorry if this may appear twice on the mailing list, the first time I've
>>> submitted this the attachment prevented this to show up. This time I'll
>>> provide a link to my DropBox for the data.
>>>
>>> Anyway, we encountered an issue with vtkContourFilter. Specifically, we
>>> have an unstructured grid containing the data to contour with and an
>>> additional vtkIntArray. The data in the int array is supposed to be of all
>>> the same value, in this case 3. However, after running the field through
>>> the contour filter, the int data contains values of 2 as well.
>>>
>>> After digging deeper into the VTK code, we discovered that the data for
>>> the int array is cast to double, interpolated and then cast back to int.
>>>
>>> So my question is:
>>> Why aren't the data values being rounded after interpolation, but
>>> instead just cast back to int? Is there a way to enable rounding or copying
>>> for only a single array in the grid?
>>>
>>> Thank you!
>>>
>>> Here's the data: https://dl.dropboxusercontent.com/u/16837761/test1.7z
>>> And here's the code we used for testing:
>>>
>>> #include <cstdint>
>>> #include <vtkSmartPointer.h>
>>> #include <vtkUnstructuredGrid.h>
>>> #include <vtkCellArray.h>
>>> #include <vtkContourFilter.h>
>>> #include <vtkPointData.h>
>>> #include <vtkXMLUnstructuredGridReader.h>
>>>
>>> // Checks whether the geo_layer values are all constant
>>> void check_geo_layer(vtkDataSet* ds)
>>> {
>>>     if (ds->GetNumberOfCells() == 0)
>>>     {
>>>         assert(false);
>>>     }
>>>
>>>     vtkDataArray* geoLayerData =
>>> ds->GetPointData()->GetArray("Geo_Layer");
>>>
>>>     double geoLayer =
>>> geoLayerData->GetComponent(ds->GetCell(0)->GetPointId(0), 0);
>>>     for (uint32_t i = 0; i < ds->GetNumberOfCells(); ++i)
>>>     {
>>>         vtkCell* c = ds->GetCell(i);
>>>         for (uint32_t j = 0; j < c->GetNumberOfPoints(); ++j)
>>>         {
>>>             double gl = geoLayerData->GetComponent(c->GetPointId(j), 0);
>>>             if (gl != geoLayer)
>>>             {
>>>                 assert(false);
>>>             }
>>>         }
>>>     }
>>> }
>>>
>>> template <class T>
>>> void check_geo_layer(vtkSmartPointer<T> vtkFilter)
>>> {
>>>     vtkFilter->Update();
>>>     auto result = vtkFilter->GetOutput();
>>>     vtkDataSet* ds = vtkDataSet::SafeDownCast(result);
>>>
>>>     check_geo_layer(ds);
>>> }
>>>
>>> int _tmain(int argc, _TCHAR* argv[])
>>> {
>>>     // Read the test field
>>>     // Field has one vtkDoubleArray ("TOTHC") and one vtkIntArray
>>> ("Geo_Layer")
>>>     vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
>>> vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
>>>     reader->SetFileName("C:\\test1.vtk");
>>>     reader->Update();
>>>     vtkDataSet* g1 = reader->GetOutputAsDataSet();
>>>
>>>     // Values in Geo_Layer are supposed to be all 3's. Check that now
>>> for the input.
>>>     check_geo_layer(g1);
>>>
>>>     const double isoValue = 0.98712348937988281;
>>>
>>>     vtkSmartPointer<vtkContourFilter> contour = vtkContourFilter::New();
>>>     contour->SetInputData(g1);
>>>     contour->ReleaseDataFlagOn();
>>>     contour->GenerateTrianglesOn();
>>>     contour->ComputeGradientsOff();
>>>     contour->ComputeNormalsOff();
>>>     contour->ComputeScalarsOff();
>>>     contour->SetNumberOfContours(1);
>>>     contour->SetValue(0, isoValue);
>>>     contour->SetInputArrayToProcess(0, 0, 0,
>>> vtkDataObject::FIELD_ASSOCIATION_POINTS, "TOTHC");
>>>
>>>     // Check again whether values in Geo_Layer are only 3's.
>>>     // PROBLEM: They are not. There are 2's in there as well, which
>>> occur from interpolating integers and not rounding in the end.
>>>     //            You can see the problem if you convert Geo_Layer to a
>>> vtkDoubleArray, in which case values of 2.9999999998 may appear after
>>> contouring.
>>>     check_geo_layer(contour);
>>>
>>>     std::cin.ignore();
>>> }
>>> _______________________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Please keep messages on-topic and check the VTK FAQ at:
>>> http://www.vtk.org/Wiki/VTK_FAQ
>>>
>>> Search the list archives at: <http://markmail.org/search/?q=vtkusers>
>>> http://markmail.org/search/?q=vtkusers
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://public.kitware.com/mailman/listinfo/vtkusers
>>>
>>
>>
>>
>>  --
>> Cory Quammen
>> R&D Engineer
>> Kitware, Inc.
>>
>>
>>
>>
>> _______________________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
>>
>> Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
>>
>> Search the list archives at: http://markmail.org/search/?q=vtkusers
>>
>> Follow this link to subscribe/unsubscribe:http://public.kitware.com/mailman/listinfo/vtkusers
>>
>>
>>
>
>
>  --
> Cory Quammen
> R&D Engineer
> Kitware, Inc.
>
>
>


-- 
Cory Quammen
R&D Engineer
Kitware, Inc.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20150723/7b62b845/attachment.html>


More information about the vtkusers mailing list