//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2016 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 <MeshProcessSignal.hpp>
#include <Progress.hpp>
#include <GraphUtils.hpp> // shiftMap

#include <QFileDialog>
#include "ui_LoadHeatMap.h"
#include <omp.h>

namespace mgx 
{
  bool ViewMeshProcess::run()
  {
    QStringList warning;
  
    Mesh* m = currentMesh();
  
    if(not m)
      throw QString("%1:run() No active mesh").arg(name());
  
    if(not parm("Show Surface").isEmpty())
      m->setShowSurface(stringToBool(parm("Show Surface")));
    if(not parm("Use Parents").isEmpty()) {
      m->setUseParents(stringToBool(parm("Use Parents")));
      m->updateTriangles();
    }
 
    if(not parm("Surface Type").isEmpty()) {
      if(parm("Surface Type") == "Vertex")
        m->setShowVertex();
      else if(parm("Surface Type") == "Tris")
        m->setShowTriangle();
      else if(parm("Surface Type") == "Cells")
        m->setShowLabel();
      else
        warning << QString("Invalid surface type '%1'").arg(parm("Surface Type"));
    }
    if(not parm("Vertex Signal").isEmpty()) {
      if(parm("Vertex Signal") == "Signal")
        m->setShowVertex("Signal");
      else if(parm("Vertex Signal") == "Stack Texture")
        m->setShowVertex("Stack Texture");
      else if(parm("Vertex Signal") == "Image Texture")
        m->setShowVertex("Image Texture");
      else
        warning << QString("Invalid signal '%1'").arg(parm("Vertex Signal"));
    }
    if(not parm("Cells Signal").isEmpty()) {
      if(parm("Cells Signal") == "Label")
        m->setShowLabel("Label");
      else if(parm("Cells Signal") == "Label Heat")
        m->setShowLabel("Label Heat");
      else if(parm("Cells Signal") == "Wall Heat")
        m->setShowLabel("Wall Heat");
      else if(parm("Cells Signal") == "Cell Color")
        m->setShowLabel("Cell Color");
      else
        warning << QString("Invalid signal '%1'").arg(parm("Cells Signal"));
    }
    if(not parm("Blend").isEmpty())
      m->setBlending(stringToBool(parm("Blend")));
    if(not parm("Cull").isEmpty())
      m->setCulling(stringToBool(parm("Cull")));
    if(not parm("Show Mesh").isEmpty()) {
      m->setShowMesh(stringToBool(parm("Show Mesh")));
    }
    if(not parm("Mesh View").isEmpty())
      m->setMeshView(parm("Mesh View"));
    if(not parm("Show Lines").isEmpty())
      m->setShowMeshLines(stringToBool(parm("Show Lines")));
    if(not parm("Show Points").isEmpty())
      m->setShowMeshPoints(stringToBool(parm("Show Points")));
    if(not parm("Show Map").isEmpty())
      m->setShowMeshCellMap(stringToBool(parm("Show Map")));
    if(not parm("Scale").isEmpty())
      m->setScaled(stringToBool(parm("Scale")));
    if(not parm("Transform").isEmpty())
      m->setTransformed(stringToBool(parm("Transform")));
    if(not parm("BBox").isEmpty())
      m->setShowBBox(stringToBool(parm("BBox")));
    if(parm("Brightness").toFloat() >= 0)
      m->setBrightness(parm("Brightness").toFloat());
    if(parm("Opacity").toFloat() >= 0)
      m->setOpacity(parm("Opacity").toFloat());
    if(not warning.isEmpty())
      setWarningMessage(warning.join("\n"));

    return true;
  }
  REGISTER_PROCESS(ViewMeshProcess);
  
  Point2f calcSignalBounds(const vvGraph &S)
  {
    if(S.size() == 0)
      return  Point2f(0,0);
  
    float minSig = std::numeric_limits<float>::max();
    float maxSig = std::numeric_limits<float>::min();
    forall(const vertex &v, S) {
      if(minSig > v->signal)
        minSig = v->signal;
      if(maxSig < v->signal)
        maxSig = v->signal;
    }
    return Point2f(minSig, maxSig);
  }
  
