def Convolve(image1, image2, MinPad=True, pad=True):
    """
    Convolves image1 with image2.

    :param image1: 2D image array
    :param image2: 2D image array
    :param MinPad: whether to use minimal padding
    :param pad: whether to pad the array
    """
    #The size of the images:
    r1, c1 = image1.shape
    r2, c2 = image2.shape

    if MinPad:
        r = r1 + r2
        c = c1 + c2
    else:
        r = 2*max(r1,r2)
        c = 2*max(c1,c2)
    
    #or in power of two
    if pad:
        pr2 = int(m.log(r)/m.log(2.) + 1.)
        pc2 = int(m.log(c)/m.log(2.) + 1.)
        rOrig = r
        cOrig = c
        r = 2**pr2
        c = 2**pc2
    
    fftimage = fft2(image1, s=(r,c))*fft2(image2[::-1,::-1],s=(r,c))

    if pad:
        return (ifft2(fftimage))[:rOrig,:cOrig].real
    else:
        return (ifft2(fftimage)).real

def InitVelField(_N, _M, _h, h, dt, rho=1.0, mu=1.0, DeltaType=0):
    WideLambda = zeros((_N, _M), float64)
    ShortLambda = zeros((_N, _M), float64)
    IB_c.InitWideLaplacian(_N, _M, _h, WideLambda)
    IB_c.InitShortLaplacian(_N, _M, _h, ShortLambda)
    DxSymbol = InitDxSymbol(_N, _M, _h)
    DySymbol = InitDySymbol(_N, _M, _h)

    r = int(ceil(3.0 * h / _h))

    fx = zeros((_N, _M), float64)
    for j in range(-r, r + 1):
        deltx = Delta(h, j * _h, DeltaType)
        for k in range(-r, r + 1):
            delt = deltx * Delta(h, k * _h, DeltaType) * 1.0
            fx[j % _N][k % _M] = fx[j % _N][k % _M] + delt
    #       print j%_N, k%_M, fx[j%_N][k%_M]

    fx, fy = fft2(dt * fx), zeros((_N, _M), float64)

    P = Solve_P_Hat(dt, WideLambda, DxSymbol, DySymbol, fx, fy)
    P[0, 0] = 0.0

    u, v = Solve_uv_Hat(dt, ShortLambda, DxSymbol, DySymbol, P, fx, fy, rho, mu)
    u = 1.0 * ifft2(u).real
    v = 1.0 * ifft2(v).real
    #    P = ifft2(P).real

    Fx1 = array(zeros((_N, _M), float64))
    Fy1 = array(zeros((_N, _M), float64))

    IB_c.WholeGridSpread(u, float(h), float(_h), int(r), Fx1, DeltaType)
    IB_c.WholeGridSpread(v, float(h), float(_h), int(r), Fy1, DeltaType)

    fy = zeros((_N, _M), float64)
    for j in range(-r, r + 1):
        deltx = Delta(h, j * _h, DeltaType)
        for k in range(-r, r + 1):
            delt = deltx * Delta(h, k * _h, DeltaType) * 1.0
            fy[j % _N][k % _M] = fy[j % _N][k % _M] + delt
    #       print j%_N, k%_M, fx[j%_N][k%_M]

    fx, fy = zeros((_N, _M), float64), fft2(dt * fy)

    P = Solve_P_Hat(dt, WideLambda, DxSymbol, DySymbol, fx, fy)
    P[0, 0] = 0.0

    u, v = Solve_uv_Hat(dt, ShortLambda, DxSymbol, DySymbol, P, fx, fy, rho, mu)
    u = 1.0 * ifft2(u).real
    v = 1.0 * ifft2(v).real

    Fx2 = array(zeros((_N, _M), float64))
    Fy2 = array(zeros((_N, _M), float64))

    IB_c.WholeGridSpread(u, float(h), float(_h), int(r), Fx2, DeltaType)
    IB_c.WholeGridSpread(v, float(h), float(_h), int(r), Fy2, DeltaType)

    return Fx1, Fy1, Fx2, Fy2

