forward_function DI_overclock, DI_stat, DI_over, DI_qmask

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                                                                    ;
; This IDL function gives  - Mask of good/bad background pixels      ;
;                                                                    ;
; CALLING SEQUENCE:                                                  ;
; function DICAL_vismodulemask, inimage2, darkimage, imgCols         ;
;                                                                    ;
; REQUIRED INPUTS:                                                   ;
; A ncol*nrow science (inimage2) and dark (darkimage) images         ;
; The dimension of the frames (ncol, nrow via imgCols)               ;
;                                                                    ;
; OUTPUTS :                                                          ;
; The Mask to use to find the stripe features after dark subtraction ;
;                                                                    ;
; PROCEDURES USED (i.e. called directly!):                           ;
; DI_overclock:gives the frame dimensions and location of the OC pix ;
; DI_stat:gives statistics (overclock, in frame pixels)              ;
; DI_over:subtract the overclock pixel median value                  ;
; DI_qmask:Computes the mask to estimate the stripes in each quadrant;
;                                                                    ;
; Frederic F. Merlin                                                 ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function DI_overclock, ncol

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This function gives the location of the overclock pixels;
; and a test on the availability of the overclock pixels  ; 
;                                                         ;
; CALLING SEQUENCE:                                       ;
;DI_overclock, ncol                                       ;
;                                                         ;
; REQUIRED INPUTS:                                        ;
;ncol:size of the image                                   ;
;                                                         ;
; OUTPUTS:                                                ;
;A vector with the location of the limits over-clock - on ;
;frame pixels and a value that determines if there is an  ;
;over-clock area (1 if yes,0 if not)                      ;
;nbrpoints: used in DI_Stripes, number of pixels need to  ;
;compute the stripe value vs frame size                   ; 
;                                                         ;
;PROCEDURES USED (i.e. called directly!):                 ;
;                                                         ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

   overvector=dblarr(9)

   IF ncol EQ 1024 THEN begin
    xleft=[1,1,1016,1016]
    xright=[7,7,1022,1022]
    test=1 
    nbrpoints=32
   ENDIF

   IF ncol EQ 512 THEN begin
    xleft=[1,1,508,508]
    xright=[3,3,510,510]
    test=1
    nbrpoints=32
   ENDIF

   IF ncol EQ 256 THEN begin
    xleft=[0,0,252,252]
    xright=[3,3,255,255]
    test=1
    nbrpoints=16
   ENDIF

   IF ncol EQ 128 THEN begin
    xleft=[0,0,126,126]
    xright=[1,1,127,127]
    test=0
    nbrpoints=8
   ENDIF

   IF ncol EQ 64 THEN begin
    xleft=[0,0,63,63]
    xright=[0,0,63,63]
    test=0
    nbrpoints=8
   ENDIF

   IF ncol EQ 1008 THEN begin
    xleft=[0,0,1007,1007]
    xright=[0,0,1007,1007]
    test=0
    nbrpoints=32
   ENDIF

   overvector=[xleft[0],xleft[1],xleft[2],xleft[3],xright[0],xright[1],xright[2],xright[3],test,nbrpoints]
   return, overvector
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function DI_stat, image, overvector, q, lengthx, lengthy

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;this function gives statistics of the images: STANDEV,   ;
;MEAN and MEDIAN of the Overclock pixels                  ;
;                                                         ;
; CALLING SEQUENCE:                                       ;
;DI_stat, image, overvector, q, lengthx, lengthy          ;
;                                                         ;
; REQUIRED INPUTS:                                        ;
;image:The Science or Dark image                          ;
;lengthx, lengthy: Dimension of the quadrant area         ;
;overvector: Vector with the edges location of the frame  ;
;q:information on the quadrant used in the image          ;
;                                                         ;
; OUTPUTS:                                                ;
;A vector with the statistic of each quadrant (in the O-C ;
;areas)                                                   ;
;                                                         ;
;PROCEDURES USED (i.e. called directly!):                 ;
;                                                         ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; PIXEL POSITION
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Defines the number of column in the over-clock pixels area
   xleft=[overvector[0],overvector[1],overvector[2],overvector[3]]
   xright=[overvector[4],overvector[5],overvector[6],overvector[7]]
   limxinf=[xright[0]+3,xright[0]+3,lengthx+xright[0]+3,lengthx+xright[0]+3]
   limyinf=[xright[0]+3,lengthy+xright[0]+3,xright[0]+3,lengthy+xright[0]+3]
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; STATISTIC
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   A=image[xleft[q]:xright[q],limyinf[q]:limyinf[q]+lengthy-1]
   medianoverclockC=MEDIAN(A)
   meanoverclockC=MEAN(A)
   stddevC=STDDEV(A)
   pos=where(A GT medianoverclockC-stddevC*4.0 AND A LT medianoverclockC+stddevC*4.0)
   IF pos[0] GT 0 THEN meanoverclockC=MEAN(A(pos))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

   IF stddevC LT 0.5 THEN stddevC=1.0
   IF stddevC GT 3.0 THEN stddevC=3.0
   IF meanoverclockC GT medianoverclockC+2.0 THEN meanoverclockC=medianoverclockC

   C=[medianoverclockC,meanoverclockC,stddevC]
   return, C
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function DI_over, inimage, overvector, ncol, nrow

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Subtract the over-clock mean value for each quadrant    ;
;                                                        ;
; CALLING SEQUENCE:                                      ;
;DI_over, inimage, overvector, ncol, nrow                ;
;                                                        ;
; REQUIRED INPUTS:                                       ;
;inimage:The Science image                               ;
;ncol/nrow: Dimension of the image                       ;
;overvector: Vector with the edges location of the frame ;
;                                                        ;
; OUTPUTS:                                               ;
;A new image with median value of the over-clock pixels  ;
;subtracted. This new image is only used to create a mask;
;                                                        ;
;PROCEDURES USED (i.e. called directly!):                ;
;DI_stat.pro                                             ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; PIXEL POSITION
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Defines the size of the quadrants, removes part of the image external-edges
   xleft=[overvector[0],overvector[1],overvector[2],overvector[3]];left limit of OVpix
   xright=[overvector[4],overvector[5],overvector[6],overvector[7]];rigth limit of OVpix
   lengthx=ncol/2.0-(xright[0]+3)
   lengthy=nrow/2.0-(xright[0]+3)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; The procedure is running with a different quadrant order
