Interpolating a hemispherical shape leads to ZeroDivisionError

[chvillanuevap]

Describe the bug
I'm trying to interpolate a hemispherical point cloud, but I get a division by zero error. I'm not sure what I'm doing wrong. This is the point cloud:

To Reproduce
Steps to reproduce the behavior:

from geomdl import fitting
from geomdl.visualization import VisMPL as vis

import numpy as np

theta = np.linspace(0, np.pi/2, num=50)
phi = np.linspace(0, 2*np.pi, num=100)
theta, phi = np.meshgrid(theta, phi)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
points = np.column_stack((x.ravel(), y.ravel(), z.ravel()))

# Set interpolation parameters
size_u = 50
size_v = 100
degree_u = 3
degree_v = 3

# Perform surface interpolation
surf = fitting.interpolate_surface(points, size_u, size_v, degree_u, degree_v,centripetal=True)

ZeroDivisionError                         Traceback (most recent call last)
Input In [32], in 
----> 1 surf = fitting.interpolate_surface(points, size_u, size_v, degree_u, degree_v,centripetal=True)

File [~/.local/lib/python3.9/site-packages/geomdl/fitting.py:82] in interpolate_surface(points, size_u, size_v, degree_u, degree_v, **kwargs)
     [79](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=78) use_centripetal = kwargs.get('centripetal', False)
     [81](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=80) # Get uk and vl
---> [82](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=81) uk, vl = compute_params_surface(points, size_u, size_v, use_centripetal)
     [84](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=83) # Compute knot vectors
     [85](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=84) kv_u = compute_knot_vector(degree_u, size_u, uk)

File [~/.local/lib/python3.9/site-packages/geomdl/fitting.py:485] in compute_params_surface(points, size_u, size_v, centripetal)
    [483](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=482) for v in range(size_v):
    [484](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=483)     pts_u = [points[v + (size_v * u)] for u in range(size_u)]
--> [485](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=484)     uk_temp += compute_params_curve(pts_u, centripetal)
    [487](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=486) # Do averaging on the u-direction
    [488](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=487) for u in range(size_u):

File [~/.local/lib/python3.9/site-packages/geomdl/fitting.py:452] in compute_params_curve(points, centripetal)
    [450](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=449) uk = [0.0 for _ in range(num_points)]
    [451](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=450) for i in range(num_points):
--> [452](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=451)     uk[i] = sum(cds[0:i + 1]) [/] d
    [454](file://~/.local/lib/python3.9/site-packages/geomdl/fitting.py?line=453) return uk

ZeroDivisionError: float division by zero

Expected Behavior
I get the control points. Also, what does centripetal do exactly? I have tried with both True and False.

Configuration:

• OS: Linux
• Python distribution: OS default
• Python version: 3.9.12
• geomdl install source: PyPI
• geomdl version/branch: 5.3.1

[chvillanuevap]

I tried the solution in #132. I no longer get the ZeroDivisionError, but I now get a closed surface as the interpolation.

(Red is data points, blue is evaluation points). Here is the code I use:

from geomdl import fitting, linalg
from geomdl.visualization import VisMPL as vis

import numpy as np
import matplotlib.pyplot as plt
import math
import h5py

np.set_printoptions(threshold=np.inf, linewidth=np.inf)

# From Github #132 
def my_compute_params(points, centripetal=False):
    if not isinstance(points, (list, tuple)):
        raise TypeError("Data points must be a list or a tuple")

    # Length of the points array
    num_points = len(points)

    # Calculate chord lengths
    cds = [0.0 for _ in range(num_points + 1)]
    cds[-1] = 1.0
    for i in range(1, num_points):
        distance = linalg.point_distance(points[i], points[i - 1])
        cds[i] = math.sqrt(distance) if centripetal else distance

    # Find the total chord length
    d = sum(cds[1:-1])

    # Divide individual chord lengths by the total chord length
    uk = [0.0 for _ in range(num_points)]
    for i in range(num_points):
        s = sum(cds[0:i + 1])
        if s == 0:
            uk[i] = i / (num_points - 1)
        else:
            try:
                uk[i] = s / d
            except ZeroDivisionError:
                uk[i] = 0

    return uk

theta = np.linspace(0, np.pi/2, num=50)
phi = np.linspace(0, 2*np.pi, num=100)
theta, phi = np.meshgrid(theta, phi)
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)
points = np.column_stack((x.ravel(), y.ravel(), z.ravel()))

# Set interpolation parameters
size_u = 50
size_v = 100
degree_u = 3
degree_v = 3

fitting.compute_params_curve = my_compute_params

# Perform surface interpolation
surf = fitting.interpolate_surface(points, size_u, size_v, degree_u, degree_v,centripetal=False)

surf.sample_size_u = 50
surf.sample_size_v = 100

# Evaluate the surface at the new resolution
surf.evaluate()

[orbingol]

The implemented algorithms, as they are, may not be able to handle surfaces like the hemisphere. Depending on how you generate the surface, you could play with the knot vectors or modify the fitting algorithm to use custom knot vectors. I am not sure custom knot vectors would be successful, but it might be worth a try.

[orbingol]

Just as a future reference to ZeroDivisionError: #168 (comment)