00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "constants.hh"
00046 #include "ioroutines.hh"
00047 #include "utilities.hh"
00048 #include "coregistration.hh"
00049 #include "conversion.hh"
00050 #include "refsystems.hh"
00051 #include "orbitbk.hh"
00052 #include "slcimage.hh"
00053 #include "productinfo.hh"
00054 #include "exceptions.hh"
00055 #include "bk_baseline.hh"
00056
00057
00058 #include <iomanip>
00059 #include <cstdlib>
00060 #include <cmath>
00061 #include <algorithm>
00062 #include <cstdio>
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 void coarseporbit(
00083 const input_ell &ell,
00084 const slcimage &master,
00085 const slcimage &slave,
00086 orbit &masterorbit,
00087 orbit &slaveorbit,
00088 const BASELINE &baseline)
00089 {
00090 TRACE_FUNCTION("coarseporbit (BK 12-Dec-1998)");
00091 const int16 MAXITER = 10;
00092 const real8 CRITERPOS = 1e-6;
00093 const real8 CRITERTIM = 1e-10;
00094
00095
00096
00097
00098
00099
00100 const uint cen_lin = (master.currentwindow.linelo+master.currentwindow.linehi)/2;
00101 const uint cen_pix = (master.currentwindow.pixlo +master.currentwindow.pixhi) /2;
00102 const real8 HEI = 0.0;
00103
00104
00105
00106 cn P;
00107 const int32 lp2xyziter =
00108 lp2xyz(cen_lin,cen_pix,ell,master,masterorbit,P,MAXITER,CRITERPOS);
00109
00110
00111 real8 lin,pix;
00112 const int32 xyz2lpiter =
00113 xyz2lp(lin,pix,slave,slaveorbit,P,MAXITER,CRITERTIM);
00114
00115
00116 const int Bt = Btemp(master.utc1,slave.utc1);
00117
00118 const real8 Bperp = baseline.get_bperp(cen_lin,cen_pix,HEI);
00119 const real8 Bpar = baseline.get_bpar(cen_lin,cen_pix,HEI);
00120 const real8 theta = rad2deg(baseline.get_theta(cen_lin,cen_pix,HEI));
00121 const real8 inc_angle = rad2deg(baseline.get_theta_inc(cen_lin,cen_pix,HEI));
00122
00123 const real8 B = baseline.get_b(cen_lin,cen_pix,HEI);
00124 const real8 alpha = rad2deg(baseline.get_alpha(cen_lin,cen_pix,HEI));
00125 const real8 Bh = baseline.get_bhor(cen_lin,cen_pix,HEI);
00126 const real8 Bv = baseline.get_bvert(cen_lin,cen_pix,HEI);
00127
00128 const real8 Hamb = baseline.get_hamb(cen_lin,cen_pix,HEI);
00129 const real8 orb_conv = rad2deg(baseline.get_orb_conv(cen_lin,cen_pix,HEI));
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 INFO << "Estimated translation (l,p): "
00156 << floor(lin-cen_lin +.5) << ", "
00157 << floor(pix-cen_pix +.5);
00158 INFO.print();
00159
00160
00161 ofstream scratchlogfile("scratchlogcoarse", ios::out | ios::trunc);
00162 bk_assert(scratchlogfile,"coarseporbit: scratchlogcoarse",__FILE__,__LINE__);
00163 scratchlogfile << "\n\n*******************************************************************"
00164 << "\n* COARSE_COREGISTRATION Orbits"
00165 << "\n*******************************************************************"
00166 << "\n(Approximate) center master (line,pixel,hei): "
00167 << cen_lin << ", " << cen_pix << ", " << HEI
00168 << "\nEllipsoid WGS84 coordinates of this pixel (x,y,z): ("
00169 << P.x << ", " << P.y << ", " << P.z << ")"
00170 << "\n(line,pixel) of these coordinates in slave: "
00171 << lin << ", " << pix
00172 << "\nEstimated translation slave w.r.t. master (l,p):"
00173 << floor(lin-cen_lin +.5)
00174 << ", "
00175 << floor(pix-cen_pix +.5)
00176 << "\nMaximum number of iterations: "
00177 << MAXITER
00178 << "\nCriterium for position (m): "
00179 << CRITERPOS
00180 << "\nCriterium for azimuth time (s): "
00181 << CRITERTIM
00182 << " (=~ " << CRITERTIM*7.e3 << "m)"
00183 << "\nNumber of iterations conversion line,pixel to xyz: "
00184 << lp2xyziter
00185 << "\nNumber of iterations conversion xyz to line,pixel: "
00186 << xyz2lpiter
00187 << "\n*******************************************************************\n";
00188 scratchlogfile.close();
00189
00190
00191 ofstream scratchresfile("scratchrescoarse", ios::out | ios::trunc);
00192 bk_assert(scratchresfile,"coarseporbit: scratchrescoarse",__FILE__,__LINE__);
00193 scratchresfile.setf(ios::right, ios::adjustfield);
00194 scratchresfile
00195 << "\n\n*******************************************************************"
00196 << "\n*_Start_" << processcontrol[pr_i_coarse]
00197 << "\n*******************************************************************"
00198 << "\nSome info for pixel: " << cen_lin << ", " << cen_pix << " (not used):"
00199 << "\n Btemp: [days]: "
00200 << setw(10)
00201 << setiosflags(ios::right)
00202 << Bt << " \t// Temporal baseline"
00203 << "\n Bperp [m]: "
00204 << setw(10)
00205 << setiosflags(ios::right)
00206 << onedecimal(Bperp) << " \t// Perpendicular baseline"
00207 << "\n Bpar [m]: "
00208 << setw(10)
00209 << setiosflags(ios::right)
00210 << onedecimal(Bpar) << " \t// Parallel baseline"
00211 << "\n Bh [m]: "
00212 << setw(10)
00213 << setiosflags(ios::right)
00214 << onedecimal(Bh) << " \t// Horizontal baseline"
00215 << "\n Bv [m]: "
00216 << setw(10)
00217 << setiosflags(ios::right)
00218 << onedecimal(Bv) << " \t// Vertical baseline"
00219 << "\n B [m]: "
00220 << setw(10)
00221 << setiosflags(ios::right)
00222 << onedecimal(B) << " \t// Baseline (distance between sensors)"
00223 << "\n alpha [deg]: "
00224 << setw(10)
00225 << setiosflags(ios::right)
00226 << onedecimal(alpha) << " \t// Baseline orientation"
00227 << "\n theta [deg]: "
00228 << setw(10)
00229 << setiosflags(ios::right)
00230 << onedecimal(theta) << " \t// look angle"
00231 << "\n inc_angle [deg]: "
00232 << setw(10)
00233 << setiosflags(ios::right)
00234 << onedecimal(inc_angle) << " \t// incidence angle"
00235 << "\n orbitconv [deg]: "
00236 << setw(10)
00237 << setiosflags(ios::right)
00238 << orb_conv << " \t// angle between orbits"
00239 << "\n Height_amb [m]: "
00240 << setw(10)
00241 << setiosflags(ios::right)
00242 << onedecimal(Hamb) << " \t// height = h_amb*phase/2pi (approximately)"
00243 << "\nEstimated translation slave w.r.t. master (slave-master):"
00244 << "\nCoarse_orbits_translation_lines: \t"
00245 << floor(lin-cen_lin +.5)
00246 << "\nCoarse_orbits_translation_pixels: \t"
00247 << floor(pix-cen_pix +.5)
00248 << "\n*******************************************************************"
00249 << "\n* End_" << processcontrol[pr_i_coarse] << "_NORMAL"
00250 << "\n*******************************************************************\n";
00251
00252
00253 scratchresfile.close();
00254 PROGRESS.print("Coarse precise orbits coregistration finished.");
00255 }
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 void coarsecorrel(
00276 const input_coarsecorr &coarsecorrinput,
00277 const slcimage &minfo,
00278 const slcimage &sinfo)
00279 {
00280 TRACE_FUNCTION("coarsecorrel (BK 12-Dec-1998)");
00281
00282 char dummyline[ONE27];
00283 const uint Mfilelines = minfo.currentwindow.lines();
00284 const uint Sfilelines = sinfo.currentwindow.lines();
00285
00286 const uint Nwin = coarsecorrinput.Nwin;
00287 const int32 initoffsetL = coarsecorrinput.initoffsetL;
00288 const int32 initoffsetP = coarsecorrinput.initoffsetP;
00289 uint MasksizeL = coarsecorrinput.MasksizeL;
00290 uint MasksizeP = coarsecorrinput.MasksizeP;
00291 const uint AccL = coarsecorrinput.AccL;
00292 const uint AccP = coarsecorrinput.AccP;
00293 bool pointsrandom = true;
00294 if (specified(coarsecorrinput.ifpositions))
00295 pointsrandom = false;
00296
00297
00298
00299
00300
00301 bool forceoddl = false;
00302 bool forceoddp = false;
00303 if (!isodd(MasksizeL))
00304 {
00305 forceoddl = true;
00306 MasksizeL+=1;
00307 }
00308 if (!isodd(MasksizeP))
00309 {
00310 forceoddp = true;
00311 MasksizeP+=1;
00312 }
00313
00314
00315
00316 const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
00317 const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
00318 const int32 sp0 = sinfo.currentwindow.pixlo - initoffsetP;
00319 const int32 spN = sinfo.currentwindow.pixhi - initoffsetP;
00320
00321
00322
00323
00324
00325
00326
00327
00328 const uint FEW=10;
00329 const uint l0 = uint(max(int32(minfo.currentwindow.linelo),sl0) + .5*MasksizeL + AccL + FEW);
00330 const uint lN = uint(min(int32(minfo.currentwindow.linehi),slN) - .5*MasksizeL - AccL - FEW);
00331 const uint p0 = uint(max(int32(minfo.currentwindow.pixlo),sp0) + .5*MasksizeP + AccP + FEW);
00332 const uint pN = uint(min(int32(minfo.currentwindow.pixhi),spN) - .5*MasksizeP - AccP - FEW);
00333 const window overlap(l0,lN,p0,pN);
00334
00335
00336
00337
00338
00339 register int32 i;
00340 matrix<uint> Centers;
00341 if (pointsrandom)
00342 {
00343 Centers = distributepoints(real4(Nwin),overlap);
00344 }
00345
00346 else
00347 {
00348 Centers.resize(Nwin,3);
00349
00350 ifstream ifpos;
00351 openfstream(ifpos,coarsecorrinput.ifpositions);
00352 bk_assert(ifpos,coarsecorrinput.ifpositions,__FILE__,__LINE__);
00353 uint ll,pp;
00354 for (i=0; i<Nwin; ++i)
00355 {
00356 ifpos >> ll >> pp;
00357 Centers(i,0) = uint(ll);
00358 Centers(i,1) = uint(pp);
00359 Centers(i,2) = uint(1);
00360 ifpos.getline(dummyline,ONE27,'\n');
00361 }
00362 ifpos.close();
00363
00364
00365 if (Centers(Nwin-1,0) == Centers(Nwin-2,0) &&
00366 Centers(Nwin-1,1) == Centers(Nwin-2,1))
00367 {
00368 Centers(Nwin-1,0) = uint(.5*(lN + l0) + 27);
00369 Centers(Nwin-1,1) = uint(.5*(pN + p0) + 37);
00370 WARNING << "CC: there should be no EOL after last point in file: "
00371 << coarsecorrinput.ifpositions;
00372 WARNING.print();
00373 }
00374
00375
00376
00377 bool troubleoverlap = false;
00378 for (i=0; i<Nwin; ++i)
00379 {
00380 if (Centers(i,0) < l0)
00381 {
00382 troubleoverlap=true;
00383 WARNING << "COARSE_CORR: point from file: "
00384 << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
00385 << " outside overlap master, slave. New position: ";
00386 Centers(i,0) = l0 + l0-Centers(i,0);
00387 WARNING << Centers(i,0) << " " << Centers(i,1);
00388 WARNING.print();
00389 }
00390 if (Centers(i,0) > lN)
00391 {
00392 troubleoverlap=true;
00393 WARNING << "COARSE_CORR: point from file: "
00394 << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
00395 << " outside overlap master, slave. New position: ";
00396 Centers(i,0) = lN + lN-Centers(i,0);
00397 WARNING << Centers(i,0) << " " << Centers(i,1);
00398 WARNING.print();
00399 }
00400 if (Centers(i,1) < p0)
00401 {
00402 troubleoverlap=true;
00403 WARNING << "COARSE_CORR: point from file: "
00404 << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
00405 << " outside overlap master, slave. New position: ";
00406 Centers(i,1) = p0 + p0-Centers(i,1);
00407 WARNING << Centers(i,0) << " " << Centers(i,1);
00408 WARNING.print();
00409 }
00410 if (Centers(i,1) > pN)
00411 {
00412 troubleoverlap=true;
00413 WARNING << "COARSE_CORR: point from file: "
00414 << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
00415 << " outside overlap master, slave. New position: ";
00416 Centers(i,1) = pN + pN-Centers(i,1);
00417 WARNING << Centers(i,0) << " " << Centers(i,1);
00418 WARNING.print();
00419 }
00420 }
00421 if (troubleoverlap)
00422 {
00423 WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
00424 << l0 << " " << lN << " " << p0 << " " << pN << ends;
00425 WARNING.print();
00426 }
00427 }
00428
00429
00430
00431
00432
00433 matrix<complr4> Mcmpl;
00434 matrix<complr4> Scmpl;
00435 matrix<real4> Master;
00436 matrix<real4> Mask;
00437 matrix<real4> Correl;
00438 matrix<real4> Result(Nwin,3);
00439
00440
00441
00442
00443 int32 percent = 0;
00444
00445 int32 tenpercent = int32((Nwin/10)+.5);
00446 if (tenpercent==0) tenpercent = 1000;
00447 for (i=0;i<Nwin;i++)
00448 {
00449
00450 if (i%tenpercent==0)
00451 {
00452
00453 PROGRESS << "COARSE_CORR: " << setw(3) << percent << "%" << ends;
00454 PROGRESS.print();
00455 percent += 10;
00456 }
00457
00458
00459 uint cenMwinL = Centers(i,0);
00460 uint cenMwinP = Centers(i,1);
00461
00462 window master;
00463 master.linelo = cenMwinL - (MasksizeL-1)/2 -AccL;
00464 master.linehi = master.linelo + MasksizeL +2*AccL - 1;
00465 master.pixlo = cenMwinP - (MasksizeP-1)/2 - AccP;
00466 master.pixhi = master.pixlo + MasksizeP +2*AccP - 1;
00467
00468
00469 window slavemask;
00470 uint cenSwinL = cenMwinL + initoffsetL;
00471 uint cenSwinP = cenMwinP + initoffsetP;
00472 slavemask.linelo = cenSwinL - (MasksizeL-1)/2;
00473 slavemask.linehi = slavemask.linelo + MasksizeL - 1;
00474 slavemask.pixlo = cenSwinP - (MasksizeP-1)/2;
00475 slavemask.pixhi = slavemask.pixlo + MasksizeP - 1;
00476
00477
00478
00479
00480 Mcmpl = minfo.readdata(master);
00481 Scmpl = sinfo.readdata(slavemask);
00482 Master = magnitude(Mcmpl);
00483 Mask = magnitude(Scmpl);
00484
00485
00486
00487
00488 Correl = correlate(Master,Mask);
00489 uint L, P;
00490 real4 corr = max(Correl, L, P);
00491
00492 uint relcenML = master.linehi - cenMwinL;
00493 uint relcenMP = master.pixhi - cenMwinP;
00494 int32 reloffsetL = relcenML - L;
00495 int32 reloffsetP = relcenMP - P;
00496 int32 offsetL = reloffsetL + initoffsetL;
00497 int32 offsetP = reloffsetP + initoffsetP;
00498
00499 Result(i,0) = offsetL;
00500 Result(i,1) = offsetP;
00501 Result(i,2) = corr;
00502 }
00503
00504
00505 int32 offsetLines = -999;
00506 int32 offsetPixels = -999;
00507 getoffset(Result,offsetLines,offsetPixels);
00508
00509
00510
00511 ofstream scratchlogfile("scratchlogcoarse2", ios::out | ios::trunc);
00512 bk_assert(scratchlogfile,"coarsecorrel: scratchlogcoarse2",__FILE__,__LINE__);
00513 scratchlogfile << "\n\n*******************************************************************"
00514 << "\n* COARSE_COREGISTRATION: Correlation"
00515 << "\n*******************************************************************"
00516 << "\nNumber of correlation windows: \t"
00517 << Nwin
00518 << "\nCorrelation window size (l,p): \t"
00519 << MasksizeL << ", " << MasksizeP;
00520 if (forceoddl) scratchlogfile << "(l forced odd) ";
00521 if (forceoddp) scratchlogfile << "(p forced odd)";
00522 scratchlogfile << "\nSearchwindow size (l,p): \t\t"
00523 << MasksizeL + 2*AccL << ", " << MasksizeP + 2*AccP
00524 << "\nNumber \tposl \tposp \toffsetl offsetp \tcorrelation\n";
00525 register int32 k;
00526 for (k=0; k<Nwin; k++)
00527 scratchlogfile << k << "\t" << Centers(k,0)
00528 << "\t" << Centers(k,1)
00529 << "\t" << Result(k,0)
00530 << "\t" << Result(k,1)
00531 << "\t" << Result(k,2) << endl;
00532 scratchlogfile << "Estimated total offset (l,p): \t"
00533 << offsetLines << ", " << offsetPixels
00534 << "\n*******************************************************************\n";
00535 scratchlogfile.close();
00536
00537 ofstream scratchresfile("scratchrescoarse2", ios::out | ios::trunc);
00538 bk_assert(scratchresfile,"coarsecorrel: scratchrescoarse2",__FILE__,__LINE__);
00539 scratchresfile << "\n\n*******************************************************************"
00540 << "\n*_Start_" << processcontrol[pr_i_coarse2]
00541 << "\n*******************************************************************"
00542 << "\nEstimated translation slave w.r.t. master:"
00543 << "\nCoarse_correlation_translation_lines: \t"
00544 << offsetLines
00545 << "\nCoarse_correlation_translation_pixels: \t"
00546 << offsetPixels
00547 << "\n*******************************************************************"
00548
00549 << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
00550 << "\n*******************************************************************\n";
00551 scratchresfile.close();
00552
00553
00554 INFO << "Individually estimated translations (#, l, p, corr, offl, offp): ";
00555 INFO.print();
00556 for (k=0; k<Nwin; k++)
00557 {
00558
00559 INFO << k << " \t"
00560 << Centers(k,0) << " \t"
00561 << Centers(k,1) << " \t"
00562 << Result(k,2) << " \t"
00563 << Result(k,0) << " \t"
00564 << Result(k,1) << " \t";
00565 INFO.print();
00566 }
00567 INFO << "Estimated translation (l,p): "
00568 << offsetLines << ", " << offsetPixels << ends;
00569 INFO.print();
00570 PROGRESS.print("Coarse coregistration based on correlation finished.");
00571 }
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592 void coarsecorrelfft(
00593 const input_coarsecorr &coarsecorrinput,
00594 const slcimage &minfo,
00595 const slcimage &sinfo)
00596 {
00597 TRACE_FUNCTION("coarsecorrelfft (BK 12-Dec-1998)");
00598 if (coarsecorrinput.method != cc_magfft)
00599 {
00600 PRINT_ERROR("unknown method, This routine is only for cc_magfft method.")
00601 throw(argument_error);
00602 }
00603
00604 char dummyline[ONE27];
00605 const uint Mfilelines = minfo.currentwindow.lines();
00606 const uint Sfilelines = sinfo.currentwindow.lines();
00607
00608 const uint Nwin = coarsecorrinput.Nwin;
00609 const int32 initoffsetL = coarsecorrinput.initoffsetL;
00610 const int32 initoffsetP = coarsecorrinput.initoffsetP;
00611 const uint MasksizeL = coarsecorrinput.MasksizeL;
00612 const uint MasksizeP = coarsecorrinput.MasksizeP;
00613
00614 bool pointsrandom = true;
00615 if (specified(coarsecorrinput.ifpositions))
00616 pointsrandom = false;
00617
00618
00619 if (!ispower2(MasksizeL))
00620 {
00621 PRINT_ERROR("coarse correl fft: MasksizeL should be 2^n")
00622 throw(input_error);
00623 }
00624 if (!ispower2(MasksizeP))
00625 {
00626 PRINT_ERROR("coarse correl fft: MasksizeP should be 2^n")
00627 throw(input_error);
00628 }
00629
00630
00631
00632 const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
00633 const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
00634 const int32 sp0 = sinfo.currentwindow.pixlo - initoffsetP;
00635 const int32 spN = sinfo.currentwindow.pixhi - initoffsetP;
00636
00637
00638
00639
00640
00641
00642
00643
00644 const uint FEW=10;
00645 const uint l0 = max(int32(minfo.currentwindow.linelo),sl0) + FEW;
00646 const uint lN = min(int32(minfo.currentwindow.linehi),slN) - MasksizeL - FEW;
00647 const uint p0 = max(int32(minfo.currentwindow.pixlo),sp0) + FEW;
00648 const uint pN = min(int32(minfo.currentwindow.pixhi),spN) - MasksizeP - FEW;
00649 const window overlap(l0,lN,p0,pN);
00650
00651
00652
00653
00654 register int32 i;
00655
00656 matrix<uint> Minlminp;
00657 if (pointsrandom)
00658 {
00659 Minlminp = distributepoints(real4(Nwin),overlap);
00660 }
00661
00662 else
00663 {
00664 Minlminp.resize(Nwin,3);
00665
00666 ifstream ifpos;
00667 openfstream(ifpos,coarsecorrinput.ifpositions);
00668 bk_assert(ifpos,coarsecorrinput.ifpositions,__FILE__,__LINE__);
00669 uint ll,pp;
00670 for (i=0; i<Nwin; ++i)
00671 {
00672 ifpos >> ll >> pp;
00673 Minlminp(i,0) = uint(ll -.5*MasksizeL);
00674 Minlminp(i,1) = uint(pp -.5*MasksizeP);
00675 Minlminp(i,2) = uint(1);
00676 ifpos.getline(dummyline,ONE27,'\n');
00677 }
00678 ifpos.close();
00679
00680
00681 if (Minlminp(Nwin-1,0) == Minlminp(Nwin-2,0) &&
00682 Minlminp(Nwin-1,1) == Minlminp(Nwin-2,1))
00683 {
00684 Minlminp(Nwin-1,0) = uint(.5*(lN + l0) + 27);
00685 Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);
00686 }
00687
00688
00689
00690 bool troubleoverlap = false;
00691 for (i=0; i<Nwin; ++i)
00692 {
00693 if (Minlminp(i,0) < l0)
00694 {
00695 troubleoverlap=true;
00696 WARNING << "COARSECORR: point from file: "
00697 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
00698 << Minlminp(i,1) +.5*MasksizeP
00699 << " outside overlap master, slave. New position: ";
00700 Minlminp(i,0) = l0 + l0-Minlminp(i,0);
00701 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
00702 WARNING.print();
00703 }
00704 if (Minlminp(i,0) > lN)
00705 {
00706 troubleoverlap=true;
00707 WARNING << "COARSECORR: point from file: "
00708 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
00709 << Minlminp(i,1) +.5*MasksizeP
00710 << " outside overlap master, slave. New position: ";
00711 Minlminp(i,0) = lN + lN-Minlminp(i,0);
00712 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
00713 WARNING.print();
00714 }
00715 if (Minlminp(i,1) < p0)
00716 {
00717 troubleoverlap=true;
00718 WARNING << "COARSECORR: point from file: "
00719 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
00720 << Minlminp(i,1) +.5*MasksizeP
00721 << " outside overlap master, slave. New position: ";
00722 Minlminp(i,1) = p0 + p0-Minlminp(i,1);
00723 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
00724 WARNING.print();
00725 }
00726 if (Minlminp(i,1) > pN)
00727 {
00728 troubleoverlap=true;
00729 WARNING << "COARSECORR: point from file: "
00730 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
00731 << Minlminp(i,1) +.5*MasksizeP
00732 << " outside overlap master, slave. New position: ";
00733 Minlminp(i,1) = pN + pN-Minlminp(i,1);
00734 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
00735 WARNING.print();
00736 }
00737 }
00738 if (troubleoverlap)
00739 {
00740 WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
00741 << l0 << " " << lN << " " << p0 << " " << pN;
00742 WARNING.print();
00743 }
00744 }
00745
00746
00747
00748
00749
00750
00751 matrix<complr4> Master;
00752 matrix<complr4> Mask;
00753 matrix<real4> Result(Nwin,3);
00754
00755
00756 int32 percent = 0;
00757
00758 int32 tenpercent = int32((Nwin/10)+.5);
00759 if (tenpercent==0) tenpercent = 1000;
00760 for (i=0; i<Nwin; ++i)
00761 {
00762 if (i%tenpercent==0)
00763 {
00764 PROGRESS << "COARSE_CORR: " << setw(3) << percent << "%";
00765 PROGRESS.print();
00766 percent += 10;
00767 }
00768
00769
00770 uint minMwinL = Minlminp(i,0);
00771 uint minMwinP = Minlminp(i,1);
00772
00773 window master(minMwinL, minMwinL+MasksizeL-1,
00774 minMwinP, minMwinP+MasksizeP-1);
00775
00776 window mask(minMwinL+initoffsetL,
00777 minMwinL+initoffsetL+MasksizeL-1,
00778 minMwinP+initoffsetP,
00779 minMwinP+initoffsetP+MasksizeP-1);
00780
00781
00782
00783
00784 Master = minfo.readdata(master);
00785 Mask = sinfo.readdata(mask);
00786
00787 matrix<real4> absMaster = magnitude(Master);
00788 matrix<real4> absMask = magnitude(Mask);
00789
00790
00791
00792 real4 offsetL, offsetP;
00793 real4 coheren;
00794
00795 coheren = corrfft(absMaster,absMask,offsetL,offsetP);
00796
00797 offsetL += initoffsetL;
00798 offsetP += initoffsetP;
00799
00800 Result(i,0) = offsetL;
00801 Result(i,1) = offsetP;
00802 Result(i,2) = coheren;
00803 }
00804
00805
00806
00807 for (i=0;i<Nwin;i++)
00808 {
00809 Minlminp(i,0) += uint(.5 * MasksizeL);
00810 Minlminp(i,1) += uint(.5 * MasksizeP);
00811 }
00812
00813
00814 int32 offsetLines = -999;
00815 int32 offsetPixels = -999;
00816 getoffset(Result,offsetLines,offsetPixels);
00817
00818
00819
00820 ofstream scratchlogfile("scratchlogcoarse2", ios::out | ios::trunc);
00821 bk_assert(scratchlogfile,"coarsecorrelfft: scratchlogcoarse2",__FILE__,__LINE__);
00822 scratchlogfile << "\n\n*******************************************************************"
00823 << "\n* COARSE_COREGISTRATION: Correlation"
00824 << "\n*******************************************************************"
00825 << "\nNumber of correlation windows: \t"
00826 << Nwin
00827 << "\nwindow size (l,p): \t"
00828 << MasksizeL << ", " << MasksizeP
00829
00830
00831 << "\n\nNumber \tposL \tposP \toffsetL offsetP\tcorrelation\n";
00832 register int32 k;
00833 for (k=0; k<Nwin; k++)
00834 scratchlogfile << k << "\t" << Minlminp(k,0)
00835 << "\t" << Minlminp(k,1)
00836 << "\t" << Result(k,0)
00837 << "\t" << Result(k,1)
00838 << "\t" << Result(k,2) << endl;
00839 scratchlogfile << "Estimated total offset (l,p): \t"
00840 << offsetLines << ", " << offsetPixels
00841 << "\n*******************************************************************\n";
00842 scratchlogfile.close();
00843
00844 ofstream scratchresfile("scratchrescoarse2", ios::out | ios::trunc);
00845 bk_assert(scratchresfile,"coarsecorrelfft: scratchrescoarse2",__FILE__,__LINE__);
00846 scratchresfile << "\n\n*******************************************************************"
00847 << "\n*_Start_" << processcontrol[pr_i_coarse2]
00848 << "\n*******************************************************************"
00849 << "\nEstimated translation slave w.r.t. master:"
00850 << "\nCoarse_correlation_translation_lines: \t"
00851 << offsetLines
00852 << "\nCoarse_correlation_translation_pixels: \t"
00853 << offsetPixels
00854 << "\n*******************************************************************"
00855
00856 << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
00857 << "\n*******************************************************************\n";
00858
00859
00860 scratchresfile.close();
00861 INFO << "Individual estimated translations (#, l, p, corr, offl, offp):";
00862 INFO.print();
00863 for (k=0; k<Nwin; k++)
00864 {
00865 INFO
00866 << k << " \t"
00867 << Minlminp(k,0) << " \t"
00868 << Minlminp(k,1) << " \t"
00869 << Result(k,2) << " \t"
00870 << Result(k,0) << " \t"
00871 << Result(k,1) << " \t";
00872 INFO.print();
00873 }
00874
00875 INFO << "Estimated translation (l,p): "
00876 << offsetLines << ", " << offsetPixels;
00877 INFO.print();
00878 PROGRESS.print("Coarse coregistration based on correlation finished.");
00879 }
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 real4 corrfft(
00907 const matrix<real4> &Master,
00908 const matrix<real4> &Mask,
00909 real4 &offsetL,
00910 real4 &offsetP)
00911 {
00912 TRACE_FUNCTION("corrfft (BK 18-Oct-1999)");
00913 const int32 L = Master.lines();
00914 const int32 P = Master.pixels();
00915 const int32 twoL = 2*L;
00916 const int32 twoP = 2*P;
00917 const int32 halfL = L/2;
00918 const int32 halfP = P/2;
00919
00920 if (L != Mask.lines() || P != Mask.pixels())
00921 {
00922 PRINT_ERROR("Mask, Master not same size.")
00923 throw(input_error);
00924 }
00925 if (!(ispower2(L) || ispower2(P)))
00926 {
00927 PRINT_ERROR("Mask, Master size not power of 2.")
00928 throw(input_error);
00929 }
00930
00931
00932 register int32 i;
00933 register int32 j;
00934
00935 const complr4 ONE(1.);
00936 matrix<complr4> Master2(twoL,twoP);
00937 matrix<complr4> Mask2(twoL,twoP);
00938 matrix<complr4> blok2(twoL,twoP);
00939
00940 const real4 meanMaster = mean(Master);
00941 const real4 meanMask = mean(Mask);
00942
00943
00944
00945 for (i=0; i<L; i++)
00946 {
00947 for (j=0; j<P; j++)
00948 {
00949 Master2(i,j) = complr4(sqr(Master(i,j)-meanMaster));
00950 Mask2(i+L,j+P) = complr4(sqr(Mask(i,j)-meanMask));
00951 blok2(i+halfL,j+halfP) = ONE;
00952 }
00953 }
00954
00955 fft2d(Master2);
00956 fft2d(Mask2);
00957 fft2d(blok2);
00958
00959
00960
00961
00962
00963
00964
00965
00966 blok2.conj();
00967 Master2 *= blok2;
00968 Master2.conj();
00969 ifft2d(Master2);
00970
00971 Mask2 *= blok2;
00972 ifft2d(Mask2);
00973
00974
00975 for (i=0; i<=L; i++)
00976 for (j=0; j<=P; j++)
00977 blok2(i,j) = complr4(sqrt(real(Master2(i,j))*real(Mask2(i,j))));
00978
00979
00980
00981 Master2.clean();
00982 Mask2.clean();
00983 for (i=0;i<L;i++)
00984 {
00985 for (j=0;j<P;j++)
00986 {
00987 Master2(i,j) = complr4(Master(i,j)-meanMaster);
00988 Mask2(i+halfL,j+halfP) = complr4(Mask(i,j)-meanMask);
00989 }
00990 }
00991
00992
00993
00994
00995 fft2d(Master2);
00996 fft2d(Mask2);
00997
00998
00999 Master2.conj();
01000 Mask2 *= Master2;
01001 ifft2d(Mask2);
01002
01003
01004
01005
01006
01007 real4 coher;
01008 real4 maxcoher = 0.;
01009 for (i=0; i<=L; i++)
01010 {
01011 for (j=0; j<=P; j++)
01012 {
01013
01014
01015
01016
01017 coher = real(Mask2(i,j)) / real(blok2(i,j));
01018 if (coher > maxcoher)
01019 {
01020 maxcoher = coher;
01021 offsetL = -halfL + i;
01022 offsetP = -halfP + j;
01023 }
01024 }
01025 }
01026 return maxcoher;
01027 }
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047 matrix<uint> distributepoints(
01048 real4 nW,
01049 const window &win)
01050 {
01051 TRACE_FUNCTION("distributepoints (BK 21-Jan-1999)")
01052 real4 lines = win.linehi - win.linelo + 1;
01053 real4 pixels = win.pixhi - win.pixlo + 1;
01054
01055 uint numw = uint(nW);
01056 matrix<uint> Result(numw,uint(3));
01057
01058 real4 wp = sqrt(nW/(lines/pixels));
01059 real4 wl = nW / wp;
01060
01061 if (wl < wp)
01062 wl = wp;
01063
01064 int32 wlint = int32(floor(wl + .5));
01065
01066 real4 deltal = (lines-1) / (real4(wlint-1));
01067 int32 totp = int32(pixels*wlint);
01068 real4 deltap = (real4(totp-1)) / (real4(nW-1));
01069
01070 real4 p = -deltap;
01071 real4 l = 0.;
01072 uint lcnt = 0;
01073
01074 register int32 i;
01075 for (i=0; i<nW; i++)
01076 {
01077 p += deltap;
01078 while (floor(p+.5)>=pixels)
01079 {
01080 p -= pixels;
01081 lcnt++;
01082 }
01083 l = lcnt * deltal;
01084 Result(i,0) = uint(floor(l+.5));
01085 Result(i,1) = uint(floor(p+.5));
01086 }
01087
01088
01089 for (i=0; i<nW; i++)
01090 {
01091 Result(i,0) += win.linelo;
01092 Result(i,1) += win.pixlo;
01093 }
01094
01095 return Result;
01096 }
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118 void getoffset(
01119 const matrix<real4> &Result,
01120 int32 &offsetLines,
01121 int32 &offsetPixels)
01122 {
01123 TRACE_FUNCTION("getoffset (BK 21-Jan-1999)")
01124 if (Result.pixels() != 3)
01125 {
01126 PRINT_ERROR("code 901: wrong input not 3 width");
01127 throw(input_error);
01128 }
01129 uint nW = Result.lines();
01130 int32 cnt;
01131 int32 valueL;
01132 int32 valueP;
01133 real4 correl;
01134 int32 highestcnt = 0;
01135 real4 highestcorrel = 0.;
01136 for (register int32 i=0;i<nW;i++)
01137 {
01138 valueL = int32(Result(i,0));
01139 valueP = int32(Result(i,1));
01140 correl = Result(i,2);
01141 if (correl > highestcorrel)
01142 highestcorrel = correl;
01143 cnt = 0;
01144 for (register int32 j=0;j<nW;j++)
01145 {
01146 if (abs(Result(j,0) - valueL) < 2 &&
01147 abs(Result(j,1) - valueP) < 2)
01148 cnt++;
01149 }
01150 if (cnt > highestcnt)
01151 {
01152 highestcnt = cnt;
01153 offsetLines = valueL;
01154 offsetPixels = valueP;
01155 }
01156 }
01157
01158
01159 real4 THRESHOLD = .4;
01160 if (nW < 6)
01161 {
01162 WARNING.print("getoffset: number of windows to estimate offset < 6");
01163 WARNING.print("(please check bottom of LOGFILE)");
01164 }
01165 if (highestcnt < .2*nW)
01166 {
01167 WARNING.print("getoffset: estimated offset not consistent with other estimates.");
01168 WARNING.print("(check bottom of LOGFILE)");
01169 }
01170 if (highestcorrel < THRESHOLD)
01171 {
01172 WARNING << "getoffset: estimated translation has correlation of: "
01173 << highestcorrel;
01174 WARNING.print();
01175 WARNING.print("(please check bottom of LOGFILE)");
01176 }
01177 }
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200 void finecoreg(
01201 const input_fine &fineinput,
01202 const slcimage &minfo,
01203 const slcimage &sinfo)
01204 {
01205 TRACE_FUNCTION("finecoreg (BK 29-Oct-99)")
01206 char dummyline[ONE27];
01207
01208 const uint Mfilelines = minfo.currentwindow.lines();
01209 const uint Sfilelines = sinfo.currentwindow.lines();
01210
01211 const uint Nwin = fineinput.Nwin;
01212 const int32 initoffsetL = fineinput.initoffsetL;
01213 const int32 initoffsetP = fineinput.initoffsetP;
01214 uint MasksizeL = fineinput.MasksizeL;
01215 uint MasksizeP = fineinput.MasksizeP;
01216
01217 bool pointsrandom = true;
01218 if (specified(fineinput.ifpositions))
01219 pointsrandom = false;
01220
01221
01222 if (fineinput.method == fc_magspace ||
01223 fineinput.method == fc_cmplxspace)
01224 {
01225 MasksizeL += 2*fineinput.AccL;
01226 MasksizeP += 2*fineinput.AccP;
01227 }
01228
01229
01230
01231
01232 const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
01233 const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
01234 const int32 sp0 = sinfo.currentwindow.pixlo - initoffsetP;
01235 const int32 spN = sinfo.currentwindow.pixhi - initoffsetP;
01236
01237
01238
01239
01240
01241
01242
01243
01244 const uint FEW=10;
01245 const uint l0 = max(int32(minfo.currentwindow.linelo),sl0) + FEW;
01246 const uint lN = min(int32(minfo.currentwindow.linehi),slN) - MasksizeL - FEW;
01247 const uint p0 = max(int32(minfo.currentwindow.pixlo),sp0) + FEW;
01248 const uint pN = min(int32(minfo.currentwindow.pixhi),spN) - MasksizeP - FEW;
01249 const window overlap(l0,lN,p0,pN);
01250
01251
01252
01253 register int32 i;
01254 matrix<uint> Minlminp;
01255 if (pointsrandom)
01256 {
01257 Minlminp = distributepoints(real4(Nwin),overlap);
01258 }
01259
01260 else
01261 {
01262
01263 Minlminp.resize(Nwin,3);
01264
01265 ifstream ifpos(fineinput.ifpositions, ios::in);
01266 bk_assert(ifpos,fineinput.ifpositions,__FILE__,__LINE__);
01267 uint ll,pp;
01268 for (i=0; i<Nwin; ++i)
01269 {
01270 ifpos >> ll >> pp;
01271 Minlminp(i,0) = uint(ll -.5*MasksizeL);
01272 Minlminp(i,1) = uint(pp -.5*MasksizeP);
01273 Minlminp(i,2) = uint(1);
01274 ifpos.getline(dummyline,ONE27,'\n');
01275 }
01276 ifpos.close();
01277
01278
01279 if (Minlminp(Nwin-1,0) == Minlminp(Nwin-2,0) &&
01280 Minlminp(Nwin-1,1) == Minlminp(Nwin-2,1))
01281 {
01282 Minlminp(Nwin-1,0) = uint(.5*(lN + l0) + 27);
01283 Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);
01284 }
01285
01286
01287
01288 bool troubleoverlap = false;
01289 for (i=0; i<Nwin; ++i)
01290 {
01291 if (Minlminp(i,0) < l0)
01292 {
01293 troubleoverlap=true;
01294 WARNING << "FINE: point from file: "
01295 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
01296 << Minlminp(i,1) +.5*MasksizeP
01297 << " outside overlap master, slave. New position: ";
01298 Minlminp(i,0) = l0 + l0-Minlminp(i,0);
01299 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
01300 WARNING.print();
01301 }
01302 if (Minlminp(i,0) > lN)
01303 {
01304 troubleoverlap=true;
01305 WARNING << "FINE: point from file: "
01306 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
01307 << Minlminp(i,1) +.5*MasksizeP
01308 << " outside overlap master, slave. New position: ";
01309 Minlminp(i,0) = lN + lN-Minlminp(i,0);
01310 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
01311 WARNING.print();
01312 }
01313 if (Minlminp(i,1) < p0)
01314 {
01315 troubleoverlap=true;
01316 WARNING << "FINE: point from file: "
01317 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
01318 << Minlminp(i,1) +.5*MasksizeP
01319 << " outside overlap master, slave. New position: ";
01320 Minlminp(i,1) = p0 + p0-Minlminp(i,1);
01321 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
01322 WARNING.print();
01323 }
01324 if (Minlminp(i,1) > pN)
01325 {
01326 troubleoverlap=true;
01327 WARNING << "FINE: point from file: "
01328 << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
01329 << Minlminp(i,1) +.5*MasksizeP
01330 << " outside overlap master, slave. New position: ";
01331 Minlminp(i,1) = pN + pN-Minlminp(i,1);
01332 WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
01333 WARNING.print();
01334 }
01335 }
01336 if (troubleoverlap)
01337 {
01338 WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
01339 << l0 << " " << lN << " " << p0 << " " << pN;
01340 WARNING.print();
01341 }
01342 }
01343
01344
01345
01346
01347
01348 matrix<complr4> Master;
01349 matrix<complr4> Mask;
01350 matrix<real4> Result(Nwin,3);
01351
01352
01353
01354 int32 tenpercent = int32((Nwin/10)+.5);
01355 if (tenpercent==0) tenpercent = 1000;
01356 int32 percent = 0;
01357
01358
01359 for (i=0;i<Nwin;i++)
01360 {
01361
01362 if (i%tenpercent==0)
01363 {
01364 PROGRESS << "FINE: " << setw(3) << percent << "%";
01365 PROGRESS.print();
01366 percent += 10;
01367 }
01368
01369
01370 uint minMwinL = Minlminp(i,0);
01371 uint minMwinP = Minlminp(i,1);
01372
01373
01374 window master(minMwinL, minMwinL+MasksizeL-1,
01375 minMwinP, minMwinP+MasksizeP-1);
01376
01377 window mask(minMwinL+initoffsetL,
01378 minMwinL+initoffsetL+MasksizeL-1,
01379 minMwinP+initoffsetP,
01380 minMwinP+initoffsetP+MasksizeP-1);
01381
01382
01383
01384
01385 Master = minfo.readdata(master);
01386 Mask = sinfo.readdata(mask);
01387
01388
01389
01390
01391 real4 offsetL, offsetP;
01392 real4 coheren;
01393 switch (fineinput.method)
01394 {
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 case fc_magfft:
01406 coheren = coherencefft(fineinput, Master, Mask, offsetL, offsetP);
01407 break;
01408 case fc_magspace:
01409 coheren = coherencespace(fineinput, Master, Mask, offsetL, offsetP);
01410 break;
01411 default:
01412 PRINT_ERROR("unknown method for fine coregistration.")
01413 throw(unhandled_case_error);
01414 }
01415
01416
01417 Result(i,0) = offsetL + initoffsetL;
01418 Result(i,1) = offsetP + initoffsetP;
01419 Result(i,2) = coheren;
01420 }
01421
01422
01423
01424 for (i=0; i<Nwin; i++)
01425 {
01426 Minlminp(i,0) += uint(.5 * MasksizeL);
01427 Minlminp(i,1) += uint(.5 * MasksizeP);
01428 }
01429
01430
01431
01432 ofstream scratchlogfile("scratchlogfine", ios::out | ios::trunc);
01433 bk_assert(scratchlogfile,"finecoreg: scratchlogfine",__FILE__,__LINE__);
01434 scratchlogfile << "\n\n*******************************************************************"
01435 << "\n* FINE_COREGISTRATION"
01436 << "\n*******************************************************************"
01437 << "\nNumber of correlation windows: \t"
01438 << Nwin
01439 << "\nwindow size (l,p): \t"
01440 << MasksizeL << ", " << MasksizeP
01441 << "\nInitial offsets: \t"
01442 << initoffsetL << ", " << initoffsetP
01443 << "\nOversampling factor: \t"
01444 << fineinput.osfactor
01445 << "\n\nNumber \tposl \tposp \toffsetl offsetp\tcorrelation\n";
01446 for (i=0;i<Nwin;i++)
01447 scratchlogfile
01448 << setiosflags(ios::fixed)
01449 << setiosflags(ios::showpoint)
01450 << setiosflags(ios::right)
01451 << setw(8)
01452 << setprecision(0)
01453 << i << " "
01454 << Minlminp(i,0) << " "
01455 << Minlminp(i,1) << " "
01456 << setprecision(3)
01457 << Result(i,0) << " "
01458 << Result(i,1) << " "
01459 << setprecision(2)
01460 << Result(i,2) << endl;
01461 scratchlogfile << "\n*******************************************************************\n";
01462 scratchlogfile.close();
01463
01464 ofstream scratchresfile("scratchresfine", ios::out | ios::trunc);
01465 bk_assert(scratchresfile,"finecoreg: scratchresfine",__FILE__,__LINE__);
01466
01467 scratchresfile
01468 << "\n\n*******************************************************************"
01469 << "\n*_Start_" << processcontrol[pr_i_fine]
01470 << "\n*******************************************************************"
01471 << "\nOversampling factor: \t\t"
01472 << fineinput.osfactor
01473 << "\nNumber_of_correlation_windows: \t"
01474 << Nwin
01475 << "\nNumber \tposL \tposP \toffsetL offsetP\tcorrelation\n";
01476 scratchresfile.close();
01477
01478 FILE *resfile;
01479 resfile=fopen("scratchresfine","a");
01480 for (i=0; i<Nwin; i++)
01481 fprintf(resfile,"%4.0f %5.0f %5.0f %# 9.2f %# 9.2f %# 6.2f\n",
01482 real4(i), real4(Minlminp(i,0)), real4(Minlminp(i,1)),
01483 Result(i,0), Result(i,1), Result(i,2));
01484 fprintf(resfile,
01485 "\n*******************************************************************");
01486 fprintf(resfile,"%s%s%s",
01487 "\n* End_", processcontrol[pr_i_fine], "_NORMAL");
01488 fprintf(resfile,
01489 "\n*******************************************************************\n");
01490
01491
01492
01493 fclose(resfile);
01494 PROGRESS.print("Fine coregistration finished.");
01495 }
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517 real4 coherencefft(
01518 const input_fine &fineinput,
01519 const matrix<complr4> &Master,
01520 const matrix<complr4> &Mask,
01521 real4 &offsetL,
01522 real4 &offsetP)
01523 {
01524 TRACE_FUNCTION("coherencefft (BK 16-Nov-1999)")
01525 const uint factor = fineinput.osfactor;
01526 const uint AccL = fineinput.AccL;
01527 const uint AccP = fineinput.AccP;
01528 const int32 L = Master.lines();
01529 const int32 P = Master.pixels();
01530 const int32 twoL = 2*L;
01531 const int32 twoP = 2*P;
01532 const int32 halfL = L/2;
01533 const int32 halfP = P/2;
01534
01535 if (fineinput.method != fc_magfft)
01536 {
01537 PRINT_ERROR("impossible error, wrong input method in fft.")
01538 throw(input_error);
01539 }
01540 if (!ispower2(factor))
01541 {
01542 PRINT_ERROR("coherence factor not power of 2")
01543 throw(input_error);
01544 }
01545 if (L != Mask.lines() || P != Mask.pixels())
01546 {
01547 PRINT_ERROR("Mask, Master not same size.")
01548 throw(input_error);
01549 }
01550 if (!(ispower2(L) || ispower2(P)))
01551 {
01552 PRINT_ERROR("Mask, Master size not power of 2.")
01553 throw(input_error);
01554 }
01555
01556
01557 matrix<real4> magMaster = magnitude(Master);
01558 matrix<real4> magMask = magnitude(Mask);
01559 magMaster -= mean(magMaster);
01560 magMask -= mean(magMask);
01561
01562
01563
01564
01565 window windef(0,0,0,0);
01566 window win1(0,L-1,0,P-1);
01567 window win2(halfL,halfL+L-1,halfP,halfP+P-1);
01568
01569 matrix<complr4> Master2(twoL,twoP);
01570 matrix<complr4> Mask2(twoL,twoP);
01571 Master2.setdata(win1,mat2cr4(magMaster),windef);
01572 Mask2.setdata(win2,mat2cr4(magMask),windef);
01573
01574 fft2d(Master2);
01575 fft2d(Mask2);
01576
01577
01578
01579 Master2.conj();
01580 Mask2 *= Master2;
01581 ifft2d(Mask2);
01582
01583
01584 window wintmp(halfL-AccL,halfL+AccL-1,halfP-AccP,halfP+AccP-1);
01585 matrix<complr4> TMP(wintmp,Mask2);
01586 matrix<real4> Covar = real(TMP);
01587
01588
01589
01590 matrix<complr4> blok(L,P);
01591 const complr4 ONE(1.);
01592 blok.setdata(ONE);
01593
01594 Master2.clean();
01595 Mask2.clean();
01596 Master2.setdata(win1,mat2cr4(sqr(magMaster)),windef);
01597 Mask2.setdata(win2,blok,windef);
01598
01599 fft2d(Master2);
01600 fft2d(Mask2);
01601
01602 Mask2.conj();
01603 Master2 *= Mask2;
01604 Master2.conj();
01605 ifft2d(Master2);
01606
01607
01608 TMP.setdata(Master2,wintmp);
01609 matrix<real4> pmaster = real(TMP);
01610
01611
01612 Master2.clean();
01613 window win5(L,twoL-1,P,twoP-1);
01614 Master2.setdata(win5,mat2cr4(sqr(magMask)),windef);
01615 fft2d(Master2);
01616 Master2 *= Mask2;
01617 ifft2d(Master2);
01618
01619
01620 TMP.setdata(Master2,wintmp);
01621 matrix<real4> pmask = real(TMP);
01622
01623
01624 pmaster *= pmask;
01625 Covar /= sqrt(pmaster);
01626
01627
01628
01629 uint offL;
01630 uint offP;
01631
01632 matrix<real4> coher8 = oversample(Covar,factor,factor);
01633 real4 maxcor = max(coher8,offL,offP);
01634 offsetL = -real4(AccL) + offL/real4(factor);
01635 offsetP = -real4(AccP) + offP/real4(factor);
01636
01637 return maxcor;
01638 }
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657 real4 coherencespace(
01658 const input_fine &fineinput,
01659 const matrix<complr4> &Master,
01660 const matrix<complr4> &Mask,
01661 real4 &offsetL,
01662 real4 &offsetP)
01663 {
01664 TRACE_FUNCTION("coherencespace (BK 03-Feb-1999)")
01665 const int32 L = Master.lines();
01666 const int32 P = Master.pixels();
01667
01668 const int32 AccL = fineinput.AccL;
01669 const int32 AccP = fineinput.AccP;
01670 const uint factor = fineinput.osfactor;
01671
01672
01673 const int32 MasksizeL = L - 2*AccL;
01674 const int32 MasksizeP = P - 2*AccP;
01675
01676 if (!ispower2(AccL) || !ispower2(AccP))
01677 {
01678 PRINT_ERROR("AccL should be power of 2 for oversampling.")
01679 throw(input_error);
01680 }
01681 if (MasksizeL < 4 || MasksizeP < 4)
01682 {
01683 PRINT_ERROR("Correlationwindow size too small (<4; size= FC_winsize-2*FC_Acc).")
01684 throw(input_error);
01685 }
01686
01687
01688 window winmask(AccL, AccL+MasksizeL-1,
01689 AccP, AccP+MasksizeP-1);
01690
01691 matrix<real4> coher(2*AccL,2*AccP);
01692
01693 window windef(0, 0, 0, 0);
01694
01695 switch (fineinput.method)
01696 {
01697 case fc_cmplxspace:
01698 {
01699 PRINT_ERROR("not implemented in v1.0")
01700 throw(unhandled_case_error);
01701 break;
01702 }
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740 case fc_magspace:
01741 {
01742 matrix<real4> magMask = magnitude(Mask);
01743 magMask -= mean(magMask);
01744
01745 matrix<real4> Mask2(winmask,magMask);
01746 real4 normmask = norm2(Mask2);
01747
01748
01749 matrix<real4> Master2(MasksizeL, MasksizeP);
01750 matrix<real4> magMaster=magnitude(Master);
01751 magMaster -= mean(magMaster);
01752 window winmaster;
01753
01754
01755 for (register int32 i=0;i<2*AccL;i++)
01756 {
01757 winmaster.linelo = i;
01758 winmaster.linehi = i+MasksizeL-1;
01759 for (register int32 j=0;j<2*AccP;j++)
01760 {
01761 winmaster.pixlo = j;
01762 winmaster.pixhi = j+MasksizeP-1;
01763 Master2.setdata(windef,magMaster,winmaster);
01764
01765
01766 real4 cohs1s2 = 0.;
01767 real4 cohs1s1 = 0.;
01768 for (register int32 k=0;k<MasksizeL;k++)
01769 {
01770 for (register int32 l=0;l<MasksizeP;l++)
01771 {
01772 cohs1s2 += (Master2(k,l) * Mask2(k,l));
01773 cohs1s1 += sqr(Master2(k,l));
01774 }
01775 }
01776 coher(i,j) = cohs1s2 / sqrt(cohs1s1 * normmask);
01777
01778 }
01779 }
01780 break;
01781 }
01782
01783 default:
01784 PRINT_ERROR("unknown method")
01785 throw(unhandled_case_error);
01786 }
01787
01788
01789
01790
01791 matrix<real4> coher8 = oversample(coher,factor,factor);
01792
01793 uint offL;
01794 uint offP;
01795 real4 maxcor = max(coher8,offL,offP);
01796 offsetL = AccL - offL/real4(factor);
01797 offsetP = AccP - offP/real4(factor);
01798
01799 return maxcor;
01800 }
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822 void coregpm(
01823 const window &originalmaster,
01824 const char* i_resfile,
01825 const input_coregpm &coregpminput,
01826 uint oversamplingsfactorfine)
01827 {
01828 TRACE_FUNCTION("coregpm (BK 26-Oct-1999)")
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839 const real4 THRESHOLD = coregpminput.threshold;
01840 const int32 DEGREE = coregpminput.degree;
01841 const int32 MAX_ITERATIONS = coregpminput.maxiter;
01842 const real4 CRIT_VALUE = coregpminput.k_alpha;
01843 const int32 Nunk = Ncoeffs(DEGREE);
01844
01845
01846 const real8 minL = originalmaster.linelo;
01847 const real8 maxL = originalmaster.linehi;
01848 const real8 minP = originalmaster.pixlo;
01849 const real8 maxP = originalmaster.pixhi;
01850
01851
01852
01853
01854
01855 const real4 ACCURACY = .5 * (1. / (real4(oversamplingsfactorfine)));
01856
01857
01858
01859
01860
01861 const real4 SIGMAL = 0.15;
01862 const real4 SIGMAP = SIGMAL/1.5;
01863
01864
01865
01866
01867 DEBUG.print("using a better sigma in range, since seems that can be estimated better");
01868 INFO << "a priori std.dev offset vectors line direction [samples]: " << SIGMAL;
01869 INFO.print();
01870 INFO << "a priori std.dev offset vectors pixel direction [samples]: " << SIGMAP;
01871 INFO.print();
01872
01873
01874 matrix<real4> Data = getofffile(i_resfile, THRESHOLD);
01875
01876
01877
01878
01879
01880 int32 ITERATION = 0;
01881 int32 DONE = 0;
01882
01883
01884 INFO << "Critical value for outlier test: " << CRIT_VALUE;
01885 INFO.print();
01886 uint winL = 0;
01887 uint winP = 0;
01888 matrix<real8> eL_hat;
01889 matrix<real8> eP_hat;
01890 matrix<real8> wtestL;
01891 matrix<real8> wtestP;
01892 matrix<real8> rhsL;
01893 matrix<real8> rhsP;
01894 matrix<real8> Qx_hat;
01895 real8 maxdev = 0.0;
01896 real8 overallmodeltestL = 0.0;
01897 real8 overallmodeltestP = 0.0;
01898 real8 maxwL;
01899 real8 maxwP;
01900 register int32 i,j,k,index;
01901 while (DONE != 1)
01902 {
01903 PROGRESS << "Start iteration " << ITERATION;
01904 PROGRESS.print();
01905
01906 if (ITERATION != 0)
01907 {
01908 matrix<real4> tmp_DATA = Data;
01909 Data.resize(Data.lines()-1, Data.pixels());
01910 j = 0;
01911 for (i=0; i<tmp_DATA.lines(); i++)
01912 {
01913 if (i != winL)
01914 {
01915 Data.setrow(j,tmp_DATA.getrow(i));
01916 j++;
01917 }
01918 else
01919 {
01920 INFO << "Removing observation " << i << " from observation vector.";
01921 INFO.print();
01922 }
01923 }
01924 }
01925
01926
01927 int32 Nobs = Data.lines();
01928 if (Nobs < Nunk)
01929 {
01930 PRINT_ERROR("coregpm: Number of windows > threshold is smaller than parameters solved for.")
01931 throw(input_error);
01932 }
01933
01934
01935
01936 matrix<real8> yL(Nobs,1);
01937 matrix<real8> yP(Nobs,1);
01938 matrix<real8> A(Nobs,Nunk);
01939 matrix<real8> Qy_1(Nobs,1);
01940
01941
01942 INFO << "coregpm: polynomial normalized by factors: "
01943 << minL << " " << maxL << " " << minP << " " << maxP << " to [-2,2]";
01944 INFO.print();
01945
01946
01947 DEBUG.print("Setting up design matrix for LS adjustment");
01948 for (i=0; i<Nobs; i++)
01949 {
01950 real8 posL = normalize(real8(Data(i,1)),minL,maxL);
01951 real8 posP = normalize(real8(Data(i,2)),minP,maxP);
01952 yL(i,0) = real8(Data(i,3));
01953 yP(i,0) = real8(Data(i,4));
01954 DEBUG << "coregpm: (" << posL << ", "<< posP << "): yL="
01955 << yL(i,0) << " yP=" << yP(i,0);
01956 DEBUG.print();
01957
01958 index = 0;
01959 for (j=0; j<=DEGREE; j++)
01960 {
01961 for (k=0; k<=j; k++)
01962 {
01963 A(i,index) = pow(posL,real8(j-k)) * pow(posP,real8(k));
01964 index++;
01965 }
01966 }
01967 }
01968
01969
01970
01971 DEBUG.print("Setting up covariance matrix for LS adjustment");
01972 switch(coregpminput.weightflag)
01973 {
01974 case 0:
01975
01976 for (i=0; i<Nobs; i++)
01977 Qy_1(i,0) = real8(1.0);
01978 break;
01979 case 1:
01980
01981 DEBUG.print("Using sqrt(coherence) as weights.");
01982 for (i=0; i<Nobs; i++)
01983 Qy_1(i,0) = real8(Data(i,5));
01984 break;
01985 case 2:
01986 DEBUG.print("Using coherence as weights.");
01987 for (i=0; i<Nobs; i++)
01988 Qy_1(i,0) = real8(Data(i,5))*real8(Data(i,5));
01989 break;
01990 default:
01991 PRINT_ERROR("Panic, NOT possible with checked input.")
01992 throw(unhandled_case_error);
01993 }
01994
01995 INFO.print("Normalizing covariance matrix for LS estimation.");
01996 Qy_1 = Qy_1 / mean(Qy_1);
01997
01998
01999
02000 matrix<real8> N = matTxmat(A,diagxmat(Qy_1,A));
02001
02002
02003
02004 rhsL = matTxmat(A,diagxmat(Qy_1,yL));
02005 rhsP = matTxmat(A,diagxmat(Qy_1,yP));
02006 Qx_hat = N;
02007
02008 choles(Qx_hat);
02009 solvechol(Qx_hat,rhsL);
02010 solvechol(Qx_hat,rhsP);
02011 invertchol(Qx_hat);
02012
02013 for (i=0; i<Qx_hat.lines(); i++)
02014 for (j=0; j<i; j++)
02015 Qx_hat(j,i) = Qx_hat(i,j);
02016 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
02017 INFO << "coregpm: max(abs(N*inv(N)-I)) = " << maxdev;
02018 INFO.print();
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036 if (maxdev > .01)
02037 {
02038 ERROR << "coregpm: maximum deviation N*inv(N) from unity = " << maxdev
02039 << ". This is larger than .01";
02040 ERROR.print(ERROR.get_str());
02041 throw(some_error);
02042 }
02043 else if (maxdev > .001)
02044 {
02045 WARNING << "coregpm: maximum deviation N*inv(N) from unity = " << maxdev
02046 << ". This is between 0.01 and 0.001";
02047 WARNING.print();
02048 }
02049
02050
02051
02052 matrix<real8> Qy_hat = A * (matxmatT(Qx_hat,A));
02053 matrix<real8> yL_hat = A * rhsL;
02054 matrix<real8> yP_hat = A * rhsP;
02055
02056
02057 eL_hat = yL - yL_hat;
02058 eP_hat = yP - yP_hat;
02059
02060 matrix<real8> Qe_hat = -Qy_hat;
02061 for (i=0; i<Nobs; i++)
02062 Qe_hat(i,i) += (1. / Qy_1(i,0));
02063
02064
02065
02066 overallmodeltestL = 0.;
02067 overallmodeltestP = 0.;
02068 for (i=0; i<Nobs; i++)
02069 {
02070 overallmodeltestL += sqr(eL_hat(i,0))*Qy_1(i,0);
02071 overallmodeltestP += sqr(eP_hat(i,0))*Qy_1(i,0);
02072 }
02073 overallmodeltestL = (overallmodeltestL/sqr(SIGMAL)) /(Nobs-Nunk);
02074 overallmodeltestP = (overallmodeltestP/sqr(SIGMAP)) /(Nobs-Nunk);
02075
02076
02077 INFO << "coregpm: overallmodeltest Lines = " << overallmodeltestL;
02078 INFO.print();
02079 INFO << "coregpm: overallmodeltest Pixels = " << overallmodeltestP;
02080 INFO.print();
02081
02082
02083
02084
02085
02086 wtestL.resize(Nobs,1);
02087 wtestP.resize(Nobs,1);
02088
02089 for (i=0; i<Nobs; i++)
02090 {
02091 wtestL(i,0) = eL_hat(i,0) / (sqrt(Qe_hat(i,i))*SIGMAL);
02092 wtestP(i,0) = eP_hat(i,0) / (sqrt(Qe_hat(i,i))*SIGMAP);
02093 }
02094
02095 uint dumm = 0;
02096 maxwL = max(abs(wtestL),winL,dumm);
02097 maxwP = max(abs(wtestP),winP,dumm);
02098 INFO << "maximum wtest statistic azimuth = " << maxwL
02099 << " for window number: "
02100 << Data(winL,0);
02101 INFO.print();
02102 INFO << "maximum wtest statistic range = " << maxwP
02103 << " for window number: "
02104 << Data(winP,0);
02105 INFO.print();
02106
02107
02108
02109
02110
02111
02112
02113 matrix<real8> wtestsum = sqr(wtestL)+sqr(wtestP);
02114 real8 maxwsum = max(wtestsum,winL,dumm);
02115 INFO << "Detected outlier: summed sqr.wtest = " << maxwsum
02116 << "; observation: " << winL
02117 << "; window number: "
02118 << Data(winL,0);
02119 INFO.print();
02120
02121
02122
02123 if (Nobs <= Nunk)
02124 {
02125 WARNING.print("NO redundancy! Exiting iterations.");
02126 DONE = 1;
02127 }
02128
02129
02130
02131
02132
02133
02134 if (max(maxwL,maxwP) <= CRIT_VALUE)
02135 {
02136 INFO.print("All outlier tests accepted! (final solution computed)");
02137 DONE = 1;
02138 }
02139 if (ITERATION >= MAX_ITERATIONS)
02140 {
02141 INFO.print("max. number of iterations reached (exiting loop).");
02142 DONE = 1;
02143 }
02144
02145
02146 if (DONE == 1)
02147 {
02148
02149 if (overallmodeltestL > 10)
02150 {
02151 WARNING << "coregpm: overallmodeltest Lines = " << overallmodeltestL << ends;
02152 WARNING.print();
02153 WARNING << " is larger than 10. (Suggest model or a priori sigma not correct.)";
02154 WARNING.print();
02155 }
02156
02157 if (overallmodeltestP > 10)
02158 {
02159 WARNING << "coregpm: overallmodeltest Pixels = " << overallmodeltestP;
02160 WARNING.print();
02161 WARNING << " is larger than 10.\n(suggests a priori sigma not correct.)";
02162 WARNING.print();
02163 }
02164
02165 if (max(maxwL,maxwP)>200.0)
02166 {
02167 WARNING << "Recommendation: remove window number: " << Data(winL,0)
02168 << " and re-run step COREGPM. max. wtest is: "
02169 << max(maxwL,maxwP) << ".";
02170 WARNING.print();
02171 }
02172
02173 real8 expected = 1. / sqr(28);
02174
02175 expected /= sqr(SIGMAL);
02176 for (i=0; i<Nobs; i++)
02177 if (Qy_hat(i,i) > expected)
02178 {
02179 WARNING << "coregpm: Qy_hat larger than it should for window number: "
02180 << Data(i,0);
02181 WARNING.print();
02182 }
02183 }
02184
02185
02186 ITERATION++;
02187 }
02188
02189
02190
02191
02192 ofstream cpmdata("CPM_Data", ios::out | ios::trunc);
02193 bk_assert(cpmdata,"coregpm: CPM_DATA",__FILE__,__LINE__);
02194 cpmdata << "File: CPM_Data"
02195 << "\nThis file contains information on the least squares"
02196 << "\n estimation of the coregistration parameters."
02197 << "\nThis info is used in the plotcmp script."
02198 << "\nThere are 10 columns with:"
02199 << "\nWindow number, position L, position P, "
02200 << "\n offsetL (observation), offsetP (observation), correlation,"
02201 << "\n estimated errorL, errorP, w-test statistics for L, P."
02202 << "\nwin posL posP offL offP corr eL eP wtstL wtstP"
02203 << "\n------------------------------------------------------------\n";
02204 cpmdata.close();
02205
02206
02207 FILE *cpm;
02208 cpm=fopen("CPM_Data","a");
02209
02210 for (i=0; i<Data.lines(); i++)
02211 fprintf(cpm,
02212 "%4.0f %5.0f %5.0f %# 9.2f %# 9.2f %# 6.2f %6.2f %6.2f %6.2f %6.2f\n",
02213 Data(i,0), Data(i,1), Data(i,2),
02214 Data(i,3), Data(i,4), Data(i,5),
02215 eL_hat(i,0), eP_hat(i,0),
02216 abs(wtestL(i,0)), abs(wtestP(i,0)));
02217 fclose(cpm);
02218
02219
02220
02221 ofstream scratchlogfile("scratchlogcpm", ios::out | ios::trunc);
02222 bk_assert(scratchlogfile,"coregpm: scratchlogcpm",__FILE__,__LINE__);
02223 scratchlogfile
02224 << "\n\n*******************************************************************"
02225 << "\n*_Start_" << processcontrol[pr_i_coregpm]
02226 << "\n*******************************************************************"
02227 << "\nA polynomial model is weighted least squares estimated"
02228 << "\nfor azimuth and range through the FINE offset vectors."
02229 << "\nThe number of coefficient are the unknowns, the number of"
02230 << "\nobservations are the offset vectors above the THRESHOLD"
02231 << "\nspecified in the input file. To estimate the unknowns, at"
02232 << "\nleast the number of observations must equal the number of unknowns."
02233 << "\nIf there are more observations, we can statistically test"
02234 << "\nwhether the observations fit the model, and whether there are"
02235 << "\noutliers in the observations, which we like to remove."
02236 << "\nThe overall model test does the first. Wtest the second."
02237 << "\nWe advice to remove some bad estimated offsets by hand based"
02238 << "\nthe largest w-test, and to iterate running this step until"
02239 << "\nno outlier is identified anymore. A great tool is plotting"
02240 << "\nthe observations and errors, which can be done with the utility"
02241 << "\nscripts provided by Doris (calls to GMT)."
02242 << "\nAlso see any book on LS methods."
02243 << "\n\nDegree of model:\t\t\t\t"
02244 << DEGREE
02245 << "\nThreshold on data (correlation):\t\t\t"
02246 << THRESHOLD
02247 << "\nOversmaplings factor used in fine: \t"
02248 << oversamplingsfactorfine
02249 << "\nThis means maximum can be found within [samples]: \t"
02250 << ACCURACY
02251 << "\nA priori sigma azimuth (based on experience): \t"
02252 << SIGMAL
02253 << "\nA priori sigma range (based on experience): \t"
02254 << SIGMAP
02255 << "\nNumber of observations: \t\t\t"
02256 << Data.lines()
02257 << "\nNumber of rejected observations: \t\t\t"
02258 << ITERATION
02259 << "\nNumber of unknowns: \t\t\t\t"
02260 << Nunk
02261 << "\nOverall model test in Azimuth direction: \t"
02262 << overallmodeltestL
02263 << "\nOverall model test in Range direction: \t\t"
02264 << overallmodeltestP
02265 << "\nLargest w test statistic in Azimuth direction: \t"
02266 << maxwL
02267 << "\n for window number: \t\t\t\t"
02268 << Data(winL,0)
02269 << "\nLargest w test statistic in Range direction: \t"
02270 << maxwP
02271 << "\n for window number: \t\t\t\t"
02272 << Data(winP,0)
02273 << "\nMaximum deviation from unity Normalmatrix*Covar(unknowns): \t"
02274 << maxdev
02275 << "\nEstimated parameters in Azimuth direction"
02276 << "\nx_hat \tstd"
02277 << "\n(a00 | a10 a01 | a20 a11 a02 | a30 a21 a12 a03 | ...)\n";
02278 for (i=0; i<Nunk; i++)
02279 scratchlogfile
02280 << setiosflags(ios::fixed)
02281 << setiosflags(ios::showpoint)
02282 << setiosflags(ios::right)
02283 << setw(8) << setprecision(4)
02284 << rhsL(i,0) << " \t" << sqrt(Qx_hat(i,i)) << endl;
02285 scratchlogfile << "\nEstimated parameters in Range direction"
02286 << "\n(b00 | b10 b01 | b20 b11 b02 | b30 b21 b12 b03 | ...)\n";
02287 for (i=0; i<Nunk; i++)
02288 scratchlogfile
02289 << setiosflags(ios::fixed)
02290 << setiosflags(ios::showpoint)
02291 << setiosflags(ios::right)
02292 << setw(8) << setprecision(4)
02293 << rhsP(i,0) << " \t" << Qx_hat(i,i) << endl;
02294
02295 scratchlogfile << "\nCovariance matrix estimated parameters:"
02296 << "\n---------------------------------------\n";
02297 for (i=0; i<Nunk; i++)
02298 {
02299 for (j=0; j<Nunk; j++)
02300 {
02301 scratchlogfile
02302 << setiosflags(ios::fixed)
02303 << setiosflags(ios::showpoint)
02304 << setiosflags(ios::right)
02305 << setw(8) << setprecision(4)
02306 << Qx_hat(i,j) << " ";
02307 }
02308 scratchlogfile << endl;
02309 }
02310
02311 scratchlogfile << "\n*******************************************************************\n";
02312 scratchlogfile.close();
02313
02314
02315 ofstream scratchresfile("scratchrescpm", ios::out | ios::trunc);
02316 bk_assert(scratchresfile,"coregpm: scratchrescpm",__FILE__,__LINE__);
02317
02318 scratchresfile.setf(ios::scientific, ios::floatfield);
02319 scratchresfile.setf(ios::right, ios::adjustfield);
02320 scratchresfile.precision(8);
02321 scratchresfile.width(18);
02322
02323 scratchresfile
02324 << "\n\n*******************************************************************"
02325 << "\n*_Start_" << processcontrol[pr_i_coregpm]
02326 << "\n*******************************************************************"
02327 << "\nDegree_cpm:\t" << DEGREE
02328 << "\nEstimated_coefficientsL:\n";
02329 int32 coeffL = 0;
02330 int32 coeffP = 0;
02331 for (i=0; i<Nunk; i++)
02332 {
02333 if (rhsL(i,0) < 0.)
02334 scratchresfile << rhsL(i,0);
02335 else
02336 scratchresfile << " " << rhsL(i,0);
02337
02338
02339 scratchresfile << " \t" << coeffL << " " << coeffP << "\n";
02340 coeffL--;
02341 coeffP++;
02342 if (coeffL == -1)
02343 {
02344 coeffL = coeffP;
02345 coeffP = 0;
02346 }
02347 }
02348
02349 coeffL = 0;
02350 coeffP = 0;
02351 scratchresfile << "\nEstimated_coefficientsP:\n";
02352 for (i=0; i<Nunk; i++)
02353 {
02354 if (rhsP(i,0) < 0.)
02355 scratchresfile << rhsP(i,0);
02356 else
02357 scratchresfile << " " << rhsP(i,0);
02358
02359
02360 scratchresfile << " \t" << coeffL << " " << coeffP << "\n";
02361 coeffL--;
02362 coeffP++;
02363 if (coeffL == -1)
02364 {
02365 coeffL = coeffP;
02366 coeffP = 0;
02367 }
02368 }
02369 scratchresfile << "*******************************************************************"
02370
02371 << "\n* End_" << processcontrol[pr_i_coregpm] << "_NORMAL"
02372 << "\n*******************************************************************\n";
02373 scratchresfile.close();
02374
02375
02376
02377
02378
02379 matrix<real8> Lcoeff = readcoeff("scratchrescpm",
02380 "Estimated_coefficientsL:",Ncoeffs(DEGREE));
02381
02382 matrix<real8> Pcoeff = readcoeff("scratchrescpm",
02383 "Estimated_coefficientsP:",Ncoeffs(DEGREE));
02384 matrix<real4> x_axis(2,1);
02385 matrix<real4> y_axis(2,1);
02386 x_axis(0,0) = minL;
02387 x_axis(1,0) = maxL;
02388 y_axis(0,0) = minP;
02389 y_axis(1,0) = maxP;
02390 normalize(x_axis,minL,maxL);
02391 normalize(y_axis,minP,maxP);
02392 matrix<real4> offsetcornersL = polyval(x_axis,y_axis,Lcoeff);
02393 matrix<real4> offsetcornersP = polyval(x_axis,y_axis,Pcoeff);
02394 INFO.print(" ");
02395 INFO.print("Modeled transformation in azimuth:");
02396 INFO.print("-------------------------------------------------");
02397 INFO << " First line: " << offsetcornersL(0,0) << " ... " << offsetcornersL(0,1);
02398 INFO.print();
02399 INFO.print(" : :");
02400 INFO << " Last line: " << offsetcornersL(1,0) << " ... " << offsetcornersL(1,1);
02401 INFO.print();
02402 INFO.print("\n");
02403 INFO.print("Modeled transformation in range:");
02404 INFO.print("-------------------------------------------------");
02405 INFO << " First line: " << offsetcornersP(0,0) << " ... " << offsetcornersP(0,1);
02406 INFO.print();
02407 INFO.print(" : :");
02408 INFO << " Last line: " << offsetcornersP(1,0) << " ... " << offsetcornersP(1,1);
02409 INFO.print();
02410 INFO.print(" ");
02411
02412
02413
02414
02415 if (coregpminput.dumpmodel)
02416 {
02417 DEBUG.print("Do evaluation of coreg model with stepsize 100 pixels or so...");
02418 DEBUG.print("And account for currentwindow, not orig window...");
02419 PROGRESS.print("Started dumping evaluated model azimuth.");
02420 TRACE.print();
02421 TRACE << "offsetazi_" << originalmaster.lines()
02422 << "_" << originalmaster.pixels()
02423 << ".r4";
02424 char fileazi[ONE27];
02425 strcpy(fileazi,TRACE.get_str());
02426 TRACE.print();
02427
02428 ofstream dumpfile;
02429 openfstream(dumpfile,fileazi,true);
02430 bk_assert(dumpfile,fileazi,__FILE__,__LINE__);
02431
02432
02433
02434 matrix<real4> l_axis(1,1);
02435 matrix<real4> p_axis(originalmaster.pixels(),1);
02436 for (i=0; i<p_axis.pixels(); ++i)
02437 p_axis(i,0) = i+originalmaster.pixlo;
02438 normalize(p_axis,minP,maxP);
02439
02440
02441 for (i =originalmaster.linelo;
02442 i<=originalmaster.linehi; ++i)
02443 {
02444
02445
02446 l_axis(0,0) = normalize(real8(i),minL,maxL);
02447
02448 matrix<real4> MODEL = polyval(l_axis,p_axis,Lcoeff);
02449 dumpfile << MODEL;
02450 }
02451 dumpfile.close();
02452 INFO << "Dumped model azimuth offset to file: " << fileazi
02453 << " format: real4; number of lines: "
02454 << originalmaster.lines()
02455 << " number of pixels: "
02456 << originalmaster.pixels();
02457 INFO.print();
02458
02459
02460 PROGRESS.print("Started dumping evaluated model range.");
02461 TRACE.print();
02462 TRACE << "offsetrange_" << originalmaster.lines()
02463 << "_" << originalmaster.pixels() << ".r4";
02464 char filerange[ONE27];
02465 strcpy(filerange,TRACE.get_str());
02466 TRACE.print();
02467 ofstream dumpfile2;
02468 openfstream(dumpfile2,filerange,true);
02469 bk_assert(dumpfile2,filerange,__FILE__,__LINE__);
02470
02471 for (i =originalmaster.linelo;
02472 i<=originalmaster.linehi; ++i)
02473 {
02474 l_axis(0,0) = normalize(real8(i),minL,maxL);
02475 matrix<real4> MODEL = polyval(l_axis,p_axis,Pcoeff);
02476 dumpfile2 << MODEL;
02477 }
02478 INFO << "Dumped model range offset to file: " << filerange
02479 << " format: real4; number of lines: "
02480 << originalmaster.lines()
02481 << " number of pixels: "
02482 << originalmaster.pixels();
02483 INFO.print();
02484 dumpfile2.close();
02485 }
02486
02487
02488 PROGRESS.print("finished computation of coregistration parameters.");
02489 }
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507 matrix<real4> getofffile(
02508 const char* file,
02509 real4 threshold)
02510 {
02511 TRACE_FUNCTION("getofffile (BK 24-Feb-1999)");
02512 char dummyline[ONE27];
02513 char word[EIGHTY];
02514 bool foundsection = false;
02515
02516 ifstream infile;
02517 openfstream(infile,file);
02518 bk_assert(infile,file,__FILE__,__LINE__);
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528 while (infile)
02529 {
02530 infile >> word;
02531 if (strcmp("Number_of_correlation_windows:",word))
02532 {
02533 infile.getline(dummyline,ONE27,'\n');
02534 }
02535 else
02536 {
02537 foundsection=true;
02538 int32 N;
02539 infile >> N;
02540 infile.getline(dummyline,ONE27,'\n');
02541 infile.getline(dummyline,ONE27,'\n');
02542 int32 pos = infile.tellg();
02543 int32 Nobs = 0;
02544 real4 winnumber, posL, posP, offL, offP, corr;
02545 register int32 i;
02546 for (i=0;i<N;i++)
02547 {
02548 infile >> winnumber >> posL >> posP >> offL >> offP >> corr;
02549 infile.getline(dummyline,ONE27,'\n');
02550 if (corr > threshold)
02551 Nobs++;
02552 }
02553
02554 if (Nobs == 0)
02555 {
02556 PRINT_ERROR("code ???: No data found > threshold.")
02557 throw(some_error);
02558 }
02559
02560 matrix<real4> Data(Nobs,6);
02561 infile.seekg(pos);
02562 int32 cnti = -1;
02563 for (i=0;i<N;i++)
02564 {
02565 infile >> winnumber >> posL >> posP >> offL >> offP >> corr;
02566 infile.getline(dummyline,ONE27,'\n');
02567 if (corr > threshold)
02568 {
02569 cnti++;
02570 Data(cnti,0) = winnumber;
02571 Data(cnti,1) = posL;
02572 Data(cnti,2) = posP;
02573 Data(cnti,3) = offL;
02574 Data(cnti,4) = offP;
02575 Data(cnti,5) = corr;
02576 }
02577 }
02578
02579 infile.close();
02580 return Data;
02581 }
02582 }
02583
02584
02585 if (!foundsection)
02586 {
02587 PRINT_ERROR("code 401: getofffile: couldn't find data section in file.");
02588 throw(some_error);
02589 }
02590 infile.close();
02591
02592
02593 return matrix<real4>(999,999);
02594 }
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610 matrix<real4> cc4(
02611 const matrix<real4> &x)
02612 {
02613 TRACE_FUNCTION("cc4 (BK 16-Mar-1999)");
02614 if (x.pixels() != 1)
02615 {
02616 PRINT_ERROR("cc4: standing vectors only.")
02617 throw(input_error);
02618 }
02619
02620 real4 alpha = -1.;
02621 matrix<real4> y(x.lines(),1);
02622 for (register int32 i=0;i<y.lines();i++)
02623 {
02624 real4 xx2 = sqr(x(i,0));
02625 real4 xx = sqrt(xx2);
02626 if (xx < 1)
02627 y(i,0) = (alpha+2)*xx2*xx - (alpha+3)*xx2 + 1;
02628 else if (xx < 2)
02629 y(i,0) = alpha*xx2*xx - 5*alpha*xx2 + 8*alpha*xx - 4*alpha;
02630 else
02631 y(i,0) = 0.;
02632 }
02633 return y;
02634 }
02635
02636
02637
02638
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652 matrix<real4> cc6(
02653 const matrix<real4> &x)
02654 {
02655 TRACE_FUNCTION("cc6 (BK 16-Mar-1999)");
02656 if (x.pixels() != 1)
02657 {
02658 PRINT_ERROR("cc6: standing vectors only.")
02659 throw(input_error);
02660 }
02661
02662 real4 alpha = -.5;
02663 real4 beta = .5;
02664 matrix<real4> y(x.lines(),1);
02665 for (register int32 i=0;i<y.lines();i++)
02666 {
02667 real4 xx2 = sqr(x(i,0));
02668 real4 xx = sqrt(xx2);
02669 if (xx < 1)
02670 y(i,0) = (alpha-beta+2)*xx2*xx - (alpha-beta+3)*xx2 + 1;
02671
02672 else if (xx < 2)
02673 y(i,0) = alpha*xx2*xx - (5*alpha-beta)*xx2
02674 + (8*alpha-3*beta)*xx - (4*alpha-2*beta);
02675 else if (xx < 3)
02676 y(i,0) = beta*xx2*xx - 8*beta*xx2 + 21*beta*xx - 18*beta;
02677 else
02678 y(i,0) = 0.;
02679 }
02680 return y;
02681 }
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696 matrix<real4> ts6(
02697 const matrix<real4> &x)
02698 {
02699 TRACE_FUNCTION("ts6 (BK 16-Mar-1999)");
02700 if (x.pixels() != 1)
02701 {
02702 PRINT_ERROR("ts6: standing vectors only.")
02703 throw(input_error);
02704 }
02705
02706 matrix<real4> y(x.lines(),1);
02707 for (register int32 i=0;i<y.lines();i++)
02708 y(i,0) = sinc(x(i,0)) * rect(x(i,0)/6.);
02709 return y;
02710 }
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725 matrix<real4> ts8(
02726 const matrix<real4> &x)
02727 {
02728 TRACE_FUNCTION("ts8 (BK 16-Mar-1999)");
02729 if (x.pixels() != 1)
02730 {
02731 PRINT_ERROR("ts8: standing vectors only.")
02732 throw(input_error);
02733 }
02734
02735 matrix<real4> y(x.lines(),1);
02736 for (register int32 i=0;i<y.lines();i++)
02737 y(i,0) = sinc(x(i,0)) * rect(x(i,0)/8.);
02738 return y;
02739 }
02740
02741
02742
02743
02744
02745
02746
02747
02748
02749
02750
02751
02752
02753
02754 matrix<real4> ts16(
02755 const matrix<real4> &x)
02756 {
02757 TRACE_FUNCTION("ts16 (BK 16-Mar-1999)");
02758 if (x.pixels() != 1)
02759 {
02760 PRINT_ERROR("ts16: standing vectors only.")
02761 throw(input_error);
02762 }
02763 matrix<real4> y(x.lines(),1);
02764 for (register int32 i=0;i<y.lines();i++)
02765 y(i,0) = sinc(x(i,0)) * rect(x(i,0)/16.);
02766 return y;
02767 }
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778
02779
02780
02781
02782 matrix<real4> rect(
02783 const matrix<real4> &x)
02784 {
02785 TRACE_FUNCTION("rect (BK 16-Mar-1999)");
02786 if (x.pixels() != 1)
02787 {
02788 PRINT_ERROR("rect: standing vectors only.");
02789 throw(input_error);
02790 }
02791
02792 matrix<real4> y(x.lines(),1);
02793 for (register int32 i=0;i<y.lines();i++)
02794 y(i,0) = rect(x(i,0));
02795 return y;
02796 }
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807
02808
02809
02810
02811 matrix<real4> tri(
02812 const matrix<real4> &x)
02813 {
02814 TRACE_FUNCTION("tri (BK 16-Mar-1999)")
02815 if (x.pixels() != 1)
02816 {
02817 PRINT_ERROR("tri: standing vectors only.")
02818 throw(input_error);
02819 }
02820 matrix<real4> y(x.lines(),1);
02821 for (register int32 i=0;i<y.lines();i++)
02822 y(i,0) = tri(x(i,0));
02823 return y;
02824 }
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838
02839
02840
02841
02842
02843
02844 matrix<real4> knab(
02845 const matrix<real4> &x,
02846 const real4 CHI,
02847 const int32 N)
02848 {
02849 TRACE_FUNCTION("knab (BK 22-Dec-2003)");
02850 if (x.pixels() != 1)
02851 {
02852 PRINT_ERROR("knab: standing vectors only.")
02853 throw(input_error);
02854 }
02855 matrix<real4> y(x.lines(),1);
02856 real4 v = 1.0-1.0/CHI;
02857 real4 vv = PI*v*real4(N)/2.0;
02858 for (register int32 i=0;i<y.lines();i++)
02859 y(i,0) = sinc(x(i,0))*cosh(vv*sqrt(1.0-sqr(2.0*x(i,0)/real4(N))))/cosh(vv);
02860 return y;
02861 }
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882
02883
02884
02885
02886
02887
02888
02889 void resample(
02890 const input_gen &generalinput,
02891 const input_resample &resampleinput,
02892 const slcimage &master,
02893 const slcimage &slave,
02894 const matrix<real8> &cpmL,
02895 const matrix<real8> &cpmP
02896 )
02897 {
02898 TRACE_FUNCTION("resample (BK 16-Mar-1999; BK 09-Nov-2000)")
02899 if (resampleinput.shiftazi==true)
02900 DEBUG.print("shifting kernelL to data fDC BK 26-Oct-2002");
02901
02902 const uint BUFFERMEMSIZE = generalinput.memory;
02903 const int32 Npoints = resampleinput.method%100;
02904 if (isodd(Npoints))
02905 {
02906 PRINT_ERROR("resample only even point interpolators.")
02907 throw(input_error);
02908 }
02909 const int32 Npointsd2 = Npoints/2;
02910 const int32 Npointsd2m1 = Npointsd2-1;
02911 const int32 Npointsm1 = Npoints-1;
02912 const uint Sfilelines = slave.currentwindow.lines();
02913 const uint sizeofci16 = sizeof(compli16);
02914 const uint sizeofcr4 = sizeof(complr4);
02915
02916
02917 const real8 minL = master.originalwindow.linelo;
02918 const real8 maxL = master.originalwindow.linehi;
02919 const real8 minP = master.originalwindow.pixlo;
02920 const real8 maxP = master.originalwindow.pixhi;
02921 INFO << "resample: polynomial normalized by factors: "
02922 << minL << " " << maxL << " " << minP << " " << maxP << " to [-2,2]";
02923 INFO.print();
02924
02925
02926 real4 CHI1 = slave.prf/slave.abw;
02927 real4 CHI2 = (slave.rsr2x/2.0)/slave.rbw;
02928 real4 CHI = (CHI1+CHI2)/2.0;
02929 DEBUG << "Oversampling factor azimuth: " << CHI1;
02930 DEBUG.print();
02931 DEBUG << "Oversampling factor azimuth: " << CHI2;
02932 DEBUG.print();
02933 DEBUG << "Mean oversampling factor data: " << CHI;
02934 DEBUG.print();
02935 if (CHI < 1.1)
02936 {
02937 WARNING << "mean oversampling factor data " << CHI << " not optimal for KNAB";
02938 WARNING.print();
02939 }
02940
02941
02942
02943
02944
02945
02946
02947
02948
02949 register int32 i;
02950
02951
02952
02953 const int32 INTERVAL = 100;
02954 const real8 dx = 1./INTERVAL;
02955 const int32 Ninterval = INTERVAL + 1;
02956 INFO << "resample: lookup table size: " << Ninterval;
02957 INFO.print();
02958
02959 matrix<real4> x_axis(Npoints,1);
02960 for (i=0; i<Npoints; ++i)
02961 x_axis(i,0) = 1. - Npointsd2 + i;
02962
02963
02964
02965
02966
02967 matrix<complr4> *pntM[Ninterval];
02968
02969
02970 matrix<real4> *pntAxis[Ninterval];
02971 for (i=0; i<Ninterval; ++i)
02972 {
02973 pntM[i] = new matrix<complr4> (Npoints,1);
02974 pntAxis[i] = new matrix<real4> (Npoints,1);
02975 switch(resampleinput.method)
02976 {
02977 case rs_rect:
02978 (*pntM[i]) = mat2cr4(rect(x_axis));
02979 break;
02980 case rs_tri:
02981 (*pntM[i]) = mat2cr4(tri(x_axis));
02982 break;
02983 case rs_cc4p:
02984 (*pntM[i]) = mat2cr4(cc4(x_axis));
02985 break;
02986 case rs_cc6p:
02987 (*pntM[i]) = mat2cr4(cc6(x_axis));
02988 break;
02989 case rs_ts6p:
02990 (*pntM[i]) = mat2cr4(ts6(x_axis));
02991 break;
02992 case rs_ts8p:
02993 (*pntM[i]) = mat2cr4(ts8(x_axis));
02994 break;
02995 case rs_ts16p:
02996 (*pntM[i]) = mat2cr4(ts16(x_axis));
02997 break;
02998 case rs_knab4p:
02999 (*pntM[i]) = mat2cr4(knab(x_axis,CHI,4));
03000 case rs_knab6p:
03001 (*pntM[i]) = mat2cr4(knab(x_axis,CHI,6));
03002 case rs_knab8p:
03003 (*pntM[i]) = mat2cr4(knab(x_axis,CHI,8));
03004 case rs_knab10p:
03005 (*pntM[i]) = mat2cr4(knab(x_axis,CHI,10));
03006 case rs_knab16p:
03007 (*pntM[i]) = mat2cr4(knab(x_axis,CHI,16));
03008 break;
03009 default:
03010 PRINT_ERROR("impossible.")
03011 throw(unhandled_case_error);
03012 }
03013 (*pntAxis[i]) = x_axis;
03014 x_axis -= dx;
03015 }
03016
03017
03018 PROGRESS.print("Resample: lookup table created (kernel and axis).");
03019
03020
03021 int32 degree_cpmL = degree(cpmL.size());
03022 int32 degree_cpmP = degree(cpmP.size());
03023
03024
03025
03026
03027
03028 window approxoverlap = getoverlap(master,slave,cpmL,cpmP);
03029
03030
03031 real8 o00L = approxoverlap.linelo
03032 + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
03033 normalize(real8(approxoverlap.pixlo),minP,maxP),
03034 cpmL,degree_cpmL);
03035 real8 o00P = approxoverlap.pixlo
03036 + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
03037 normalize(real8(approxoverlap.pixlo),minP,maxP),
03038 cpmP,degree_cpmP);
03039 real8 o0NL = approxoverlap.linelo
03040 + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
03041 normalize(real8(approxoverlap.pixhi),minP,maxP),
03042 cpmL,degree_cpmL);
03043 real8 o0NP = approxoverlap.pixhi
03044 + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
03045 normalize(real8(approxoverlap.pixhi),minP,maxP),
03046 cpmP,degree_cpmP);
03047 real8 oN0L = approxoverlap.linehi
03048 + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
03049 normalize(real8(approxoverlap.pixlo),minP,maxP),
03050 cpmL,degree_cpmL);
03051 real8 oN0P = approxoverlap.pixlo
03052 + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
03053 normalize(real8(approxoverlap.pixlo),minP,maxP),
03054 cpmP,degree_cpmP);
03055 real8 oNNL = approxoverlap.linehi
03056 + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
03057 normalize(real8(approxoverlap.pixhi),minP,maxP),
03058 cpmL,degree_cpmL);
03059 real8 oNNP = approxoverlap.pixhi
03060 + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
03061 normalize(real8(approxoverlap.pixhi),minP,maxP),
03062 cpmP,degree_cpmP);
03063
03064
03065
03066
03067
03068
03069
03070 int32 correctlinelo = slave.currentwindow.linelo
03071 - (int32(ceil(min(o00L,o0NL)) - Npointsd2));
03072 if (correctlinelo < 0) correctlinelo = 0;
03073
03074 int32 correctpixlo = slave.currentwindow.pixlo
03075 - (int32(ceil(min(o00P,oN0P)) - Npointsd2));
03076 if (correctpixlo < 0) correctpixlo = 0;
03077
03078 int32 correctlinehi = slave.currentwindow.linehi
03079 - (int32(floor(max(oN0L,oNNL)) + Npointsd2));
03080 if (correctlinehi > 0) correctlinehi = 0;
03081
03082 int32 correctpixhi = slave.currentwindow.pixhi
03083 - (int32(floor(max(o0NP,oNNP)) + Npointsd2));
03084 if (correctpixhi > 0) correctpixhi = 0;
03085
03086
03087
03088 window overlap(approxoverlap.linelo + correctlinelo,
03089 approxoverlap.linehi + correctlinehi,
03090 approxoverlap.pixlo + correctpixlo,
03091 approxoverlap.pixhi + correctpixhi);
03092
03093
03094 int32 write0lines1 = 0;
03095 int32 write0linesN = 0;
03096 int32 write0pixels1 = 0;
03097 int32 write0pixelsN = 0;
03098 if (!(resampleinput.dbow.linelo == 0 &&
03099 resampleinput.dbow.linehi == 0 &&
03100 resampleinput.dbow.pixlo == 0 &&
03101 resampleinput.dbow.pixhi == 0 ))
03102 {
03103
03104 if (resampleinput.dbow.linelo > overlap.linehi)
03105 {
03106 PRINT_ERROR("RS_DBOW: specified min. line larger than max. line of overlap.")
03107 throw(input_error);
03108 }
03109 if (resampleinput.dbow.linehi < overlap.linelo)
03110 {
03111 PRINT_ERROR("RS_DBOW: specified max. line smaller than min. line of overlap.")
03112 throw(input_error);
03113 }
03114 if (resampleinput.dbow.pixlo > overlap.pixhi)
03115 {
03116 PRINT_ERROR("RS_DBOW: specified min. pixel larger than max. pixel of overlap.")
03117 throw(input_error);
03118 }
03119 if (resampleinput.dbow.pixhi < overlap.pixlo)
03120 {
03121 PRINT_ERROR("RS_DBOW: specified max. pixel smaller than min. pixel of overlap.")
03122 throw(input_error);
03123 }
03124
03125 write0lines1 = overlap.linelo - resampleinput.dbow.linelo;
03126 if ( write0lines1 < 0 ) write0lines1 = 0;
03127 write0linesN = -overlap.linehi + resampleinput.dbow.linehi;
03128 if ( write0linesN < 0 ) write0linesN = 0;
03129 write0pixels1 = overlap.pixlo - resampleinput.dbow.pixlo;
03130 if ( write0pixels1 < 0 ) write0pixels1 = 0;
03131 write0pixelsN = -overlap.pixhi + resampleinput.dbow.pixhi;
03132 if ( write0pixelsN < 0 ) write0pixelsN = 0;
03133 if (resampleinput.dbow.linelo < overlap.linelo)
03134 {
03135 WARNING << "RS_DBOW: min. line < overlap (writing: " << write0lines1
03136 << " lines with zeros before first resampled line).";
03137 WARNING.print();
03138 }
03139 else
03140 overlap.linelo = resampleinput.dbow.linelo;
03141 if (resampleinput.dbow.linehi > overlap.linehi)
03142 {
03143 WARNING << "RS_DBOW: max. line > overlap (writing: " << write0linesN
03144 << " lines with zeros after last resampled line).";
03145 WARNING.print();
03146 }
03147 else
03148 overlap.linehi = resampleinput.dbow.linehi;
03149 if (resampleinput.dbow.pixlo < overlap.pixlo)
03150 {
03151 WARNING << "RS_DBOW: min. pixel < overlap (writing: " << write0pixels1
03152 << " columns with zeros before first resampled column).";
03153 WARNING.print();
03154 }
03155 else
03156 overlap.pixlo = resampleinput.dbow.pixlo;
03157 if (resampleinput.dbow.pixhi > overlap.pixhi)
03158 {
03159 WARNING << "RS_DBOW: max. pixel > overlap (writing: " << write0pixelsN
03160 << " columns with zeros after last resampled column).";
03161 WARNING.print();
03162 }
03163 else
03164 overlap.pixhi = resampleinput.dbow.pixhi;
03165 }
03166
03167
03168
03169 const int32 Npointsxsize = Npoints*sizeofcr4;
03170 const int32 npixels = slave.currentwindow.pixels();
03171 const real8 bytesperline = sizeofcr4 * npixels;
03172
03173
03174 const real8 bigmatrices = 2.5;
03175 const int32 nlines = int32((BUFFERMEMSIZE/bigmatrices)/bytesperline);
03176
03177
03178 matrix<complr4> BUFFER;
03179 matrix<complr4> RESULT(nlines,overlap.pixhi-overlap.pixlo+1);
03180 matrix<complr4> PART(Npoints,Npoints);
03181
03182 #ifdef __USE_VECLIB_LIBRARY__
03183 matrix<complr4> TMPRES(Npoints,1);
03184 int32 Np = Npoints;
03185 int32 ONEint = 1;
03186 complr4 c4alpha(1.,0.);
03187 complr4 c4beta(0.,0.);
03188 STUPID_cr4 ANS;
03189 #endif
03190
03191
03192
03193 ofstream ofile;
03194 openfstream(ofile,resampleinput.fileout,generalinput.overwrit);
03195 bk_assert(ofile,resampleinput.fileout,__FILE__,__LINE__);
03196
03197
03198 register int32 thisline, thispixel;
03199 switch (resampleinput.oformatflag)
03200 {
03201 case FORMATCR4:
03202 {
03203 const complr4 zerocr4(0,0);
03204 for (thisline=0; thisline<write0lines1; ++thisline)
03205 for (thispixel=0;
03206 thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
03207 ++thispixel)
03208 ofile.write((char*)&zerocr4,sizeofcr4);
03209 break;
03210 }
03211 case FORMATCI2:
03212 {
03213 const compli16 zeroci16(0,0);
03214 for (thisline=0; thisline<write0lines1; ++thisline)
03215 for (thispixel=0;
03216 thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
03217 ++thispixel)
03218 ofile.write((char*)&zeroci16,sizeofci16);
03219 break;
03220 }
03221 default:
03222 PRINT_ERROR("impossible format")
03223 throw(unhandled_case_error);
03224 }
03225
03226
03227 INFO << "Overlap window: "
03228 << overlap.linelo << ":" << overlap.linehi << ", "
03229 << overlap.pixlo << ":" << overlap.pixhi;
03230 INFO.print();
03231
03232
03233 int32 percent = 0;
03234 int32 tenpercent = int32((overlap.lines()/10)+.5);
03235 if (tenpercent==0) tenpercent = 1000;
03236
03237
03238 bool newbufferrequired = true;
03239 register int32 linecnt = -1;
03240 int32 firstline;
03241 register int32 line;
03242 register int32 pixel;
03243 for (line=overlap.linelo; line<=overlap.linehi; line++)
03244 {
03245
03246 if (((line-overlap.linelo)%tenpercent)==0)
03247 {
03248 PROGRESS << "RESAMPLE: " << setw(3) << percent << "%";
03249 PROGRESS.print();
03250 percent += 10;
03251 }
03252
03253
03254 if (linecnt==RESULT.lines()-1)
03255 {
03256 newbufferrequired = true;
03257 DEBUG << "Writing slave: ["
03258 << line-RESULT.lines() << ":" << line-1 << ", "
03259 << overlap.pixlo << ":" << overlap.pixhi
03260 << "] (master coord. system)";
03261 DEBUG.print();
03262 linecnt = 0;
03263
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278 switch (resampleinput.oformatflag)
03279 {
03280 case FORMATCR4:
03281 {
03282
03283 const complr4 zerocr4(0,0);
03284 for (thisline=0; thisline<RESULT.lines(); ++thisline)
03285 {
03286
03287 for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03288 {
03289 ofile.write((char*)&zerocr4,sizeofcr4);
03290 }
03291
03292 ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeof(RESULT(0,0)));
03293
03294 for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
03295 {
03296 ofile.write((char*)&zerocr4,sizeofcr4);
03297 }
03298 }
03299 break;
03300 }
03301 case FORMATCI2:
03302 {
03303 const compli16 zeroci16(0,0);
03304 compli16 castedresult;
03305 for (thisline=0; thisline<RESULT.lines(); ++thisline)
03306 {
03307
03308 for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03309 {
03310 ofile.write((char*)&zeroci16,sizeofci16);
03311 }
03312
03313 for (thispixel=0; thispixel<RESULT.pixels(); ++thispixel)
03314 {
03315
03316 castedresult = cr4toci2(RESULT(thisline,thispixel));
03317 ofile.write((char*)&castedresult,sizeofci16);
03318 }
03319
03320 for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
03321 {
03322 ofile.write((char*)&zeroci16,sizeofci16);
03323 }
03324 }
03325 break;
03326 }
03327 default:
03328 PRINT_ERROR("impossible format")
03329 throw(unhandled_case_error);
03330 }
03331 }
03332
03333 else
03334 {
03335 linecnt++;
03336 }
03337
03338
03339 if (newbufferrequired==true)
03340 {
03341 newbufferrequired = false;
03342 firstline = int32(ceil(min(line +
03343 polyval(normalize(real4(line),minL,maxL),
03344 normalize(real4(overlap.pixlo),minP,maxP),
03345 cpmL,degree_cpmL),
03346 line +
03347 polyval(normalize(real4(line),minL,maxL),
03348 normalize(real4(overlap.pixhi),minP,maxP),
03349 cpmL,degree_cpmL))))
03350 - Npoints;
03351 int32 line2 = line + nlines - 1;
03352 int32 lastline = int32(ceil(min(line2 +
03353 polyval(normalize(real4(line2),minL,maxL),
03354 normalize(real4(overlap.pixlo),minP,maxP),
03355 cpmL,degree_cpmL),
03356 line2 +
03357 polyval(normalize(real4(line2),minL,maxL),
03358 normalize(real4(overlap.pixhi),minP,maxP),
03359 cpmL,degree_cpmL))))
03360 + Npoints;
03361
03362 const int32 FORSURE = 25;
03363 firstline -= FORSURE;
03364 lastline += FORSURE;
03365
03366 if (firstline < int32(slave.currentwindow.linelo))
03367 firstline = slave.currentwindow.linelo;
03368 if (lastline > int32(slave.currentwindow.linehi))
03369 lastline = slave.currentwindow.linehi;
03370
03371 window winslavefile(firstline, lastline,
03372 slave.currentwindow.pixlo,
03373 slave.currentwindow.pixhi);
03374 DEBUG << "Reading slave: ["
03375 << winslavefile.linelo << ":" << winslavefile.linehi << ", "
03376 << winslavefile.pixlo << ":" << winslavefile.pixhi << "]";
03377 DEBUG.print();
03378 BUFFER = slave.readdata(winslavefile);
03379
03380
03381
03382
03383
03384 }
03385
03386
03387
03388 for (pixel=overlap.pixlo; pixel<=overlap.pixhi; pixel++)
03389 {
03390
03391
03392
03393
03394
03395 const real4 interpL = line +
03396 polyval(normalize(real4(line),minL,maxL),
03397 normalize(real4(pixel),minP,maxP),
03398 cpmL,degree_cpmL);
03399 const real4 interpP = pixel +
03400 polyval(normalize(real4(line),minL,maxL),
03401 normalize(real4(pixel),minP,maxP),
03402 cpmP,degree_cpmP);
03403
03404
03405
03406
03407
03408
03409
03410
03411 const int32 fl_interpL = int32(interpL);
03412 const int32 fl_interpP = int32(interpP);
03413 const int32 firstL = fl_interpL - Npointsd2m1;
03414 const int32 firstP = fl_interpP - Npointsd2m1;
03415 const real4 interpLdec = interpL - fl_interpL;
03416 const real4 interpPdec = interpP - fl_interpP;
03417
03418
03419
03420 int32 kernelnoL = int32(interpLdec*INTERVAL + .5);
03421 int32 kernelnoP = int32(interpPdec*INTERVAL + .5);
03422 matrix<complr4> kernelL = (*pntM[kernelnoL]);
03423 matrix<complr4> kernelP = (*pntM[kernelnoP]);
03424
03425
03426 #ifdef __DEBUG
03427
03428 if (firstL < slave.currentwindow.linelo)
03429 {
03430 WARNING.print("firstL smaller then on disk (required for interpolation). continuing");
03431 RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
03432 continue;
03433 }
03434 if (firstL+Npointsm1 > slave.currentwindow.linehi)
03435 {
03436 WARNING.print("lastL larger then on disk (required for interpolation). continuing");
03437 RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
03438 continue;
03439 }
03440 if (firstP < slave.currentwindow.pixlo)
03441 {
03442 WARNING.print("firstP smaller then on disk (required for interpolation). continuing");
03443 RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
03444 continue;
03445 }
03446 if (firstP+Npointsm1 > slave.currentwindow.pixhi)
03447 {
03448 WARNING.print("lastP larger then on disk (required for interpolation). continuing");
03449 RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
03450 continue;
03451 }
03452 #endif
03453
03454
03455
03456 if (resampleinput.shiftazi==true)
03457 {
03458
03459
03460
03461
03462
03463
03464
03465 const real4 tmp = 2.0*PI*slave.pix2fdc(interpP)/slave.prf;
03466
03467
03468
03469 for (i=0; i<Npoints; ++i)
03470 {
03471
03472
03473
03474
03475 real4 t = ((*pntAxis[kernelnoL])(i,0))*tmp;
03476 kernelL(i,0) *= complr4(cos(t),-sin(t));
03477 }
03478 }
03479
03480
03481
03482
03483
03484 for (i=0; i<Npoints; i++)
03485 memcpy(PART[i],
03486 BUFFER[i+firstL-firstline]+
03487 firstP-slave.currentwindow.pixlo,Npointsxsize);
03488
03489
03490 #ifdef __USE_VECLIB_LIBRARY__
03491
03492
03493
03494
03495 cgemv("T",&Np,&Np,&c4alpha, PART[0],&Np,
03496 kernelP[0],&ONEint,
03497 &c4beta,TMPRES[0],&ONEint,1);
03498
03499
03500
03501
03502 ANS = cdotu(&Np,TMPRES[0],&ONEint, kernelL[0],&ONEint);
03503 RESULT(linecnt,pixel-overlap.pixlo) = complr4(ANS.re,ANS.im);
03504 #else // do use VECLIB
03505
03506
03507
03508 RESULT(linecnt,pixel-overlap.pixlo) =
03509 ((matTxmat(PART*kernelP, kernelL))(0,0));
03510 #endif // VECLIB
03511 }
03512 }
03513
03514
03515
03516 DEBUG << "Writing slave: ["
03517
03518 << overlap.linehi-linecnt << ":" << overlap.linehi << ", "
03519 << overlap.pixlo << ":" << overlap.pixhi
03520 << "] (master coord. system)";
03521 DEBUG.print();
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531
03532
03533
03534
03535
03536
03537
03538
03539
03540 switch (resampleinput.oformatflag)
03541 {
03542 case FORMATCR4:
03543 {
03544 const complr4 zerocr4(0,0);
03545 for (thisline=0; thisline<=linecnt; thisline++)
03546 {
03547
03548 for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03549 {
03550 ofile.write((char*)&zerocr4,sizeofcr4);
03551 }
03552
03553 ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeofcr4);
03554
03555 for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
03556 {
03557 ofile.write((char*)&zerocr4,sizeofcr4);
03558 }
03559 }
03560 break;
03561 }
03562 case FORMATCI2:
03563 {
03564 const compli16 zeroci16(0,0);
03565 compli16 castedresult;
03566 for (thisline=0; thisline<=linecnt; thisline++)
03567 {
03568
03569 for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03570 {
03571 ofile.write((char*)&zeroci16,sizeofci16);
03572 }
03573
03574 for (thispixel=0; thispixel<RESULT.pixels(); thispixel++)
03575 {
03576 castedresult = cr4toci2(RESULT(thisline,thispixel));
03577
03578
03579
03580 ofile.write((char*)&castedresult,sizeofci16);
03581 }
03582
03583 for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
03584 {
03585 ofile.write((char*)&zeroci16,sizeofci16);
03586 }
03587 }
03588 break;
03589 }
03590 default:
03591 PRINT_ERROR("impossible format")
03592 throw(unhandled_case_error);
03593 }
03594
03595
03596
03597 switch (resampleinput.oformatflag)
03598 {
03599 case FORMATCR4:
03600 {
03601 complr4 zerocr4(0,0);
03602 for (thisline=0; thisline<write0linesN; ++thisline)
03603 for (thispixel=0;
03604 thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
03605 ++thispixel)
03606 ofile.write((char*)&zerocr4,sizeofcr4);
03607 break;
03608 }
03609 case FORMATCI2:
03610 {
03611 compli16 zeroci16(0,0);
03612 for (thisline=0; thisline<write0linesN; ++thisline)
03613 for (thispixel=0;
03614 thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
03615 ++thispixel)
03616 ofile.write((char*)&zeroci16,sizeofci16);
03617 break;
03618 }
03619 default:
03620 PRINT_ERROR("impossible format")
03621 throw(unhandled_case_error);
03622 }
03623 ofile.close();
03624
03625
03626
03627
03628 char rsmethod[EIGHTY];
03629 switch(resampleinput.method)
03630 {
03631 case rs_rect:
03632 strcpy(rsmethod,"nearest neighbour");
03633 break;
03634 case rs_tri:
03635 strcpy(rsmethod,"piecewise linear");
03636 break;
03637 case rs_cc4p:
03638 strcpy(rsmethod,"4 point cubic convolution");
03639 break;
03640 case rs_cc6p:
03641 strcpy(rsmethod,"6 point cubic convolution");
03642 break;
03643 case rs_ts6p:
03644 strcpy(rsmethod,"6 point truncated sinc");
03645 break;
03646 case rs_ts8p:
03647 strcpy(rsmethod,"8 point truncated sinc");
03648 break;
03649 case rs_ts16p:
03650 strcpy(rsmethod,"16 point truncated sinc");
03651 break;
03652 case rs_knab4p:
03653 strcpy(rsmethod,"4 point knab window");
03654 break;
03655 case rs_knab6p:
03656 strcpy(rsmethod,"6 point knab window");
03657 break;
03658 case rs_knab8p:
03659 strcpy(rsmethod,"8 point knab window");
03660 break;
03661 case rs_knab10p:
03662 strcpy(rsmethod,"10 point knab window");
03663 break;
03664 case rs_knab16p:
03665 strcpy(rsmethod,"16 point knab window");
03666 break;
03667 default:
03668 PRINT_ERROR("impossible.")
03669 throw(unhandled_case_error);
03670 }
03671
03672 char rsoformat[EIGHTY];
03673 switch(resampleinput.oformatflag)
03674 {
03675 case FORMATCR4:
03676 strcpy(rsoformat,"complex_real4");
03677 break;
03678 case FORMATCI2:
03679 strcpy(rsoformat,"complex_short");
03680 break;
03681 default:
03682 PRINT_ERROR("impossible.")
03683 throw(unhandled_case_error);
03684 }
03685
03686
03687 ofstream scratchlogfile("scratchlogresample", ios::out | ios::trunc);
03688 bk_assert(scratchlogfile,"resample: scratchlogresample",__FILE__,__LINE__);
03689
03690 scratchlogfile
03691 << "\n\n*******************************************************************"
03692 << "\n*_Start_resample"
03693 << "\nData_output_file: \t\t\t"
03694 << resampleinput.fileout
03695 << "\nData_output_format: \t\t\t"
03696 << rsoformat
03697 << "\nInterpolation kernel: \t\t\t"
03698 << rsmethod
03699 << "\nResampled slave size in master system: \t"
03700 << overlap.linelo - write0lines1 << ", "
03701 << overlap.linehi + write0linesN << ", "
03702 << overlap.pixlo - write0pixels1 << ", "
03703 << overlap.pixhi + write0pixelsN
03704
03705
03706 << "\n*******************************************************************\n";
03707 scratchlogfile.close();
03708
03709 ofstream scratchresfile("scratchresresample", ios::out | ios::trunc);
03710 bk_assert(scratchresfile,"resample: scratchresresample",__FILE__,__LINE__);
03711 scratchresfile
03712 << "\n\n*******************************************************************"
03713 << "\n*_Start_" << processcontrol[pr_s_resample]
03714 << "\n*******************************************************************"
03715 << "\nShifted azimuth spectrum: \t\t"
03716 << resampleinput.shiftazi
03717 << "\nData_output_file: \t\t"
03718 << resampleinput.fileout
03719 << "\nData_output_format: \t\t"
03720 << rsoformat
03721 << "\nInterpolation kernel: \t\t"
03722 << rsmethod
03723 << "\nFirst_line (w.r.t. original_master): \t\t"
03724 << overlap.linelo - write0lines1
03725 << "\nLast_line (w.r.t. original_master): \t\t"
03726 << overlap.linehi + write0linesN
03727 << "\nFirst_pixel (w.r.t. original_master): \t\t"
03728 << overlap.pixlo - write0pixels1
03729 << "\nLast_pixel (w.r.t. original_master): \t\t"
03730 << overlap.pixhi + write0pixelsN
03731 << "\n*******************************************************************"
03732 << "\n* End_" << processcontrol[pr_s_resample] << "_NORMAL"
03733 << "\n*******************************************************************\n";
03734 scratchresfile.close();
03735
03736
03737
03738
03739
03740
03741 DEBUG.print("Exiting resample.");
03742
03743 }
03744