[vtkusers] vtkContourFilter rounding issue with vtkIntArray
Cory Quammen
cory.quammen at kitware.com
Thu Jul 23 13:10:51 EDT 2015
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> 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>
>> 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>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 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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20150723/2667260d/attachment.html>
More information about the vtkusers
mailing list