;;; initial position=['LeftDown','LeftUp','RightDown','RightUp']

   FOR g=0,3, 1 DO BEGIN
   val=[1,0,3,2]
   q=val[g]

;;; Define the limits of the overclock pixels for each quadrant
   limxinf=[xright[0]+3,xright[0]+3,lengthx+xright[0]+3,lengthx+xright[0]+3]
   limyinf=[xright[0]+3,lengthy+xright[0]+3,xright[0]+3,lengthy+xright[0]+3]
   limx=[0,0,ncol/2.0,ncol/2.0]
   limy=[0,nrow/2.0,0,nrow/2.0]
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; STATISTIC
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Analyzes the science image
   statima=DI_stat(inimage,overvector,q,lengthx,lengthy)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; COSMETIC - BAD PIXELS REMOVAL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Remove the bad pixels in the OC
   FOR i=xleft[q],xright[q],1 DO BEGIN
    FOR j=limyinf[q],limyinf[q]+lengthy-1,1 DO BEGIN
     IF ABS(inimage[i,j]-statima[1]) GT statima[2]*4.0 THEN inimage[i,j]=statima[1]
    ENDFOR
   ENDFOR
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; DARK REMOVAL, ROW BY ROW
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;remove the dark and the stripe features from the OC pixels - only used for the mask computation
;;;We need to take into account the y-shift between FOV and OC pixels
   shiftoc=[-1,1,-1,1]

;;;(A) Computation of the Y-profile/Test on the lenght of the OC area (if there is)! 
   biasvec=fltarr(1,nrow/2.0)
   IF xright[q]-xleft[q]-1 GE 1 THEN BEGIN
     FOR j=0, nrow/2.0-1, 1 DO BEGIN
      biasvec[0,j]=(TOTAL(inimage[xleft[q]:xright[q],j+limy[q]])-MIN(inimage[xleft[q]:xright[q],j+limy[q]])-MAX(inimage[xleft[q]:xright[q],j+limy[q]]))/(xright[q]-xleft[q]-1.0)
     ENDFOR
   ENDIF ELSE BEGIN
     FOR j=0, nrow/2.0-1, 1 DO BEGIN
      biasvec[0,j]=(TOTAL(inimage[xleft[q]:xright[q],j+limy[q]]))/(xright[q]-xleft[q]+1.0)
     ENDFOR
   ENDELSE
;;;Smooth the profile in order to reduce the high frequency variation and adjust the OC level (+0.8) compared to the FOV level/
   biasvecs=biasvec
   FOR j=3, nrow/2.0-4, 1 DO BEGIN
    biasvecs[0,j]=TOTAL(biasvec[0,j-2+shiftoc[q]:j+2+shiftoc[q]])/5.0+1.0
   ENDFOR
;;;Add median value for the edges and for bad row
    biasvecs[0,0:2]=MEDIAN(biasvecs[0,*])
    biasvecs[0,nrow/2.0-3:nrow/2.0-1]=MEDIAN(biasvecs[0,*])

;;;Remove the stripe features at the first order
    FOR j=0, nrow/2.0-1, 1 DO BEGIN
     inimage[limx[q]:limx[q]+ncol/2.0-1,j+limy[q]]=inimage[limx[q]:limx[q]+ncol/2.0-1,j+limy[q]]-biasvecs[0,j]
    ENDFOR

   ENDFOR
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   return, inimage

