python类splprep()的实例源码

def resample_nd(curve, npoints, smooth = 0, periodic = False, derivative = 0):
  """Resample n points using n equidistant points along a curve

  Arguments:
    points (mxd array): coordinate of the reference points for the curve
    npoints (int): number of resamples equidistant points
    smooth (number): smoothness factor
    periodic (bool): if True assumes the curve is a closed curve

  Returns:
    (nxd array): resampled equidistant points
  """

  cinterp, u = splprep(curve.T, u = None, s = smooth, per = periodic);
  if npoints is all:
    npoints = curve.shape[0];
  us = np.linspace(u.min(), u.max(), npoints)
  curve = splev(us, cinterp, der = derivative);
  return np.vstack(curve).T;

def resample_spline(points, smooth=.001, count=None):
    from scipy.interpolate import splprep, splev
    if count is None: 
        count = len(points)    
    points = np.asanyarray(points)
    closed = np.linalg.norm(points[0] - points[-1]) < tol.merge 

    tpl = splprep(points.T, s=smooth)[0]
    i = np.linspace(0.0, 1.0, count)
    resampled = np.column_stack(splev(i, tpl))

    if closed:
        shared = resampled[[0,-1]].mean(axis=0)
        resampled[0]  = shared
        resampled[-1] = shared

    return resampled

def points_to_spline_entity(points, smooth=.0005, count=None):
    from scipy.interpolate import splprep

    if count is None: 
        count = len(points)
    points = np.asanyarray(points)
    closed = np.linalg.norm(points[0] - points[-1]) < tol.merge 

    knots, control, degree =  splprep(points.T, s=smooth)[0]
    control = np.transpose(control)
    index = np.arange(len(control))

    if closed:
        control[0] = control[[0,-1]].mean(axis=0)
        control = control[:-1]
        index[-1] = index[0]

    entity = BSpline(points = index, 
                     knots  = knots,
                     closed = closed)

    return entity, control

def interpolate_path(list_of_xys):
    x = list_of_xys[0]
    y = list_of_xys[1]
    #print x, y

    #x = np.array([0,1,2,3,4,5,6,7,6,5,4,3,2,1,0])
    #y = np.array([0,1,2,3,4,4,5,6,7,8,7,8,9,9,9])

    tck, u = interpolate.splprep([x, y], s = 0.03)
    unew = np.arange(0,9,0.03)
    #print unew
    interpolated_path = interpolate.splev(unew, tck, ext = 1)
    return interpolated_path

#a = interpolate_path(array)
#plot_path(array, a)

def simple_poly_fit(x , y , number_of_uni_point, smoothness):

    number_of_point = len(x)
    t = np.zeros_like(x)
    for i in range(0 , number_of_point):
        if i > 0:
            t[i] = t[i - 1] + np.sqrt((x[i] - x[i - 1]) ** 2 + (y[i] - y[i - 1]) ** 2)
        else:
            t[i] = 0

    k = min(3 , number_of_point - 1) # spline order
    nest = -1 # estimate of number of knots needed (-1 = maximal)
    tckp , u = splprep([x , y , t] , s = smoothness , k = k , nest = -1)

    x_new , y_new, t_new = splev(linspace(0,1,number_of_uni_point),tckp)

    return x_new , y_new

########################################################################
########################################################################

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 _interpolate(xy, num_points):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = linspace(0, 1, num_points)
  out = splev(unew, tck)
  return column_stack(out)

def _rnd_interpolate(xy, num_points, ordered=False):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = random(num_points)
  if ordered:
    unew = sort(unew)
  out = splev(unew, tck)
  return column_stack(out)

def getBoundaryPoints(x , y):
    tck,u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 10000)
    xnew,ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int),ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return coord[:,0],coord[:,1]

def getBoundaryPoints(x = [], y = []):
    tck,u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 1000)
    xnew,ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int),ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return coord[:,0],coord[:,1]

