-
Notifications
You must be signed in to change notification settings - Fork 0
/
lstdrv_segsrc.pro
128 lines (92 loc) · 2.82 KB
/
lstdrv_segsrc.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
PRO lstdrv_segsrc, sx,sy,pfwhm250,pfwhm350,pfwhm500,hdr250, map250,sig250 $
,hdr350, map350,sig350 $
,hdr500, map500,sig500,segsrc,bsn=bsn
; set SN threshold below which we will segment map
if(not keyword_set(bsn)) then bsn=0.5
prf250=psf_gaussian(fwhm=pfwhm250,npix=[11,11])
if(keyword_set(realin)) then prf250=prf250/max(prf250)
; build blend map
naxis1=sxpar(hdr250,'naxis1')
naxis2=sxpar(hdr250,'naxis2')
bim=fltarr(naxis1,naxis2)
print,'build blend map'
; 250 micron SN map
bim[sx,sy]=1.
flat=fltarr(49,49)+1.
cb250=convol(bim,flat) ; only interested in regions where there are sources (need this for speed)
bad=where(cb250 gt 1)
cb250[bad]=1.
sm250=cb250*map250/sig250
; blank is anything below some arbitrary SN threshold
blank250=where(sm250 le bsn, nblank250,comp=source250)
sm250[blank250]=-100.
smap250=intarr(naxis1,naxis2)
smap250[blank250]=-1
smap250[source250]=0
;350 micron SN map
sm350=map350/sig350
hastrom,sm350,hdr350,smap350,chdr350,hdr250,interp=1
smap350=smap350*cb250
blank350=where(smap350 le bsn, nblank350,comp=source350)
smap350[blank350]=-1.
smap350[source350]=0
;500 micron SN map
sm500=map500/sig500
hastrom,sm500,hdr500,smap500,chdr500,hdr250,interp=1
smap500=smap500*cb250
blank500=where(smap500 le bsn, nblank500,comp=source500)
smap500[blank500]=-1
smap500[source500]=0
smap=smap250*smap350*smap500
lim=where(smap ne 0)
smap[lim]=-1.
tacc=1
snum=0L
; flood fill each segment
while tacc ne 0 do begin
; find first unallocated pixel
tp=where(smap eq 0,tacc)
; fill that pixel with current segment
; number
snum=snum+1L
; try and fill surrounding pixels
sp=[tp[0]]
fill=where(smap[sp] eq 0.,fm)
k=0
stnum=snum
while fm ne 0 do begin
; check for mergers
merge=where(smap[sp] ne stnum and smap[sp] ne 0 and smap[sp] ne -1)
if(merge[0] ne -1) then begin
remark=where(smap eq stnum)
stnum=smap[sp[merge[0]]]
smap[remark]=stnum
endif
smap[sp[fill]]=stnum
if(smap[sp[0]] ne -1)then sp=[sp[0]-1,sp]
if(smap[sp[n_elements(sp)-1]] ne -1) then sp=[sp,sp[n_elements(sp)-1]+1]
fill=where(smap[sp] eq 0 ,fm)
if fm eq 0 then begin
good=where(smap[sp] eq stnum)
sp=sp[good]+naxis1
fill=where(smap[sp] ge 0,fm)
endif
endwhile
endwhile
; find segment for each source
n_src=n_elements(sx)
segsrc=lonarr(n_src)
for i=0L,n_src-1 do begin
segsrc[i]=smap[sx[i],sy[i]]
rad=1
while(segsrc[i] le 0) do begin
xrad=sx[i]+lindgen(2*rad+1)-rad
yrad=sy[i]+lindgen(2*rad+1)-rad
make_2d,xrad,yrad,dx,dy
tmap=smap[dx,dy]
tst=where(tmap gt 0)
if(tst[0] ne -1) then segsrc[i]=tmap[tst[0]]
rad=rad+1
endwhile
endfor
end