fslmeants.cc

/*  fslmeants.cc

    Mark Jenkinson and Matthew Webster, FMRIB Image Analysis Group

    Copyright (C) 2007 University of Oxford  */

/*  CCOPYRIGHT  */


// Creates a mean time series (ignoring zeros) from the input 4D volume
// Saves the results as a column in a text file

#include "newimage/newimageall.h"
#include "miscmaths/miscmaths.h"
#include "utils/options.h"

using namespace NEWIMAGE;
using namespace MISCMATHS;


using namespace Utilities;

// The two strings below specify the title and example usage that is
//  printed out as the help or usage message

string title="fslmeants (Version 1.2)\nCopyright(c) 2004-2009, University of Oxford (Mark Jenkinson, Christian F. Beckmann)\nPrints average timeseries (intensities) to the screen (or saves to a file).\nThe average is taken over all voxels in the mask (or all voxels in the image if no mask is specified).\n";
string examples="fslmeants -i filtered_func_data -o meants.txt -m my_mask\nfslmeants -i filtered_func_data -m my_mask\nfslmeants -i filtered_func_data -c 24 19 10";

// Each (global) object below specificies as option and can be accessed
//  anywhere in this file (since they are global).  The order of the
//  arguments needed is: name(s) of option, default value, help message,
//       whether it is compulsory, whether it requires arguments
// Note that they must also be included in the main() function or they
//  will not be active.

Option<bool> verbose(string("-v,--verbose"), false, 
		     string("switch on diagnostic messages"), 
		     false, no_argument);
Option<bool> help(string("-h,--help"), false,
		  string("display this message"),
		  false, no_argument);
Option<bool> usemm(string("--usemm"), false,
		  string("\tuse mm instead of voxel coordinates (for -c option)"),
		  false, no_argument);
Option<bool> showall(string("--showall"), false,
		  string("show all voxel time series (within mask) instead of averaging"),
		  false, no_argument);
Option<string> inname(string("-i"), string(""),
		  string("~<filename>\tinput 4D image"),
		  true, requires_argument);
Option<string> maskname(string("-m"), string(""),
		  string("~<filename>\tinput 3D mask"),
		  false, requires_argument);
Option<string> outmat(string("-o"), string(""),
		  string("~<filename>\toutput text matrix"),
		  false, requires_argument);
Option<float> coordval(string("-c"), 0.0,
		  string("~<x y z>\trequested spatial coordinate (instead of mask)"), 
		  false, requires_3_arguments);
Option<bool> eig(string("--eig"), false,
		  string("        calculate Eigenvariate(s) instead of mean (output will have 0 mean)"),
		  false, no_argument);
Option<bool> bin_mask(string("--no_bin"), true,
			  string("        do not binarise the mask for calculation of Eigenvariates"),
			  false, no_argument);
Option<int> order(string("--order"), 1,
		  string("        select number of Eigenvariates (default 1)"),
		  false, requires_argument);
Option<bool> transpose(string("--transpose"), false,
		  string("        output results in transpose format (one row per voxel/mean)"),
		  false, no_argument);

int nonoptarg;



