#! /usr/bin/perl
#
# Simple Quality checks on the white and gray surfaces.
#   - Count how many classified white voxels outside white surface
#   - Count how many classified gray voxels outside gray surface
#
# Copyright Alan C. Evans
# Professor of Neurology
# McGill University
#

use strict;
use warnings "all";

use Getopt::Tabular;
use MNI::Startup;
use MNI::FileUtilities;
use MNI::Spawn;

MNI::Spawn::RegisterPrograms
  ( [qw/  rm
     mv
     minccalc
     mincstats
     mincresample
     surface_mask2/ ] )
  or exit 1;

# --- set the help & usage strings ---
my $help = <<HELP;
Simple Quality checks on the white and gray surfaces.

HELP

my $usage = <<USAGE;
usage: $ProgramName cls_correct.mnc white_mask_left.mnc white_mask_right.mnc \
white_surf_left.obj white_surf_right.obj gray_surf_left.obj gray_surf_right.obj
USAGE

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

# --- initialization ---
my @options = 
  ( @DefaultArgs,     # from MNI::Startup
  );

GetOptions( \@options, \@ARGV )
  or exit 1;
die "$usage\n" unless @ARGV == 9;

my $cls_correct = shift;
my $white_mask_left = shift;
my $white_mask_right = shift;
my $white_left = shift;
my $white_right = shift;
my $gray_left = shift;
my $gray_right = shift;
my $brain_mask = shift;
my $info_file = shift;

# These are temporary files.
MNI::FileUtilities::check_output_path("${TmpDir}")
or exit 1;

#  Quality checks on the white and gray surfaces.

my ($ret,$bla,$sum_left,$sum_right,$sum_white,$sum_gray,$sum);

# White surface:

# fill up the white_masks if there are small holes inside where
# the white surface cannot reach (for a better quality check).

&run( "mincdefrag", $white_mask_left, "${TmpDir}/white_mask_left_filled.mnc", 0, 6 );
&run( "mincdefrag", $white_mask_right, "${TmpDir}/white_mask_right_filled.mnc", 0, 6 );

# Resample at half resolution for more accurate calculation. Do
# the resampling by hand because of a bug in autocrop with start.

my $dx = `mincinfo -attvalue xspace:step $white_mask_left`; chomp( $dx );
my $dy = `mincinfo -attvalue yspace:step $white_mask_left`; chomp( $dy );
my $dz = `mincinfo -attvalue zspace:step $white_mask_left`; chomp( $dz );
my $nx = `mincinfo -attvalue xspace:length $white_mask_left`; chomp( $nx );
my $ny = `mincinfo -attvalue yspace:length $white_mask_left`; chomp( $ny );
my $nz = `mincinfo -attvalue zspace:length $white_mask_left`; chomp( $nz );
my $sx = `mincinfo -attvalue xspace:start $white_mask_left`; chomp( $sx );
my $sy = `mincinfo -attvalue yspace:start $white_mask_left`; chomp( $sy );
my $sz = `mincinfo -attvalue zspace:start $white_mask_left`; chomp( $sz );
$sx -= 0.25 * $dx;  # This shift is missing in autocrop.
$sy -= 0.25 * $dy;
$sz -= 0.25 * $dz;
$dx *= 0.5;
$dy *= 0.5;
$dz *= 0.5;
$nx *= 2;
$ny *= 2;
$nz *= 2;
&run( 'mincresample', '-clob', '-nearest', '-start', $sx, $sy, $sz, 
      '-step', $dx, $dy, $dz, '-nelements', $nx, $ny, $nz, 
      "${TmpDir}/white_mask_left_filled.mnc",
      "${TmpDir}/white_mask_left_filled_iso.mnc" );
unlink( "${TmpDir}/white_mask_left_filled.mnc" );
&run( 'mincresample', '-clob', '-nearest', '-start', $sx, $sy, $sz, 
      '-step', $dx, $dy, $dz, '-nelements', $nx, $ny, $nz, 
      "${TmpDir}/white_mask_right_filled.mnc",
      "${TmpDir}/white_mask_right_filled_iso.mnc" );
unlink( "${TmpDir}/white_mask_right_filled.mnc" );

# number of voxels classified as white
$sum_left = `mincstats -quiet -sum "${TmpDir}/white_mask_left_filled_iso.mnc"`;
$sum_right = `mincstats -quiet -sum "${TmpDir}/white_mask_right_filled_iso.mnc"`;
chomp( $sum_left );
chomp( $sum_right );
$sum_white = 1 + $sum_left + $sum_right;

# Get the voxels inside the surface.

&run( "surface_mask2", "-binary_mask", 
      "${TmpDir}/white_mask_left_filled_iso.mnc", $white_left, 
      "${TmpDir}/wqc_tmp.mnc" );

# Unfortunately, surface_mask2 puts the outlining voxels outside the surface
# so try to exclude the border line of the surface as it will otherwise
# count as a layer of voxel. We could use scan_object_to_volume, but the
# latter creates a thick outline (about 2 voxels thick) so removing this
# thick outline would remove too much. Instead, simply do a dilation of
# the volume into the mask (this should be a decent compromise).

&run( "dilate_volume", "${TmpDir}/wqc_tmp.mnc", "${TmpDir}/wqc.mnc", 0, 6, 1 );

# Compute difference between the surface mask and the target mask. This is not
# the most accurate because there are errors "voxelizing" the surface.

&run( "minccalc", "-clobber", 
      "-expression", 'if(abs(A[0]-A[1])>0.5){out=1;}else{out=0;}',
      "${TmpDir}/white_mask_left_filled_iso.mnc", "${TmpDir}/wqc.mnc", 
      "${TmpDir}/diff.mnc" );

