#!/usr/bin/perl

###############################################
# This script is used to Warp subject image
# to model image and for normalization and group
# analysis. The output is the deformation
# field from Model to Subject space, called
# "Field.1.img.Mdl2Subj", and a warped subject
# image(Warped.SubjSeg.img, already in
# model space)
#
# 12/05/2003
#
# Modified by lbloy 11/28/08
#   added a flag to NOT use fsl...
#   Default is to use fsl
###############################################

use Getopt::Std;

#get arguments list
$argc=@ARGV;

$useFSL=1;

if ($argc < 1){
   &showUsage ;
}

my $opt_string = 'hTFS:M:i:l:R:r:m:c:X:E:D:V:v:';

getopts("$opt_string")  || &showUsage ;

if($opt_h){
  &showUsage;
  }

#define argument list and get the arguments


# get Input files:
if($opt_S){
   $FileSegmented=$opt_S;
}else{
   print "\n\n Must input a subject DT image\n\n";
   &showUsage;
}

if($opt_M){
   $MdlFile =$opt_M;
}else{
   print "\n\n Must input a model DT image\n\n";
   &showUsage;
}


#iteration number
$iteration = 50;
if($opt_i){
  $iteration=$opt_i;
}

#iteration number for orientation features
$iter_orient = 0;
if($opt_l){
  $iter_orient=$opt_l;
}

# step ratio
$step_ration = 0.020;
if($opt_R){
  $step_ration=$opt_R;
}

# neighborhood search size, default: 12 for Low resolution,
# 10 for Middl resolution 8 for High resolutions
@search_size = (12,10,8);
if($opt_r){
  $sizes = $opt_r;
  @search_size = split(",", $sizes);
}

# option: -m
$m_degree=0.65;
if($opt_m){
   $m_degree=$opt_m;
}

# option: -c
# paramter used for the first parameter of option -c
$smooth = 0.5;
if($opt_c){
   $smooth=$opt_c;
}

#option: -D
$mdl_path = "." ;
if($opt_D){
   $mdl_path=$opt_D;
}
$mdl_path=$mdl_path."/";

# option: -X
@image_size=(256,256,70);
if($opt_X){
   $img_size=$opt_X;
   @image_size=split(",",$img_size);
}

# option: -E
# parameter used for the second paramter of option -E
$edge_smooth=1;
if($opt_E){
   $edge_smooth=$opt_E;
}

# option: -V
# voxel resolution for image model
@model_voxel=(0.9375,0.9375,2.5);
if($opt_V){
   $voxels = $opt_V;
   @model_voxel = split(",", $voxels);
}

# option: -v
# voxel resolution for image subject
@subject_voxel=(0.9375,0.9375,2.5);
if($opt_v){
   $voxels = $opt_v;
   @subject_voxel = split(",", $voxels);
}

# option: -F
# use FSL for affine registration
$useFSL = 1;
if($opt_F){
  print "turning off useFSl\n";
  $useFSL = 0;
}


print "\n\nWarp sujbect to Model by HAMMER\n\n";
print "Input information: \n";
print "\tinput subject segmented image is $FileSegmented\n";
print "\tinput subject original image is $FileSegmented\n";
print "\tinput Model image is $MdlFile\n";
print "\tdiretory where subject image stored: $mdl_path\n";
print "\timage size on XYZ dimension: @image_size\n";
print "\titeration number: $iteration\n";
print "\titeration number for orientation features: $iter_orient\n";
print "\tstep ratio: $step_ration\n";
print "\tneighborhood search size: ($search_size[0], $search_size[1], $search_size[2])\n";
print "\tsimilarity degree: $m_degree\n";
print "\tsmoothness factor for deformation field: $smooth\n";
print "\tTimes for smoothing deformation fields: $edge_smooth\n";
if( $useFSL == 1){
  print "\tUsing FSL Flirt for afine registration\n";
}
else{
  print "\tNot using FSL Flirt for afine registration\n";
}

print "\n\n\n";

