Python scipy.ndimage.map_coordinates() Examples

def bicubic_interpolation(LR, HR_shape):

    LR_padded = np.zeros((LR.shape[0]+1,LR.shape[1]+1,LR.shape[2]+1))
    LR_padded[:-1,:-1,:-1] = LR[:,:,:]
    LR_padded[-1,:-1,:-1] = LR[0,:,:]
    LR_padded[:-1,-1,:-1] = LR[:,0,:]
    LR_padded[:-1,:-1,-1] = LR[:,:,0]
    LR_padded[:-1,-1,-1] = LR[:,0,0]
    LR_padded[-1,:-1,-1] = LR[0,:,0]
    LR_padded[-1,-1,:-1] = LR[0,0,:]
    LR_padded[-1,-1,-1] = LR[0,0,0]
    
    x_HR = np.linspace(0, LR.shape[0], num=HR_shape[0]+1)[:-1]
    y_HR = np.linspace(0, LR.shape[1], num=HR_shape[1]+1)[:-1]
    z_HR = np.linspace(0, LR.shape[2], num=HR_shape[2]+1)[:-1]
    
    xx, yy, zz = np.meshgrid(x_HR, y_HR, z_HR, indexing='ij')
    
    xx = xx.reshape((-1))
    yy = yy.reshape((-1))
    zz = zz.reshape((-1))
    out_BC = ndimage.map_coordinates(LR_padded, [xx, yy, zz], order=3, mode='wrap').reshape(HR_shape)
    
    return out_BC 

def test_map_coordinates03(self):
        data = numpy.array([[4, 1, 3, 2],
                               [7, 6, 8, 5],
                               [3, 5, 3, 6]], order='F')
        idx = numpy.indices(data.shape) - 1
        out = ndimage.map_coordinates(data, idx)
        assert_array_almost_equal(out, [[0, 0, 0, 0],
                                   [0, 4, 1, 3],
                                   [0, 7, 6, 8]])
        assert_array_almost_equal(out, ndimage.shift(data, (1, 1)))
        idx = numpy.indices(data[::2].shape) - 1
        out = ndimage.map_coordinates(data[::2], idx)
        assert_array_almost_equal(out, [[0, 0, 0, 0],
                                   [0, 4, 1, 3]])
        assert_array_almost_equal(out, ndimage.shift(data[::2], (1, 1)))
        idx = numpy.indices(data[:,::2].shape) - 1
        out = ndimage.map_coordinates(data[:,::2], idx)
        assert_array_almost_equal(out, [[0, 0], [0, 4], [0, 7]])
        assert_array_almost_equal(out, ndimage.shift(data[:,::2], (1, 1))) 

def query(self, coords, order=1):
        """
        Returns the map value at the specified location(s) on the sky.

        Args:
            coords (`astropy.coordinates.SkyCoord`): The coordinates to query.
            order (Optional[int]): Interpolation order to use. Defaults to `1`,
                for linear interpolation.

        Returns:
            A float array containing the map value at every input coordinate.
            The shape of the output will be the same as the shape of the
            coordinates stored by `coords`.
        """
        out = np.full(len(coords.l.deg), np.nan, dtype='f4')

        for pole in self.poles:
            m = (coords.b.deg >= 0) if pole == 'ngp' else (coords.b.deg < 0)

            if np.any(m):
                data, w = self._data[pole]
                x, y = w.wcs_world2pix(coords.l.deg[m], coords.b.deg[m], 0)
                out[m] = map_coordinates(data, [y, x], order=order, mode='nearest')

        return out 

def __call__(self, img, steps=2):
        n, m = np.shape(img)
        def_imgs = []


        xn, yn = np.meshgrid(np.arange(m), np.arange(n))
        xn_mapped, yn_mapped = xn, yn

        for i in range(steps):
            if i == 0:
                Ic = img
            else:
                if self.multiplicative:
                    xn_mapped, yn_mapped = self.coodinate_mapper(xn_mapped, yn_mapped)
                else:
                    xn_mapped, yn_mapped = self.coodinate_mapper(xn, yn,frame=i)

                Ic = nd.map_coordinates(img, np.array([yn_mapped, xn_mapped]), order=self.order, cval=0)

            def_imgs.append(Ic)

        return def_imgs 

