marching_cubes.pl.in

#! @PERL@
#
# Author: Claude Lepage <claude@bic.mni.mcgill.ca>
#
# May 2011
#
# Copyright Alan C. Evans
# Professor of Neurology
# McGill University
#

# Extract a white matter surface using a marching-cubes
# algorithm, in a first step, then run ASP to complete
# the convergence.

use strict;
use warnings "all";
use File::Basename;
use File::Spec;
use File::Temp qw/ tempdir /;

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

# --- set the help & usage strings ---
my $help = <<HELP;
Required parameters:
  wm_mask.mnc  : white matter mask
  white.obj    : white matter surface (output)
HELP

my $license = <<LICENSE;
Copyright Alan C. Evans
Professor of Neurology
McGill University

LICENSE

my $usage = <<USAGE;
Usage: $ProgramName t1.mnc cls.mnc skull_mask.mnc wm_mask.mnc white.obj 
       surf_reg_model white_surf_mask.txt 
       [-left] [-right] [-refine] [-subsample] [-calibrate]
       $ProgramName -help to list options

$license
USAGE

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

my $side = undef;
my $refine = 0;
my $subsample = 0;
my $calibrate = 0;
my $mask = undef;
my $mask_label = 0;

my @options = ( 
  ['-left', 'const', "Left", \$side, "Extract left surface"],
  ['-right', 'const', "Right", \$side, "Extract right surface"],
  ['-refine', 'boolean', 0, \$refine,
   "Create a high-resolution surface at 327680 polygons"],
  ['-subsample', 'boolean', 0, \$subsample,
   "Subsample white matter mask at half voxel size"],
  ['-calibrate', 'boolean', 0, \$calibrate,
   "Calibrate white surface to gray-white gradient"],
  ['-mask', 'string', 1, \$mask,
   "Surface mask where to ignore gray-white gradient correction"],
  ['-mask-label', 'string', 1, \$mask_label,
   "Surface mask label where to ignore gray-white gradient correction"],
  );

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

my $t1 = shift;
my $cls = shift;
my $skull_mask = shift;
my $original_white_matter_mask = shift;
my $white_surface = shift;
my $surfreg_model = shift;
my $white_surf_mask = shift;

if( !( defined $side ) ) {
  print "$usage\n";
  die "You must specify left or right hemisphere.\n";
}
if( !( -e $t1 ) ) {
  print "$usage\n";
  die "t1 image must exist.\n";
}
if( !( -e $cls ) ) {
  print "$usage\n";
  die "Tissue classification image must exist.\n";
}
if( !( -e $skull_mask ) ) {
  print "$usage\n";
  die "Brain mask must exist.\n";
}
if( !( -e $original_white_matter_mask ) ) {
  print "$usage\n";
  die "White matter mask must exist.\n";
}
if( $surfreg_model ne "none" && !( -e $surfreg_model ) ) {
  print "$usage\n";
  die "You must specify a model for surface registration.\n";
}

#my $surfreg = "geom_surfreg.pl";
my $surfreg = "bestsurfreg.pl";

my $tmpdir = &tempdir( "mcubes-XXXXXX", TMPDIR => 1, CLEANUP => 1 );

my $ICBM_white_model = MNI::DataDir::dir("surface-extraction") .
                       "/white_model_320.obj";
my $initial_model = "${tmpdir}/initial_white_model.obj";
if( $side eq "Left" ) {
  &run( "param2xfm", "-translation", -25, 0, 0,
        "${tmpdir}/slide_left.xfm" );
  &run( "transform_objects", $ICBM_white_model,
        "${tmpdir}/slide_left.xfm", $initial_model );
  unlink( "${tmpdir}/slide_left.xfm" );
} else {
  &run( "param2xfm", "-scales", -1, 1, 1,
        "${tmpdir}/flip.xfm" );
  &run( "transform_objects", $ICBM_white_model,
        "${tmpdir}/flip.xfm", $initial_model );
  unlink( "${tmpdir}/flip.xfm" );
  &run( "param2xfm", "-translation", 25, 0, 0,
        "${tmpdir}/slide_right.xfm" );
  &run( "transform_objects", $initial_model,
        "${tmpdir}/slide_right.xfm", $initial_model );
  unlink( "${tmpdir}/slide_right.xfm" );
}

