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

$mode = 1;
$mcycle = 100.;
$scale = .7;	#display scale factor
$exp = 0.35;	#display exponent
#command line options parsing, the rest of the arguments are passed to @ARGV
GetOptions ('e=f' => \$exp, 's=f' => \$scale);

if (($#ARGV + 1) < 5){die <<EOS ;}
*** $0
*** Copyright 2013, Gamma Remote Sensing, v1.4 19-Jul-2013 clw ***
*** Calculate height maps from a stack of unwrapped differential interferograms ***

usage: $0 <RSLC_tab> <itab> <mli> <diff_dir> [type] [m/cycle]

    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)
    mli       (input) MLI image file with same range and azimuth looks as the interferogram (float)
    diff_dir  differential interferogram directory     
    type      unwrapped differential interferogram type:
                0: unfiltered differential interferogram (*.unw)
	        1: (default) adf filtered differential interferogram (*.adf.unw)
    m/cycle   meters per color cycle (default=100.)
    -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_dir  = $ARGV[3];

if($#ARGV >= 4){
  $mode = $ARGV[4];
  ($mode >= 0  && $mode <= 1) or die "\nERROR: invalid processing mode: $mode\n";
}

if($#ARGV >= 5){
  $mcycle = $ARGV[5];
  ($mcycle > 0) or die "\nERROR: invalid value for meters/cycle: $mcycle\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 reference image file does not exist: $rmli\n";
-d $diff_dir or die "ERROR $0: differential interferogram directory does not exist\n";

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

$time = localtime;
print "\n$0  @ARGV\n";
print "start: $time  log file: $log\n\n";
print "unwrapping mode:        $mode\n";
print "RSLC_tab:               $rslc_tab\n";
print "itab:                   $itab\n";
print "diff. directory:        $diff_dir\n\n";

print LOG "\n$0  @ARGV\n";
print LOG "start: $time  log file: $log\n\n";
print LOG "unwrapping mode:    $mode\n";
print LOG "RSLC_tab:           $rslc_tab\n";
print LOG "itab:               $itab\n";
print LOG "diff. directory:    $diff_dir\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];
  printf "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  printf LOG "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  $i++;
} 
print LOG"\n";
print "\n";
close 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 unwrap interferograms with the layer flag set to 0
  $slc1 = $slc[$fields[0]];
  $slc1_par = $slc_par[$fields[0]];
#  print " slc1: $slc1\n";
  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]];
#  print " slc2: $slc2\n";
  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;
  $diff_par = $fn.".diff_par";
  $off =  $fn.".off";
  $sim_unw = $fn.".sim_unw";
  $hgt = $fn.".hgt";
  $gr = $fn.".gr";
  $int_unw = $fn.".int_unw";
  $base = $fn.".base";

  $time = localtime;
  printf "\n*** diff. interferogram: %d  $time ***\n",$n;
  print "slc1: $slc1   slc1_par: $slc1_par\n";
  print "slc2: $slc2   slc2_par: $slc2_par\n";
  print "diff: $diff   diff_par: $diff_par\n";
  print "diff. interf. parameters: $off\n";
  print "reference SLC parameter file: $ref_par\n";
  print "height map: $hgt\n";
  print "raster height image: $hgt.ras  meters/color cycle: $mcycle\n";
  print "ground range map: $gr\n\n";
  
  @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";
  
  open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";
  print LOG "\n*** diff. interferogram:%3d  $time ***\n",$n;
  print LOG "slc1: $slc1   slc1_par: $slc1_par\n";
  print LOG "slc2: $slc2   slc2_par: $slc2_par\n";
  print LOG "diff: $diff   diff_par: $diff_par\n";
  print LOG "diff. interf. parameters: $off\n";
  print LOG "reference SLC parameter file: $ref_par\n";
  print LOG "height map: $hgt\n";
  print LOG "raster height image: $hgt.ras  meters/color cycle: $mcycle\n";
  print LOG "ground range map: $gr\n\n";
  print LOG "SLC-1 radar frequency: @slc1_freq\n";
  print LOG "SLC-2 radar frequency: @slc2_freq\n";
  close LOG;
  
  $f2flg = 0;
  if($slc1_freq[1] != $slc2_freq[1]){
    printf "SLC scenes have different radar carrier frequencies!\n"; 
    $f2flg = 1
  }
  
#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\n";
  
  if ($mode == 0){
    $diff = $fn.".diff";
    $unw = $fn.".unw";
  }
  if ($mode == 1){
    $diff = $fn.".adf.diff";
    $unw =  $fn.".adf.unw";
  }
      
#add unwrapped interferogram to simulated interferogram
  execute("lin_comb 2 $unw $sim_unw 0. 1. 1. $int_unw $width[1]",$log);
  if($f2flg == 1){
    execute("hgt_map $int_unw $slc1_par $off $base $hgt $gr 0 0 - $slc2_par", $log);
  }
  else {
    execute("hgt_map $int_unw $slc1_par $off $base $hgt $gr 0 ", $log);  
  }
  execute("rashgt $hgt $rmli $width[1] 1 1 0 1 1 $mcycle $scale $exp",$log);
}
$time = localtime;
print "\nprocessing end: $time\n";
open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";
print LOG "\nprocessing 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";
}