#!/usr/bin/perl
use FileHandle;
use File::Basename;
use Getopt::Long;

$model = 0;
$dr = 1;
$daz = 1;
$mode = 0;
$scale = .7;    #display intensity relative scale factor
$exp = .35;	#display intensity exponent
$rmli_dir = "";
$diff_tab = '';
$pcycle = 6.28318530718;  #phase cycle = 28PI
$mask = "-";

#command line options parsing, the rest of the arguments are passed to @ARGV
GetOptions ('p=f' => \$pcycle,'d=s' => \$diff_tab, 'b=s' => \$rmli_dir, 'e=f' => \$exp, 's=f' => \$scale);

if (($#ARGV + 1) < 10){die <<EOS ;}
*** $0
*** Copyright 2015, Gamma Remote Sensing, v1.8 10-Apr-2015 clw ***
*** Estimate and subtract quadratic polynomial phase model from a stack of differential interferograms ***

usage: $0 <SLC_tab> <itab> <mli> <diff_dir1> <diff_dir2> <unw_type> <model> <mode> <roff> <loff> [nr] [nl] [mask] [dr] [daz] 

    SLC_tab    (input) two column list of co-registered SLC filenames and SLC parameter filenames
             	 1. SLC filename  (includes path)
	         2. SLC parameter filename (includes path)
    itab       (input) table associating interferogram stack records with pairs of RSLC_tab records (text)		
             	 1. row number in SLC_tab of the reference SLC 
             	 2. row number in SLC_tab of SLC-2 of the interferogram
             	 3. line number in the itab
             	 4. flag used to indicate if this interferogram is used for IPTA processing (0:not used  1:used)
    mli        (input) MLI image file with same range and azimuth looks as the interferogram (float)
    diff_dir1  input directory containing input differential interferograms
    diff_dir2  output directory containing detrended differential interferograms
    unw_type   unwrapped phase data used for L.S. model fit:
                 0: unwrapped phase (*.diff.unw)
	         1: unwrapped phase with adf filter (*.adf.diff.unw)
    model      polynomial phase model:
                 0: a0 + a1*y + a2*x + a3*x*y + a4*x^2 + a5*y^2 (default)
                 1: a0 + a4*x^2 + a5*y^2
                 2: a0 + a1*y + a2*x + a3*x*y
                 3: a0 + a1*y + a2*x
                 4: a0 + a2*x + a4*x^2
                 5: a0 + a2*x
    mode       processing mode:
                 0: subtract the phase model from the fcomplex differential interferogram (default)
                 1: subtract the phase model from the unwrapped interferogram
    roff       offset to starting sample for the reference region (samples)
    loff       offset to starting line (lines)
    nr         number of range samples for the reference region (enter - for default: 3)
    nl         number of lines to copy (enter - for default: 3)
    mask       raster mask file (Sun raster or BMP), 0 valued pixels are excluded from the L.S. fit (enter - for none)
    dr         range sample increment for fit values (default: $dr)
    daz        azimuth line increment (default: $daz)

    -b rmli_dir (option) use MLI2 as the background image for display rather than MLI image specified on the command line
    -d diff_tab (option) output a DIFF_tab file containing 2 column list of unwrapped diff. interferograms and delta_T values in decimal days
    -p pcycle   (option) set phase cycle for output display,  (default: $pcycle)
    -s scale  (option) set image display scale factor (default: $scale)   
    -e exp    (option) set image display exponent (default: $exp)                     
EOS

$rslc_tab  = $ARGV[0];
$itab      = $ARGV[1];
$rmli      = $ARGV[2];
$diff_dir1 = $ARGV[3];
$diff_dir2 = $ARGV[4];
$unw_type  = $ARGV[5];
$model     = $ARGV[6];
$mode      = $ARGV[7];
$roff      = $ARGV[8];
$loff      = $ARGV[9];

($unw_type >= 0  && $unw_type <= 1) or die "\nERROR: invalid unwrapped data type: $unw_type\n";
($model <= 5 && $model >= 0) or die "\nERROR: invalid polynomial phase model: $model\n";
($mode <=1  && $mode >= 0) or die "\nERROR $0: invalid processing mode: $mode\n";

if($#ARGV >= 10){
  $nr = $ARGV[10];
}

if($#ARGV >= 11){
  $nl= $ARGV[11];
}

if($#ARGV >= 12){
  $mask = $ARGV[12];
  -e $mask || $mask =~ /-/ or die "ERROR $0: quad_fit mask file does not exist: $mask\n";
}

if($#ARGV >= 13){
  $dr = $ARGV[13];
  ($dr > 0) or die "\nERROR: invalid value for quad_fit range step parameter dr: $dr\n";
}

if($#ARGV >= 14){
  $daz = $ARGV[14];
  ($daz > 0) or die "\nERROR: invalid value for quad_fit azimuth step parameter daz: $daz\n";
}

if($diff_tab ne ''){
  open(DIFF_TAB,">$diff_tab") or die "ERROR $0: cannot open DIFF_tab file: $log\n";
}

open(RSLC_TAB, "<$rslc_tab") or die "\nERROR $0: table of resampled SLC images and parameter files does not exist: $rslc_tab\n\n";
print "resampled SLC images and parameter files: $rslc_tab\n"; 

open(ITAB, "<$itab") or die "\nERROR $0: table of interferometric pairs does not exist: $itab\n\n";
print "table of interferometric pairs: $itab\n"; 

-e $rmli or die "ERROR $0: RMLI intensity image file does not exist: $rmli\n";
-d $diff_dir1 or die "ERROR $0: input differential interferogram directory does not exist: $diff_dir1\n";
($diff_dir1 ne $diff_dir2) or die "\nERROR $0: input and output directories are the same, they must be different!\n";

unless (-d $diff_dir2) {
  print "creating output directory: $diff_dir2\n";
  $exit = system("mkdir $diff_dir2");
  $exit == 0 or die "ERROR $0: non-zero exit status: mkdir $diff_dir2\n"
}

$log = "$diff_dir2/mk_quad_2d.log";  
open(LOG,">$log") or die "ERROR $0: cannot open mk_quad_2d log file: $log\n";
LOG->autoflush;
print LOG "$0 @ARGV\n\n";	#print command line that was used

$time = localtime;
print "start: $time  log file: $log\n";
print "RSLC_tab: $rslc_tab\n";
print "itab:     $itab\n";
print "input diff. directory:  $diff_dir1\n";
print "output diff. directory: $diff_dir2\n";
print "unw. phase file type:   $unw_type\n";
print "phase model:   $model\n";
print "mask file:     $mask\n";
print "sample range step: $dr    azimuth step: $daz\n";
print "processing mode mode:    $mode\n";
print "phase reference region range: $roff  azimuth: $loff\n";
print "phase reference region size range: $nr   azimuth: $naz\n";
print "quad_fit phase cycle:   $pcycle\n";
print "display scale factor: $scale  exponent: $exp\n\n";

print LOG "start: $time  log file: $log\n";
print LOG "RSLC_tab: $rslc_tab\n";
print LOG "itab:     $itab\n";
print LOG "input diff. directory:  $diff_dir1\n";
print LOG "output diff. directory: $diff_dir2\n";
print LOG "unw. phase file type:   $unw_type\n";
print LOG "phase model:   $model\n";
print LOG "mask file:     $mask\n";
print LOG "sample range step: $dr    azimuth step: $daz\n";
print LOG "processing mode mode:    $mode\n";
print LOG "phase reference region range: $roff  azimuth: $loff\n";
print LOG "phase reference region size range: $nr   azimuth: $naz\n";
print LOG "quad_fit phase cycle:   $pcycle\n";
print LOG "display scale factor: $scale  exponent: $exp\n\n";

$i = 1;
#read in the contents of the RSLC_TAB file and store in arrays
LINE: while (<RSLC_TAB>) {	#read lines of processing list file
  chomp $_;			#remove new line from record
  next LINE if /^$/; 		#skip blank lines in processing list
  next LINE if /^#/; 		#skip comments in processing list
  @fields = split;		#extract the scene identifier and other parameters on the line if present
  $slc[$i] = $fields[0];
  $slc_par[$i] = $fields[1];
  -e $slc[$i] or die "ERROR $0: SLC image file does not exist: $slc[$i]\n";     
  -e $slc_par[$i] or die "ERROR $0: ISP image parameter file does not exist for SLC: $slc_par[$i]\n";
  printf LOG "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  $i++;
} 

$n = 1;
LINE: while (<ITAB>) {		#read lines of processing list file
  chomp $_;			#remove new line from record
  next LINE if /^$/; 		#skip blank lines in processing list
  next LINE if /^#/; 		#skip comments in processing list
  @fields = split;		#extract the scene identifier and other parameters on the line if present
  next LINE if ($fields[3] == 0); #ignore interferograms with the layer flag set to 0
  $slc1 = $slc[$fields[0]];
  $slc1_par = $slc_par[$fields[0]];
  my ($s1, $dir1, $ext1) = fileparse($slc1, qr/\.[^.]*/); #extension is the last bit after the last period   

  $slc2 = $slc[$fields[1]];
  $slc2_par = $slc_par[$fields[1]];
  my ($s2, $dir2, $ext2) = fileparse($slc2, qr/\.[^.]*/); #extension is the last bit after the last period   

  @slc1_date = extract_param($slc1_par, "date:");
  @slc2_date = extract_param($slc2_par, "date:");
  @stime1 = extract_param($slc1_par, "start_time:"); #determine start time of SLC 1
  @stime2 = extract_param($slc2_par, "start_time:"); #determine start time of SLC 2
  $st1 = $stime1[1];
  $st2 = $stime2[1];

  $jd1 = julday($slc1_date[1], $slc1_date[2], $slc1_date[3]) + $st1/86400.;
  $jd2 = julday($slc2_date[1], $slc2_date[2], $slc2_date[3]) + $st2/86400.;
  $year = $slc2_date[1];
  $month = $slc2_date[2];
  $day = $slc2_date[3];
  $hrs = int($st2/3600);
  $min = int(($st2 - $hrs*3600.)/60.);
  $sec = $st2 - $hrs*3600. - $min*60.;
  $tstr = sprintf("%4d-%02d-%02dT%02d:%02d:%06.3f",$year,$month,$day,$hrs,$min,$sec);
  print "SLC-1 date: @slc1_date   julian day: $jd1\n";
  print "SLC-2 date: @slc2_date   julian day: $jd2   ISO8601: $tstr\n";
  $delta_t = $jd2 - $jd1;
  printf "interferogram time interval (days): %.8f\n", $delta_t;

  $fn1 = $diff_dir1."/".$s1."_".$s2;
  $fn2 = $diff_dir2."/".$s1."_".$s2;
  $diff_par1 = $fn1.".diff_par";
  $diff_par2 = $fn2.".diff_par";
  $off1 =  $fn1.".off";
  $off2 =  $fn2.".off";
  $plt = "-";

  if ($rmli_dir ne ""){		#background intensity image is rmli2 and changes from image to image
    $rmli2 = $rmli_dir."/".$s2.".rmli";
    unless (-e $rmli2){
      $rmli2 = $rmli_dir."/".$s2.".mli";
      -e $rmli2 or die "MLI-2 does not exist: $rmli2\n";
    }
  }
  else {
    $rmli2 = $rmli;
  }
 
  $qtmp = $fn2.".tmp";	#temporary output file
  if ($unw_type == 0){
    $diff = $fn1.".diff";
    $unw = $fn1.".unw";
    if ($mode == 0) {
      $qdiff = $fn2.".diff";
    }
    else{
      $qdiff = $fn2.".unw";
    }
  }
  
  if ($unw_type == 1){
    $diff = $fn1.".adf.diff";
    $unw = $fn1.".adf.unw";
    if ($mode == 0) {
      $qdiff = $fn2.".adf.diff";
    }
    else{
      $qdiff = $fn2.".adf.unw";
    }
  } 
  
  $time = localtime;
  printf "\n*** Diff. interf.:%5d  $time ***\n",$n;
  print "SLC-1: $slc1   SLC1_par: $slc1_par\n";
  print "SLC-2: $slc2   SLC2_par: $slc2_par\n";
  print "DIFF_par: $diff_par1\n";
  print "OFF_par:  $off1\n";
  print "Diff. interf. after subtraction of phase model: $qdiff\n"; 

  printf LOG "\n*** Diff. interf.:%5d  $time ***\n",$n;
  print LOG "SLC-1: $slc1   SLC1_par: $slc1_par\n";
  print LOG "SLC-2: $slc2   SLC2_par: $slc2_par\n";
  print LOG "DIFF_par: $diff_par1\n";
  print LOG "OFF_par:  $off1\n";
  print LOG "Diff. interf. after subtraction of phase model: $qdiff\n"; 
  @width = extract_param($off1, "interferogram_width:");
  execute("quad_fit $unw $diff_par1 $dr $daz $mask $plt $model", $log);

  if ($mode == 0){
    if ($s1 =~ m/$s2/) {
      execute("cp $diff $qdiff",$log);
    }
    else {
      execute("quad_sub $diff $diff_par1 $qtmp 1 0", $log);
      execute("cpx_math $qtmp - $qdiff $width[1] 0 $roff $loff $nr $nl - - - 1", $log);	#remove phase bias at reference
      unlink $qtmp; 
    }
    execute("rasmph_pwr $qdiff $rmli2 $width[1] 1 1 0 1 1 $scale $exp", $log);
  }
  else {
    if ($s1 =~ m/$s2/) {
      execute("cp $unw $qdiff",$log);
    }
    else {
      execute("quad_sub $unw $diff_par1 $qtmp 0 0", $log);
      execute("float_math $qtmp - $qdiff $width[1] 0 $roff $loff $nr $nl - ", $log);	#remove phase bias at reference
      unlink $qtmp;
    }
    execute("rasdt_pwr24 $qdiff $rmli2 $width[1] 1 1 0 1 1 $pcycle $scale $exp", $log);
  }
  execute("cp $diff_par1 $diff_par2", $log);
  execute("cp $off1 $off2", $log);

  if($diff_tab ne ''){
    printf DIFF_TAB "$qdiff   %13.7lf   %s\n",$delta_t, $tstr;
  }
  $n++;
}

print "\n";
$time = localtime;
print "\n$0 processing end: $time\n";
print LOG "\n$0 processing end: $time\n";
exit 0;

sub execute{
  my ($command, $log) = @_;  
  if (-e $log){open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";}
  else {open(LOG,">$log") or die "ERROR $0 : cannot open log file: $log\n";} 
  LOG->autoflush; 
  print "$command\n";
  print LOG ("\n${command}\n");
  close LOG;  
  $exit = system("$command 1>> $log");
  $exit == 0 or die "ERROR $0: non-zero exit status: $command\n"
}

sub extract_param{
  my ($infile,$keyword) = @_;
  open(PAR_IN,$infile) || die "ERROR $0: cannot open parameter file: $infile\n";

  while(<PAR_IN>){
    chomp;
    @tokens = split;
    if($tokens[0] eq $keyword){close PAR_IN; return @tokens;}
  }
  close PAR_IN;
  die "ERROR $0: keyword $keyword not found in file: $infile\n";
}

sub julday {
##################################
# Julian Day from Gregorian Date #
##################################
use integer;
my ( $y, $m, $d) = @_;
#print "y,m,d: $y $m $d\n";
my $jd = (1461 * ($y + 4800 + 
($m - 14) / 12)) / 4 + (367 * 
($m - 2 - 12 * (($m - 14) / 12))) 
/ 12 - (3 * (($y + 4900 + ($m - 14)
/ 12) / 100)) / 4 + $d - 32075;
return $jd;
}