# Make sure that the mask is not empty.

my $sum = `mincstats -quiet -sum $original_white_matter_mask`;
chomp( $sum );
die "Empty white matter mask.\n" if( $sum == 0 );

# Remove loosely connected voxels and fill-in tightly connected voxels 
# to have a smoother surface.

open KRNL, "> ${tmpdir}/ngh_count.kernel";
print KRNL "MNI Morphology Kernel File\n";
print KRNL "Kernel_Type = Normal_Kernel;\n";
print KRNL "Kernel =\n";
print KRNL " 1.0  0.0  0.0  0.0  0.0  1.0\n";
print KRNL "-1.0  0.0  0.0  0.0  0.0  1.0\n";
print KRNL " 0.0  1.0  0.0  0.0  0.0  1.0\n";
print KRNL " 0.0 -1.0  0.0  0.0  0.0  1.0\n";
print KRNL " 0.0  0.0  1.0  0.0  0.0  1.0\n";
print KRNL " 0.0  0.0 -1.0  0.0  0.0  1.0;\n";
close KRNL;

my $ngh_count = "${tmpdir}/ngh_count.mnc";
my $white_matter_mask = "${tmpdir}/wm_mask_clean.mnc";
my $tmp_wm_mask = "${tmpdir}/tmp_wm_mask.mnc";

&run( 'cp', '-f', $original_white_matter_mask, $white_matter_mask );
my $ret = `mincstats -sum $white_matter_mask`; chomp( $ret );
print "Connectivity iteration 0: $ret\n";

# We will use a white matter mask with simplified connectivity
# for the raw marching-cubes surface, but we will later fit the
# surface to the original white matter mask.

my $lower = 2.5;   ## 1.5;
for( my $iter = 1; $iter <= 5; $iter++ ) {
  &run( 'mincmorph', '-clobber', '-unsigned', '-byte', '-convolve', '-kernel', 
        "${tmpdir}/ngh_count.kernel", $white_matter_mask, $ngh_count );
  &run( 'minccalc', '-clob', '-quiet', '-expr', 
        "if(A[1]<$lower){0}else{if(A[1]>4.5){1}else{A[0]}}",
        $white_matter_mask, $ngh_count, $tmp_wm_mask );
  &run( 'mv', '-f', $tmp_wm_mask, $white_matter_mask );
  my $ret = `mincstats -sum $white_matter_mask`; chomp( $ret );
  print "Connectivity iteration $iter: $ret\n";
  ## $lower = 0.5;
}

unlink( $ngh_count );
unlink( "${tmpdir}/ngh_count.kernel" );

# Preparation of the white matter mask for extraction of the surface
# using marching cubes.

&run( 'minccalc', '-clobber', '-quiet', '-expression', 'out=1',
      '-unsigned', '-byte', $white_matter_mask, "${tmpdir}/filled.mnc" );
&run( 'surface_mask2', '-binary', "${tmpdir}/filled.mnc",
      $initial_model, "${tmpdir}/s0.mnc" );
&run( 'mincresample', '-clobber', '-quiet', '-like',
      "${tmpdir}/filled.mnc", "${tmpdir}/s0.mnc",
      "${tmpdir}/s1.mnc" );
&run( 'minccalc', '-clobber', '-quiet', '-expression',
      'if(A[0]>0.5||A[1]>0.5){1}else{0}',
      $white_matter_mask, "${tmpdir}/s1.mnc", "${tmpdir}/s0.mnc" );
&run( 'dilate_volume', "${tmpdir}/s0.mnc", "${tmpdir}/s1.mnc", 1, 26, 1 );
&run( 'minccalc', '-clobber', '-quiet', '-expression', 'A[0]+A[1]', 
      $white_matter_mask, "${tmpdir}/s1.mnc", "${tmpdir}/s0.mnc" );
&run( 'mincreshape', '-quiet', '-clobber', '-unsigned', '-byte',
      '-image_range', 0, 255, '-valid_range', 0, 255, "${tmpdir}/s0.mnc",
      "${tmpdir}/s1.mnc" );