def get_stat_spin_struct(filenames,nsamp):
    """
    Gets the static structure factor flatened.
    The q-vectors are given by:
    q=arange(L)
    qx,qy=meshgrid(q,q)
    qx=qx.flatten()
    qy=qy.flatten()
    """
    if type(filenames)!=list:
        filenames=[filenames]
    Sq=load.get_quantity(filenames,nsamp)
    params=load.get_attr(filenames[0])
    Lx=int(params['L'])
    Ly=int(params['L'])
    hLx=int(Lx/2)
    hLy=int(Ly/2)
    N=Lx*Ly
    if Sq.shape[2]==N:
        # old file format, struct stored in Fourier components
        Sqxx=np.reshape(0.25*(Sq[:,1,:]+Sq[:,2,:]+Sq[:,3,:]+Sq[:,4,:]),(Sq.shape[0],Lx,Ly))
        Sqyy=np.reshape(0.25*(Sq[:,1,:]+Sq[:,2,:]-Sq[:,3,:]-Sq[:,4,:]),(Sq.shape[0],Lx,Ly))
        Sqzz=np.reshape(Sq[:,0,:],(Sq.shape[0],Lx.Ly))
        Srxx=fft.fftshift(fft.fft2(Sqxx,axes=(1,2)),axes=(1,2))/N
        Sryy=fft.fftshift(fft.fft2(Sqyy,axes=(1,2)),axes=(1,2))/N
        Srzz=fft.fftshift(fft.fft2(Sqzz,axes=(1,2)),axes=(1,2))/N
    else :
        # new file format, struct stored as real space site pairs.
        rx,ry=np.meshgrid(np.arange(Lx,dtype=int),np.arange(Ly,dtype=int))
        rx=rx.ravel()
        ry=ry.ravel()
        rix,rjx=np.meshgrid(rx,rx)
        riy,rjy=np.meshgrid(ry,ry)
        rijx=rjx-rix
        rijy=rjy-riy
        rijx[rijx>=hLx]-=Lx
        rijx[rijx<-hLx]+=Lx
        rijy[rijy>=hLy]-=Ly
        rijy[rijy<-hLy]+=Ly
        rijx=rijx.ravel()
        rijy=rijy.ravel()
        Sr=np.zeros((Sq.shape[0],5,N),dtype=complex)
        for samp in range(Sq.shape[0]):
            for t in range(N):
                Sr[samp,:,t]=np.sum(Sq[samp,:,np.where((rijy+hLy)*Lx+rijx+hLx==t)[0]],axis=0)/N
        Srxx=np.zeros((Sq.shape[0],Lx,Ly),dtype=complex)
        Sryy=np.zeros((Sq.shape[0],Lx,Ly),dtype=complex)
        Srzz=np.zeros((Sq.shape[0],Lx,Ly),dtype=complex)
        for samp in range(Sq.shape[0]):
            Srxx[samp,:,:]=np.reshape(0.25*np.sum(Sr[samp,1:,:],axis=0),(Lx,Ly))
            Sryy[samp,:,:]=np.reshape(0.25*(np.sum(Sr[samp,1:3,:],axis=0)-np.sum(Sr[samp,3:,:],axis=0)),(Lx,Ly))
            Srzz[samp,:,:]=np.reshape(Sr[samp,0,:],(Lx,Ly))
        Sqxx=fft.ifft2(fft.fftshift(Srxx,axes=(1,2)),axes=(1,2))*N
        Sqyy=fft.ifft2(fft.fftshift(Sryy,axes=(1,2)),axes=(1,2))*N
        Sqzz=fft.ifft2(fft.fftshift(Srzz,axes=(1,2)),axes=(1,2))*N
    return (Sqxx,Sqyy,Sqzz),(Srxx,Sryy,Srzz)

