toNurbs() broken for circular arcs?

Need help, or want to share a macro? Post here!
Forum rules
Be nice to others! Respect the FreeCAD code of conduct!
edwilliams16
Veteran
Posts: 3107
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

toNurbs() broken for circular arcs?

Post by edwilliams16 »

Maybe I'm doing it wrong - but if I create a circular arc:

Code: Select all

from math import degrees
doc = App.ActiveDocument
arc = doc.addObject("Part::Circle","Arc")
arc.Radius=10.0000
arc.Angle1 =25.0000
arc.Angle2=55.0000
doc.recompute()
then convert it to nurbs (which should be exact for conics), the start point is correct but the end point wrong

Code: Select all

edge = arc.Shape.Edge1
curve = edge.Curve.trim(*arc.Shape.ParameterRange)
nurbs = curve.toNurbs() # or nurbs = curve.toNurbs(*arc.Shape.ParameterRange)
print(' Arc      Nurbs')
print(f'Start Points {curve.StartPoint}  {nurbs.StartPoint}')
print(f'End Points {curve.EndPoint}   {nurbs.EndPoint}')
arcangle = degrees(curve.StartPoint.getAngle(curve.EndPoint))
nurbsangle = degrees(nurbs.StartPoint.getAngle(nurbs.EndPoint))
print(f'Subtend {arcangle}  {nurbsangle}')
Output:

Code: Select all

>>> print(' Arc      Nurbs')
 Arc      Nurbs
>>> print(f'Start Points {curve.StartPoint}  {nurbs.StartPoint}')
Start Points Vector (9.063077870366499, 4.2261826174069945, 0.0)  Vector (9.063077870366499, 4.2261826174069945, 0.0)
>>> print(f'End Points {curve.EndPoint}   {nurbs.EndPoint}')
End Points Vector (5.735764363510461, 8.191520442889917, 0.0)   Vector (6.740463181335613, 7.386890827747417, 0.0)
>>> arcangle = degrees(curve.StartPoint.getAngle(curve.EndPoint))
>>> nurbsangle = degrees(nurbs.StartPoint.getAngle(nurbs.EndPoint))
>>> print(f'Subtend {arcangle}  {nurbsangle}')
Subtend 30.000000000000004  22.619864948040433
>>> 
Is this a bug?

EDIT:
FWIW I also tried this with a spiral arc - no problem.
It also works for ellipses, elliptic arcs and parabolic arcs. So far only circular arcs are a problem.
wmayer
Founder
Posts: 20243
Joined: Thu Feb 19, 2009 10:32 am
Contact:

Re: toNurbs() broken for circular arcs?

Post by wmayer »

Yes, something looks quite wrong. If you display the nurbs curve with

Code: Select all

Part.show(nurbs.toShape())
then the arc is too short.

If you create another arc with e.g. 26 deg as start and 337 deg as end angles then the converted nurbs has the right length but the start and end points are displaced.
User avatar
jnxd
Posts: 951
Joined: Mon Mar 30, 2015 2:30 pm
Contact:

Re: toNurbs() broken for circular arcs?

Post by jnxd »

@edwilliams16 could you share the version you tested with? Could you also test with 0.20/0.19 for good measure in case something changed?
My latest (or last) project: B-spline Construction Project.
wmayer
Founder
Posts: 20243
Joined: Thu Feb 19, 2009 10:32 am
Contact:

Re: toNurbs() broken for circular arcs?

Post by wmayer »

Could you also test with 0.20/0.19 for good measure in case something changed?
They all suffer from the same problem.

Code: Select all

OS: Ubuntu 22.04.1 LTS (XFCE/xubuntu)
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.19.24416 (Git)
Build type: Unknown
Branch: releases/FreeCAD-0-19
Hash: 0a49fd05a46447aba70b6f18b5bb8210175dd471
Python version: 3.10.6
Qt version: 5.15.3
Coin version: 4.0.0
OCC version: 7.5.1
Locale: German/Germany (de_DE)

Code: Select all

OS: Ubuntu 22.04.1 LTS (XFCE/xubuntu)
Word size of FreeCAD: 64-bit
Version: 0.20.2.29614 (Git)
Build type: Unknown
Branch: releases/FreeCAD-0-20
Hash: 74e5a7675929a5a058b36a05b0a0f01a19017bf4
Python 3.10.6, Qt 5.15.3, Coin 4.0.0, Vtk 7.1.1, OCC 7.5.1
Locale: German/Germany (de_DE)