&run( 'mincdefrag', "${tmpdir}/s1.mnc", "${tmpdir}/s0.mnc", 2, 19 );  ## was 6 before
unlink( "${tmpdir}/filled.mnc" );
unlink( "${tmpdir}/s1.mnc" );

# Do not use sub-sampling if voxel resolution is already below 1mm.
# Too slow.

my $dx = `mincinfo -attvalue xspace:step $white_matter_mask`; chomp($dx);
$subsample = 0 if( $dx < 0.90 );
unlink( $white_matter_mask );

# Crop volume to smallest extents to speed up marching-cubes.
my @crop = split( ' ', `mincbbox ${tmpdir}/s0.mnc -mincreshape` );
&run( 'mincreshape', '-quiet', '-clobber', @crop, "${tmpdir}/s0.mnc",
      "${tmpdir}/s1.mnc" );
unlink( "${tmpdir}/s0.mnc" );

# Run the marching-cubes algorithm on the mask.

&run( 'sphere_mesh', "${tmpdir}/s1.mnc", $white_surface, 
      ( $subsample ) ? "-subsample" : "" );
unlink( "${tmpdir}/s1.mnc" );

# Coarsen and smooth the original marching-cubes surface.

my $white_surface_sm = "${tmpdir}/white_surf_sm.obj";
&run( 'adapt_object_mesh', $white_surface, $white_surface_sm, 120000, 1, 50, 1 );

# Inflate the white surface onto a unit sphere.

my $white_sphere_sm = "${tmpdir}/white_sphere_sm.obj";
&run( 'inflate_to_sphere', $white_surface_sm, $white_sphere_sm );

# Interpolate from sphere-to-sphere to resample the white surface 
# using the 40962 vertices on the standard ICBM surface average 
# template. This unit sphere is the one used for surface registration.

my $unit_sphere = "${tmpdir}/unit_sphere.obj";
&run( 'create_tetra', $unit_sphere, 0, 0, 0, 1, 1, 1, 81920 );
if( $side eq "Right" ) {
  &run( "param2xfm", "-scales", -1, 1, 1,
        "${tmpdir}/flip.xfm" );
  &run( "transform_objects", $unit_sphere,
        "${tmpdir}/flip.xfm", $unit_sphere );
  unlink( "${tmpdir}/flip.xfm" );
}

# Evaluate the white surface from the marching-cubes surface.

&resample_white_surface( $white_surface_sm, $white_sphere_sm,
                         $unit_sphere, $white_surface );

# Do surface registration on to the average white population model.
# This will be a first alignment that will ensure better isotropic mesh.