# affine transformation: for upsample MdlFile
$MdlFile_withPath = $mdl_path."".$MdlFile;


open (IN, "$FileSegmented") || die "\nError open: $FileSegmented\n\n";
close(IN);
open (IN, "$MdlFile_withPath") || die "\nError open: $MdlFile_withPath\n\n";
close(IN);


$File_Rigid=$FileSegmented.".rig";
#print "XXXXXXXXXXX\n";
#print "/sbia/home/xwu/bin/RigidTransformImgA2ImgB $FileSegmented $MdlFile_withPath $File_Rigid -x$image_size,$image_size -T$FileSegmented\n";
#print "XXXXXXXXXXX\n";
#system("/sbia/home/xwu/bin/RigidTransformImgA2ImgB $FileSegmented $MdlFile_withPath $File_Rigid -x$image_size,$image_size -T$FileSegmented");

# affine transformation: for upsample MdlFile
$command="rigidlyregisterdt.sh";


print "\n\nLinear align two images\n";
if ($useFSL){
  $fcmd = "$command  $FileSegmented $MdlFile_withPath $File_Rigid 12 $image_size[0] $image_size[1] $image_size[2] $subject_voxel[0] $subject_voxel[1] $subject_voxel[2] $image_size[0] $image_size[1] $image_size[2] $model_voxel[0] $model_voxel[1] $model_voxel[2] Y";
}
else{
  $fcmd = "$command  $FileSegmented $MdlFile_withPath $File_Rigid 12 $image_size[0] $image_size[1] $image_size[2] $subject_voxel[0] $subject_voxel[1] $subject_voxel[2] $image_size[0] $image_size[1] $image_size[2] $model_voxel[0] $model_voxel[1] $model_voxel[2] N";
}
print "\n".$fcmd;

#this system call is retarded lbloy...
system("$fcmd") && die "Error: Linear align two images, check the FSL package  was installed properly and the input images\n";
#"$command  $FileSegmented $MdlFile_withPath $File_Rigid 12 $image_size[0] $image_size[1] $image_size[2] $subject_voxel[0] $subject_voxel[1] $subject_voxel[2] $image_size[0] $image_size[1] $image_size[2] $model_voxel[0] $model_voxel[1] $model_voxel[2]") && die "Error: Linear align two images, check the FSL package  was installed properly and the input images\n";
print "Done rigid transformation\n\n";

$HGH_FnlID=$iteration -1;

#remove previous results
if ( -e "*sample?.*gz") {
  system("rm *sample?.*gz");
}


#Low resoultion: 4
print "\n";
print "*******************************************************************\n";
print "Low resoultion:\n";
print "*******************************************************************\n";
print "Warping '$File_Rigid' in low-resolution\n";

#print "XXXXXXXXXXX\n";
#print "LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration  -r$search_size[0] -m$m_degree -s4 -c0.5 -B -E1 -p$smooth,$edge_smooth -X$image_size[0] \n";
print "TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration  -r$search_size[0] -m$m_degree -s4 -c0.5 -B -E1 -p$smooth,$edge_smooth -X$image_size[0] \n";
#print "XXXXXXXXXXX\n";
#system("LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[0] -m$m_degree -s4 -c0.5 -B -E1 -p$smooth,$edge_smooth -X$image_size[0]  > NoUse");
system("TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[0] -m$m_degree -s4 -c0.5 -B -E1 -p$smooth,$edge_smooth -X$image_size[0]  > NoUse");

#system("rm Warped.OriginalImg.sample4.*");
system("rm Res.img.sample4.*");
system("splitfile DeformationField.float.img NoUse 8 Field.4.img");

#Mid resoultion: 2
print "\n";
print "*******************************************************************\n";
print "Mid resoultion:\n";
print "*******************************************************************\n";
print "Warping '$File_Rigid' in mid-resolution\n";

