source: Sophya/trunk/Cosmo/RadioBeam/fgndsub.cc@ 3975

Last change on this file since 3975 was 3830, checked in by ansari, 15 years ago

Amelioration soustraction des avant-plans par l'introduction d'un fit polynome deg 2 sur ln(Temp)=f(ln(freq)), Reza 04/08/2010
File size: 9.0 KB
Line 
1/* ------------------------ Projet BAORadio --------------------
2 Classe ForegroundCleaner
3 R. Ansari , C. Magneville - Juin 2010
4--------------------------------------------------------------- */
5
6#include "fgndsub.h"
7#include "lobe.h"
8
9#include "cubedef.h"
10#include "matharr.h"
11#include "poly.h"
12
13#include "ctimer.h"
14
15/* --Methode-- */
16ForegroundCleaner::ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube)
17 : arrep_(arrep) , tbeam_(tbeam), skycube_(skycube)
18{
19 double dxdeg = ThetaSizeDegre/(double)NTheta;
20 double dydeg = PhiSizeDegre/(double)NPhi;
21 dx_ = DegreeToRadian(dxdeg);
22 dy_ = DegreeToRadian(dydeg);
23 freq0_ = Freq0MHz;
24 dfreq_ = FreqSizeMHz/(double)NFreq;
25 cout << " ForegroundCleaner: " << " dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx_ << " dy_rad=" << dy_ << ")"
26 << " Freq0=" << freq0_ << " deltaFreq=" << dfreq_ << " MHz" << endl;
27 skycube.Show();
28}
29
30/* --Methode-- */
31void ForegroundCleaner::BeamCorrections()
32{
33 BeamEffect beam(arrep_);
34 beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_);
35 cout << " ForegroundCleaner::BeamCorrections() done " << endl;
36}
37
38/* --Methode-- */
39int ForegroundCleaner::CleanNegatives(TF seuil)
40{
41 sa_size_t nneg = 0.;
42 for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++)
43 for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
44 for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++)
45 if (skycube_(kx, ky, kz) < seuil) {
46 nneg++; skycube_(kx, ky, kz)=seuil;
47 }
48 cout << " ForegroundCleaner::CleanNegatives " << nneg << " sky-pixels <" << seuil << " changed to" << seuil << endl;
49 return (int)nneg;
50}
51
52/* --Methode-- */
53int ForegroundCleaner::CleanPointSources(double nsigmas)
54{
55 Timer tm("ForegroundCleaner::CleanPointSources");
56 TArray< TF > sky2d(skycube_.SizeX(), skycube_.SizeY());
57 for(sa_size_t ky=0; ky<sky2d.SizeY(); ky++)
58 for(sa_size_t kx=0; kx<sky2d.SizeX(); kx++)
59 sky2d(kx, ky) = skycube_(Range(kx), Range(ky), Range::all()).Sum();
60
61
62 double mean, sigma;
63
64
65 TArray< TF > amz(1,1,skycube_.SizeZ()), asz(1,1,skycube_.SizeZ());
66 for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++) {
67 TArray< TF > slice = skycube_(Range::all() , Range::all(), Range(kz));
68 MeanSigma(slice, mean, sigma);
69 amz(0,0,kz)=mean; asz(0,0,kz)=sigma;
70 }
71
72 MeanSigma(sky2d, mean, sigma);
73 cout << " ForegroundCleaner::CleanPointSources 2D Sky projection, mean=" << mean << " sigma=" << sigma << endl;
74
75 TF seuil = (TF)(mean+nsigmas*sigma);
76
77 sa_size_t srccnt=0;
78 for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
79 for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) {
80 if (sky2d(kx,ky)>seuil) {
81 srccnt++;
82 skycube_(Range(kx), Range(ky), Range::all()) = amz;
83 }
84 }
85
86 cout << " Cleaned NSrc= " << srccnt << " 2D source/pixels (TotNPix=" << sky2d.Size()
87 << ")-> " << 100.*srccnt/sky2d.Size() << "% with S>" << seuil << " NSigmas=" << nsigmas << endl;
88 return (int)srccnt;
89}
90
91
92/* --Methode-- */
93TArray< TF > ForegroundCleaner::extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx)
94{
95 Timer tm("ForegroundCleaner::extractLSSCubeP1");
96// Inputs : maplss, mapsyc, freq0, dfreq
97// Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
98// Return_Array : foreground subtracted LSS signal
99 sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
100 synctemp.SetSize(2, sz);
101 specidx.SetSize(2, sz);
102 TArray<r_4> omap;
103 omap.SetSize(skycube_, true);
104 Vector vlnf(skycube_.SizeZ());
105 int nprt = 0;
106 // double freq0 : Frequence premier index en k (MHz)
107 // double dfreq : // largeur en frequence de chaque plan (Mhz)
108 for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
109 vlnf(k)=log((double)k*dfreq_+freq0_);
110
111 sa_size_t nbinfini=0;
112 sa_size_t nbbad=0;
113 sa_size_t imodprt=skycube_.SizeX()/6;
114 sa_size_t jmodprt=skycube_.SizeY()/6;
115 for(sa_size_t i=0; i<skycube_.SizeX(); i++)
116 for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
117 r_8 s1, sx, sx2, sy, sxy;
118 s1=sx=sx2=sy=sxy=0.;
119 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
120 double lnf = vlnf(k);
121 double ttot=(r_8)(skycube_(i,j,k));
122 if (ttot < 1.e-5) continue;
123 double lntt=log(ttot);
124 s1+=1.; sx+=lnf; sx2+=(lnf*lnf);
125 sy+=lntt; sxy+=(lnf*lntt);
126 }
127 double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
128 double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
129 double T0 = exp(beta+alpha*vlnf(0));
130
131 bool fgnan = false;
132 if (!isfinite(alpha)||(!isfinite(beta))) {
133 cout << "extractLSSCubeP1[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
134 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
135 alpha=beta=-999.;
136 fgnan = true; nbinfini++;
137 }
138 else {
139 double axp1 = beta+alpha*vlnf(0);
140 double axp2 = beta+alpha*vlnf(vlnf.Size()-1);
141
142 if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
143 cout << "extractLSSCubeP1[" << i << "," << j << "] BAD alpha=" << alpha
144 << " beta=" << beta << " T0=" << T0 << " - (mapsync="
145 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
146 fgnan = true; nbbad++;
147 }
148 }
149 if ((i%imodprt==0)&&(j%jmodprt==0))
150 cout << "extractLSSCubeP1[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha
151 << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
152 synctemp(i,j) = T0;
153 specidx(i,j) = alpha;
154 if (fgnan) {
155 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
156 }
157 else {
158 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
159 r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
160 omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
161 }
162 }
163 }
164 cout << " ForegroundCleaner::extractLSSCubeP1() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
165 return omap;
166}
167
168/*
169static inline val_polyn2(double alpha, double beta, double gamma, double x)
170{
171 return (beta+alpha*x+gamma*x*x);
172}
173*/
174
175/* --Methode-- */
176TArray< TF > ForegroundCleaner::extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx)
177{
178 Timer tm("ForegroundCleaner::extractLSSCubeP2");
179// Inputs : maplss, mapsyc, freq0, dfreq
180// Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
181// Return_Array : foreground subtracted LSS signal
182 sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
183 synctemp.SetSize(2, sz);
184 specidx.SetSize(2, sz);
185 TArray<r_4> omap;
186 omap.SetSize(skycube_, true);
187 Vector vlnf(skycube_.SizeZ());
188 Vector vlnT(skycube_.SizeZ());
189 int nprt = 0;
190 // double freq0 : Frequence premier index en k (MHz)
191 // double dfreq : // largeur en frequence de chaque plan (Mhz)
192 for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
193 vlnf(k)=log((double)k*dfreq_+freq0_);
194 vlnf -= vlnf(0);
195 // cout << " DBG*extractLSSCubeP2 vlnf(0)=" << vlnf(0) << " vlnf(1)=" << vlnf(1)
196 // << "vlnf(last)=" << vlnf(vlnf.Size()-1) << endl;
197 sa_size_t nbinfini=0;
198 sa_size_t nbbad=0;
199 sa_size_t imodprt=skycube_.SizeX()/6;
200 sa_size_t jmodprt=skycube_.SizeY()/6;
201 for(sa_size_t i=0; i<skycube_.SizeX(); i++)
202 for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
203 vlnT = -12.;
204 Poly polyn;
205 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
206 double lnf = vlnf(k);
207 double ttot=(r_8)(skycube_(i,j,k));
208 if (ttot < 1.e-5) continue;
209 vlnT(k)=log(ttot);
210 }
211 polyn.Fit(vlnf,vlnT,2);
212 double beta = polyn[0];
213 double alpha = polyn[1];
214 double gamma = polyn[2];
215 double T0 = exp(polyn(vlnf(0))); // exp( val_polyn2(alpha, beta, gamma, vlnf(0)) );
216
217 bool fgnan = false;
218 if (!isfinite(alpha)||(!isfinite(beta))||(!isfinite(gamma))) {
219 cout << "extractLSSCubeP2[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
220 << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
221 alpha=beta=gamma=-999.;
222 fgnan = true; nbinfini++;
223 }
224 else {
225 double axp1 = polyn(vlnf(0));
226 double axp2 = polyn(vlnf(vlnf.Size()-1));
227
228 if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
229 cout << "extractLSSCubeP2[" << i << "," << j << "] BAD alpha=" << alpha
230 // << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " axp1=" << axp1 << " axp2=" << axp2
231 << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " - (mapsync=" << skycube_(i,j,0)
232 << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
233 fgnan = true; nbbad++;
234 }
235 }
236 if ((i%imodprt==0)&&(j%jmodprt==0))
237 cout << "extractLSSCubeP2[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha << " gamma=" << gamma
238 << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
239 synctemp(i,j) = T0;
240 specidx(i,j) = alpha;
241 if (fgnan) {
242 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
243 }
244 else {
245 for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
246 r_4 fittedtemp = (r_4)( exp(polyn(vlnf(k))) );
247 omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
248 }
249 }
250 }
251 cout << " ForegroundCleaner::extractLSSCubeP2() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
252 return omap;
253}
254
Note: See TracBrowser for help on using the repository browser.