#! /usr/bin/env perl

############################################################################
############################################################################
### Post-processing module to the CIVET pipeline to perform basic
### quality control.
###                                       
###   Authors: Claude Lepage
###            August, 2009
###
###         Copyright Alan C. Evans
###         Professor of Neurology
###         McGill University
###
### For more information, have a look at our documentation at:
###    http://wiki.bic.mni.mcgill.ca/index.php/CIVET
###                                 
############################################################################
############################################################################

use strict;
use FindBin;
use Cwd qw( abs_path );
use Env qw( PATH );
use List::Util qw( min max );
use File::Temp qw/ tempdir /;
use File::Basename;

# All modules that will be used in either case are declared here
use Getopt::Tabular;
use MNI::Startup;
use MNI::PathUtilities qw(split_path);
use MNI::FileUtilities qw(check_output_dirs check_output_path);
use MNI::DataDir;

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

$PATH = "$FindBin::Bin/progs:${PATH}";

# Set interrupt handler (for cleaning of lock files)
$SIG{'INT'} = 'CLEANUP';
$SIG{'TERM'} = 'CLEANUP';

# Create tmp dir.
my $tmpdir = tempdir( CLEANUP => 1 );

my $version = "1.1.12";
my $versionDate= "February, 2013";
my $authors= "Claude Lepage";
my $usage = "\nUSAGE:\n$ProgramName -sourcedir <dir> -targetdir <dir> -prefix <prefix> [options] id1 id2 ... idn > <logfile> &\n 
ALTERNATIVE USAGE:\n
$ProgramName -sourcedir <dir> -targetdir <dir> -prefix <prefix> -id-file <idfile> [options]  > <logfile> &\n\n";

my $whatsnew = <<NEW;
    *************************************************************************
    New in $versionDate quarantine, version $version

    1 - Initial implementation

    More to come soon:

    *************************************************************************
NEW

my $help = <<HELP;

$ProgramName, version $version, released $versionDate.
    Released by $authors.

    Basic quality control for CIVET output.

$whatsnew
HELP
Getopt::Tabular::SetHelp($help, $usage);


############# The status report filename:
my $statusReportFile = "CIVET_status_report.csv";

############# Some directories that the user will need to specify
my $sourceDir = undef;
my $target = undef;
my $prefix = undef;
my $sourceSubDir = "noIdSubDir";
my $idTextFile = undef;

############# Options table
my @leftOverArgs;
my @argTbl = (
          ["File options", "section"],
          ["-sourcedir", "string", 1, \$sourceDir, "Directory containing the source files.", "<dir>"],
          ["-targetdir", "string", 1, \$target, "Directory where processed data will be placed.", "<dir>"],
          ["-prefix", "string", 1, \$prefix, "File prefix to be used in naming output files.", "<prefix>"],
          ["-id-subdir", "const", "IdSubDir", \$sourceSubDir, 
           "Indicate that the source directory contains sub-directories for each id"],
          ["-id-file", "string", 1, \$idTextFile, "A text file that contains all the subject id\'s (separated by space, tab, return or comma) that CIVET will run on.", "<file>"],

#         ["CIVET options", "section"],  
#         ["-thickness", "string", 2, \@thickness,
#          "compute cortical thickness and blur [tlink|tlaplace|tnear|tnormal] [kernel size in mm]"],

          );

GetOptions(\@argTbl, \@ARGV, \@leftOverArgs) or die "\n";


############# Basic usage
my @dsids;

if ($idTextFile) {
    open (IDTEXTFILE, "$idTextFile") or die ("Cannot open '$idTextFile': $!");
    # read the whole text file into one string
    my $idstext = "";
    while (my $idline = <IDTEXTFILE>) {
    $idstext .= $idline;
    }
    close (IDTEXTFILE) or die ("Cannot close '$idTextFile': $!");
    # split the string on whitespace (\s) or comma
    @dsids = split(/[\s,]+/, $idstext);
} 
else {
    @dsids = @leftOverArgs or die $usage;
}

unless ($prefix && $target && $sourceDir) {
    die "\n\n" .
        "*************************** ERROR ****************************: \n" .
        "     You must specify -prefix, -targetdir, and -sourcedir \n" .
        "**************************************************************\n\n\n";
}
$target =~ s#/+$##;      # remove trailing / at end of directory name, if any
$target = abs_path( $target );
$sourceDir = abs_path( $sourceDir );

############# Set no file buffering for stdout (buffer is printed every 1 line)
$| = 1;

############# Print the CIVET options list and related error messages
my $DATE = `date`;
chomp( $DATE );
my $UNAME = `uname -s -n -r`;
chomp( $UNAME );
print "\n\n* Pipeline started at $DATE on $UNAME \n";
print "\n* This is CIVET $version, $versionDate \n";
print "\n* CIVET Command line is:\n       $0 @ARGV \n";
print "\n* The source directory is: '$sourceDir' \n";
print "* The target directory is: '$target' \n";
print "* The prefix is: '$prefix' \n";

