def getCircularBounds(fitCloud=None,width=64,height=64,smoothing=0.01):
    circumference = 2*(width+height)
    
    if not fitCloud is None:
        cx = np.mean(fitCloud[:,0])
        cy = np.mean(fitCloud[:,1])
        r = 0.5* max( np.max(fitCloud[:,0])- np.min(fitCloud[:,0]),np.max(fitCloud[:,1])- np.min(fitCloud[:,1]))
    else:
        r = circumference /(2.0*math.pi)
        cx = cy = r
    perimeterPoints = np.zeros((circumference,2),dtype=float)
    for i in range(circumference):
        angle = (2.0*math.pi)*float(i) / circumference - math.pi * 0.5 
        perimeterPoints[i][0] = cx + r * math.cos(angle)
        perimeterPoints[i][1] = cy + r * math.sin(angle)
        
        
    bounds = {'top':perimeterPoints[0:width],
              'right':perimeterPoints[width-1:width+height-1],
              'bottom':perimeterPoints[width+height-2:2*width+height-2],
              'left':perimeterPoints[2*width+height-3:]}
    
    bounds['s_top'],u = interpolate.splprep([bounds['top'][:,0], bounds['top'][:,1]],s=smoothing)
    bounds['s_right'],u = interpolate.splprep([bounds['right'][:,0],bounds['right'][:,1]],s=smoothing)
    bounds['s_bottom'],u = interpolate.splprep([bounds['bottom'][:,0],bounds['bottom'][:,1]],s=smoothing)
    bounds['s_left'],u = interpolate.splprep([bounds['left'][:,0],bounds['left'][:,1]],s=smoothing)
   
    
    return bounds

def interpolation_polynom(path,grad):
#    data=np.ndarray(shape=(len(path),3),dtype=float)   #create an array of float type for the input points
#    #fill the array with the Pathdata
#    a=path[0]
#    b=path[1]
#    c=path[2]
#    for i in range(len(a)):
#        data[i,0]=a[i]
#        data[i,1]=b[i]
#        data[i,2]=c[i]
#    #arrange the data to use the function
#    data = data.transpose()
    #interpolate polynom degree 1
    if grad==1:
        tck, u= interpolate.splprep(path,k=1,s=10)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    #interpolate polynom degree 2
    if grad==2:
        tck, u= interpolate.splprep(path,k=2,s=10)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    #interpolate polynom degree 3
    if grad==3:
        tck, u= interpolate.splprep(path, w=None, u=None, ub=None, ue=None, k=3, task=0, s=0.3, t=None, full_output=0, nest=None, per=0, quiet=1)
        path = interpolate.splev(np.linspace(0,1,200), tck)
    return path

def joinJoints(j1, j2):
    assert(j1.dtype==j2.dtype)
    atype = j1.dtype
    param  = np.array([0,0.33, 0.66, 1])



    j1_r  = np.concatenate([j1['r_arm'], [j1['r_gripper']]])
    j2_r  = np.concatenate([j2['r_arm'], [j2['r_gripper']]])
    j1_l  = np.concatenate([j1['l_arm'], [j1['l_gripper']]])
    j2_l  = np.concatenate([j2['l_arm'], [j2['l_gripper']]])
    

    combined_r  = np.concatenate([[j1_r], [0.67*j1_r + 0.33*j2_r], [0.33*j1_r + 0.67*j2_r], [j2_r]])
    combined_l  = np.concatenate([[j1_l], [0.67*j1_l + 0.33*j2_l], [0.33*j1_l + 0.67*j2_l], [j2_l]])
    
    N = 100
    
    (r_tck, _) = si.splprep(combined_r.T, s=0.3, u=param, k=3)                                 
    smooth_r = np.r_[si.splev(np.linspace(0,1,N), r_tck, der=0)].T                               
    
    (l_tck, _) = si.splprep(combined_l.T, s=0.3, u=param, k=3)                                 
    smooth_l   = np.r_[si.splev(np.linspace(0,1,N), l_tck, der=0)].T                               
    

    smooth  = np.zeros(N, dtype=atype)
    for i in xrange(0,N):
        smooth[i]['r_arm'][:]     = smooth_r[i][:-1]
        smooth[i]['r_gripper']    = smooth_r[i][-1]
        smooth[i]['l_arm'][:]     = smooth_l[i][:-1]
        smooth[i]['l_gripper']    = smooth_l[i][-1]
    
    return smooth

