coregistration.cc

Go to the documentation of this file.

/*
 * Copyright (c) 1999-2003 Bert Kampes
 * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
 *
 * This file is part of Doris, the Delft o-o radar interferometric software.
 *
 * Doris program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Doris is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Publications that contain results produced by the Doris software should
 * contain an acknowledgment. (For example: The interferometric processing
 * was performed using the freely available Doris software package developed
 * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
 * University of Technology, or include a reference to: Bert Kampes and
 * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
 * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
 *
 */
/****************************************************************
 * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/coregistration.cc,v $     *
 * $Revision: 3.30 $                                            *
 * $Date: 2005/04/11 13:47:45 $                                 *
 * $Author: kampes $                                            *
 *                                                              *
 * -coarse based on orbits.                                     *
 * -coarse based on correlation with fft/in space domain.       *
 * -fine coregistration offset vector computation.              *
 * -computation coregistration parameters.                      *
 * -computation flat earth correction.                          *
 * -resampling of slave to master grid.                         *
 ****************************************************************/


#include "constants.hh"         // typedefs etc.
#include "ioroutines.hh"        // error function etc.
#include "utilities.hh"         // isodd
#include "coregistration.hh"    // header file
#include "conversion.hh"        // lp2xyz etc.
#include "refsystems.hh"        // for dsk WGS84
#include "orbitbk.hh"           // my orbit class
#include "slcimage.hh"          // my slc image class
#include "productinfo.hh"       // my 'products' class
#include "exceptions.hh"        // my exceptions class
#include "bk_baseline.hh"       // my exceptions class


#include <iomanip>              // setw
#include <cstdlib>              // system
#include <cmath>                // sqrt
#include <algorithm>            // max
#include <cstdio>               // some compilers, remove function



/****************************************************************
 *    coarseporbit                                              *
 *                                                              *
 * computes translation of slave w.r.t. master                  *
 * slave(some point) = master(same point) + trans(l,p) =>       *
 *  trans = slavecoordinates - mastercoordinates                *
 * uses orbits to find coordinates of center of master          *
 *  then solves for lin,pix of these cn. for slave.             *
 *                                                              *
 * input:                                                       *
 *  -                                                           *
 * output:                                                      *
 *  -                                                           *
 *                                                              *
 *    Bert Kampes, 12-Dec-1998                                  *
 ****************************************************************/
void coarseporbit(
        const input_ell &ell,
        const slcimage  &master,
        const slcimage  &slave, 
        orbit           &masterorbit,  // cannot be const for spline
        orbit           &slaveorbit,   // cannot be const for spline
        const BASELINE  &baseline)
  {
  TRACE_FUNCTION("coarseporbit (BK 12-Dec-1998)");
  const int16   MAXITER   = 10;                 // maximum number of iterations
  const real8   CRITERPOS = 1e-6;                       // 1micrometer
  const real8   CRITERTIM = 1e-10;              // seconds (~10-6 m)

  // ______Initial values______         :master.approxcentreoriginal.x .y .z
  // ______Window______                 :master.currentwindow.linelo/hi , pixlo/hi
  // ______Time______                   :master.t_azi0/N , t_range0/N

  // ______Get (approx) center pixel of current window master______
  const uint cen_lin = (master.currentwindow.linelo+master.currentwindow.linehi)/2;
  const uint cen_pix = (master.currentwindow.pixlo +master.currentwindow.pixhi) /2;
  const real8 HEI = 0.0;

  // ______ Compute x,y,z (fill P) ______
  // ______ P.x/y/z contains (converged) solution ______
  cn P;
  const int32 lp2xyziter =
    lp2xyz(cen_lin,cen_pix,ell,master,masterorbit,P,MAXITER,CRITERPOS);

  // ______Compute line,pixel for slave of this xyz______
  real8 lin,pix;
  const int32 xyz2lpiter = 
    xyz2lp(lin,pix,slave,slaveorbit,P,MAXITER,CRITERTIM);

  // ______ Some extra parameters (not used, just info) ______ // BK 19-Oct-2000
  const int Bt          = Btemp(master.utc1,slave.utc1);
  // ______ Modeled quantities ______
  const real8 Bperp     = baseline.get_bperp(cen_lin,cen_pix,HEI);
  const real8 Bpar      = baseline.get_bpar(cen_lin,cen_pix,HEI);
  const real8 theta     = rad2deg(baseline.get_theta(cen_lin,cen_pix,HEI));
  const real8 inc_angle = rad2deg(baseline.get_theta_inc(cen_lin,cen_pix,HEI));
  // ______ Derived quantities ______
  const real8 B         = baseline.get_b(cen_lin,cen_pix,HEI);
  const real8 alpha     = rad2deg(baseline.get_alpha(cen_lin,cen_pix,HEI));
  const real8 Bh        = baseline.get_bhor(cen_lin,cen_pix,HEI);
  const real8 Bv        = baseline.get_bvert(cen_lin,cen_pix,HEI);

  const real8 Hamb      = baseline.get_hamb(cen_lin,cen_pix,HEI);
  const real8 orb_conv  = rad2deg(baseline.get_orb_conv(cen_lin,cen_pix,HEI));


/*   removed: #%// Bert Kampes, 06-Apr-2005
//  cn M,S;
//  xyz2orb(M,master,masterorbit,P); 
//  xyz2orb(S,slave,slaveorbit,P);  
//  real8 B, Bperp, Bpar, Bh, Bv, alpha, theta, Hamb;
//  BalphaBhBvBparBperpTheta(B,alpha,Bh,Bv,Bpar,Bperp,theta,M,P,S);
//  const cn r1 = P.min(M);
//  Hamb = (master.wavelength*r1.norm()*sin(theta)) / (2*Bperp);
//  theta = rad2deg(theta);
//  alpha = rad2deg(alpha);
//  // BK 28-Nov-2000:
//  real8 inc_angle = P.angle(M.min(P));
//  inc_angle = rad2deg(inc_angle);
//  const real8 m_tazi = master.line2ta(cen_lin);
//  const real8 s_tazi = slave.line2ta(lin);
//  cn Mdot = masterorbit.getxyzdot(m_tazi);
//  cn Sdot = slaveorbit.getxyzdot(s_tazi);
//  real8 orb_conv = Mdot.angle(Sdot);
//  orb_conv = rad2deg(orb_conv);
//  //const real8 heading = angle(Mdot,[1 0 0])?
*/

  // ______ offset = P_slave - P_master = lin - cen_lin ______
  INFO << "Estimated translation (l,p): "
       << floor(lin-cen_lin +.5) << ", "
       << floor(pix-cen_pix +.5);
  INFO.print();

  // ______ Write to tmp files ______
  ofstream scratchlogfile("scratchlogcoarse", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"coarseporbit: scratchlogcoarse",__FILE__,__LINE__);
  scratchlogfile << "\n\n*******************************************************************"
                 << "\n* COARSE_COREGISTRATION Orbits"
                 << "\n*******************************************************************"
                 << "\n(Approximate) center master (line,pixel,hei): "
                 <<  cen_lin << ", " << cen_pix << ", " << HEI
                 << "\nEllipsoid WGS84 coordinates of this pixel (x,y,z): ("
                 <<  P.x << ", " << P.y << ", " << P.z << ")"
                 << "\n(line,pixel) of these coordinates in slave: "
                 <<  lin << ", " << pix
                 << "\nEstimated translation slave w.r.t. master (l,p):"
                 <<  floor(lin-cen_lin +.5)                             // round
                 << ", "
                 <<  floor(pix-cen_pix +.5)                             // round
                 << "\nMaximum number of iterations: "
                 <<  MAXITER
                 << "\nCriterium for position (m): "
                 <<  CRITERPOS
                 << "\nCriterium for azimuth time (s): "
                 <<  CRITERTIM
                 << " (=~ " << CRITERTIM*7.e3 << "m)"
                 << "\nNumber of iterations conversion line,pixel to xyz: "
                 <<  lp2xyziter
                 << "\nNumber of iterations conversion xyz to line,pixel: "
                 <<  xyz2lpiter
                 << "\n*******************************************************************\n";
  scratchlogfile.close();

  // ______ give some extra info in resfile: Bperp, Bpar, Bh, Bv, Btemp ______
  ofstream scratchresfile("scratchrescoarse", ios::out | ios::trunc);
  bk_assert(scratchresfile,"coarseporbit: scratchrescoarse",__FILE__,__LINE__);
  scratchresfile.setf(ios::right, ios::adjustfield);
  scratchresfile
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_i_coarse]
    << "\n*******************************************************************"
    << "\nSome info for pixel: " << cen_lin << ", " << cen_pix << " (not used):"
    << "\n  Btemp:     [days]:  "
    << setw(10)
    << setiosflags(ios::right)
    << Bt                << "      \t// Temporal baseline"
    << "\n  Bperp      [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(Bperp) << "      \t// Perpendicular baseline"
    << "\n  Bpar       [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(Bpar)  << "      \t// Parallel baseline"
    << "\n  Bh         [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(Bh)    << "      \t// Horizontal baseline"
    << "\n  Bv         [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(Bv)    << "      \t// Vertical baseline"
    << "\n  B          [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(B)     << "      \t// Baseline (distance between sensors)"
    << "\n  alpha      [deg]:   "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(alpha) << "      \t// Baseline orientation"
    << "\n  theta      [deg]:   "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(theta) << "      \t// look angle"
    << "\n  inc_angle  [deg]:   "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(inc_angle) << "      \t// incidence angle"
    << "\n  orbitconv  [deg]:   "
    << setw(10)
    << setiosflags(ios::right)
    << orb_conv               << "      \t// angle between orbits"
    << "\n  Height_amb [m]:     "
    << setw(10)
    << setiosflags(ios::right)
    << onedecimal(Hamb)       << "      \t// height = h_amb*phase/2pi (approximately)"
    << "\nEstimated translation slave w.r.t. master (slave-master):"
    << "\nCoarse_orbits_translation_lines:  \t"
    <<  floor(lin-cen_lin +.5)                          // round
    << "\nCoarse_orbits_translation_pixels: \t"
    <<  floor(pix-cen_pix +.5)                          // round
    << "\n*******************************************************************"
    << "\n* End_" << processcontrol[pr_i_coarse] << "_NORMAL"
    << "\n*******************************************************************\n";

  // ______Tidy up______
  scratchresfile.close();
  PROGRESS.print("Coarse precise orbits coregistration finished.");
  } // END coarseporbit



/****************************************************************
 *    coarsecorrel                                              *
 *                                                              *
 * computes translation of slave w.r.t. master                  *
 * slave(some point) = master(same point) + trans(l,p) =>       *
 *  trans = slavecoordinates - mastercoordinates                *
 * uses correlation between magnitude of slave/master image     *
 *                                                              *
 * requires thingsa on disk, input                              *
 * input:                                                       *
 *  -                                                           *
 * output:                                                      *
 *  -                                                           *
 *                                                              *
 *    Bert Kampes, 12-Dec-1998                                  *
 ****************************************************************/
void coarsecorrel(
        const input_coarsecorr &coarsecorrinput,
        const slcimage         &minfo,
        const slcimage         &sinfo)
  {
  TRACE_FUNCTION("coarsecorrel (BK 12-Dec-1998)");

  char          dummyline[ONE27];                               // for errormessages
  const uint Mfilelines = minfo.currentwindow.lines();
  const uint Sfilelines = sinfo.currentwindow.lines();

  const uint Nwin         = coarsecorrinput.Nwin;               // number of windows
  const int32 initoffsetL = coarsecorrinput.initoffsetL;        // initila offset
  const int32 initoffsetP = coarsecorrinput.initoffsetP;        // initila offset
  uint MasksizeL          = coarsecorrinput.MasksizeL;          // size of correlation window
  uint MasksizeP          = coarsecorrinput.MasksizeP;          // size of correlation window
  const uint AccL         = coarsecorrinput.AccL;               // accuracy of initial offset
  const uint AccP         = coarsecorrinput.AccP;               // accuracy of initial offset
  bool pointsrandom = true;
  if (specified(coarsecorrinput.ifpositions))   // filename specified
    pointsrandom = false;                       // only use those points


//  INFO("Masksize ...

// ______Only odd Masksize possible_____
  bool forceoddl = false;
  bool forceoddp = false;
  if (!isodd(MasksizeL)) 
    {
    forceoddl = true; 
    MasksizeL+=1;                       // force oddness
    }
  if (!isodd(MasksizeP))
    {
    forceoddp = true; 
    MasksizeP+=1;                       // force oddness
    }

// ______Corners of slave in master system______
// ______offset = A(slave system) - A(master system)______
  const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
  const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
  const int32 sp0 = sinfo.currentwindow.pixlo  - initoffsetP;
  const int32 spN = sinfo.currentwindow.pixhi  - initoffsetP;

// ______Corners of overlap master,slave in master system______
//  const int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
//  const int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
//  const int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
//  const int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
 
// ______Corners of useful overlap master,slave in master system______
  const uint FEW=10;
  const uint l0 = uint(max(int32(minfo.currentwindow.linelo),sl0) + .5*MasksizeL + AccL + FEW);
  const uint lN = uint(min(int32(minfo.currentwindow.linehi),slN) - .5*MasksizeL - AccL - FEW);
  const uint p0 = uint(max(int32(minfo.currentwindow.pixlo),sp0)  + .5*MasksizeP + AccP + FEW);
  const uint pN = uint(min(int32(minfo.currentwindow.pixhi),spN)  - .5*MasksizeP - AccP - FEW);
  const window overlap(l0,lN,p0,pN);

// ______Distribute Nwin points over window______
// ______Centers(i,0): line, (i,1): pixel, (i,2) flagfromdisk______
// old  matrix<uint> Centers = distributepoints(Nwin,overlap);

  register int32 i;
  matrix<uint> Centers;
  if (pointsrandom)                             // no filename specified
    {
    Centers = distributepoints(real4(Nwin),overlap);
    }

  else  // read in points (center of windows) from file
    {
    Centers.resize(Nwin,3);
    //ifstream ifpos(coarsecorrinput.ifpositions, ios::in);
    ifstream ifpos;
    openfstream(ifpos,coarsecorrinput.ifpositions);
    bk_assert(ifpos,coarsecorrinput.ifpositions,__FILE__,__LINE__);
    uint ll,pp;
    for (i=0; i<Nwin; ++i)
      {
      ifpos >> ll >> pp;
      Centers(i,0) = uint(ll);                  // correct for lower left corner
      Centers(i,1) = uint(pp);                  // correct for lower left corner
      Centers(i,2) = uint(1);                   // flag from file
      ifpos.getline(dummyline,ONE27,'\n');              // goto next line.
      }
    ifpos.close();

// ______ Check last point ivm. EOL after last position in file ______
    if (Centers(Nwin-1,0) == Centers(Nwin-2,0) &&
        Centers(Nwin-1,1) == Centers(Nwin-2,1))
      {
      Centers(Nwin-1,0) = uint(.5*(lN + l0) + 27);      // random
      Centers(Nwin-1,1) = uint(.5*(pN + p0) + 37);      // random
      WARNING << "CC: there should be no EOL after last point in file: "
           << coarsecorrinput.ifpositions;
      WARNING.print();
      }

// ______ Check if points are in overlap ______
// ______ no check for uniqueness of points ______
    bool troubleoverlap = false;
    for (i=0; i<Nwin; ++i)
      {
      if (Centers(i,0) < l0)
        {
        troubleoverlap=true;
        WARNING << "COARSE_CORR: point from file: "
             << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
             << " outside overlap master, slave. New position: ";
        Centers(i,0) = l0 + l0-Centers(i,0);
        WARNING << Centers(i,0) << " " << Centers(i,1);
        WARNING.print();
        }
      if (Centers(i,0) > lN)
        {
        troubleoverlap=true;
        WARNING << "COARSE_CORR: point from file: "
             << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
             << " outside overlap master, slave. New position: ";
        Centers(i,0) = lN + lN-Centers(i,0);
        WARNING << Centers(i,0) << " " << Centers(i,1);
        WARNING.print();
        }
      if (Centers(i,1) < p0)
        {
        troubleoverlap=true;
        WARNING << "COARSE_CORR: point from file: "
             << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
             << " outside overlap master, slave. New position: ";
        Centers(i,1) = p0 + p0-Centers(i,1);
        WARNING << Centers(i,0) << " " << Centers(i,1);
        WARNING.print();
        }
      if (Centers(i,1) > pN)
        {
        troubleoverlap=true;
        WARNING << "COARSE_CORR: point from file: "
             << i+1 << " " << Centers(i,0) << " " << Centers(i,1)
             << " outside overlap master, slave. New position: ";
        Centers(i,1) = pN + pN-Centers(i,1);
        WARNING << Centers(i,0) << " " << Centers(i,1);
        WARNING.print();
        }
      }
    if (troubleoverlap) // give some additional info
      {
      WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
           << l0 << " " << lN << " " << p0 << " " << pN << ends;
      WARNING.print();
      }
    }


// ______Compute correlation of these points______
  //matrix<complr4> Mcmpl(2*AccL+MasksizeL,2*AccP+MasksizeP);   // on file
  //matrix<complr4> Scmpl(MasksizeL,MasksizeP);                         // on file
  matrix<complr4> Mcmpl;
  matrix<complr4> Scmpl;
  matrix<real4> Master;                         // amplitude master
  matrix<real4> Mask;                           // amplitude slave
  matrix<real4> Correl;                         // matrix with correlations
  matrix<real4> Result(Nwin,3);                 // R(i,0):delta l; (i,1):delta p; (i,2):correl

  //if (Nwin<10)
  //  PROGRESS("COARSE_CORR:  0%");
  // ______ Progress messages ______
  int32 percent    = 0;
  //int32 tenpercent = int32(floor(overlap.lines()/10)+.5);    // round
  int32 tenpercent = int32((Nwin/10)+.5);       // round
  if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
  for (i=0;i<Nwin;i++)
    {
    //if (Nwin>=10) if (!(i%((Nwin)/10)))               // indicate 10% progress
    if (i%tenpercent==0)
      {
      //int32 percentage = int32(100*(real4(i)/real4(Nwin)));
      PROGRESS << "COARSE_CORR: " << setw(3) << percent << "%" << ends;
      PROGRESS.print();
      percent += 10;
      }

// ______Center of window in master system______
    uint cenMwinL = Centers(i,0);
    uint cenMwinP = Centers(i,1);

    window master;                                      // size=masksize+2*acc.
    master.linelo = cenMwinL - (MasksizeL-1)/2 -AccL;   // ML is forced odd
    master.linehi = master.linelo + MasksizeL +2*AccL - 1;
    master.pixlo  = cenMwinP - (MasksizeP-1)/2 - AccP;  // MP is forced odd
    master.pixhi  = master.pixlo + MasksizeP +2*AccP - 1;

// ______Same points in slave system (disk)______
    window slavemask;                                   // size=masksize
    uint cenSwinL = cenMwinL + initoffsetL;             // adjust initoffset
    uint cenSwinP = cenMwinP + initoffsetP;             // adjust initoffset
    slavemask.linelo = cenSwinL - (MasksizeL-1)/2;      // ML is forced odd
    slavemask.linehi = slavemask.linelo + MasksizeL - 1;
    slavemask.pixlo  = cenSwinP - (MasksizeP-1)/2;      // MP is forced odd
    slavemask.pixhi  = slavemask.pixlo + MasksizeP - 1;

// ______Read windows from files, compute magnitude______
    //minfo.readdata(Mcmpl,master);
    //sinfo.readdata(Scmpl,slavemask);
    Mcmpl = minfo.readdata(master);
    Scmpl = sinfo.readdata(slavemask);
    Master = magnitude(Mcmpl); 
    Mask   = magnitude(Scmpl); 


// ______Compute correlation matrix and find maximum______
//    DEBUG("returning matrix is not necessary, could return max,l,p"); 
    Correl = correlate(Master,Mask); 
    uint L, P;
    real4 corr = max(Correl, L, P);             // returns also L,P

    uint relcenML = master.linehi - cenMwinL;   // system of matrix
    uint relcenMP = master.pixhi  - cenMwinP;   // system of matrix
    int32 reloffsetL = relcenML - L;
    int32 reloffsetP = relcenMP - P;
    int32 offsetL = reloffsetL + initoffsetL;   // estimated offset lines
    int32 offsetP = reloffsetP + initoffsetP;   // estimated offset pixels

    Result(i,0) = offsetL;
    Result(i,1) = offsetP;
    Result(i,2) = corr;
    }

// ______Get correct offsetL, offsetP______
  int32 offsetLines = -999;
  int32 offsetPixels = -999;
  getoffset(Result,offsetLines,offsetPixels);


// ______Write to files______
  ofstream scratchlogfile("scratchlogcoarse2", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"coarsecorrel: scratchlogcoarse2",__FILE__,__LINE__);
  scratchlogfile << "\n\n*******************************************************************"
                 << "\n* COARSE_COREGISTRATION: Correlation"
                 << "\n*******************************************************************"
                 << "\nNumber of correlation windows: \t"
                 <<  Nwin
                 << "\nCorrelation window size (l,p): \t"
                 <<  MasksizeL << ", " << MasksizeP;
    if (forceoddl) scratchlogfile << "(l forced odd) ";
    if (forceoddp) scratchlogfile << "(p forced odd)";
  scratchlogfile << "\nSearchwindow size (l,p): \t\t"
                 <<  MasksizeL + 2*AccL << ", " << MasksizeP + 2*AccP
                 << "\nNumber \tposl \tposp \toffsetl offsetp \tcorrelation\n";
  register int32 k;
  for (k=0; k<Nwin; k++)
    scratchlogfile << k << "\t" << Centers(k,0) 
                        << "\t" << Centers(k,1) 
                        << "\t" << Result(k,0) 
                        << "\t" << Result(k,1) 
                        << "\t" << Result(k,2) << endl;
  scratchlogfile << "Estimated total offset (l,p): \t"
                 <<  offsetLines << ", " << offsetPixels
                 << "\n*******************************************************************\n";
  scratchlogfile.close();

  ofstream scratchresfile("scratchrescoarse2", ios::out | ios::trunc);
  bk_assert(scratchresfile,"coarsecorrel: scratchrescoarse2",__FILE__,__LINE__);
  scratchresfile << "\n\n*******************************************************************"
                 << "\n*_Start_" << processcontrol[pr_i_coarse2]
                 << "\n*******************************************************************"
                 << "\nEstimated translation slave w.r.t. master:"
                 << "\nCoarse_correlation_translation_lines: \t"
                 <<  offsetLines                                // 1 digit after point?
                 << "\nCoarse_correlation_translation_pixels: \t"
                 <<  offsetPixels                               // 1 digit after point?
                 << "\n*******************************************************************"
                 //<< "\n* End_coarse_correlation:_NORMAL"
                 << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
                 << "\n*******************************************************************\n";
  scratchresfile.close();

// ______Tidy up______
  INFO << "Individually estimated translations (#, l, p, corr, offl, offp): ";
  INFO.print();
  for (k=0; k<Nwin; k++)
    {
    // "Estimate (#, l, p, corr, offl, offp): "
    INFO << k            << " \t"
         << Centers(k,0) << " \t"
         << Centers(k,1) << " \t"
         << Result(k,2)  << " \t"
         << Result(k,0)  << " \t"
         << Result(k,1)  << " \t";
    INFO.print();
    }
  INFO << "Estimated translation (l,p): "
       << offsetLines << ", " << offsetPixels << ends;
  INFO.print();
  PROGRESS.print("Coarse coregistration based on correlation finished.");
  } // END coarsecorrel