def resample_inplace(self, nbins):
        """
        Change the shape of the histogram using an n-dimensional
        interpolation function (via `ndimage.map_coordinates`).

        Parameters
        ----------
        nbins: tuple of int
            a tuple of the new number of bins to resample_inplace
            this histogram over (e.g. if the original histogram was
            (100,100), setting bins to (200,200) would provide double
            the resolution, with the interviening bins interpolated.
        """

        oldbins = self._nbins
        # iold = np.indices(oldbins)
        inew = np.indices(nbins)

        coords = np.array(
            [inew[X] * (oldbins[X]) / float(nbins[X]) for X in range(len(nbins))]
        )

        self._nbins = nbins
        self.data = ndimage.map_coordinates(self.data, coords)
        self._bin_lower_edges = None  # need to be recalculated 

def test_map_coordinates03(self):
        data = numpy.array([[4, 1, 3, 2],
                            [7, 6, 8, 5],
                            [3, 5, 3, 6]], order='F')
        idx = numpy.indices(data.shape) - 1
        out = ndimage.map_coordinates(data, idx)
        assert_array_almost_equal(out, [[0, 0, 0, 0],
                                        [0, 4, 1, 3],
                                        [0, 7, 6, 8]])
        assert_array_almost_equal(out, ndimage.shift(data, (1, 1)))
        idx = numpy.indices(data[::2].shape) - 1
        out = ndimage.map_coordinates(data[::2], idx)
        assert_array_almost_equal(out, [[0, 0, 0, 0],
                                        [0, 4, 1, 3]])
        assert_array_almost_equal(out, ndimage.shift(data[::2], (1, 1)))
        idx = numpy.indices(data[:, ::2].shape) - 1
        out = ndimage.map_coordinates(data[:, ::2], idx)
        assert_array_almost_equal(out, [[0, 0], [0, 4], [0, 7]])
        assert_array_almost_equal(out, ndimage.shift(data[:, ::2], (1, 1))) 

def test_uint64_max():
    # Test interpolation respects uint64 max.  Reported to fail at least on
    # win32 (due to the 32 bit visual C compiler using signed int64 when
    # converting between uint64 to double) and Debian on s390x.
    # Interpolation is always done in double precision floating point, so we
    # use the largest uint64 value for which int(float(big)) still fits in
    # a uint64.
    big = 2**64-1025
    arr = np.array([big, big, big], dtype=np.uint64)
    # Tests geometric transform (map_coordinates, affine_transform)
    inds = np.indices(arr.shape) - 0.1
    x = ndimage.map_coordinates(arr, inds)
    assert_equal(x[1], int(float(big)))
    assert_equal(x[2], int(float(big)))
    # Tests zoom / shift
    x = ndimage.shift(arr, 0.1)
    assert_equal(x[1], int(float(big)))
    assert_equal(x[2], int(float(big))) 

def interpolate_slice(f3d, rot, pfac=2, size=None):
    nhalf = f3d.shape[0] / 2
    if size is None:
        phalf = nhalf
    else:
        phalf = size / 2
    qot = rot * pfac  # Scaling!
    px, py, pz = np.meshgrid(np.arange(-phalf, phalf), np.arange(-phalf, phalf), 0)
    pr = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
    pcoords = np.vstack([px.reshape(-1), py.reshape(-1), pz.reshape(-1)])
    mcoords = qot.T.dot(pcoords)
    mcoords = mcoords[:, pr.reshape(-1) < nhalf]
    pvals = map_coordinates(np.real(f3d), mcoords, order=1, mode="wrap") + \
             1j * map_coordinates(np.imag(f3d), mcoords, order=1, mode="wrap")
    pslice = np.zeros(pr.shape, dtype=np.complex)
    pslice[pr < nhalf] = pvals
    return pslice 