print "\n\n\n* Data-set Subject ID(s) is/are: '@dsids'\n\n\n";

my $qc_dir = "${target}/QC";
system("mkdir -p $qc_dir") if (! -d $qc_dir );

# Default parameters to create image object using MRI_Image

my $inputType = "t1only";
my $correctPVE = 0;
my $maskType = "t1only";
my $interpMethod = "trilinear";
my $headheight = 0;
my $MaskBloodVessels = 0;
my $N3Distance = 200;
my $N3Damping = "1.0e-06";
my $lsqtype = "-lsq9";
my $surface = "SURFACE";
my @thickness = ("tlink","30");
my $ResampleSurfaces = 1;
my $MeanCurvature = 0;
my $Area_fwhm = 40;
my $Volume_fwhm = 40;
my $CombineSurfaces = 0;
my $VBM = "noVBM";
my $VBM_fwhm = 8;
my $VBM_symmetry = "noSymmetry";
my $VBM_cerebellum = "Cerebellum";
my $animal = "noANIMAL";
my $Template = MNI::DataDir::dir("ICBM") . "icbm_template_1.00mm.mnc";

my %CIVETmodels = ( 'icbm152nl' => {
                      'RegLinDir'        => MNI::DataDir::dir("mni-models"),
                      'RegLinModel'      => "icbm_avg_152_t1_tal_nlin_symmetric_VI",
                      'RegNLDir'         => MNI::DataDir::dir("mni-models"),
                      'RegNLModel'       => "icbm_avg_152_t1_tal_nlin_symmetric_VI",
                      'TemplateDir'      => MNI::DataDir::dir("ICBM"),
                      'TemplateModel'    => "icbm_template",
                      'SurfaceMaskDir'   => MNI::DataDir::dir("mni-models"),
                      'SurfaceMask'      => "icbm_avg_152_t1_tal_nlin_symmetric_VI_mask.obj",
                      'SurfRegModelDir'  => "$FindBin::Bin/models",
                      'SurfRegModel'     => "surf_reg_model_left.obj",
                      'SurfRegDataTerm'  => "surf_reg_model_left.txt",
                      'SurfAtlas'        => "surface_atlas.txt",
                      'TagFileDir'       => MNI::DataDir::dir("classify"),
                      'TagFile'          => "n/a",
                      'AnimalAtlas'      => undef,
                      'AnimalAtlasDir'   => undef,
                      'AnimalNLRegDir'   => undef,
                      'AnimalNLRegModel' => undef
                    }
                  ); 

#########################################################
# Collect data for each subject.
#########################################################

my $good = "lightgreen";
my $warning = "orange";
my $danger = "red";
my $incomplete = "yellow";
my $good_level = 0;
my $warning_level = 1;
my $danger_level = 5;
my $incomplete_level = 3;

my $html = "<table cellspacing=\"2\" bgcolor=\"white\">\n";
$html .= "<tr bgcolor=\"\#cccccc\">" .
         "<th>No.</th>" .
         "<th>ID</th>" .
         "<th>x-step</th>" .
         "<th>y-step</th>" .
         "<th>z-step</th>" .
         "<th>StxMaskErr</th>" .
         "<th>CSFcls</th>" .
         "<th>GMcls</th>" .
         "<th>WMcls</th>" .
         "<th>GMCortical</th>" .
         "<th>WMsurf</th>" .
         "<th>GMsurf</th>" .
         "<th>SRLeft</th>" .
         "<th>SRRight</th>" .
         "<th>SSLeft</th>" .
         "<th>SSRight</th>" .
         "<th>MeanCTLeft</th>" .
         "<th>MeanCTRight</th>" .
         "<th>MeanWM-T1</th>" .
         "<th>StdDevWM-T1</th>" .
         "<th>MeanGM-T1</th>" .
         "<th>StdDevGM-T1</th>" .
         "<th>GIleft</th>" .
         "<th>GIright</th>" .
         "<th>GIfull</th>" .
         "</tr>\n";
$html .= "<tr> </tr>\n";

my $info = "# CIVET QC pipeline\n";

my $avg_stx_mask_error = 0.0;
my $stddev_stx_mask_error = 0.0;
my $count_stx_mask_error = 0;
my $avg_csf_cls = 0.0;
my $stddev_csf_cls = 0.0;
my $avg_gm_cls = 0.0;
my $stddev_gm_cls = 0.0;
my $avg_wm_cls = 0.0;
my $stddev_wm_cls = 0.0;
my $count_cls = 0.0;
my $avg_gm_cortical = 0.0;
my $stddev_gm_cortical = 0.0;
my $count_gm_cortical = 0.0;
my $avg_wm_surf = 0.0;
my $stddev_wm_surf = 0.0;
my $avg_gm_surf = 0.0;
my $stddev_gm_surf = 0.0;
my $count_gm_wm_surf = 0.0;
my $avg_surfreg_left = 0.0;
my $stddev_surfreg_left = 0.0;
my $count_surfreg_left = 0.0;
my $avg_surfreg_right = 0.0;
my $stddev_surfreg_right = 0.0;
my $count_surfreg_right = 0.0;
my $avg_surfsurf_left = 0.0;
my $stddev_surfsurf_left = 0.0;
my $count_surfsurf_left = 0.0;
my $avg_surfsurf_right = 0.0;
my $stddev_surfsurf_right = 0.0;
my $count_surfsurf_right = 0.0;
my $avg_meanct = 0.0;
my $stddev_meanct = 0.0;
my $count_meanct = 0.0;

