Division by zero in backward_substitution while fitting spline

[GarthSnyder]

Reproducible crash via division by zero. Should be reproducible from the code below. I took a look at the library code, but the exact issue isn't clear to me. Most likely it's something peculiar about the input, but I'm not sure exactly what.

Looks like matrix_u[q - 1][q - 1] is probably zero on line 616 of linalg.py.

There are 88 points, and accepting the default number of control points from fit.approximate_curve (which I assume is 87) causes the crash, as does any number of control points above this. The points fit just fine if fewer control points are requested.

I also notice that rounding all the floating point values avoids the problem, even if you round them to 15 digits. Very strange.

Environment: NURBS-Python 5.3.1, macOS Big Sur. Python 3.7 or 3.8.3. Regular Python is installed via Anaconda, but I'm working on a plugin for Fusion 360 which has its own Python environment; either reproduces the issue.

import geomdl
import geomdl.fitting as fit

points = [
    [2.844226598739624, -11.780293452192772, 0.09999999336045856],
    [2.7841949334192995, -11.766323767032048, 0.09999967816303046],
    [2.726926206709484, -11.752261321710279, 0.09999834367677496],
    [2.672275643786814, -11.73818458975435, 0.09999921155914912],
    [2.620158046052624, -11.72397519981413, 0.09999996191039656],
    [2.570553374917445, -11.709268999929673, 0.09999962409195731],
    [2.5232627391815186, -11.694355010986328, 0.10000000149011612],
    [2.478258938913115, -11.678929052642333, 0.09999946225906724],
    [2.435358747058011, -11.66315338641665, 0.09999991570513912],
    [2.394573211669922, -11.646668434143066, 0.10000000149011612],
    [2.355729818344116, -11.629592895507812, 0.10000000149011612],
    [2.3187353452599866, -11.611833538564408, 0.09999994560062432],
    [2.283545474902697, -11.59321500316845, 0.0999996723968506],
    [2.2499759197235107, -11.573817253112793, 0.10000000149011612],
    [2.217386549203785, -11.552999730025487, 0.09999874565806666],
    [2.1852133705257297, -11.530521281819617, 0.09999987781828826],
    [2.1537365913391113, -11.505938529968262, 0.10000000149011612],
    [2.122917890548706, -11.479263305664062, 0.10000000149011612],
    [2.092893351165359, -11.450263706561202, 0.09999973866089655],
    [2.0636894982676828, -11.418829252146876, 0.09999995528148392],
    [2.0354888439178467, -11.384699821472168, 0.10000000149011612],
    [2.0084729194641113, -11.34798812866211, 0.10000000149011612],
    [1.9837802648544312, -11.310267448425293, 0.10000000149011612],
    [1.9615306854248047, -11.272104263305664, 0.10000000149011612],
    [1.9416178464889526, -11.233770370483398, 0.10000000149011612],
    [1.9239351749420166, -11.195538520812988, 0.10000000149011612],
    [1.9083763360977173, -11.157683372497559, 0.10000000149011612],
    [1.894835114479065, -11.120477676391602, 0.10000000149011612],
    [1.8829901218414307, -11.083454132080078, 0.10000000149011612],
    [1.8725316644939427, -11.045799074770361, 0.09999994854771904],
    [1.8636144155532972, -11.007607214944535, 0.09999973105743598],
    [1.8561846017837524, -10.969038009643555, 0.10000000149011612],
    [1.8503996133333516, -10.930196091253926, 0.10000389344817495],
    [1.8471901202260397, -10.90117315661032, 0.1021130635318992],
    [1.8449852466583252, -10.872684478759766, 0.10844724625349045],
    [1.8437870740890503, -10.845902442932129, 0.11862903833389282],
    [1.843483992081289, -10.822776128482914, 0.13129674342530556],
    [1.8497409859620302, -10.79790687621756, 0.15000001355427042],
    [1.855139136314392, -10.781584739685059, 0.16294921934604645],
    [1.8570542447382592, -10.768168074449688, 0.17169737639867397],
    [1.8596694886032168, -10.754084788851447, 0.17936588915690657],
    [1.8629006147384644, -10.739490509033203, 0.18585821986198425],
    [1.866892695426941, -10.724578857421875, 0.19108659029006958],
    [1.8769696950912476, -10.694561004638672, 0.19784031808376312],
    [1.8894060850143433, -10.665621757507324, 0.20000000298023224],
    [1.9103977680206299, -10.627848625183105, 0.20000000298023224],
    [1.9356015920639038, -10.592745780944824, 0.20000000298023224],
    [1.9646819829940796, -10.56078052520752, 0.20000000298023224],
    [1.9972517490386963, -10.532379150390625, 0.20000000298023224],
    [2.032876968383789, -10.507918357849121, 0.20000000298023224],
    [2.071082830429077, -10.487727165222168, 0.20000000298023224],
    [2.111361265182495, -10.472070693969727, 0.20000000298023224],
    [2.153174877166748, -10.461159706115723, 0.20000000298023224],
    [2.195967197418213, -10.45513916015625, 0.20000000298023224],
    [2.239168167114258, -10.454089164733887, 0.20000000298023224],
    [2.2822024822235107, -10.458023071289062, 0.20000000298023224],
    [2.3244969844818115, -10.466890335083008, 0.20000000298023224],
    [2.365474916391236, -10.4806007670158, 0.1999998747902143],
    [2.4046443724843556, -10.498852355657805, 0.19999999233961327],
    [2.461829900741577, -10.526902198791504, 0.20000000298023224],
    [2.5210773944854736, -10.550283432006836, 0.20000000298023224],
    [2.5819969177246094, -10.568877220153809, 0.20000000298023224],
    [2.644202709197998, -10.582564353942871, 0.20000000298023224],
    [2.7073001861572266, -10.59126091003418, 0.20000000298023224],
    [2.7708895206451416, -10.594908714294434, 0.20000000298023224],
    [2.836341619491577, -10.596054077148438, 0.20000000298023224],
    [2.87705397605896, -10.59399700164795, 0.20000000298023224],
    [2.9171111583709717, -10.58643913269043, 0.20000000298023224],
    [2.9557740688323975, -10.573519706726074, 0.20000000298023224],
    [2.992328643798828, -10.555480003356934, 0.20000000298023224],
    [3.0261001586914062, -10.532649993896484, 0.20000000298023224],
    [3.0564651489257812, -10.505453109741211, 0.20000000298023224],
    [3.082862615585327, -10.474390983581543, 0.20000000298023224],
    [3.1048054695129395, -10.440035820007324, 0.20000000298023224],
    [3.1218883991241455, -10.403024673461914, 0.20000000298023224],
    [3.133796215057373, -10.364038467407227, 0.20000000298023224],
    [3.1403088569641113, -10.323798179626465, 0.20000000298023224],
    [3.141305923461914, -10.283045768737793, 0.20000000298023224],
    [3.1367692947387695, -10.242535591125488, 0.20000000298023224],
    [3.1267828941345215, -10.20301342010498, 0.20000000298023224],
    [3.1115307807922363, -10.165209770202637, 0.20000000298023224],
    [3.091294527053833, -10.1298246383667, 0.20000000298023224],
    [3.0664479732513428, -10.09750747680664, 0.20000000298023224],
    [3.037449598312378, -10.06885814666748, 0.20000000298023224],
    [3.0048351287841797, -10.044404029846191, 0.20000000298023224],
    [2.9692065715789795, -10.024596214294434, 0.20000000298023224],
    [2.9312217235565186, -10.00980281829834, 0.20000000298023224],
    [2.8915822505950928, -10.00029468536377, 0.20000000298023224]
]

