pro es_img_sek, bs_img, es_img, img ; ; ; ; to run this program you need to do first: ; IDL> .r es_img_2phil ; IDL> es_img_sek, bs_img, es_img, img ; once you get img, make a simple program ; yourself to write img into a jpeg file ; as shown in es_img_sek4.jpg, adjust the ; contrast whatsoever to your taste - I just do not ; know where is my program!! and then you ; can use this jpeg file as input for cover.pro ; KEEP IN MIND that this image is from 00feb29!! ; ; get these two following files from the backup CDROM in my ; office (if not there, snag at Jeff!!) ; bs_img = readfits('l00feb29036.fts', /silent) es_img = readfits('l00feb29037.fts', /silent) help, bs_img help, es_img bs_img(0:20, *) = min(bs_img) bs_img(491:511, *) = min(bs_img) bs_img(*, 0:20) = min(bs_img) bs_img(*, 491:511) = min(bs_img) es_img(0:20, *) = min(es_img) es_img(491:511, *) = min(es_img) es_img(*, 0:20) = min(es_img) es_img(*, 491:511) = min(es_img) bs_cx = 276. bs_cy = 245. bs_r = 190.2128 bs_tilt = 2.549692 bs_lib_beta = 0.008751445 bs_lib_alpha = -.04280026 es_cx = 276. es_cy = 248. es_r = 189.9373 es_tilt = 2.538734 es_lib_beta = 0.08716932 es_lib_alpha = -0.04277911 bs_img = rotate(bs_img, 7) es_img = rotate(es_img, 7) bs_img = shift(bs_img, -bs_cx+256, -bs_cy+256) es_img = shift(es_img, -es_cx+256, -es_cy+256) ratio = es_r/bs_r es_img = congrid(es_img, 512/ratio, 512/ratio, /interp) ndim = n_elements(es_img(*, 0)) print, ndim, (ndim-512)/2 es_img = es_img((ndim-512)/2:(ndim-512)/2+511, (ndim-512)/2:(ndim-512)/2+511) bs_img = rot(bs_img, (es_tilt-bs_tilt)/!dtor/1., cubic = -0.5) bs_img = rotate(bs_img, 7) es_img = rotate(es_img, 7) img = bs_img + es_img*1.4 img = enhance_mem(img) tvscl, img<45000 ; ;img1 = filt_img(bs_img/1000., 0.01, /high) ;img2 = filt_img(es_img/1000., 0.01, /high) ;img = filt_img(img/1000., 0.01, /high);<15>(-20) ;tvscl, img<20>(-4) end ;*********************************************************************** function filt_img, img, fac, high=high, low=low, exp = exp ;using FFT to change the contrast of image. ;*********************************************************************** sz=size(img) nx=sz(1) & ny=sz(2) fimg=fft(img, 1) psf=dist(nx, ny) psf=psf/max(psf) if n_elements(fac) ne 0 then fac=float(fac) key=0 if keyword_set(high) then key=1 if keyword_set(low) then key=-1 if keyword_set(high) AND keyword_set(low) then key=0 case key of 1: filt=1.-1./(1.+(temporary(psf)/fac)) -1: filt=1./(1.+(temporary(psf)/fac)) 0: filt=1. endcase if keyword_set(exp) then filt = 1-exp(-psf*fac) nimg=fft(temporary(fimg)*temporary(filt), -1) nimg=float(nimg) ;hlp, nimg ;tvscl, nimg<1000>(-1000) return, nimg bail: end ;************************************************************** function enhance_mem, image, fwhm=fwhm, $ magnification=magnification, maxit=maxit, reg=reg if n_elements(magnification) ne 1 then magnification=1 if n_elements(maxit) eq 0 then maxit=8 if n_elements(reg) eq 0 then reg=0.05 s=size(image) x=(findgen(s(1))-s(1)/2)#replicate(1, s(2)) y=replicate(1, s(1))#(findgen(s(2))-s(2)/2) if n_elements(fwhm) eq 0 then fwhm=1.5 a=fwhm/1.55 k = shift(2*!pi*sqrt((x/s(1))^2+(y/s(2))^2), -s(1)/2, -s(2)/2) mtf = exp(-a*k) nx=2^fix(alog(s(1))/alog(2.)) nx=nx*(1+(s(1) gt nx)) ny=2^fix(alog(s(1))/alog(2.)) ny =ny*(1+(s(2) gt ny)) mem96, congrid(image, nx, ny, cubic=-0.5), x, $ otf=congrid(mtf, nx, ny, cubic=-0.5), $ /quiet, maxit=maxit, reg=reg return, congrid(x, s(1)*magnification, s(2)*magnification, $ cubic=-0.5) end ;*********************************************************************************** pro mem96, y, x, psf = psf, otf = otf, scale=a, bg=b, reg=reg, error = error, $ maxit=maxit, minit=minit, quiet=quiet ;+ ; ; NAME: ; mem96 ; PURPOSE: ; Deconvolve an image whose point spread function or ; optical transfer function is known by use of maximum ; entropy algorithm. ; CALLING SEQUENCE: ; mem96, y, x, psf = psf, otf = otf, scale=scale, bg=bg, $ ; reg=reg, error = error, maxit=maxit, minit=minit ; INPUT: ; y input image to be deconvolved (!-D or 2-D) ; ; KEWORD INPUT: ; psf point spread function ; otf optical transfer function or Fourier transform ; of point spread function. ; Either psf or otf should be given as input. ; Note that both psf and otf have their orgin of ; coordinate at lower left corner. ; scale scale coefficient a in Eq 3.17 ; of Chae's thesis (Default = 1.) ; bg background coefficient b in Eq 3.17 ; of Chae's thesis (Default = 0.) ; reg dimentionless regularizing parameter given ; by Eq 3.36 (Default=1.E-2) ; error error criterion for stopping iteration (Default=1.E-4) ; maxit maximum iteration number (Default=50) ; minit minimum iteration number (Default=5) ; quiet if set, no display of intermediate results ; ; OUTPUT: ; x deconvolved image ; REQUIRED ROUTINES: NONE ; REFERENCE : Chae's Ph.D. Thesis ; History: ; witten by Jong Chul Chae, March 1996 ; ;- ;=--------------------- KEYWORD PARAMETERS SETTING ------------------ if not keyword_set(maxit) then maxit=50 if not keyword_set(minit) then minit=5 if not keyword_set( error) then error=1.0e-3 if not keyword_set(reg) then reg=1.0e-2 if not keyword_set(a) then a = 1. if not keyword_set(b) then b = 0. if not keyword_set(quiet) then show=1 else show=0 if keyword_set(otf) then ftpsf=otf if keyword_set(psf) then ftpsf=fft(psf, 1) if show then begin print, ' ' print, ' ==============================================================' print, ' Welcome to Maximum Entropy Method Program for Deconvolution ' print, ' ==============================================================' endif ;=--------------------- INITIALIZATION --------------------------------- s=size(y) s1=s(1) s2=s(2) if show then begin window,11, xs=s1, ys=s2,r=2 window, 12, xs=s1, ys=s2, r=2 endif npixels = s1*s2 fty=fft(y, 1) tau2 = abs(ftpsf)^2 cftpsf =(temporary(ftpsf)) ftypsf = fty*temporary(cftpsf) m=a*y+b maxm=max(m) m=m>(maxm*0.01) x=(m-b)/a lambda00=total(1./m)/n_elements(m) ; Eq 3.34 reg2=reg/lambda00 ;------------------------------------------------------------------- if show then begin items = " ITER # BETA 1/LAMBDA REG1 ERROR " print, ' ' print, ' Now MEM iteration begins.' print, ' ' endif fx=fft(x, 1) iter=0 lambda=lambda00 ; ===================== ITERATION BLOCK BEGINNING========================== repeat begin lambda0=lambda reg1 = reg2*lambda q = fft(a*alog((a*x+b)/m), 1) ; Eq 3.23 and 3.27 dfx = -((fx*tau2-ftypsf)+reg2*q)/(tau2+reg1) ; Eq 3.27 and 3.33 dx = float(fft(dfx, -1)) beta1 = min((a*x+b)/((-a*dx)>0.00001*maxm)) beta = beta1*0.9 < 1. x0=x+dx*beta fx=temporary(fx)+dfx*beta dx=x0-x if iter le 100 then $ lambda= 0.5*(total((a*dx)^2/(a*x+b))/total(dx^2)+lambda0) ; Eq 3.31 x=temporary(x0) iter=iter+1 ;--------------------------------------------------------------------- error1=max((dx)/maxm, min=mm) if error1 lt abs(mm) then error1=mm converge = (abs(error1) le error) converge = converge and (beta ge 0.8) exceed = iter ge maxit ; -------------------------------------------------------------------- if show then begin wset, 12 tv, bytscl(dx, min=-0.005*maxm,max=0.005*maxm) wset, 11 tvscl, x if iter mod 5 eq 0 then begin print, '' print, items print, '' endif print, '' check = string(iter,format='(" ", i3)') +string(beta, format='(" ",g10.3)') $ +string(1./lambda, format='(" ",g10.3)')+string(reg1, format='(" ",e10.3)') $ +string(error1, format='(" " ,g10.3)') print, check if converge then print, ' Iteration has converged.' if exceed then print, ' Iteration # has excceeded maximum.' endif endrep until (iter gt minit) and (converge or exceed) ;========================== ITERATION BLOCK END ======================== if show then begin tvscl, x wait, 3 print, '' print, ' ============================================================ ' print, ' End of Maximum Entropy Method for Deconvolution ' print, ' ============================================================ ' endif return end