my $count = 0;
foreach my $dsid (@dsids) {

    $count++;
    my $level = 0;

    print "Processing subject $dsid...\n";

    ##### Create image object. #####

    # depending on the two options your files can be in on target
    # or there can be subdirs for every subject
    my $Source_Base;
    if ($sourceSubDir eq "noIdSubDir") {
      $Source_Base = "${sourceDir}/";
    } else {
      $Source_Base = "${sourceDir}/${dsid}/";
    }

    my $image = MRI_Image->new( $version,
                                $Source_Base,
                                $target,
                                $prefix,
                                $dsid,
                                $inputType,
                                $correctPVE,
                                $maskType,
                                $interpMethod,
                                $headheight,
                                $MaskBloodVessels,
                                $N3Distance,
                                $N3Damping,
                                $lsqtype,
                                $surface,
                                \@thickness,
                                $ResampleSurfaces,
                                $MeanCurvature,
                                $Area_fwhm,
                                $Volume_fwhm,
                                $CombineSurfaces,
                                $VBM,
                                $VBM_fwhm,
                                $VBM_symmetry,
                                $VBM_cerebellum,
                                $animal,
                                $Template,
                                \$CIVETmodels{icbm152nl} );

    # Column 1: subject id.
    my $infoline = "$dsid ";
    my $htmlline = "";

    # Columns 2,3,4 : voxel sizes in native space
    if( -e $image->{t1}{native} ) {
      my $dx = `mincinfo -attvalue xspace:step $image->{t1}{native}`;
      my $dy = `mincinfo -attvalue yspace:step $image->{t1}{native}`;
      my $dz = `mincinfo -attvalue zspace:step $image->{t1}{native}`;
      chomp( $dx ); chomp( $dy ); chomp( $dz );
      $dx = abs($dx); $dy = abs($dy); $dz = abs($dz);
      $infoline .= sprintf( "%5.2f %5.2f %5.2f ", abs($dx), abs($dy), abs($dz) );
      
      if( abs($dx) >= 2.5 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( abs($dx) >= 1.75 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= sprintf( "%5.2f", abs($dx) ) . "</td>";
      if( abs($dy) >= 2.5 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( abs($dy) >= 1.75 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= sprintf( "%5.2f", abs($dy) ) . "</td>";
      if( abs($dz) >= 2.5 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( abs($dx) >= 1.75 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= sprintf( "%5.2f", abs($dz) ) . "</td> ";
    } else {
      $infoline .= "  0  0  0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> <td bgcolor=\"$incomplete\">0</td> <td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Column 5: Error in stereotaxic brain mask.
    my $skull_mask = $image->{skull_mask_tal};
    if( -e $skull_mask ) {

      my $model_mask = "$CIVETmodels{icbm152nl}{RegLinDir}/$CIVETmodels{icbm152nl}{RegLinModel}_mask.mnc";

      my $skull_mask_rsl = "${tmpdir}/skull_mask_rsl.mnc";
      `mincresample -clobber -quiet -like $model_mask -nearest_neighbour $skull_mask $skull_mask_rsl`;

      my $diff = "${tmpdir}/diff_mask_stx.mnc";
      `minccalc -clobber -quiet -expression 'A[1]!=A[0]' $model_mask $skull_mask_rsl $diff`;

      my $diff_volume;
      chomp( $diff_volume = `mincstats -quiet -sum $diff` );
      my $mask2_volume;
      chomp( $mask2_volume = `mincstats -quiet -sum $model_mask` );
      my $error = 100.0 * $diff_volume / $mask2_volume;

      unlink( $skull_mask_rsl );
      unlink( $diff );

      if( $error < 20.0 ) {   # exclude errors from avg and stddev
        $avg_stx_mask_error += $error;
        $stddev_stx_mask_error += $error*$error;
        $count_stx_mask_error++; 
      }
      $error = sprintf( "%5.2f", $error );
      $infoline .= "$error ";
      if( $error >= 20.0 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( $error >= 10.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $error . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 6,7,8: percentages of CSF, GM, WM tissue in brain.
    my $classify_qc = $image->{classify_qc};
    if( -e $classify_qc ) {
      open (IDFILE, $classify_qc);
      my $idline = <IDFILE>;
      $idline =~ /classified image CSF (.*)\%  GM (.*)\%  WM (.*)\%/;
      my $csf = $1;
      my $gm = $2;
      my $wm = $3;
      $infoline .= "$csf $gm $wm ";
      close (IDFILE);
      $avg_csf_cls += $csf;
      $stddev_csf_cls += $csf*$csf;
      $avg_gm_cls += $gm;
      $stddev_gm_cls += $gm*$gm;
      $avg_wm_cls += $wm;
      $stddev_wm_cls += $wm*$wm;
      $count_cls++;
      if( $csf >= 50.0 || $csf <= 5.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $csf . "</td> ";
      if( $gm >= 80.0 || $gm <= 15.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $gm . "</td> ";
      if( $wm >= 80.0 || $wm <= 15.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $wm . "</td> ";
    } else {
      $infoline .= " 0   0   0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> <td bgcolor=\"$incomplete\">0</td> <td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Column 9: Volume of gray matter inside the gray surface (cortical gray
    #           without sub-cortical).
    my $cerebral_volume_qc = $image->{cerebral_volume};
    if( -e $cerebral_volume_qc ) {
      open (IDFILE, $cerebral_volume_qc);
      my $idline = <IDFILE>;
      my $idline = <IDFILE>;  # need second line for cortical gray
      close (IDFILE);
      $idline =~ s/^\s+//;
      my ($lbl, $gm) = split(/\s+/, $idline); chomp( $gm );
      $gm /= 1000.0;    # convert from mm^3 to cc
      $infoline .= "$gm ";
      $htmlline .= "<td>" . $gm . "</td>";
      $count_gm_cortical++;
      $avg_gm_cortical += $gm;
      $stddev_gm_cortical += $gm*$gm;
    } else {
      $infoline .= "0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 10,11: Error in wm surface and gm surface.
    my $surface_qc = $image->{surface_qc};
    if( -e $surface_qc ) {
      open (IDFILE, $surface_qc);
      my $idline = <IDFILE>;
      $idline =~ /white surface \((.*)\%\), gray surface \((.*)\%\)/;
      my $wm = $1;
      my $gm = $2;
      $infoline .= "$wm $gm ";
      close (IDFILE);
      $avg_gm_surf += $gm;
      $stddev_gm_surf += $gm*$gm;
      $avg_wm_surf += $wm;
      $stddev_wm_surf += $wm*$wm;
      $count_gm_wm_surf++;
      if( $wm >= 20.0 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( $wm >= 10.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td> ";
      }
      $htmlline .= $wm . "</td>";
      if( $gm >= 20.0 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( $gm >= 10.0 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $gm . "</td> ";

    } else {
      $infoline .= " 100 100 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">100</td> <td bgcolor=\"$incomplete\">100</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 12: number of self-intersections on resampled left mid surface
    my $mid_rsl_left = $image->{mid_surface_rsl}{left};
    if( -e $mid_rsl_left ) {
      my @ret = `check_self_intersect $mid_rsl_left`;
      my $left_inter = 0;
      if( $ret[1] =~ m/self_intersecting/ ) {      # little legacy typo
        $ret[1] =~ /Number of self_intersecting triangles = (\d+)/;
        $left_inter = $1;
      } else {
        $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
        $left_inter = $1;
      }
      $infoline .= "$left_inter ";
      if( $left_inter < 500 ) {
        $avg_surfreg_left += $left_inter;
        $stddev_surfreg_left += $left_inter*$left_inter;
        $count_surfreg_left++;
      }
      if( $left_inter >= 500 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( $left_inter >= 250 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $left_inter . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 13: number of self-intersections on resampled right mid surface
    my $mid_rsl_right = $image->{mid_surface_rsl}{right};
    if( -e $mid_rsl_right ) {
      my @ret = `check_self_intersect $mid_rsl_right`;
      my $right_inter = 0;
      if( $ret[1] =~ m/self_intersecting/ ) {      # little legacy typo
        $ret[1] =~ /Number of self_intersecting triangles = (\d+)/;
        $right_inter = $1;
      } else {
        $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
        $right_inter = $1;
      }
      $infoline .= "$right_inter ";
      if( $right_inter < 500 ) {
        $avg_surfreg_right += $right_inter;
        $stddev_surfreg_right += $right_inter*$right_inter;
        $count_surfreg_right++;
      }
      if( $right_inter >= 500 ) {
        $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
      } elsif( $right_inter >= 250 ) {
        $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
      } else {
        $htmlline .= "<td>";
      }
      $htmlline .= $right_inter . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 14,15: number of surface-surface intersections between left/right wm-gm surfaces
    foreach my $side ("left","right") {
      my $gm = $image->{gray}{$side};
      my $wm = $image->{cal_white}{$side};
      if( -e $gm && -e $wm ) {
        my $tmpobj = "${tmpdir}/surfsurf_${side}.obj";
        my $tmptxt = "${tmpdir}/surfsurf_${side}.txt";
        `objconcat $wm $gm none none $tmpobj none`;
        `check_self_intersect $tmpobj $tmptxt`;
        my @ret = split( / /, `wc -l $tmptxt` );
        my $npoints = $ret[0]; chomp( $npoints ); $npoints /= 2;
        `head -${npoints} $tmptxt > $tmpobj`;
        `mv -f $tmpobj $tmptxt`;
        `vertstats_math -const2 -0.001 0.001 -seg -old_style_file $tmptxt $tmptxt`;
        @ret = `vertstats_stats $tmptxt |grep Sum`;
        unlink( $tmptxt );
        $ret[0] =~ / Sum: (.*)/;
        my $num_inter = $1;
        $infoline .= "$num_inter ";
        if( $num_inter < 500 ) {
          if( $side eq "left" ) {
            $avg_surfsurf_left += $num_inter;
            $stddev_surfsurf_left += $num_inter*$num_inter;
            $count_surfsurf_left++;
          } else {
            $avg_surfsurf_right += $num_inter;
            $stddev_surfsurf_right += $num_inter*$num_inter;
            $count_surfsurf_right++;
          }
        }
        if( $num_inter >= 500 ) {
          $htmlline .= "<td bgcolor=\"$danger\">"; $level = max($level,$danger_level);
        } elsif( $num_inter >= 250 ) {
          $htmlline .= "<td bgcolor=\"$warning\">"; $level = max($level,$warning_level);
        } else {
          $htmlline .= "<td>";
        }
        $htmlline .= $num_inter . "</td> ";
      } else {
        $infoline .= " 0 ";
        $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
        $level = max($level,$incomplete_level);
      }
    }

    # Colums 16,17: mean cortical thickness for left,right hemispheres
    foreach my $side ("left","right") {
      my $meanct_file = $image->{lobe_thickness}{$side};
      if( -e $meanct_file ) {
        my $idline = `grep Total $meanct_file`;
        chomp( $idline );
        $idline =~ / Total  (.*)/;
        my $meanct = $1;
        $infoline .= "$meanct ";
        $htmlline .= "<td>" . $meanct . "</td> ";
        $avg_meanct += $meanct;
        $stddev_meanct += $meanct*$meanct;
        $count_meanct++;
      } else {
        $infoline .= " 0 ";
        $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
        $level = max($level,$incomplete_level);
      }
    }

    # Columns 18,19,20,21: T1 mean and stddev for WM and GM.
    my $t1 = $image->{t1}{final};
    # my $cls = $image->{cls_correct};
    my $cls = "$image->{pve_prefix}\_disc.mnc";
    my $mask = $image->{skull_mask_tal};
    if( -e $t1 && -e $cls && -e $mask ) {
      my $cls_masked = "${tmpdir}/cls_masked.mnc";
      `minccalc -quiet -clobber -expression 'if(A[1]>0.5){A[0]}else{0}' $cls $mask $cls_masked`;
      my $WMmean = `mincstats -quiet -mask $cls_masked -mask_binvalue 3 -mean $t1`;
      my $WMsdev = `mincstats -quiet -mask $cls_masked -mask_binvalue 3 -stddev $t1`;
      my $GMmean = `mincstats -quiet -mask $cls_masked -mask_binvalue 2 -mean $t1`;
      my $GMsdev = `mincstats -quiet -mask $cls_masked -mask_binvalue 2 -stddev $t1`;
      chomp( $WMmean );
      chomp( $WMsdev );
      chomp( $GMmean );
      chomp( $GMsdev );
      unlink( $cls_masked );
      $infoline .= "$WMmean $WMsdev $GMmean $GMsdev ";
      $htmlline .= "<td>" . $WMmean . "</td> ";
      $htmlline .= "<td>" . $WMsdev . "</td> ";
      $htmlline .= "<td>" . $GMmean . "</td> ";
      $htmlline .= "<td>" . $GMsdev . "</td> ";
    } else {
      $infoline .= " 0 0 0 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 22: gyrification index for left hemisphere
    my $gi_left = $image->{gyrification_index}{left};
    if( -e $gi_left ) {
      open (IDFILE, $gi_left);
      my $idline = <IDFILE>;
      $idline =~ /gyrification index gray: (.*)/;
      my $gi = $1;
      $infoline .= "$gi ";
      close (IDFILE);
      $htmlline .= "<td>" . $gi . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 23: gyrification index for right hemisphere
    my $gi_right = $image->{gyrification_index}{right};
    if( -e $gi_right ) {
      open (IDFILE, $gi_right);
      my $idline = <IDFILE>;
      $idline =~ /gyrification index gray: (.*)/;
      my $gi = $1;
      $infoline .= "$gi ";
      close (IDFILE);
      $htmlline .= "<td>" . $gi . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    # Columns 24: gyrification index for left+right hemispheres combined
    my $gi_full = $image->{gyrification_index}{full};
    if( -e $gi_full ) {
      open (IDFILE, $gi_full);
      my $idline = <IDFILE>;
      $idline =~ /gyrification index gray: (.*)/;
      my $gi = $1;
      $infoline .= "$gi ";
      close (IDFILE);
      $htmlline .= "<td>" . $gi . "</td> ";
    } else {
      $infoline .= " 0 ";
      $htmlline .= "<td bgcolor=\"$incomplete\">0</td> "; 
      $level = max($level,$incomplete_level);
    }

    $infoline .= $dsid . "/thickness/" . basename( $image->{rms_rsl}{left} ) . " ";
    $infoline .= $dsid . "/thickness/" . basename( $image->{rms_rsl}{right} ) . " ";
    $infoline .= $dsid . "/surfaces/" . basename( $image->{surface_area_rsl}{left} ) . " ";
    $infoline .= $dsid . "/surfaces/" . basename( $image->{surface_area_rsl}{right} ) . " ";
    $infoline .= $dsid . "/surfaces/" . basename( $image->{surface_volume_rsl}{left} ) . " ";
    $infoline .= $dsid . "/surfaces/" . basename( $image->{surface_volume_rsl}{right} ) . " ";
    $infoline .= $dsid . "/surfaces/" . basename( $image->{rms_rsl}{asym_hemi} ) . " ";

    $infoline .= "\n";
    $info .= $infoline;

    # Complete html line for global error good.

    my $idcolourcode;
    if( $level == $warning_level ) {
      $idcolourcode = $warning;
    } elsif( $level == $danger_level ) {
      $idcolourcode = $danger;
    } elsif( $level == $incomplete_level ) {
      $idcolourcode = $incomplete;
    } else {
      $idcolourcode = $good;
    }

    $htmlline = "<th bgcolor=\"${idcolourcode}\">" . 
                "<a target=\"_blank\" href=\"${dsid}.html\">$dsid</a></th>" . $htmlline;
    $htmlline = "<tr bgcolor=\"#dddddd\">" . "<td>" . sprintf( "%4d", $count ) . "</td>" .
                $htmlline . "</tr>\n";
    $html .= $htmlline;

    open FILE, ">${qc_dir}/${dsid}.html";
    print FILE "<table>\n";
    print FILE "<tr><td><image src=\"../${dsid}/verify/${prefix}_${dsid}_verify.png\"></td></tr>\n";
    print FILE "<tr><td><image src=\"../${dsid}/verify/${prefix}_${dsid}_clasp.png\"></td></tr>\n";
    print FILE "<tr><td><image src=\"../${dsid}/verify/${prefix}_${dsid}_atlas.png\"></td></tr>\n";
    print FILE "<tr><td><image src=\"../${dsid}/verify/${prefix}_${dsid}_surfsurf.png\"></td></tr>\n";
    print FILE "<tr><td><image src=\"../${dsid}/verify/${prefix}_${dsid}_laplace.png\"></td></tr>\n";
    print FILE "</table>\n";
    close FILE;

}

# Line of averages and stddev on some fields.

$avg_stx_mask_error /= $count_stx_mask_error;
$avg_csf_cls /= $count_cls;
$avg_gm_cls /= $count_cls;
$avg_wm_cls /= $count_cls;
$avg_gm_cortical /= $count_gm_cortical;
$avg_gm_surf /= $count_gm_wm_surf if( $count_gm_wm_surf );
$avg_wm_surf /= $count_gm_wm_surf if( $count_gm_wm_surf );
$avg_surfreg_left /= $count_surfreg_left if( $count_surfreg_left );
$avg_surfreg_right /= $count_surfreg_right if( $count_surfreg_right );
$avg_surfsurf_left /= $count_surfsurf_left if( $count_surfsurf_left );
$avg_surfsurf_right /= $count_surfsurf_right if( $count_surfsurf_right );
$avg_meanct /= $count_meanct if( $count_meanct );

$html .= "<tr> </tr> ";
$html .= "<tr bgcolor=\"#dddddd\">" .
         "<th>Mean</th> <td></td> <td></td> <td></td> <td></td> " .
         "<td>" . sprintf( "%8.3f", $avg_stx_mask_error ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_csf_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_gm_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_wm_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_gm_cortical ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_wm_surf ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_gm_surf ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $avg_surfreg_left ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $avg_surfreg_right ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $avg_surfsurf_left ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $avg_surfsurf_right ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $avg_meanct ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $avg_meanct ) . "</td>" .
         "<td></td> <td></td> <td></td> <td></td> " .
         "<td></td> <td></td> <td></td>" .
         "</tr>\n";

$stddev_stx_mask_error = sqrt( $stddev_stx_mask_error/$count_stx_mask_error -
                               $avg_stx_mask_error*$avg_stx_mask_error );
$stddev_csf_cls = sqrt( $stddev_csf_cls/$count_cls - $avg_csf_cls*$avg_csf_cls );
$stddev_gm_cls = sqrt( $stddev_gm_cls/$count_cls - $avg_gm_cls*$avg_gm_cls );
$stddev_wm_cls = sqrt( $stddev_wm_cls/$count_cls - $avg_wm_cls*$avg_wm_cls );
$stddev_gm_cortical = sqrt( $stddev_gm_cortical/$count_gm_cortical - 
                            $avg_gm_cortical * $avg_gm_cortical );
$stddev_gm_surf = sqrt( $stddev_gm_surf/$count_gm_wm_surf - 
                        $avg_gm_surf*$avg_gm_surf ) if( $count_gm_wm_surf );
$stddev_wm_surf = sqrt( $stddev_wm_surf/$count_gm_wm_surf - 
                        $avg_wm_surf*$avg_wm_surf ) if( $count_gm_wm_surf );
$stddev_surfreg_left = sqrt( $stddev_surfreg_left/$count_surfreg_left -
                             $avg_surfreg_left*$avg_surfreg_left ) if( $count_surfreg_left );
$stddev_surfreg_right = sqrt( $stddev_surfreg_right/$count_surfreg_right -
                              $avg_surfreg_right*$avg_surfreg_right ) if( $count_surfreg_right );
$stddev_surfsurf_left = sqrt( $stddev_surfsurf_left/$count_surfsurf_left -
                              $avg_surfsurf_left*$avg_surfsurf_left ) if( $count_surfsurf_left );
$stddev_surfsurf_right = sqrt( $stddev_surfsurf_right/$count_surfsurf_right -
                               $avg_surfsurf_right*$avg_surfsurf_right ) if( $count_surfsurf_right );
$stddev_meanct = sqrt( $stddev_meanct/$count_meanct -
                       $avg_meanct*$avg_meanct ) if( $count_meanct );

$html .= "<tr bgcolor=\"#dddddd\">" .
         "<th>StdDev</th> <td></td> <td></td> <td></td> <td></td> " .
         "<td>" . sprintf( "%8.3f", $stddev_stx_mask_error ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_csf_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_gm_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_wm_cls ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_gm_cortical ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_wm_surf ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_gm_surf ) . "</td> " .
         "<td>" . sprintf( "%5.2f", $stddev_surfreg_left ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $stddev_surfreg_right ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $stddev_surfsurf_left ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $stddev_surfsurf_right ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $stddev_meanct ) . "</td>" .
         "<td>" . sprintf( "%5.2f", $stddev_meanct ) . "</td>" .
         "<td></td> <td></td> <td></td> <td></td>" .
         "<td></td> <td></td> <td></td>" .
         "</tr>\n";
$html .= "</table>";

# Do some scatter plots of the various fields.

open FILE, ">${qc_dir}/civet_${prefix}.glm";
$info =~ s/ +/ /g;
print FILE $info;
close FILE;

# ----------------------------------------------------------------------------------------------
open FILE, ">${qc_dir}/civet_${prefix}.gnu";
my $figure = 0;

$figure++;
print FILE "# Stereotaxic space brain masking error\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Stereotaxic space brain mask error\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Error %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/StxMaskErr.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 5 w po, $avg_stx_mask_error w li 2, " .
           sprintf("%f",$avg_stx_mask_error-$stddev_stx_mask_error) . " w li 3, " .
           sprintf("%f",$avg_stx_mask_error+$stddev_stx_mask_error) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"StxMaskErr.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Error between stereotaxic space brain mask " .
         "and the brain mask of the stereotaxic model. </td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Percentage of the classified CSF in the brain (masked)\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Percentage of the classified CSF in the brain (masked)\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"CSF %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/CSFclsPct.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 6 w po, $avg_csf_cls w li 2, " .
           sprintf("%f",$avg_csf_cls-$stddev_csf_cls) . " w li 3, " .
           sprintf("%f",$avg_csf_cls+$stddev_csf_cls) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"CSFclsPct.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Percentage of the classified CSF in the brain (masked)." .
         "</td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Percentage of the classified GM in the brain (masked)\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Percentage of the classified GM in the brain (masked)\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"GM %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GMclsPct.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 7 w po, $avg_gm_cls w li 2, " .
           sprintf("%f",$avg_gm_cls-$stddev_gm_cls) . " w li 3, " .
           sprintf("%f",$avg_gm_cls+$stddev_gm_cls) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GMclsPct.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Percentage of the classified GM in the brain (masked)." .
         "</td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Percentage of the classified WM in the brain (masked)\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Percentage of the classified WM in the brain (masked)\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"WM %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/WMclsPct.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 8 w po, $avg_wm_cls w li 2, " .
           sprintf("%f",$avg_wm_cls-$stddev_wm_cls) . " w li 3, " .
           sprintf("%f",$avg_wm_cls+$stddev_wm_cls) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"WMclsPct.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Percentage of the classified WM in the brain (masked)." .
         "</td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------

$figure++;
print FILE "# Volume of total cortical gray matter (without sub-cortical gray)\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Volume of total cortical gray matter (without sub-cortical gray)\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"GM cc\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GMCortical.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 9 w po, $avg_gm_cortical w li 2, " .
           sprintf("%f",$avg_gm_cortical-$stddev_gm_cortical) . " w li 3, " .
           sprintf("%f",$avg_gm_cortical+$stddev_gm_cortical) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GMCortical.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Volume of total cortical gray matter (without sub-cortical gray)." .
         "</td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Percentage error in the extraction of the white surface\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Percentage error in extraction of white surface\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"White surface error %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/WMSurfPct.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 10 w po, $avg_wm_surf w li 2, " .
           sprintf("%f",$avg_wm_surf-$stddev_wm_surf) . " w li 3, " .
           sprintf("%f",$avg_wm_surf+$stddev_wm_surf) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"WMSurfPct.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Percentage of classified white voxels outside " .
         "the white surface (large values mean underexpanded surface). </td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Percentage error in the expansion of the gray surface\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Percentage error in expansion of gray surface\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Gray surface error %\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GMSurfPct.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 11 w po, $avg_gm_surf w li 2, " .
           sprintf("%f",$avg_gm_surf-$stddev_gm_surf) . " w li 3, " .
           sprintf("%f",$avg_gm_surf+$stddev_gm_surf) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GMSurfPct.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Percentage of classified gray voxels outside " .
         "the gray surface (large values mean underexpanded surface). </td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Number of self-intersections in resampled mid surface\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "set key\n";
print FILE "set title \"Number of self-intersections in resampled mid surface\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Number of self-intersections\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/SurfRegInter.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 12 t \"Left\" w po 1, " .
           "\"${qc_dir}/civet_${prefix}.glm\" u 13 t \"Right\" w po 2\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"SurfRegInter.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Number of self-intersections in resampled " .
         "mid surface as a result of surface registration. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Number of surface-surface intersections between wm-gm surfaces\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "set key\n";
print FILE "set title \"Number of surface-surface intersections between wm-gm surfaces\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Number of intersections\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/SurfSurfInter.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 14 t \"Left\" w po 1, " .
           "\"${qc_dir}/civet_${prefix}.glm\" u 15 t \"Right\" w po 2\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"SurfSurfInter.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Number of surface-surface intersections between " .
         "wm-gm surfaces. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Mean cortical thickness for left/right hemispheres\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Mean cortical thickness for left/right hemispheres\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Mean CT%\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/MeanCT.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 16 t \"Left\" w po, " .
           "\"${qc_dir}/civet_${prefix}.glm\" u 17 t \"Right\" w po," .
           " $avg_meanct w li 2, " .
           sprintf("%f",$avg_meanct-$stddev_meanct) . " w li 3, " .
           sprintf("%f",$avg_meanct+$stddev_meanct) . " w li 3\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"MeanCT.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Mean cortical thickness </td></tr>\n" .
         "<tr> <td></td> <td>Green line is the mean; blue lines are mean +/- 1 stddev. </td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Mean T1 intensity for classified WM and GM\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "set key\n";
print FILE "set title \"Mean T1 intensity for classified WM and GM\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"Mean T1\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/MeanT1WMGM.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 18 t \"WM\" w po, " .
           "\"${qc_dir}/civet_${prefix}.glm\" u 20 t \"GM\" w po\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"MeanT1WMGM.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Mean T1 intensity for classified WM and GM." .
         "</td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# StdDev T1 intensity for classified WM and GM\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "set key\n";
print FILE "set title \"StdDev T1 intensity for classified WM and GM\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"StdDev T1\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/StdDevT1WMGM.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 19 t \"WM\" w po, " .
           "\"${qc_dir}/civet_${prefix}.glm\" u 21 t \"GM\" w po\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"StdDevT1WMGM.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> StdDev T1 intensity for classified WM and GM." .
         "</td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Gyrification index for the left hemisphere\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Gyrification index for left hemisphere\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"GI\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GILeft.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 22 w po\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GILeft.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Gyrification index for the left hemisphere." .
         "</td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Gyrification index for the right hemisphere\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Gyrification index for right hemisphere\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"GI\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GIRight.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 23 w po\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GIRight.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Gyrification index for the right hemisphere." .
         "</td></tr>\n" .
         "</table>\n";

# ----------------------------------------------------------------------------------------------
$figure++;
print FILE "# Gyrification index for the full brain\n";
print FILE "clear\n";
print FILE "set grid\n";
print FILE "unset key\n";
print FILE "set title \"Gyrification index for full brain\"\n";
print FILE "set xlabel \"Subjects\"\n";
print FILE "set ylabel \"GI\"\n";
print FILE "set terminal png\n";
print FILE "set output \"${qc_dir}/GIFull.png\"\n";
print FILE "plot \"${qc_dir}/civet_${prefix}.glm\" u 24 w po\n";

$html .= "<table>\n" .
         "<tr> <td></td> <td> <img src=\"GIFull.png\"> </td></tr>\n" .
         "<tr> <td> Figure $figure: </td> <td> Gyrification index for the full brain." .
         "</td></tr>\n" .
         "</table>\n";

close FILE;

`gnuplot \< "${qc_dir}/civet_${prefix}.gnu"`;


##
## <img src="toto.png">
## 

open FILE, ">${qc_dir}/civet_${prefix}.html";
print FILE $html;
close FILE;

############# Voila!! #############

# Set interrupt handler (for cleaning of lock files)

sub CLEANUP {
  exit(1);
}