00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046 #ifndef BK_BASELINE_H
00047 #define BK_BASELINE_H
00048 using namespace std;
00049
00050 #include <iostream>
00051 #include <cstdlib>
00052
00053 #include "constants.hh"
00054 #include "bk_messages.hh"
00055 #include "slcimage.hh"
00056 #include "productinfo.hh"
00057 #include "orbitbk.hh"
00058 #include "conversion.hh"
00059 #include "utilities.hh"
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 class BASELINE
00076 {
00077
00078 private:
00079 bool initialized;
00080 real8 master_wavelength;
00081 real8 nearrange;
00082 real8 drange_dp;
00083 real8 orbit_convergence;
00084 real8 orbit_heading;
00085 real8 H_min;
00086 real8 H_max;
00087
00088
00089
00090
00091 int32 N_coeffs;
00092
00093 matrix<real8> BPERP_cf;
00094 matrix<real8> BPAR_cf;
00095 matrix<real8> THETA_cf;
00096 matrix<real8> THETA_INC_cf;
00097
00098
00099
00100 BASELINE(const BASELINE& A)
00101 {};
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 real8 polyval(
00112 const matrix<real8> &C,
00113 const real8 line,
00114 const real8 pixel,
00115 const real8 height) const
00116 {
00117 #ifdef __DEBUG
00118 TRACE_FUNCTION("baseline::polyval (BK 05-Mar-2005)")
00119 #endif
00120 if (C.size()!=10) throw("error");
00121 return
00122 C(0,0) +
00123 C(1,0)*line + C(2,0)*pixel + C(3,0)*height +
00124 C(4,0)*line*pixel + C(5,0)*line*height + C(6,0)*pixel*height +
00125 C(7,0)*sqr(line) + C(8,0)*sqr(pixel) + C(9,0)*sqr(height);
00126 };
00127
00128
00129
00130
00131 void BBparBperpTheta(real8 &B, real8 &Bpar, real8 &Bperp, real8 &theta,
00132 const cn Master, const cn Point, const cn Slave) const
00133 {
00134 #ifdef __DEBUG
00135 TRACE_FUNCTION("BBparBperp (BK 05-Mar-2005)")
00136 #endif
00137 B = Master.dist(Slave);
00138 const real8 range1 = Master.dist(Point);
00139 const real8 range2 = Slave.dist(Point);
00140 Bpar = range1-range2;
00141 const cn r1 = Master.min(Point);
00142 const cn r2 = Slave.min(Point);
00143 theta = Master.angle(r1);
00144 Bperp = sqr(B)-sqr(Bpar);
00145 if (Bperp < 0.0) Bperp=0.0;
00146 else Bperp = (theta > Master.angle(r2)) ?
00147 sqrt(Bperp) : -sqrt(Bperp);
00148 };
00149
00150
00151
00152
00153
00154
00155
00156
00157 real8 IncidenceAngle(const cn Master, const cn Point) const
00158 {
00159 #ifdef __DEBUG
00160 TRACE_FUNCTION("IncidenceAngle (BK 05-Mar-2005)")
00161 #endif
00162 const cn r1 = Master.min(Point);
00163 const real8 inc = Point.angle(r1);
00164 return inc;
00165 };
00166
00167
00168
00169 public:
00170
00171 BASELINE()
00172 {
00173 TRACE_FUNCTION("BASELINE() (BK 06-Mar-2005)");
00174 initialized = false;
00175 N_coeffs = 10;
00176 H_min = 0.0;
00177 H_max = 5000.0;
00178 master_wavelength = 0.0;
00179 orbit_convergence = 0.0;
00180 orbit_heading = 0.0;
00181 }
00182
00183
00184 ~BASELINE()
00185 {
00186 TRACE_FUNCTION("~BASELINE() (BK 06-Mar-2005)");
00187 ;
00188 }
00189
00190
00191
00192 void model_parameters(
00193 const slcimage &master,
00194 const slcimage &slave,
00195 orbit &masterorbit,
00196 orbit &slaveorbit,
00197
00198 const input_ell &ellips)
00199 {
00200 TRACE_FUNCTION("model_parameters (BK 06-Mar-2005)");
00201 if (masterorbit.is_initialized() == false)
00202 {
00203 DEBUG.print("Baseline cannot be computed, master orbit not initialized.");
00204 return;
00205 }
00206 if (slaveorbit.is_initialized() == false)
00207 {
00208 DEBUG.print("Baseline cannot be computed, slave orbit not initialized.");
00209 return;
00210 }
00211
00212
00213 initialized = true;
00214 master_wavelength = master.wavelength;
00215
00216 nearrange = master.pix2range(1.0);
00217 drange_dp = master.pix2range(2.0)-master.pix2range(1.0);
00218 nearrange -= drange_dp;
00219
00220
00221
00222 register int32 cnt = 0;
00223 const int N_pointsL = 10;
00224 const int N_pointsP = 10;
00225 const int N_heights = 4;
00226 const real8 deltapixels = master.currentwindow.pixels() / N_pointsP;
00227 const real8 deltalines = master.currentwindow.lines() / N_pointsL;
00228 const real8 deltaheight = (H_max-H_min) / N_heights;
00229
00230
00231 matrix<real8> BPERP(N_pointsL*N_pointsP*N_heights,1);
00232 matrix<real8> BPAR(N_pointsL*N_pointsP*N_heights,1);
00233 matrix<real8> THETA(N_pointsL*N_pointsP*N_heights,1);
00234 matrix<real8> THETA_INC(N_pointsL*N_pointsP*N_heights,1);
00235 matrix<real8> AMATRIX(N_pointsL*N_pointsP*N_heights,N_coeffs);
00236
00237
00238 for (register int32 k=0; k<N_heights; ++k)
00239 {
00240 const real8 HEIGHT = H_min + k*deltaheight;
00241 input_ell ELLIPS;
00242 ELLIPS.a = ellips.a + HEIGHT;
00243 ELLIPS.b = ellips.b + HEIGHT;
00244 ELLIPS.ecc1st_sqr();
00245 ELLIPS.ecc2nd_sqr();
00246
00247 for (register int32 i=0; i<N_pointsL; ++i)
00248 {
00249 const real8 line = master.currentwindow.linelo + i*deltalines;
00250 cn P;
00251 real8 s_tazi;
00252 real8 s_trange;
00253 const int32 MAXITER = 10;
00254 const real8 CRITERPOS = 1e-6;
00255 const real8 CRITERTIM = 1e-10;
00256
00257 const real8 m_tazi = master.line2ta(line);
00258
00259 const cn M = masterorbit.getxyz(m_tazi);
00260
00261 for (register int32 j=0; j<N_pointsP; ++j)
00262 {
00263 const real8 pixel = master.currentwindow.pixlo + j*deltapixels;
00264
00265
00266 lp2xyz(line,pixel,ELLIPS,master,masterorbit,P,MAXITER,CRITERPOS);
00267
00268 xyz2t(s_tazi,s_trange,slave,slaveorbit,P,MAXITER,CRITERTIM);
00269
00270 const cn S = slaveorbit.getxyz(s_tazi);
00271
00272 const cn Mdot = masterorbit.getxyzdot(m_tazi);
00273 const cn Sdot = slaveorbit.getxyzdot(s_tazi);
00274 const real8 angleorbits = Mdot.angle(Sdot);
00275 DEBUG << "Angle between orbits master-slave (at l,p="
00276 << line << "," << pixel << ") = "
00277 << rad2deg(angleorbits) << " [deg]";
00278 DEBUG.print();
00279 orbit_convergence = angleorbits;
00280
00281
00282
00283
00284
00285
00286 real8 B,Bpar,Bperp,theta;
00287 BBparBperpTheta(B,Bpar,Bperp,theta,M,P,S);
00288 const real8 theta_inc = IncidenceAngle(M,P);
00289
00290
00291 BPERP(cnt,0) = Bperp;
00292 BPAR(cnt,0) = Bpar;
00293 THETA(cnt,0) = theta;
00294 THETA_INC(cnt,0) = theta_inc;
00295
00296
00297
00298
00299 AMATRIX(cnt,0) = 1.0;
00300 AMATRIX(cnt,1) = line;
00301 AMATRIX(cnt,2) = pixel;
00302 AMATRIX(cnt,3) = HEIGHT;
00303 AMATRIX(cnt,4) = line*pixel;
00304 AMATRIX(cnt,5) = line*HEIGHT;
00305 AMATRIX(cnt,6) = pixel*HEIGHT;
00306 AMATRIX(cnt,7) = sqr(line);
00307 AMATRIX(cnt,8) = sqr(pixel);
00308 AMATRIX(cnt,9) = sqr(HEIGHT);
00309 cnt++;
00310
00311 const real8 alpha = (Bpar==0 && Bperp==0) ?
00312 NaN :
00313 theta - atan2(Bpar,Bperp);
00314
00315 const real8 Bh = B * cos(alpha);
00316 const real8 Bv = B * sin(alpha);
00317
00318 const real8 hambiguity = (Bperp==0) ?
00319 Inf :
00320 -master.wavelength*(M.min(P)).norm()*sin(theta)/(2.0*Bperp);
00321
00322
00323
00324 DEBUG << "The baseline parameters for (l,p,h) = "
00325 << line << ", " << pixel << ", " << HEIGHT;
00326 DEBUG.print();
00327 DEBUG << "\talpha (deg), baseline: \t" << rad2deg(alpha) << " \t" << B;
00328 DEBUG.print();
00329 DEBUG << "\tBpar, Bperp: \t" << Bpar << " \t" << Bperp;
00330 DEBUG.print();
00331 DEBUG << "\tBh, Bv: \t" << Bh << " \t" << Bv;
00332 DEBUG.print();
00333 DEBUG << "\tHeight ambiguity: \t" << hambiguity;
00334 DEBUG.print();
00335 DEBUG << "\ttheta (deg): \t" << rad2deg(theta);
00336 DEBUG.print();
00337 DEBUG << "\ttheta_inc (deg): \t" << rad2deg(theta_inc);
00338 DEBUG.print();
00339 DEBUG.precision(10);
00340 DEBUG.width(11);
00341 DEBUG.rewind();
00342 DEBUG << "\tM (x,y,z) = " << M.x << ", " << M.y << ", " << M.z;
00343 DEBUG.print();
00344 DEBUG << "\tS (x,y,z) = " << S.x << ", " << S.y << ", " << S.z;
00345 DEBUG.print();
00346 DEBUG << "\tP (x,y,z) = " << P.x << ", " << P.y << ", " << P.z;
00347 DEBUG.print();
00348 DEBUG.reset();
00349 }
00350 }
00351 }
00352
00353
00354 matrix<real8> N = matTxmat(AMATRIX,AMATRIX);
00355 matrix<real8> rhsBperp = matTxmat(AMATRIX,BPERP);
00356 matrix<real8> rhsBpar = matTxmat(AMATRIX,BPAR);
00357 matrix<real8> rhsT = matTxmat(AMATRIX,THETA);
00358 matrix<real8> rhsT_INC = matTxmat(AMATRIX,THETA_INC);
00359 matrix<real8> Qx_hat = N;
00360 choles(Qx_hat);
00361 solvechol(Qx_hat,rhsBperp);
00362 solvechol(Qx_hat,rhsBpar);
00363 solvechol(Qx_hat,rhsT);
00364 solvechol(Qx_hat,rhsT_INC);
00365 invertchol(Qx_hat);
00366
00367 ;
00368 BPERP_cf = rhsBperp;
00369 BPAR_cf = rhsBpar;
00370 THETA_cf = rhsT;
00371 THETA_INC_cf = rhsT_INC;
00372
00373
00374 for (int32 i=0; i<Qx_hat.lines(); i++)
00375 for (int32 j=0; j<i; j++)
00376 Qx_hat(j,i) = Qx_hat(i,j);
00377 const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00378 DEBUG << "baseline: max(abs(N*inv(N)-I)) = " << maxdev;
00379 DEBUG.print();
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396 if (maxdev > .01)
00397 {
00398 WARNING << "baseline: max. deviation N*inv(N) from unity = "
00399 << maxdev << ". This is larger than .01: do not use this.";
00400 WARNING.print();
00401 }
00402 else if (maxdev > .001)
00403 {
00404 WARNING << "baseline: max. deviation N*inv(N) from unity = "
00405 << maxdev << ". This is between 0.01 and 0.001 (do not use it)";
00406 WARNING.print();
00407 }
00408
00409
00410 matrix<real8> y_hatBperp = AMATRIX * rhsBperp;
00411 matrix<real8> e_hatBperp = BPERP - y_hatBperp;
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422 DEBUG.print("--------------------");
00423 DEBUG.print("Result of modeling: Bperp(l,p) = a000 + a100*l + a010*p + a001*h + ");
00424 DEBUG.print(" a110*l*p + a101*l*h + a011*p*h + a200*l^2 + a020*p^2 + a002*h^2");
00425 DEBUG.print("l,p,h in unnormalized, absolute, coordinates (1:N).");
00426 DEBUG << "Bperp_a000 = " << rhsBperp(0,0); DEBUG.print();
00427 DEBUG << "Bperp_a100 = " << rhsBperp(1,0); DEBUG.print();
00428 DEBUG << "Bperp_a010 = " << rhsBperp(2,0); DEBUG.print();
00429 DEBUG << "Bperp_a001 = " << rhsBperp(3,0); DEBUG.print();
00430 DEBUG << "Bperp_a110 = " << rhsBperp(4,0); DEBUG.print();
00431 DEBUG << "Bperp_a101 = " << rhsBperp(5,0); DEBUG.print();
00432 DEBUG << "Bperp_a011 = " << rhsBperp(6,0); DEBUG.print();
00433 DEBUG << "Bperp_a200 = " << rhsBperp(7,0); DEBUG.print();
00434 DEBUG << "Bperp_a020 = " << rhsBperp(8,0); DEBUG.print();
00435 DEBUG << "Bperp_a002 = " << rhsBperp(9,0); DEBUG.print();
00436 real8 maxerr = max(abs(e_hatBperp));
00437 if (maxerr > 2.00)
00438 {
00439 WARNING << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
00440 WARNING.print();
00441 }
00442 else
00443 {
00444 INFO << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
00445 INFO.print();
00446 }
00447 DEBUG.print();
00448 DEBUG.print("--------------------");
00449 DEBUG.print("Range: r(p) = r0 + dr*p");
00450 DEBUG.print("l and p in unnormalized, absolute, coordinates (1:N).");
00451 const real8 range1 = master.pix2range(1.0);
00452 const real8 range5000 = master.pix2range(5000.0);
00453 const real8 drange = (range5000-range1)/5000;
00454 DEBUG << "range = " << range1-drange << " + " << drange << "*p";
00455 DEBUG.print();
00456
00457 };
00458
00459
00460 inline real8 get_range(const real8 pixel) const
00461 {
00462 return nearrange + drange_dp*pixel;
00463 };
00464
00465
00466
00467
00468
00469
00470
00471 inline real8 get_bperp(const real8 line, const real8 pixel, const real8 height=0.0) const
00472 {
00473 return polyval(BPERP_cf, line, pixel, height);
00474 };
00475
00476 inline real8 get_bpar(const real8 line, const real8 pixel, const real8 height=0.0) const
00477 {
00478 return polyval(BPAR_cf, line, pixel, height);
00479 };
00480
00481 inline real8 get_theta(const real8 line, const real8 pixel, const real8 height=0.0) const
00482 {
00483 return polyval(THETA_cf, line, pixel, height);
00484 };
00485
00486 inline real8 get_theta_inc(const real8 line, const real8 pixel, const real8 height=0.0) const
00487 {
00488 return polyval(THETA_INC_cf, line, pixel, height);
00489 };
00490
00491
00492
00493 inline real8 get_b(const real8 line, const real8 pixel, const real8 height=0.0) const
00494 {
00495 return sqrt(sqr(get_bpar(line,pixel,height))+sqr(get_bperp(line,pixel,height)));
00496 };
00497
00498
00499 inline real8 get_alpha(const real8 line, const real8 pixel, const real8 height=0.0) const
00500 {
00501 const real8 Bperp = get_bperp(line,pixel,height);
00502 const real8 Bpar = get_bpar(line,pixel,height);
00503 const real8 B = get_b(line,pixel,height);
00504 const real8 theta = get_theta(line,pixel,height);
00505 const real8 alpha = (Bpar==0 && Bperp==0) ?
00506 NaN :
00507 theta-atan2(Bpar,Bperp);
00508 return alpha;
00509 };
00510
00511
00512 inline real8 get_bhor(const real8 line, const real8 pixel, const real8 height=0.0) const
00513 {
00514 const real8 Bperp = get_bperp(line,pixel,height);
00515 const real8 Bpar = get_bpar(line,pixel,height);
00516 const real8 B = get_b(line,pixel,height);
00517 const real8 theta = get_theta(line,pixel,height);
00518 const real8 alpha = get_alpha(line,pixel,height);
00519 return B*cos(alpha);
00520 };
00521
00522
00523 inline real8 get_bvert(const real8 line, const real8 pixel, const real8 height=0.0) const
00524 {
00525 const real8 Bperp = get_bperp(line,pixel,height);
00526 const real8 Bpar = get_bpar(line,pixel,height);
00527 const real8 B = get_b(line,pixel,height);
00528 const real8 theta = get_theta(line,pixel,height);
00529 const real8 alpha = get_alpha(line,pixel,height);
00530 return B*sin(alpha);
00531 };
00532
00533
00534 inline real8 get_hamb(const real8 line, const real8 pixel, const real8 height=0.0) const
00535 {
00536
00537 const real8 theta_inc = get_theta_inc(line,pixel,height);
00538 const real8 Bperp = get_bperp(line,pixel,height);
00539 const real8 range_MP = get_range(pixel);
00540 const real8 h_amb = (Bperp==0) ?
00541 Inf :
00542 -master_wavelength*range_MP*sin(theta_inc)/(2.0*Bperp);
00543
00544 return h_amb;
00545 };
00546
00547
00548 inline real8 get_orb_conv(const real8 line, const real8 pixel, const real8 height=0.0) const
00549 {
00550
00551 return orbit_convergence;
00552 };
00553
00554
00555
00556 void dump(const real8 line, const real8 pixel, const real8 height=0.0)
00557 {
00558 if (initialized==false)
00559 {
00560 DEBUG.print("Exiting dumpbaseline, not initialized.");
00561 return;
00562 }
00563
00564 const real8 Bperp = get_bperp(line,pixel,height);
00565 const real8 Bpar = get_bpar(line,pixel,height);
00566 const real8 theta = get_theta(line,pixel,height);
00567 const real8 theta_inc = get_theta_inc(line,pixel,height);
00568
00569 const real8 B = get_b(line,pixel,height);
00570 const real8 alpha = get_alpha(line,pixel,height);
00571 const real8 Bh = get_bhor(line,pixel,height);
00572 const real8 Bv = get_bvert(line,pixel,height);
00573 const real8 h_amb = get_hamb(line,pixel,height);
00574
00575
00576 INFO << "The baseline parameters for (l,p,h) = "
00577 << line << ", " << pixel << ", " << height;
00578 INFO.print();
00579 INFO << "\tBpar, Bperp: \t" << Bpar << " \t" << Bperp;
00580 INFO.print();
00581 DEBUG << "\tB, alpha (deg): \t" << B << " \t" << rad2deg(alpha);
00582 DEBUG.print();
00583 DEBUG << "\tBh, Bv: \t" << Bh << " \t" << Bv;
00584 DEBUG.print();
00585 INFO << "\tHeight ambiguity: \t" << h_amb;
00586 INFO.print();
00587 INFO << "\tLook angle (deg): \t" << rad2deg(theta);
00588 INFO.print();
00589 DEBUG << "\tIncidence angle (deg): \t" << rad2deg(theta_inc);
00590 DEBUG.print();
00591 }
00592
00593 };
00594
00595 #endif // BK_BASELINE_H
00596