END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function DI_qmask, quadrant, lengthx, lengthy, valuemax

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;This function computes the Mask with 'valuemax'          ;
;given in the main program                                ;
;                                                         ;
; CALLING SEQUENCE:                                       ;
;DI_qmask, quadrant, lengthx, lengthy, valuemax           ;
;                                                         ;
; REQUIRED INPUTS:                                        ;
;quadrant:The part of the image where the mask is computed;
;lengthx/lengthy: Dimension of the quadrant               ;
;valuemax: Limit value to compute the mask                ;
;                                                         ;
; OUTPUTS:                                                ;
;Mask frame for a quadrant                                ;
;                                                         ;
;PROCEDURES USED (i.e. called directly!):                 ;
;                                                         ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; The mask is initialized at 1
   mask=fltarr(lengthx, lengthy)+1.0

;;; Test on the value of quadrant[i,j] compared to the 'valuemax'
;;; Valuemax is an estimation of the upper limit (in DN) of the real expected background value.
   FOR j=0, lengthy-1, 1 DO BEGIN
   minimaline=MIN(quadrant[0:lengthx-1,j])
   IF minimaline GE -2 AND minimaline LE 2 THEN valuemax2=valuemax+minimaline
   IF minimaline LT -2 THEN valuemax2=valuemax-2.0; If flux or wrong line subtraction
   IF minimaline GT 2 THEN valuemax2=valuemax+1.0; IF flux or wrong line subtraction
    FOR i=0, lengthx-1, 1 DO BEGIN
     IF quadrant[i,j] GT valuemax2 or quadrant[i,j] LT minimaline-10.0 THEN mask[i,j] = 0.0
    ENDFOR
   ENDFOR

   return, mask
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;
; Main function ;
;;;;;;;;;;;;;;;;;

function DICAL_vismodulemask, inimage2, darkimage, imgCols

;;; LOCAL PARAMETERS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   ncol=imgCols[1]+imgCols[0]+1 
   nrow=ncol
   inimage2test = inimage2
;;; Defines the output images - initialization
   mask_final=dblarr(ncol,nrow); mask 
;;; Test on the image type size to compute the median overclock pixels (OVpix) value
   overvector=dblarr(9)
   overvector=DI_overclock(ncol)
   xleft=[overvector[0],overvector[1],overvector[2],overvector[3]];left limit of OVpix
   xright=[overvector[4],overvector[5],overvector[6],overvector[7]];rigth limit of OVpix

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; REMOVES THE DARK LEVEL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Removes the over-clock pixels median value if possible (overclock pixels exist i.e:testover=1)
   IF overvector[8] EQ 1 THEN BEGIN
    inimage2test=DI_over(inimage2, overvector, ncol, nrow)
    mask_final[0,*]=inimage2test[0,*]
    mask_final[ncol-1,*]=inimage2test[ncol-1,*]
   ENDIF ELSE BEGIN
;;; Else, use the value of the dark image
    inimage2test=inimage2-darkimage
   ENDELSE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

   BADPOS=where(inimage2test LT -10.0)
   IF BADPOS[0] GT 0 THEN inimage2test(BADPOS)=MAX(inimage2test)

;;; CREATES THE MASK IMAGE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Defines the size of the quadrants, takes into account the external edges of the FOV
   lengthx=ncol/2.0-(xright[0]+3)
   lengthy=nrow/2.0-(xright[0]+3)
;;; Defines the limits of the over-clock pixels areas for each quadrant
   limxinf=[xright[0]+3,xright[0]+3,lengthx+xright[0]+3,lengthx+xright[0]+3]
   limyinf=[xright[0]+3,lengthy+xright[0]+3,xright[0]+3,lengthy+xright[0]+3]
;;; Position of the quadrants=['LeftDown','LeftUp','RightDown','RightUp']

;;; The procedure is performed for each quadrant
   FOR g=0,3,1 DO BEGIN

;;; The procedure is running with a different quadrant order
   val=[1,0,3,2]
   q=val[g]

;;; Creates the Quadrants where the procedure will compute the mask
   quadrant_init=fltarr(lengthx,lengthy)+0.0
   quadrant_init[0:lengthx-1,0:lengthy-1]=1.0*inimage2test[limxinf[q]:limxinf[q]+lengthx-1,limyinf[q]:limyinf[q]+lengthy-1]

;;; Changes the initial mask (flag=0 or 1)
   mask_dn=DI_qmask(quadrant_init,lengthx,lengthy,5)
   mask_final[limxinf[q]:limxinf[q]+lengthx-1,limyinf[q]:limyinf[q]+lengthy-1]=mask_dn[0:lengthx-1,0:lengthy-1]

  ENDFOR

   mask_final[0,0]=overvector[9]

;;; Make some memory space, vectors to unit
inimage2test=0
mask_dn=0
quadrant_init=0

   return, mask_final
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

END