[MINC-users] program to compute determinants of xfm files?

[Claude LEPAGE]

>
> Does such a thing exist?  It'd be easy to write one for the case of a file containing a single linear transformation ... just wondering if there's possibly already something around that might handle more general cases?
>

I wrote some little script a long time ago to compute the magnitude of
the Jacobian from a non-linear xfm transformation. I think it's correct.
Here it is, cut-and-paste (in perl):

#! /usr/bin/env perl
#
# Evaluate the magnitude of the Jacobian from a non-linear transformation.

#
# Note: All of this can be replaced by mincblob -determinant on the grid
#       file after resampling it in z,y,x order. You must add 1 to the
#       mincblob result.
#

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

# --- set the help & usage strings ---
my $help = <<HELP;
Required parameters:
  nl_grid_0.mnc      : displacement file from a non-linear transform
  out.mnc            : output for magnitude of Jacobian transform
HELP

my $usage = <<USAGE;
Usage: jacobian.pl grid.mnc out.mnc [options]
       jacobian.pl -help to list options

USAGE

die "$usage\n" unless @ARGV == 2;

my $grid = shift;
my $out = shift;

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

# Extract x-, y-, z- displacements from grid file.

my $gridX = "${tmpdir}/gridX.mnc";
my $gridY = "${tmpdir}/gridY.mnc";
my $gridZ = "${tmpdir}/gridZ.mnc";
&run( 'mincreshape', '-quiet', '-clobber', '-dimrange', 'vector_dimension=0',
      $grid, $gridX );
&run( 'mincreshape', '-quiet', '-clobber', '-dimrange', 'vector_dimension=1',
      $grid, $gridY );
&run( 'mincreshape', '-quiet', '-clobber', '-dimrange', 'vector_dimension=2',
      $grid, $gridZ );

# Define kernels for gradients.

my $kerndx="${tmpdir}/dx.kern";
my $kerndy="${tmpdir}/dy.kern";
my $kerndz="${tmpdir}/dz.kern";

&CreateKernelFiles;

my $dudx = "${tmpdir}/dudx.mnc";
my $dudy = "${tmpdir}/dudy.mnc";
my $dudz = "${tmpdir}/dudz.mnc";
my $dvdx = "${tmpdir}/dvdx.mnc";
my $dvdy = "${tmpdir}/dvdy.mnc";
my $dvdz = "${tmpdir}/dvdz.mnc";
my $dwdx = "${tmpdir}/dwdx.mnc";
my $dwdy = "${tmpdir}/dwdy.mnc";
my $dwdz = "${tmpdir}/dwdz.mnc";

# compute gradient volume

&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndx, $gridX, $dudx );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndy, $gridX, $dudy );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndz, $gridX, $dudz );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndx, $gridY, $dvdx );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndy, $gridY, $dvdy );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndz, $gridY, $dvdz );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndx, $gridZ, $dwdx );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndy, $gridZ, $dwdy );
&run( "mincmorph", "-clobber", "-convolve", "-kernel", $kerndz, $gridZ, $dwdz );

unlink( $kerndx );
unlink( $kerndy );
unlink( $kerndz );
unlink( $gridX );
unlink( $gridY );
unlink( $gridZ );

my $dx = `mincinfo -attvalue xspace:step $grid`; chomp( $dx ); $dx += 0;
my $dy = `mincinfo -attvalue yspace:step $grid`; chomp( $dy ); $dy += 0;
my $dz = `mincinfo -attvalue zspace:step $grid`; chomp( $dz ); $dz += 0;

my $expr = "(($dx+A[0])*(($dy+A[4])*($dz+A[8])-A[5]*A[7])+" .
            "A[1]*(A[5]*A[6]-A[3]*($dz+A[8]))+" .
            "A[2]*(A[3]*A[7]-A[6]*($dy+A[4])))/($dx*$dy*$dz)";

&run( 'minccalc', '-quiet', '-clobber', '-expression', $expr,
      $dudx, $dudy, $dudz, $dvdx, $dvdy, $dvdz, $dwdx, $dwdy, $dwdz,
      $out );

unlink( $dudx );
unlink( $dudy );
unlink( $dudz );
unlink( $dvdx );
unlink( $dvdy );
unlink( $dvdz );
unlink( $dwdx );
unlink( $dwdy );
unlink( $dwdz );

# -- done --

#Create the kernel files used in taking the image derivative

sub CreateKernelFiles {

  open DXFILE, "> $kerndx";
  print DXFILE "MNI Morphology Kernel File\n";
  print DXFILE "Kernel_Type = Normal_Kernel;\n";
  print DXFILE "Kernel =\n";
  print DXFILE "  -1.0  0.0  0.0  0.0  0.0     -0.5\n";
  print DXFILE "   1.0  0.0  0.0  0.0  0.0      0.5;\n";
  close DXFILE;

  open DYFILE, "> $kerndy";
  print DYFILE "MNI Morphology Kernel File\n";
  print DYFILE "Kernel_Type = Normal_Kernel;\n";
  print DYFILE "Kernel =\n";
  print DYFILE "   0.0 -1.0  0.0  0.0  0.0     -0.5\n";
  print DYFILE "   0.0  1.0  0.0  0.0  0.0      0.5;\n";
  close DYFILE;

  open DZFILE, "> $kerndz";
  print DZFILE "MNI Morphology Kernel File\n";
  print DZFILE "Kernel_Type = Normal_Kernel;\n";
  print DZFILE "Kernel =\n";
  print DZFILE "   0.0  0.0 -1.0  0.0  0.0     -0.5\n";
  print DZFILE "   0.0  0.0  1.0  0.0  0.0      0.5;\n";
  close DZFILE;

}


# Execute a system call.

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