def update_g(F, G_arg, C, iterations=1, eps=1e-9):
    G = G_arg.copy()
    g = ifft2(G)
    f_rev = np.rot90(ifft2(F), k=2)
    for k in range(int(iterations)):
        blur = (C / (G * F + eps)) * fft2(f_rev)
        G = fftconvolve(G, blur, "same")
        g = ifft2(blur) * g
        G = fft(g)
    return G, ifft2(G)

def do_n_extended_fits(nfits, xsh, ysh, imsize,  gaussfit=False,
        maxoff=None, return_error=False, powerlaw=2.0, noise=1.0,
        unsharp_mask=False, smoothfactor=5, zeropad=0,
        shift_func=cross_correlation_shifts, sfkwargs={},
        doplot=False,
        **kwargs):

    try: 
        import progressbar
        widgets = [progressbar.FormatLabel('Processed: %(value)d offsets in %(elapsed)s)'), progressbar.Percentage()]
        progress = progressbar.ProgressBar(widgets=widgets)
    except ImportError:
        def progress(x):
            yield x

    image = make_extended(imsize, powerlaw=powerlaw)
    if zeropad > 0:
        newsize = [s+zeropad for s in image.shape]
        ylen,xlen = newsize
        xcen = xlen/2-(1-xlen%2) 
        ycen = ylen/2-(1-ylen%2) 
        newim = np.zeros(newsize)
        newim[ycen-image.shape[0]/2:ycen+image.shape[0]/2, xcen-image.shape[1]/2:xcen+image.shape[1]/2] = image
        image = newim


    if unsharp_mask:
        from AG_fft_tools import smooth
        offsets = []
        for ii in progress(xrange(nfits)):
            inim = image-smooth(image,smoothfactor)
            offim = make_offset_extended(image, xsh, ysh, noise=noise, **kwargs)
            offim -= smooth(offim,smoothfactor)
            offsets.append( shift_func( inim, offim,  return_error=return_error, **sfkwargs) )
    else:
        offsets = []
        if doplot:
            import pylab
            pylab.figure(3); pylab.subplot(221); pylab.imshow(image-image.mean()); pylab.subplot(222); pylab.imshow(offim-offim.mean())
            #subplot(223); pylab.imshow((abs(fft2(image-image.mean())*conj(fft2(offim-offim.mean())))))
            pylab.subplot(223); pylab.imshow(abs(ifft2((fft2(image)*conj(fft2(offim))))))
            pylab.subplot(224); pylab.imshow(abs(ifft2((fft2(image-image.mean())*conj(fft2(offim-offim.mean()))))))
            draw()
        for ii in progress(xrange(nfits)):
            offim = make_offset_extended(image, xsh, ysh, noise=noise, **kwargs)
            offsets.append( shift_func( 
                image,
                offim,
                return_error=return_error, **sfkwargs)
                )

    return offsets