def makeSpline(pointList,smPnts):   
    x = [p[0] for p in pointList]
    y = [p[1] for p in pointList] 
        
    xRed = [p[0] for p in smPnts]
    yRed = [p[1] for p in smPnts]

#    print xRed
#    print yRed                                     
    tck,uout = splprep([xRed,yRed],s=0.,k=2,per=False)
    tckOri, uout = splprep([x,y],s=0.,k=2,per=False)
                                
    N=300
                
    uout = list((float(i) / N for i in xrange(N + 1)))
                            
    xOri, yOri = splev(uout,tckOri)                        
    xSp,ySp = splev(uout,tck)         
                
    import dtw
    diff = dtw.dynamicTimeWarp(zip(xOri,yOri), zip(xSp,ySp))
                
    err =  diff/len(xSp)
        
    return tck,err

def muDDot(time,mu):
    tck,uout = interpolate.splprep([time,mu],s=0.,k=2,per=False)
    dx,dy = interpolate.splev(uout,tck,der=1)
    mudot = dy/dx
    tck,uout = interpolate.splprep([time,mudot],s=0.,k=2,per=False)
    ddx,ddy = interpolate.splev(uout,tck,der=1)
    muddot = ddy/ddx
    return muddot

def midcrossings(a, b=None, thresh=1e-3, k=5):
    """usage res = midcrossings([x,] y)

    returns fwhm of y (a function with a single maximum value)
    from spline-interpolated midpoint crossings"""
    if b == None:
        y = array(a)
        x = arange(y)
    else:
        x = array(a)
        y = array(b)
    try:
        assert x.shape[0] == y.shape[0]
    except AssertionError:
        print "x and y must be same length"
        return None
    maxind = where(y == y.max())[0].flatten()[0]  # uses only first max pt

    (y1, y2) = y[:maxind], y[maxind:]
    (x1, x2) = x[:maxind], x[maxind:]

    s = 1.0
    nest = -1
    interpvals = linspace(0, 1, 251)
    # 251 simply gives enough points to give one point close to 0.5

    print "thresholding to:", thresh
    # lower half
    nob1 = where(y1 > thresh)
    # ^ need to ignore baseline values when fitting splines
    y1tofit = y1[nob1] / y1[nob1].max()

    tckp, u = splprep([y1tofit, x1[nob1]], s=s, k=k, nest=nest)
    y1p, x1p = splev(interpvals, tckp)
    dtohalf = abs(y1p - 0.5)
    # 0.5 because want width at _half_ max
    closest = where(dtohalf == dtohalf.min())
    lowval = x1p[closest]

    # upper half
    nob2 = where(y2 > thresh)
    y2tofit = y2[nob2] / y2[nob2].max()
    tckp, u = splprep([y2tofit, x2[nob2]], s=s, k=k, nest=nest)
    y2p, x2p = splev(interpvals, tckp)
    dtohalf = abs(y2p - 0.5)
    closest = where(dtohalf == dtohalf.min())
    hival = x2p[closest]

    if graphic:
        fwhm = hival - lowval
        figure(2)
        clf()
        plot(x, y / y.max(), "bx-", label="original")
        plot(x1p, y1p, "r-", x2p, y2p, "r-", label="spline fit")
        plot((lowval, hival), (0.5, 0.5), "k-", label="width")

    return fwhm

    def interpolate(self, fit_x, fit_y, smoothing = None, newu = None):
        if smoothing is not None:
            tck, u = interpolate.splprep([fit_x, fit_y], k=4, s=smoothing)
        else:
            tck, u = interpolate.splprep([fit_x, fit_y], k=4)

        if (newu is not None):
            u = newu

        out = interpolate.splev(u, tck)
        return out, u, tck

 def __init__(self, pts, V, dV=None, V_spline_samples=100,
              extend_to_minima=False, reeval_distances=True):
     assert len(pts) > 1
     # 1. Find derivs
     dpts = _pathDeriv(pts)
     # 2. Extend the path
     if extend_to_minima:
         def V_lin(x, p0, dp0, V): return V(p0+x*dp0)
         # extend at the front of the path
         xmin = optimize.fmin(V_lin, 0.0, args=(pts[0], dpts[0], V), 
                              xtol=1e-6, disp=0)[0]
         if xmin > 0.0: xmin = 0.0
         nx = np.ceil(abs(xmin)-.5) + 1
         x = np.linspace(xmin, 0, nx)[:, np.newaxis]
         pt_ext = pts[0] + x*dpts[0]
         pts = np.append(pt_ext, pts[1:], axis=0)
         # extend at the end of the path
         xmin = optimize.fmin(V_lin, 0.0, args=(pts[-1], dpts[-1], V), 
                              xtol=1e-6, disp=0)[0]
         if xmin < 0.0: xmin = 0.0
         nx = np.ceil(abs(xmin)-.5) + 1
         x = np.linspace(xmin, 0, nx)[::-1, np.newaxis]
         pt_ext = pts[-1] + x*dpts[-1]
         pts = np.append(pts[:-1], pt_ext, axis=0)
         # Recalculate the derivative
         dpts = _pathDeriv(pts)
     # 3. Find knot positions and fit the spline.
     pdist = integrate.cumtrapz(np.sqrt(np.sum(dpts*dpts, axis=1)), 
                                initial=0.0)
     self.L = pdist[-1]
     k = min(len(pts)-1, 3) # degree of the spline
     self._path_tck = interpolate.splprep(pts.T, u=pdist, s=0, k=k)[0]
     # 4. Re-evaluate the distance to each point.
     if reeval_distances:
         def dpdx(_, x): 
             dp = np.array(interpolate.splev(x, self._path_tck, der=1))
             return np.sqrt(np.sum(dp*dp))
         pdist = integrate.odeint(dpdx, 0., pdist,
                                  rtol=0, atol=pdist[-1]*1e-8)[:,0]
         self.L = pdist[-1]
         self._path_tck = interpolate.splprep(pts.T, u=pdist, s=0, k=k)[0]
     # Now make the potential spline.
     self._V = V
     self._dV = dV
     self._V_tck = None
     if V_spline_samples is not None:
         x = np.linspace(0,self.L,V_spline_samples)
         # extend 20% beyond this so that we more accurately model the
         # path end points
         x_ext = np.arange(x[1], self.L*.2, x[1])
         x = np.append(-x_ext[::-1], x)
         x = np.append(x, self.L+x_ext)
         y = self.V(x)
         self._V_tck = interpolate.splrep(x,y,s=0)