def _profile(self) -> np.ndarray:
        """The actual profile array; private attr that is passed to MultiProfile."""
        return ndimage.map_coordinates(self.image_array, [self.y_locations, self.x_locations], order=0) 

def advection_correction(R, T=5, t=1):
    """
    R = np.array([qpe_previous, qpe_current])
    T = time between two observations (5 min)
    t = interpolation timestep (1 min)
    """

    # Evaluate advection
    oflow_method = motion.get_method("LK")
    fd_kwargs = {"buffer_mask": 10}  # avoid edge effects
    V = oflow_method(np.log(R), fd_kwargs=fd_kwargs)

    # Perform temporal interpolation
    Rd = np.zeros((R[0].shape))
    x, y = np.meshgrid(
        np.arange(R[0].shape[1], dtype=float), np.arange(R[0].shape[0], dtype=float),
    )
    for i in range(t, T + t, t):

        pos1 = (y - i / T * V[1], x - i / T * V[0])
        R1 = map_coordinates(R[0], pos1, order=1)

        pos2 = (y + (T - i) / T * V[1], x + (T - i) / T * V[0])
        R2 = map_coordinates(R[1], pos2, order=1)

        Rd += (T - i) * R1 + i * R2

    return t / T ** 2 * Rd


###############################################################################
# Finally, we apply the advection correction to the whole sequence of radar
# images and produce the rainfall accumulation map. 

def fuv2img(fuv, coorW=1024, floorW=1024, floorH=512):
    '''
    Project 1d signal in uv space to 2d floor plane image
    '''
    floor_plane_x, floor_plane_y = np.meshgrid(range(floorW), range(floorH))
    floor_plane_x, floor_plane_y = -(floor_plane_y - floorH / 2), floor_plane_x - floorW / 2
    floor_plane_coridx = (np.arctan2(floor_plane_y, floor_plane_x) / (2 * PI) + 0.5) * coorW - 0.5
    floor_plane = map_coordinates(fuv, floor_plane_coridx.reshape(1, -1), order=1, mode='wrap')
    floor_plane = floor_plane.reshape(floorH, floorW)
    return floor_plane 

def pano_stretch(img, mask, corners, kx, ky, order=1):
    '''
    img:     [H, W, C]
    corners: [N, 2] in image coordinate (x, y) format
    kx:      Stretching along front-back direction
    ky:      Stretching along left-right direction
    order:   Interpolation order. 0 for nearest-neighbor. 1 for bilinear.
    '''

    # Process image
    sin_u, cos_u, tan_v = uv_tri(img.shape[1], img.shape[0])
    u0 = np.arctan2(sin_u * kx / ky, cos_u)
    v0 = np.arctan(tan_v * np.sin(u0) / sin_u * ky)

    refx = (u0 / (2 * np.pi) + 0.5) * img.shape[1] - 0.5
    refy = (v0 / np.pi + 0.5) * img.shape[0] - 0.5

    # [TODO]: using opencv remap could probably speedup the process a little
    stretched_img = np.stack([
        map_coordinates(img[..., i], [refy, refx], order=order, mode='wrap')
        for i in range(img.shape[-1])
    ], axis=-1)
    stretched_mask = np.stack([
        map_coordinates(mask[..., i], [refy, refx], order=order, mode='wrap')
        for i in range(mask.shape[-1])
    ], axis=-1)
    #stretched_label = np.stack([
    #    map_coordinates(label[..., i], [refy, refx], order=order, mode='wrap')
    #    for i in range(label.shape[-1])
    #], axis=-1)

    # Process corners
    corners_u0 = coorx2u(corners[:, 0], img.shape[1])
    corners_v0 = coory2v(corners[:, 1], img.shape[0])
    corners_u = np.arctan2(np.sin(corners_u0) * ky / kx, np.cos(corners_u0))
    corners_v = np.arctan(np.tan(corners_v0) * np.sin(corners_u) / np.sin(corners_u0) / ky)
    cornersX = u2coorx(corners_u, img.shape[1])
    cornersY = v2coory(corners_v, img.shape[0])
    stretched_corners = np.stack([cornersX, cornersY], axis=-1)

    return stretched_img, stretched_mask, stretched_corners 