if( $surfreg_model ne "none" ) {

  &run( 'adapt_object_mesh', $white_surface, $white_surface, 0, 10, 0, 0 );

  &run( $surfreg, '-clobber', '-min_control_mesh', 320, '-max_control_mesh',
        20480, '-blur_coef', '1.0', '-neighbourhood_radius', '2.8', 
        '-mode', 'stiff', $surfreg_model, $white_surface, "${tmpdir}/white.sm" );

  # Resample the unit sphere.

  my $unit_sphere_rsl = "${tmpdir}/unit_sphere_rsl.obj";
  &run( 'sphere_resample_obj', '-clobber', $unit_sphere, "${tmpdir}/white.sm",
        $unit_sphere_rsl );
  unlink( "${tmpdir}/white.sm" );

  # Remove self-intersections in unit_sphere_rsl. Yes, there can be some due
  # to surface registration itself. I know, there shouldn't be. First, scale
  # the unit sphere to radius 100 because the hard-coded parameters in 
  # check_self_intersect are based on the size of a human head. 

  my $big_sphere_rsl = "${tmpdir}/big_sphere_rsl.obj";
  &run( 'param2xfm', '-scales', 100, 100, 100, "${tmpdir}/scale100.xfm" );
  &run( 'transform_objects', $unit_sphere_rsl, "${tmpdir}/scale100.xfm",
        $big_sphere_rsl );
  my @ret = `check_self_intersect $big_sphere_rsl -fix $big_sphere_rsl`;
  $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
  my $num_inter = $1;
  if( $num_inter > 0 ) {
    my $prev_num_inter = $num_inter;
    my $increasing = 0;
    for( my $i = 1; $i <= 500; $i++ ) {
      my @ret = `check_self_intersect $big_sphere_rsl -fix $big_sphere_rsl`;
      $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
      $num_inter = $1;
      printf "Iter = $i  Self-Intersections = $num_inter\n";
      last if( $num_inter == 0 );
      if( $num_inter >= $prev_num_inter ) {
        if( $increasing >= 0 ) {
          $increasing++;
        } else {
          $increasing = 1;
        }
      } else {
        $increasing = 0;
      }
      $prev_num_inter = $num_inter;
      last if( $increasing >= 5 && $i > 50 );
    }
    if( $num_inter > 0 ) {
      print "Warning: There are still $num_inter self-intersections in resampled sphere.\n";
    }
    &run( 'flatten_to_sphere', $big_sphere_rsl, $unit_sphere_rsl );
  
    # no need to normalize sphere to unit radius. It's done in interpolate_sphere.
    # &run( 'xfminvert', "${tmpdir}/scale100.xfm", "${tmpdir}/scale100.xfm" );
    # &run( 'transform_objects', $big_sphere_rsl, "${tmpdir}/scale100.xfm",
    #       $unit_sphere_rsl );
  }
  unlink( $big_sphere_rsl );
  unlink( "${tmpdir}/scale100.xfm" );

  # New interpolation of surface based on points on the resampled sphere.
  # Note: no self-inter in unit_sphere_rsl. Can have self-inter in white_sphere_sm
  #       if inflation to sphere is difficult.
  
  &run( 'interpolate_sphere', $white_surface_sm, $white_sphere_sm, $unit_sphere_rsl,
        $white_surface );

  unlink( $unit_sphere_rsl );
  unlink( $white_surface_sm );
  unlink( $white_sphere_sm );

}
unlink( $unit_sphere );

# Check for self-intersections in the marching-cubes surface.

my $prev_num_inter = 999999;
my $increasing = 0;
for( my $i = 1; $i <= 500; $i++ ) {
  my @ret = `check_self_intersect $white_surface -fix $white_surface`;
  $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
  my $num_inter = $1;
  printf "Iter = $i  Self-Intersections = $num_inter\n";
  if( $num_inter == 0 ) {
    $prev_num_inter = 0;
    last;
  }
  if( $num_inter >= $prev_num_inter ) {
    if( $increasing >= 0 ) {
      $increasing++;
    } else {
      $increasing = 1;
    }
  } else {
    $increasing = 0;
  }
  $prev_num_inter = $num_inter;
  last if( $increasing >= 5 && $i > 50 );
}
if( $prev_num_inter > 0 ) {
  my $failed_surface = $white_surface;
  $failed_surface =~ s/\.obj$/-failed\.obj/;
  `mv -f $white_surface $failed_surface`;
  # unlink( $white_surface );   ## is this desirable???
  die "Failed interpolation of marching-cubes surface with $prev_num_inter self-intersections.\n";
}

# Run ASP on the resampled white surface to converge it fully
# to the white matter mask. The ICBM model is used by ASP to
# define the distribution of the edge lengths for the stretch
# constraint. Here, we use the original white matter mask 
# without the changes for dangling white voxels.

&run_asp( $white_surface, $original_white_matter_mask, 0.5, $initial_model,
          ($refine && !$calibrate), $tmpdir );

# Run ASP on the resampled white surface to move the surface to
# the maximum gm-wm gradient. cls must be pve_classify, without
# ventricles filled.

if( $calibrate ) {
  &calibrate_white( $t1, $cls, $skull_mask, $white_surface, 
                    $white_surf_mask, $initial_model, 
                    $refine, $tmpdir );
}

# Resample the white surface from a standard sphere from the 
# hi-res marching-cubes white surface. The distribution of
# vertices on the sphere is adapted such as to produce an
# interpolated surface with triangles of nearly the same size.

