python类RectBivariateSpline()的实例源码

def __init__(self, scale_list, values, extrapolate=False):
        # error checking
        if len(values.shape) != 2:
            raise ValueError('This class only works for 2D data.')
        elif len(scale_list) != 2:
            raise ValueError('input and output dimension mismatch.')

        nx, ny = values.shape
        offset, scale = scale_list[0]
        x = np.linspace(offset, (nx - 1) * scale + offset, nx)  # type: np.multiarray.ndarray
        offset, scale = scale_list[1]
        y = np.linspace(offset, (ny - 1) * scale + offset, ny)  # type: np.multiarray.ndarray

        self._min = x[0], y[0]
        self._max = x[-1], y[-1]

        DiffFunction.__init__(self, [(x[0], x[-1]), (y[0], y[-1])], delta_list=None)

        self.fun = interp.RectBivariateSpline(x, y, values)
        self._extrapolate = extrapolate

def get_interpolated_pixel_color_rbspline(pts, s_im, size):
    """given pts in floats, linear interpolate pixel values nearby to get a good colour"""
    pts = clamp(pts, size)

    s_im = np.atleast_3d(s_im)
    ys,xs = size
    ycoords, xcoords = np.arange(ys), np.arange(xs)
    out = np.empty(pts.shape[1:] + (s_im.shape[-1],),dtype=s_im.dtype)

    pts_vec = pts.reshape((2,-1))
    out_vec = out.reshape((-1,s_im.shape[-1])) #flatten for easier vectorization
    for i in range(s_im.shape[-1]): #loop over color channels
        rbspline = RectBivariateSpline(ycoords, xcoords, s_im[...,i])
        out_vec[:,i] = rbspline.ev(pts_vec[0],pts_vec[1])
    return out


### Functions generating SL(2,C) matrices ###
# Do not need to be vectorized #

def resample_2d(array, sample_pts, query_pts):
    ''' Resamples 2D array to be sampled along queried points.

    Args:
        array (numpy.ndarray): 2D array.

        sample_pts (tuple): pair of numpy.ndarray objects that contain the x and y sample locations,
            each array should be 1D.

        query_pts (tuple): points to interpolate onto, also 1D for each array.

    Returns:
        numpy.ndarray.  array resampled onto query_pts via bivariate spline.

    '''
    xq, yq = np.meshgrid(*query_pts)
    interpf = interpolate.RectBivariateSpline(*sample_pts, array)
    return interpf.ev(yq, xq)

def __get_2d_function(x, y, z, kind):
        """
        Train 2-D interpolator
        :param x: x-coordinates of data points
        :param y: y-coordinates of data points
        :param z: z-coordinates of data points
        :param kind: 'linear' or 'spline' interpolation type
        :return Interpolator callable object
        """
        try:
            from scipy.interpolate import RectBivariateSpline
        except ImportError as e:
            logger.error("Please install scipy on your platform to be able to use spline-based interpolation")
            raise e
        if len(x) == 2:  # fall-back to linear interpolation
            kx = 1
        elif len(x) == 3:  # fall-back to 2nd degree spline
            kx = 2
        else:
            kx = 3
        if len(y) == 2:  # fall-back to linear interpolation
            ky = 1
        elif len(y) == 3:  # fall-back to 2nd degree spline
            ky = 2
        else:
            ky = 3
        if kind == 'linear':
            kx, ky = 1, 1
        x_array, y_array, z_array = x, y, z
        if not isinstance(x, np.ndarray):
            x_array = np.asarray(x)
        if not isinstance(y, np.ndarray):
            y_array = np.asarray(y)
        if not isinstance(z, np.ndarray):
            z_array = np.asarray(z)
        result = RectBivariateSpline(x_array, y_array, z_array, kx=kx, ky=ky, s=0)
        return result

