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 #include "constants.hh"
00045 #include "matrixbk.hh"
00046 #include "slcimage.hh"
00047 #include "productinfo.hh"
00048 #include "geocode.hh"
00049 #include "utilities.hh"
00050 #include "ioroutines.hh"
00051 #include "coregistration.hh"
00052 #include "conversion.hh"
00053 #include "exceptions.hh"
00054
00055 #include <strstream>
00056 #include <cctype>
00057 #include <algorithm>
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 void slant2hschwabisch(
00085 const input_gen &generalinput,
00086 const input_slant2h &slant2hinput,
00087 const input_ell &ellips,
00088 const slcimage &master,
00089 const slcimage &slave,
00090 const productinfo &unwrappedinterf,
00091 orbit &masterorbit,
00092 orbit &slaveorbit)
00093 {
00094 TRACE_FUNCTION("slant2hschwabisch (BK 02-Jun-1999)")
00095
00096 const int32 MAXITER = 10;
00097 const real8 CRITERPOS = 1e-6;
00098 const real8 CRITERTIM = 1e-10;
00099 const real8 m_minpi4cdivlambda = (-4.*PI*SOL)/master.wavelength;
00100 const real8 s_minpi4cdivlambda = (-4.*PI*SOL)/slave.wavelength;
00101
00102 const int32 Npoints = slant2hinput.Npoints;
00103 const int32 DEGREE1D = slant2hinput.degree1d;
00104 const int32 DEGREE2D = slant2hinput.degree2d;
00105 const int32 Nheights = slant2hinput.Nheights;
00106
00107 const int32 MAXHEIGHT = 5000;
00108 const int32 TEN = 10;
00109 if (DEGREE1D - 1 > TEN)
00110 {
00111 PRINT_ERROR("panic, programmers problem: increase TEN.")
00112 throw(some_error);
00113 }
00114
00115 const int32 HEIGHTSTEP = MAXHEIGHT / (Nheights - 1);
00116
00117
00118
00119
00120
00121
00122
00123 matrix<real8> PHASE(Npoints,Nheights);
00124
00125
00126
00127
00128 matrix<uint> Position = distributepoints(Npoints,unwrappedinterf.win);
00129
00130 cn pospoint;
00131 input_ell ELLIPS;
00132
00133
00134
00135
00136 PROGRESS.print("S2H: schwabisch: STEP1: compute reference phase for Nheights.");
00137
00138 register int32 numheight;
00139 register int32 i,j,k,l,index;
00140 register int32 alfa;
00141 for (numheight=0; numheight<Nheights; numheight++)
00142 {
00143 const int32 HEIGHT = numheight * HEIGHTSTEP;
00144 ELLIPS.a = ellips.a + HEIGHT;
00145 ELLIPS.b = ellips.b + HEIGHT;
00146 ELLIPS.ecc1st_sqr();
00147 ELLIPS.ecc2nd_sqr();
00148
00149
00150 for (i=0;i<Npoints;i++)
00151 {
00152 const real8 line = Position(i,0);
00153 const real8 pixel = Position(i,1);
00154 const real8 m_trange = master.pix2tr(pixel);
00155
00156
00157 lp2xyz(line,pixel,ELLIPS,master,masterorbit,
00158 pospoint,MAXITER,CRITERPOS);
00159
00160
00161 real8 s_tazi;
00162 real8 s_trange;
00163 xyz2t(s_tazi,s_trange,slave, slaveorbit,
00164 pospoint,MAXITER,CRITERTIM);
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 PHASE(i,numheight) = m_minpi4cdivlambda*m_trange-
00175 s_minpi4cdivlambda*s_trange;
00176 }
00177 }
00178
00179
00180
00181
00182
00183
00184 for (i=0; i<Npoints; ++i)
00185 {
00186 real8 offset = PHASE(i,0);
00187 PHASE(i,0) = 0.;
00188 for (j=1; j<numheight; ++j)
00189 {
00190 PHASE(i,j) -= offset;
00191 }
00192 }
00193
00194
00195
00196
00197 PROGRESS.print("S2H: schwabisch: STEP2: estimate coefficients 1d polynomial.");
00198
00199
00200
00201 matrix<real8> DESIGN(Nheights,DEGREE1D+1);
00202 matrix<real8> ALPHAS(Npoints,DEGREE1D+1);
00203 matrix<real8> HEI(Nheights,1);
00204 for (i=0; i<Nheights; i++)
00205 HEI(i,0) = i*HEIGHTSTEP;
00206
00207
00208 const real8 minphi = min(PHASE);
00209 const real8 maxphi = max(PHASE);
00210 normalize(PHASE,minphi,maxphi);
00211
00212 for (i=0; i<Npoints; i++)
00213 {
00214
00215
00216 for (j=0; j<Nheights; j++)
00217 for (k=0; k<=DEGREE1D; k++)
00218 DESIGN(j,k) = pow(PHASE(i,j),real8(k));
00219
00220
00221
00222 matrix<real8> N = matTxmat(DESIGN,DESIGN);
00223 matrix<real8> rhs = matTxmat(DESIGN,HEI);
00224 matrix<real8> Qx_hat = N;
00225 choles(Qx_hat);
00226 solvechol(Qx_hat,rhs);
00227
00228
00229
00230
00231 invertchol(Qx_hat);
00232 for (k=0; k<Qx_hat.lines(); k++)
00233 for (j=0; j<k; j++)
00234 Qx_hat(j,k) = Qx_hat(k,j);
00235 const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00236 INFO << "s2h schwaebisch: max(abs(N*inv(N)-I)) = " << maxdev;
00237 INFO.print();
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 if (maxdev > .01)
00256 WARNING.print("wrong solution for 1d polynomial? (decrease d1d or nhei)");
00257
00258
00259
00260 for (alfa=0; alfa<=DEGREE1D; alfa++)
00261 ALPHAS(i,alfa) = rhs(alfa,0);
00262 }
00263
00264
00265
00266
00267 PROGRESS.print("S2H: schwabisch: STEP3: estimate coefficients for 2d polynomial.");
00268
00269
00270
00271
00272 const int32 Nunk = Ncoeffs(DEGREE2D);
00273
00274
00275 if (Npoints < Nunk)
00276 {
00277 PRINT_ERROR("slant2hschwabisch: N_observations<N_unknowns (increase S2H_NPOINTS or decrease S2H_DEGREE2D.")
00278 throw(input_error);
00279 }
00280
00281 matrix<real8> A(Npoints,Nunk);
00282
00283
00284
00285
00286 const real8 minL = min(Position.getcolumn(0));
00287 const real8 maxL = max(Position.getcolumn(0));
00288 const real8 minP = min(Position.getcolumn(1));
00289 const real8 maxP = max(Position.getcolumn(1));
00290
00291 for (i=0; i<Npoints; i++)
00292 {
00293
00294 const real8 posL = normalize(real8(Position(i,0)),minL,maxL);
00295 const real8 posP = normalize(real8(Position(i,1)),minP,maxP);
00296
00297 index = 0;
00298 for (j=0; j<=DEGREE2D; j++)
00299 {
00300 for (k=0; k<=j; k++)
00301 {
00302 A(i,index) = pow(posL,real8(j-k)) * pow(posP,real8(k));
00303 index++;
00304 }
00305 }
00306 }
00307
00308
00309 matrix<real8> N = matTxmat(A,A);
00310 matrix<real8> rhs = matTxmat(A,ALPHAS);
00311 matrix<real8> Qx_hat = N;
00312 choles(Qx_hat);
00313
00314
00315
00316
00317
00318 for (i=0; i<rhs.pixels(); ++i)
00319 {
00320 matrix<real8> rhs_alphai = rhs.getcolumn(i);
00321 solvechol(Qx_hat,rhs_alphai);
00322 rhs.setcolumn(i,rhs_alphai);
00323 }
00324
00325
00326
00327 invertchol(Qx_hat);
00328 for (i=0; i<Qx_hat.lines(); i++)
00329 for (j=0; j<i; j++)
00330 Qx_hat(j,i) = Qx_hat(i,j);
00331 const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00332 INFO << "s2h schwaebisch: max(abs(N*inv(N)-I)) = " << maxdev;
00333 INFO.print();
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 if (maxdev > .01)
00352 {
00353 WARNING << "slant2h: possibly wrong solution. deviation from unity AtA*inv(AtA) = "
00354 << maxdev << " >.01";
00355 WARNING.print();
00356 }
00357
00358
00359
00360
00361 PROGRESS.print("S2H: schwabisch: STEP4: compute height for all pixels.");
00362
00363
00364
00365 const real8 multiL = unwrappedinterf.multilookL;
00366 const real8 multiP = unwrappedinterf.multilookP;
00367
00368
00369
00370 const int32 mllines = int32(floor(real8(
00371 unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
00372 const int32 mlpixels = int32(floor(real8(
00373 unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1) / multiP));
00374
00375
00376
00377 const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
00378 (real8(multiL) - 1.) / 2.;
00379 const real8 firstpixel = real8(unwrappedinterf.win.pixlo) +
00380 (real8(multiP) - 1.) / 2.;
00381
00382
00383
00384 matrix<real4> p_axis(mlpixels,1);
00385 for (i=0; i<mlpixels; i++)
00386
00387 p_axis(i,0) = firstpixel + i*multiP;
00388 normalize(p_axis,minP,maxP);
00389
00390 const int32 NUMMAT = 1 + (1 + DEGREE1D);
00391 int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
00392 if (bufferlines > mllines)
00393 bufferlines = mllines;
00394
00395 const int32 FULLBUFFERS = mllines / bufferlines;
00396 const int32 RESTLINES = mllines % bufferlines;
00397 const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
00398
00399
00400 const uint dummy = 999999;
00401 window bufferwin(1, bufferlines, 1, mlpixels);
00402 window offsetbuffer(1,dummy,1,dummy);
00403
00404
00405 ofstream ofile;
00406 openfstream(ofile,slant2hinput.fohei,generalinput.overwrit);
00407 bk_assert(ofile,slant2hinput.fohei,__FILE__,__LINE__);
00408
00409
00410
00411 for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; buffer++)
00412 {
00413
00414
00415 PROGRESS << "SLANT2H: "
00416 << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
00417 PROGRESS.print();
00418
00419
00420 const real8 firstline = veryfirstline + (buffer-1) * bufferlines * multiL;
00421
00422
00423 bufferwin.linelo = 1 + (buffer-1) * bufferlines;
00424 if (buffer == FULLBUFFERS+1)
00425 {
00426 bufferlines = RESTLINES;
00427
00428 }
00429 bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
00430
00431
00432 matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
00433
00434
00435
00436
00437 matrix<real4> l_axis(bufferlines,1);
00438 for (k=0; k<bufferlines; k++)
00439 l_axis(k,0) = firstline + k*multiL;
00440 normalize(l_axis,minL,maxL);
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452 matrix<real4> *pntALPHA[TEN];
00453
00454 for (k=0; k<=DEGREE1D; k++)
00455 {
00456 matrix<real8> beta(Ncoeffs(DEGREE2D),1);
00457 for (l=0; l<Ncoeffs(DEGREE2D); l++)
00458 {
00459 beta(l,0) = rhs(l,k);
00460 }
00461 pntALPHA[k] = new matrix<real4> (l_axis.lines(),p_axis.lines());
00462 (*pntALPHA[k]) = polyval(l_axis, p_axis, beta, DEGREE2D);
00463 }
00464
00465
00466
00467
00468
00469 matrix<real8> coeff_thispoint(DEGREE1D+1,1);
00470 for (register int32 line=0; line<BUFFER.lines(); line++)
00471 {
00472 for (register int32 pixel=0; pixel<BUFFER.pixels(); pixel++)
00473 {
00474
00475
00476 if (BUFFER(line,pixel) != NaN)
00477 {
00478 for (k=0; k<DEGREE1D+1; k++)
00479 {
00480 coeff_thispoint(k,0) = (*pntALPHA[k])[line][pixel];
00481 }
00482 BUFFER(line,pixel) = polyval1d(normalize(
00483 BUFFER(line,pixel),minphi,maxphi),
00484 coeff_thispoint);
00485 }
00486 }
00487 }
00488
00489
00490 ofile << BUFFER;
00491 }
00492 ofile.close();
00493
00494
00495
00496 ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
00497 bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
00498 scratchlogfile << "\n\n*******************************************************************"
00499 << "\n*_Start_" << processcontrol[pr_i_slant2h]
00500 << "\n*******************************************************************"
00501 << "\nMethod: \t\t\tschwabisch"
00502 << "\nData_output_file: \t\t"
00503 << slant2hinput.fohei
00504 << "\nData_output_format: \t\t"
00505 << "real4"
00506 << "\nEllipsoid (name,a,b): \t\t"
00507 << ellips.name << " "
00508 << ellips.a << " "
00509 << ellips.b
00510 << endl;
00511 scratchlogfile.close();
00512
00513
00514 ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
00515 bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
00516 scratchresfile << "\n\n*******************************************************************"
00517 << "\n*_Start_" << processcontrol[pr_i_slant2h]
00518 << "\n*******************************************************************"
00519 << "\nMethod: \t\t\tschwabisch"
00520 << "\nData_output_file: \t"
00521 << slant2hinput.fohei
00522 << "\nData_output_format: \t"
00523 << "real4"
00524 << "\nFirst_line (w.r.t. original_master): \t"
00525 << unwrappedinterf.win.linelo
00526 << "\nLast_line (w.r.t. original_master): \t"
00527 << unwrappedinterf.win.linehi
00528 << "\nFirst_pixel (w.r.t. original_master): \t"
00529 << unwrappedinterf.win.pixlo
00530 << "\nLast_pixel (w.r.t. original_master): \t"
00531 << unwrappedinterf.win.pixhi
00532 << "\nMultilookfactor_azimuth_direction: \t"
00533 << unwrappedinterf.multilookL
00534 << "\nMultilookfactor_range_direction: \t"
00535 << unwrappedinterf.multilookP
00536
00537 << "\nEllipsoid (name,a,b): \t"
00538 << ellips.name << " "
00539 << ellips.a << " "
00540 << ellips.b
00541 << "\n*******************************************************************"
00542
00543 << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
00544 << "\n*******************************************************************\n";
00545
00546 scratchresfile.close();
00547 PROGRESS.print("finished slant2hschwabisch.");
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567 void slant2hambiguity(
00568 const input_gen &generalinput,
00569 const input_slant2h &slant2hinput,
00570 const input_ell &ellips,
00571 const slcimage &master,
00572 const slcimage &slave,
00573 const productinfo &unwrappedinterf,
00574 orbit &masterorbit,
00575 orbit &slaveorbit,
00576 const BASELINE &baseline)
00577 {
00578 TRACE_FUNCTION("slant2hambiguity (BK 07-Jul-1999)")
00579
00580 const int32 MAXITER = 10;
00581 const real8 CRITERPOS = 1e-6;
00582 const real8 CRITERTIM = 1e-10;
00583
00584 const int32 MAXITERHERE = 10;
00585 const real8 CRITERHERE = 0.05;
00586
00587
00588
00589 const real8 multiL = unwrappedinterf.multilookL;
00590 const real8 multiP = unwrappedinterf.multilookP;
00591
00592
00593 const int32 mllines = int32(floor(real8(
00594 unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
00595 const int32 mlpixels = int32(floor(real8(
00596 unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1) / multiP));
00597
00598
00599 const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
00600 (real8(multiL) - 1.) / 2.;
00601 const real8 firstpixel = real8(unwrappedinterf.win.pixlo) +
00602 (real8(multiP) - 1.) / 2.;
00603
00604
00605
00606
00607 const int32 NUMMAT = 3;
00608 int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
00609 if (bufferlines > mllines)
00610 bufferlines = mllines;
00611
00612 const int32 FULLBUFFERS = mllines / bufferlines;
00613 const int32 RESTLINES = mllines % bufferlines;
00614 const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
00615
00616
00617
00618
00619
00620
00621
00622 window bufferwin(1, bufferlines, 1, mlpixels);
00623
00624
00625
00626 ofstream fohei;
00627 openfstream(fohei,slant2hinput.fohei,generalinput.overwrit);
00628 bk_assert(fohei,slant2hinput.fohei,__FILE__,__LINE__);
00629
00630 ofstream fophi;
00631 openfstream(fophi,slant2hinput.fophi,generalinput.overwrit);
00632 bk_assert(fophi,slant2hinput.fophi,__FILE__,__LINE__);
00633
00634 ofstream folambda;
00635 openfstream(folambda,slant2hinput.folam,generalinput.overwrit);
00636 bk_assert(folambda,slant2hinput.folam,__FILE__,__LINE__);
00637
00638
00639 cn M;
00640 cn Mdot;
00641 cn S;
00642 cn P;
00643 real8 sintheta = 0.0;
00644 real8 costheta = 0.0;
00645 real8 inc_angle = 0.0;
00646 matrix<real8> observations(3,1);
00647 matrix<real8> partials(3,3);
00648 matrix<real8> solution(3,1);
00649
00650
00651
00652 for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; ++buffer)
00653 {
00654
00655 PROGRESS << "SLANT2H: "
00656 << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
00657 PROGRESS.print();
00658
00659
00660
00661 real8 firstline = veryfirstline + real8((buffer-1) * bufferlines * multiL);
00662
00663
00664 bufferwin.linelo = 1 + (buffer-1) * bufferlines;
00665 if (buffer == FULLBUFFERS+1)
00666 {
00667 bufferlines = RESTLINES;
00668 }
00669 bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
00670
00671
00672 matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
00673 matrix<real4> PHI(BUFFER.lines(),BUFFER.pixels());
00674 matrix<real4> LAMBDA(BUFFER.lines(),BUFFER.pixels());
00675
00676
00677
00678
00679 input_ell ELLIPS;
00680 register int32 i,j;
00681 register real8 line = firstline - multiL;
00682 for (i=0; i<BUFFER.lines(); i++)
00683 {
00684 line += multiL;
00685 real8 currentheight = 0.;
00686 real8 lastheight = 0.;
00687 register real8 pixel = firstpixel - multiP;
00688
00689
00690
00691
00692 const real8 m_tazi = master.line2ta(line);
00693 M = masterorbit.getxyz(m_tazi);
00694 Mdot = masterorbit.getxyzdot(m_tazi);
00695 const real8 norm2M = M.norm2();
00696
00697
00698
00699
00700 real8 tempphase = 0.0;
00701 real8 temppixel = 0.0;
00702 real8 Bpar = 0.0;
00703 real8 Bperp = 0.0;
00704
00705 bool lineokunwrapped = false;
00706 for (int t1=0; t1<BUFFER.pixels(); ++t1)
00707 {
00708 if (BUFFER(i,t1) != NaN)
00709 lineokunwrapped = true;
00710 }
00711
00712 if (lineokunwrapped)
00713 {
00714
00715 int32 middle = BUFFER.pixels()/2;
00716 for (j=0; j<middle+1; ++j)
00717 {
00718 if (BUFFER(i,middle-j) != NaN)
00719 {
00720 tempphase = BUFFER(i,middle-j);
00721 temppixel = firstpixel + (middle-j)*multiP;
00722 break;
00723 }
00724 else if (BUFFER(i,middle+j) != NaN)
00725 {
00726 tempphase = BUFFER(i,middle+j);
00727 temppixel = firstpixel + (middle+j)*multiP;
00728 break;
00729 }
00730 }
00731
00732
00733
00734
00735 for (j=0; j<=MAXITERHERE; ++j)
00736 {
00737
00738 bool DO_NEW_METHOD=false;
00739 if (DO_NEW_METHOD==false)
00740 {
00741 real8 s_tazi;
00742 real8 s_trange;
00743 ELLIPS.a = ellips.a + currentheight;
00744 ELLIPS.b = ellips.b + currentheight;
00745 lp2xyz(line,temppixel,
00746 ELLIPS,master,masterorbit,
00747 P,MAXITER,CRITERPOS);
00748
00749
00750 xyz2t(s_tazi,s_trange,slave,slaveorbit,
00751 P,MAXITER,CRITERTIM);
00752 S = slaveorbit.getxyz(s_tazi);
00753
00754
00755
00756
00757 const real8 B = M.dist(S);
00758 Bpar = P.dist(M) - P.dist(S);
00759
00760
00761
00762
00763
00764
00765 costheta = ((-(P.norm2()) + norm2M +
00766 sqr(master.pix2range(temppixel))) /
00767 (2*sqrt(norm2M)*master.pix2range(temppixel)));
00768 sintheta = sqrt(1-sqr(costheta));
00769
00770 const cn r2 = S.min(P);
00771
00772
00773
00774
00775
00776
00777
00778 const real8 costheta2 = M.in(r2) / (M.norm()*r2.norm());
00779 Bperp = (costheta < costheta2) ?
00780 sqrt(sqr(B)-sqr(Bpar)) :
00781 -sqrt(sqr(B)-sqr(Bpar));
00782
00783 }
00784 else
00785 {
00786
00787 Bpar = baseline.get_bpar(line,temppixel,currentheight);
00788 Bperp = baseline.get_bperp(line,temppixel,currentheight);
00789 costheta = cos(baseline.get_theta_inc(line,temppixel,currentheight));
00790 sintheta = sqrt(1-sqr(costheta));
00791 inc_angle = baseline.get_theta_inc(line,temppixel,currentheight);
00792 }
00793
00794
00795 lastheight = currentheight;
00796 const real8 m_tr = master.pix2tr(temppixel);
00797 const real8 tempr1 = SOL*m_tr;
00798
00799
00800 currentheight = (-master.wavelength*tempr1*sin(inc_angle)*tempphase)/
00801 (4.*PI*Bperp);
00802
00803
00804
00805 if (abs(currentheight-lastheight) < CRITERHERE)
00806 break;
00807 }
00808
00809
00810 if (j >= MAXITERHERE)
00811 {
00812 WARNING << "slant2hambiguity: maxiter reached. "
00813 << "MAXITER: " << MAXITERHERE
00814 << "CRITER: " << CRITERHERE << "m "
00815 << "last correction: " << currentheight-lastheight;
00816 WARNING.print();
00817 }
00818 }
00819
00820
00821
00822
00823 const real8 Bhor = Bperp*costheta + Bpar*sintheta;
00824 const real8 Bver = Bperp*sintheta - Bpar*costheta;
00825 currentheight = 0.;
00826 lastheight = 0.;
00827
00828
00829
00830 for (j=0; j<BUFFER.pixels(); j++)
00831 {
00832 pixel += multiP;
00833
00834
00835 if (BUFFER(i,j) == NaN)
00836 {
00837 PHI(i,j) = NaN;
00838 LAMBDA(i,j) = NaN;
00839 }
00840
00841 else
00842 {
00843
00844
00845 const real8 m_trange = master.pix2tr(pixel);
00846 const real8 normr1 = SOL*m_trange;
00847 const real8 partnumerator = norm2M + sqr(normr1);
00848 const real8 denominator = 2.0*sqrt(norm2M)*normr1;
00849 const real8 m_lamr1phidiv4pi = master.wavelength*normr1*BUFFER(i,j) / (4.0*PI);
00850
00851
00852 register int32 iteration;
00853 for (iteration=0; iteration<=MAXITERHERE; ++iteration)
00854 {
00855 ELLIPS.a = ellips.a + currentheight;
00856 ELLIPS.b = ellips.b + currentheight;
00857
00858
00859
00860
00861
00862
00863
00864
00865 lp2xyz(line,pixel,
00866 ELLIPS,master,masterorbit,
00867 P,MAXITER,CRITERPOS);
00868
00869
00870
00871
00872 costheta = (-(P.norm2()) + partnumerator) / denominator;
00873 sintheta = sqrt(1-sqr(costheta));
00874 lastheight = currentheight;
00875
00876
00877 inc_angle = baseline.get_theta_inc(line,pixel,currentheight);
00878 currentheight = (-m_lamr1phidiv4pi * sin(inc_angle)) /
00879 (Bhor*costheta + Bver*sintheta);
00880
00881
00882
00883
00884 if (abs(currentheight-lastheight) < CRITERHERE)
00885 break;
00886 }
00887
00888
00889 if (iteration >= MAXITERHERE)
00890 {
00891 WARNING << "slant2hambiguity: maxiter reached. "
00892 << "MAXITER: " << MAXITERHERE
00893 << "CRITER: " << CRITERHERE << "m "
00894 << "last correction: " << currentheight-lastheight;
00895 WARNING.print();
00896 }
00897
00898
00899
00900 BUFFER(i,j) = currentheight;
00901
00902
00903 ELLIPS.a = ellips.a + currentheight;
00904 ELLIPS.b = ellips.b + currentheight;
00905 ELLIPS.ecc1st_sqr();
00906 ELLIPS.ecc2nd_sqr();
00907 const real8 r = sqrt(sqr(P.x)+sqr(P.y));
00908 const real8 nu = atan2((P.z*ELLIPS.a),(r*ELLIPS.b));
00909 const real8 sin3 = pow(sin(nu),3);
00910 const real8 cos3 = pow(cos(nu),3);
00911 PHI(i,j) = rad2deg(atan2((P.z+ELLIPS.e2b*ELLIPS.b*sin3),
00912 (r-ELLIPS.e2*ELLIPS.a*cos3)));
00913 LAMBDA(i,j) = rad2deg(atan2(P.y,P.x));
00914
00915 }
00916 }
00917 }
00918 fohei << BUFFER;
00919 fophi << PHI;
00920 folambda << LAMBDA;
00921 }
00922
00923
00924
00925 ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
00926 bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
00927 scratchlogfile << "\n\n*******************************************************************"
00928 << "\n*_Start_" << processcontrol[pr_i_slant2h]
00929 << "\n*******************************************************************"
00930 << "\nMethod: \t\t\tambiguity"
00931 << "\nData_output_file: \t\t"
00932 << slant2hinput.fohei
00933 << "\nData_output_format: \t\t"
00934 << "real4"
00935 << "\nData_output_file_phi: \t\t"
00936 << slant2hinput.fophi
00937 << "\nData_output_format: \t\t"
00938 << "real4"
00939 << "\nData_output_file_lam: \t\t"
00940 << slant2hinput.folam
00941 << "\nData_output_format: \t\t"
00942 << "real4"
00943 << "\nEllipsoid (name,a,b): \t\t"
00944 << ellips.name << " "
00945 << ellips.a << " "
00946 << ellips.b
00947 << endl;
00948 scratchlogfile.close();
00949
00950
00951 ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
00952 bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
00953 scratchresfile << "\n\n*******************************************************************"
00954 << "\n*_Start_" << processcontrol[pr_i_slant2h]
00955 << "\n*******************************************************************"
00956 << "\nMethod: \t"
00957 << "ambiguity"
00958 << "\nData_output_file: \t"
00959 << slant2hinput.fohei
00960 << "\nData_output_format: \t"
00961 << "real4"
00962 << "\nData_output_file_phi: \t"
00963 << slant2hinput.fophi
00964 << "\nData_output_format: \t"
00965 << "real4"
00966 << "\nData_output_file_lam: \t"
00967 << slant2hinput.folam
00968 << "\nData_output_format: \t"
00969 << "real4"
00970 << "\nFirst_line (w.r.t. original_master): \t"
00971 << unwrappedinterf.win.linelo
00972 << "\nLast_line (w.r.t. original_master): \t"
00973 << unwrappedinterf.win.linehi
00974 << "\nFirst_pixel (w.r.t. original_master): \t"
00975 << unwrappedinterf.win.pixlo
00976 << "\nLast_pixel (w.r.t. original_master): \t"
00977 << unwrappedinterf.win.pixhi
00978 << "\nMultilookfactor_azimuth_direction: \t"
00979 << unwrappedinterf.multilookL
00980 << "\nMultilookfactor_range_direction: \t"
00981 << unwrappedinterf.multilookP
00982
00983 << "\nEllipsoid (name,a,b): \t"
00984 << ellips.name << " "
00985 << ellips.a << " "
00986 << ellips.b
00987 << "\n*******************************************************************"
00988
00989 << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
00990 << "\n*******************************************************************\n";
00991
00992
00993 scratchresfile.close();
00994 PROGRESS.print("finished slant2hambiguity.");
00995 }
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 void slant2hrodriguez(
01013 const input_gen &generalinput,
01014 const input_slant2h &slant2hinput,
01015 const input_ell &ellips,
01016 const slcimage &master,
01017 const slcimage &slave,
01018 const productinfo &unwrappedinterf,
01019 const matrix<real8> &coeff_flatearth,
01020 orbit &masterorbit,
01021 orbit &slaveorbit,
01022 const BASELINE &baseline)
01023 {
01024 TRACE_FUNCTION("slant2hrodriguez (BK 30-Sep-1999)")
01025 WARNING.print("this method is based on wrong approximations. (?)");
01026
01027 const int32 MAXITER = 10;
01028 const real8 CRITERPOS = 1e-6;
01029 const real8 CRITERTIM = 1e-10;
01030
01031 const int32 degreecfe = degree(coeff_flatearth.size());
01032
01033
01034 const real8 minL = master.originalwindow.linelo;
01035 const real8 maxL = master.originalwindow.linehi;
01036 const real8 minP = master.originalwindow.pixlo;
01037 const real8 maxP = master.originalwindow.pixhi;
01038
01039
01040
01041 const real8 multiL = unwrappedinterf.multilookL;
01042 const real8 multiP = unwrappedinterf.multilookP;
01043 const real8 m_pi4divlambda = (4.*PI)/master.wavelength;
01044
01045
01046 const int32 mllines = int32(floor(real8(
01047 unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
01048 const int32 mlpixels = int32(floor(real8(
01049 unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1) / multiP));
01050
01051
01052 const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
01053 (real8(multiL) - 1.) / 2.;
01054 const real8 firstpixel = real8(unwrappedinterf.win.pixlo) +
01055 (real8(multiP) - 1.) / 2.;
01056
01057
01058
01059 const int32 NUMMAT = 1;
01060 int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
01061 if (bufferlines > mllines)
01062 bufferlines = mllines;
01063
01064 const int32 FULLBUFFERS = mllines / bufferlines;
01065 const int32 RESTLINES = mllines % bufferlines;
01066 const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
01067
01068
01069
01070 const uint dummy = 999999;
01071 window bufferwin(1, bufferlines, 1, mlpixels);
01072 window offsetbuffer(1,dummy,1,dummy);
01073
01074
01075
01076 ofstream fohei;
01077 openfstream(fohei,slant2hinput.fohei,generalinput.overwrit);
01078 bk_assert(fohei,slant2hinput.fohei,__FILE__,__LINE__);
01079
01080
01081
01082
01083
01084 for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; ++buffer)
01085 {
01086
01087 PROGRESS << "SLANT2H: "
01088 << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
01089 PROGRESS.print();
01090
01091
01092
01093 real8 firstline = veryfirstline + real8((buffer-1) * bufferlines * multiL);
01094
01095
01096 bufferwin.linelo = 1 + (buffer-1) * bufferlines;
01097 if (buffer == FULLBUFFERS+1)
01098 {
01099 bufferlines = RESTLINES;
01100
01101 }
01102 bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
01103
01104
01105 matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
01106
01107
01108
01109 register int32 i,j;
01110 register real8 line = firstline - multiL;
01111 for (i=0; i<BUFFER.lines(); i++)
01112 {
01113 line += multiL;
01114 register real8 pixel = firstpixel - multiP;
01115
01116
01117
01118
01119 const real8 m_tazi = master.line2ta(line);
01120 cn M = masterorbit.getxyz(m_tazi);
01121 cn P;
01122
01123
01124
01125
01126 const real8 middlepixel = pixel+(.5*BUFFER.pixels()*multiP);
01127 lp2xyz(line,middlepixel,
01128 ellips,master,masterorbit,
01129 P,MAXITER,CRITERPOS);
01130
01131
01132
01133
01134
01135 real8 s_tazi;
01136 real8 s_trange;
01137 xyz2t(s_tazi,s_trange,slave, slaveorbit,
01138 P,MAXITER,CRITERTIM);
01139 cn S = slaveorbit.getxyz(s_tazi);
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153 const real8 B = M.dist(S);
01154 const real8 Bpartemp = P.dist(M) - P.dist(S);
01155
01156
01157 const real8 rho1sqr = M.norm2();
01158 const real8 costhetatemp = ((-(P.norm2()) + rho1sqr +
01159 sqr(master.pix2range(middlepixel))) /
01160 (2*sqrt(rho1sqr)*master.pix2range(middlepixel)));
01161
01162
01163 const cn r2 = S.min(P);
01164 const real8 Bperp = (acos(costhetatemp) > M.angle(r2)) ?
01165 sqrt(sqr(B)-sqr(Bpartemp)) :
01166 -sqrt(sqr(B)-sqr(Bpartemp));
01167
01168 const real8 alpha = acos(costhetatemp) - atan2(Bpartemp,Bperp);
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179 const real8 rho1 = sqrt(rho1sqr);
01180 real8 satphi,satlambda,satheight;
01181 xyz2ell(ellips,M,satphi,satlambda,satheight);
01182 const real8 Radius = rho1 - satheight;
01183
01184
01185
01186 for (j=0; j<BUFFER.pixels(); j++)
01187 {
01188 pixel += multiP;
01189
01190
01191 if (BUFFER(i,j) != NaN)
01192 {
01193
01194
01195 const real8 r1 = master.pix2range(pixel);
01196 const real8 r1sqr = sqr(r1);
01197
01198
01199 const real8 phiref= polyval(normalize(line,minL,maxL),
01200 normalize(pixel,minP,maxP),
01201 coeff_flatearth,degreecfe);
01202
01203
01204 const real8 Bpar = -(BUFFER(i,j)+phiref)/m_pi4divlambda;
01205
01206 const real8 sinthetaminalpha =
01207 (sqr(r1)+sqr(B)-sqr(r1-Bpar))/(2*r1*B);
01208
01209
01210
01211
01212 real8 theta = asin(sinthetaminalpha)+alpha;
01213 if (theta<0.0 || theta>1.0)
01214 theta = PI-asin(sinthetaminalpha)+alpha;
01215
01216 const real8 costheta = cos(theta);
01217 const real8 psqr = rho1sqr+r1sqr - 2*rho1*r1*costheta;
01218 const real8 p = sqrt(psqr);
01219 const real8 cosmu = ((-r1sqr+psqr+rho1sqr) / (2*p*rho1));
01220 const real8 H = rho1 - Radius*cosmu;
01221
01222
01223
01224 BUFFER(i,j) = H - r1*costheta;
01225 }
01226 }
01227 }
01228 fohei << BUFFER;
01229 }
01230
01231
01232
01233 ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
01234 bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
01235 scratchlogfile << "\n\n*******************************************************************"
01236 << "\n*_Start_" << processcontrol[pr_i_slant2h]
01237 << "\n*******************************************************************"
01238 << "\nMethod: \t\t\trodriguez"
01239 << "\nData_output_file: \t\t"
01240 << slant2hinput.fohei
01241 << "\nData_output_format: \t\t"
01242 << "real4"
01243 << "\nEllipsoid (name,a,b): \t\t"
01244 << ellips.name << " "
01245 << ellips.a << " "
01246 << ellips.b
01247 << endl;
01248 scratchlogfile.close();
01249
01250
01251 ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
01252 bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
01253 scratchresfile << "\n\n*******************************************************************"
01254 << "\n*_Start_" << processcontrol[pr_i_slant2h]
01255 << "\n*******************************************************************"
01256 << "\nMethod: \t"
01257 << "rodriguez"
01258 << "\nData_output_file: \t"
01259 << slant2hinput.fohei
01260 << "\nData_output_format: \t"
01261 << "real4"
01262 << "\nFirst_line (w.r.t. original_master): \t"
01263 << unwrappedinterf.win.linelo
01264 << "\nLast_line (w.r.t. original_master): \t"
01265 << unwrappedinterf.win.linehi
01266 << "\nFirst_pixel (w.r.t. original_master): \t"
01267 << unwrappedinterf.win.pixlo
01268 << "\nLast_pixel (w.r.t. original_master): \t"
01269 << unwrappedinterf.win.pixhi
01270 << "\nMultilookfactor_azimuth_direction: \t"
01271 << unwrappedinterf.multilookL
01272 << "\nMultilookfactor_range_direction: \t"
01273 << unwrappedinterf.multilookP
01274
01275 << "\nEllipsoid (name,a,b): \t"
01276 << ellips.name << " "
01277 << ellips.a << " "
01278 << ellips.b
01279 << "\n*******************************************************************"
01280
01281 << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
01282 << "\n*******************************************************************\n";
01283
01284
01285
01286 scratchresfile.close();
01287 PROGRESS.print("finished slant2hrodriguez.");
01288 }
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309 void geocode(
01310 const input_gen &generalinput,
01311 const input_geocode &geocodeinput,
01312 const input_ell &ellips,
01313 const slcimage &master,
01314 const productinfo &heightinradarsystem,
01315 orbit &masterorbit)
01316 {
01317 TRACE_FUNCTION("geocode (BK 02-Jun-1999)")
01318 const int32 MAXITER = 10;
01319 const real8 CRITERPOS = 1e-6;
01320 const real8 CRITERTIM = 1e-10;
01321
01322
01323
01324
01325
01326 const real8 multiL = heightinradarsystem.multilookL;
01327 const real8 multiP = heightinradarsystem.multilookP;
01328
01329
01330 const int32 mllines = int32(floor(real8(
01331 heightinradarsystem.win.linehi-heightinradarsystem.win.linelo+1) / multiL));
01332 const int32 mlpixels = int32(floor(real8(
01333 heightinradarsystem.win.pixhi-heightinradarsystem.win.pixlo+1) / multiP));
01334
01335
01336 const real8 veryfirstline = real8(heightinradarsystem.win.linelo) +
01337 (multiL - 1.) / 2.;
01338 const real8 firstpixel = real8(heightinradarsystem.win.pixlo) +
01339 (multiP - 1.) / 2.;
01340
01341
01342
01343 const int32 NUMMAT = 3;
01344 int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
01345 if (bufferlines > mllines)
01346 bufferlines = mllines;
01347
01348 const int32 FULLBUFFERS = mllines / bufferlines;
01349 const int32 RESTLINES = mllines % bufferlines;
01350 const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
01351
01352
01353
01354 matrix<real4> PHI(bufferlines,mlpixels);
01355 matrix<real4> LAMBDA(bufferlines,mlpixels);
01356 matrix<real4> HEIGHT(bufferlines,mlpixels);
01357 const uint dummy = 999999;
01358 window bufferwin(1, bufferlines, 1, mlpixels);
01359 window offsetbuffer(1,dummy,1,dummy);
01360
01361
01362
01363 ofstream fophi;
01364 openfstream(fophi,geocodeinput.fophi,generalinput.overwrit);
01365 bk_assert(fophi,geocodeinput.fophi,__FILE__,__LINE__);
01366
01367 ofstream folam;
01368 openfstream(folam,geocodeinput.folam,generalinput.overwrit);
01369 bk_assert(folam,geocodeinput.folam,__FILE__,__LINE__);
01370
01371
01372
01373 for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; buffer++)
01374 {
01375
01376 PROGRESS << "GEOCODE: "
01377 << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
01378 PROGRESS.print();
01379
01380
01381 real8 firstline = veryfirstline + (buffer-1) * bufferlines * multiL;
01382
01383
01384 bufferwin.linelo = 1 + (buffer-1) * bufferlines;
01385 if (buffer == FULLBUFFERS+1)
01386 {
01387 bufferlines = RESTLINES;
01388 PHI.resize(bufferlines,mlpixels);
01389 LAMBDA.resize(bufferlines,mlpixels);
01390 HEIGHT.resize(bufferlines,mlpixels);
01391 }
01392 bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
01393
01394
01395
01396 switch (heightinradarsystem.formatflag)
01397 {
01398 case FORMATR4:
01399 readfile (HEIGHT, heightinradarsystem.file, mllines, bufferwin,
01400 offsetbuffer);
01401 break;
01402
01403 default:
01404 PRINT_ERROR("geocode format flag on file heights (s2h output) only real4 possible.");
01405 throw(unhandled_case_error);
01406 }
01407
01408
01409 register int32 i,j;
01410 register real8 line = firstline - multiL;
01411 cn pospoint;
01412 input_ell ELLIPS;
01413 real8 r;
01414 real8 nu;
01415 real8 sin3;
01416 real8 cos3;
01417
01418
01419 for (i=0; i<HEIGHT.lines(); i++)
01420 {
01421 line += multiL;
01422 register real8 pixel = firstpixel - multiP;
01423
01424 for (j=0; j<HEIGHT.pixels(); j++)
01425 {
01426 pixel += multiP;
01427
01428
01429 if (HEIGHT(i,j) == NaN)
01430 {
01431 PHI(i,j) = NaN;
01432 LAMBDA(i,j) = NaN;
01433 }
01434
01435 else
01436 {
01437
01438 ELLIPS.a = ellips.a + HEIGHT(i,j);
01439 ELLIPS.b = ellips.b + HEIGHT(i,j);
01440 ELLIPS.ecc1st_sqr();
01441 ELLIPS.ecc2nd_sqr();
01442
01443
01444
01445
01446 lp2xyz(line,pixel,ELLIPS,master,masterorbit,
01447 pospoint,MAXITER,CRITERPOS);
01448
01449
01450 r = sqrt(sqr(pospoint.x)+sqr(pospoint.y));
01451 nu = atan2((pospoint.z*ELLIPS.a),(r*ELLIPS.b));
01452 sin3 = pow(sin(nu),3);
01453 cos3 = pow(cos(nu),3);
01454 PHI(i,j) = rad2deg(atan2((pospoint.z+ELLIPS.e2b*ELLIPS.b*sin3),
01455 (r-ELLIPS.e2*ELLIPS.a*cos3)));
01456 LAMBDA(i,j) = rad2deg(atan2(pospoint.y,pospoint.x));
01457 }
01458 }
01459 }
01460 fophi << PHI;
01461 folam << LAMBDA;
01462 }
01463
01464
01465 fophi.close();
01466 folam.close();
01467
01468
01469 ofstream scratchlogfile("scratchloggeocode", ios::out | ios::trunc);
01470 bk_assert(scratchlogfile,"geocode: scratchloggeocode",__FILE__,__LINE__);
01471 scratchlogfile << "\n\n*******************************************************************"
01472 << "\n*_Start_" << processcontrol[pr_i_geocoding]
01473 << "\n*******************************************************************"
01474 << "\nData_output_file_hei (slant2h): "
01475 << heightinradarsystem.file
01476 << "\nData_output_file_phi: \t\t"
01477 << geocodeinput.fophi
01478 << "\nData_output_file_lambda: \t"
01479 << geocodeinput.folam
01480 << endl;
01481 scratchlogfile.close();
01482
01483
01484 ofstream scratchresfile("scratchresgeocode", ios::out | ios::trunc);
01485 bk_assert(scratchresfile,"geocode: scratchresgeocode",__FILE__,__LINE__);
01486 scratchresfile << "\n\n*******************************************************************"
01487 << "\n*_Start_" << processcontrol[pr_i_geocoding]
01488 << "\n*******************************************************************"
01489 << "\nData_output_file_hei (slant2h): "
01490 << heightinradarsystem.file
01491 << "\nData_output_file_phi: \t"
01492 << geocodeinput.fophi
01493 << "\nData_output_file_lamda: \t"
01494 << geocodeinput.folam
01495 << "\nData_output_format: \t"
01496 << "real4"
01497 << "\nFirst_line (w.r.t. original_master): \t"
01498 << heightinradarsystem.win.linelo
01499 << "\nLast_line (w.r.t. original_master): \t"
01500 << heightinradarsystem.win.linehi
01501 << "\nFirst_pixel (w.r.t. original_master): \t"
01502 << heightinradarsystem.win.pixlo
01503 << "\nLast_pixel (w.r.t. original_master): \t"
01504 << heightinradarsystem.win.pixhi
01505 << "\nMultilookfactor_azimuth_direction: \t"
01506 << heightinradarsystem.multilookL
01507 << "\nMultilookfactor_range_direction: \t"
01508 << heightinradarsystem.multilookP
01509 << "\n*******************************************************************"
01510 << "\n* End_" << processcontrol[pr_i_geocoding] << "_NORMAL"
01511 << "\n*******************************************************************\n";
01512 scratchresfile.close();
01513 PROGRESS.print("finished geocode.");
01514 }
01515