sub resample_white_surface {

  my $white_mc = shift;     # hi-res raw marching-cubes surface
  my $sphere_mc = shift;    # sphere corresponding to white_mc
  my $unit_sphere = shift;  # standard sphere
  my $output = shift;       # output white surface with uniform triangles

  my @conf = ( { 'size' => 320,     # most of the motion occurs early
                 'fwhm' => 20.0,
                 'niter' => 500 },
               { 'size' => 1280,
                 'fwhm' => 10.0,
                 'niter' => 500 },
               { 'size' => 5120,
                 'fwhm' => 5.0,
                 'niter' => 300 },
               { 'size' => 20480,
                 'fwhm' => 2.0,
                 'niter' => 150 } );

  my $start = 320;
  my $end = 20480;

  my $npolys = `print_n_polygons $unit_sphere`;
  chomp( $npolys );

  my $current_sphere = "${tmpdir}/current_sphere.obj";
  &run( 'cp', '-f', $unit_sphere, $current_sphere );

  # obtain initial white surface

  &run( 'interpolate_sphere', $white_mc, $sphere_mc,
        $current_sphere, $output );

  # Multi-resolution approach from coarse to fine mesh.

  for( my $idx = 0; $idx <= $#conf; $idx++ ) {

    my $size = $conf[$idx]{size};

    next if( $size < $start );
    last if( $size > $end );

    print "Sphere adaptation at $size vertices...\n";

    # Obtain the triangle areas from current white surface to
    # the current sphere at size npolys.

    my $white_area = "${tmpdir}/white_area.txt";
    my $sphere_area = "${tmpdir}/sphere_area.txt";
    &run( 'depth_potential', '-area_simple', $output, $white_area );
    &run( 'depth_potential', '-area_simple', $current_sphere, $sphere_area );
    &run( 'vertstats_math', '-old_style_file', '-div', $white_area,
          $sphere_area, $white_area );
    unlink( $sphere_area );
    if( $conf[$idx]{fwhm} > 0 ) {
      &run( 'depth_potential', '-smooth', $conf[$idx]{fwhm},
            $white_area, $output, $white_area );
    }

    # adapt the current_sphere at this size based on the areas.
  
    &subdivide_mesh( $current_sphere, $size, $current_sphere );
    &run( 'adapt_metric', $current_sphere, $white_area,
          $current_sphere, $conf[$idx]{niter} );
    unlink( $white_area );

    # interpolate relative to the original white surface at npolys.

    &subdivide_mesh( $current_sphere, $npolys, $current_sphere );
    &run( 'interpolate_sphere', $white_mc, $sphere_mc,
          $current_sphere, $output );

  }

  # Create a new hi-res background mesh with uniform triangles.
  # This new background mesh will be used for interpolating the
  # resampled white surface after surface registration. This way,
  # we can associate the standard sphere to this new bg mesh 
  # since the standard sphere is used during surface registration.

  $npolys *= 4;
  &subdivide_mesh( $current_sphere, $npolys, $current_sphere );

  &run( 'interpolate_sphere', $white_mc, $sphere_mc,
        $current_sphere, $white_mc );

  &subdivide_mesh( $unit_sphere, $npolys, $sphere_mc );

  unlink( $current_sphere );

}

# Run ASP on the resampled white surface to converge it fully
# to the white matter mask. This is a simplified version of 
# extract_white_surface without the coarse steps.

