#!/usr/bin/env perl

use strict;

use FindBin;
use lib "$FindBin::Bin";

use BXHPerlUtils;

my $starttime = time();

# This script takes an input diffusion image and runs FSL's dtifit on it.
# Thanks to Chris Petty for the workflow.

my $progdiffdirs;
my $progbxhreorient;
my $progbxh2analyze;
my $proganalyze2bxh;
my $progbet;
my $progeddy;
my $progfslsplit;
my $progdtifit;

my %exechash =
  (
   'extractdiffdirs' => \$progdiffdirs,
   'bxhreorient' => \$progbxhreorient,
   'bxh2analyze' => \$progbxh2analyze,
   'analyze2bxh' => \$proganalyze2bxh,
   'bet' => \$progbet,
   'eddy_correct' => \$progeddy,
   'fslsplit' => \$progfslsplit,
   'dtifit' => \$progdtifit,
  );
foreach my $execname (keys %exechash) {
  ${$exechash{$execname}} = findexecutable($execname);
}
foreach my $execname (keys %exechash) {
  if (!defined(${$exechash{$execname}})) {
    print STDERR "Can't find required executable \"$execname\"!\n";
    exit -1;
  }
}

if (scalar(@ARGV) != 2) {
  print STDERR <<EOM;
ERROR: two arguments required: inputfile outputprefix"
inputfile is a diffusion image in .bxh or XCEDE format."
All output filenames will start with outputprefix.  outputprefix may"
contain directory names.  If outputprefix is an existing directory,"
files will be written to that directory."
EOM
  exit -1;
}

my $inputfile = shift;
my $outputprefix = shift;
if (-d $outputprefix || $outputprefix =~ m%/$%) {
  mkdir $outputprefix;
  $outputprefix = "${outputprefix}/";
} else {
  $outputprefix = "${outputprefix}_";
}

my $logfh = undef;
open($logfh, '>', "${outputprefix}LOG.txt") || die "Error opening '${outputprefix}LOG.txt' for writing: $!\n";

print $logfh "START: ", scalar(localtime($starttime)), "\n";

#### converts dti and orient as LAS ####

unlink glob("${outputprefix}tmp.*");
run_cmd([$logfh], $progbxhreorient, '--orientation=LAS', $inputfile, "${outputprefix}tmp.bxh");
unlink glob("${outputprefix}dti*");
run_cmd([$logfh], $progbxh2analyze, '-s', '--niigz', "${outputprefix}tmp.bxh", "${outputprefix}dti");
print $logfh "Removing '${outputprefix}tmp.*'\n";
unlink glob("${outputprefix}tmp.*");

#### grab correct gradient directions from *reoriented* volumes (since
#### FSL expects gradient directions to be in image space, doing it before
#### this won't work).
run_cmd([$logfh], $progdiffdirs, '--overwrite', '--fsl', "${outputprefix}dti.bxh", "${outputprefix}bvecs", "${outputprefix}bvals");

#### skull strip the dti
run_cmd([$logfh], $progbet, "${outputprefix}dti", "${outputprefix}dti_brain", '-F', '-f', '.25');

#### eddy current correction
run_cmd([$logfh], $progeddy, "${outputprefix}dti_brain", "${outputprefix}data", '0');
print $logfh "Removing '${outputprefix}dti_brain.nii.gz'\n";
unlink "${outputprefix}dti_brain.nii.gz";
print $logfh "Renaming '${outputprefix}dti_brain_mask.nii.gz' to '${outputprefix}nodif_brain_mask.nii.gz'\n";
rename "${outputprefix}dti_brain_mask.nii.gz", "${outputprefix}nodif_brain_mask.nii.gz";
run_cmd([$logfh], $progfslsplit, "${outputprefix}data", "${outputprefix}tmp", '-t');
print $logfh "Renaming '${outputprefix}tmp0000.nii.gz' to '${outputprefix}nodif_brain.nii.gz'\n";
rename "${outputprefix}tmp0000.nii.gz", "${outputprefix}nodif_brain.nii.gz";
print $logfh "Removing '${outputprefix}tmp.*'\n";
unlink glob("${outputprefix}tmp*");

##run dtifit
run_cmd([$logfh], $progdtifit, '-k', "${outputprefix}data", '-m', "${outputprefix}nodif_brain_mask", '-o', "${outputprefix}dti", '-r', "${outputprefix}bvecs", '-b', "${outputprefix}bvals");

for my $img ('FA', 'S0', 'MO', 'V1', 'V2', 'V3', 'L1', 'L2', 'L3', 'MD') {
  unlink "${outputprefix}dti_${img}.bxh";
  run_cmd([$logfh], $proganalyze2bxh, "${outputprefix}dti_${img}.nii.gz", "${outputprefix}dti_${img}.bxh");
}
for my $img ('nodif_brain_mask', 'nodif_brain', 'data') {
  unlink "${outputprefix}${img}.bxh";
  run_cmd([$logfh], $proganalyze2bxh, "${outputprefix}${img}.nii.gz", "${outputprefix}${img}.bxh");
}

my $endtime = time();
print $logfh "END: ", scalar(localtime($endtime)), "\n";

my $runtime = $endtime - $starttime;
my $runsecs = $runtime % 60;
my $runmins = int($runtime / 60) % 60;
my $runhours = int($runtime / (60 * 60));
print $logfh "Run time: ${runhours}h ${runmins}m ${runsecs}s\n";

close $logfh;

# $Log: In-line log eliminated on transition to SVN; use svn log instead. $
# Revision 1.2  2009/04/09 15:51:34  gadde
# Move extraction of diffusion directions to *after* reorient!
#
# Revision 1.1  2009/04/08 16:34:37  gadde
# Initial import.
#