def correct_te(tck, k):
    """Corrects the trailing edge of a flatback airfoil.

    This corrections will make the trailing edge of the normalized flatback
    airfoil align with the y-axis.

    Args:
        tck (tuple): A tuple (t,c,k) containing the vector of knots, the
            B-spline coefficients, and the degree of the spline.
        k (int): The degree of the returned bspline

    Return:
        tuple: A tuple (t,c,k) containing the vector of knots, the
            B-spline coefficients, and the degree of the spline.

    """
    try:
        u0_x = bspl_find_x(x_loc=1.0, start=0.0, end=0.1, tck=tck)
    except ValueError:
        u0_x = None
    try:
        u1_x = bspl_find_x(x_loc=1.0, start=0.9, end=1.0, tck=tck)
    except ValueError:
        u1_x = None

    if u0_x is not None and u1_x is not None:
        u = np.linspace(u0_x, u1_x, 1000)
        points = interpolate.splev(u, tck, der=0)
        tck_norm_mod = interpolate.splprep(points, s=0.0, k=k)
    elif u0_x is None and u1_x is not None:
        u = np.linspace(0.0, u1_x, 1000)
        points = interpolate.splev(u, tck, der=0)
        p_u0 = [points[0][0], points[1][0]]
        u0_grad = interpolate.splev(0.0, tck, der=1)
        dx = 1.0 - p_u0[0]
        dy = dx * u0_grad[1] / u0_grad[0]
        p_new = [1.0, p_u0[1] + dy]
        x_pts = np.insert(points[0], 0, p_new[0])
        y_pts = np.insert(points[1], 0, p_new[1])
        tck_norm_mod, _ = interpolate.splprep([x_pts, y_pts], s=0.0, k=k)
    elif u0_x is not None and u1_x is None:
        u = np.linspace(u0_x, 1.0, 1000)
        points = interpolate.splev(u, tck, der=0)
        p_u1 = [points[0][-1], points[1][-1]]
        u1_grad = interpolate.splev(1.0, tck, der=1)
        dx = 1.0 - p_u1[0]
        dy = dx * u1_grad[1] / u1_grad[0]
        p_new = [1.0, p_u1[1] + dy]
        x_pts = np.append(points[0], p_new[0])
        y_pts = np.append(points[1], p_new[1])
        tck_norm_mod, _ = interpolate.splprep([x_pts, y_pts], s=0.0, k=k)
    else:
        raise ValueError('Something is wrong with the bspline!')
    return tck_norm_mod

    def test_splprep(self):
        x = np.arange(15).reshape((3, 5))
        b, u = splprep(x)
        tck, u1 = _impl.splprep(x)

        # test the roundtrip with splev for both "old" and "new" output
        assert_allclose(u, u1, atol=1e-15)
        assert_allclose(splev(u, b), x, atol=1e-15)
        assert_allclose(splev(u, tck), x, atol=1e-15)

        # cover the ``full_output=True`` branch
        (b_f, u_f), _, _, _ = splprep(x, s=0, full_output=True)
        assert_allclose(u, u_f, atol=1e-15)
        assert_allclose(splev(u_f, b_f), x, atol=1e-15)

    def test_splprep_errors(self):
        # test that both "old" and "new" code paths raise for x.ndim > 2
        x = np.arange(3*4*5).reshape((3, 4, 5))
        with assert_raises(ValueError, message="too many values to unpack"):
            splprep(x)
        with assert_raises(ValueError, message="too many values to unpack"):
            _impl.splprep(x)

        # input below minimum size
        x = np.linspace(0, 40, num=3)
        with assert_raises(TypeError, message="m > k must hold"):
            splprep([x])
        with assert_raises(TypeError, message="m > k must hold"):
            _impl.splprep([x])

        # automatically calculated parameters are non-increasing
        # see gh-7589
        x = [-50.49072266, -50.49072266, -54.49072266, -54.49072266]
        with assert_raises(ValueError, message="Invalid inputs"):
            splprep([x])
        with assert_raises(ValueError, message="Invalid inputs"):
            _impl.splprep([x])

        # given non-increasing parameter values u
        x = [1, 3, 2, 4]
        u = [0, 0.3, 0.2, 1]
        with assert_raises(ValueError, message="Invalid inputs"):
            splprep(*[[x], None, u])

 def _smooth_segs(self,seg,kind,avg=False):
     useg,per=unique_vec(seg)
     if len(useg)>self.k: #check if long enough to smooth
         tck,u=splprep(useg.T,s=self.smooth,k=self.k,per=per)
         seg=np.vstack(splev(self.unew,tck)).T
         if (self.smooth>0) and (avg):
             #when smooth>0 you get different results if you reverse the order of the points
             #do it both ways and take the average of the two solutions
             #and return that instead (needed for filled contours to line up correctly)
             tck2,u2=splprep(useg[::-1].T,s=self.smooth,k=self.k,per=per)
             seg2=np.vstack(splev(self.unew[::-1],tck2)).T
             seg=np.dstack([seg,seg2]).mean(axis=2)
         if kind is not None:
             kind=self.kinds_fill
     return seg,kind

