-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlstdrv_initsolvesp.pro
144 lines (109 loc) · 4.29 KB
/
lstdrv_initsolvesp.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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
pro lstdrv_initsolveSP,sx,sy,segsrc,prf,fwhm,paxis,noisy_map,sig_map,x_pix,y_pix,seg,bl,bu,psrc,psig,pmat,amat2,sn_map,bvec,xvec,iBIC,solve=solve,sx_pix=sx_pix,sy_pix=sy_pix
snsrc=n_elements(segsrc)
paxis1=paxis[0]
paxis2=paxis[1]
prf=psf_gaussian(fwhm=fwhm,npixel=[paxis1,paxis2])
ssx=sx[segsrc]
ssy=sy[segsrc]
mux=max(ssx,min=mlx)
muy=max(ssy,min=mly)
npix=where(x_pix lt mux+20 and y_pix lt muy+20 and y_pix ge mly-20 and x_pix ge mlx-20,snpix)
sx_pix=x_pix[npix]
sy_pix=y_pix[npix]
snoisy_map=noisy_map[npix]
ssig_map=sig_map[npix]
pmat=dblarr(snsrc,snsrc)
bvec=dblarr(snsrc)
amat=dblarr(snpix,snsrc)
amat2=dblarr(snpix,snsrc)
for s=0L,snsrc-1 do begin
dx = -round(ssx[s])+(paxis1-1.)/2.+sx_pix
dy = -round(ssy[s])+(paxis2-1.)/2.+sy_pix
pindx=findgen(paxis1)-ssx[s]+round(ssx[s])
pindy=findgen(paxis2)-ssy[s]+round(ssy[s])
px=ssx[s]-round(ssx[s])+(paxis1-1.)/2.
py=ssy[s]-round(ssy[s])+(paxis2-1.)/2.
good = where(dx GE 0 AND dx LT paxis1 AND dy GE 0 AND dy LT paxis2, ngood, comp=bad, ncomp=nbad)
if(ngood gt 0.5*n_elements(prf)) then begin
make_2d,pindx,pindy,ipx2,ipy2
nprf=interpolate(prf,ipx2,ipy2,cubic=-0.5)
amat[good,s]=nprf[round(dx[good]),round(dy[good])]/sqrt(ssig_map[good])
amat2[good,s]=nprf[round(dx[good]),round(dy[good])]/(ssig_map[good])
if(s eq 0) then begin
arow=fltarr(ngood)+s
acol=good
endif else begin
arow=[arow,fltarr(ngood)+s]
acol=[acol,good]
endelse
endif
endfor
; do sparse multiplication
t1=systime(/sec)
pmat=dblarr(snsrc,snsrc)
sindex=sort(acol)
ncnt=n_elements(acol)
; create pmat
cnt=0L
while(cnt lt ncnt) do begin
wpix=acol[sindex[cnt]]
wind=cnt
nipix=1L
cnt=cnt+1
if(cnt lt ncnt) then begin
while(acol[sindex[cnt]] eq wpix ) do begin
nipix=nipix+1L
wind=[wind,cnt]
cnt=cnt+1L
if(cnt eq ncnt) then break
endwhile
endif
ipix=sindex[wind]
irow=arow[ipix]
make_2d,irow,irow,tPx,tPy
make_2d,ipix,ipix,tIx,tIy
tPx=reform(tPx,nipix^2)
tPy=reform(tPy,nipix^2)
tIx=reform(tIx,nipix^2)
tIy=reform(tIy,nipix^2)
pmat[tPx,tPy]=pmat[tPx,tPy]+amat[acol[tIx],arow[tIx]]*amat[acol[tIy],arow[tIy]]
endwhile
;psig
psig=dblarr(snsrc,snsrc)
sindex=sort(acol)
ncnt=n_elements(acol)
; create pmat
cnt=0L
while(cnt lt ncnt) do begin
wpix=acol[sindex[cnt]]
wind=cnt
nipix=1L
cnt=cnt+1
if(cnt lt ncnt) then begin
while(acol[sindex[cnt]] eq wpix ) do begin
nipix=nipix+1L
wind=[wind,cnt]
cnt=cnt+1L
if(cnt eq ncnt) then break
endwhile
endif
ipix=sindex[wind]
irow=arow[ipix]
make_2d,irow,irow,tPx,tPy
make_2d,ipix,ipix,tIx,tIy
tPx=reform(tPx,nipix^2)
tPy=reform(tPy,nipix^2)
tIx=reform(tIx,nipix^2)
tIy=reform(tIy,nipix^2)
psig[tPx,tPy]=psig[tPx,tPy]+amat2[acol[tIx],arow[tIx]]*amat2[acol[tIy],arow[tIy]]
endwhile
; create bvec
bvec=dblarr(snsrc)
for i=0L,snsrc-1 do begin
ipix=where(arow eq i)
if(ipix[0] ne -1) then bvec[i]=total(amat2[acol[ipix],arow[ipix]]*snoisy_map[acol[ipix]]/ssig_map[acol[ipix]])
endfor
print,'finished constructing arrays ',systime(/sec)-t1
sn_map=snoisy_map/(ssig_map)
xvec=dblarr(snsrc)
end