avwconv.cc

/*  fslconv.cc

    Mark Jenkinson, FMRIB Image Analysis Group

    Copyright (C) 2003-2004 University of Oxford  */

/*  CCOPYRIGHT  */

// General convolution and filtering utility

#define _GNU_SOURCE 1
#define POSIX_SOURCE 1

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

using namespace MISCMATHS;
using namespace NEWIMAGE;
using namespace Utilities;


string title="fslconv (Version 1.0)\nCopyright(c) 2003, University of Oxford (Mark Jenkinson)";
string examples="fslconv [options] -i <input image> -k <kernel image> -o <output image>";


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> forcesum(string("--forcesum"), false,
		  string("forces convolution to be done by summation"),
		  false, no_argument);
Option<bool> forcefft(string("--forcefft"), false,
		  string("forces convolution to be done by fft"),
		  false, no_argument);
Option<bool> dilate(string("--dilate"), false,
		  string("perform gray-value dilation (max filter) using kernel as mask"),
		  false, no_argument);
Option<bool> erode(string("--erode"), false,
		  string("perform gray-value erosion (min filter) using kernel as mask"),
		  false, no_argument);
Option<bool> medianfilt(string("--median"), false,
		  string("perform median filtering using kernel as mask"),
		  false, no_argument);
Option<float> boxsize(string("-b,--box"), 0,
		  string("size of box kernel (in mm)"),
		  false, requires_argument);
Option<float> sphere(string("-s,--sphere"), 0,
		  string("radius of spherical kernel (in mm)"),
		  false, requires_argument);
Option<float> gaussian(string("-g,--gaussian"), 0,
		  string("size of Gaussian kernel (sigma in mm)"),
		  false, requires_argument);
Option<string> inname(string("-i,--in"), string(""),
		  string("input filename for base image"),
		  true, requires_argument);
Option<string> kername(string("-k,--kernel"), string(""),
		  string("input filename for kernel"),
		  false, requires_argument);
Option<string> outname(string("-o,--out"), string(""),
		  string("output filename for convolved image"),
		  true, requires_argument);
int nonoptarg;

////////////////////////////////////////////////////////////////////////////

// Convolution Kernel Calculations

volume<float> box_kernel(float boxsize, const volume<float>& vin) {
  float b = boxsize;
  int sx = ((int) floor(b/vin.xdim()/2))*2 + 1;
  int sy = ((int) floor(b/vin.ydim()/2))*2 + 1;
  int sz = ((int) floor(b/vin.zdim()/2))*2 + 1;
  volume<float> vker(sx,sy,sz);
  vker.copyproperties(vin);
  vker = 1.0;
  return vker;
}


volume<float> spherical_kernel(float radius, const volume<float>& vin) {
  int sx = MISCMATHS::round(radius/vin.xdim())*2 + 1;
  int sy = MISCMATHS::round(radius/vin.ydim())*2 + 1;
  int sz = MISCMATHS::round(radius/vin.zdim())*2 + 1;
  volume<float> vker(sx,sy,sz);
  vker.copyproperties(vin);
  vker = 0.0;
  float dx2=Sqr(vin.xdim());
  float dy2=Sqr(vin.ydim());
  float dz2=Sqr(vin.zdim());
  for (int z=-sz/2; z<=sz/2; z++) {
    for (int y=-sy/2; y<=sy/2; y++) {
      for (int x=-sx/2; x<=sx/2; x++) {
	if ((x*x*dx2+y*y*dy2+z*z*dz2)<=Sqr(radius)) { 
	  vker(x+sx/2,y+sy/2,z+sz/2)=1.0; 
	}
      }
    }
  }
  return vker;
}


volume<float> gaussian_kernel(float sigma, const volume<float>& vin) {
  int sx = ((int) ceil(sigma*3/vin.xdim()))*2 + 1;
  int sy = ((int) ceil(sigma*3/vin.ydim()))*2 + 1;
  int sz = ((int) ceil(sigma*3/vin.zdim()))*2 + 1;
  volume<float> vker(sx,sy,sz);
  vker.copyproperties(vin);
  float dx2=Sqr(vin.xdim());
  float dy2=Sqr(vin.ydim());
  float dz2=Sqr(vin.zdim());
  for (int z=-sz/2; z<=sz/2; z++) {
    for (int y=-sy/2; y<=sy/2; y++) {
      for (int x=-sx/2; x<=sx/2; x++) {
	vker(x+sx/2,y+sy/2,z+sz/2)=exp(-(x*x*dx2+y*y*dy2+z*z*dz2)/(2*sigma*sigma));
      }
    }
  }
  return vker;
}

