function gaincalib, images, x, y, object=object, maxiter=maxiter, additive=additive ;+ ; NAME: GAINCALIB ; PURPOSE: Produce a gain table from a set of images with relative offsets ; CALIING SEQUENCE: ; Result = gaincalib(images, object=object, additive=additive) ; INPUT: ; images a three-dimensional array representing ; a sequence of two-dimensional images ; images(*,*,k) ( k=0, 1,.., N-1). ; It is strongly recommended that the number ; of elements in each image is a power of 2 ; for the efficent FFT processing. ; x an array of x-shift ; y an array of y-shift ; ; OUTPUT: ; Result the gain table if the keyword ADDITIVE is not set ; or the offset table if the keyword is set. ; INPUT KEYWORDS: ; maxiter maximum # of iternation (default=5) ; OUTPUT KEYWORD: ; object flat-field corrected object ; History: ; 1999 May, Jongchul Chae ;- if n_elements(maxiter) eq 0 then maxiter=5 s=size(images) i = indgen(s(1))#replicate(1, s(2)) j = replicate(1, s(1))#indgen(s(2)) x0 = x - median(x) y0 = y - median(y) ; Initial Estimates if keyword_set(additive) then Logimages = images $ else Logimages = alog(images>(0.01*median(images(*,*,0)))) Flat = total(Logimages, 3)/s(3) & Flat = Flat-median(Flat) C = fltarr(s(3)) for k=0, s(3)-1 do begin C(k) = median(logimages(*,*,k)-flat) endfor C = C -C(0) ; Start Iteration Object = 0. for iter=0, maxiter do begin aa=0.0 & bb=0.0 for k=0, s(3)-1 do begin weight = (i ge -x0(k)) and (i le -x0(k)+s(1)-1) $ and (j ge -y0(k)) and (j le -y0(k)+s(2)-1) aa = aa + (-C(k) -Object $ +shift_sub(Logimages(*,*,k)-Flat, -x0(k), -y0(k)) )*weight bb = bb+weight endfor DelObject = aa/(bb+1.e-5) Object = Object + DelObject aa=0. & bb=0.0 for k=0, s(3)-1 do begin weight = (i ge (x0(k)) ) and (i le x0(k)+s(1)-1) $ and (j ge (y0(k))) and (j le y0(k)+s(2)-1) ob = (-C(k)-shift_sub(Object, x0(k), y0(k))-Flat+ $ Logimages(*,*,k))*weight DelC = total(ob)/total(weight) C(k) = C(k)+DelC aa = aa + ob bb = bb+weight endfor DelFlat = aa/(bb+1.e-5) Flat = Flat+DelFlat error = max(abs(Delflat)) print, 'iteration # =', iter , ' error=', error if error le 1.0e-6 then goto, final endfor final: if not keyword_set(additive) then begin Object = exp(Object) Flat = exp(Flat) C = exp(C) endif return, flat end