//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//
#include "StackProcessAnalyses.hpp"

#include <QFile>
#include <QFileDialog>
#include <QTextStream>
#include "Misc.hpp"
#include "Triangulate.hpp" // calcNearestPointOnBezierLine
#include "GraphUtils.hpp" // shiftMap

namespace mgx
{
  bool ComputeVolume::run(const Store* input, QString filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
      return setErrorMessage(QString("Error, cannot open file '%1' for writing").arg(filename));
    std::vector<unsigned long> volumes(1 << 16);   // List of volumes per label
    const Stack* stk = input->stack();
    float dv = stk->step().x() * stk->step().y() * stk->step().z();   // volume of a voxel
    const HVecUS& data = input->data();
    for(uint i = 0; i < data.size(); ++i) {
      ushort label = data[i];
      ++volumes[label];
    }

    QTextStream ts(&file);
    ts << QString("Label,Volume (%1)").arg(UM3) << endl;
    for(uint lab = 1; lab < 1 << 16; ++lab) {
      if(volumes[lab] > 0)
        ts << lab << "," << (volumes[lab] * dv) << endl;
    }
    file.close();
    return true;
  }

  REGISTER_PROCESS(ComputeVolume);


   bool StackVoxelCount::run(Stack* s1, Store* store1, QString filename)
  {

    const HVecUS& data = store1->data();
    Point3i imgSize(s1->size());

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    std::map<int, std::map<int, int> > occuranceMap;
    typedef std::pair<int, std::map<int, int> > IntMapP;
    std::set<int> voxValues;

    //#pragma omp parallel for schedule(guided)
    for(int z=0; z<zDim; z++){
      std::cout << "z " << z << std::endl;
      for(int x=0; x<xDim; x++){
        for(int y=0; y<yDim; y++){
          Point3i pOut (x,y,z);
          occuranceMap[z][data[s1->offset(pOut)]]++;
         // #pragma omp critical
          voxValues.insert(data[s1->offset(pOut)]);
        }
      }
    }

    std::map<int, int> voxValueIdxMap;
    std::vector<int> voxValueVec;

    forall(int vox, voxValues){
      voxValueVec.push_back(vox);
      voxValueIdxMap[vox] = voxValueVec.size()-1;
    }

     // file for root cell data
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);

    out << "z-layer";
    forall(int vox, voxValueVec){
      out << "," << vox;
    }
    out << endl;

    forall(IntMapP p, occuranceMap){
      out << p.first;

      forall(int vox, voxValueVec){
        out << "," << p.second[vox];
      }

      out << endl;
    }