////////////////////////////////////////////////////////////////////////////

// Convolution code

// template <class T, class S>
// int insertpart(volume<T>& v1, const volume<S>& v2) 
// {
//   for (int z=v2.minz(); z<=v2.maxz(); z++) {
//     for (int y=v2.miny(); y<=v2.maxy(); y++) {
//       for (int x=v2.minx(); x<=v2.maxx(); x++) {
// 	v1(x,y,z)=(T) v2(x,y,z);
//       }
//     }
//   }
//   return 0;
// }


// template <class T, class S>
// volume<S> extractpart(const volume<T>& v1, const volume<S>& v2, 
// 		      const volume<S>& kernel) 
// {
//   volume<S> vout=v2;
//   vout = (S) 0.0;
//   int kxoff = (kernel.xsize()-1)/2;
//   int kyoff = (kernel.ysize()-1)/2;
//   int kzoff = (kernel.zsize()-1)/2;
//   for (int z=v2.minz(); z<=v2.maxz(); z++) {
//     for (int y=v2.miny(); y<=v2.maxy(); y++) {
//       for (int x=v2.minx(); x<=v2.maxx(); x++) {
// 	vout(x,y,z)=(S) v1(x+kxoff,y+kyoff,z+kzoff);
//       }
//     }
//   }
//   return vout;
// }


// float fsllog2(float x)
// {
//   // a cygwin annoyance!
//   return log(x)/log(2);
// }


volume<float> efficient_convolve(const volume<float>& vin,
				 const volume<float>& vker)
{
  bool usefft=true;
  if (forcesum.set()) { usefft=false; }
  else if (forcefft.set()) { usefft=true; }
  else {
    // estimate calculation time for the two methods and pick the best
    float offt = 2 * vin.nvoxels() * fsllog2(2 * vin.nvoxels());
    float osum = vin.nvoxels() * vker.nvoxels();
    float scalefactor = 44;  // relative unit operation cost for fft vs sum
    usefft = (osum > offt * scalefactor);
  }
  if (usefft) {
    int sx = Max(vin.xsize(),vker.xsize())*2;
    int sy = Max(vin.ysize(),vker.ysize())*2;
    int sz = Max(vin.zsize(),vker.zsize())*2;
    complexvolume vif, vkf;
    vif.re().reinitialize(sx,sy,sz);
    //vif.re().copyproperties(vin);
    vif.re() = 0.0;
    vif.im() = vif.re();
    vkf = vif;
    insertpart(vif.re(),vin);
    insertpart(vkf.re(),vker);
    fft3(vif);
    fft3(vkf);
    vif *= vkf;
    ifft3(vif);
    return extractpart(vif.re(),vin,vker);
  } else {
    // Direct intensity-based summation method
    return convolve(vin,vker);
  }
}

////////////////////////////////////////////////////////////////////////////

// Main execution stuff

int do_work(int argc, char* argv[]) 
{
  volume<float> vin, vout, vker;
  read_volume(vin,inname.value());

  if (boxsize.set()) { vker = box_kernel(boxsize.value(),vin); }

  if (sphere.set()) { vker = spherical_kernel(sphere.value(),vin); }

  if (gaussian.set()) { vker = gaussian_kernel(gaussian.value(),vin); }

  if (kername.set()) { read_volume(vker,kername.value()); }

  if (verbose.value()) { print_volume_info(vker,"kernel"); }

  if (dilate.value()) {
    vout = morphfilter(vin,vker,"dilate");
  } else if (erode.value()) {
    vout = morphfilter(vin,vker,"erode");
  } else if (medianfilt.value()) {
    vout = morphfilter(vin,vker,"median");
  } else {
    if (kername.unset()) { vker /= vker.sum(); }
    vout = efficient_convolve(vin,vker);
  }

  save_volume(vout,outname.value());
  return 0;
}

////////////////////////////////////////////////////////////////////////////

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

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

  try {
    options.add(inname);
    options.add(outname);
    options.add(kername);
    options.add(gaussian);
    options.add(boxsize);
    options.add(sphere);
    options.add(dilate);
    options.add(erode);
    options.add(medianfilt);
    options.add(forcesum);
    options.add(forcefft);
    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)) ||
	 (kername.unset() && boxsize.unset() && sphere.unset() 
	  && gaussian.unset()) )
      {
	options.usage();
	exit(EXIT_FAILURE);
      }
    
  }  catch(X_OptionError& e) {
    options.usage();
    cerr << endl << e.what() << endl;
    exit(EXIT_FAILURE);
  } catch(std::exception &e) {
    cerr << e.what() << endl;
  } 

  // Call the local functions

  return do_work(argc,argv);
}