sub run_asp {

  my $white = shift;
  my $wm_mask = shift;
  my $isovalue = shift;
  my $white_model = shift;
  my $refine = shift;
  my $tmpdir = shift;

  my $self_dist2 = 0.01;
  my $self_weight2 = 1e06;
  my $n_selfs = 9;

  my $stop_threshold = 1e-3;
  my $stop_iters = 1000;
  my $n_per = 5;
  my $tolerance = 1.0e-03;
  my $f_tolerance = 1.0e-06;
  my $oo_scale = 0.5;

  my @schedule = (
  #  size    sw   n_itr  inc   offo   offi   si  in  out  over   sw   self
  #  -----  ----  -----  ----  ----   ----   --  --  ---  ----  ----  ----
    81920,  100,   200,   50,   50,    1,    1,  20,  5,    2,  1e0,  .25,
    81920,   20,   200,   50,   50,    1,    1,  20,  3,    2,  1e0,  .25,
    81920,    5,   200,   50,   50,    5,    1,  10,  2,    2,  1e0,  .25,
    81920,    2,   200,   50,   50,    5,  0.5,  10,  2,    2,  1e0,  .25,
    81920,    1,   200,   50,   20,   10,  0.5,   5,  1,    2,  1e0,  .25,
    # these 2 stages are to "plug" holes
    81920,    1,   200,   50,    2,    2,  0.5,   1,  1,    2,  1e0,  .25,
    81920,    1,   100,   50,   10,    5,  0.5,   2,  1,    2,  1e0,  .25,
   327680,  5e4,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
   327680,  2e4,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
   327680,  1e4,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
   327680,  5e3,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
   327680,  2e3,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
   327680,  1e3,    50,   50,   10,    2, 0.25,   5,  1,    1,  1e0, .125,
  );
  my $sched_size =  12;
  my $num_steps = @schedule / $sched_size;

  &run( 'cp', '-f', $white_model, "${tmpdir}/white_model_tmp.obj" );
  $white_model = "${tmpdir}/white_model_tmp.obj";
  subdivide_mesh( $white_model, 81920, $white_model );

  for( my $i = 0;  $i < @schedule;  $i += $sched_size ) {
    my ( $size, $sw, $n_iters, $iter_inc, $offo, $offi,
         $si_step, $in_dist, $out_dist, $oversample, $self_weight,
         $self_dist ) = @schedule[$i..$i+$sched_size-1];
    if( !( $refine ) ) {
      last if( $size == 327680 );
    }

    my $prev_size = `print_n_polygons $white`;
    chomp( $prev_size );
    if( $prev_size != $size ) {
      my @objsuffix = ( ".obj", "_320", "_1280", "_5120", "_20480",
                        "_81920", "_327680" );
      my $white_dir = dirname( $white );
      my $white_prefix = basename( $white, @objsuffix ) ;
      $white_prefix = "${white_dir}/${white_prefix}";
      subdivide_mesh( $white, $size, "${white_prefix}_$size.obj" );
      $white = "${white_prefix}_$size.obj";
      subdivide_mesh( $white_model, $size, $white_model );
    }

    $oversample *= $oo_scale;
    my $self2 = get_self_intersect( $self_weight, $self_weight2, $n_selfs, 
                                    $self_dist, $self_dist2 );
  
    my $b2 = " -boundary $offo $offi $wm_mask " .
             " $isovalue - $out_dist $in_dist 0 0 $oversample ";
  
    for( my $iter = 0;  $iter < $n_iters;  $iter += $iter_inc ) {
      print "echo Step ${size}: $iter / $n_iters    $sw\n";

      my $ni = $n_iters - $iter;
      if( $ni > $iter_inc )  { $ni = $iter_inc; }

      # Add a little bit of Taubin smoothing between cycles.
      &taubinize_surface( $white, 3 );
    
      my $command = "surface_fit -mode two -surface $white $white " .
                    " -stretch $sw $white_model -.9 0 0 0 " .
                    " $b2 $self2 -step $si_step " .
                    " -fitting $ni $n_per $tolerance " .
                    " -ftol $f_tolerance " .
                    " -stop $stop_threshold $stop_iters ";
      system( $command ) == 0 or die "Command $command failed with status: $?";
    }
  }
  unlink( $white_model );
}

# Run ASP on the resampled white surface to move the surface to
# the maximum gm-wm gradient.