def plot_contour(self, xv, yv, cost_grid):
        """
        Function constructs contour lines
        """
        contour = Contour(xv, yv, cost_grid)
        contour_lines = contour.contours(
            np.linspace(np.min(cost_grid), np.max(cost_grid), 20))

        series = []
        count = 0
        for key, value in contour_lines.items():
            for line in value:
                if len(line) > 3:
                    tck, u = splprep(np.array(line).T, u=None, s=0.0, per=0)
                    u_new = np.linspace(u.min(), u.max(), 100)
                    x_new, y_new = splev(u_new, tck, der=0)
                    interpol_line = np.c_[x_new, y_new]
                else:
                    interpol_line = line

                series.append(dict(data=interpol_line,
                                   color=self.contour_color,
                                   type="spline",
                                   lineWidth=0.5,
                                   marker=dict(enabled=False),
                                   name="%g" % round(key, 2),
                                   enableMouseTracking=False
                                   ))
                count += 1
        return series

def callable_from_trajectory(t, curves):
    """
    Use scipy.interpolate splprep to build cubic b-spline interpolating
    functions over a set of curves.

    Parameters
    ----------
    t : 1D array_like
        Array of m time indices of trajectory
    curves : 2D array_like
        Array of m x n vector samples at the time indices. First dimension
        indexes time, second dimension indexes vector components

    Returns
    -------
    interpolated_callable : callable
        Callable which interpolates the given curve/trajectories
    """
    tck_splprep = interpolate.splprep(
        x=[curves[:, i] for i in range(curves.shape[1])], u=t, s=0)

    def interpolated_callable(t, *args):
        return np.array(
            interpolate.splev(t, tck_splprep[0], der=0)
        ).T.squeeze()

    return interpolated_callable

def get_boundary_points(x, y):
    tck, u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 1000)
    xnew, ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int), ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return np.array(coord[:, 0], dtype=np.int32), np.array(coord[:, 1], dtype=np.int32)

def get_boundary_points(x, y):
    tck, u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 1000)
    xnew, ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int), ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return np.array(coord[:, 0], dtype=np.int32), np.array(coord[:, 1], dtype=np.int32)

def _interpolate(xy, num_points):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = linspace(0, 1, num_points)
  out = splev(unew, tck)
  return column_stack(out)

def _rnd_interpolate(xy, num_points, ordered=False):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = random(num_points)
  if sort:
    unew = sort(unew)
  out = splev(unew, tck)
  return column_stack(out)

def set_points_uniform_length(self, npoints = None):
    """Resample the curves base points such that the distances between the curve vectors are uniform

    Arguments:
      npoints (int or None): number of sample points    
    """
    if npoints is None:
      if self._npoints is None:
        raise ValueError('cannot determine number of points for uniform sampling!');
      npoints = self._npoints; 

    tck,u = splprep(self.values.T, s = 0);
    points = np.linspace(0,1,npoints)
    values = splev(points, tck, der = 0, ext = 1);  
    self.set_values(np.vstack(values).T, points);

def tck(self, dimension = None):
    """Returns a tck tuple for use with spline functions

    Arguments:
      dimension (None, list or int): dimension(s) for which to return tck tuple, if None return for all

    Returns:
      tuple: tck couple as returned by splprep
    """
    if dimension is None:
      return (self._knots_all, self._parameter.T, self.degree);
    else:
      if isinstance(dimension,int):
        dimension = [dimension];
      return (self._knots_all, self._parameter[:,dimension].T, self.degree);