  void ProjectSignal::projSignal(const Stack* stack, const HVecUS& data, vertex v, float mindist, float maxdist)
  {
    Point3d vpos = v->pos;
    Point3d vnrml = v->nrml;
    if(stack->showScale() and norm(stack->scale()) != 1.0) {
      vpos = divide(vpos, Point3d(stack->scale()));
      vnrml = divide(vnrml, Point3d(stack->scale()));
    }
    Point3d pmin = (vpos - vnrml * mindist);
    Point3d pmax = (vpos - vnrml * maxdist);
    Point3i pimin(stack->worldToImagei(Point3f(pmin)));
    Point3i pimax(stack->worldToImagei(Point3f(pmax)));
    int steps = int((pimax - pimin).norm() * 2 + 3);
  
    float sum = 0;
    for(int i = 0; i <= steps; i++) {
      Point3i oi(stack->worldToImagei(Point3f(pmax + v->nrml * (maxdist - mindist) * float(i) / float(steps))));
      if(stack->boundsOK(oi.x(), oi.y(), oi.z()))
        sum += data[stack->offset(oi.x(), oi.y(), oi.z())];
    }
    v->signal = float(sum) / float(steps);
  }
  
  bool ProjectSignal::run(const Store* store, Mesh* mesh, bool useAbsSignal, float mindist, 
    float maxdist, float absSignalMin, float absSignalMax)
  {
    const VtxVec& vs = mesh->activeVertices();
    setNormals(mesh->graph());
    const HVecUS& data = store->data();
    const Stack* stack = store->stack();
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i)
      projSignal(stack, data, vs[i], mindist, maxdist);
  
  
    if(useAbsSignal)
      mesh->signalBounds() = Point2f(absSignalMin, absSignalMax);
    else
      mesh->signalBounds() = calcSignalBounds(mesh->graph());
  
    SETSTATUS("Signal projection, min:" << mesh->signalBounds().x() << " max:" << mesh->signalBounds().y());
    mesh->signalUnit() = "";
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(ProjectSignal);
  

  void ProjectLabel::projLabel(const Stack* stack, const HVecUS& data, vertex v, float mindist, float maxdist, int zeroLabel, 
    QString method, bool ignoreZero)
  {
    Point3d vpos = v->pos;
    Point3d vnrml = v->nrml;
    if(stack->showScale() and norm(stack->scale()) != 1.0) {
      vpos = divide(vpos, Point3d(stack->scale()));
      vnrml = divide(vnrml, Point3d(stack->scale()));
    }
    Point3d pmin = (vpos - vnrml * mindist);
    Point3d pmax = (vpos - vnrml * maxdist);
    Point3i pimin(stack->worldToImagei(Point3f(pmin)));
    Point3i pimax(stack->worldToImagei(Point3f(pmax)));
    int steps = int((pimax - pimin).norm() * 2 + 3);
  
    std::map<int, int> labelCounterMap;
    for(int i = 0; i <= steps; i++) {
      Point3i oi(stack->worldToImagei(Point3f(pmax + v->nrml * (maxdist - mindist) * float(i) / float(steps))));
      if(stack->boundsOK(oi.x(), oi.y(), oi.z()))
        labelCounterMap[data[stack->offset(oi.x(), oi.y(), oi.z())]]++;
    }

    int maxC = -1;
    int winningLabel = -1;
    int totalCount = 0;
    forall(auto p, labelCounterMap){
      totalCount+= p.second;
      if(maxC < p.second){
        if(ignoreZero and p.first == 0) continue;
        maxC = p.second;
        winningLabel = p.first;
      }
    }

    v->label = winningLabel;

    if(method == "Absolute Majority"){
      if(maxC < totalCount/2) v->label = 0;
    }

    if(winningLabel == 0)
      v->label = zeroLabel;
    
  }
  