sub calibrate_white {

  my $t1 = shift;
  my $cls = shift;
  my $skull_mask = shift;
  my $white = shift;
  my $white_surf_mask = shift;
  my $white_model = shift;
  my $refine = shift;
  my $tmpdir = shift;

  my $self_dist2 = 0.01;
  my $self_weight2 = 1e06;
  my $n_selfs = 9;

  my $stop_threshold = 1e-3;
  my $stop_iters = 1000;

  my $n_per = 5;
  my $tolerance = 1.0e-03;
  my $f_tolerance = 1.0e-06;
  my $oo_scale = 0.5;

  my @schedule = (
  #  size   blur   sw  n_itr  inc  gdist  ginc  gwgt over   si   sw   self
  #  -----  ----  ---- -----  ---- -----  ----  ---- ----   --  ----  ----
    81920,   90,  100,  150,   50,  5.0,  0.1,  1.0,  1,   1.0,  1e0,  .25,
    81920,   50,   20,  200,   50,  5.0,  0.1,  2.0,  1,   0.5,  1e0,  .25,
    81920,   20,   10,  200,   50,  5.0, 0.05,  5.0,  1,   0.5,  1e0,  .25,
    81920,   10,   10,  200,   50,  5.0, 0.05, 10.0,  1,   0.5,  1e0,  .25,
   327680,   10,  5e4,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
   327680,   10,  2e4,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
   327680,   10,  1e4,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
   327680,   10,  5e3,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
   327680,   10,  2e3,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
   327680,   10,  1e3,   50,   50,  5.0, 0.05, 10.0,  1,   0.25, 1e0, .125,
  );

  my $sched_size =  12;
  my $num_steps = @schedule / $sched_size;

  ##
  ## TO DO: clean-up t1 image to remove local min/max
  ##

  # Apply skull mask on the cls image. We want to preserve the original
  # classification with CSF in the ventricles, not filled with WM, so
  # that the gradient correction can compute good tissue thresholds.

  my $cls_rsl = "${tmpdir}/cls_rsl.mnc";
  my $cls_mc = "${tmpdir}/cls_rsl_mc.mnc";
  &run( 'mincresample', '-clobber', '-quiet', '-like', $skull_mask,
        $cls, $cls_rsl );
  &run( 'minccalc', '-clobber', '-quiet', '-expression',
        "if(A[1]>0.5){out=A[0]}else{0}", $cls_rsl,
        $skull_mask, $cls_mc );
  unlink( $cls_rsl );

  # Create surface mask for medial plane and brainstem.
  my $white_surf_mask_mc = "${tmpdir}/white_surf_mask_mc.txt";
  &run( 'vertstats_math', '-old_style_file', '-const2', -0.5, 0.5, '-seg',
        $white_surf_mask, $white_surf_mask_mc );

  &run( 'cp', '-f', $white_model, "${tmpdir}/white_model_tmp.obj" );
  $white_model = "${tmpdir}/white_model_tmp.obj";
  subdivide_mesh( $white_model, 81920, $white_model );

  for( my $i = 0;  $i < @schedule;  $i += $sched_size ) {
    my ( $size, $blur, $sw, $n_iters, $iter_inc, $grad_dist, $grad_inc, 
         $grad_wgt, $oversample, $si_step, $self_weight, $self_dist ) = 
       @schedule[$i..$i+$sched_size-1];
    if( !( $refine ) ) {
      last if( $size == 327680 );
    }

    my $t1_blur = $t1;
    if( $blur > 0 ) {
      $t1_blur = &basename( $t1, ".mnc" );
      $t1_blur = "${tmpdir}/${t1_blur}" . "_blur.mnc";
      &run( 'geo_smooth', 0.004, $blur, $t1, $t1_blur );
    }

    my $prev_size = `print_n_polygons $white`;
    chomp( $prev_size );
    if( $prev_size != $size ) {
      my @objsuffix = ( ".obj", "_320", "_1280", "_5120", "_20480",
                        "_81920", "_327680" );
      my $white_dir = dirname( $white );
      my $white_prefix = basename( $white, @objsuffix ) ;
      $white_prefix = "${white_dir}/${white_prefix}";
      subdivide_mesh( $white, $size, "${white_prefix}_$size.obj" );
      $white = "${white_prefix}_$size.obj";
      subdivide_mesh( $white_model, $size, $white_model );
    }

    my $self2 = get_self_intersect( $self_weight, $self_weight2, $n_selfs,
                                    $self_dist, $self_dist2 );

    for( my $iter = 0;  $iter < $n_iters;  $iter += $iter_inc ) {
      print "echo Step ${size}: $iter / $n_iters    $sw\n";

      my $ni = $n_iters - $iter;
      if( $ni > $iter_inc )  { $ni = $iter_inc; }

      # Add a little bit of Taubin smoothing between cycles.
      &taubinize_surface( $white, 3 );
    
      my $command = "surface_fit -mode two -surface $white $white " .
                    " -stretch $sw $white_model -.9 0 0 0 " .
                    " $self2 -step $si_step " .
                    " -gw_gradient $grad_wgt $oversample $grad_dist $grad_inc $t1_blur " .
                    " $cls_mc $white_surf_mask_mc " .
                    " -fitting $ni $n_per $tolerance " .
                    " -ftol $f_tolerance " .
                    " -stop $stop_threshold $stop_iters ";
      system( $command ) == 0 or die "Command $command failed with status: $?";
    }
  }
  unlink( $white_model );
  unlink( $white_surf_mask_mc );
  unlink( $cls_mc );
}