def move_forward_center_discrete(distance, center, straight = True):
  """Move worm forward peristaltically

  Arguments:
    distance (float): distance to move forward
    center (nx2 array): center points
    length (float or None): length to use for position update
    straight (bool): if True extrapolated points move straight

  Note:
    The head is first point in center line and postive distances will move the
    worm in this direction.
  """  
  cinterp, u = splprep(center.T, u = None, s = 0, per = 0)
  us = u - distance;
  x, y = splev(us, cinterp, der = 0); 
  cline2 = np.array([x,y]).T;   

  if straight:
     l = length_from_center_discrete(center);
     if distance > 0:
       idx = np.where(us < 0)[0];
       if len(idx) > 0:
         d = center[0,:] - center[1,:];
         d = d / np.linalg.norm(d) * l;
         for i in idx:
           cline2[i,:] = center[0,:] - d * us[i];
     elif distance < 0:
       idx = np.where(us > 1)[0];
       if len(idx) > 0:
         d = center[-1,:] - center[-2,:];
         d = d / np.linalg.norm(d) * l;
         for i in idx:
           cline2[i,:] = center[-1,:] + d * (us[i]-1);

  return cline2;

def _interpolate(xy, num_points):
    tck, u = splprep([
        xy[:, 0],
        xy[:, 1]],
        s=0
    )
    unew = linspace(0, 1, num_points)
    out = splev(unew, tck)
    return column_stack(out)

def _rnd_interpolate(xy, num_points, ordered=False):
    tck, u = splprep([
        xy[:, 0],
        xy[:, 1]],
        s=0
    )
    unew = random(num_points)
    if sort:
        unew = sort(unew)
    out = splev(unew, tck)
    return column_stack(out)

def _interpolate_loops(self, ds):
        """
        Interpolate all loops to a resolution (`ds`) below the minimum bin width of all of the
        instruments. This ensures that the image isn't 'patchy' when it is binned.
        """
        if ds is None:
            ds = 0.1*np.min([min(instr.resolution.x.value, instr.resolution.y.value)
                             for instr in self.instruments])*self.instruments[0].resolution.x.unit
        ds = self.field._convert_angle_to_length(ds)
        # FIXME: memory requirements for this list will grow with number of loops, consider saving 
        # it to the instrument files, both the interpolated s and total_coordinates
        total_coordinates = []
        interpolated_loop_coordinates = []
        for loop in self.field.loops:
            self.logger.debug('Interpolating loop {}'.format(loop.name))
            n_interp = int(np.ceil(loop.full_length/ds))
            interpolated_s = np.linspace(loop.field_aligned_coordinate.value[0],
                                         loop.field_aligned_coordinate.value[-1], n_interp)
            interpolated_loop_coordinates.append(interpolated_s)
            nots, _ = splprep(loop.coordinates.value.T)
            _tmp = splev(np.linspace(0, 1, n_interp), nots)
            total_coordinates += [(x, y, z) for x, y, z in zip(_tmp[0], _tmp[1], _tmp[2])]

        total_coordinates = np.array(total_coordinates)*loop.coordinates.unit

        return total_coordinates, interpolated_loop_coordinates

def _interpolate(xy, num_points):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = linspace(0, 1, num_points)
  out = splev(unew, tck)
  return column_stack(out)

def _rnd_interpolate(xy, num_points, ordered=False):
  tck,u = splprep([
    xy[:,0],
    xy[:,1]],
    s=0
  )
  unew = random(num_points)
  if ordered:
    unew = sort(unew)
  out = splev(unew, tck)
  return column_stack(out)