#print "XXXXXXXXXXX\n";
#print "LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[1] -m$m_degree -s2 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] \n";
print "TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[1] -m$m_degree -s2 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] \n";
#print "XXXXXXXXXXX\n";
#system("LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[1] -m$m_degree -s2 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] > NoUse");
system("TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[1] -m$m_degree -s2 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] > NoUse");

#system("rm Warped.OriginalImg.sample2.*");
system("rm Res.img.sample2.*");
system("splitfile DeformationField.float.img NoUse 8 Field.2.img");

#High resoultion: 1
print "\n";
print "*******************************************************************\n";
print "high resoultion:\n";
print "*******************************************************************\n";
print "Warping '$File_Rigid' in high-resolution\n";

#print "XXXXXXXXXXX\n";
#print "LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration  -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] \n";
print "TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration  -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0] \n";
#print "XXXXXXXXXXX\n";
#system("LinuxTensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0]  > NoUse");
system("TensorStructRegistration_V2.0 $MdlFile $File_Rigid -R$step_ration  -D$mdl_path -i$iteration -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E1 -p$smooth,$edge_smooth -F -X$image_size[0]  > NoUse");

#****************************************************
# generate deformation field from model to subject
#****************************************************
print "splitfile DeformationField.float.img NoUse 8 Field.1.img\n";
system("splitfile DeformationField.float.img NoUse 8 Field.1.img");

print "transform3dvectorfield Field.1.img Field.1.img.trans GlobalTransformation.Affine -v$image_size[0]\n";
system("transform3dvectorfield Field.1.img Field.1.img.trans GlobalTransformation.Affine -v$image_size[0]");

system("mv Field.1.img.trans Field.1.img.Mdl2Subj");

#############begin the registration with orientation features################

for ($count=0; $count<$iter_orient; $count++)
{
  #*****************************************************
  # generate warped subject image to model space
  #*****************************************************

  #transform vector field for warpDT
  print "reversedeformationfield Field.1.img.Mdl2Subj Field.1.img.Subj2Mdl -X $image_size[0] -z $image_size[2] \n";
  system("reversedeformationfield Field.1.img.Mdl2Subj Field.1.img.Subj2Mdl -X $image_size[0] -z $image_size[2]");

  print "hammer2xdr Field.1.img.Mdl2Subj Field.1.img.Mdl2Subj.xdr $image_size[0] $image_size[1] $image_size[2] \n";
  system("hammer2xdr Field.1.img.Mdl2Subj Field.1.img.Mdl2Subj.xdr $image_size[0] $image_size[1] $image_size[2]");
  print "hammer2xdr Field.1.img.Subj2Mdl Field.1.img.Subj2Mdl.xdr $image_size[0] $image_size[1] $image_size[2] \n";
  system("hammer2xdr Field.1.img.Subj2Mdl Field.1.img.Subj2Mdl.xdr $image_size[0] $image_size[1] $image_size[2]");

  print "dtiPD 1 $image_size[0] $image_size[1] $image_size[2] $FileSegmented Subj.pd1.img 1 1 0 \n";
  system("dtiPD 1 $image_size[0] $image_size[1] $image_size[2] $FileSegmented Subj.pd1.img 1 1 0");

  #warpDT...
  print "warpDT.V5 $FileSegmented Field.1.img.Subj2Mdl.xdr Field.1.img.Mdl2Subj.xdr Subj.pd1.img -O Warped.Subj.dt -X $image_size[0] -Z $image_size[2] -R $subject_voxel[0],$subject_voxel[1],$subject_voxel[2] -N 0 \n";
  system("warpDT.V5 $FileSegmented Field.1.img.Subj2Mdl.xdr Field.1.img.Mdl2Subj.xdr Subj.pd1.img -O Warped.Subj.dt -X $image_size[0] -Z $image_size[2] -R $subject_voxel[0],$subject_voxel[1],$subject_voxel[2] -N 0");

  # Registration....
  #print "XXXXXXXXXXX\n";
#  print "LinuxTensorStructRegistration_V2.0 $MdlFile Warped.Subj.dt -R$step_ration  -D$mdl_path -i1 -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E3 -p$smooth,$edge_smooth -X$image_size[0] -o 0.4\n";
  print "TensorStructRegistration_V2.0 $MdlFile Warped.Subj.dt -R$step_ration  -D$mdl_path -i1 -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E3 -p$smooth,$edge_smooth -X$image_size[0] -o 0.4\n";
  #print "XXXXXXXXXXX\n";
#  system("LinuxTensorStructRegistration_V2.0 $MdlFile Warped.Subj.dt -R$step_ration  -D$mdl_path -i1 -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E3 -p$smooth,$edge_smooth -X$image_size[0] -o 0.4 > NoUse");
  system("TensorStructRegistration_V2.0 $MdlFile Warped.Subj.dt -R$step_ration  -D$mdl_path -i1 -r$search_size[2] -m$m_degree -s1 -c0.5 -B -E3 -p$smooth,$edge_smooth -X$image_size[0] -o 0.4 > NoUse");

  print "splitfile DeformationField.float.img NoUse 8 Field.1.img\n";
  system("splitfile DeformationField.float.img NoUse 8 Field.1.img");

  print "combinetwodeformationfields Field.1.img.Mdl2Subj Field.1.img Field.1.img.Mdl2Subj -d$image_size[0],$image_size[1],$image_size[2] -e$image_size[0],$image_size[1],$image_size[2] -S \n";
  system("combinetwodeformationfields Field.1.img.Mdl2Subj Field.1.img Field.1.img.Mdl2Subj -d$image_size[0],$image_size[1],$image_size[2] -e$image_size[0],$image_size[1],$image_size[2] -S");

}

