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

$type = 1;
$nr = 32;
$naz = 32;
$gcp_win = 3;
$mask = "-";

if (($#ARGV + 1) < 5){die <<EOS ;}
*** $0
*** Copyright 2015, Gamma Remote Sensing, v2.5 13-Apr-2015 clw ***
*** Estimate precision baselines from unwrapped phase and DEM in range-Doppler coordinates  ***

usage: $0 <RSLC_tab> <itab> <DEM_rdc> <diff_dir> <pbase> [mask] [nr] [naz] [gcp_win] [type] 

    RSLC_tab  (input) two column list of resampled SLC filenames and SLC parameter filenames (including paths) (ascii)
                1. SLC filename  (includes path)
                2. SLC parameter filename (includes path)
    itab      (input) table associating interferogram stack records with pairs of SLC stack records (ascii)
                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)
    DEM_rdc   (input) terrain height in radar range-Doppler coordinates (meters, float format)   
    diff_dir  directory containing unwrapped differential interferograms and baselines  
    pbase     (output) baseline parameter stack (enter - for none)
    mask      (input) mask for selection of valid GCPs (Sun raster or BMP format, enter - for none)
    nr        number of GCP selection points in range (default: $nr)
    naz       number of GCP selection points in azimuth (default: $naz)
    gcp_win   window size for averaging unwrapped phase, must be odd (default: $gcp_win)
    type      differential interferogram type:
             	0: unfiltered differential interferogram (*.unw)
	     	1: (default) adf filtered differential interferogram (*.adf.unw) 
                           
EOS

$rslc_tab  = $ARGV[0];
$itab      = $ARGV[1];
$dem       = $ARGV[2];
$diff_dir  = $ARGV[3];
$pbase     = $ARGV[4];

if($#ARGV >= 5){
  $mask = $ARGV[5];
  (-e $mask || $mask =~ m/-/ ) or die "\nERROR: raster 2D image mask does not exist: $mask\n";
}

if ($mask eq "-"){
  $mask = " ";			#set mask to blank, if input is - 
}

if($#ARGV >= 6){
  $nr = $ARGV[6];
  ($nr >= 0) or die "\nERROR: invalid number of GCP points in range: $nr\n";
}

if($#ARGV >= 7){
  $naz = $ARGV[7];
  ($nr >= 0) or die "\nERROR: invalid number of GCP points in azimuth: $naz\n";
}

if($#ARGV >= 8){
  $gcp_win = $ARGV[8];
  ($gcp_win >= 1 && $gcp_win%2 != 0) or die "\nERROR: invalid GCP window size: $gcp_win\n";
}

if($#ARGV >= 9){
  $type = $ARGV[9];
  ($type >= 0 && $type <= 1 ) or die "\nERROR: differential interferogram type: $type\n";
}
print "$0  @ARGV\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"; 
-d $diff_dir or die "ERROR $0: differential interferogram directory does not exist\n";
-e $dem or die "ERROR $0: DEM does not exist: $dem\n";

$log = "$diff_dir/mk_base_2d.log";  
open(LOG,">$log") or die "ERROR $0: cannot open log file: $log\n";
LOG->autoflush;

$time = localtime;
print "start: $time  log file: $log\n";
print "RSLC_tab:           $rslc_tab\n";
print "itab:               $itab\n";
print "DIFF. directory:    $diff_dir\n";
print "DIFF. interf. type: $type\n";
print "GCP mask file:      $mask\n";

print LOG "\n$0  @ARGV\n";
print LOG "start: $time  log file: $log\n";
print LOG "RSLC_tab:           $rslc_tab\n";
print LOG "itab:               $itab\n";
print LOG "DIFF. directory:    $diff_dir\n";
print LOG "DIFF. interf. type: $type\n";
print LOG "GCP mask file:      $mask\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_par[$i] or die "ERROR $0: ISP image parameter file does not exist for slc: $slc2_par[$i]\n";
  printf LOG "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  $i++;
} 
print LOG"\n";
close LOG;