def plot_contour(self):
        """
        Function constructs contour lines
        """
        self.scatter.remove_contours()
        if self.contours_enabled:
            is_tree = type(self.learner) in [RandomForestLearner, TreeLearner]
            # tree does not need smoothing
            contour = Contour(
                self.xv, self.yv, self.probabilities_grid
                if is_tree else self.blur_grid(self.probabilities_grid))
            contour_lines = contour.contours(
                np.hstack(
                    (np.arange(0.5, 0, - self.contour_step)[::-1],
                     np.arange(0.5 + self.contour_step, 1, self.contour_step))))
            # we want to have contour for 0.5

            series = []
            count = 0
            for key, value in contour_lines.items():
                for line in value:
                    if (len(line) > self.degree and
                                type(self.learner) not in
                                [RandomForestLearner, TreeLearner]):
                        # if less than degree interpolation fails
                        tck, u = splprep(
                            [list(x) for x in zip(*reversed(line))],
                            s=0.001, k=self.degree,
                            per=(len(line)
                                 if np.allclose(line[0], line[-1])
                                 else 0))
                        new_int = np.arange(0, 1.01, 0.01)
                        interpol_line = np.array(splev(new_int, tck)).T.tolist()
                    else:
                        interpol_line = line

                    series.append(dict(data=self.labeled(interpol_line, count),
                                       color=self.contour_color,
                                       type="spline",
                                       lineWidth=0.5,
                                       showInLegend=False,
                                       marker=dict(enabled=False),
                                       name="%g" % round(key, 2),
                                       enableMouseTracking=False
                                       ))
                    count += 1
            self.scatter.add_series(series)
        self.scatter.redraw_series()

def _get_spline(cls, points, pix_err=2, need_sort=True, **kwargs):
        '''
        Returns a closed spline for the points handed in.
        Input is assumed to be a (2xN) array

        =====
        input
        =====

        :param points: the points to fit the spline to
        :type points: a 2xN ndarray or a list of len =2 tuples

        :param pix_err: the error is finding the spline in pixels
        :param need_sort: if the points need to be sorted
            or should be processed as-is

        =====
        output
        =====
        tck
           The return data from the spline fitting
        '''
        if type(points) is np.ndarray:
            # make into a list
            pt_lst = zip(*points)
            # get center
            center = np.mean(points, axis=1).reshape(2, 1)
        else:
            # make a copy of the list
            pt_lst = list(points)

            # compute center
            tmp_fun = lambda x, y: (x[0] + y[0], x[1] + y[1])

            center = np.array(reduce(tmp_fun, pt_lst)).reshape(2, 1)
            center /= len(pt_lst)
        if len(pt_lst) < 5:
            raise TooFewPointsException("not enough points")

        if need_sort:
            # sort the list by angle around center
            pt_lst.sort(key=lambda x: np.arctan2(x[1] - center[1],
                                                 x[0] - center[0]))
        # add first point to end because it is periodic (makes the
        # interpolation code happy)
        pt_lst.append(pt_lst[0])
        # make array for handing in to spline fitting
        pt_array = np.vstack(pt_lst).T
        # do spline fitting

        tck, u = si.splprep(pt_array, s=len(pt_lst) * (pix_err ** 2), per=True)
        return tck

def measure(self):
    """Measure some properties at the same time"""

    if self.valid:
      cl = self.center_line();
      n2 = self.center_index();

      # positions
      pos_head = cl[0];
      pos_tail = cl[1];
      pos_center = cl[n2];


      # head tail distance:
      dist_head_tail = np.linalg.norm(pos_head-pos_tail)
      dist_head_center = np.linalg.norm(pos_head-pos_center)
      dist_tail_center = np.linalg.norm(pos_tail-pos_center)

      #average curvature
      xymintp, u = splprep(cl.T, u = None, s = 1.0, per = 0);
      dcx, dcy = splev(u, xymintp, der = 1)
      d2cx, d2cy = splev(u, xymintp, der = 2)
      ck = (dcx * d2cy - dcy * d2cx)/np.power(dcx**2 + dcy**2, 1.5);
      curvature_mean = np.mean(ck);
      curvature_variation = np.sum(np.abs(ck))

      curled = False;
    else:
      pos_head = np.array([0,0]);
      pos_tail = pos_head;
      pos_center = pos_head;


      # head tail distance:
      dist_head_tail = 0.0
      dist_head_center = 0.0
      dist_tail_center = 0.0

      #average curvature
      curvature_mean = 0.0;
      curvature_variation = 0.0

      curled = True;

    data = np.hstack([pos_head, pos_tail, pos_center, dist_head_tail, dist_head_center, dist_tail_center, curvature_mean, curvature_variation, curled]);
    return data


  ############################################################################
  ### Worm shape deformations, Worm motions

