/*  surface_norm.cc

    Mark Jenkinson, FMRIB Image Analysis Group

    Copyright (C) 2004 University of Oxford  */

/*  CCOPYRIGHT  */

// Calculates the surface normals for a mask, using a smoothed
//  gradient calculation (all non-surface points get zero ouput)

#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;

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

string title="surface_norm (Version 1.0)\nCopyright(c) 2004, University of Oxford (Mark Jenkinson)";
string examples="surface_norm [options] -i <input mask image> -s <smoothing in mm> -o <output surface normal image>";

// 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<float>  smoothing(string("-s"),3.0,
			 string("smoothing in mm (default is 3.0)"),
			 false, requires_argument);
Option<string> inname(string("-i"), string(""),
		      string("input mask filename"),
		      true, requires_argument);
Option<string> outname(string("-o"), string(""),
		       string("output surface normal filename"),
		       true, requires_argument);
int nonoptarg;

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

// Local functions

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

  // force input to be a binary mask
  vin.binarise(0.5);
  if (verbose.value()) print_info(vin,"vin");

  // smooth the mask volume and calculate the gradient
  if (verbose.value()) { cerr << "Performing smoothing" << endl; }
  { 
    volume<float> smoothv;
    smoothv = smooth(vin,smoothing.value());
    if (verbose.value()) print_info(smoothv,"smoothv");
    if (verbose.value()) { cerr << "Calculating the gradient" << endl; }
    gradient(smoothv,snorm);
  }

  if (verbose.value()) print_info(snorm,"snorm");

  bool atedge;
  if (verbose.value()) { cerr << "Normalising gradients" << endl; }
  // zero all non-surface values and normalise vectors at the surface
  for (int z=vin.minz(); z<=vin.maxz(); z++) {
    for (int y=vin.miny(); y<=vin.maxy(); y++) {
      for (int x=vin.minx(); x<=vin.maxx(); x++) {
	atedge = false;
	if ( (vin(x,y,z)>0.5) ) {
	  if (vin(x,y,z-1)<0.5) atedge=true;
	  else { 
	    if (vin(x,y-1,z)<0.5) atedge=true;
	    else { 
	      if (vin(x-1,y,z)<0.5) atedge=true;
	      else {
		if (vin(x+1,y,z)<0.5) atedge=true;
		else {
		  if (vin(x,y+1,z)<0.5) atedge=true;
		  else {
		    if (vin(x,y,z+1)<0.5) atedge=true;
		  }
		}
	      }
	    }
	  }
	}
	if (atedge) {
	  float norm;
	  norm = 1.0 / sqrt(Sqr(snorm(x,y,z,0)) + Sqr(snorm(x,y,z,1)) 
			    + Sqr(snorm(x,y,z,2)) );
	  snorm(x,y,z,0) *= norm;
	  snorm(x,y,z,1) *= norm;
	  snorm(x,y,z,2) *= norm;
	} else {
	  snorm(x,y,z,0) = 0.0;
	  snorm(x,y,z,1) = 0.0;
	  snorm(x,y,z,2) = 0.0;
	}
      }
    }
    if (verbose.value()) { cerr << "."; }
  }
  if (verbose.value()) { cerr << endl; }

  // save the results
  save_volume4D(snorm,outname.value());

  return 0;
}

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

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

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

  try {
    // must include all wanted options here (the order determines how
    //  the help message is printed)
    options.add(inname);
    options.add(outname);
    options.add(smoothing);
    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);
      }
    
  }  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);
}