peak_detection.py 文件源码

def gaussian2D(window):
    """Real 2D Gaussian interpolation for sub pixel shift"""
    #ref on paper
    w = np.ones((3, 3))*(1./9)
    rhs = np.zeros(6)
    M = np.zeros((6,6))
    for i in [-1, 0, 1]:
        for j in [-1, 0, 1]:
            rhs = rhs     +     np.array([i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                          j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        i*i*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                        j*j*w[i+1, j+1]*np.log(np.abs(window[i+1, j+1])),
                                            w[i+1, j+1]*np.log(np.abs(window[i+1, j+1]))], dtype='float')

            M = M + w[i+1, j+1]*np.array([[  i*i,   i*j,   i*i*j,   i*i*i,   i*j*j,   i],
                                          [  i*j,   j*j,   i*j*j,   i*i*j,   j*j*j,   j],
                                          [i*i*j, i*j*j, i*i*j*j, i*i*i*j, i*j*j*j, i*j],
                                          [i*i*i, i*i*j, i*i*i*j, i*i*i*i, i*i*j*j, i*i],
                                          [i*j*j, j*j*j, i*j*j*j, i*i*j*j, j*j*j*j, j*j],
                                          [    i,     j,     i*j,     i*i,     j*j,   1]], dtype='float')
    solution = nl.solve(M, rhs)

    dx = (    solution[2]*solution[1] - 2.0*solution[0]*solution[4])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    dy = (    solution[2]*solution[0] - 2.0*solution[1]*solution[3])/ \
         (4.0*solution[3]*solution[4] -     solution[2]*solution[2])

    return dx, dy