curve = fit.approximate_curve(points, 5)

The error output is:

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-28-132e6a53193c> in <module>
----> 1 fit.approximate_curve(points, 5)

~/bin/anaconda3/lib/python3.8/site-packages/geomdl/fitting.py in approximate_curve(points, degree, **kwargs)
    196         b = [pt[i] for pt in vector_r]
    197         y = linalg.forward_substitution(matrix_l, b)
--> 198         x = linalg.backward_substitution(matrix_u, y)
    199         for j in range(1, num_cpts - 1):
    200             ctrlpts[j][i] = x[j - 1]

~/bin/anaconda3/lib/python3.8/site-packages/geomdl/linalg.py in backward_substitution(matrix_u, matrix_y)
    614     q = len(matrix_y)
    615     matrix_x = [0.0 for _ in range(q)]
--> 616     matrix_x[q - 1] = float(matrix_y[q - 1]) / float(matrix_u[q - 1][q - 1])
    617     for i in range(q - 2, -1, -1):
    618         matrix_x[i] = float(matrix_y[i]) - sum([matrix_u[i][j] * matrix_x[j] for j in range(i, q)])

ZeroDivisionError: float division by zero

[orbingol]

@GarthSnyder thanks for the bug report. I am guessing it would be a good idea to add a try-except block for ZeroDivisionError on line 616. Line 615 initializes the array to 0.0, so something like this should be good enough:

try:
  matrix_x[q - 1] = float(matrix_y[q - 1]) / float(matrix_u[q - 1][q - 1])
except ZeroDivisionError:
  pass

[orbingol]

It looks like changing the input degree also fixes the issue.

curve = fit.approximate_curve(points, 4)  # tested with 2, 3, 4, 6

Is there a specific reason for setting degree to 5 @GarthSnyder? Usually quadratic or cubic should be enough.

[GarthSnyder]

Fixed by #121

[neevparikh]

Hi I seem to have this issue when fitting a surface. Here is my reproducible snippet.

from geomdl.fitting import approximate_surface

pts = [
    (5.29724, 2.39077, -4.18169),
    (5.29212, 2.31915, -4.09873),
    (5.27242, 2.33416, -4.07725),
    (5.28571, 2.2567, -4.03315),
    (5.2707, 2.27565, -4.01087),
    (5.26706, 2.29867, -3.9834),
    (5.2648, 2.25254, -3.93714),
    (5.26704, 2.22952, -3.9618),
    (5.27635, 2.2083, -3.97793),
    (5.32034, 2.61671, -4.44898),
    (5.31134, 2.54153, -4.35878),
    (5.28584, 2.55673, -4.34687),
    (5.3047, 2.46589, -4.27014),
    (5.28338, 2.48316, -4.25543),
    (5.27473, 2.50775, -4.23376),
    (5.27542, 2.4077, -4.16565),
    (5.26814, 2.43144, -4.14071),
    (5.33211, 2.86337, -4.70251),
    (5.26721, 2.35746, -4.0507),
    (5.34108, 2.93704, -4.792),
    (5.30763, 2.8091, -4.59725),
    (5.36651, 2.92402, -4.80067),
    (5.32066, 2.88619, -4.68443),
    (5.33322, 2.9603, -4.77297),
    (5.27569, 2.65359, -4.42211),
    (5.2725, 2.57996, -4.32821),
    (5.29434, 2.73155, -4.51016),
    (5.30647, 2.70779, -4.52962),
    (5.29296, 2.63058, -4.43997),
    (5.32025, 2.78582, -4.61585),
    (5.34551, 2.77244, -4.62424),
    (5.33331, 2.69442, -4.53731),
    (5.35631, 2.85074, -4.71051),
]
surface_1 = approximate_surface(pts, 3, 11, 3, 3)
surface_2 = approximate_surface(pts, 3, 11, 3, 3, ctrlpts_size_u=4, ctrlpts_size_v=4)

Surface_1 results in a index error.

Stack trace:

Traceback (most recent call last):
  File "datagen/seam.py", line 123, in <module>
    surface_1 = approximate_surface(pts, 3, 11, 3, 3)
  File "/home/neev/.miniconda/envs/csm3.7/lib/python3.7/site-packages/geomdl/fitting.py", line 261, in approximate_surface
    matrix_ntnu = linalg.matrix_multiply(matrix_ntu, matrix_nu)
  File "/home/neev/.miniconda/envs/csm3.7/lib/python3.7/site-packages/geomdl/linalg.py", line 457, in matrix_multiply
    p1 = len(mat1[0])
IndexError: list index out of range

Trying to set the num_ctrl_pts to be 3 seems to have an issue. Surface_2 has this error:

Traceback (most recent call last):
  File "datagen/seam.py", line 124, in <module>
    surface_2 = approximate_surface(pts, 3, 11, 3, 3, ctrlpts_size_u=4, ctrlpts_size_v=4)
  File "/home/neev/.miniconda/envs/csm3.7/lib/python3.7/site-packages/geomdl/fitting.py", line 294, in approximate_surface
    x = linalg.backward_substitution(matrix_ntnuu, y)
  File "/home/neev/.miniconda/envs/csm3.7/lib/python3.7/site-packages/geomdl/linalg.py", line 616, in backward_substitution
    matrix_x[q - 1] = float(matrix_y[q - 1]) / float(matrix_u[q - 1][q - 1])
ZeroDivisionError: float division by zero

According to the code comments, the size of u or v needs to be such that (size_u + 1) > num_ctrl_pts, so in theory, for a 3x11 patch, the 3 pt direction should still work with 3 ctrlpts right?

Environment is: Python 3.7, geomdl 5.3.1, Linux.

[GarthSnyder]

@neevparikh, did you try this on the #121 fork? It would be interesting to know whether that fix addresses the 2D case. The fix is in a utility routine, so it might or might not be in the path for surfaces.

[neevparikh]

Oh okay I thought the #121 was merged in. I'll test and report back

[neevparikh]

No fix from that sadly. The 2nd surface error changes to the same error as the first one now.

[GarthSnyder]

It's been a while since I looked at NURBS math in detail, but I think you may be running afoul of some intrinsic requirements.

I don't think the standard approximation method implemented by NURBS-Python (from The NURBS Book) will extrapolate additional points in a dimension. If you have an input dimension of N, you can't ask for N + M control points as an output in that dimension; it has to be at most N and is N - 1 by default.

Similarly, you need at least degree + 1 control points in a dimension. This is a basic requirement for NURBS curves/surfaces. Otherwise, there's nothing to interpolate - all your control points are active at both the start and end of the interval.

Your input data is 3 x 11, so the highest number of control points you can obtain in the first dimension is 3. And that limits your output degree for that dimension to 2.

Rewriting your example in light of these limits and also using the #121 patch works fine. I'm not sure if it works without the #121 patch. (More specifically, it may run without errors without #121, although I believe the approximation is still potentially biased. You may see ringing in the output surface, but then again you might not. It depends on the input data.)

surface_1 = fit.approximate_surface(pts, 3, 11, 2, 3, ctrlpts_size_u=3)
surface_2 = fit.approximate_surface(pts, 3, 11, 2, 3, ctrlpts_size_u=3, ctrlpts_size_v=4)

Note that you have to specify the number of control points in U explicitly for both cases. Otherwise it defaults to 3 - 1 = 2, which is not sufficient to support a degree-2 NURBS.

[neevparikh]

Oh I see! I think I'd assumed 3 meant degree 3, but in hindsight that's obviously under determined, because you need n+1 pts to fit an n-dim polynomial. Thanks!

[neevparikh]

I wonder if the default behavior should warn the user of this edge case? Or maybe the default number of control pts should be max(size_u - 1, degree_u + 1)?

[GarthSnyder]

I wonder if the default behavior should warn the user of this edge case? Or maybe the default number of control pts should be max(size_u - 1, degree_u + 1)?

That seems like a good idea, but you should probably file it as a separate issue.