def interpolate(p, axis_values, pixelgrid, order=1, mode='constant', cval=0.0):
    """
    Interpolates in a grid prepared by create_pixeltypegrid().
    
    p is an array of parameter arrays
    
    @param p: Npar x Ninterpolate array
    @type p: array
    @return: Ndata x Ninterpolate array
    @rtype: array
    """
    # convert requested parameter combination into a coordinate
    p_ = np.array([np.searchsorted(av_,val) for av_, val in zip(axis_values,p)])
    lowervals_stepsize = np.array([[av_[p__-1], av_[p__]-av_[p__-1]] \
                         for av_, p__ in zip(axis_values,p_)])
    p_coord = (p-lowervals_stepsize[:,0])/lowervals_stepsize[:,1] + p_-1
    
    # interpolate
    if order > 1:
        prefilter = False
    else:
        prefilter = False

    return [ndimage.map_coordinates(pixelgrid[...,i],p_coord, order=order,
                                    prefilter=prefilter, mode=mode, cval=cval) \
                for i in range(np.shape(pixelgrid)[-1])] 

def fit_avg_z(cor_id, cor_id_xy, cor_img):
    score = map_coordinates(cor_img, [cor_id[:, 1], cor_id[:, 0]])
    c = np.sqrt((cor_id_xy ** 2).sum(1))
    cor_id_v = ((cor_id[:, 1] + 0.5) / cor_img.shape[0] - 0.5) * np.pi
    z = c * np.tan(cor_id_v)
    fit_z = (z * score).sum() / score.sum()
    return fit_z 

def permute_image(sy_image, mask, multiplier=1):
    shifter = Shifter(mask, multiplier)
    coords_in_input = shifter.get_new_coords()
    sy_permuted = ndi.map_coordinates(sy_image, coords_in_input)

    # sy_permuted = ndimage.geometric_transform(sy_image, shifter.geometric_shift, mode='nearest')
    return sy_permuted 

def _profile(self) -> np.ndarray:
        """The actual profile array; private attr that is passed to MultiProfile."""
        profile = np.zeros(len(self._multi_x_locations[0]))
        for radius, x, y in zip(self._radii, self._multi_x_locations, self._multi_y_locations):
            profile += ndimage.map_coordinates(self.image_array, [y, x], order=0)
        profile /= self.num_profiles
        return profile 

def test_map_coordinates(x):
    coordinates = np.stack([np.arange(10) ** 2] * x.ndim)
    coordinates = np.random.randint(0,min(x.shape),(x.ndim,100))
    print(coordinates.shape, x.shape)
    out1 = map_coordinates(x, coordinates, interpolation="linear")
    out2 = ndimage.map_coordinates(x, coordinates, order=1, prefilter=False)
    return out1, out2 

def get_values(self, coordi=None, interpolation_mode=0, border='constant', cval=0.0):
        """
        This function returns the intensity value of the image at the position coordi (can be a list of coordinates).
        :param coordi: continuouspix
        :param interpolation_mode: 0=nearest neighbor, 1= linear, 2= 2nd-order spline, 3= 2nd-order spline, 4= 2nd-order spline, 5= 5th-order spline
        :return: intensity values at continuouspix with interpolation_mode
        """
        return map_coordinates(self.data, coordi, output=np.float32, order=interpolation_mode, mode=border, cval=cval) 

