#include <odindata/data.h>
#include <odindata/fileio.h>
#include <odindata/utils.h>
#include <odindata/statistics.h>
#include <odindata/correlation.h>
#include <odindata/filter.h>

/**
  * \page micalc Basic calculations with medical image data
  * \verbinclude cmdline-utils/micalc.usage
  *
  * \section micalc_examples Some examples how to use micalc:
  *
  * Calculate the difference between two DICOM files 'image1.dcm' 'and image2.dcm' and store it in
  * float-type raw data file 'diff.float':
  *   \verbatim
       micalc -if1 image1.dcm -op '-' -if2 image2.dcm -of diff.float
      \endverbatim
  *
  *
  * Multiply image magnitude of Nifti file 'image.nii' by a factor 10:
  *   \verbatim
       micalc -if1 image.nii -op '*' -in2 10 -of image.nii
      \endverbatim
  *
  *
  * Accumulate all values in file 'image.float' and write result to console:
  *   \verbatim
       micalc -if image.float -op '+'
      \endverbatim
  *
  *
  * Negate values (invert sign) of file image.hdr and write result to neg.hdr
  *   \verbatim
       micalc -if image.hdr -op '-' neg.hdr
      \endverbatim
  *
  *
  * Print statistics about file 'data.asc' to console:
  *   \verbatim
       micalc -if data.asc
      \endverbatim
  *
  *
  * Write time-course mean and fluctuation (relative standard deviation over time) of
  * file 'tcourse.v' to file 'mean.v' and  'stdev.v', respectively:
  *   \verbatim
       micalc -if tcourse.v -mean mean.v -fluct stdev.v
      \endverbatim
  *
  *
  * Calculate histogram in interval (0.0,100.0) with 64 steps of DICOMs in directory cbv:
  *   \verbatim
       micalc -if cbv -hist histogram.asc -histmin 0.0 -histmax 100.0 -histslots 64
      \endverbatim
  *
  *
  * Please note that it might be necessary to quote the operation argument (as shown above)
  * in order to prevent the shell from expanding it.
  */

void usage(const Protocol& prot) {
  FileReadOpts ropts;
  FileWriteOpts wopts;
  FilterChain filter;
  STD_cout << "micalc:  Performs basic mathematics with data sets" << STD_endl;
  STD_cout << "         File formats are automatically identified by their file extension." << STD_endl;
  STD_cout << "         For binary operations with two input files, the protocol of the first input file will be used for the output file." << STD_endl;
  STD_cout << configInfo() << STD_endl;
  STD_cout << STD_endl;
  STD_cout << "micalc can be used in one of the following modes:" << STD_endl;
  STD_cout << "  Binary operation with data sets and/or scalar numbers:" << STD_endl;
  STD_cout << "       micalc [-if1 <input-file1> | -in1 <input-number1>] -op <operation(+,-,*,/,min,max,lcorr,kcorr)> [-if2 <input-file2> | -in2 <input-number1>] -of <output-file>" << STD_endl;
  STD_cout << "  Accumulation/reduction of one data set:" << STD_endl;
  STD_cout << "       micalc -if <input-file> -op <operation(+,*,mean,median,stdev,min,max)>" << STD_endl;
  STD_cout << "  Transformation (magnitude, logarithm, ...) of one data set:" << STD_endl;
  STD_cout << "       micalc -if <input-file> -op <operation(abs,sqrt,log,exp,-,/,acos)> -of <output-file>" << STD_endl;
  STD_cout << "  Statistics of one data set:" << STD_endl;
  STD_cout << "       micalc -if <input-file>" << STD_endl;
  STD_cout << "         [-mean <time-mean-file>] [-stdev <time-stdev-file>] [-fluct <relative-time-stdev-file>] [-tcourse <average time course of all voxels>]" << STD_endl;
  STD_cout << "         [-hist <histogram> -histslots <numofslots> [-histmin <minval>] [-histmax <maxval>] [-rightstairs] [-histfract] [-loghist] [-histslotmap <file>]]" << STD_endl;
  STD_cout << "         [-cog <center-of-gravity-file>]" << STD_endl;

  STD_cout << "Extra options:" << STD_endl;
  STD_cout << "\t-mask <Binary mask file: voxels with zeroes are discarded>" << STD_endl;
  STD_cout << "Protocol options (overrides parameters in input data):" << STD_endl;
  STD_cout << prot.get_cmdline_usage("\t");
  STD_cout << "File read options (will be applied to input data before processing):" << STD_endl;
  STD_cout << ropts.get_cmdline_usage("\t");
  STD_cout << "File write options (will be applied to output data after processing):" << STD_endl;
  STD_cout << wopts.get_cmdline_usage("\t");
  STD_cout << "Filters/processing applied to input file(s) and mask, applied in the order they appear on the command line:" << STD_endl;
  STD_cout << filter.get_cmdline_usage("\t");
  STD_cout << "Other options:" << STD_endl;
  STD_cout << "\t" << LogBase::get_usage() << STD_endl;
  STD_cout << "\t" << helpUsage() << STD_endl;
  STD_cout << "Supported file extensions(formats):" << STD_endl;
  STD_cout << FileIO::autoformats_str("\t") << STD_endl;
}


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


