Improved PGA method

[isaacgerg]

In talking with Hayden Callow, I was able to improve your existing autofocus2 method.

In it, I use the ML estimator which is unbiased. Incidentally, I find this method performs better than the existing method on my SAS data of the SASSED dataset.

Here is my updated version of autofocus2 method.

def autoFocus2(img, win = 'auto', win_params = [100,0.5]):
#  Originally written by Douglas Macdonald. 
#  Based upon:
#  D. Wahl, P. Eichel, D. Ghiglia, and J. Jakowatz, C.V., \Phase gradient 
#     autofocus-a robust tool for high resolution sar phase correction,"
#     Aerospace and Electronic Systems, IEEE Transactions on, vol. 30,
#     pp. 827{835, Jul 1994.
# Assumes image azimuth is vertical dimension and range increases in horizontal dim.

    #Derive parameters
    npulses = int(img.shape[0])
    nsamples = int(img.shape[1])

    #Initialize loop variables
    img_af = 1.0*img
    max_iter = 30
    af_ph = 0
    rms = []

    #Compute phase error and apply correction
    for iii in range(max_iter):

        #Find brightest azimuth sample in each range bin
        index = np.argsort(np.abs(img_af), axis=0)[-1]

        #Circularly shift image so max values line up
        f = np.zeros(img.shape)+0j
        for i in range(nsamples):
            f[:,i] = np.roll(img_af[:,i], int(npulses/2-index[i]))

        if win == 'auto':
            #Compute window width
            s = np.sum(f*np.conj(f), axis = -1)
            s = 10*np.log10(s/s.max())
            #For first two iterations use all azimuth data
            #and half of azimuth data, respectively
            if iii == 0:
                width = npulses
            #For all other iterations, use twice the 10 dB threshold
            else:
                width = np.sum(s>-10) 
            window = np.arange(npulses/2-width/2,npulses/2+width/2)
        else:
            #Compute window width using win_params if win not set to 'auto'
            width = int(win_params[0]*win_params[1]**iii)
            window = np.arange(npulses/2-width/2,npulses/2+width/2)
            if width<5:
                break

        window = window.astype('int')
        #Window image
        g = np.zeros(img.shape)+0j
        g[window] = f[window]

        #Fourier Transform
        G = sig.ift(g, ax=0)

        # ML method
        phi_dot = np.angle(np.sum(np.conj(G[:-1, :]) * G[1:, :], axis=1)) 
        phi = np.concatenate([[0],  np.cumsum(phi_dot)])

        #Remove linear trend
        t = np.arange(0,nsamples)
        slope, intercept, r_value, p_value, std_err = linregress(t,phi)
        line = slope*t+intercept
        phi = phi-line
        rms.append(np.sqrt(np.mean(phi**2)))

        if win == 'auto':
            if rms[iii]<0.01:
                break

        #Apply correction
        phi2 = np.tile(np.array([phi]).T,(1,nsamples))
        IMG_af = sig.ift(img_af, ax=0)
        IMG_af = IMG_af*np.exp(-1j*phi2)
        img_af = sig.ft(IMG_af, ax=0)

        #Store phase
        af_ph += phi

    print('number of iterations: %i'%(iii+1))

    return(img_af, np.flip(af_ph), rms)