# Inaccurate principal curvature near peaks
**URL:** https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603
**Category:** Support
**Created:** [April 21, 2021, 4:30pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603 "2021-04-21T16:30:37Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 4:30pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/1 "2021-04-21T16:30:37Z")
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](https://vtk.org/doc/nightly/html/classvtkCurvatures.html#details). The numpy ones are saved as point-data and looks more reasonable and continuous.
Any suggestions? Explanations?
**Min Curvature**
**Max Curvature**
Code snippet to compute curvature:
```auto
>>> vtk.VTK_VERSION
'9.0.1'
```
```auto
def _vtk_curvature(pd, curvature_type_vtk, curvature_name):
c1 = vtk.vtkCurvatures()
c1.SetInputData(pd)
c1.SetCurvatureType(curvature_type_vtk)
c1.Update()
return c1.GetOutput(), numpy_support.vtk_to_numpy(c1.GetOutput().GetPointData().GetArray(curvature_name))
cv2vtk = {
'Maximum_Curvature': vtk.VTK_CURVATURE_MAXIMUM,
'Minimum_Curvature': vtk.VTK_CURVATURE_MINIMUM,
'Gauss_Curvature': vtk.VTK_CURVATURE_GAUSS,
'Mean_Curvature': vtk.VTK_CURVATURE_MEAN,
}
curvature_np = {}
for cv_name, cv_vtk in cv2vtk.items():
polydata, curvature_np[cv_name] = _vtk_curvature(polydata, cv_vtk, cv_name)
```
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 4:58pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/2 "2021-04-21T16:58:46Z")
Could you share your data set?
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 5:03pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/3 "2021-04-21T17:03:44Z")
I tried and got “Sorry, new users can not upload attachments.”
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 5:09pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/4 "2021-04-21T17:09:09Z")
> [@Yohann\_Bearzi](#):
>
> Could you share you data set?
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:
```auto
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
```
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 5:26pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/5 "2021-04-21T17:26:00Z")
FYI @cbhushan you should be now able to upload files. You profile just got upgraded. 🙂
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 5:31pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/6 "2021-04-21T17:31:18Z")
Thanks. The original file is attached.
[random\_hill1.vtk](https://discourse.vtk.org/uploads/short-url/kvHMp9zUBlS8I4A4uxrx2xRWANL.vtk) (365.4 KB)
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 6:13pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/7 "2021-04-21T18:13:58Z")
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.
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 6:16pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/8 "2021-04-21T18:16:55Z")
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.
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 6:40pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/9 "2021-04-21T18:40:35Z")
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.
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 6:42pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/10 "2021-04-21T18:42:09Z")
> [@Yohann\_Bearzi](#):
>
> 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.
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
```auto
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
```
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 6:43pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/11 "2021-04-21T18:43:34Z")
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.
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 6:53pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/12 "2021-04-21T18:53:47Z")
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?
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 7:00pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/13 "2021-04-21T19:00:02Z")
It is in `Filters/General/vtkCurvatures.cxx`, in `GetMinimumCurvature` and `GetMaximumCurvature`. It works for an error of 0.1, which is kind of large…
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 7:06pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/14 "2021-04-21T19:06:59Z")
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.
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 7:27pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/15 "2021-04-21T19:27:11Z")
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 🙂 ):
```auto
if std::abs(tmp) < this->errorWindow
{
tmp = 0.0;
}
```
---
### Author:  [@dgobbi](https://discourse.vtk.org/u/dgobbi)
#### Post date: [April 21, 2021, 7:27pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/16 "2021-04-21T19:27:38Z")
This topic attracted my attention, because I’m all in for fixing numerical instabilities in VTK.
I noticed this code in `ComputeGaussCurvature()`:
```auto
// 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.
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 7:45pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/17 "2021-04-21T19:45:02Z")
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:
```auto
radius = 1
sph = vtk.vtkSphereSource()
sph.SetRadius(radius)
sph.SetPhiResolution(360)
sph.SetThetaResolution(360)
sph.LatLongTessellationOff()
sph.Update()
poly = sph.GetOutput()
```
---
### Author:  [@Yohann\_Bearzi](https://discourse.vtk.org/u/Yohann_Bearzi)
#### Post date: [April 21, 2021, 7:52pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/18 "2021-04-21T19:52:34Z")
Can you display the Gaussian curvature on your examples?
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 7:57pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/19 "2021-04-21T19:57:24Z")
Just to be clear: numpy is not doing anything magical. Numpy calculations in [my script](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/10) 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).
---
### Author:  [@cbhushan](https://discourse.vtk.org/u/cbhushan)
#### Post date: [April 21, 2021, 7:58pm UTC](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603/20 "2021-04-21T19:58:20Z")
> [@Yohann\_Bearzi](#):
>
> Can you display the Gaussian curvature on your examples?
**Top: Gaussian; Bottom: Mean**
[Next page](https://discourse.vtk.org/t/inaccurate-principal-curvature-near-peaks/5603.md?page=2)