FUNCTION Rotscalefind, pic1, pic2 ;+ ; NAME: ; ; PURPOSE: ; ; CALLING SEQUENCE: ; ; INPUTS: ; ; OPTIONAL PARAMETERS: ; ; KEYWORDS: ; ; OUTPUTS: ; ; COMMON BLOCKS: ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; PROCEDURE: ; ; MODIFICATION HISTORY: ; - -200 P.Suetterlin, SIU ; 15-Oct-2004 Change algo to find lower left shift ; increase angle range to get enough overlap for ; lager rotation values ; 09-Oct-2005 more efficient way to compute x/y matrix ;- s1 = size(pic1, /dim) s2 = size(pic2, /dim) ;;; Find the shift in the lower left corner. With rotation, shifts ;;; might be severe, so starting with a small window gives bad ;;; results. Start with a 256 window and narrow down to 64 sh = [[0, 0], [-shc(pic1(0:255, 0:255), pic2(0:255, 0:255), /filt)]] sh(1, *) = sh(1, *)-min(sh(1, *)) sh(0, *) = sh(0, *)-min(sh(0, *)) sh(*, 1)=sh(*, 1)-shc(pic1(sh(0, 0):sh(0, 0)+127, sh(1, 0):sh(1, 0)+127), $ pic2(sh(0, 1):sh(0, 1)+127, sh(1, 1):sh(1, 1)+127), /filt) sh(1, *) = sh(1, *)-min(sh(1, *)) sh(0, *) = sh(0, *)-min(sh(0, *)) sh(*, 1)=sh(*, 1)-shc(pic1(sh(0, 0):sh(0, 0)+63, sh(1, 0):sh(1, 0)+63), $ pic2(sh(0, 1):sh(0, 1)+63, sh(1, 1):sh(1, 1)+63), /filt) sh(1, *) = sh(1, *)-min(sh(1, *)) sh(0, *) = sh(0, *)-min(sh(0, *)) ;;; remaining image size newx = (s1(0)-sh(0, 0)) < (s2(0)-sh(0, 1)) newy = (s1(1)-sh(1, 0)) < (s2(1)-sh(1, 1)) ;;; set up rotation index: ~40 degrees in steps of 1/20 phi = findgen(768)/20*!dtor ;;; set up (log) distance index. mx = newx-32-1 rm = fix(1000*alog(mx))/1000. ;step = alog(mx)-alog(mx-1) ;;; wanted acuracy: 1 pixel from scaling over ~1000 pixels = 1e-4 step = 1e-4 nr = 1024 r = findgen(nr)*step r = exp(r-r(nr-1)+rm) x = fltarr(nr, 768) y = fltarr(nr, 768) FOR j=0, 767 DO BEGIN cp = cos(phi(j)) sp = sin(phi(j)) FOR i=0, nr-1 DO BEGIN x(i, j) = r(i)*cp y(i, j) = r(i)*sp ENDFOR ENDFOR ;FOR i=0, nr-1 DO FOR j=0, 767 DO x(i, j) = r(i)*cos(phi(j)) ;FOR i=0, nr-1 DO FOR j=0, 767 DO y(i, j) = r(i)*sin(phi(j)) p1 = interpolate(pic1, x+sh(0, 0)+32, y+sh(1, 0)+32) p2 = interpolate(pic2, x+sh(0, 1)+32, y+sh(1, 1)+32) sh2 = shc(p1, p2, /int, /filt) return, [-sh2(1)/20, exp(sh2(0)*step)] END