PRO AUTOMUXSTRIPE,infile,outfile,thresh=thresh,VERBOSE=verbose,$$
                  CHECKIMAGE=checkimage

;
; NAME:
;        AUTOMUXSTRIPE
;
; PURPOSE:
;        Locate and correct muxstriping in IRAC images.
;
; EXPLANATION
;         "muxstriping" is related to muxbleed. Bright objects can
;         can cause a readout channel (column) to be depressed for
;         some number of rows, creating a horizontal "stripe" in the 
;         image. This routine works by medianing each readout channel
;         and looking to see if any have large numbers of statistically low
;         pixels. If they do, it scans that readout channel looking for
;         consecutive pixels that are low, then high. Once such a strip
;         is found, it is corrected.
;  
;
;
; REQUIREMENTS:
;         IDL 5.2 or higher, and the IDLastro package
;
; INPUTS:
;        infile - string containing the name of an IRAC FITS file.
;        outfile - string containing the name of a FITS file to
;                  which the output will be written.
;
;
; OPTIONAL INPUTS:
;        thresh - threshold in sigma for the stripe-finder. 3 is 
;                 recommended, and the default.
;        VERBOSE - if set, output diagnostic text.
;        CHECKIMAGE - if set, write a 4x256 diagnostic image called
;                     "check.fits".
;
;
;
;  EXAMPLE:
;       automuxstripe,'in.fits','out.fits',thresh=3,/VERBOSE,/CHECKIMAGE
;
;
;  NOTE:
;        The image needs to be scrubbed first for muxbleed and 
;        column-pulldown.
;
;        Also, the code was designed specifically for bright objects on
;        a nearly uniform background, i.e. the simplest case. Other cases
;        will require adapatation.
;
;        One refinement developed by the S-COSMOS teams and used by
;        S-COSMOS and SWIRE utilizes a wrapper that iterates the code
;        based on the number of 2MASS stars present in the frame.
;
;
;  HISTORY:
;       1.0 - 11/27/06, first release version, JAS
;
;


; set the defaults
if (keyword_set(thresh) NE 1) then thresh=3



bcd=readfits(infile,bcd_hd,/SILENT)

; define the arrays for the column-wise breakdown
collapse=fltarr(4,256)
hold1=fltarr(64)
hold2=fltarr(64)
hold3=fltarr(64)
hold4=fltarr(64)

; fill the columnwise arrays
for i=0, 255 do begin
  for j=0, 63 do begin

       hold1(j)=bcd(j*4,i)
       hold2(j)=bcd(j*4+1,i)
       hold3(j)=bcd(j*4+2,i)
       hold4(j)=bcd(j*4+3,i)

   endfor

    ; create a median column for each readout by taking median of just
    ;   just that readout
    collapse(0,i)=median(hold1)
    collapse(1,i)=median(hold2)
    collapse(2,i)=median(hold3)
    collapse(3,i)=median(hold4)

endfor 

; Removing any residual baseline and tilt in the columns
;  by fitting and subtracting a line from each one
;
if keyword_set(VERBOSE) then print,'Flattening collapsed rows.'
if keyword_set(VERBOSE) then print,'     Offset    Slope'
linex=findgen(256)
for i=0,3 do begin
   line=ladfit(linex,collapse(i,*))  ;  fit a line to each column in turn
   collapse(i,*)=collapse(i,*)-line(0)-line(1)*linex ; subtract the line fit
   if keyword_set(VERBOSE) then print,line(0),line(1)
endfor

; if set, write a 4x256 image containing the individual 
;  median readout columns
if keyword_set(CHECKIMAGE) then writefits,'check.fits',collapse

; find the overall noise level
;   note the factor of 8, to account for the averageing of the data
sig=robust_sigma(bcd)/8


; searching for which column has the stripe
;
;  this is done by measuring for each median readout the number of 
;    pixels _below_ the threshold, in sigma
;    the column witht he most low pixels is considered the "bad"
;    column for further analysis
badcolumn=intarr(4)
for i=0,255 do begin
   badcolumn(0)=badcolumn(0)+(collapse(0,i) LT (-thresh*sig))
   badcolumn(1)=badcolumn(1)+(collapse(1,i) LT (-thresh*sig))  
   badcolumn(2)=badcolumn(2)+(collapse(2,i) LT (-thresh*sig))  
   badcolumn(3)=badcolumn(3)+(collapse(3,i) LT (-thresh*sig))  
endfor

badcol=max(badcolumn)
badcol=!C
if keyword_set(VERBOSE) then begin
     print,'Column scores:'
     print,badcolumn
     print,'Bad Column is',badcol+1
endif


; now that the column containing the stripe has been identified
; begin search for it's location
; this is done by searching for three consecutive points 3-sigma low
stripestart=999
for i=255,10,-1 do begin

loscore=(collapse(badcol,i) LT (-thresh*sig))+$$
        (collapse(badcol,i-1) LT (-thresh*sig))+$$
        (collapse(badcol,i-2) LT (-thresh*sig))

if (loscore EQ 3)then begin
   stripestart=i
   if keyword_set(VERBOSE) then print,'Stripe starts at',stripestart
   i=1
   endif

endfor

; if the a starting point was found, now look for the endpoint
if (stripestart LT 256) then begin

; find end of stripe
thresh=1
stripestop=1
for i=(stripestart-5),2,-1 do begin

  hiscore=(collapse(badcol,i) GT (-thresh*sig))+$$
        (collapse(badcol,i-1) GT (-thresh*sig))+$$
        (collapse(badcol,i-2) GT (-thresh*sig))

if (hiscore EQ 3)then begin
   stripestop=i
   if keyword_set(VERBOSE) then print,'Stripe stops at',stripestop
   i=1
   endif

endfor

endif

out=bcd
if (stripestart GT 900) then begin
   if keyword_set(VERBOSE) then print,"No stripe detected."
   endif else begin

; if the stripe was found,execute the following


; determine mean offset
collapse2=collapse
collapse2(badcol,stripestop:stripestart)=-999
collapse3=collapse(where(collapse2 GT -900))
offset=median(collapse3)-median(collapse(badcol,stripestop:stripestart))
if keyword_set(VERBOSE) then print,'Stripe offset is',offset

; apply the offset
for i=0,63 do begin
     out(i*4+badcol,stripestop:stripestart)=bcd(i*4+badcol,stripestop:stripestart)+offset
endfor

; prepare header info
sxaddhist,'Automuxstripe applied.',bcd_hd
output='Bad column was'+string(badcol)
sxaddhist,output,bcd_hd
output='Offset was'+string(offset)
sxaddhist,output,bcd_hd
output='Started at row'+string(stripestart)
sxaddhist,output,bcd_hd
output='Stopped at row'+string(stripestop)
sxaddhist,output,bcd_hd


endelse

; finally, write the output
; always written, no matter what happens
writefits,outfile,out,bcd_hd



END