# Inaccurate principal curvature near peaks

Could you share your data set?

I tried and got “Sorry, new users can not upload attachments.”

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
``````

FYI @cbhushan you should be now able to upload files. You profile just got upgraded.

Thanks. The original file is attached.

random_hill1.vtk (365.4 KB)

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.

I am setting those points to zero. Code snippet below. I do not think those points belong to the case of H^2 - K < 0

``````    AA = np.square(curvature_np['Mean_Curvature']) - curvature_np['Gauss_Curvature']
if not np.alltrue(AA>=0):
logging.warning(f'{np.count_nonzero(AA<0)} values below-sqrt is not non-negative!')
BB = np.sqrt(AA)
BB[AA<0] = 0
curvature_np['Max_Curv_NP'] = curvature_np['Mean_Curvature'] + BB
curvature_np['Min_Curv_NP'] = curvature_np['Mean_Curvature'] - BB
``````

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.

I do think your initial guess is in the right direction, i.e. kmin ~ kmax ~ mean curvature seems to hit some precision or error-window issues.

Can you point me to the corresponding C++ code where principal curvatures are computed?

It is in `Filters/General/vtkCurvatures.cxx`, in `GetMinimumCurvature` and `GetMaximumCurvature`. It works for an error of 0.1, which is kind of large…

What I propose is writing, when H^2 - K < 0:

`k_max = std::exp(-0.5 * tmp * tmp / SIGMA * SIGMA) * h;`

Where SIGMA equals 0.001 for example.
This will allow a smooth transition from undefined to small.

1 Like

I like above approach.

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.

2 Likes

Here is another example which shows clear issues when kmin, kmax are identical or very very close. This is a sphere with radius 1.

Min curvature. Top VTK. Bottom numpy

Max curvature. Top VTK. Bottom numpy

Sphere was generated using this code snippet:

``````radius = 1
sph = vtk.vtkSphereSource()
sph.SetPhiResolution(360)
sph.SetThetaResolution(360)
sph.LatLongTessellationOff()

sph.Update()
poly = sph.GetOutput()
``````

Can you display the Gaussian curvature on your examples?

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

Top: Gaussian; Bottom: Mean

Same as above, but rescaled in range [0.999, 1.001] to show minor differences (numerical or discretization?) that appear.