$sum_left = `mincstats -quiet -sum "${TmpDir}/diff.mnc"`;
chomp( $sum_left );

# Repeat for right surface.

&run( "surface_mask2", "-binary_mask", 
      "${TmpDir}/white_mask_right_filled_iso.mnc", $white_right, 
      "${TmpDir}/wqc_tmp.mnc" );
&run( "dilate_volume", "${TmpDir}/wqc_tmp.mnc", "${TmpDir}/wqc.mnc", 0, 6, 1 );
&run( "minccalc", "-clobber", 
      "-expression", 'if(abs(A[0]-A[1])>0.5){out=1;}else{out=0;}',
      "${TmpDir}/white_mask_right_filled_iso.mnc", "${TmpDir}/wqc.mnc", 
      "${TmpDir}/diff.mnc" );
$sum_right = `mincstats -quiet -sum "${TmpDir}/diff.mnc"`;
chomp( $sum_right );

unlink( "${TmpDir}/white_mask_left_filled_iso.mnc" );
unlink( "${TmpDir}/white_mask_right_filled_iso.mnc" );
unlink( "${TmpDir}/wqc.mnc" );
unlink( "${TmpDir}/diff.mnc" );

# Total for left+right.

my $info = sprintf( "white surface (%5.2f%%), ", 100.0*(${sum_left}+${sum_right})/${sum_white} );

# Gray surface:

# number of voxels classified as gray (label 2).

&run( "minccalc", "-clobber", "-expression", 
      'if(A[0]>1.5&&A[1]>0.5){out=A[0];}else{out=0;}',
      $cls_correct, $brain_mask, "${TmpDir}/diff.mnc" );
&run( "mincdefrag", "${TmpDir}/diff.mnc", "${TmpDir}/diff.mnc", 2, 6 );

# Resample at half resolution for more accurate calculation. Do
# the resampling by hand because of a bug in autocrop with start.

$dx = `mincinfo -attvalue xspace:step $cls_correct`; chomp( $dx );
$dy = `mincinfo -attvalue yspace:step $cls_correct`; chomp( $dy );
$dz = `mincinfo -attvalue zspace:step $cls_correct`; chomp( $dz );
$nx = `mincinfo -attvalue xspace:length $cls_correct`; chomp( $nx );
$ny = `mincinfo -attvalue yspace:length $cls_correct`; chomp( $ny );
$nz = `mincinfo -attvalue zspace:length $cls_correct`; chomp( $nz );
$sx = `mincinfo -attvalue xspace:start $cls_correct`; chomp( $sx );
$sy = `mincinfo -attvalue yspace:start $cls_correct`; chomp( $sy );
$sz = `mincinfo -attvalue zspace:start $cls_correct`; chomp( $sz );
$sx -= 0.25 * $dx;  # This shift is missing in autocrop.
$sy -= 0.25 * $dy;
$sz -= 0.25 * $dz;
$dx *= 0.5;
$dy *= 0.5;
$dz *= 0.5;
$nx *= 2;
$ny *= 2;
$nz *= 2;
&run( 'mincresample', '-clob', '-nearest', '-start', $sx, $sy, $sz, 
      '-step', $dx, $dy, $dz, '-nelements', $nx, $ny, $nz, 
      "${TmpDir}/diff.mnc", "${TmpDir}/diff_iso.mnc" );
unlink( "${TmpDir}/diff.mnc" );

$sum_gray = `mincstats -quiet -count -mask "${TmpDir}/diff_iso.mnc" -mask_binvalue 2 "${TmpDir}/diff_iso.mnc"`;
chomp( $sum_gray );
$sum_gray += 1;

# voxels outside gray surface that are classified as gray matter (label 2)
&run( "surface_mask2", "-binary_mask", "${TmpDir}/diff_iso.mnc",
       $gray_left, "${TmpDir}/gqc_left.mnc" );
&run( "surface_mask2", "-binary_mask", "${TmpDir}/diff_iso.mnc",
       $gray_right, "${TmpDir}/gqc_right.mnc" );

# dilate before combining (we may lose centerline, but this avoids outer 
# ouline of gray surfaces to cross - will not be able to dilate if touching).
 
&run( "dilate_volume", "${TmpDir}/gqc_left.mnc", 
      "${TmpDir}/gqc_left.mnc", 0, 6, 1 );
&run( "dilate_volume", "${TmpDir}/gqc_right.mnc", 
      "${TmpDir}/gqc_right.mnc", 0, 6, 1 );

# Combine left + right, 
&run( "minccalc", "-clobber", "-expression", 'if( A[0]>0.5 || A[1]>0.5 ) {out=1;} else {out=0;};', 
      "${TmpDir}/gqc_left.mnc", "${TmpDir}/gqc_right.mnc", "${TmpDir}/gqc.mnc" );
unlink( "${TmpDir}/gqc_left.mnc" );
unlink( "${TmpDir}/gqc_right.mnc" );

&run( "minccalc", "-clobber", "-expression", 
      'if((A[0]>0.5&&A[1]<0.5)||(A[0]<0.5&&A[1]>0.5)){1}else{0}',
      "${TmpDir}/diff_iso.mnc", "${TmpDir}/gqc.mnc", 
      "${TmpDir}/diff.mnc" );
$sum = `mincstats -quiet -sum "${TmpDir}/diff.mnc"`;
chomp( $sum );

unlink( "${TmpDir}/diff.mnc" );
unlink( "${TmpDir}/diff_iso.mnc" );
unlink( "${TmpDir}/gqc.mnc" );

$info .= sprintf( "gray surface (%5.2f%%)\n", 100.0*${sum}/${sum_gray} );

open PIPE, ">$info_file";
print PIPE $info;
close PIPE;

#Execute a system call.

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