Code: Select all

OS: Ubuntu 22.04.1 LTS (XFCE/xubuntu)
Word size of FreeCAD: 64-bit
Version: 0.21.0.31757 +6 (Git)
Build type: debug
Branch: refactoring_settings
Hash: 2e3feb156c4822bad2116c0009b675d813384ee4
Python 3.10.6, Qt 5.15.3, Coin 4.0.0, Vtk 7.1.1, OCC 7.5.1
Locale: German/Germany (de_DE)
edwilliams16
Veteran
Posts: 3107
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: toNurbs() broken for circular arcs?

Post by edwilliams16 »

Code: Select all

OS: macOS 10.16
Word size of FreeCAD: 64-bit
Version: 0.21.0.31709 (Git)
Build type: Release
Branch: master
Hash: e188802ca6997d2564e7570ab648462e6a059f87
Python 3.10.8, Qt 5.15.6, Coin 4.0.0, Vtk 9.1.0, OCC 7.6.3
Locale: C/Default (C)
Installed mods: 
  * MeshRemodel 1.8919.0
  * Trails 2022.1.0
  * offline-documentation.backup1670387333.286728 1.0.0-alpha (Disabled)
  * DynamicData 2.46.0
  * offline-documentation.backup1668887250.997013 1.0.0-alpha (Disabled)
  * Curves.backup1667610597.862062 0.5.12 (Disabled)
  * fcgear 1.0.0
  * workfeature
  * offline-documentation.backup1673656435.692312 1.0.0-alpha (Disabled)
  * Pyramids-and-Polyhedrons
  * QuickMeasure-main
  * GDML 1.8.0
  * offline-documentation
  * Manipulator 1.4.9
  * fasteners 0.4.53
  * lattice2 1.0.0
  * CurvedShapes 1.0.4
  * toSketch 1.0.1
  * Help 1.0.3
  * Curves 0.6.5
and

Code: Select all

OS: macOS 10.16
Word size of OS: 64-bit
Word size of FreeCAD: 64-bit
Version: 0.19.24276 (Git)
Build type: Release
Branch: (HEAD detached at 0.19.1)
Hash: a88db11e0a908f6e38f92bfc5187b13ebe470438
Python version: 3.8.8
Qt version: 5.12.9
Coin version: 4.0.0
OCC version: 7.4.0
Locale: C/Default (C)


and

Code: Select all

OS: macOS 10.16
Word size of FreeCAD: 64-bit
Version: 0.20.29177 (Git)
Build type: Release
Branch: (HEAD detached at 0.20)
Hash: 68e337670e227889217652ddac593c93b5e8dc94
Python 3.9.13, Qt 5.12.9, Coin 4.0.0, Vtk 9.1.0, OCC 7.5.3
Locale: C/Default (C)
Installed mods: 
  * MeshRemodel 1.8919.0
  * Trails 2022.1.0
  * offline-documentation.backup1670387333.286728 1.0.0-alpha (Disabled)
  * DynamicData 2.46.0
  * offline-documentation.backup1668887250.997013 1.0.0-alpha (Disabled)
  * Curves.backup1667610597.862062 0.5.12 (Disabled)
  * fcgear 1.0.0
  * workfeature
  * offline-documentation.backup1673656435.692312 1.0.0-alpha (Disabled)
  * Pyramids-and-Polyhedrons
  * QuickMeasure-main
  * GDML 1.8.0
  * offline-documentation
  * Manipulator 1.4.9
  * fasteners 0.4.53
  * lattice2 1.0.0
  * CurvedShapes 1.0.4
  * toSketch 1.0.1
  * Help 1.0.3
  * Curves 0.6.5
all fail. However the (wrong) endpoint of the arcs differ with version. This suggests a possibility that the circle is being computed as some limiting case of a conic, where one ends up using near-zeros or -infinities with drastic loss of precision.
wmayer
Founder
Posts: 20243
Joined: Thu Feb 19, 2009 10:32 am
Contact:

Re: toNurbs() broken for circular arcs?

Post by wmayer »

What works is to create the nurbs shape directly with

Code: Select all

nurbs = edge.toNurbs()
With

