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)
python类splev()的实例源码
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 circ(self):
'''returns a rough estimate of the circumference'''
if self._circ is None:
new_pts = si.splev(np.linspace(0, 1, 1000), self.tck, ext=2)
self._circ = np.sum(np.sqrt(np.sum(np.diff(new_pts, axis=1) ** 2,
axis=0)))
return self._circ
def cntr(self):
'''returns a rough estimate of the circumference'''
if self._cntr is None:
new_pts = si.splev(np.linspace(0, 1, 1000), self.tck, ext=2)
self._cntr = np.mean(new_pts, 1)
return self._cntr
def ts():
return splev(np.linspace(0,1,s.npoints), tck);
def get_values(self, points = None, parameter = None, dimension = None, derivative = 0, extrapolation = 1):
"""Calculates the values of the curve along the sample points
Arguments:
points (array or None): the sample points for the values, if None use internal samples points
parameter (array or None): the bspline parameter, if None use internal parameter
dimensions (None, list or int): the dimension(s) for which to return values
derivative (int): the order of the derivative
extrapolation (int): 0=extrapolated value, 1=return 0, 2=raise a ValueError, 3=boundary value
Returns:
array: the values of the spline at the sample points
"""
if dimension is None:
dimension = range(self.ndim);
if isinstance(dimension,int):
dimension = [dimension];
points = self.get_points(points, error = 'sample points need to be specified for the calculation of the values of this curve!')
if parameter is None:
parameter = self._parameter[:,dimension];
if points is self._points and derivative == 0:
if parameter is self._parameter and self._values is not None: #cached version
return self._values[:,dimension];
if self.with_projection:
return self.projection.dot(parameter);
# full interpolation
tck = (self._knots_all, parameter.T, self.degree);
values = splev(points, tck, der = derivative, ext = extrapolation);
return np.vstack(values).T
def get_points_and_values(self, points = None, parameter = None, dimension = None, derivative = 0, extrapolation = 1):
"""Calculates the values of the curve along the sample points
Arguments:
points (array or None): the sample points for the values, if None use internal samples points
parameter (array or None): the bspline parameter, if None use internal parameter
dimensions (None, list or int): the dimension(s) for which to return values
derivative (int): the order of the derivative
extrapolation (int): 0=extrapolated value, 1=return 0, 2=raise a ValueError, 3=boundary value
Returns:
array: tha sample points
array: the values of the spline at the sample points
"""
if dimension is None:
dimension = range(self.ndim);
if isinstance(dimension,int):
dimension = [dimension];
points = self.get_points(points, error = 'sample points need to be specified for the calculation of the values of this curve!')
if parameter is None:
parameter = self._parameter[:,dimension];
if points is self._points and derivative == 0:
if parameter is self._parameter and self._values is not None: #cached version
return points, self._values[:,dimension];
if self.with_projection:
return points, self.projection.dot(parameter);
# full interpolation
tck = (self._knots_all, parameter.T, self.degree);
values = splev(points, tck, der = derivative, ext = extrapolation);
return points, np.vstack(values).T
def theta(self, points = None, with_xy = True, with_length = True, with_orientation = True, reference = 0.5):
"""Returns a Spline representing the derivative of the tangent angle along a 2d curve
Arguments:
points (int, array or None): sample points used to determine theta
with_lenth (bool): if True also return length of the curve
with_position (bool): if True also return absolute position
with_orientation (bool): if True also return absolute orientation of the curve
reference (float): reference point for absolute position and orientation
Returns:
Spline: spline of theta
Note:
To fully reconstruct the curve, the center point, length and orientation is needed.
"""
if self.ndim != 2:
raise RuntimeError('theta angle can only be computed for 2d curves');
points = self.get_points(points, error = 'cannot determine sample points needed for the calculation of theta');
#get the tangents and tangent angles
tgs = splev(points, self.tck(), der = 1);
tgs = np.vstack(tgs).T;
phi = np.arctan2(tgs[:,1],tgs[:,0]);
phi = np.mod(phi + np.pi, 2 * np.pi) - np.pi;
#phi = Spline(phi, points = self.points, knots = self.knots, degree = self.degree + 1);
tck = splrep(points, phi, s = 0.0, k = self.degree + 1);
phi = Spline(tck = tck);
orientation = phi(reference);
theta = phi.derivative();
rtr = [theta];
if with_xy:
rtr.append(self(reference));
if with_length:
rtr.append(self.length());
if with_orientation:
rtr.append(orientation);
return tuple(rtr);
def get_values(self, points = None, parameter = None, derivative = 0, extrapolation = 1):
"""Calculates the values of the spline along the sample points
Arguments:
points (array or None): the sample points for the values, if None use internal samples points
parameter (array or None): the bspline parameter, if None use internal parameter
derivative (int): the order of the derivative
extrapolation (int): 0=extrapolated value, 1=return 0, 2=raise a ValueError, 3=boundary value
Returns:
array: the values of the spline at the sample points
"""
points = self.get_points(points, error = 'sample points need to be specified for the calculation of the values of this spline!')
if parameter is None:
parameter = self._parameter;
if points is self._points and derivative == 0:
if parameter is self._parameter and self._values is not None: #cached version
return self._values;
if self.with_projection:
return self.projection.dot(parameter);
# full interpolation
pp = np.pad(parameter, (0,self.degree+1), 'constant');
tck = (self._knots_all, pp, self.degree);
return splev(points, tck, der = derivative, ext = extrapolation);
def get_values(self, points = None, parameter = None, derivative = 0, extrapolation = 1):
"""Calculates the values of the spline along the sample points
Arguments:
points (array or None): the sample points for the values, if None use internal samples points
parameter (array or None): the bspline parameter, if None use internal parameter
derivative (int): the order of the derivative
extrapolation (int): 0=extrapolated value, 1=return 0, 2=raise a ValueError, 3=boundary value
Returns:
array: the values of the spline at the sample points
"""
points = self.get_points(points, error = 'sample points need to be specified for the calculation of the values of this spline!')
if parameter is None:
parameter = self._parameter;
if points is self._points and derivative == 0:
if parameter is self._parameter and self._values is not None: #cached version
return self._values;
if self.with_projection:
return self.projection.dot(parameter);
# full interpolation
pp = np.pad(parameter, (0,self.degree+1), 'constant');
tck = (self._knots_all, pp, self.degree);
return splev(points, tck, der = derivative, ext = extrapolation);
def get_points_and_values(self, points = None, parameter = None, derivative = 0, extrapolation = 1, error = None):
"""Calculates the values of the spline along the sample points
Arguments:
points (array or None): the sample points for the values, if None use internal samples points
parameter (array or None): the bspline parameter, if None use internal parameter
derivative (int): the order of the derivative
extrapolation (int): 0=extrapolated value, 1=return 0, 2=raise a ValueError, 3=boundary value
Returns:
array: tha sample points
array: the values of the spline at the sample points
"""
points = self.get_points(points, error = error)
if parameter is None:
parameter = self._parameter;
if points is self._points and derivative == 0:
if parameter is self._parameter and self._values is not None: #cached version
return points, self._values;
if self.with_projection:
return points, self.projection.dot(parameter);
# full interpolation
pp = np.pad(parameter, (0,self.degree+1), 'constant');
tck = (self._knots_all, pp, self.degree);
return points, splev(points, tck, der = derivative, ext = extrapolation);
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 convex_hull_removal(pixel, wvl):
"""
Remove the convex-hull of the signal by hull quotient.
Parameters:
pixel: `list`
1D HSI data (p), a pixel.
wvl: `list`
Wavelength of each band (p x 1).
Results: `list`
Data with convex hull removed (p).
Reference:
Clark, R.N. and T.L. Roush (1984) Reflectance Spectroscopy: Quantitative
Analysis Techniques for Remote Sensing Applications, J. Geophys. Res., 89,
6329-6340.
"""
points = list(zip(wvl, pixel))
# close the polygone
poly = [(points[0][0],0)]+points+[(points[-1][0],0)]
hull = _jarvis.convex_hull(poly)
# the last two points are on the x axis, remove it
hull = hull[:-2]
x_hull = [u for u,v in hull]
y_hull = [v for u,v in hull]
tck = interpolate.splrep(x_hull, y_hull, k=1)
iy_hull = interpolate.splev(wvl, tck, der=0)
norm = []
for ysig, yhull in zip(pixel, iy_hull):
if yhull != 0:
norm.append(ysig/yhull)
else:
norm.append(1)
return norm, x_hull, y_hull
def smooth_plot(figure,X,Y,color1,color2,xlabel="",ylabel="",filled=False,n_points=400,smooth_factor=1.0,spline_order=3,linewidth=3,alpha=1.0,label=""):
"""
"""
X_smooth = np.linspace(X.min(),X.max(),n_points)
tck = splrep(X,Y,s=smooth_factor,k=spline_order)
Y_smooth = splev(X_smooth,tck,der=0)
font = fm.FontProperties(family = 'Trebuchet', weight ='light')
#font = fm.FontProperties(family = 'CenturyGothic',fname = '/Library/Fonts/Microsoft/Century Gothic', weight ='light')
figure.patch.set_facecolor('white')
axes = figure.add_subplot(111)
axes.plot(X,Y,linewidth=1,color=tuple(color2),alpha=0.2)
if filled:
axes.fill_between(X_smooth,Y_smooth,0,color=color2,alpha=0.1)
for i in xrange(100):
color = tuple(color1*(1.0-i/100.0) + color2*(i/100.0))
if i == 0:
axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha,label=label)
else:
axes.plot(X_smooth[(i*n_points/100):((i+1)*n_points)/100+1],Y_smooth[(i*n_points/100):((i+1)*n_points)/100+1],linewidth=linewidth,color=color,alpha=alpha)
axes.set_xlim(X.min(),X.max())
axes.set_xlabel(xlabel,fontproperties=font, size=10, style='italic')
axes.set_xticklabels(axes.get_xticks(),fontproperties=font, size=12)
if '%' in ylabel:
axes.set_ylim(0,np.minimum(2*Y.max(),100))
axes.set_ylabel(ylabel, fontproperties=font, size=10, style='italic')
axes.set_yticklabels(axes.get_yticks(),fontproperties=font, size=12)
def _interpolate(self, xa, xb, data):
from scipy import interpolate
f = interpolate.splrep(xa, data)
return interpolate.splev(xb, f)
def discretize_bspline(control, knots, count=None, scale=1.0):
'''
Given a B-Splines control points and knot vector, return
a sampled version of the curve.
Arguments
----------
control: (o,d) list of control points of the b- spline.
knots: (j) list of knots
count: int, number of sections to discretize the spline in to.
If not specified, RES_LENGTH will be used to inform this.
Returns
----------
discrete: (count,d) list of points, a polyline of the B-spline.
'''
# evaluate the b-spline using scipy/fitpack
from scipy.interpolate import splev
# (n, d) control points where d is the dimension of vertices
control = np.array(control)
degree = len(knots) - len(control) - 1
if count is None:
norm = np.linalg.norm(np.diff(control, axis=0), axis=1).sum()
count = int(np.clip(norm / (res.seg_frac*scale),
res.min_sections*len(control),
res.max_sections*len(control)))
ipl = np.linspace(knots[0], knots[-1], count)
discrete = splev(ipl, [knots, control.T, degree])
discrete = np.column_stack(discrete)
return discrete
def formPEC(R, Rmin, Rmax, E, De, limb):
evcm = const.e/(const.c*const.h*100) # converts cm-1 to eV
# combine Rmin with Rmax to form PEC
Re = (Rmin[0] + Rmax[0])/2
print(u"RKR: Re = {:g}".format(Re))
RTP = np.append(Rmin[::-1], Rmax, 0) # radial positions of turning-points
PTP = np.append(E[::-1], E, 0) # potential energy at turning point
# interpolate
psp = splrep(RTP, PTP, s=0)
# Interpolate RKR curve to this grid
PEC = splev(R, psp, der=0)
# extrapolate using analytical function
inner_limb_Morse(R, PEC, RTP, PTP, Re, De)
if limb=='L':
outer_limb_LeRoy(R, PEC, RTP, PTP, De)
else:
outer_limb_Morse(R, PEC, RTP, PTP, De)
PTP /= evcm
PEC /= evcm # convert to eV
return PEC, RTP, PTP
# analytical functions
def p2eRatio(self):
"""
Calculates the proton density to electron density ratio using Eq. 7 of [1]_.
Notes
------
Uses the abundance and ionization equilibrium.
References
----------
.. [1] Young, P. R. et al., 2003, ApJS, `144, 135 <http://adsabs.harvard.edu/abs/2003ApJS..144..135Y>`_
"""
if hasattr(self, 'Temperature'):
temperature = self.Temperature
else:
temperature = self.IoneqAll['ioneqTemperature']
if not hasattr(self, 'AbundanceName'):
AbundanceName = self.Defaults['abundfile']
else:
AbundanceName = self.AbundanceName
tmp_abundance = io.abundanceRead(abundancename=AbundanceName)
abundance = tmp_abundance['abundance'][tmp_abundance['abundance']>0]
denominator = np.zeros(len(self.IoneqAll['ioneqTemperature']))
for i in range(len(abundance)):
for z in range(1,i+2):
denominator += z*self.IoneqAll['ioneqAll'][i,z,:]*abundance[i]
p2eratio = abundance[0]*self.IoneqAll['ioneqAll'][0,1,:]/denominator
nots = interpolate.splrep(np.log10(self.IoneqAll['ioneqTemperature']),p2eratio,s=0)
self.ProtonDensityRatio = interpolate.splev(np.log10(temperature),nots,der=0,ext=1)