/****************************************************************
 *    coarsecorrelfft                                           *
 *                                                              *
 * computes translation of slave w.r.t. master                  *
 * slave(some point) = master(same point) + trans(l,p) =>       *
 *  trans = slavecoordinates - mastercoordinates                *
 * uses correlation between magnitude of slave/master image     *
 * uses fft to compute coherence, no subtraction of mean        *
 *                                                              *
 * requires thingsa on disk, input                              *
 * input:                                                       *
 *  -                                                           *
 * output:                                                      *
 *  -                                                           *
 *                                                              *
 *    Bert Kampes, 12-Dec-1998                                  *
 ****************************************************************/
void coarsecorrelfft(
    const input_coarsecorr &coarsecorrinput,
    const slcimage         &minfo,
    const slcimage         &sinfo)
  {
  TRACE_FUNCTION("coarsecorrelfft (BK 12-Dec-1998)");
  if (coarsecorrinput.method != cc_magfft)
    {
    PRINT_ERROR("unknown method, This routine is only for cc_magfft method.")
    throw(argument_error);
    }

  char          dummyline[ONE27];                       // for errormessages
  const uint Mfilelines   = minfo.currentwindow.lines();
  const uint Sfilelines   = sinfo.currentwindow.lines();

  const uint Nwin           = coarsecorrinput.Nwin;             // number of windows
  const int32 initoffsetL = coarsecorrinput.initoffsetL;        // initila offset
  const int32 initoffsetP = coarsecorrinput.initoffsetP;        // initila offset
  const uint MasksizeL    = coarsecorrinput.MasksizeL;  // size of correlation window
  const uint MasksizeP    = coarsecorrinput.MasksizeP;  // size of correlation window

  bool pointsrandom = true;
  if (specified(coarsecorrinput.ifpositions))   // filename specified
    pointsrandom = false;                       // only use these points

  // ______Only pow2 Masksize possible_____
  if (!ispower2(MasksizeL)) 
    {
    PRINT_ERROR("coarse correl fft: MasksizeL should be 2^n")
    throw(input_error);
    }
  if (!ispower2(MasksizeP))
    {
    PRINT_ERROR("coarse correl fft: MasksizeP should be 2^n")
    throw(input_error);
    }

  // ______Corners of slave in master system______
  // ______offset = [A](slave system) - [A](master system)______
  const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
  const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
  const int32 sp0 = sinfo.currentwindow.pixlo  - initoffsetP;
  const int32 spN = sinfo.currentwindow.pixhi  - initoffsetP;

// ______Corners of overlap master,slave in master system______
//  int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
//  int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
//  int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
//  int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
 
  // ______Corners of useful overlap master,slave in master system______
  const uint FEW=10;
  const uint l0 = max(int32(minfo.currentwindow.linelo),sl0) + FEW;
  const uint lN = min(int32(minfo.currentwindow.linehi),slN) - MasksizeL - FEW;
  const uint p0 = max(int32(minfo.currentwindow.pixlo),sp0)  + FEW;
  const uint pN = min(int32(minfo.currentwindow.pixhi),spN)  - MasksizeP - FEW;
  const window overlap(l0,lN,p0,pN);

  // ______Distribute Nwin points over window______
  // ______Minlminp(i,0): line, (i,1): pixel, (i,2) flagfromdisk______
  // old  matrix<uint> Minlminp = distributepoints(Nwin,overlap);
  register int32 i;

  matrix<uint> Minlminp;
  if (pointsrandom)                             // no filename specified
    {
    Minlminp = distributepoints(real4(Nwin),overlap);
    }

  else  // read in points (center of windows) from file
    {
    Minlminp.resize(Nwin,3);
    //ifstream ifpos(coarsecorrinput.ifpositions, ios::in);
    ifstream ifpos;
    openfstream(ifpos,coarsecorrinput.ifpositions);
    bk_assert(ifpos,coarsecorrinput.ifpositions,__FILE__,__LINE__);
    uint ll,pp;
    for (i=0; i<Nwin; ++i)
      {
      ifpos >> ll >> pp;
      Minlminp(i,0) = uint(ll -.5*MasksizeL);   // correct for lower left corner
      Minlminp(i,1) = uint(pp -.5*MasksizeP);   // correct for lower left corner
      Minlminp(i,2) = uint(1);                  // flag from file
      ifpos.getline(dummyline,ONE27,'\n');      // goto next line.
      }
    ifpos.close();

    // ______ Check last point ivm. EOL after last position in file ______
    if (Minlminp(Nwin-1,0) == Minlminp(Nwin-2,0) &&
        Minlminp(Nwin-1,1) == Minlminp(Nwin-2,1))
      {
      Minlminp(Nwin-1,0) = uint(.5*(lN + l0) + 27);     // random
      Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);     // random
      }

    // ______ Check if points are in overlap ______
    // ______ no check for uniqueness of points ______
    bool troubleoverlap = false;
    for (i=0; i<Nwin; ++i)
      {
      if (Minlminp(i,0) < l0)
        {
        troubleoverlap=true;
        WARNING << "COARSECORR: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,0) = l0 + l0-Minlminp(i,0);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,0) > lN)
        {
        troubleoverlap=true;
        WARNING << "COARSECORR: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,0) = lN + lN-Minlminp(i,0);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,1) < p0)
        {
        troubleoverlap=true;
        WARNING << "COARSECORR: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,1) = p0 + p0-Minlminp(i,1);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,1) > pN)
        {
        troubleoverlap=true;
        WARNING << "COARSECORR: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,1) = pN + pN-Minlminp(i,1);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      }
    if (troubleoverlap) // give some additional info
      {
      WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
           << l0 << " " << lN << " " << p0 << " " << pN;
      WARNING.print();
      }
    }



  // ______Compute coherence of these points______
  //matrix<complr4> Master(MasksizeL,MasksizeP);        // on file
  //matrix<complr4> Mask(MasksizeL,MasksizeP);  // on file
  matrix<complr4> Master;
  matrix<complr4> Mask;
  matrix<real4>   Result(Nwin,3);               // R(i,0):delta l; 
                                                //  R(i,1):delta p; R(i,2):correl
  // ______ Progress messages ______
  int32 percent    = 0;
  //int32 tenpercent = int32((overlap.lines()/10)+.5);    // round
  int32 tenpercent = int32((Nwin/10)+.5);       // round
  if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
  for (i=0; i<Nwin; ++i)
    {
    if (i%tenpercent==0)
      {
      PROGRESS << "COARSE_CORR: " << setw(3) << percent << "%";
      PROGRESS.print();
      percent += 10;
      }

    // ______Minlminp (lower left corners) of window in master system______
    uint minMwinL = Minlminp(i,0);
    uint minMwinP = Minlminp(i,1);
                                                        // size=masksize
    window master(minMwinL, minMwinL+MasksizeL-1,
                  minMwinP, minMwinP+MasksizeP-1);
    // ______Same points in slave system (disk)______
    window mask(minMwinL+initoffsetL, 
                minMwinL+initoffsetL+MasksizeL-1,
                minMwinP+initoffsetP,
                minMwinP+initoffsetP+MasksizeP-1);

    // ______Read windows from files______
    //minfo.readdata(Master,master);
    //sinfo.readdata(Mask,mask);
    Master = minfo.readdata(master);
    Mask   = sinfo.readdata(mask);

    matrix<real4> absMaster = magnitude(Master);
    matrix<real4> absMask = magnitude(Mask);

    // ______Coherence______
    // ______update offsetL/P______
    real4 offsetL, offsetP;
    real4 coheren;

    coheren = corrfft(absMaster,absMask,offsetL,offsetP);

    offsetL += initoffsetL;                     // estimated offset lines
    offsetP += initoffsetP;                     // estimated offset pixels

    Result(i,0) = offsetL;
    Result(i,1) = offsetP;
    Result(i,2) = coheren;
    } // for nwin

  // ______Position approx. with respect to center of window______
  // ______correct position array for center instead of lower left______
  for (i=0;i<Nwin;i++)
    {
    Minlminp(i,0) += uint(.5 * MasksizeL);
    Minlminp(i,1) += uint(.5 * MasksizeP);
    }

  // ______Get correct offsetL, offsetP______
  int32 offsetLines = -999;
  int32 offsetPixels = -999;
  getoffset(Result,offsetLines,offsetPixels);


  // ______Write to files______
  ofstream scratchlogfile("scratchlogcoarse2", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"coarsecorrelfft: scratchlogcoarse2",__FILE__,__LINE__);
  scratchlogfile << "\n\n*******************************************************************"
                 << "\n* COARSE_COREGISTRATION: Correlation"
                 << "\n*******************************************************************"
                 << "\nNumber of correlation windows: \t"
                 <<  Nwin
                 << "\nwindow size (l,p):             \t"
                 <<  MasksizeL << ", " << MasksizeP
//               << "\nSearchwindow size (l,p):       \t"
//               <<  MasksizeL + 2*AccL << ", " << MasksizeP + 2*AccP
                 << "\n\nNumber \tposL \tposP \toffsetL offsetP\tcorrelation\n";
  register int32 k;
  for (k=0; k<Nwin; k++)
    scratchlogfile << k << "\t" << Minlminp(k,0) 
                        << "\t" << Minlminp(k,1) 
                        << "\t" << Result(k,0) 
                        << "\t" << Result(k,1) 
                        << "\t" << Result(k,2) << endl;
  scratchlogfile << "Estimated total offset (l,p): \t"
                 <<  offsetLines << ", " << offsetPixels
                 << "\n*******************************************************************\n";
  scratchlogfile.close();

  ofstream scratchresfile("scratchrescoarse2", ios::out | ios::trunc);
  bk_assert(scratchresfile,"coarsecorrelfft: scratchrescoarse2",__FILE__,__LINE__);
  scratchresfile << "\n\n*******************************************************************"
                 << "\n*_Start_" << processcontrol[pr_i_coarse2]
                 << "\n*******************************************************************"
                 << "\nEstimated translation slave w.r.t. master:"
                 << "\nCoarse_correlation_translation_lines: \t"
                 <<  offsetLines                        // 1 digit after point?
                 << "\nCoarse_correlation_translation_pixels: \t"
                 <<  offsetPixels                       // 1 digit after point?
                 << "\n*******************************************************************"
                 //<< "\n* End_coarse_correlation:_NORMAL"
                 << "\n* End_" << processcontrol[pr_i_coarse2] << "_NORMAL"
                 << "\n*******************************************************************\n";

// ______Tidy up______
  scratchresfile.close();
  INFO << "Individual estimated translations (#, l, p, corr, offl, offp):";
  INFO.print();
  for (k=0; k<Nwin; k++)
    {
    INFO // << "Estimate (#, l, p, corr, offl, offp): "
         << k             << " \t"
         << Minlminp(k,0) << " \t"
         << Minlminp(k,1) << " \t"
         << Result(k,2)   << " \t"
         << Result(k,0)   << " \t"
         << Result(k,1)   << " \t";
    INFO.print();
    }

  INFO << "Estimated translation (l,p): "
       << offsetLines << ", " << offsetPixels;
  INFO.print();
  PROGRESS.print("Coarse coregistration based on correlation finished.");
  } // END coarsecorrelfft



/****************************************************************
 * corrfft                                                      *
 *                                                              *
 * coherence in spectral domain by fft's                        *
 *  uses extension with zeros                                   *
 *  pixel level                                                 *
 *                                                              *
 * input:                                                       *
 *  - Master                                                    *
 *  - Mask (size Master)                                        *
 * output:                                                      *
 *  - coherence value                                           *
 *  - updated offsetL, P                                        *
 *                                                              *
 *    Bert Kampes, 16-Feb-1999                                  *
 *                                                              *
 * note: this routine can be speeded up by removing matrices    *
 * powerma* and by using *= instead of dotmult                  *
 * for now this is not done because it requires only little time*
 * note also that for coarse coregistration division by powers  *
 * is not really required, cross products are good enough.      *
 *    Bert Kampes, 18-Oct-1999                                  *
 ****************************************************************/
real4 corrfft(
         const matrix<real4> &Master,                   // magnitude image
         const matrix<real4> &Mask,                     // magnitude image
         real4 &offsetL,                                // updated
         real4 &offsetP)                                // updated 
  {
  TRACE_FUNCTION("corrfft (BK 18-Oct-1999)");
  const int32 L     = Master.lines();
  const int32 P     = Master.pixels();
  const int32 twoL  = 2*L;
  const int32 twoP  = 2*P;
  const int32 halfL = L/2;
  const int32 halfP = P/2;

  if (L != Mask.lines() || P != Mask.pixels())
    {
    PRINT_ERROR("Mask, Master not same size.")
    throw(input_error);
    }
  if (!(ispower2(L) || ispower2(P)))
    {
    PRINT_ERROR("Mask, Master size not power of 2.")
    throw(input_error);
    }

  // ======Compute powers for submatrices======
  register int32 i;
  register int32 j;

  const complr4 ONE(1.);
  matrix<complr4> Master2(twoL,twoP);                   // init 0
  matrix<complr4> Mask2(twoL,twoP);                     // init 0
  matrix<complr4> blok2(twoL,twoP);                     // init 0

  const real4 meanMaster = mean(Master);
  const real4 meanMask   = mean(Mask);

// ====== Powers, misuse Master2, blok2 ======
// ______ First powers to save a matrix, 3 is minimum ______
  for (i=0; i<L; i++)
    {
    for (j=0; j<P; j++)
      {
      Master2(i,j)   = complr4(sqr(Master(i,j)-meanMaster));    // intensity image
      Mask2(i+L,j+P) = complr4(sqr(Mask(i,j)-meanMask));// only real part
      blok2(i+halfL,j+halfP) = ONE;                     // only real part
      }
    }

  fft2d(Master2);
  fft2d(Mask2);
  fft2d(blok2);

// old way:
//  matrix<complr4> Powermaster = dotmult(conj(Master2),blok2); // twice conj(M2)
//  ifft2d(Powermaster);
//  matrix<complr4> Powermask = dotmult(conj(blok2),Mask2);
//  ifft2d(Powermask);

// new way: test this
  blok2.conj();                                         // conjugated in blok2
  Master2 *= blok2;
  Master2.conj();                                       // M2=original(b2)*conj(m2)
  ifft2d(Master2);                                      // norms master in Master2

  Mask2 *= blok2;                                       // powers in spectral domain
  ifft2d(Mask2);                                        // norms slave in Mask2

// ______ Use real(block2) to store sqrt(norms1*norms2) ______
  for (i=0; i<=L; i++)                          // all shifts
    for (j=0; j<=P; j++)                        // all shifts
      blok2(i,j) = complr4(sqrt(real(Master2(i,j))*real(Mask2(i,j))));


// ====== Cross products covariance master/slave ======
  Master2.clean();                                      // init 0
  Mask2.clean();                                        // init 0
  for (i=0;i<L;i++)
    {
    for (j=0;j<P;j++)
      {
      Master2(i,j)           = complr4(Master(i,j)-meanMaster); // only real
      Mask2(i+halfL,j+halfP) = complr4(Mask(i,j)-meanMask);     // part mag. image
      }
    }


// ======FFT's of master/mask======
// padded with N zeros to prevent periodical convolution
  fft2d(Master2);
  fft2d(Mask2);

// ______ Store in Mask2 crossproducts in spectral/space domain ______
  Master2.conj();
  Mask2 *= Master2;                                     // corr. by zero padding
  ifft2d(Mask2);                                        // space domain (real only)
//  matrix<complr4> Cor = dotmult(conj(Master2),Mask2); // correlation with zero ext.
//  ifft2d(Cor);                                        // space


// ====== Correlation in space domain for all shifts [-N/2,N/2] ======
  real4 coher;
  real4 maxcoher = 0.;
  for (i=0; i<=L; i++)                          // all shifts
    {
    for (j=0; j<=P; j++)                        // all shifts
      {
      //coher = (sqr(real(Cor(i,j)))) /                         // old way
//            ((real(Powermaster(i,j))*real(Powermask(i,j))));
//      coher = (sqr(real(Cor(i,j)))) /                                 // new way
//            ((real(Master2(i,j))*real(Mask2(i,j))));
      coher = real(Mask2(i,j)) / real(blok2(i,j));
      if (coher > maxcoher)
        {
        maxcoher = coher;
        offsetL = -halfL + i;                   // update by reference
        offsetP = -halfP + j;                   // update by reference
        }
      }
    }
  return maxcoher;
  } // END corrfft



/****************************************************************
 *    distributepoints                                          *
 *                                                              *
 * Returns matrix with distributed points in of input.          *
 * First window at (win.linelo,win.pixlo),                      *
 *  divided over wl lines,                                      *
 *  with dp distance in pixel direction.                        *
 *                                                              *
 * input:                                                       *
 *  - number of windows                                         *
 *  - window which should be divided                            *
 * output:                                                      *
 *  - matrix <uint> (NW,3) =(l,p, flagfromdisk==0)              *
 *                                                              *
 *    Bert Kampes, 21-Jan-1999                                  *
 ****************************************************************/
matrix<uint> distributepoints(
        real4 nW,
        const window &win)
  {
  TRACE_FUNCTION("distributepoints (BK 21-Jan-1999)")
  real4 lines  = win.linehi - win.linelo + 1;
  real4 pixels = win.pixhi  - win.pixlo  + 1;

  uint numw = uint(nW);
  matrix<uint> Result(numw,uint(3));
  // ______Distribution for dl=dp______
  real4 wp = sqrt(nW/(lines/pixels));   // wl: #windows in line direction
  real4 wl = nW / wp;                   // wp: #windows in pixel direction

  if (wl < wp)                          // switch wl,wp : later back
    wl = wp; 

  int32 wlint  = int32(floor(wl + .5)); // round largest

  real4 deltal = (lines-1) / (real4(wlint-1));
  int32 totp   = int32(pixels*wlint);
  real4 deltap = (real4(totp-1)) / (real4(nW-1)); 

  real4 p      = -deltap;
  real4 l      = 0.;
  uint lcnt    = 0;

  register int32 i;
  for (i=0; i<nW; i++)
    {
    p += deltap;
    while (floor(p+.5)>=pixels)         // round
      {
      p -= pixels;
      lcnt++;
      }
    l = lcnt * deltal;
    Result(i,0) = uint(floor(l+.5));
    Result(i,1) = uint(floor(p+.5));
    }

// ______Correct distribution to window______
  for (i=0; i<nW; i++)
    {
    Result(i,0) += win.linelo;
    Result(i,1) += win.pixlo;
    }

  return Result;
  } // END distributepoints



/****************************************************************
 *    getoffset                                                 *
 *                                                              *
 * Returns offset in line and pixel direction                   *
 *  based on matrix with estimated offests                      *
 *  by correlation                                              *
 * Checks on consistency, THRESHOLD 0.4 for correlation         *
 *                                                              *
 * input:                                                       *
 *  - matrix<real4> with offsets(l,p),correlation               *
 *  - (offL,offP)                                               *
 * output:                                                      *
 *  - updated (offL,offP)                                       *
 *                                                              *
 * See also Documentation page 6.                               *
 *                                                              *
 *    Bert Kampes, 21-Jan-1999                                  *
 ****************************************************************/
void getoffset(
        const matrix<real4> &Result,
        int32 &offsetLines,
        int32 &offsetPixels)
  {
  TRACE_FUNCTION("getoffset (BK 21-Jan-1999)")
  if (Result.pixels() != 3)
    {
    PRINT_ERROR("code 901: wrong input not 3 width");
    throw(input_error);
    }
  uint nW = Result.lines();
  int32 cnt;
  int32 valueL;
  int32 valueP;
  real4 correl;
  int32 highestcnt    = 0;
  real4 highestcorrel = 0.;
  for (register int32 i=0;i<nW;i++)
    {
    valueL = int32(Result(i,0));
    valueP = int32(Result(i,1));
    correl = Result(i,2);
    if (correl > highestcorrel) 
      highestcorrel = correl;
    cnt = 0;
    for (register int32 j=0;j<nW;j++)
      {
      if (abs(Result(j,0) - valueL) < 2  &&  
          abs(Result(j,1) - valueP) < 2)
        cnt++;
      }
    if (cnt > highestcnt)
      {
      highestcnt = cnt;
      offsetLines  = valueL;                    // Return these
      offsetPixels = valueP;                    // Return these
      }
    }

  // ______ Check result ______
  real4 THRESHOLD = .4;
  if (nW < 6)
    {
    WARNING.print("getoffset: number of windows to estimate offset < 6");
    WARNING.print("(please check bottom of LOGFILE)");
    }
  if (highestcnt < .2*nW)
    {
    WARNING.print("getoffset: estimated offset not consistent with other estimates.");
    WARNING.print("(check bottom of LOGFILE)");
    }
  if (highestcorrel < THRESHOLD)
    {
    WARNING << "getoffset: estimated translation has correlation of: "
         << highestcorrel;
    WARNING.print();
    WARNING.print("(please check bottom of LOGFILE)");
    }
  } // END getoffset



/****************************************************************
 *    finecoreg                                                 *
 *                                                              *
 * computes translation of slave w.r.t. master                  *
 * slave(some point) = master(same point) + trans(l,p) =>       *
 *  trans = slavecoordinates - mastercoordinates                *
 * in NWIN windows.                                             *
 * Then solves polynomial for best transformation to master     *
 * with coregpm routine/step                                    *
 *                                                              *
 * input:                                                       *
 *  -                                                           *
 * output:                                                      *
 *  -                                                           *
 *    Bert Kampes, 12-Dec-1998                                  *
 * distribute points can be with input file as well, besides    *
 * letting Doris randomly distribute npoints.                   *
 *    BK 29-Oct-99                                              *
 ****************************************************************/
void finecoreg(
        const input_fine &fineinput,
        const slcimage   &minfo,
        const slcimage   &sinfo)
  {
  TRACE_FUNCTION("finecoreg (BK 29-Oct-99)")
  char dummyline[ONE27];

  const uint Mfilelines = minfo.currentwindow.lines();
  const uint Sfilelines = sinfo.currentwindow.lines();

  const uint Nwin   = fineinput.Nwin;           // n windows, from file or random
  const int32 initoffsetL = fineinput.initoffsetL;      // initial offset
  const int32 initoffsetP = fineinput.initoffsetP;      // initial offset
  uint MasksizeL    = fineinput.MasksizeL;      // size of correlation window
  uint MasksizeP    = fineinput.MasksizeP;      // size of correlation window

  bool pointsrandom = true;
  if (specified(fineinput.ifpositions))         // filename specified
    pointsrandom = false;                       // only use these points

// ______Correct sizes if in space domain______
  if (fineinput.method == fc_magspace ||
      fineinput.method == fc_cmplxspace)
    {
    MasksizeL += 2*fineinput.AccL;
    MasksizeP += 2*fineinput.AccP;
    }


// ______Corners of slave in master system______
// ______offset = [A](slave system) - [A](master system)______
  const int32 sl0 = sinfo.currentwindow.linelo - initoffsetL;
  const int32 slN = sinfo.currentwindow.linehi - initoffsetL;
  const int32 sp0 = sinfo.currentwindow.pixlo  - initoffsetP;
  const int32 spN = sinfo.currentwindow.pixhi  - initoffsetP;

// ______Corners of overlap master,slave in master system______
//  int32 overl0 = max(int32(minfo.currentwindow.linelo),sl0);
//  int32 overlN = min(int32(minfo.currentwindow.linehi),slN);
//  int32 overp0 = max(int32(minfo.currentwindow.pixlo),sp0);
//  int32 overp0 = min(int32(minfo.currentwindow.pixhi),spN);
 
// ______Corners of useful overlap master,slave in master system______
  const uint FEW=10;
  const uint l0 = max(int32(minfo.currentwindow.linelo),sl0) + FEW;
  const uint lN = min(int32(minfo.currentwindow.linehi),slN) - MasksizeL - FEW;
  const uint p0 = max(int32(minfo.currentwindow.pixlo),sp0)  + FEW;
  const uint pN = min(int32(minfo.currentwindow.pixhi),spN)  - MasksizeP - FEW;
  const window overlap(l0,lN,p0,pN);

// ______ Distribute Nwin points over window, or read from file ______
// ______ Minlminp(i,0): line, (i,1): pixel, (i,2) flagfromdisk ______
  register int32 i;
  matrix<uint> Minlminp;
  if (pointsrandom)                             // no filename specified
    {
    Minlminp = distributepoints(real4(Nwin),overlap);
    }

  else  // read in points (center of windows) from file
    {
    // const int32 numwin = filelines(fineinput.ifpositions); // done in readinput
    Minlminp.resize(Nwin,3);
    //ifstream ifpos(fineinput.ifpositions, ios::in | ios::nocreate);
    ifstream ifpos(fineinput.ifpositions, ios::in);
    bk_assert(ifpos,fineinput.ifpositions,__FILE__,__LINE__);
    uint ll,pp;
    for (i=0; i<Nwin; ++i)
      {
      ifpos >> ll >> pp;
      Minlminp(i,0) = uint(ll -.5*MasksizeL);   // correct for lower left corner
      Minlminp(i,1) = uint(pp -.5*MasksizeP);   // correct for lower left corner
      Minlminp(i,2) = uint(1);                  // flag from file
      ifpos.getline(dummyline,ONE27,'\n');      // goto next line.
      }
    ifpos.close();

// ______ Check last point ivm. EOL after last position in file ______
    if (Minlminp(Nwin-1,0) == Minlminp(Nwin-2,0) &&
        Minlminp(Nwin-1,1) == Minlminp(Nwin-2,1))
      {
      Minlminp(Nwin-1,0) = uint(.5*(lN + l0) + 27);     // random
      Minlminp(Nwin-1,1) = uint(.5*(pN + p0) + 37);     // random
      }

// ______ Check if points are in overlap ______
// ______ no check for uniqueness of points ______
    bool troubleoverlap = false;
    for (i=0; i<Nwin; ++i)
      {
      if (Minlminp(i,0) < l0)
        {
        troubleoverlap=true;
        WARNING << "FINE: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,0) = l0 + l0-Minlminp(i,0);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,0) > lN)
        {
        troubleoverlap=true;
        WARNING << "FINE: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,0) = lN + lN-Minlminp(i,0);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,1) < p0)
        {
        troubleoverlap=true;
        WARNING << "FINE: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,1) = p0 + p0-Minlminp(i,1);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      if (Minlminp(i,1) > pN)
        {
        troubleoverlap=true;
        WARNING << "FINE: point from file: "
             << i+1 << " " << Minlminp(i,0) +.5*MasksizeL << " "
             << Minlminp(i,1) +.5*MasksizeP
             << " outside overlap master, slave. New position: ";
        Minlminp(i,1) = pN + pN-Minlminp(i,1);
        WARNING << Minlminp(i,0) << " " << Minlminp(i,1);
        WARNING.print();
        }
      }
    if (troubleoverlap) // give some additional info
      {
      WARNING << "FINE: there were points from file outside overlap (l0,lN,p0,pN): "
           << l0 << " " << lN << " " << p0 << " " << pN;
      WARNING.print();
      }
    }


  // ______Compute coherence of these points______
  //matrix<complr4> Master(MasksizeL,MasksizeP);        // on file
  //matrix<complr4> Mask(MasksizeL,MasksizeP);  // on file
  matrix<complr4> Master;
  matrix<complr4> Mask;
  matrix<real4>   Result(Nwin,3);               // R(i,0):delta l; 
                                                //  R(i,1):delta p; R(i,2):correl

  // ______ Progress message ______
  int32 tenpercent = int32((Nwin/10)+.5);    // round
  if (tenpercent==0) tenpercent = 1000;
  int32 percent = 0;

  // ====== Compute for all locations ======
  for (i=0;i<Nwin;i++)
    {
    // ______ Give progress message ______
    if (i%tenpercent==0)
      {
      PROGRESS << "FINE: " << setw(3) << percent << "%";
      PROGRESS.print();
      percent += 10;
      }

    // ______Minlminp (lower left corners) of window in master system______
    uint minMwinL = Minlminp(i,0);
    uint minMwinP = Minlminp(i,1);

                                                        // size=masksize
    window master(minMwinL, minMwinL+MasksizeL-1,
                  minMwinP, minMwinP+MasksizeP-1);
    // ______Same points in slave system (disk)______
    window mask(minMwinL+initoffsetL, 
                minMwinL+initoffsetL+MasksizeL-1,
                minMwinP+initoffsetP,
                minMwinP+initoffsetP+MasksizeP-1);

    // ______Read windows from files______
    //minfo.readdata(Master,master);
    //sinfo.readdata(Mask,mask);
    Master = minfo.readdata(master);
    Mask   = sinfo.readdata(mask);


    // ______Coherence______
    // ______update offsetL/P______
    real4 offsetL, offsetP;
    real4 coheren;
    switch (fineinput.method)
      {

//        case fc_cmplxfft:
//      WARNING("THIS METHOD IS NOT OK YET, I RECOMMEND MAGNITUDE.");
//          coheren = coherencefft(fineinput, Master, Mask, offsetL, offsetP);
//      break;
//        case fc_cmplxspace:
//      WARNING("THIS METHOD IS NOT OK YET, I RECOMMEND MAGNITUDE.");
//          coheren = coherencespace(fineinput, Master, Mask, offsetL, offsetP);
//      break;

      case fc_magfft:
        coheren = coherencefft(fineinput, Master, Mask, offsetL, offsetP);
        break;
      case fc_magspace:
        coheren = coherencespace(fineinput, Master, Mask, offsetL, offsetP);
        break;
      default:
        PRINT_ERROR("unknown method for fine coregistration.")
        throw(unhandled_case_error);
      } // switch method


    Result(i,0) = offsetL + initoffsetL;
    Result(i,1) = offsetP + initoffsetP;
    Result(i,2) = coheren;
    } // for nwin

// ______Position approx. with respect to center of window______
// ______correct position array for center instead of lower left______
  for (i=0; i<Nwin; i++)
    {
    Minlminp(i,0) += uint(.5 * MasksizeL);
    Minlminp(i,1) += uint(.5 * MasksizeP);
    }


// ______Write to files______
  ofstream scratchlogfile("scratchlogfine", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"finecoreg: scratchlogfine",__FILE__,__LINE__);
  scratchlogfile << "\n\n*******************************************************************"
                 << "\n* FINE_COREGISTRATION"
                 << "\n*******************************************************************"
                 << "\nNumber of correlation windows: \t"
                 <<  Nwin
                 << "\nwindow size (l,p):             \t"
                 <<  MasksizeL << ", " << MasksizeP
                 << "\nInitial offsets:               \t"
                 <<  initoffsetL << ", " << initoffsetP
                 << "\nOversampling factor:           \t"
                 <<  fineinput.osfactor
                 << "\n\nNumber \tposl \tposp \toffsetl offsetp\tcorrelation\n";
  for (i=0;i<Nwin;i++)
    scratchlogfile
      << setiosflags(ios::fixed)
      << setiosflags(ios::showpoint)
      << setiosflags(ios::right)
      << setw(8)
      << setprecision(0)
      << i << " "
      << Minlminp(i,0) << " "
      << Minlminp(i,1) << " "
      << setprecision(3)
      << Result(i,0) << " "
      << Result(i,1) << " "
      << setprecision(2)
      << Result(i,2) << endl;
  scratchlogfile << "\n*******************************************************************\n";
  scratchlogfile.close();

  ofstream scratchresfile("scratchresfine", ios::out | ios::trunc);
  bk_assert(scratchresfile,"finecoreg: scratchresfine",__FILE__,__LINE__);

  scratchresfile 
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_i_fine]
    << "\n*******************************************************************"
    << "\nOversampling factor: \t\t"
    <<  fineinput.osfactor
    << "\nNumber_of_correlation_windows: \t"
    <<  Nwin
    << "\nNumber \tposL \tposP \toffsetL offsetP\tcorrelation\n";
  scratchresfile.close();

  FILE *resfile;
  resfile=fopen("scratchresfine","a");
  for (i=0; i<Nwin; i++)
    fprintf(resfile,"%4.0f %5.0f %5.0f %# 9.2f %# 9.2f %# 6.2f\n",
            real4(i), real4(Minlminp(i,0)), real4(Minlminp(i,1)),
            Result(i,0), Result(i,1), Result(i,2));
  fprintf(resfile,
  "\n*******************************************************************");
  fprintf(resfile,"%s%s%s",
  "\n* End_", processcontrol[pr_i_fine], "_NORMAL");
  fprintf(resfile,
  "\n*******************************************************************\n");


// ______Tidy up______
  fclose(resfile);
  PROGRESS.print("Fine coregistration finished.");
  } // END finecoreg




/****************************************************************
 * coherencefft                                                 *
 *                                                              *
 * coherence in spectral domain by fft's based on magnitude     *
 *  uses extension with zeros                                   *
 *                                                              *
 * input:                                                       *
 *  - Master                                                    *
 *  - Mask (size Master)                                        *
 * output:                                                      *
 *  - coherence value [-1 1]                                    *
 *  - updated offsetL, P                                        *
 *                                                              *
 *    Bert Kampes, 03-Feb-1999                                  *
 * bugfix? streamlined, based on magnitude forced               *
 *    Bert Kampes, 16-Nov-1999                                  *
 ****************************************************************/
real4 coherencefft(
        const input_fine &fineinput,
        const matrix<complr4> &Master,          // data
        const matrix<complr4> &Mask,            // data
        real4 &offsetL,                         // returned
        real4 &offsetP)                         // returned
  {
  TRACE_FUNCTION("coherencefft (BK 16-Nov-1999)")
  const uint factor = fineinput.osfactor;
  const uint AccL   = fineinput.AccL;
  const uint AccP   = fineinput.AccP;
  const int32 L     = Master.lines();
  const int32 P     = Master.pixels();
  const int32 twoL  = 2*L;
  const int32 twoP  = 2*P;
  const int32 halfL = L/2;
  const int32 halfP = P/2;

  if (fineinput.method != fc_magfft)
    {
    PRINT_ERROR("impossible error, wrong input method in fft.")
    throw(input_error);
    }
  if (!ispower2(factor))
    {
    PRINT_ERROR("coherence factor not power of 2")
    throw(input_error);
    }
  if (L != Mask.lines() || P != Mask.pixels())
    {
    PRINT_ERROR("Mask, Master not same size.")
    throw(input_error);
    }
  if (!(ispower2(L) || ispower2(P)))
    {
    PRINT_ERROR("Mask, Master size not power of 2.")
    throw(input_error);
    }

// ______ Zero mean magnitude images ______
  matrix<real4> magMaster = magnitude(Master);
  matrix<real4> magMask   = magnitude(Mask);
  magMaster              -= mean(magMaster);
  magMask                -= mean(magMask);


// ====== FFT's of master/mask ======
// ______ Pad with N zeros to prevent periodical convolution ______
  window windef(0,0,0,0);                       // defaults to total matrix
  window win1(0,L-1,0,P-1);     
  window win2(halfL,halfL+L-1,halfP,halfP+P-1);

  matrix<complr4> Master2(twoL,twoP);                   // initial 0
  matrix<complr4> Mask2(twoL,twoP);                     // initial 0
  Master2.setdata(win1,mat2cr4(magMaster),windef);      // zero mean magnitude
  Mask2.setdata(win2,mat2cr4(magMask),windef);          // zero mean magnitude

  fft2d(Master2);
  fft2d(Mask2);

// ______ Crossproducts in spectral/space domain ______
// ______ Use Mask2 to store cross products temporarly ______
  Master2.conj();
  Mask2 *= Master2;                                     // corr = conj(M).*S
  ifft2d(Mask2);                                        // cross prod. in space

// ______ Resize to shifts [-AccL,+AccL) ______
  window wintmp(halfL-AccL,halfL+AccL-1,halfP-AccP,halfP+AccP-1);
  matrix<complr4> TMP(wintmp,Mask2);
  matrix<real4> Covar = real(TMP);              // imag == 0 (?)


// ====== Compute norms, zero padded matrices ======
  matrix<complr4> blok(L,P);
  const complr4 ONE(1.);
  blok.setdata(ONE);                                    // only real part

  Master2.clean();                                      // reset to zeros
  Mask2.clean();                                        // reset to zeros
  Master2.setdata(win1,mat2cr4(sqr(magMaster)),windef);// use Master2 for intensity
  Mask2.setdata(win2,blok,windef);                      // use Mask2 for padded Block

  fft2d(Master2);                                       // (intensity of master)
  fft2d(Mask2);                                         // (block)

  Mask2.conj();                                         // conj(block)
  Master2 *= Mask2;
  Master2.conj();                               // Master2 == conj(Master)*block
  ifft2d(Master2);
// ______ Master2 now contains norms of master image in space domain ______
// ______ Resize to shifts [-AccL,+AccL) ______
  TMP.setdata(Master2,wintmp);                          // fill TMP
  matrix<real4> pmaster = real(TMP);                    // norms in pmaster

// ====== Now norms for slave image ======
  Master2.clean();                                      // reset to zeros
  window win5(L,twoL-1,P,twoP-1);
  Master2.setdata(win5,mat2cr4(sqr(magMask)),windef);
  fft2d(Master2);                                       // (intensity of slave)
  Master2 *= Mask2;                             // Master2 == conj(block)*Slave
  ifft2d(Master2);
// ______ Master2 now contains norms of slave image in space domain ______
// ______Resize to shifts [-AccL,+AccL)______
  TMP.setdata(Master2,wintmp);                          // fill TMP
  matrix<real4> pmask = real(TMP);                      // norms in pmask

// ====== Correlation in space domain ======
  pmaster *= pmask;
  Covar /= sqrt(pmaster);


// ====== Estimate shift by oversampling ======
  uint offL;
  uint offP;
  // old: matrix<real4> coher8 = oversample(Covar,factor);
  matrix<real4> coher8 = oversample(Covar,factor,factor);
  real4 maxcor = max(coher8,offL,offP);
  offsetL = -real4(AccL) + offL/real4(factor);           // update by reference
  offsetP = -real4(AccP) + offP/real4(factor);           // update by reference

  return maxcor;
  } // END coherence fft



/****************************************************************
 * coherencespace                                               *
 *                                                              *
 * coherence in space domain based on magnitude                 *
 *  uses extension with zeros                                   *
 *                                                              *
 * input:                                                       *
 *  - Master                                                    *
 *  - Mask (size Master)                                        *
 * output:                                                      *
 *  - coherence value                                           *
 *  - updated offsetL, P                                        *
 *                                                              *
 *    Bert Kampes, 03-Feb-1999                                  *
 ****************************************************************/
real4 coherencespace(
        const input_fine &fineinput,
        const matrix<complr4> &Master,          // complex data
        const matrix<complr4> &Mask,            // complex data
        real4 &offsetL,                         // returned
        real4 &offsetP)                         // returned
  {
  TRACE_FUNCTION("coherencespace (BK 03-Feb-1999)")
  const int32 L         = Master.lines();
  const int32 P         = Master.pixels();

  const int32 AccL      = fineinput.AccL;
  const int32 AccP      = fineinput.AccP;
  const uint factor     = fineinput.osfactor;

  // ______Select parts of Master/slave______
  const int32 MasksizeL = L - 2*AccL;
  const int32 MasksizeP = P - 2*AccP;

  if (!ispower2(AccL) || !ispower2(AccP))
    {
    PRINT_ERROR("AccL should be power of 2 for oversampling.")
    throw(input_error);
    }
  if (MasksizeL < 4 || MasksizeP < 4)
    {
    PRINT_ERROR("Correlationwindow size too small (<4; size= FC_winsize-2*FC_Acc).")
    throw(input_error);
    }

  // ______Shift center of Slave over Master______
  window winmask(AccL, AccL+MasksizeL-1,
                 AccP, AccP+MasksizeP-1);

  matrix<real4> coher(2*AccL,2*AccP);           // store result
                                                // 1st element: shift==AccL
  window windef(0, 0, 0, 0);                    // defaults to total

  switch (fineinput.method)
    {
    case fc_cmplxspace:
      {
      PRINT_ERROR("not implemented in v1.0")
      throw(unhandled_case_error);
      break;
      }

//    if (fineinput.method == fc_cmplxspace)
//      {
//  DEBUG("mean should be subtracted.        e");
//      matrix<complr4> Mask2(winmask,Mask);            // construct as part 
//      real4 normmask = norm2(Mask2);
//    
//      matrix<complr4> Master2(MasksizeL, MasksizeP); 
//      window winmaster;
//      
//      for (register int32 i=0;i<2*AccL;i++)
//        {
//        winmaster.linelo = i;
//        winmaster.linehi = i+MasksizeL-1;
//        for (register int32 j=0;j<2*AccP;j++)
//          {
//          winmaster.pixlo = j;
//          winmaster.pixhi = j+MasksizeP-1;
//          Master2.setdata(windef,Master,winmaster);
//       
//  // ______Compute coherence for this position______
//          complr4 cohs1s2(0.,0.);
//          real4 cohs1s1 = 0.;
//  DEBUG("nog multilooked coherence-ref.phase");
//          for (register int32 k=0;k<MasksizeL;k++)
//            {
//            for (register int32 l=0;l<MasksizeP;l++)
//              {
//          cohs1s2 += (Master2(k,l) * conj(Mask2(k,l)));
//          cohs1s1 += norm(Master2(k,l));
//              }
//            }
//          coher(i,j) = norm(cohs1s2) / (cohs1s1 * normmask);
//          }
//        }
//      }

    case fc_magspace:
      {
      matrix<real4> magMask = magnitude(Mask);          // magnitude
      magMask -= mean(magMask);                         // subtract mean
      //matrix<real4> Mask2(winmask,magnitude(Mask));   // construct as part 
      matrix<real4> Mask2(winmask,magMask);             // construct as part 
      real4 normmask = norm2(Mask2);


      matrix<real4> Master2(MasksizeL, MasksizeP); 
      matrix<real4> magMaster=magnitude(Master);
      magMaster -= mean(magMaster);
      window winmaster;


      for (register int32 i=0;i<2*AccL;i++)
        {
        winmaster.linelo = i;
        winmaster.linehi = i+MasksizeL-1;
        for (register int32 j=0;j<2*AccP;j++)
          {
          winmaster.pixlo = j;
          winmaster.pixhi = j+MasksizeP-1;
          Master2.setdata(windef,magMaster,winmaster);
  
// ______Coherence for this position______
          real4 cohs1s2 = 0.;
          real4 cohs1s1 = 0.;
          for (register int32 k=0;k<MasksizeL;k++)
            {
            for (register int32 l=0;l<MasksizeP;l++)
              {
              cohs1s2 += (Master2(k,l) * Mask2(k,l));
              cohs1s1 += sqr(Master2(k,l));
              }
            }
          coher(i,j) = cohs1s2 / sqrt(cohs1s1 * normmask);      // [-1 1]
          //coher(i,j) = sqr(cohs1s2) / (cohs1s1 * normmask);
          }
        }
      break;
      }
   
    default:
      PRINT_ERROR("unknown method")
      throw(unhandled_case_error);
    } // switch method


// ______ Correlation in space domain ______
  // old matrix<real4> coher8 = oversample(coher,factor);
  matrix<real4> coher8 = oversample(coher,factor,factor);

  uint offL;
  uint offP;
  real4 maxcor = max(coher8,offL,offP);
  offsetL = AccL - offL/real4(factor);           // update by reference
  offsetP = AccP - offP/real4(factor);           // update by reference

  return maxcor;
  } // END coherencespace



/****************************************************************
 * coregpm                                                      *
 *                                                              *
 * Compute coregistration parameters (least squares)            *
 *                                                              *
 * input:                                                       *
 *  - Position of windows                                       *
 *  - Computed offsets                                          *
 *                                                              *
 * output:                                                      *
 *  - coregistration parameters to file                         *
 *    (wrt. normalized master grid)                             *
 *                                                              *
 *    Bert Kampes, 22-Feb-1999                                  *
 *    Bert Kampes, 26-Oct-1999 normalized coordinates           *
 * changed matrices from real4 to real8,                        *
 #%// BK 22-Mar-2001                                            *
 ****************************************************************/
void coregpm(
        const window &originalmaster,           // normalization factors
        const char* i_resfile,
        const input_coregpm &coregpminput,
        uint oversamplingsfactorfine)
  {
  TRACE_FUNCTION("coregpm (BK 26-Oct-1999)")
  // ______Names of variables in this routine______
  // unknowns:          x
  // solution unknowns: placed in rhsL/P
  // observations:      y
  // covariance obs.:   Qy_1
  // designmatrix:      A
  // normalmatrix:      N = At. Qy-1. A
  // covariance unkn.:  N_1 (Qx_hat, inverse normalmatrix)
  // estimates:         *_hat

  const real4 THRESHOLD = coregpminput.threshold;       // threshold ...
  const int32 DEGREE    = coregpminput.degree;          // degree of polynomial
  const int32 MAX_ITERATIONS = coregpminput.maxiter;    // max. of pnts to remove
  const real4 CRIT_VALUE = coregpminput.k_alpha;        // crit. value outlier removal
  const int32 Nunk      = Ncoeffs(DEGREE);              // Number of unknowns/direction

  // ______ Normalize data for polynomial ______
  const real8 minL = originalmaster.linelo;
  const real8 maxL = originalmaster.linehi;
  const real8 minP = originalmaster.pixlo;
  const real8 maxP = originalmaster.pixhi;

  // ______ A priori sigma of  offset ______
  // _____ oversampling factor is bin in which maximum can be found _____
  // _____ ovsf=16-->apriorisigma=0.03
  //const real4 SIGMA   = .5 * (1. / (real4(oversamplingsfactorfine)));
  const real4 ACCURACY = .5 * (1. / (real4(oversamplingsfactorfine)));
  // but we need coreg accuracy of 0.1 pixel about.  therefore use a priori
  // based on experience here, and different for azimuth and range
  // this also helps our automated outlier detection and testing hopefully.
  // BK 15-Apr-2003
  //const real4 SIGMAL   = 0.075;// hmmm. ?? sigma in pixels
  const real4 SIGMAL   = 0.15;// hmmm. ?? sigma in pixels
  const real4 SIGMAP   = SIGMAL/1.5;// seems pixel is 2x better ???
  //const real4 SIGMA_1L = 1./SIGMAL;
  //const real4 SIGMA_1P = 1./SIGMAP;
  //const real4 SIGMA2L  = sqr(SIGMAL);
  //const real4 SIGMA2P  = sqr(SIGMAP);
  DEBUG.print("using a better sigma in range, since seems that can be estimated better");
  INFO << "a priori std.dev offset vectors line direction [samples]:  " << SIGMAL;
  INFO.print();
  INFO << "a priori std.dev offset vectors pixel direction [samples]: " << SIGMAP;
  INFO.print();

  // ______ Find #points > threshold ______
  matrix<real4> Data = getofffile(i_resfile, THRESHOLD);
  // ______Data contains:______
  // Data(i,0) = winnumber; Data(i,1) = posL; Data(i,2) = posP; 
  // Data(i,3) = offL;      Data(i,4) = offP; Data(i,5) = corr;


  int32 ITERATION  = 0;
  int32 DONE       = 0;
  // sqr: level significance: alpha=0.001; power of test: gamma=0.80
  //real4 CRIT_VALUE = sqrt(3.29);
  INFO << "Critical value for outlier test: " << CRIT_VALUE;
  INFO.print();
  uint winL = 0;// window number to be removed
  uint winP = 0;// window number of largest w -test in range
  matrix<real8> eL_hat;
  matrix<real8> eP_hat;
  matrix<real8> wtestL;
  matrix<real8> wtestP;
  matrix<real8> rhsL;
  matrix<real8> rhsP;
  matrix<real8> Qx_hat;
  real8 maxdev = 0.0;
  real8 overallmodeltestL = 0.0;
  real8 overallmodeltestP = 0.0;
  real8 maxwL;
  real8 maxwP;
  register int32 i,j,k,index;
  while (DONE != 1)
    {
    PROGRESS << "Start iteration " << ITERATION;
    PROGRESS.print();
    // ______ Remove identified outlier from previous estimation ______
    if (ITERATION != 0)
      {
      matrix<real4> tmp_DATA = Data; //(remove_observation_i,*);
      Data.resize(Data.lines()-1, Data.pixels());
      j = 0;// counter over reduced obs.vector
      for (i=0; i<tmp_DATA.lines(); i++)// counter over original window numbers
        {
        if (i != winL)// do not copy the one to be removed.
          {
          Data.setrow(j,tmp_DATA.getrow(i));// copy back without removed obs.
          j++;// fill next row of Data
          }
        else
          {
          INFO << "Removing observation " << i << " from observation vector.";
          INFO.print();
          }
        }
      }

    // ______Check redundancy______
    int32 Nobs = Data.lines();                          // Number of points > threshold
    if (Nobs < Nunk)
      {
      PRINT_ERROR("coregpm: Number of windows > threshold is smaller than parameters solved for.")
      throw(input_error);
      }
  
    // ______Set up system of equations______
    // ______Order unknowns: A00 A10 A01 A20 A11 A02 A30 A21 A12 A03 for degree=3______
    matrix<real8> yL(Nobs,1);                   // observation
    matrix<real8> yP(Nobs,1);                   // observation
    matrix<real8> A(Nobs,Nunk);                 // designmatrix
    matrix<real8> Qy_1(Nobs,1);                 // a priori covariance matrix (diag)
  
    // ______ Normalize data for polynomial ______
    INFO << "coregpm: polynomial normalized by factors: "
         << minL << " " << maxL << " " << minP << " " << maxP << " to [-2,2]";
    INFO.print();
  
    // ______Fill matrices______
    DEBUG.print("Setting up design matrix for LS adjustment");
    for (i=0; i<Nobs; i++)
      {
      real8 posL = normalize(real8(Data(i,1)),minL,maxL);
      real8 posP = normalize(real8(Data(i,2)),minP,maxP);
      yL(i,0)    = real8(Data(i,3));
      yP(i,0)    = real8(Data(i,4));
      DEBUG << "coregpm: (" << posL << ", "<< posP << "): yL=" 
            << yL(i,0) << " yP=" << yP(i,0);
      DEBUG.print();
      // ______Set up designmatrix______
      index = 0;
      for (j=0; j<=DEGREE; j++)
        {
        for (k=0; k<=j; k++)
          {
          A(i,index) = pow(posL,real8(j-k)) * pow(posP,real8(k));
          index++;
          }
        }
      }
  

    // ______Weight matrix data______
    DEBUG.print("Setting up covariance matrix for LS adjustment");
    switch(coregpminput.weightflag)
      {
      case 0:
        //Qy_1 = 1.0;
        for (i=0; i<Nobs; i++)
          Qy_1(i,0)     = real8(1.0); 
        break;
      case 1:
        //Qy_1 = Data(*,5)
        DEBUG.print("Using sqrt(coherence) as weights.");
        for (i=0; i<Nobs; i++)
          Qy_1(i,0)       = real8(Data(i,5));     // more weight to higher correlation
        break;
      case 2:
        DEBUG.print("Using coherence as weights.");
        for (i=0; i<Nobs; i++)
          Qy_1(i,0)       = real8(Data(i,5))*real8(Data(i,5));// more weight to higher correlation
        break;
      default:
        PRINT_ERROR("Panic, NOT possible with checked input.")
        throw(unhandled_case_error);
      }
    // ______ Normalize weights to avoid influence on estimated var.factor ______
    INFO.print("Normalizing covariance matrix for LS estimation.");
    Qy_1 = Qy_1 / mean(Qy_1);// normalize weights (for tests!)

 
    // ______Compute Normalmatrix, rghthandside______
    matrix<real8> N    = matTxmat(A,diagxmat(Qy_1,A));
    //matrix<real8> rhsL = matTxmat(A,diagxmat(Qy_1,yL));
    //matrix<real8> rhsP = matTxmat(A,diagxmat(Qy_1,yP));
    //matrix<real8> Qx_hat = N;
    rhsL = matTxmat(A,diagxmat(Qy_1,yL));
    rhsP = matTxmat(A,diagxmat(Qy_1,yP));
    Qx_hat = N;
    // ______Compute solution______
    choles(Qx_hat);             // Cholesky factorisation normalmatrix
    solvechol(Qx_hat,rhsL);     // Solution unknowns in rhs
    solvechol(Qx_hat,rhsP);     // Solution unknowns in rhs
    invertchol(Qx_hat);         // Covariance matrix of unknowns
    // ______Test inverse______
    for (i=0; i<Qx_hat.lines(); i++)
      for (j=0; j<i; j++)
        Qx_hat(j,i) = Qx_hat(i,j);// repair Qx
    maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
    INFO << "coregpm: max(abs(N*inv(N)-I)) = " << maxdev;
    INFO.print();
    /*
    matrix<real8> unity = N * Qx_hat;
    maxdev = abs(unity(0,0)-1);
    for (i=1; i<unity.lines(); i++)
      {
      if (abs(unity(i,i)-1) > maxdev)
        maxdev = abs(unity(i,i)-1);
      for (j=0; j<i-1; j++)
        {
        if (abs(unity(i,j)) > maxdev)
          maxdev = abs(unity(i,j));
        if (abs(unity(j,i)) > maxdev)
          maxdev = abs(unity(j,i));
        }
      }
    */
    // ___ use trace buffer to store string, remember to rewind it ___
    if (maxdev > .01) 
      {
      ERROR << "coregpm: maximum deviation N*inv(N) from unity = " << maxdev
            << ". This is larger than .01";
      ERROR.print(ERROR.get_str());
      throw(some_error);
      }
    else if (maxdev > .001) 
      {
      WARNING << "coregpm: maximum deviation N*inv(N) from unity = " << maxdev
              << ". This is between 0.01 and 0.001";
      WARNING.print();
      }


    // ______Some other stuff, scale is ok______
    matrix<real8> Qy_hat        = A * (matxmatT(Qx_hat,A));
    matrix<real8> yL_hat        = A * rhsL;
    matrix<real8> yP_hat        = A * rhsP;
    //matrix<real8> eL_hat      = yL - yL_hat;
    //matrix<real8> eP_hat      = yP - yP_hat;
    eL_hat      = yL - yL_hat;
    eP_hat      = yP - yP_hat;
    //  matrix<real4> Qe_hat    = Qy - Qy_hat;
    matrix<real8> Qe_hat = -Qy_hat;
    for (i=0; i<Nobs; i++)
      Qe_hat(i,i) += (1. / Qy_1(i,0));
  
  
    // ______Overall model test (variance factor)______
    overallmodeltestL = 0.;
    overallmodeltestP = 0.;
    for (i=0; i<Nobs; i++)
      {
      overallmodeltestL += sqr(eL_hat(i,0))*Qy_1(i,0);
      overallmodeltestP += sqr(eP_hat(i,0))*Qy_1(i,0);
      }
    overallmodeltestL = (overallmodeltestL/sqr(SIGMAL)) /(Nobs-Nunk);// this is sigma hat!
    overallmodeltestP = (overallmodeltestP/sqr(SIGMAP)) /(Nobs-Nunk);// not OMT!
  

    INFO << "coregpm: overallmodeltest Lines = " << overallmodeltestL;
    INFO.print();
    INFO << "coregpm: overallmodeltest Pixels = " << overallmodeltestP;
    INFO.print();
  
  
    // ______Datasnooping, assume Qy diag______
    //matrix<real8> wtestL(Nobs,1);
    //matrix<real8> wtestP(Nobs,1);
    wtestL.resize(Nobs,1);
    wtestP.resize(Nobs,1);
  
    for (i=0; i<Nobs; i++)
      {
      wtestL(i,0) = eL_hat(i,0) / (sqrt(Qe_hat(i,i))*SIGMAL);// computed excl.var.factor
      wtestP(i,0) = eP_hat(i,0) / (sqrt(Qe_hat(i,i))*SIGMAP);
      }

    uint dumm = 0;
    maxwL     = max(abs(wtestL),winL,dumm);     // returns winL
    maxwP     = max(abs(wtestP),winP,dumm);     // returns winP
    INFO << "maximum wtest statistic azimuth = " << maxwL
         << " for window number: " 
         <<  Data(winL,0);
    INFO.print();
    INFO << "maximum wtest statistic range   = " << maxwP
         << " for window number: " 
         <<  Data(winP,0);
    INFO.print();
    //if (maxwP>maxwL)
    //  {
    //  winL  = winP;// winL points to rejected window number
    //  //maxwL = maxwP;
    //  }
    // --- use summed wtest for outlier detection ---
    // #%// BK 21-Oct-2003
    matrix<real8> wtestsum = sqr(wtestL)+sqr(wtestP);// (Nobs,1)
    real8 maxwsum = max(wtestsum,winL,dumm);// idx to remove
    INFO << "Detected outlier:  summed sqr.wtest = " << maxwsum
         << "; observation: " << winL
         << "; window number: " 
         <<  Data(winL,0);
    INFO.print();


    // ______ Test if we are done yet ______
    if (Nobs <= Nunk)
      {
      WARNING.print("NO redundancy!  Exiting iterations.");
      DONE = 1;// cannot remove more than this
      }
// seems something fishy here...
//    if (max(overallmodeltestL,overallmodeltestP) < 1.0)
//      {
//      INFO.print("OMTs accepted, not iterating anymore (final solution reached).");
//      DONE = 1;// ok (?).
//      }
    if (max(maxwL,maxwP) <= CRIT_VALUE)// all tests accepted?
      {
      INFO.print("All outlier tests accepted! (final solution computed)");
      DONE = 1;// yeah!
      }
    if (ITERATION >= MAX_ITERATIONS)
      {
      INFO.print("max. number of iterations reached (exiting loop).");
      DONE = 1;// we reached max. (or no max_iter specified)
      }

    // ______ Only warn if last iteration has been done ______
    if (DONE == 1) 
      {
      // ___ use trace buffer to store string, remember to rewind it ___
      if (overallmodeltestL > 10)
        {
        WARNING << "coregpm: overallmodeltest Lines = " << overallmodeltestL << ends;
        WARNING.print();
        WARNING << " is larger than 10. (Suggest model or a priori sigma not correct.)";
        WARNING.print();
        }
      // ___ use trace buffer to store string, remember to rewind it ___
      if (overallmodeltestP > 10)
        {
        WARNING << "coregpm: overallmodeltest Pixels = " << overallmodeltestP;
        WARNING.print();
        WARNING << " is larger than 10.\n(suggests a priori sigma not correct.)";
        WARNING.print();
        }
      // if a priori sigma is correct, max wtest should be something like 1.96
      if (max(maxwL,maxwP)>200.0)
        {
        WARNING << "Recommendation: remove window number: " << Data(winL,0)
               << " and re-run step COREGPM.  max. wtest is: "
               <<  max(maxwL,maxwP) << ".";
        WARNING.print();
        }
      // ______Test of Jaron Samson's thesis: not ok/ depends on SIGMA ...______
      real8 expected = 1. / sqr(28);    //784.;
      //real4 expected = 1. / 400.;             // WRONG, 1./sqr(28) ???
      expected /= sqr(SIGMAL);                  // correct for variance factor
      for (i=0; i<Nobs; i++)
        if (Qy_hat(i,i) > expected)
          {
          WARNING << "coregpm: Qy_hat larger than it should for window number: "
               << Data(i,0);
          WARNING.print();
          }
      }// Only warn when done iterating.


    ITERATION++;// update counter here!
    }// iterations remove outliers
  


  // ______ Create dump file for making plots ______
  ofstream cpmdata("CPM_Data", ios::out | ios::trunc);
  bk_assert(cpmdata,"coregpm: CPM_DATA",__FILE__,__LINE__);
  cpmdata        << "File: CPM_Data"
                 << "\nThis file contains information on the least squares"
                 << "\n estimation of the coregistration parameters."
                 << "\nThis info is used in the plotcmp script."
                 << "\nThere are 10 columns with:"
                 << "\nWindow number, position L, position P, "
                 << "\n offsetL (observation), offsetP (observation), correlation,"
                 << "\n estimated errorL, errorP, w-test statistics for L, P."
                 << "\nwin   posL  posP      offL      offP  corr      eL     eP  wtstL  wtstP"
                 << "\n------------------------------------------------------------\n";
  cpmdata.close();

  // ______ Only way to format in c++ since stupid iomanip dont work? ______
  FILE *cpm;
  cpm=fopen("CPM_Data","a");
  //for (i=0; i<Nobs; i++)
  for (i=0; i<Data.lines(); i++)
    fprintf(cpm,
    "%4.0f %5.0f %5.0f %# 9.2f %# 9.2f %# 6.2f %6.2f %6.2f %6.2f %6.2f\n",
            Data(i,0), Data(i,1), Data(i,2), 
            Data(i,3), Data(i,4), Data(i,5),
            eL_hat(i,0), eP_hat(i,0), 
            abs(wtestL(i,0)), abs(wtestP(i,0)));
  fclose(cpm);


  // ====== Write results to scratch files ======
  ofstream scratchlogfile("scratchlogcpm", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"coregpm: scratchlogcpm",__FILE__,__LINE__);
  scratchlogfile
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_i_coregpm]
    << "\n*******************************************************************"
    << "\nA polynomial model is weighted least squares estimated"
    << "\nfor azimuth and range through the FINE offset vectors."  
    << "\nThe number of coefficient are the unknowns, the number of"
    << "\nobservations are the offset vectors above the THRESHOLD"
    << "\nspecified in the input file.  To estimate the unknowns, at"
    << "\nleast the number of observations must equal the number of unknowns."
    << "\nIf there are more observations, we can statistically test"
    << "\nwhether the observations fit the model, and whether there are"
    << "\noutliers in the observations, which we like to remove."
    << "\nThe overall model test does the first.  Wtest the second."
    << "\nWe advice to remove some bad estimated offsets by hand based"
    << "\nthe largest w-test, and to iterate running this step until"
    << "\nno outlier is identified anymore.  A great tool is plotting"
    << "\nthe observations and errors, which can be done with the utility"
    << "\nscripts provided by Doris (calls to GMT)."
    << "\nAlso see any book on LS methods."
    << "\n\nDegree of model:\t\t\t\t" 
    << DEGREE
    << "\nThreshold on data (correlation):\t\t\t" 
    << THRESHOLD
    << "\nOversmaplings factor used in fine:           \t"
    << oversamplingsfactorfine
    << "\nThis means maximum can be found within [samples]: \t"
    << ACCURACY
    << "\nA priori sigma azimuth (based on experience): \t"
    << SIGMAL
    << "\nA priori sigma range (based on experience): \t"
    << SIGMAP
    << "\nNumber of observations: \t\t\t" 
    << Data.lines()
    << "\nNumber of rejected observations: \t\t\t" 
    << ITERATION
    << "\nNumber of unknowns: \t\t\t\t" 
    << Nunk
    << "\nOverall model test in Azimuth direction: \t" 
    << overallmodeltestL
    << "\nOverall model test in Range direction: \t\t"
    << overallmodeltestP
    << "\nLargest w test statistic in Azimuth direction: \t"
    << maxwL
    << "\n  for window number: \t\t\t\t"
    <<  Data(winL,0)
    << "\nLargest w test statistic in Range direction: \t" 
    << maxwP
    << "\n  for window number: \t\t\t\t"
    <<  Data(winP,0)
    << "\nMaximum deviation from unity Normalmatrix*Covar(unknowns): \t"
    << maxdev
    << "\nEstimated parameters in Azimuth direction"
    << "\nx_hat \tstd"
    << "\n(a00 | a10 a01 | a20 a11 a02 | a30 a21 a12 a03 | ...)\n";
  for (i=0; i<Nunk; i++)
    scratchlogfile
      << setiosflags(ios::fixed)
      << setiosflags(ios::showpoint)
      << setiosflags(ios::right)
      << setw(8) << setprecision(4)
      <<  rhsL(i,0) << " \t" << sqrt(Qx_hat(i,i)) << endl;
  scratchlogfile << "\nEstimated parameters in Range direction"
                 << "\n(b00 | b10 b01 | b20 b11 b02 | b30 b21 b12 b03 | ...)\n";
  for (i=0; i<Nunk; i++)
    scratchlogfile
      << setiosflags(ios::fixed)
      << setiosflags(ios::showpoint)
      << setiosflags(ios::right)
      << setw(8) << setprecision(4)
      <<  rhsP(i,0) << " \t" << Qx_hat(i,i) << endl;

  scratchlogfile << "\nCovariance matrix estimated parameters:"
                 << "\n---------------------------------------\n";
  for (i=0; i<Nunk; i++)
    {
    for (j=0; j<Nunk; j++)
      {
      scratchlogfile
      << setiosflags(ios::fixed)
      << setiosflags(ios::showpoint)
      << setiosflags(ios::right)
      << setw(8) << setprecision(4)
      << Qx_hat(i,j) << " ";
      }
    scratchlogfile << endl;
    }

  scratchlogfile << "\n*******************************************************************\n";
  scratchlogfile.close();


  ofstream scratchresfile("scratchrescpm", ios::out | ios::trunc);
  bk_assert(scratchresfile,"coregpm: scratchrescpm",__FILE__,__LINE__);

  scratchresfile.setf(ios::scientific, ios::floatfield);
  scratchresfile.setf(ios::right, ios::adjustfield);
  scratchresfile.precision(8);
  scratchresfile.width(18);

  scratchresfile
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_i_coregpm]
    << "\n*******************************************************************"
    << "\nDegree_cpm:\t" << DEGREE
    << "\nEstimated_coefficientsL:\n";
  int32 coeffL = 0;
  int32 coeffP = 0;
  for (i=0; i<Nunk; i++)
    {
    if (rhsL(i,0) < 0.)
      scratchresfile <<         rhsL(i,0);
    else
      scratchresfile << " " <<  rhsL(i,0);

// ______ Add coefficient number behind value ______
    scratchresfile << " \t" <<  coeffL << " " << coeffP << "\n";
    coeffL--;
    coeffP++;
    if (coeffL == -1)
      {
      coeffL = coeffP;
      coeffP = 0;
      }
    }

  coeffL = 0;
  coeffP = 0;
  scratchresfile << "\nEstimated_coefficientsP:\n";
  for (i=0; i<Nunk; i++)
    {
    if (rhsP(i,0) < 0.)
      scratchresfile <<         rhsP(i,0);
    else
      scratchresfile << " " <<  rhsP(i,0);

// ______ Add coefficient number behind value ______
    scratchresfile << " \t" <<  coeffL << " " << coeffP << "\n";
    coeffL--;
    coeffP++;
    if (coeffL == -1)
      {
      coeffL = coeffP;
      coeffP = 0;
      }
    }
  scratchresfile << "*******************************************************************"
                 //<< "\n* End_coregpm:_NORMAL"
                 << "\n* End_" << processcontrol[pr_i_coregpm] << "_NORMAL"
                 << "\n*******************************************************************\n";
  scratchresfile.close();


// ====== Compute offsets for corners ======
// BK 18-May-2000
  // read rhsL from file due top format... double
  matrix<real8> Lcoeff = readcoeff("scratchrescpm",
                     "Estimated_coefficientsL:",Ncoeffs(DEGREE));
  // read rhsP from file due top format... double
  matrix<real8> Pcoeff = readcoeff("scratchrescpm",
                     "Estimated_coefficientsP:",Ncoeffs(DEGREE));
  matrix<real4> x_axis(2,1);
  matrix<real4> y_axis(2,1);
  x_axis(0,0) = minL;
  x_axis(1,0) = maxL;
  y_axis(0,0) = minP;
  y_axis(1,0) = maxP;
  normalize(x_axis,minL,maxL);
  normalize(y_axis,minP,maxP);
  matrix<real4> offsetcornersL = polyval(x_axis,y_axis,Lcoeff);
  matrix<real4> offsetcornersP = polyval(x_axis,y_axis,Pcoeff);
  INFO.print(" ");
  INFO.print("Modeled transformation in azimuth:");
  INFO.print("-------------------------------------------------");
  INFO << "  First line:    " << offsetcornersL(0,0)  << " ... " << offsetcornersL(0,1);
  INFO.print();
  INFO.print("                    :           :");
  INFO << "  Last line:     " << offsetcornersL(1,0)  << " ... " << offsetcornersL(1,1);
  INFO.print();
  INFO.print("\n");
  INFO.print("Modeled transformation in range:");
  INFO.print("-------------------------------------------------");
  INFO << "  First line:    " << offsetcornersP(0,0)  << " ... " << offsetcornersP(0,1);
  INFO.print();
  INFO.print("                    :           :");
  INFO << "  Last line:     " << offsetcornersP(1,0)  << " ... " << offsetcornersP(1,1);
  INFO.print();
  INFO.print(" ");


  // ====== Dump evaluated polynomial if requested ======
  // BK 17-May-2000
  if (coregpminput.dumpmodel)
    {
    DEBUG.print("Do evaluation of coreg model with stepsize 100 pixels or so...");
    DEBUG.print("And account for currentwindow, not orig window...");
    PROGRESS.print("Started dumping evaluated model azimuth.");
    TRACE.print();// empty buffer to be sure
    TRACE << "offsetazi_" << originalmaster.lines() 
               <<          "_" << originalmaster.pixels() 
               << ".r4";
    char fileazi[ONE27];
    strcpy(fileazi,TRACE.get_str());
    TRACE.print();// empty buffer to be sure

    ofstream dumpfile;
    openfstream(dumpfile,fileazi,true);
    bk_assert(dumpfile,fileazi,__FILE__,__LINE__);

    // polyval both standing x,y... (?)
    // matrix<real4> p_axis(1,originalmaster.pixels());
    matrix<real4> l_axis(1,1);                  // ...
    matrix<real4> p_axis(originalmaster.pixels(),1);
    for (i=0; i<p_axis.pixels(); ++i)
      p_axis(i,0) = i+originalmaster.pixlo;                  // multilook==1 ?
    normalize(p_axis,minP,maxP);

    // azimuth
    for (i =originalmaster.linelo;
         i<=originalmaster.linehi; ++i) // all lines
      {
      //l_axis(0,0) = i;
      //normalize(l_axis,minL,maxL);
      l_axis(0,0) = normalize(real8(i),minL,maxL);
      //matrix<real4> MODEL = polyval(l_axis,p_axis,rhsL);
      matrix<real4> MODEL = polyval(l_axis,p_axis,Lcoeff);
      dumpfile << MODEL;
      }
    dumpfile.close();
    INFO << "Dumped model azimuth offset to file: " << fileazi
         << " format: real4; number of lines: "
         << originalmaster.lines()
         << " number of pixels: "
         << originalmaster.pixels();
    INFO.print();
    
    // ______ same for range ______
    PROGRESS.print("Started dumping evaluated model range.");
    TRACE.print();// empty buffer to be sure
    TRACE << "offsetrange_" << originalmaster.lines() 
             << "_" << originalmaster.pixels() << ".r4";
    char filerange[ONE27];
    strcpy(filerange,TRACE.get_str());
    TRACE.print();// empty buffer to be sure
    ofstream dumpfile2;
    openfstream(dumpfile2,filerange,true);
    bk_assert(dumpfile2,filerange,__FILE__,__LINE__);

    for (i =originalmaster.linelo;
         i<=originalmaster.linehi; ++i) // all lines
      {
      l_axis(0,0) = normalize(real8(i),minL,maxL);
      matrix<real4> MODEL = polyval(l_axis,p_axis,Pcoeff);
      dumpfile2 << MODEL;
      }
    INFO << "Dumped model range offset to file: " << filerange
         << " format: real4; number of lines: "
         << originalmaster.lines()
         << " number of pixels: "
         << originalmaster.pixels();
    INFO.print();
    dumpfile2.close();
    }

  // ====== Tidy up ======
  PROGRESS.print("finished computation of coregistration parameters.");
  } // END coregpm



/****************************************************************
 *    getofffile                                                *
 *                                                              *
 * Returns matrix (real4) with data of fine coreg from file     *
 *   mat(i,0)=window number                                     *
 *   mat(i,1)=position: line coordinate                         *
 *   mat(i,2)=position: pixels coordinate                       *
 *   mat(i,3)=offset: line direction                            *
 *   mat(i,4)=offset: pixle direction                           *
 *   mat(i,5)=correlation:                                      *
 * searches for "Number_of_correlation_windows:"                *
 *                                                              *
 *    Bert Kampes, 24-Feb-1999                                  *
 ****************************************************************/
matrix<real4> getofffile(
        const char* file,
        real4 threshold)
  {
  TRACE_FUNCTION("getofffile (BK 24-Feb-1999)");
  char                  dummyline[ONE27];
  char                  word[EIGHTY];
  bool                  foundsection = false;

  ifstream infile;
  openfstream(infile,file);
  bk_assert(infile,file,__FILE__,__LINE__);

  //ifstream infile(file, ios::in);
  //if (!infile)
  //  {
  //  DEBUG << "getofffile: problem opening file: " << file;
  //  DEBUG.print();
  //  }

// ======Search file for data section======
  while (infile)
    {
    infile >> word;
    if (strcmp("Number_of_correlation_windows:",word))  // no pattern match.
      {
      infile.getline(dummyline,ONE27,'\n');             // goto next line.
      }
    else                                                // in data section
      {
      foundsection=true;
      int32 N;                                          // number of points
      infile >> N;
      infile.getline(dummyline,ONE27,'\n');             // next line
      infile.getline(dummyline,ONE27,'\n');             // skip line with info
      int32 pos  = infile.tellg();                      // position of start data
      int32 Nobs = 0;                                   // number points > threshold
      real4 winnumber, posL, posP, offL, offP, corr;    // on file
      register int32 i;
      for (i=0;i<N;i++)
        {
        infile >> winnumber >> posL >> posP >> offL >> offP >> corr;
        infile.getline(dummyline,ONE27,'\n');           // goto next data record
        if (corr > threshold)
          Nobs++; 
        }

      if (Nobs == 0)
        {
        PRINT_ERROR("code ???: No data found > threshold.")
        throw(some_error);
        }

      matrix<real4> Data(Nobs,6);
      infile.seekg(pos);                                // return to start data
      int32 cnti = -1;
      for (i=0;i<N;i++)
        {
        infile >> winnumber >> posL >> posP >> offL >> offP >> corr;
        infile.getline(dummyline,ONE27,'\n');           // goto next data record
        if (corr > threshold)
          {
          cnti++;
          Data(cnti,0) = winnumber;
          Data(cnti,1) = posL;
          Data(cnti,2) = posP;
          Data(cnti,3) = offL;
          Data(cnti,4) = offP;
          Data(cnti,5) = corr;
          }
        }

      infile.close();
      return Data;
      } // else
    } // file

// ______Tidy up______
  if (!foundsection)
    {
    PRINT_ERROR("code 401: getofffile: couldn't find data section in file.");
    throw(some_error);
    }
  infile.close();

  // --- return a dummy here since some compiler like that ---
  return matrix<real4>(999,999);// BK 07-Apr-2003
  } // END getofffile



/****************************************************************
 *    cc4                                                       *
 *                                                              *
 * cubic convolution 4 points                                   *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> cc4(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("cc4 (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("cc4: standing vectors only.")
    throw(input_error);
    }

  real4 alpha = -1.;
  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    {
    real4 xx2 = sqr(x(i,0));
    real4 xx  = sqrt(xx2);
    if      (xx < 1)
      y(i,0) = (alpha+2)*xx2*xx - (alpha+3)*xx2 + 1;
    else if (xx < 2)
      y(i,0) = alpha*xx2*xx - 5*alpha*xx2 + 8*alpha*xx - 4*alpha;
    else 
      y(i,0) = 0.;
    }
  return y;
  } // END cc4


/****************************************************************
 *    cc6                                                       *
 *                                                              *
 * cubic convolution 6 points                                   *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 * corrected (alfa+beta)->(alfa-beta) after correction in paper *
 * by Ramon Hanssen                                             *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> cc6(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("cc6 (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("cc6: standing vectors only.")
    throw(input_error);
    }

  real4 alpha = -.5;
  real4 beta  =  .5;
  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    {
    real4 xx2 = sqr(x(i,0));
    real4 xx  = sqrt(xx2);
    if      (xx < 1)
      y(i,0) = (alpha-beta+2)*xx2*xx - (alpha-beta+3)*xx2 + 1;
//      y(i,0) = (alpha+beta+2)*xx2*xx - (alpha+beta+3)*xx2 + 1;
    else if (xx < 2)
      y(i,0) =   alpha*xx2*xx - (5*alpha-beta)*xx2 
               + (8*alpha-3*beta)*xx - (4*alpha-2*beta);
    else if (xx < 3)
      y(i,0) = beta*xx2*xx - 8*beta*xx2 + 21*beta*xx - 18*beta;
    else 
      y(i,0) = 0.;
    }
  return y;
  } // END cc6


/****************************************************************
 *    ts6                                                       *
 *                                                              *
 * truncated sinc 6 points                                      *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> ts6(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("ts6 (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("ts6: standing vectors only.")
    throw(input_error);
    }

  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = sinc(x(i,0)) * rect(x(i,0)/6.);
  return y;
  } // END ts6


/****************************************************************
 *    ts8                                                       *
 *                                                              *
 * truncated sinc 8 points                                      *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> ts8(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("ts8 (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("ts8: standing vectors only.")
    throw(input_error);
    }

  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = sinc(x(i,0)) * rect(x(i,0)/8.);
  return y;
  } // END ts8


/****************************************************************
 *    ts16                                                      *
 *                                                              *
 * truncated sinc 16 points                                     *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> ts16(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("ts16 (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("ts16: standing vectors only.")
    throw(input_error);
    }
  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = sinc(x(i,0)) * rect(x(i,0)/16.);
  return y;
  } // END ts16


/****************************************************************
 *    rect                                                      *
 *                                                              *
 * rect function for matrix (stepping function?)                *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> rect(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("rect (BK 16-Mar-1999)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("rect: standing vectors only.");
    throw(input_error);
    }

  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = rect(x(i,0));
  return y;
  } // END rect


/****************************************************************
 *    tri                                                       *
 *                                                              *
 * tri function for matrix (piecewize linear?, triangle)        *
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 ****************************************************************/
matrix<real4> tri(
        const matrix<real4> &x)
  {
  TRACE_FUNCTION("tri (BK 16-Mar-1999)")
  if (x.pixels() != 1)
    {
    PRINT_ERROR("tri: standing vectors only.")
    throw(input_error);
    }
  matrix<real4> y(x.lines(),1);
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = tri(x(i,0));
  return y;
  } // END tri


/****************************************************************
 *    knab                                                      *
 *                                                              *
 * KNAB window of N points, oversampling factor CHI             *
 *                                                              *
 * defined by: Migliaccio IEEE letters vol41,no5, pp1105,1110, 2003 *
 * k = sinc(x).*(cosh((pi*v*L/2)*sqrt(1-(2.*x./L).^2))/cosh(pi*v*L/2));
 *                                                              *
 * input:                                                       *
 *  - x-axis                                                    *
 *  - oversampling factor of bandlimited sigal CHI              *
 *  - N points of kernel size                                   *
 * output:                                                      *
 *  - y=f(x); function evaluated at x                           *
 *                                                              *
 *    Bert Kampes, 22-DEC-2003                                  *
 ****************************************************************/
matrix<real4> knab(
        const matrix<real4> &x,
        const real4 CHI,
        const int32 N)
  {
  TRACE_FUNCTION("knab (BK 22-Dec-2003)");
  if (x.pixels() != 1)
    {
    PRINT_ERROR("knab: standing vectors only.")
    throw(input_error);
    }
  matrix<real4> y(x.lines(),1);
  real4 v  = 1.0-1.0/CHI;
  real4 vv = PI*v*real4(N)/2.0;
  for (register int32 i=0;i<y.lines();i++)
    y(i,0) = sinc(x(i,0))*cosh(vv*sqrt(1.0-sqr(2.0*x(i,0)/real4(N))))/cosh(vv);
  return y;
  } // END knab


/****************************************************************
 *    resample                                                  *
 *                                                              *
 * Resample slave to master grid based on coregistration        *
 *  parameters.                                                 *
 * if dbow==0 then default to overlap, else dbow,               *
 *  write 0's where it does not overlap                         *
 * (later at interf.comp. if master<slave then doris exits! bug)*
 *                                                              *
 * input:                                                       *
 *  - inputoptions                                              *
 * output:                                                      *
 *  - void (file)                                               *
 *                                                              *
 *    Bert Kampes, 16-Mar-1999                                  *
 * added DBOW master add zeros.                                 *
 #%// BK 21-Aug-2000BOW master add zeros.                       *
 * shift data to center of spectrum before resampling, shift    *
 * back afterwards. (see e.g. thesis Geudtner)                  *
 #%// BK 09-Nov-2000                                            *
 * Seems to be a bug in shifting the data spectrum if more      *
 * buffers are used, working on it.                             *
 * (Increase FORSURE variable if crash)                         *
 #%// BK 19-Nov-2000                                            *
 ****************************************************************/
void resample(
        const input_gen         &generalinput,
        const input_resample    &resampleinput,
        const slcimage          &master,
        const slcimage          &slave,
        const matrix<real8>     &cpmL,          // coregistration parameters
        const matrix<real8>     &cpmP           // coregistration parameters
)
  {
  TRACE_FUNCTION("resample (BK 16-Mar-1999; BK 09-Nov-2000)")
  if (resampleinput.shiftazi==true)
    DEBUG.print("shifting kernelL to data fDC BK 26-Oct-2002");
  // ___ Handle input ___
  const uint BUFFERMEMSIZE = generalinput.memory;       // Bytes
  const int32 Npoints      = resampleinput.method%100;  // #pnts interpolator
  if (isodd(Npoints))
    {
    PRINT_ERROR("resample only even point interpolators.")
    throw(input_error);
    }
  const int32 Npointsd2    = Npoints/2;
  const int32 Npointsd2m1  = Npointsd2-1;
  const int32 Npointsm1    = Npoints-1;
  const uint  Sfilelines   = slave.currentwindow.lines();
  const uint sizeofci16    = sizeof(compli16);
  const uint sizeofcr4     = sizeof(complr4);

  // ______ Normalize data for polynomial ______
  const real8 minL         = master.originalwindow.linelo;
  const real8 maxL         = master.originalwindow.linehi;
  const real8 minP         = master.originalwindow.pixlo;
  const real8 maxP         = master.originalwindow.pixhi;
  INFO << "resample: polynomial normalized by factors: "
       << minL << " " << maxL << " " << minP << " " << maxP << " to [-2,2]";
  INFO.print();

  // ______ For KNAB window only if requested ______
  real4 CHI1 = slave.prf/slave.abw;// oversampling factor az
  real4 CHI2 = (slave.rsr2x/2.0)/slave.rbw;// oversampling factor rg
  real4 CHI  = (CHI1+CHI2)/2.0;// mean oversampling factor
  DEBUG << "Oversampling factor azimuth:   " << CHI1;
  DEBUG.print();
  DEBUG << "Oversampling factor azimuth:   " << CHI2;
  DEBUG.print();
  DEBUG << "Mean oversampling factor data: " << CHI;
  DEBUG.print();
  if (CHI < 1.1)
    {
    WARNING << "mean oversampling factor data " << CHI << " not optimal for KNAB";
    WARNING.print();
    }


  // ======Create lookup table======
  // ______e.g. four point interpolator
  // ______interpolating point: p=6.4925
  // ______required points: 5, 6, 7, 8
  // ______kernel number from lookup table: floor(.4925*INTERVAL+.5)
  // ______ table[0]= 0 1 0 0 ;table[INTERVAL]= 0 0 1 0
  // ______intervals in lookup table: dx
  register int32 i;
  // ______ it seems 20 is enough, but as fast as 40... ______
  //const int32 INTERVAL = 20;                          // precision 1./INTERVAL
  //const int32 INTERVAL = 40;                          // precision: 1./INTERVAL
  const int32 INTERVAL = 100;                           // precision: 1./INTERVAL
  const real8 dx = 1./INTERVAL;                         // interval look up table
  const int32 Ninterval = INTERVAL + 1;                 // size of lookup table
  INFO << "resample: lookup table size: " << Ninterval;
  INFO.print();

  matrix<real4> x_axis(Npoints,1);
  for (i=0; i<Npoints; ++i)
    x_axis(i,0) = 1. - Npointsd2 + i;                   // start at [-1 0 1 2]

  // ______ Lookup table complex because of multiplication with complex ______
  // ______ This is faster if veclib is used, slower if for loop... ______
  // ______ Could make a loopkup table for azimuth and range and ______
  // ______ shift spectrum of azi kernel with doppler centroid (not done) ______
  matrix<complr4> *pntM[Ninterval];
  // ______ same axis required for shift azimuth spectrum as used ______
  // ______ for kernel to avoid phase shift (Raffaele Nutricato) ______
  matrix<real4>   *pntAxis[Ninterval];
  for (i=0; i<Ninterval; ++i)
    {
    pntM[i]    = new matrix<complr4> (Npoints,1);
    pntAxis[i] = new matrix<real4>   (Npoints,1);// used only for azishift
    switch(resampleinput.method)
      {
      case rs_rect:
        (*pntM[i]) = mat2cr4(rect(x_axis)); 
        break;
      case rs_tri:
        (*pntM[i]) = mat2cr4(tri(x_axis));  
        break;
      case rs_cc4p:
        (*pntM[i]) = mat2cr4(cc4(x_axis));  
        break;
      case rs_cc6p:
        (*pntM[i]) = mat2cr4(cc6(x_axis));  
        break;
      case rs_ts6p:
        (*pntM[i]) = mat2cr4(ts6(x_axis));  
        break;
      case rs_ts8p:
        (*pntM[i]) = mat2cr4(ts8(x_axis));  
        break;
      case rs_ts16p:
        (*pntM[i]) = mat2cr4(ts16(x_axis)); 
        break;
      case rs_knab4p:
        (*pntM[i]) = mat2cr4(knab(x_axis,CHI,4)); 
      case rs_knab6p:
        (*pntM[i]) = mat2cr4(knab(x_axis,CHI,6)); 
      case rs_knab8p:
        (*pntM[i]) = mat2cr4(knab(x_axis,CHI,8)); 
      case rs_knab10p:
        (*pntM[i]) = mat2cr4(knab(x_axis,CHI,10)); 
      case rs_knab16p:
        (*pntM[i]) = mat2cr4(knab(x_axis,CHI,16)); 
        break;
      default:
        PRINT_ERROR("impossible.")
        throw(unhandled_case_error);
      }
    (*pntAxis[i]) = x_axis;// to shift kernelL use: k*=exp(-i*2pi*axis*fdc/prf) 
    x_axis -= dx;               // Note: 'wrong' way (mirrored)
    }
  // Bert: usage: 
  //  pntM[0]->showdata(); or (*pntM[0][0]).showdata(); 
  PROGRESS.print("Resample: lookup table created (kernel and axis).");

  // ______Save some time by computing degree here______
  int32 degree_cpmL = degree(cpmL.size());
  int32 degree_cpmP = degree(cpmP.size());


  // ______Corners of overlap in slave system______
  // ______offset = A(slave system) - A(master system)______
  // BK 21-09-00, get approx overlap here, what is diff with below???
  window approxoverlap = getoverlap(master,slave,cpmL,cpmP);

  // compute it here as well??? (doesn't hurt)
  real8 o00L = approxoverlap.linelo 
                + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
                          normalize(real8(approxoverlap.pixlo),minP,maxP),
                          cpmL,degree_cpmL);
  real8 o00P = approxoverlap.pixlo
                + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
                          normalize(real8(approxoverlap.pixlo),minP,maxP),
                          cpmP,degree_cpmP);
  real8 o0NL = approxoverlap.linelo
                + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
                          normalize(real8(approxoverlap.pixhi),minP,maxP),
                          cpmL,degree_cpmL);
  real8 o0NP = approxoverlap.pixhi
                + polyval(normalize(real8(approxoverlap.linelo),minL,maxL),
                          normalize(real8(approxoverlap.pixhi),minP,maxP),
                          cpmP,degree_cpmP);
  real8 oN0L = approxoverlap.linehi
                + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
                          normalize(real8(approxoverlap.pixlo),minP,maxP),
                          cpmL,degree_cpmL);
  real8 oN0P = approxoverlap.pixlo
                + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
                          normalize(real8(approxoverlap.pixlo),minP,maxP),
                          cpmP,degree_cpmP);
  real8 oNNL = approxoverlap.linehi
                + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
                          normalize(real8(approxoverlap.pixhi),minP,maxP),
                          cpmL,degree_cpmL);
  real8 oNNP = approxoverlap.pixhi
                + polyval(normalize(real8(approxoverlap.linehi),minL,maxL),
                          normalize(real8(approxoverlap.pixhi),minP,maxP),
                          cpmP,degree_cpmP);


  // ______Compute if it is possible to resample entire overlap______
  // ______Assume +- linear transformation model?______
  //  min required line coord: ceil(min(oL00,oL0N)) - Npointsd2
  //  max required line coord: floor(max(oLN0,oLNN)) + Npointsd2
  //  difference with present: correction = slave.linelo - min.required
  int32 correctlinelo = slave.currentwindow.linelo 
                        - (int32(ceil(min(o00L,o0NL))  - Npointsd2));
  if (correctlinelo < 0) correctlinelo = 0;

  int32 correctpixlo  = slave.currentwindow.pixlo
                        - (int32(ceil(min(o00P,oN0P))  - Npointsd2));
  if (correctpixlo  < 0) correctpixlo = 0;

  int32 correctlinehi = slave.currentwindow.linehi
                        - (int32(floor(max(oN0L,oNNL)) + Npointsd2));
  if (correctlinehi > 0) correctlinehi = 0;

  int32 correctpixhi  = slave.currentwindow.pixhi
                        - (int32(floor(max(o0NP,oNNP)) + Npointsd2));
  if (correctpixhi  > 0) correctpixhi = 0;

  //  window overlap;   // in master system -> res file for update productinfo
  // no problem with uint initialization
  window overlap(approxoverlap.linelo + correctlinelo,
                 approxoverlap.linehi + correctlinehi,
                 approxoverlap.pixlo  + correctpixlo,
                 approxoverlap.pixhi  + correctpixhi);

  // ====== Adjust overlap possibly for RS_DBOW card ======
  int32 write0lines1  = 0;                      // DBOW card, 0's at start
  int32 write0linesN  = 0;
  int32 write0pixels1 = 0;
  int32 write0pixelsN = 0;
  if (!(resampleinput.dbow.linelo == 0 &&       // as such initialized by readinput
        resampleinput.dbow.linehi == 0 &&
        resampleinput.dbow.pixlo  == 0 &&
        resampleinput.dbow.pixhi  == 0    ))
    {
    // ______ Check if overlap is large enough to contain DBOW ______
    if (resampleinput.dbow.linelo > overlap.linehi)
      {
      PRINT_ERROR("RS_DBOW: specified min. line larger than max. line of overlap.")
      throw(input_error);
      }
    if (resampleinput.dbow.linehi < overlap.linelo)
      {
      PRINT_ERROR("RS_DBOW: specified max. line smaller than min. line of overlap.")
      throw(input_error);
      }
    if (resampleinput.dbow.pixlo > overlap.pixhi)
      {
      PRINT_ERROR("RS_DBOW: specified min. pixel larger than max. pixel of overlap.")
      throw(input_error);
      }
    if (resampleinput.dbow.pixhi < overlap.pixlo)
      {
      PRINT_ERROR("RS_DBOW: specified max. pixel smaller than min. pixel of overlap.")
      throw(input_error);
      }

    write0lines1  =  overlap.linelo - resampleinput.dbow.linelo;
    if ( write0lines1 < 0 ) write0lines1 = 0;   // smaller window selected
    write0linesN  = -overlap.linehi + resampleinput.dbow.linehi;
    if ( write0linesN < 0 ) write0linesN = 0;   // smaller window selected
    write0pixels1  =  overlap.pixlo - resampleinput.dbow.pixlo;
    if ( write0pixels1 < 0 ) write0pixels1 = 0; // smaller window selected
    write0pixelsN  = -overlap.pixhi + resampleinput.dbow.pixhi;
    if ( write0pixelsN < 0 ) write0pixelsN = 0; // smaller window selected
    if (resampleinput.dbow.linelo < overlap.linelo)
      {
      WARNING << "RS_DBOW: min. line < overlap (writing: " << write0lines1
           << " lines with zeros before first resampled line).";
      WARNING.print();
      }
    else
      overlap.linelo = resampleinput.dbow.linelo;       // correct it
    if (resampleinput.dbow.linehi > overlap.linehi)
      {
      WARNING << "RS_DBOW: max. line > overlap (writing: " << write0linesN
           << " lines with zeros after last resampled line).";
      WARNING.print();
      }
    else
      overlap.linehi = resampleinput.dbow.linehi;       // correct it
    if (resampleinput.dbow.pixlo < overlap.pixlo)
      {
      WARNING << "RS_DBOW: min. pixel < overlap (writing: " << write0pixels1
           << " columns with zeros before first resampled column).";
      WARNING.print();
      }
    else
      overlap.pixlo = resampleinput.dbow.pixlo;         // correct it
    if (resampleinput.dbow.pixhi > overlap.pixhi)
      {
      WARNING << "RS_DBOW: max. pixel > overlap (writing: " << write0pixelsN
           << " columns with zeros after last resampled column).";
      WARNING.print();
      }
    else
      overlap.pixhi = resampleinput.dbow.pixhi;         // correct it
    } // adjust overlap


  // ______ Buffersize output matrix ______
  const int32 Npointsxsize = Npoints*sizeofcr4; // size for memcpy (fill PART)
  const int32 npixels      = slave.currentwindow.pixels();
  const real8 bytesperline = sizeofcr4 * npixels;
  // ___ COMMENTED OUT, OLD WAY SHIFT DATA, now shiftkernel ___
  //const real8 bigmatrices  = (resampleinput.shiftazi==true) ? 6:2.5;
  const real8 bigmatrices  = 2.5;
  const int32 nlines       = int32((BUFFERMEMSIZE/bigmatrices)/bytesperline);

  // ______ Declare/allocate matrices ______
  matrix<complr4> BUFFER;                       // load after output is written
  matrix<complr4> RESULT(nlines,overlap.pixhi-overlap.pixlo+1);
  matrix<complr4> PART(Npoints,Npoints);

