Bounds errors in curvature calculation

[bragostin]

I am using the following code (CQ master + OCP 7.5.1beta) to integrate the curvature over the surface of a quarter of a cylindrical shell:

import cadquery as cq
import numpy as np
from scipy.integrate import dblquad

from OCP.GeomLProp import GeomLProp_SLProps, GeomLProp_SurfaceTool
from OCP.BRep import BRep_Tool

def curvature(cq_obj):

    def C(u, v, curvature):
        curvature.SetParameters(u, v)
        if curvature.IsCurvatureDefined():
            crv_max = curvature.MaxCurvature() # (1/mm)
        else:
            crv_max = 0
        return crv_max

    faces = cq_obj.val().Faces()
    h_srf = [BRep_Tool().Surface_s(f.wrapped) for f in faces]

    Suv = []
    res = []
    for h in h_srf:
        uv = GeomLProp_SurfaceTool().Bounds_s(h)
        umin, umax, vmin, vmax = uv[0], uv[2], uv[1], uv[3]
        uv_s = (umax - umin) * (vmax - vmin) 
        print("umin, umax, vmin, vmax = ", umin, umax, vmin, vmax)
        curvature = GeomLProp_SLProps(h, 2, 1e-6)
        res.append(dblquad(C, umin, umax, vmin, vmax, args=(curvature,), epsabs=1.49e-08, epsrel=1.49e-08))
        Suv.append(uv_s)

    print("res = ", res)
    print("Suv = ", Suv)

    crv = np.array([r[0]/S for r, S in zip(res, Suv)])

    return crv

r = 5
FIN01 = cq.Workplane('XZ').rect(10,0.01).revolve(angleDegrees=90, axisStart=(0, r, 0), axisEnd=(r, r, 0))
cq.exporters.export(FIN01, "curv.step")
cv = curvature(FIN01)
print("cv = ", cv)

which gives the output:

umin, umax, vmin, vmax =  0.0 6.283185307179586 -2e+100 2e+100
umin, umax, vmin, vmax =  -2e+100 2e+100 -2e+100 2e+100
umin, umax, vmin, vmax =  0.0 6.283185307179586 -2e+100 2e+100
umin, umax, vmin, vmax =  -2e+100 2e+100 -2e+100 2e+100
umin, umax, vmin, vmax =  -2e+100 2e+100 -2e+100 2e+100
umin, umax, vmin, vmax =  -2e+100 2e+100 -2e+100 2e+100
res =  [(0.0, 0), (0.0, 0), (5.031579825569239e+100, 5.586175772586398e+86), (0.0, 0), (0.0, 0), (0.0, 0)]
Suv =  [2.5132741228718345e+101, 1.6e+201, 2.5132741228718345e+101, 1.6e+201, 1.6e+201, 1.6e+201]
cv =  [0.        0.        0.2002002 0.        0.        0.       ]

The surface is indeed made of 6 faces and the final curvature value of 0.2 = 1/5 seems correct. There are however 2 issues in my opinion: 1 - the bounds do not seem correct : +-2e100 is way too large 2 - there is only 1 face with a curvature, while there are 2 cylindrical surfaces

So I am wondering if : 1 - this is normal 2 - it is a CQ issue 3 - it is an OCC issue

Any clue?

[adam-urbanczyk]

You are using wrong bounds - underlying geometry as opposed to the trimmed face. There is a CQ method for this: cq.Face._uvBounds(). Let me know if it solves the curvature issue too.

[bragostin]

Now I get:

umin, umax, vmin, vmax =  0.0 1.5707963267948966 0.0 10.0
umin, umax, vmin, vmax =  3.0585553808705147e-16 5.005 0.0 5.005
umin, umax, vmin, vmax =  0.0 1.5707963267948966 0.0 10.0
umin, umax, vmin, vmax =  -5.005 -3.0585553808705147e-16 0.0 5.005
umin, umax, vmin, vmax =  -0.005 0.005 -5.0 5.0
umin, umax, vmin, vmax =  -0.005 0.005 -5.0 5.0
res =  [(0.0, 0), (0.0, 0), (3.144737390980774, 3.491359857866499e-14), (0.0, 0), (0.0, 0), (0.0, 0)]
Suv =  [15.707963267948966, 25.050024999999998, 15.707963267948966, 25.050024999999998, 0.1, 0.1]
cv =  [0.        0.        0.2002002 0.        0.        0.       ]

So the 2e100 have indeed disappeared, however:

• still only one face yields a non zero curvature
• the 15.71 mm2 (1/4 circle area) area does not correspond to any existing face

That said in the end the curvature value is correct, but it is strange.

Updated code:

import cadquery as cq
import numpy as np
from scipy.integrate import dblquad