def stretch_mat_creation(refcc, str_range=0.01, nstr=1001):
    """ Matrix of stretched instance of a reference trace.

    From the MIIC Development Team (eraldo.pomponi@uni-leipzig.de)
    The reference trace is stretched using a cubic spline interpolation
    algorithm form ``-str_range`` to ``str_range`` (in %) for totally
    ``nstr`` steps.
    The output of this function is a matrix containing the stretched version
    of the reference trace (one each row) (``strrefmat``) and the corresponding
    stretching amount (`strvec```).
    :type refcc: :class:`~numpy.ndarray`
    :param refcc: 1d ndarray. The reference trace that will be stretched
    :type str_range: float
    :param str_range: Amount, in percent, of the desired stretching (one side)
    :type nstr: int
    :param nstr: Number of stretching steps (one side)
    :rtype: :class:`~numpy.ndarray` and float
    :return: **strrefmat**: 2d ndarray of stretched version of the reference trace.
    :rtype: float
    :return: **strvec**: List of float, stretch amount for each row of ``strrefmat``
    """

    n = len(refcc)
    samples_idx = np.arange(-n // 2 + 1, n // 2 + 1)
    strvec = np.linspace(1 - str_range, 1 + str_range, nstr)
    str_timemat = samples_idx / strvec[:,None] + n // 2
    tiled_ref = np.tile(refcc,(nstr,1))
    coord = np.vstack([(np.ones(tiled_ref.shape) * np.arange(tiled_ref.shape[0])[:,None]).flatten(),str_timemat.flatten()])
    strrefmat = map_coordinates(tiled_ref, coord)
    strrefmat = strrefmat.reshape(tiled_ref.shape)
    return strrefmat, strvec 

def resample_images(images, target_shape, ratio="same", order=1, mode="nearest"):
    output_images = np.zeros(
        (len(images), images[0].shape[0]) + target_shape, dtype=images[0].dtype
    )

    for i, image in enumerate(images):

        # the first step is to resample to fit it into the target shape
        if ratio == "same":
            width_change = target_shape[1] / image.shape[2]
            height_change = target_shape[0] / image.shape[1]
            change = min(width_change, height_change)

        x = np.linspace(0, image.shape[1] - 1, int(np.ceil(image.shape[1] * change)))
        y = np.linspace(0, image.shape[2] - 1, int(np.ceil(image.shape[2] * change)))
        coordinates = np.stack(np.meshgrid(x, y))
        coordinates = np.stack([coordinates[0].reshape(-1), coordinates[1].reshape(-1)])
        print(coordinates)
        new_image = np.stack(
            [
                ndimage.map_coordinates(channel, coordinates, order=order, mode=mode)
                .reshape((len(y), len(x)))
                .T
                for channel in image
            ]
        )
        # now we position this image into the output
        middle_height = (target_shape[0] - len(x)) // 2
        middle_width = (target_shape[1] - len(y)) // 2
        output_images[
            i,
            :,
            middle_height : middle_height + len(x),
            middle_width : middle_width + len(y),
        ] = new_image
    return output_images 

def local_z_distorted(self, x, y):
        if get_distorted_surface is None:
            return 0
        coords = np.array(
            [(x-self.limPhysX[0]) /
             (self.limPhysX[1]-self.limPhysX[0]) * (self.warpNX-1),
             (y-self.limPhysY[0]) /
             (self.limPhysY[1]-self.limPhysY[0]) * (self.warpNY-1)])
# coords.shape = (2, self.nrays)
        z = ndimage.map_coordinates(self.warpSplineZ, coords, prefilter=True)
        return z 

def local_n_distorted(self, x, y):
        if get_distorted_surface is None:
            return
#        a = np.zeros_like(x)
#        b = np.ones_like(x)
        coords = np.array(
            [(x-self.limPhysX[0]) /
             (self.limPhysX[1]-self.limPhysX[0]) * (self.warpNX-1),
             (y-self.limPhysY[0]) /
             (self.limPhysY[1]-self.limPhysY[0]) * (self.warpNY-1)])
        a = ndimage.map_coordinates(self.warpSplineA, coords, prefilter=True)
        b = ndimage.map_coordinates(self.warpSplineB, coords, prefilter=True)
        return b, -a 

def fit_avg_z(cor_id, cor_id_xy, cor_img):
    score = map_coordinates(cor_img, [cor_id[:, 1], cor_id[:, 0]])
    c = np.sqrt((cor_id_xy ** 2).sum(1))
    cor_id_v = ((cor_id[:, 1] + 0.5) / cor_img.shape[0] - 0.5) * np.pi
    z = c * np.tan(cor_id_v)
    fit_z = (z * score).sum() / score.sum()
    return fit_z 

def restore_spline_arrays(
            self, pickleName, findNewImax=True, IminCutOff=1e-50):
        if sys.version_info < (3, 1):
            kw = {}
        else:
            kw = dict(encoding='latin1')
        with open(pickleName, 'rb') as f:
            Imax, savedSplines = pickle.load(f, **kw)
        try:
            if findNewImax:
                ind = [i for i in range(self.Es.shape[0]) if
                       self.eMinRays <= self.Es[i] <= self.eMaxRays]
                if len(ind) == 0:
                    fact = self.eN / (self.eMax-self.eMin)
                    ind = [(self.eMinRays-self.eMin) * fact,
                           (self.eMaxRays-self.eMin) * fact]
                elif len(ind) == 1:
                    ind = [ind[0], ind[0]]
                coords = np.mgrid[ind[0]:ind[-1],
                                  self.extraRows:self.nx+1+self.extraRows,
                                  self.extraRows:self.nz+1+self.extraRows]

                I = ndimage.map_coordinates(
                    savedSplines[0], coords, order=self.order,
                    prefilter=self.prefilter)
                Imax = I.max()
                _eMinRays = self.Es[min(np.nonzero(I > Imax*IminCutOff)[0])]
                if _eMinRays > self.eMinRays:
                    self.eMinRays = _eMinRays
                    print('eMinRays has been corrected up to {0}'.format(
                        _eMinRays))
                _eMaxRays = self.Es[max(np.nonzero(I > Imax*IminCutOff)[0])]
                if _eMaxRays < self.eMaxRays:
                    self.eMaxRays = _eMaxRays
                    print('eMaxRays has been corrected down to {0}'.format(
                        _eMaxRays))
        except ValueError:
            pass
        return savedSplines, Imax 

def intensities_on_mesh(self):
        Is = []
        coords = np.mgrid[0:self.Es.shape[0],
                          self.extraRows:self.nx+1+self.extraRows,
                          self.extraRows:self.nz+1+self.extraRows]
        for a in self.splines:
            aM = ndimage.map_coordinates(a, coords, order=self.order,
                                         prefilter=self.prefilter)
            Is.append(aM)
        return Is 

def plot_NOM_3D(fname):
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator, FormatStrFormatter

    xL, yL, zL = np.loadtxt(fname+'.dat', unpack=True)
    nX = (yL == yL[0]).sum()
    nY = (xL == xL[0]).sum()
    x = xL.reshape((nY, nX))
    y = yL.reshape((nY, nX))
    z = zL.reshape((nY, nX))
    x1D = xL[:nX]
    y1D = yL[::nX]
#    z += z[::-1, :]
    zmax = abs(z).max()

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False, alpha=0.5)
    ax.set_zlim(-zmax, zmax)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

    fig.colorbar(surf, shrink=0.5, aspect=5)

    splineZ = ndimage.spline_filter(z.T)
    nrays = 1e3
    xnew = np.random.uniform(x1D[0], x1D[-1], nrays)
    ynew = np.random.uniform(y1D[0], y1D[-1], nrays)
    coords = np.array([(xnew-x1D[0]) / (x1D[-1]-x1D[0]) * (nX-1),
                       (ynew-y1D[0]) / (y1D[-1]-y1D[0]) * (nY-1)])
    znew = ndimage.map_coordinates(splineZ, coords, prefilter=True)
    ax.scatter(xnew, ynew, znew, c=znew, marker='o', color='gray', s=50,
               cmap=cm.coolwarm)

    fig.savefig(fname+'_3d.png')
    plt.show() 