int main(int argc,char *argv[])
{

  Tracer tr("main");
  OptionParser options(title, examples);

  options.add(inname);
  options.add(outmat);
  options.add(maskname);
  options.add(coordval);
  options.add(usemm);
  options.add(showall);
  options.add(eig);
  options.add(order);
  options.add(bin_mask);
  options.add(transpose);
  options.add(verbose);
  options.add(help);
  
  nonoptarg = options.parse_command_line(argc, argv);
  
  // line below stops the program if the help was requested or 
  //  a compulsory option was not set
  if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
    {
      options.usage();
      exit(EXIT_FAILURE);
    }
  

  // OK, now do the job ...

  volume4D<float> vin;
  read_volume4D(vin,inname.value());

  volume<float> mask;
  volume4D<float> mask_nonbin;
  if (maskname.set()) {
    read_volume(mask,maskname.value());
  } else {
    mask = vin[0];
    mask = 1.0;
  }


  // override the mask if a coordinate is specified
  if (coordval.set()) {
    mask = vin[0];
    mask = 0.0;
    float x, y, z;
    x = coordval.value(0);
    y = coordval.value(1);
    z = coordval.value(2);
    ColumnVector v(4);
    v << x << y << z << 1.0;
    if (usemm.value()) {
      // convert from mm to newimage voxels
      v = (vin[0].newimagevox2mm_mat()).i() * v;
    } else {
      // convert from nifti voxels (input) to newimage voxels (internal)
      v = vin[0].niftivox2newimagevox_mat() * v;
    }
    x = v(1);  y = v(2);  z = v(3);
    mask(MISCMATHS::round(x),MISCMATHS::round(y),MISCMATHS::round(z)) = 1.0;
  }

  if (!samesize(vin[0],mask)) {
    cerr << "ERROR: Mask and Input volumes have different (x,y,z) size." 
	 << endl;
    return -1;
  }

  mask_nonbin.addvolume(mask);
  mask.binarise(1e-8);  // arbitrary "0" threshold

  if(eig.value()){
    Matrix dat, evecs, scales;
    scales = mask_nonbin.matrix(mask);
    dat = vin.matrix(mask);
    if(!bin_mask.value())
      dat = SP (dat, ones(dat.Nrows(),1) * scales);
    dat = remmean(dat,1);
    
    if (verbose.value()) {
      cout << "Number of voxels used = " << dat.Ncols() << endl;
    }
    
    SymmetricMatrix Corr;
    Corr << dat * dat.t() / dat.Ncols();
    DiagonalMatrix tmpD;
    EigenValues(Corr,tmpD,evecs);	
    evecs = fliplr(evecs.Columns(evecs.Ncols()-order.value()+1 , evecs.Ncols())) * std::sqrt(dat.Nrows());
    
    Matrix res2 = mean(dat,2);
    res2 = res2.Column(1).t() * evecs.Column(1);
    
    if((float)res2.AsScalar()<0)  
      evecs = -1.0 * evecs;
    
    if (transpose.value()) { evecs=evecs.t(); }
    if (outmat.set()) {
      write_ascii_matrix(evecs,outmat.value());
    } else {
      cout << evecs << endl;
    }
    if (transpose.value()) { evecs=evecs.t(); }
  }
  else{	
    int nt = vin.tsize();
    int nvox = 1;
    if (showall.value()) {
      nvox = MISCMATHS::round(mask.sum());
      nt += 3;
    }
    Matrix meants(nt,nvox);
    meants = 0;
    long int num=0;
    
    for (int z=mask.minz(); z<=mask.maxz(); z++) {
      for (int y=mask.miny(); y<=mask.maxy(); y++) {
        for (int x=mask.minx(); x<=mask.maxx(); x++) {
	  if (mask(x,y,z)>0.5) {
	    num++;
	    if (showall.value()) {
	      ColumnVector v(4);
	      v << x << y << z << 1.0;
	      v = (vin[0].niftivox2newimagevox_mat()).i() * v;
	      meants(1,num) = v(1);
	      meants(2,num) = v(2);
	      meants(3,num) = v(3);
	      meants.SubMatrix(4,nt,num,num) = vin.voxelts(x,y,z);
	    } else {
	      meants += vin.voxelts(x,y,z);
	    }
	  }
        }
      }
    }
    
    if (verbose.value()) {
      cout << "Number of voxels used = " << num << endl;
    }
    
    // normalise for number of valid entries if averaging
    if (!showall.value()) {
      if (num>0) meants /= (float) num;
    }
    
    // save the result
    if (transpose.value()) { meants=meants.t(); }
    if (outmat.set()) {
      write_ascii_matrix(meants,outmat.value());
    } else {
      cout << meants << endl;
    }
    if (transpose.value()) { meants=meants.t(); }
  }
  return 0;
}