[vtkusers] vtkContourFilter rounding issue with vtkIntArray

[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> 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
>
> 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/bce06750/attachment.html>