<filename>scratch_AKS/DHC.py<gh_stars>1-10
def self_conj(a):
    return np.real(np.multiply(a,np.conj(a)))

def complex_rotate_image(image, angle):
    image_center = tuple(np.array(image.shape[1::-1]) / 2)
    rot_mat = cv2.getRotationMatrix2D(image_center, angle, 1.0)
    result_real = cv2.warpAffine(np.real(image), rot_mat, image.shape[1::-1], flags=cv2.INTER_LANCZOS4)
    result_imag = cv2.warpAffine(np.imag(image), rot_mat, image.shape[1::-1], flags=cv2.INTER_LANCZOS4)
    return result_real+result_imag*1j

def phase_stretch(im_rd_0_1,stretch,p0=0):
    im_stretch = np.abs(im_rd_0_1)*np.exp(stretch*np.angle(im_rd_0_1)*1j+np.pi*p0*1j)
    return np.fft.fft2(im_stretch)

def phase_stretch_f(U_1_c_f,stretch,p0=0):
    im_rs = np.fft.ifft2(U_1_c_f)
    im_stretch = np.abs(im_rs)*np.exp(stretch*np.angle(im_rs)*1j+np.pi*p0*1j)
    return np.fft.fft2(im_stretch)

def dual_conj(a,b):
    return np.abs(np.multiply(a,np.conj(b)))

def DHC(image,filter_set,phase_on = 1, rotate_on = 1, coeff_12_on = 1, norm_on = 1):
    Nx, Ny = image.shape

    out_coeff = []

    ##Coeff at (1,0) level
    #takes care of image normalization
    S0 = np.zeros((2,1),dtype=np.float64)
    if norm_on == 1:
        S0[0] = np.mean(image)
        norm_im = image-S0[0]
        S0[1] = np.sum(np.square(norm_im))/(Nx*Ny)
        norm_im = norm_im/np.sqrt(Nx*Ny*S0[1])
    else:
        norm_im = image

    out_coeff.append(S0)

    ##Coeff at (1,1) level
    #store the fft of images and im*filters for later
    im_fd_0 = np.fft.fft2(norm_im)

    S1 = np.zeros((8,8),dtype=np.float64)
    im_fd_0_1 = np.zeros((8,8,256,256),dtype=np.complex128)
    for j in range(J):
        for l in range(L):
            im_fd_0_1[j,l] = np.multiply(im_fd_0,filter_set[j,l]) #wavelet already in fft domain not shifted
            S1[j,l]=np.sum(self_conj(im_fd_0_1[j,l])) #normalization choice arb to make order unity

    out_coeff.append(S1)

    ##Coeff at (2,0) level
    S20 = np.zeros((J,J,L,L),dtype=np.float64)

    im_rd_0_1 = np.zeros((8,8,256,256),dtype=np.complex128)
    for j1 in range(J):
        for l1 in range(L):
            im_rd_0_1[j1,l1] = np.fft.ifft2(im_fd_0_1[j1,l1])

    for j1 in range(J):
        for l1 in range(L):
            for l2 in range(L):
                for j2 in range(J):
                    S20[j1,j2,l1,l2] = np.sum(
                        dual_conj(
                            im_rd_0_1[j1,l1],
                            im_rd_0_1[j2,l2]
                        )
                    )

    out_coeff.append(S20)

    ##Coeff at (1,2) level
    if coeff_12_on == 1:
        if np.logical_and(rotate_on == 0,phase_on == 0):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    not_rot_im = np.fft.fftshift(im_fd_0_1[j1,l1,:,:])
                    for l2 in range(L):
                        for j2 in range(J):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    not_rot_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 1,phase_on == 0):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    for l2 in range(L):
                        rot_im = complex_rotate_image(np.fft.fftshift(im_fd_0_1[j1,l1,:,:]),(l1-l2)*180/L)
                        for j2 in range(J):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    rot_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 0,phase_on == 1):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for j2 in range(J):
                    for l1 in range(L):
                        stretch_im = np.fft.fftshift(phase_stretch(im_rd_0_1[j1,l1,:,:],2**(j1-j2)))
                        for l2 in range(L):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    stretch_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 1,phase_on == 1):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    for l2 in range(L):
                        rot_im = complex_rotate_image(np.fft.fftshift(im_fd_0_1[j1,l1,:,:]),(l1-l2)*180/L)
                        for j2 in range(J):
                            stretch_im = np.fft.fftshift(phase_stretch_f(np.fft.fftshift(rot_im),2**(j1-j2)))
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    stretch_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

    return out_coeff

<filename>scratch_AKS/DHC.py<gh_stars>1-10
def self_conj(a):
    return np.real(np.multiply(a,np.conj(a)))

