;+ ;NAME: ; test4.pro ;PURPOSE: ; Test script for coregistrate.pro, with graphical and numerical ; diagnostic outputs. ;CALLING SEQUENCE: ; .run test4 ;DEPENDENCIES: ; Requires DSCN3316.JPG or test4QuickStart.sav as a source of test ; image data. If you want a fresh randomized example, delete or ; rename the .sav file. ;KNOWN ISSUES: ; The only way currently to change between having noise and not ; having noise is to edit the source code (NoiseFlag =...). ; The choice of merit function used for coregistration is also ; hardwired in COREGISTRATE. ;MODIFICATION HISTORY: ;2006-Jun-11 CCK. Based on test3. In test3, the test image was ; warped, and COREGISTRATE was used to unwarp it. But it turns ; out that the form of the nonlinear warping is such that its ; inverse usually cannot be expressed in the same form. Now in ; test4, the initial warping is done to form the reference ; image, and COREGISTRATE will seek to replicate the warp. So ; it should be possible to attain an exact solution, even if ; the warping has a strong nonlinear term. ;2008-Jul-01 CCK. Reintroduced random shift of the reference image, ; leaving the old fixed offsets (which were designed to be ; a difficult case) in the comments for future reference. Improved ; comments throughout test4. ;- missing = 0.0 ;same convention as MOSES_PREP NaN = !values.f_nan forward_function CR_badness if file_test("test4QuickStart.sav") then begin print,"Restoring test4QuickStart.sav" restore, "test4QuickStart.sav" print,"Skipping image setup calculations, using canned example..." goto, QUICKSTART endif else begin print,"Calculating image and reference..." endelse ;This high resolution JPEG is used to create all the test data. read_jpeg, 'DSCN3316.JPG', rawimage, /grayscale isize = size(rawimage) NxRaw = isize[1] NyRaw = isize[2] ;We will downsample rawimage by integer factor mag. The idea is to ;do all the initial image processing at super high resolution; this ;represents the action of two "instruments" producing the "image" and ;the "reference" both from the same ideal "object". These two instruments ;have different sensitivity, different PSFs, and different optical distortion. mag = 4 Nx = NxRaw/mag Ny = NyRaw/mag ;Now apply PSF and downsample to create the image. foo = double(rawimage) kernel = [[0,1,1],[1,2,1],[1,1,0]]/8.0d print,"Applying image PSF..." smoothimage = double(convolf(foo, kernel, reps=160.0, /complex)) image = 15.0*rebin(smoothimage, Nx, Ny) ;And now apply a DIFFERENT PSF, distort and downsample to create the reference. foo = double(rawimage) kernel = [[1,1,0],[1,2,1],[0,1,1]]/8.0d print,"Applying different PSF to reference..." smoothimage = double(convolf(foo, kernel, reps=160.0, /complex)) ;Warp the reference, and make it interesting! ;Choose control points at image corners: imxRaw = [0.0,NxRaw,0,NxRaw] imyRaw = [0.0,0,NyRaw,NyRaw] ;Random x-y offsets of the 4 control points: dxRaw = 10*mag*randomn(seed, 4) ;formerly used fixed values: [6.0, 7, 7, 6] dyRaw = 10*mag*randomn(seed, 4) ;formerly used fixed values: [-2.0,2, 2,-2] refxRaw = imxRaw + dxRaw refyRaw = imyRaw + dyRaw polywarp, imxRaw, imyRaw, refxRaw, refyRaw, 1, KxRaw, KyRaw, /double ;Fill border of image with NaN to work around POLY_2D missing data bug smoothimage[0, *] = NaN smoothimage[NxRaw-1, *] = NaN smoothimage[*, 0] = NaN smoothimage[*, NyRaw-1] = NaN warpedimage = poly_2d(smoothimage, KxRaw, KyRaw, 2, cubic=-0.5, missing=NaN) ;Shrink to size and scale reference = 30.0*rebin(warpedimage, Nx, Ny) ;Figure out where the control points went with rebinning ;CCK 2006-Jun-10 see "Coordinates in Rebinned Images", rebin_coords.pdf refxBinned = (refxRaw / mag) + (1.0/(2.0*mag)) - 0.5 imxBinned = (imxRaw / mag) + (1.0/(2.0*mag)) - 0.5 refyBinned = (refyRaw / mag) + (1.0/(2.0*mag)) - 0.5 imyBinned = (imyRaw / mag) + (1.0/(2.0*mag)) - 0.5 print,"Saving image data to test4QuickStart.sav to speed up future test4 runs." ;but first, delete some big stuff that I don't need anymore. smoothimage = 0 warpedimage = 0 saturated = 0 rawimage = 0 foo = 0 save, file="test4QuickStart.sav" QUICKSTART: print,"Now just add noise an mark saturated pixels..." ;Restore seed so that the random number generator can pick ;up where it left off last time. if file_test("seed.sav") then restore, "seed.sav" ;****************************** ;* SOME ADJUSTABLE PARAMETERS * ;****************************** percentile = 0.95 ;fraction of 'good' (unsaturated) pixels in image. NoiseFlag = 1 ;1 to add Poisson noise; 0 to leave it out. ;Add Poisson noise if NoiseFlag then begin print,"Applying Poisson noise." pnoise, reference, seed=seed pnoise, image, seed=seed endif else begin print,"Noise is disabled. Using perfect fake data." endelse ;Mark 'saturated' or NaN data missing in image saturated = image( (sort(image))[percentile*Nx*Ny] ) ;saturated pixels. ss = where( image ge saturated or (~ finite(image)) ) image[ss] = missing ;mark as missing ;Mark missing data in reference ss = where(~ finite(reference)) reference[ss] = missing ;mark as missing ;Display image and reference. window, 1, title='reference', xsize=Nx, ysize=Ny tv,displayscale(reference, missing=missing, /ct) window, 2, title='image', xsize=Nx, ysize=Ny tv, displayscale(image, missing=missing, /ct) window,3, title='residuals', xsize=Nx, ysize=Ny ;for use by COREGISTRATE. ;Parameters for coregistrate: refx = refxBinned refy = refyBinned ;Produce an initial guess that's in the ballpark, but not too close!w imx = imxBinned + 10*randomn(seed,4) imy = imyBinned + 10*randomn(seed,4) ;****************** C O R E G I S T R A T I O N ************************* coregistrate, reference, image, 1, refx, refy, imx, imy, /verbose, $ control_pts=control_pts, Kx=Kx, Ky=Ky, image2=image2, missing=missing ;************************************************************************ ;Explore the result! print,"Here is where COREGISTRATE thinks the optimum alignment is:" fu = CR_badness([control_pts.imx,control_pts.imy],warped=coregistrated,$ ss=ss_cor) print,"Here is where I think the optimum alignment should be:" imx_expect = imxBinned imy_expect = imyBinned window,5,xsize=Nx, ysize=Ny, title='residuals for expected alignment' fu = CR_badness([imx_expect,imy_expect],warped=expected, ss=ss_exp) ;Numerical comparison x_err = control_pts.imx - imx_expect y_err = control_pts.imy - imy_expect rms_xerr = sqrt(mean(x_err^2)) rms_yerr = sqrt(mean(y_err^2)) rms_err = sqrt( mean(x_err^2 + y_err^2) ) print,"RMS error in x: ",rms_xerr,"; in y: ",rms_yerr,"; in r: ",rms_err ;Print comments regarding noise if NoiseFlag then begin print,"Test included Poisson noise." endif else begin print,"Noise was disabled for this run." endelse ;Report comments regarding bad pixels print,"The fraction of unsaturated pixels in image was: ",percentile ;Save seed so that next time, the random number generator can pick ;up where it left off. save, seed, file="seed.sav" end