#*****************************************************
# generate warped subject image to model space
#*****************************************************
#print "\nLinuxPerformDeformationOnImgUsingVectorField Field.1.img.Mdl2Subj $FileSegmented Warped.SubjSeg.img -v$image_size -s$image_size -I\n\n";
#system("LinuxPerformDeformationOnImgUsingVectorField Field.1.img.Mdl2Subj $FileSegmented Warped.SubjSeg.img -v$image_size -s$image_size -I");

#transform vector field for warpDT
print "reversedeformationfield Field.1.img.Mdl2Subj Field.1.img.Subj2Mdl -X $image_size[0] -z $image_size[2] \n";
system("reversedeformationfield Field.1.img.Mdl2Subj Field.1.img.Subj2Mdl -X $image_size[0] -z $image_size[2]");

print "hammer2xdr Field.1.img.Mdl2Subj Field.1.img.Mdl2Subj.xdr $image_size[0] $image_size[1] $image_size[2] \n";
system("hammer2xdr Field.1.img.Mdl2Subj Field.1.img.Mdl2Subj.xdr $image_size[0] $image_size[1] $image_size[2]");
print "hammer2xdr Field.1.img.Subj2Mdl Field.1.img.Subj2Mdl.xdr $image_size[0] $image_size[1] $image_size[2] \n";
system("hammer2xdr Field.1.img.Subj2Mdl Field.1.img.Subj2Mdl.xdr $image_size[0] $image_size[1] $image_size[2]");

print "dtiPD 1 $image_size[0] $image_size[1] $image_size[2] $FileSegmented Subj.pd1.img 1 1 0 \n";
system("dtiPD 1 $image_size[0] $image_size[1] $image_size[2] $FileSegmented Subj.pd1.img 1 1 0");

#warpDT...
print "warpDT.V5 $FileSegmented Field.1.img.Subj2Mdl.xdr Field.1.img.Mdl2Subj.xdr Subj.pd1.img -O Warped.Subj.dt -X $image_size[0] -Z $image_size[2] -R $subject_voxel[0],$subject_voxel[1],$subject_voxel[2] -N 0 \n";
system("warpDT.V5 $FileSegmented Field.1.img.Subj2Mdl.xdr Field.1.img.Mdl2Subj.xdr Subj.pd1.img -O Warped.Subj.dt -X $image_size[0] -Z $image_size[2] -R $subject_voxel[0],$subject_voxel[1],$subject_voxel[2] -N 0");



