[vtkusers] vtkContourFilter rounding issue with vtkIntArray

[Bengt Rosenberger]

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();
}