#! /usr/bin/env perl
#
# Make an asymmetry map for cortical thickness.
#
# Claude Lepage - claude@bic.mni.mcgill.ca
#
# Copyright Alan C. Evans
# Professor of Neurology
# McGill University
#

use strict;
use warnings "all";
use Getopt::Tabular;
use File::Basename;
use File::Temp qw/ tempdir /;

my($Help, $Usage, $me);
$me = &basename($0);

$Help = <<HELP;
| Make an asymmetry map for cortical thickness. Use the resampled
| thickness maps for a meaningful comparison. The map is defined
| as
|           (L-R)/(0.5*(L+R))
| 
| Problems or comments should be sent to: claude\@bic.mni.mcgill.ca
HELP

$Usage = "Usage: $me [options] left_thickness.txt right_thickness.txt" .
         " diff_thickness.txt [diff_thickness_full.txt]\n".
         "       $me -help to list options\n\n";

my $verbose = 0;
my $clobber = 0;

my @opt_table = (
   ["-verbose", "boolean", 0, \$verbose,
      "be verbose" ],
   ["-clobber", "boolean", 0, \$clobber,
      "clobber existing check files" ],
   );

# Check arguments
&Getopt::Tabular::SetHelp($Help, $Usage);
&GetOptions (\@opt_table, \@ARGV) || exit 1;

die $Usage if( !( $#ARGV == 2 || $#ARGV == 3 ) );

my $left_thickness = shift(@ARGV);
my $right_thickness = shift(@ARGV);
my $diff_thickness = shift(@ARGV);
my $diff_thickness_full = shift(@ARGV);

# check for files
die "$me: Couldn't find input file: $left_thickness\n\n" if (!-e $left_thickness);
die "$me: Couldn't find input file: $right_thickness\n\n" if (!-e $right_thickness);
if(-e $diff_thickness && !$clobber){
   die "$me: $diff_thickness exists, -clobber to overwrite\n\n";
}

# make tmpdir
my $tmpdir = &tempdir( "$me-XXXXXXXX", TMPDIR => 1, CLEANUP => 1 );

# Left/right differences by hemisphere.
#
#     (L-R)/(0.5*(L+R))
#

&do_cmd( 'vertstats_math', '-old_style_file', '-sub', $left_thickness, $right_thickness,
         "${tmpdir}/l_minus_r.txt" );
&do_cmd( 'vertstats_math', '-old_style_file', '-add', $left_thickness, $right_thickness,
         "${tmpdir}/l_plus_r.txt" );
&do_cmd( 'vertstats_math', '-old_style_file', '-mult', "${tmpdir}/l_plus_r.txt", 
         '-const', '0.5', "${tmpdir}/l_plus_r.txt" );
# add very small coeff to avoid division by zero.
&do_cmd( 'vertstats_math', '-old_style_file', '-add', "${tmpdir}/l_plus_r.txt", 
         '-const', '1.e-8', "${tmpdir}/l_plus_r.txt" );
&do_cmd( 'vertstats_math', '-old_style_file', '-div', "${tmpdir}/l_minus_r.txt",
         "${tmpdir}/l_plus_r.txt", $diff_thickness );

unlink( "${tmpdir}/l_minus_r.txt" );
unlink( "${tmpdir}/l_plus_r.txt" );

# Left/right differences by hemisphere, on reconstructed full surface
# (left+right surfaces combined).

if( defined $diff_thickness_full ) {

  if(-e $diff_thickness_full && !$clobber){
     die "$me: $diff_thickness_full exists, -clobber to overwrite\n\n";
  }

  &do_cmd( 'vertstats_math', '-old_style_file', '-mult', $diff_thickness,
           '-const', '-1', "${tmpdir}/negative_diff.txt" );
  &do_cmd( 'cat', $diff_thickness, "${tmpdir}/negative_diff.txt", ' > ',
           $diff_thickness_full );
  unlink( "${tmpdir}/negative_diff.txt" );
}

sub do_cmd { 
   print STDOUT "@_\n" if $verbose;
   system("@_") == 0 or die;
}