def move_forward(self, distance, smooth = 1.0, straight = True):
    """Move worm peristaltically forward

    Arguments:
      distance (number): distance to move forward
      smooth (number): smoothness of the interpolation
      straight (bool): if True extrapolated points move straight

    Note:
      The head is first point in center line and postive distances will move the
      worm in this direction.
    """

    length = self.length;
    cline = self.center_line();

    cinterp, u = splprep(cline.T, u = None, s = smooth, per = 0)
    us = u - distance / length;
    x, y = splev(us, cinterp, der = 0); 
    cline2 = np.array([x,y]).T;   
    if straight:
       if distance > 0:
         idx = np.nonzero(us < 0)[0];
         if len(idx) > 0:
           d = cline[0,:] - cline[1,:];
           #l = np.linalg.norm(d);
           l = self.length / (self.npoints -1);
           d = d / l;
           m = idx.max();
           for i in idx:
             cline2[i,:] = cline[0,:] + distance * d * (m + 1.0 - i)/(m + 1.0);
       elif distance < 0:
         idx = np.nonzero(us > 1)[0];
         if len(idx) > 0:
           d = cline[-1,:] - cline[-2,:];
           #l = np.linalg.norm(d);
           l = self.length / (self.npoints -1);
           d = d / l;
           m = idx.min(); mx = idx.max();
           for i in idx:
             cline2[i,:] = cline[-1,:] - distance * d * (i - m + 1.0)/(mx + 1.0 - m);

    self.from_center_line(cline2, self.width);
    self.stretch(length / self.length);

def set_values(self, values, points = None, dimension = None):
    """Calcualte the bspline parameter for the data points y

    Arguments:
      values (array): values of data points
      points (array or None): sample points for the data values, if None use internal points or linspace(0,1,values.shape[0])
      dimension (int, list or None): the dimension(s) at which to change the curve, if None change dimension to values.shape[0]
    """
    values = np.array(values, dtype = float);
    if values.ndim == 1:
        values = values[:,np.newaxis];
    vdims = range(values.shape[1]);

    # determine the dimensions at which to change curve
    if dimension is None:
      pdims = range(values.shape[1]);
      self.ndim = values.shape[1];
    else:
      pdims = np.array(dimension, dtype = int);

    if len(vdims) != len(pdims) or len(pdims) != values.shape[1] or max(pdims) > self.ndim:
      raise RuntimeError('inconsistent number of dimensions %d, values %d and parameter %d and curve %d' % (values.shape[1], len(vdims), len(pdims), self.ndim));

    #set points
    if points is None:
      if self._points is None:
        self.set_points(values.shape[0]);
    else:
      self.set_points(points);

    if values.shape[0] != self._points.shape[0]:
      raise ValueError('number of values %d mismatch number of points %d' % (values.shape[0], self._points.shape[0]));


    #set parameter from values
    if self.with_projection and self.projection_inverse is not False:
      self._parameter[:,pdims] = self._projection_inverse.dot(values);
      #self._values[:,pdims] = values;
    else:
      #tck,u = splprep(values, u = self.points, t = self.knots, task = -1, k = self.degree, s = 0); # splprep crashes due to erros in fitpack
      #for d in range(self.ndim):
      #  self.parameter[d] = tck[1][d];
      #  self.values[d] = self.basis.dot(self.parameter[d]);
      for v,p in zip(vdims, pdims):
        tck = splrep(self.points, values[:,v], t = self.knots, task = -1, k = self.degree);
        self.parameter[:,p] = tck[1][:self.nparameter];

      # values will change
      self._values = None;
      #self._values = values; #fast but imprecise as values of spline due approximation might differ!