| 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-- */
|
|---|
| 16 | PowerLawChecker::PowerLawChecker(TArray< TF >& skycube)
|
|---|
| 17 | : skycube_(skycube)
|
|---|
| 18 | {
|
|---|
| 19 | }
|
|---|
| 20 | /* --Methode-- */
|
|---|
| 21 | void PowerLawChecker::CheckXYMean()
|
|---|
| 22 | {
|
|---|
| 23 | // double freq0 : Frequence premier index en k (MHz)
|
|---|
| 24 | // double dfreq : // largeur en frequence de chaque plan (Mhz)
|
|---|
| 25 | double freq0_ = Freq0MHz;
|
|---|
| 26 | double dfreq_ = FreqSizeMHz/(double)NFreq;
|
|---|
| 27 |
|
|---|
| 28 | double tempfirst,templast;
|
|---|
| 29 | r_8 s1, sx, sx2, sy, sxy;
|
|---|
| 30 | s1=sx=sx2=sy=sxy=0.;
|
|---|
| 31 | double lnf0=0.;
|
|---|
| 32 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 33 | double lnf=log((double)k*dfreq_+freq0_);
|
|---|
| 34 | double ttot = Mean(skycube_(Range::all(), Range::all(), Range(k)));
|
|---|
| 35 | if (k==0) { tempfirst=ttot; lnf0=lnf; }
|
|---|
| 36 | if (k==skycube_.SizeZ()-1) templast=ttot;
|
|---|
| 37 | if (ttot < 1.e-5) continue;
|
|---|
| 38 | double lntt=log(ttot);
|
|---|
| 39 | s1+=1.; sx+=lnf; sx2+=(lnf*lnf);
|
|---|
| 40 | sy+=lntt; sxy+=(lnf*lntt);
|
|---|
| 41 | }
|
|---|
| 42 | double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
|
|---|
| 43 | double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
|
|---|
| 44 | double T0 = exp(beta+alpha*lnf0);
|
|---|
| 45 |
|
|---|
| 46 | cout << " PowerLawChecker::CheckMean() meanTemp(0 ...last) " << tempfirst << " ... "
|
|---|
| 47 | << templast << endl;
|
|---|
| 48 | bool fgnan = false;
|
|---|
| 49 | if (!isfinite(alpha)||(!isfinite(beta))) {
|
|---|
| 50 | cout << "ePowerLawChecker::CheckMean() Not finite alpha, beta " << endl;
|
|---|
| 51 | fgnan = true;
|
|---|
| 52 | }
|
|---|
| 53 | cout << "PowerLawChecker::CheckMean() - T0=" << T0 << " alpha=" << alpha << "(beta=" << beta << ")" << endl;
|
|---|
| 54 | return;
|
|---|
| 55 | }
|
|---|
| 56 |
|
|---|
| 57 | /* --Methode-- */
|
|---|
| 58 | ForegroundCleaner::ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube, double maxratio)
|
|---|
| 59 | : arrep_(arrep) , tbeam_(tbeam), skycube_(skycube), maxratio_(maxratio)
|
|---|
| 60 | {
|
|---|
| 61 | double dxdeg = ThetaSizeDegre/(double)NTheta;
|
|---|
| 62 | double dydeg = PhiSizeDegre/(double)NPhi;
|
|---|
| 63 | dx_ = DegreeToRadian(dxdeg);
|
|---|
| 64 | dy_ = DegreeToRadian(dydeg);
|
|---|
| 65 | freq0_ = Freq0MHz;
|
|---|
| 66 | dfreq_ = FreqSizeMHz/(double)NFreq;
|
|---|
| 67 | cout << " ForegroundCleaner: " << " dx=" << dxdeg << " dy=" << dydeg << " degres ( dx_rad=" << dx_ << " dy_rad=" << dy_ << ")"
|
|---|
| 68 | << " Freq0=" << freq0_ << " deltaFreq=" << dfreq_ << " MHz" << endl;
|
|---|
| 69 | skycube.Show();
|
|---|
| 70 | }
|
|---|
| 71 |
|
|---|
| 72 | /* --Methode-- */
|
|---|
| 73 | void ForegroundCleaner::BeamCorrections()
|
|---|
| 74 | {
|
|---|
| 75 | BeamEffect beam(arrep_);
|
|---|
| 76 | beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_, maxratio_);
|
|---|
| 77 | cout << " ForegroundCleaner::BeamCorrections() done " << endl;
|
|---|
| 78 | }
|
|---|
| 79 |
|
|---|
| 80 | /* --Methode-- */
|
|---|
| 81 | int ForegroundCleaner::FixMeanXYTemp(double T0, double alpha)
|
|---|
| 82 | {
|
|---|
| 83 | cout << "ForegroundCleaner::FixMeanXYTemp(T0=" << T0 << ",alpha=" << alpha << ")" << endl;
|
|---|
| 84 | double lnf0=log(freq0_);
|
|---|
| 85 | sa_size_t modprt=skycube_.SizeZ()/12;
|
|---|
| 86 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 87 | double lnf=log((double)k*dfreq_+freq0_);
|
|---|
| 88 | double fittedtemp = T0*exp(alpha*(lnf-lnf0));
|
|---|
| 89 | TArray<TF> slice = skycube_(Range::all(), Range::all(), Range(k));
|
|---|
| 90 | TF deltatemp = (TF)(fittedtemp-(double)Mean(slice));
|
|---|
| 91 | slice += deltatemp;
|
|---|
| 92 | if (k%modprt == 0)
|
|---|
| 93 | cout << "FixMeanXYTemp[k=" << k << " MeanXYTemp=" << fittedtemp-deltatemp << " -> " << fittedtemp
|
|---|
| 94 | << " (DeltaTemp=" << deltatemp << ")" << endl;
|
|---|
| 95 | }
|
|---|
| 96 | cout << "ForegroundCleaner::FixMeanXYTemp done" << endl;
|
|---|
| 97 | return 0;
|
|---|
| 98 | }
|
|---|
| 99 |
|
|---|
| 100 | /* --Methode-- */
|
|---|
| 101 | int ForegroundCleaner::CleanNegatives(TF seuil)
|
|---|
| 102 | {
|
|---|
| 103 | sa_size_t nneg = 0.;
|
|---|
| 104 | for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++)
|
|---|
| 105 | for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
|
|---|
| 106 | for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++)
|
|---|
| 107 | if (skycube_(kx, ky, kz) < seuil) {
|
|---|
| 108 | nneg++; skycube_(kx, ky, kz)=seuil;
|
|---|
| 109 | }
|
|---|
| 110 | cout << " ForegroundCleaner::CleanNegatives " << nneg << " sky-pixels <" << seuil << " changed to" << seuil << endl;
|
|---|
| 111 | return (int)nneg;
|
|---|
| 112 | }
|
|---|
| 113 |
|
|---|
| 114 | /* --Methode-- */
|
|---|
| 115 | int ForegroundCleaner::CleanPointSources(double nsigmas)
|
|---|
| 116 | {
|
|---|
| 117 | Timer tm("ForegroundCleaner::CleanPointSources");
|
|---|
| 118 | TArray< TF > sky2d(skycube_.SizeX(), skycube_.SizeY());
|
|---|
| 119 | for(sa_size_t ky=0; ky<sky2d.SizeY(); ky++)
|
|---|
| 120 | for(sa_size_t kx=0; kx<sky2d.SizeX(); kx++)
|
|---|
| 121 | sky2d(kx, ky) = skycube_(Range(kx), Range(ky), Range::all()).Sum();
|
|---|
| 122 |
|
|---|
| 123 |
|
|---|
| 124 | double mean, sigma;
|
|---|
| 125 |
|
|---|
| 126 |
|
|---|
| 127 | TArray< TF > amz(1,1,skycube_.SizeZ()), asz(1,1,skycube_.SizeZ());
|
|---|
| 128 | for(sa_size_t kz=0; kz<skycube_.SizeZ(); kz++) {
|
|---|
| 129 | TArray< TF > slice = skycube_(Range::all() , Range::all(), Range(kz));
|
|---|
| 130 | MeanSigma(slice, mean, sigma);
|
|---|
| 131 | amz(0,0,kz)=mean; asz(0,0,kz)=sigma;
|
|---|
| 132 | }
|
|---|
| 133 |
|
|---|
| 134 | MeanSigma(sky2d, mean, sigma);
|
|---|
| 135 | cout << " ForegroundCleaner::CleanPointSources 2D Sky projection, mean=" << mean << " sigma=" << sigma << endl;
|
|---|
| 136 |
|
|---|
| 137 | TF seuil = (TF)(mean+nsigmas*sigma);
|
|---|
| 138 |
|
|---|
| 139 | sa_size_t srccnt=0;
|
|---|
| 140 | for(sa_size_t ky=0; ky<skycube_.SizeY(); ky++)
|
|---|
| 141 | for(sa_size_t kx=0; kx<skycube_.SizeX(); kx++) {
|
|---|
| 142 | if (sky2d(kx,ky)>seuil) {
|
|---|
| 143 | srccnt++;
|
|---|
| 144 | skycube_(Range(kx), Range(ky), Range::all()) = amz;
|
|---|
| 145 | }
|
|---|
| 146 | }
|
|---|
| 147 |
|
|---|
| 148 | cout << " Cleaned NSrc= " << srccnt << " 2D source/pixels (TotNPix=" << sky2d.Size()
|
|---|
| 149 | << ")-> " << 100.*srccnt/sky2d.Size() << "% with S>" << seuil << " NSigmas=" << nsigmas << endl;
|
|---|
| 150 | return (int)srccnt;
|
|---|
| 151 | }
|
|---|
| 152 |
|
|---|
| 153 |
|
|---|
| 154 | /* --Methode-- */
|
|---|
| 155 | TArray< TF > ForegroundCleaner::extractLSSCubeP1(TArray< TF >& synctemp, TArray< TF >& specidx)
|
|---|
| 156 | {
|
|---|
| 157 | Timer tm("ForegroundCleaner::extractLSSCubeP1");
|
|---|
| 158 | // Inputs : maplss, mapsyc, freq0, dfreq
|
|---|
| 159 | // Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
|
|---|
| 160 | // Return_Array : foreground subtracted LSS signal
|
|---|
| 161 | sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
|
|---|
| 162 | synctemp.SetSize(2, sz);
|
|---|
| 163 | specidx.SetSize(2, sz);
|
|---|
| 164 | TArray<r_4> omap;
|
|---|
| 165 | omap.SetSize(skycube_, true);
|
|---|
| 166 | Vector vlnf(skycube_.SizeZ());
|
|---|
| 167 | int nprt = 0;
|
|---|
| 168 | // double freq0 : Frequence premier index en k (MHz)
|
|---|
| 169 | // double dfreq : // largeur en frequence de chaque plan (Mhz)
|
|---|
| 170 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
|
|---|
| 171 | vlnf(k)=log((double)k*dfreq_+freq0_);
|
|---|
| 172 |
|
|---|
| 173 | sa_size_t nbinfini=0;
|
|---|
| 174 | sa_size_t nbbad=0;
|
|---|
| 175 | sa_size_t imodprt=skycube_.SizeX()/6;
|
|---|
| 176 | sa_size_t jmodprt=skycube_.SizeY()/6;
|
|---|
| 177 | for(sa_size_t i=0; i<skycube_.SizeX(); i++)
|
|---|
| 178 | for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
|
|---|
| 179 | r_8 s1, sx, sx2, sy, sxy;
|
|---|
| 180 | s1=sx=sx2=sy=sxy=0.;
|
|---|
| 181 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 182 | double lnf = vlnf(k);
|
|---|
| 183 | double ttot=(r_8)(skycube_(i,j,k));
|
|---|
| 184 | if (ttot < 1.e-5) continue;
|
|---|
| 185 | double lntt=log(ttot);
|
|---|
| 186 | s1+=1.; sx+=lnf; sx2+=(lnf*lnf);
|
|---|
| 187 | sy+=lntt; sxy+=(lnf*lntt);
|
|---|
| 188 | }
|
|---|
| 189 | double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
|
|---|
| 190 | double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
|
|---|
| 191 | double T0 = exp(beta+alpha*vlnf(0));
|
|---|
| 192 |
|
|---|
| 193 | bool fgnan = false;
|
|---|
| 194 | if (!isfinite(alpha)||(!isfinite(beta))) {
|
|---|
| 195 | cout << "extractLSSCubeP1[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
|
|---|
| 196 | << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 197 | alpha=beta=-999.;
|
|---|
| 198 | fgnan = true; nbinfini++;
|
|---|
| 199 | }
|
|---|
| 200 | else {
|
|---|
| 201 | double axp1 = beta+alpha*vlnf(0);
|
|---|
| 202 | double axp2 = beta+alpha*vlnf(vlnf.Size()-1);
|
|---|
| 203 |
|
|---|
| 204 | if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
|
|---|
| 205 | cout << "extractLSSCubeP1[" << i << "," << j << "] BAD alpha=" << alpha
|
|---|
| 206 | << " beta=" << beta << " T0=" << T0 << " - (mapsync="
|
|---|
| 207 | << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 208 | fgnan = true; nbbad++;
|
|---|
| 209 | }
|
|---|
| 210 | }
|
|---|
| 211 | if ((i%imodprt==0)&&(j%jmodprt==0))
|
|---|
| 212 | cout << "extractLSSCubeP1[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha
|
|---|
| 213 | << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 214 | synctemp(i,j) = T0;
|
|---|
| 215 | specidx(i,j) = alpha;
|
|---|
| 216 | if (fgnan) {
|
|---|
| 217 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
|
|---|
| 218 | }
|
|---|
| 219 | else {
|
|---|
| 220 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 221 | r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
|
|---|
| 222 | omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
|
|---|
| 223 | }
|
|---|
| 224 | }
|
|---|
| 225 | }
|
|---|
| 226 | cout << " ForegroundCleaner::extractLSSCubeP1() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
|
|---|
| 227 | return omap;
|
|---|
| 228 | }
|
|---|
| 229 |
|
|---|
| 230 | /*
|
|---|
| 231 | static inline val_polyn2(double alpha, double beta, double gamma, double x)
|
|---|
| 232 | {
|
|---|
| 233 | return (beta+alpha*x+gamma*x*x);
|
|---|
| 234 | }
|
|---|
| 235 | */
|
|---|
| 236 |
|
|---|
| 237 | /* --Methode-- */
|
|---|
| 238 | TArray< TF > ForegroundCleaner::extractLSSCubeP2(TArray< TF >& synctemp, TArray< TF >& specidx)
|
|---|
| 239 | {
|
|---|
| 240 | Timer tm("ForegroundCleaner::extractLSSCubeP2");
|
|---|
| 241 | // Inputs : maplss, mapsyc, freq0, dfreq
|
|---|
| 242 | // Outputs : synctemp, specidx (reconstructed foreground temperature and spectral index
|
|---|
| 243 | // Return_Array : foreground subtracted LSS signal
|
|---|
| 244 | sa_size_t sz[5]; sz[0]=skycube_.SizeX(); sz[1]=skycube_.SizeY();
|
|---|
| 245 | synctemp.SetSize(2, sz);
|
|---|
| 246 | specidx.SetSize(2, sz);
|
|---|
| 247 | TArray<r_4> omap;
|
|---|
| 248 | omap.SetSize(skycube_, true);
|
|---|
| 249 | Vector vlnf(skycube_.SizeZ());
|
|---|
| 250 | Vector vlnT(skycube_.SizeZ());
|
|---|
| 251 | int nprt = 0;
|
|---|
| 252 | // double freq0 : Frequence premier index en k (MHz)
|
|---|
| 253 | // double dfreq : // largeur en frequence de chaque plan (Mhz)
|
|---|
| 254 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++)
|
|---|
| 255 | vlnf(k)=log((double)k*dfreq_+freq0_);
|
|---|
| 256 | vlnf -= vlnf(0);
|
|---|
| 257 | // cout << " DBG*extractLSSCubeP2 vlnf(0)=" << vlnf(0) << " vlnf(1)=" << vlnf(1)
|
|---|
| 258 | // << "vlnf(last)=" << vlnf(vlnf.Size()-1) << endl;
|
|---|
| 259 | sa_size_t nbinfini=0;
|
|---|
| 260 | sa_size_t nbbad=0;
|
|---|
| 261 | sa_size_t imodprt=skycube_.SizeX()/6;
|
|---|
| 262 | sa_size_t jmodprt=skycube_.SizeY()/6;
|
|---|
| 263 | for(sa_size_t i=0; i<skycube_.SizeX(); i++)
|
|---|
| 264 | for(sa_size_t j=0; j<skycube_.SizeY(); j++) {
|
|---|
| 265 | vlnT = -12.;
|
|---|
| 266 | Poly polyn;
|
|---|
| 267 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 268 | double lnf = vlnf(k);
|
|---|
| 269 | double ttot=(r_8)(skycube_(i,j,k));
|
|---|
| 270 | if (ttot < 1.e-5) continue;
|
|---|
| 271 | vlnT(k)=log(ttot);
|
|---|
| 272 | }
|
|---|
| 273 | polyn.Fit(vlnf,vlnT,2);
|
|---|
| 274 | double beta = polyn[0];
|
|---|
| 275 | double alpha = polyn[1];
|
|---|
| 276 | double gamma = polyn[2];
|
|---|
| 277 | double T0 = exp(polyn(vlnf(0))); // exp( val_polyn2(alpha, beta, gamma, vlnf(0)) );
|
|---|
| 278 |
|
|---|
| 279 | bool fgnan = false;
|
|---|
| 280 | if (!isfinite(alpha)||(!isfinite(beta))||(!isfinite(gamma))) {
|
|---|
| 281 | cout << "extractLSSCubeP2[" << i << "," << j << "]/ Not finite alpha, beta - (mapsync="
|
|---|
| 282 | << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 283 | alpha=beta=gamma=-999.;
|
|---|
| 284 | fgnan = true; nbinfini++;
|
|---|
| 285 | }
|
|---|
| 286 | else {
|
|---|
| 287 | double axp1 = polyn(vlnf(0));
|
|---|
| 288 | double axp2 = polyn(vlnf(vlnf.Size()-1));
|
|---|
| 289 |
|
|---|
| 290 | if ((axp1<-70.)||(axp1>70.)||(axp2<-70.)||(axp2>70.)) {
|
|---|
| 291 | cout << "extractLSSCubeP2[" << i << "," << j << "] BAD alpha=" << alpha
|
|---|
| 292 | // << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " axp1=" << axp1 << " axp2=" << axp2
|
|---|
| 293 | << " beta=" << beta << " gamma=" << gamma << " T0=" << T0 << " - (mapsync=" << skycube_(i,j,0)
|
|---|
| 294 | << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 295 | fgnan = true; nbbad++;
|
|---|
| 296 | }
|
|---|
| 297 | }
|
|---|
| 298 | if ((i%imodprt==0)&&(j%jmodprt==0))
|
|---|
| 299 | cout << "extractLSSCubeP2[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha << " gamma=" << gamma
|
|---|
| 300 | << " (mapsync=" << skycube_(i,j,0) << " ... " << skycube_(i,j,skycube_.SizeZ()-1) << ")" << endl;
|
|---|
| 301 | synctemp(i,j) = T0;
|
|---|
| 302 | specidx(i,j) = alpha;
|
|---|
| 303 | if (fgnan) {
|
|---|
| 304 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) omap(i,j,k) = 0.;
|
|---|
| 305 | }
|
|---|
| 306 | else {
|
|---|
| 307 | for(sa_size_t k=0; k<skycube_.SizeZ(); k++) {
|
|---|
| 308 | r_4 fittedtemp = (r_4)( exp(polyn(vlnf(k))) );
|
|---|
| 309 | omap(i,j,k) = skycube_(i,j,k)-fittedtemp;
|
|---|
| 310 | }
|
|---|
| 311 | }
|
|---|
| 312 | }
|
|---|
| 313 | cout << " ForegroundCleaner::extractLSSCubeP2() - NbNan alpha/beta=" << nbinfini << " NbBAD =" << nbbad << endl;
|
|---|
| 314 | return omap;
|
|---|
| 315 | }
|
|---|
| 316 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
