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

filtering.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/filtering.cc,v $  *
00032  * $Revision: 3.15 $                                            *
00033  * $Date: 2005/04/22 11:30:15 $                                 *
00034  * $Author: kampes $                                            *
00035  *                                                              *
00036  * Implementation of filtering routines:                        *
00037  * -rangefilter (adaptive, after resampling)                    *
00038  * -rangefilter (based on baseline, points on ellips)           *
00039  * -phasefilter: goldstein                                      *
00040  * -phasefilter: spatial convolution                            *
00041  * -phasefilter: spectral method                                *
00042  * -azimuthfilter: based on fDC polynomials                     *
00043  ****************************************************************/
00044 
00045 // #define _LIB_VERSION _IEEE_
00046 
00047 #include "matrixbk.hh"
00048 #include "constants.hh"                 // global constants
00049 #include "filtering.hh"                 // header file
00050 #include "utilities.hh"                 // myhamming
00051 #include "conversion.hh"                // cr4toci2
00052 #include "slcimage.hh"                  // my slc image class
00053 #include "productinfo.hh"               // my 'products' class
00054 #include "orbitbk.hh"                   // my orbit class
00055 #include "exceptions.hh"                 // my exceptions class
00056 
00057 #include <strstream>                    // for memory stream
00058 #include <iomanip>                      // setp
00059 
00060 
00061 
00062 
00063 
00064 
00065 /****************************************************************
00066  *    rangefilter                                               *
00067  * adaptive range filtering. 2b performed after resampling.     *
00068  * peak estimation of local fringe frequency in spectral        *
00069  * domain of oversampled([2] or 4) complex interferogram.       *
00070  * then master and slave spectra are filtered, optionally       *
00071  * de- and re-weight with hamming.                              *
00072  *                                                              *
00073  * Input:                                                       *
00074  *  - files for master/slave, input options                     *
00075  * Output:                                                      *
00076  *  - filtered files for master/slave                           *
00077  *                                                              *
00078  *    Bert Kampes, 29-Mar-2000                                  *
00079  ****************************************************************/
00080 void rangefilter(
00081         const input_gen       &generalinput,
00082         const slcimage        &master,
00083         const slcimage        &slave,
00084         const productinfo     &interferogram,
00085         const input_filtrange &filtrangeinput)
00086   {
00087   TRACE_FUNCTION("rangefilter (BK 29-Mar-2000)");
00088   INFO << "RANGEFILT: Master input file: " << master.file;
00089   INFO.print();
00090   INFO << "RANGEFILT: Slave input file:  " << slave.file;
00091   INFO.print();
00092 
00093   // ====== Handle input ======
00094   const real8 hamming      = filtrangeinput.hammingalpha;// alpha=1: const
00095   const int32 nlmean       = filtrangeinput.nlmean;      // odd number range lines to average
00096   const int32 fftlength    = filtrangeinput.fftlength;   // length of adaptive
00097   const int32 OVERLAP      = filtrangeinput.overlap;     // overlap input blocks
00098   const real8 SNRthreshold = filtrangeinput.SNRthreshold;// criterium to filter
00099   const bool doweightcorrel= filtrangeinput.doweightcorrel;
00100   const int32 oversamplefactor = filtrangeinput.oversample;
00101 
00102   // RSR,RBW in Hz
00103   const real8 RSR = 0.5*master.rsr2x;           // range sampling rate fr 18.96MHz
00104   const real8 RBW = master.rbw;                 // range band width fs 15.55MHz
00105 
00106 
00107 // ______ Number of lines on file ______
00108 // ______ Assume resampling has been done and interferowin contains master coord.
00109   const int32 Mfilelines        = master.currentwindow.lines();
00110   const int32 Ifilelines        = interferogram.win.lines();
00111   const int32 numpixels         = interferogram.win.pixels();
00112 
00113   // ______ Memory buffers ______
00114   const uint  BUFFERMEMSIZE  = generalinput.memory;
00115   const int32 bytesperline   = numpixels * sizeof(complr4);
00116   const real4 numbigmatrices = 4.2;                     // in/out Master/Slave +rest
00117   int32 bufferlines =                                   // numlines in buffer
00118     int32((BUFFERMEMSIZE/numbigmatrices)/bytesperline); 
00119 
00120 
00121 // ====== Open output files ======
00122   ofstream ofilemaster;
00123   openfstream(ofilemaster,filtrangeinput.fomaster,generalinput.overwrit);
00124   bk_assert(ofilemaster,filtrangeinput.fomaster,__FILE__,__LINE__);
00125 
00126   ofstream ofileslave;
00127   openfstream(ofileslave, filtrangeinput.foslave,generalinput.overwrit);
00128   bk_assert(ofileslave, filtrangeinput.foslave, __FILE__,__LINE__);
00129 
00130 
00131 // ====== Compute indices, counters, etc. ======
00132   const int32 numrangeblocks  = int32(numpixels/(fftlength-2*OVERLAP)); // approx
00133   real8 overallSNR         = 0.;                        // statistics
00134   real8 overallnotfiltered = 0.;                        // statistics
00135 
00136   // ______ Buffer counters ______
00137   int32 azimuthbuffer      = 1;                         // loop counter
00138   bool azimuthdone         = false;
00139   int32 startlinethisblock = interferogram.win.linelo;
00140   window filteredpart      (0, bufferlines-uint((nlmean+1)/2),
00141                               0, numpixels-1);                  // output to disk
00142 
00143   // ______ Give info ______
00144   INFO << "overlap master slave: ["
00145        << interferogram.win.linelo << ":" << interferogram.win.linehi
00146        << "; " << interferogram.win.pixlo << ":" << interferogram.win.pixhi
00147        << "] (master system)";
00148   INFO.print();
00149 
00150   // ====== Loop over azimuth buffers over whole width ======
00151   // ______ loop forever, break when finished ______
00152   for (azimuthbuffer=1; azimuthbuffer<1e5; ++azimuthbuffer)
00153     {
00154     PROGRESS << "azimuthbuffer (" << bufferlines << " lines): " 
00155          << azimuthbuffer << ".";
00156     PROGRESS.print();
00157 
00158     // ====== Read data master slave for this block ======
00159     int32 endlinethisblock = startlinethisblock+bufferlines-1;
00160     // ______ check slightly larger last block ______
00161     if (endlinethisblock>=interferogram.win.linehi-nlmean-int32(bufferlines/10))
00162       {
00163       azimuthdone         = true;                               // finish this one
00164       endlinethisblock    = interferogram.win.linehi;           // last line on file
00165       bufferlines         = endlinethisblock-startlinethisblock+1;
00166       filteredpart.linehi = bufferlines - 1;    // write all at end
00167       }
00168 
00169     // ______ Allocate matrices and read window from file ______
00170     const window winfile(startlinethisblock, endlinethisblock,
00171                          interferogram.win.pixlo,
00172                          interferogram.win.pixhi);
00173     // ______ Input buffer ______
00174     const matrix<complr4> MASTER_ORIG = master.readdata(winfile);
00175     const matrix<complr4> SLAVE_ORIG  = slave.readdata(winfile);
00176     // ______ Output buffer ______
00177     matrix<complr4> MASTER_FILT(MASTER_ORIG.lines(),MASTER_ORIG.pixels());
00178     matrix<complr4> SLAVE_FILT (SLAVE_ORIG.lines(), SLAVE_ORIG.pixels());
00179 
00180     // ====== Loop over range buffers dividing width ======
00181     // ______ Indices are in inbuffer local coordinates [0:numpixs] ______
00182     int32 startpixel    = 0;                    // valid for first range block
00183     int32 filteredpixlo = 0;                    // valid for first range block
00184     int32 partpixlo     = 0;                    // local system in PART
00185     int32 partpixhi     = fftlength-1-OVERLAP;  // local system in PART
00186     bool  rangedone     = false;                // not finished yet
00187     for (;;)    // forever
00188       {
00189       int32 endpixel      = startpixel+fftlength-1;     // BUFFER -> PART
00190       int32 filteredpixhi = endpixel - OVERLAP;         // PART   -> BUFFER
00191       if (endpixel>=MASTER_ORIG.pixels()-1)     // last block back-shifted
00192         {
00193         rangedone     = true;
00194         endpixel      = MASTER_ORIG.pixels()-1;
00195         startpixel    = endpixel - fftlength + 1;
00196         filteredpixhi = endpixel;               // write all for last block
00197         partpixhi     = fftlength-1;            // write all for last block
00198         partpixlo     = partpixhi - (filteredpixhi-filteredpixlo);
00199         }
00200 
00201       // ______ Some info ______
00202       TRACE << "range block [" << startlinethisblock << ":" << endlinethisblock
00203            << "; " << startpixel+interferogram.win.pixlo << ":"
00204                    << endpixel+interferogram.win.pixlo << "]";
00205       TRACE.print();
00206 
00207       // ______ Cut out master/slave for adaptive range filtering ______
00208       const window win(0,MASTER_ORIG.lines()-1,startpixel,endpixel);
00209       matrix<complr4> PARTMASTER(win,MASTER_ORIG);              // construct as part
00210       matrix<complr4> PARTSLAVE (win,SLAVE_ORIG);               // construct as part
00211 
00212       // ______ Do the actual range filtering ______
00213       // ______ filtered lines: index[0:L-1]: [(nlmean-1)/2:L-1-(nlmean-1)/2] ______
00214       real8 SNRmean, percentnotfiltered;
00215       rfilterblock(PARTMASTER,PARTSLAVE,                        // returned
00216                    nlmean,SNRthreshold,RSR,RBW,                 // parameters
00217                    hamming,oversamplefactor,doweightcorrel,     // control algorithm
00218                    SNRmean, percentnotfiltered);                // statistics returned
00219       overallSNR         += SNRmean;
00220       overallnotfiltered += percentnotfiltered;
00221 
00222       // ______ Put filtered part in outputbuffer ______
00223       const window winpart(win.linelo,win.linehi, partpixlo,partpixhi);
00224       const window wintoset(win.linelo,win.linehi, filteredpixlo,filteredpixhi);
00225       MASTER_FILT.setdata(wintoset,PARTMASTER,winpart);
00226       SLAVE_FILT.setdata (wintoset,PARTSLAVE, winpart);
00227 
00228       // ______ Update counters and check if it is already time to go ______
00229       if (rangedone==true) break;
00230       partpixlo     = OVERLAP;                  // valid for all next blocks
00231       startpixel    = endpixel - 2*OVERLAP + 1; // get PART from inputbuffer
00232       filteredpixlo = startpixel + OVERLAP;     // valid for all blocks except first
00233       } // loop over range blocks
00234 
00235     // ====== Write filtered parts master,slave to output files ======
00236     switch (filtrangeinput.oformatflag)
00237       {
00238       case FORMATCR4:
00239         {
00240         writefile(ofilemaster,MASTER_FILT,filteredpart);
00241         writefile(ofileslave, SLAVE_FILT, filteredpart);
00242         break;
00243         }
00244       // ______ Convert first to ci2 before writing to file ______
00245       case FORMATCI2:
00246         {
00247         register int32 ii,jj;
00248         matrix <compli16> TMP(MASTER_FILT.lines(),MASTER_FILT.pixels());
00249         for (ii=filteredpart.linelo; ii<=filteredpart.linehi; ++ii)
00250           for (jj=filteredpart.pixlo; jj<=filteredpart.pixhi; ++jj)
00251             TMP(ii,jj) = cr4toci2(MASTER_FILT(ii,jj));
00252         writefile(ofilemaster,TMP,filteredpart);
00253         for (ii=filteredpart.linelo; ii<=filteredpart.linehi; ++ii)
00254           for (jj=filteredpart.pixlo; jj<=filteredpart.pixhi; ++jj)
00255             TMP(ii,jj) = cr4toci2(SLAVE_FILT(ii,jj));
00256         writefile(ofileslave, TMP, filteredpart);
00257         break;
00258         }
00259       default:
00260         {
00261         PRINT_ERROR("Totally impossible, checked input.")
00262         throw(unhandled_case_error);
00263         }
00264       }
00265     // ______ Update counters, check finished ______
00266     if (azimuthdone==true) break;
00267     filteredpart.linelo = uint((nlmean-1.)/2.); // same for all buffers except first
00268     startlinethisblock  = endlinethisblock - nlmean + 2; 
00269     } // loop over azimuth buffers
00270 
00271   // _______ Some stats, better dump SUN rasterfile all shifts ______
00272   overallSNR         /= (azimuthbuffer*numrangeblocks);
00273   overallnotfiltered /= (azimuthbuffer*numrangeblocks);
00274   INFO << "Mean SNR for all blocks (approx): " << overallSNR;
00275   INFO.print();
00276   INFO << "Filtered (approx): " << setprecision(3) 
00277        << 100.00-overallnotfiltered << "%";
00278   if (overallnotfiltered<60.) 
00279     {
00280     INFO.print();
00281     }
00282   else
00283     {
00284     WARNING.print(INFO.get_str());
00285     INFO.rewind();
00286     }
00287 
00288 
00289 // ====== Write results to file ======
00290   ofstream scratchlogfile("scratchlogfiltrange", ios::out | ios::trunc);
00291   bk_assert(scratchlogfile,"filtrange: scratchlogfiltrange",__FILE__,__LINE__);
00292   scratchlogfile
00293     << "\n\n*******************************************************************"
00294     << "\nRange filtering adaptive for master and slave..."
00295     << "\nInput file master (format): \t\t\t"
00296     <<  master.file << " " << master.formatflag
00297     << "\nInput file slave (format): \t\t\t"
00298     <<  slave.file << " " << slave.formatflag
00299     << "\nOutput file filtered master (format): \t\t\t"
00300     <<  filtrangeinput.fomaster << " ";
00301     if (filtrangeinput.oformatflag==FORMATCR4)
00302       scratchlogfile << "complex r4";
00303     if (filtrangeinput.oformatflag==FORMATCI2)
00304       scratchlogfile << "complex short int";
00305   scratchlogfile
00306     << "\nOutput file filtered slave (format): \t\t\t"
00307     <<  filtrangeinput.foslave << " " << FORMATCR4
00308     << "\n..."
00309     << "\n*******************************************************************\n";
00310   scratchlogfile.close();
00311 
00312   for (int32 fileID=1; fileID<=2; ++fileID)
00313     {
00314     char oresfile[EIGHTY];
00315     char odatafile[EIGHTY];
00316     char odataformat[EIGHTY];
00317     char processcf[EIGHTY];
00318     if (filtrangeinput.oformatflag==FORMATCR4)
00319       strcpy(odataformat,"complex_real4");
00320     else if (filtrangeinput.oformatflag==FORMATCI2)
00321       strcpy(odataformat,"complex_short");
00322     else
00323       {
00324       PRINT_ERROR("case problem.")
00325       throw(unhandled_case_error);
00326       }
00327     switch (fileID)
00328       {
00329       case 1:                                                   // master
00330         strcpy(oresfile,"scratchresMfiltrange");
00331         strcpy(odatafile,filtrangeinput.fomaster);
00332         strcpy(processcf,processcontrol[pr_m_filtrange]);       // control flag
00333         break;
00334       case 2:                                                   // slave
00335         strcpy(oresfile,"scratchresSfiltrange");
00336         strcpy(odatafile,filtrangeinput.foslave);
00337         strcpy(processcf,processcontrol[pr_s_filtrange]);       // control flag
00338         break;
00339       default:
00340         {
00341         PRINT_ERROR("panic: ID!={1,2}")
00342         throw(unhandled_case_error);
00343         }
00344       }
00345 
00346     // ______ updateproductinfo greps info from this file ______
00347     ofstream scratchresfile(oresfile, ios::out | ios::trunc);
00348     bk_assert(scratchresfile,oresfile,__FILE__,__LINE__);
00349     scratchresfile
00350       << "\n\n*******************************************************************"
00351       << "\n*_Start_" << processcf
00352       << "\n*******************************************************************"
00353       << "\nMethod_rangefilt:                     \t"
00354       << "adaptive"
00355       << "\nData_output_file:                     \t"
00356       <<  odatafile
00357       << "\nData_output_format:                   \t"
00358       <<  odataformat
00359       << "\nFirst_line (w.r.t. original_master):  \t"
00360       <<  interferogram.win.linelo
00361       << "\nLast_line (w.r.t. original_master):   \t"
00362       <<  interferogram.win.linehi
00363       << "\nFirst_pixel (w.r.t. original_master): \t"
00364       <<  interferogram.win.pixlo
00365       << "\nLast_pixel (w.r.t. original_master):  \t"
00366       <<  interferogram.win.pixhi
00367       << "\n*******************************************************************"
00368       << "\n* End_" << processcf << "_NORMAL"
00369       << "\n*******************************************************************\n";
00370     scratchresfile.close();
00371     } // for master and slave
00372   // ______Tidy up______
00373   } // END rangefilter
00374 
00375 
00376 
00377 /****************************************************************
00378  *    rfilterblock                                              *
00379  * Computes powerspectrum of complex interferogram.             *
00380  * (product oversampled master. conj oversampled slave in range)*
00381  * A peak in this spectrum corresponds to frequency shift.      *
00382  * The master and slave are LPF filtered for this shift.        *
00383  *                                                              *
00384  * Optionally the oversampling can be turned off, since no use  *
00385  * if only small baseline, and flat terrain.                    *
00386  * The powerspectrum can be weighted to give more influence to  *
00387  * higher frequencies (conv. of 2 blocks, should be hamming).   *
00388  * The peak is detected by taking a mean of nlmean lines (odd). *
00389  * Filtering is applied if the SNR (N*power peak / power rest)  *
00390  * is above a user supplied threshold.                          *
00391  * At LPF filtering of the master/slave a hamming window may    *
00392  * be applied first to deweight, then to re-weight the spectrum *
00393  *                                                              *
00394  * Should filter based on zero terrain slope if below SNR,      *
00395  * but this requried knowledge of orbits, pixel,line coordinate *
00396  #%// BK 13-Nov-2000                                            *
00397  *                                                              *
00398  * Input:                                                       *
00399  *  - MASTER: block of master, that will be filtered            *
00400  *  - SLAVE:  block of slave, that will be filtered             *
00401  * Output:                                                      *
00402  *  - MASTER (SLAVE): filtered from indeces[0:numl-1]           *
00403  *    (nlmean-1)/2 to numlines-(nlmean-1)/2-1                   *
00404  *                                                              *
00405  *    Bert Kampes, 29-Mar-2000                                  *
00406  ****************************************************************/
00407 void rfilterblock(
00408         matrix<complr4> &master,                // updated
00409         matrix<complr4> &slave,                 // updated
00410         int32 nlmean,                           // number of lines to take mean over
00411         real8 SNRthreshold,                     // 
00412         real8 RSR,                              // in MHz
00413         real8 RBW,                              // in MHz
00414         real8 alphahamming,                     // parameter hamming filtering [0,1]
00415         int32  osfactor,                        // oversampling factor interf. gen.
00416         bool  doweightcorrel,                   // correct correl for #elem.
00417         real8 &meanSNR,                         // returned
00418         real8 &percentnotfiltered)              // returned
00419   {
00420   const int32 numlines    = master.lines();
00421   const int32 numpixs     = master.pixels();
00422   const int32 outputlines = numlines - nlmean + 1;
00423   const int32 firstline   = int32((nlmean-1)/2);        // indices in matrix system
00424   const int32 lastline    = firstline + outputlines - 1;
00425   const bool dohamming    = (alphahamming<0.9999) ? true : false;
00426   // use oversampling before int. gen.
00427   const bool dooversample = (osfactor!=1) ? true : false;
00428   int32 notfiltered=0;                                  // counter
00429 
00430 
00431 #ifdef __DEBUG
00432   TRACE_FUNCTION("rfilterblock (BK 29-Mar-2000)")
00433   if (!isodd(nlmean))          
00434     {
00435     PRINT_ERROR("nlmean has to be odd.")
00436     throw(argument_error);
00437     }
00438   if (!ispower2(numpixs))
00439     {
00440     PRINT_ERROR("numpixels (FFT) has to be power of 2.")
00441     throw(argument_error);
00442     }
00443   if (!ispower2(osfactor))
00444     {
00445     PRINT_ERROR("oversample factor (FFT) has to be power of 2.")
00446     throw(argument_error);
00447     }
00448   if (slave.lines()!=numlines)
00449     {
00450     PRINT_ERROR("slave not same size as master.")
00451     throw(argument_error);
00452     }
00453   if (slave.pixels()!=numpixs)
00454     {
00455     PRINT_ERROR("slave not same size as master.")
00456     throw(argument_error);
00457     }
00458 #endif
00459   if (outputlines < 1)
00460     {
00461     WARNING.print("no outputlines, continuing."); 
00462     return;
00463     }
00464 
00465 
00466   // ______ Shift parameters ______
00467   register int32 i,j;
00468   const real4 deltaf =  RSR/real4(numpixs);
00469   const real4 fr     = -RSR/2.;
00470   matrix<real4> freqaxis(1,numpixs);
00471   for (i=0; i<numpixs; ++i)
00472     freqaxis(0,i) = fr+real4(i)*deltaf;
00473  
00474   matrix<real4> inversehamming;
00475   if (dohamming)
00476     {
00477     inversehamming = myhamming(freqaxis,RBW,RSR,alphahamming);
00478     for (i=0; i<numpixs; ++i)
00479       if (inversehamming(0,i)!=0.)
00480         inversehamming(0,i)= 1./inversehamming(0,i);
00481     }
00482 
00483 // ______ Compute complex interferogram -> power ______
00484   matrix<complr4> cint;
00485   if (dooversample)
00486     cint = dotmult(oversample(master,1,osfactor),oversample(conj(slave),1,osfactor));
00487   else
00488     cint = dotmult(master,conj(slave));
00489   const int32 fftlength = cint.pixels();
00490 
00491   DEBUG.print("is real4 accurate enough?");// seems so
00492   fft(cint,2);                                  // cint=fft over rows
00493   matrix<real4> power   = intensity(cint);      // power=cint.*conj(cint);
00494 
00495 
00496   // ______ Use weighted correlation due to bias in normal definition ______
00497   // ______ Actually better deweight with autoconvoluted hamming.
00498   // ______ No use a triangle for #points used for correlation estimation
00499   // ______ not in combination with dooversample...
00500   if (doweightcorrel)
00501     {
00502     INFO.print("De-weighting power spectrum of convolution.");
00503     //matrix<real4> weighting(1,fftlength);
00504     //for (i=0; i<fftlength; ++i)
00505     //  weighting(0,i) = abs(i)...;
00506     //power *= weightingfunction...== inv.triangle
00507 
00508     // weigth = numpoints in spectral convolution for fft squared for power...
00509     const int32 indexnopeak = int32((1.-(RBW/RSR))*real4(numpixs));
00510     for (j=0; j<fftlength; ++j)
00511       {
00512       const int32 npnts = abs(numpixs-j);
00513       const real4 weight =                      // oversample: numpixs==fftlength/2
00514         (npnts<indexnopeak) ? sqr(numpixs) : sqr(npnts);        // ==zero
00515       for (i=0; i<numlines; ++i)
00516         {
00517         power(i,j) /= weight;
00518         }
00519       }
00520     }
00521 
00522 
00523   // ______ Average power to reduce noise ______
00524   fft(master,2);                                // master=fft over rows
00525   fft(slave,2);                                 // slave=fft over rows
00526   TRACE.print("Took FFT over rows of master, slave.");
00527   matrix<real4> nlmeanpower = sum(power(0,nlmean-1, 0,fftlength-1),1);
00528 
00529   uint shift      = 0;                          // returned by max
00530   uint dummy      = 0;                          // returned by max
00531   meanSNR         = 0.;
00532   real8 meanSHIFT = 0.;
00533 
00534   // ______ Start actual filtering ______
00535   // ______ oline is index in matrix system ______
00536   for (register int32 oline=firstline; oline<=lastline; ++oline)
00537     {
00538     matrix<real4> totalp   = sum(nlmeanpower,2);        // 1x1 matrix ...
00539     const real4 totalpower = totalp(0,0);
00540     const real4 maxvalue   = max(nlmeanpower,dummy,shift);      // shift returned
00541     uint lastshift         = shift;     // use this if current shift not ok.
00542     const real8 SNR        = fftlength*(maxvalue/(totalpower-maxvalue));
00543     meanSNR               += SNR;
00544 
00545     // ______ Check for negative shift ______
00546     bool negshift = false;
00547     if (shift > fftlength/2)
00548       {
00549       // ERR: uint! shift    = uint(abs(real8(shift-fftlength)));
00550       shift     = uint(fftlength)-shift;
00551       lastshift = shift;// use this if current shift not OK.
00552       negshift  = true;
00553       }
00554 
00555     // ______ Do actual filtering ______
00556     if (SNR<SNRthreshold)
00557       {
00558       notfiltered++;                                    // update counter
00559       shift = lastshift;// but use this anyway.
00560       //WARNING.print("using last shift for filter");
00561       }
00562     meanSHIFT += shift;
00563     matrix<real4> filter;
00564     if (dohamming)
00565         {
00566         // ______ Newhamming is scaled and centered around new mean ______
00567         filter  = myhamming(freqaxis-real4(.5*shift*deltaf),
00568                             RBW-real8(shift*deltaf),
00569                             RSR,alphahamming);          // fftshifted
00570         filter *= inversehamming;
00571         }
00572       else                                              // no weighting of spectra
00573         {
00574         filter = myrect((freqaxis-real4(.5*shift*deltaf)) /
00575                           (real4(RBW-shift*deltaf)));   // fftshifted
00576         }
00577       // ______ Use freq. as returned by fft ______
00578       // ______ Note that filter_s = fliplr(filter_m) ______
00579       // ______ and that this is also valid after ifftshift ______
00580       ifftshift(filter);                                // fftsh works on data!
00581 
00582       // ====== Actual spectral filtering ======
00583       // ______ Decide which side to filter, may be dependent on ______
00584       // ______ definition of FFT, this is ok for VECLIB ______
00585       // ______ if you have trouble with this step, either check your FFT
00586       // ______ (or use uinternal one, or add a card for changing false to true below
00587       if (negshift==false)
00588         {
00589         dotmult(master[oline],filter,1); 
00590         filter.fliplr();                                // works on data!
00591         dotmult(slave[oline],filter,1);
00592         }
00593       else
00594         {
00595         dotmult(slave[oline],filter,1);
00596         filter.fliplr();                                // works on data!
00597         dotmult(master[oline],filter,1);
00598         }
00599 // following is removed, we now always filter with last know spectral offset
00600 // Bert Kampes, 13-Sep-2004
00601 //      } // SNR>threshold
00602 //    else
00603 //      {
00604 //      notfiltered++;                                  // update counter
00605 //      }
00606 
00607     // ______ Update 'walking' mean ______
00608     if (oline!=lastline)                        // then breaks
00609       {
00610       matrix<real4> line1 = power.getrow(oline-firstline);
00611       matrix<real4> lineN = power.getrow(oline-firstline+nlmean);
00612       nlmeanpower += (lineN-line1);
00613       }
00614     } // loop over outputlines
00615 
00616   // ______ IFFT of spectrally filtered data, and return these ______
00617   ifft(master,2);                               // master=ifft over rows
00618   ifft(slave,2);                                // slave=ifft over rows
00619 
00620   // ______ Return these to main ______
00621   meanSHIFT          /= (outputlines-notfiltered);
00622   meanSNR            /= outputlines;
00623   percentnotfiltered  = 100. * (real4(notfiltered)/real4(outputlines));
00624 
00625 
00626   // ______ Some info for this block ______
00627   const real8 meanfrfreq = meanSHIFT*deltaf;    // Hz?
00628   DEBUG << "mean SHIFT for block"
00629        << ": "  << meanSHIFT
00630        << " = " << meanfrfreq/1e6 << " MHz (fringe freq.).";
00631   DEBUG.print();
00632   DEBUG << "mean SNR for block"
00633        << ": " << meanSNR;
00634   DEBUG.print();
00635   DEBUG << "filtered for block"
00636        << ": " << setprecision(3) << 100.00-percentnotfiltered << "%" << ends;
00637   DEBUG.print();
00638   if (percentnotfiltered>60.0) 
00639     {
00640     WARNING.print(DEBUG.get_str());
00641     DEBUG.reset();
00642     }
00643   } // END rfilterblock
00644 
00645 
00646 
00647 /****************************************************************
00648  *    phasefilter                                               *
00649  * goldsteins method, see routine goldstein and smooth.         *
00650  * After Goldstein and Werner, Radar interferogram filtering    *
00651  * for geophysical applications. GRL 25-21 pp 4035-4038, 1998.  *
00652  * and: ESA Florence 1997, vol2, pp969-972, Goldstein & Werner  *
00653  * "Radar ice motion interferometry".                           *
00654  #%// BK 24-Oct-2000                                            *
00655  ****************************************************************/
00656 void phasefilter(
00657         const input_gen       &generalinput,
00658         const productinfo     &interferogram,
00659         const input_filtphase &filtphaseinput)
00660   {
00661   TRACE_FUNCTION("phasefilter (BK 24-Oct-2000)")
00662 
00663   // ====== Handle input ======
00664   bool doexternalfile = false;
00665   if (specified(filtphaseinput.fifiltphase))
00666     doexternalfile = true;
00667   char infile[EIGHTY];                          // file 2b filtered
00668   strcpy(infile,interferogram.file);
00669   int32 numlinesinput  = int32(interferogram.win.lines()/interferogram.multilookL);
00670   int32 numpixelsinput = int32(interferogram.win.pixels()/interferogram.multilookP);
00671   if (doexternalfile)                           // overwrite default file
00672     {
00673     numlinesinput  = filtphaseinput.finumlines;
00674     strcpy(infile,filtphaseinput.fifiltphase);
00675     ifstream ifile;
00676     openfstream(ifile,infile);
00677     bk_assert(ifile,infile,__FILE__,__LINE__);
00678     ifile.seekg(0,ios::end);                    // pointer to end...
00679     const int32 sizefile     = ifile.tellg();   // opened ate
00680     ifile.close();
00681     numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;// floor
00682     if (numlinesinput*numpixelsinput*sizeof(complr4) != sizefile)
00683       WARNING.print("Format infile not CR4, or numlines not correct.");
00684     INFO.print("Using input file for phase filtering, not default.");
00685     }
00686   INFO << "phasefilter: inputfile: (#l,#p): " << infile << " ("
00687        << numlinesinput << "," << numpixelsinput << ")";
00688   INFO.print();
00689 
00690   // ______ Set variables ______
00691   const real8 ALPHA    = filtphaseinput.alpha;          // 0<a<1;
00692   const int32 SIZE     = filtphaseinput.blocksize;      // power of 2
00693   const int32 OVERLAP  = filtphaseinput.overlap;        // half of the overlap
00694 
00695   // ______ Open output file ______
00696   ofstream ofile;
00697   openfstream(ofile,filtphaseinput.fofiltphase,generalinput.overwrit);
00698   bk_assert(ofile,filtphaseinput.fofiltphase,__FILE__,__LINE__);
00699 
00700 
00701   // ______ initialize indices ______
00702   const int32 numout   = SIZE-(2*OVERLAP);      // number of output lines
00703   int32 cintlinelo     = 0;                     // index in CINT to get 1st buffer
00704   int32 cintlinehi     = SIZE-1;                // index in CINT to get 1st buffer
00705   int32 filteredlinelo = 0;                     // index in FILTERED to write (1st)
00706   int32 filteredlinehi = SIZE-OVERLAP-1;        // index in FILTERED to write (1st)
00707   bool  lastbuffer     = false;                 // not yet, just starting...
00708 
00709   if (!doexternalfile)
00710     if (interferogram.formatflag!=FORMATCR4)
00711       {
00712       PRINT_ERROR("Sorry, for phasefilter, format of interferogram must be CR4")
00713       throw(argument_error);
00714       }
00715   matrix<complr4> CINT(SIZE,numpixelsinput);    // type must be cr4!
00716 
00717   int32 lineswritten     = OVERLAP;             // first buffer extra OVERLAP lines
00718   const int32 numbuffers = numlinesinput/numout;        // approx?
00719   int32 tenpercent = int32((numbuffers/10)+.5); // round
00720   if (tenpercent==0) tenpercent = 1000;
00721   PROGRESS.print("FILTPHASE:  0%");
00722   int32 percent          = 10;
00723 
00724   // ====== Start filter loop buffers of BLOCKSIZE ======
00725   for (int32 buffer=1; buffer<1e6; ++buffer)
00726     {
00727     // ______ Check if this will be the last buffer ______
00728     if (cintlinehi >= numlinesinput-1)          // -1 since in matrix system
00729       {
00730       lastbuffer     = true;
00731       cintlinehi     = numlinesinput-1;
00732       cintlinelo     = cintlinehi-SIZE+1;       // make sure SIZE lines are read
00733       filteredlinehi = SIZE-1;                  // write upto lastline for last buffer
00734       const int32 lines2bwritten = numlinesinput-lineswritten;
00735       filteredlinelo = filteredlinehi - lines2bwritten + 1;
00736       if (lines2bwritten<1) 
00737         WARNING.print("PANIC: this will crash, lines2bwritten<1?");
00738       }
00739 
00740     // ______ Read in buffers of BLOCKSIZE lines complex interferogram ______
00741     const window windummy (0,0,0,0);            // no offset in readfile
00742     const window wincint  (cintlinelo,cintlinehi,0,numpixelsinput-1);
00743     readfile(CINT,infile,numlinesinput,wincint,windummy);
00744 
00745     // ______ Filter data in buffer of BLOCKSIZE lines ______
00746     // ______ output has same size, but write only part to disk ______
00747     matrix<complr4> FILTERED = goldstein(CINT,ALPHA,OVERLAP,filtphaseinput.kernel);
00748 
00749     // ______ Write filtered data ______
00750     const window winfiltered(filteredlinelo,filteredlinehi,0,numpixelsinput-1);
00751     writefile(ofile,FILTERED,winfiltered);
00752 
00753     // ______ Check if all done ______
00754     if (lastbuffer)
00755       break;
00756 
00757     // ______ Update indexes in matrices, will be corrected for last block ______
00758     lineswritten  += numout;            // first block overlap added
00759     cintlinelo    += numout;            // next buffer
00760     cintlinehi    += numout;            // next buffer
00761     filteredlinelo = OVERLAP;           // index in buffer,
00762                                         // valid for all middle buffers, but not last
00763     // ______ Give progress message ______
00764     if (buffer%tenpercent==0)
00765       {
00766       PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
00767       PROGRESS.print();
00768       percent += 10;
00769       }
00770     } // end loop buffers
00771   ofile.close();
00772 
00773   // ______ exit if only external file was desired ______
00774   if (doexternalfile)
00775     {
00776     cerr  << "\n Finished external file phasefilter, Exiting\n";
00777     PROGRESS.print("\n Finished external file phasefilter, Exiting\n");
00778     exit(0);
00779     }
00780 
00781   // ====== Write results to file ======
00782   ofstream scratchlogfile("scratchlogfiltphase", ios::out | ios::trunc);
00783   bk_assert(scratchlogfile,"filtphase: scratchlogfiltphase",__FILE__,__LINE__);
00784   scratchlogfile
00785     << "\n\n*******************************************************************"
00786     << "\nPhase filter complex interferogram: \t"
00787     <<  infile
00788     << "\nOutput file filtered master (format): \t\t\t"
00789     <<  filtphaseinput.fofiltphase << " "
00790     << "\n..."
00791     << "\n*******************************************************************\n";
00792   scratchlogfile.close();
00793 
00794   // ______ updateproductinfo greps info from this file ______
00795   ofstream scratchresfile("scratchresfiltphase", ios::out | ios::trunc);
00796   bk_assert(scratchresfile,"scratchresfiltphase",__FILE__,__LINE__);
00797   scratchresfile
00798     << "\n\n*******************************************************************"
00799     << "\n*_Start_" << processcontrol[pr_i_filtphase]
00800     << "\n*******************************************************************"
00801     << "\nMethod_phasefilt: goldstein: size, alpha, overlap: \t"
00802     <<  SIZE << " " << ALPHA << " " << OVERLAP
00803     << "\n1D Smoothing kernel for |spectrum|:   \t";
00804   for (int32 ii=0; ii<filtphaseinput.kernel.pixels(); ++ii)
00805     scratchresfile << " " << filtphaseinput.kernel(0,ii);
00806   scratchresfile
00807     << "\nInput_file:                           \t"
00808     <<  infile
00809     << "\nData_output_file:                     \t"
00810     <<  filtphaseinput.fofiltphase
00811     << "\nData_output_format:                   \t"
00812     <<  "complex_real4"
00813     << "\nFirst_line (w.r.t. original_master):  \t"
00814     <<  interferogram.win.linelo
00815     << "\nLast_line (w.r.t. original_master):   \t"
00816     <<  interferogram.win.linehi
00817     << "\nFirst_pixel (w.r.t. original_master): \t"
00818     <<  interferogram.win.pixlo
00819     << "\nLast_pixel (w.r.t. original_master):  \t"
00820     <<  interferogram.win.pixhi
00821     << "\nMultilookfactor_azimuth_direction:    \t"
00822     <<  interferogram.multilookL
00823     << "\nMultilookfactor_range_direction:      \t"
00824     <<  interferogram.multilookP
00825     << "\nNumber of lines (multilooked):        \t"
00826     <<  numlinesinput
00827     << "\nNumber of pixels (multilooked):       \t"
00828     <<  numpixelsinput
00829     << "\n*******************************************************************"
00830     << "\n* End_" << processcontrol[pr_i_filtphase] << "_NORMAL"
00831     << "\n*******************************************************************\n";
00832   scratchresfile.close();
00833   // ====== Tidy up ======
00834   } // END phasefilter
00835 
00836 
00837 
00838 /****************************************************************
00839  *    phasefilter goldstein                                     *
00840  * Input is matrix of SIZE (e.g. 32) lines, and N range pixels. *
00841  * Filtered OUTPUT is same size as input block.                 *
00842  * Because of overlap, only write to disk in calling routine    *
00843  * part (in matrix coord.) [OVERLAP:SIZE-OVERLAP-1]             *
00844  *                                                              *
00845  * Smoothing of the amplitude of the spectrum is performed by   *
00846  * spatial convolution with a block kernel of size 2*SMOOTH+1.  *
00847  * (Which is done by FFT's). e.g. a spatial moving average with *
00848  * kernel (1d) k=[1 1 1 1 1]/5; kernel2d = transpose(k)*k.      *
00849  * Blocks in range direction.                                   * 
00850  *                                                              *
00851  * After Goldstein and Werner, Radar interferogram filtering    *
00852  * for geophysical applications. GRL 25-21 pp 4035-4038, 1998.  *
00853  * and: ESA Florence 1997, vol2, pp969-972, Goldstein & Werner  *
00854  * "Radar ice motion interferometry".                           *
00855  #%// BK 25-Oct-2000                                            *
00856  ****************************************************************/
00857 matrix<complr4> goldstein(
00858         const matrix<complr4> &CINT,
00859         const real4            ALPHA,
00860         const int32            OVERLAP,
00861         const matrix<real4>   &smoothkernel)    // lying down
00862   {
00863   TRACE_FUNCTION("ROUTINE: goldstein (BK 25-Oct-2000)")
00864   //#define CHECKINDICESCALLING
00865   #ifdef CHECKINDICESCALLING
00866     return CINT;
00867   #else
00868   // ______ Allocate output matrix ______
00869   const int32 SIZE = CINT.lines();
00870   const int32 NPIX = CINT.pixels();
00871   matrix<complr4> FILTERED(SIZE,NPIX);          // output
00872 
00873   // ______ Get block from buffer ______
00874   const int32 numout  = SIZE-(2*OVERLAP);       // number of output pixels
00875   int32 cintpixlo     = 0;                      // index in CINT to get 1st block
00876   int32 cintpixhi     = SIZE-1;                 // index in CINT to get 1st block
00877   int32 outblockpixlo = 0;                      // index in BLOCK (only 1st block)
00878   int32 outblockpixhi = SIZE-1-OVERLAP;         // index in BLOCK (except last block)
00879   int32 outpixlo      = outblockpixlo;          // index in FILTERED (1st block)
00880   int32 outpixhi      = outblockpixhi;          // index in FILTERED
00881   bool  lastblockdone = false;                  // only just started...
00882   // note that int32() floors division 
00883   const int32 SMOOTH  = int32(smoothkernel.pixels())/2; // half block size, odd kernel
00884   const bool dosmooth = (SMOOTH==0) ? false : true;
00885   DEBUG << "SMOOTH: " << SMOOTH;// problem with uint<0 index in smoothkernel
00886   DEBUG.print();
00887 
00888   // ______ use FFT's for convolution with smoothkernel ______
00889   // ______ this could also be done static, or in the calling routine ______
00890   // ______ KERNEL2D is FFT2 of even kernel (no imag part after fft!) ______
00891   matrix<complr4> KERNEL2D;
00892   if (dosmooth==true)
00893     {
00894     matrix<complr4> kernel(1,SIZE);             // init to zeros
00895     for (register int32 ii=-SMOOTH; ii<=SMOOTH; ++ii)// 1d kernel function of block
00896       {
00897       //kernel(0,(ii+SIZE)%SIZE) = smoothkernel(0,ii-SMOOTH);
00898       // BK 07-Apr-2003
00899       // Cygwin fails on index very large number, i.e., uint problem somewhere
00900       // e.g.: [30,31,0,1,2] <--> [0,1,2,3,4]
00901       int32 tmp1 = (ii+SIZE)%SIZE;
00902       int32 tmp2 = ii+SMOOTH;// used to be ii-SMOOTH: wrong
00903       DEBUG << "tmp1: " << tmp1 << "; tmp2: " << tmp2;
00904       DEBUG.print();
00905       kernel(0,tmp1) = complr4(smoothkernel(0,tmp2),real4(0.0));
00906       }
00907     KERNEL2D = matTxmat(kernel,kernel);
00908     fft2d(KERNEL2D);                            // should be real sinc
00909     }
00910   DEBUG.print("kernel created for smoothing spectrum");
00911 
00912   // ====== Loop forever, stop after lastblockdone ======
00913   for (;;)      //forever
00914     {
00915     if (cintpixhi>=NPIX-1)                      // check if we are doing the last block
00916       {
00917       lastblockdone = true;
00918       cintpixhi     = NPIX-1;                   // prevent reading after file
00919       cintpixlo     = cintpixhi-SIZE+1;         // but make sure SIZE pixels are read
00920       outpixhi      = cintpixhi;                // index in FILTERED 2b written
00921       outblockpixhi = SIZE-1;                   // write all to the end
00922       outblockpixlo = outblockpixhi - (outpixhi-outpixlo+1) + 1;
00923       }
00924     const window wincint     (0,SIZE-1,cintpixlo,cintpixhi);
00925     const window winblock    (0,SIZE-1,outblockpixlo,outblockpixhi);
00926     const window winfiltered (0,SIZE-1,outpixlo,outpixhi);
00927 
00928     // ______ Construct BLOCK as part of CINT ______
00929     matrix<complr4> BLOCK(wincint,CINT);
00930 
00931     //#define CHECKINDEXONLY
00932     #ifndef CHECKINDEXONLY
00933     // ______ Get spectrum/amplitude/smooth/filter ______
00934     fft2d(BLOCK);
00935     matrix<real4> AMPLITUDE = magnitude(BLOCK);
00936 
00937     // ______ use FFT's for convolution with rect ______
00938     if (dosmooth==true)
00939       AMPLITUDE = smooth(AMPLITUDE,KERNEL2D);
00940 
00941     //dumpasc("A",AMPLITUDE);
00942     //AMPLITUDE = smooth(AMPLITUDE,SMOOTH);
00943     //dumpasc("As",AMPLITUDE); exit(1);
00944     const real4 maxamplitude = max(AMPLITUDE);
00945     if (maxamplitude>1e-20) //?
00946       {
00947       AMPLITUDE /= maxamplitude;
00948       AMPLITUDE.mypow(ALPHA);
00949       BLOCK     *= AMPLITUDE;           // weight spectrum
00950       }
00951     else
00952       {
00953       WARNING.print("no filtering, maxamplitude<1e-20, zeros in this block?");
00954       }
00955     ifft2d(BLOCK);
00956     #endif // check index blocks
00957     // ______ Set correct part that is filtered in output matrix ______
00958     FILTERED.setdata(winfiltered,BLOCK,winblock);
00959 
00960     // ______ Exit if finished ______
00961     if (lastblockdone)
00962       return FILTERED;                  // return
00963 
00964     // ______ Update indexes in matrices, will be corrected for last block ______
00965     cintpixlo    += numout;             // next block
00966     cintpixhi    += numout;             // next block
00967     outblockpixlo = OVERLAP;            // index in block, valid for all middle blocks
00968     outpixlo      = outpixhi+1;         // index in FILTERED, next range line
00969     outpixhi      = outpixlo+numout-1;  // index in FILTERED
00970     } // for all blocks in this buffer
00971   #endif // check index calling routine
00972   } // END goldstein phase filter
00973 
00974 
00975 
00976 /****************************************************************
00977  * B = smooth(A,KERNEL)                                         *
00978  * (circular) spatial moving average with a (2N+1,2N+1) block.  *
00979  * See also matlab script smooth.m for some tests.              *
00980  * implementation as convolution with FFT's                     *
00981  * input: KERNEL is the FFT of the kernel (block)               *
00982  #%// BK 26-Oct-2000                                            *
00983  ****************************************************************/
00984 matrix<real4> smooth(
00985         const matrix<real4> &A,
00986         const matrix<complr4> &KERNEL2D)
00987   {
00988   #ifdef __DEBUGMAT2
00989   matDEBUG.print("smooth (BK 26-Oct-2000)");
00990   #endif
00991   const int32 L = A.lines();
00992   const int32 P = A.pixels();
00993  
00994   matrix<complr4> DATA = mat2cr4(A);            // or define fft(R4)
00995   fft2d(DATA);                                  // or define fft(R4)
00996   // ______ create kernel in calling routine, e.g., like ______
00997   // ______ Kernel has to be even! ______
00998   //matrix<complr4> kernel(1,L);                        // init to zeros
00999   //for (register int32 ii=-N; ii<=N; ++ii)     // 1d kernel function of block
01000   //  kernel(0,(ii+L)%L) = 1./(2*N+1);
01001   //matrix<complr4> KERNEL2D = matTxmat(kernel,kernel);
01002   //fft2d(KERNEL2D);                            // should be real sinc
01003   DATA *= KERNEL2D;                             // no need for conj. with real fft...
01004   ifft2d(DATA);                                 // convolution, but still complex...
01005   return real(DATA);                            // you know it is real only...
01006   } // END smooth
01007 
01008 
01009 
01010 /****************************************************************
01011  * B = smooth(A,blocksize)                                      *
01012  * (circular) spatial moving average with a (2N+1,2N+1) block.  *
01013  * See also matlab script smooth.m for some tests.              *
01014  #%// BK 26-Oct-2000                                            *
01015  ****************************************************************/
01016 matrix<real4> smooth(
01017         const matrix<real4> &A,
01018         int32 N)
01019   {
01020   #ifdef __DEBUGMAT2
01021   matDEBUG.print("ROUTINE: smooth (BK 26-Oct-2000)");
01022   #endif
01023   // ______ Check if smoothing desired ______
01024   if (N==0) return A;
01025 
01026   const int32 L = A.lines();
01027   const int32 P = A.pixels();
01028 
01029   // define one of SPACEDOMAIN and SPECTRALDOMAIN to select algorithm
01030   // #define SPACEDOMAIN
01031   #define SPECTRALDOMAIN
01032 
01033   // ====== In space domain ======
01034   #ifdef SPACEDOMAIN
01035   matrix<real4> SMOOTH(L,P);            // init to zero...
01036   register real8 sum = 0.;
01037   register int32 indexii;
01038   const real8 Nsmooth = (2*N+1)*(2*N+1);
01039   for (register int32 i=0; i<L; ++i)
01040     {
01041     for (register int32 j=0; j<P; ++j)
01042       {
01043       // ______ Smooth this pixel ______
01044       for (register int32 ii=-N; ii<=N; ++ii)
01045         {
01046         indexii = (i+ii+L)%L;
01047         for (register int32 jj=-N; jj<=N; ++jj)
01048           {
01049           sum += A(indexii,(j+jj+P)%P);
01050           }
01051         }
01052       SMOOTH(i,j) = real4(sum/Nsmooth);
01053       sum = 0.;
01054       }
01055     }
01056   return SMOOTH;
01057   #endif
01058 
01059   // ====== Do the same but faster ======
01060   // ______ FFT's, but overhead due to conversion r4<->cr4 ______
01061   #ifdef SPECTRALDOMAIN
01062   matrix<complr4> DATA = mat2cr4(A);            // or define fft(R4)
01063   fft2d(DATA);                                  // or define fft(R4)
01064   matrix<complr4> kernel(1,L);                  // init to zeros
01065   for (register int32 ii=-N; ii<=N; ++ii)       // 1d kernel function of block
01066     kernel(0,(ii+L)%L) = complr4(1.0/(2*N+1),0.0);// BK 07-Apr-2003
01067   matrix<complr4> KERNEL2D = matTxmat(kernel,kernel);
01068   fft2d(KERNEL2D);                              // should be real sinc
01069   DATA *= KERNEL2D;                             // no need for conj. with real fft...
01070   ifft2d(DATA);                                 // convolution, but still complex...
01071   return real(DATA);                            // you know it is real only...
01072   #endif
01073   } // END smooth
01074 
01075 
01076 
01077 /****************************************************************
01078  *    spatialphasefilt                                          *
01079  * For the first block the part [0:OVERLAP-1] is set to 0.      *
01080  * For the last block the part [NPIX-1-OVERLAP:NPIX-1] is 0.    *
01081  #%// BK 30-Oct-2000                                            *
01082  ****************************************************************/
01083 void spatialphasefilt(
01084         const input_gen       &generalinput,
01085         const productinfo     &interferogram,
01086         const input_filtphase &filtphaseinput)
01087   {
01088   TRACE_FUNCTION("spatialphasefilt (BK 30-Oct-2000)")
01089 
01090   // ====== Handle input ======
01091   bool doexternalfile = false;
01092   char infile[EIGHTY];                          // file 2b filtered
01093   strcpy(infile,interferogram.file);
01094   int32 numlinesinput  = int32(interferogram.win.lines()/interferogram.multilookL);
01095   int32 numpixelsinput = int32(interferogram.win.pixels()/interferogram.multilookP);
01096 
01097   if (specified(filtphaseinput.fifiltphase))
01098     {
01099     doexternalfile = true;
01100     numlinesinput  = filtphaseinput.finumlines;
01101     strcpy(infile,filtphaseinput.fifiltphase);
01102     ifstream ifile;
01103     openfstream(ifile,infile);
01104     bk_assert(ifile,infile,__FILE__,__LINE__);
01105     ifile.seekg(0,ios::end);                    // pointer to end...
01106     const int32 sizefile     = ifile.tellg();   // opened ate
01107     ifile.close();
01108     numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;// floor
01109     if (numlinesinput*numpixelsinput*sizeof(complr4) != sizefile)
01110       WARNING.print("Format infile not CR4, or numlines not correct.");
01111     INFO.print("Using input file for phase filtering, not default.");
01112     }
01113   INFO << "phasefilter: inputfile: (#l,#p): " << infile << " ("
01114        << numlinesinput << "," << numpixelsinput << ")";
01115   INFO.print();
01116 
01117   // ______ Set variables ______
01118   // ______ edge of image not filtered ______
01119   //const int32 SIZE     = 512;                         // the larger the better
01120   int32 SIZE       = 1024;                              // blocksize power of 2
01121   if (generalinput.memory<100e6) SIZE/=2;               // check memory size 512
01122   if (generalinput.memory<20e6)  SIZE/=2;               // check memory size 256
01123   while (numlinesinput  < SIZE ) SIZE/=2;
01124   while (numpixelsinput < SIZE ) SIZE/=2;
01125   int32 KERNELSIZE_LINES;
01126   int32 KERNELSIZE_PIXELS;                              // ==L for 1D
01127   int32 OVERLAP_LINES;
01128   int32 OVERLAP_PIXELS;                                 // ==L for 1D
01129 
01130   // ______ use FFT's for convolution with rect ______
01131   // ______ this could also be done static, or in the calling routine ______
01132   // ______ KERNEL2D is FFT2 of even kernel (no imag part after fft!) ______
01133   // ______ one could also use ascii input file for 2d kernel ______
01134   matrix<complr4> KERNEL2D;
01135   if (!specified(filtphaseinput.fikernel2d))
01136     {
01137     KERNELSIZE_LINES  = filtphaseinput.kernel.pixels();
01138     KERNELSIZE_PIXELS = KERNELSIZE_LINES;
01139     OVERLAP_LINES     = int32(floor(KERNELSIZE_LINES/2.));
01140     OVERLAP_PIXELS    = OVERLAP_LINES;
01141 
01142     // ______ 1d kernel function ______
01143     matrix<complr4> kernel(1,SIZE);             // init to zeros
01144     for (register int32 ii=-OVERLAP_LINES; ii<=OVERLAP_LINES; ++ii)
01145       kernel(0,(ii+SIZE)%SIZE) = 
01146                           complr4(filtphaseinput.kernel(0,ii+OVERLAP_LINES),0.0);// BK 07-Apr-2003
01147       //kernel(0,(ii+SIZE)%SIZE) = filtphaseinput.kernel(0,ii+OVERLAP_LINES);
01148     KERNEL2D = matTxmat(kernel,kernel);
01149     }
01150   else // use filename and blocksize
01151     {
01152     KERNEL2D.resize(SIZE,SIZE);
01153     INFO.print("Reading 2d kernel from file (PF_IN_KERNEL2D card).");
01154     ifstream kernel2d(filtphaseinput.fikernel2d, ios::in);
01155     bk_assert(kernel2d,filtphaseinput.fikernel2d,__FILE__,__LINE__);
01156     real4 scalefactor;
01157     kernel2d >> KERNELSIZE_LINES >> KERNELSIZE_PIXELS >> scalefactor;
01158     if (abs(scalefactor)<=EPS) scalefactor=1.;
01159     INFO << "kernel2d: file: lines: rows: scale: "
01160          << filtphaseinput.fikernel2d << " " << KERNELSIZE_LINES << " "
01161          << KERNELSIZE_PIXELS << " " << scalefactor;
01162     INFO.print();
01163     if (!isodd(KERNELSIZE_LINES))
01164       {
01165       PRINT_ERROR("2D kernel must have odd number of rows!")
01166       throw(argument_error);
01167       }
01168     if (!isodd(KERNELSIZE_PIXELS))
01169       {
01170       PRINT_ERROR("2D kernel must have odd number of columns!")
01171       throw(argument_error);
01172       }
01173     if (KERNELSIZE_LINES >SIZE)
01174       {
01175       PRINT_ERROR("2D kernel has more rows than blocksize")
01176       throw(argument_error);
01177       }
01178     if (KERNELSIZE_PIXELS>SIZE)
01179       {
01180       PRINT_ERROR("2D kernel has more columns than blocksize")
01181       throw(argument_error);
01182       }
01183     // ______ DO NOT normalize to 1! ______
01184     OVERLAP_LINES  = int32((KERNELSIZE_LINES/2));// i.e., int32() floors
01185     OVERLAP_PIXELS = int32((KERNELSIZE_PIXELS/2));// i.e., int32() floors
01186     // ______ Shift kernel so centered in space domain around pixels ______
01187     for (int32 ii=-OVERLAP_LINES; ii<=OVERLAP_LINES; ++ii)
01188       {
01189       real4 tmpvalue;
01190       char dummyline[10*ONE27];                 // take care of very large kernel
01191       kernel2d.getline(dummyline,10*ONE27,'\n');
01192       const int32 indexii = (ii+SIZE)%SIZE;
01193       for (int32 jj=-OVERLAP_PIXELS; jj<=OVERLAP_PIXELS; ++jj)
01194         {
01195         const int32 indexjj = (jj+SIZE)%SIZE;
01196         kernel2d >> tmpvalue;
01197         KERNEL2D(indexii,indexjj) = complr4(scalefactor*tmpvalue);
01198         }
01199       }
01200     } // 2Dkernel
01201 
01202   // ______ Give some info ______
01203   INFO << "buffersize: (" << SIZE << ", " << SIZE << ")";
01204   INFO.print();
01205   INFO << "kernelsize: ("
01206        << KERNELSIZE_LINES << ", " << KERNELSIZE_PIXELS << ")";
01207   INFO.print();
01208 
01209   // ______ Prepare kernel for routine convolution ______
01210   fft2d(KERNEL2D);      // input to myconv function
01211   KERNEL2D.conj();      // input to myconv function, not required here, but correct.
01212                         // since fft of even function is real only...
01213 
01214   // ______ Open output file ______
01215   ofstream ofile;
01216   openfstream(ofile,filtphaseinput.fofiltphase,generalinput.overwrit);
01217   bk_assert(ofile,filtphaseinput.fofiltphase,__FILE__,__LINE__);
01218 
01219 
01220   // ______ initialize indices ______
01221   const int32 numout   = SIZE-(2*OVERLAP_LINES);// number of output lines per buffer
01222   int32 cintlinelo     = 0;                     // index in CINT to get 1st buffer
01223   int32 cintlinehi     = SIZE-1;                // index in CINT to get 1st buffer
01224   int32 filteredlinelo = 0;                     // index in FILTERED to write (1st)
01225   int32 filteredlinehi = SIZE-OVERLAP_LINES-1;  // index in FILTERED to write
01226   bool  lastbuffer     = false;                 // not yet, just starting...
01227 
01228   if (!doexternalfile)
01229     if (interferogram.formatflag!=FORMATCR4)
01230       {
01231       PRINT_ERROR("Sorry, for phasefilter, format of interferogram must be CR4")
01232       throw(argument_error);
01233       }
01234   matrix<complr4> CINT(SIZE,numpixelsinput);    // type must be cr4!
01235 
01236   int32 lineswritten     = OVERLAP_LINES;               // first buffer extra OVERLAP lines
01237   const int32 numbuffers = numlinesinput/numout;        // approx?
01238   int32 tenpercent       = int32((numbuffers/10)+.5); // round
01239   int32 percent          = 10;
01240   if (tenpercent==0) tenpercent = 1000;                 // avoid error x%0;
01241   PROGRESS.print("filtphase:  0%");
01242 
01243   // ====== Start filter loop buffers of BLOCKSIZE ======
01244   for (int32 buffer=1; buffer<1e6; ++buffer)
01245     {
01246     // ______ Check if this will be the last buffer ______
01247     if (cintlinehi >= numlinesinput-1)          // -1 since in matrix system
01248       {
01249       lastbuffer     = true;
01250       cintlinehi     = numlinesinput-1;
01251       cintlinelo     = cintlinehi-SIZE+1;       // make sure SIZE lines are read
01252       filteredlinehi = SIZE-1;                  // write 0's OVERLAP lines
01253       const int32 lines2bwritten = numlinesinput-lineswritten;
01254       filteredlinelo = filteredlinehi - lines2bwritten + 1;
01255       if (lines2bwritten<1) 
01256         WARNING.print("panic, this will crash, lines2bwritten<1?");
01257       }
01258 
01259     // ______ Read in buffers of BLOCKSIZE lines complex interferogram ______
01260     const window windummy (0,0,0,0);            // no offset in readfile
01261     const window wincint  (cintlinelo,cintlinehi,0,numpixelsinput-1);
01262     readfile(CINT,infile,numlinesinput,wincint,windummy);
01263 
01264     // ______ Filter data in buffer of BLOCKSIZE lines ______
01265     // ______ output has same size, but write only part to disk ______
01266     matrix<complr4> FILTERED = convbuffer(CINT,KERNEL2D,OVERLAP_PIXELS);
01267 
01268     // ______ Write filtered data (i.e. part of FILTERED) ______
01269     // ______ Write overlap zero lines for first buffer and lastbuffer (circular) ______
01270     if (filteredlinelo==0)                      // firstblock
01271       for (int ii=0; ii<OVERLAP_LINES; ++ii)
01272         FILTERED.setrow(ii,complr4(0.,0.));
01273     if (lastbuffer==true)
01274       for (int ii=SIZE-OVERLAP_LINES; ii<SIZE; ++ii)
01275         FILTERED.setrow(ii,complr4(0.,0.));
01276     const window winfiltered(filteredlinelo,filteredlinehi,0,numpixelsinput-1);
01277     writefile(ofile,FILTERED,winfiltered);
01278 
01279     // ______ Check if all done ______
01280     if (lastbuffer)
01281       break;
01282 
01283     // ______ Update indexes in matrices for all buffers except last ______
01284     lineswritten  += numout;            // first block overlap added
01285     cintlinelo    += numout;            // next buffer
01286     cintlinehi    += numout;            // next buffer
01287     filteredlinelo = OVERLAP_LINES;     // index in buffer for all restbuffers
01288 
01289     // ______ Give progress message ______
01290     if (buffer%tenpercent==0)
01291       {
01292       PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
01293       PROGRESS.print();
01294       percent += 10;
01295       }
01296     } // end loop buffers
01297   ofile.close();
01298 
01299   // ______ exit if only external file was desired ______
01300   if (doexternalfile)
01301     {
01302     cerr  << "\n Finished external file phasefilter, Exiting\n";
01303     PROGRESS.print("\n Finished external file phasefilter, Exiting\n");
01304     exit(0);
01305     }
01306 
01307 
01308   // ====== Write results to file ======
01309   ofstream scratchlogfile("scratchlogfiltphase", ios::out | ios::trunc);
01310   bk_assert(scratchlogfile,"filtphase: scratchlogfiltphase",__FILE__,__LINE__);
01311   scratchlogfile
01312     << "\n\n*******************************************************************"
01313     << "\nPhase filter complex interferogram: \t"
01314     <<  infile
01315     << "\nOutput file filtered master (format): \t\t\t"
01316     <<  filtphaseinput.fofiltphase << " "
01317     << "\n..."
01318     << "\n*******************************************************************\n";
01319   scratchlogfile.close();
01320 
01321   // ______ updateproductinfo greps info from this file ______
01322   ofstream scratchresfile("scratchresfiltphase", ios::out | ios::trunc);
01323   bk_assert(scratchresfile,"scratchresfiltphase",__FILE__,__LINE__);
01324   scratchresfile
01325     << "\n\n*******************************************************************"
01326     << "\n*_Start_" << processcontrol[pr_i_filtphase]
01327     << "\n*******************************************************************"
01328     << "\nMethod_phasefilt: spatial convolution: \t"
01329 //    <<  kernel...
01330     << "\nInput_file:                           \t"
01331     <<  infile
01332     << "\nData_output_file:                     \t"
01333     <<  filtphaseinput.fofiltphase
01334     << "\nData_output_format:                   \t"
01335     <<  "complex_real4"
01336     << "\nFirst_line (w.r.t. original_master):  \t"
01337     <<  interferogram.win.linelo
01338     << "\nLast_line (w.r.t. original_master):   \t"
01339     <<  interferogram.win.linehi
01340     << "\nFirst_pixel (w.r.t. original_master): \t"
01341     <<  interferogram.win.pixlo
01342     << "\nLast_pixel (w.r.t. original_master):  \t"
01343     <<  interferogram.win.pixhi
01344     << "\nMultilookfactor_azimuth_direction:    \t"
01345     <<  interferogram.multilookL
01346     << "\nMultilookfactor_range_direction:      \t"
01347     <<  interferogram.multilookP
01348     << "\nNumber of lines (multilooked):        \t"
01349     <<  numlinesinput
01350     << "\nNumber of pixels (multilooked):       \t"
01351     <<  numpixelsinput
01352     << "\n*******************************************************************"
01353     << "\n* End_" << processcontrol[pr_i_filtphase] << "_NORMAL"
01354     << "\n*******************************************************************\n";
01355   scratchresfile.close();
01356   // ====== Tidy up ======
01357   } // END spatialphasefilt
01358 
01359 
01360 
01361 /****************************************************************
01362  *    phasefilter buffer by spatial conv. with kernel.          *
01363  * Input is matrix of SIZE (e.g. 256) lines, and N range pixels.*
01364  * Filtered OUTPUT is same size as input block.                 *
01365  * Because of overlap, only write to disk in calling routine    *
01366  * part (in matrix coord.) [OVERLAP:SIZE-OVERLAP-1]             *
01367  * (in line direction)                                          *
01368  * spatial convolution with a kernel function, such as a block  *
01369  * function 111 (1D) (By FFT's).                                *
01370  * Processing is done in blocks in range direction.             * 
01371  * For the first block the part [0:OVERLAP-1] is set to 0.      *
01372  * For the last block the part [NPIX-1-OVERLAP:NPIX-1] is 0.    *
01373  *                                                              *
01374  * Input:                                                       *
01375  *  - matrix to be filtered of blocklines * numpixs             *
01376  *  - kernel2d: fft2 of 2d spatial kernel.                      *
01377  *  - overlap: half of the kernel size, e.g., 1 for 111.        *
01378  * Output:                                                      *
01379  *  - filtered matrix.                                          *
01380  *    ifft2d(BLOCK .* KERNEL2D) is returned, so if required for *
01381  *    non symmetrical kernel, offer the conj(KERNEL2D)!         *
01382  *                                                              *
01383  #%// BK 30-Oct-2000                                            *
01384  ****************************************************************/
01385 matrix<complr4> convbuffer(
01386         const matrix<complr4> &CINT,
01387         const matrix<complr4> &KERNEL2D,
01388         const int32            OVERLAP)         // overlap in column direction
01389   {
01390   TRACE_FUNCTION("convbuffer (BK 30-Oct-2000)")
01391 
01392   //#define CHECKINDICESCALLING
01393   #ifdef CHECKINDICESCALLING
01394     return CINT;
01395   #else
01396   // ______ Allocate output matrix ______
01397   const int32 SIZE = CINT.lines();
01398   const int32 NPIX = CINT.pixels();
01399   matrix<complr4> FILTERED(SIZE,NPIX);          // allocate output (==0)
01400 
01401   // ______ Get block from buffer ______
01402   const int32 numout  = SIZE-(2*OVERLAP);       // number of output pixels per block
01403   int32 cintpixlo     = 0;                      // index in CINT to get 1st block
01404   int32 cintpixhi     = SIZE-1;                 // index in CINT to get 1st block
01405   //int32 outblockpixlo = 0;                    // index in BLOCK (only 1st block)
01406   int32 outblockpixlo = OVERLAP;                // index in block
01407   int32 outblockpixhi = SIZE-1-OVERLAP;         // index in BLOCK (except last block)
01408   int32 outpixlo      = outblockpixlo;          // index in FILTERED (1st block)
01409   int32 outpixhi      = outblockpixhi;          // index in FILTERED
01410   bool  lastblockdone = false;                  // only just started...
01411 
01412 
01413   // ====== Loop forever, stop after lastblockdone ======
01414   for (;;)      //forever
01415     {
01416     if (cintpixhi>=NPIX-1)                      // check if we are doing the last block
01417       {
01418       lastblockdone = true;
01419       cintpixhi     = NPIX-1;                   // prevent reading after file
01420       cintpixlo     = cintpixhi-SIZE+1;         // but make sure SIZE pixels are read
01421       // leave last few==0
01422       outpixhi      = NPIX-1-OVERLAP;           // index in FILTERED 2b written
01423       //outblockpixhi = SIZE-1;                 // write all to the end
01424       outblockpixlo = outblockpixhi - (outpixhi-outpixlo+1) + 1;
01425       }
01426     const window wincint     (0,SIZE-1,cintpixlo,cintpixhi);
01427     const window winblock    (0,SIZE-1,outblockpixlo,outblockpixhi);
01428     const window winfiltered (0,SIZE-1,outpixlo,outpixhi);
01429 
01430     // ______ Construct BLOCK as part of CINT ______
01431     matrix<complr4> BLOCK(wincint,CINT);
01432 
01433     //#define CHECKINDEXONLY
01434     #ifndef CHECKINDEXONLY
01435     // ______ use FFT's for convolution with kernel ______
01436     fft2d(BLOCK);
01437     BLOCK *= KERNEL2D;                  // no need for conj. cause kernel is even?
01438     ifft2d(BLOCK);                      // conv. in space domain
01439     #endif // check index blocks
01440     // ______ Set correct part that is filtered in output matrix ______
01441     FILTERED.setdata(winfiltered,BLOCK,winblock);
01442 
01443     // ______ Exit if finished ______
01444     if (lastblockdone)
01445       return FILTERED;                  // return
01446 
01447     // ______ Update indexes in matrices, will be corrected for last block ______
01448     cintpixlo    += numout;             // next block
01449     cintpixhi    += numout;             // next block
01450     outpixlo      = outpixhi+1;         // index in FILTERED, next range line
01451     outpixhi      = outpixlo+numout-1;  // index in FILTERED
01452     } // for all blocks in this buffer
01453   #endif // check index calling routine
01454   } // END convbuffer
01455 
01456 
01457 
01458 /****************************************************************
01459  *    phasefilterspectral                                       *
01460  * loop over whole file and multiply spectrum of interferogram  *
01461  * with kernel specified in input file.                         *
01462  #%// BK 31-Oct-2000                                            *
01463  ****************************************************************/
01464 void phasefilterspectral(
01465         const input_gen       &generalinput,
01466         const productinfo     &interferogram,
01467         const input_filtphase &filtphaseinput)
01468   {
01469   TRACE_FUNCTION("phasefilterspectral (BK 31-Oct-2000)")
01470 
01471   // ====== Handle input ======
01472   bool doexternalfile = false;
01473   if (specified(filtphaseinput.fifiltphase))
01474     doexternalfile = true;
01475   char infile[EIGHTY];                          // file 2b filtered
01476   strcpy(infile,interferogram.file);
01477   int32 numlinesinput  = int32(interferogram.win.lines()/interferogram.multilookL);
01478   int32 numpixelsinput = int32(interferogram.win.pixels()/interferogram.multilookP);
01479   if (doexternalfile)                           // overwrite default file
01480     {
01481     numlinesinput  = filtphaseinput.finumlines;
01482     strcpy(infile,filtphaseinput.fifiltphase);
01483     ifstream ifile;
01484     openfstream(ifile,infile);
01485     bk_assert(ifile,infile,__FILE__,__LINE__);
01486     ifile.seekg(0,ios::end);                    // pointer to end...
01487     const int32 sizefile     = ifile.tellg();   // opened ate
01488     ifile.close();
01489     numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;// floor
01490     if (numlinesinput*numpixelsinput*sizeof(complr4) != sizefile)
01491       WARNING.print("Format infile not CR4, or numlines not correct.");
01492     INFO.print("Using input file for phase filtering, not default.");
01493     }
01494   INFO << "phasefilter: inputfile: (#l,#p): " << infile << " ("
01495        << numlinesinput << "," << numpixelsinput << ")";
01496   INFO.print();
01497 
01498   // ______ Set variables ______
01499   const int32 SIZE     = filtphaseinput.blocksize;      // power of 2
01500   const int32 OVERLAP  = filtphaseinput.overlap;        // half of the overlap
01501 
01502   // ______ Open output file ______
01503   ofstream ofile;
01504   openfstream(ofile,filtphaseinput.fofiltphase,generalinput.overwrit);
01505   bk_assert(ofile,filtphaseinput.fofiltphase,__FILE__,__LINE__);
01506 
01507 
01508   if (!doexternalfile)
01509     if (interferogram.formatflag!=FORMATCR4)
01510       {
01511       PRINT_ERROR("Sorry, for phasefilter, format of interferogram must be CR4")
01512       throw(argument_error);
01513       }
01514   matrix<complr4> CINT(SIZE,numpixelsinput);    // type must be cr4!
01515 
01516   // ====== Obtain KERNEL2D to multiply spectrum by from input file ======
01517   INFO.print("Reading 2d kernel from file (PF_IN_KERNEL2D card).");
01518   ifstream kernel2d(filtphaseinput.fikernel2d, ios::in);
01519   bk_assert(kernel2d,filtphaseinput.fikernel2d,__FILE__,__LINE__);
01520   // ______ First read header of kernel2d file ______
01521   int32 KERNELSIZE_LINES  = 0;
01522   int32 KERNELSIZE_PIXELS = 0;
01523   real4 scalefactor       = 0.;
01524   kernel2d >> KERNELSIZE_LINES >> KERNELSIZE_PIXELS >> scalefactor;
01525   if (abs(scalefactor)<=EPS) scalefactor=1.;
01526   INFO << "kernel2d: file: lines: rows: scale: "
01527        << filtphaseinput.fikernel2d << " " << KERNELSIZE_LINES << " "
01528        << KERNELSIZE_PIXELS << " " << scalefactor;
01529   INFO.print();
01530   // why should it be odd for spectral filter...
01531   //if (!isodd(KERNELSIZE_LINES))  PRINT_ERROR("2D kernel must have odd number of rows!")
01532   //if (!isodd(KERNELSIZE_PIXELS)) PRINT_ERROR("2D kernel must have odd number of columns!");
01533   if (KERNELSIZE_LINES >SIZE)
01534     {
01535     PRINT_ERROR("2D kernel has more rows than blocksize");
01536     throw(argument_error);
01537     }
01538   if (KERNELSIZE_PIXELS>SIZE) 
01539     {
01540     PRINT_ERROR("2D kernel has more columns than blocksize");
01541     throw(argument_error);
01542     }
01543 
01544   // ______ Then read values of kernel2d ______
01545   // ______ Shift kernel so centered in spectral domain around zero freq. ______
01546   // ______ DO NOT normalize to 1 ______
01547   //const int32 HBS_L = int32(floor(KERNELSIZE_LINES/2));
01548   //const int32 HBS_P = int32(floor(KERNELSIZE_PIXELS/2));
01549   const int32 HBS_L = int32((KERNELSIZE_LINES/2));
01550   const int32 HBS_P = int32((KERNELSIZE_PIXELS/2));
01551   const int32 EXTRA_L = (iseven(KERNELSIZE_LINES))  ? 1:0;      // 1 less to fill
01552   const int32 EXTRA_P = (iseven(KERNELSIZE_PIXELS)) ? 1:0;      // 1 less to fill
01553   matrix<real4> KERNEL2D(SIZE,SIZE);                    // allocate THE matrix
01554   for (int32 ii=-HBS_L+EXTRA_L; ii<=HBS_L; ++ii)
01555     {
01556     char dummyline[10*ONE27];                           // prevent very large kernels
01557     kernel2d.getline(dummyline,10*ONE27,'\n');
01558     const int32 indexii = (ii+SIZE)%SIZE;
01559     for (int32 jj=-HBS_P+EXTRA_P; jj<=HBS_P; ++jj)
01560       {
01561       const int32 indexjj = (jj+SIZE)%SIZE;
01562       kernel2d >> KERNEL2D(indexii,indexjj);
01563       }
01564     }
01565   if (scalefactor!=1) KERNEL2D *= scalefactor;
01566 
01567 
01568   // ______ Initialize indices ______
01569   const int32 numout   = SIZE-(2*OVERLAP);      // number of output lines
01570   int32 cintlinelo     = 0;                     // index in CINT to get 1st buffer
01571   int32 cintlinehi     = SIZE-1;                // index in CINT to get 1st buffer
01572   int32 filteredlinelo = 0;                     // index in FILTERED to write (1st)
01573   int32 filteredlinehi = SIZE-OVERLAP-1;        // index in FILTERED to write (1st)
01574   bool  lastbuffer     = false;                 // not yet, just starting...
01575   int32 lineswritten   = OVERLAP;               // first buffer extra OVERLAP lines
01576 
01577   // ______ Progress messages ______
01578   int32 percent          = 10;
01579   //const int32 numbuffers = numlinesinput/numout;      // approx?
01580   const int32 numbuffers = numlinesinput/SIZE;  // approx for large mem?
01581   int32 tenpercent = int32((numbuffers/10)+.5); // round
01582   if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
01583   PROGRESS.print("filtphase:  0%");
01584 
01585   // ====== Start filter loop buffers of BLOCKSIZE ======
01586   for (int32 buffer=1; buffer<1e6; ++buffer)
01587     {
01588     // ______ Check if this will be the last buffer ______
01589     if (cintlinehi >= numlinesinput-1)          // -1 since in matrix system
01590       {
01591       lastbuffer     = true;
01592       cintlinehi     = numlinesinput-1;
01593       cintlinelo     = cintlinehi-SIZE+1;       // make sure SIZE lines are read
01594       filteredlinehi = SIZE-1;                  // write upto lastline for last buffer
01595       const int32 lines2bwritten = numlinesinput-lineswritten;
01596       filteredlinelo = filteredlinehi - lines2bwritten + 1;
01597       if (lines2bwritten<1) 
01598         WARNING.print("panic, this will crash, lines2bwritten<1?");
01599       }
01600 
01601     // ______ Read in buffers of BLOCKSIZE lines complex interferogram ______
01602     const window windummy (0,0,0,0);            // no offset in readfile
01603     const window wincint  (cintlinelo,cintlinehi,0,numpixelsinput-1);
01604     readfile(CINT,infile,numlinesinput,wincint,windummy);
01605 
01606     // ______ Filter data in buffer of BLOCKSIZE lines ______
01607     // ______ output has same size, but write only part to disk ______
01608     const matrix<complr4> FILTERED = spectralfilt(CINT,KERNEL2D,OVERLAP);
01609 
01610     // ______ Write filtered data ______
01611     const window winfiltered (filteredlinelo,filteredlinehi,0,numpixelsinput-1);
01612     writefile(ofile,FILTERED,winfiltered);
01613 
01614     // ______ Check if all done ______
01615     if (lastbuffer)
01616       break;
01617 
01618     // ______ Update indexes in matrices, will be corrected for last block ______
01619     lineswritten  += numout;            // first block overlap added
01620     cintlinelo    += numout;            // next buffer
01621     cintlinehi    += numout;            // next buffer
01622     filteredlinelo = OVERLAP;           // index in buffer,
01623                                         // valid for all middle buffers, but not last
01624     // ______ Give progress message ______
01625     if (buffer%tenpercent==0)
01626       {
01627       PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
01628       PROGRESS.print();
01629       percent += 10;
01630       }
01631     } // end loop buffers
01632   ofile.close();
01633 
01634   // ______ exit if only external file was desired ______
01635   if (doexternalfile)
01636     {
01637     cerr << "\n  Finished external file phasefilter, Exiting\n";
01638     PROGRESS.print("\n Finished external file phasefilter, Exiting\n");
01639     exit(0);
01640     }
01641 
01642   // ====== Write results to file ======
01643   ofstream scratchlogfile("scratchlogfiltphase", ios::out | ios::trunc);
01644   bk_assert(scratchlogfile,"filtphase: scratchlogfiltphase",__FILE__,__LINE__);
01645   scratchlogfile
01646     << "\n\n*******************************************************************"
01647     << "\nPhase filter complex interferogram: \t"
01648     <<  infile
01649     << "\nOutput file filtered master (format): \t\t\t"
01650     <<  filtphaseinput.fofiltphase << " "
01651 // TODO
01652     << "\n..."
01653     << "\n*******************************************************************\n";
01654   scratchlogfile.close();
01655 
01656   // ______ updateproductinfo greps info from this file ______
01657   ofstream scratchresfile("scratchresfiltphase", ios::out | ios::trunc);
01658   bk_assert(scratchresfile,"scratchresfiltphase",__FILE__,__LINE__);
01659   scratchresfile
01660     << "\n\n*******************************************************************"
01661     << "\n*_Start_" << processcontrol[pr_i_filtphase]
01662     << "\n*******************************************************************"
01663     << "\nMethod_phasefilt: spectral: size, overlap, kernelfile: \t"
01664 // TODO
01665     <<  SIZE << " " << OVERLAP << " " << filtphaseinput.fikernel2d
01666     << "\nInput_file:                           \t"
01667     <<  infile
01668     << "\nData_output_file:                     \t"
01669     <<  filtphaseinput.fofiltphase
01670     << "\nData_output_format:                   \t"
01671     <<  "complex_real4"
01672     << "\nFirst_line (w.r.t. original_master):  \t"
01673     <<  interferogram.win.linelo
01674     << "\nLast_line (w.r.t. original_master):   \t"
01675     <<  interferogram.win.linehi
01676     << "\nFirst_pixel (w.r.t. original_master): \t"
01677     <<  interferogram.win.pixlo
01678     << "\nLast_pixel (w.r.t. original_master):  \t"
01679     <<  interferogram.win.pixhi
01680     << "\nMultilookfactor_azimuth_direction:    \t"
01681     <<  interferogram.multilookL
01682     << "\nMultilookfactor_range_direction:      \t"
01683     <<  interferogram.multilookP
01684     << "\nNumber of lines (multilooked):        \t"
01685     <<  numlinesinput
01686     << "\nNumber of pixels (multilooked):       \t"
01687     <<  numpixelsinput
01688     << "\n*******************************************************************"
01689     << "\n* End_" << processcontrol[pr_i_filtphase] << "_NORMAL"
01690     << "\n*******************************************************************\n";
01691   scratchresfile.close();
01692   // ====== Tidy up ======
01693   } // END phasefilterspectral
01694 
01695 
01696 
01697 /****************************************************************
01698  *    phasefilter spectral                                      *
01699  * Input is matrix of SIZE (e.g. 32) lines, and N range pixels. *
01700  * Filtered OUTPUT is same size as input block.                 *
01701  * Because of overlap, only write to disk in calling routine    *
01702  * part (in matrix coord.) [OVERLAP:SIZE-OVERLAP-1]             *
01703  *                                                              *
01704  * Filtering is performed by pointwise multiplication of the    *
01705  * spectrum per block by the KERNEL2D (input).                  *
01706  * Blocks in range direction,                                   * 
01707  *                                                              *
01708  #%// BK 31-Oct-2000                                            *
01709  ****************************************************************/
01710 matrix<complr4> spectralfilt(
01711         const matrix<complr4> &CINT,
01712         const matrix<real4>   &KERNEL2D,
01713         const int32            OVERLAP)
01714   {
01715   TRACE_FUNCTION("spectralfilt (BK 31-Oct-2000)")
01716   // ______ Allocate output matrix ______
01717   const int32 SIZE = CINT.lines();
01718   const int32 NPIX = CINT.pixels();
01719   matrix<complr4> FILTERED(SIZE,NPIX);          // output
01720 
01721   // ______ Get block from buffer ______
01722   const int32 numout  = SIZE-(2*OVERLAP);       // number of output pixels
01723   int32 cintpixlo     = 0;                      // index in CINT to get 1st block
01724   int32 cintpixhi     = SIZE-1;                 // index in CINT to get 1st block
01725   int32 outblockpixlo = 0;                      // index in BLOCK (only 1st block)
01726   int32 outblockpixhi = SIZE-1-OVERLAP;         // index in BLOCK (except last block)
01727   int32 outpixlo      = outblockpixlo;          // index in FILTERED (1st block)
01728   int32 outpixhi      = outblockpixhi;          // index in FILTERED
01729   bool  lastblockdone = false;                  // only just started...
01730 
01731   // ====== Loop forever, stop after lastblockdone ======
01732   for (;;)      //forever
01733     {
01734     if (cintpixhi>=NPIX-1)                      // check if we are doing the last block
01735       {
01736       lastblockdone = true;
01737       cintpixhi     = NPIX-1;                   // prevent reading after file
01738       cintpixlo     = cintpixhi-SIZE+1;         // but make sure SIZE pixels are read
01739       outpixhi      = cintpixhi;                // index in FILTERED 2b written
01740       outblockpixhi = SIZE-1;                   // write all to the end
01741       outblockpixlo = outblockpixhi - (outpixhi-outpixlo+1) + 1;
01742       }
01743     const window wincint     (0,SIZE-1,cintpixlo,cintpixhi);
01744     const window winblock    (0,SIZE-1,outblockpixlo,outblockpixhi);
01745     const window winfiltered (0,SIZE-1,outpixlo,outpixhi);
01746 
01747     // ______ Construct BLOCK as part of CINT ______
01748     matrix<complr4> BLOCK(wincint,CINT);
01749 
01750     // ______ Get spectrum/filter/ifft ______
01751     fft2d(BLOCK);
01752     BLOCK *= KERNEL2D;                  // the filter...
01753     ifft2d(BLOCK);
01754     // ______ Set correct part that is filtered in output matrix ______
01755     FILTERED.setdata(winfiltered,BLOCK,winblock);
01756 
01757     // ______ Exit if finished ______
01758     if (lastblockdone)
01759       return FILTERED;                  // return
01760 
01761     // ______ Update indexes in matrices, will be corrected for last block ______
01762     cintpixlo    += numout;             // next block
01763     cintpixhi    += numout;             // next block
01764     outblockpixlo = OVERLAP;            // index in block, valid for all middle blocks
01765     outpixlo      = outpixhi+1;         // index in FILTERED, next range line
01766     outpixhi      = outpixlo+numout-1;  // index in FILTERED
01767     } // for all blocks in this buffer
01768   } // END spectralfilt
01769 
01770 
01771 
01772 /****************************************************************
01773  *    azimuthfilter                                             *
01774  * Loop over whole master and slave image and filter out        *
01775  * part of the spectrum that is not common.                     *
01776  * Only do zero doppler freq. offset.                           *
01777  * do not use a polynomial from header for now.                 *
01778  * (next we will, but assume image are almost coreg. in range,  *
01779  *  so f_dc polynomial can be eval. same)                       * 
01780  * Per block in azimuth [1024] use a certain overlap with the   *
01781  * next block so that same data is partially used for spectrum  *
01782  * (not sure if this is requried).                              *
01783  * Filter is composed of: DE-hamming, RE-hamming (for correct   *
01784  * new size and center of the spectrum).                        *
01785  * Trick in processor.c: First call routine as:                 *
01786  *  (generalinput,filtaziinput,master,slave)                    *
01787  * in order to process the master, and then as:                 *
01788  *  (generalinput,filtaziinput,slave,master)                    *
01789  * to filter the slave slc image.                               *
01790  #%// BK 01-Nov-2000                                            *
01791  ****************************************************************/
01792 void azimuthfilter(
01793         const input_gen       &generalinput,
01794         const input_filtazi   &filtaziinput,
01795         slcimage        &master,  // not const, fdc reset possibly
01796         slcimage        &slave)   // not const, fdc reset possibly
01797   {
01798   TRACE_FUNCTION("azimuthfilter (BK 01-Nov-2000)")
01799 
01800   // ====== Handle input ======
01801   char infile[EIGHTY];                          // file 2b filtered
01802   strcpy(infile,master.file);
01803   int32 numlinesinput  = master.currentwindow.lines();
01804   int32 numpixelsinput = master.currentwindow.pixels();
01805   INFO << "azimuthfilter: inputfile: (#l,#p): " << infile << " ("
01806        << numlinesinput << "," << numpixelsinput << ")";
01807   INFO.print();
01808 
01809   // ______ Set variables ______
01810   const int32 SIZE     = filtaziinput.fftlength;// power of 2
01811   const int32 OVERLAP  = filtaziinput.overlap;  // half of the overlap
01812   const real8 HAMMING  = filtaziinput.hammingalpha;
01813   matrix<real4> FILTER;                         // i.e., THE filter
01814 
01815 
01816   // ====== If master and slave both have almost constant Doppler, ======
01817   // ====== then reset polynomial to easy case ======
01818   // --- Check variation of Doppler for this crop; zeros in spectrum ---
01819   const real8 max_fdc_change = abs(master.prf - master.abw);// 150 Hz or so for ERS
01820   const real8 master_fdc_p0  = master.pix2fdc(master.currentwindow.pixlo);
01821   const real8 master_fdc_pN  = master.pix2fdc(master.currentwindow.pixhi);
01822   const real8 slave_fdc_p0   = slave.pix2fdc(slave.currentwindow.pixlo);
01823   const real8 slave_fdc_pN   = slave.pix2fdc(slave.currentwindow.pixhi);
01824   INFO << "Master: fDC Variation [Hz] = [" << master_fdc_p0 << ":" << master_fdc_pN << "]"; 
01825   INFO.print();
01826   INFO << "Slave: fDC Variation [Hz]  = [" << slave_fdc_p0  << ":" << slave_fdc_pN  << "]"; 
01827   INFO.print();
01828   if (abs(master_fdc_pN-master_fdc_p0)<max_fdc_change &&
01829       abs(slave_fdc_pN-slave_fdc_p0)<max_fdc_change) 
01830     {
01831     INFO << "Variation of fDC for master and slave < " << max_fdc_change;
01832     INFO.print();
01833     INFO.print("Reseting Doppler polynomial to constant of center crop:");
01834     master.f_DC_a0 = master.pix2fdc(0.5*(master.currentwindow.pixlo+master.currentwindow.pixhi));
01835     master.f_DC_a1 = 0.0;
01836     master.f_DC_a2 = 0.0;
01837     INFO << "Master: constant fDC [Hz] = " << master.f_DC_a0;
01838     INFO.print();
01839     slave.f_DC_a0  = slave.pix2fdc(0.5*(slave.currentwindow.pixlo+slave.currentwindow.pixhi));
01840     slave.f_DC_a1  = 0.0;
01841     slave.f_DC_a2  = 0.0;
01842     INFO << "Slave: constant fDC [Hz] = " << slave.f_DC_a0;
01843     INFO.print();
01844     }
01845 
01846   // ====== Find out if same filter for all columns can be used ======
01847   bool samefd0 = false;                         // assume NOT same filter
01848   if (abs(master.f_DC_a1)<EPS &&
01849       abs(master.f_DC_a2)<EPS &&
01850       abs( slave.f_DC_a1)<EPS &&
01851       abs( slave.f_DC_a2)<EPS)          // true for eSARp'ed processed raw data.
01852     samefd0 = true;
01853   else
01854     DEBUG.print("Filtering data by computing fDC for each column."); 
01855 
01856   // ______ Compute filter here if same for all columns ______
01857   // ______ Else call routine filtazibuffer to compute filt per column ______
01858   // ______ based on polynomial for doppler centroid. Then assume ______
01859   // ______ images are almost alligned for 1st column in original system ______
01860   if (samefd0==true)
01861     {
01862     DEBUG.print("Filtering data by same fDC for each column."); 
01863     const bool dohamming = (HAMMING<0.9999) ? true : false;
01864     const real8 PRF      = master.prf;          // pulse repetition freq. [Hz]
01865     const real8 ABW      = master.abw;          // azimuth band width [Hz]
01866     const real8 fDC_m    = master.f_DC_a0;      // zero doppler freq. [Hz]
01867     const real8 fDC_s    = slave.f_DC_a0;       // zero doppler freq. [Hz]
01868     const real8 fDC_mean = 0.5*(fDC_m+fDC_s);   // mean doppler centroid freq.
01869     const real8 ABW_new  = max(1.0, 
01870       2.0*(0.5*ABW-abs(fDC_m-fDC_mean)));       // new bandwidth>1.0
01871     INFO << "New Azimuth Bandwidth: " << ABW_new << " [Hz]";
01872     INFO.print();
01873     INFO << "New central frequency: " << fDC_mean << " [Hz]";
01874     INFO.print();
01875 
01876     const real4 deltaf   =  PRF/real4(SIZE);
01877     const real4 fr       = -PRF/2.0;
01878     matrix<real4> freqaxis(1,SIZE);
01879     for (int32 i=0; i<SIZE; ++i)
01880       freqaxis(0,i) = fr+real4(i)*deltaf;       // [-fr:df:fr-df]
01881     if (dohamming)
01882       {
01883       // ______ NOT a good implementation for per col., cause wshift *AND* fftshift.
01884       // ______ DE-weight spectrum at centered at fDC_m ______
01885       // ______ spectrum should be periodic! (use wshift) ______
01886       matrix<real4> inversehamming = myhamming(freqaxis,ABW,PRF,HAMMING);
01887       for (int32 i=0; i<SIZE; ++i)
01888         if (inversehamming(0,i)!=0.0)// check no zero division
01889           inversehamming(0,i) = 1.0/inversehamming(0,i);
01890       // ______ Shift this circular by myshift pixels ______
01891       int32 myshift = int32(((SIZE*fDC_m/PRF) + 0.5));  // round
01892       wshift(inversehamming,-myshift);          // center at fDC_m
01893 
01894       // ______ Newhamming is scaled and centered around new mean ______
01895       myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));             // round
01896       FILTER  = myhamming(freqaxis,
01897                           ABW_new,
01898                           PRF,HAMMING);         // fftshifted
01899       wshift(FILTER,-myshift);                  // center at fDC_mean
01900       //#define REALLYDEBUG
01901       #ifdef REALLYDEBUG
01902       WARNING.print("really debug switched on, dumping filters to ascii file.");
01903       dumpasc("INVHAM",inversehamming);
01904       dumpasc("NEWHAM",FILTER);
01905       #endif
01906       FILTER *= inversehamming; // filterslve=1/filter, or flipped in meanfdc...
01907       }
01908     else        // no weighting, but center at fDC_mean, size ABW_new
01909       {
01910       int32 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));       // round
01911       FILTER = myrect(freqaxis/real4(ABW_new)); // fftshifted
01912       wshift(FILTER,-myshift);                  // center at fDC_mean
01913       #ifdef REALLYDEBUG
01914       dumpasc("NEWHAM",FILTER);
01915       #endif
01916       }
01917     ifftshift(FILTER);                          // fftsh works on data!
01918     #ifdef REALLYDEBUG
01919     dumpasc("TOTshifted",FILTER);
01920     #endif
01921     }
01922 
01923 
01924   // ______ Open output file ______
01925   ofstream ofile;
01926   openfstream(ofile,filtaziinput.foname,generalinput.overwrit);
01927   bk_assert(ofile,filtaziinput.foname,__FILE__,__LINE__);
01928 
01929   // ______ Initialize indices ______
01930   const int32 numout   = SIZE-(2*OVERLAP);      // number of output lines
01931   int32 slclinelo      = 0;                     // index in SLCIMAGE for 1st buffer
01932   int32 slclinehi      = SIZE-1;                // index in SLCIMAGE for 1st buffer
01933   int32 filteredlinelo = 0;                     // index in FILTERED to write (1st)
01934   int32 filteredlinehi = SIZE-OVERLAP-1;        // index in FILTERED to write (1st)
01935   bool  lastbuffer     = false;                 // not yet, just starting...
01936   int32 lineswritten   = OVERLAP;               // first buffer extra OVERLAP lines
01937 
01938   // ______ Progress messages (bug report Mohanty Kamini Kanta 12/06/03) ______
01939   int32 percent          = 10;
01940   const int32 numbuffers = (numlinesinput+OVERLAP)/numout;      // approx?
01941   int32 tenpercent       = int32((numbuffers/10.0)+0.5);        // round
01942   if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
01943   PROGRESS.print("filtazi:  0%");
01944 
01945   // ====== Start filter loop buffers of BLOCKSIZE ======
01946   matrix<complr4> SLCIMAGE(SIZE,numpixelsinput);
01947   for (int32 buffer=1; buffer<1e6; ++buffer)
01948     {
01949     // ______ Check if this will be the last buffer ______
01950     if (slclinehi >= numlinesinput-1)           // -1 since in matrix system
01951       {
01952       lastbuffer     = true;
01953       slclinehi      = numlinesinput-1;
01954       slclinelo      = slclinehi-SIZE+1;        // make sure SIZE lines are read
01955       filteredlinehi = SIZE-1;                  // write upto lastline for last buffer
01956       // if there is 1 buffer but smaller than size=fftlength, it cannot work
01957       //if (slclinelo < 0) slclinelo = 0;     //KKM (not OK, BK: have to take fft)
01958       if (slclinelo < 0) 
01959         {
01960         WARNING.print("I think fftlength is too big (>size image).");       
01961         WARNING.print("...but i will continue and probably crash.");       
01962         WARNING.print("PLease run again with smaller fftlength.");       
01963         }
01964 
01965       const int32 lines2bwritten = numlinesinput-lineswritten;
01966       filteredlinelo = filteredlinehi - lines2bwritten + 1;
01967       if (lines2bwritten<1) 
01968         WARNING.print("panic, this will crash, lines2bwritten<1?");
01969       }
01970 
01971     // ______ Read buffer of BLOCKSIZE lines in SLCIMAGE ______
01972     const window windummy (0,0,0,0);            // no offset in readfile
01973     const window winslc   (slclinelo,slclinehi,0,numpixelsinput-1);
01974     //readfile(SLCIMAGE,infile,numlinesinput,winslc,windummy);
01975     // req. correct startline#:
01976     //matrix<complr4> SLCIMAGE = master.readdata(SLCIMAGE);
01977     switch (master.formatflag)
01978       {
01979       case FORMATCI2:
01980         {
01981         fileci2tomatcr4(SLCIMAGE,infile,numlinesinput,winslc,windummy);
01982         break;
01983         }
01984       case FORMATCR4:
01985         {
01986         readfile(SLCIMAGE,infile,numlinesinput,winslc,windummy);
01987         break;
01988         }
01989       default:
01990         {
01991         PRINT_ERROR("readdata::not correct format on file.")
01992         throw(unhandled_case_error);
01993         }
01994       }
01995 
01996     // ______ Filter data in buffer of BLOCKSIZE lines ______
01997     // ______ output has same size, but write only part to disk ______
01998     matrix<complr4> FILTERED;
01999     if (samefd0==true)
02000       {
02001       // ______ Same filter for all columns ______
02002       fft(SLCIMAGE,1);                          // fft foreach column
02003       FILTERED = diagxmat(FILTER,SLCIMAGE);     // filter each column
02004       ifft(FILTERED,1);                         // ifft foreach column
02005       }
02006     else
02007       {
02008       // ______ Filter depends on columns number ______
02009       FILTERED = blockazifilt(SLCIMAGE,master,slave,HAMMING);
02010       }
02011 
02012     // ______ Write filtered data ______
02013     const window winfiltered(filteredlinelo,filteredlinehi,0,numpixelsinput-1);
02014     switch (filtaziinput.oformatflag)
02015       {
02016       case FORMATCR4:
02017         {
02018         writefile(ofile,FILTERED,winfiltered);
02019         break;
02020         }
02021       // ______ Convert first to ci2 before writing to file ______
02022       case FORMATCI2:
02023         {
02024         matrix <compli16> TMP(FILTERED.lines(),FILTERED.pixels());
02025         for (int32 ii=winfiltered.linelo; ii<=winfiltered.linehi; ++ii)
02026           for (int32 jj=winfiltered.pixlo; jj<=winfiltered.pixhi; ++jj)
02027             TMP(ii,jj) = cr4toci2(FILTERED(ii,jj));
02028         writefile(ofile,TMP,winfiltered);
02029         break;
02030         }
02031       default:
02032         {
02033         PRINT_ERROR("Totally impossible, checked input.")
02034         throw(unhandled_case_error);
02035         }
02036       }
02037 
02038     // ______ Check if all done ______
02039     if (lastbuffer)
02040       break;
02041 
02042     // ______ Update indexes in matrices, will be corrected for last block ______
02043     lineswritten  += numout;            // first block overlap added
02044     slclinelo     += numout;            // next buffer
02045     slclinehi     += numout;            // next buffer
02046     filteredlinelo = OVERLAP;           // index in buffer,
02047                                         // valid for all middle buffers, but not last
02048     // ______ Give progress message ______
02049     if (buffer%tenpercent==0)
02050       {
02051       PROGRESS << "FILTAZI: " << setw(3) << percent << "%";
02052       PROGRESS.print();
02053       percent += 10;
02054       }
02055     } // end loop buffers
02056   ofile.close();
02057 
02058 
02059   // ====== Write results to file ======
02060   ofstream scratchlogfile("scratchlogfiltazi", ios::out | ios::trunc);
02061   bk_assert(scratchlogfile,"filtazi: scratchlogfiltazi",__FILE__,__LINE__);
02062   scratchlogfile
02063     << "\n\n*******************************************************************"
02064     << "\nPhase filter complex interferogram: \t"
02065     <<  infile
02066     << "\nOutput file filtered : \t\t\t"
02067     << "\n..."
02068     << "\n*******************************************************************\n";
02069   scratchlogfile.close();
02070 
02071   // ______ updateproductinfo greps info from this file ______
02072   ofstream scratchresfile("scratchresfiltazi", ios::out | ios::trunc);
02073   bk_assert(scratchresfile,"scratchresfiltazi",__FILE__,__LINE__);
02074   scratchresfile
02075     << "\n\n*******************************************************************"
02076     << "\n*_Start_" << processcontrol[pr_m_filtazi]
02077     << "\n*******************************************************************"
02078     << "\nInput_file:                           \t"
02079     <<  infile
02080     << "\nData_output_file:                     \t"
02081     <<  filtaziinput.foname
02082     << "\nData_output_format:                   \t";
02083   // ___ report of Mohanty Kamini Kanta i forget this jun12/03___
02084   if (filtaziinput.oformatflag==FORMATCR4) scratchresfile <<  "complex_real4";
02085   if (filtaziinput.oformatflag==FORMATCI2) scratchresfile <<  "complex_short";
02086   scratchresfile
02087     << "\nFirst_line (w.r.t. original_master):  \t"
02088     <<  master.currentwindow.linelo
02089     << "\nLast_line (w.r.t. original_master):   \t"
02090     <<  master.currentwindow.linehi
02091     << "\nFirst_pixel (w.r.t. original_master): \t"
02092     <<  master.currentwindow.pixlo
02093     << "\nLast_pixel (w.r.t. original_master):  \t"
02094     <<  master.currentwindow.pixhi
02095     << "\n*******************************************************************"
02096     << "\n* End_" << processcontrol[pr_m_filtazi] << "_NORMAL"
02097     << "\n*******************************************************************\n";
02098   scratchresfile.close();
02099   // ====== Tidy up ======
02100   } // END azimuthfilter
02101 
02102 
02103 
02104 /****************************************************************
02105  *    azimuth filter per block                                  *
02106  * Input is matrix of SIZE (e.g. 1024) lines, and N range pixs. *
02107  * Input is SLC of master. slave_info gives fDC polynomial      *
02108  * for slave + coarse offset. HAMMING is alpha for myhamming f. *
02109  * Filtered OUTPUT is same size as input block.                 *
02110  * Because of overlap (azimuth), only write to disk in calling  *
02111  * routine part (in matrix coord.) [OVERLAP:SIZE-OVERLAP-1]     *
02112  * = SIZE-(2*OVERLAP);  // number of output pixels              *
02113  *                                                              *
02114  * Filtering is performed in the spectral domain                *
02115  * (1DFFT over azimuth for all columns at once)                 *
02116  * Filter is different for each column due to shift in fd_c     *
02117  * doppler centroid frequency.                                  *
02118  *                                                              *
02119  * ! It should still be studied if real4 matrices are accurate  *
02120  * enough, but I guess it is (BK).                              *
02121  *                                                              *
02122  #%// BK 06-Nov-2000                                            *
02123  ****************************************************************/
02124 matrix<complr4> blockazifilt(
02125         const matrix<complr4> &SLCIMAGE,
02126         const slcimage        &master,          // PRF, BW, fd0
02127         const slcimage        &slave,           // PRF, BW, fd0
02128         const real8            HAMMING)
02129   {
02130   TRACE_FUNCTION("blockazifilt (BK 06-Nov-2000)")
02131   const int32 SIZE     = SLCIMAGE.lines();      // fftlength
02132   const int32 NCOLS    = SLCIMAGE.pixels();     // width
02133   if (NCOLS != master.currentwindow.pixels())
02134     WARNING.print("this will crash, size input matrix not ok...");
02135 
02136   // ______ Compute fDC_master, fDC_slave for all columns ______
02137   // ______ Create axis to evaluate fDC polynomial for master/slave ______
02138   // ______ fDC(column) = fdc_a0 + fDC_a1*(col/RSR) + fDC_a2*(col/RSR)^2 ______
02139   // ______ fDC = y = Ax ______
02140   // ______ Capitals indicate matrices (FDC_M <-> fDC_m) ______
02141   register int32 i;
02142   DEBUG.print("Filtering data by evaluated polynomial fDC for each column."); 
02143   matrix<real8> xaxis(1,master.currentwindow.pixels());         // lying
02144   for (i=master.currentwindow.pixlo; i<=master.currentwindow.pixhi; ++i)
02145     xaxis(0,i-master.currentwindow.pixlo) = real8(i-1.0);
02146   xaxis /= (master.rsr2x/2.0);
02147   matrix<real8> FDC_M = (master.f_DC_a1 * xaxis);
02148   FDC_M += master.f_DC_a0;
02149   FDC_M += (master.f_DC_a2 * sqr(xaxis));// TODO: better use master.pix2fdc()
02150   // Bert Kampes, 20-Apr-2005
02151 
02152 
02153   // ______ fDC_slave for same(!) columns (coarse offset). ______
02154   // ______ offset defined as: cols=colm+offsetP ______
02155   for (i=master.currentwindow.pixlo; i<=master.currentwindow.pixhi; ++i)
02156     xaxis(0,i-master.currentwindow.pixlo) = real8(i-1.0) + real8(slave.coarseoffsetP);
02157   xaxis /= (slave.rsr2x/2.0);
02158   matrix<real8> FDC_S = (slave.f_DC_a1 * xaxis);
02159   FDC_S += slave.f_DC_a0;
02160   FDC_S += (slave.f_DC_a2 * sqr(xaxis));
02161 
02162   #ifdef __DEBUG
02163   DEBUG.print("Dumping matrices fDC_m, fDC_s (__DEBUG defined)");
02164   dumpasc("fDC_m",FDC_M);
02165   dumpasc("fDC_s",FDC_S);
02166   #endif
02167   
02168   // ______ Axis for filter in frequencies ______
02169 // TODO check, rather shift, test matlab... or wshift,1D over dim1
02170 // use fft properties to shift...
02171 
02172   const bool dohamming = (HAMMING<0.9999) ? true : false;
02173   const real8 PRF      = master.prf;            // pulse repetition freq. [Hz]
02174   const real8 ABW      = master.abw;            // azimuth band width [Hz]
02175 
02176   const real4 deltaf   =  PRF/real4(SIZE);
02177   const real4 fr       = -PRF/2.0;
02178   matrix<real4> freqaxis(1,SIZE);
02179   for (int32 i=0; i<SIZE; ++i)
02180     freqaxis(0,i) = fr+real4(i)*deltaf; // [-fr:df:fr-df]
02181 
02182   matrix<real4> FILTER;                         // i.e., the filter per column
02183   matrix<real4> FILTERMAT(SIZE,NCOLS);          // i.e., THE filter
02184   for (i=0; i<NCOLS; ++i)
02185     {
02186     const real8 fDC_m    = FDC_M(0,i);          // zero doppler freq. [Hz]
02187     const real8 fDC_s    = FDC_S(0,i);          // zero doppler freq. [Hz]
02188     const real8 fDC_mean = 0.5*(fDC_m+fDC_s);   // mean doppler centroid freq.
02189     const real8 ABW_new  = max(1.0,
02190       2.0*(0.5*ABW-abs(fDC_m-fDC_mean)));       // new bandwidth > 1.0
02191     if (dohamming)
02192       {
02193       // ______ NOT a good implementation for per col., cause wshift AND fftshift.
02194       // ______ DE-weight spectrum at centered at fDC_m ______
02195       // ______ spectrum should be periodic! (use wshift) ______
02196       matrix<real4> inversehamming = myhamming(freqaxis,ABW,PRF,HAMMING);
02197       for (int32 ii=0; ii<SIZE; ++ii)
02198         if (inversehamming(0,ii)!=0.0)// check zero division
02199               inversehamming(0,ii) = 1.0/inversehamming(0,ii);
02200       // ______ Shift this circular by myshift pixels ______
02201       int32 myshift = int32(((SIZE*fDC_m/PRF) + 0.5));  // round
02202       wshift(inversehamming,-myshift);          // center at fDC_m
02203   
02204       // ______ Newhamming is scaled and centered around new mean ______
02205       myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));             // round
02206       FILTER  = myhamming(freqaxis,
02207                           ABW_new,
02208                           PRF,HAMMING);         // fftshifted
02209       wshift(FILTER,-myshift);                  // center at fDC_mean
02210       FILTER *= inversehamming;
02211       }
02212     else        // no weighting, but center at fDC_mean, size ABW_new
02213       {
02214       int32 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));       // round
02215       FILTER = myrect(freqaxis/real4(ABW_new)); // fftshifted
02216       wshift(FILTER,-myshift);                  // center at fDC_mean
02217       }
02218     ifftshift(FILTER);                          // fftsh works on data!
02219     FILTERMAT.setcolumn(i,FILTER);
02220     } // foreach column
02221 
02222   // ______ Filter slcdata ______
02223   matrix<complr4> FILTERED = SLCIMAGE;
02224   fft(FILTERED,1);                              // fft foreach column
02225   FILTERED *= FILTERMAT;                        // filter each column
02226   ifft(FILTERED,1);                             // ifft foreach column
02227   return FILTERED;
02228   } // END blockazifilt
02229 
02230 
02231 
02232 
02233 /****************************************************************
02234  *    rangefilterorbits                                         *
02235  * -getoverlap approx between master/slave.                     *
02236  * -get next line master/slave.                                 *
02237  * -compute freq. shift due to different incidence angle.       *
02238  * -filter master and slave.                                    *
02239  * -write output.                                               *
02240  #%// BK 13-Nov-2000                                            *
02241  * BUG detected by Rens Swart (26-Apr-2001):                    *
02242  *  if fftlength < width/2 then output consists of replicas     *
02243  *  number of replicas = int(width/fftlength) - 1               *
02244  ****************************************************************/
02245 void rangefiltporbits(
02246          const input_gen       &generalinput,
02247          const input_filtrange &filtrangeinput,
02248          const input_ell       &ellips,
02249          const slcimage        &master,
02250          const slcimage        &slave,
02251                orbit           &masterorbit,
02252                orbit           &slaveorbit)
02253   {
02254   TRACE_FUNCTION("rangefiltporbits (BK 13-Nov-2000)")
02255 
02256   // ______ RSR,RBW in Hz, SI ______
02257   const real8 HAMMINGA  = filtrangeinput.hammingalpha;
02258   const bool dohamming  = (HAMMINGA<0.9999) ? true : false;
02259   const real8 RSR       = 0.5*master.rsr2x;     // range sampling rate fr 18.96MHz
02260   const real8 RBW       = master.rbw;           // range band width fs 15.55MHz
02261 
02262   // ______ Use approximate overlap master slave ______
02263   // ______ s=m+offset --> m=s-offset ______
02264   // ______ prevent a-b uint (< 0) == very large ______ 
02265   window wintmp;
02266   int32 tmp     = slave.currentwindow.linelo - slave.coarseoffsetL;
02267   wintmp.linelo = (tmp<1) ? 1 : tmp;
02268   tmp           = slave.currentwindow.linehi - slave.coarseoffsetL;
02269   wintmp.linehi = (tmp<1) ? 1 : tmp;
02270   tmp           = slave.currentwindow.pixlo  - slave.coarseoffsetP;
02271   wintmp.pixlo  = (tmp<1) ? 1 : tmp;
02272   tmp           = slave.currentwindow.pixhi  - slave.coarseoffsetP;
02273   wintmp.pixhi  = (tmp<1) ? 1 : tmp;
02274   // ______ master.currentwin and wintmp are now in same system ______
02275   window overlap = getoverlap(master.currentwindow,wintmp);
02276 
02277   TRACE << "overlap window: " << overlap.linelo << " " << overlap.linehi
02278        << " " << overlap.pixlo << " " << overlap.pixhi;
02279   TRACE.print();
02280 
02281   // ______ Initialize loop parameters etc. ______
02282   int32 FFTLENGTH = filtrangeinput.fftlength;
02283   const int32 NCOLS     = overlap.pixels();             // width 2b processed
02284   while (FFTLENGTH>NCOLS)
02285     {
02286     WARNING << "FFTLENGTH>overlap: New length: "
02287          << FFTLENGTH << " --> " << FFTLENGTH/2;
02288     WARNING.print();
02289     FFTLENGTH /= 2;
02290     }
02291   const int32 NBLOCKS   = int32((NCOLS/FFTLENGTH));
02292   const int32 EXTRA     = ((NCOLS%FFTLENGTH)==0) ? 0:1; // last overlaps previous
02293 
02294   cn M;                 // master position
02295   cn S;                 // slave position
02296   cn P;                 // P on ellipsoid, perpendicular to orbits;
02297   const int32 MAXITER   = 10;
02298   const real8 CRITERPOS = 1e-6;
02299   const real8 CRITERTIM = 1e-10;
02300   INFO << "rangefiltporbits: MAXITER: "   << MAXITER   << "; "
02301        << "CRITERPOS: " << CRITERPOS << " m; "
02302        << "CRITERTIM: " << CRITERTIM << " s";
02303   INFO.print();
02304 
02305   // ====== Open output files ======
02306   ofstream of_m;
02307   openfstream(of_m,filtrangeinput.fomaster,generalinput.overwrit);
02308   bk_assert(of_m,filtrangeinput.fomaster,__FILE__,__LINE__);
02309 
02310   ofstream of_s;
02311   openfstream(of_s,filtrangeinput.foslave,generalinput.overwrit);
02312   bk_assert(of_s,filtrangeinput.foslave, __FILE__,__LINE__);
02313 
02314   // ______ Shift parameters ______
02315   register int32 i;
02316   const real4 df =  RSR/real4(FFTLENGTH);
02317   const real4 fr = -RSR/2.;
02318   matrix<real4> freqaxis(1,FFTLENGTH);
02319   for (i=0; i<freqaxis.pixels(); ++i)
02320     freqaxis(0,i) = fr+real4(i)*df;
02321   matrix<real4> inversehamming;
02322   if (dohamming)
02323     {
02324     inversehamming = myhamming(freqaxis,RBW,RSR,HAMMINGA);
02325     for (i=0; i<inversehamming.pixels(); ++i)
02326       if (inversehamming(0,i)!=0.)
02327         inversehamming(0,i)= 1./inversehamming(0,i);
02328     }
02329 
02330   // ______ Some extra info for debuggers ______
02331   DEBUG << "parameters used: NCOLS " << NCOLS   << " "
02332        << "FFTLENGTH: " << FFTLENGTH << " "
02333        << "NBLOCKS: " << NBLOCKS << " "
02334        << "EXTRA: " << EXTRA;
02335   DEBUG.print();
02336   DEBUG << "overlap window: " << overlap.linelo << " " << overlap.linehi
02337        << " " << overlap.pixlo << " " << overlap.pixhi;
02338   DEBUG.print();
02339 
02340   // ______ Progress messages ______
02341   int32 percent    = 10;
02342   int32 tenpercent = int32((overlap.lines()/10)+.5);    // round
02343   if (tenpercent==0) tenpercent = 1000;                 // avoid error: x%0
02344   PROGRESS.print("filtrange:  0%");
02345 
02346   real8 totalblocks = 0;        // to give some info at end
02347   real8 totaldeltaf = 0;        // to give some info at end
02348 
02349   // ====== Foreach line, read block, filter, write ======
02350   for (int32 line_m=overlap.linelo; line_m<=overlap.linehi; ++line_m)   // in master coord.
02351     {
02352     //const real8 m_tazi = line2ta(line_m,master.t_azi1,master.prf);
02353     const real8 m_tazi = master.line2ta(line_m);
02354     M              = masterorbit.getxyz(m_tazi);
02355     int32 pixlo_m  = overlap.pixlo;
02356     int32 line_s   = line_m + slave.coarseoffsetL;              // in slave coord. system
02357     bool lastblock = false;
02358     // ______ For this line, filter all blocks ______
02359     for (int32 rangeblock=1; rangeblock<=NBLOCKS+EXTRA; ++rangeblock)
02360       {
02361       totalblocks++;
02362       int32 pixhi_m = pixlo_m + FFTLENGTH - 1;                  // i.e., fftlength pixels
02363       if (rangeblock==NBLOCKS+EXTRA)                            // only last one
02364         {
02365         pixhi_m   = overlap.pixhi;
02366         pixlo_m   = pixhi_m - FFTLENGTH + 1;
02367         lastblock = true;
02368         }
02369 
02370 //cerr << "BERT: block, pixlo, hi: "
02371 //<< rangeblock << " " << pixlo_m << " " << pixhi_m << endl;
02372 
02373       // ______ Compute frequency shift for this block ______
02374       const real8 middlepix = real8(pixlo_m+pixhi_m)/2.; // theta,Bperp here
02375       // better compute range with m_trange...
02376       lp2xyz(line_m,middlepix,ellips,master,masterorbit,P,MAXITER,CRITERPOS);
02377       // ______ Compute xyz slave satellite from P ______
02378       real8 s_trange,s_tazi;
02379       xyz2t(s_tazi,s_trange,slave,slaveorbit,P,MAXITER,CRITERTIM);
02380       // ______ Slave position ______
02381       S = slaveorbit.getxyz(s_tazi);
02382       cn MP = M.min(P);
02383       cn SP = S.min(P);
02384       const real8 incidence1 = MP.angle(P);             // theta1
02385       const real8 incidence2 = SP.angle(P);             // theta2
02386       const real8 deltatheta = incidence1-incidence2;   // counterclockwise
02387       // ______ (approx. in Hz, see also manual) ______
02388       const real4 deltaf = real4((-SOL*deltatheta) /
02389             (master.wavelength*tan(incidence1-filtrangeinput.terrainslope)));
02390       totaldeltaf += deltaf;
02391 
02392 //cerr << "deltaf: " << deltaf << " Hz\n";
02393       // ______ Read block from disk ______
02394       const int32 pixlo_s    = pixlo_m + slave.coarseoffsetP;
02395       const int32 pixhi_s    = pixhi_m + slave.coarseoffsetP;
02396       const window win_m     (line_m,line_m,pixlo_m,pixhi_m);
02397       const window win_s     (line_s,line_s,pixlo_s,pixhi_s);
02398 
02399 //  cout << "win_m:\n";
02400 //  win_m.disp();
02401 //  cout << "win_s:\n";
02402 //  win_s.disp();
02403 
02404       matrix<complr4> MASTER = master.readdata(win_m);
02405       matrix<complr4> SLAVE  = slave.readdata(win_s);
02406       fft(MASTER,2);
02407       fft(SLAVE,2);
02408 
02409       // ______ Compose filter for master and slave, and filter both ______
02410       matrix<real4> FILTER;                             // for master, slave=fliplr(m)
02411       if (dohamming)
02412         {
02413         // ______ Newhamming is scaled and centered around new mean ______
02414         // Plus or min?, master or slave?
02415         FILTER  = myhamming(freqaxis-real4(.5*deltaf),// new center
02416                             RBW-deltaf,                 // new width
02417                             RSR,HAMMINGA);              // fftshifted
02418         FILTER *= inversehamming;                       // rect f. included
02419         }
02420       else                                              // no weighting of spectra
02421         {
02422         FILTER  = myrect((freqaxis-real4(.5*deltaf)) /
02423                           real4(RBW-deltaf));           // fftshifted
02424         }
02425       ifftshift(FILTER);        // fftsh works on data!
02426 
02427       // ______ Note that filter_s = fliplr(filter_m) ______
02428       // ______ and that this is also valid after ifftshift ______
02429       // ______ I tested what side had to be filtered: master first ______
02430       // ______ or slave first. Now it should always work. ______
02431       // ______ But comment out define if you think wrong side is filtered ______
02432       #define THISside
02433       #ifdef THISside
02434       MASTER *= FILTER;         // filtered master range spectral domain
02435       FILTER.fliplr();
02436       SLAVE *= FILTER;          // filtered slave range spectral domain
02437       #else
02438       WARNING.print("This is not as theory, sign Bperp defined otherwise?");
02439       SLAVE *= FILTER;          // filtered slave range spectral domain
02440       FILTER.fliplr();
02441       MASTER *= FILTER;         // filtered master range spectral domain
02442       #endif
02443       ifft(MASTER,2);           // back to space domain
02444       ifft(SLAVE,2);            // back to space domain
02445 
02446       // ______ Write to new output file for master and slave ______
02447       if (lastblock==false)
02448         {
02449         switch (filtrangeinput.oformatflag)
02450           {
02451           case FORMATCR4:
02452             {
02453             of_m << MASTER;
02454             of_s << SLAVE;
02455             break;
02456             }
02457           // ______ Convert first to ci2 before writing to file ______
02458           case FORMATCI2:
02459             {
02460             register int32 ii;
02461             matrix <compli16> TMP(1,MASTER.pixels());
02462             for (ii=0; ii<TMP.pixels(); ++ii)
02463               TMP(0,ii) = cr4toci2(MASTER(0,ii));
02464             of_m << TMP;
02465             for (ii=0; ii<TMP.pixels(); ++ii)
02466               TMP(0,ii) = cr4toci2(SLAVE(0,ii));
02467             of_s << TMP;
02468             break;
02469             }
02470           default:
02471             {
02472             PRINT_ERROR("Totally impossible, checked input.")
02473             throw(unhandled_case_error);
02474             }
02475           } // switch
02476         }
02477       // ______ Filtered part F(NUMBLOCKS*FFTLENGTH:pixels()-1) ______
02478       else // only write extra part, which less than fftlength
02479         {
02480         const int32 lastwritten = NBLOCKS*FFTLENGTH;
02481         const int32 numtowrite  = overlap.pixels() - lastwritten;
02482         const window wintmp
02483           (0,0,MASTER.pixels()-numtowrite,MASTER.pixels()-1);
02484         matrix<complr4> TMP_M(wintmp,MASTER);   // construct as part
02485         matrix<complr4> TMP_S(wintmp,SLAVE);    // construct as part
02486         switch (filtrangeinput.oformatflag)
02487           {
02488           case FORMATCR4:
02489             {
02490             of_m << TMP_M;
02491             of_s << TMP_S;
02492             break;
02493             }
02494           // ______ Convert first to ci2 before writing to file ______
02495           case FORMATCI2:
02496             {
02497             register int32 ii;
02498             matrix <compli16> TMP(1,TMP_M.pixels());
02499             for (ii=0; ii<TMP.pixels(); ++ii)
02500               TMP(0,ii) = cr4toci2(TMP_M(0,ii));
02501             of_m << TMP;
02502             for (ii=0; ii<TMP.pixels(); ++ii)
02503               TMP(0,ii) = cr4toci2(TMP_S(0,ii));
02504             of_s << TMP;
02505             break;
02506             }
02507           default:
02508             {
02509             PRINT_ERROR("Totally impossible, checked input.")
02510             throw(unhandled_case_error);
02511             }
02512           } // switch format
02513         } // write output
02514 
02515         // ______ Update block index ______
02516         // pixlo_m   = pixhi_m - FFTLENGTH + 1;
02517         pixlo_m   = pixhi_m + 1; // BUGFIX
02518       } // all blocks
02519 
02520     // ______ Give progress message ______
02521     if ((line_m-overlap.linelo+1%tenpercent)==0)
02522       {
02523       PROGRESS << "FILTRANGE: " << setw(3) << percent << "%";
02524       PROGRESS.print();
02525       percent += 10;
02526       }
02527     } // all lines of overlap
02528 
02529   // ====== Close files and finish up ======
02530   of_m.close();
02531   of_s.close();
02532   INFO << "Mean frequency shift over " << totalblocks << "blocks = "
02533        <<  (totaldeltaf/1e6)/totalblocks << " MHz.";
02534   INFO.print();
02535 
02536 
02537   // ====== Write results to file ======
02538   ofstream scratchlogfile("scratchlogfiltrange", ios::out | ios::trunc);
02539   bk_assert(scratchlogfile,"filtrange: scratchlogfiltrange",__FILE__,__LINE__);
02540   scratchlogfile
02541     << "\n\n*******************************************************************"
02542     << "\nRange filtering based on orbits for master and slave:"
02543     << "\nInput file master (format): \t\t\t"
02544     <<  master.file << " " << master.formatflag
02545     << "\nInput file slave (format): \t\t\t"
02546     <<  slave.file << " " << slave.formatflag
02547     << "\nOutput file filtered master (format): \t\t\t"
02548     <<  filtrangeinput.fomaster << " ";
02549     if (filtrangeinput.oformatflag==FORMATCR4)
02550       scratchlogfile << "complex r4";
02551     if (filtrangeinput.oformatflag==FORMATCI2)
02552       scratchlogfile << "complex short int";
02553   scratchlogfile
02554     << "\nOutput file filtered slave (format): \t\t\t"
02555     <<  filtrangeinput.foslave << " ";
02556     if (filtrangeinput.oformatflag==FORMATCR4)
02557       scratchlogfile << "complex r4";
02558     if (filtrangeinput.oformatflag==FORMATCI2)
02559       scratchlogfile << "complex short int";
02560   scratchlogfile
02561     << "\nmean shift: \t\t\t"
02562     <<  (totaldeltaf/1e6)/totalblocks << " MHz."
02563     << "\n*******************************************************************\n";
02564   scratchlogfile.close();
02565 
02566   for (int32 fileID=1; fileID<=2; ++fileID)
02567     {
02568     char oresfile[EIGHTY];
02569     char odatafile[EIGHTY];
02570     char odataformat[EIGHTY];
02571     char processcf[EIGHTY];
02572     int offsetL = 0;
02573     int offsetP = 0;
02574     if (filtrangeinput.oformatflag==FORMATCR4)
02575       strcpy(odataformat,"complex_real4");
02576     else if (filtrangeinput.oformatflag==FORMATCI2)
02577       strcpy(odataformat,"complex_short");
02578     else
02579       {
02580       PRINT_ERROR("case error")
02581       throw(unhandled_case_error);
02582       }
02583     switch (fileID)
02584       {
02585       case 1:                                                   // master
02586         strcpy(oresfile,"scratchresMfiltrange");
02587         strcpy(odatafile,filtrangeinput.fomaster);
02588         strcpy(processcf,processcontrol[pr_m_filtrange]);       // control flag
02589         break;
02590       case 2:                                                   // slave
02591         strcpy(oresfile,"scratchresSfiltrange");
02592         strcpy(odatafile,filtrangeinput.foslave);
02593         strcpy(processcf,processcontrol[pr_s_filtrange]);       // control flag
02594         offsetL = slave.coarseoffsetL;
02595         offsetP = slave.coarseoffsetP;
02596         break;
02597       default:
02598         {
02599         PRINT_ERROR("panic: ID!={1,2}")
02600         throw(unhandled_case_error);
02601         }
02602       }
02603 
02604     // ______ updateproductinfo greps info from this file ______
02605     ofstream scratchresfile(oresfile, ios::out | ios::trunc);
02606     bk_assert(scratchresfile,oresfile,__FILE__,__LINE__);
02607     scratchresfile
02608       << "\n\n*******************************************************************"
02609       << "\n*_Start_" << processcf
02610       << "\n*******************************************************************"
02611       << "\nMethod_rangefilt:                     \t"
02612       << "porbits"
02613       << "\nData_output_file:                     \t"
02614       <<  odatafile
02615       << "\nData_output_format:                   \t"
02616       <<  odataformat
02617       // ______ s=m+offset --> m=s-offset ______
02618       // ______ still in own system! ______
02619       << "\nFirst_line (w.r.t. original_image):   \t"
02620       <<  overlap.linelo + offsetL
02621       << "\nLast_line (w.r.t. original_image):    \t"
02622       <<  overlap.linehi + offsetL
02623       << "\nFirst_pixel (w.r.t. original_image):  \t"
02624       <<  overlap.pixlo  + offsetP
02625       << "\nLast_pixel (w.r.t. original_image):   \t"
02626       <<  overlap.pixhi  + offsetP
02627       << "\n*******************************************************************"
02628       << "\n* End_" << processcf << "_NORMAL"
02629       << "\n*******************************************************************\n";
02630     scratchresfile.close();
02631     } // for master and slave
02632   } // END rangefiltporbits
02633 
02634 
02635 
02636 

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