;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;pro combine dithers ; ;PURPOSE: ; To form a single image from a stack of images taken with the HxRG detectors. ; The file name should have the location of the dither as well as the time ; at which the exposure was taken. For each pixel, both a median and a ; mean after rejection of outliers is calculated. ; ; In some cases it might be advisable to star doing the centroids with a ; large centroid box and run the script once. Then the FILELIST will ; be updated with the calculated offsets. Running again with a smaller ; centroid box should then increase the accuracy. ; ;INPUTS: ; FITSDIR: (string) ; The name of the directory where BOTH the FILELIST and the fits files ; are located. ; FILELIST: (string) ; The list of files that are to be stacked. The format of this text ; file is the following: ; ; FitsFileName (\t) XOffset (\t) YOffset ; ; where FitsFileName is the fully qualified path of the file and XOffset ; and YOffset are the offsets of object positiions in that image relative to the ; object position in the first file of the stack. ; ;KEYWORDS: ; Keywords taken from KeywordStruct.pro ; ; WAITPERSIST: (int) ; The number of dithers to wait until the pixels that were saturated ; are used in the mean or median of the mosaic. ; REJECTPERSIST: (int) ; If this keyword is set, the pixels that have been saturated in ; a frame will be set to zero for WAITPERSIST frames afterwards. ; RADMAX: (int) ; The radius (in pixels) from the center of a star for which the ; pixels will be zeroed if REJECTPERSIST is set. ; DITHERCHECK: (int) ; 0 - Don't tv the image that will be written ; 1 - Tv the image that will be written ; STARREGIONSIZE: (int) ; The box that will be used to contain the star and the sky surrounding ; it measured in pixels. The box will be STARREGIONSIZE x STARREGIONSIZE ; SKIPSTAR: (int) ; The number of stars to skip over before the plotting begins in the ; case that TVFLAGLOCAL != 0 ; RADMAX: (int) ; The radius of pixels to include in the rejection of persistence. ; The centers will be determined from a centroid search of the SATMASK. ; VERBOSE: (int) ; If 1, there will be information printed to the screen. ; ;CALLING SEQUENCE: ; Combine_Dither_xRG, FitsDir, FileList ; ;EXAMPLES: ; Combine_Dithers_xRG, '/nfs/slac/g/ki/ki04/lances/H2RG-32-147/ASIC/Reduced/07Dec12/NGC2683/', 'NGC2683_I_Slopes.lst' ; Combine_Dithers_xRG, '/nfs/slac/g/ki/ki04/lances/H2RG-32-147/ASIC/Reduced/07Dec14/M77/','M77_I_Slopes.lst', ; Combine_Dithers_xRG, '/nfs/slac/g/ki/ki04/lances/H1RG-022/LEACH/Reduced/2007Nov18/Orion/','OrionTrapeziumGIY_Slopes.lst' ; pro Combine_Dithers_xRG, FitsDir, FileList, $ TvFlag = TvFlag, DitherCheck = DitherCheck, $ SkipStar = SkipStar,RejectPersist = RejectPersist, $ WaitPersist = WaitPersist, StarRegionSize=StarRegionSize, $ Verbose = Verbose, RadMax = RadMax, ObjectText=ObjectText ;KEYWORD Parameters If N_Elements(RejectPersist) eq 0 then RejectPersist = 0 If N_Elements(WaitPersist) eq 0 then WaitPersist = 3 If N_Elements(SkipStar) eq 0 then SkipStar = 0 If N_Elements(DitherCheck) eq 0 then DitherCheck = 0 If N_Elements(RadMax) eq 0 then RadMax = 20 If N_Elements(StarRegionSize) eq 0 then StarRegionSize=20 If N_Elements(ObjectText) eq 0 then ObjectText='' ;Make some common blocks Common KeyParams, KeyStr Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri ;Get defaults for keywords @KeywordStruct_xRG.pro @PlotSettings_xRG ;Move to the directory cd, FitsDir ;Check to see whether persistence rejection is going to happen If KeyStr.RejectPersist then begin KeyStr.RejectPersistStr = 'PersistReject'+Strtrim(WaitPersist,2)+'_' EndIf ;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 KeyStr.NumDithers = N_Elements(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 If KeyStr.TvFlag then window, xsize=700,ysize=700 ;Get a few kewyords, pointings, and offsets FileCheck = Return_File_Params_xRG(RawFiles(0), /slope) ;Try out a first guess at the offsets from our 20 arcsecond dithers Return_Dither_Offsets_xRG, XStart, YStart, TvFlagLocal=TvFlag, $ StarRegionSize=StarRegionSize, 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) - KeyStr.Ref NewArrYSize = KeyStr.Naxis2 + (MaxYOff - MinYOff) - KeyStr.Ref ;Form the Array for Airmasses and UTs AirMasses = [0.] UTs = [0.] ;Form an array big enough for all the dithers FullIms = DblArr(NewArrXSize, NewArrYSize, NumRawFiles) ColorArr = LonArr(NewArrXSize, NewArrYSize, 3) PersistArr = BytArr(KeyStr.Naxis1, KeyStr.Naxis2) ;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_xRG(MeanHeader) MkHdr, MedHeader, MedIm MedHeader = Update_Header_xRG(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 and the saturation mask Fits_Read, RawFiles(FileNum), Frame, FitsHeader, exten_no = 0 Fits_Read, RawFiles(FileNum), SatMask, exten_no = 1 ;Don't include the reference pixels If KeyStr.ElecStr eq 'ASIC' then begin Frame = Frame(KeyStr.Ref/2:KeyStr.Naxis1-KeyStr.Ref/2-1, $ KeyStr.Ref/2:KeyStr.Naxis2-KeyStr.Ref/2-1) SatMask = SatMask(KeyStr.Ref/2:KeyStr.Naxis1-KeyStr.Ref/2-1, $ KeyStr.Ref/2:KeyStr.Naxis2-KeyStr.Ref/2-1) EndIf Else Begin Frame = Frame(KeyStr.Ref/2+1:KeyStr.Naxis1-KeyStr.Ref/2-1, $ KeyStr.Ref/2:KeyStr.Naxis2-KeyStr.Ref/2-1) SatMask = SatMask(KeyStr.Ref/2+1:KeyStr.Naxis1-KeyStr.Ref/2-1, $ KeyStr.Ref/2:KeyStr.Naxis2-KeyStr.Ref/2-1) EndElse ;Get the Airmass and UTs KeyStr.AirMasses(FileNum) = Double(SxPar(FitsHeader,'AIRMASS')) KeyStr.UTs(FileNum) = Double(SxPar(FitsHeader,'UT')) print,min(frame), max(frame) SxAddPar, MedHeader, 'FILE'+Strtrim(FileNum,2), $ File_BaseName(RawFiles(FileNum)) SxAddPar, MeanHeader, 'FILE'+Strtrim(FileNum,2), $ File_BaseName(RawFiles(FileNum)) ;Zero the pixels that were previously saturated and If KeyStr.RejectPersist then begin PersistPix = where(PersistArr gt 0, NumPersist) If NumPersist gt 0 then begin Frame(PersistPix) = 0 PersistArr(PersistPix) -= 1 EndIf EndIf ;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 If KeyStr.ElecStr eq 'ASIC' then begin ;Put it in the right spot in the dithered large cube FullIms(XOff : XOff +KeyStr.Naxis1-KeyStr.Ref-1, $ YOff : YOff +KeyStr.Naxis2-KeyStr.Ref-1, $ FileNum) = Frame EndIf Else Begin ;Put it in the right spot, avoiding bad column FullIms(XOff : XOff +KeyStr.Naxis1-KeyStr.Ref-2, $ YOff : YOff +KeyStr.Naxis2-KeyStr.Ref-1, $ FileNum) = Frame EndElse ;Tag the saturated pixels in the previous array and set them to zero ;Wait at least 2 dithers to use these pixels If KeyStr.RejectPersist then begin PersistenceMask = Return_Persistence_Mask_xRG(SatMask, $ RadMax=RadMax, TvFlag=KeyStr.TvFlag) SatMask = SatMask + PersistenceMask PersistPix = where(SatMask gt 0,NumPersist) If NumPersist gt 0 then PersistArr(PersistPix) += WaitPersist EndIf 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 StarNum = 0, N_Elements(XStart)-1 do begin printf, OffsetsFile, Format= '(%"%I\t%I")', $ XStart(StarNum), YStart(StarNum) 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' ;Compute the Average AirMass using Stetson's formula with Simpson avg ValidAirMasses = KeyStr.AirMasses(0:NumRawFiles-1) ValidUTs = KeyStr.UTs(0:NumRawFiles-1) ValidUTsSort = ValidUTs(sort(ValidUTs)) ValidAirMasses = ValidAirMasses(sort(ValidUTs)) MiddleUT = ValidUTs(NumRawFiles/2) ;AirMass = (AM_beg+4*AM_middle+AM_last)/6 AvgAirMass = (ValidAirMasses(0) + 4*ValidAirMasses(NumRawFiles/2) +$ ValidAirMasses(NumRawFiles-1))/6 FxAddPar, MedHeader, 'AIRMASS', AvgAirMass, 'Average Airmass for mosaic' FxAddPar, MeanHeader, 'AIRMASS', AvgAirMass, 'Average Airmass for mosaic' FxAddPar, MedHeader, 'ITIME', KeyStr.ITIME, 'Integration time (slope = 1 sec.)' FxAddPar, MeanHeader, 'ITIME', KeyStr.ITIME, 'Integration time (slope = 1 sec.)' FxAddPar, MedHeader, 'UTMID', MiddleUT, 'Universal time at middle of mosaic' FxAddPar, MeanHeader, 'UTMID', MiddleUT, 'Universal time at middle of mosaic' ;Form the Filenames and write the data OutFileBase = KeyStr.RedObjectDir+KeyStr.ObjectName+'_'+$ Strtrim(KeyStr.ThisFilter,2)+'_'+$ 'SlopeFlatFielded'+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr OutMedFile = OutFileBase+ 'DitheredMedian.fits' OutMeanFile = OutFileBase+ 'DitheredMean.fits' Writefits, OutMedFile, MedIm, MedHeader Writefits, OutMeanFile, MeanIm, MeanHeader ;Form the image of the center file that has Naxis1 columns and Naxis2 rows XCenOff = -(Median(KeyStr.XOffs(0:NumRawFiles-1)) - MaxXOff) YCenOff = -(Median(KeyStr.YOffs(0:NumRawFiles-1)) - 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 ;Write the Centered Images OutCFileBase = KeyStr.RedObjectDir+KeyStr.ObjectName+'_'+$ Strtrim(KeyStr.ThisFilter,2)+'_'+$ 'SlopeFlatFielded'+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr OutCMedFile = OutFileBase+ 'CenteredDitheredMedian.fits' OutCMeanFile = OutFileBase+ 'CenteredDitheredMean.fits' Writefits, OutCMedFile, MedIm, MedHeader Writefits, OutCMeanFile, MeanIm, MeanHeader stop end