def get_affine_subwindow(self,img, pos,sc, rot, window_sz):
        def simiparam2mat(tx,ty,rot,s):
            sn,cs=np.sin(rot),np.cos(rot)
            p=[tx,ty,s[0]*cs,-s[1]*sn,s[0]*sn,s[1]*cs]
            return p

        def interp2(img, Xq, Yq):
            wimg = map_coordinates(img, [Xq.ravel(), Yq.ravel()], order=1, mode='wrap')
            wimg = wimg.reshape(Xq.shape)
            return wimg

        def mwarpimg(img,p,sz):
            imsz=img.shape
            w,h=sz
            x,y=np.meshgrid(np.arange(1,w+1)-w//2,np.arange(1,h+1)-h//2)
            tmp1=np.zeros((h*w,3))
            tmp1[:,0]=1
            tmp1[:,1]=x.flatten()
            tmp1[:,2]=y.flatten()
            tmp2=np.array([[p[0],p[1]],[p[2],p[4]],[p[3],p[5]]])
            tmp3=tmp1.dot(tmp2)
            tmp3=np.clip(tmp3,a_min=1,a_max=None)
            tmp3[:, 0] = np.clip(tmp3[:, 0], a_min=None,a_max=imsz[1])
            tmp3[:, 1] = np.clip(tmp3[:, 1], a_min=None,a_max=imsz[0])
            pos=np.reshape(tmp3,(h,w,2))
            c=img.shape[2]
            wimg=np.zeros((sz[1],sz[0],c))
            pos=pos-1
            for i in range(c):
                wimg[:,:,i]=interp2(img[:,:,i],pos[:,:,1],pos[:,:,0])
            return wimg
        x,y=pos
        param0=simiparam2mat(x,y,rot,sc)
        out=mwarpimg(img.astype(np.float32),param0,window_sz).astype(np.uint8)
        #cv2.imshow('affine_window',out.astype(np.uint8))
        #cv2.waitKey(1)
        return out 

def stretch_mat_creation(refcc, str_range=0.01, nstr=1001):
    """ Matrix of stretched instance of a reference trace.

    From the MIIC Development Team (eraldo.pomponi@uni-leipzig.de)
    The reference trace is stretched using a cubic spline interpolation
    algorithm form ``-str_range`` to ``str_range`` (in %) for totally
    ``nstr`` steps.
    The output of this function is a matrix containing the stretched version
    of the reference trace (one each row) (``strrefmat``) and the corresponding
    stretching amount (`strvec```).
    :type refcc: :class:`~numpy.ndarray`
    :param refcc: 1d ndarray. The reference trace that will be stretched
    :type str_range: float
    :param str_range: Amount, in percent, of the desired stretching (one side)
    :type nstr: int
    :param nstr: Number of stretching steps (one side)
    :rtype: :class:`~numpy.ndarray` and float
    :return: **strrefmat**: 2d ndarray of stretched version of the reference trace.
    :rtype: float
    :return: **strvec**: List of float, stretch amount for each row of ``strrefmat``
    """

    n = len(refcc)
    samples_idx = np.arange(-n // 2 + 1, n // 2 + 1)
    strvec = np.linspace(1 - str_range, 1 + str_range, nstr)
    str_timemat = samples_idx / strvec[:,None] + n // 2
    tiled_ref = np.tile(refcc,(nstr,1))
    coord = np.vstack([(np.ones(tiled_ref.shape) * np.arange(tiled_ref.shape[0])[:,None]).flatten(),str_timemat.flatten()])
    strrefmat = map_coordinates(tiled_ref, coord)
    strrefmat = strrefmat.reshape(tiled_ref.shape)
    return strrefmat, strvec 

def test_map_coordinates02(self):
        data = numpy.array([[4, 1, 3, 2],
                               [7, 6, 8, 5],
                               [3, 5, 3, 6]])
        idx = numpy.indices(data.shape, numpy.float64)
        idx -= 0.5
        for order in range(0, 6):
            out1 = ndimage.shift(data, 0.5, order=order)
            out2 = ndimage.map_coordinates(data, idx,
                                                     order=order)
            assert_array_almost_equal(out1, out2) 

def sample_equirec(e_img, coor_xy, order):
    w = e_img.shape[1]
    coor_x, coor_y = np.split(coor_xy, 2, axis=-1)
    pad_u = np.roll(e_img[[0]], w // 2, 1)
    pad_d = np.roll(e_img[[-1]], w // 2, 1)
    e_img = np.concatenate([e_img, pad_d, pad_u], 0)
    return map_coordinates(e_img, [coor_y, coor_x],
                           order=order, mode='wrap')[..., 0]