Code: Select all

nurbs.Edge1.Curve.getWeights()
you can check that it's a real nurbs, not just a B-spline.

The function to blame is https://github.com/FreeCAD/FreeCAD/blob ... .cpp#L2411
edwilliams16
Veteran
Posts: 3107
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: toNurbs() broken for circular arcs?

Post by edwilliams16 »

wmayer wrote: Mon Feb 06, 2023 1:30 pm
The function to blame is https://github.com/FreeCAD/FreeCAD/blob ... .cpp#L2411
I believe this function leads one to https://github.com/FreeCAD/FreeCAD/blob ... .cpp#L2262
which it seems creates the segment of the curve of parameter length in u of last-first.
But the u-parameter is not angle (theta). If I plot uniformly spaced u's I get
Screenshot 2023-02-06 at 7.35.19 AM.png
Screenshot 2023-02-06 at 7.35.19 AM.png (16.67 KiB) Viewed 693 times
The exact representation of a circle as a Nurbs parameterizes the circle with the parameter t = tan(theta/2). Then sin(theta) = 2t/(1+t^2) and cos(theta) = (1 - t^2)/(1+ t^2), which are rational functions, representable by NURBS. But it isn't possible to exactly represent a circle as a NURBS directly parameterized by theta as that would imply that sin(theta) and cos(theta) are equal to some ratio of polynomials - which they aren't.

I'll try and figure out what edge.toNurbs() is doing.
edwilliams16
Veteran
Posts: 3107
Joined: Thu Sep 24, 2020 10:31 pm
Location: Hawaii
Contact:

Re: toNurbs() broken for circular arcs?

Post by edwilliams16 »

@wmayer

There's a fix to the inaccurate trimming - solve for the u-parameter as a function of angle. This is a straightforward Newton solve. My C++ is none existent, but here's the fix in Python

Code: Select all

import math
def findUCircle(circ, angle,  maxiter = 100):
    ''' NURBS fit to a circle: circ
    with parameter range (0, 2*pi)
    angle CCW from x-axis (radians): angle
    return the corresponding u-parameter and success/failure
    '''
    if abs(angle - 2*math.pi) < 1e-15: #problem when angle = 2*pi
        angle = 0
    u = angle
    for i in range(maxiter):
        v = circ.value(u)
        anglenew = math.atan2(v.y, v.x)
        if anglenew < 0: # need  0<= anglenew <= 2*pi
            anglenew += 2 * math.pi
        unew = u + (angle - anglenew) #approx Newton solve - unit derivative
        #print(f'{angle}  {anglenew}  {u} {unew}')
        if abs(unew - u) < 1e-15:
            solve = True
            break
        u = unew

    return (unew, solve)

def trimNurbsCirc(circ, ang1, ang2):
    u1, solve1 = findUCircle(circ, ang1)
    u2, solve2 = findUCircle(circ, ang2)
    if solve1 and solve2:
        return circ.trim(u1,u2)
    else:
        print('findUCircle failed to converge') #have not been able to make it fail - but...
        return None


doc = App.ActiveDocument
arc = doc.addObject("Part::Circle","Arc")
arc.Radius=10.0000
arc.Angle1 =0.0000
arc.Angle2=60.0000
doc.recompute()
circle = doc.addObject("Part::Circle","Circle")
circle.Radius = 10
edge = arc.Shape.Edge1
curve = edge.Curve
nurbs = trimNurbsCirc(curve, *arc.Shape.ParameterRange)
Part.show(nurbs.toShape(), "Nurbs")
wmayer
Founder
Posts: 20243
Joined: Thu Feb 19, 2009 10:32 am
Contact:

Re: toNurbs() broken for circular arcs?

Post by wmayer »

It also works for ellipses, elliptic arcs and parabolic arcs. So far only circular arcs are a problem.
There it works because for them GeomCurve::toNurbs() is used that uses OCC's ShapeConstruct_Curve::ConvertToBSpline(). This function takes care of the correct parameter range but only creates non-rational splines.
What works is to create the nurbs shape directly with
Looking at the OCC code I found the relevant function to create real and trimmed NURBS. It's GeomConvert::CurveToBSplineCurve
Here you see the way how it's used: https://git.dev.opencascade.org/gitweb/ ... =HEAD#l473
Post Reply