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

geocode.cc

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /****************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/geocode.cc,v $    *
00032  * $Revision: 3.13 $                                            *
00033  * $Date: 2005/04/11 13:47:45 $                                 *
00034  * $Author: kampes $                                            *
00035  *                                                              *
00036  * implementation of computation of endproducts (DEM, defo.map) *
00037  * -slant range 2 height (schwabisch)                           *
00038  * -slant range 2 height (rodriguez, exact, some problems)      *
00039  * -slant range 2 height (ambiguity)                            *
00040  * -geocode heightmatrix to phi,lambda                          *
00041  ****************************************************************/
00042 
00043 
00044 #include "constants.hh"                 // global constants
00045 #include "matrixbk.hh"                  // my matrix class
00046 #include "slcimage.hh"                  // my slc image class
00047 #include "productinfo.hh"               // my 'products' class
00048 #include "geocode.hh"                   // header file
00049 #include "utilities.hh"                 // utils
00050 #include "ioroutines.hh"                // error messages
00051 #include "coregistration.hh"            // distributepoints
00052 #include "conversion.hh"                // 
00053 #include "exceptions.hh"                 // my exceptions class
00054 
00055 #include <strstream>                    // for memory stream
00056 #include <cctype>                       // isspace
00057 #include <algorithm>                    // max
00058 
00059 
00060 
00061 /****************************************************************
00062  *    slant2h eight (schwabisch)                                *
00063  *                                                              *
00064  * compute height in radar coded system (master):               *
00065  *                                                              *
00066  * Input:                                                       *
00067  *  -                                                           *
00068  * Output:                                                      *
00069  *  -                                                           *
00070  *                                                              *
00071  * See thesis swabisch for method                               *
00072  * 1.  compute reference phase for h=0,2000,4000 in Npoints     *
00073  * 2.  solve system: h(phi) = a_0 + a_1*phi + a_2*phi*phi       *
00074  *      (2nd degree 1D polynomial) for all Npoints              *
00075  * 3.  compute a_i (l,p) = DEGREE2D 2D polynomial               *
00076  * 4.0 set offset to the one of first pixel , add this to all   *
00077  *     this step is skipped, phase is w.r.t. h=0, ref. is subtracted *
00078  * 4.1 evaluate polynomial of 3. for all points (l,p) of        *
00079  *      (multilooked) unwrapped interferogram                   *
00080  * solution to system for betas seems not very stable?          *
00081  *                                                              *
00082  *    Bert Kampes, 02-Jun-1999                                  *
00083  ****************************************************************/
00084 void slant2hschwabisch(
00085         const input_gen     &generalinput,
00086         const input_slant2h &slant2hinput,
00087         const input_ell     &ellips,
00088         const slcimage      &master,
00089         const slcimage      &slave,
00090         const productinfo   &unwrappedinterf,
00091         orbit               &masterorbit,
00092         orbit               &slaveorbit)
00093   {
00094   TRACE_FUNCTION("slant2hschwabisch (BK 02-Jun-1999)")
00095 
00096   const int32 MAXITER   = 10;
00097   const real8 CRITERPOS = 1e-6;
00098   const real8 CRITERTIM = 1e-10;
00099   const real8 m_minpi4cdivlambda = (-4.*PI*SOL)/master.wavelength;
00100   const real8 s_minpi4cdivlambda = (-4.*PI*SOL)/slave.wavelength;
00101 
00102   const int32 Npoints   = slant2hinput.Npoints;         // where ref.phase is evaluated
00103   const int32 DEGREE1D  = slant2hinput.degree1d;        // only possible now.
00104   const int32 DEGREE2D  = slant2hinput.degree2d;
00105   const int32 Nheights  = slant2hinput.Nheights;
00106 
00107   const int32 MAXHEIGHT = 5000;                         // max hei for ref.phase
00108   const int32 TEN       = 10;                           // used in pointer
00109   if (DEGREE1D - 1 > TEN)
00110     {
00111     PRINT_ERROR("panic, programmers problem: increase TEN.")
00112     throw(some_error);
00113     }
00114 
00115   const int32 HEIGHTSTEP = MAXHEIGHT / (Nheights - 1);  // heights to eval ref.phase
00116 
00117 
00118 
00119 // ______ Matrices for storing phase for all ref. ellipsoids ______
00120 // ______ PHASE(i,0)  phase for height 0
00121 // ______ PHASE(i,1)  phase for height Heigthsep * 1
00122 // ______ PHASE(i,Nh) phase for height 4000
00123   matrix<real8>         PHASE(Npoints,Nheights);        // pseudo-observation
00124 
00125 
00126 // ______ Distribute points in original master system (not multilooked) ______
00127 // ______ (i,0): line, (i,1): pixel, (i,2) flagfromdisk (not used here) ______
00128   matrix<uint> Position = distributepoints(Npoints,unwrappedinterf.win);
00129 
00130   cn pospoint;                          // point, returned by lp2xyz
00131   input_ell ELLIPS;
00132 
00133 
00134 
00135 // ====== STEP 1 ======
00136   PROGRESS.print("S2H: schwabisch: STEP1: compute reference phase for Nheights.");
00137 // ====== Compute reference phase in N points for height (numheight) ======
00138   register int32 numheight;
00139   register int32 i,j,k,l,index;
00140   register int32 alfa;
00141   for (numheight=0; numheight<Nheights; numheight++)
00142     {
00143     const int32 HEIGHT = numheight * HEIGHTSTEP;
00144     ELLIPS.a   = ellips.a + HEIGHT;
00145     ELLIPS.b   = ellips.b + HEIGHT;
00146     ELLIPS.ecc1st_sqr();
00147     ELLIPS.ecc2nd_sqr();
00148     
00149     // ______ Compute delta r for all points ______
00150     for (i=0;i<Npoints;i++)
00151       {
00152       const real8 line  = Position(i,0);
00153       const real8 pixel = Position(i,1);
00154       const real8 m_trange = master.pix2tr(pixel);
00155   
00156       // ______ Compute xyz of point P on ELLIPS for this line,pixel ______
00157       lp2xyz(line,pixel,ELLIPS,master,masterorbit,
00158              pospoint,MAXITER,CRITERPOS);
00159   
00160       // ______ Compute xyz of slave satelite in orbit_slave from P ______
00161       real8 s_tazi;                             // returned
00162       real8 s_trange;                           // returned, unused?
00163       xyz2t(s_tazi,s_trange,slave, slaveorbit,
00164             pospoint,MAXITER,CRITERTIM);
00165   
00166       // ______ Compute delta range ~= phase, store in matrix ______
00167       // ______ real8 dr = dist(m_possat,pospoint) - dist(s_possat,pospoint);
00168       // ______ real8 phase = -pi4*(dr/LAMBDA);
00169       // ______  dr    == M-S   cause if no flatearth M-S - flatearth = M-S-(M-S)=0
00170       // ______  phase == -4pi*dr/lambda == 4pi*(S-M)/lambda
00171       //  PHASE(i,numheight) = pi4divlambda*
00172       //    (s_possat.dist(pospoint)-m_possat.dist(pospoint));
00173       //PHASE(i,numheight) = minpi4cdivlambda*(m_trange-s_trange);
00174       PHASE(i,numheight) = m_minpi4cdivlambda*m_trange-
00175                            s_minpi4cdivlambda*s_trange;
00176       }
00177     }
00178 
00179 
00180   // ______ Subtract ref. phase at h=0 for all point ______
00181   // ______ this is the same as adding reference phase for all in uint ______
00182   // ______ BK tested 4/oct/99 ______
00183   // ______ This possibly needs to be CHANGED later for other situations ______
00184   for (i=0; i<Npoints; ++i)
00185     {
00186     real8 offset = PHASE(i,0);
00187     PHASE(i,0) = 0.;
00188     for (j=1; j<numheight; ++j)
00189       {
00190       PHASE(i,j) -= offset;
00191       }
00192     }
00193 
00194 
00195 
00196 // ====== STEP 2 ======
00197   PROGRESS.print("S2H: schwabisch: STEP2: estimate coefficients 1d polynomial.");
00198   // ====== Compute alpha coefficients of polynomials for these points ======
00199   // ______ h = sum (alpha_i phi^i); ______
00200   // ______ bk 26/10/99 rescale phi -> phi[0,1] ______
00201   matrix<real8>         DESIGN(Nheights,DEGREE1D+1);    // design matrix
00202   matrix<real8>         ALPHAS(Npoints,DEGREE1D+1);     // pseudo-observation
00203   matrix<real8>         HEI(Nheights,1);
00204   for (i=0; i<Nheights; i++)
00205     HEI(i,0) = i*HEIGHTSTEP;                            // 0, .., 5000
00206 
00207   // ______ normalize data to [0,1] ______
00208   const real8 minphi = min(PHASE);
00209   const real8 maxphi = max(PHASE);
00210   normalize(PHASE,minphi,maxphi);                       // regrid data
00211 
00212   for (i=0; i<Npoints; i++)     // solve system for all points
00213     {
00214 
00215     // ______ Set up design matrix ______
00216     for (j=0; j<Nheights; j++)
00217       for (k=0; k<=DEGREE1D; k++)
00218         DESIGN(j,k) = pow(PHASE(i,j),real8(k));                 // PHASE is normalized
00219 
00220 
00221     // ______ Solve by cholesky (even the exactly determined case) ______
00222     matrix<real8> N   = matTxmat(DESIGN,DESIGN);
00223     matrix<real8> rhs = matTxmat(DESIGN,HEI);
00224     matrix<real8> Qx_hat = N;
00225     choles(Qx_hat);             // Cholesky factorisation normalmatrix
00226     solvechol(Qx_hat,rhs);      // Estimate of unknowns (alphas) in rhs
00227 
00228 
00229 //#ifdef __DEBUG        // else trust solution, errors from lib.
00230 // ______ Test inverse ______
00231     invertchol(Qx_hat);         // Covariance matrix (lower)
00232     for (k=0; k<Qx_hat.lines(); k++)
00233       for (j=0; j<k; j++)
00234         Qx_hat(j,k) = Qx_hat(k,j);              // was only stored in lower
00235     const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00236     INFO << "s2h schwaebisch: max(abs(N*inv(N)-I)) = " << maxdev;
00237     INFO.print();
00238 /*
00239     matrix<real8> unity = N * Qx_hat;
00240     real8 maxdev = abs(unity(0,0)-1);
00241     for (k=1; k<unity.lines(); k++)
00242       {
00243       if (abs(unity(k,k)-1) > maxdev)
00244         maxdev = abs(unity(k,k)-1);
00245       A
00246       for (j=0; j<k-1; j++)
00247         {
00248         if (abs(unity(k,j)) > maxdev)
00249           maxdev = abs(unity(k,j));
00250         if (abs(unity(j,k)) > maxdev)
00251           maxdev = abs(unity(j,k));
00252         }
00253       }
00254 */
00255     if (maxdev > .01)
00256       WARNING.print("wrong solution for 1d polynomial? (decrease d1d or nhei)");
00257 //#endif
00258     // ______ Scale back unknowns: alpha_i <= alpha_i * (scale)^i ______
00259     // ______ Store solution in ALPHAS ______
00260     for (alfa=0; alfa<=DEGREE1D; alfa++)
00261       ALPHAS(i,alfa) = rhs(alfa,0);
00262     } // loop over all points
00263 
00264 
00265 
00266 // ====== STEP 3 ======
00267   PROGRESS.print("S2H: schwabisch: STEP3: estimate coefficients for 2d polynomial.");
00268   // ====== Compute alpha_i coefficients of polynomials as function of (l,p) ======
00269   // ______ alpha_i = sum(k,l) beta_kl l^k p^l; ______
00270   // ______ Solve simultaneous for all betas ______
00271   // ______ this does not seem to be possibly with my routine, so do per alfa_i ______
00272   const int32 Nunk     = Ncoeffs(DEGREE2D);              // Number of unknowns
00273 
00274   // ______ Check redundancy is done before? ______
00275   if (Npoints < Nunk)
00276     {
00277     PRINT_ERROR("slant2hschwabisch: N_observations<N_unknowns (increase S2H_NPOINTS or decrease S2H_DEGREE2D.")
00278     throw(input_error);
00279     }
00280 
00281   matrix<real8>         A(Npoints,Nunk);                // designmatrix
00282 
00283 
00284   // ______ Set up system of equations ______
00285   // ______ Order unknowns: B00 B10 B01 B20 B11 B02 B30 B21 B12 B03 for degree=3 ______
00286   const real8 minL = min(Position.getcolumn(0)); 
00287   const real8 maxL = max(Position.getcolumn(0)); 
00288   const real8 minP = min(Position.getcolumn(1)); 
00289   const real8 maxP = max(Position.getcolumn(1)); 
00290 
00291   for (i=0; i<Npoints; i++)
00292     {
00293     // ______ normalize coordinates ______
00294     const real8 posL = normalize(real8(Position(i,0)),minL,maxL);
00295     const real8 posP = normalize(real8(Position(i,1)),minP,maxP);
00296    
00297     index = 0;
00298     for (j=0; j<=DEGREE2D; j++)
00299       {
00300       for (k=0; k<=j; k++)
00301         {
00302         A(i,index) = pow(posL,real8(j-k)) * pow(posP,real8(k));
00303         index++;
00304         }
00305       }
00306     }
00307 
00308   // ______ Solve 2d polynomial system for alfas at these points _____
00309   matrix<real8> N   = matTxmat(A,A);
00310   matrix<real8> rhs = matTxmat(A,ALPHAS);
00311   matrix<real8> Qx_hat = N;
00312   choles(Qx_hat);                       // Cholesky factorisation normalmatrix
00313 /////  solvechol(Qx_hat,rhs);        // Estimate of unknowns (betas) in rhs, NOT OK!
00314 
00315 
00316   // ______ Solve the normal equations for all alpha_i ______
00317   // ______ Simultaneous solution doesn't work somehow ______
00318   for (i=0; i<rhs.pixels(); ++i)
00319     {
00320     matrix<real8> rhs_alphai = rhs.getcolumn(i);
00321     solvechol(Qx_hat,rhs_alphai);       // Solution in rhs_alphai
00322     rhs.setcolumn(i,rhs_alphai);        // place solution back
00323     }
00324 
00325 
00326   // ______ Test solution by inverse ______
00327   invertchol(Qx_hat);                   // Covariance matrix (lower)
00328   for (i=0; i<Qx_hat.lines(); i++)
00329     for (j=0; j<i; j++)
00330       Qx_hat(j,i) = Qx_hat(i,j);        // was only stored in lower
00331   const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00332   INFO << "s2h schwaebisch: max(abs(N*inv(N)-I)) = " << maxdev;
00333   INFO.print();
00334 /*
00335   matrix<real8> unity = N * Qx_hat;
00336   real8 maxdev = abs(unity(0,0)-1);
00337   for (i=1; i<unity.lines(); i++)
00338     {
00339     if (abs(unity(i,i)-1) > maxdev)
00340       maxdev = abs(unity(i,i)-1);
00341     for (j=0; j<i-1; j++)
00342       {
00343       if (abs(unity(i,j)) > maxdev)
00344         maxdev = abs(unity(i,j));
00345       if (abs(unity(j,i)) > maxdev)
00346         maxdev = abs(unity(j,i));
00347       }
00348     }
00349 */
00350 
00351   if (maxdev > .01)
00352     {
00353     WARNING << "slant2h: possibly wrong solution. deviation from unity AtA*inv(AtA) = "
00354          << maxdev << " >.01";
00355     WARNING.print();
00356     }
00357 
00358 
00359 
00360 // ====== STEP 4 ======
00361   PROGRESS.print("S2H: schwabisch: STEP4: compute height for all pixels.");
00362   // ====== Evaluate for all points interferogram h=f(l,p,phase) ======
00363   // ______ recon with multilook, degree1D, degree2D free
00364   // ______ Multilook factors ______
00365   const real8 multiL = unwrappedinterf.multilookL;
00366   const real8 multiP = unwrappedinterf.multilookP;
00367 
00368 
00369   // ______ Number of lines/pixels of multilooked unwrapped interferogram ______
00370   const int32 mllines = int32(floor(real8(
00371     unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
00372   const int32 mlpixels = int32(floor(real8(
00373     unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1)   / multiP));
00374 
00375 
00376   // ______ Line/pixel of first point in original master coordinates ______
00377   const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
00378                                 (real8(multiL) - 1.) / 2.;
00379   const real8 firstpixel    = real8(unwrappedinterf.win.pixlo)  +
00380                                 (real8(multiP) - 1.) / 2.;
00381 
00382 
00383   // ______ Constant axis of pixel coordinates ______
00384   matrix<real4> p_axis(mlpixels,1);
00385   for (i=0; i<mlpixels; i++)
00386   //    p_axis(i,0) = real8((firstpixel + i*multiP) - minP) / (maxP - minP);
00387     p_axis(i,0) = firstpixel + i*multiP;
00388   normalize(p_axis,minP,maxP);
00389 
00390   const int32 NUMMAT = 1 + (1 + DEGREE1D);              // number of heavy matrices
00391   int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
00392   if (bufferlines > mllines)                            // whole image fits in BUFFER
00393     bufferlines = mllines;
00394 
00395   const int32 FULLBUFFERS = mllines / bufferlines;
00396   const int32 RESTLINES   = mllines % bufferlines;
00397   const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
00398   
00399   // ______ Window to be read into BUFFER from file in multilooked system ______
00400   const uint dummy = 999999;                    // large to force error if not ok
00401   window bufferwin(1, bufferlines, 1, mlpixels);        // initial
00402   window offsetbuffer(1,dummy,1,dummy); // dummy not used in readfile, no offset
00403 
00404   // ______ Open output file ______
00405   ofstream ofile;
00406   openfstream(ofile,slant2hinput.fohei,generalinput.overwrit);
00407   bk_assert(ofile,slant2hinput.fohei,__FILE__,__LINE__); 
00408 
00409 
00410   // ====== Process BUFFERS ======
00411   for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; buffer++)
00412     {
00413 
00414     // ______ Give progress ______
00415     PROGRESS << "SLANT2H: " 
00416          << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
00417     PROGRESS.print();
00418 
00419     // ______ In original master coordinate sytem ______
00420     const real8 firstline = veryfirstline + (buffer-1) * bufferlines * multiL;
00421 
00422     // ______ Set indices for loading / check last buffer ______
00423     bufferwin.linelo = 1 + (buffer-1) * bufferlines;    // Update window to be read from file
00424     if (buffer == FULLBUFFERS+1)
00425       {
00426       bufferlines = RESTLINES;
00427       //BUFFER.resize(bufferlines,mlpixels);
00428       }
00429     bufferwin.linehi = bufferwin.linelo + bufferlines - 1; // window 2b read from file
00430 
00431     // ______ Read in phase buffer of unwrapped interferogram ______
00432     matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
00433 
00434 
00435     // ______ Evaluate polynomial coefficients for these points ______
00436     // ______ Compute first line of current buffer in master coordinates ______
00437     matrix<real4> l_axis(bufferlines,1);
00438     for (k=0; k<bufferlines; k++)
00439       l_axis(k,0) = firstline + k*multiL;
00440     normalize(l_axis,minL,maxL);
00441 
00442     
00443 
00444 // ______ Lookup table because not known in advance what DEGREE1D is ______
00445 // BK: usage:
00446 //  pntALPHA[0]->showdata();
00447 //  (*pntALPHA[0][0]).showdata();
00448 // containing a grid of 1D coefficients
00449 
00450     // not possible??? (constant integral): matrix<real4> *pntALPHA[DEGREE1D+1];
00451     // but this is ansi c++ ?
00452     matrix<real4> *pntALPHA[TEN];
00453 
00454     for (k=0; k<=DEGREE1D; k++)
00455       {
00456       matrix<real8> beta(Ncoeffs(DEGREE2D),1);
00457       for (l=0; l<Ncoeffs(DEGREE2D); l++)
00458         {
00459         beta(l,0) = rhs(l,k);           // solution stored in rhs
00460         }
00461       pntALPHA[k] = new matrix<real4> (l_axis.lines(),p_axis.lines());
00462       (*pntALPHA[k]) = polyval(l_axis, p_axis, beta, DEGREE2D);
00463       }
00464 
00465 
00466     // ______ Matrix should be deleted ??? (later) ______
00467 
00468     // ______ Evaluate h=f(l,p,phi) for all points in grid in BUFFER ______
00469     matrix<real8> coeff_thispoint(DEGREE1D+1,1);
00470     for (register int32 line=0; line<BUFFER.lines(); line++)
00471       {
00472       for (register int32 pixel=0; pixel<BUFFER.pixels(); pixel++)
00473         {
00474 
00475         // ______ Check if unwrapped ok, else compute h ______
00476         if (BUFFER(line,pixel) != NaN)          // else leave NaN
00477           {
00478           for (k=0; k<DEGREE1D+1; k++)
00479             {
00480             coeff_thispoint(k,0) = (*pntALPHA[k])[line][pixel];
00481             }
00482           BUFFER(line,pixel) = polyval1d(normalize(
00483                     BUFFER(line,pixel),minphi,maxphi),  // regrid phi [0,1]
00484                     coeff_thispoint);
00485           }
00486         }
00487       }
00488 
00489     // ______ Write computed heights to file ______
00490     ofile << BUFFER;
00491     }                                           // loop over BUFFERS
00492   ofile.close();
00493 
00494 
00495 // ====== Write to result files ======
00496   ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
00497   bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
00498   scratchlogfile << "\n\n*******************************************************************"
00499                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
00500                  << "\n*******************************************************************"
00501                  << "\nMethod: \t\t\tschwabisch"
00502                  << "\nData_output_file: \t\t"
00503                  <<  slant2hinput.fohei
00504                  << "\nData_output_format: \t\t"
00505                  << "real4"
00506                  << "\nEllipsoid (name,a,b): \t\t"
00507                  <<  ellips.name << " " 
00508                  <<  ellips.a    << " " 
00509                  <<  ellips.b 
00510                  << endl;
00511   scratchlogfile.close();
00512 
00513 
00514   ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
00515   bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
00516   scratchresfile << "\n\n*******************************************************************"
00517                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
00518                  << "\n*******************************************************************"
00519                  << "\nMethod: \t\t\tschwabisch"
00520                  << "\nData_output_file:                     \t"
00521                  <<  slant2hinput.fohei
00522                  << "\nData_output_format:                   \t"
00523                  << "real4"
00524                  << "\nFirst_line (w.r.t. original_master):  \t"
00525                  <<  unwrappedinterf.win.linelo
00526                  << "\nLast_line (w.r.t. original_master):   \t"
00527                  <<  unwrappedinterf.win.linehi
00528                  << "\nFirst_pixel (w.r.t. original_master): \t"
00529                  <<  unwrappedinterf.win.pixlo
00530                  << "\nLast_pixel (w.r.t. original_master):  \t"
00531                  <<  unwrappedinterf.win.pixhi
00532                  << "\nMultilookfactor_azimuth_direction:    \t"
00533                  <<  unwrappedinterf.multilookL
00534                  << "\nMultilookfactor_range_direction:      \t"
00535                  <<  unwrappedinterf.multilookP
00536 
00537                  << "\nEllipsoid (name,a,b):                 \t"
00538                  <<  ellips.name << " " 
00539                  <<  ellips.a    << " " 
00540                  <<  ellips.b
00541                  << "\n*******************************************************************"
00542                  //<< "\n* End_slant2h:_NORMAL"
00543                  << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
00544                  << "\n*******************************************************************\n";
00545 
00546   scratchresfile.close();
00547   PROGRESS.print("finished slant2hschwabisch.");
00548   } // END slant2h
00549 
00550 
00551 
00552 /****************************************************************
00553  *    slant2h ambiguity method                                  *
00554  *                                                              *
00555  * compute height in radar coded system (master):               *
00556  * use Bhor Bver                                                *
00557  * (constant per line for parallel orbits, otherwize nearly)    *
00558  * and transformation model to get slave position               *
00559  *                                                              *
00560  * Input:                                                       *
00561  *  -                                                           *
00562  * Output:                                                      *
00563  *  -                                                           *
00564  *                                                              *
00565  *    Bert Kampes, 07-Jul-1999                                  *
00566  ****************************************************************/
00567 void slant2hambiguity(
00568         const input_gen     &generalinput,
00569         const input_slant2h &slant2hinput,
00570         const input_ell     &ellips,
00571         const slcimage      &master,
00572         const slcimage      &slave,
00573         const productinfo   &unwrappedinterf,
00574         orbit               &masterorbit,
00575         orbit               &slaveorbit,
00576         const BASELINE      &baseline)
00577   {
00578   TRACE_FUNCTION("slant2hambiguity (BK 07-Jul-1999)")
00579 
00580   const int32 MAXITER   = 10;                   // iterations for lp2xyz
00581   const real8 CRITERPOS = 1e-6;                 // stop criterium for lp2xyz
00582   const real8 CRITERTIM = 1e-10;                // stop criterium for lp2xyz
00583 
00584   const int32 MAXITERHERE = 10;                 // iterations for h
00585   const real8 CRITERHERE  = 0.05;               // m (delta h in iterations)
00586 
00587 
00588   // ______ Multilook factors ______
00589   const real8 multiL = unwrappedinterf.multilookL;
00590   const real8 multiP = unwrappedinterf.multilookP;
00591 
00592   // ______ Number of lines/pixels of multilooked unwrapped interferogram ______
00593   const int32 mllines = int32(floor(real8(
00594     unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
00595   const int32 mlpixels = int32(floor(real8(
00596     unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1)   / multiP));
00597 
00598   // ______ Line/pixel of first point in original master coordinates ______
00599   const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
00600                                 (real8(multiL) - 1.) / 2.;
00601   const real8 firstpixel    = real8(unwrappedinterf.win.pixlo)  +
00602                                 (real8(multiP) - 1.) / 2.;
00603 
00604 
00605   // ====== Compute number of buffers required ======
00606   // BK 8/10/99: why use buffer? just read in line by line is as efficient...
00607   const int32 NUMMAT = 3;       // number of large matrices BUFFER PHI LAMBDA
00608   int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
00609   if (bufferlines > mllines)                            // whole image fits in BUFFER
00610     bufferlines = mllines;
00611 
00612   const int32 FULLBUFFERS = mllines / bufferlines;
00613   const int32 RESTLINES   = mllines % bufferlines;
00614   const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
00615   
00616 
00617 // ______ Window to be read into BUFFER from file in multilooked system ______
00618 //  matrix<real4> BUFFER(bufferlines,mlpixels);         // unwrapped phase and height
00619 //  matrix<real4> PHI(BUFFER.lines(),BUFFER.pixels());  // geocoded
00620 //  matrix<real4> LAMBDA(BUFFER.lines(),BUFFER.pixels());       // geocoded
00621 
00622   window bufferwin(1, bufferlines, 1, mlpixels);        // initial
00623 
00624 
00625   // ______ Open output file ______
00626   ofstream fohei;
00627   openfstream(fohei,slant2hinput.fohei,generalinput.overwrit);
00628   bk_assert(fohei,slant2hinput.fohei,__FILE__,__LINE__); 
00629 
00630   ofstream fophi;
00631   openfstream(fophi,slant2hinput.fophi,generalinput.overwrit);
00632   bk_assert(fophi,slant2hinput.fophi,__FILE__,__LINE__); 
00633 
00634   ofstream folambda;
00635   openfstream(folambda,slant2hinput.folam,generalinput.overwrit);
00636   bk_assert(folambda,slant2hinput.folam,__FILE__,__LINE__); 
00637 
00638   // ______ Local variables ______
00639   cn M;                                 // coordinates of master in orbit
00640   cn Mdot;                              // velocity of master in orbit
00641   cn S;                                 // coordinates of slave in orbit
00642   cn P;                                 // coordinates of point on earth
00643   real8 sintheta = 0.0;                 // sine looking angle
00644   real8 costheta = 0.0;                 // cosine looking angle
00645   real8 inc_angle = 0.0;                // incidence angle by KKM
00646   matrix<real8> observations(3,1);      // setofeq
00647   matrix<real8> partials(3,3);          // partial derivatives
00648   matrix<real8> solution(3,1);          // solution of setofeq.
00649 
00650   
00651   // ====== Process BUFFERS ======
00652   for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; ++buffer)
00653     {
00654     // ______ Give progress ______
00655     PROGRESS << "SLANT2H: " 
00656          << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
00657     PROGRESS.print();
00658 
00659     // ====== First read in data ======
00660     // ______ firstline of this buffer in master coordinate system ______
00661     real8 firstline = veryfirstline + real8((buffer-1) * bufferlines * multiL); 
00662 
00663     // ______ Set indices to be read from file / check if last buffer ______
00664     bufferwin.linelo = 1 + (buffer-1) * bufferlines;
00665     if (buffer == FULLBUFFERS+1)
00666       {
00667       bufferlines = RESTLINES;
00668       }
00669     bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
00670 
00671     // ______ Read in buffer of unwrapped interferogram ______
00672     matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
00673     matrix<real4> PHI(BUFFER.lines(),BUFFER.pixels());          // geocoded
00674     matrix<real4> LAMBDA(BUFFER.lines(),BUFFER.pixels());       // geocoded
00675 
00676 
00677     // ====== Actually compute h for all points ======
00678     // ====== (better use inverse function in baseline, ie., schwabisch method.)
00679     input_ell ELLIPS;                                   // to correct height of ellips
00680     register int32 i,j;
00681     register real8 line = firstline - multiL;           // in master coordinate system
00682     for (i=0; i<BUFFER.lines(); i++)
00683       {
00684       line                += multiL;                    // in master coordinate system
00685       real8 currentheight  = 0.;                        // this iteration
00686       real8 lastheight     = 0.;                        // last iteration
00687       register real8 pixel = firstpixel - multiP;       // in master coordinate system
00688 
00689 
00690       // ______ Evaluate position M,S,P(ell(h)) for this pixel(l,p) ______
00691       // ______ only dependent on line ______
00692       const real8 m_tazi = master.line2ta(line);
00693       M    = masterorbit.getxyz(m_tazi);
00694       Mdot = masterorbit.getxyzdot(m_tazi);
00695       const real8 norm2M = M.norm2();                   // (sqr) constant per line
00696 
00697       // ______ compute baseline for point on ellips (h) ______
00698       // ______ Compute P(x,y,z), also iteratively ______
00699       //      real8 middlepixel = pixel+(.5*real8(BUFFER.pixels())*multiP);
00700       real8 tempphase = 0.0;    // to find out baseline parameters
00701       real8 temppixel = 0.0;  // to find out baseline parameters
00702       real8 Bpar      = 0.0;
00703       real8 Bperp     = 0.0;
00704 
00705       bool lineokunwrapped = false;
00706       for (int t1=0; t1<BUFFER.pixels(); ++t1)
00707         {
00708         if (BUFFER(i,t1) != NaN)
00709           lineokunwrapped = true;
00710         }
00711 
00712       if (lineokunwrapped)
00713         {
00714         // ______ Get phase for a pixel to compute h to obtain baseline parameters ______
00715         int32 middle = BUFFER.pixels()/2;       // floor
00716         for (j=0; j<middle+1; ++j)
00717           {
00718           if (BUFFER(i,middle-j) != NaN)
00719             {
00720             tempphase = BUFFER(i,middle-j);
00721             temppixel = firstpixel + (middle-j)*multiP;
00722             break; // for loop
00723             }
00724           else if (BUFFER(i,middle+j) != NaN)
00725             {
00726             tempphase = BUFFER(i,middle+j);
00727             temppixel = firstpixel + (middle+j)*multiP;
00728             break; // for loop
00729             }
00730           }
00731 
00732 //for loop.... to find out h of point for correct baseline computation
00733 //might be done by transformation model as well.(better?)
00734 // ______ Compute h iteratively to get ok baseline parameters ______
00735         for (j=0; j<=MAXITERHERE; ++j)
00736           {
00737 
00738 bool DO_NEW_METHOD=false;// no time to test this, methdo seems wrong??
00739 if (DO_NEW_METHOD==false)// old method, not using BASELINE class:
00740 {
00741           real8 s_tazi;                                 // returned
00742           real8 s_trange;                               // returned, unused?
00743           ELLIPS.a = ellips.a + currentheight;          //next height
00744           ELLIPS.b = ellips.b + currentheight;          //next height
00745           lp2xyz(line,temppixel,
00746                  ELLIPS,master,masterorbit,
00747                  P,MAXITER,CRITERPOS);                  // P returned
00748 
00749 // ______ Compute S(x,y,z) ______ 
00750           xyz2t(s_tazi,s_trange,slave,slaveorbit,
00751                 P,MAXITER,CRITERTIM);
00752           S = slaveorbit.getxyz(s_tazi);
00753 
00754 // ====== The baseline parameters, derived from the positions (x,y,z) ======
00755 // ====== Compute Bhor Bver (assumed contant per line) ======
00756 // ______ theta is angle (M,M-P) ______
00757           const real8 B = M.dist(S);                    // abs. value
00758           Bpar = P.dist(M) - P.dist(S); // sign ok
00759 
00760 
00761 // ______ if (MP>SP) then S is to the right of slant line, then B perp is positive.
00762 //          costheta = ((-(P.norm2()) + norm2M +                // cosine law
00763 //            sqr(pix2range(temppixel,master.t_range1,master.rsr2x))) /
00764 //           (2*sqrt(norm2M)*pix2range(temppixel,master.t_range1,master.rsr2x)));
00765           costheta = ((-(P.norm2()) + norm2M +          // cosine law
00766               sqr(master.pix2range(temppixel))) /
00767              (2*sqrt(norm2M)*master.pix2range(temppixel)));
00768           sintheta = sqrt(1-sqr(costheta));             // cos^2 + sin^2 = 1
00769  
00770           const cn r2 = S.min(P);                       // vector; to find sign
00771           //Bperp = (acos(costheta) > M.angle(r2)) ?    // sign ok
00772           //   sqrt(sqr(B)-sqr(Bpar)) :
00773           //  -sqrt(sqr(B)-sqr(Bpar));
00774 
00775           // ___ following is less comp. intens.
00776           // (cost1>cost2) then t1<t2 then Bperp<0
00777           // BK 23-Oct-2000
00778           const real8 costheta2 = M.in(r2) / (M.norm()*r2.norm());
00779           Bperp = (costheta < costheta2) ?              // sign ok
00780              sqrt(sqr(B)-sqr(Bpar)) :
00781             -sqrt(sqr(B)-sqr(Bpar));
00782           // ____
00783 }
00784 else
00785 {
00786 // --- New method with BASELINE class and incidence angle ---
00787 Bpar     = baseline.get_bpar(line,temppixel,currentheight);
00788 Bperp    = baseline.get_bperp(line,temppixel,currentheight);
00789 costheta = cos(baseline.get_theta_inc(line,temppixel,currentheight));
00790 sintheta = sqrt(1-sqr(costheta));
00791 inc_angle = baseline.get_theta_inc(line,temppixel,currentheight);//KKM
00792 }
00793 
00794           // --- Update height ---
00795           lastheight         = currentheight;
00796           const real8 m_tr   = master.pix2tr(temppixel);
00797           const real8 tempr1 = SOL*m_tr;
00798           //currentheight      = (-master.wavelength*tempr1*sintheta*tempphase)/
00799                                //(4.*PI*Bperp); KKM commented the old code
00800           currentheight      = (-master.wavelength*tempr1*sin(inc_angle)*tempphase)/
00801                                (4.*PI*Bperp);//KKM added this
00802 
00803 
00804           // ______ Check convergence ______
00805           if (abs(currentheight-lastheight) < CRITERHERE)
00806             break; // iterate to get h
00807           }                             // loop iterations (j)
00808 
00809         // ______ Check number of iterations ______
00810         if (j >= MAXITERHERE)
00811           {
00812           WARNING << "slant2hambiguity: maxiter reached. "
00813                << "MAXITER: " << MAXITERHERE
00814                << "CRITER: "  << CRITERHERE << "m "
00815                << "last correction: " << currentheight-lastheight;
00816           WARNING.print();
00817           }
00818         } // check lineokunwrapped (all NaNs)
00819 
00820 // ______ The baseline parameters that are used foreach pixel ______
00821 // this is not correct, Bhor not same for each pixel? 
00822 // BK 09-Aug-2000
00823       const real8 Bhor  = Bperp*costheta + Bpar*sintheta;
00824       const real8 Bver  = Bperp*sintheta - Bpar*costheta;
00825       currentheight     = 0.;                   // this iteration
00826       lastheight        = 0.;                   // last iteration
00827 
00828 
00829 // ====== Start loop over pixels ======
00830       for (j=0; j<BUFFER.pixels(); j++)
00831         {
00832         pixel += multiP;                        // in master coordinate system
00833 
00834         // ______ Check if conversion is necessary _______
00835         if (BUFFER(i,j) == NaN)                 // leave NaN in buffer
00836           {
00837           PHI(i,j)    = NaN;                    // not geocoded
00838           LAMBDA(i,j) = NaN;                    // not geocoded
00839           }
00840 
00841         else
00842           {
00843 
00844           // ______ Compute some constants (per pixel) ______
00845           const real8 m_trange       = master.pix2tr(pixel);
00846           const real8 normr1         = SOL*m_trange;
00847           const real8 partnumerator  = norm2M + sqr(normr1);
00848           const real8 denominator    = 2.0*sqrt(norm2M)*normr1;
00849           const real8 m_lamr1phidiv4pi = master.wavelength*normr1*BUFFER(i,j) / (4.0*PI);
00850 
00851           // ______ Iteratively solve for height ______
00852           register int32 iteration;
00853           for (iteration=0; iteration<=MAXITERHERE; ++iteration)
00854             {
00855             ELLIPS.a = ellips.a + currentheight;//next height
00856             ELLIPS.b = ellips.b + currentheight;//next height
00857             //ELLIPS.ecc1st_sqr();              // not used in lp2xyz
00858             //ELLIPS.ecc2nd_sqr();              // not used in lp2xyz
00859       
00860 
00861 
00862             // ______ Evaluate position P(for l,p) at h ______
00863             // use function in orbitbk.cc cause eq1 there private f
00864             // BK 09-Aug-2000
00865             lp2xyz(line,pixel,
00866                    ELLIPS,master,masterorbit,
00867                    P,MAXITER,CRITERPOS);                        // P returned
00868 
00869 
00870             // ______ Compute h(theta,B,phi) ______
00871             // ______ dh=-lambdaover4pi*R1*(sintheta/(Bhor*costheta+Bver*sintheta))*(phi1-phi2)
00872             costheta      = (-(P.norm2()) + partnumerator) / denominator;       // cosine law
00873             sintheta      = sqrt(1-sqr(costheta));              // cos^2 + sin^2 = 1
00874             lastheight    = currentheight;
00875             //currentheight = (-m_lamr1phidiv4pi * sintheta) / //KKM commented this existing line
00876                             //(Bhor*costheta + Bver*sintheta); //KKM commented this existing line
00877             inc_angle = baseline.get_theta_inc(line,pixel,currentheight);//KKM
00878             currentheight = (-m_lamr1phidiv4pi * sin(inc_angle)) /
00879                             (Bhor*costheta + Bver*sintheta);//KKM
00880 
00881 
00882 
00883             // ______ Check convergence ______
00884             if (abs(currentheight-lastheight) < CRITERHERE)
00885               break;
00886             }                           // loop iterations
00887 
00888           // ______ Check number of iterations ______
00889           if (iteration >= MAXITERHERE)
00890             {
00891             WARNING << "slant2hambiguity: maxiter reached. "
00892                  << "MAXITER: " << MAXITERHERE
00893                  << "CRITER: "  << CRITERHERE << "m "
00894                  << "last correction: " << currentheight-lastheight;
00895             WARNING.print();
00896             }
00897 
00898 
00899           // ______ Put computed height in buffer ______
00900           BUFFER(i,j) = currentheight;
00901 
00902           // ______ Geocode P(x,y,z) ______
00903           ELLIPS.a = ellips.a + currentheight;
00904           ELLIPS.b = ellips.b + currentheight;
00905           ELLIPS.ecc1st_sqr();
00906           ELLIPS.ecc2nd_sqr();
00907           const real8 r    = sqrt(sqr(P.x)+sqr(P.y));
00908           const real8 nu   = atan2((P.z*ELLIPS.a),(r*ELLIPS.b));
00909           const real8 sin3 = pow(sin(nu),3);
00910           const real8 cos3 = pow(cos(nu),3);
00911           PHI(i,j)         = rad2deg(atan2((P.z+ELLIPS.e2b*ELLIPS.b*sin3),
00912                                            (r-ELLIPS.e2*ELLIPS.a*cos3)));
00913           LAMBDA(i,j)      = rad2deg(atan2(P.y,P.x));
00914 
00915           }                             // unwrapped ok
00916         }                               // loop over pixels
00917       }                                 // loop over lines
00918     fohei << BUFFER;                    // write output
00919     fophi << PHI;                       // write output
00920     folambda << LAMBDA;                 // write output
00921     }                                   // loop over buffers
00922 
00923 
00924 // ====== Write to result files ======
00925   ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
00926   bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
00927   scratchlogfile << "\n\n*******************************************************************"
00928                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
00929                  << "\n*******************************************************************"
00930                  << "\nMethod: \t\t\tambiguity"
00931                  << "\nData_output_file: \t\t"
00932                  <<  slant2hinput.fohei
00933                  << "\nData_output_format: \t\t"
00934                  << "real4"
00935                  << "\nData_output_file_phi: \t\t"
00936                  <<  slant2hinput.fophi
00937                  << "\nData_output_format: \t\t"
00938                  << "real4"
00939                  << "\nData_output_file_lam: \t\t"
00940                  <<  slant2hinput.folam
00941                  << "\nData_output_format: \t\t"
00942                  << "real4"
00943                  << "\nEllipsoid (name,a,b): \t\t"
00944                  <<  ellips.name << " " 
00945                  <<  ellips.a    << " " 
00946                  <<  ellips.b 
00947                  << endl;
00948   scratchlogfile.close();
00949 
00950 
00951   ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
00952   bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
00953   scratchresfile << "\n\n*******************************************************************"
00954                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
00955                  << "\n*******************************************************************"
00956                  << "\nMethod:                               \t"
00957                  << "ambiguity"
00958                  << "\nData_output_file:                     \t"
00959                  <<  slant2hinput.fohei
00960                  << "\nData_output_format:                   \t"
00961                  << "real4"
00962                  << "\nData_output_file_phi:                 \t"
00963                  <<  slant2hinput.fophi
00964                  << "\nData_output_format:                   \t"
00965                  << "real4"
00966                  << "\nData_output_file_lam:                 \t"
00967                  <<  slant2hinput.folam
00968                  << "\nData_output_format:                   \t"
00969                  << "real4"
00970                  << "\nFirst_line (w.r.t. original_master):  \t"
00971                  <<  unwrappedinterf.win.linelo
00972                  << "\nLast_line (w.r.t. original_master):   \t"
00973                  <<  unwrappedinterf.win.linehi
00974                  << "\nFirst_pixel (w.r.t. original_master): \t"
00975                  <<  unwrappedinterf.win.pixlo
00976                  << "\nLast_pixel (w.r.t. original_master):  \t"
00977                  <<  unwrappedinterf.win.pixhi
00978                  << "\nMultilookfactor_azimuth_direction:    \t"
00979                  <<  unwrappedinterf.multilookL
00980                  << "\nMultilookfactor_range_direction:      \t"
00981                  <<  unwrappedinterf.multilookP
00982 
00983                  << "\nEllipsoid (name,a,b):                 \t"
00984                  <<  ellips.name << " " 
00985                  <<  ellips.a    << " " 
00986                  <<  ellips.b
00987                  << "\n*******************************************************************"
00988                  //<< "\n* End_slant2h:_NORMAL"
00989                  << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
00990                  << "\n*******************************************************************\n";
00991 
00992 // ====== Tidy up ======
00993   scratchresfile.close();
00994   PROGRESS.print("finished slant2hambiguity.");
00995   } // END slant2hambiguity
00996 
00997 
00998 
00999 /****************************************************************
01000  *    slant2h rodriguez92 method                                *
01001  *                                                              *
01002  * compute height in radar coded system (master):               *
01003  * use own derivations, check carefully                         *
01004  *                                                              *
01005  * Input:                                                       *
01006  *  -                                                           *
01007  * Output:                                                      *
01008  *  -                                                           *
01009  *                                                              *
01010  *    Bert Kampes, 30-Sep-1999                                  *
01011  ****************************************************************/
01012 void slant2hrodriguez(
01013         const input_gen     &generalinput,
01014         const input_slant2h &slant2hinput,
01015         const input_ell     &ellips,
01016         const slcimage      &master,
01017         const slcimage      &slave,
01018         const productinfo   &unwrappedinterf,
01019         const matrix<real8> &coeff_flatearth,
01020         orbit               &masterorbit,
01021         orbit               &slaveorbit,
01022         const BASELINE      &baseline)
01023   {
01024   TRACE_FUNCTION("slant2hrodriguez (BK 30-Sep-1999)")
01025   WARNING.print("this method is based on wrong approximations. (?)");
01026 
01027   const int32 MAXITER   = 10;                   // iterations for lp2xyz
01028   const real8 CRITERPOS = 1e-6;                 // stop criterium for lp2xyz
01029   const real8 CRITERTIM = 1e-10;                // stop criterium for lp2xyz
01030 
01031   const int32 degreecfe = degree(coeff_flatearth.size());
01032 
01033 // ______ Normalize data for polynomial ______
01034   const real8 minL = master.originalwindow.linelo;
01035   const real8 maxL = master.originalwindow.linehi;
01036   const real8 minP = master.originalwindow.pixlo;
01037   const real8 maxP = master.originalwindow.pixhi;
01038 
01039 
01040   // ______ Multilook factors ______
01041   const real8 multiL = unwrappedinterf.multilookL;
01042   const real8 multiP = unwrappedinterf.multilookP;
01043   const real8 m_pi4divlambda = (4.*PI)/master.wavelength;
01044 
01045   // ______ Number of lines/pixels of multilooked unwrapped interferogram ______
01046   const int32 mllines = int32(floor(real8(
01047     unwrappedinterf.win.linehi-unwrappedinterf.win.linelo+1) / multiL));
01048   const int32 mlpixels = int32(floor(real8(
01049     unwrappedinterf.win.pixhi-unwrappedinterf.win.pixlo+1)   / multiP));
01050 
01051   // ______ Line/pixel of first point in original master coordinates ______
01052   const real8 veryfirstline = real8(unwrappedinterf.win.linelo) +
01053                                 (real8(multiL) - 1.) / 2.;
01054   const real8 firstpixel    = real8(unwrappedinterf.win.pixlo)  +
01055                                 (real8(multiP) - 1.) / 2.;
01056 
01057 
01058   // ====== Compute number of buffers required ======
01059   const int32 NUMMAT = 1;       // number of large matrices to determine size of matrices
01060   int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
01061   if (bufferlines > mllines)                            // whole image fits in BUFFER
01062     bufferlines = mllines;
01063 
01064   const int32 FULLBUFFERS = mllines / bufferlines;
01065   const int32 RESTLINES   = mllines % bufferlines;
01066   const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
01067   
01068 
01069   // ______ Window to be read into BUFFER from file in multilooked system ______
01070   const uint dummy = 999999;                    // large to force error if not ok
01071   window bufferwin(1, bufferlines, 1, mlpixels);        // initial
01072   window offsetbuffer(1,dummy,1,dummy); // dummy not used in readfile, no offset
01073 
01074 
01075   // ______ Open output file ______
01076   ofstream fohei;
01077   openfstream(fohei,slant2hinput.fohei,generalinput.overwrit);
01078   bk_assert(fohei,slant2hinput.fohei,__FILE__,__LINE__); 
01079 
01080 
01081 
01082   
01083 // ====== Process BUFFERS ======
01084   for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; ++buffer)
01085     {
01086     // ______ Give progress ______
01087     PROGRESS << "SLANT2H: " 
01088          << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
01089     PROGRESS.print();
01090 
01091     // ====== First read in data ======
01092     // ______ firstline of this buffer in master coordinate system ______
01093     real8 firstline = veryfirstline + real8((buffer-1) * bufferlines * multiL); 
01094 
01095     // ______ Set indices to be read from file / check if last buffer ______
01096     bufferwin.linelo = 1 + (buffer-1) * bufferlines;
01097     if (buffer == FULLBUFFERS+1)
01098       {
01099       bufferlines = RESTLINES;
01100       //BUFFER.resize(bufferlines,mlpixels);
01101       }
01102     bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
01103 
01104     // ______ Read in buffer of unwrapped interferogram ______
01105     matrix<real4> BUFFER = unwrappedinterf.readphase(bufferwin);
01106 
01107 
01108     // ====== Actually compute h for all points ======
01109     register int32 i,j;
01110     register real8 line = firstline - multiL;           // in master coordinate system
01111     for (i=0; i<BUFFER.lines(); i++)
01112       {
01113       line += multiL;                                   // in master coordinate system
01114       register real8 pixel = firstpixel - multiP;       // in master coordinate system
01115 
01116 
01117       // ______ Evaluate position M,S,P(ell(h)) for this pixel(l,p) ______
01118       // ______ only dependent on line ______
01119       const real8 m_tazi = master.line2ta(line);
01120       cn M = masterorbit.getxyz(m_tazi);
01121       cn P;                                     // coordinates of point on earth
01122 
01123 
01124       // ______ compute baseline for point on ellips (h) ______
01125       // ______ Compute P(x,y,z) only to get baseline and alpha for this line ______
01126       const real8 middlepixel = pixel+(.5*BUFFER.pixels()*multiP);
01127       lp2xyz(line,middlepixel,
01128              ellips,master,masterorbit,
01129              P,MAXITER,CRITERPOS);              // P returned
01130 
01131 
01132       // ====== Compute which azimuth time this is for slave ======
01133       // ______ and compute S(x,y,z) for this time ______
01134       // ______ Compute S(x,y,z) could be by transf. model coregistration ______ 
01135       real8 s_tazi;                             // returned
01136       real8 s_trange;                           // returned, unused?
01137       xyz2t(s_tazi,s_trange,slave, slaveorbit,
01138             P,MAXITER,CRITERTIM);
01139       cn S = slaveorbit.getxyz(s_tazi);
01140 
01141 
01142 // !! this position is not correct, depends on P(h) if orbits are not parallel.
01143 // you should compute S by transformation model?
01144 // error is probably very small and since H is not compute exact I will leave it for now.
01145 // M -> P(h0) -> S -> B -> h1 = H - r costheta;
01146 //   -> P(h1) -> S -> B -> h2
01147 //   -> P(h2) -> S -> B -> h3 etc.
01148 
01149 
01150       // ====== The baseline parameters, derived from the positions (x,y,z) ======
01151       // ====== Compute B and alpha (contant per line) ======
01152       // ______ theta is angle (M,M-P) ______
01153       const real8 B        = M.dist(S);                         // abs. value
01154       const real8 Bpartemp = P.dist(M) - P.dist(S);             // sign ok
01155 
01156       // ______ if (MP>SP) then S is to the right of slant line, then B perp is positive.
01157       const real8 rho1sqr = M.norm2();                          // (sqr) constant per line
01158       const real8 costhetatemp = ((-(P.norm2()) + rho1sqr +             // cosine law
01159                    sqr(master.pix2range(middlepixel))) /
01160                    (2*sqrt(rho1sqr)*master.pix2range(middlepixel)));
01161       //sintheta = sqrt(1-sqr(costhetatemp));           // cos^2 + sin^2 = 1
01162 
01163       const cn r2 = S.min(P);                                   // vector; to find out sign
01164       const real8 Bperp = (acos(costhetatemp) > M.angle(r2)) ?  // sign ok
01165          sqrt(sqr(B)-sqr(Bpartemp)) :
01166         -sqrt(sqr(B)-sqr(Bpartemp));
01167 
01168       const real8 alpha  = acos(costhetatemp) - atan2(Bpartemp,Bperp);
01169 
01170       // ______ not used here ______
01171       //  const real8 Bhor  = Bperp*costhetatemp + Bpartemp*sintheta;
01172       //  const real8 Bver  = Bperp*sintheta - Bpartemp*costhetatemp;
01173       //  const real8 Bhor  = B * cos(alfa);
01174       //  const real8 Bver  = B * sin(alfa);
01175 
01176 
01177       // ______ Find out height of satellite ______
01178       // ______ assume radius to P is equal to this Radius. ______
01179       const real8 rho1  = sqrt(rho1sqr);
01180       real8 satphi,satlambda,satheight;
01181       xyz2ell(ellips,M,satphi,satlambda,satheight);
01182       const real8 Radius = rho1 - satheight;
01183 
01184 
01185       // ====== Compute per pixel: phi->Bpar, r, ->theta,p,mu,H,h ======
01186       for (j=0; j<BUFFER.pixels(); j++)
01187         {
01188         pixel += multiP;                        // in master coordinate system
01189 
01190         // ______ Check if conversion is necessary _______
01191         if (BUFFER(i,j) != NaN)                 // leave NaN in buffer
01192           {
01193 
01194           // ______ Compute some constants (per pixel) ______
01195           const real8 r1    = master.pix2range(pixel);
01196           const real8 r1sqr = sqr(r1);
01197           // ______ Compute reference phase to obtain Bpar ______
01198           // const real8 phiref= polyval(line,pixel,coeff_flatearth,degreecfe);
01199           const real8 phiref= polyval(normalize(line,minL,maxL),
01200                                       normalize(pixel,minP,maxP),
01201                                       coeff_flatearth,degreecfe);
01202           //const real8 Bpar  = -(BUFFER(i,j)+phiref)/pi4divlambda;
01203           // master [or slave] wavelength???
01204           const real8 Bpar  = -(BUFFER(i,j)+phiref)/m_pi4divlambda;
01205 
01206           const real8 sinthetaminalpha =
01207             (sqr(r1)+sqr(B)-sqr(r1-Bpar))/(2*r1*B);
01208 
01209 // ______ two possible solutions, find out later how to do this more efficient ______
01210 // ______ probably if Bperp is <0 then choose pi-theta
01211 // ______ for now use theta is approx 20 degrees., (theta=1)=57 degrees
01212           real8 theta = asin(sinthetaminalpha)+alpha;
01213           if (theta<0.0 || theta>1.0)
01214             theta = PI-asin(sinthetaminalpha)+alpha;
01215 
01216           const real8 costheta = cos(theta);
01217           const real8 psqr  = rho1sqr+r1sqr - 2*rho1*r1*costheta;       // cosine law
01218           const real8 p     = sqrt(psqr);
01219           const real8 cosmu = ((-r1sqr+psqr+rho1sqr) / (2*p*rho1));     // cosine law
01220           const real8 H = rho1 - Radius*cosmu;
01221 
01222 
01223           // ______ Put computed height in buffer ______
01224           BUFFER(i,j) = H - r1*costheta;
01225           }                             // unwrapped ok
01226         }                               // loop over pixels
01227       }                                 // loop over lines
01228     fohei << BUFFER;                    // write output
01229     }                                   // loop over buffers
01230 
01231 
01232 // ====== Write to result files ======
01233   ofstream scratchlogfile("scratchlogslant2h", ios::out | ios::trunc);
01234   bk_assert(scratchlogfile,"slant2h: scratchlogslant2h",__FILE__,__LINE__);
01235   scratchlogfile << "\n\n*******************************************************************"
01236                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
01237                  << "\n*******************************************************************"
01238                  << "\nMethod: \t\t\trodriguez"
01239                  << "\nData_output_file: \t\t"
01240                  <<  slant2hinput.fohei
01241                  << "\nData_output_format: \t\t"
01242                  << "real4"
01243                  << "\nEllipsoid (name,a,b): \t\t"
01244                  <<  ellips.name << " " 
01245                  <<  ellips.a    << " " 
01246                  <<  ellips.b 
01247                  << endl;
01248   scratchlogfile.close();
01249 
01250 
01251   ofstream scratchresfile("scratchresslant2h", ios::out | ios::trunc);
01252   bk_assert(scratchresfile,"slant2h: scratchresslant2h",__FILE__,__LINE__);
01253   scratchresfile << "\n\n*******************************************************************"
01254                  << "\n*_Start_" << processcontrol[pr_i_slant2h]
01255                  << "\n*******************************************************************"
01256                  << "\nMethod:                               \t"
01257                  << "rodriguez"
01258                  << "\nData_output_file:                     \t"
01259                  <<  slant2hinput.fohei
01260                  << "\nData_output_format:                   \t"
01261                  << "real4"
01262                  << "\nFirst_line (w.r.t. original_master):  \t"
01263                  <<  unwrappedinterf.win.linelo
01264                  << "\nLast_line (w.r.t. original_master):   \t"
01265                  <<  unwrappedinterf.win.linehi
01266                  << "\nFirst_pixel (w.r.t. original_master): \t"
01267                  <<  unwrappedinterf.win.pixlo
01268                  << "\nLast_pixel (w.r.t. original_master):  \t"
01269                  <<  unwrappedinterf.win.pixhi
01270                  << "\nMultilookfactor_azimuth_direction:    \t"
01271                  <<  unwrappedinterf.multilookL
01272                  << "\nMultilookfactor_range_direction:      \t"
01273                  <<  unwrappedinterf.multilookP
01274 
01275                  << "\nEllipsoid (name,a,b):                 \t"
01276                  <<  ellips.name << " " 
01277                  <<  ellips.a    << " " 
01278                  <<  ellips.b
01279                  << "\n*******************************************************************"
01280                  //<< "\n* End_slant2h:_NORMAL"
01281                  << "\n* End_" << processcontrol[pr_i_slant2h] << "_NORMAL"
01282                  << "\n*******************************************************************\n";
01283 
01284 
01285 // ====== Tidy up ======
01286   scratchresfile.close();
01287   PROGRESS.print("finished slant2hrodriguez.");
01288   } // END slant2rodriguez
01289 
01290 
01291 
01292 /****************************************************************
01293  *    geocode                                                   *
01294  *                                                              *
01295  * compute phi, lambda, (height is input)                       *
01296  *                                                              *
01297  * Input:                                                       *
01298  *  - slant2h done, h filename in interferogram struct          *
01299  * Output:                                                      *
01300  *  - geocoded image                                            *
01301  *                                                              *
01302  * See thesis swabisch for method 3eq.                          *
01303  *  ellips is added known height,                               *
01304  *  point(x,y,z) is evaluated at that ellips.                   *
01305  * This is converted to ell. coord. (bowrings method)           *
01306  *                                                              *
01307  *    Bert Kampes, 02-Jun-1999                                  *
01308  ****************************************************************/
01309 void geocode(
01310         const input_gen     &generalinput,
01311         const input_geocode &geocodeinput,
01312         const input_ell     &ellips,
01313         const slcimage      &master,
01314         const productinfo   &heightinradarsystem,
01315         orbit               &masterorbit)
01316   {
01317   TRACE_FUNCTION("geocode (BK 02-Jun-1999)")
01318   const int32 MAXITER   = 10;
01319   const real8 CRITERPOS = 1e-6;
01320   const real8 CRITERTIM = 1e-10;
01321 
01322 
01323 // ====== Evaluate for all points interferogram pos=f(l,p,height) ======
01324 // ______ recon with multilook, buffers
01325 // ______ Multilook factors ______
01326   const real8 multiL = heightinradarsystem.multilookL;
01327   const real8 multiP = heightinradarsystem.multilookP;
01328 
01329 // ______ Number of lines/pixels of multilooked unwrapped interferogram ______
01330   const int32 mllines = int32(floor(real8(
01331     heightinradarsystem.win.linehi-heightinradarsystem.win.linelo+1) / multiL));
01332   const int32 mlpixels = int32(floor(real8(
01333     heightinradarsystem.win.pixhi-heightinradarsystem.win.pixlo+1)   / multiP));
01334 
01335 // ______ Line/pixel of first point in original master coordinates ______
01336   const real8 veryfirstline = real8(heightinradarsystem.win.linelo) +
01337                                 (multiL - 1.) / 2.;
01338   const real8 firstpixel    = real8(heightinradarsystem.win.pixlo)  +
01339                                 (multiP - 1.) / 2.;
01340 
01341 
01342 // ====== Compute number of buffers required ======
01343   const int32 NUMMAT = 3;       // number of large matrices to determine size of matrices
01344   int32 bufferlines = generalinput.memory / (NUMMAT * (mlpixels * sizeof(real4)));
01345   if (bufferlines > mllines)                            // whole image fits in BUFFER
01346     bufferlines = mllines;
01347 
01348   const int32 FULLBUFFERS = mllines / bufferlines;
01349   const int32 RESTLINES   = mllines % bufferlines;
01350   const int32 EXTRABUFFER = RESTLINES ? 1 : 0;
01351   
01352 
01353 // ______ Window to be read into BUFFER from file in multilooked system ______
01354   matrix<real4> PHI(bufferlines,mlpixels);              // 
01355   matrix<real4> LAMBDA(bufferlines,mlpixels);           // 
01356   matrix<real4> HEIGHT(bufferlines,mlpixels);           // also output
01357   const uint dummy = 999999;                    // large to force error if not ok
01358   window bufferwin(1, bufferlines, 1, mlpixels);        // initial
01359   window offsetbuffer(1,dummy,1,dummy); // dummy not used in readfile, no offset
01360 
01361 
01362 // ______ Open output files ______
01363   ofstream fophi;
01364   openfstream(fophi,geocodeinput.fophi,generalinput.overwrit);
01365   bk_assert(fophi,geocodeinput.fophi,__FILE__,__LINE__); 
01366 
01367   ofstream folam;
01368   openfstream(folam,geocodeinput.folam,generalinput.overwrit);
01369   bk_assert(folam,geocodeinput.folam,__FILE__,__LINE__); 
01370 
01371   
01372 // ====== Process BUFFERS ======
01373   for (register int32 buffer=1; buffer<=FULLBUFFERS+EXTRABUFFER; buffer++)
01374     {
01375     // ______ Give progress ______
01376     PROGRESS << "GEOCODE: " 
01377          << 100*(buffer-1)/(FULLBUFFERS+EXTRABUFFER) << "%";
01378     PROGRESS.print();
01379 
01380     // ______ firstline of this buffer in master coordinate system ______
01381     real8 firstline = veryfirstline + (buffer-1) * bufferlines * multiL; 
01382 
01383     // ______ Set indices to be read from file / check if last buffer ______
01384     bufferwin.linelo = 1 + (buffer-1) * bufferlines;
01385     if (buffer == FULLBUFFERS+1)
01386       {
01387       bufferlines = RESTLINES;
01388       PHI.resize(bufferlines,mlpixels);
01389       LAMBDA.resize(bufferlines,mlpixels);
01390       HEIGHT.resize(bufferlines,mlpixels);
01391       }
01392     bufferwin.linehi = bufferwin.linelo + bufferlines - 1;
01393   
01394     // ______ Read in buffer of s2h ______
01395     // only real4 logical..., no member f. necessary.
01396     switch (heightinradarsystem.formatflag)
01397       {
01398       case FORMATR4:
01399         readfile (HEIGHT, heightinradarsystem.file, mllines, bufferwin,
01400         offsetbuffer);
01401         break;
01402   
01403       default:
01404         PRINT_ERROR("geocode format flag on file heights (s2h output) only real4 possible.");
01405         throw(unhandled_case_error);
01406       } // switch formatflag
01407   
01408 // ====== Compute xyz for all points on their height ======
01409     register int32 i,j;
01410     register real8 line = firstline - multiL;           // in master coordinate system
01411     cn pospoint;                                        // returned by lp2xyz
01412     input_ell ELLIPS;                                   // to correct height of ellips
01413     real8 r;                                            // for conversion xyz2philambda
01414     real8 nu;                                           // for conversion xyz2philambda
01415     real8 sin3;                                         // for conversion xyz2philambda
01416     real8 cos3;                                         // for conversion xyz2philambda
01417   
01418   
01419     for (i=0; i<HEIGHT.lines(); i++)
01420       {
01421       line += multiL;                                   // in master coordinate system
01422       register real8 pixel = firstpixel - multiP;       // in master coordinate system
01423   
01424       for (j=0; j<HEIGHT.pixels(); j++)
01425         {
01426         pixel += multiP;                                // in master coordinate system
01427   
01428 // ______ Check if conversion is necessary _______
01429         if (HEIGHT(i,j) == NaN)
01430           {
01431           PHI(i,j) = NaN;
01432           LAMBDA(i,j) = NaN;
01433           }
01434 
01435         else
01436           {
01437 // ______ Adjust height of ellips _______
01438           ELLIPS.a    = ellips.a + HEIGHT(i,j);
01439           ELLIPS.b    = ellips.b + HEIGHT(i,j);
01440           ELLIPS.ecc1st_sqr();
01441           ELLIPS.ecc2nd_sqr();
01442       
01443 // ______ Compute point on this ellips (Bowring) ______
01444 // THIS SHOULD BE DONE INTERNALLY TO SPEED UP ???
01445 // FIRST PROFILE HOW MUCH TIME THIS TAKES !!!
01446           lp2xyz(line,pixel,ELLIPS,master,masterorbit,
01447                  pospoint,MAXITER,CRITERPOS);
01448        
01449 // why real4??
01450           r           = sqrt(sqr(pospoint.x)+sqr(pospoint.y));
01451           nu          = atan2((pospoint.z*ELLIPS.a),(r*ELLIPS.b));
01452           sin3        = pow(sin(nu),3);
01453           cos3        = pow(cos(nu),3);
01454           PHI(i,j)    = rad2deg(atan2((pospoint.z+ELLIPS.e2b*ELLIPS.b*sin3),
01455                          (r-ELLIPS.e2*ELLIPS.a*cos3)));
01456           LAMBDA(i,j) = rad2deg(atan2(pospoint.y,pospoint.x));
01457           } //else
01458         } // loop j
01459       } // loop i
01460     fophi << PHI;
01461     folam << LAMBDA;
01462     } // loop buffers
01463   
01464 
01465   fophi.close();
01466   folam.close();
01467   
01468 // ====== Write to result files ======
01469   ofstream scratchlogfile("scratchloggeocode", ios::out | ios::trunc);
01470   bk_assert(scratchlogfile,"geocode: scratchloggeocode",__FILE__,__LINE__);
01471   scratchlogfile << "\n\n*******************************************************************"
01472                  << "\n*_Start_" << processcontrol[pr_i_geocoding]
01473                  << "\n*******************************************************************"
01474                  << "\nData_output_file_hei (slant2h): "
01475                  <<  heightinradarsystem.file
01476                  << "\nData_output_file_phi: \t\t"
01477                  <<  geocodeinput.fophi
01478                  << "\nData_output_file_lambda: \t"
01479                  <<  geocodeinput.folam
01480                  << endl;
01481   scratchlogfile.close();
01482 
01483 
01484   ofstream scratchresfile("scratchresgeocode", ios::out | ios::trunc);
01485   bk_assert(scratchresfile,"geocode: scratchresgeocode",__FILE__,__LINE__);
01486   scratchresfile << "\n\n*******************************************************************"
01487                  << "\n*_Start_" << processcontrol[pr_i_geocoding]
01488                  << "\n*******************************************************************"
01489                  << "\nData_output_file_hei (slant2h): "
01490                  <<  heightinradarsystem.file
01491                  << "\nData_output_file_phi:                 \t"
01492                  <<  geocodeinput.fophi
01493                  << "\nData_output_file_lamda:               \t"
01494                  <<  geocodeinput.folam
01495                  << "\nData_output_format:                   \t"
01496                  << "real4"
01497                  << "\nFirst_line (w.r.t. original_master):  \t"
01498                  <<  heightinradarsystem.win.linelo
01499                  << "\nLast_line (w.r.t. original_master):   \t"
01500                  <<  heightinradarsystem.win.linehi
01501                  << "\nFirst_pixel (w.r.t. original_master): \t"
01502                  <<  heightinradarsystem.win.pixlo
01503                  << "\nLast_pixel (w.r.t. original_master):  \t"
01504                  <<  heightinradarsystem.win.pixhi
01505                  << "\nMultilookfactor_azimuth_direction:    \t"
01506                  <<  heightinradarsystem.multilookL
01507                  << "\nMultilookfactor_range_direction:      \t"
01508                  <<  heightinradarsystem.multilookP
01509                  << "\n*******************************************************************"
01510                  << "\n* End_" << processcontrol[pr_i_geocoding] << "_NORMAL"
01511                  << "\n*******************************************************************\n";
01512   scratchresfile.close();
01513   PROGRESS.print("finished geocode.");
01514   } // END geocode
01515 

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