  bool ProjectLabel::run(const Store* store, Mesh* mesh, float mindist, 
    float maxdist, QString method, bool ignoreZero)
  {
    const VtxVec& vs = mesh->activeVertices();
    const vvGraph& S = mesh->graph();
    setNormals(S);
    const HVecUS& data = store->data();
    const Stack* stack = store->stack();

    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i)
      projLabel(stack, data, vs[i], mindist, maxdist, -1, method, ignoreZero);
  
    int maxLabel = 0;
    forall(const vertex& v, S){
      if(v->label > maxLabel) maxLabel = v->label;
    }

    forall(const vertex& v, S){
      if(v->label == -1 and !ignoreZero) v->label = maxLabel+1;
    }

    //SETSTATUS("Signal projection, min:" << mesh->signalBounds().x() << " max:" << mesh->signalBounds().y());
    //mesh->signalUnit() = "";
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(ProjectLabel);


   bool ProjectLabelsOnMesh::run(Stack* s1, Store* store1, Mesh* m, bool splitBorder, int searchNeighborhood)
  {
    // find nearest voxel for each mesh vertex and label it accordingly

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

    vvGraph& S = m->graph();

    //if(zDim > 1) return setErrorMessage("z dimension larger 1");

    forall(const vertex& v, S){
      Point3f p = s1->worldToImagef(Point3f(v->pos));
      p.x() = std::round(p.x());
      p.y() = std::round(p.y());
      p.z() = std::round(p.z());
      int label = data[s1->offset(Point3i(p))];
      if(label == 0){
        double minDis = 1E20;
        for(int sX = -searchNeighborhood; sX <= searchNeighborhood; sX++){
          for(int sY = -searchNeighborhood; sY <= searchNeighborhood; sY++){
            for(int sZ = -searchNeighborhood; sZ <= searchNeighborhood; sZ++){
              if(data[s1->offset(Point3i(p)+Point3i(sX,sY,sZ))] > 0 and norm(Point3i(sX,sY,sZ)) < minDis){
                label = data[s1->offset(Point3i(p)+Point3i(sX,sY,sZ))];
                minDis = norm(Point3i(sX,sY,sZ));
              }
            }
          }
        }
      }
      v->label = label;
    }
    if(splitBorder){
      Subdivide *sDiv = NULL;

      typedef std::pair<vertex,vertex> VtxVtx;
      std::set<VtxVtx> edgesToSplit;

      forall(const vertex &v, S){
        forall(const vertex &n, S.neighbors(v)){
          if(n->label != v->label and (n->label > 0 and v->label > 0)){
            VtxVtx vv;
            if(v < n) vv = std::make_pair(v,n);
            else vv = std::make_pair(n,v);
            edgesToSplit.insert(vv);
          }
        }
      }

      forall(const VtxVtx& vv, edgesToSplit){
        splitEdge(S, vv.first, vv.second, sDiv);
        vertex newV;
        forall(const vertex& v, S.neighbors(vv.first)){
          if(norm((vv.first->pos + vv.second->pos)/2.0 - v->pos) < 0.01){ // find the added vertex and set it to -1
            newV = v;
          }
        }
        newV->label = -1;
      }
    }

    m->updateAll();

    return true;
  }
  REGISTER_PROCESS(ProjectLabelsOnMesh);

