;pro dither_median
;
;PURPOSE:
;	To form a single image from a 3x3 set of dithered images.  
;	For each image, the median sky at a given pixel is removed.
;
;KEYWORDS:
;	Keywords taken from KeywordStruct.pro file in same directory as this file
;
;Star for M57: 3100, 3560

pro Dither_Slope_Median, FileList, XStart = XStart, YStart=YStart, $
		   RootDir=RootDir, Date=Date, $
		   ReducedDir = ReducedDir, $
		   SkySubtract= SkySubtract, ObjectName=ObjectName, $
		   NReads = NReads, ThisFilter = ThisFilter, $
		   FlatFielded = FlatFielded, SkipFileNum = SkipFileNum,$
		   TvFlag = TvFlag, DitherCheck = DitherCheck, SkipStar = SkipStar

  ;Make some common blocks
  Common KeyParams, KeyStr
  Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

  ;Get defaults for keywords
  @KeywordStruct.pro
  @PlotSettings
  cd, ReducedDir

  If N_Elements(SkipFileNum) eq 0 then SkipFileNum = 0
  If FlatFielded then FlatFieldStr = 'FlatFielded' else FlatFieldStr = ''
  If N_Elements(SkipStar) eq 0 then SkipStar = 0
  If N_Elements(DitherCheck) eq 0 then DitherCheck = 0

  ;Determine whether files should be taken from a list or searched for
  If N_Elements(FileList) eq 0 then begin
     ;Get the files from a directory search
     RawFiles = File_Search(KeyStr.ObjectName+'*'+strtrim(KeyStr.NReads, 2) +$
                 '_Reads*'+FlatFieldStr+'*.fits')
     RawFiles = Return_File_List(RawFiles, ThisFilter=ThisFilter)
     RawFiles = RawFiles(SkipFileNum:N_Elements(RawFiles)-1)
     print, 'Dithering with Files : '
     print, RawFiles
     ObjectDir = ReducedDir+KeyStr.ObjectName+PathDelim
     File_MkDir, ObjectDir
  EndIf else begin
     ;Get the files from a list in a file and the offsets
     ReadCol, FileList, RawFiles, XOff, YOff, Format = 'A,I,I'
     ;Put the info in the structure
     KeyStr.DitherFiles(0:N_Elements(RawFiles)-1) = RawFiles

     ;Get the bright star coordinates from the same file
     ReadCol, FileList, XStart, YStart, Format = 'I,I', $
       SkipLine = N_Elements(RawFiles)
     For i = 0, N_Elements(XStart)-1 do begin 
       KeyStr.RAs(0:N_Elements(RawFiles)-1,i)  = XOff
       KeyStr.DECs(0:N_Elements(RawFiles)-1,i) = YOff
     EndFor
  EndElse 

  If KeyStr.TvFlag then window, xsize=700,ysize=700
  ;Get a few kewyords, pointings, and offsets
  FileCheck = Return_File_Params(RawFiles(0), /slope)	
  ;Try out a first guess at the offsets from our 20 arcsecond dithers
  Return_Dither_Offsets, XStart, YStart, TvFlag=TvFlag, SkipStar = SkipStar
  NumRawFiles = N_elements(RawFiles)

  ;Find the Maximums for the dimensions of the array
  MinXOff  = Min(KeyStr.XOffs)  & MaxXOff  = Max(KeyStr.XOffs)
  MinYOff  = Min(KeyStr.YOffs)  & MaxYOff = Max(KeyStr.YOffs)
  NewArrXSize = KeyStr.Naxis1 + (MaxXOff - MinXOff)
  NewArrYSize = KeyStr.Naxis2 + (MaxYOff - MinYOff)

  ;Form an array big enough for all the dithers
  FullIms  = DblArr(NewArrXSize, NewArrYSize, NumRawFiles)
  ColorArr = LonArr(NewArrXSize, NewArrYSize, 3)

  ;Start arrays and header that will contain the names of the files used in the
  ;dither. Form Array for Median and Mean of Non-Outliers of the stack
  MedIm  = FltArr(NewArrXSize, NewArrYSize)
  MeanIm = FltArr(NewArrXSize, NewArrYSize)
  MkHdr, MeanHeader, MeanIm
  MeanHeader = Update_Header(MeanHeader)
  MkHdr, MedHeader, MedIm
  MedHeader = Update_Header(MedHeader)

  ; ****************************************************DITHER
  ; Form Dithered Image with Desired Filter
  For FileNum= 0, NumRawFiles - 1 do begin
      print, 'Dithering with File : ' + RawFiles(FileNum)

      ;Read in the frame
      Fits_Read,  RawFiles(FileNum), Frame, FitsHeader
      print,min(frame), max(frame)
      FxAddPar, MedHeader, 'FILE'+Strtrim(FileNum,2), $
	  RawFiles(FileNum)
      FxAddPar, MeanHeader, 'FILE'+Strtrim(FileNum,2), $
	  RawFiles(FileNum)

      ;tvscl, congrid(Frame,512,512)

      ;Shift the images in the opposite direction to stack
      XOff = -(KeyStr.XOffs(FileNum)-MaxXOff)
      YOff = -(KeyStr.YOffs(FileNum)-MaxYOff)
 
      If DitherCheck eq 1 then begin
         ColorArr(XOff:XOff+KeyStr.Naxis1-1, $
		  YOff : YOff +KeyStr.Naxis2-1, $
		  *) = $
	          [[[replicate(ColNumStr.ColNums(0,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]],$ 
		   [[replicate(ColNumStr.ColNums(1,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]],$
		   [[replicate(ColNumStr.ColNums(2,FileNum),$
		   KeyStr.Naxis1,KeyStr.Naxis2)]]]
	 !p.noerase = 1
         tvimage, congrid(colorarr,700,700,3), true=3,/tv
      EndIf
      ;Put it in the right spot in the dithered large cube
      FullIms(XOff : XOff +KeyStr.Naxis1-1, $
	      YOff : YOff +KeyStr.Naxis2-1, $
	      FileNum) = Frame

  Endfor

  ; **********************************************************
  ; Run through the stacks 
  For Row = 0, NewArrYSize - 1 do begin 
    If Row mod (KeyStr.Naxis2/4) eq 0 then print, 'Got to Row: '+Strtrim(Row,2)

    For Col = 0, NewArrXSize - 1 do begin 
        PixelStack = FullIms(Col,Row,*)
        ;First off: remove Nans
        OutRange = where(finite(PixelStack) eq 0, NumOut,Complement =InRange )
	If NumOut gt 0 then PixelStack = PixelStack(InRange)

	;Next off: only deal with positive count values
        GoodPix    = where(PixelStack gt 0, NumGood)

	If NumGood eq 0 then begin
           ;No good values - set to zero
           MedIm(Col,Row)  = 0
           MeanIm(Col,Row) = 0
	EndIf Else If NumGood eq 1 then begin
           ;Only one good value - go with it
           MedIm(Col,Row)  = PixelStack(GoodPix)
           MeanIm(Col,Row) = PixelStack(GoodPix)
        EndIf Else Begin
	   ;More than one value - do some cases
           PixelStackKeep = PixelStack(GoodPix)
           MeanPixStack = Mean(PixelStackKeep)
           StdDevPixStack = StdDev(PixelStackKeep)
           If NumGood eq 2 then begin
	     ;Two Good Values, take mean for both since median chooses higher
             ;Check to see if one is very high - most likely semi-hot pixel
             ;so choose lower value
             If Max(PixelStackKeep) gt 3*Min(PixelStackKeep) then begin
		MedIm(Col,Row)  = Min(PixelStackKeep)
		MeanIm(Col,Row) = Min(PixelStackKeep)
             EndIf Else Begin 
	        MedIm(Col,Row)  = Mean(PixelStack(GoodPix))
                MeanIm(Col,Row) = Mean(PixelStack(GoodPix))
	     EndElse
           EndIf Else Begin
             PixelStackKeep = PixelStack(GoodPix)
             ;Take the median - should reject outliers by definition
             MedIm(Col,Row) = Median(PixelStackKeep)

	     ;Do sigma clipping for the mean value
     	     Rej = 1
             While Rej eq 1 do begin

	       ;Look for an obvious outlier hot pixel 
               If Max(PixelStackKeep) gt 3* Min(PixelStackKeep) then $
		  PixelStackKeep = PixelStackKeep($
			where(PixelStackKeep ne max(PixelStackKeep)))
	       If N_Elements(PixelStackKeep) lt 2 then begin
		  MeanIm(Col,Row) = Mean(PixelStackKeep)
		  break
	       EndIf
	       ;Look for the mean - Reject values outside the mean
	       MedianPixStack = Median(PixelStackKeep)
	       MeanPixStack = Mean(PixelStackKeep)
               StdDevPixStack = StdDev(PixelStackKeep)
	     
	       SigThresh = 1.75
	       ;Do a few simple clips
	       If MeanPixStack gt 2*MedianPixStack then SigThresh = 1
	     
	       ClipInds=where($
	       abs(PixelStackKeep-MeanPixStack) gt SigThresh*StdDevPixStack,$
		  NumClip,complement = KeepInds)
               If NumClip gt 0 then begin
	          PixelStackKeep = PixelStackKeep(KeepInds)
                  Rej = 1
	       EndIf Else begin
                  MeanIm(Col,Row) = Mean(PixelStackKeep)
		  Rej = 0
               EndElse
	     EndWhile

	   EndElse
	EndElse
     EndFor
  EndFor

  ;**************************************************************************
  ;Update the file that holds the list of filenames and offsets with the 
  ;current ones. Don't change the star coordinates
  get_lun, OffsetsFile
  openw, OffsetsFile, FileList
  For FileNum = 0, NumRawFiles-1 do begin
     printf, OffsetsFile, Format= '(%"%s\t%I\t%I")', $
	RawFiles(FileNum), KeyStr.XOffs(FileNum), KeyStr.YOffs(FileNum)
  EndFor
  printf, OffsetsFile, Format = '(%"%s")', $
	';Bright Star Coordinates'
  For FileNum = 0, NumRawFiles-1 do begin
     printf, OffsetsFile, Format= '(%"%I\t%I")', $
        XStart(FileNum), YStart(FileNum)
  EndFor
  close, OffsetsFile
  free_lun, OffsetsFile 


  ;Write the full fits files and find the offset of the most populated 
  ;area of the image
  XCenOff = -(KeyStr.XOffs(0) - MaxXOff)
  YCenOff = -(KeyStr.YOffs(0) - MaxYOff)
  FxAddPar, MedHeader, 'DITHXOFF', XCenOff , 'Offset in X of Central Dither'
  FxAddPar, MedHeader, 'DITHYOFF', YCenOff , 'Offset in Y of Central Dither'
  FxAddPar, MeanHeader, 'DITHXOFF', XCenOff , 'Offset in X of Central Dither'
  FxAddPar, MeanHeader, 'DITHYOFF', YCenOff , 'Offset in Y of Central Dither'

  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
	'Slope'+FlatFieldStr+'DitheredMedian.fits', MedIm, MedHeader
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
	'Slope'+FlatFieldStr+'DitheredMean.fits', MeanIm, MeanHeader

  ;Form the image of the center file that has Naxis1 columns and Naxis2 rows
  XCenOff = -(KeyStr.XOffs(0) - MaxXOff) 
  YCenOff = -(KeyStr.YOffs(0) - MaxYOff)
  MeanCenIm = MeanIm(XCenOff : XCenOff + KeyStr.Naxis1-1, $
		     YCenOff : YCenOff + KeyStr.Naxis2-1 )
  MedCenIm  =  MedIm(XCenOff : XCenOff + KeyStr.Naxis1-1, $
                     YCenOff : YCenOff + KeyStr.Naxis2-1 )

  ;Modify the headers
  FxAddPar, MedHeader,  'NAXIS1', KeyStr.Naxis1
  FxAddPar, MedHeader,  'NAXIS2', KeyStr.Naxis2
  FxAddPar, MeanHeader, 'NAXIS1', KeyStr.Naxis1
  FxAddPar, MeanHeader, 'NAXIS2', KeyStr.Naxis2
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
        'Slope'+FlatFieldStr+'CenteredDitheredMedian.fits', MedCenIm, MedHeader
  Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$
	Strtrim(KeyStr.ThisFilter,2)+'_'+$
        'Slope'+FlatFieldStr+'CenteredDitheredMean.fits', MeanCenIm, MeanHeader

  stop


end