def fitCurve(points, maxError):
    leftTangent = normalize(points[1] - points[0])
    rightTangent = normalize(points[-2] - points[-1])
    return fitCubic(points, leftTangent, rightTangent, maxError)

>>> import bezier
>>> from fitCurves import fitCurve
>>> import numpy as np
>>> import sys
>>>
>>> points = np.array([[0, 0], [5,50], [50,40], [75,50]], dtype=np.float32)
>>> samples_t = np.linspace(0, 1, 100, dtype=np.float32)
>>> samples_pt = [bezier.q(points, t) for t in samples_t]
>>> calculated_bezier = fitCurve(samples_pt, sys.float_info.max)
>>>
>>> print(calculated_bezier)
[[array([0., 0.], dtype=float32), array([ 4.945499, 45.22236 ], dtype=float32),
array([42.795097, 37.474834], dtype=float32), array([75., 50.], dtype=float32)]]

>>> print(np.array(calculated_bezier) - points)
[[[ 0.          0.        ]
  [-0.05450106 -4.7776413 ]
  [-7.2049026  -2.5251656 ]
  [ 0.          0.        ]]]