#***************************************************
# remove unwanted files
#****************************************************
#system("rm Warped.OriginalImg.sample1.*");
system("rm Res.img.sample1.*");

system("rm GlobalTransformation.Affine");
system("rm DeformationField.float.img");

system("rm MdlImg_*");
system("rm ObjImg_*");
system("rm Field.2.img");
system("rm Field.4.img");
#system("rm NoUse");
system("rm Field.1.img");


system("rm $File_Rigid");

system("rm -f optimized.d1");
system("rm -f Field.1.img.*.xdr");
system("rm -f Field.1.img.Subj2Mdl");
system("rm -f Subj.pd1.img");

system("rm -f NoUse");

# *************************************************
# function:  showUsage
# arguments: None
# *************************************************
sub showUsage{
  print "\n\n";
  print "This program is used to warp subject DT image to modle space. \n";
  print "The output is the deformation field from Model to Subject(Field.1.img.Mdl2Subj),\n";
  print "and a warped subject image (Warped.Subj.dt), already in model space.\n";
  print "\n";
  print "\n";
  print "Usage: $0 [options]\n";

  print "\t-h                          : display this (help) message\n";
  print "\t-M <string>                 : model image file name \n";
  print "\t-D <string>                 : directory path where model file is stored, default: current working directory\n";
  print "\t-S <string>                 : subject image file name, including path if not at current working direcotory\n";
  print "\t-X <int>,<int>,<int>        : image size on X,Y,Z dimension, default: <256,256,70>\n";
  print "\t-i <int>                    : iteration number, default: 50\n";
  print "\t-l <int>                    : iteration number for orientation features, default: 0 (without this feature)\n";
  print "\t                              keep the number small to save time for registration.\n";
  print "\t-R <float>                  : step ratio, default: 0.02\n";
  print "\t-r <int, int, int>          : neighborhood search size, different at each resolution.\n";
  print "\t                              input as < value_at_low, value_at_middle, value_at_high>\n";
  print "\t                              default: <12, 10, 8> mm/voxels\n";
  print "\t-V <float, float, float>    : voxel resolution for image model, default: <0.9375, 0.9375, 2.5>\n";
  print "\t-v <float, float, float>    : voxel resolution for subject, default: <0.9375, 0.9375, 2.5>\n";
  print "\t-m <float>                  : the similarity degree, default: 0.75\n";
  print "\t-c <float>                  : smoothness factor for deformation field, [0,1], 1: very smooth; default: 0.5\n";
  print "\t-E <int>                    : times for smoothing deformation fields, default: 1\n";
  print "\t-F                          : set if you would like to NOT use FSL flirt for affine Registations\n";
  print "                                if not set fsl flirt will be used for affine registration\n";

  print "\nNotes:\n";
  print "\t    This program only takes voxel-wise DT images as input. Each voxel contains 6 components [Dxx Dyy Dzz Dxy Dxz Dyz] in float format.\n\n";

  print "Example:\n";
  print "\tusing default value: \n";
  print "\t\t$0 -M modle.dt -S subject.dt\n\n";
  print "\twith options:\n";
  print "\t\t$0 -M modle.dt -S subject.dt -D/sbia/model -X256,256,70 -i40 -r10,8,6 -R0.07 -m0.8 -c0.3 -E2 \n\n";

    #print "References:\n";
  #print "\t1. Shen, D. and C. Davatzikos, \"Very high-resolution morphometry using mass-preserving\n";
  #print "\t   deformations and HAMMER elastic registration\". Neuroimage, 2003. 18(1): p. 28-41.\n";
  #print "\t2. Shen, D. and C. Davatzikos, \"HAMMER: Hierarchical attribute matching mechanism\n";
  #print "\t   for elastic registration\". IEEE Trans on Med. Imaging, 2002. 21(11): p. 1421-1439.\n\n";
    print "\n\n";

  exit(1);
}
