I have been noticing some inaccurate principal curvature estimates near peaks of some of manifolds I was working with. Example included below. Areas with issue are pointed by arrows. These inaccuracies are affecting calculation of other curvature-based metrics like curvedness and shape-index.
I am using vtk.vtkCurvatures() in python and visualizing values in paraview. The mean and Gaussian curvatures seems to be okay. To compare, I also computed the principal curvatures (i.e. min, max) from VTK-computed mean and Gaussian curvatures using numpy, using the same formula as that described in documentation. The numpy ones are saved as point-data and looks more reasonable and continuous.
As new discourse user, I am not allowed to upload attachments.
This workaround might work: I used random hill class. Here is the code snipped with fixed random seed that should generate the same surface, albeit not sure:
rh = vtk.vtkParametricRandomHills()
rh.SetRandomSeed(1256)
rhFnSrc = vtk.vtkParametricFunctionSource()
rhFnSrc.SetParametricFunction(rh)
rhFnSrc.Update()
poly_new = rhFnSrc.GetOutput() # all curvature computation done on this polydata
How do you treat the case H^2 - K < 0 with your numpy script? Because there is a square root involved. VTK catches a H^2 - K < 0 at the locations of the discontinuity points you have on your data.
My guts tell me that at these point locations, the gaussian curvature almost equals the squared mean curvature, which means that kmin ~ kmax ~ mean curvature. But in reality there should be a complex number looking like kmin = mean curvature + i * something_very_small.
It is actually worse. At those points the gaussian curvature is very big and likely holds a large error because of the discretization. What happens when you have H^2 - K < 0 is that the estimated surface is not anymore a nice manifold with points with real coordinates.
You can convince yourself by writing down those curvatures using complex polynomials. Let’s say you do the taylor expansion of a real heightfield where heights are defined by projecting the neighborhood to the local normal using complex polynomials, you can write:
f(r, theta) = a r^2 + b r^2 exp(2 i theta) + c r^ 2 exp(- 2 i theta) + o(r^2).
This f(r, theta) is real because we said so (so b and c are the conjugate of each other in this instance).
You can show by making those coefficients correspond with what you would have in cartesian coordinates, that
K = 4(a^2 - 4 bc) and H = 2a.
When you write H^2 - K, you end up with 16bc. If b and c are conjugate, as we said they were, this is a positive real number. If not, it means that your heightfield is complex, and you’re in trouble.
So basically you cannot really trust the computations of H and K in your example.
In the data set you shared, those points are not set to zero, they are equal to the mean curvature. You can see it by looking with a spreadsheet view in ParaView.
EDIT:
My bad, I was pinging the wrong point. In those instances, you have H^2 - K almost equalling zero. So I suppose numpy did the computation differently. We could fix that in VTK by adding a larger error window than there is currently.
However, in certain scenario it may be desirable to allow custom error-window (at class level) that should be treated as zero (please ignore any syntax errors ):
if std::abs(tmp) < this->errorWindow
{
tmp = 0.0;
}
This topic attracted my attention, because I’m all in for fixing numerical instabilities in VTK.
I noticed this code in ComputeGaussCurvature():
// angles
// I get lots of acos domain errors so clamp the value to +/-1 as the
// normalize function can return 1.000000001 etc (I think)
double ac1 = vtkMath::Dot(e1, e2);
double ac2 = vtkMath::Dot(e2, e0);
double ac3 = vtkMath::Dot(e0, e1);
alpha0 = acos(-CLAMP_MACRO(ac1));
alpha1 = acos(-CLAMP_MACRO(ac2));
alpha2 = acos(-CLAMP_MACRO(ac3));
The math in vtkMath::AngleBetweenVectors() will give a more precise answer than acos(), and it might be faster, too, because it doesn’t require normalized vectors.
Just to be clear: numpy is not doing anything magical. Numpy calculations in my script is treating the H^2 - K ~ 0 case differently than how VTK c++ code is treating. VTK is setting kmax/kmin to zero, my script is setting them to H (mean curvature).
