#! /usr/bin/env perl
#
# Surface Feature Extractor
#
# Oliver Lyttelton oliver@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);
my(@opt_table, $source, $target, $outxfm, @args, $tmpdir);

$me = &basename($0);

my $feature_choice = 0;
my $measurement_choice = 0;
my $extraction_transform = undef;
my $blurfwhm = undef;
my $alignment_transform = undef;
my $model_mid_surface = undef;
my $verbose = 0;
my $clobber = 0;
my $template_paint = undef;
my $template_labels = undef;
my $flip_transform = undef;
my $volume_feature = undef;
my $volume_feature_intersect = undef;

   
$Help = <<HELP;
|    $me extracts a surface feature set with options to 
|    choose the feature,  the type of measurement (roi vs vbc)
|    vary the extraction space (native or tal) the blurring kernel,
|    the registration space (tal or nl tal),
| 
| Problems or comments should be sent to: oliver\@bic.mni.mcgill.ca
HELP

$Usage = "Usage: $me [options] classified.mnc hemi_model_dir output_left_hemi.mnc output_right_hemi.mnc\n".
         "       $me -help to list options\n\n";

@opt_table = (
         
        ["-verbose","const","1", \$verbose,
        "be verbose" ],
        ["-clobber", "const","1", \$clobber,
        "clobber existing check features" ],  
   );

# Check arguments
&Getopt::Tabular::SetHelp($Help, $Usage);
&GetOptions (\@opt_table, \@ARGV) || exit 1;
die $Usage if($#ARGV != 3);
my $classify = shift(@ARGV);
my $hemi_model_dir= shift(@ARGV);
my $output_left_hemi = shift(@ARGV);
my $output_right_hemi = shift(@ARGV);

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



 my $wm_crop_left = "${tmpdir}/wm_crop_left.mnc";
 my $wm_crop_right  = "${tmpdir}/wm_crop_right.mnc";
 my $wm_defrag_left  = "${tmpdir}/wm_defrag_left.mnc";
 my $wm_defrag_right  = "${tmpdir}/wm_defrag_right.mnc";
 my $wm_add_back_left  = "${tmpdir}/wm_add_back_left.mnc";
 my $wm_add_back_right    = "${tmpdir}/wm_add_back_right.mnc";
 my $hemi_left_mask = "${hemi_model_dir}/hemi_left_mask.mnc"; 
my  $hemi_right_mask = "${hemi_model_dir}/hemi_right_mask.mnc";
my $slide_left_xfm =  "${hemi_model_dir}/slide_left.xfm";
my $slide_right_xfm =  "${hemi_model_dir}/slide_right.xfm";
                
#set up internal variables


 &do_cmd("minccalc","-clobber","-expr","out=(abs(A[0]-3)<0.45)&&A[1];",$classify,$hemi_left_mask,$wm_crop_left);

&do_cmd("minccalc","-clobber","-expr","out=(abs(A[0]-3)<0.45)&&A[1];",$classify,$hemi_right_mask,$wm_crop_right);
 	
&do_cmd("mincdefrag",$wm_crop_left,$wm_defrag_left,"1","27","100000");

&do_cmd("mincdefrag",$wm_crop_right,$wm_defrag_right,"1","27","100000");

&do_cmd("minccalc","-clobber","-expr","out=A[0]||(A[1]&&!A[2]);",$wm_defrag_left,$wm_crop_right,$wm_defrag_right,$wm_add_back_left);

&do_cmd("minccalc","-clobber","-expr","out=A[0]||(A[1]&&!A[2]);",$wm_defrag_right,$wm_crop_left,$wm_defrag_left,$wm_add_back_right);

&do_cmd("mincresample","-clobber","-like",$wm_add_back_left,"-transform",$slide_right_xfm,$wm_add_back_left,$output_left_hemi);

&do_cmd("mincresample","-clobber","-like",$wm_add_back_right,"-transform",$slide_left_xfm,$wm_add_back_right,$output_right_hemi);

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  if ($feature_choice ==0){
          #thickness
          &do_cmd("dump_rms",$source_white,$source_gray,$source_feature);
  }
 if ($feature_choice ==1){
          #surface_area                               
               &do_cmd("/data/aces/aces1/oliver/cvs/pipelines/civet/progs/cortex_area","-surface", $mid_surface,"-output",$source_feature);
  }
   if ($feature_choice ==2){
          #cortical_complexity
          &do_cmd("cortical_complexity",$mid_surface,$source_feature) ;
  }
     if ($feature_choice ==3){
          #cortical volume
           &do_cmd("cortical_volume",$source_white,$source_gray,$source_feature);
          
 }
 if ($feature_choice ==4){
          #volume_feature
           &do_cmd("volume_feature_on_surface",$source_white,$source_gray,$volume_feature,$volume_feature_intersect,$source_feature);
          
 }
   if ($feature_choice ==5){
          #callosal distance
          &do_cmd("callosal_distance",$source_white,$source_feature) ;
  }
 
 
 #ROI Analysis
 if ($measurement_choice==1) {
       #then we are doing ROI analysis 
       if ($alignment_transform){   
                 &do_cmd("surface-resample2","-clobber",$alignment_transform, $template_paint,$paint_on_target);
                 $template_paint =$paint_on_target ;
       }
          &do_cmd("roi_process","-clobber",$template_paint,$source_feature,$output);    
	}
   
 
 #VBC Cortical complexity, callosal distance, and thickness	
 if (($measurement_choice!=1)&&!($feature_choice==1 || $feature_choice==4)){ 
         if ($blurfwhm){
                &do_cmd( "diffuse", "-kernel", $blurfwhm, "-parametric", "0", $mid_surface, $source_feature, $blurred_feature);
                $source_feature = $blurred_feature; 
         }
          if ($alignment_transform){  
                  &do_cmd("surface-resample2","-clobber",$alignment_transform, $source_feature,$resampled_feature ); 
          }
 }
                                                            
#VBC Surface area and volume   
    if (($measurement_choice!=1) &&($feature_choice==1 || $feature_choice==4)) {               
           #create a node paint feature
           open(myFile,">${node_paint}");
           my $i = 1;
           while($i<40963){
           print myFile "${i}\n";
           $i=$i+1;
           }
           close myFile;
          &do_cmd("surface-resample2","-clobber",$alignment_transform, $node_paint,$paint_on_target ); 
          &do_cmd("zone_to_vertex_areas",$paint_on_target,$source_feature,$resampled_feature);
          $source_feature = $resampled_feature;
          if ($blurfwhm){
          &do_cmd( "diffuse", "-kernel", $blurfwhm, "-parametric", "0", $model_mid_surface, $source_feature, $blurred_feature);
         $source_feature = $blurred_feature;               
           }                            
}

#Optional Asymmetry stages                                                                      
  	if ($measurement_choice==2){
	    #perform_symmetry step
		&do_cmd("surface_resample2",$flip_transform,$source_feature, $flip_feature);
        #perform_symmetry step
        &do_cmd("surfmath","-clobber","-sub",$source_feature,$flip_feature,$diff_feature);
		$source_feature = $diff_feature;	
  }
  
  #copy to output
  &do_cmd("cp", $source_feature ,$output);
  


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