search

orbingol / NURBS-Python

Object-oriented pure Python B-Spline and NURBS library
https://onurraufbingol.com/NURBS-Python/
MIT License
640 stars 156 forks source link

Help wanted: unknown ZeroDivisionError #168

Open chenboshuo opened 1 year ago

chenboshuo commented 1 year ago

this is my code

from geomdl import BSpline

# Create a BSpline surface instance (Bezier surface)
surf = BSpline.Surface()

surf.degree_u=2
surf.degree_v=2
surf.ctrlpts_size_u=5
surf.ctrlpts_size_v=9
surf.ctrlpts=[
                (52.5, 60, 0),
                (52.5, 60, 60),
                (52.5, 0, 60),
                (52.5, -59.999999999999986, 60.000000000000007),
                (52.5, -60, 0),
                (52.5, -60.000000000000007, -59.999999999999986),
                (52.5, 0, -60),
                (52.5, 60, -60.000000000000021),
                (52.5, 60, 0),
                (17.500000000000004, 59.999999999999986, 0),
                (17.5, 59.999999999999986, 59.999999999999986),
                (17.500000000000004, 0, 59.999999999999986),
                (17.5, -59.999999999999979, 59.999999999999993),
                (17.500000000000004, -59.999999999999986, 0),
                (17.5, -59.999999999999993, -59.999999999999979),
                (17.500000000000004, 0, -59.999999999999986),
                (17.5, 59.999999999999986, -60),
                (17.500000000000004, 59.999999999999986, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
]
surf.weights=[
1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 6.95192e-310, 52.5, 60, 0, 52.5, 60, 60, 52.5, 3.67e-15, 60, 52.5, -60, 60, 52.5, -60, 7.35e-15, 52.5, -60, ]
surf.knotvector_u=[0, 0, 0, 1, 1, 1, 1, 1, ]
surf.knotvector_v=[0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1, ]
# Evaluate surface points
surf.evaluate()

and I get below error

ZeroDivisionError                         Traceback (most recent call last)
Input In [83], in <cell line: 1>()
     60 surf.knotvector_v=[0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1, ]
     61 # Evaluate surface points
---> 62 surf.evaluate()
     64 # Plot the control points grid and the evaluated surface
     65 surf.vis = VisMPL.VisSurface(vis_config)

File ~/anaconda3/lib/python3.9/site-packages/geomdl/BSpline.py:616, in Surface.evaluate(self, **kwargs)
    613 self.reset(evalpts=True)
    615 # Evaluate and cache
--> 616 self._eval_points = self._evaluator.evaluate(self.data,
    617                                              start=(start_u, start_v),
    618                                              stop=(stop_u, stop_v))

File ~/anaconda3/lib/python3.9/site-packages/geomdl/evaluators.py:291, in SurfaceEvaluator.evaluate(self, datadict, **kwargs)
    289     knots = linalg.linspace(start[idx], stop[idx], sample_size[idx], decimals=precision)
    290     spans[idx] = helpers.find_spans(degree[idx], knotvector[idx], size[idx], knots, self._span_func)
--> 291     basis[idx] = helpers.basis_functions(degree[idx], knotvector[idx], spans[idx], knots)
    293 eval_points = []
    294 for i in range(len(spans[0])):

File ~/anaconda3/lib/python3.9/site-packages/geomdl/helpers.py:250, in basis_functions(degree, knot_vector, spans, knots)
    248 basis = []
    249 for span, knot in zip(spans, knots):
--> 250     basis.append(basis_function(degree, knot_vector, span, knot))
    251 return basis

File ~/anaconda3/lib/python3.9/site-packages/geomdl/helpers.py:167, in basis_function(degree, knot_vector, span, knot)
    165 saved = 0.0
    166 for r in range(0, j):
--> 167     temp = N[r] / (right[r + 1] + left[j - r])
    168     N[r] = saved + right[r + 1] * temp
    169     saved = left[j - r] * temp

ZeroDivisionError: float division by zero

The equivalent ACIS code

// spline
    int degree_of_u = 2;
    logical is_rational_in_u = TRUE;
    int OPEN = 0;
    int form_u = OPEN;
    int NO_SINGULARITY_AT_U_BOUNDARY = 0;
    int pole_u = NO_SINGULARITY_AT_U_BOUNDARY;
    int control_points_in_u = 5;
    int degree_v = 2;
    logical is_rational_in_v = TRUE;
    int form_v = 2;
    int pole_v = NO_SINGULARITY_AT_U_BOUNDARY;
    int control_points_in_v = 9;
    SPAposition control_points[] = {
      {52.50000000000000000, 60.00000000000000000,  0.00000000000000000  },
      {52.50000000000000000, 60.00000000000000000,  60.00000000000000000 },
      {52.50000000000000000, 0.00000000000000367,   60.00000000000000000 },
      {52.50000000000000000, -59.99999999999998579, 60.00000000000000711 },
      {52.50000000000000000, -60.00000000000000000, 0.00000000000000735  },
      {52.50000000000000000, -60.00000000000000711, -59.99999999999998579},
      {52.50000000000000000, -0.00000000000001102,  -60.00000000000000000},
      {52.50000000000000000, 60.00000000000000000,  -60.00000000000002132},
      {52.50000000000000000, 60.00000000000000000,  0.00000000000000000  },
      {17.50000000000000355, 59.99999999999998579,  0.00000000000000000  },
      {17.50000000000000000, 59.99999999999998579,  59.99999999999998579 },
      {17.50000000000000355, 0.00000000000000367,   59.99999999999998579 },
      {17.50000000000000000, -59.99999999999997868, 59.99999999999999289 },
      {17.50000000000000355, -59.99999999999998579, 0.00000000000000735  },
      {17.50000000000000000, -59.99999999999999289, -59.99999999999997868},
      {17.50000000000000355, -0.00000000000001102,  -59.99999999999998579},
      {17.50000000000000000, 59.99999999999998579,  -60.00000000000000000},
      {17.50000000000000355, 59.99999999999998579,  0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
      {17.50000000000000000, 0.00000000000000000,   0.00000000000000000  },
    };
    double weights[] = {
      1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000,
      0.70710678118654779, 0.50000000000000022, 0.70710678118654779, 0.50000000000000022, 0.70710678118654779, 0.50000000000000022, 0.70710678118654779, 0.50000000000000022, 0.70710678118654779,
      1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000, 0.70710678118654757, 1.00000000000000000,
    };
    double control_points_tolerance = SPAresabs;
    int num_knots_u = 6;
    double knots_u[] = {
      0.00000000000000000, 0.00000000000000000, 0.00000000000000000, 1.00000000000000000, 1.00000000000000000, 1.00000000000000000,
    };
    int num_knots_v = 12;
    double knots_v[] = {
      0.00000000000000000, 0.00000000000000000, 0.00000000000000000, 0.25000000000000000, 0.25000000000000000, 0.50000000000000000,
      0.50000000000000000, 0.75000000000000000, 0.75000000000000000, 1.00000000000000000, 1.00000000000000000, 1.00000000000000000,
    };
    double knot_tol = SPAresabs;
    //
    bs3_surface bezier_surface = bs3_surface_from_ctrlpts(degree_of_u, is_rational_in_u, form_u, pole_u, control_points_in_u, degree_v, is_rational_in_v, form_v, pole_v, control_points_in_v, control_points, weights, control_points_tolerance, num_knots_u,
                                                          knots_u, num_knots_v, knots_v, knot_tol);
    spline spl = spline(bezier_surface);

Could you please take a look at this bug ? I would greatly appreciate your assistance. Thank you in advance!

orbingol commented 1 year ago

I don't have a valid ACIS license anymore, so I cannot test it. ZeroDivisionError may happen, because the algorithms are using the user input directly. There is a frange implementation, but I doubt that's the problem. There were several other problems like that and converting values like -59.99999999999998579 to -60 worked fine at that point. You could try that.

I was planning to implement a more precise floating point control to eliminate such issues, but never had time to do it without adding any additional dependencies.

chenboshuo commented 1 year ago

I don't have a valid ACIS license anymore, so I cannot test it. ZeroDivisionError may happen, because the algorithms are using the user input directly. There is a frange implementation, but I doubt that's the problem. There were several other problems like that and converting values like -59.99999999999998579 to -60 worked fine at that point. You could try that.

I was planning to implement a more precise floating point control to eliminate such issues, but never had time to do it without adding any additional dependencies.

I changed this snippets to

from geomdl import BSpline

# Create a BSpline surface instance (Bezier surface)
surf = BSpline.Surface()

surf.degree_u=2
surf.degree_v=2
surf.ctrlpts_size_u=5
surf.ctrlpts_size_v=9
surf.ctrlpts=[
                (52.5, 60, 0),
                (52.5, 60, 60),
                (52.5, 0, 60),
                (52.5, -60, 60.0),
                (52.5, -60, 0),
                (52.5, -60.0, -60),
                (52.5, 0, -60),
                (52.5, 60, -60.0),
                (52.5, 60, 0),
                (17.5, 60, 0),
                (17.5, 60, 60),
                (17.5, 0, 60),
                (17.5, -60, 60),
                (17.5, -60, 0),
                (17.5, -60, -60),
                (17.5, 0, -60),
                (17.5, 60, -60),
                (17.5, 60, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (17.5, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
                (0, 0, 0),
]
surf.weights=[
1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 0.5, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0.707107, 1, 0, 52.5, 60, 0, 52.5, 60, 60, 52.5, 0, 60, 52.5, -60, 60, 52.5, -60, 0, 52.5, -60, ]
surf.knotvector_u=[0, 0, 0, 1, 1, 1, 1, 1, ]
surf.knotvector_v=[0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1, ]
# Evaluate surface points
surf.evaluate()

# Plot the control points grid and the evaluated surface
surf.vis = VisMPL.VisSurface(vis_config)

I am still encountering the same error.

orbingol commented 1 year ago

Could you also update the weights like you did for the control points? I'd try converting 0.707107 to 0.70711 and remove one significant if it doesn't work.