def computeInterpolatedICImg(self):

        """Computes interpolation of initial condition image.

        Interpolation is done as in :py:func:`pyfrp.modules.pyfrp_sim_module.applyInterpolatedICs`.

        Returns:
            tuple: Tuple containing:

                * xInt (numpy.ndarray): Meshgrid x-coordinates.
                * yInt (numpy.ndarray): Meshgrid y-coordinates.
                * f (numpy.ndarray): Interpolated image.

        """

        #Get image resolution and center of geometry
        res=self.ICimg.shape[0]
        center=self.embryo.geometry.getCenter()

        #Define x/y coordinates of interpolation
        if 'quad' in self.embryo.analysis.process.keys():
            #Shift everything by center to fit with the mesh
            xInt = np.arange(center[0]+1, center[0]+res+1, 1)
            yInt = np.arange(center[1]+1, center[1]+res+1, 1)       
        else:
            xInt = np.arange(1, res+1, 1)
            yInt = np.arange(1, res+1, 1)

        #Generate interpolation function
        f=interp.RectBivariateSpline(xInt, yInt, self.ICimg, bbox=[None, None, None, None], kx=3, ky=3, s=0)

        return xInt, yInt, f

def __init__(self):
        aftmap_filename = os.path.join(DATA_DIR, 's1_aft_rz_02Mar2017.json')
        with open(aftmap_filename) as data_file:
            data = json.load(data_file)
        r_pts = np.array(data['r_pts'])
        z_pts = np.array(data['z_pts'])
        aft_vals = np.array(data['map']).reshape(len(r_pts), len(z_pts))
        self.aft_map = RectBivariateSpline(r_pts, z_pts, aft_vals)

def fref(rad,cosb,height,t_setx,r_setx,s_setx,sang):
    ref=reflectance(rad,cosb,t_setx,height,r_setx,s_setx,sang)
    return interpolate.RectBivariateSpline(t_setx,r_setx,ref)

def omega_spline(self, parameters, s=0):
        omega = self.omega(parameters)
        return spinterp.RectBivariateSpline(self.kx[:,0], self.ky[0,:], omega, s=s)

def rotate_ground(original, theta, horizon=60, half_height=360 / 2, focal=1.0):
        height, width, channel = original.shape
        # the target grids
        yp = range(height - horizon, height)
        xp = range(0, width)

        # from pixel to coordinates
        y0 = (np.array(yp) - half_height) * 1.0 / half_height
        x0 = (np.array(xp) - width / 2) / (width / 2.0)

        # form the mesh
        mesh = MyDataset.generate_meshlist(x0, y0)
        # compute the source coordinates
        st = math.sin(theta)
        ct = math.cos(theta)
        deno = ct * focal + st * mesh[:, 0]
        out = np.array([(-st * focal + ct * mesh[:, 0]) / deno, mesh[:, 1] / deno])

        # interpolate
        vout = []
        for i in range(3):
            f = interpolate.RectBivariateSpline(y0, x0, original[- horizon:, :, i])
            values = f(out[1, :], out[0, :], grid=False)
            vout.append(values)

        lower = np.reshape(vout, (3, width, horizon)).transpose((2, 1, 0)).astype("uint8")

        # compute the upper part
        out = np.reshape(out[0, :], (width, horizon))
        out = out[:, 0]
        f = interpolate.interp1d(x0, original[:-horizon, :, :], axis=1,
                                 fill_value=(original[:-horizon, 0, :], original[:-horizon, -1, :]),
                                 bounds_error=False)
        ans = f(out)
        ans = ans.astype("uint8")

        return np.concatenate((ans, lower), axis=0)

def __call__(self, *args, **kwargs):
        # TODO: Should be more defensive about this
        if len(args) == 1:
            # We have a dataframe
            df = args[0]
        else:
            # We have parameters for a single invocation
            df = pd.DataFrame(kwargs)

        if self.key_columns:
            sub_tables = df.groupby(self.key_columns)
        else:
            sub_tables = [(None, df)]

        result = pd.DataFrame(index=df.index)
        for key, sub_table in sub_tables:
            if sub_table.empty:
                continue
            funcs = self.interpolations[key]
            parameters = tuple(sub_table[k] for k in self.parameter_columns)
            for value_column, func in funcs.items():
                out = func(*parameters)
                # This reshape is necessary because RectBivariateSpline and InterpolatedUnivariateSpline return results
                # in slightly different shapes and we need them to be consistent
                if out.shape:
                    result.loc[sub_table.index, value_column] = out.reshape((out.shape[0],))
                else:
                    result.loc[sub_table.index, value_column] = out

        if self.func:
            return self.func(result)

        if len(result.columns) == 1:
            return result[result.columns[0]]

        return result