bool read_infile(const STD_string& filename, Data<float,4>& data, const FileReadOpts& ropts, Protocol& prot, const FilterChain& filterchain, ProgressMeter* progmeter) {
  Log<OdinData> odinlog("micalc","read_infile");

  Data<float,4> datain;

  if(datain.autoread(filename,ropts,&prot, progmeter)<0) return false;
  if(!filterchain.apply(prot,datain)) return false;
  data.resize(datain.shape());

  data=datain; // In case storage order is non-default

  return true;
}


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


int main(int argc, char* argv[]) {
  LogBase::set_log_levels(argc,argv);

  Log<OdinData> odinlog("micalc","main");

  Range all=Range::all();

  FileReadOpts  ropts;
  FileWriteOpts wopts;
  Protocol prot;

  if(hasHelpOption(argc,argv)) {usage(prot); return 0;}

  // set defaults
  prot.seqpars.set_MatrixSize(readDirection,1);
  prot.seqpars.set_MatrixSize(phaseDirection,1);
  prot.seqpars.set_MatrixSize(sliceDirection,1);

  char optval[ODIN_MAXCHAR];
  STD_string operation;
  STD_string infile1;
  STD_string infile2;
  STD_string outfile;
  STD_string maskfile;

  float innumber1=0.0;
  float innumber2=0.0;


  bool i1flag=false;
  bool i2flag=false;
  bool offlag=false;

  bool ifflag=false;
  bool opflag=false;

  bool binaryop=false;
  bool accumulate=false;
  bool transform=false;
  bool stateval=false;

  // command line parsing, parse before filters to remove options
  if(getCommandlineOption(argc,argv,"-op",optval,ODIN_MAXCHAR)) {operation=optval; opflag=true;}

  if(getCommandlineOption(argc,argv,"-if1",optval,ODIN_MAXCHAR)) {infile1=optval; i1flag=true;}
  if(getCommandlineOption(argc,argv,"-if2",optval,ODIN_MAXCHAR)) {infile2=optval; i2flag=true;}
  if(getCommandlineOption(argc,argv,"-if",optval,ODIN_MAXCHAR)) {infile1=optval; ifflag=true;}
  if(getCommandlineOption(argc,argv,"-in1",optval,ODIN_MAXCHAR)) {innumber1=atof(optval); i1flag=true;}
  if(getCommandlineOption(argc,argv,"-in2",optval,ODIN_MAXCHAR)) {innumber2=atof(optval); i2flag=true;}

  if(getCommandlineOption(argc,argv,"-of",optval,ODIN_MAXCHAR)) {outfile=optval; offlag=true;}

  if(getCommandlineOption(argc,argv,"-mask",optval,ODIN_MAXCHAR)) {maskfile=optval;}

  STD_string histfile;
  if(getCommandlineOption(argc,argv,"-hist",optval,ODIN_MAXCHAR)) histfile=optval;

  STD_string cogfile;
  if(getCommandlineOption(argc,argv,"-cog",optval,ODIN_MAXCHAR)) cogfile=optval;

  STD_string meanfile;
  STD_string fluctfile;
  STD_string stdevfile;
  STD_string tcoursefile;
  STD_string histslots;
  STD_string histmin;
  STD_string histmax;
  STD_string histslotmapfile;
  if(getCommandlineOption(argc,argv,"-mean",optval,ODIN_MAXCHAR))  meanfile=optval;
  if(getCommandlineOption(argc,argv,"-fluct",optval,ODIN_MAXCHAR)) fluctfile=optval;
  if(getCommandlineOption(argc,argv,"-stdev",optval,ODIN_MAXCHAR)) stdevfile=optval;
  if(getCommandlineOption(argc,argv,"-tcourse",optval,ODIN_MAXCHAR)) tcoursefile=optval;
  if(getCommandlineOption(argc,argv,"-histslots",optval,ODIN_MAXCHAR)) histslots=optval;
  if(getCommandlineOption(argc,argv,"-histmin",optval,ODIN_MAXCHAR)) histmin=optval;
  if(getCommandlineOption(argc,argv,"-histmax",optval,ODIN_MAXCHAR)) histmax=optval;
  bool rightstairs=isCommandlineOption(argc,argv,"-rightstairs");
  bool histfract=isCommandlineOption(argc,argv,"-histfract");
  bool loghist=isCommandlineOption(argc,argv,"-loghist");
  if(getCommandlineOption(argc,argv,"-histslotmap",optval,ODIN_MAXCHAR)) histslotmapfile=optval;



  // Parse before filters to remove options
  ropts.parse_cmdline_options(argc,argv);
  wopts.parse_cmdline_options(argc,argv);
  prot.parse_cmdline_options(argc,argv);


  ProgressDisplayConsole display;
  ProgressMeter progmeter(display);


  FilterChain filterchain(argc,argv);


  if(i1flag && i2flag && opflag) binaryop=true;
  if(ifflag && opflag) {
    if(offlag) transform=true;
    else accumulate=true;
  }
  if(ifflag && (!opflag)) stateval=true;

  if( !( binaryop || accumulate || transform || stateval ) ) {usage(prot); return 0;}

  if(accumulate || stateval) FileIO::set_trace_status(false);

  if(binaryop && infile1=="" && infile2=="") {
    ODINLOG(odinlog,errorLog) << "Specify at least one file" << STD_endl;
    return -1;
  }

  Data<float,4> data1;
  Data<float,4> data2;

  if(infile1!="") {
    if(!read_infile(infile1, data1, ropts, prot, filterchain,&progmeter)) return -1;
  }

  if(binaryop && infile2!="") {
    Protocol protcopy(prot); // dispensable protocol
    Protocol* prot2p=&protcopy;
    if(infile1=="") prot2p=&prot; // use only if 1st operand is number
    if(!read_infile(infile2, data2, ropts, *prot2p, filterchain,&progmeter)) return -1;
  }

  TinyVector<int,4> outshape;


  if(binaryop || transform) {

    bool has_shape1=sum(data1.shape());
    bool has_shape2=sum(data2.shape());

    if( has_shape1 && has_shape2 && !same_shape(data1,data2) ) {
      ODINLOG(odinlog,errorLog) << "Shape mismatch:  data1shape/data2shape=" << data1.shape() << "/" << data2.shape() << STD_endl;
      return -1;
    } else outshape=data1.shape();

    if(has_shape1 && !has_shape2) {
      outshape=data1.shape();
      data2.resize(outshape);
      data2=innumber2;
    }
    if(!has_shape1 && has_shape2) {
      outshape=data2.shape();
      data1.resize(outshape);
      data1=innumber1;
    }
    if(!has_shape1 && !has_shape2) {
      ODINLOG(odinlog,errorLog) << "Specify at least one valid file" << STD_endl;
      return -1;
    }

  }


  Data<float,4>* maskptr=0;
  Data<float,4> maskdata;
  if(maskfile!="") {
    Data<float,4> inmaskdata;
    Protocol maskprot;
    if(!read_infile(maskfile, inmaskdata, ropts, maskprot, filterchain,&progmeter)) return -1;
//    if(inmaskdata.autoread(maskfile,ropts)<0) return -1;
    if( inmaskdata.extent(0) && same_shape(data1, inmaskdata, TinyVector<int,4>(0,1,1,1) ) ) {
      maskdata.resize(data1.shape());
      maskdata=0.0;
      for(unsigned int i=0; i<data1.numElements(); i++) {
        TinyVector<int,4> index=data1.create_index(i);
        TinyVector<int,4> maskindex=index; maskindex(0)=0;
        if(inmaskdata(maskindex)!=0.0) maskdata(index)=1.0;
      }
      maskptr=&maskdata;
    } else {
      ODINLOG(odinlog,errorLog) << "Shape mismatch:  datashape/maskshape=" << data1.shape() << "/" << inmaskdata.shape() << STD_endl;
      return -1;
    }
  }


  prot.system.set_data_type(TypeTraits::type2label(float(0))); // Store results in float

  if(binaryop) {
    if(maskptr) {
      ODINLOG(odinlog,normalDebug) << "data1/data2/maskdata" << data1.shape() << "/" << data2.shape() << "/" << maskdata.shape() << STD_endl;
      data1=data1*maskdata;
      data2=data2*maskdata;
    }

    if(operation=="lcorr" || operation=="kcorr") {
      int npts=data1.size();
      if(maskptr) npts=int(sum(maskdata)+0.5);
      Array<float,1> vec1(npts);
      Array<float,1> vec2(npts);
      int vecindex=0;
      for(unsigned int i=0; i<data1.numElements(); i++) {
        TinyVector<int,4> index=data1.create_index(i);
        bool include=true;
        if( maskptr && maskdata(index)==0.0 ) include=false;
        if(include) {
          vec1(vecindex)=data1(index);
          vec2(vecindex)=data2(index);
          vecindex++;
        }
      }
      correlationResult coresult;
      if(operation=="lcorr") coresult=correlation(vec1, vec2);
      if(operation=="kcorr") coresult=kendall(vec1, vec2);
      STD_cout << "r/p/z(" << npts << ")=" << coresult.r << "/" << coresult.p << "/"  << coresult.z << STD_endl;
    } else {
      Data<float,4> result(outshape);
      ODINLOG(odinlog,infoLog) << "calculating " << data1.shape() << " " << operation << " " << data2.shape() << " ..." << STD_endl;
      bool validop=false;
      if(operation=="+") {validop=true; result=data1+data2;}
      if(operation=="-") {validop=true; result=data1-data2;}
      if(operation=="*") {validop=true; result=data1*data2;}
      if(operation=="/") {validop=true; result=secureDivision(data1,data2);}
      if(operation=="min" || operation=="max") {
        Data<float,4> minval(outshape);
        Data<float,4> maxval(outshape);
        for(unsigned int i=0; i<data1.numElements(); i++) {
          TinyVector<int,4> index=data1.create_index(i);
          minval(index)=STD_min(data1(index), data2(index));
          maxval(index)=STD_max(data1(index), data2(index));
        }
        if(operation=="min") result=minval;
        else                 result=maxval;
        validop=true;
      }
      if(!validop) {
        ODINLOG(odinlog,errorLog) << "invalid binary operation: " << operation << STD_endl;
        return -1;
      }
      if(result.autowrite(outfile,wopts,&prot)<0) return -1;
    }
    return 0;
  }

  if(accumulate) {
    if(maskptr) {
      data1=data1*maskdata;
    }
    double accuresult=0.0;
    bool validop=false;
    if(operation=="+") {validop=true; accuresult=sum(data1);}
    if(operation=="*") {validop=true; accuresult=product(data1);}
    if(operation=="mean") {validop=true; accuresult=statistics(data1).mean;}
    if(operation=="median") {validop=true; accuresult=median(data1);}
    if(operation=="stdev") {validop=true; accuresult=statistics(data1).stdev;}
    if(operation=="min") {validop=true; accuresult=min(data1);}
    if(operation=="max") {validop=true; accuresult=max(data1);}
    if(!validop) {
      ODINLOG(odinlog,errorLog) << "invalid accumulation operator: " << operation << STD_endl;
      return -1;
    }
    STD_cout << accuresult << STD_endl;
  }

  if(transform) {
    if(maskptr) {
      data1=data1*maskdata;
    }
    Data<float,4> result(outshape);
    ODINLOG(odinlog,infoLog) << "Calculating " << operation << " " << data1.shape() << " ..." << STD_endl;
    bool validop=false;
    if(operation=="abs")  {validop=true; result=fabs(data1);}
    if(operation=="sqrt") {validop=true; result=sqrt(data1);}
    if(operation=="log")  {validop=true; result=where(Array<float,4>(data1)>0.0, Array<float,4>(log(data1)), float(0.0));}
    if(operation=="exp")  {validop=true; result=exp(data1);}
    if(operation=="-")    {validop=true; result=-data1;}
    if(operation=="/")    {validop=true; result=where(Array<float,4>(data1)!=0.0, Array<float,4>(1.0/data1), float(0.0));}
    if(operation=="acos") {validop=true; result=where(Array<float,4>(data1)>=-1.0 && Array<float,4>(data1)<=1.0, Array<float,4>(acos(data1)), float(0.0));}
    if(!validop) {
      ODINLOG(odinlog,errorLog) << "invalid unary operation: " << operation << STD_endl;
      return -1;
    }
    if(result.autowrite(outfile,wopts,&prot)<0) return -1;
    return 0;
  }

  if(stateval) {

    if(histfile!="") {

      int nvals=data1.numElements();
      if(maskptr) {
        nvals=0;
        for(unsigned int i=0; i<maskdata.numElements(); i++) {
          TinyVector<int,4> index=maskdata.create_index(i);
          if( maskdata(index)>0.0 ) {
            nvals++;
          }
        }
      }


      if(histslots=="") {usage(prot); return -1;}
      int nslots=atoi(histslots.c_str());


      statisticResult stats=statistics(data1,maskptr);
      float minval=stats.min;
      if(histmin!="") minval=atof(histmin.c_str());
      float maxval=stats.max;
      if(histmax!="") maxval=atof(histmax.c_str());
      ODINLOG(odinlog,infoLog) << "calculating histogram with nslots/minval/maxval=" << nslots <<  "/" << minval << "/" << maxval << STD_endl;

      if(minval>=maxval) {
        ODINLOG(odinlog,errorLog) << "histmax must be larger than histmin" << STD_endl;
        return -1;
      }

      float stairoffset=0.5; // centered stairs
      if(rightstairs) stairoffset=0.0; // assume rightstairs when plotting (xmgrace)

      Data<float,1> xvals(nslots);
      Data<float,1> yvals(nslots); yvals=0.0;


      Data<float,4> histslotmap(data1.shape()); histslotmap=0.0;


//      if(minval>0.0 && maxval>0.0) { // Use log histogram
      if(loghist) { // Use log histogram

        if( minval<=0.0 || maxval<=0.0 ) {
          ODINLOG(odinlog,errorLog) << "positive minval/maxval required for log histogram, current minval/maxval=" << minval << "/" << maxval << STD_endl;
          return -1;
        }

        float minpow=log10(minval);
        float maxpow=log10(maxval);
        float powrange=maxpow-minpow;
        float step=secureDivision(powrange,nslots);
        ODINLOG(odinlog,normalDebug) << "minpow/maxpow/powrange/step=" << minpow << "/" << maxpow << "/" << powrange << "/" << step << STD_endl; 

        for(int i=0; i<nslots; i++) xvals(i)=pow(float(10.0),float(minpow+(float(i)+stairoffset)*step));

        for(unsigned int i=0; i<data1.numElements(); i++) {
          TinyVector<int,4> index=data1.create_index(i);
          if( !(maskptr && maskdata(index)==0.0) ) {
            float pow=log10(data1(index));
            int slot=int(secureDivision(pow-minpow,powrange)*nslots);
            ODINLOG(odinlog,normalDebug) << "pow/slot=" << pow << "/" << slot << STD_endl; 
            if(slot>=0 && slot<nslots) {
              yvals(slot)++;
              histslotmap(index)=slot+stairoffset;
            }
          }
        }


      } else { // no log histogram

        float range=maxval-minval;
        float step=secureDivision(range,nslots);
        ODINLOG(odinlog,normalDebug) << "range/step=" << range << "/" << step << STD_endl; 

        for(int i=0; i<nslots; i++) xvals(i)=minval+(float(i)+stairoffset)*step;

        for(unsigned int i=0; i<data1.numElements(); i++) {
          TinyVector<int,4> index=data1.create_index(i);
          if( !(maskptr && maskdata(index)==0.0) ) {
            float val=data1(index);
            int slot=int(secureDivision(val-minval,range)*nslots);
            ODINLOG(odinlog,normalDebug) << "val/slot=" << val << "/" << slot << STD_endl; 
            if(slot>=0 && slot<nslots) {
              yvals(slot)++;
              histslotmap(index)=slot+stairoffset;
            }
          }
        }
      }

      if(histfract) yvals/=sum(yvals); // Each slot reflects the fraction of total points
      yvals.write_asc_file(histfile, xvals);

      if(histslotmapfile!="") if(histslotmap.autowrite(histslotmapfile,wopts,&prot)<0) return -1;

    }


    if(meanfile!="" || fluctfile!="" || stdevfile!="" || tcoursefile!="" || cogfile!="") {
      TinyVector<int,4> inshape=data1.shape();
      outshape=inshape;
      ODINLOG(odinlog,infoLog) << "calculating general statistics for shape " << outshape <<  " ..." << STD_endl;
      outshape(0)=1;
      Data<float,4> meanresult(outshape); meanresult=0.0;
      Data<float,4> fluctresult(outshape); fluctresult=0.0;
      Data<float,4> stdevresult(outshape); stdevresult=0.0;
      Data<float,1> tcourseresult(data1.extent(0)); tcourseresult=0.0;

      farray cogsum(inshape(0), 3);
      fvector masssum(inshape(0));

      int maskpoints=0;
      for(unsigned int i=0; i<meanresult.numElements(); i++) {
        TinyVector<int,4> index=meanresult.create_index(i);
        if( !(maskptr && maskdata(index)==0.0) ) {
          statisticResult stats=statistics(data1(all,index(1),index(2),index(3)));
          meanresult(index)= stats.mean;
          fluctresult(index)=secureDivision(stats.stdev,stats.mean);
          stdevresult(index)=stats.stdev;
          tcourseresult(all)+=data1(all,index(1),index(2),index(3));
          ODINLOG(odinlog,normalDebug) << "mean/stdev(" << index << ")=" << stats.mean << "/" << stats.stdev << STD_endl;

          for(int irep=0; irep<inshape(0); irep++) {
            float mass=data1(irep,index(1),index(2),index(3));
            masssum[irep]+=mass;
            for(int idim=0; idim<3; idim++) {
              cogsum(irep,idim)+=mass*index(1+idim);
            }
          }

          maskpoints++;
        }
      }
      ODINLOG(odinlog,infoLog) << "maskpoints=" << maskpoints << STD_endl;
      tcourseresult/=float(maskpoints);

      if(meanfile!="")  if(meanresult.autowrite(meanfile,wopts,&prot)<0) return -1;
      if(fluctfile!="") if(fluctresult.autowrite(fluctfile,wopts,&prot)<0) return -1;
      if(stdevfile!="") if(stdevresult.autowrite(stdevfile,wopts,&prot)<0) return -1;
      if(tcoursefile!="") if(tcourseresult.autowrite(tcoursefile,wopts,&prot)<0) return -1;

      if(cogfile!="") {
        STD_string cogstr;
        for(int irep=0; irep<inshape(timeDim); irep++) {
          cogstr+=ftos(cogsum(irep, 0)/masssum[irep])+" "+ftos(cogsum(irep, 1)/masssum[irep])+" "+ftos(cogsum(irep, 2)/masssum[irep])+"\n";
        }
        write(cogstr, cogfile);
      }


    } else {
      STD_cout << statistics(data1,maskptr) << STD_endl;
    }
//    STD_cout << "median = " << median(data1) << STD_endl;



  }

  return 0;
}