def tract_prototype_mean(tractography, smooth_order, file_output=None):
    from .tract_obb import prototype_tract

    tracts = tractography.tracts()
    prototype_ix, leave_centers = prototype_tract(tracts, return_leave_centers=True)

    median_tract = tracts[prototype_ix]

    mean_tract = numpy.empty_like(median_tract)
    centers_used = set()
    for point in median_tract:
        closest_leave_center_ix = (
            ((leave_centers - point[None, :]) ** 2).sum(1)
        ).argmin()

        if closest_leave_center_ix in centers_used:
            continue

        mean_tract[len(centers_used)] = leave_centers[closest_leave_center_ix]
        centers_used.add(closest_leave_center_ix)

    mean_tract = mean_tract[:len(centers_used)]

    if smooth_order > 0:
        try:
            from scipy import interpolate

            tck, u = interpolate.splprep(mean_tract.T)
            mean_tract = numpy.transpose(interpolate.splev(u, tck))
        except ImportError:
            warn("A smooth order larger than 0 needs scipy installed")

    return Tractography([mean_tract], {}, **tractography.extra_args)

def unif_resample(x,n=None,tol=0,deg=None, seg_len = .02):
 
    if deg is None: deg = min(3, len(x) - 1)

    x = np.atleast_2d(x)
    x = remove_duplicate_rows(x)
    
    (tck,_) = si.splprep(x.T,k=deg,s = tol**2*len(x),u=np.linspace(0,1,len(x)))
    xup = np.array(si.splev(np.linspace(0,1, 10*len(x),.1),tck)).T
    dl = norms(xup[1:] - xup[:-1],1)
    l = np.cumsum(np.r_[0,dl])
    (tck,_) = si.splprep(xup.T,k=deg,s = tol**2*len(xup),u=l)

    if n is not None: newu = np.linspace(0,l[-1],n)
    else: newu = np.linspace(0, l[-1], l[-1]//seg_len)
    return np.array(si.splev(newu,tck)).T    

    def get_measurement_lines(self, tail):
        """
        determines the measurement segments that are used for the line scan
        """
        f_c = self.params['measurement/line_offset']
        f_o = 1 - f_c

        centerline = tail.centerline
        result = []
        for side in tail.sides:
            # find the line between the centerline and the ventral line
            points = []
            for p_c in centerline:
                p_o = curves.get_projection_point(side, p_c) #< outer line
                points.append((f_c*p_c[0] + f_o*p_o[0],
                               f_c*p_c[1] + f_o*p_o[1]))
                
            # do spline fitting to smooth the line
            smoothing = self.params['measurement/spline_smoothing']*len(points)
            tck, _ = interpolate.splprep(np.transpose(points),
                                         k=2, s=smoothing)
            
            points = interpolate.splev(np.linspace(-0.5, .8, 100), tck)
            points = zip(*points) #< transpose list
    
            # restrict centerline to object
            mline = geometry.LineString(points).intersection(tail.polygon)
            
            # pick longest line if there are many due to strange geometries
            if isinstance(mline, geometry.MultiLineString):
                mline = mline[np.argmax([l.length for l in mline])]
                
            result.append(np.array(mline.coords))
            
        return result   

def spline_fitted_magnitudes_brute(times, magnitudes, errors, requested_times):
    # Spline parameters:
    s=len(times)/100. # smoothness parameter
    k=5 # spline order
    nest=-1 # estimate of number of knots needed (-1 = maximal)

    # Find the knot points.
    tckp, u = splprep([times, magnitudes],s=s,k=k,nest=-1)

    # Evaluate spline, including interpolated points.
    new_times, new_magnitudes = splev(linspace(0,1,len(times)),tckp)

    # Define an interpolating function along the spline fit.
    interpolating_function = interp1d(new_times, new_magnitudes, kind = "linear")

    # Interpolate linerarly along the spline at the requested times.
    fitted_magnitudes = interpolating_function(requested_times)

    fitted_errors = array([])
    for n in range(len(requested_times)):
        for m in range(len(times)-1):
            if (requested_times[n] > times[m]) and (requested_times[n] < times[m+1]):
                error = (errors[m] + errors[m+1])*0.5
                fitted_errors = append(fitted_errors, error)
    return fitted_magnitudes, fitted_errors

def second_derivative(xdata, inds, gt=False, s=0):
    '''
    The second derivative of d^2 xdata / d inds^2

    why inds for interpolation, not log l?
    if not using something like model number instead of log l,
    the tmin will get hidden by data with t < tmin but different
    log l. This is only a problem for very low Z.
    If I find the arg min of teff to be very close to MS_BEG it
    probably means the MS_BEG is at a lower Teff than Tmin.
    '''
    tckp, _ = splprep([inds, xdata], s=s, k=3)
    arb_arr = np.arange(0, 1, 1e-2)
    xnew, ynew = splev(arb_arr, tckp)
    # second derivative, bitches.
    ddxnew, ddynew = splev(arb_arr, tckp, der=2)
    ddyddx = ddynew / ddxnew
    # not just argmin, but must be actual min...
    try:
        if gt:
            aind = [a for a in np.argsort(ddyddx) if ddyddx[a-1] < 0][0]
        else:
            aind = [a for a in np.argsort(ddyddx) if ddyddx[a-1] > 0][0]
    except IndexError:
        return -1
    tmin_ind, _ = closest_match2d(aind, inds, xdata, xnew, ynew)
    return inds[tmin_ind]

def interpolation_polynom(path,grad):
    (x,y)=path.shape
    anzahl=y*40
    #interpolate polynom degree 1
    if grad==1:
        tck, u= interpolate.splprep(path,k=1,s=0.2)
        path = interpolate.splev(np.linspace(0,1,anzahl), tck)
    #interpolate polynom degree 2
    if grad==2:
        tck, u= interpolate.splprep(path,k=2,s=0.2)
        path = interpolate.splev(np.linspace(0,1,anzahl), tck)
    #interpolate polynom degree 3
    if grad==3:
        tck, u= interpolate.splprep(path, w=None, u=None, ub=None, ue=None, k=3, task=0, s=0.2, t=None, full_output=0, nest=None, per=0, quiet=1)
        path = interpolate.splev(np.linspace(0,1,anzahl), tck)
    return path

def make_joint_trajectory_with_limits(positions, vel_limits, acc_limits):
    positions = np.asarray(positions)
    vel_limits = np.asarray(vel_limits)
    acc_limits = np.asarray(acc_limits)
    n_waypoints, n_joints = positions.shape

    # estimate the time, to pick a reasonable number of samples
    t_est = (abs(positions[1:] - positions[:-1]).sum(axis=0) / vel_limits).max()
    # samples is an integer multiple of n_waypoints, so that we hit all waypoints
    upsample_ratio = max(1, int(np.ceil(t_est * 10 / n_waypoints)))
    n_samples = n_waypoints * upsample_ratio

    # upsample and smooth a little bit
    k = min(3, n_waypoints - 1)     
    (tck, _) = si.splprep(positions.T, s = .001**2*n_waypoints, u=linspace(0,1,n_waypoints), k=k) # todo: is s total or per waypoint?
    sampled_positions = r_[si.splev(linspace(0,1,n_samples),tck)].T


    velocities = np.zeros((n_samples, n_joints))
    times = np.zeros(n_samples)
    
    for i in xrange(1,n_samples):
        
        dpos = (sampled_positions[i] - sampled_positions[i-1])
        # amount of time if we're at velocity speed limit:
        dt_vel = norm(dpos / vel_limits, inf)

        # search for minimal dt that satisfies velocity and acceleration limits
        f = lambda dt: (abs(dpos/dt - velocities[i-1])/dt < acc_limits).all()
        dt = line_search(f, dt_vel)
        
        times[i] = times[i-1]+dt
        velocities[i] = dpos/dt

    return sampled_positions, velocities, times

 def rotate_te(self, alpha, nsamples, te_smooth=1, smoothing=0.0, degree=5,
               ins_pt=None, plot=False):
     """Splits the airfoil in two parts at the leading edge, rotates both
     parts by alpha/2.0 around the trailing edge and than reconnects both
     parts again."""
     te_point = self.get_te_point()
     u_le, le_point = self.get_le_point()
     te_smoothing = 1.0*te_smooth/100.0
     u0 = 0.0
     u1 = u_le - te_smoothing
     u2 = u_le + te_smoothing
     u3 = 1.0
     ss_pts = self._rotate_around_point(alpha=-alpha/2.0, u0=u0, u1=u1,
                                        rot_pt=te_point, nsamples=nsamples)
     ps_pts = self._rotate_around_point(alpha=alpha/2.0, u0=u2, u1=u3,
                                        rot_pt=te_point, nsamples=nsamples)
     if ins_pt is None:
         new_pts = np.vstack((ss_pts, ps_pts))
     else:
         new_pts = np.vstack((ss_pts, ins_pt, ps_pts))
     if plot:
         plt.figure('Leading edge rotation')
         plt.plot(new_pts[:, 0], new_pts[:, 1], 'or', label='fit points')
         plt.axis('equal')
         plt.grid(True)
         plt.legend()
         plt.show()
     x = [new_pts[:, 0], new_pts[:, 1]]
     self.tck, _ = interpolate.splprep(x, s=smoothing, k=degree)
     return ss_pts, ps_pts