unless ($pbase =~ m/-/) {  	#generate blank pbase stack
  $btmp = "/tmp/btmp";
  execute("rm -f $pbase",$log) if -e $pbase;
  open(BTMP,">/tmp/btmp") or die "ERROR $0: cannot open temporary baseline file: $BTMP\n";
  print BTMP <<EOSTR; 
initial_baseline(TCN):          0.0000000        0.0000000        0.0000000   m   m   m
initial_baseline_rate:          0.0000000        0.0000000        0.0000000   m/s m/s m/s
precision_baseline(TCN):        0.0000000        0.0000000        0.0000000   m   m   m
precision_baseline_rate:        0.0000000        0.0000000        0.0000000   m/s m/s m/s
unwrap_phase_constant:          0.00000     radians
EOSTR
  close BTMP;
  LINE: while (<ITAB>) {	#read lines of ITAB 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
    $n > 0 or $rmax = $fields[2]; #determine maximum record number 
    if($rmax < $fields[2]){$rmax =  $fields[2];}
  }

  print "maximum record number in baseline stack: $rmax\n\n";
  seek ITAB, 0, 0;			#rewind the ITAB file to the beginning
  for($i = 1; $i <= $rmax; $i++){	#created blank filled baseline stack
    execute("base_par_pt $btmp $pbase $i 1",$log);
  }	
  execute("rm -f $btmp",$log);
}

$n = 0;
LINE: while (<ITAB>) {		#read lines of processing list file
  $n++;
  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); #do not process interferograms with the layer flag set to 0
  $time = localtime;
  printf "\n*** Diff. interferogram: %4d  $time ***\n",$n;

  $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   
  next LINE if ($s1 =~ m/$s2/);
  $fn = $diff_dir."/".$s1."_".$s2;

  $off =  $fn.".off";
  $base = $fn.".base";
  $diff_par = $fn.".diff_par";
  $gcp = $diff_dir."/gcp";
  $gcp_ph = $fn.".gcp_ph";
  $unw_int = $diff_dir."/unw_int";
  $sim = $fn.".sim";
  $sim_orb = $fn.".sim_orb";

  @slc1_freq = extract_param($slc1_par,"radar_frequency:");
  @slc2_freq = extract_param($slc2_par,"radar_frequency:");
  print "SLC-1 radar frequency: @slc1_freq\n";
  print "SLC-2 radar frequency: @slc2_freq\n";
  $f2flg = 0;
  if($slc1_freq[1] != $slc2_freq[1]){
    printf "SLC scenes have different radar carrier frequencies!\n"; 
    $f2flg = 1
  }
  
  print "SLC1: $slc1   slc1_par: $slc1_par\n";
  print "SLC2: $slc2   slc2_par: $slc2_par\n";
  print "DIFF: $fn   diff_par: $diff_par\n";
  print "DIFF. interferogram parameters: $off\n";
  print "GCP values: $gcp\n";
  print "GCP and phase data values: $gcp_ph\n";
  print "baseline file: $base\n";

#read interferogram size from *.off interferogram parameter file
  @width = extract_param($off, "interferogram_width:");
  @lines = extract_param($off, "interferogram_azimuth_lines:");
  @rstart = extract_param($off, "first_nonzero_range_pixel:");
  @rpix = extract_param($off, "number_of_nonzero_range_pixels:");
  print "interferogram width: $width[1]  lines: $lines[1]\n";
  
  if ($type == 0){
    $unw = $fn.".unw";
  }
  if ($type == 1){
    $unw = $fn.".adf.unw";
  }
  
  -e $unw or die "ERROR $0: unwrapped differential interferogram does not exist: $unw\n";
  if (-e $sim_orb) {
    print "simulated interferogram calculated using orbit state vectors: $sim_orb\n\n";
  }
  else {
    if (-e $sim) {
      print "simulated interferogram calculated using baseline model: $sim\n\n";
    }
    else {
      die "ERROR $0: simulated unwrapped interferogram does not exist: $sim\n";
    }
  }
  -e $off or die "ERROR $0: interferogram parameter file does not exist: $off\n";
  -e $base or die "ERROR $0: baseline file does not exist: $base\n";

  execute("extract_gcp $dem $off $gcp $nr $naz $mask", $log);

  if(-e $sim_orb){     #estimate only residual baseline
    execute("gcp_phase $unw $off $gcp $gcp_ph $gcp_win", $log);
  }
  else {
    execute("sub_phase $unw $sim $diff_par $unw_int 0 1",$log);   #add back simulated unwrapped interferogram back to the phase if baseline was used originally
    execute("gcp_phase $unw_int $off $gcp $gcp_ph $gcp_win", $log);
  }
  
  if($f2flg == 1){
    execute("base_ls $slc1_par $off $gcp_ph $base 0 1 1 1 1 5. $slc2_par", $log);
  }
  else{
    execute("base_ls $slc1_par $off $gcp_ph $base 0 1 1 1 1 5.", $log); 
  }

#copy baseline into pbase stack if given
  execute("base_par_pt $base $pbase $fields[2] 1",$log) unless $pbase =~ m/-/;
  print "deleting $unw_int\n";
  unlink $unw_int;
}

$time = localtime;
print "\nprocessing end: $time\n";
open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";
print LOG "\n*** mk_base_2d 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";
}