Singularity in fitting.interpolate_surface()

[sphh]

I try to fit a surface through a field of points, where one row of points collapse to a singularity, similar to this:
https://developer.rhino3d.com/guides/general/essential-mathematics/parametric-curves-surfaces/#singularity-in-nurbs-surfaces

At the moment I get a ZeroDivisionError:

...
  File "/usr/local/lib/python3.8/dist-packages/geomdl/fitting.py", line 38, in interpolate_curve
    uk = compute_params_curve(points, use_centripetal)
  File "/usr/local/lib/python3.8/dist-packages/geomdl/fitting.py", line 452, in compute_params_curve
    uk[i] = sum(cds[0:i + 1]) / d
ZeroDivisionError: float division by zero

Is there a workaround to get such a surface? Is it possible at all?

[orbingol]

compute_params_curve needs to compute the chord length and when there is a singularity, the chord length becomes zero. I think there are 2 things to do:

1. It might not work perfectly but there could be a chance to fix the issue by catching ZeroDivisionError and filling the uk list separately
2. Put a very small distance between the data points that are causing the singularity. Therefore, the division will be a very small number but not zero

[sphh]

Currently I use 2, but would like to have a proper singularity.

If 1. is possible, that would be great. What would uk look like if d→0? Is this limes defined?

[orbingol]

I'll need to test for Option 1. I am not sure if it will even work. The concept of chord length breaks the idea of singularity. However, it looks possible when we look at it from the programming perspective.

[portnov]

Usually in NURBS algorithms, if you're encountering 0/0, you just say it is 0, and algorithm works fine. Does this trick help here?

[orbingol]

Usually in NURBS algorithms, if you're encountering 0/0, you just say it is 0, and algorithm works fine. Does this trick help here?

I was thinking the same thing but I'd prefer testing before coming to the conclusion. Sometimes, it ends up being a completely different bug.

[sphh]

I replaced line 452

        uk[i] = sum(cds[0:i + 1]) / d

with

        try:
            uk[i] = sum(cds[0:i + 1]) / d
        except ZeroDivisionError:
            if sum(cds[0:i + 1]) == 0:
                uk[i] = 0
            else:
                raise

(maybe this is faster:

        if (s := sum(cds[0:i + 1])) == 0 and d == 0:
            uk[i] = 0
        else:
            uk[i] = s / d

is faster? But it needs Python ≥3.8)

Both my test and real cases are working perfectly with both changes.

[orbingol]

ZeroDivisionError occurs when d == 0. If sum(cds[0:i + 1]) == 0 and d != 0, it will evaluate to 0. I think it would be enough to update the line 452 like this:

try:
   uk[i] = sum(cds[0:i + 1]) / d
except ZeroDivisionError:
   uk[i] = 0

Note: I haven't tested it.

Elimination of if-else will automatically remove some of the overhead and it will become just a bit faster.

[sphh]

I have not looked into the code more closely, but I wanted to implement it strictly according to @portnov! If sum(cds[0:1 + 1]) is always 0 whenever d==0, my approach is certainly overkill.

I have different experiences when I timed try ... except ... and if ... else ... ages ago. I remembered, that try ... took more time. I just had another look into it again and noticed, that it depends, how often the exception is triggered (https://stackoverflow.com/questions/1835756/using-try-vs-if-in-python, https://www.geeksforgeeks.org/try-except-vs-if-in-python/). Based on these findings (and that singularities do not appear often), it looks like your approach is the faster one. Still, is it assured, that the sum() is always 0 when d is 0?

[sphh]

Oh, well I managed to construct a test case, which fails:

import numpy as np

from geomdl import fitting
from geomdl.visualization import VisMPL

def ellipse(x, a, b):
    """Elliptic function."""
    return b * np.sqrt(1 - (x/a)**2)

def profile(x, a, b, t_te):
    """Profile with open trailing edge."""
    return ellipse(x - (a-t_te), a+t_te, b)

X = np.linspace(0, 1, 5, endpoint=True)
Y = np.linspace(0, 1, 5, endpoint=True)

sections = []
for xi in X:
    chord = ellipse(xi, 1, 0.5)
    if chord == 0:
        chord = 0.000001
    thickness = ellipse(xi, 1, 0.1)
    # Upper face
    section = [
        [xi, (yi-0.5)*chord, profile(yi, 0.5, thickness, 0.01)]
        for yi in Y]
    # Lower face
    section.extend([
        [xi, (yi-0.5)*chord, -profile(yi, 0.5, thickness, 0.01)]
        for yi in reversed(Y[:-1])])
    sections.append(section)

surface = fitting.interpolate_surface(
    [list(point)
     for section in sections
     for point in section],
    len(sections), len(sections[0]),
    3, 3)

vis = VisMPL.VisSurface()
surface.vis = vis
surface.render()

With this code I get the following wing, which is how I expect it to look like:

If I comment out these lines (which removed the singularity a the tip) to have a singularity at the tip x=0:

#     if chord == 0:
#         chord = 0.000001

I get the following wing:

[portnov]

@sphh as far as I understood, you want your profile to be elliptical. For ellipses and other conic sections there are known algorithms that define control points / weights / knot vectors of NURBS curves directly, without interpolation. Then you can apply lofting algorithm to generate surface. Such surface can be much "better" in several aspects, comparing to one which you can generate by interpolation.
Just an example:

[sphh]

Well, it's true, that I used elliptical sections, section widths and thickness here. But this was only used to construct the simple test case and in real life different sections and distributions will be used. But thanks for pointing that out. It's highly appreciated. (@orbingol: Do you want to implement them in the geomdl.shapes package??)

BTW which nice program are you using here?

[portnov]

Blender (blender.org) + Sverchok (https://github.com/nortikin/sverchok) :)

[portnov]

as far as I see, geomdl.shapes.curve2d module already has methods for generating circles. To generate an ellipse, it is enough to generate a circle and then scale it's control points in one direction.

[sphh]

@portnov: As mentioned earlier, this is just a testcase, the actual case is much more complicated. So let us not get side-tracked by figuring out, how the test case can be created in a different/better/faster/easier way, ok?