[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