; KOREG.PRO : 
;
; 1996  by Andreas Reigber
;------------------------------------------------------------------
; Automatic coregistration for two complex SAR images with
; subpixel precision. Algorithm based on coherenece optimisation.
; Useful only for spaceborne images where the offset is spatilly
; invariant.
;------------------------------------------------------------------
; This software has been released under the terms of the GNU Public
; license. See http://www.gnu.org/copyleft/gpl.html for details.
;------------------------------------------------------------------


pro koreg,bild1,bild2
        s=size(bild1)
        xpos=intarr(4)
        ypos=intarr(4)
        xw=[ s(1)/4 , s(1)*3/4, s(1)/4 , s(1)*3/4]
        yw=[ s(2)/4 , s(1)/4, s(1)*3/4 , s(1)*3/4]  
        for i=0,3 do begin
                auxcor = fft( fft(abs(bild1(xw(i)-256:xw(i)+255,yw(i)-256:yw(i)+255)),-1) $
                       * conj( fft(abs(bild2(xw(i)-256:xw(i)+255,yw(i)-256:yw(i)+255)),-1)) , 1)       
                aux = max( abs( auxcor ), auxpos )      
                posx = auxpos mod 512
                posy = auxpos/512
                if (posx gt 256) then posx = -512 + posx 
                if (posy gt 256) then posy = -512 + posy 
                xpos(i)=posx
                ypos(i)=posy
        endfor 
        shg_x=median(xpos)
        shg_y=median(ypos)
        bild2=shift(bild2,shg_x,shg_y)
        xpos=fltarr(4)
        ypos=fltarr(4)
        kpos=fltarr(4)
        for k=0,3 do begin
                sbild1=bild1(xw(k)-32:xw(k)+31,yw(k)-32:yw(k)+31)
                sbild2=bild2(xw(k)-32:xw(k)+31,yw(k)-32:yw(k)+31)

                koh=fltarr(17,17)
                fbild = fft(sbild2,-1)

                freq=shift((findgen(64)-32),32)
                zero_a = fltarr(64)+1

                fbild2=complexarr(64,64)
                for i=0,16 do begin
                        phase_ran = (i-8)/10.0  *2*!pi/64
                        for j=0,16 do begin
                                phase_azi = (j-8)/10.0 *2*!pi/64      
                                shift_matrix = ((phase_ran*freq)#zero_a)+(zero_a#(phase_azi*freq))
                                fbild2 = fbild*exp(complex(0,shift_matrix))
                                fbild2=fft(fbild2,1)
                                koh(i,j) = total(abs( smooth(sbild1*conj(fbild2),7,/edge_truncate) /  $
                                sqrt( smooth(abs(sbild1)^2,7,/edge_truncate) * smooth(abs(fbild2)^2,7,/edge_truncate) ) ))  
                        endfor
                endfor
                aux=max(koh,auxpos)
                posx=(auxpos mod 17)
                posy=(auxpos  /  17)
                xpos(k)=(posx-8.0)/10.0
                ypos(k)=(posy-8.0)/10.0
                kpos(k)=aux/64/64
        endfor 

        shf_x = total(kpos*xpos)/total(kpos)
        shf_y = total(kpos*ypos)/total(kpos)
        bild2 = fft(bild2,-1)
        freq1=shift((findgen(s(1))-s(1)/2),s(1)/2)
        freq2=shift((findgen(s(2))-s(2)/2),s(2)/2)
        zero_a1=fltarr(s(2))+1
        zero_a2=fltarr(s(1))+1
        phase_ran = shf_x  *2*!pi/s(1)
        phase_azi = shf_y  *2*!pi/s(2)   
        shift_matrix = ((phase_ran*freq1)# zero_a1)+(zero_a2#(phase_azi*freq2))
        bild2 = bild2*exp(complex(0,shift_matrix))
        bild2=fft(bild2,1)	

        print,'X-Offset :',shg_x+shf_x 
        print,'Y-Offset :',shg_y+shf_y 
        return
end