[vtkusers] vtkContourFilter rounding issue with vtkIntArray

[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.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20150723/a4e03e44/attachment.html>