;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 ; ; WAITPERSIST: ; The number of dithers to wait until the pixels that were saturated ; are used in the mean or median of the mosaic. ; REJECTPERSIST: ; If this keyword is set, the pixels that have been saturated in ; a frame will be set to zero for WAITPERSIST frames afterwards. ; RADMAX: ; The radius (in pixels) from the center of a star for which the ; pixels will be zeroed if REJECTPERSIST is set. ;Star for M57: 3100, 3560 ;CALLING SEQUENCE: ; ; dither_slope_median,'/nfs/slac/g/ki/ki09/lances/Reduced/07Apr29/M16/M16GSlopes.lst',/rejectpersist 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, Rem60Hz = Rem60Hz, $ RejectPersist = RejectPersist, WaitPersist = WaitPersist, $ Verbose = Verbose, RadMax = RadMax If N_Elements(WaitPersist) eq 0 then WaitPersist = 3 ;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 If N_Elements(RadMax) eq 0 then RadMax = 50 ;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(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 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(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 and the saturation mask Fits_Read, RawFiles(FileNum), Frame, FitsHeader, exten_no = 0 Fits_Read, RawFiles(FileNum), SatMask, exten_no = 1 ;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 ;Allow subtraction of commond mode noise before dithering If KeyStr.Rem60Hz then $ Frame = Subtract_CM(Frame) ;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 ;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(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' Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$ Strtrim(KeyStr.ThisFilter,2)+'_'+$ 'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr+$ 'DitheredMedian.fits', MedIm, MedHeader Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$ Strtrim(KeyStr.ThisFilter,2)+'_'+$ 'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr+$ '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+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr+$ 'CenteredDitheredMedian.fits', MedCenIm, MedHeader Writefits, KeyStr.ObjectDir+KeyStr.ObjectName+'_'+$ Strtrim(KeyStr.ThisFilter,2)+'_'+$ 'Slope'+FlatFieldStr+KeyStr.Rem60HzStr+$ KeyStr.RejectPersistStr+$ 'CenteredDitheredMean.fits', MeanCenIm, MeanHeader stop end