Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

coregistration.cc

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /****************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/coregistration.cc,v $     *
00032  * $Revision: 3.30 $                                            *
00033  * $Date: 2005/04/11 13:47:45 $                                 *
00034  * $Author: kampes $                                            *
00035  *                                                              *
00036  * -coarse based on orbits.                                     *
00037  * -coarse based on correlation with fft/in space domain.       *
00038  * -fine coregistration offset vector computation.              *
00039  * -computation coregistration parameters.                      *
00040  * -computation flat earth correction.                          *
00041  * -resampling of slave to master grid.                         *
00042  ****************************************************************/
00043 
00044 
00045 #include "constants.hh"         // typedefs etc.
00046 #include "ioroutines.hh"        // error function etc.
00047 #include "utilities.hh"         // isodd
00048 #include "coregistration.hh"    // header file
00049 #include "conversion.hh"        // lp2xyz etc.
00050 #include "refsystems.hh"        // for dsk WGS84
00051 #include "orbitbk.hh"           // my orbit class
00052 #include "slcimage.hh"          // my slc image class
00053 #include "productinfo.hh"       // my 'products' class
00054 #include "exceptions.hh"        // my exceptions class
00055 #include "bk_baseline.hh"       // my exceptions class
00056 
00057 
00058 #include <iomanip>              // setw
00059 #include <cstdlib>              // system
00060 #include <cmath>                // sqrt
00061 #include <algorithm>            // max
00062 #include <cstdio>               // some compilers, remove function
00063 
00064 
00065 
00066 /****************************************************************
00067  *    coarseporbit                                              *
00068  *                                                              *
00069  * computes translation of slave w.r.t. master                  *
00070  * slave(some point) = master(same point) + trans(l,p) =>       *
00071  *  trans = slavecoordinates - mastercoordinates                *
00072  * uses orbits to find coordinates of center of master          *
00073  *  then solves for lin,pix of these cn. for slave.             *
00074  *                                                              *
00075  * input:                                                       *
00076  *  -                                                           *
00077  * output:                                                      *
00078  *  -                                                           *
00079  *                                                              *
00080  *    Bert Kampes, 12-Dec-1998                                  *
00081  ****************************************************************/
00082 void coarseporbit(
00083         const input_ell &ell,
00084         const slcimage  &master,
00085         const slcimage  &slave, 
00086         orbit           &masterorbit,  // cannot be const for spline
00087         orbit           &slaveorbit,   // cannot be const for spline
00088         const BASELINE  &baseline)
00089   {
00090   TRACE_FUNCTION("coarseporbit (BK 12-Dec-1998)");
00091   const int16   MAXITER   = 10;                 // maximum number of iterations
00092   const real8   CRITERPOS = 1e-6;                       // 1micrometer
00093   const real8   CRITERTIM = 1e-10;              // seconds (~10-6 m)
00094 
00095   // ______Initial values______         :master.approxcentreoriginal.x .y .z
00096   // ______Window______                 :master.currentwindow.linelo/hi , pixlo/hi
00097   // ______Time______                   :master.t_azi0/N , t_range0/N
00098 
00099   // ______Get (approx) center pixel of current window master______
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   // ______ Compute x,y,z (fill P) ______
00105   // ______ P.x/y/z contains (converged) solution ______
00106   cn P;
00107   const int32 lp2xyziter =
00108     lp2xyz(cen_lin,cen_pix,ell,master,masterorbit,P,MAXITER,CRITERPOS);
00109 
00110   // ______Compute line,pixel for slave of this xyz______
00111   real8 lin,pix;
00112   const int32 xyz2lpiter = 
00113     xyz2lp(lin,pix,slave,slaveorbit,P,MAXITER,CRITERTIM);
00114 
00115   // ______ Some extra parameters (not used, just info) ______ // BK 19-Oct-2000
00116   const int Bt          = Btemp(master.utc1,slave.utc1);
00117   // ______ Modeled quantities ______
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   // ______ Derived quantities ______
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 /*   removed: #%// Bert Kampes, 06-Apr-2005
00133 //  cn M,S;
00134 //  xyz2orb(M,master,masterorbit,P); 
00135 //  xyz2orb(S,slave,slaveorbit,P);  
00136 //  real8 B, Bperp, Bpar, Bh, Bv, alpha, theta, Hamb;
00137 //  BalphaBhBvBparBperpTheta(B,alpha,Bh,Bv,Bpar,Bperp,theta,M,P,S);
00138 //  const cn r1 = P.min(M);
00139 //  Hamb = (master.wavelength*r1.norm()*sin(theta)) / (2*Bperp);
00140 //  theta = rad2deg(theta);
00141 //  alpha = rad2deg(alpha);
00142 //  // BK 28-Nov-2000:
00143 //  real8 inc_angle = P.angle(M.min(P));
00144 //  inc_angle = rad2deg(inc_angle);
00145 //  const real8 m_tazi = master.line2ta(cen_lin);
00146 //  const real8 s_tazi = slave.line2ta(lin);
00147 //  cn Mdot = masterorbit.getxyzdot(m_tazi);
00148 //  cn Sdot = slaveorbit.getxyzdot(s_tazi);
00149 //  real8 orb_conv = Mdot.angle(Sdot);
00150 //  orb_conv = rad2deg(orb_conv);
00151 //  //const real8 heading = angle(Mdot,[1 0 0])?
00152 */
00153 
00154   // ______ offset = P_slave - P_master = lin - cen_lin ______
00155   INFO << "Estimated translation (l,p): "
00156        << floor(lin-cen_lin +.5) << ", "
00157        << floor(pix-cen_pix +.5);
00158   INFO.print();
00159 
00160   // ______ Write to tmp files ______
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)                             // round
00174                  << ", "
00175                  <<  floor(pix-cen_pix +.5)                             // round
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   // ______ give some extra info in resfile: Bperp, Bpar, Bh, Bv, Btemp ______
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)                          // round
00246     << "\nCoarse_orbits_translation_pixels: \t"
00247     <<  floor(pix-cen_pix +.5)                          // round
00248     << "\n*******************************************************************"
00249     << "\n* End_" << processcontrol[pr_i_coarse] << "_NORMAL"
00250     << "\n*******************************************************************\n";
00251 
00252   // ______Tidy up______
00253   scratchresfile.close();
00254   PROGRESS.print("Coarse precise orbits coregistration finished.");
00255   } // END coarseporbit
00256 
00257 
00258 
00259 /****************************************************************
00260  *    coarsecorrel                                              *
00261  *                                                              *
00262  * computes translation of slave w.r.t. master                  *
00263  * slave(some point) = master(same point) + trans(l,p) =>       *
00264  *  trans = slavecoordinates - mastercoordinates                *
00265  * uses correlation between magnitude of slave/master image     *
00266  *                                                              *
00267  * requires thingsa on disk, input                              *
00268  * input:                                                       *
00269  *  -                                                           *
00270  * output:                                                      *
00271  *  -                                                           *
00272  *                                                              *
00273  *    Bert Kampes, 12-Dec-1998                                  *
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];                               // for errormessages
00283   const uint Mfilelines = minfo.currentwindow.lines();
00284   const uint Sfilelines = sinfo.currentwindow.lines();
00285 
00286   const uint Nwin         = coarsecorrinput.Nwin;               // number of windows
00287   const int32 initoffsetL = coarsecorrinput.initoffsetL;        // initila offset
00288   const int32 initoffsetP = coarsecorrinput.initoffsetP;        // initila offset
00289   uint MasksizeL          = coarsecorrinput.MasksizeL;          // size of correlation window
00290   uint MasksizeP          = coarsecorrinput.MasksizeP;          // size of correlation window
00291   const uint AccL         = coarsecorrinput.AccL;               // accuracy of initial offset
00292   const uint AccP         = coarsecorrinput.AccP;               // accuracy of initial offset
00293   bool pointsrandom = true;
00294   if (specified(coarsecorrinput.ifpositions))   // filename specified
00295     pointsrandom = false;                       // only use those points
00296 
00297 
00298 //  INFO("Masksize ...
00299 
00300 // ______Only odd Masksize possible_____
00301   bool forceoddl = false;
00302   bool forceoddp = false;
00303   if (!isodd(MasksizeL)) 
00304     {
00305     forceoddl = true; 
00306     MasksizeL+=1;                       // force oddness
00307     }
00308   if (!isodd(MasksizeP))
00309     {
00310     forceoddp = true; 
00311     MasksizeP+=1;                       // force oddness
00312     }
00313 
00314 // ______Corners of slave in master system______
00315 // ______offset = A(slave system) - A(master system)______
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 // ______Corners of overlap master,slave in master system______
00322 //  const int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
00323 //  const int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
00324 //  const int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
00325 //  const int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
00326  
00327 // ______Corners of useful overlap master,slave in master system______
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 // ______Distribute Nwin points over window______
00336 // ______Centers(i,0): line, (i,1): pixel, (i,2) flagfromdisk______
00337 // old  matrix<uint> Centers = distributepoints(Nwin,overlap);
00338 
00339   register int32 i;
00340   matrix<uint> Centers;
00341   if (pointsrandom)                             // no filename specified
00342     {
00343     Centers = distributepoints(real4(Nwin),overlap);
00344     }
00345 
00346   else  // read in points (center of windows) from file
00347     {
00348     Centers.resize(Nwin,3);
00349     //ifstream ifpos(coarsecorrinput.ifpositions, ios::in);
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);                  // correct for lower left corner
00358       Centers(i,1) = uint(pp);                  // correct for lower left corner
00359       Centers(i,2) = uint(1);                   // flag from file
00360       ifpos.getline(dummyline,ONE27,'\n');              // goto next line.
00361       }
00362     ifpos.close();
00363 
00364 // ______ Check last point ivm. EOL after last position in file ______
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);      // random
00369       Centers(Nwin-1,1) = uint(.5*(pN + p0) + 37);      // random
00370       WARNING << "CC: there should be no EOL after last point in file: "
00371            << coarsecorrinput.ifpositions;
00372       WARNING.print();
00373       }
00374 
00375 // ______ Check if points are in overlap ______
00376 // ______ no check for uniqueness of points ______
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) // give some additional info
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 // ______Compute correlation of these points______
00431   //matrix<complr4> Mcmpl(2*AccL+MasksizeL,2*AccP+MasksizeP);   // on file
00432   //matrix<complr4> Scmpl(MasksizeL,MasksizeP);                         // on file
00433   matrix<complr4> Mcmpl;
00434   matrix<complr4> Scmpl;
00435   matrix<real4> Master;                         // amplitude master
00436   matrix<real4> Mask;                           // amplitude slave
00437   matrix<real4> Correl;                         // matrix with correlations
00438   matrix<real4> Result(Nwin,3);                 // R(i,0):delta l; (i,1):delta p; (i,2):correl
00439 
00440   //if (Nwin<10)
00441   //  PROGRESS("COARSE_CORR:  0%");
00442   // ______ Progress messages ______
00443   int32 percent    = 0;
00444   //int32 tenpercent = int32(floor(overlap.lines()/10)+.5);    // round
00445   int32 tenpercent = int32((Nwin/10)+.5);       // round
00446   if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
00447   for (i=0;i<Nwin;i++)
00448     {
00449     //if (Nwin>=10) if (!(i%((Nwin)/10)))               // indicate 10% progress
00450     if (i%tenpercent==0)
00451       {
00452       //int32 percentage = int32(100*(real4(i)/real4(Nwin)));
00453       PROGRESS << "COARSE_CORR: " << setw(3) << percent << "%" << ends;
00454       PROGRESS.print();
00455       percent += 10;
00456       }
00457 
00458 // ______Center of window in master system______
00459     uint cenMwinL = Centers(i,0);
00460     uint cenMwinP = Centers(i,1);
00461 
00462     window master;                                      // size=masksize+2*acc.
00463     master.linelo = cenMwinL - (MasksizeL-1)/2 -AccL;   // ML is forced odd
00464     master.linehi = master.linelo + MasksizeL +2*AccL - 1;
00465     master.pixlo  = cenMwinP - (MasksizeP-1)/2 - AccP;  // MP is forced odd
00466     master.pixhi  = master.pixlo + MasksizeP +2*AccP - 1;
00467 
00468 // ______Same points in slave system (disk)______
00469     window slavemask;                                   // size=masksize
00470     uint cenSwinL = cenMwinL + initoffsetL;             // adjust initoffset
00471     uint cenSwinP = cenMwinP + initoffsetP;             // adjust initoffset
00472     slavemask.linelo = cenSwinL - (MasksizeL-1)/2;      // ML is forced odd
00473     slavemask.linehi = slavemask.linelo + MasksizeL - 1;
00474     slavemask.pixlo  = cenSwinP - (MasksizeP-1)/2;      // MP is forced odd
00475     slavemask.pixhi  = slavemask.pixlo + MasksizeP - 1;
00476 
00477 // ______Read windows from files, compute magnitude______
00478     //minfo.readdata(Mcmpl,master);
00479     //sinfo.readdata(Scmpl,slavemask);
00480     Mcmpl = minfo.readdata(master);
00481     Scmpl = sinfo.readdata(slavemask);
00482     Master = magnitude(Mcmpl); 
00483     Mask   = magnitude(Scmpl); 
00484 
00485 
00486 // ______Compute correlation matrix and find maximum______
00487 //    DEBUG("returning matrix is not necessary, could return max,l,p"); 
00488     Correl = correlate(Master,Mask); 
00489     uint L, P;
00490     real4 corr = max(Correl, L, P);             // returns also L,P
00491 
00492     uint relcenML = master.linehi - cenMwinL;   // system of matrix
00493     uint relcenMP = master.pixhi  - cenMwinP;   // system of matrix
00494     int32 reloffsetL = relcenML - L;
00495     int32 reloffsetP = relcenMP - P;
00496     int32 offsetL = reloffsetL + initoffsetL;   // estimated offset lines
00497     int32 offsetP = reloffsetP + initoffsetP;   // estimated offset pixels
00498 
00499     Result(i,0) = offsetL;
00500     Result(i,1) = offsetP;
00501     Result(i,2) = corr;
00502     }
00503 
00504 // ______Get correct offsetL, offsetP______
00505   int32 offsetLines = -999;
00506   int32 offsetPixels = -999;
00507   getoffset(Result,offsetLines,offsetPixels);
00508 
00509 
00510 // ______Write to files______
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                                // 1 digit after point?
00545                  << "\nCoarse_correlation_translation_pixels: \t"
00546                  <<  offsetPixels                               // 1 digit after point?
00547                  << "\n*******************************************************************"
00548                  //<< "\n* End_coarse_correlation:_NORMAL"
00549                  << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
00550                  << "\n*******************************************************************\n";
00551   scratchresfile.close();
00552 
00553 // ______Tidy up______
00554   INFO << "Individually estimated translations (#, l, p, corr, offl, offp): ";
00555   INFO.print();
00556   for (k=0; k<Nwin; k++)
00557     {
00558     // "Estimate (#, l, p, corr, offl, offp): "
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   } // END coarsecorrel
00572 
00573 
00574 
00575 /****************************************************************
00576  *    coarsecorrelfft                                           *
00577  *                                                              *
00578  * computes translation of slave w.r.t. master                  *
00579  * slave(some point) = master(same point) + trans(l,p) =>       *
00580  *  trans = slavecoordinates - mastercoordinates                *
00581  * uses correlation between magnitude of slave/master image     *
00582  * uses fft to compute coherence, no subtraction of mean        *
00583  *                                                              *
00584  * requires thingsa on disk, input                              *
00585  * input:                                                       *
00586  *  -                                                           *
00587  * output:                                                      *
00588  *  -                                                           *
00589  *                                                              *
00590  *    Bert Kampes, 12-Dec-1998                                  *
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];                       // for errormessages
00605   const uint Mfilelines   = minfo.currentwindow.lines();
00606   const uint Sfilelines   = sinfo.currentwindow.lines();
00607 
00608   const uint Nwin           = coarsecorrinput.Nwin;             // number of windows
00609   const int32 initoffsetL = coarsecorrinput.initoffsetL;        // initila offset
00610   const int32 initoffsetP = coarsecorrinput.initoffsetP;        // initila offset
00611   const uint MasksizeL    = coarsecorrinput.MasksizeL;  // size of correlation window
00612   const uint MasksizeP    = coarsecorrinput.MasksizeP;  // size of correlation window
00613 
00614   bool pointsrandom = true;
00615   if (specified(coarsecorrinput.ifpositions))   // filename specified
00616     pointsrandom = false;                       // only use these points
00617 
00618   // ______Only pow2 Masksize possible_____
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   // ______Corners of slave in master system______
00631   // ______offset = [A](slave system) - [A](master system)______
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 // ______Corners of overlap master,slave in master system______
00638 //  int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
00639 //  int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
00640 //  int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
00641 //  int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
00642  
00643   // ______Corners of useful overlap master,slave in master system______
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   // ______Distribute Nwin points over window______
00652   // ______Minlminp(i,0): line, (i,1): pixel, (i,2) flagfromdisk______
00653   // old  matrix<uint> Minlminp = distributepoints(Nwin,overlap);
00654   register int32 i;
00655 
00656   matrix<uint> Minlminp;
00657   if (pointsrandom)                             // no filename specified
00658     {
00659     Minlminp = distributepoints(real4(Nwin),overlap);
00660     }
00661 
00662   else  // read in points (center of windows) from file
00663     {
00664     Minlminp.resize(Nwin,3);
00665     //ifstream ifpos(coarsecorrinput.ifpositions, ios::in);
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);   // correct for lower left corner
00674       Minlminp(i,1) = uint(pp -.5*MasksizeP);   // correct for lower left corner
00675       Minlminp(i,2) = uint(1);                  // flag from file
00676       ifpos.getline(dummyline,ONE27,'\n');      // goto next line.
00677       }
00678     ifpos.close();
00679 
00680     // ______ Check last point ivm. EOL after last position in file ______
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);     // random
00685       Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);     // random
00686       }
00687 
00688     // ______ Check if points are in overlap ______
00689     // ______ no check for uniqueness of points ______
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) // give some additional info
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   // ______Compute coherence of these points______
00749   //matrix<complr4> Master(MasksizeL,MasksizeP);        // on file
00750   //matrix<complr4> Mask(MasksizeL,MasksizeP);  // on file
00751   matrix<complr4> Master;
00752   matrix<complr4> Mask;
00753   matrix<real4>   Result(Nwin,3);               // R(i,0):delta l; 
00754                                                 //  R(i,1):delta p; R(i,2):correl
00755   // ______ Progress messages ______
00756   int32 percent    = 0;
00757   //int32 tenpercent = int32((overlap.lines()/10)+.5);    // round
00758   int32 tenpercent = int32((Nwin/10)+.5);       // round
00759   if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
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     // ______Minlminp (lower left corners) of window in master system______
00770     uint minMwinL = Minlminp(i,0);
00771     uint minMwinP = Minlminp(i,1);
00772                                                         // size=masksize
00773     window master(minMwinL, minMwinL+MasksizeL-1,
00774                   minMwinP, minMwinP+MasksizeP-1);
00775     // ______Same points in slave system (disk)______
00776     window mask(minMwinL+initoffsetL, 
00777                 minMwinL+initoffsetL+MasksizeL-1,
00778                 minMwinP+initoffsetP,
00779                 minMwinP+initoffsetP+MasksizeP-1);
00780 
00781     // ______Read windows from files______
00782     //minfo.readdata(Master,master);
00783     //sinfo.readdata(Mask,mask);
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     // ______Coherence______
00791     // ______update offsetL/P______
00792     real4 offsetL, offsetP;
00793     real4 coheren;
00794 
00795     coheren = corrfft(absMaster,absMask,offsetL,offsetP);
00796 
00797     offsetL += initoffsetL;                     // estimated offset lines
00798     offsetP += initoffsetP;                     // estimated offset pixels
00799 
00800     Result(i,0) = offsetL;
00801     Result(i,1) = offsetP;
00802     Result(i,2) = coheren;
00803     } // for nwin
00804 
00805   // ______Position approx. with respect to center of window______
00806   // ______correct position array for center instead of lower left______
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   // ______Get correct offsetL, offsetP______
00814   int32 offsetLines = -999;
00815   int32 offsetPixels = -999;
00816   getoffset(Result,offsetLines,offsetPixels);
00817 
00818 
00819   // ______Write to files______
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 //               << "\nSearchwindow size (l,p):       \t"
00830 //               <<  MasksizeL + 2*AccL << ", " << MasksizeP + 2*AccP
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                        // 1 digit after point?
00852                  << "\nCoarse_correlation_translation_pixels: \t"
00853                  <<  offsetPixels                       // 1 digit after point?
00854                  << "\n*******************************************************************"
00855                  //<< "\n* End_coarse_correlation:_NORMAL"
00856                  << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
00857                  << "\n*******************************************************************\n";
00858 
00859 // ______Tidy up______
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 // << "Estimate (#, l, p, corr, offl, offp): "
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   } // END coarsecorrelfft
00880 
00881 
00882 
00883 /****************************************************************
00884  * corrfft                                                      *
00885  *                                                              *
00886  * coherence in spectral domain by fft's                        *
00887  *  uses extension with zeros                                   *
00888  *  pixel level                                                 *
00889  *                                                              *
00890  * input:                                                       *
00891  *  - Master                                                    *
00892  *  - Mask (size Master)                                        *
00893  * output:                                                      *
00894  *  - coherence value                                           *
00895  *  - updated offsetL, P                                        *
00896  *                                                              *
00897  *    Bert Kampes, 16-Feb-1999                                  *
00898  *                                                              *
00899  * note: this routine can be speeded up by removing matrices    *
00900  * powerma* and by using *= instead of dotmult                  *
00901  * for now this is not done because it requires only little time*
00902  * note also that for coarse coregistration division by powers  *
00903  * is not really required, cross products are good enough.      *
00904  *    Bert Kampes, 18-Oct-1999                                  *
00905  ****************************************************************/
00906 real4 corrfft(
00907          const matrix<real4> &Master,                   // magnitude image
00908          const matrix<real4> &Mask,                     // magnitude image
00909          real4 &offsetL,                                // updated
00910          real4 &offsetP)                                // updated 
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   // ======Compute powers for submatrices======
00932   register int32 i;
00933   register int32 j;
00934 
00935   const complr4 ONE(1.);
00936   matrix<complr4> Master2(twoL,twoP);                   // init 0
00937   matrix<complr4> Mask2(twoL,twoP);                     // init 0
00938   matrix<complr4> blok2(twoL,twoP);                     // init 0
00939 
00940   const real4 meanMaster = mean(Master);
00941   const real4 meanMask   = mean(Mask);
00942 
00943 // ====== Powers, misuse Master2, blok2 ======
00944 // ______ First powers to save a matrix, 3 is minimum ______
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));    // intensity image
00950       Mask2(i+L,j+P) = complr4(sqr(Mask(i,j)-meanMask));// only real part
00951       blok2(i+halfL,j+halfP) = ONE;                     // only real part
00952       }
00953     }
00954 
00955   fft2d(Master2);
00956   fft2d(Mask2);
00957   fft2d(blok2);
00958 
00959 // old way:
00960 //  matrix<complr4> Powermaster = dotmult(conj(Master2),blok2); // twice conj(M2)
00961 //  ifft2d(Powermaster);
00962 //  matrix<complr4> Powermask = dotmult(conj(blok2),Mask2);
00963 //  ifft2d(Powermask);
00964 
00965 // new way: test this
00966   blok2.conj();                                         // conjugated in blok2
00967   Master2 *= blok2;
00968   Master2.conj();                                       // M2=original(b2)*conj(m2)
00969   ifft2d(Master2);                                      // norms master in Master2
00970 
00971   Mask2 *= blok2;                                       // powers in spectral domain
00972   ifft2d(Mask2);                                        // norms slave in Mask2
00973 
00974 // ______ Use real(block2) to store sqrt(norms1*norms2) ______
00975   for (i=0; i<=L; i++)                          // all shifts
00976     for (j=0; j<=P; j++)                        // all shifts
00977       blok2(i,j) = complr4(sqrt(real(Master2(i,j))*real(Mask2(i,j))));
00978 
00979 
00980 // ====== Cross products covariance master/slave ======
00981   Master2.clean();                                      // init 0
00982   Mask2.clean();                                        // init 0
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); // only real
00988       Mask2(i+halfL,j+halfP) = complr4(Mask(i,j)-meanMask);     // part mag. image
00989       }
00990     }
00991 
00992 
00993 // ======FFT's of master/mask======
00994 // padded with N zeros to prevent periodical convolution
00995   fft2d(Master2);
00996   fft2d(Mask2);
00997 
00998 // ______ Store in Mask2 crossproducts in spectral/space domain ______
00999   Master2.conj();
01000   Mask2 *= Master2;                                     // corr. by zero padding
01001   ifft2d(Mask2);                                        // space domain (real only)
01002 //  matrix<complr4> Cor = dotmult(conj(Master2),Mask2); // correlation with zero ext.
01003 //  ifft2d(Cor);                                        // space
01004 
01005 
01006 // ====== Correlation in space domain for all shifts [-N/2,N/2] ======
01007   real4 coher;
01008   real4 maxcoher = 0.;
01009   for (i=0; i<=L; i++)                          // all shifts
01010     {
01011     for (j=0; j<=P; j++)                        // all shifts
01012       {
01013       //coher = (sqr(real(Cor(i,j)))) /                         // old way
01014 //            ((real(Powermaster(i,j))*real(Powermask(i,j))));
01015 //      coher = (sqr(real(Cor(i,j)))) /                                 // new way
01016 //            ((real(Master2(i,j))*real(Mask2(i,j))));
01017       coher = real(Mask2(i,j)) / real(blok2(i,j));
01018       if (coher > maxcoher)
01019         {
01020         maxcoher = coher;
01021         offsetL = -halfL + i;                   // update by reference
01022         offsetP = -halfP + j;                   // update by reference
01023         }
01024       }
01025     }
01026   return maxcoher;
01027   } // END corrfft
01028 
01029 
01030 
01031 /****************************************************************
01032  *    distributepoints                                          *
01033  *                                                              *
01034  * Returns matrix with distributed points in of input.          *
01035  * First window at (win.linelo,win.pixlo),                      *
01036  *  divided over wl lines,                                      *
01037  *  with dp distance in pixel direction.                        *
01038  *                                                              *
01039  * input:                                                       *
01040  *  - number of windows                                         *
01041  *  - window which should be divided                            *
01042  * output:                                                      *
01043  *  - matrix <uint> (NW,3) =(l,p, flagfromdisk==0)              *
01044  *                                                              *
01045  *    Bert Kampes, 21-Jan-1999                                  *
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   // ______Distribution for dl=dp______
01058   real4 wp = sqrt(nW/(lines/pixels));   // wl: #windows in line direction
01059   real4 wl = nW / wp;                   // wp: #windows in pixel direction
01060 
01061   if (wl < wp)                          // switch wl,wp : later back
01062     wl = wp; 
01063 
01064   int32 wlint  = int32(floor(wl + .5)); // round largest
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)         // round
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 // ______Correct distribution to window______
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   } // END distributepoints
01097 
01098 
01099 
01100 /****************************************************************
01101  *    getoffset                                                 *
01102  *                                                              *
01103  * Returns offset in line and pixel direction                   *
01104  *  based on matrix with estimated offests                      *
01105  *  by correlation                                              *
01106  * Checks on consistency, THRESHOLD 0.4 for correlation         *
01107  *                                                              *
01108  * input:                                                       *
01109  *  - matrix<real4> with offsets(l,p),correlation               *
01110  *  - (offL,offP)                                               *
01111  * output:                                                      *
01112  *  - updated (offL,offP)                                       *
01113  *                                                              *
01114  * See also Documentation page 6.                               *
01115  *                                                              *
01116  *    Bert Kampes, 21-Jan-1999                                  *
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;                    // Return these
01154       offsetPixels = valueP;                    // Return these
01155       }
01156     }
01157 
01158   // ______ Check result ______
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   } // END getoffset
01178 
01179 
01180 
01181 /****************************************************************
01182  *    finecoreg                                                 *
01183  *                                                              *
01184  * computes translation of slave w.r.t. master                  *
01185  * slave(some point) = master(same point) + trans(l,p) =>       *
01186  *  trans = slavecoordinates - mastercoordinates                *
01187  * in NWIN windows.                                             *
01188  * Then solves polynomial for best transformation to master     *
01189  * with coregpm routine/step                                    *
01190  *                                                              *
01191  * input:                                                       *
01192  *  -                                                           *
01193  * output:                                                      *
01194  *  -                                                           *
01195  *    Bert Kampes, 12-Dec-1998                                  *
01196  * distribute points can be with input file as well, besides    *
01197  * letting Doris randomly distribute npoints.                   *
01198  *    BK 29-Oct-99                                              *
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;           // n windows, from file or random
01212   const int32 initoffsetL = fineinput.initoffsetL;      // initial offset
01213   const int32 initoffsetP = fineinput.initoffsetP;      // initial offset
01214   uint MasksizeL    = fineinput.MasksizeL;      // size of correlation window
01215   uint MasksizeP    = fineinput.MasksizeP;      // size of correlation window
01216 
01217   bool pointsrandom = true;
01218   if (specified(fineinput.ifpositions))         // filename specified
01219     pointsrandom = false;                       // only use these points
01220 
01221 // ______Correct sizes if in space domain______
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 // ______Corners of slave in master system______
01231 // ______offset = [A](slave system) - [A](master system)______
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 // ______Corners of overlap master,slave in master system______
01238 //  int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
01239 //  int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
01240 //  int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
01241 //  int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
01242  
01243 // ______Corners of useful overlap master,slave in master system______
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 // ______ Distribute Nwin points over window, or read from file ______
01252 // ______ Minlminp(i,0): line, (i,1): pixel, (i,2) flagfromdisk ______
01253   register int32 i;
01254   matrix<uint> Minlminp;
01255   if (pointsrandom)                             // no filename specified
01256     {
01257     Minlminp = distributepoints(real4(Nwin),overlap);
01258     }
01259 
01260   else  // read in points (center of windows) from file
01261     {
01262     // const int32 numwin = filelines(fineinput.ifpositions); // done in readinput
01263     Minlminp.resize(Nwin,3);
01264     //ifstream ifpos(fineinput.ifpositions, ios::in | ios::nocreate);
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);   // correct for lower left corner
01272       Minlminp(i,1) = uint(pp -.5*MasksizeP);   // correct for lower left corner
01273       Minlminp(i,2) = uint(1);                  // flag from file
01274       ifpos.getline(dummyline,ONE27,'\n');      // goto next line.
01275       }
01276     ifpos.close();
01277 
01278 // ______ Check last point ivm. EOL after last position in file ______
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);     // random
01283       Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);     // random
01284       }
01285 
01286 // ______ Check if points are in overlap ______
01287 // ______ no check for uniqueness of points ______
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) // give some additional info
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   // ______Compute coherence of these points______
01346   //matrix<complr4> Master(MasksizeL,MasksizeP);        // on file
01347   //matrix<complr4> Mask(MasksizeL,MasksizeP);  // on file
01348   matrix<complr4> Master;
01349   matrix<complr4> Mask;
01350   matrix<real4>   Result(Nwin,3);               // R(i,0):delta l; 
01351                                                 //  R(i,1):delta p; R(i,2):correl
01352 
01353   // ______ Progress message ______
01354   int32 tenpercent = int32((Nwin/10)+.5);    // round
01355   if (tenpercent==0) tenpercent = 1000;
01356   int32 percent = 0;
01357 
01358   // ====== Compute for all locations ======
01359   for (i=0;i<Nwin;i++)
01360     {
01361     // ______ Give progress message ______
01362     if (i%tenpercent==0)
01363       {
01364       PROGRESS << "FINE: " << setw(3) << percent << "%";
01365       PROGRESS.print();
01366       percent += 10;
01367       }
01368 
01369     // ______Minlminp (lower left corners) of window in master system______
01370     uint minMwinL = Minlminp(i,0);
01371     uint minMwinP = Minlminp(i,1);
01372 
01373                                                         // size=masksize
01374     window master(minMwinL, minMwinL+MasksizeL-1,
01375                   minMwinP, minMwinP+MasksizeP-1);
01376     // ______Same points in slave system (disk)______
01377     window mask(minMwinL+initoffsetL, 
01378                 minMwinL+initoffsetL+MasksizeL-1,
01379                 minMwinP+initoffsetP,
01380                 minMwinP+initoffsetP+MasksizeP-1);
01381 
01382     // ______Read windows from files______
01383     //minfo.readdata(Master,master);
01384     //sinfo.readdata(Mask,mask);
01385     Master = minfo.readdata(master);
01386     Mask   = sinfo.readdata(mask);
01387 
01388 
01389     // ______Coherence______
01390     // ______update offsetL/P______
01391     real4 offsetL, offsetP;
01392     real4 coheren;
01393     switch (fineinput.method)
01394       {
01395 
01396 //        case fc_cmplxfft:
01397 //      WARNING("THIS METHOD IS NOT OK YET, I RECOMMEND MAGNITUDE.");
01398 //          coheren = coherencefft(fineinput, Master, Mask, offsetL, offsetP);
01399 //      break;
01400 //        case fc_cmplxspace:
01401 //      WARNING("THIS METHOD IS NOT OK YET, I RECOMMEND MAGNITUDE.");
01402 //          coheren = coherencespace(fineinput, Master, Mask, offsetL, offsetP);
01403 //      break;
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       } // switch method
01415 
01416 
01417     Result(i,0) = offsetL + initoffsetL;
01418     Result(i,1) = offsetP + initoffsetP;
01419     Result(i,2) = coheren;
01420     } // for nwin
01421 
01422 // ______Position approx. with respect to center of window______
01423 // ______correct position array for center instead of lower left______
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 // ______Write to files______
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 // ______Tidy up______
01493   fclose(resfile);
01494   PROGRESS.print("Fine coregistration finished.");
01495   } // END finecoreg
01496 
01497 
01498 
01499 
01500 /****************************************************************
01501  * coherencefft                                                 *
01502  *                                                              *
01503  * coherence in spectral domain by fft's based on magnitude     *
01504  *  uses extension with zeros                                   *
01505  *                                                              *
01506  * input:                                                       *
01507  *  - Master                                                    *
01508  *  - Mask (size Master)                                        *
01509  * output:                                                      *
01510  *  - coherence value [-1 1]                                    *
01511  *  - updated offsetL, P                                        *
01512  *                                                              *
01513  *    Bert Kampes, 03-Feb-1999                                  *
01514  * bugfix? streamlined, based on magnitude forced               *
01515  *    Bert Kampes, 16-Nov-1999                                  *
01516  ****************************************************************/
01517 real4 coherencefft(
01518         const input_fine &fineinput,
01519         const matrix<complr4> &Master,          // data
01520         const matrix<complr4> &Mask,            // data
01521         real4 &offsetL,                         // returned
01522         real4 &offsetP)                         // returned
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 // ______ Zero mean magnitude images ______
01557   matrix<real4> magMaster = magnitude(Master);
01558   matrix<real4> magMask   = magnitude(Mask);
01559   magMaster              -= mean(magMaster);
01560   magMask                -= mean(magMask);
01561 
01562 
01563 // ====== FFT's of master/mask ======
01564 // ______ Pad with N zeros to prevent periodical convolution ______
01565   window windef(0,0,0,0);                       // defaults to total matrix
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);                   // initial 0
01570   matrix<complr4> Mask2(twoL,twoP);                     // initial 0
01571   Master2.setdata(win1,mat2cr4(magMaster),windef);      // zero mean magnitude
01572   Mask2.setdata(win2,mat2cr4(magMask),windef);          // zero mean magnitude
01573 
01574   fft2d(Master2);
01575   fft2d(Mask2);
01576 
01577 // ______ Crossproducts in spectral/space domain ______
01578 // ______ Use Mask2 to store cross products temporarly ______
01579   Master2.conj();
01580   Mask2 *= Master2;                                     // corr = conj(M).*S
01581   ifft2d(Mask2);                                        // cross prod. in space
01582 
01583 // ______ Resize to shifts [-AccL,+AccL) ______
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);              // imag == 0 (?)
01587 
01588 
01589 // ====== Compute norms, zero padded matrices ======
01590   matrix<complr4> blok(L,P);
01591   const complr4 ONE(1.);
01592   blok.setdata(ONE);                                    // only real part
01593 
01594   Master2.clean();                                      // reset to zeros
01595   Mask2.clean();                                        // reset to zeros
01596   Master2.setdata(win1,mat2cr4(sqr(magMaster)),windef);// use Master2 for intensity
01597   Mask2.setdata(win2,blok,windef);                      // use Mask2 for padded Block
01598 
01599   fft2d(Master2);                                       // (intensity of master)
01600   fft2d(Mask2);                                         // (block)
01601 
01602   Mask2.conj();                                         // conj(block)
01603   Master2 *= Mask2;
01604   Master2.conj();                               // Master2 == conj(Master)*block
01605   ifft2d(Master2);
01606 // ______ Master2 now contains norms of master image in space domain ______
01607 // ______ Resize to shifts [-AccL,+AccL) ______
01608   TMP.setdata(Master2,wintmp);                          // fill TMP
01609   matrix<real4> pmaster = real(TMP);                    // norms in pmaster
01610 
01611 // ====== Now norms for slave image ======
01612   Master2.clean();                                      // reset to zeros
01613   window win5(L,twoL-1,P,twoP-1);
01614   Master2.setdata(win5,mat2cr4(sqr(magMask)),windef);
01615   fft2d(Master2);                                       // (intensity of slave)
01616   Master2 *= Mask2;                             // Master2 == conj(block)*Slave
01617   ifft2d(Master2);
01618 // ______ Master2 now contains norms of slave image in space domain ______
01619 // ______Resize to shifts [-AccL,+AccL)______
01620   TMP.setdata(Master2,wintmp);                          // fill TMP
01621   matrix<real4> pmask = real(TMP);                      // norms in pmask
01622 
01623 // ====== Correlation in space domain ======
01624   pmaster *= pmask;
01625   Covar /= sqrt(pmaster);
01626 
01627 
01628 // ====== Estimate shift by oversampling ======
01629   uint offL;
01630   uint offP;
01631   // old: matrix<real4> coher8 = oversample(Covar,factor);
01632   matrix<real4> coher8 = oversample(Covar,factor,factor);
01633   real4 maxcor = max(coher8,offL,offP);
01634   offsetL = -real4(AccL) + offL/real4(factor);           // update by reference
01635   offsetP = -real4(AccP) + offP/real4(factor);           // update by reference
01636 
01637   return maxcor;
01638   } // END coherence fft
01639 
01640 
01641 
01642 /****************************************************************
01643  * coherencespace                                               *
01644  *                                                              *
01645  * coherence in space domain based on magnitude                 *
01646  *  uses extension with zeros                                   *
01647  *                                                              *
01648  * input:                                                       *
01649  *  - Master                                                    *
01650  *  - Mask (size Master)                                        *
01651  * output:                                                      *
01652  *  - coherence value                                           *
01653  *  - updated offsetL, P                                        *
01654  *                                                              *
01655  *    Bert Kampes, 03-Feb-1999                                  *
01656  ****************************************************************/
01657 real4 coherencespace(
01658         const input_fine &fineinput,
01659         const matrix<complr4> &Master,          // complex data
01660         const matrix<complr4> &Mask,            // complex data
01661         real4 &offsetL,                         // returned
01662         real4 &offsetP)                         // returned
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   // ______Select parts of Master/slave______
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   // ______Shift center of Slave over Master______
01688   window winmask(AccL, AccL+MasksizeL-1,
01689                  AccP, AccP+MasksizeP-1);
01690 
01691   matrix<real4> coher(2*AccL,2*AccP);           // store result
01692                                                 // 1st element: shift==AccL
01693   window windef(0, 0, 0, 0);                    // defaults to total
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 //    if (fineinput.method == fc_cmplxspace)
01705 //      {
01706 //  DEBUG("mean should be subtracted.        e");
01707 //      matrix<complr4> Mask2(winmask,Mask);            // construct as part 
01708 //      real4 normmask = norm2(Mask2);
01709 //    
01710 //      matrix<complr4> Master2(MasksizeL, MasksizeP); 
01711 //      window winmaster;
01712 //      
01713 //      for (register int32 i=0;i<2*AccL;i++)
01714 //        {
01715 //        winmaster.linelo = i;
01716 //        winmaster.linehi = i+MasksizeL-1;
01717 //        for (register int32 j=0;j<2*AccP;j++)
01718 //          {
01719 //          winmaster.pixlo = j;
01720 //          winmaster.pixhi = j+MasksizeP-1;
01721 //          Master2.setdata(windef,Master,winmaster);
01722 //       
01723 //  // ______Compute coherence for this position______
01724 //          complr4 cohs1s2(0.,0.);
01725 //          real4 cohs1s1 = 0.;
01726 //  DEBUG("nog multilooked coherence-ref.phase");
01727 //          for (register int32 k=0;k<MasksizeL;k++)
01728 //            {
01729 //            for (register int32 l=0;l<MasksizeP;l++)
01730 //              {
01731 //          cohs1s2 += (Master2(k,l) * conj(Mask2(k,l)));
01732 //          cohs1s1 += norm(Master2(k,l));
01733 //              }
01734 //            }
01735 //          coher(i,j) = norm(cohs1s2) / (cohs1s1 * normmask);
01736 //          }
01737 //        }
01738 //      }
01739 
01740     case fc_magspace:
01741       {
01742       matrix<real4> magMask = magnitude(Mask);          // magnitude
01743       magMask -= mean(magMask);                         // subtract mean
01744       //matrix<real4> Mask2(winmask,magnitude(Mask));   // construct as part 
01745       matrix<real4> Mask2(winmask,magMask);             // construct as part 
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 // ______Coherence for this position______
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);      // [-1 1]
01777           //coher(i,j) = sqr(cohs1s2) / (cohs1s1 * normmask);
01778           }
01779         }
01780       break;
01781       }
01782    
01783     default:
01784       PRINT_ERROR("unknown method")
01785       throw(unhandled_case_error);
01786     } // switch method
01787 
01788 
01789 // ______ Correlation in space domain ______
01790   // old matrix<real4> coher8 = oversample(coher,factor);
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);           // update by reference
01797   offsetP = AccP - offP/real4(factor);           // update by reference
01798 
01799   return maxcor;
01800   } // END coherencespace
01801 
01802 
01803 
01804 /****************************************************************
01805  * coregpm                                                      *
01806  *                                                              *
01807  * Compute coregistration parameters (least squares)            *
01808  *                                                              *
01809  * input:                                                       *
01810  *  - Position of windows                                       *
01811  *  - Computed offsets                                          *
01812  *                                                              *
01813  * output:                                                      *
01814  *  - coregistration parameters to file                         *
01815  *    (wrt. normalized master grid)                             *
01816  *                                                              *
01817  *    Bert Kampes, 22-Feb-1999                                  *
01818  *    Bert Kampes, 26-Oct-1999 normalized coordinates           *
01819  * changed matrices from real4 to real8,                        *
01820  #%// BK 22-Mar-2001                                            *
01821  ****************************************************************/
01822 void coregpm(
01823         const window &originalmaster,           // normalization factors
01824         const char* i_resfile,
01825         const input_coregpm &coregpminput,
01826         uint oversamplingsfactorfine)
01827   {
01828   TRACE_FUNCTION("coregpm (BK 26-Oct-1999)")
01829   // ______Names of variables in this routine______
01830   // unknowns:          x
01831   // solution unknowns: placed in rhsL/P
01832   // observations:      y
01833   // covariance obs.:   Qy_1
01834   // designmatrix:      A
01835   // normalmatrix:      N = At. Qy-1. A
01836   // covariance unkn.:  N_1 (Qx_hat, inverse normalmatrix)
01837   // estimates:         *_hat
01838 
01839   const real4 THRESHOLD = coregpminput.threshold;       // threshold ...
01840   const int32 DEGREE    = coregpminput.degree;          // degree of polynomial
01841   const int32 MAX_ITERATIONS = coregpminput.maxiter;    // max. of pnts to remove
01842   const real4 CRIT_VALUE = coregpminput.k_alpha;        // crit. value outlier removal
01843   const int32 Nunk      = Ncoeffs(DEGREE);              // Number of unknowns/direction
01844 
01845   // ______ Normalize data for polynomial ______
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   // ______ A priori sigma of  offset ______
01852   // _____ oversampling factor is bin in which maximum can be found _____
01853   // _____ ovsf=16-->apriorisigma=0.03
01854   //const real4 SIGMA   = .5 * (1. / (real4(oversamplingsfactorfine)));
01855   const real4 ACCURACY = .5 * (1. / (real4(oversamplingsfactorfine)));
01856   // but we need coreg accuracy of 0.1 pixel about.  therefore use a priori
01857   // based on experience here, and different for azimuth and range
01858   // this also helps our automated outlier detection and testing hopefully.
01859   // BK 15-Apr-2003
01860   //const real4 SIGMAL   = 0.075;// hmmm. ?? sigma in pixels
01861   const real4 SIGMAL   = 0.15;// hmmm. ?? sigma in pixels
01862   const real4 SIGMAP   = SIGMAL/1.5;// seems pixel is 2x better ???
01863   //const real4 SIGMA_1L = 1./SIGMAL;
01864   //const real4 SIGMA_1P = 1./SIGMAP;
01865   //const real4 SIGMA2L  = sqr(SIGMAL);
01866   //const real4 SIGMA2P  = sqr(SIGMAP);
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   // ______ Find #points > threshold ______
01874   matrix<real4> Data = getofffile(i_resfile, THRESHOLD);
01875   // ______Data contains:______
01876   // Data(i,0) = winnumber; Data(i,1) = posL; Data(i,2) = posP; 
01877   // Data(i,3) = offL;      Data(i,4) = offP; Data(i,5) = corr;
01878 
01879 
01880   int32 ITERATION  = 0;
01881   int32 DONE       = 0;
01882   // sqr: level significance: alpha=0.001; power of test: gamma=0.80
01883   //real4 CRIT_VALUE = sqrt(3.29);
01884   INFO << "Critical value for outlier test: " << CRIT_VALUE;
01885   INFO.print();
01886   uint winL = 0;// window number to be removed
01887   uint winP = 0;// window number of largest w -test in range
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     // ______ Remove identified outlier from previous estimation ______
01906     if (ITERATION != 0)
01907       {
01908       matrix<real4> tmp_DATA = Data; //(remove_observation_i,*);
01909       Data.resize(Data.lines()-1, Data.pixels());
01910       j = 0;// counter over reduced obs.vector
01911       for (i=0; i<tmp_DATA.lines(); i++)// counter over original window numbers
01912         {
01913         if (i != winL)// do not copy the one to be removed.
01914           {
01915           Data.setrow(j,tmp_DATA.getrow(i));// copy back without removed obs.
01916           j++;// fill next row of Data
01917           }
01918         else
01919           {
01920           INFO << "Removing observation " << i << " from observation vector.";
01921           INFO.print();
01922           }
01923         }
01924       }
01925 
01926     // ______Check redundancy______
01927     int32 Nobs = Data.lines();                          // Number of points > threshold
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     // ______Set up system of equations______
01935     // ______Order unknowns: A00 A10 A01 A20 A11 A02 A30 A21 A12 A03 for degree=3______
01936     matrix<real8> yL(Nobs,1);                   // observation
01937     matrix<real8> yP(Nobs,1);                   // observation
01938     matrix<real8> A(Nobs,Nunk);                 // designmatrix
01939     matrix<real8> Qy_1(Nobs,1);                 // a priori covariance matrix (diag)
01940   
01941     // ______ Normalize data for polynomial ______
01942     INFO << "coregpm: polynomial normalized by factors: "
01943          << minL << " " << maxL << " " << minP << " " << maxP << " to [-2,2]";
01944     INFO.print();
01945   
01946     // ______Fill matrices______
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       // ______Set up designmatrix______
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     // ______Weight matrix data______
01971     DEBUG.print("Setting up covariance matrix for LS adjustment");
01972     switch(coregpminput.weightflag)
01973       {
01974       case 0:
01975         //Qy_1 = 1.0;
01976         for (i=0; i<Nobs; i++)
01977           Qy_1(i,0)     = real8(1.0); 
01978         break;
01979       case 1:
01980         //Qy_1 = Data(*,5)
01981         DEBUG.print("Using sqrt(coherence) as weights.");
01982         for (i=0; i<Nobs; i++)
01983           Qy_1(i,0)       = real8(Data(i,5));     // more weight to higher correlation
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));// more weight to higher correlation
01989         break;
01990       default:
01991         PRINT_ERROR("Panic, NOT possible with checked input.")
01992         throw(unhandled_case_error);
01993       }
01994     // ______ Normalize weights to avoid influence on estimated var.factor ______
01995     INFO.print("Normalizing covariance matrix for LS estimation.");
01996     Qy_1 = Qy_1 / mean(Qy_1);// normalize weights (for tests!)
01997 
01998  
01999     // ______Compute Normalmatrix, rghthandside______
02000     matrix<real8> N    = matTxmat(A,diagxmat(Qy_1,A));
02001     //matrix<real8> rhsL = matTxmat(A,diagxmat(Qy_1,yL));
02002     //matrix<real8> rhsP = matTxmat(A,diagxmat(Qy_1,yP));
02003     //matrix<real8> Qx_hat = N;
02004     rhsL = matTxmat(A,diagxmat(Qy_1,yL));
02005     rhsP = matTxmat(A,diagxmat(Qy_1,yP));
02006     Qx_hat = N;
02007     // ______Compute solution______
02008     choles(Qx_hat);             // Cholesky factorisation normalmatrix
02009     solvechol(Qx_hat,rhsL);     // Solution unknowns in rhs
02010     solvechol(Qx_hat,rhsP);     // Solution unknowns in rhs
02011     invertchol(Qx_hat);         // Covariance matrix of unknowns
02012     // ______Test inverse______
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);// repair Qx
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     matrix<real8> unity = N * Qx_hat;
02021     maxdev = abs(unity(0,0)-1);
02022     for (i=1; i<unity.lines(); i++)
02023       {
02024       if (abs(unity(i,i)-1) > maxdev)
02025         maxdev = abs(unity(i,i)-1);
02026       for (j=0; j<i-1; j++)
02027         {
02028         if (abs(unity(i,j)) > maxdev)
02029           maxdev = abs(unity(i,j));
02030         if (abs(unity(j,i)) > maxdev)
02031           maxdev = abs(unity(j,i));
02032         }
02033       }
02034     */
02035     // ___ use trace buffer to store string, remember to rewind it ___
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     // ______Some other stuff, scale is ok______
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     //matrix<real8> eL_hat      = yL - yL_hat;
02056     //matrix<real8> eP_hat      = yP - yP_hat;
02057     eL_hat      = yL - yL_hat;
02058     eP_hat      = yP - yP_hat;
02059     //  matrix<real4> Qe_hat    = Qy - Qy_hat;
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     // ______Overall model test (variance factor)______
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);// this is sigma hat!
02074     overallmodeltestP = (overallmodeltestP/sqr(SIGMAP)) /(Nobs-Nunk);// not OMT!
02075   
02076 
02077     INFO << "coregpm: overallmodeltest Lines = " << overallmodeltestL;
02078     INFO.print();
02079     INFO << "coregpm: overallmodeltest Pixels = " << overallmodeltestP;
02080     INFO.print();
02081   
02082   
02083     // ______Datasnooping, assume Qy diag______
02084     //matrix<real8> wtestL(Nobs,1);
02085     //matrix<real8> wtestP(Nobs,1);
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);// computed excl.var.factor
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);     // returns winL
02097     maxwP     = max(abs(wtestP),winP,dumm);     // returns winP
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     //if (maxwP>maxwL)
02107     //  {
02108     //  winL  = winP;// winL points to rejected window number
02109     //  //maxwL = maxwP;
02110     //  }
02111     // --- use summed wtest for outlier detection ---
02112     // #%// BK 21-Oct-2003
02113     matrix<real8> wtestsum = sqr(wtestL)+sqr(wtestP);// (Nobs,1)
02114     real8 maxwsum = max(wtestsum,winL,dumm);// idx to remove
02115     INFO << "Detected outlier:  summed sqr.wtest = " << maxwsum
02116          << "; observation: " << winL
02117          << "; window number: " 
02118          <<  Data(winL,0);
02119     INFO.print();
02120 
02121 
02122     // ______ Test if we are done yet ______
02123     if (Nobs <= Nunk)
02124       {
02125       WARNING.print("NO redundancy!  Exiting iterations.");
02126       DONE = 1;// cannot remove more than this
02127       }
02128 // seems something fishy here...
02129 //    if (max(overallmodeltestL,overallmodeltestP) < 1.0)
02130 //      {
02131 //      INFO.print("OMTs accepted, not iterating anymore (final solution reached).");
02132 //      DONE = 1;// ok (?).
02133 //      }
02134     if (max(maxwL,maxwP) <= CRIT_VALUE)// all tests accepted?
02135       {
02136       INFO.print("All outlier tests accepted! (final solution computed)");
02137       DONE = 1;// yeah!
02138       }
02139     if (ITERATION >= MAX_ITERATIONS)
02140       {
02141       INFO.print("max. number of iterations reached (exiting loop).");
02142       DONE = 1;// we reached max. (or no max_iter specified)
02143       }
02144 
02145     // ______ Only warn if last iteration has been done ______
02146     if (DONE == 1) 
02147       {
02148       // ___ use trace buffer to store string, remember to rewind it ___
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       // ___ use trace buffer to store string, remember to rewind it ___
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       // if a priori sigma is correct, max wtest should be something like 1.96
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       // ______Test of Jaron Samson's thesis: not ok/ depends on SIGMA ...______
02173       real8 expected = 1. / sqr(28);    //784.;
02174       //real4 expected = 1. / 400.;             // WRONG, 1./sqr(28) ???
02175       expected /= sqr(SIGMAL);                  // correct for variance factor
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       }// Only warn when done iterating.
02184 
02185 
02186     ITERATION++;// update counter here!
02187     }// iterations remove outliers
02188   
02189 
02190 
02191   // ______ Create dump file for making plots ______
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   // ______ Only way to format in c++ since stupid iomanip dont work? ______
02207   FILE *cpm;
02208   cpm=fopen("CPM_Data","a");
02209   //for (i=0; i<Nobs; i++)
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   // ====== Write results to scratch files ======
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 // ______ Add coefficient number behind value ______
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 // ______ Add coefficient number behind value ______
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                  //<< "\n* End_coregpm:_NORMAL"
02371                  << "\n* End_" << processcontrol[pr_i_coregpm] << "_NORMAL"
02372                  << "\n*******************************************************************\n";
02373   scratchresfile.close();
02374 
02375 
02376 // ====== Compute offsets for corners ======
02377 // BK 18-May-2000
02378   // read rhsL from file due top format... double
02379   matrix<real8> Lcoeff = readcoeff("scratchrescpm",
02380                      "Estimated_coefficientsL:",Ncoeffs(DEGREE));
02381   // read rhsP from file due top format... double
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   // ====== Dump evaluated polynomial if requested ======
02414   // BK 17-May-2000
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();// empty buffer to be sure
02421     TRACE << "offsetazi_" << originalmaster.lines() 
02422                <<          "_" << originalmaster.pixels() 
02423                << ".r4";
02424     char fileazi[ONE27];
02425     strcpy(fileazi,TRACE.get_str());
02426     TRACE.print();// empty buffer to be sure
02427 
02428     ofstream dumpfile;
02429     openfstream(dumpfile,fileazi,true);
02430     bk_assert(dumpfile,fileazi,__FILE__,__LINE__);
02431 
02432     // polyval both standing x,y... (?)
02433     // matrix<real4> p_axis(1,originalmaster.pixels());
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;                  // multilook==1 ?
02438     normalize(p_axis,minP,maxP);
02439 
02440     // azimuth
02441     for (i =originalmaster.linelo;
02442          i<=originalmaster.linehi; ++i) // all lines
02443       {
02444       //l_axis(0,0) = i;
02445       //normalize(l_axis,minL,maxL);
02446       l_axis(0,0) = normalize(real8(i),minL,maxL);
02447       //matrix<real4> MODEL = polyval(l_axis,p_axis,rhsL);
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     // ______ same for range ______
02460     PROGRESS.print("Started dumping evaluated model range.");
02461     TRACE.print();// empty buffer to be sure
02462     TRACE << "offsetrange_" << originalmaster.lines() 
02463              << "_" << originalmaster.pixels() << ".r4";
02464     char filerange[ONE27];
02465     strcpy(filerange,TRACE.get_str());
02466     TRACE.print();// empty buffer to be sure
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) // all lines
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   // ====== Tidy up ======
02488   PROGRESS.print("finished computation of coregistration parameters.");
02489   } // END coregpm
02490 
02491 
02492 
02493 /****************************************************************
02494  *    getofffile                                                *
02495  *                                                              *
02496  * Returns matrix (real4) with data of fine coreg from file     *
02497  *   mat(i,0)=window number                                     *
02498  *   mat(i,1)=position: line coordinate                         *
02499  *   mat(i,2)=position: pixels coordinate                       *
02500  *   mat(i,3)=offset: line direction                            *
02501  *   mat(i,4)=offset: pixle direction                           *
02502  *   mat(i,5)=correlation:                                      *
02503  * searches for "Number_of_correlation_windows:"                *
02504  *                                                              *
02505  *    Bert Kampes, 24-Feb-1999                                  *
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   //ifstream infile(file, ios::in);
02521   //if (!infile)
02522   //  {
02523   //  DEBUG << "getofffile: problem opening file: " << file;
02524   //  DEBUG.print();
02525   //  }
02526 
02527 // ======Search file for data section======
02528   while (infile)
02529     {
02530     infile >> word;
02531     if (strcmp("Number_of_correlation_windows:",word))  // no pattern match.
02532       {
02533       infile.getline(dummyline,ONE27,'\n');             // goto next line.
02534       }
02535     else                                                // in data section
02536       {
02537       foundsection=true;
02538       int32 N;                                          // number of points
02539       infile >> N;
02540       infile.getline(dummyline,ONE27,'\n');             // next line
02541       infile.getline(dummyline,ONE27,'\n');             // skip line with info
02542       int32 pos  = infile.tellg();                      // position of start data
02543       int32 Nobs = 0;                                   // number points > threshold
02544       real4 winnumber, posL, posP, offL, offP, corr;    // on file
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');           // goto next data record
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);                                // return to start data
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');           // goto next data record
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       } // else
02582     } // file
02583 
02584 // ______Tidy up______
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   // --- return a dummy here since some compiler like that ---
02593   return matrix<real4>(999,999);// BK 07-Apr-2003
02594   } // END getofffile
02595 
02596 
02597 
02598 /****************************************************************
02599  *    cc4                                                       *
02600  *                                                              *
02601  * cubic convolution 4 points                                   *
02602  *                                                              *
02603  * input:                                                       *
02604  *  - x-axis                                                    *
02605  * output:                                                      *
02606  *  - y=f(x); function evaluated at x                           *
02607  *                                                              *
02608  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END cc4
02635 
02636 
02637 /****************************************************************
02638  *    cc6                                                       *
02639  *                                                              *
02640  * cubic convolution 6 points                                   *
02641  *                                                              *
02642  * input:                                                       *
02643  *  - x-axis                                                    *
02644  * output:                                                      *
02645  *  - y=f(x); function evaluated at x                           *
02646  *                                                              *
02647  *    Bert Kampes, 16-Mar-1999                                  *
02648  * corrected (alfa+beta)->(alfa-beta) after correction in paper *
02649  * by Ramon Hanssen                                             *
02650  *    Bert Kampes, 16-Mar-1999                                  *
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 //      y(i,0) = (alpha+beta+2)*xx2*xx - (alpha+beta+3)*xx2 + 1;
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   } // END cc6
02682 
02683 
02684 /****************************************************************
02685  *    ts6                                                       *
02686  *                                                              *
02687  * truncated sinc 6 points                                      *
02688  *                                                              *
02689  * input:                                                       *
02690  *  - x-axis                                                    *
02691  * output:                                                      *
02692  *  - y=f(x); function evaluated at x                           *
02693  *                                                              *
02694  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END ts6
02711 
02712 
02713 /****************************************************************
02714  *    ts8                                                       *
02715  *                                                              *
02716  * truncated sinc 8 points                                      *
02717  *                                                              *
02718  * input:                                                       *
02719  *  - x-axis                                                    *
02720  * output:                                                      *
02721  *  - y=f(x); function evaluated at x                           *
02722  *                                                              *
02723  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END ts8
02740 
02741 
02742 /****************************************************************
02743  *    ts16                                                      *
02744  *                                                              *
02745  * truncated sinc 16 points                                     *
02746  *                                                              *
02747  * input:                                                       *
02748  *  - x-axis                                                    *
02749  * output:                                                      *
02750  *  - y=f(x); function evaluated at x                           *
02751  *                                                              *
02752  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END ts16
02768 
02769 
02770 /****************************************************************
02771  *    rect                                                      *
02772  *                                                              *
02773  * rect function for matrix (stepping function?)                *
02774  *                                                              *
02775  * input:                                                       *
02776  *  - x-axis                                                    *
02777  * output:                                                      *
02778  *  - y=f(x); function evaluated at x                           *
02779  *                                                              *
02780  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END rect
02797 
02798 
02799 /****************************************************************
02800  *    tri                                                       *
02801  *                                                              *
02802  * tri function for matrix (piecewize linear?, triangle)        *
02803  *                                                              *
02804  * input:                                                       *
02805  *  - x-axis                                                    *
02806  * output:                                                      *
02807  *  - y=f(x); function evaluated at x                           *
02808  *                                                              *
02809  *    Bert Kampes, 16-Mar-1999                                  *
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   } // END tri
02825 
02826 
02827 /****************************************************************
02828  *    knab                                                      *
02829  *                                                              *
02830  * KNAB window of N points, oversampling factor CHI             *
02831  *                                                              *
02832  * defined by: Migliaccio IEEE letters vol41,no5, pp1105,1110, 2003 *
02833  * k = sinc(x).*(cosh((pi*v*L/2)*sqrt(1-(2.*x./L).^2))/cosh(pi*v*L/2));
02834  *                                                              *
02835  * input:                                                       *
02836  *  - x-axis                                                    *
02837  *  - oversampling factor of bandlimited sigal CHI              *
02838  *  - N points of kernel size                                   *
02839  * output:                                                      *
02840  *  - y=f(x); function evaluated at x                           *
02841  *                                                              *
02842  *    Bert Kampes, 22-DEC-2003                                  *
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   } // END knab
02862 
02863 
02864 /****************************************************************
02865  *    resample                                                  *
02866  *                                                              *
02867  * Resample slave to master grid based on coregistration        *
02868  *  parameters.                                                 *
02869  * if dbow==0 then default to overlap, else dbow,               *
02870  *  write 0's where it does not overlap                         *
02871  * (later at interf.comp. if master<slave then doris exits! bug)*
02872  *                                                              *
02873  * input:                                                       *
02874  *  - inputoptions                                              *
02875  * output:                                                      *
02876  *  - void (file)                                               *
02877  *                                                              *
02878  *    Bert Kampes, 16-Mar-1999                                  *
02879  * added DBOW master add zeros.                                 *
02880  #%// BK 21-Aug-2000BOW master add zeros.                       *
02881  * shift data to center of spectrum before resampling, shift    *
02882  * back afterwards. (see e.g. thesis Geudtner)                  *
02883  #%// BK 09-Nov-2000                                            *
02884  * Seems to be a bug in shifting the data spectrum if more      *
02885  * buffers are used, working on it.                             *
02886  * (Increase FORSURE variable if crash)                         *
02887  #%// BK 19-Nov-2000                                            *
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,          // coregistration parameters
02895         const matrix<real8>     &cpmP           // coregistration parameters
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   // ___ Handle input ___
02902   const uint BUFFERMEMSIZE = generalinput.memory;       // Bytes
02903   const int32 Npoints      = resampleinput.method%100;  // #pnts interpolator
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   // ______ Normalize data for polynomial ______
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   // ______ For KNAB window only if requested ______
02926   real4 CHI1 = slave.prf/slave.abw;// oversampling factor az
02927   real4 CHI2 = (slave.rsr2x/2.0)/slave.rbw;// oversampling factor rg
02928   real4 CHI  = (CHI1+CHI2)/2.0;// mean oversampling factor
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   // ======Create lookup table======
02943   // ______e.g. four point interpolator
02944   // ______interpolating point: p=6.4925
02945   // ______required points: 5, 6, 7, 8
02946   // ______kernel number from lookup table: floor(.4925*INTERVAL+.5)
02947   // ______ table[0]= 0 1 0 0 ;table[INTERVAL]= 0 0 1 0
02948   // ______intervals in lookup table: dx
02949   register int32 i;
02950   // ______ it seems 20 is enough, but as fast as 40... ______
02951   //const int32 INTERVAL = 20;                          // precision 1./INTERVAL
02952   //const int32 INTERVAL = 40;                          // precision: 1./INTERVAL
02953   const int32 INTERVAL = 100;                           // precision: 1./INTERVAL
02954   const real8 dx = 1./INTERVAL;                         // interval look up table
02955   const int32 Ninterval = INTERVAL + 1;                 // size of lookup table
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;                   // start at [-1 0 1 2]
02962 
02963   // ______ Lookup table complex because of multiplication with complex ______
02964   // ______ This is faster if veclib is used, slower if for loop... ______
02965   // ______ Could make a loopkup table for azimuth and range and ______
02966   // ______ shift spectrum of azi kernel with doppler centroid (not done) ______
02967   matrix<complr4> *pntM[Ninterval];
02968   // ______ same axis required for shift azimuth spectrum as used ______
02969   // ______ for kernel to avoid phase shift (Raffaele Nutricato) ______
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);// used only for azishift
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;// to shift kernelL use: k*=exp(-i*2pi*axis*fdc/prf) 
03014     x_axis -= dx;               // Note: 'wrong' way (mirrored)
03015     }
03016   // Bert: usage: 
03017   //  pntM[0]->showdata(); or (*pntM[0][0]).showdata(); 
03018   PROGRESS.print("Resample: lookup table created (kernel and axis).");
03019 
03020   // ______Save some time by computing degree here______
03021   int32 degree_cpmL = degree(cpmL.size());
03022   int32 degree_cpmP = degree(cpmP.size());
03023 
03024 
03025   // ______Corners of overlap in slave system______
03026   // ______offset = A(slave system) - A(master system)______
03027   // BK 21-09-00, get approx overlap here, what is diff with below???
03028   window approxoverlap = getoverlap(master,slave,cpmL,cpmP);
03029 
03030   // compute it here as well??? (doesn't hurt)
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   // ______Compute if it is possible to resample entire overlap______
03066   // ______Assume +- linear transformation model?______
03067   //  min required line coord: ceil(min(oL00,oL0N)) - Npointsd2
03068   //  max required line coord: floor(max(oLN0,oLNN)) + Npointsd2
03069   //  difference with present: correction = slave.linelo - min.required
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   //  window overlap;   // in master system -> res file for update productinfo
03087   // no problem with uint initialization
03088   window overlap(approxoverlap.linelo + correctlinelo,
03089                  approxoverlap.linehi + correctlinehi,
03090                  approxoverlap.pixlo  + correctpixlo,
03091                  approxoverlap.pixhi  + correctpixhi);
03092 
03093   // ====== Adjust overlap possibly for RS_DBOW card ======
03094   int32 write0lines1  = 0;                      // DBOW card, 0's at start
03095   int32 write0linesN  = 0;
03096   int32 write0pixels1 = 0;
03097   int32 write0pixelsN = 0;
03098   if (!(resampleinput.dbow.linelo == 0 &&       // as such initialized by readinput
03099         resampleinput.dbow.linehi == 0 &&
03100         resampleinput.dbow.pixlo  == 0 &&
03101         resampleinput.dbow.pixhi  == 0    ))
03102     {
03103     // ______ Check if overlap is large enough to contain DBOW ______
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;   // smaller window selected
03127     write0linesN  = -overlap.linehi + resampleinput.dbow.linehi;
03128     if ( write0linesN < 0 ) write0linesN = 0;   // smaller window selected
03129     write0pixels1  =  overlap.pixlo - resampleinput.dbow.pixlo;
03130     if ( write0pixels1 < 0 ) write0pixels1 = 0; // smaller window selected
03131     write0pixelsN  = -overlap.pixhi + resampleinput.dbow.pixhi;
03132     if ( write0pixelsN < 0 ) write0pixelsN = 0; // smaller window selected
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;       // correct it
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;       // correct it
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;         // correct it
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;         // correct it
03165     } // adjust overlap
03166 
03167 
03168   // ______ Buffersize output matrix ______
03169   const int32 Npointsxsize = Npoints*sizeofcr4; // size for memcpy (fill PART)
03170   const int32 npixels      = slave.currentwindow.pixels();
03171   const real8 bytesperline = sizeofcr4 * npixels;
03172   // ___ COMMENTED OUT, OLD WAY SHIFT DATA, now shiftkernel ___
03173   //const real8 bigmatrices  = (resampleinput.shiftazi==true) ? 6:2.5;
03174   const real8 bigmatrices  = 2.5;
03175   const int32 nlines       = int32((BUFFERMEMSIZE/bigmatrices)/bytesperline);
03176 
03177   // ______ Declare/allocate matrices ______
03178   matrix<complr4> BUFFER;                       // load after output is written
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;                           // must be non-constant
03185   int32 ONEint = 1;                             // must have pointer to 1
03186   complr4 c4alpha(1.,0.);
03187   complr4 c4beta(0.,0.);
03188   STUPID_cr4 ANS;                               // VECLIB struct return type
03189 #endif
03190 
03191 
03192   // ====== Open output file ======
03193   ofstream ofile;
03194   openfstream(ofile,resampleinput.fileout,generalinput.overwrit);
03195   bk_assert(ofile,resampleinput.fileout,__FILE__,__LINE__);
03196 
03197   // ________ First write zero lines if appropriate (DBOW) ______
03198   register int32 thisline, thispixel;           // local loop variable
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   // ______ Info ______
03227   INFO << "Overlap window: "
03228        << overlap.linelo << ":" << overlap.linehi << ", "
03229        << overlap.pixlo  << ":" << overlap.pixhi;
03230   INFO.print();
03231 
03232   // ______ Progress messages ______
03233   int32 percent    = 0;
03234   int32 tenpercent = int32((overlap.lines()/10)+.5);  // round
03235   if (tenpercent==0) tenpercent = 1000;                    // avoid error: x%0
03236 
03237   // ====== Resample all lines that are requested ======
03238   bool newbufferrequired = true;                // read initial slave buffer
03239   register int32 linecnt = -1;                  // indicate output buffer full
03240   int32 firstline;                              // slave system
03241   register int32 line;                          // loop counter master system
03242   register int32 pixel;                         // loop counter master system
03243   for (line=overlap.linelo; line<=overlap.linehi; line++)
03244     {
03245     // ______ Progress messages ______
03246     if (((line-overlap.linelo)%tenpercent)==0)
03247       {
03248       PROGRESS << "RESAMPLE: " << setw(3) << percent << "%";
03249       PROGRESS.print();
03250       percent += 10;
03251       }
03252 
03253     // ====== Write RESULT to disk if it is full (write last bit at end) ======
03254     if (linecnt==RESULT.lines()-1)              // ==nlines
03255       {
03256       newbufferrequired = true;                 // do load slave from file
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       // ______ Shift azimuth spectrum of resampled data back to fDC ______
03264       // ______ Require pixelnumber in slave system == pix_master + offsetP ______
03265       // ______ +shift for forward shifting ______
03266       // ___ COMMENTED OUT, OLD WAY SHIFT DATA ___
03267 //      if (resampleinput.shiftazi==true)
03268 //      {
03269 //        real4 shift = overlap.pixlo +
03270 //        polyval(normalize(real4(line-(linecnt/2.)),minL,maxL),// middle line
03271 //                normalize(real4(overlap.pixlo),minP,maxP),    // first pixel
03272 //                cpmP,degree_cpmP);                            // e.g. 2.5232
03273 //        if (shift < 0)
03274 //          ERROR("sorry, but this really is impossible.");
03275 //        shiftazispectrum(RESULT,slave,shift);
03276 //      }
03277       // ______ Actually write ______
03278       switch (resampleinput.oformatflag)
03279         {
03280         case FORMATCR4:
03281           {
03282           // old, now first write zeropixels...: ofile << RESULT;
03283           const complr4 zerocr4(0,0);
03284           for (thisline=0; thisline<RESULT.lines(); ++thisline)
03285             {
03286             // ______ Write zero pixels at start ______
03287             for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03288               {
03289               ofile.write((char*)&zerocr4,sizeofcr4);
03290               }
03291             // ______ WRITE the interpolated data per row ______
03292             ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeof(RESULT(0,0)));
03293             // ______ Write zero pixels at end ______
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             // ______ Write zero pixels at start ______
03308             for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03309               {
03310               ofile.write((char*)&zeroci16,sizeofci16);
03311               }
03312             // ______ Write the interpolated data per row ______
03313             for (thispixel=0; thispixel<RESULT.pixels(); ++thispixel)
03314               {
03315               // no default conversion, this seems slow, test this (BK)
03316               castedresult = cr4toci2(RESULT(thisline,thispixel));
03317               ofile.write((char*)&castedresult,sizeofci16);
03318               }
03319             // ______ Write zero pixels at end ______
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 // output buffer not full yet
03334       {
03335       linecnt++;
03336       }
03337 
03338     // ====== Read slave buffer if justwritten || firstblock ======
03339     if (newbufferrequired==true)
03340       {
03341       newbufferrequired = false;        // only load after output written
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       //const int32 FORSURE = 5;                // buffer larger 2*FORSURE start/end
03362       const int32 FORSURE = 25;         // buffer larger 2*FORSURE start/end
03363       firstline -= FORSURE;
03364       lastline  += FORSURE;
03365       // ______ Don't compare apples with pears, uint<->int! ______
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       // ______ Fill slave BUFFER from disk ______
03371       window winslavefile(firstline, lastline,  // part of slave loaded
03372                            slave.currentwindow.pixlo,   // from file in BUFFER.
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       // ______ Shift azimuthspectrum from fDC to center ______
03380       // ______ - for inverse shift, pixlo for evaluating fDC ______
03381       // ___ COMMENTED OUT, OLD WAY SHIFT DATA ___
03382 //      if (resampleinput.shiftazi==true)
03383 //        shiftazispectrum(BUFFER,slave,-real4(slave.currentwindow.pixlo));
03384       } // ___end: Read new slave buffer to resample outputbuffer
03385 
03386 
03387     // ====== Actual resample all pixels this output line ======
03388     for (pixel=overlap.pixlo; pixel<=overlap.pixhi; pixel++)
03389       {
03390       // ______ Compute coregistration ______
03391       // bk 25-10-99 why don't i do this per buffer, that's faster. (but more mem)
03392       //interpL = line  + polyval(line,pixel,cpmL,degree_cpmL); // e.g. 255.35432
03393       //interpP = pixel + polyval(line,pixel,cpmP,degree_cpmP); // e.g. 2.5232
03394       // ______ BK USE normalized coordinates, do this smarter .... !!!!
03395       const real4 interpL = line  + 
03396         polyval(normalize(real4(line),minL,maxL),
03397                 normalize(real4(pixel),minP,maxP),
03398                 cpmL,degree_cpmL);                              // e.g. 255.35432
03399       const real4 interpP = pixel +
03400         polyval(normalize(real4(line),minL,maxL),
03401                 normalize(real4(pixel),minP,maxP),
03402                 cpmP,degree_cpmP);                              // e.g. 2.5232
03403 
03404       // ______ Get correct lines for interpolation ______
03405       // ______ Note: following 4 not ok, floor(x) != ceil(x)-1 ______
03406       // ______ correct to int32 for speed, floor is slower, all positive indices.
03407       //int32 firstL     = ceil(interpL) - Npointsd2;           // e.g. 254 (5 6 7)
03408       //int32 firstP     = ceil(interpP) - Npointsd2;           // e.g. 1   (2 3 4)
03409       //int32 kernelnoL  = floor(.5+interpLdec*INTERVAL) - 1;   // lookup table index
03410       //int32 kernelnoP  = floor(.5+interpPdec*INTERVAL) - 1;   // lookup table index
03411       const int32 fl_interpL = int32(interpL);
03412       const int32 fl_interpP = int32(interpP);
03413       const int32 firstL     = fl_interpL - Npointsd2m1;        // e.g. 254 (5 6 7)
03414       const int32 firstP     = fl_interpP - Npointsd2m1;        // e.g. 1   (2 3 4)
03415       const real4 interpLdec = interpL - fl_interpL;            // e.g. .35432
03416       const real4 interpPdec = interpP - fl_interpP;            // e.g. .5232
03417 
03418       // ___ copy kernels here, change kernelL if required ___
03419       // BK 26-Oct-2002
03420       int32 kernelnoL  = int32(interpLdec*INTERVAL + .5);       // lookup table index
03421       int32 kernelnoP  = int32(interpPdec*INTERVAL + .5);       // lookup table index
03422       matrix<complr4> kernelL = (*pntM[kernelnoL]);// local copy
03423       matrix<complr4> kernelP = (*pntM[kernelnoP]);// local copy
03424 
03425 
03426 #ifdef __DEBUG
03427       // ______This shouldn't be possible...______
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;               // with next pixel
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;               // with next pixel
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;               // with next pixel
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;               // with next pixel
03451         }
03452 #endif
03453 
03454 
03455       // ______ Shift azimuth kernel with fDC before interpolation ______
03456       if (resampleinput.shiftazi==true)
03457         {
03458         // ___ Doppler centroid is function of range only ____
03459         //const real8 tau = interpP/(slave.rsr2x/2.0);
03460         //const real8 fDC = slave.f_DC_a0 + 
03461         //                  slave.f_DC_a1*tau +
03462         //                  slave.f_DC_a2*sqr(tau);
03463         //const real4 tmp = 2.0*PI*fDC/slave.prf;
03464         //const real4 tmp = 2.0*PI*pix2fdc(interpP,slave)/slave.prf;
03465         const real4 tmp = 2.0*PI*slave.pix2fdc(interpP)/slave.prf;
03466         // ___ to shift spectrum of convolution kernel to fDC of data, multiply
03467         // ___ in the space domain with a phase trend of -2pi*t*fdc/prf
03468         // ___ (to shift back (no need) use +fdc)
03469         for (i=0; i<Npoints; ++i)
03470           {
03471           //complr4 t = complr4(cos(i*2.0*PI*fDC/slave.prf),
03472           //                   -sin(i*2.0*PI*fDC/slave.prf)); // exp(i*trend)
03473           //kernelL(i,0) *= t;
03474           // ___ Modify kernel, shift spectrum to fDC ___
03475           real4 t       = ((*pntAxis[kernelnoL])(i,0))*tmp;
03476           kernelL(i,0) *= complr4(cos(t),-sin(t));// note '-' (see manual)
03477           }
03478         }
03479 
03480 
03481         // ______ For speed: define setdata internally (memcpy) ______
03482         //linelo = firstL - firstline;
03483         //pixlo  = firstP - 1;          // !! -1 should be -currentwin.pixlo
03484         for (i=0; i<Npoints; i++)
03485           memcpy(PART[i],
03486           BUFFER[i+firstL-firstline]+
03487           firstP-slave.currentwindow.pixlo,Npointsxsize);
03488 
03489         // ====== Some speed considerations ======
03490 #ifdef __USE_VECLIB_LIBRARY__
03491         // ______Compute PART * kernelP______
03492         //cgemv("T",&Np,&Np,&c4alpha, PART[0],&Np,
03493         //      (*pntM[int32(interpPdec*INTERVAL+.5)])[0],&ONEint,
03494         //      &c4beta,TMPRES[0],&ONEint,1);
03495         cgemv("T",&Np,&Np,&c4alpha, PART[0],&Np,
03496               kernelP[0],&ONEint,
03497               &c4beta,TMPRES[0],&ONEint,1);
03498 
03499         // ______Compute Result * kernelL; put in matrix RESULT______
03500         //ANS = cdotu(&Np,TMPRES[0],&ONEint,
03501         //     (*pntM[int32(interpLdec*INTERVAL+.5)])[0],&ONEint);
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         // ______ NOvec slower, but ok: better use CR4*R4 ______
03506         //RESULT(linecnt,pixel-overlap.pixlo) = 
03507         //((matTxmat(PART * (*pntM[kernelnoP]), (*pntM[kernelnoL])))(0,0));
03508         RESULT(linecnt,pixel-overlap.pixlo) = 
03509              ((matTxmat(PART*kernelP, kernelL))(0,0));
03510 #endif // VECLIB 
03511       } // for all pixels in overlap
03512     } // for all lines in overlap
03513 
03514 
03515   // ______ Write last lines of Result to disk (filled upto linecnt) ______
03516   DEBUG << "Writing slave: ["
03517        //<< overlap.linehi-linecnt+1 << ":" << overlap.linehi << ", "
03518        << overlap.linehi-linecnt << ":" << overlap.linehi << ", "
03519        << overlap.pixlo << ":" << overlap.pixhi
03520        << "] (master coord. system)";
03521   DEBUG.print();
03522   // ______ Shift azimuth spectrum of resampled data back to fDC ______
03523   // ______ coordinate of first column in slave system equals: ______
03524   // ______ offset defined as: cols=colm+offsetP ______
03525   // ______ This is not totally correct... (e.g., if slave is scaled back) ______
03526   // ______ +shift for forward shifting ______
03527   // ______ We shift to much here, but this should not mather. ______
03528 //  if (resampleinput.shiftazi==true)
03529 //    {
03530 //    real4 shift = overlap.pixlo + polyval(normalize(real4(
03531 //              overlap.linehi-(RESULT.lines()/2.)),minL,maxL), // middle line
03532 //                  normalize(real4(overlap.pixlo),minP,maxP),  // first pixel
03533 //                  cpmP,degree_cpmP);                          // e.g. 2.5232
03534 //    if (shift < 0)
03535 //      ERROR("sorry, this is impossible, use noshiftazi.");
03536 //    shiftazispectrum(RESULT,slave,shift);
03537 //    }
03538 
03539   // ______ Actually write ______
03540   switch (resampleinput.oformatflag)
03541     {
03542     case FORMATCR4:
03543       {
03544       const complr4 zerocr4(0,0);
03545       for (thisline=0; thisline<=linecnt; thisline++)
03546         {
03547         // ______ Write zero pixels at start ______
03548         for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03549           {
03550           ofile.write((char*)&zerocr4,sizeofcr4);
03551           }
03552         // ______ WRITE the interpolated data per row ______
03553         ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeofcr4);
03554         // ______ Write zero pixels at end ______
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         // ______ Write zero pixels at start ______
03569         for (thispixel=0; thispixel<write0pixels1; ++thispixel)
03570           {
03571           ofile.write((char*)&zeroci16,sizeofci16);
03572           }
03573         // ______ Write the interpolated data per row ______
03574         for (thispixel=0; thispixel<RESULT.pixels(); thispixel++)
03575           {
03576           castedresult = cr4toci2(RESULT(thisline,thispixel));
03577           //?? BK 6/1 castedresult = cr4toci2(RESULT[thisline][thispixel].real());
03578           // castedresult = compli16(RESULT[thisline][thispixel].real(),
03579           //                       RESULT[thisline][thispixel].imag());
03580           ofile.write((char*)&castedresult,sizeofci16);
03581           }
03582         // ______ Write zero pixels at end ______
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   // ====== Write last zero lines if appropriate (DBOW card) ======
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   // ====== Write results to slave resfile ======
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     //<< "\nData first shifted to center azimuth spectrum by: "
03705     //<< fDC...
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 // ______Tidy up______
03738 //  for (i=0;i<Ninterval;i++)   // like this ???
03739 //    delete    pntM[i];
03740 //    delete [] pntM[i];
03741   DEBUG.print("Exiting resample.");
03742 
03743   } // END resample
03744 

Generated on Fri Apr 22 15:57:50 2005 for Doris by doxygen 1.3.6