# Add a little bit of Taubin smoothing between cycles. This
# can introduce self-intersections, so try to fix those as
# well, in any. If the surface cannot be improved, return the
# original.

sub taubinize_surface {

  my $surf = shift;
  my $iter = shift;

  my $tmp_surf = "${tmpdir}/surface_taubin.obj";

  &run( 'adapt_object_mesh', $surf, $tmp_surf, 0, $iter, 0, 0 );
  my $tries = 0;
  do {
    &run( 'check_self_intersect', $tmp_surf, '-fix', $tmp_surf );
    my @ret = `check_self_intersect $tmp_surf`;
    $ret[1] =~ /Number of self-intersecting triangles = (\d+)/;
    if( $1 == 0 ) {
      `mv -f $tmp_surf $surf`;
      $tries = 9999;
    }
    $tries++;
  } while( $tries < 10 );
  unlink( $tmp_surf ) if( -e $tmp_surf );
}

# from surface-extraction/deform_utils.pl

sub  get_self_intersect( $$$$$ ) {

    my( $self_weight, $self_weight2, $n_selfs, $self_dist, $self_dist2 ) = @_;
    my( $self, $weight, $weight_factor, $s, $dist );

    if( $self_weight > 0.0 ) {
        $self = "";
        $weight = $self_weight;
        $weight_factor = ( $self_weight2 / $self_weight ) **
                         ( 1.0 / $n_selfs );

        for( $s = 0;  $s < $n_selfs;  ++$s ) {
            $dist = $self_dist + ($self_dist2 - $self_dist) *
                    $s / $n_selfs;
            $self = $self . " -self_intersect $weight $dist ";
            $weight *= $weight_factor;
        }
        $self = $self . " -self_intersect $self_weight2 $self_dist2 ";
    } else {
        $self = "";
    }
    $self;
}

# Check if the input surface has the same side orientation (left)
# as the default template model.

sub CheckFlipOrientation {

  my $obj = shift;

  my $npoly = `print_n_polygons $obj`;
  chomp( $npoly );

  my $ret = `tail -5 $obj`;
  my @verts = split( ' ', $ret );
  my @last3 = ( $verts[$#verts-2], $verts[$#verts-1], $verts[$#verts] );

  my $dummy_sphere = "${tmpdir}/dummy_sphere.obj";
  &run('create_tetra',$dummy_sphere,0,0,0,1,1,1,$npoly);
  $ret = `tail -5 $dummy_sphere`;
  unlink( $dummy_sphere );
  @verts = split( ' ', $ret );
  my @sphere3 = ( $verts[$#verts-2], $verts[$#verts-1], $verts[$#verts] );
  if( $last3[0] == $verts[$#verts-2] &&
      $last3[1] == $verts[$#verts-1] &&
      $last3[2] == $verts[$#verts-0] ) {
    return 0;
  } else {
    return 1;
  }
}

# subdivide a surface taking into account if it's a left or right hemisphere.

sub subdivide_mesh {

  my $input = shift;
  my $npoly = shift;
  my $output = shift;

  my $npoly_input = `print_n_polygons $input`;
  chomp( $npoly_input );
  if( !CheckFlipOrientation( $input ) ) {
    &run( "subdivide_polygons", $input, $output, $npoly );
  } else {
    # flip right as left first before subdividing, then flip back.
    &run( "param2xfm", '-clobber', '-scales', -1, 1, 1,
          "${tmpdir}/flip.xfm" );
    my $input_flipped = "${tmpdir}/right_flipped.obj";
    &run( "transform_objects", $input,
          "${tmpdir}/flip.xfm", $input_flipped );
    &run( "subdivide_polygons", $input_flipped, $output, $npoly );
    &run( "transform_objects", $output,
          "${tmpdir}/flip.xfm", $output );  # flip.xfm is its own inverse
    unlink( $input_flipped );
    unlink( "${tmpdir}/flip.xfm" );
  }
}

# Execute a system call.

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