#!/usr/bin/env perl
#

# This step eliminates the skull and meninges in stereotaxic
# space using mincbet. For better performance by mincbet, use 
# uniform steps for (x,y,z). Also, resample in the same way 
# as the original t1 image (the mincbet mask is z-y-x so if 
# the original image is different, do a mincresample in order 
# to be able to apply the mask to its image).
#
# Copyright Alan C. Evans
# Professor of Neurology
# McGill University
#

use strict;
use warnings "all";

use MNI::Startup;
use MNI::FileUtilities qw(check_output_dirs);

# define input variables:

my $maskType=$ARGV[0];
my $t1_input=$ARGV[1];
my $t2_input=$ARGV[2];
my $pd_input=$ARGV[3];
my $headmask=$ARGV[4];
my $skull_mask=$ARGV[5];

# Directory for temporary files.

MNI::FileUtilities::check_output_dirs($TmpDir)
or exit 1;

# define internal variables:

my $bet_prefix = "${TmpDir}/bet";
my $bet_mask = "${bet_prefix}_mask.mnc";
my $bet_mask_temp = "${bet_prefix}_mask_temp.mnc";
my $input_resampled = "${TmpDir}/bet_resampled.mnc";

my $mincbet = "mincbet";

# In stereotaxic space, clear the background outside the registered
# head mask for t1. This should be only a mild attenuation of the 
# background such that we can distinguish it from CSF (usually for 
# 3-T scans). For t2/pd, CSF is inverted, so no need to do this.

my $mask_tmp = "${TmpDir}/nuc_mask.mnc";
&run( "mincresample", "-like", $t1_input, "-nearest_neighbour",
      $headmask, $mask_tmp );
my $tmp_input = "${TmpDir}/t1_nobg.mnc";
my $thresh = `mincstats -quiet -pctT 50 $t1_input`;
chomp( $thresh );
&run( "minccalc", "-clobber", "-quiet", "-expression",
      "if(A[1]<0.5&&A[0]<$thresh){A[0]/2.0}else{A[0]}", $t1_input,
      $mask_tmp, $tmp_input );
unlink( $mask_tmp );

# mask non-cortical tissues using mincbet

# run FSL bet on the classified t1 image (output is "${bet_prefix}_mask.mnc")
# (Use -h switch only for t1.)
my $bet_t1_mask = "${bet_prefix}_t1_mask.mnc";

&run( $mincbet, $tmp_input, "${bet_prefix}_t1", 
      "-n", "-m", "-h", "1.15", "-f", "0.50" );

&run( "mv", "-f", $bet_t1_mask, $bet_mask );

# Use t2 and pd to obtain a better mask.

if( $maskType eq "multispectral" ) {

  my $bet_expr = 'out=A[0]*A[1];';

  if( -e $t2_input ) {
    my $bet_t2_mask = "${bet_prefix}_t2_mask.mnc";
    &run( $mincbet, $t2_input, "${bet_prefix}_t2", 
          "-n", "-m", "-f", "0.50", "-r" );
    &run( "mincresample", "-clobber", "-unsigned", "-byte", "-nearest_neighbour",
          "-like", $bet_mask, $bet_t2_mask, $input_resampled );
    &run( "mv", "-f", $bet_mask, $bet_mask_temp );
    &run( "minccalc", "-clobber", "-expression", $bet_expr,
          $bet_mask_temp, $input_resampled, $bet_mask );
  }

  if( -e $pd_input ) {
    my $bet_pd_mask = "${bet_prefix}_pd_mask.mnc";
    &run( $mincbet, $pd_input, "${bet_prefix}_pd", 
          "-n", "-m", "-f", "0.50", "-r" );
    &run( "mincresample", "-clobber", "-unsigned", "-byte", "-nearest_neighbour",
          "-like", $bet_mask, $bet_pd_mask, $input_resampled );
    &run( "mv", "-f", $bet_mask, $bet_mask_temp );
    &run( "minccalc", "-clobber", "-expression", $bet_expr,
          $bet_mask_temp, $input_resampled, $bet_mask );
  }
}

# fill up holes when mincbet fails (rare case of inverted surface).

my $bet_filled = "${bet_prefix}_filled.mnc";
&run( "dilate_volume", $bet_mask, $bet_filled, 1, 6, 1 );
&run( "mincdefrag", $bet_filled, $bet_filled, 0, 6 );
&run( "dilate_volume", $bet_filled, $bet_filled, 0, 6, 2 );
&run( "minccalc", "-clobber", "-unsigned", "-byte", "-expression", 
      "if(A[0]+A[1]>0.5){out=1;}else{out=0;}",
      $bet_mask, $bet_filled, $bet_mask_temp );

# getting mincbet data in correct form -like $t1_input. This also
# applies when we have t2 and pd since the latter are coregistered
# to t1.

&run( "mincresample", "-clobber", "-unsigned", "-byte", "-nearest_neighbour",
      "-like", $t1_input, $bet_mask_temp, $skull_mask );

#Execute a system call.

sub run {
  print "@_\n";
  system(@_)==0 or die "Command @_ failed with status: $?";
}