00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 #include "matrixbk.hh"
00048 #include "constants.hh"
00049 #include "filtering.hh"
00050 #include "utilities.hh"
00051 #include "conversion.hh"
00052 #include "slcimage.hh"
00053 #include "productinfo.hh"
00054 #include "orbitbk.hh"
00055 #include "exceptions.hh"
00056
00057 #include <strstream>
00058 #include <iomanip>
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
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
00094 const real8 hamming = filtrangeinput.hammingalpha;
00095 const int32 nlmean = filtrangeinput.nlmean;
00096 const int32 fftlength = filtrangeinput.fftlength;
00097 const int32 OVERLAP = filtrangeinput.overlap;
00098 const real8 SNRthreshold = filtrangeinput.SNRthreshold;
00099 const bool doweightcorrel= filtrangeinput.doweightcorrel;
00100 const int32 oversamplefactor = filtrangeinput.oversample;
00101
00102
00103 const real8 RSR = 0.5*master.rsr2x;
00104 const real8 RBW = master.rbw;
00105
00106
00107
00108
00109 const int32 Mfilelines = master.currentwindow.lines();
00110 const int32 Ifilelines = interferogram.win.lines();
00111 const int32 numpixels = interferogram.win.pixels();
00112
00113
00114 const uint BUFFERMEMSIZE = generalinput.memory;
00115 const int32 bytesperline = numpixels * sizeof(complr4);
00116 const real4 numbigmatrices = 4.2;
00117 int32 bufferlines =
00118 int32((BUFFERMEMSIZE/numbigmatrices)/bytesperline);
00119
00120
00121
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
00132 const int32 numrangeblocks = int32(numpixels/(fftlength-2*OVERLAP));
00133 real8 overallSNR = 0.;
00134 real8 overallnotfiltered = 0.;
00135
00136
00137 int32 azimuthbuffer = 1;
00138 bool azimuthdone = false;
00139 int32 startlinethisblock = interferogram.win.linelo;
00140 window filteredpart (0, bufferlines-uint((nlmean+1)/2),
00141 0, numpixels-1);
00142
00143
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
00151
00152 for (azimuthbuffer=1; azimuthbuffer<1e5; ++azimuthbuffer)
00153 {
00154 PROGRESS << "azimuthbuffer (" << bufferlines << " lines): "
00155 << azimuthbuffer << ".";
00156 PROGRESS.print();
00157
00158
00159 int32 endlinethisblock = startlinethisblock+bufferlines-1;
00160
00161 if (endlinethisblock>=interferogram.win.linehi-nlmean-int32(bufferlines/10))
00162 {
00163 azimuthdone = true;
00164 endlinethisblock = interferogram.win.linehi;
00165 bufferlines = endlinethisblock-startlinethisblock+1;
00166 filteredpart.linehi = bufferlines - 1;
00167 }
00168
00169
00170 const window winfile(startlinethisblock, endlinethisblock,
00171 interferogram.win.pixlo,
00172 interferogram.win.pixhi);
00173
00174 const matrix<complr4> MASTER_ORIG = master.readdata(winfile);
00175 const matrix<complr4> SLAVE_ORIG = slave.readdata(winfile);
00176
00177 matrix<complr4> MASTER_FILT(MASTER_ORIG.lines(),MASTER_ORIG.pixels());
00178 matrix<complr4> SLAVE_FILT (SLAVE_ORIG.lines(), SLAVE_ORIG.pixels());
00179
00180
00181
00182 int32 startpixel = 0;
00183 int32 filteredpixlo = 0;
00184 int32 partpixlo = 0;
00185 int32 partpixhi = fftlength-1-OVERLAP;
00186 bool rangedone = false;
00187 for (;;)
00188 {
00189 int32 endpixel = startpixel+fftlength-1;
00190 int32 filteredpixhi = endpixel - OVERLAP;
00191 if (endpixel>=MASTER_ORIG.pixels()-1)
00192 {
00193 rangedone = true;
00194 endpixel = MASTER_ORIG.pixels()-1;
00195 startpixel = endpixel - fftlength + 1;
00196 filteredpixhi = endpixel;
00197 partpixhi = fftlength-1;
00198 partpixlo = partpixhi - (filteredpixhi-filteredpixlo);
00199 }
00200
00201
00202 TRACE << "range block [" << startlinethisblock << ":" << endlinethisblock
00203 << "; " << startpixel+interferogram.win.pixlo << ":"
00204 << endpixel+interferogram.win.pixlo << "]";
00205 TRACE.print();
00206
00207
00208 const window win(0,MASTER_ORIG.lines()-1,startpixel,endpixel);
00209 matrix<complr4> PARTMASTER(win,MASTER_ORIG);
00210 matrix<complr4> PARTSLAVE (win,SLAVE_ORIG);
00211
00212
00213
00214 real8 SNRmean, percentnotfiltered;
00215 rfilterblock(PARTMASTER,PARTSLAVE,
00216 nlmean,SNRthreshold,RSR,RBW,
00217 hamming,oversamplefactor,doweightcorrel,
00218 SNRmean, percentnotfiltered);
00219 overallSNR += SNRmean;
00220 overallnotfiltered += percentnotfiltered;
00221
00222
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
00229 if (rangedone==true) break;
00230 partpixlo = OVERLAP;
00231 startpixel = endpixel - 2*OVERLAP + 1;
00232 filteredpixlo = startpixel + OVERLAP;
00233 }
00234
00235
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
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
00266 if (azimuthdone==true) break;
00267 filteredpart.linelo = uint((nlmean-1.)/2.);
00268 startlinethisblock = endlinethisblock - nlmean + 2;
00269 }
00270
00271
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
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:
00330 strcpy(oresfile,"scratchresMfiltrange");
00331 strcpy(odatafile,filtrangeinput.fomaster);
00332 strcpy(processcf,processcontrol[pr_m_filtrange]);
00333 break;
00334 case 2:
00335 strcpy(oresfile,"scratchresSfiltrange");
00336 strcpy(odatafile,filtrangeinput.foslave);
00337 strcpy(processcf,processcontrol[pr_s_filtrange]);
00338 break;
00339 default:
00340 {
00341 PRINT_ERROR("panic: ID!={1,2}")
00342 throw(unhandled_case_error);
00343 }
00344 }
00345
00346
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 }
00372
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 void rfilterblock(
00408 matrix<complr4> &master,
00409 matrix<complr4> &slave,
00410 int32 nlmean,
00411 real8 SNRthreshold,
00412 real8 RSR,
00413 real8 RBW,
00414 real8 alphahamming,
00415 int32 osfactor,
00416 bool doweightcorrel,
00417 real8 &meanSNR,
00418 real8 &percentnotfiltered)
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);
00424 const int32 lastline = firstline + outputlines - 1;
00425 const bool dohamming = (alphahamming<0.9999) ? true : false;
00426
00427 const bool dooversample = (osfactor!=1) ? true : false;
00428 int32 notfiltered=0;
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
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
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?");
00492 fft(cint,2);
00493 matrix<real4> power = intensity(cint);
00494
00495
00496
00497
00498
00499
00500 if (doweightcorrel)
00501 {
00502 INFO.print("De-weighting power spectrum of convolution.");
00503
00504
00505
00506
00507
00508
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 =
00514 (npnts<indexnopeak) ? sqr(numpixs) : sqr(npnts);
00515 for (i=0; i<numlines; ++i)
00516 {
00517 power(i,j) /= weight;
00518 }
00519 }
00520 }
00521
00522
00523
00524 fft(master,2);
00525 fft(slave,2);
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;
00530 uint dummy = 0;
00531 meanSNR = 0.;
00532 real8 meanSHIFT = 0.;
00533
00534
00535
00536 for (register int32 oline=firstline; oline<=lastline; ++oline)
00537 {
00538 matrix<real4> totalp = sum(nlmeanpower,2);
00539 const real4 totalpower = totalp(0,0);
00540 const real4 maxvalue = max(nlmeanpower,dummy,shift);
00541 uint lastshift = shift;
00542 const real8 SNR = fftlength*(maxvalue/(totalpower-maxvalue));
00543 meanSNR += SNR;
00544
00545
00546 bool negshift = false;
00547 if (shift > fftlength/2)
00548 {
00549
00550 shift = uint(fftlength)-shift;
00551 lastshift = shift;
00552 negshift = true;
00553 }
00554
00555
00556 if (SNR<SNRthreshold)
00557 {
00558 notfiltered++;
00559 shift = lastshift;
00560
00561 }
00562 meanSHIFT += shift;
00563 matrix<real4> filter;
00564 if (dohamming)
00565 {
00566
00567 filter = myhamming(freqaxis-real4(.5*shift*deltaf),
00568 RBW-real8(shift*deltaf),
00569 RSR,alphahamming);
00570 filter *= inversehamming;
00571 }
00572 else
00573 {
00574 filter = myrect((freqaxis-real4(.5*shift*deltaf)) /
00575 (real4(RBW-shift*deltaf)));
00576 }
00577
00578
00579
00580 ifftshift(filter);
00581
00582
00583
00584
00585
00586
00587 if (negshift==false)
00588 {
00589 dotmult(master[oline],filter,1);
00590 filter.fliplr();
00591 dotmult(slave[oline],filter,1);
00592 }
00593 else
00594 {
00595 dotmult(slave[oline],filter,1);
00596 filter.fliplr();
00597 dotmult(master[oline],filter,1);
00598 }
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 if (oline!=lastline)
00609 {
00610 matrix<real4> line1 = power.getrow(oline-firstline);
00611 matrix<real4> lineN = power.getrow(oline-firstline+nlmean);
00612 nlmeanpower += (lineN-line1);
00613 }
00614 }
00615
00616
00617 ifft(master,2);
00618 ifft(slave,2);
00619
00620
00621 meanSHIFT /= (outputlines-notfiltered);
00622 meanSNR /= outputlines;
00623 percentnotfiltered = 100. * (real4(notfiltered)/real4(outputlines));
00624
00625
00626
00627 const real8 meanfrfreq = meanSHIFT*deltaf;
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 }
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
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
00664 bool doexternalfile = false;
00665 if (specified(filtphaseinput.fifiltphase))
00666 doexternalfile = true;
00667 char infile[EIGHTY];
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)
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);
00679 const int32 sizefile = ifile.tellg();
00680 ifile.close();
00681 numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;
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
00691 const real8 ALPHA = filtphaseinput.alpha;
00692 const int32 SIZE = filtphaseinput.blocksize;
00693 const int32 OVERLAP = filtphaseinput.overlap;
00694
00695
00696 ofstream ofile;
00697 openfstream(ofile,filtphaseinput.fofiltphase,generalinput.overwrit);
00698 bk_assert(ofile,filtphaseinput.fofiltphase,__FILE__,__LINE__);
00699
00700
00701
00702 const int32 numout = SIZE-(2*OVERLAP);
00703 int32 cintlinelo = 0;
00704 int32 cintlinehi = SIZE-1;
00705 int32 filteredlinelo = 0;
00706 int32 filteredlinehi = SIZE-OVERLAP-1;
00707 bool lastbuffer = false;
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);
00716
00717 int32 lineswritten = OVERLAP;
00718 const int32 numbuffers = numlinesinput/numout;
00719 int32 tenpercent = int32((numbuffers/10)+.5);
00720 if (tenpercent==0) tenpercent = 1000;
00721 PROGRESS.print("FILTPHASE: 0%");
00722 int32 percent = 10;
00723
00724
00725 for (int32 buffer=1; buffer<1e6; ++buffer)
00726 {
00727
00728 if (cintlinehi >= numlinesinput-1)
00729 {
00730 lastbuffer = true;
00731 cintlinehi = numlinesinput-1;
00732 cintlinelo = cintlinehi-SIZE+1;
00733 filteredlinehi = SIZE-1;
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
00741 const window windummy (0,0,0,0);
00742 const window wincint (cintlinelo,cintlinehi,0,numpixelsinput-1);
00743 readfile(CINT,infile,numlinesinput,wincint,windummy);
00744
00745
00746
00747 matrix<complr4> FILTERED = goldstein(CINT,ALPHA,OVERLAP,filtphaseinput.kernel);
00748
00749
00750 const window winfiltered(filteredlinelo,filteredlinehi,0,numpixelsinput-1);
00751 writefile(ofile,FILTERED,winfiltered);
00752
00753
00754 if (lastbuffer)
00755 break;
00756
00757
00758 lineswritten += numout;
00759 cintlinelo += numout;
00760 cintlinehi += numout;
00761 filteredlinelo = OVERLAP;
00762
00763
00764 if (buffer%tenpercent==0)
00765 {
00766 PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
00767 PROGRESS.print();
00768 percent += 10;
00769 }
00770 }
00771 ofile.close();
00772
00773
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
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
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
00834 }
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857 matrix<complr4> goldstein(
00858 const matrix<complr4> &CINT,
00859 const real4 ALPHA,
00860 const int32 OVERLAP,
00861 const matrix<real4> &smoothkernel)
00862 {
00863 TRACE_FUNCTION("ROUTINE: goldstein (BK 25-Oct-2000)")
00864
00865 #ifdef CHECKINDICESCALLING
00866 return CINT;
00867 #else
00868
00869 const int32 SIZE = CINT.lines();
00870 const int32 NPIX = CINT.pixels();
00871 matrix<complr4> FILTERED(SIZE,NPIX);
00872
00873
00874 const int32 numout = SIZE-(2*OVERLAP);
00875 int32 cintpixlo = 0;
00876 int32 cintpixhi = SIZE-1;
00877 int32 outblockpixlo = 0;
00878 int32 outblockpixhi = SIZE-1-OVERLAP;
00879 int32 outpixlo = outblockpixlo;
00880 int32 outpixhi = outblockpixhi;
00881 bool lastblockdone = false;
00882
00883 const int32 SMOOTH = int32(smoothkernel.pixels())/2;
00884 const bool dosmooth = (SMOOTH==0) ? false : true;
00885 DEBUG << "SMOOTH: " << SMOOTH;
00886 DEBUG.print();
00887
00888
00889
00890
00891 matrix<complr4> KERNEL2D;
00892 if (dosmooth==true)
00893 {
00894 matrix<complr4> kernel(1,SIZE);
00895 for (register int32 ii=-SMOOTH; ii<=SMOOTH; ++ii)
00896 {
00897
00898
00899
00900
00901 int32 tmp1 = (ii+SIZE)%SIZE;
00902 int32 tmp2 = ii+SMOOTH;
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);
00909 }
00910 DEBUG.print("kernel created for smoothing spectrum");
00911
00912
00913 for (;;)
00914 {
00915 if (cintpixhi>=NPIX-1)
00916 {
00917 lastblockdone = true;
00918 cintpixhi = NPIX-1;
00919 cintpixlo = cintpixhi-SIZE+1;
00920 outpixhi = cintpixhi;
00921 outblockpixhi = SIZE-1;
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
00929 matrix<complr4> BLOCK(wincint,CINT);
00930
00931
00932 #ifndef CHECKINDEXONLY
00933
00934 fft2d(BLOCK);
00935 matrix<real4> AMPLITUDE = magnitude(BLOCK);
00936
00937
00938 if (dosmooth==true)
00939 AMPLITUDE = smooth(AMPLITUDE,KERNEL2D);
00940
00941
00942
00943
00944 const real4 maxamplitude = max(AMPLITUDE);
00945 if (maxamplitude>1e-20)
00946 {
00947 AMPLITUDE /= maxamplitude;
00948 AMPLITUDE.mypow(ALPHA);
00949 BLOCK *= AMPLITUDE;
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
00958 FILTERED.setdata(winfiltered,BLOCK,winblock);
00959
00960
00961 if (lastblockdone)
00962 return FILTERED;
00963
00964
00965 cintpixlo += numout;
00966 cintpixhi += numout;
00967 outblockpixlo = OVERLAP;
00968 outpixlo = outpixhi+1;
00969 outpixhi = outpixlo+numout-1;
00970 }
00971 #endif // check index calling routine
00972 }
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
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);
00995 fft2d(DATA);
00996
00997
00998
00999
01000
01001
01002
01003 DATA *= KERNEL2D;
01004 ifft2d(DATA);
01005 return real(DATA);
01006 }
01007
01008
01009
01010
01011
01012
01013
01014
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
01024 if (N==0) return A;
01025
01026 const int32 L = A.lines();
01027 const int32 P = A.pixels();
01028
01029
01030
01031 #define SPECTRALDOMAIN
01032
01033
01034 #ifdef SPACEDOMAIN
01035 matrix<real4> SMOOTH(L,P);
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
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
01060
01061 #ifdef SPECTRALDOMAIN
01062 matrix<complr4> DATA = mat2cr4(A);
01063 fft2d(DATA);
01064 matrix<complr4> kernel(1,L);
01065 for (register int32 ii=-N; ii<=N; ++ii)
01066 kernel(0,(ii+L)%L) = complr4(1.0/(2*N+1),0.0);
01067 matrix<complr4> KERNEL2D = matTxmat(kernel,kernel);
01068 fft2d(KERNEL2D);
01069 DATA *= KERNEL2D;
01070 ifft2d(DATA);
01071 return real(DATA);
01072 #endif
01073 }
01074
01075
01076
01077
01078
01079
01080
01081
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
01091 bool doexternalfile = false;
01092 char infile[EIGHTY];
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);
01106 const int32 sizefile = ifile.tellg();
01107 ifile.close();
01108 numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;
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
01118
01119
01120 int32 SIZE = 1024;
01121 if (generalinput.memory<100e6) SIZE/=2;
01122 if (generalinput.memory<20e6) SIZE/=2;
01123 while (numlinesinput < SIZE ) SIZE/=2;
01124 while (numpixelsinput < SIZE ) SIZE/=2;
01125 int32 KERNELSIZE_LINES;
01126 int32 KERNELSIZE_PIXELS;
01127 int32 OVERLAP_LINES;
01128 int32 OVERLAP_PIXELS;
01129
01130
01131
01132
01133
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
01143 matrix<complr4> kernel(1,SIZE);
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);
01147
01148 KERNEL2D = matTxmat(kernel,kernel);
01149 }
01150 else
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
01184 OVERLAP_LINES = int32((KERNELSIZE_LINES/2));
01185 OVERLAP_PIXELS = int32((KERNELSIZE_PIXELS/2));
01186
01187 for (int32 ii=-OVERLAP_LINES; ii<=OVERLAP_LINES; ++ii)
01188 {
01189 real4 tmpvalue;
01190 char dummyline[10*ONE27];
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 }
01201
01202
01203 INFO << "buffersize: (" << SIZE << ", " << SIZE << ")";
01204 INFO.print();
01205 INFO << "kernelsize: ("
01206 << KERNELSIZE_LINES << ", " << KERNELSIZE_PIXELS << ")";
01207 INFO.print();
01208
01209
01210 fft2d(KERNEL2D);
01211 KERNEL2D.conj();
01212
01213
01214
01215 ofstream ofile;
01216 openfstream(ofile,filtphaseinput.fofiltphase,generalinput.overwrit);
01217 bk_assert(ofile,filtphaseinput.fofiltphase,__FILE__,__LINE__);
01218
01219
01220
01221 const int32 numout = SIZE-(2*OVERLAP_LINES);
01222 int32 cintlinelo = 0;
01223 int32 cintlinehi = SIZE-1;
01224 int32 filteredlinelo = 0;
01225 int32 filteredlinehi = SIZE-OVERLAP_LINES-1;
01226 bool lastbuffer = false;
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);
01235
01236 int32 lineswritten = OVERLAP_LINES;
01237 const int32 numbuffers = numlinesinput/numout;
01238 int32 tenpercent = int32((numbuffers/10)+.5);
01239 int32 percent = 10;
01240 if (tenpercent==0) tenpercent = 1000;
01241 PROGRESS.print("filtphase: 0%");
01242
01243
01244 for (int32 buffer=1; buffer<1e6; ++buffer)
01245 {
01246
01247 if (cintlinehi >= numlinesinput-1)
01248 {
01249 lastbuffer = true;
01250 cintlinehi = numlinesinput-1;
01251 cintlinelo = cintlinehi-SIZE+1;
01252 filteredlinehi = SIZE-1;
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
01260 const window windummy (0,0,0,0);
01261 const window wincint (cintlinelo,cintlinehi,0,numpixelsinput-1);
01262 readfile(CINT,infile,numlinesinput,wincint,windummy);
01263
01264
01265
01266 matrix<complr4> FILTERED = convbuffer(CINT,KERNEL2D,OVERLAP_PIXELS);
01267
01268
01269
01270 if (filteredlinelo==0)
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
01280 if (lastbuffer)
01281 break;
01282
01283
01284 lineswritten += numout;
01285 cintlinelo += numout;
01286 cintlinehi += numout;
01287 filteredlinelo = OVERLAP_LINES;
01288
01289
01290 if (buffer%tenpercent==0)
01291 {
01292 PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
01293 PROGRESS.print();
01294 percent += 10;
01295 }
01296 }
01297 ofile.close();
01298
01299
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
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
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
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
01357 }
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385 matrix<complr4> convbuffer(
01386 const matrix<complr4> &CINT,
01387 const matrix<complr4> &KERNEL2D,
01388 const int32 OVERLAP)
01389 {
01390 TRACE_FUNCTION("convbuffer (BK 30-Oct-2000)")
01391
01392
01393 #ifdef CHECKINDICESCALLING
01394 return CINT;
01395 #else
01396
01397 const int32 SIZE = CINT.lines();
01398 const int32 NPIX = CINT.pixels();
01399 matrix<complr4> FILTERED(SIZE,NPIX);
01400
01401
01402 const int32 numout = SIZE-(2*OVERLAP);
01403 int32 cintpixlo = 0;
01404 int32 cintpixhi = SIZE-1;
01405
01406 int32 outblockpixlo = OVERLAP;
01407 int32 outblockpixhi = SIZE-1-OVERLAP;
01408 int32 outpixlo = outblockpixlo;
01409 int32 outpixhi = outblockpixhi;
01410 bool lastblockdone = false;
01411
01412
01413
01414 for (;;)
01415 {
01416 if (cintpixhi>=NPIX-1)
01417 {
01418 lastblockdone = true;
01419 cintpixhi = NPIX-1;
01420 cintpixlo = cintpixhi-SIZE+1;
01421
01422 outpixhi = NPIX-1-OVERLAP;
01423
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
01431 matrix<complr4> BLOCK(wincint,CINT);
01432
01433
01434 #ifndef CHECKINDEXONLY
01435
01436 fft2d(BLOCK);
01437 BLOCK *= KERNEL2D;
01438 ifft2d(BLOCK);
01439 #endif // check index blocks
01440
01441 FILTERED.setdata(winfiltered,BLOCK,winblock);
01442
01443
01444 if (lastblockdone)
01445 return FILTERED;
01446
01447
01448 cintpixlo += numout;
01449 cintpixhi += numout;
01450 outpixlo = outpixhi+1;
01451 outpixhi = outpixlo+numout-1;
01452 }
01453 #endif // check index calling routine
01454 }
01455
01456
01457
01458
01459
01460
01461
01462
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
01472 bool doexternalfile = false;
01473 if (specified(filtphaseinput.fifiltphase))
01474 doexternalfile = true;
01475 char infile[EIGHTY];
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)
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);
01487 const int32 sizefile = ifile.tellg();
01488 ifile.close();
01489 numpixelsinput = sizefile/sizeof(complr4)/numlinesinput;
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
01499 const int32 SIZE = filtphaseinput.blocksize;
01500 const int32 OVERLAP = filtphaseinput.overlap;
01501
01502
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);
01515
01516
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
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
01531
01532
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
01545
01546
01547
01548
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;
01552 const int32 EXTRA_P = (iseven(KERNELSIZE_PIXELS)) ? 1:0;
01553 matrix<real4> KERNEL2D(SIZE,SIZE);
01554 for (int32 ii=-HBS_L+EXTRA_L; ii<=HBS_L; ++ii)
01555 {
01556 char dummyline[10*ONE27];
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
01569 const int32 numout = SIZE-(2*OVERLAP);
01570 int32 cintlinelo = 0;
01571 int32 cintlinehi = SIZE-1;
01572 int32 filteredlinelo = 0;
01573 int32 filteredlinehi = SIZE-OVERLAP-1;
01574 bool lastbuffer = false;
01575 int32 lineswritten = OVERLAP;
01576
01577
01578 int32 percent = 10;
01579
01580 const int32 numbuffers = numlinesinput/SIZE;
01581 int32 tenpercent = int32((numbuffers/10)+.5);
01582 if (tenpercent==0) tenpercent = 1000;
01583 PROGRESS.print("filtphase: 0%");
01584
01585
01586 for (int32 buffer=1; buffer<1e6; ++buffer)
01587 {
01588
01589 if (cintlinehi >= numlinesinput-1)
01590 {
01591 lastbuffer = true;
01592 cintlinehi = numlinesinput-1;
01593 cintlinelo = cintlinehi-SIZE+1;
01594 filteredlinehi = SIZE-1;
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
01602 const window windummy (0,0,0,0);
01603 const window wincint (cintlinelo,cintlinehi,0,numpixelsinput-1);
01604 readfile(CINT,infile,numlinesinput,wincint,windummy);
01605
01606
01607
01608 const matrix<complr4> FILTERED = spectralfilt(CINT,KERNEL2D,OVERLAP);
01609
01610
01611 const window winfiltered (filteredlinelo,filteredlinehi,0,numpixelsinput-1);
01612 writefile(ofile,FILTERED,winfiltered);
01613
01614
01615 if (lastbuffer)
01616 break;
01617
01618
01619 lineswritten += numout;
01620 cintlinelo += numout;
01621 cintlinehi += numout;
01622 filteredlinelo = OVERLAP;
01623
01624
01625 if (buffer%tenpercent==0)
01626 {
01627 PROGRESS << "FILTPHASE: " << setw(3) << percent << "%";
01628 PROGRESS.print();
01629 percent += 10;
01630 }
01631 }
01632 ofile.close();
01633
01634
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
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
01652 << "\n..."
01653 << "\n*******************************************************************\n";
01654 scratchlogfile.close();
01655
01656
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
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
01693 }
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
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
01717 const int32 SIZE = CINT.lines();
01718 const int32 NPIX = CINT.pixels();
01719 matrix<complr4> FILTERED(SIZE,NPIX);
01720
01721
01722 const int32 numout = SIZE-(2*OVERLAP);
01723 int32 cintpixlo = 0;
01724 int32 cintpixhi = SIZE-1;
01725 int32 outblockpixlo = 0;
01726 int32 outblockpixhi = SIZE-1-OVERLAP;
01727 int32 outpixlo = outblockpixlo;
01728 int32 outpixhi = outblockpixhi;
01729 bool lastblockdone = false;
01730
01731
01732 for (;;)
01733 {
01734 if (cintpixhi>=NPIX-1)
01735 {
01736 lastblockdone = true;
01737 cintpixhi = NPIX-1;
01738 cintpixlo = cintpixhi-SIZE+1;
01739 outpixhi = cintpixhi;
01740 outblockpixhi = SIZE-1;
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
01748 matrix<complr4> BLOCK(wincint,CINT);
01749
01750
01751 fft2d(BLOCK);
01752 BLOCK *= KERNEL2D;
01753 ifft2d(BLOCK);
01754
01755 FILTERED.setdata(winfiltered,BLOCK,winblock);
01756
01757
01758 if (lastblockdone)
01759 return FILTERED;
01760
01761
01762 cintpixlo += numout;
01763 cintpixhi += numout;
01764 outblockpixlo = OVERLAP;
01765 outpixlo = outpixhi+1;
01766 outpixhi = outpixlo+numout-1;
01767 }
01768 }
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792 void azimuthfilter(
01793 const input_gen &generalinput,
01794 const input_filtazi &filtaziinput,
01795 slcimage &master,
01796 slcimage &slave)
01797 {
01798 TRACE_FUNCTION("azimuthfilter (BK 01-Nov-2000)")
01799
01800
01801 char infile[EIGHTY];
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
01810 const int32 SIZE = filtaziinput.fftlength;
01811 const int32 OVERLAP = filtaziinput.overlap;
01812 const real8 HAMMING = filtaziinput.hammingalpha;
01813 matrix<real4> FILTER;
01814
01815
01816
01817
01818
01819 const real8 max_fdc_change = abs(master.prf - master.abw);
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
01847 bool samefd0 = false;
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)
01852 samefd0 = true;
01853 else
01854 DEBUG.print("Filtering data by computing fDC for each column.");
01855
01856
01857
01858
01859
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;
01865 const real8 ABW = master.abw;
01866 const real8 fDC_m = master.f_DC_a0;
01867 const real8 fDC_s = slave.f_DC_a0;
01868 const real8 fDC_mean = 0.5*(fDC_m+fDC_s);
01869 const real8 ABW_new = max(1.0,
01870 2.0*(0.5*ABW-abs(fDC_m-fDC_mean)));
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;
01881 if (dohamming)
01882 {
01883
01884
01885
01886 matrix<real4> inversehamming = myhamming(freqaxis,ABW,PRF,HAMMING);
01887 for (int32 i=0; i<SIZE; ++i)
01888 if (inversehamming(0,i)!=0.0)
01889 inversehamming(0,i) = 1.0/inversehamming(0,i);
01890
01891 int32 myshift = int32(((SIZE*fDC_m/PRF) + 0.5));
01892 wshift(inversehamming,-myshift);
01893
01894
01895 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));
01896 FILTER = myhamming(freqaxis,
01897 ABW_new,
01898 PRF,HAMMING);
01899 wshift(FILTER,-myshift);
01900
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;
01907 }
01908 else
01909 {
01910 int32 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));
01911 FILTER = myrect(freqaxis/real4(ABW_new));
01912 wshift(FILTER,-myshift);
01913 #ifdef REALLYDEBUG
01914 dumpasc("NEWHAM",FILTER);
01915 #endif
01916 }
01917 ifftshift(FILTER);
01918 #ifdef REALLYDEBUG
01919 dumpasc("TOTshifted",FILTER);
01920 #endif
01921 }
01922
01923
01924
01925 ofstream ofile;
01926 openfstream(ofile,filtaziinput.foname,generalinput.overwrit);
01927 bk_assert(ofile,filtaziinput.foname,__FILE__,__LINE__);
01928
01929
01930 const int32 numout = SIZE-(2*OVERLAP);
01931 int32 slclinelo = 0;
01932 int32 slclinehi = SIZE-1;
01933 int32 filteredlinelo = 0;
01934 int32 filteredlinehi = SIZE-OVERLAP-1;
01935 bool lastbuffer = false;
01936 int32 lineswritten = OVERLAP;
01937
01938
01939 int32 percent = 10;
01940 const int32 numbuffers = (numlinesinput+OVERLAP)/numout;
01941 int32 tenpercent = int32((numbuffers/10.0)+0.5);
01942 if (tenpercent==0) tenpercent = 1000;
01943 PROGRESS.print("filtazi: 0%");
01944
01945
01946 matrix<complr4> SLCIMAGE(SIZE,numpixelsinput);
01947 for (int32 buffer=1; buffer<1e6; ++buffer)
01948 {
01949
01950 if (slclinehi >= numlinesinput-1)
01951 {
01952 lastbuffer = true;
01953 slclinehi = numlinesinput-1;
01954 slclinelo = slclinehi-SIZE+1;
01955 filteredlinehi = SIZE-1;
01956
01957
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
01972 const window windummy (0,0,0,0);
01973 const window winslc (slclinelo,slclinehi,0,numpixelsinput-1);
01974
01975
01976
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
01997
01998 matrix<complr4> FILTERED;
01999 if (samefd0==true)
02000 {
02001
02002 fft(SLCIMAGE,1);
02003 FILTERED = diagxmat(FILTER,SLCIMAGE);
02004 ifft(FILTERED,1);
02005 }
02006 else
02007 {
02008
02009 FILTERED = blockazifilt(SLCIMAGE,master,slave,HAMMING);
02010 }
02011
02012
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
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
02039 if (lastbuffer)
02040 break;
02041
02042
02043 lineswritten += numout;
02044 slclinelo += numout;
02045 slclinehi += numout;
02046 filteredlinelo = OVERLAP;
02047
02048
02049 if (buffer%tenpercent==0)
02050 {
02051 PROGRESS << "FILTAZI: " << setw(3) << percent << "%";
02052 PROGRESS.print();
02053 percent += 10;
02054 }
02055 }
02056 ofile.close();
02057
02058
02059
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
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
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
02100 }
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124 matrix<complr4> blockazifilt(
02125 const matrix<complr4> &SLCIMAGE,
02126 const slcimage &master,
02127 const slcimage &slave,
02128 const real8 HAMMING)
02129 {
02130 TRACE_FUNCTION("blockazifilt (BK 06-Nov-2000)")
02131 const int32 SIZE = SLCIMAGE.lines();
02132 const int32 NCOLS = SLCIMAGE.pixels();
02133 if (NCOLS != master.currentwindow.pixels())
02134 WARNING.print("this will crash, size input matrix not ok...");
02135
02136
02137
02138
02139
02140
02141 register int32 i;
02142 DEBUG.print("Filtering data by evaluated polynomial fDC for each column.");
02143 matrix<real8> xaxis(1,master.currentwindow.pixels());
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));
02150
02151
02152
02153
02154
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
02169
02170
02171
02172 const bool dohamming = (HAMMING<0.9999) ? true : false;
02173 const real8 PRF = master.prf;
02174 const real8 ABW = master.abw;
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;
02181
02182 matrix<real4> FILTER;
02183 matrix<real4> FILTERMAT(SIZE,NCOLS);
02184 for (i=0; i<NCOLS; ++i)
02185 {
02186 const real8 fDC_m = FDC_M(0,i);
02187 const real8 fDC_s = FDC_S(0,i);
02188 const real8 fDC_mean = 0.5*(fDC_m+fDC_s);
02189 const real8 ABW_new = max(1.0,
02190 2.0*(0.5*ABW-abs(fDC_m-fDC_mean)));
02191 if (dohamming)
02192 {
02193
02194
02195
02196 matrix<real4> inversehamming = myhamming(freqaxis,ABW,PRF,HAMMING);
02197 for (int32 ii=0; ii<SIZE; ++ii)
02198 if (inversehamming(0,ii)!=0.0)
02199 inversehamming(0,ii) = 1.0/inversehamming(0,ii);
02200
02201 int32 myshift = int32(((SIZE*fDC_m/PRF) + 0.5));
02202 wshift(inversehamming,-myshift);
02203
02204
02205 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));
02206 FILTER = myhamming(freqaxis,
02207 ABW_new,
02208 PRF,HAMMING);
02209 wshift(FILTER,-myshift);
02210 FILTER *= inversehamming;
02211 }
02212 else
02213 {
02214 int32 myshift = int32(((SIZE*fDC_mean/PRF) + 0.5));
02215 FILTER = myrect(freqaxis/real4(ABW_new));
02216 wshift(FILTER,-myshift);
02217 }
02218 ifftshift(FILTER);
02219 FILTERMAT.setcolumn(i,FILTER);
02220 }
02221
02222
02223 matrix<complr4> FILTERED = SLCIMAGE;
02224 fft(FILTERED,1);
02225 FILTERED *= FILTERMAT;
02226 ifft(FILTERED,1);
02227 return FILTERED;
02228 }
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239
02240
02241
02242
02243
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
02257 const real8 HAMMINGA = filtrangeinput.hammingalpha;
02258 const bool dohamming = (HAMMINGA<0.9999) ? true : false;
02259 const real8 RSR = 0.5*master.rsr2x;
02260 const real8 RBW = master.rbw;
02261
02262
02263
02264
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
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
02282 int32 FFTLENGTH = filtrangeinput.fftlength;
02283 const int32 NCOLS = overlap.pixels();
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;
02293
02294 cn M;
02295 cn S;
02296 cn P;
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
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
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
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
02341 int32 percent = 10;
02342 int32 tenpercent = int32((overlap.lines()/10)+.5);
02343 if (tenpercent==0) tenpercent = 1000;
02344 PROGRESS.print("filtrange: 0%");
02345
02346 real8 totalblocks = 0;
02347 real8 totaldeltaf = 0;
02348
02349
02350 for (int32 line_m=overlap.linelo; line_m<=overlap.linehi; ++line_m)
02351 {
02352
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;
02357 bool lastblock = false;
02358
02359 for (int32 rangeblock=1; rangeblock<=NBLOCKS+EXTRA; ++rangeblock)
02360 {
02361 totalblocks++;
02362 int32 pixhi_m = pixlo_m + FFTLENGTH - 1;
02363 if (rangeblock==NBLOCKS+EXTRA)
02364 {
02365 pixhi_m = overlap.pixhi;
02366 pixlo_m = pixhi_m - FFTLENGTH + 1;
02367 lastblock = true;
02368 }
02369
02370
02371
02372
02373
02374 const real8 middlepix = real8(pixlo_m+pixhi_m)/2.;
02375
02376 lp2xyz(line_m,middlepix,ellips,master,masterorbit,P,MAXITER,CRITERPOS);
02377
02378 real8 s_trange,s_tazi;
02379 xyz2t(s_tazi,s_trange,slave,slaveorbit,P,MAXITER,CRITERTIM);
02380
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);
02385 const real8 incidence2 = SP.angle(P);
02386 const real8 deltatheta = incidence1-incidence2;
02387
02388 const real4 deltaf = real4((-SOL*deltatheta) /
02389 (master.wavelength*tan(incidence1-filtrangeinput.terrainslope)));
02390 totaldeltaf += deltaf;
02391
02392
02393
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
02400
02401
02402
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
02410 matrix<real4> FILTER;
02411 if (dohamming)
02412 {
02413
02414
02415 FILTER = myhamming(freqaxis-real4(.5*deltaf),
02416 RBW-deltaf,
02417 RSR,HAMMINGA);
02418 FILTER *= inversehamming;
02419 }
02420 else
02421 {
02422 FILTER = myrect((freqaxis-real4(.5*deltaf)) /
02423 real4(RBW-deltaf));
02424 }
02425 ifftshift(FILTER);
02426
02427
02428
02429
02430
02431
02432 #define THISside
02433 #ifdef THISside
02434 MASTER *= FILTER;
02435 FILTER.fliplr();
02436 SLAVE *= FILTER;
02437 #else
02438 WARNING.print("This is not as theory, sign Bperp defined otherwise?");
02439 SLAVE *= FILTER;
02440 FILTER.fliplr();
02441 MASTER *= FILTER;
02442 #endif
02443 ifft(MASTER,2);
02444 ifft(SLAVE,2);
02445
02446
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
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 }
02476 }
02477
02478 else
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);
02485 matrix<complr4> TMP_S(wintmp,SLAVE);
02486 switch (filtrangeinput.oformatflag)
02487 {
02488 case FORMATCR4:
02489 {
02490 of_m << TMP_M;
02491 of_s << TMP_S;
02492 break;
02493 }
02494
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 }
02513 }
02514
02515
02516
02517 pixlo_m = pixhi_m + 1;
02518 }
02519
02520
02521 if ((line_m-overlap.linelo+1%tenpercent)==0)
02522 {
02523 PROGRESS << "FILTRANGE: " << setw(3) << percent << "%";
02524 PROGRESS.print();
02525 percent += 10;
02526 }
02527 }
02528
02529
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
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:
02586 strcpy(oresfile,"scratchresMfiltrange");
02587 strcpy(odatafile,filtrangeinput.fomaster);
02588 strcpy(processcf,processcontrol[pr_m_filtrange]);
02589 break;
02590 case 2:
02591 strcpy(oresfile,"scratchresSfiltrange");
02592 strcpy(odatafile,filtrangeinput.foslave);
02593 strcpy(processcf,processcontrol[pr_s_filtrange]);
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
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
02618
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 }
02632 }
02633
02634
02635
02636