def complex_rotate_image(image, angle):
    image_center = tuple(np.array(image.shape[1::-1]) / 2)
    rot_mat = cv2.getRotationMatrix2D(image_center, angle, 1.0)
    result_real = cv2.warpAffine(np.real(image), rot_mat, image.shape[1::-1], flags=cv2.INTER_LANCZOS4)
    result_imag = cv2.warpAffine(np.imag(image), rot_mat, image.shape[1::-1], flags=cv2.INTER_LANCZOS4)
    return result_real+result_imag*1j

def phase_stretch(im_rd_0_1,stretch,p0=0):
    im_stretch = np.abs(im_rd_0_1)*np.exp(stretch*np.angle(im_rd_0_1)*1j+np.pi*p0*1j)
    return np.fft.fft2(im_stretch)

def phase_stretch_f(U_1_c_f,stretch,p0=0):
    im_rs = np.fft.ifft2(U_1_c_f)
    im_stretch = np.abs(im_rs)*np.exp(stretch*np.angle(im_rs)*1j+np.pi*p0*1j)
    return np.fft.fft2(im_stretch)

def dual_conj(a,b):
    return np.abs(np.multiply(a,np.conj(b)))

def DHC(image,filter_set,phase_on = 1, rotate_on = 1, coeff_12_on = 1, norm_on = 1):
    Nx, Ny = image.shape

    out_coeff = []

    ##Coeff at (1,0) level
    #takes care of image normalization
    S0 = np.zeros((2,1),dtype=np.float64)
    if norm_on == 1:
        S0[0] = np.mean(image)
        norm_im = image-S0[0]
        S0[1] = np.sum(np.square(norm_im))/(Nx*Ny)
        norm_im = norm_im/np.sqrt(Nx*Ny*S0[1])
    else:
        norm_im = image

    out_coeff.append(S0)

    ##Coeff at (1,1) level
    #store the fft of images and im*filters for later
    im_fd_0 = np.fft.fft2(norm_im)

    S1 = np.zeros((8,8),dtype=np.float64)
    im_fd_0_1 = np.zeros((8,8,256,256),dtype=np.complex128)
    for j in range(J):
        for l in range(L):
            im_fd_0_1[j,l] = np.multiply(im_fd_0,filter_set[j,l]) #wavelet already in fft domain not shifted
            S1[j,l]=np.sum(self_conj(im_fd_0_1[j,l])) #normalization choice arb to make order unity

    out_coeff.append(S1)

    ##Coeff at (2,0) level
    S20 = np.zeros((J,J,L,L),dtype=np.float64)

    im_rd_0_1 = np.zeros((8,8,256,256),dtype=np.complex128)
    for j1 in range(J):
        for l1 in range(L):
            im_rd_0_1[j1,l1] = np.fft.ifft2(im_fd_0_1[j1,l1])

    for j1 in range(J):
        for l1 in range(L):
            for l2 in range(L):
                for j2 in range(J):
                    S20[j1,j2,l1,l2] = np.sum(
                        dual_conj(
                            im_rd_0_1[j1,l1],
                            im_rd_0_1[j2,l2]
                        )
                    )

    out_coeff.append(S20)

    ##Coeff at (1,2) level
    if coeff_12_on == 1:
        if np.logical_and(rotate_on == 0,phase_on == 0):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    not_rot_im = np.fft.fftshift(im_fd_0_1[j1,l1,:,:])
                    for l2 in range(L):
                        for j2 in range(J):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    not_rot_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 1,phase_on == 0):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    for l2 in range(L):
                        rot_im = complex_rotate_image(np.fft.fftshift(im_fd_0_1[j1,l1,:,:]),(l1-l2)*180/L)
                        for j2 in range(J):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    rot_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 0,phase_on == 1):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for j2 in range(J):
                    for l1 in range(L):
                        stretch_im = np.fft.fftshift(phase_stretch(im_rd_0_1[j1,l1,:,:],2**(j1-j2)))
                        for l2 in range(L):
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    stretch_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

        if np.logical_and(rotate_on == 1,phase_on == 1):
            S12 = np.zeros((J,J,L,L),dtype=np.float64)
            for j1 in range(J):
                for l1 in range(L):
                    for l2 in range(L):
                        rot_im = complex_rotate_image(np.fft.fftshift(im_fd_0_1[j1,l1,:,:]),(l1-l2)*180/L)
                        for j2 in range(J):
                            stretch_im = np.fft.fftshift(phase_stretch_f(np.fft.fftshift(rot_im),2**(j1-j2)))
                            S12[j1,j2,l1,l2] = np.sum(
                                dual_conj(
                                    stretch_im,
                                    np.fft.fftshift(im_fd_0_1[j2,l2,:,:])
                                )
                            )

            out_coeff.append(S12)

    return out_coeff

##Coeff at (1,0) level #takes care of image normalization ##Coeff at (1,1) level #store the fft of images and im*filters for later #wavelet already in fft domain not shifted #normalization choice arb to make order unity ##Coeff at (2,0) level ##Coeff at (1,2) level