def approximate_curve(points, degree, **kwargs):
    """ Curve approximation using least squares method with fixed number of control points.

    Please refer to The NURBS Book (2nd Edition), pp.410-413 for details.

    Keyword Arguments:
        * ``centripetal``: activates centripetal parametrization method. *Default: False*
        * ``ctrlpts_size``: number of control points. *Default: len(points) - 1*

    :param points: data points
    :type points: list, tuple
    :param degree: degree of the output parametric curve
    :type degree: int
    :return: approximated B-Spline curve
    :rtype: BSpline.Curve
    """
    # Number of data points
    num_dpts = len(points)  # corresponds to variable "r" in the algorithm

    # Get keyword arguments
    use_centripetal = kwargs.get('centripetal', False)
    num_cpts = kwargs.get('ctrlpts_size', num_dpts - 1)

    # Dimension
    dim = len(points[0])

    # Get uk
    uk = compute_params_curve(points, use_centripetal)

    # Compute knot vector
    kv = compute_knot_vector2(degree, num_dpts, num_cpts, uk)

    # Compute matrix N
    matrix_n = []
    for i in range(1, num_dpts - 1):
        m_temp = []
        for j in range(1, num_cpts - 1):
            m_temp.append(helpers.basis_function_one(degree, kv, j, uk[i]))
        matrix_n.append(m_temp)

    # Compute NT
    matrix_nt = linalg.matrix_transpose(matrix_n)

    # Compute NTN matrix
    matrix_ntn = linalg.matrix_multiply(matrix_nt, matrix_n)

    # LU-factorization
    matrix_l, matrix_u = linalg.lu_decomposition(matrix_ntn)

    # Initialize control points array
    ctrlpts = [[0.0 for _ in range(dim)] for _ in range(num_cpts)]

    # Fix start and end points
    ctrlpts[0] = list(points[0])
    ctrlpts[-1] = list(points[-1])

    # Compute Rk - Eqn 9.63
    pt0 = points[0]  # Qzero
    ptm = points[-1]  # Qm
    rk = []
    for i in range(1, num_dpts - 1):
        ptk = points[i]
        n0p = helpers.basis_function_one(degree, kv, 0, uk[i])
        nnp = helpers.basis_function_one(degree, kv, num_cpts - 1, uk[i])
        elem2 = [c * n0p for c in pt0]
        elem3 = [c * nnp for c in ptm]
        rk.append([a - b - c for a, b, c in zip(ptk, elem2, elem3)])

    # Compute R - Eqn. 9.67
    vector_r = [[0.0 for _ in range(dim)] for _ in range(num_cpts - 2)]
    for i in range(1, num_cpts - 1):
        ru_tmp = []
        for idx, pt in enumerate(rk):
            ru_tmp.append([p * helpers.basis_function_one(degree, kv, i, uk[idx + 1]) for p in pt])
        for d in range(dim):
            for idx in range(len(ru_tmp)):
                vector_r[i - 1][d] += ru_tmp[idx][d]

    # Computer control points
    for i in range(dim):
        b = [pt[i] for pt in vector_r]
        y = linalg.forward_substitution(matrix_l, b)
        x = linalg.backward_substitution(matrix_u, y)
        for j in range(1, num_cpts - 1):
            ctrlpts[j][i] = x[j - 1]

    # Generate B-spline curve
    curve = BSpline.Curve()
    curve.degree = degree
    curve.ctrlpts = ctrlpts
    curve.knotvector = kv

    return curve
Curve approximation using least squares method with fixed number of control points.

    Please refer to The NURBS Book (2nd Edition), pp.410-413 for details.

    Keyword Arguments:
        * ``centripetal``: activates centripetal parametrization method. *Default: False*
        * ``ctrlpts_size``: number of control points. *Default: len(points) - 1*

    :param points: data points
    :type points: list, tuple
    :param degree: degree of the output parametric curve
    :type degree: int
    :return: approximated B-Spline curve
    :rtype: BSpline.Curve