Skip to content

Fix wrong cell area / convexity computed by vtkPolygon

Tomas requested to merge tomj/vtk:Fix-vktPolygon.ComputeArea/IsConvex into master

As was pointed out by this stackoverflow question: https://stackoverflow.com/questions/49187430/strange-vtk-results-on-the-computation-of-surface-cell-areas-using-python/49194724#49194724 , the vtkPolygon::ComputeArea() method can give wrong results and I believe it may in some cases actually cause an index out of array bounds error. The reason is that if the cell is part of some larger polydata, the method for computing area uses the "global" point indices (i.e. pointing into the polydata's points array) to index the "local" points array (i.e. array with points only of the cell). This patch fixes that. A more detailed description, if required, follows:

Here is a small example to reproduce the issue, using the cone from vtkConeSource as test data (but anything else that contains at least one vtkPolygon should give similar results):

#include <string>
#include <iostream>
#include "vtkConeSource.h"
#include "vtkPolygon.h"
#include "vtkTriangle.h"

int main()
{
	vtkSmartPointer<vtkConeSource> coneSource = vtkSmartPointer<vtkConeSource>::New();
	//coneSource->SetResolution(40);
	coneSource->Update();
	vtkSmartPointer<vtkPolyData> cone = coneSource->GetOutput();

	std::cout << "# of points total: " << cone->GetPoints()->GetNumberOfPoints() << std::endl;
	double normal[3], point1[3], point2[3];
	for (vtkIdType i = 0; i < cone->GetNumberOfCells(); i++)
	{
		if (std::string(cone->GetCell(i)->GetClassName()).compare(std::string("vtkPolygon")) == 0) // only vtkPolygon seems to be wrong, vtkTriangle is OK
		{
			std::cout << "# of points in cell: " << cone->GetCell(i)->GetNumberOfPoints() << std::endl;
			for (vtkIdType j = 0; j < cone->GetCell(i)->GetPoints()->GetNumberOfPoints(); j++)
			{
				vtkIdType id = cone->GetCell(i)->GetPointIds()->GetId(j);
				cone->GetPoints()->GetPoint(id, point1);
				cone->GetCell(i)->GetPoints()->GetPoint(id, point2);
				std::cout << "ID " << id << " in global: " << point1[0] << " " << point1[1] << " " << point1[2] << std::endl;
				std::cout << "ID " << id << " in local: " << point2[0] << " " << point2[1] << " " << point2[2] << std::endl;
			}
			std::cout << "using cone->GetPoints(): " << vtkPolygon::ComputeArea(cone->GetPoints(),
				cone->GetCell(i)->GetNumberOfPoints(), cone->GetCell(i)->GetPointIds()->GetPointer(0), normal);
			std::cout << "; using ComputeArea(): " << vtkPolygon::SafeDownCast(cone->GetCell(i))->ComputeArea() << std::endl << std::endl;
		}
	}
	return 0;
}

Only the first cell of the cone is vtkPolygon, so there will be just one printout. The last line contains areas computed in two different ways: first computed using the static vtkPolygon::ComputeArea(vtkPoints *points, vtkIdType noOfPoints, vtkIdType *pointIds, double *normal) the second is computed by casting the cell to vtkPolygon and using the vtkPolygon::ComputeArea() member method. You should see different results (with the first one I believe being the correct one). The reason is that in the first approach, I am providing cone->GetPoints() as the first parameter, i.e. an array that contains all the points in the polydata, whereas in the second, the implementation of vtkPolygon::ComputeArea() currently looks like this:

double vtkPolygon::ComputeArea()
{
  double normal[3]; //not used, but required for the
                    //following ComputeArea call
  return vtkPolygon::ComputeArea(this->GetPoints(),
                             this->GetNumberOfPoints(),
                             this->GetPointIds()->GetPointer(0),
                             normal);

}

So it calls the same method, but using this->GetPoints() as the first parameter, which would be equivalent to replacing the cone->GetPoints() in the first approach by cone->GetCell(i)->GetPoints(), i.e. an array that contains only points that are in the cell. The problem is, the PointIds array of the cell contains indices into the "global" point array of the polydata, not the "local" of the cell. So when the local array is indexed by the "global" indices, the results are obviously wrong (the sample program prints out the point coordinates that you get by using the indices to index both arrays to see the difference).

I believe that if the "local" array is used, the used indices should simply be 0 to N-1 (number of points in the cell) - note that for example vtkTriangle does it this way - and if the "global" array is used, the PointIds array of the cell should be used. However, current implementation of vtkPolygon::ComputeArea() mangles the two together - using "local" array and "global" indices.

Hence I propose changing the PointIds parameter to nullptr to ensure a 0 to N-1 indexing is used.

This solution does not seem to fail any of the current tests, though maybe the TestPolygon etc. tests could be modified to test for this possibility. I made a search for ComputeArea through the whole code base, but vtkPolygon seems to be the only class doing it this "wrong" way.

Also, I noticed that the vtkPolygon::IsConvex() method, is doing the same mistake, so I made the change there too.

Merge request reports