#ifdef __USE_VECLIB_LIBRARY__
  matrix<complr4> TMPRES(Npoints,1);
  int32 Np = Npoints;                           // must be non-constant
  int32 ONEint = 1;                             // must have pointer to 1
  complr4 c4alpha(1.,0.);
  complr4 c4beta(0.,0.);
  STUPID_cr4 ANS;                               // VECLIB struct return type
#endif


  // ====== Open output file ======
  ofstream ofile;
  openfstream(ofile,resampleinput.fileout,generalinput.overwrit);
  bk_assert(ofile,resampleinput.fileout,__FILE__,__LINE__);

  // ________ First write zero lines if appropriate (DBOW) ______
  register int32 thisline, thispixel;           // local loop variable
  switch (resampleinput.oformatflag)
    {
    case FORMATCR4:
      {
      const complr4 zerocr4(0,0);
      for (thisline=0; thisline<write0lines1; ++thisline)
        for (thispixel=0;
             thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
             ++thispixel)
          ofile.write((char*)&zerocr4,sizeofcr4);
      break;
      }
    case FORMATCI2:
      {
      const compli16 zeroci16(0,0);
      for (thisline=0; thisline<write0lines1; ++thisline)
        for (thispixel=0; 
             thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
             ++thispixel)
          ofile.write((char*)&zeroci16,sizeofci16);
      break;
      }
    default:
      PRINT_ERROR("impossible format")
      throw(unhandled_case_error);
    }

  // ______ Info ______
  INFO << "Overlap window: "
       << overlap.linelo << ":" << overlap.linehi << ", "
       << overlap.pixlo  << ":" << overlap.pixhi;
  INFO.print();

  // ______ Progress messages ______
  int32 percent    = 0;
  int32 tenpercent = int32((overlap.lines()/10)+.5);  // round
  if (tenpercent==0) tenpercent = 1000;                    // avoid error: x%0

  // ====== Resample all lines that are requested ======
  bool newbufferrequired = true;                // read initial slave buffer
  register int32 linecnt = -1;                  // indicate output buffer full
  int32 firstline;                              // slave system
  register int32 line;                          // loop counter master system
  register int32 pixel;                         // loop counter master system
  for (line=overlap.linelo; line<=overlap.linehi; line++)
    {
    // ______ Progress messages ______
    if (((line-overlap.linelo)%tenpercent)==0)
      {
      PROGRESS << "RESAMPLE: " << setw(3) << percent << "%";
      PROGRESS.print();
      percent += 10;
      }

    // ====== Write RESULT to disk if it is full (write last bit at end) ======
    if (linecnt==RESULT.lines()-1)              // ==nlines
      {
      newbufferrequired = true;                 // do load slave from file
      DEBUG << "Writing slave: ["
           << line-RESULT.lines() << ":" << line-1 << ", "
           << overlap.pixlo << ":" << overlap.pixhi
           << "] (master coord. system)";
      DEBUG.print();
      linecnt = 0;
      // ______ Shift azimuth spectrum of resampled data back to fDC ______
      // ______ Require pixelnumber in slave system == pix_master + offsetP ______
      // ______ +shift for forward shifting ______
      // ___ COMMENTED OUT, OLD WAY SHIFT DATA ___
//      if (resampleinput.shiftazi==true)
//      {
//        real4 shift = overlap.pixlo +
//        polyval(normalize(real4(line-(linecnt/2.)),minL,maxL),// middle line
//                normalize(real4(overlap.pixlo),minP,maxP),    // first pixel
//                cpmP,degree_cpmP);                            // e.g. 2.5232
//        if (shift < 0)
//          ERROR("sorry, but this really is impossible.");
//        shiftazispectrum(RESULT,slave,shift);
//      }
      // ______ Actually write ______
      switch (resampleinput.oformatflag)
        {
        case FORMATCR4:
          {
          // old, now first write zeropixels...: ofile << RESULT;
          const complr4 zerocr4(0,0);
          for (thisline=0; thisline<RESULT.lines(); ++thisline)
            {
            // ______ Write zero pixels at start ______
            for (thispixel=0; thispixel<write0pixels1; ++thispixel)
              {
              ofile.write((char*)&zerocr4,sizeofcr4);
              }
            // ______ WRITE the interpolated data per row ______
            ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeof(RESULT(0,0)));
            // ______ Write zero pixels at end ______
            for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
              {
              ofile.write((char*)&zerocr4,sizeofcr4);
              }
            }
          break;
          }
        case FORMATCI2:
          {
          const compli16 zeroci16(0,0);
          compli16 castedresult;
          for (thisline=0; thisline<RESULT.lines(); ++thisline)
            {
            // ______ Write zero pixels at start ______
            for (thispixel=0; thispixel<write0pixels1; ++thispixel)
              {
              ofile.write((char*)&zeroci16,sizeofci16);
              }
            // ______ Write the interpolated data per row ______
            for (thispixel=0; thispixel<RESULT.pixels(); ++thispixel)
              {
              // no default conversion, this seems slow, test this (BK)
              castedresult = cr4toci2(RESULT(thisline,thispixel));
              ofile.write((char*)&castedresult,sizeofci16);
              }
            // ______ Write zero pixels at end ______
            for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
              {
              ofile.write((char*)&zeroci16,sizeofci16);
              }
            }
          break;
          }
        default:
          PRINT_ERROR("impossible format")
          throw(unhandled_case_error);
        }
      }

    else // output buffer not full yet
      {
      linecnt++;
      }

    // ====== Read slave buffer if justwritten || firstblock ======
    if (newbufferrequired==true)
      {
      newbufferrequired = false;        // only load after output written
      firstline = int32(ceil(min(line + 
                        polyval(normalize(real4(line),minL,maxL),
                                normalize(real4(overlap.pixlo),minP,maxP),
                                cpmL,degree_cpmL),
                                 line + 
                        polyval(normalize(real4(line),minL,maxL),
                                normalize(real4(overlap.pixhi),minP,maxP),
                                cpmL,degree_cpmL))))
                             - Npoints;
      int32 line2 = line + nlines - 1;
      int32 lastline  = int32(ceil(min(line2 + 
                        polyval(normalize(real4(line2),minL,maxL),
                                normalize(real4(overlap.pixlo),minP,maxP),
                                cpmL,degree_cpmL),
                                 line2 + 
                        polyval(normalize(real4(line2),minL,maxL),
                                normalize(real4(overlap.pixhi),minP,maxP),
                                cpmL,degree_cpmL))))
                             + Npoints;
      //const int32 FORSURE = 5;                // buffer larger 2*FORSURE start/end
      const int32 FORSURE = 25;         // buffer larger 2*FORSURE start/end
      firstline -= FORSURE;
      lastline  += FORSURE;
      // ______ Don't compare apples with pears, uint<->int! ______
      if (firstline < int32(slave.currentwindow.linelo))
        firstline = slave.currentwindow.linelo;
      if (lastline > int32(slave.currentwindow.linehi)) 
        lastline = slave.currentwindow.linehi;
      // ______ Fill slave BUFFER from disk ______
      window winslavefile(firstline, lastline,  // part of slave loaded
                           slave.currentwindow.pixlo,   // from file in BUFFER.
                           slave.currentwindow.pixhi);
      DEBUG << "Reading slave: ["
           << winslavefile.linelo << ":" << winslavefile.linehi << ", "
           << winslavefile.pixlo  << ":" << winslavefile.pixhi  << "]";
      DEBUG.print();
      BUFFER = slave.readdata(winslavefile);
      // ______ Shift azimuthspectrum from fDC to center ______
      // ______ - for inverse shift, pixlo for evaluating fDC ______
      // ___ COMMENTED OUT, OLD WAY SHIFT DATA ___
//      if (resampleinput.shiftazi==true)
//        shiftazispectrum(BUFFER,slave,-real4(slave.currentwindow.pixlo));
      } // ___end: Read new slave buffer to resample outputbuffer


    // ====== Actual resample all pixels this output line ======
    for (pixel=overlap.pixlo; pixel<=overlap.pixhi; pixel++)
      {
      // ______ Compute coregistration ______
      // bk 25-10-99 why don't i do this per buffer, that's faster. (but more mem)
      //interpL = line  + polyval(line,pixel,cpmL,degree_cpmL); // e.g. 255.35432
      //interpP = pixel + polyval(line,pixel,cpmP,degree_cpmP); // e.g. 2.5232
      // ______ BK USE normalized coordinates, do this smarter .... !!!!
      const real4 interpL = line  + 
        polyval(normalize(real4(line),minL,maxL),
                normalize(real4(pixel),minP,maxP),
                cpmL,degree_cpmL);                              // e.g. 255.35432
      const real4 interpP = pixel +
        polyval(normalize(real4(line),minL,maxL),
                normalize(real4(pixel),minP,maxP),
                cpmP,degree_cpmP);                              // e.g. 2.5232

      // ______ Get correct lines for interpolation ______
      // ______ Note: following 4 not ok, floor(x) != ceil(x)-1 ______
      // ______ correct to int32 for speed, floor is slower, all positive indices.
      //int32 firstL     = ceil(interpL) - Npointsd2;           // e.g. 254 (5 6 7)
      //int32 firstP     = ceil(interpP) - Npointsd2;           // e.g. 1   (2 3 4)
      //int32 kernelnoL  = floor(.5+interpLdec*INTERVAL) - 1;   // lookup table index
      //int32 kernelnoP  = floor(.5+interpPdec*INTERVAL) - 1;   // lookup table index
      const int32 fl_interpL = int32(interpL);
      const int32 fl_interpP = int32(interpP);
      const int32 firstL     = fl_interpL - Npointsd2m1;        // e.g. 254 (5 6 7)
      const int32 firstP     = fl_interpP - Npointsd2m1;        // e.g. 1   (2 3 4)
      const real4 interpLdec = interpL - fl_interpL;            // e.g. .35432
      const real4 interpPdec = interpP - fl_interpP;            // e.g. .5232

      // ___ copy kernels here, change kernelL if required ___
      // BK 26-Oct-2002
      int32 kernelnoL  = int32(interpLdec*INTERVAL + .5);       // lookup table index
      int32 kernelnoP  = int32(interpPdec*INTERVAL + .5);       // lookup table index
      matrix<complr4> kernelL = (*pntM[kernelnoL]);// local copy
      matrix<complr4> kernelP = (*pntM[kernelnoP]);// local copy


#ifdef __DEBUG
      // ______This shouldn't be possible...______
      if (firstL < slave.currentwindow.linelo)
        {
        WARNING.print("firstL smaller then on disk (required for interpolation). continuing");
        RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
        continue;               // with next pixel
        }
      if (firstL+Npointsm1 > slave.currentwindow.linehi)
        {
        WARNING.print("lastL larger then on disk (required for interpolation). continuing");
        RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
        continue;               // with next pixel
        }
      if (firstP < slave.currentwindow.pixlo)
        {
        WARNING.print("firstP smaller then on disk (required for interpolation). continuing");
        RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
        continue;               // with next pixel
        }
      if (firstP+Npointsm1 > slave.currentwindow.pixhi)
        {
        WARNING.print("lastP larger then on disk (required for interpolation). continuing");
        RESULT(linecnt,pixel-overlap.pixlo) = complr4(0.,0.);
        continue;               // with next pixel
        }
#endif


      // ______ Shift azimuth kernel with fDC before interpolation ______
      if (resampleinput.shiftazi==true)
        {
        // ___ Doppler centroid is function of range only ____
        //const real8 tau = interpP/(slave.rsr2x/2.0);
        //const real8 fDC = slave.f_DC_a0 + 
        //                  slave.f_DC_a1*tau +
        //                  slave.f_DC_a2*sqr(tau);
        //const real4 tmp = 2.0*PI*fDC/slave.prf;
        //const real4 tmp = 2.0*PI*pix2fdc(interpP,slave)/slave.prf;
        const real4 tmp = 2.0*PI*slave.pix2fdc(interpP)/slave.prf;
        // ___ to shift spectrum of convolution kernel to fDC of data, multiply
        // ___ in the space domain with a phase trend of -2pi*t*fdc/prf
        // ___ (to shift back (no need) use +fdc)
        for (i=0; i<Npoints; ++i)
          {
          //complr4 t = complr4(cos(i*2.0*PI*fDC/slave.prf),
          //                   -sin(i*2.0*PI*fDC/slave.prf)); // exp(i*trend)
          //kernelL(i,0) *= t;
          // ___ Modify kernel, shift spectrum to fDC ___
          real4 t       = ((*pntAxis[kernelnoL])(i,0))*tmp;
          kernelL(i,0) *= complr4(cos(t),-sin(t));// note '-' (see manual)
          }
        }


        // ______ For speed: define setdata internally (memcpy) ______
        //linelo = firstL - firstline;
        //pixlo  = firstP - 1;          // !! -1 should be -currentwin.pixlo
        for (i=0; i<Npoints; i++)
          memcpy(PART[i],
          BUFFER[i+firstL-firstline]+
          firstP-slave.currentwindow.pixlo,Npointsxsize);

        // ====== Some speed considerations ======
#ifdef __USE_VECLIB_LIBRARY__
        // ______Compute PART * kernelP______
        //cgemv("T",&Np,&Np,&c4alpha, PART[0],&Np,
        //      (*pntM[int32(interpPdec*INTERVAL+.5)])[0],&ONEint,
        //      &c4beta,TMPRES[0],&ONEint,1);
        cgemv("T",&Np,&Np,&c4alpha, PART[0],&Np,
              kernelP[0],&ONEint,
              &c4beta,TMPRES[0],&ONEint,1);

        // ______Compute Result * kernelL; put in matrix RESULT______
        //ANS = cdotu(&Np,TMPRES[0],&ONEint,
        //     (*pntM[int32(interpLdec*INTERVAL+.5)])[0],&ONEint);
        ANS = cdotu(&Np,TMPRES[0],&ONEint, kernelL[0],&ONEint);
        RESULT(linecnt,pixel-overlap.pixlo) = complr4(ANS.re,ANS.im);
#else // do use VECLIB
        // ______ NOvec slower, but ok: better use CR4*R4 ______
        //RESULT(linecnt,pixel-overlap.pixlo) = 
        //((matTxmat(PART * (*pntM[kernelnoP]), (*pntM[kernelnoL])))(0,0));
        RESULT(linecnt,pixel-overlap.pixlo) = 
             ((matTxmat(PART*kernelP, kernelL))(0,0));
#endif // VECLIB 
      } // for all pixels in overlap
    } // for all lines in overlap


  // ______ Write last lines of Result to disk (filled upto linecnt) ______
  DEBUG << "Writing slave: ["
       //<< overlap.linehi-linecnt+1 << ":" << overlap.linehi << ", "
       << overlap.linehi-linecnt << ":" << overlap.linehi << ", "
       << overlap.pixlo << ":" << overlap.pixhi
       << "] (master coord. system)";
  DEBUG.print();
  // ______ Shift azimuth spectrum of resampled data back to fDC ______
  // ______ coordinate of first column in slave system equals: ______
  // ______ offset defined as: cols=colm+offsetP ______
  // ______ This is not totally correct... (e.g., if slave is scaled back) ______
  // ______ +shift for forward shifting ______
  // ______ We shift to much here, but this should not mather. ______