def F(wt, t, nu, KX, KY, K):

    # calculate psi in fourier space
    psit = -wt/K

    # calculate the derivatives of psi and w
    # in real space
    psi_x = np.real(fft.ifft2(1.j*KX*psit))
    psi_y = np.real(fft.ifft2(1.j*KY*psit))
    w_x   = np.real(fft.ifft2(1.j*KX*wt))
    w_y   = np.real(fft.ifft2(1.j*KY*wt))

    return fft.fft2(w_x*psi_y - w_y*psi_x) - nu*K*wt

    def setUp(self):
        # Input
        self.data = hdf5.read_data("example_input.h5")
        C, T, Z, Y, X = self.data[0].shape
        self.data[0].shape = Y, X
        self.small_data = self.data[0][350:400, 325:375]

        # Input drifted by known value
        self.data_drifted = hdf5.read_data("example_drifted.h5")
        C, T, Z, Y, X = self.data_drifted[0].shape
        self.data_drifted[0].shape = Y, X

        # Input drifted by random value
        z = 1j  # imaginary unit
        self.deltar = numpy.random.uniform(-100, 100)
        self.deltac = numpy.random.uniform(-100, 100)
        nr, nc = self.data[0].shape
        array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
        array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
        Nr = fft.ifftshift(array_nr)
        Nc = fft.ifftshift(array_nc)
        [Nc, Nr] = numpy.meshgrid(Nc, Nr)
        self.data_random_drifted = fft.ifft2(fft.fft2(self.data[0]) * numpy.power(math.e,
        				z * 2 * math.pi * (self.deltar * Nr / nr + self.deltac * Nc / nc)))

        # Noisy inputs
        noise = random.normal(0, 3000, self.data[0].size)
        noise_array = noise.reshape(self.data[0].shape[0], self.data[0].shape[1])

        self.data_noisy = self.data[0] + noise_array
        self.data_drifted_noisy = self.data_drifted[0] + noise_array
        self.data_random_drifted_noisy = self.data_random_drifted + noise_array

        # Small input drifted by random value
        self.small_deltar = numpy.random.uniform(-10, 10)
        self.small_deltac = numpy.random.uniform(-10, 10)
        nr, nc = self.small_data.shape
        array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
        array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
        Nr = fft.ifftshift(array_nr)
        Nc = fft.ifftshift(array_nc)
        [Nc, Nr] = numpy.meshgrid(Nc, Nr)
        self.small_data_random_drifted = fft.ifft2(fft.fft2(self.small_data) * numpy.power(math.e,
        				z * 2 * math.pi * (self.small_deltar * Nr / nr + self.small_deltac * Nc / nc)))

        # Small noisy inputs
        small_noise = random.normal(0, 3000, self.small_data.size)
        small_noise_array = small_noise.reshape(self.small_data.shape[0], self.small_data.shape[1])

        self.small_data_noisy = self.small_data + small_noise_array
        self.small_data_random_drifted_noisy = self.small_data_random_drifted + small_noise_array

def task3_2():
    print('3.2')
    low_pass = normalize_intensity(imread(LOW_PASS))
    high_pass = normalize_intensity(imread(HIGH_PASS))

    img_freq_dom = fft2(a=normalize_intensity(imread(BRICKS_2)))
    apply_low_pass = img_freq_dom * low_pass
    ifft_res = ifft2(a=apply_low_pass)
    output_path = os.path.join(OUTPUT_DIR, "3_2_low_pass_" + os.path.split(BRICKS_2)[-1])
    imsave(output_path, abs(ifft_res))

    apply_high_pass = img_freq_dom * high_pass
    ifft_res = ifft2(a=apply_high_pass)
    output_path = os.path.join(OUTPUT_DIR, "3_2_high_pass_" + os.path.split(BRICKS_2)[-1])
    imsave(output_path, abs(ifft_res))

