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 "orbitbk.hh"
00047 #include "slcimage.hh"
00048 #include "productinfo.hh"
00049 #include "products.hh"
00050 #include "utilities.hh"
00051 #include "ioroutines.hh"
00052 #include "conversion.hh"
00053 #include "exceptions.hh"
00054 #include "bk_baseline.hh"
00055
00056 #include <algorithm>
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079 void compinterfero(
00080 const slcimage &master,
00081 const slcimage &slave,
00082 const input_gen &input_general,
00083 const input_interfero &input_i_interfero,
00084 const matrix<real8> &coeff_flatearth)
00085 {
00086 TRACE_FUNCTION("compinterfero (BK 06-Oct-1999)");
00087 INFO << "INTERFERO: master input file: " << master.file;
00088 INFO.print();
00089 INFO << "INTERFERO: slave input file: " << slave.file;
00090 INFO.print();
00091
00092 #ifdef __DEBUG
00093
00094 if (slave.currentwindow.linelo < master.currentwindow.linelo ||
00095 slave.currentwindow.linehi > master.currentwindow.linehi ||
00096 slave.currentwindow.pixlo < master.currentwindow.pixlo ||
00097 slave.currentwindow.pixhi > master.currentwindow.pixhi )
00098 {
00099 slave.currentwindow.disp();
00100 master.currentwindow.disp();
00101 slave.currentwindow.disp();
00102 PRINT_ERROR("Panic, interferowin smaller than master win.")
00103 throw(file_error);
00104 }
00105 #endif
00106
00107 INFO << "compinterfero: inteferogram for (master coord.): "
00108 << slave.currentwindow.linelo << ":" << slave.currentwindow.linehi << "; "
00109 << slave.currentwindow.pixlo << ":" << slave.currentwindow.pixhi;
00110 INFO.print();
00111
00112
00113 bool nocint = true;
00114 bool noint = true;
00115 bool noflatearthcorrection = false;
00116 if (specified(input_i_interfero.focint))
00117 nocint = false;
00118 if (specified(input_i_interfero.foint))
00119 noint = false;
00120 if (coeff_flatearth.size() == 0)
00121 noflatearthcorrection = true;
00122
00123 const uint BUFFERMEMSIZE = input_general.memory;
00124 const int32 multiL = input_i_interfero.multilookL;
00125 const int32 multiP = input_i_interfero.multilookP;
00126
00127
00128
00129 const real8 minL = master.originalwindow.linelo;
00130 const real8 maxL = master.originalwindow.linehi;
00131 const real8 minP = master.originalwindow.pixlo;
00132 const real8 maxP = master.originalwindow.pixhi;
00133 INFO << "compinterfero: polynomial normalized by factors: "
00134 << minL << " " << maxL << " " << minP << " " << maxP
00135 << " to [-2,2]";
00136 INFO.print();
00137
00138
00139 ofstream ofilecint;
00140 if (!nocint)
00141 {
00142 openfstream(ofilecint,input_i_interfero.focint,input_general.overwrit);
00143 bk_assert(ofilecint,input_i_interfero.focint,__FILE__,__LINE__);
00144 }
00145 ofstream ofileint;
00146 if (!noint)
00147 {
00148 openfstream(ofileint,input_i_interfero.foint,input_general.overwrit);
00149 bk_assert(ofileint,input_i_interfero.foint,__FILE__,__LINE__);
00150 }
00151
00152
00153
00154 const int32 numpixels = (slave.currentwindow.pixels());
00155 const int32 bytesperline = numpixels * sizeof(complr4);
00156
00157 const real4 numbigmatrices = (noflatearthcorrection) ? 3.2 : 4.2;
00158 int32 numlines = int32((BUFFERMEMSIZE / numbigmatrices) / bytesperline);
00159
00160 while (numlines%multiL)
00161 numlines -= 1;
00162 if (numlines%multiL)
00163 {
00164 PRINT_ERROR("panic : totally impossible on HP aCC compiler.")
00165 throw(input_error);
00166 }
00167 int32 nummllines = numlines/multiL;
00168 int32 nummlpixels = numpixels/multiP;
00169 if (nummllines < 1)
00170 {
00171 PRINT_ERROR("Please increase memory (MEMORY card) or decrease multiL (INT_MULTILOOK card).")
00172 throw(input_error);
00173 }
00174
00175
00176
00177
00178 const uint Mfilelines = master.currentwindow.lines();
00179 const uint Ifilelines = slave.currentwindow.lines();
00180
00181
00182 window winfile(slave.currentwindow.linelo,
00183 slave.currentwindow.linelo + numlines - 1,
00184 slave.currentwindow.pixlo,
00185 slave.currentwindow.pixhi);
00186
00187
00188
00189
00190
00191 register int32 i,blocks;
00192 const int32 numfullbuffers = Ifilelines/numlines;
00193 const int32 numrestlines = Ifilelines%numlines;
00194 const int32 nummlrestlines = numrestlines/multiL;
00195 const int32 EXTRABUFFER = nummlrestlines ? 1 : 0;
00196
00197 matrix<real4> p_axis(numpixels,1);
00198 for (i=0; i<numpixels; i++)
00199 p_axis(i,0) = winfile.pixlo + i;
00200 normalize(p_axis,minP,maxP);
00201
00202 for (blocks=1; blocks<=numfullbuffers+EXTRABUFFER; blocks++)
00203 {
00204
00205 PROGRESS << "INTERFERO: "
00206 << int32(100*real4(blocks-1)/(real4(Ifilelines)/real4(numlines))+.5)
00207 << "%";
00208 PROGRESS.print();
00209
00210
00211 if (blocks == (numfullbuffers+1))
00212 {
00213 numlines = multiL * nummlrestlines;
00214 winfile.linehi = numlines + winfile.linelo - 1;
00215 }
00216
00217
00218 matrix<complr4> MASTER = master.readdata(winfile);
00219 matrix<complr4> SLAVE = slave.readdata(winfile);
00220
00221
00222
00223
00224 if (!noflatearthcorrection)
00225 {
00226 matrix<real4> l_axis(numlines,1);
00227 for (i=0; i<numlines; i++)
00228 l_axis(i,0) = winfile.linelo + i;
00229
00230 normalize(l_axis,minL,maxL);
00231
00232 matrix<real4> REFPHASE = polyval(l_axis, p_axis, coeff_flatearth);
00233 SLAVE *= angle2cmplx(REFPHASE);
00234 }
00235
00236
00237 MASTER *= conj(SLAVE);
00238
00239
00240 matrix<complr4> ML = multilook(MASTER,multiL,multiP);
00241
00242
00243 if (!nocint)
00244 ofilecint << ML;
00245 if (!noint)
00246 ofileint << angle(ML);
00247
00248
00249
00250 winfile.linelo = winfile.linehi + 1;
00251 winfile.linehi += numlines;
00252 }
00253
00254
00255
00256
00257 ofstream scratchlogfile("scratchloginterfero", ios::out | ios::trunc);
00258 bk_assert(scratchlogfile,"compinterfero: scratchloginterfero",__FILE__,__LINE__);
00259
00260 scratchlogfile
00261 << "\n\n*******************************************************************"
00262 << "\n Computation of interferogram"
00263 << "\nInput file master (format): \t\t\t"
00264 << master.file << " " << master.formatflag
00265 << "\nInput file slave (format): \t\t\t"
00266 << slave.file << " " << slave.formatflag
00267 << "\ncomplex interferogram: \t\t"
00268 << input_i_interfero.focint
00269 << "\ninterferogram: \t\t"
00270 << input_i_interfero.foint
00271 << "\nNumber of lines (multilooked): \t"
00272 << Ifilelines / multiL
00273 << "\nNumber of pixels (multilooked): \t"
00274 << nummlpixels
00275 << "\nMultilookfactor in line direction: \t"
00276 << multiL
00277 << "\nMultilookfactor in pixel direction: \t"
00278 << multiP
00279 << "\nTip: dismph " << input_i_interfero.focint
00280 << " " << nummlpixels << " 1 500"
00281 << "\n*******************************************************************\n";
00282 scratchlogfile.close();
00283
00284 ofstream scratchresfile("scratchresinterfero", ios::out | ios::trunc);
00285 bk_assert(scratchresfile,"compinterfero: scratchresinterfero",__FILE__,__LINE__);
00286
00287 scratchresfile << "\n\n*******************************************************************"
00288 << "\n*_Start_" << processcontrol[pr_i_interfero]
00289 << "\n*******************************************************************";
00290 if (!nocint)
00291 {
00292 scratchresfile
00293 << "\nData_output_file: \t\t\t"
00294 << input_i_interfero.focint
00295 << "\nData_output_format: \t\t\t"
00296 << "complex_real4";
00297 }
00298 if (!noint)
00299 {
00300 scratchresfile
00301 << "\nData_output_file_real_interferogram: \t"
00302 << input_i_interfero.foint
00303 << "\nData_output_format_real_interferogram: \t\t"
00304 << "real4";
00305 }
00306 scratchresfile << "\nFlatearth correction subtracted: \t";
00307 if (!noflatearthcorrection)
00308 scratchresfile << "yes";
00309 else
00310 scratchresfile << "no";
00311 scratchresfile
00312 << "\nFirst_line (w.r.t. original_master): \t"
00313 << slave.currentwindow.linelo
00314 << "\nLast_line (w.r.t. original_master): \t"
00315 << slave.currentwindow.linehi
00316 << "\nFirst_pixel (w.r.t. original_master): \t"
00317 << slave.currentwindow.pixlo
00318 << "\nLast_pixel (w.r.t. original_master): \t"
00319 << slave.currentwindow.pixhi
00320 << "\nMultilookfactor_azimuth_direction: \t"
00321 << multiL
00322 << "\nMultilookfactor_range_direction: \t"
00323 << multiP
00324 << "\nNumber of lines (multilooked): \t\t"
00325 << Ifilelines / multiL
00326 << "\nNumber of pixels (multilooked): \t"
00327 << nummlpixels
00328 << "\n*******************************************************************"
00329 << "\n* End_" << processcontrol[pr_i_interfero] << "_NORMAL"
00330 << "\n*******************************************************************\n";
00331 scratchresfile.close();
00332
00333
00334
00335 PROGRESS.print("finished compinterfero.");
00336 if (!nocint)
00337 ofilecint.close();
00338 if (!noint)
00339 ofileint.close();
00340 }
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 void compcoherence(
00360 const slcimage &master,
00361 const slcimage &slave,
00362 const input_gen &input_general,
00363 const input_coherence &input_i_coherence,
00364 const matrix<real8> &coeff_flatearth)
00365 {
00366 TRACE_FUNCTION("compcoherence (BK 16-Apr-1999)")
00367 #ifdef __DEBUG
00368
00369 if (slave.currentwindow.linelo < master.currentwindow.linelo ||
00370 slave.currentwindow.linehi > master.currentwindow.linehi ||
00371 slave.currentwindow.pixlo < master.currentwindow.pixlo ||
00372 slave.currentwindow.pixhi > master.currentwindow.pixhi )
00373 {
00374 PRINT_ERROR("Panic, impossible 3333.")
00375 throw(input_error);
00376 }
00377 #endif
00378
00379 bool nocoh = true;
00380 bool noccoh = true;
00381 bool noflatearthcorrection = false;
00382 if (specified(input_i_coherence.focoh))
00383 nocoh = false;
00384 if (specified(input_i_coherence.foccoh))
00385 noccoh = false;
00386 if (coeff_flatearth.size() == 0)
00387 {
00388 WARNING.print("Computation of coherence without subtracting of reference image.");
00389 noflatearthcorrection = true;
00390 }
00391
00392
00393 const real8 minL = master.originalwindow.linelo;
00394 const real8 maxL = master.originalwindow.linehi;
00395 const real8 minP = master.originalwindow.pixlo;
00396 const real8 maxP = master.originalwindow.pixhi;
00397 INFO << "compcoherence: polynomial normalized by factors: "
00398 << minL << " " << maxL << " " << minP << " " << maxP
00399 << " to [-2,2]";
00400 INFO.print();
00401
00402
00403 const uint BUFFERMEMSIZE = input_general.memory;
00404 const int32 multiL = input_i_coherence.multilookL;
00405 const int32 multiP = input_i_coherence.multilookP;
00406
00407 const int32 winsizeL = input_i_coherence.cohsizeL;
00408 const int32 winsizeP = input_i_coherence.cohsizeP;
00409
00410
00411 ofstream ofileccoh;
00412 if (!noccoh)
00413 {
00414 openfstream(ofileccoh,input_i_coherence.foccoh,input_general.overwrit);
00415 bk_assert(ofileccoh,input_i_coherence.foccoh,__FILE__,__LINE__);
00416 }
00417 ofstream ofilecoh;
00418 if (!nocoh)
00419 {
00420 openfstream(ofilecoh,input_i_coherence.focoh,input_general.overwrit);
00421 bk_assert(ofilecoh,input_i_coherence.focoh,__FILE__,__LINE__);
00422 }
00423
00424
00425
00426 const int32 numpixels = (slave.currentwindow.pixhi-slave.currentwindow.pixlo+1);
00427 const int32 bytesperline = numpixels * sizeof(complr4);
00428 int32 numlines;
00429
00430
00431 if (noflatearthcorrection)
00432 numlines = int32((BUFFERMEMSIZE / 3.0) / bytesperline);
00433 else
00434 numlines = int32((BUFFERMEMSIZE / 3.5) / bytesperline);
00435
00436 while (numlines%multiL)
00437 numlines -= 1;
00438 if (numlines%multiL)
00439 {
00440 PRINT_ERROR("PANIC: multilookwindow not exactly in numlines.")
00441 throw(input_error);
00442 }
00443
00444 int32 nummllines = numlines/multiL;
00445 int32 nummlpixels = numpixels/multiP;
00446 if (nummllines < 1)
00447 {
00448 PRINT_ERROR("increase memory (MEMORY card) or reduce multiL (COH_MULTILOOK card).")
00449 throw(input_error);
00450 }
00451
00452
00453 int32 numlines2 = numlines;
00454 numlines += (winsizeL - 1);
00455
00456
00457 const uint Mfilelines = master.currentwindow.lines();
00458 const uint Ifilelines = slave.currentwindow.lines();
00459
00460
00461
00462
00463
00464
00465 register int32 i,j;
00466 uint zerolinesstart = (winsizeL-1)/2;
00467 uint trailinglines = (winsizeL)/2;
00468 int32 extrazerolines = 0;
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499 const int32 BUF= 1+winsizeL/multiL;
00500 window winfile(slave.currentwindow.linelo,
00501 slave.currentwindow.linelo + BUF*multiL + trailinglines - 1,
00502 slave.currentwindow.pixlo,
00503 slave.currentwindow.pixhi);
00504
00505
00506 matrix<complr4> SLAVE(winfile.linehi-winfile.linelo+1,numpixels);
00507 matrix<complr4> MASTER;
00508
00509
00510 matrix<real4> p_axis;
00511 if (!noflatearthcorrection)
00512 {
00513 p_axis.resize(numpixels,1);
00514 for (i=0; i<numpixels; i++)
00515 p_axis(i,0) = winfile.pixlo + i;
00516
00517 normalize(p_axis,minP,maxP);
00518 }
00519
00520
00521
00522
00523 bool firstblock = true;
00524 bool lastblock = false;
00525 for (register int32 blocks=0; blocks<100000000; blocks++)
00526 {
00527 INFO << "coherence block: " << blocks;
00528 INFO.print();
00529
00530 if (lastblock)
00531 break;
00532 if (winfile.linehi > slave.currentwindow.linehi)
00533 {
00534 lastblock = true;
00535
00536 int32 restlines = Ifilelines%multiL;
00537 if (winfile.linelo+zerolinesstart >= slave.currentwindow.linehi-restlines+1)
00538 break;
00539
00540
00541
00542 int32 lastlinerequired = slave.currentwindow.linehi - restlines + trailinglines;
00543 if (lastlinerequired > slave.currentwindow.linehi)
00544 {
00545 winfile.linehi = slave.currentwindow.linehi;
00546 extrazerolines = lastlinerequired - slave.currentwindow.linehi;
00547 }
00548 else
00549 {
00550 winfile.linehi = lastlinerequired;
00551 }
00552 INFO << "coherence: extra zero lines last buffer: " << extrazerolines;
00553 INFO.print();
00554 MASTER.resize(winfile.linehi-winfile.linelo+1,numpixels);
00555 SLAVE.resize (winfile.linehi-winfile.linelo+1,numpixels);
00556 }
00557
00558 else
00559 {
00560
00561 PROGRESS << "COHERENCE: "
00562 << int32(100*real4(blocks)/(real4(Ifilelines)/real4(numlines2))+.5)
00563 << "%";
00564 PROGRESS.print();
00565 }
00566
00567
00568 INFO << "coherence winfile: [" << winfile.linelo << ":" << winfile.linehi
00569 << ", " << winfile.pixlo << ":" << winfile.pixhi << "]";
00570 INFO.print();
00571 MASTER = master.readdata(winfile);
00572 SLAVE = slave.readdata(winfile);
00573
00574
00575 if (firstblock == true)
00576 {
00577 matrix<complr4> TMP = SLAVE;
00578 SLAVE.resize(BUF*multiL+winsizeL-1,TMP.pixels());
00579 const window wintmp(SLAVE.lines()-TMP.lines(), SLAVE.lines()-1,
00580 0, SLAVE.pixels()-1);
00581 const window windef(0, 0, 0, 0);
00582 SLAVE.setdata(wintmp,TMP,windef);
00583 TMP = MASTER;
00584 MASTER.resize(BUF*multiL+winsizeL-1,TMP.pixels());
00585 MASTER.setdata(wintmp,TMP,windef);
00586 INFO << "coherence: increase lines for first block: "
00587 << SLAVE.lines()-TMP.lines();
00588 INFO.print();
00589 }
00590
00591
00592
00593 if (extrazerolines != 0)
00594 {
00595 const window wintmp(0, MASTER.lines()-1, 0, MASTER.pixels()-1);
00596 const window windef(0, 0, 0, 0);
00597 matrix<complr4> TMP = SLAVE;
00598 SLAVE.resize(TMP.lines()+extrazerolines,TMP.pixels());
00599 SLAVE.setdata(wintmp,TMP,windef);
00600 TMP = MASTER;
00601 MASTER.resize(TMP.lines()+extrazerolines,TMP.pixels());
00602 MASTER.setdata(wintmp,TMP,windef);
00603 INFO << "coherence: increase lines for last block: "
00604 << SLAVE.lines()-TMP.lines();
00605 INFO.print();
00606 }
00607
00608
00609 uint sizeL = SLAVE.lines();
00610 uint sizeP = SLAVE.pixels();
00611 if (!noflatearthcorrection)
00612 {
00613 matrix<real4> l_axis(sizeL,1);
00614 for (i=0; i<sizeL; i++)
00615 l_axis(i,0) = winfile.linelo + i;
00616
00617 normalize(l_axis,minL,maxL);
00618 matrix<real4> REFPHASE = polyval(l_axis, p_axis, coeff_flatearth);
00619
00620 for (i=0; i<sizeL; i++)
00621 {
00622 for (j=0; j<sizeP; j++)
00623 {
00624 real4 tmp = norm(MASTER(i,j));
00625 MASTER(i,j) *= (conj(SLAVE(i,j)) *
00626 complr4(cos(REFPHASE(i,j)),-sin(REFPHASE(i,j))));
00627 SLAVE(i,j) = complr4(norm(SLAVE(i,j)),tmp);
00628 }
00629 }
00630 }
00631
00632
00633 else
00634 {
00635 for (i=0; i<sizeL; i++)
00636 {
00637 for (j=0; j<sizeP; j++)
00638 {
00639 real4 tmp = norm(MASTER(i,j));
00640 MASTER(i,j) *= conj(SLAVE(i,j));
00641 SLAVE(i,j) = complr4(norm(SLAVE(i,j)),tmp);
00642 }
00643 }
00644 }
00645
00646
00647
00648 INFO << "block: " << blocks << "; num output lines: "
00649 << real8(SLAVE.lines()-winsizeL+1.0)/real8(multiL);
00650 INFO.print();
00651 INFO << "SLAVE number of lines this buffer block: " << SLAVE.lines();
00652 INFO.print();
00653 if (!noccoh)
00654 {
00655 if (!nocoh)
00656 {
00657 matrix<complr4> CCOHERENCE = multilook(
00658 coherence(MASTER,SLAVE,winsizeL,winsizeP),multiL,multiP);
00659 ofileccoh << CCOHERENCE;
00660 ofilecoh << magnitude(CCOHERENCE);
00661 }
00662 else
00663 {
00664 ofileccoh << multilook(
00665 coherence(MASTER,SLAVE,winsizeL,winsizeP),multiL,multiP);
00666 }
00667 }
00668
00669 else if (!nocoh)
00670 {
00671
00672
00673
00674 matrix<real4> COHERENCE =
00675 coherence2(MASTER,SLAVE,winsizeL,winsizeP);
00676 INFO << "block: " << blocks << "; size COHERENCE matrix: "
00677 << COHERENCE.lines();
00678 INFO.print();
00679 ofilecoh << multilook(COHERENCE, multiL,multiP);
00680 }
00681
00682
00683 else
00684 {
00685 PRINT_ERROR("2212 panic impossible.")
00686 throw(unhandled_case_error);
00687 }
00688
00689
00690
00691 winfile.linelo = winfile.linehi - trailinglines + 1 - zerolinesstart;
00692 winfile.linehi = winfile.linelo + numlines - 1;
00693
00694 if (firstblock)
00695 {
00696 firstblock = false;
00697 MASTER.resize(winfile.linehi-winfile.linelo+1,numpixels);
00698 SLAVE.resize(winfile.linehi-winfile.linelo+1,numpixels);
00699 }
00700 }
00701
00702
00703
00704
00705 ofstream scratchlogfile("scratchlogcoherence", ios::out | ios::trunc);
00706 bk_assert(scratchlogfile,"compcoherence: scratchlogcoherence",__FILE__,__LINE__);
00707 scratchlogfile
00708 << "\n\n*******************************************************************"
00709 << "\n Computation of coherence"
00710 << "\nInput file master (format): \t"
00711 << master.file << " " << master.formatflag
00712 << "\nInput file slave (format): \t"
00713 << slave.file << " " << slave.formatflag
00714 << "\ncomplex coherence: \t"
00715 << input_i_coherence.foccoh
00716 << "\ncoherence: \t"
00717 << input_i_coherence.focoh
00718 << "\nNumber of lines (multilooked): \t"
00719 << Ifilelines / multiL
00720 << "\nNumber of pixels (multilooked): \t"
00721 << nummlpixels
00722 << "\nMultilookfactor in line direction: \t"
00723 << multiL
00724 << "\nMultilookfactor in pixel direction: \t"
00725 << multiP
00726 << "\nNumber of lines window for coherence estimation: "
00727 << winsizeL
00728 << "\nNumber of pixels window for coherence estimation: "
00729 << winsizeP
00730 << "\nNumber of ml lines per buffer during computation: "
00731 << BUF
00732 << "\nTip: disfloat " << input_i_coherence.focoh
00733 << " " << nummlpixels << " 1 500"
00734 << "\n*******************************************************************\n";
00735 scratchlogfile.close();
00736
00737 ofstream scratchresfile("scratchrescoherence", ios::out | ios::trunc);
00738 bk_assert(scratchresfile,"compcoherence: scratchrescoherence",__FILE__,__LINE__);
00739 scratchresfile << "\n\n*******************************************************************"
00740
00741 << "\n*_Start_" << processcontrol[pr_i_coherence]
00742 << "\n*******************************************************************";
00743 if (!nocoh)
00744 {
00745 scratchresfile
00746 << "\nData_output_file: \t\t\t"
00747 << input_i_coherence.focoh
00748 << "\nData_output_format: \t\t\t"
00749 << "real4";
00750 }
00751 if (!noccoh)
00752 {
00753 scratchresfile
00754 << "\nData_output_file_complex_coherence: "
00755 << input_i_coherence.foccoh
00756 << "\nData_output_format_complex_coherence: \t\t"
00757 << "complex_real4";
00758 }
00759 scratchresfile
00760 << "\nFirst_line (w.r.t. original_master): \t"
00761 << slave.currentwindow.linelo
00762 << "\nLast_line (w.r.t. original_master): \t"
00763 << slave.currentwindow.linehi
00764 << "\nFirst_pixel (w.r.t. original_master): \t"
00765 << slave.currentwindow.pixlo
00766 << "\nLast_pixel (w.r.t. original_master): \t"
00767 << slave.currentwindow.pixhi
00768 << "\nMultilookfactor_azimuth_direction: \t"
00769 << multiL
00770 << "\nMultilookfactor_range_direction: \t"
00771 << multiP
00772 << "\nNumber of lines (multilooked): \t\t"
00773 << Ifilelines / multiL
00774 << "\nNumber of pixels (multilooked): \t"
00775 << nummlpixels
00776 << "\n*******************************************************************"
00777 << "\n* End_" << processcontrol[pr_i_coherence] << "_NORMAL"
00778 << "\n*******************************************************************\n";
00779 scratchresfile.close();
00780
00781
00782 PROGRESS.print("finished compcoherence.");
00783 if (!noccoh)
00784 ofileccoh.close();
00785 if (!nocoh)
00786 ofilecoh.close();
00787
00788 }
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 void subtrrefpha(
00806 const slcimage &master,
00807 const productinfo &interferogram,
00808 const input_gen &input_general,
00809 const input_subtrrefpha &input_i_subtrrefpha,
00810 const matrix<real8> &coeff_flatearth)
00811 {
00812 TRACE_FUNCTION("subtrrefpha (polynomial) (BK 09-Feb-2000)")
00813
00814 const int32 multiL = input_i_subtrrefpha.multilookL;
00815 const int32 multiP = input_i_subtrrefpha.multilookP;
00816
00817
00818 const int32 mldiskL = interferogram.multilookL;
00819 const int32 mldiskP = interferogram.multilookP;
00820
00821 const int32 numlinesdisk = interferogram.win.lines()/mldiskL;
00822 const int32 numpixelsdisk = interferogram.win.pixels()/mldiskP;
00823
00824 const real4 veryfirstline = real4(interferogram.win.linelo) +
00825 (real4(mldiskL) - 1.) / 2.;
00826 const real4 firstpixel = real4(interferogram.win.pixlo) +
00827 (real4(mldiskP) - 1.) / 2.;
00828
00829 const int32 totalmlL = mldiskL*multiL;
00830 const int32 totalmlP = mldiskP*multiP;
00831 const int32 numlinesoutput = numlinesdisk/multiL;
00832 const int32 numpixelsoutput = numpixelsdisk/multiP;
00833 const int32 lastlineoutput =
00834 interferogram.win.linelo + totalmlL*numlinesoutput - 1;
00835 const int32 lastpixeloutput =
00836 interferogram.win.pixlo + totalmlP*numpixelsoutput - 1;
00837
00838
00839 const real8 minL = master.originalwindow.linelo;
00840 const real8 maxL = master.originalwindow.linehi;
00841 const real8 minP = master.originalwindow.pixlo;
00842 const real8 maxP = master.originalwindow.pixhi;
00843 INFO << "subtrrefpha: polynomial ref.phase normalized by factors: "
00844 << minL << " " << maxL << " " << minP << " " << maxP
00845 << " to [-2,2]";
00846 INFO.print();
00847
00848
00849
00850 ofstream ofilecint;
00851 ofstream ofrefpha;
00852 if (input_i_subtrrefpha.dumponlyrefpha)
00853 {
00854 openfstream(ofrefpha,input_i_subtrrefpha.forefpha,input_general.overwrit);
00855 bk_assert(ofrefpha,input_i_subtrrefpha.forefpha,__FILE__,__LINE__);
00856 }
00857 else
00858 {
00859 openfstream(ofilecint,input_i_subtrrefpha.focint,input_general.overwrit);
00860 bk_assert(ofilecint,input_i_subtrrefpha.focint,__FILE__,__LINE__);
00861 }
00862
00863
00864 const uint BUFFERMEMSIZE = input_general.memory;
00865 const int32 bytesperline = 2*sizeof(real4)*numpixelsdisk;
00866
00867 ifstream ifcint;
00868 openfstream(ifcint,interferogram.file);
00869 bk_assert(ifcint,interferogram.file,__FILE__,__LINE__);
00870 if (interferogram.formatflag != FORMATCR4)
00871 {
00872 PRINT_ERROR("code: .. Complex interferogram on disk assumed to be complex real4.")
00873 throw(unhandled_case_error);
00874 }
00875
00876 const real4 nummatrices = 2;
00877 int32 bufferlines =
00878 int32((BUFFERMEMSIZE / nummatrices) / bytesperline);
00879 while (bufferlines%multiL)
00880 bufferlines -= 1;
00881 DEBUG << "bufferlines per buffer: " << bufferlines;
00882 DEBUG.print();
00883 if (bufferlines%multiL != 0)
00884 {
00885 PRINT_ERROR("panic: totally impossible on HP.")
00886 throw(some_error);
00887 }
00888
00889 const int32 numfullbuffers = numlinesdisk/bufferlines;
00890 int32 restlines = numlinesdisk%bufferlines;
00891 while (restlines%multiL)
00892 restlines -= 1;
00893 DEBUG << "restlines last buffer: " << restlines
00894 << "; and: " << numfullbuffers
00895 << " buffers of #lines: " << bufferlines;
00896 DEBUG.print();
00897
00898 if (restlines%multiL != 0)
00899 {
00900 PRINT_ERROR("panic: i thought restlines is always exactly fitted this way??.")
00901 throw(some_error);
00902 }
00903 const int32 EXTRABUFFER = (restlines) ? 1 : 0;
00904
00905
00906 register int32 i;
00907 matrix<real4> p_axis(numpixelsdisk,1);
00908 for (i=0; i<numpixelsdisk; i++)
00909 p_axis(i,0) = firstpixel + i*mldiskP;
00910 normalize(p_axis,minP,maxP);
00911
00912
00913 for (register int32 buffer=0; buffer<numfullbuffers+EXTRABUFFER; ++buffer)
00914 {
00915 const real4 firstline = veryfirstline + buffer*bufferlines*mldiskL;
00916 if (buffer == numfullbuffers)
00917 bufferlines = restlines;
00918
00919
00920 PROGRESS << "SUBTRREFPHA: buffer: " << buffer+1 << "; "
00921 << "line: " << firstline << " : "
00922 << firstline+bufferlines-1;
00923 PROGRESS.print();
00924
00925
00926
00927 matrix<real4> l_axis(bufferlines,1);
00928 for (i=0; i<bufferlines; i++)
00929 l_axis(i,0) = firstline + i*mldiskL;
00930 normalize(l_axis,minL,maxL);
00931
00932 matrix<real4> REFPHASE = polyval(l_axis, p_axis, coeff_flatearth);
00933
00934 if (input_i_subtrrefpha.dumponlyrefpha)
00935 {
00936 matrix<complr4> REFPHASECR4 = angle2cmplx(REFPHASE);
00937 ofrefpha << REFPHASECR4;
00938 continue;
00939 }
00940
00941
00942 matrix<complr4> CINT(bufferlines,numpixelsdisk);
00943 ifcint >> CINT;
00944
00945
00946
00947
00948 dotmultconjphase(CINT,REFPHASE);
00949
00950
00951 ofilecint << multilook(CINT, multiL, multiP);
00952 }
00953 ifcint.close();
00954 ofilecint.close();
00955
00956 if (input_i_subtrrefpha.dumponlyrefpha)
00957 ofrefpha.close();
00958
00959
00960
00961 ofstream scratchlogfile("scratchlogsubtrrefpha", ios::out | ios::trunc);
00962 bk_assert(scratchlogfile,"subtrrefpha: scratchlogsubtrrefpha",__FILE__,__LINE__);
00963 scratchlogfile
00964 << "\n\n* SUBTRREFPHA *****************************************************";
00965 if (input_i_subtrrefpha.dumponlyrefpha)
00966 {
00967 scratchlogfile
00968 << "\nOnly dump of reference phase, no subtraction."
00969 << "\nData_output_file_ref.phase: "
00970 << input_i_subtrrefpha.forefpha;
00971 }
00972 else
00973 {
00974 scratchlogfile
00975 << "\nInput file complex interferogram: \t\t\t"
00976 << interferogram.file
00977 << "\nData_output_file_complex_interferogram: "
00978 << input_i_subtrrefpha.focint
00979 << "\nData_output_format_cint: \t\t"
00980 << "complex_real4";
00981 }
00982 scratchlogfile
00983 << "\nFirst_line (w.r.t. original_master): \t"
00984 << interferogram.win.linelo
00985 << "\nLast_line (w.r.t. original_master): \t"
00986 << lastlineoutput
00987 << "\nFirst_pixel (w.r.t. original_master): \t"
00988 << interferogram.win.pixlo
00989 << "\nLast_pixel (w.r.t. original_master): \t"
00990 << lastpixeloutput
00991 << "\nMultilookfactor_azimuth_direction: \t"
00992 << totalmlL
00993 << "\nMultilookfactor_range_direction: \t"
00994 << totalmlP
00995 << "\nNumber of lines (multilooked): \t\t"
00996 << numlinesoutput
00997 << "\nNumber of pixels (multilooked): \t"
00998 << numpixelsoutput
00999 << "\n\nMultilookfactors input complex interferogram: \t"
01000 << mldiskL << " " << mldiskP
01001 << "\nMultilookfactors requested in this step: \t"
01002 << multiL << " " << multiP
01003 << "\nNumber of buffers used (size): \t"
01004 << numfullbuffers << "(" << bufferlines << "," << numpixelsdisk << ")"
01005 << "\n*******************************************************************\n\n";
01006 scratchlogfile.close();
01007
01008
01009 ofstream scratchresfile("scratchressubtrrefpha", ios::out | ios::trunc);
01010 bk_assert(scratchresfile,"subtrrefpha: scratchressubtrrefpha",__FILE__,__LINE__);
01011 scratchresfile
01012
01013 << "\n\n*******************************************************************"
01014 << "\n*_Start_" << processcontrol[pr_i_subtrrefpha]
01015 << "\n* Method: polynomial"
01016 << "\n*******************************************************************"
01017 << "\nData_output_file: \t\t\t";
01018 if (input_i_subtrrefpha.dumponlyrefpha)
01019 {
01020 scratchresfile
01021 << "NO_OUTPUT_ONLY_DUMPING_REF_PHA"
01022 << "\nFile_name of ref.phase: \t\t"
01023 << input_i_subtrrefpha.forefpha
01024 << "\nData_output_format: \t\t\t"
01025 << "complex_real4";
01026 }
01027 else
01028 {
01029 scratchresfile
01030 << input_i_subtrrefpha.focint
01031 << "\nData_output_format: \t\t\t"
01032 << "complex_real4";
01033 }
01034 scratchresfile
01035 << "\nFirst_line (w.r.t. original_master): \t"
01036 << interferogram.win.linelo
01037 << "\nLast_line (w.r.t. original_master): \t"
01038 << lastlineoutput
01039 << "\nFirst_pixel (w.r.t. original_master): \t"
01040 << interferogram.win.pixlo
01041 << "\nLast_pixel (w.r.t. original_master): \t"
01042 << lastpixeloutput
01043 << "\nMultilookfactor_azimuth_direction: \t"
01044 << totalmlL
01045 << "\nMultilookfactor_range_direction: \t"
01046 << totalmlP
01047 << "\nNumber of lines (multilooked): \t\t"
01048 << numlinesoutput
01049 << "\nNumber of pixels (multilooked): \t"
01050 << numpixelsoutput
01051 << "\n*******************************************************************"
01052 << "\n* End_" << processcontrol[pr_i_subtrrefpha] << "_NORMAL"
01053 << "\n*******************************************************************\n";
01054 scratchresfile.close();
01055
01056
01057
01058 }
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 void subtrrefpha(
01076 const input_ell &ellips,
01077 const slcimage &master,
01078 const slcimage &slave,
01079 const productinfo &interferogram,
01080 const input_gen &input_general,
01081 const input_subtrrefpha &input_i_subtrrefpha,
01082 orbit &masterorbit,
01083 orbit &slaveorbit)
01084 {
01085 TRACE_FUNCTION("subtrrefpha (BK 03-Jul-2000)")
01086 const int32 MAXITER = 10;
01087 const real8 CRITERPOS = 1e-6;
01088 const real8 CRITERTIM = 1e-10;
01089 INFO << "SUBTRREFPHA: MAXITER: " << MAXITER << "; "
01090 << "CRITERPOS: " << CRITERPOS << " m; "
01091 << "CRITERTIM: " << CRITERTIM << " s";
01092 INFO.print();
01093
01094
01095 const int32 multiL = input_i_subtrrefpha.multilookL;
01096 const int32 multiP = input_i_subtrrefpha.multilookP;
01097
01098
01099 const int32 mldiskL = interferogram.multilookL;
01100 const int32 mldiskP = interferogram.multilookP;
01101
01102 const int32 numlinesdisk = interferogram.win.lines()/mldiskL;
01103 const int32 numpixelsdisk = interferogram.win.pixels()/mldiskP;
01104
01105 const real4 veryfirstline = real4(interferogram.win.linelo) +
01106 (real4(mldiskL) - 1.) / 2.;
01107 const real4 firstpixel = real4(interferogram.win.pixlo) +
01108 (real4(mldiskP) - 1.) / 2.;
01109
01110 const int32 totalmlL = mldiskL*multiL;
01111 const int32 totalmlP = mldiskP*multiP;
01112 const int32 numlinesoutput = numlinesdisk/multiL;
01113 const int32 numpixelsoutput = numpixelsdisk/multiP;
01114 const int32 lastlineoutput =
01115 interferogram.win.linelo + totalmlL*numlinesoutput - 1;
01116 const int32 lastpixeloutput =
01117 interferogram.win.pixlo + totalmlP*numpixelsoutput - 1;
01118
01119
01120
01121 ifstream ifcint;
01122 ofstream ofilecint;
01123 ofstream ofrefpha;
01124 if (input_i_subtrrefpha.dumponlyrefpha)
01125 {
01126 openfstream(ofrefpha,input_i_subtrrefpha.forefpha,input_general.overwrit);
01127 bk_assert(ofrefpha,input_i_subtrrefpha.forefpha,__FILE__,__LINE__);
01128 }
01129 else
01130 {
01131 openfstream(ofilecint,input_i_subtrrefpha.focint,input_general.overwrit);
01132 bk_assert(ofilecint,input_i_subtrrefpha.focint,__FILE__,__LINE__);
01133
01134 if (interferogram.formatflag != FORMATCR4)
01135 {
01136 PRINT_ERROR("Complex interferogram on disk must be complex real4.")
01137 throw(unhandled_case_error);
01138 }
01139 openfstream(ifcint,interferogram.file);
01140 bk_assert(ifcint,interferogram.file,__FILE__,__LINE__);
01141 }
01142
01143
01144
01145 matrix <real4> REFPHA(multiL,numpixelsdisk);
01146 int32 tenpercent = int32(floor(numlinesoutput/10.));
01147 if (tenpercent==0) tenpercent = 1000;
01148 int32 percentage = 0;
01149 real8 line = veryfirstline;
01150 const real8 m_minpi4cdivlam = (-4*PI*SOL)/master.wavelength;
01151 const real8 s_minpi4cdivlam = (-4*PI*SOL)/slave.wavelength;
01152
01153
01154
01155 for (register int32 buffer=0; buffer<numlinesoutput; ++buffer)
01156 {
01157
01158 if (buffer%tenpercent==0)
01159 {
01160 PROGRESS << "SUBTRREFPHA: " << setw(3) << percentage << "%";
01161 PROGRESS.print();
01162 percentage += 10;
01163 }
01164
01165
01166 for (register int32 i=0; i<multiL; ++i)
01167 {
01168 for (register int32 j=0; j<numpixelsdisk; ++j)
01169 {
01170 real8 pixel = firstpixel + j*mldiskP;
01171
01172
01173
01174 const real8 m_trange = master.pix2tr(pixel);
01175
01176
01177 cn P;
01178 lp2xyz(line,pixel,ellips,master,masterorbit,
01179 P,MAXITER,CRITERPOS);
01180
01181
01182 real8 s_tazi;
01183 real8 s_trange;
01184 xyz2t(s_tazi,s_trange,slave,
01185 slaveorbit,
01186 P,MAXITER,CRITERTIM);
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199 REFPHA(i,j) = m_minpi4cdivlam*m_trange -
01200 s_minpi4cdivlam*s_trange;
01201 }
01202
01203 line += mldiskL;
01204 }
01205
01206
01207 if (input_i_subtrrefpha.dumponlyrefpha)
01208 {
01209 matrix<complr4> REFPHASECR4 = angle2cmplx(REFPHA);
01210 ofrefpha << multilook(REFPHASECR4,multiL,multiP);
01211 }
01212 else
01213 {
01214 matrix <complr4> CINT(multiL,numpixelsdisk);
01215 ifcint >> CINT;
01216
01217 dotmultconjphase(CINT,REFPHA);
01218 ofilecint << multilook(CINT,multiL,multiP);
01219 }
01220 }
01221
01222
01223
01224 if (input_i_subtrrefpha.dumponlyrefpha)
01225 ofrefpha.close();
01226 else
01227 {
01228 ofilecint.close();
01229 ifcint.close();
01230 }
01231
01232
01233
01234 ofstream scratchlogfile("scratchlogsubtrrefpha", ios::out | ios::trunc);
01235 bk_assert(scratchlogfile,"subtrrefpha: scratchlogsubtrrefpha",__FILE__,__LINE__);
01236 scratchlogfile
01237 << "\n\n* SUBTRREFPHA *****************************************************";
01238 if (input_i_subtrrefpha.dumponlyrefpha)
01239 {
01240 scratchlogfile
01241 << "\nOnly dump of reference phase, no subtraction."
01242 << "\nData_output_file_ref.phase: "
01243 << input_i_subtrrefpha.forefpha;
01244 }
01245 else
01246 {
01247 scratchlogfile
01248 << "\nInput file complex interferogram: \t\t\t"
01249 << interferogram.file
01250 << "\nData_output_file_complex_interferogram: "
01251 << input_i_subtrrefpha.focint
01252 << "\nData_output_format_complex_interferogram: "
01253 << "complex_real4";
01254 }
01255 scratchlogfile
01256 << "\nFirst_line (w.r.t. original_master): \t"
01257 << interferogram.win.linelo
01258 << "\nLast_line (w.r.t. original_master): \t"
01259 << lastlineoutput
01260 << "\nFirst_pixel (w.r.t. original_master): \t"
01261 << interferogram.win.pixlo
01262 << "\nLast_pixel (w.r.t. original_master): \t"
01263 << lastpixeloutput
01264 << "\nMultilookfactor_azimuth_direction: \t"
01265 << totalmlL
01266 << "\nMultilookfactor_range_direction: \t"
01267 << totalmlP
01268 << "\nNumber of lines (multilooked): \t\t"
01269 << numlinesoutput
01270 << "\nNumber of pixels (multilooked): \t"
01271 << numpixelsoutput
01272 << "\n\nMultilookfactors input complex interferogram: \t"
01273 << mldiskL << " " << mldiskP
01274 << "\nMultilookfactors requested in this step: \t"
01275 << multiL << " " << multiP
01276 << "\n*******************************************************************\n\n";
01277 scratchlogfile.close();
01278
01279
01280 ofstream scratchresfile("scratchressubtrrefpha", ios::out | ios::trunc);
01281 bk_assert(scratchresfile,"subtrrefpha: scratchressubtrrefpha",__FILE__,__LINE__);
01282 scratchresfile
01283
01284 << "\n\n*******************************************************************"
01285 << "\n*_Start_" << processcontrol[pr_i_subtrrefpha]
01286 << "\n* Method: exact"
01287 << "\n*******************************************************************"
01288 << "\nData_output_file: \t\t\t";
01289 if (input_i_subtrrefpha.dumponlyrefpha)
01290 {
01291 scratchresfile
01292 << "NO_OUTPUT_ONLY_DUMPING_REF_PHA"
01293 << "\nFile_name of ref.phase: \t\t"
01294 << input_i_subtrrefpha.forefpha
01295 << "\nData_output_format: \t\t\t"
01296 << "complex_real4";
01297 }
01298 else
01299 {
01300 scratchresfile
01301 << input_i_subtrrefpha.focint
01302 << "\nData_output_format: \t\t\t"
01303 << "complex_real4";
01304 }
01305 scratchresfile
01306 << "\nFirst_line (w.r.t. original_master): \t"
01307 << interferogram.win.linelo
01308 << "\nLast_line (w.r.t. original_master): \t"
01309 << lastlineoutput
01310 << "\nFirst_pixel (w.r.t. original_master): \t"
01311 << interferogram.win.pixlo
01312 << "\nLast_pixel (w.r.t. original_master): \t"
01313 << lastpixeloutput
01314 << "\nMultilookfactor_azimuth_direction: \t"
01315 << totalmlL
01316 << "\nMultilookfactor_range_direction: \t"
01317 << totalmlP
01318 << "\nNumber of lines (multilooked): \t\t"
01319 << numlinesoutput
01320 << "\nNumber of pixels (multilooked): \t"
01321 << numpixelsoutput
01322 << "\n*******************************************************************"
01323 << "\n* End_" << processcontrol[pr_i_subtrrefpha] << "_NORMAL"
01324 << "\n*******************************************************************\n";
01325 scratchresfile.close();
01326
01327
01328
01329
01330 }
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 void subtrrefdem(
01354 const productinfo &interferogram,
01355 const productinfo &radarcodedrefdem,
01356 const input_gen &input_general,
01357 const input_subtrrefdem &input_i_subtrrefdem)
01358 {
01359 TRACE_FUNCTION("subtrrefdem (BK 10-Feb-2000)")
01360 const int32 additional_offsetL = input_i_subtrrefdem.offsetL;
01361 const int32 additional_offsetP = input_i_subtrrefdem.offsetP;
01362
01363
01364
01365
01366
01367
01368
01369
01370 window refdem = radarcodedrefdem.win;
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380 refdem.linelo += (additional_offsetL * radarcodedrefdem.multilookL);
01381 refdem.linehi += (additional_offsetL * radarcodedrefdem.multilookL);
01382 refdem.pixlo += (additional_offsetP * radarcodedrefdem.multilookP);
01383 refdem.pixhi += (additional_offsetP * radarcodedrefdem.multilookP);
01384
01385
01386
01387
01388
01389 if (interferogram.multilookL != radarcodedrefdem.multilookL)
01390 {
01391 PRINT_ERROR("MultilookfactorL complex interferogram, ref. dem not equal.")
01392 throw(file_error);
01393 }
01394 if (interferogram.multilookP != radarcodedrefdem.multilookP)
01395 {
01396 PRINT_ERROR("MultilookfactorP complex interferogram, ref. dem not equal.")
01397 throw(file_error);
01398 }
01399 if ((interferogram.win.linelo-refdem.linelo)%radarcodedrefdem.multilookL != 0)
01400 WARNING.print("Seems reference phase DEM does not lie at same grid as complex interferogram.");
01401 if ((interferogram.win.pixlo-refdem.pixlo)%radarcodedrefdem.multilookP != 0)
01402 WARNING.print("Seems reference phase DEM does not lie at same grid as complex interferogram.");
01403
01404 const int32 cintfilelines =
01405 interferogram.win.lines() /interferogram.multilookL;
01406 const int32 cintfilepixels =
01407 interferogram.win.pixels()/interferogram.multilookP;
01408 const int32 demfilelines =
01409 radarcodedrefdem.win.lines() /radarcodedrefdem.multilookL;
01410 const int32 demfilepixels =
01411 radarcodedrefdem.win.pixels()/radarcodedrefdem.multilookP;
01412
01413
01414
01415 const int32 offsetL =
01416 (int32(refdem.linelo)-int32(interferogram.win.linelo))/
01417 int32(interferogram.multilookL);
01418 const int32 skiplinesstart = max(0,offsetL);
01419 const int32 skiplinesend = max(0,cintfilelines-demfilelines-offsetL);
01420 const int32 numpixoverlap =
01421 int32(min(int32(interferogram.win.pixhi),int32(refdem.pixhi)) -
01422 max(int32(interferogram.win.pixlo),int32(refdem.pixlo)) + 1 ) /
01423 int32(radarcodedrefdem.multilookP);
01424 const int32 startcintP =
01425 max(0,(int32(refdem.pixlo-interferogram.win.pixlo))/
01426 int32(radarcodedrefdem.multilookP));
01427 const int32 startrefdemP =
01428 max(0,(int32(interferogram.win.pixlo-refdem.pixlo))/
01429 int32(radarcodedrefdem.multilookP));
01430
01431 DEBUG << " skiplinesstart: " << skiplinesstart << " skiplinesend: " << skiplinesend
01432 << " offsetL: " << offsetL << " numpixoverlap: " << numpixoverlap
01433 << " startcintP: " << startcintP << " startrefdemP: " << startrefdemP;
01434 DEBUG.print();
01435
01436
01437
01438 if (interferogram.formatflag != FORMATCR4)
01439 {
01440 PRINT_ERROR("code ..: Complex interferogram on disk assumed to be complex real4.")
01441 throw(file_error);
01442 }
01443 if (radarcodedrefdem.formatflag != FORMATR4)
01444 {
01445 PRINT_ERROR("code ..: Reference phase DEM on disk assumed to be real4.")
01446 throw(file_error);
01447 }
01448
01449 ifstream ifcint;
01450 openfstream(ifcint,interferogram.file);
01451 bk_assert(ifcint,interferogram.file,__FILE__,__LINE__);
01452
01453 ifstream ifrefdem;
01454 openfstream(ifrefdem,radarcodedrefdem.file);
01455 bk_assert(ifrefdem,radarcodedrefdem.file,__FILE__,__LINE__);
01456
01457 ofstream ofilecint;
01458 openfstream(ofilecint,input_i_subtrrefdem.focint,input_general.overwrit);
01459 bk_assert(ofilecint,input_i_subtrrefdem.focint,__FILE__,__LINE__);
01460
01461
01462
01463 matrix<complr4> CINT (1,cintfilepixels);
01464 matrix<real4> REFDEM (1,demfilepixels);
01465 register int32 i,line;
01466
01467 #define COMP_COHER
01468 #ifdef COMP_COHER
01469 real8 coher = 0;
01470 #endif
01471
01472 for (line=0; line<skiplinesstart; ++line)
01473 {
01474 DEBUG.print("Copying line of interferogram since no overlap at start.");
01475 ifcint >> CINT;
01476 ofilecint << CINT;
01477 }
01478 for (line=0; line<-offsetL; ++line)
01479 ifrefdem >> REFDEM;
01480 for (line=0; line<cintfilelines-skiplinesstart-skiplinesend; ++line)
01481 {
01482
01483
01484 ifcint >> CINT;
01485 ifrefdem >> REFDEM;
01486
01487
01488
01489 for (i=0; i<numpixoverlap; ++i)
01490 CINT[0][i+startcintP] *= complr4(cos(REFDEM[0][i+startrefdemP]),
01491 -sin(REFDEM[0][i+startrefdemP]));
01492 #ifdef COMP_COHER
01493
01494
01495
01496 complr4 coher_line = 0.;
01497 real4 sum_line = 0.;
01498 for (i=0; i<numpixoverlap; ++i)
01499 {
01500 coher_line += CINT[0][i+startcintP];
01501 sum_line += abs(CINT[0][i+startcintP]);
01502 }
01503 coher += abs(coher_line)/sum_line;
01504 #endif
01505
01506
01507 ofilecint << CINT;
01508 }
01509 ifrefdem.close();
01510 #ifdef COMP_COHER
01511 coher /= (cintfilelines-skiplinesstart-skiplinesend);
01512 INFO << "coherence interferogram_synthetic phase = " << coher;
01513 INFO.print();
01514 #endif
01515 for (line=0; line<skiplinesend; ++line)
01516 {
01517 DEBUG.print("Copying line of interferogram since no overlap at end.");
01518 ifcint >> CINT;
01519 ofilecint << CINT;
01520 }
01521 ifcint.close();
01522 ofilecint.close();
01523
01524
01525
01526 ofstream scratchlogfile("scratchlogsubtrrefdem", ios::out | ios::trunc);
01527 bk_assert(scratchlogfile,"subtrrefdem: scratchlogsubtrrefdem",__FILE__,__LINE__);
01528 scratchlogfile
01529 << "\n\n*******************************************************************"
01530 << "\n*_Start_subtrrefdem"
01531 << "\nInput file complex interferogram: \t\t\t"
01532 << interferogram.file
01533 << "\nInput file reference dem: \t\t\t"
01534 << radarcodedrefdem.file
01535 << "\nAdditional_azimuth_shift: \t"
01536 << additional_offsetL
01537 << "\nAdditional_range_shift: \t"
01538 << additional_offsetP
01539 << "\nCoherence IFG DEMPHASE: \t"
01540 << coher
01541 << "\n*******************************************************************\n\n";
01542 scratchlogfile.close();
01543
01544 ofstream scratchresfile("scratchressubtrrefdem", ios::out | ios::trunc);
01545 bk_assert(scratchresfile,"subtrrefdem: scratchressubtrrefdem",__FILE__,__LINE__);
01546 scratchresfile
01547 << "\n\n*******************************************************************"
01548 << "\n*_Start_" << processcontrol[pr_i_subtrrefdem]
01549 << "\n*******************************************************************"
01550 << "\nMethod: \t"
01551 << "NOT_USED"
01552 << "\nAdditional_azimuth_shift: \t"
01553 << additional_offsetL
01554 << "\nAdditional_range_shift: \t"
01555 << additional_offsetP
01556 << "\nData_output_file: \t"
01557 << input_i_subtrrefdem.focint
01558 << "\nData_output_format: \t"
01559 << "complex_real4"
01560 << "\nFirst_line (w.r.t. original_master): \t"
01561 << interferogram.win.linelo
01562 << "\nLast_line (w.r.t. original_master): \t"
01563 << interferogram.win.linehi
01564 << "\nFirst_pixel (w.r.t. original_master): \t"
01565 << interferogram.win.pixlo
01566 << "\nLast_pixel (w.r.t. original_master): \t"
01567 << interferogram.win.pixhi
01568 << "\nMultilookfactor_azimuth_direction: \t"
01569 << interferogram.multilookL
01570 << "\nMultilookfactor_range_direction: \t"
01571 << interferogram.multilookP
01572 << "\nNumber of lines (multilooked): \t"
01573 << cintfilelines
01574 << "\nNumber of pixels (multilooked): \t"
01575 << cintfilepixels
01576 << "\n*******************************************************************"
01577 << "\n* End_" << processcontrol[pr_i_subtrrefdem] << "_NORMAL"
01578 << "\n*******************************************************************\n";
01579 scratchresfile.close();
01580 }
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619 void dinsar(
01620 const input_gen &input_general,
01621 const input_dinsar &dinsarinput,
01622 const input_ell &ellips,
01623 const slcimage &master,
01624 orbit &masterorbit,
01625 const slcimage &defoslave,
01626 orbit &defoorbit,
01627 const productinfo &defointerferogram
01628 )
01629 {
01630 TRACE_FUNCTION("dinsar (BK 22-Sep-2000)")
01631
01632 char difffile[ONE27];
01633 strcpy(difffile,dinsarinput.fodinsar);
01634
01635
01636
01637 slcimage topomaster;
01638 slcimage toposlave;
01639 productinfo topounwrappedinterf;
01640 orbit topomasterorbit;
01641 orbit toposlaveorbit;
01642
01643
01644 bool FOURPASS = true;
01645 if (!specified(dinsarinput.topomasterresfile))
01646 FOURPASS = false;
01647 if (!strcmp(dinsarinput.topomasterresfile,input_general.m_resfile))
01648 FOURPASS = false;
01649
01650 if (FOURPASS==true)
01651 {
01652 INFO.print("Using 4 pass differential interferometry (card IN_TOPOMASTER).");
01653 INFO.print("Reading primary (topography pair) master parameters.");
01654 topomaster.fillslcimage(dinsarinput.topomasterresfile);
01655 INFO.print("Modelling primary (topography pair) master orbit.");
01656 topomasterorbit.initialize(dinsarinput.topomasterresfile);
01657 }
01658
01659 INFO.print("Reading primary (topography pair) slave parameters.");
01660 toposlave.fillslcimage(dinsarinput.toposlaveresfile);
01661 INFO.print("Reading (topography pair) unwrapped section.");
01662 char SECTIONID[ONE27];
01663 strcpy(SECTIONID,"*_Start_");
01664 strcat(SECTIONID,processcontrol[pr_i_unwrap]);
01665 topounwrappedinterf.fillproductinfo(dinsarinput.topointresfile,SECTIONID);
01666 INFO.print("Modelling primary (topography pair) slave orbit.");
01667 toposlaveorbit.initialize(dinsarinput.toposlaveresfile);
01668
01669
01670
01671 if (defointerferogram.multilookL != topounwrappedinterf.multilookL)
01672 {
01673 WARNING << "multilookfactor defo != factor topo: " << defointerferogram.multilookL
01674 << " != " << topounwrappedinterf.multilookL;
01675 WARNING.print();
01676 }
01677 if (defointerferogram.multilookP != topounwrappedinterf.multilookP)
01678 {
01679 WARNING << "multilookfactor defo != factor topo: " << defointerferogram.multilookP
01680 << " != " << topounwrappedinterf.multilookP;
01681 WARNING.print();
01682 }
01683 if (defointerferogram.win != topounwrappedinterf.win)
01684 WARNING.print("window defo pair != window topo pair");
01685
01686 if (defointerferogram.formatflag != FORMATCR4)
01687 {
01688 PRINT_ERROR("format defo interferogram not complr4.")
01689 throw(file_error);
01690 }
01691
01692
01693
01694
01695
01696
01697
01698 const real8 minL = master.originalwindow.linelo;
01699 const real8 maxL = master.originalwindow.linehi;
01700 const real8 minP = master.originalwindow.pixlo;
01701 const real8 maxP = master.originalwindow.pixhi;
01702 INFO << "dinsar: polynomial for ratio normalized by: "
01703 << minL << " " << maxL << " " << minP << " " << maxP
01704 << " to [-2,2]";
01705 INFO.print();
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727 const int32 numpointsL = 20;
01728 const int32 numpointsP = 10;
01729 real8 dlines = (topounwrappedinterf.win.linehi-topounwrappedinterf.win.linelo) /
01730 (numpointsL-1);
01731 real8 dpixels = (topounwrappedinterf.win.pixhi -topounwrappedinterf.win.pixlo) /
01732 (numpointsP-1);
01733 INFO << "Computing baseline on grid ("
01734 << topounwrappedinterf.win.linelo << ":" << dlines << ":"
01735 << topounwrappedinterf.win.linehi << ","
01736 << topounwrappedinterf.win.pixlo << ":" << dpixels << ":"
01737 << topounwrappedinterf.win.pixhi
01738 << ") = (" << numpointsL << "x" << numpointsP << ")";
01739 INFO.print();
01740
01741 matrix<real8> LINENUMBER(numpointsL*numpointsP,1);
01742 matrix<real8> PIXELNUMBER(numpointsL*numpointsP,1);
01743 int32 i,j;
01744 int32 k=0;
01745 for (i=0; i<numpointsL; ++i)
01746 {
01747 for (j=0; j<numpointsP; ++j)
01748 {
01749 LINENUMBER(k,0) = topounwrappedinterf.win.linelo + i*dlines;
01750 PIXELNUMBER(k,0) = topounwrappedinterf.win.pixlo + j*dpixels;
01751 ++k;
01752 }
01753 }
01754
01755
01756 cn masterpos[numpointsL*numpointsP];
01757 cn pointpos[numpointsL*numpointsP];
01758 cn defoslavepos[numpointsL*numpointsP];
01759 cn topomasterpos[numpointsL*numpointsP];
01760 cn toposlavepos[numpointsL*numpointsP];
01761
01762 real8 lastline = -1.;
01763 real8 B,Bpar,Bperp;
01764 matrix<real8> Bperptopo(LINENUMBER.lines(),1);
01765 matrix<real8> Bperpdefo(LINENUMBER.lines(),1);
01766 for (i=0; i<LINENUMBER.lines(); ++i)
01767 {
01768 real8 line = LINENUMBER(i,0);
01769 real8 pixel = PIXELNUMBER(i,0);
01770 if (line == lastline)
01771 {
01772 masterpos[i] = masterpos[i-1];
01773 }
01774 else
01775 {
01776 lastline = line;
01777 const real8 m_tazi = master.line2ta(line);
01778 masterpos[i] = masterorbit.getxyz(m_tazi);
01779 }
01780
01781 lp2xyz(line,pixel,ellips,master,masterorbit,pointpos[i]);
01782 xyz2orb(defoslavepos[i],defoslave,defoorbit,pointpos[i]);
01783 xyz2orb(toposlavepos[i],toposlave,toposlaveorbit,pointpos[i]);
01784 if (FOURPASS==true)
01785 xyz2orb(topomasterpos[i],topomaster,topomasterorbit,pointpos[i]);
01786 else
01787 topomasterpos[i] = masterpos[i];
01788
01789
01790 BBparBperp(B,Bpar,Bperp,topomasterpos[i],pointpos[i],toposlavepos[i]);
01791 Bperptopo(i,0) = Bperp;
01792 BBparBperp(B,Bpar,Bperp,masterpos[i],pointpos[i],defoslavepos[i]);
01793 Bperpdefo(i,0) = Bperp;
01794 }
01795
01796
01797
01798
01799 matrix <real8> Ratio = Bperpdefo/Bperptopo;
01800
01801 #ifdef __DEBUG
01802 DEBUG.print("data dump to files: LINENUMBER, PIXELNUMBER, Bperptopo, Bperpdefo, Ratio");
01803 dumpasc("LINENUMBER",LINENUMBER);
01804 dumpasc("PIXELNUMBER",PIXELNUMBER);
01805 dumpasc("Bperptopo",Bperptopo);
01806 dumpasc("Bperpdefo",Bperpdefo);
01807 dumpasc("Ratio",Ratio);
01808 #endif
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818 matrix<real8> A(Ratio.lines(),3);
01819 for (i=0; i<A.lines(); ++i)
01820 {
01821 A(i,0) = 1.;
01822 A(i,1) = normalize(LINENUMBER(i,0),minL,maxL);
01823 A(i,2) = normalize(PIXELNUMBER(i,0),minP,maxP);
01824 }
01825 matrix<real8> N = matTxmat(A,A);
01826 matrix<real8> rhs = matTxmat(A,Ratio);
01827
01828
01829 matrix<real8> Qx_hat = N;
01830 choles(Qx_hat);
01831 solvechol(Qx_hat,rhs);
01832 invertchol(Qx_hat);
01833
01834
01835
01836 for (i=0; i<Qx_hat.lines(); i++)
01837 for (j=0; j<i; j++)
01838 Qx_hat(j,i) = Qx_hat(i,j);
01839 const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
01840 INFO << "dinsar: max(abs(N*inv(N)-I)) = " << maxdev;
01841 INFO.print();
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858 if (maxdev > .01)
01859 {
01860 ERROR << ". Too large, normalization factors <-> crop?";
01861 PRINT_ERROR(ERROR.get_str())
01862 throw(some_error);
01863 }
01864 else if (maxdev > .001)
01865 {
01866 WARNING.print("Deviation quite large. careful!");
01867 }
01868
01869
01870
01871
01872 matrix<real8> y_hat = A * rhs;
01873 matrix<real8> e_hat = Ratio - y_hat;
01874
01875 uint pos,dummy;
01876 const real8 maxerrorratio = max(abs(e_hat),pos,dummy);
01877 const real8 maxrelerror = 100. * maxerrorratio / Ratio(pos,0);
01878 INFO << "maximum error for l,p: " << LINENUMBER(pos,0) << ","
01879 << PIXELNUMBER(pos,0) << "; Ratio=" << Ratio(pos,0)
01880 << " estimate=" << y_hat(pos,0) << "; rel. err=" << maxrelerror << "%. ";
01881 INFO.print();
01882 if (maxrelerror < 5 )
01883 {
01884 INFO.print("max err OK");
01885 }
01886 else
01887 {
01888 WARNING.print("max err quite large");
01889 WARNING.print("Error in deformation vector larger than 5% due to mismodeling baseline!");
01890 }
01891
01892
01893
01894
01895 ofstream ofcint,ofscaledtopo;
01896 openfstream(ofcint,dinsarinput.fodinsar,input_general.overwrit);
01897 bk_assert(ofcint,dinsarinput.fodinsar,__FILE__,__LINE__);
01898
01899 bool writescaledtopo=false;
01900 if (specified(dinsarinput.foscaleduint))
01901 {
01902 INFO.print("writing scaled version of unwrapped topo interferogram in real4 format.");
01903 writescaledtopo=true;
01904 openfstream(ofscaledtopo,dinsarinput.foscaleduint,input_general.overwrit);
01905 bk_assert(ofscaledtopo,dinsarinput.foscaleduint,__FILE__,__LINE__);
01906 }
01907
01908 const int32 numlines = defointerferogram.win.lines() / defointerferogram.multilookL;
01909 const int32 numpixels = defointerferogram.win.pixels() / defointerferogram.multilookP;
01910 const real4 firstline = defointerferogram.win.linelo +
01911 (real8(defointerferogram.multilookL-1.)/2.);
01912 const real4 firstpixel= defointerferogram.win.pixlo +
01913 (real8(defointerferogram.multilookP-1.)/2.);
01914
01915 matrix<real4> ratioline(1,numpixels);
01916 for (i=0; i<numpixels; ++i)
01917 ratioline(0,i) = firstpixel+i*defointerferogram.multilookP;
01918 normalize(ratioline,real4(minP),real4(maxP));
01919 ratioline *= real4(rhs(2,0));
01920 ratioline += real4(rhs(0,0));
01921
01922
01923 ifstream ifdefocint;
01924 openfstream(ifdefocint,defointerferogram.file);
01925 bk_assert(ifdefocint,defointerferogram.file,__FILE__,__LINE__);
01926 matrix<complr4> DEFO(1,numpixels);
01927
01928
01929
01930
01931
01932
01933 int32 tenpercent = int32(floor(numlines/10.));
01934 if (tenpercent==0) tenpercent = 1000;
01935 int32 percent = 0;
01936 const complr4 CR4ZERO = complr4(0.,0.);
01937 for (i=0; i<numlines; ++i)
01938 {
01939 if (i%tenpercent==0)
01940 {
01941 PROGRESS << "DINSAR: " << setw(3) << percent << "%";
01942 PROGRESS.print();
01943 percent += 10;
01944 }
01945
01946 const real4 line = firstline + i*defointerferogram.multilookL;
01947 matrix<real4> ratio = ratioline +
01948 real4(rhs(1,0))*normalize(line,real4(minL),real4(maxL));
01949
01950
01951 ifdefocint >> DEFO;
01952 const window filewin(i+1,i+1,1,numpixels);
01953 matrix<real4> TOPO = topounwrappedinterf.readphase(filewin);
01954
01955
01956
01957
01958
01959
01960
01961
01962 for (j=0; j<numpixels; ++j)
01963 {
01964 (TOPO(0,j)==NaN) ?
01965 DEFO(0,j) = CR4ZERO :
01966 DEFO(0,j) *= complr4(cos(ratio(0,j)*TOPO(0,j)),-sin(ratio(0,j)*TOPO(0,j)));
01967 }
01968 ofcint << DEFO;
01969
01970
01971 if (writescaledtopo)
01972 {
01973 if (!(i%1000))
01974 PROGRESS.print("Writing scaled topo interferogram to file.");
01975 TOPO *= ratio;
01976 ofscaledtopo << TOPO;
01977 }
01978 }
01979 ifdefocint.close();
01980 ofcint.close();
01981 if (writescaledtopo)
01982 ofscaledtopo.close();
01983
01984
01985
01986 ofstream scratchlogfile("scratchlogdinsar", ios::out | ios::trunc);
01987 bk_assert(scratchlogfile,"dinsar: scratchlogdinsar",__FILE__,__LINE__);
01988 scratchlogfile
01989 << "\n\n*******************************************************************"
01990 << "\n*_Start_" << processcontrol[pr_i_dinsar]
01991 << "\nDegree 2d polynomial for baseline modeling: "
01992 << "1"
01993
01994 << "\nEstimated coefficients:\n"
01995 << "\nx_hat \tstd:\n";
01996 for (i=0; i<rhs.lines(); ++i)
01997 scratchlogfile << rhs(i,0) << " \t" << sqrt(Qx_hat(i,i)) << endl;
01998
01999 scratchlogfile << "\nQx_hat:\n";
02000 for (i=0; i<Qx_hat.lines(); i++)
02001 {
02002 for (j=0; j<Qx_hat.pixels(); j++)
02003 {
02004 scratchlogfile << Qx_hat(i,j) << " ";
02005 }
02006 scratchlogfile << endl;
02007 }
02008 scratchlogfile << "\nMaximum deviation N*inv(N):"
02009 << maxdev
02010 << "\nBaseline for pixel: "
02011 << LINENUMBER(0,0) << ", " << PIXELNUMBER(0,0) << ":"
02012 << "\n Bperp_defo: " << Bperpdefo(0,0)
02013 << "\n Bperp_topo: " << Bperptopo(0,0) << endl
02014 << "\nSome more info for each observation:"
02015 << "\nline \t pixel \t obs \t obs_hat \t err_hat\n";
02016 for (i=0; i<LINENUMBER.lines(); i++)
02017 scratchlogfile << LINENUMBER(i,0) << "\t" << PIXELNUMBER(i,0) << "\t"
02018 << Ratio(i,0) << "\t" << y_hat(i,0) << "\t" << e_hat(i,0) << "\n";
02019 scratchlogfile
02020 << "\nMaximum absolute error: \t"
02021 << maxerrorratio << " = " << maxrelerror << "%"
02022 << "\ninput filename topo interferogram: "
02023 << defointerferogram.file
02024 << "\ninput filename defo interferogram: "
02025 << topounwrappedinterf.file
02026 << "\noutput filename topocorrected defo: "
02027 << dinsarinput.fodinsar;
02028 if (writescaledtopo)
02029 scratchlogfile
02030 << "\noutput filename scaled topo interferogram: "
02031 << dinsarinput.foscaleduint;
02032 scratchlogfile
02033 << "\n*******************************************************************\n";
02034 scratchlogfile.close();
02035
02036 ofstream scratchresfile("scratchresdinsar", ios::out | ios::trunc);
02037 bk_assert(scratchresfile,"dinsar: scratchresdinsar",__FILE__,__LINE__);
02038 scratchresfile.setf(ios::scientific, ios::floatfield);
02039 scratchresfile.setf(ios::right, ios::adjustfield);
02040 scratchresfile.precision(8);
02041 scratchresfile.width(18);
02042 scratchresfile
02043 << "\n\n*******************************************************************"
02044 << "\n*_Start_" << processcontrol[pr_i_dinsar]
02045 << "\n*******************************************************************"
02046 << "\nData_output_file: \t"
02047 << dinsarinput.fodinsar
02048 << "\nData_output_format: \t"
02049 << "complex_real4"
02050 << "\nFirst_line (w.r.t. original_master): \t"
02051 << defointerferogram.win.linelo
02052 << "\nLast_line (w.r.t. original_master): \t"
02053 << defointerferogram.win.linehi
02054 << "\nFirst_pixel (w.r.t. original_master): \t"
02055 << defointerferogram.win.pixlo
02056 << "\nLast_pixel (w.r.t. original_master): \t"
02057 << defointerferogram.win.pixhi
02058 << "\nMultilookfactor_azimuth_direction: \t"
02059 << defointerferogram.multilookL
02060 << "\nMultilookfactor_range_direction: \t"
02061 << defointerferogram.multilookP
02062 << "\nNumber of lines (multilooked): \t"
02063 << numlines
02064 << "\nNumber of pixels (multilooked): \t"
02065 << numpixels
02066 << "\n*******************************************************************"
02067 << "\n* End_" << processcontrol[pr_i_dinsar] << "_NORMAL"
02068 << "\n*******************************************************************\n";
02069 scratchresfile.close();
02070 }
02071
02072