from OCP.GeomLProp import GeomLProp_SLProps, GeomLProp_SurfaceTool
from OCP.BRep import BRep_Tool

def curvature(cq_obj):

    def C(u, v, curvature):
        curvature.SetParameters(u, v)
        if curvature.IsCurvatureDefined():
            crv_max = curvature.MaxCurvature() # (1/mm)
        else:
            crv_max = 0
        return crv_max

    faces = cq_obj.val().Faces()
    h_srf = [BRep_Tool().Surface_s(f.wrapped) for f in faces]

    Suv = []
    res = []
    for h, f in zip(h_srf, faces):
        uv = cq.Face(f.wrapped)._uvBounds()
        umin, umax, vmin, vmax = uv[0], uv[1], uv[2], uv[3]       
        uv_s = (umax - umin) * (vmax - vmin) 
        print("umin, umax, vmin, vmax = ", umin, umax, vmin, vmax)
        curvature = GeomLProp_SLProps(h, 2, 1e-6)
        res.append(dblquad(C, umin, umax, vmin, vmax, args=(curvature,), epsabs=1.49e-08, epsrel=1.49e-08))
        Suv.append(uv_s)

    print("res = ", res)
    print("Suv = ", Suv)

    crv = np.array([r[0]/S for r, S in zip(res, Suv)])

    return crv

r = 5
FIN01 = cq.Workplane('XZ').rect(10,0.01).revolve(angleDegrees=90, axisStart=(0, r, 0), axisEnd=(r, r, 0))
cq.exporters.export(FIN01, "curv.step")
cv = curvature(FIN01)
print("cv = ", cv)

[bragostin]

I think I got a grip of what is going on:

• one face has 0.2 curvature and the other -0.2, so with MaxCurvature one yields 0.2 and the other one 0. Using both MinCurvature and MaxCurvature yields the expected results of 2 faces with curvatures.
• The 2 faces yielding curvatures have 0<u<pi and 0<v<10 as parameters, not cartesian coordinates, thus u*v does not yield a metric area, just the integration domain.

Updated code:

import cadquery as cq
import numpy as np
from scipy.integrate import dblquad

from OCP.GeomLProp import GeomLProp_SLProps
from OCP.BRep import BRep_Tool

def curvature(cq_obj):

    def C(u, v, curvature):
        curvature.SetParameters(u, v)
        if curvature.IsCurvatureDefined():
            crv_max = curvature.MaxCurvature() # (1/mm)
            crv_min = curvature.MinCurvature() # (1/mm)
        else:
            crv_max = 0
            crv_min = 0
        cv = [crv_min, crv_max]
        return sum(cv) 

    faces = cq_obj.val().Faces()

    Suv = []
    res = []
    curvature = GeomLProp_SLProps(N=2, Resolution=1e-6)
    for f in faces:
        h = BRep_Tool().Surface_s(f.wrapped)
        uv = cq.Face(f.wrapped)._uvBounds()
        umin, umax, vmin, vmax = uv[0], uv[1], uv[2], uv[3]       
        uv_s = (umax - umin) * (vmax - vmin) 
        print("umin, umax, vmin, vmax = ", umin, umax, vmin, vmax)
        curvature.SetSurface(h)
        res.append(dblquad(C, umin, umax, vmin, vmax, args=(curvature,), epsabs=1.49e-08, epsrel=1.49e-08))
        Suv.append(uv_s)

    print("res = ", res)
    print("Suv = ", Suv)

    crv = np.array([r[0]/S for r, S in zip(res, Suv)])

    return crv

r = 5
FIN01 = cq.Workplane('XZ').rect(10,0.01).revolve(angleDegrees=180, axisStart=(0, r, 0), axisEnd=(r, r, 0))
cq.exporters.export(FIN01, "curv.step")
cv = curvature(FIN01)
print("cv = ", cv)

yielding

umin, umax, vmin, vmax =  -5.005000000000001 5.005 0.0 5.005
umin, umax, vmin, vmax =  0.0 3.141592653589793 0.0 10.0
umin, umax, vmin, vmax =  -5.005 5.005 0.0 5.005
umin, umax, vmin, vmax =  -0.005 0.005 -5.0 5.0
umin, umax, vmin, vmax =  -0.005 0.005 -5.0 5.0
res =  [(-6.276908398780805, 6.968768227789474e-14), (0.0, 0), (6.289474781961548, 6.982719715732997e-14), (0.0, 0), (0.0, 0), (0.0, 0)]
Suv =  [31.41592653589793, 50.10005000000001, 31.41592653589793, 50.100049999999996, 0.1, 0.1]
cv =  [-0.1998002  0.         0.2002002  0.         0.         0.       ]
``