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

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

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