def whiten_olsh_lee_inner(image, f_0=None, central_clip=(None, None), normalize_pre=True, normalize_post=True,
                          no_filter=False):
    height, width = image.shape
    assert height % 2 == 0 and width % 2 == 0, "image must have even size!"
    if normalize_pre:
        image = image - image.mean()  # I personally think this is useless, since rho will make (0,0) freq compoenent 0.
        std_im = image.std(ddof=1)
        assert std_im != 0, "constant image unsupported!"
        image /= std_im

    fx, fy = np.meshgrid(np.arange(-height / 2, height / 2), np.arange(-width / 2, width / 2), indexing='ij')
    rho = np.sqrt(fx * fx + fy * fy)
    if f_0 is None:
        f_0 = 0.4 * (height + width) / 2
    filt = rho * np.exp(-((rho / f_0) ** 4))

    im_f = fft2(image)
    if not no_filter:
        fft_filtered_old = im_f * ifftshift(filt)
    else:  # hack to only lower frequency response.
        print('no real filtering!')
        fft_filtered_old = im_f
    fft_filtered_old = fftshift(fft_filtered_old)
    if central_clip != (None, None):
        fft_filtered_old = fft_filtered_old[height // 2 - central_clip[0] // 2:height // 2 + central_clip[0] // 2,
                           width // 2 - central_clip[1] // 2:width // 2 + central_clip[1] // 2]
    im_out = np.real(ifft2(ifftshift(fft_filtered_old)))
    # I believe since the rho at the (0,0) frequency part is zero, then the whole image should be zero as well.
    # so explicit DC removing is useless.
    if normalize_post:
        assert abs(im_out.mean()) < 1e-6  # should be extremely small.
        std_im_out = im_out.std(ddof=1)
    else:
        std_im_out = 1
    return im_out / std_im_out

def high_pass_filter(img, filtersize=10):
    """
    A FFT implmentation of high pass filter from pyKLIP.
    Args:
        img: a 2D image
        filtersize: size in Fourier space of the size of the space. In image space, size=img_size/filtersize
    Returns:
        filtered: the filtered image
    """
    if filtersize == 0:
        return img
    # mask NaNs
    nan_index = np.where(np.isnan(img))
    img[nan_index] = 0

    transform = fft.fft2(img)

    # coordinate system in FFT image
    u,v = np.meshgrid(fft.fftfreq(transform.shape[1]), fft.fftfreq(transform.shape[0]))
    # scale u,v so it has units of pixels in FFT space
    rho = np.sqrt((u*transform.shape[1])**2 + (v*transform.shape[0])**2)
    # scale rho up so that it has units of pixels in FFT space
    # rho *= transform.shape[0]
    # create the filter
    filt = 1. - np.exp(-(rho**2/filtersize**2))

    filtered = np.real(fft.ifft2(transform*filt))

    # restore NaNs
    filtered[nan_index] = np.nan
    img[nan_index] = np.nan

    return filtered

    def _simulate_image(self):
        """
        Generates the fake output.
        """
        with self._acquisition_init_lock:
            pos = self.align.position.value
            logging.debug("Simulating image shift by %s", pos)
            ac, bc = pos.get("a"), pos.get("b")
            ang = math.radians(135)
            # AB->XY
            xc = -(ac * math.sin(ang) + bc * math.cos(ang))
            yc = -(ac * math.cos(ang) - bc * math.sin(ang))
            pixelSize = self.fake_img.metadata[model.MD_PIXEL_SIZE]
            self.fake_img.metadata[model.MD_ACQ_DATE] = time.time()
            x_pxs = xc / pixelSize[0]
            y_pxs = yc / pixelSize[1]

            # Image shifted based on LensAligner position
            z = 1j  # imaginary unit
            self.deltar = x_pxs
            self.deltac = y_pxs
            nr, nc = self.fake_img.shape
            array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
            array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
            Nr = fft.ifftshift(array_nr)
            Nc = fft.ifftshift(array_nc)
            [Nc, Nr] = numpy.meshgrid(Nc, Nr)
            sim_img = fft.ifft2(fft.fft2(self.fake_img) * numpy.power(math.e,
                            z * 2 * math.pi * (self.deltar * Nr / nr + self.deltac * Nc / nc)))
            output = model.DataArray(abs(sim_img), self.fake_img.metadata)
            return output

def acf(x, axes=(0,1)):
    """
    2D ACF using fft

    Inputs:

    x is ndarray with shape (nt, nx)
    """
    from numpy.fft import fft2, ifft2, fftshift

    if x.ndim == 1:
        x = x[:,None]
    elif x.ndim == 2:
        pass
    else:
        raise NotImplementedError

    nt, nx = x.shape

    padding = np.zeros((nt, nx))

    x = np.concatenate((x, padding))

    fx = fft2(x, axes=axes)
    ac=  np.real(ifft2(fx * np.conj(fx), axes=axes))[:(nt-10),:] / nx / np.arange(nt, 10, -1)[:,None]

    ac = ac[:nt/2, :nx/2]

    return ac

def calculate(imageFFT, phaseSymmetry, symmetryTotal, amplitudeTotal):
	for a in range(1,orientations + 1):
		#print(a)
		symmetryAtOrientation = np.zeros((image.shape[0], image.shape[1]))
		amplitudeAtOrientation = np.zeros((image.shape[0], image.shape[1]))
		for n in range(scales):		
			kernel = bank[a - 1][n]		
			convolved = imageFFT * kernel
			s1 = np.array(kernel.shape)
			s2 = np.array(imageFFT.shape)
			convolved = fft.ifft2(convolved)
			
			evens = np.real(convolved)
			odds = np.imag(convolved)
			amplitude = np.sqrt(np.power(evens, 2) + np.power(odds, 2)) 
			
			amplitudeAtOrientation += amplitude
			symmetryAtOrientation += np.maximum((np.abs(evens) - np.abs(odds))  , floor)
		
		amplitudeTotal += np.add(amplitudeAtOrientation, 0.00001)
		symmetryTotal += symmetryAtOrientation

	phaseSymmetry = np.divide(symmetryTotal , amplitudeTotal)
	phaseSymmetry[mask < 0.2] = 0.0
	return phaseSymmetry

    def compute(self, scene: Scene):
        """ Compute optical irradiance map
        Computation proccedure:
            1) convert radiance to irradiance
            2) apply lens and macular transmittance
            3) apply off-axis fall-off (cos4th)
            4) apply optical transfert function

        Args:
            scene (pyEyeBall.Scene): instance of Scene class, containing the radiance and other scene information

        Examples:
            >>> oi = Optics()
            >>> oi.compute(Scene())
        """
        # set field of view and wavelength samples
        self.fov = scene.fov
        scene.wave = self._wave
        self.dist = scene.dist

        # compute irradiance
        self.photons = pi / (1 + 4 * self.f_number**2 * (1 + abs(self.magnification))**2) * scene.photons

        # apply ocular transmittance
        self.photons *= self.ocular_transmittance

        # apply the relative illuminant (off-axis) fall-off: cos4th function
        x, y = self.spatial_support
        s_factor = np.sqrt(self.image_distance**2 + x**2 + y**2)
        self.photons *= (self.image_distance / s_factor[:, :, None])**4

        # apply optical transfer function of the optics
        for ii in range(self.wave.size):
            otf = fftshift(self.otf(self._wave[ii], self.frequency_support_x, self.frequency_support_y))
            self.photons[:, :, ii] = np.abs(ifftshift(ifft2(otf * fft2(fftshift(self.photons[:, :, ii])))))

def remove_lowest(shape, dim=1):
    x_min = amin(real(fft.ifft2(change_n(shape.x_hat, 1E2 * ones(FFT_NDIM)), 
                 axes=FFT_AXES))[...,dim])

    g = 1/(5*absolute(shape.x[...,dim] - x_min) + 0.5) - 0.5
    g[g<0] = 0
    return g

def circulantPinkNoiseIterator(p, powerspectrum_sample_size, patchSampler, orientation='F'):
    """

    Samples pxp patches from circulant pink noise with power spectrum from patchSampler.

    The images are vectorized in FORTRAN/MATLAB style by default.

    :param p: patch size
    :type p: int
    :param powerspectrum_sample_size: number of patches to sample from sampler for power spectrum
    :type powerspectrum_sample_size: int
    :param patchSampler: Iterator to sample patches from
    :type patchSampler: Iterator
    :param orientation: 'C' (C/Python, row-major) or 'F' (FORTRAN/MATLAB, column-major) vectorized patches
    :type orientation: string
    :returns: Iterator that samples from all files
    :rtype: Iterator

    """

    patches = zeros((p**2, powerspectrum_sample_size))
    stdout.write('Initialize power spectrum of circulant pink noise generator. This may take a moment...\n')
    stdout.flush()
    for ii in xrange(powerspectrum_sample_size):
        patches[:,ii] = patchSampler.next()
    PATCHES = fft.fft2(patches.reshape(p,p,powerspectrum_sample_size), axes=(0,1))
    powerspec = (PATCHES*PATCHES.conj()).mean(2).real
    stdout.write('Powerspectrum completed.\n')

    while True:
        phi = rand(p,p)*2*pi - pi
        X = fft.ifft2(sqrt(powerspec)*exp(1J*phi), axes=(0,1)).real.flatten(orientation)
        yield X

    return

def cross_corr(img1,img2,mask=None):
    '''Compute the autocorrelation of two images.
        Right now does not take mask into account.
        todo: take mask into account (requires extra calculations)
        input: 
            img1: first image
            img2: second image
            mask: a mask array
        output:
            the autocorrelation of the two images (same shape as the correlated images)
        
    '''
    #if(mask is not None):
    #   img1 *= mask
    #  img2 *= mask
    
    #img1_mean = np.mean( img1.flat )
    #img2_mean = np.mean( img2.flat )
    
    # imgc = fftshift( ifft2(     
    #        fft2(img1/img1_mean -1.0 )*np.conj(fft2( img2/img2_mean -1.0 ))).real )
    
    #imgc = fftshift( ifft2(     
    #        fft2(  img1/img1_mean  )*np.conj(fft2(  img2/img2_mean   ))).real )
    
    imgc = fftshift( ifft2(     
            fft2(  img1  )*np.conj(fft2(  img2  ))).real )
    
    #imgc /= (img1.shape[0]*img1.shape[1])**2
    if(mask is not None):
        maskc = cross_corr(mask,mask)        
        imgc /= np.maximum( 1, maskc )
            
            
    return imgc

    def gradient_helper(self, x, shape, ctr):
        """Not sure exactly what this does yet.

        Parameters
        ----------
        i_t : int
            Epoch index.
        x : np.ndarray (3-d)
            Same shape as *data* for single epoch (nw, ny, nx).
        xcoords : np.ndarray (1-d)
        ycoords : np.ndarray (1-d)

        Returns
        -------
        x : np.ndarray (3-d)
            Shape is (nw, len(ycoords), len(xcoords)).
        """

        # shift necessary to put model onto data coordinates
        offset = yxoffset((self.ny, self.nx), shape, ctr)
        fshift = fft_shift_phasor_2d((self.ny, self.nx), (-offset[0], -offset[1]))
        fshift = np.asarray(fshift, dtype=self.complex_dtype)

        # create output array
        out = np.zeros((self.nw, self.ny, self.nx), dtype=self.dtype)
        out[:, : x.shape[1], : x.shape[2]] = x

        for i in range(self.nw):
            tmp = ifft2(np.conj(self.fftconv[i, :, :] * fshift) * fft2(out[i, :, :]))
            out[i, :, :] = tmp.real

        return out

def delay_fringe_rate(tf_plane,padding=1):
    if isnan(tf_plane).sum() > 0:
        ValueError('*tf_plane* contains NaN values. Please sanitize it before plotting using, e.g. tf_plane[isnan(tf_plane)] == 0.0, or pyautoplot.utilities.set_nan_zero(tf_plane)')
    nt,nf = tf_plane.shape
    padded_plane=zeros((nt,padding*nf),dtype=complex64)
    padded_plane[:,(padding//2):(padding//2+nf)] = tf_plane.data*logical_not(tf_plane.mask)
    return fftshift(ifft2(padded_plane))

def subpixel_shift(img_arr, dx=0.5, inplace=True):
    '''
    Shifts an image by a fraction of a pixel in a random direction,
    with circular boundary conditions.

    Arguments:
        img_arr (2D ndarray): array to subpixel shift in a random direction
        dx (float): distance in pixel to shift by
        inplace (bool): if True, will modify the array inplace

    '''

    f = fft.fft2(img_arr, axes=(0, 1))
    g = np.meshgrid(
        fft.fftfreq(img_arr.shape[-3]),
        fft.fftfreq(img_arr.shape[-2])
    )

    u = npr.normal(0, 1, size=2)
    fk = g[0] * u[0] + g[1] * u[1]
    phi = np.exp(2j * np.pi * dx * fk)

    out = np.real(fft.ifft2(phi[..., np.newaxis] * f, axes=(0, 1)))
    if inplace:
        img_arr[:] = out[:]
    else:
        return out