;function returnclassifyevents ;PURPOSE: ; This function takes in an Naxis1xNaxis2 image, Im1, and returns ; all of the sets of contiguous pixels to the user so that he or she ; may classify the event amongst 5 categories. The categories are ; ; Spot An event that is small (less than 3x3) ; Fake Star An event that looks and has the properties of a star ; Delta An event that splits into two separte tracks ; Worm An event that snakes around into wiggles ; Streak An event that looks like a straight trajectory ; ; For each set of contiugous pixels, it returns the second moments, ; ellipticity, length (as calculated by sqrt((MaxX-MinX)^2-(MaxY-MinY)^2), ; and a number of other paramters ;INPUTS: ; im1 An Naxis1xNaxis2 array with data (i.e. fits array) ; Naxis1 Number of Columns ; Naxis2 Number of Rows ; Threshold Value above which pixel is considered part of event ; Pitch The width/legnth of pixel in microns ; Thickness The thickness of CCD in microns ; ;RETURN VALUES: An Array of (11,Number of Events) quantities ; will be returned. To extract the values, use these ; indices. (t=thickness, L =track length) ; ; Lambda1Arr=ReturnArr(0,*) Max Second Moment ; Lambda2Arr=ReturnArr(1,*) Min Second Moment ; EArr=ReturnArr(2,*) Ellipticity ; FluxArr=ReturnArr(3,*) Total Counts ; TrackLength2DArr=ReturnArr(4,*) 2D Length in pixels ; TrackLength3DArr=ReturnArr(5,*) 3D Length in microns ; PerpCountsArr=ReturnArr(6,*) PerpCounts=Flux*t/L ; HitSizeArr=ReturnArr(7,*) Number of Pixels affected in Hit ; HitWidthArr=ReturnArr(8,*) Width = HitSize/L ; XExtentArr=ReturnArr(9,*) Max X value - Min X Value ; YExtentArr=ReturnArr(10,*) Max Y value - Min Y Value function classifyevents,Im1,Naxis1,Naxis2,MinThreshold,ThresholdSlope,$ Sigma,Pitch,Thickness,TvFlagOriginal,ClassFlag,$ OutDir,ParamStr TvFlag=TvFlagOriginal ImV=Im1 ;The indices of the Array NoPixels=N_elements(Im1) PixInd=lindgen(NoPixels) PixX=(PixInd mod Naxis1) PixY=(PixInd/Naxis1) ;///////////////////////////////////////////////////////////////////// ;DATA CONTAINERS NumHits=0. Lambda1Arr=[0.] ;The long second moment Lambda2Arr=[0.] ;The short second moment EArr=[0.] ;Scalar Ellipticity FluxArr=[0.] ;Total Counts in a region TrackLength2DArr=[0.] ;Length of the track in Pixels 2D TrackLength3DArr=[0.] ;Lenght of the track in physical units PerpCountsArr=[0.] ;Counts scaled by length HitSizeArr=[0.] ;Total number of pixels affected HitWidthArr=[0.] ;HitWidth as judged by Length XExtentArr=[0.] ;Extent of Hit in X pixels YExtentArr=[0.] ;Extent of Hit in Y pixels EventColor=[0.] ;Color of Event to Plot Deltas=[0.] ;An attempt at Discerning a delta Worms=[0.] ;An attempt at Discrening a worm NumQuantities=15 NoSpots=0 ;Number of "Spots" NoStreaks=0 ;Number of "Streaks" NoWorms=0 ;Number of "Worms" NoDeltas=0 ;Number of "Delta Rays" NoFakeStars=0 ;Number of Objects that look like stars NoUnidentified=0 ;Number of objects not identified as above ;////////////////////////////////////////////////////////////////////// CurX=0 CurY=0 NoHits=0. for Pix=0.,NoPixels-1 do begin ;Dark Images Have significant variation in dark current from top to bottom Threshold=MinThreshold+ThresholdSlope*float(CurY)/Naxis2 if Im1(Pix) gt Sigma*Threshold then begin ;Arrays to Hold Pixels ContX=[PixX(Pix)] ContY=[PixY(Pix)] ContInd=[PixInd(Pix)] ContCou=[Im1(Pix)] PossibleDelta=-1 PossibleWorm=-1 Im1(Pix)=0 CurSubInd=0 ;SubIndex NumPixelsSearched=0 ;NumPixels already searched NumPixelsTotal=1 ;Total NumPixels while NumPixelsSearched lt NumPixelsTotal do begin CurInd=ContInd(CurSubInd) CurX=ContX(CurSubInd) CurY=ContY(CurSubInd) ;A very slow, but exact way to find CR ;Being Careful About Edges and Corners if CurX eq 0 then begin LeftEdge=1 RightEdge=0 endif else if CurX eq Naxis1-1 then begin RightEdge=1 LeftEdge=0 endif else begin RightEdge=0 LeftEdge=0 endelse if CurY eq 0 then begin LowerEdge=1 UpperEdge=0 endif else if CurY eq Naxis2-1 then begin UpperEdge=1 LowerEdge=0 endif else begin LowerEdge=0 UpperEdge=0 endelse ;Look Through Adjacent Pixels. Set them to zero after they've been tacked If (RightEdge) ne 1 then begin ;RIGHT SRight=Im1(CurX+1,CurY) if SRight gt Threshold Then begin ContCou=[ContCou,SRight] ContInd=[ContInd,CurInd+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY)=0 endif Endif If (RightEdge ne 1) and (UpperEdge ne 1) then begin ;UP,RIGHT URight=Im1(CurX+1,CurY+1) if URight gt Threshold then begin ContCou=[ContCou,URight] ContInd=[ContInd,CurInd+Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY+1)=0 endif Endif If (UpperEdge ne 1) then begin ;UP SUp=Im1(CurX,CurY+1) if SUp gt Threshold then begin ContCou=[ContCou,SUp] ContInd=[ContInd,CurInd+Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX,CurY+1)=0 endif Endif If (LeftEdge ne 1) and (UpperEdge ne 1) then begin ;UP,LEFT ULeft=Im1(CurX-1,CurY+1) if ULeft gt Threshold then begin ContCou=[ContCou,ULeft] ContInd=[ContInd,CurInd+Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY+1)=0 endif Endif If (LeftEdge ne 1) then begin ;LEFT SLeft=Im1(CurX-1,CurY) if SLeft gt Threshold then begin ContCou=[ContCou,SLeft] ContInd=[ContInd,CurInd-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY)=0 endif Endif If (LeftEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,LEFT DLeft=Im1(CurX-1,CurY-1) if DLeft gt Threshold then begin ContCou=[ContCou,DLeft] ContInd=[ContInd,CurInd-Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY-1)=0 endif Endif If (LowerEdge ne 1) then begin ;DOWN SDown=Im1(CurX,CurY-1) if SDown gt Threshold then begin ContCou=[ContCou,SDown] ContInd=[ContInd,CurInd-Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX,CurY-1)=0 endif Endif If (RightEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,RIGHT DRight=Im1(CurX+1,CurY-1) if DRight gt Threshold then begin ContCou=[ContCou,DRight] ContInd=[ContInd,CurInd-Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY-1)=0 endif Endif NumPixelsSearched=NumPixelsSearched+1 CurSubInd=CurSubInd+1 endwhile If N_elements(ContX) gt 1 then begin ;Event Has been Found Now Classify it ////////////////////////////////// NoHits=NoHits+1. Flux=total(ContCou) ContXAvg=round(total(ContCou*ContX)/Flux) ;First Moment in X ContYAvg=round(total(ContCou*ContY)/Flux) ;First Moment in Y AbsX=ContX ;Maintain Absloute Coordinates AbsY=ContY ContX=ContX-ContXAvg ContY=ContY-ContYAvg FluxArr=[FluxArr,Flux] ContCountNorm=float(ContCou)/Flux ;Calculate Second Moments and find them in rotated frame I11=total(ContX^2*ContCou)/Flux I22=total(ContY^2*ContCou)/Flux I12=total(ContX*ContY*ContCou)/Flux E1=(I11-I22)/(I11+I22) E2=(2*I12)/(I11+I22) PA=.5*atan(E2,E1) XMat=[[I11,I12],[I12,I22]] Evals=Hqr(Elmhes(XMat),/Double) Lambda1Arr=[Lambda1Arr,Real_part(max(Evals))] Lambda2Arr=[Lambda2Arr,Real_part(min(Evals))] EArr=[EArr,E1^2+E2^2] ;Get Track Length and Number of Perpendicular Counts XExtentArr=[XExtentArr,Max(ContX)-Min(ContX)] YExtentArr=[YExtentArr,Max(ContY)-Min(ContY)] TrackLength3DArr=[TrackLength3DArr,$ float(sqrt((float(Pitch)*XExtentArr(NoHits))^2+$ (float(Pitch)*YExtentArr(NoHits))^2+$ float(Thickness)^2))] TrackLength2DArr=[TrackLength2DArr,sqrt((XExtentArr(NoHits)+1)^2+$ (YExtentArr(NoHits)+1)^2)] PerpCountsArr=[PerpCountsArr,FluxArr(NoHits)*Thickness/TrackLength3DArr(NoHits)] HitSizeArr=[HitSizeArr,N_elements(ContCou)] HitWidthArr=[HitWidthArr,HitSizeArr(NoHits)/TrackLength2DArr(NoHits)] ;Find Breaks in the Event to see if it is a delta ray or worm if HitSizeArr(NoHits) gt 5 then begin IndH=sort(ContInd) IndV=AbsX*Naxis1+AbsY XordH=ContX(IndH) YordH=ContY(IndH) XordV=ContX(sort(IndV)) YordV=ContY(sort(IndV)) XBreakH=XordH(1:HitSizeArr(NoHits)-1)-XordH(0:HitSizeArr(NoHits)-2) YBreakH=YordH(1:HitSizeArr(NoHits)-1)-YordH(0:HitSizeArr(NoHits)-2) XBreakV=XordV(1:HitSizeArr(NoHits)-1)-XordV(0:HitSizeArr(NoHits)-2) YBreakV=YordV(1:HitSizeArr(NoHits)-1)-YordV(0:HitSizeArr(NoHits)-2) PossibleDelta=where((((abs(XBreakH) gt 1) and (abs(YBreakH) eq 0)) or $ ((abs(XBreakH) eq 0) and (abs(YBreakH) gt 1))) or $ (((abs(XBreakV) gt 1) and (abs(YBreakV) eq 0)) or $ ((abs(XBreakV) eq 0) and (abs(YBreakV) gt 1)))) If PossibleDelta(0) eq -1 then begin Deltas=[Deltas,0] PossibleDelta=0 endif else begin Deltas=[Deltas,N_elements(PossibleDelta)] PossibleDelta=N_elements(PossibleDelta) endelse PossibleWorm=where((((abs(XBreakH) gt 3) and (abs(YBreakH) eq 0)) or $ ((abs(XBreakH) eq 0) and (abs(YBreakH) gt 3))) or $ (((abs(XBreakV) gt 3) and (abs(YBreakV) eq 0)) or $ ((abs(XBreakV) eq 0) and (abs(YBreakV) gt 3)))) If PossibleWorm(0) eq -1 then begin Worms=[Worms,0] PossibleWorm=0 endif else begin Worms=[Worms,N_elements(PossibleWorm)] PossibleWorm=N_elements(PossibleWorm) endelse endif else begin Deltas=[Deltas,0] PossibleDelta=0 Worms=[Worms,0] PossibleWorm=0 endelse if (TVFlag ge 0) and (TVFlag le 1) then begin erase !p.noerase=1 print,'Event Number = ',NoHits print,'Coordinates = ',CurX,' ',CurY print,'Ellipticity = ',Earr(noHits) print,'E Angle = ',PA print,'XExtent = ',XExtentArr(NoHits) print,'YExtent = ',YExtentArr(NoHits) print,'Lambda1 = ',Lambda1Arr(NoHits) print,'Lambda2 = ',Lambda2Arr(NoHits) print,'Lambda Ratio = ',Lambda2Arr(NoHits)/Lambda1Arr(NoHits) print,'Hit Size = ',HitSizeArr(NoHits) print,'TrackLength2D = ',TrackLength2DArr(NoHits) print,'Hit Width = ',HitWidthArr(NoHits) print,'PerpCounts = ',PerpCountsArr(NoHits) print,'Worm Index = ',PossibleWorm print,'Delta Index = ',PossibleDelta if TvFlag eq 1 then begin Device, Get_Screen_Size=ScreenSize window,xsize=.5*ScreenSize(0),ysize=.5*ScreenSize(1) tvscl,congrid($ (ImV(Min(AbsX):Max(AbsX),Min(AbsY):Max(AbsY))),$ 0.25*ScreenSize(0),0.45*ScreenSize(1),/center) colorbar,range=[min(ContCou),max(ContCou)],$ format='(f10.2)',$ position=[0.1,0.92,0.5,0.97],color=255,$ title='intensity' endif else if TvFlag eq 2 then begin tvscl,congrid($ (ImV(Min(AbsX):Max(AbsX),Min(AbsY):Max(AbsY))),$ 512,512,/center),$ 0.2,0.4,$ xsize=2.5,ysize=2.5,/inches !p.charsize=0.6 xyouts,0.05,0.94,'Ellipticity = '+$ strtrim(Earr(noHits),2),/norm xyouts,0.05,0.90,'E Angle = '+$ strtrim(PA,2),/norm xyouts,0.05,0.86,'XExtent = '+$ strtrim(XExtentArr(NoHits),2),/norm xyouts,0.05,0.82,'YExtent = '+$ strtrim(YExtentArr(NoHits),2),/norm xyouts,0.05,0.78,'Lambda1 = '+$ strtrim(Lambda1Arr(NoHits),2),/norm xyouts,0.05,0.74,'Lambda2 = '+$ strtrim(Lambda2Arr(NoHits),2),/norm xyouts,0.05,0.70,'Lambda Ratio = '+$ strtrim(Lambda2Arr(NoHits)/Lambda1Arr(NoHits),2),/norm xyouts,0.05,0.66,'Hit Size = '+$ strtrim(HitSizeArr(NoHits),2),/norm xyouts,0.05,0.62,'TrackLength2D = '+$ strtrim(TrackLength2DArr(NoHits),2),/norm xyouts,0.05,0.58,'Hit Width = '+$ strtrim(HitWidthArr(NoHits),2),/norm xyouts,0.05,0.54,'PerpCounts = '+$ strtrim(PerpCountsArr(NoHits),2),/norm endif endif HLineInd=where(ContY eq 0) VLineInd=where(ContX eq 0) HLineContX=ContX(HLineInd) HLineContCou=ContCou(HLineInd) VLineContY=ContY(VLineInd) VLineContCou=ContCou(VLineInd) VSort=sort(VLineContY) VLineContY=VLineContY(VSort) VLineContCou=VLineContCou(VSort) if TvFlag lt 2 then begin plot,VLineContY,VLineContCou-Threshold,$ psym=2,title='Vertical Profile at x=0',$ position=[0.55,0.05,0.95,.5],/norm plot,HLineContX,HLineContCou-Threshold,$ psym=2,title='Horizontal Profile at y=0',$ position=[0.55,0.55,0.95,.95],/norm endif ;Manual Classification if ClassFlag eq 1 then begin EventSuc=Select_Event(ScreenSize) ;Create Buttons If EventSuc eq 'Cursor Value' then begin curval while EventSuc eq 'Cursor Value' do begin EventSuc=Select_Event(ScreenSize) If EventSuc eq 'Cursor Value' then begin print, 'Pick Something Different' endif endwhile endif else if EventSuc eq 'Spot' then begin ;SPOT NoSpots=NoSpots+1 EventColor=[EventColor,0] print,"Classification: Spot" endif else if EventSuc eq 'Fake Star' then begin ;FAKE STAR NoFakeStars=NoFakeStars+1 EventColor=[EventColor,1] print,"Classification: Fake Star" endif else if EventSuc eq 'Worm' then begin ;WORM NoWorms=NoWorms+1 EventColor=[EventColor,2] print,"Classification: Worm" endif else if EventSuc eq 'Delta' then begin ;DELTA NoDeltas=NoDeltas+1 EventColor=[EventColor,3] print,"Classification: Delta" endif else if EventSuc eq 'Streak' then begin ;STREAK NoStreaks=NoStreaks+1 EventColor=[EventColor,4] print,"Classification: Streak" endif Endif else if ClassFlag eq 2 then begin if (XExtentArr(NoHits) le 3) and (YExtentArr(NoHits) le 3) then begin NoSpots=NoSpots+1 EventColor=[EventColor,0] ;Tiny Spot print,"Classification: Spot" endif else if (EArr(noHits) lt 0.10) and $ (Lambda1Arr(noHits) gt 0.9) and (Lambda2Arr(noHits) gt 0.9) and $ (Lambda2Arr(noHits)/Lambda1Arr(noHits) gt 0.8) and $ (Lambda2Arr(noHits)/Lambda1Arr(noHits) lt 1.2) $ then begin NoFakeStars=NoFakeStars+1 EventColor=[EventColor,1] print,"Classification: Fake Star" Tvflag=TvFlagOriginal endif else if Lambda1Arr(Nohits) gt 2 and $ Lambda2Arr(NoHits)/Lambda1Arr(NoHits) lt 0.45 then begin NoStreaks=NoStreaks+1 EventColor=[EventColor,4] print,"Classification: Streak" endif else if PossibleWorm(0) ne -1 then begin NoWorms=NoWorms+1 EventColor=[EventColor,2] print,"Classification: Worm" endif else if PossibleDelta(0) ne -1 then begin NoDeltas=NoDeltas+1 EventColor=[EventColor,3] print,"Classification: Delta" endif else begin NoUnidentified=NoUnidentified+1 EventColor=[EventColor,5] print,"Classification: Unidentified" endelse EndIf endif endif endfor ReturnArr=fltarr(NumQuantities,NoHits) ReturnArr(0,*)=Lambda1Arr(1:NoHits) ReturnArr(1,*)=Lambda2Arr(1:NoHits) ReturnArr(2,*)=EArr(1:NoHits) ReturnArr(3,*)=FluxArr(1:NoHits) ReturnArr(4,*)=TrackLength2DArr(1:NoHits) ReturnArr(5,*)=TrackLength3DArr(1:NoHits) ReturnArr(6,*)=PerpCountsArr(1:NoHits) ReturnArr(7,*)=HitSizeArr(1:NoHits) ReturnArr(8,*)=HitWidthArr(1:NoHits) ReturnArr(9,*)=XExtentArr(1:NoHits) ReturnArr(10,*)=YExtentArr(1:NoHits) ReturnArr(11,*)=EventColor(1:NoHits) ReturnArr(12,*)=Deltas(1:NoHits) ReturnArr(13,*)=Worms(1:NoHits) ReturnArr(14,0)=NoSpots ReturnArr(14,1)=NoDeltas ReturnArr(14,2)=NoWorms ReturnArr(14,3)=NoFakeStars ReturnArr(14,4)=NoStreaks ReturnArr(14,5)=NoUnidentified get_lun,EventClassFile openw,EventClassFile,OutDir+ParamStr+'.class' printf,EventClassFile,$ format='(%"%s\t\t%s\t\t%s\t\t%s\t\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s")',$ 'Lambda1','Lambda2','E','Flux',$ 'TrackL 2D','TrackL 3D ','Perp Counts','Hit Size ',$ 'HitWidth ','X Extent ','YExtent ','EventColor',$ 'Pos. Worm','Poss Delta' for EventNo=0, NoHits-1 do begin printf,EventClassFile,$ format='(%"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%d\t\t%d\t\t%d\t\t%d\t\t%d")',$ ReturnArr(0,EventNo),ReturnArr(1,EventNo),ReturnArr(2,EventNo),$ ReturnArr(3,EventNo),ReturnArr(4,EventNo),ReturnArr(5,EventNo),$ ReturnArr(6,EventNo),ReturnArr(7,EventNo),ReturnArr(8,EventNo),$ ReturnArr(9,EventNo),ReturnArr(10,EventNo),ReturnArr(11,EventNo),$ ReturnArr(12,EventNo),ReturnArr(13,EventNo) endfor close,EventClassFile free_lun,EventClassFile !p.noerase=0 return,ReturnArr end