;+
;NAME: ivm_geomreflat
;
;PURPOSE: correct alternating-frame residual fringes in geometry-camera
;	   data w/o over-correcting photometry, so fringes not imparted later
;	   on data camera, but photometry changes are left for the deblur algorithm.
;
;CATEGORY: image processing
;
;CALLING SEQUENCE: 
;	ivm_geomreflat,geom,oblk
;
;INPUTS: geomin:  input geometry data = float(nx,ny, nimages)  [default input]
;	   oblk:  variable holding the IVM header structure, i.e. oblk
;        filein:  if set to 1, read input file ('geofin'); else, use array
;       fileout:  if set to 1, write output file ('geofout'); else, use array
;
;OUTPUTS: de-fringed geom+"_fl2" data.
;       geomout:  output de-fringed geom data (nx, ny, nimages)
;
;OPTIONAL KEYWORDS:
;	 geofin:  input geometry filename, flat-fielded. [use if filein]
;       geofout:  output filename if set, written if 'fileout' is true
;      output4D:  if set, reform the output array to (nx,ny, nf, nw*ns)
;    smoothsize: 
;
;SIDE EFFECTS:
;PROCEDURE: similar to ivm_desee in building reference image, then determine
;		solely the "dc offset" of tiles as destretchdeblur does, and only
; 		those (hopefully) in fringe pattern w/in some limit of a mean.
;
;MODIFICATION HISTORY: Written based on ivm_desee and ivm_destretchdeblur, to
;     correct fringes observed in *RAW* geometry-camera data which alternate
;     between even/odd frame number, and thus are not corrected for with the
;     default flat-field approach of the fringes from averaging multiple frames.
;     unknown how widespread these residual fringes are; difference images
;     show them clearly, and not in the data camera, from 2000 epoch data.
;
;     20091120 ELW - Modified to use array I/O, with optional file output
;     20100108 ELW - Floats propagated by 'ntiles=nn/16.' gone.  Assuming at
;			this stage that all nx or ny are divisible by 16.
;     20100526 ELW - file input option
;-

PRO ivm_geomreflat,geomin,geomout, oblk, filein, fileout, verbose=verbose, $
			geofin=geofin, geofout=geofout, $
			smoothsize=smoothsize, output4D=output4D

if (keyword_set(verbose)) then vb=verbose else vb=0
if (vb ge 1) then PRINT,'IVM_GEOMREFLAT:  ', SYSTIME()

nn = oblk.PPIX
nw = oblk.N_IMAGE_SETS
ns = oblk.NTIMES
nf = oblk.NFRAMES
nimages=nw*ns*nf

;nsmooth = 2 * nn/256 + 1 ;basic smoothing to counter pixellation noise.
; smoothing found to be bad in this application.  kept option below, and 
; as a keyword.
nsmooth = 1.0
if (keyword_set(smoothsize)) then nsmooth = smoothsize

nws = nw*ns

if filein && keyword_set(geofin) then begin
    geomin = fltarr(nn, nn, nimages)
    OPENR, lun1, geofin, /get_lun    		; Open geom input file
    readu, lun1, geomin			; 20100419 ELW - 'read'->'readu'
    CLOSE,lun1 & FREE_LUN,lun1
    ; geomin = ASSOC(lun1,FLTARR(nn,nn))
endif

geomout = fltarr(nn,nn, nimages)

;first, construct reference image from all images:
refim=fltarr(nn,nn)
FOR j=0,nimages-1 DO BEGIN
    refim = temporary(refim) + smooth(geomin(*,*,j),nsmooth)
ENDFOR
refim=temporary(refim)/float(nimages)
ntiles=nn / 16   ; number of tiles in x,y direction, for fringe pattern from DC offsets
tilesize=nn/ntiles

nt=3 ; number of derivative terms, will be augmented by '1' for dc-type offset.

coefs=fltarr(nimages,nn,nn,nt+1)
xgr = (findgen(nn) - tilesize/2) / tilesize
ygr=xgr ; assume square.
av=fltarr(nimages)
FOR j=0,nimages-1 DO BEGIN
      geoin=geomin(*,*,j)
      ; get rough mean lightlevel
      av(j)=rebin(geoin(ntiles:nn-ntiles-1,ntiles:nn-ntiles-1),1,1)
      IMAGE_DIF,geoin,refim,tilesize,cco
	; returns ntile x ntile coefficient arrays.
      ;; now make it big
      FOR k=0,nt DO BEGIN
         coefs(j,*,*,k) = INTERPOLATE(cco(*,*,k),xgr,ygr,/grid, cubic=-0.5)
      ENDFOR
ENDFOR

ff0=fltarr(nn,nn) ; this will be even-image fring pattern
ff1=fltarr(nn,nn) ; odd-image
FOR i=0,nimages/2 - 1 DO BEGIN
   ff0=ff0 + coefs(2*i,*,*,0)
   ff1=ff1 + coefs(2*i+1,*,*,0)
ENDFOR
ff0=ff0 / (nimages/2.) ; average it all out
ff1=ff1 / (nimages/2.) 

if (vb ge 1) then begin
   print,'  Magnitudes of even/odd residual correction:  '
   pmm,ff0(ntiles:nn-ntiles-1,ntiles:nn-ntiles-1),ff1(ntiles:nn-ntiles-1,ntiles:nn-ntiles-1)
endif

avmean=mean(av)

FOR i = 0, nimages/2 - 1 DO BEGIN
   twoi = 2*i
   geomout(*,*,twoi)   = geomin(*,*,twoi)   - ff0*av(twoi)/avmean
   geomout(*,*,twoi+1) = geomin(*,*,twoi+1) - ff1*av(twoi+1)/avmean
ENDFOR

if (fileout) && keyword_set(geofout) then begin
   OPENW, lun2, geofout, /get_lun
   WRITEU,lun2, geomout
   CLOSE, lun2 & FREE_LUN,lun2
endif

if (keyword_set(output4D)) then $
    geomout = reform(geomout, nn, nn, nf, nws, /overwrite)

END