def interpolate_2d_features(features):
    out_size = feature_size
    x = np.arange(features.shape[0])
    y = np.arange(features.shape[1])
    z = features
    xx = np.linspace(x.min(), x.max(), out_size)
    yy = np.linspace(y.min(), y.max(), out_size)
    new_kernel = interpolate.RectBivariateSpline(x, y, z, kx=1, ky=1)
    kernel_out = new_kernel(xx, yy)
    return kernel_out


# Interpolation 2d of each channel, so we obtain 3d interpolated feature maps

def py_sampleImage(reference_image, dxy, udat, vdat, dRA=0., dDec=0., PA=0.):
    """
    Python implementation of sampleImage.

    """
    nxy = reference_image.shape[0]

    dRA *= 2.*np.pi
    dDec *= 2.*np.pi

    du = 1. / (nxy*dxy)

    # Real to Complex transform
    fft_r2c_shifted = np.fft.fftshift(
                        np.fft.rfft2(
                            np.fft.fftshift(reference_image)), axes=0)

    # apply rotation
    cos_PA = np.cos(PA)
    sin_PA = np.sin(PA)

    urot = udat * cos_PA - vdat * sin_PA
    vrot = udat * sin_PA + vdat * cos_PA

    dRArot = dRA * cos_PA - dDec * sin_PA
    dDecrot = dRA * sin_PA + dDec * cos_PA

    # interpolation indices
    uroti = np.abs(urot)/du
    vroti = nxy/2. + vrot/du
    uneg = urot < 0.
    vroti[uneg] = nxy/2 - vrot[uneg]/du

    # coordinates of FT
    u_axis = np.linspace(0., nxy // 2, nxy // 2 + 1)
    v_axis = np.linspace(0., nxy - 1, nxy)

    # We use RectBivariateSpline to do only linear interpolation, which is faster
    # than interp2d for our case of a regular grid.
    # RectBivariateSpline does not work for complex input, so we need to run it twice.
    f_re = RectBivariateSpline(v_axis, u_axis, fft_r2c_shifted.real, kx=1, ky=1, s=0)
    ReInt = f_re.ev(vroti, uroti)
    f_im = RectBivariateSpline(v_axis, u_axis, fft_r2c_shifted.imag, kx=1, ky=1, s=0)
    ImInt = f_im.ev(vroti, uroti)

    # correct for Real to Complex frequency mapping
    uneg = urot < 0.
    ImInt[uneg] *= -1.

    # apply the phase change
    theta = urot*dRArot + vrot*dDecrot
    vis = (ReInt + 1j*ImInt) * (np.cos(theta) + 1j*np.sin(theta))

    return vis

def ionex_stec_map(ionex_fname,
                   augmented_arc_map):
    """
    """
    logger.info('computing interpolated and mapped STEC from {}'.format(ionex_fname))
    # compute temporally interpolated IONEX VTEC maps
    dt_list = get_dt_list(augmented_arc_map)
    (grid_lon, grid_lat, vtec,
     _, sat_biases, _) = interpolate2D_temporal(ionex_fname,
                                                dt_list)
    # compute spline interpolators
    bbox = [-180, 180, -90, 90]
    # check that latitude grid is in decreasing order
    assert grid_lat[1] < grid_lat[0]
    interp_map = {dt: RectBivariateSpline(grid_lon,
                                          grid_lat[::-1],
                                          vtec[:, ::-1, i],
                                          bbox=bbox) for i, dt in enumerate(dt_list)}
    # compute interpolated stec
    stec_map = defaultdict(list)
    for key, arc_list in augmented_arc_map.iteritems():
        for arc in arc_list:
            ionex_stec = []
            for (dt_i, ipp_lat_i, ipp_lon_i, el_map_i) in zip(arc.dt,
                                                              arc.ipp_lat,
                                                              arc.ipp_lon,
                                                              arc.el_map):
                i, _ = find_le(dt_list, dt_i)
                interpolator = interp_map[dt_list[i]]
                vtec_i = float(interpolator.ev(ipp_lon_i,
                                               ipp_lat_i))
                ionex_stec.append(vtec_i * el_map_i)
            stec_map[key].append(ionex_stec)
    return stec_map, sat_biases

def scale_field(self, f):
        """
        Scaling of the field according to calculated meshgrid.

        For the interpolation rectangular bivariate Splines are used.
        These are implemented from the scipy function `RectBivariateSpline <https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>`_.

        :param f: field to be interpolated
        :returns: field after interpolation
        """
        rbs = RBS(self._tx, self._ty, f)
        return rbs.ev(self._out_x, self._out_y)

def correlate(image, pattern,
              include_direct_beam=False,
              sim_threshold=1e-5,
              interpolate=False,
              **kwargs):
    """The correlation between a diffraction pattern and a simulation.
    Calculated using
        .. math::
            \frac{\sum_{j=1}^m P(x_j, y_j) T(x_j, y_j)}{\sqrt{\sum_{j=1}^m P^2(x_j, y_j)} \sqrt{\sum_{j=1}^m T^2(x_j, y_j)}}
    Parameters
    ----------
    image : :class:`np.ndarray`
        A single electron diffraction signal. Should be appropriately scaled
        and centered.
    pattern : :class:`DiffractionSimulation`
        The pattern to compare to.
    sim_threshold : float
        The threshold simulation intensity to consider for correlation
    interpolate : bool
        If True, perform sub-pixel interpolation of the image.
    **kwargs
        Arguments to pass to scipy.interpolate.RectBivariateSpline
    Returns
    -------
    float
        The correlation coefficient.
    References
    ----------
    E. F. Rauch and L. Dupuy, “Rapid Diffraction Patterns identification through
        template matching,” vol. 50, no. 1, pp. 87–99, 2005.
    """
    shape = image.shape
    half_shape = tuple(i // 2 for i in shape)

    pixel_coordinates = pattern.calibrated_coordinates.astype(int)[
        :, :2] + half_shape
    in_bounds = np.product((pixel_coordinates > 0) *
                           (pixel_coordinates < shape[0]), axis=1).astype(bool)
    pattern_intensities = pattern.intensities
    large_intensities = pattern_intensities > sim_threshold
    mask = np.logical_and(in_bounds, large_intensities)

    if interpolate:
        x = np.arange(shape[0], dtype='float') - half_shape[0]
        y = np.arange(shape[1], dtype='float') - half_shape[1]
        for ar, i in zip([x, y], shape):
            if not i % 2:
                ar += 0.5
        x = x * pattern.calibration[0]
        y = y * pattern.calibration[1]
        ip = RectBivariateSpline(x, y, image.T, **kwargs)
        image_intensities = ip.ev(pattern.coordinates[:, 0][mask],
                                  pattern.coordinates[:, 1][mask])
    else:
        image_intensities = image.T[pixel_coordinates[:, 0][mask], pixel_coordinates[:, 1][mask]]
    pattern_intensities = pattern_intensities[mask]
    return np.nan_to_num(_correlate(image_intensities, pattern_intensities))

def interp_spline(I, xn, yn, compute_derivs=True):
    """ Perform spline interpolation of I.

    I is evaluated at xn, yn.

    Returns
    -------
    I_warped : array_like
        Warped image
    dI_warped_dx : array_like
        Derivative of warped image in x direction
    dI_warped_dy : array_like
        Derivative of warped image in y direction
    """

    h,w = I.shape[:2]
    xar = np.arange(w)
    yar = np.arange(h)

    if I.ndim > 2:
        I_warped = np.zeros_like(I)
        dI_warped_dx = np.zeros_like(I)
        dI_warped_dy = np.zeros_like(I)

        for d in range(I.shape[2]):
            result = interp_spline(I[:,:,d], xn, yn, compute_derivs=compute_derivs)
            if compute_derivs:
                I_warped[:,:,d] = result[0]
                dI_warped_dx[:,:,d] = result[1]
                dI_warped_dy[:,:,d] = result[2]
            else:
                I_warped[:,:,d] = result
    else:
        S = interpolate.RectBivariateSpline(yar,xar,I,kx=2,ky=2)
        I_warped = np.clip(S.ev(yn,xn),0,1.0)
        #I_warped = S.ev(yn,xn)
        if compute_derivs:
            dI_warped_dx = S.ev(yn,xn,dy=1) # Note that y and x are swapped!
            dI_warped_dy = S.ev(yn,xn,dx=1)

    if compute_derivs:
        return I_warped, dI_warped_dx, dI_warped_dy
    else:
        return I_warped

def interpolate_curvature(curvatures, times, kind='cubic'):
    """ Interpolate the curvatures.

    A curvature is a parametrisation `s`, d(angle)/ ds and a parameter between [0,1].
    Return a surface f(s,t).
    """
    import numpy as np

    from scipy import interpolate
    from scipy.ndimage import measurements

    if len(curvatures) == 3 and not isinstance(curvatures[2], tuple):
        curvatures= [curvatures]

    curv_abs = [c[0] for c in curvatures]
    curves = [c[1] for c in curvatures]
    params = [c[2] for c in curvatures]

    n = len(params)
    if n==1:
        curv_abs*=2
        curves*=2
        curves[0]=np.zeros(len(curves[0]))
        params = [0.,1.]
    elif False:
        if params[0] != 0.:
            params.insert(0,0.)
            curves.insert(0,curves[0])
            curv_abs.insert(0, curv_abs[0])
        if params[-1] != 0:
            params.append(1.)
            curves.append(curves[-1])
            curv_abs.append(curv_abs[-1])

    # compute a common parametrisation s
    # We conserve the last parametrisation because the result is very sensitive
    if True:
        min_s = min(np.diff(s).min() for s in curv_abs)
        s_new = np.unique(np.array(curv_abs).flatten())
        ds = np.diff(s_new)
        k = np.cumsum(ds >= min_s)
        labels = range(k.min(), k.max()+1)
        s = np.zeros(len(labels)+1)
        s[1:] = measurements.mean(s_new[1:], k, labels)
        s[-1]=1.
        s = s_new
    else:
        s = np.array(curv_abs[-1])

    # renormalise all the curves
    curves = [np.interp(s[1:-1], old_s[1:-1], old_crv) for old_s, old_crv in zip(curv_abs,curves)]

    # interpolate correctly the curvatures
    x = s[1:-1]
    y = np.array(params)
    z = np.array(curves)

    #f = interpolate.interp2d(y, x, z, kind=kind)
    f = interpolate.RectBivariateSpline(x, y, z.T, kx=1, ky=1)
    return s, f(x,times)

def __init__(self, data, categorical_parameters, continuous_parameters, func=None, order=1):
        self.key_columns = categorical_parameters
        self.parameter_columns = continuous_parameters
        self.func = func

        if len(self.parameter_columns) not in [1, 2]:
            raise ValueError("Only interpolation over 1 or 2 variables is supported")
        if len(self.parameter_columns) == 1 and order == 0:
            raise ValueError("Order 0 only supported for 2d interpolation")

        # These are the columns which the interpolation function will approximate
        value_columns = sorted(data.columns.difference(set(self.key_columns)|set(self.parameter_columns)))

        if self.key_columns:
            # Since there are key_columns we need to group the table by those
            # columns to get the sub-tables to fit
            sub_tables = data.groupby(self.key_columns)
        else:
            # There are no key columns so we will fit the whole table
            sub_tables = {None: data}.items()

        self.interpolations = {}

        for key, base_table in sub_tables:
            if base_table.empty:
                continue
            # For each permutation of the key columns build interpolations
            self.interpolations[key] = {}
            for value_column in value_columns:
                # For each value in the table build an interpolation function
                if len(self.parameter_columns) == 2:
                    # 2 variable interpolation
                    if order == 0:
                        x = base_table[list(self.parameter_columns)]
                        y = base_table[value_column]
                        func = interpolate.NearestNDInterpolator(x=x.values, y=y.values)
                    else:
                        index, column = self.parameter_columns
                        table = base_table.pivot(index=index, columns=column, values=value_column)
                        x = table.index.values
                        y = table.columns.values
                        z = table.values
                        func = interpolate.RectBivariateSpline(x=x, y=y, z=z, ky=order, kx=order).ev
                else:
                    # 1 variable interpolation
                    base_table = base_table.sort_values(by=self.parameter_columns[0])
                    x = base_table[self.parameter_columns[0]]
                    y = base_table[value_column]
                    func = interpolate.InterpolatedUnivariateSpline(x, y, k=order)
                self.interpolations[key][value_column] = func