    return true;
  }
  REGISTER_PROCESS(StackVoxelCount);

 bool StackCircularHistogram::initialize(QWidget* parent)
  {

    QString fileName = parm("Filename");
    if(parent)
      fileName = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(fileName.isEmpty())
      return false;
    if(!fileName.endsWith(".csv", Qt::CaseInsensitive))
      fileName += ".csv";
    setParm("Filename", fileName);;
    return true;
  }


   bool StackCircularHistogram::run(Stack* s1, Store* store1, Store* store2, QString centralAxis, int binNumber, double thresholdValue, double minDis, double maxDis, QString align, bool weightVol, QString filename, bool writeStack)
  {
    const HVecUS& data = store1->data();
    HVecUS& data2 = store2->data();
    Point3i imgSize(s1->size());

    Matrix4d rotMatrixS1;
    s1->getFrame().getMatrix(rotMatrixS1.data());

    // fix the rotations
    Matrix4d mGLRotMat = transpose(rotMatrixS1);

    Point3d cAxis(0,0,1);
    Point3d refDir(1,0,0);
    bool zAxis = true, yAxis = false, xAxis = false;
    if(centralAxis == "X"){
      cAxis = Point3d(1,0,0);
      refDir = Point3d(0,0,1);
      zAxis = false;
      yAxis = false;
      xAxis = true;
    }
    else if(centralAxis == "Y"){
      cAxis = Point3d(0,1,0);
      refDir = Point3d(1,0,0);
      zAxis = false;
      yAxis = true;
      xAxis = false;
    }

    double binSize = 360. / (double)binNumber;

    std::map<int, double> histogramSumMap;
    std::map<int, double> histogramMaxMap;

    std::map<int,int> newKeyOldKeyMap;

    for(int i = 0; i < binNumber; i++){
      histogramSumMap[i] = 0;
      histogramMaxMap[i] = 0;
      newKeyOldKeyMap[i] = i;
    }

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    Point3f voxelSizes = s1->step();

    double voxelVolume = voxelSizes.x() * voxelSizes.y() * voxelSizes.z();

    Point3i lower(0,0,0);
    Point3i upper(xDim-1, yDim-1, zDim-1);

    Point3d midpoint = Point3d(s1->imageToWorld((lower+upper)/2.));

    Point3d refDir2 = cAxis ^ refDir;

    double totalSignal = 0;

    double radConversion = 180./3.14159;
    double valueThreshold = thresholdValue*655.35; // * 16bit / 100

    //#pragma omp parallel for schedule(guided)
    for(int z=0; z<zDim; z++){
      for(int x=0; x<xDim; x++){
        for(int y=0; y<yDim; y++){
          Point3i ps(x,y,z);
          data2[s1->offset(ps)]=0;
          Point3d p = Point3d(s1->imageToWorld(ps));
          if(data[s1->offset(ps)] < valueThreshold) continue;
          totalSignal += data[s1->offset(ps)];

          p = multMatrix4Point3(mGLRotMat,p); // correct rotations
          // calculate angle
          Point3d pointPlane = projectPointOnPlane(p, midpoint, cAxis);

          if(maxDis > 0 or minDis>0){
            double dis;
            if(zAxis) dis = pointPlane.x()*pointPlane.x() + pointPlane.y()*pointPlane.y();
            else if(yAxis) dis = pointPlane.x()*pointPlane.x() + pointPlane.z()*pointPlane.z();
            else if(xAxis) dis = pointPlane.z()*pointPlane.z() + pointPlane.y()*pointPlane.y();
            dis = sqrt(dis);
            if(dis > maxDis and maxDis > 0) continue;
            if(dis < minDis and minDis > 0) continue;
          }



          Point3d v1 = pointPlane - midpoint;
          double scalar = v1 * refDir;
          double angle = acos(scalar / norm(v1))*radConversion;
          if(v1 * refDir2 < 0) angle = 360-angle;

          //if(angle < 0 or angle > 360) std::cout << "weird angle? " << angle << "/" << scalar << std::endl;
          // find bin
          int bin = angle / binSize;
          if(bin>=binNumber) bin = binNumber-1;
          if(bin<0) bin = 0;
          // add signal to map
          double weightedValue = voxelVolume*data[s1->offset(ps)];
          if(!weightVol) weightedValue = data[s1->offset(ps)];
          histogramSumMap[bin] += weightedValue;
          histogramMaxMap[bin] = std::max(histogramMaxMap[bin],weightedValue);
          if(writeStack){
            data2[s1->offset(ps)]=bin;
          }
        }
      }
    }

    int maxKey = 0;

    if(align == "Align at Max of Signal Sum"){
      maxKey = getKeyOfMaxValue(histogramSumMap);
      shiftMap(histogramSumMap, maxKey);
      newKeyOldKeyMap = shiftMap(histogramMaxMap, maxKey);
      std::cout << "alignsum max at bin:" << maxKey << std::endl;
    } else if(align == "Align at Signal Max"){
      maxKey = getKeyOfMaxValue(histogramMaxMap);
      shiftMap(histogramSumMap, maxKey);
      newKeyOldKeyMap = shiftMap(histogramMaxMap, maxKey);
      std::cout << "alignmax max at bin:" << maxKey << std::endl;
    }

     // file for root cell data
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);

    out << "bin aligned,bin,angle start,angle end,signal max,Signal Sum (Signal*um3)" << endl;

    forall(auto p, histogramSumMap){
      out << p.first << "," << newKeyOldKeyMap[p.first] << "," << binSize * newKeyOldKeyMap[p.first] << "," << binSize * (newKeyOldKeyMap[p.first]+1) << "," << histogramMaxMap[p.first] << "," << p.second << endl;
    }
    out << endl;

    std::cout << "Total Signal: " << totalSignal << std::endl;
    if(writeStack) store2->changed();
    return true;
  }
  REGISTER_PROCESS(StackCircularHistogram);


 bool StackHistogramBezier::initialize(QWidget* parent)
  {

    QString fileName = parm("Filename");
    if(parent)
      fileName = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(fileName.isEmpty())
      return false;
    if(!fileName.endsWith(".csv", Qt::CaseInsensitive))
      fileName += ".csv";
    setParm("Filename", fileName);;
    return true;
  }

   bool StackHistogramBezier::run(Stack* s1, Store* store1, Store* store2, Mesh* m, double binSize, double ignoreLow, double minDis, double maxDis,
    int bezPoints, QString filename, bool writeStack)
  {

    if(ignoreLow >= 100.0) return setErrorMessage("No voxels included. Lower the Voxel Value Threshold.");

    const HVecUS& dataStore = store1->data();
    HVecUS& dataStore2 = store2->data();
    Point3i imgSize(s1->size());

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    std::map<int, double> histogram, histcount;


    Point3f voxelSizes = s1->step();
    double voxelVolume = voxelSizes.x() * voxelSizes.y() * voxelSizes.z();

    // calculate length of the bezier line (for converting longitudinal coord into um)
    CuttingSurface* cutSurf = cuttingSurface();
    double totalBezLength;
    int dataPointsBezier = bezPoints;
    std::vector<Point3d> bezierVec, bDiff;

    Matrix4d rotMatrixS1, rotMatrixCS;
    s1->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());

    // fix the rotations
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    Bezier b = cutSurf->bezier();

    b.discretizeLineEqualWeight(dataPointsBezier, 0.0, mGLTot, bezierVec, bDiff, totalBezLength);

    int binNr = totalBezLength / binSize;
    double bezPointDis = totalBezLength / (double)dataPointsBezier;


    for(int i = 0; i <= binNr+1; i++){
      histogram[i] = 0;
      histcount[i] = 0;
    }

    double valueThreshold = ignoreLow*655.35; // * 16bit / 100

    // go through every voxel, check for thresholds, add it to the map
    #pragma omp parallel for schedule(guided)
    for(int x=0; x<xDim;x++){
      for(int y=0; y<yDim;y++){
        for(int z=0; z<zDim;z++){
          Point3i pPos (x,y,z);
          if(writeStack){
            dataStore2[s1->offset(pPos)]=0;
          }
          if(dataStore[s1->offset(pPos)] < valueThreshold) continue;

          Point3d pOutWorld(s1->imageToWorld(pPos));
          int idx;
          double idxWeight;
          Point3d closest = calcNearestPointOnBezierLine(pOutWorld, bezierVec, idx, idxWeight);
          //if(idxWeight < 0.5) idx--;

          if(maxDis > 0 or minDis>0){
            double dis = norm(closest - pOutWorld);
            if(dis > maxDis and maxDis > 0) continue;
            if(dis < minDis and minDis > 0) continue;
          }

          double value = (double)dataStore[s1->offset(pPos)]*voxelVolume;
          double bezDis = ((idx-1) + idxWeight) * bezPointDis;
          int bin = bezDis / binSize +1;
          if(idx == 1 and idxWeight < 0.01) bin = 0;
          if(idx == dataPointsBezier-1 and idxWeight > 0.99) bin = binNr+1;

          #pragma omp critical
          {
            histogram[bin] += value;
            histcount[bin]++;
            if(writeStack){
              if(bin == binNr+1) bin = 0;
              dataStore2[s1->offset(pPos)]=bin;
            }
          }
        }
      }
    }

    // write output to file
    int totalVoxels = 0;
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out(&file);

    out << "Bin,Sum Signal (Signal*um3),Nr Voxels With Signal above ";
    out << ignoreLow;
    out << "%";
    out << endl;

    forall(auto p, histogram){
      if(histcount[p.first] > 0){
        totalVoxels+=histcount[p.first];
      }
      if(p.first == 0) out << "-" << "," << p.second << ",";
      else if(p.first == binNr+1) out << "+" << "," << p.second << ",";
      else out << p.first << "," << p.second << ",";
      out << histcount[p.first];
      out << endl;
    }

    std::cout << "Considered " << totalVoxels << " Voxels (" << (double)totalVoxels/(double)(xDim*yDim*zDim)*100 << "% of all Voxels)" << std::endl;

    if(writeStack) store2->changed();

    return true;
  }
  REGISTER_PROCESS(StackHistogramBezier);



   bool StackBorderFromSegmented::run(Stack* s1, Store* store1, int width, bool ignX, bool ignY, bool ignZ)
  {

    Store* store2 = s1->work();

    const HVecUS& data = store1->data();
    HVecUS& data2 = store2->data();
    Point3i imgSize(s1->size());

    int xDim = imgSize.x();
    int yDim = imgSize.y();
    int zDim = imgSize.z();

    Point3f voxelSizes = s1->step();

    double voxelVolume = voxelSizes.x() * voxelSizes.y() * voxelSizes.z();

    Point3i lower(0,0,0);
    Point3i upper(xDim-1, yDim-1, zDim-1);

    int fillValue = 30000;

    // clear the work stack
    for(int z=0; z<zDim; z++){
      for(int x=0; x<xDim; x++){
        for(int y=0; y<yDim; y++){
          Point3i p(x,y,z);
          data2[s1->offset(p)] = 0;
        }
      }
    }

    //#pragma omp parallel for schedule(guided)
    // check for borders in y
    if(!ignY){
      for(int z=0; z<zDim; z++){
        for(int x=0; x<xDim; x++){
          Point3i last(x,0,z);
          for(int y=0; y<yDim; y++){
            Point3i current(x,y,z);

            if(data[s1->offset(last)] != data[s1->offset(current)]){
              // border found
              data2[s1->offset(last)] = fillValue;
              data2[s1->offset(current)] = fillValue;
            }
            last = current;
          }
        }
      }
    }


    // check for borders in x
    if(!ignX){
      for(int z=0; z<zDim; z++){
        for(int y=0; y<yDim; y++){
          Point3i last(0,y,z);
          for(int x=0; x<xDim; x++){
            Point3i current(x,y,z);

            if(data[s1->offset(last)] != data[s1->offset(current)]){
              // border found
              data2[s1->offset(last)] = fillValue;
              data2[s1->offset(current)] = fillValue;
            }
            last = current;
          }
        }
      }
    }

    // check for borders in z
    if(!ignZ){
      for(int x=0; x<xDim; x++){
        for(int y=0; y<yDim; y++){
          Point3i last(x,y,0);
          for(int z=0; z<zDim; z++){
            Point3i current(x,y,z);

            if(data[s1->offset(last)] != data[s1->offset(current)]){
              // border found
              data2[s1->offset(last)] = fillValue;
              data2[s1->offset(current)] = fillValue;
            }
            last = current;
          }
        }
      }
    }

    store2->copyMetaData(store1);
    store2->changed();

    return true;
  }
  REGISTER_PROCESS(StackBorderFromSegmented);

}