  bool SmoothMeshSignal::run(Mesh* mesh, uint passes)
  {
    VtxFloatAttr tSignal;
    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
    if(vs.empty())
      return true;
    for(uint i = 0; i < passes; i++) {
      forall(const vertex& v, vs) {
        float tarea = 0;
        v->*tSignal = 0;
        forall(const vertex& n, S.neighbors(v)) {
          vertex m = S.nextTo(v, n);
          if(!S.edge(m, n))
            continue;
          float a = triangleArea(v->pos, n->pos, m->pos);
          tarea += a;
          v->*tSignal += a * (v->signal + n->signal + m->signal) / 3.0;
        }
        if(tarea > 0)
          v->*tSignal /= tarea;
        else
          v->*tSignal = v->signal;
      }
      forall(const vertex& v, vs)
        v->signal = v->*tSignal;
    }
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(SmoothMeshSignal);

  bool MeshBrightness::run(Mesh* mesh, float amount)
  {
    const std::vector<vertex>& vs = mesh->activeVertices();
    if(vs.empty())
      return true;

    forall(const vertex& v, vs)
      v->signal *= amount;

    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshBrightness);
 
  bool ClearMeshSignal::run(Mesh* mesh, uint value)
  {
    const std::vector<vertex>& vs = mesh->activeVertices();
    if(vs.empty())
      return true;
    forall(const vertex& v, vs) {
      v->signal = value;
    }
    SETSTATUS("Mesh signal set to:" << value);
    mesh->signalUnit() = "";
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(ClearMeshSignal);
  
  ProjectCurvatureOp::ProjectCurvatureOp(const vvGraph& _S, const VtxVec& _vs, std::string _type, float _radius,
                                         FloatVec& _values, std::vector<Curvature>& _curvs)
    : S(_S), vs(_vs), type(_type), radius(_radius), values(_values) , curvs(_curvs) {}
  
  void ProjectCurvatureOp::operator()(int i)
  {
    vertex v = vs[i];
    VtxSet nbs;
    Mesh::getNhbd(S, v, radius, nbs);
  
    std::vector<Point3f> plist, nlist;
    plist.push_back(Point3f(v->pos));
    nlist.push_back(Point3f(v->nrml));
    forall(const vertex& n, nbs) {
      if(n == v)
        continue;
      plist.push_back(Point3f(n->pos));
      nlist.push_back(Point3f(n->nrml));
    }
    if(!curvs.empty()) {
      curvs[i].update(plist, nlist);
      v->signal = curvs[i].get_curvature();
    } else {
      Curvature curv(type);
      curv.update(plist, nlist);
      v->signal = curv.get_curvature();
    }
    if(!values.empty())
      values[i] = v->signal;
  }
  
  bool ProjectCurvature::run(Mesh* mesh, QString output, QString type, float radius, bool auto_scale,
                                    float mincurv, float maxcurv, int percentile)
  {
    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
    if(vs.empty())
      return true;
    type = type.toLower();
    if((type == "sumsquare" or type == "anisotropy") and mincurv < 0)
      mincurv = 0;
    if(!auto_scale and mincurv >= maxcurv)
      return true;
    if(auto_scale) {
      mincurv = FLT_MAX;
      maxcurv = -FLT_MAX;
    }
  
    std::vector<float> values(auto_scale ? vs.size() : 0);
    std::vector<Curvature> curvatures;
    if(!output.isEmpty()) {
      Curvature curv(type.toStdString());
      for(uint i = 0; i < vs.size(); i++)
        curvatures.push_back(curv);
    }
  
    progressStart(QString("Computing curvature for Mesh %1").arg(mesh->userId()),
                      vs.size() / omp_get_max_threads());
    ProjectCurvatureOp op(S, vs, type.toStdString(), radius, values, curvatures);
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    if(auto_scale) {
      std::sort(values.begin(), values.end());
      int rankmax = values.size() - 1;
      int rankmin = 0;
      if(percentile < 100) {
        rankmax = (values.size() * percentile) / 100;
        rankmin = values.size() - rankmax;
      }
      mincurv = values[rankmin];
      maxcurv = values[rankmax];
      if(maxcurv * mincurv < 0) {
        maxcurv = std::max(fabs(maxcurv), fabs(mincurv));
        mincurv = -maxcurv;
      } else if(maxcurv > 0)
        mincurv = 0;
      else if(mincurv < 0)
        maxcurv = 0;
      else {
        mincurv = 0;
        maxcurv = 1;
      }
    }
  
    if(!output.isEmpty()) {
      QFile f(output);
      if(f.open(QIODevice::WriteOnly)) {
        QTextStream ts(&f);
        ts << QString("X,Y,Z,Label,X1,Y1,Z1,X2,Y2,Z2,k1 (1/%1),k2 (1/%1)\n").arg(UM);
        size_t i = 0;
        forall(const vertex &v, S) {
          Point3f pos(v->pos);
          Curvature& curv = curvatures[i];
          Point3f e1, e2;
          float c1, c2;
          curv.get_eigenVectors(e1, e2);
          curv.get_eigenValues(c1, c2);
          QStringList fields;
          fields << QString::number(pos.x(), 'g', 10) << QString::number(pos.y(), 'g', 10)
                 << QString::number(pos.z(), 'g', 10) << QString::number(v->label)
                 << QString::number(e1.x(), 'g', 10) << QString::number(e1.y(), 'g', 10)
                 << QString::number(e1.z(), 'g', 10) << QString::number(e2.x(), 'g', 10)
                 << QString::number(e2.y(), 'g', 10) << QString::number(e2.z(), 'g', 10)
                 << QString::number(c1, 'g', 10) << QString::number(c2, 'g', 10);
          ts << fields.join(",") << "\n";
          ++i;
        }
      } else
        setWarningMessage(QString("Warning, could not open file '%1' for writing. Not output written").arg(output));
    }
  
    if(type == "gaussian" or type == "sumsquare")
      mesh->signalUnit() = UM_2;
    else if(type == "anisotropy")
      mesh->signalUnit() = "";
    else
      mesh->signalUnit() = UM_1;
  
    mesh->signalBounds() = Point2f(mincurv, maxcurv);
  
    mesh->setShowVertex("Signal");
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(ProjectCurvature);
 


  MeshSignalGradOp::MeshSignalGradOp(const vvGraph& _S, const VtxVec& _vs, VtxFloatAttr & _result)
    : S(_S), vs(_vs), result(_result) {}

  void MeshSignalGradOp::operator()(int i)
  {
    vertex v = vs[i];
    float num = 0,  wt = 0;
    Point3d gradient(0,0,0);

    // Gaussian blur over the cell
    forall(const vertex& r, S.neighbors(v)) {
	    const vertex &n = S.nextTo(v,r);
	    num++;  
	    wt++;//norm((r->pos-v->pos)^(n->pos-v->pos));
	    Point3d tgrad = triangleGradient(v->pos, r->pos, n->pos, double(v->signal), double(r->signal), double(n->signal)) /
                                                                   2.0 / norm((r->pos-v->pos)^(n->pos-v->pos));
	    gradient+= tgrad;
    }

	  //std::cout<<"In op: finished"<<std::endl;
    v->*result = norm(gradient)/wt;   // Blurring the signal
  }
  
  bool MeshSignalGrad::run(Mesh* mesh)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    std::cout<<"running sig grad"<<std::endl;
    progressStart(QString("Gaussian Blur on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxFloatAttr result;

    std::cout<<"Initializing op (sig grad)"<<std::endl;

    MeshSignalGradOp op(S, vs, result);

    std::cout<<"Running op (sig grad)"<<std::endl;

    //#pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();

    // Update the Signal
    forall(const vertex& v, vs)
      v->signal = v->*result;

    std::cout<<"Updating tri (sig grad)"<<std::endl;
 
    mesh->updateTriangles();

    std::cout<<"Completed sig grad"<<std::endl;


    return true;
  }
  REGISTER_PROCESS(MeshSignalGrad);


 
  MeshGaussianBlurOp::MeshGaussianBlurOp(
       const vvGraph& _S, const VtxVec& _vs, float _radius, VtxFloatAttr & _result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshGaussianBlurOp::operator()(int i)
  {
    vertex v = vs[i];
    float sigma = radius / 2;
    VtxSet VSet;
    float den = 0, num = 0, g = 0, ar = 0, wt = 0;
    Mesh::getNhbd(S, v, radius, VSet);   // Neighbourhood processing
  
    // Gaussian blur over the cell
    forall(const vertex& r, VSet) {
      if(norm(r->pos - v->pos) < radius) {
        den = sqrt(6.28 * (pow(sigma, 2)));
        num = exp(-(pow((norm(r->pos - v->pos)), 2) / (2 * pow(sigma, 2))));
        g = num / den;
        ar = ar + (g * (r->signal));
        wt = wt + g;
      }
    }
    v->*result = ar / wt;   // Blurring the signal
  }
  
  bool MeshGaussianBlur::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
  
    progressStart(QString("Gaussian Blur on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxFloatAttr result;
    MeshGaussianBlurOp op(S, vs, radius, result);
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the Signal
    forall(const vertex& v, vs)
      v->signal = v->*result;
  
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshGaussianBlur);
  
  bool MeshDiffGaussians::run(Mesh* mesh, float radius1, float radius2)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
  
    // Find first Gaussian, allocate space for result, need to initialize or not thread safe
    VtxFloatAttr firstGauss;
    progressStart(QString("Calculating first Gaussian on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    if(radius1 > 0) {
      MeshGaussianBlurOp op1(S, vs, radius1, firstGauss);
      #pragma omp parallel for
      for(size_t i = 0; i < vs.size(); ++i) {
        if(!progressAdvance(i))
          continue;
        op1(i);
      }
      if(!progressAdvance(vs.size()))
        userCancel();
    }
  
    progressStart(QString("Calculating second Gaussian on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Find second Gaussian, allocate space for result, need to initialize or not thread safe
    VtxFloatAttr secondGauss;
    MeshGaussianBlurOp op2(S, vs, radius2, secondGauss);
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op2(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the Signal
    forall(const vertex& v, vs)
      v->signal = v->*firstGauss - v->*secondGauss;
  
    mesh->signalBounds() = calcSignalBounds(mesh->graph());
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshDiffGaussians);
  
  MeshLocalMinimaOp::MeshLocalMinimaOp(const vvGraph& _S, const VtxVec& _vs, float _radius, int _label, bool _overwrite)
    : S(_S) , vs(_vs) , radius(_radius) , label(_label), overwrite(_overwrite) {}
  
  void MeshLocalMinimaOp::operator()(int i)
  {
    vertex v = vs[i];
    VtxSet VSet;
    Mesh::getNhbd(S, v, radius, VSet);   // Neighbourhood processing
  
    // Finding the lowest v->signal in the set of vertices respective to a cell
    float lowest = HUGE_VAL;
    bool otherLabel = false;;
    forall(const vertex& r, VSet){
      if(r->label > 0)
        otherLabel = true;
      if(r->signal < lowest)
        lowest = r->signal;
    }

    // Label the vertex at the local minima
    if(v->signal == lowest){
      if(!otherLabel or !overwrite)
        v->label = label;
    }
  }
  
  bool MeshLocalMinima::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
  
    progressStart(QString("Finding Local Minima for Mesh %1").arg(mesh->userId()),
                      vs.size() / omp_get_max_threads());
    int startLabel = mesh->nextLabel();
    bool overwrite = stringToBool(parm("Consider existing Labels"));
    MeshLocalMinimaOp op(S, vs, radius, startLabel, overwrite);
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
    for(size_t i = 0; i < vs.size(); ++i)
      if(vs[i]->label == startLabel) {
        int label = mesh->nextLabel();
        vs[i]->label = label;
        forall(const vertex& n, S.neighbors(vs[i]))
          n->label = label;
      }
  
    mesh->setShowLabel();
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshLocalMinima);
  
  MeshNormalizeOp::MeshNormalizeOp(const vvGraph& _S, const VtxVec& _vs, 
        float _radius, VtxFloatAttr & _result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshNormalizeOp::operator()(int i)
  {
    VtxSet VSet;
    vertex v = vs[i];
    float min = FLT_MAX;
    float max = -FLT_MAX;
    Mesh::getNhbd(S, v, radius, VSet);
  
    // Normalizing the signal of mesh
    forall(const vertex& n, VSet) {
      if(min > n->signal)
        min = n->signal;
      if(max < n->signal)
        max = n->signal;
    }
    v->*result = (v->signal - min) / (max - min);
  }
  
  bool MeshNormalize::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    progressStart(QString("Normalizing Signal on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    VtxFloatAttr result;
    MeshNormalizeOp op(S, vs, radius, result);
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the Signal
    forall(const vertex& v, vs)
      v->signal = v->*result * 65535.0;
  
    mesh->signalBounds() = calcSignalBounds(S);
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshNormalize);
  
  MeshDilationOp::MeshDilationOp(const vvGraph& _S, const VtxVec& _vs, float _radius, VtxFloatAttr &_result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshDilationOp::operator()(int i)
  {
    VtxSet VSet;
    vertex v = vs[i];
    float max = -FLT_MAX;
    Mesh::getNhbd(S, v, radius, VSet);
  
    // Dilation
    forall(const vertex& n, VSet)
      if(max < n->signal)
        max = n->signal;
  
    v->*result = max;
  };
  
  bool MeshDilation::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    VtxSet VSet;
    progressStart(QString("Dilating Signal on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxFloatAttr result;
    MeshDilationOp op(S, vs, radius, result);
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the signal
    forall(const vertex& v, vs)
      v->signal = v->*result;
  
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshDilation);
  

  MeshLabelDilationOp::MeshLabelDilationOp(const vvGraph& _S, const VtxVec& _vs, float _radius, VtxIntAttr &_result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshLabelDilationOp::operator()(int i)
  {
    VtxSet VSet;
    vertex v = vs[i];
    std::set<int> labels;
    Mesh::getNhbd(S, v, radius, VSet);
  
    v->*result = v->label;

    if(v->label !=0)
      return;

    // Dilation
    forall(const vertex& n, VSet){
      if(n->label > 0 and v->label != n->label)
        labels.insert(n->label);
    }
  
    if(labels.size() == 1)
      v->*result = *(labels.begin());

  };
  
  bool MeshLabelDilation::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    VtxSet VSet;
    progressStart(QString("Dilating Signal on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxIntAttr result;
    MeshLabelDilationOp op(S, vs, radius, result);
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the signal
    forall(const vertex& v, vs)
      v->label = v->*result;
  
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshLabelDilation);

  MeshErosionOp::MeshErosionOp(const vvGraph& _S, const VtxVec& _vs, float _radius, VtxFloatAttr &_result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshErosionOp::operator()(int i)
  {
    VtxSet VSet;
    vertex v = vs[i];
    float min = FLT_MAX;
    Mesh::getNhbd(S, v, radius, VSet);
  
    // Erosion
    forall(const vertex& n, VSet)
      if(min > n->signal)
        min = n->signal;
  
    v->*result = min;
  };
  
  bool MeshErosion::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    VtxSet VSet;
    progressStart(QString("Eroding Signal on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxFloatAttr result;
    MeshErosionOp op(S, vs, radius, result);
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the signal
    forall(const vertex& v, vs)
      v->signal = v->*result;
  
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshErosion);
  
  MeshLabelErosionOp::MeshLabelErosionOp(const vvGraph& _S, const VtxVec& _vs, float _radius, VtxIntAttr &_result)
    : S(_S), vs(_vs), radius(_radius), result(_result) {}
  
  void MeshLabelErosionOp::operator()(int i)
  {
    VtxSet VSet;
    vertex v = vs[i];
    int min = v->label;
    Mesh::getNhbd(S, v, radius, VSet);
  
    // Erosion
    forall(const vertex& n, VSet)
      if(n->label != v->label)
        min = 0;
  
    v->*result = min;
  };

  bool MeshLabelErosion::run(Mesh* mesh, float radius)
  {
    vvGraph& S = mesh->graph();
    const VtxVec& vs = mesh->activeVertices();
    VtxSet VSet;
    progressStart(QString("Eroding Labels on Mesh %1").arg(mesh->userId()), vs.size() / omp_get_max_threads());
    // Allocate space for result, need to initialize or not thread safe
    VtxIntAttr result;
    MeshLabelErosionOp op(S, vs, radius, result);
  
    #pragma omp parallel for
    for(size_t i = 0; i < vs.size(); ++i) {
      if(!progressAdvance(i))
        continue;
      op(i);
    }
    if(!progressAdvance(vs.size()))
      userCancel();
  
    // Update the label
    forall(const vertex& v, vs)
      v->label = v->*result;
  
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshLabelErosion);

  bool MeshClosing::run(Mesh* mesh, float radius)
  {
    MeshDilation md(*this);
    md.run(mesh, radius);

    MeshErosion me(*this);
    me.run(mesh, radius);

    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshClosing);
  
  bool MeshOpening::run(Mesh* mesh, float radius)
  {
    MeshErosion me(*this);
    me.run(mesh, radius);

    MeshDilation md(*this);
    md.run(mesh, radius);

    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(MeshOpening);
  
  bool RescaleSignal::run(Mesh* mesh, bool use_zero, float percentile, float minimum, float maximum)
  {
    if(percentile > 0) {
      vvGraph& S = mesh->graph();
      std::vector<float> values(S.size());
      int i = 0;
      forall(const vertex& v, S)
        values[i++] = v->signal;
      std::sort(values.begin(), values.end());
  
      int rankmax = values.size() - 1;
      int rankmin = 0;
      if(percentile < 100) {
        rankmax = (values.size() * percentile) / 100;
        rankmin = values.size() - rankmax;
      }
      minimum = values[rankmin];
      maximum = values[rankmax];
      if(use_zero) {
        if(minimum * maximum < 0) {
          maximum = std::max(fabs(maximum), fabs(minimum));
          minimum = -maximum;
        } else if(maximum > 0)
          minimum = 0;
        else if(minimum < 0)
          maximum = 0;
        else {
          minimum = 0;
          maximum = 1;
        }
      }
    }
  
    mesh->signalBounds() = Point2f(minimum, maximum);
    return true;
  }
  REGISTER_PROCESS(RescaleSignal);
  
  bool SetSignalProcess::run(Mesh* m, float value, bool rescale, float percentile, bool use_zero)
  {
    if(rescale and (percentile <= 0 or percentile > 100))
      return setErrorMessage("Percentile must be greater than 0 and less or equal to 100.");
    const std::vector<vertex>& vs = m->activeVertices();
    forall(const vertex& v, vs)
      v->signal = value;
    if(rescale) {
      RescaleSignal rescale(*this);
      rescale.run(m, use_zero, percentile, 0, 0);
    }
    m->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(SetSignalProcess);


 bool MeshCircularHistogram::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 MeshCircularHistogram::run(Stack* s1, Mesh* m1, QString centralAxis, int binNumber, double thresholdValue,
    double minDis, double maxDis, QString align, bool weightVol, QString filename, bool writeLabels)
  {

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

    Point3d refDir2 = cAxis ^ refDir;

    double totalSignal = 0;

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

    vvGraph& S = m1->graph();
    Point3d zeroP(0,0,0);

    // //#pragma omp parallel for schedule(guided)
    forall(vertex v, S){
      forall(vertex n, S.neighbors(v)){
        vertex m = S.nextTo(v,n);
        if(!S.uniqueTri(v,n,m)) continue;

        // Total area, signal, and label center
        float area = triangleArea(v->pos, n->pos, m->pos);
        Point3d triPoint = (v->pos + n->pos + m->pos) / 3.0;
        triPoint = multMatrix4Point3(mGLRotMat,triPoint); // correct rotations

        Point3d pointPlane = projectPointOnPlane(triPoint, zeroP, 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 - zeroP;
        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 = area * (v->signal + n->signal + m->signal) / 3.0;
        if(!weightVol) weightedValue = (v->signal + n->signal + m->signal) / 3.0;
        histogramSumMap[bin] += weightedValue;
        histogramMaxMap[bin] = std::max(histogramMaxMap[bin],weightedValue);
        if(writeLabels){
          v->label = bin;
          n->label = bin;
          m->label = 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;
    }

     // save to file
    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*um2)" << 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;

    m1->updateTriangles();

    return true;
  }
  REGISTER_PROCESS(MeshCircularHistogram);


}
