Computes the control points of the rational/non-rational spline after knot removal.

    Implementation based on Algorithm A5.8 and Equation 5.28 of The NURBS Book by Piegl & Tiller

    Keyword Arguments:
        * ``num``: number of knot removals

    :param degree: degree
    :type degree: int
    :param knotvector: knot vector
    :type knotvector: list, tuple
    :param ctrlpts: control points
    :type ctrlpts: list
    :param u: knot to be removed
    :type u: float
    :return: updated control points
    :rtype: list

def knot_removal(degree, knotvector, ctrlpts, u, **kwargs):
    """ Computes the control points of the rational/non-rational spline after knot removal.

    Implementation based on Algorithm A5.8 and Equation 5.28 of The NURBS Book by Piegl & Tiller

    Keyword Arguments:
        * ``num``: number of knot removals

    :param degree: degree
    :type degree: int
    :param knotvector: knot vector
    :type knotvector: list, tuple
    :param ctrlpts: control points
    :type ctrlpts: list
    :param u: knot to be removed
    :type u: float
    :return: updated control points
    :rtype: list
    """
    tol = kwargs.get('tol', 10e-4)  # Refer to Eq 5.30 for the meaning
    num = kwargs.get('num', 1)  # number of same knot removals
    s = kwargs.get('s', find_multiplicity(u, knotvector))  # multiplicity
    r = kwargs.get('span', find_span_linear(degree, knotvector, len(ctrlpts), u))  # knot span

    # Edge case
    if num < 1:
        return ctrlpts

    # Initialize variables
    first = r - degree
    last = r - s

    # Don't change input variables, prepare new ones for updating
    ctrlpts_new = deepcopy(ctrlpts)

    is_volume = True
    if isinstance(ctrlpts_new[0][0], float):
        is_volume = False

    # Initialize temp array for storing new control points
    if is_volume:
        temp = [[[] for _ in range(len(ctrlpts_new[0]))] for _ in range((2 * degree) + 1)]
    else:
        temp = [[] for _ in range((2 * degree) + 1)]

    # Loop for Eqs 5.28 & 5.29
    for t in range(0, num):
        temp[0] = ctrlpts[first - 1]
        temp[last - first + 2] = ctrlpts[last + 1]
        i = first
        j = last
        ii = 1
        jj = last - first + 1
        remflag = False

        # Compute control points for one removal step
        while j - i >= t:
            alpha_i = knot_removal_alpha_i(u, degree, tuple(knotvector), t, i)
            alpha_j = knot_removal_alpha_j(u, degree, tuple(knotvector), t, j)
            if is_volume:
                for idx in range(len(ctrlpts[0])):
                    temp[ii][idx] = [(cpt - (1.0 - alpha_i) * ti) / alpha_i for cpt, ti
                                     in zip(ctrlpts[i][idx], temp[ii - 1][idx])]
                    temp[jj][idx] = [(cpt - alpha_j * tj) / (1.0 - alpha_j) for cpt, tj
                                     in zip(ctrlpts[j][idx], temp[jj + 1][idx])]
            else:
                temp[ii] = [(cpt - (1.0 - alpha_i) * ti) / alpha_i for cpt, ti in zip(ctrlpts[i], temp[ii - 1])]
                temp[jj] = [(cpt - alpha_j * tj) / (1.0 - alpha_j) for cpt, tj in zip(ctrlpts[j], temp[jj + 1])]
            i += 1
            j -= 1
            ii += 1
            jj -= 1

        # Check if the knot is removable
        if j - i < t:
            if is_volume:
                if linalg.point_distance(temp[ii - 1][0], temp[jj + 1][0]) <= tol:
                    remflag = True
            else:
                if linalg.point_distance(temp[ii - 1], temp[jj + 1]) <= tol:
                    remflag = True
        else:
            alpha_i = knot_removal_alpha_i(u, degree, tuple(knotvector), t, i)
            if is_volume:
                ptn = [(alpha_i * t1) + ((1.0 - alpha_i) * t2) for t1, t2 in zip(temp[ii + t + 1][0], temp[ii - 1][0])]
            else:
                ptn = [(alpha_i * t1) + ((1.0 - alpha_i) * t2) for t1, t2 in zip(temp[ii + t + 1], temp[ii - 1])]
            if linalg.point_distance(ctrlpts[i], ptn) <= tol:
                remflag = True

        # Check if we can remove the knot and update new control points array
        if remflag:
            i = first
            j = last
            while j - i > t:
                ctrlpts_new[i] = temp[i - first + 1]
                ctrlpts_new[j] = temp[j - first + 1]
                i += 1
                j -= 1

        # Update indices
        first -= 1
        last += 1

    # Fix indexing
    t += 1

    # Shift control points (refer to p.183 of The NURBS Book, 2nd Edition)
    j = int((2*r - s - degree) / 2)  # first control point out
    i = j
    for k in range(1, t):
        if k % 2 == 1:
            i += 1
        else:
            j -= 1
    for k in range(i+1, len(ctrlpts)):
        ctrlpts_new[j] = ctrlpts[k]
        j += 1

    # Slice to get the new control points
    ctrlpts_new = ctrlpts_new[0:-t]

    return ctrlpts_new