;+ ;NAME: ; COREGISTRATE ;PURPOSE: ; Coregister an image to a reference image to sub-pixel precision, allowing ; arbitrary image warping. The image is translated, rotated, ; sheared and otherwise distorted by iteratively adjusting a list of ; control points to maximize is correspondence with the reference. The key ; question in development has been how to measure that correspondence, as ; the cross-correlation seems not to yield reliable results. ;CALLING SEQUENCE: ; coregistrate, reference, image, N, refx, refy $ ; [, imx] [, imy] [, /verbose] [, image2=image2] $ ; [, missing=missing] [control_pts=control_pts] $ ; [, Kx=Kx] [, Ky=Ky] ;INPUTS: ; reference --- reference to which image should be coregistered. Must ; have the same dimensions as image. You can NaN-pad one or both ; images if needed. ; image --- a 2D image. Must have the same dimensions as reference. You ; may NaN-pad one or both images if needed. ; N --- degree of polynomial to use for warping image. See POLY_2D. To ; prevent degeneracy of the result, it is required that the number of ; control points equal (N+1)^2. ; refx --- (N+1)^2-element array of x-coordinates of control points in ; reference. ; refy --- (N+1)^2-element array of y-coordinates of control points in ; reference . ;OPTIONAL INPUTS: ; imx --- (N+1)^2-element array of x-coordinates of control points in image, ; as an initial guess for alignment to the reference. If not supplied, ; then imx defaults to refx. ; imy --- (N+1)^2-element array of y-coordinates of control points in image, ; as an initial guess for alignment to the reference. If not supplied, ; then imy defaults to refy. ;OPTIONAL KEYOWRD INPUTS: ; verbose --- if set, then diagnostic messages are printed to the screen. ; missing --- value to use for missing data on input reference, image, and ; warped image. Default = 0.0. ;OPTIONAL KEYWORD OUTPUTS: ; image2 --- coregistered version of image. ; control_pts --- structure containing final alignment of control points ; in image coordinates as control_pts.imx and control_pts.imy. Also ; contains .refx and .refy fields for completeness of description. ; Kx, Ky --- final coefficient matrices from POLYWARP for POLY_2D. ;COMMON BLOCKS: ; COREGISTRATE_CB --- used to convey context to the function of merit. ;BUGS (OR UNINTENDED FEATURES): ; 2009-06-03 CCK: Kathryn Williamson and I realized that IMAGE and ; REFERENCE must be of the same size. They can be padded (preferably ; with a bad data marker like NaN) to accomplish this. ;MODIFICATION HISTORY: ; 2006-05-09 C. Kankelborg. ; 2006-05-10 C. Kankelborg. Fixed CR_badness. A bug in POLY_2D produces ; anomalous results at the edges of the warped image, where it ; is extrapolating and really ought to be marked as missing data. ; 2006-05-10 C. Kankelborg switched CR_badness from cross-correlation to ; chi-squared, because cross-correlation is not, in general, maximized ; at perfect alignment when we allow image stretching. ; 2006-05-12 C. Kankelbog. Added input filtering on image and reference. ; This is necessary because warping the images smooths them, suppressing ; the noise. If the filtering is not done ahead of time, then the noise ; level will depend upon the phase of the image offset. This in turn ; causes a bias in the coregistration (it prefers 1/2-pixel offset, since ; that minimizes the high frequency noise). ; I have now incorporated both chi squared and cross-correlation ; badness routines within coregistrate.pro. Switch between them by ; editing function CR_badness, which is now just a front end. ; 2006-05-19 C. Kankelborg. Modified behavior of the missing keyword so ; that missing data in both the image and reference inputs is handled ; properly. Also made the verbose reporting a little better. ; 2006-05-23 C. Kankelborg. Fixed a bug that caused COREGISTRATE to crash ; when run in verbose mode on the Z-buffer. Now diagnostic displays ; operate only under X windows. ; 2006-05-24 C. Kankelborg. tightened the convergence criterion for amoeba ; from 1e-4 to 1e-5, and now allow 2000 instead of 1000 calls. ; 2006-05-31 C. Kankelborg. Applied the POLY_2D workaround (2006-05-10) to ; the output keyword IMAGE2. ; 2006-05-31 C. Kankelborg. Discovered that the smoothing that was ; introduced to beat down the noise causes another serious problem. ; If the filtering is allowed to soften the edge between the ; missing data and the good pixels, then two bad things happen: ; (1) formerly bad pixels along that edge are no longer recognized ; as missing and (2) good pixels along that edge have their contents ; influenced by proximity to bad pixels. ; 2006-06-01 C. Kankelborg. Implemented solution to yesterday's problem. ; Used new routine, CK_CONVOL. ; 2006-06-03 C. Kankelborg. The previous POLY_2D workaround was not ; general enough. New approach is to convert missing data to NaN, ; mark image border with NaN before running POLY_2D. Convert NaN ; values back to whatever is the desired 'missing' value on program ; return. Also removed CHI_EDGE_BADNESS, which was going nowhere. ; Fixed crash when running in batch mode with Z-buffer. ; 2006-06-09 C. Kankelborg. Added IDL_CORR_BADNESS, using ; a cross-correlation as defined for IDL's C_CORRELATE. ; This performs very badly. ;- ;************************************************** ;* CR_badness is the function to be minimized. * ;* It is really just a front end for several * ;* possible functions of merit, such as: * ;* chisq_badness --- chi squared statistic * ;* corr_badness --- (minus) cross-correlation * ;* * ;************************************************** function CR_badness, prams, warped=warped, terminate=terminate, ss=ss forward_function chisq_badness, corr_badness, idl_corr_badness return, chisq_badness(prams, warped=warped, terminate=terminate, ss=ss) end ;************************************************** ;* chisq_badness is one possible merit function. * ;* It is a chi-squared statistic. * ;* The keywords facilitate debugging: * ;* WARPED --- image warped according to the * ;* control points Crefx, Crefy, imx, imy. * ;* TERMINATE --- if set, exit for debugging. * ;************************************************** function chisq_badness, prams, warped=warped, terminate=terminate, ss=ss common coregistrate_cb, Creference, Cimage, CN, Crefx, Crefy, $ Cverbose, weights ;extract image coordinates of control points from input prams m = n_elements(prams) if m ne 2*(CN+1L)^2 then message,'Wrong number of parameters?!' imx = prams[0:m/2-1] imy = prams[m/2:m-1] ;Calculate warped image polywarp, imx, imy, Crefx, Crefy, CN, Kx, Ky, /double warped = poly_2d(Cimage, Kx, Ky, 2, cubic=-0.5, missing=NaN) ;Create renormalized reference image ss = where(finite(warped) and finite(Creference)) ref = Creference * total(warped[ss])/total(Creference[ss]) ;Calculate chi squared chisq = total( (warped[ss] - ref[ss])^2 /(warped[ss]>1.0) ) chisq_reduced = chisq/n_elements(ss) if (Cverbose) then begin if (!D.NAME eq 'X') then tv, displayscale(warped-ref, missing=NaN) print,'CHISQ_BADNESS: chisq = ',chisq,' chisq_reduced = ',chisq_reduced print,'imx = ',imx print,'imy = ',imy endif if keyword_set(terminate) then $ message,'Terminated with extreme prejudice. Time to debug!' return, chisq_reduced end ;*************************************************** ;* corr_badness is one possible merit function. * ;* It is one minus the cross-correlation. * ;* The keywords facilitate debugging: * ;* WARPED --- image warped according to the * ;* control points Crefx, Crefy, imx, imy. * ;* TERMINATE --- if set, exit for debugging. * ;*************************************************** function corr_badness, prams, warped=warped, terminate=terminate, ss=ss common coregistrate_cb, Creference, Cimage, CN, Crefx, Crefy, $ Cverbose, weights ;extract image coordinates of control points from input prams m = n_elements(prams) if m ne 2*(CN+1L)^2 then message,'Wrong number of parameters?!' imx = prams[0:m/2-1] imy = prams[m/2:m-1] ;Calculate warped image polywarp, imx, imy, Crefx, Crefy, CN, Kx, Ky, /double warped = poly_2d(Cimage, Kx, Ky, 2, cubic=-0.5, missing=NaN) ;Calculate cross-covariance ss = where(finite(warped) and finite(Creference)) covariance = total( ( warped[ss] - mean(warped[ss]) ) $ * ( Creference[ss] - mean(Creference[ss]) ) ) correlation = covariance / sqrt( total( ( warped[ss] - mean(warped[ss]) )^2 ) $ *total( ( Creference[ss] - mean(Creference[ss]) )^2 ) ) if (Cverbose) then begin ;Create renormalized reference image ref = Creference * total(warped[ss])/total(Creference[ss]) if (!D.NAME eq 'X') then tv, displayscale(warped-ref, missing=NaN) print,'CORR_BADNESS: covariance = ',covariance,' correlation = ',correlation print,'imx = ',imx print,'imy = ',imy endif if keyword_set(terminate) then $ message,'Terminated with extreme prejudice. Time to debug!' return, 1-correlation end ;*************************************************** ;* The following merit function is provided as a * ;* demonstration only. It is not recommended! * ;*************************************************** ;* idl_corr_badness is a possible merit function * ;* It is one minus the cross-correlation like * ;* corr_badness, but correlation is calculated in * ;* the manner of IDL's C_CORRELATE. Bad idea! * ;* The keywords facilitate debugging: * ;* WARPED --- image warped according to the * ;* control points Crefx, Crefy, imx, imy. * ;* TERMINATE --- if set, exit for debugging. * ;*************************************************** function IDL_corr_badness, prams, warped=warped, terminate=terminate, ss=ss common coregistrate_cb, Creference, Cimage, CN, Crefx, Crefy, $ Cverbose, weights ;extract image coordinates of control points from input prams m = n_elements(prams) if m ne 2*(CN+1L)^2 then message,'Wrong number of parameters?!' imx = prams[0:m/2-1] imy = prams[m/2:m-1] ;Calculate warped image polywarp, imx, imy, Crefx, Crefy, CN, Kx, Ky, /double warped = poly_2d(Cimage, Kx, Ky, 2, cubic=-0.5, missing=NaN) ;Calculate image and reference means and normalization factors based ;on the mathematical definitions given for IDL C_CORRELATE: ss = where(finite(Cimage)) mean_image = mean(Cimage[ss]) norm_image = sqrt( total( (Cimage[ss] - mean_image)^2 ) ) ss = where(finite(Creference)) mean_reference = mean(Creference[ss]) norm_reference = sqrt( total( (Creference[ss] - mean_reference)^2 ) ) ;Calculate cross-covariance ss = where(finite(warped) and finite(Creference)) covariance = total( ( warped[ss] - mean_image ) $ * ( Creference[ss] - mean_reference ) ) correlation = covariance / (norm_image * norm_reference) if (Cverbose) then begin ;Create renormalized reference image ref = Creference * total(warped[ss])/total(Creference[ss]) if (!D.NAME eq 'X') then tv, displayscale(warped-ref, missing=NaN) print,'IDL_CORR_BADNESS: covariance = ',covariance,' correlation = ',correlation print,'imx = ',imx print,'imy = ',imy endif if keyword_set(terminate) then $ message,'Terminated with extreme prejudice. Time to debug!' return, 1-correlation end ;************************************ ;* M A I N P R O G R A M * ;************************************ pro coregistrate, reference, image, N, refx, refy, imx, imy, $ verbose=verbose, missing=missing, control_pts=control_pts, $ Kx=Kx, Ky=Ky, image2=image2 NaN = !values.f_nan ;IEEE NaN ;What shall we use to signify missing data? Default = NaN. if n_elements(missing) ne 1 then missing = NaN ;Check image and reference control points if n_elements(imx) eq 0 then imx = refx if n_elements(imy) eq 0 then imy = refy if n_elements(refx) ne (N+1L)^2 then message,'REFX must have (N+1)^2 elements.' if n_elements(refy) ne (N+1L)^2 then message,'REFY must have (N+1)^2 elements.' if n_elements(imx) ne (N+1L)^2 then message, 'IMX must have (N+1)^2 elements.' if n_elements(imy) ne (N+1L)^2 then message, 'IMY must have (N+1)^2 elements.' isize = size(image) Nx = isize[1] Ny = isize[2] ;Create common block shared with the badness function common coregistrate_cb, Creference, Cimage, CN, Crefx, Crefy, $ Cverbose, weights ;Make filtered versions of image and reference to beat down the high frequency ;noise. Populate the common block. ;N.B.: Treat missing data conservatively. See comments 2006-May-31 thru Jun-1. kernel=[[1,2,1],[2,4,2],[1,2,1]]/16.0d ;make kernel double precision so that ck_convol is all double precision. Cimage = ck_convol(double(image), kernel, missing=missing, $ method="conserve_taint",/edge_truncate) Creference = ck_convol(double(reference), kernel, missing=missing, $ method="conserve_taint",/edge_truncate) ;(Use all double precision variables for computing badness.) CN = N Crefx = refx Crefy = refy weights = 0 ;Reset the weights array for this run of coregistrate. if keyword_set(verbose) then Cverbose=1.0 else Cverbose=0.0 ;Convert all 'missing' data to NaN for compatibility with POLY_2D. ss_im = where( (Cimage eq missing) or (~ finite(Cimage)) ) ;list of bad image pixels. Note that ~ is the right kind of 'not'. ss_ref= where( (Creference eq missing) or (~ finite(Creference)) ) ;list of bad ref pixels. if ( ss_im[0] ne -1) then Cimage[ss_im] = NaN if (ss_ref[0] ne -1) then Creference[ss_ref] = NaN ;Mark data around the edges of Cimage missing, to prevent POLYWARP from ;extrapolating at the edges. Cimage[0,*] = NaN Cimage[Nx-1,*] = NaN Cimage[*,0] = NaN Cimage[*,Ny-1] = NaN if keyword_set(verbose and !D.NAME eq 'X') then begin print,"Here are the two images at the starting point." disp_im = displayscale(Cimage, missing=NaN, /ct) disp_ref = displayscale(Creference, missing=NaN) for i = 1,5 do begin tv, disp_im wait,1 tv, disp_ref wait,1 endfor endif ;Time to burn some CPU cycles. if keyword_set(verbose) then $ message,'Starting alignment at '+systime(), /informational prams = amoeba( 1e-5, function_name='CR_badness', ncalls=ncalls, $ nmax = 2000, P0=[imx,imy], scale=replicate(4.0, 2*(N+1L)^2) ) if keyword_set(verbose) then begin message,'Alignment completed at '+systime(), /informational print,'ncalls = ',ncalls endif if prams[0] eq -1 then message,'Amoeba failed to converge.' ;Dish up the results, including coregistered image m = n_elements(imx) imx = prams[0:m-1] imy = prams[m:2*m-1] control_pts = {imx:imx, imy:imy, refx:refx, refy:refy} polywarp, imx, imy, refx, refy, N, Kx, Ky, /double image2 = poly_2d(image, Kx, Ky, 2, cubic=-0.5, missing=NaN) ;Convert from NaN back to whatever the user specified for missing data. ss = where( ~ finite(image2) ) if (ss[0] ne -1) then image2[ss] = missing end