[vtkusers] vtkContourFilter rounding issue with vtkIntArray
Bengt Rosenberger
bengt at ctech.com
Thu Jul 23 14:34:32 EDT 2015
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
>
> 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
> <mailto: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 <mailto: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 <http://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.
>>
>>
>>
>> _______________________________________________
>> Powered bywww.kitware.com <http://www.kitware.com>
>>
>> Visit other Kitware open-source projects athttp://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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20150723/6b8b34cb/attachment.html>
More information about the vtkusers
mailing list