//  if (resampleinput.shiftazi==true)
//    {
//    real4 shift = overlap.pixlo + polyval(normalize(real4(
//              overlap.linehi-(RESULT.lines()/2.)),minL,maxL), // middle line
//                  normalize(real4(overlap.pixlo),minP,maxP),  // first pixel
//                  cpmP,degree_cpmP);                          // e.g. 2.5232
//    if (shift < 0)
//      ERROR("sorry, this is impossible, use noshiftazi.");
//    shiftazispectrum(RESULT,slave,shift);
//    }

  // ______ Actually write ______
  switch (resampleinput.oformatflag)
    {
    case FORMATCR4:
      {
      const complr4 zerocr4(0,0);
      for (thisline=0; thisline<=linecnt; thisline++)
        {
        // ______ Write zero pixels at start ______
        for (thispixel=0; thispixel<write0pixels1; ++thispixel)
          {
          ofile.write((char*)&zerocr4,sizeofcr4);
          }
        // ______ WRITE the interpolated data per row ______
        ofile.write((char*)&RESULT[thisline][0],RESULT.pixels()*sizeofcr4);
        // ______ Write zero pixels at end ______
        for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
          {
          ofile.write((char*)&zerocr4,sizeofcr4);
          }
        }
      break;
      }
    case FORMATCI2:
      {
      const compli16 zeroci16(0,0);
      compli16 castedresult;
      for (thisline=0; thisline<=linecnt; thisline++)
        {
        // ______ Write zero pixels at start ______
        for (thispixel=0; thispixel<write0pixels1; ++thispixel)
          {
          ofile.write((char*)&zeroci16,sizeofci16);
          }
        // ______ Write the interpolated data per row ______
        for (thispixel=0; thispixel<RESULT.pixels(); thispixel++)
          {
          castedresult = cr4toci2(RESULT(thisline,thispixel));
          //?? BK 6/1 castedresult = cr4toci2(RESULT[thisline][thispixel].real());
          // castedresult = compli16(RESULT[thisline][thispixel].real(),
          //                       RESULT[thisline][thispixel].imag());
          ofile.write((char*)&castedresult,sizeofci16);
          }
        // ______ Write zero pixels at end ______
        for (thispixel=0; thispixel<write0pixelsN; ++thispixel)
          {
          ofile.write((char*)&zeroci16,sizeofci16);
          }
        }
      break;
      }
    default:
      PRINT_ERROR("impossible format")
      throw(unhandled_case_error);
    }


  // ====== Write last zero lines if appropriate (DBOW card) ======
  switch (resampleinput.oformatflag)
    {
    case FORMATCR4:
      {
      complr4 zerocr4(0,0);
      for (thisline=0; thisline<write0linesN; ++thisline)
        for (thispixel=0;
             thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
             ++thispixel)
          ofile.write((char*)&zerocr4,sizeofcr4);
      break;
      }
    case FORMATCI2:
      {
      compli16 zeroci16(0,0);
      for (thisline=0; thisline<write0linesN; ++thisline)
        for (thispixel=0;
             thispixel<RESULT.pixels()+write0pixels1+write0pixelsN;
             ++thispixel)
          ofile.write((char*)&zeroci16,sizeofci16);
      break;
      }
    default:
      PRINT_ERROR("impossible format")
      throw(unhandled_case_error);
    }
  ofile.close();



  // ====== Write results to slave resfile ======
  char rsmethod[EIGHTY];
  switch(resampleinput.method)
    {
    case rs_rect:
      strcpy(rsmethod,"nearest neighbour");
      break;
    case rs_tri:
      strcpy(rsmethod,"piecewise linear");
      break;
    case rs_cc4p:
      strcpy(rsmethod,"4 point cubic convolution");
      break;
    case rs_cc6p:
      strcpy(rsmethod,"6 point cubic convolution"); 
      break;
    case rs_ts6p:
      strcpy(rsmethod,"6 point truncated sinc"); 
      break;
    case rs_ts8p:
      strcpy(rsmethod,"8 point truncated sinc"); 
      break;
    case rs_ts16p:
      strcpy(rsmethod,"16 point truncated sinc"); 
      break;
    case rs_knab4p:
      strcpy(rsmethod,"4 point knab window"); 
      break;
    case rs_knab6p:
      strcpy(rsmethod,"6 point knab window"); 
      break;
    case rs_knab8p:
      strcpy(rsmethod,"8 point knab window"); 
      break;
    case rs_knab10p:
      strcpy(rsmethod,"10 point knab window"); 
      break;
    case rs_knab16p:
      strcpy(rsmethod,"16 point knab window"); 
      break;
    default:
      PRINT_ERROR("impossible.")
      throw(unhandled_case_error);
    }

  char rsoformat[EIGHTY];
  switch(resampleinput.oformatflag)
    {
    case FORMATCR4:
      strcpy(rsoformat,"complex_real4");
      break;
    case FORMATCI2:
      strcpy(rsoformat,"complex_short");
      break;
    default:
      PRINT_ERROR("impossible.")
      throw(unhandled_case_error);
    }


  ofstream scratchlogfile("scratchlogresample", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"resample: scratchlogresample",__FILE__,__LINE__);

  scratchlogfile
    << "\n\n*******************************************************************"
    << "\n*_Start_resample"
    << "\nData_output_file: \t\t\t"
    <<  resampleinput.fileout
    << "\nData_output_format: \t\t\t"
    << rsoformat
    << "\nInterpolation kernel: \t\t\t"
    <<  rsmethod
    << "\nResampled slave size in master system: \t"
    <<  overlap.linelo - write0lines1 << ", "
    <<  overlap.linehi + write0linesN << ", "
    <<  overlap.pixlo  - write0pixels1 << ", "
    <<  overlap.pixhi  + write0pixelsN
    //<< "\nData first shifted to center azimuth spectrum by: "
    //<< fDC...
    << "\n*******************************************************************\n";
  scratchlogfile.close();

  ofstream scratchresfile("scratchresresample", ios::out | ios::trunc);
  bk_assert(scratchresfile,"resample: scratchresresample",__FILE__,__LINE__);
  scratchresfile 
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_s_resample]
    << "\n*******************************************************************"
    << "\nShifted azimuth spectrum:             \t\t"
    <<  resampleinput.shiftazi
    << "\nData_output_file:                     \t\t"
    <<  resampleinput.fileout
    << "\nData_output_format:                   \t\t"
    << rsoformat
    << "\nInterpolation kernel:                 \t\t"
    <<  rsmethod
    << "\nFirst_line (w.r.t. original_master):  \t\t"
    <<  overlap.linelo - write0lines1
    << "\nLast_line (w.r.t. original_master):   \t\t"
    <<  overlap.linehi + write0linesN
    << "\nFirst_pixel (w.r.t. original_master): \t\t"
    <<  overlap.pixlo  - write0pixels1
    << "\nLast_pixel (w.r.t. original_master):  \t\t"
    <<  overlap.pixhi  + write0pixelsN
    << "\n*******************************************************************"
    << "\n* End_" << processcontrol[pr_s_resample] << "_NORMAL"
    << "\n*******************************************************************\n";
  scratchresfile.close();


// ______Tidy up______
//  for (i=0;i<Ninterval;i++)   // like this ???
//    delete    pntM[i];
//    delete [] pntM[i];
  DEBUG.print("Exiting resample.");

  } // END resample

---