//
// 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 <SystemProcessLoad.hpp>
#include "MeshProcessCellAxis.hpp"
#include "MeshProcessLineage.hpp"
#include "Dir.hpp"
#include "Triangulate.hpp"
#include "MeshProcessCurv.hpp"
#include "MeshProcessFibril.hpp"
#include "MeshProcessMeasures.hpp"
#include "MeshProcessPDG.hpp"
#include "MeshProcessSignalOrientation.hpp"
#include "PCAnalysis.hpp"
#include "Progress.hpp"
#include <QFileDialog>
#include "GraphUtils.hpp"
#include "GraphAlgorithms.hpp"

namespace mgx {

  bool SaveCellAxis::initialize(QWidget* parent)
  {
    QString filename = parm("Output File");
    if(filename.isEmpty() and parent)
      filename = QDir::currentPath();
    filename = QFileDialog::getSaveFileName(parent, "Choose spreadsheet file to save", filename, "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Output File", stripCurrentDir(filename));
    return true;
  }

  bool SaveCellAxis::run(Mesh* mesh, const QString& filename, const QString& type)
  {
    QFile outputFile;
    QTextStream ts;
    if(filename.isEmpty()) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }

    outputFile.setFileName(filename);
    if(!outputFile.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("Unable to write to output file"));
      return false;
    }

    ts.setDevice(&outputFile);
    if(mesh->cellAxisType() == "PDG"){
		  if(type == "Original Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,stretchRatioMax,stretchRatioMin,stretchRatioNormal") << endl;
			else if(type == "Custom Tensor")
        ts << QString("Label,customXMax,customYMax,customZMax,customXMin,customYMin,customZMin,stretchRatioMax,stretchRatioMin,shear") << endl;
      else
        ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,xNorm,yNorm,zNorm") << endl;
		}
    else if(mesh->cellAxisType() == "curvature")
      ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,1/radiusCurvMax,1/radiusCurvMin,1/radiusCurvNormal") << endl;
    else if(mesh->cellAxisType() == "fibril")
      ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,alignmentMax,alignmentMin,alignmentNormal") << endl;
		else if(mesh->cellAxisType() == "polarization"){// from the addOn "SignalOrientation" (for PIN signal)
      ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,signalMax,signalMin,signalNormal") << endl;
			}
    else if(mesh->cellAxisType() == "PCA"){
		  if(type == "Original Tensor")
        ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,pcaMax,pcaMin,pcaNormal") << endl;
      else if(type == "Custom Tensor")
        ts << QString("Label,customXMax,customYMax,customZMax,customXMin,customYMin,customZMin,pcaMax,pcaMin,shear") << endl;
      else
        ts << QString("Label,xMax,yMax,zMax,xMin,yMin,zMin,xNorm,yNorm,zNorm") << endl;
			}
		else {

		  //setErrorMessage(QString("Unknown cell axis type"));
      std::cout << "Unknwon cell axis type! Saving Default file" << std::endl;
      return false;
    }

    if(type == "Original Tensor"){

      forall(const IntSymTensorPair& p, mesh->cellAxis()) {
        int label = p.first;
        const SymmetricTensor& tensor = p.second;
        const Point3f& vMax = tensor.ev1();
        const Point3f& vMin = tensor.ev2();
        const Point3f& eig = tensor.evals();
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
           << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << endl;
      }

    } else if(type == "Custom Tensor"){

      const IntSymTensorAttr& cellAxisCustom = currentMesh()->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
      if(cellAxisCustom.empty()) return setErrorMessage("No Custom Directions found.");

      forall(const IntSymTensorPair& p, mesh->cellAxis()) {
        int label = p.first;
        Point3f customX = cellAxisCustom[label].ev1();
        Point3f customY = cellAxisCustom[label].ev2();
        const SymmetricTensor& tensor = p.second;

        SymmetricTensor tensorCustom = customDirectionTensor(tensor, customX, customY, false);

        const Point3f& vMax = tensorCustom.ev1();
        const Point3f& vMin = tensorCustom.ev2();
        const Point3f& eig = tensorCustom.evals();
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
           << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << endl;
      }

    } else { // save visualized tensor

      forall(const auto& p, mesh->cellAxisVis()) {
        int label = p.first;
        const Point3f& vMax = Point3f(p.second[0][0], p.second[0][1], p.second[0][2]);
        const Point3f& vMin = Point3f(p.second[1][0], p.second[1][1], p.second[1][2]);
        const Point3f& eig = Point3f(p.second[2][0], p.second[2][1], p.second[2][2]);
        ts << label << "," << vMax.x() << "," << vMax.y() << "," << vMax.z() << "," << vMin.x() << "," << vMin.y()
             << "," << vMin.z() << "," << eig.x() << "," << eig.y() << "," << eig.z() << endl;
      }
    }

    return true;
  }
  REGISTER_PROCESS(SaveCellAxis);

  bool LoadCellAxis::initialize(QWidget* parent)
  {
    QString filename = parm("Input File");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getOpenFileName(parent, "Choose spreadsheet to load", QDir::currentPath(),
                                              "CSV files (*.csv);;All files (*.*)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Input File", stripCurrentDir(filename));
    return true;
  }

  bool LoadCellAxis::run(Mesh* mesh, const QString& type, const QString& filename)
  {
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ts(&file);
    QString line = ts.readLine();
    QStringList fields = line.split(",");
    QHash<QString, int> fieldMap;
    int field_xMax, field_yMax, field_zMax;
    int field_xMin, field_yMin, field_zMin;
    int field_k1, field_k2, field_k3;
    int field_label;
    if(fields.size() < 10) {
      setErrorMessage(QString("File '%1' should be a CSV file with at least 10 columns").arg(filename));
      return false;
    }
    // First 7 columns of csv files are the same for all cell axis types:
		// Label, first direction (x,y,z), second direction (x,y,z)
    field_label = fields.indexOf("Label");
    if(field_label < 0)
      return setErrorMessage("Error, no field 'Label' in the CSV file");
    field_xMax = max(fields.indexOf("xMax"), fields.indexOf("x1"));
    if(field_xMax < 0)
      return setErrorMessage("Error, no field 'xMax' nor 'x1' in the CSV file");
    field_yMax = max(fields.indexOf("yMax"), fields.indexOf("y1"));
    if(field_yMax < 0)
      return setErrorMessage("Error, no field 'yMax' nor 'y1' in the CSV file");
    field_zMax = max(fields.indexOf("zMax"), fields.indexOf("z1"));
    if(field_zMax < 0)
      return setErrorMessage("Error, no field 'zMax' nor 'z1' in the CSV file");
    field_xMin = max(fields.indexOf("xMin"), fields.indexOf("x2"));
    if(field_xMin < 0)
      return setErrorMessage("Error, no field 'xMin' nor 'x2' in the CSV file");
    field_yMin = max(fields.indexOf("yMin"), fields.indexOf("y2"));
    if(field_yMin < 0)
      return setErrorMessage("Error, no field 'yMin' nor 'y2' in the CSV file");
    field_zMin = max(fields.indexOf("zMin"), fields.indexOf("z2"));
    if(field_zMin < 0)
      return setErrorMessage("Error, no field 'zMin' nor 'z2' in the CSV file");

    // Next columns of csv files are called differently depending on cell axis types.
		// norm of first vector, norm of second vector, norm of third vector (= 0 for 2D cell axis)
    if(type == QString("PDG")) {
      field_k1 = max(fields.indexOf("stretchRatioMax"), fields.indexOf("k1"));
      if(field_k1 < 0)
        return setErrorMessage("Error, no field 'stretchRatioMax' nor 'k1' in the CSV file");
      field_k2 = max(fields.indexOf("stretchRatioMin"), fields.indexOf("k2"));
      if(field_k2 < 0)
        return setErrorMessage("Error, no field 'stretchRatioMin' nor 'k2' in the CSV file");
      field_k3 = max(fields.indexOf("stretchRatioNormal"), fields.indexOf("k3"));
      if(field_k3 < 0)
        return setErrorMessage("Error, no field 'stretchRatioNormal' nor 'k3' in the CSV file");
    } else if(type == QString("curvature")) {
      field_k1 = max(fields.indexOf("1/radiusCurvMax"), fields.indexOf("k1"));
      if(field_k1 < 0)
        return setErrorMessage("Error, no field '1/radiusCurvMax' nor 'k1' in the CSV file");
      field_k2 = max(fields.indexOf("1/radiusCurvMin"), fields.indexOf("k2"));
      if(field_k2 < 0)
        return setErrorMessage("Error, no field '1/radiusCurvMin' nor 'k2' in the CSV file");
      field_k3 = max(fields.indexOf("1/radiusCurvNormal"), fields.indexOf("k3"));
      if(field_k3 < 0)
        return setErrorMessage("Error, no field '1/radiusCurvNormal' nor 'k3' in the CSV file");
    } else if(type == QString("fibril")) {
      field_k1 = max(fields.indexOf("alignmentMax"), fields.indexOf("k1"));
      if(field_k1 < 0)
        return setErrorMessage("Error, no field 'alignmentMax' nor 'k1' in the CSV file");
      field_k2 = max(fields.indexOf("alignmentMin"), fields.indexOf("k2"));
      if(field_k2 < 0)
        return setErrorMessage("Error, no field 'alignmentMin' nor 'k2' in the CSV file");
      field_k3 = max(fields.indexOf("alignmentNormal"), fields.indexOf("k3"));
      if(field_k3 < 0)
        return setErrorMessage("Error, no field 'alignmentNormal' nor 'k3' in the CSV file");
    } else if(type == QString("polarization")) {
      field_k1 = max(fields.indexOf("signalMax"), fields.indexOf("k1"));
      if(field_k1 < 0)
        return setErrorMessage("Error, no field 'signalMax' nor 'k1' in the CSV file");
      field_k2 = max(fields.indexOf("signalMin"), fields.indexOf("k2"));
      if(field_k2 < 0)
        return setErrorMessage("Error, no field 'signalMin' nor 'k2' in the CSV file");
      field_k3 = max(fields.indexOf("signalNormal"), fields.indexOf("k3"));
      if(field_k3 < 0)
        return setErrorMessage("Error, no field 'signalNormal' nor 'k3' in the CSV file");
    } else if(type == QString("PCA")) {
      field_k1 = max(fields.indexOf("pcaMax"), fields.indexOf("k1"));
      if(field_k1 < 0)
        return setErrorMessage("Error, no field 'pcaMax' nor 'k1' in the CSV file");
      field_k2 = max(fields.indexOf("pcaMin"), fields.indexOf("k2"));
      if(field_k2 < 0)
        return setErrorMessage("Error, no field 'pcaMin' nor 'k2' in the CSV file");
      field_k3 = max(fields.indexOf("pcaNormal"), fields.indexOf("k3"));
      if(field_k3 < 0)
        return setErrorMessage("Error, no field 'pcaNormal' nor 'k3' in the CSV file");
    }	else {
		  setErrorMessage(QString("Unknown cell axis type"));
      return false;
    }

    int nb_fields = fields.size();

    mesh->clearCellAxis();
    IntSymTensorAttr& cellAxis = mesh->cellAxis();

    int linenum = 1;
    while(ts.status() == QTextStream::Ok) {
      ++linenum;
      fields = ts.readLine().split(",");
      if(fields.empty() or fields.size() == 1)
        break;

      if(fields.size() < nb_fields)
        return setErrorMessage(
          QString("Error on line %1: not enough fields (Expected: %2, Read: %3)").arg(linenum).arg(nb_fields).arg(
            fields.size()));

      bool ok;
      int label = fields[field_label].toInt(&ok);
      if(!ok) {
        setErrorMessage(
          QString("Error line %1: invalid label '%2'").arg(linenum).arg(fields[field_label].trimmed()));
        return false;
      }

      Point3f ev1, ev2, evals;
      ev1.x() = fields[field_xMax].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid xMax '%2'").arg(linenum).arg(fields[field_xMax]));
      ev1.y() = fields[field_yMax].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid yMax '%2'").arg(linenum).arg(fields[field_yMax]));
      ev1.z() = fields[field_zMax].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid zMax '%2'").arg(linenum).arg(fields[field_zMax]));

      ev2.x() = fields[field_xMin].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid xMin '%2'").arg(linenum).arg(fields[field_xMin]));
      ev2.y() = fields[field_yMin].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid yMin '%2'").arg(linenum).arg(fields[field_yMin]));
      ev2.z() = fields[field_zMin].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid zMin '%2'").arg(linenum).arg(fields[field_zMin]));

      evals.x() = fields[field_k1].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid k1 '%2'").arg(linenum).arg(fields[field_k1]));
      evals.y() = fields[field_k2].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid k2 '%2'").arg(linenum).arg(fields[field_k2]));
      evals.z() = fields[field_k3].toFloat(&ok);
      if(!ok)
        return setErrorMessage(QString("Error line %1: invalid k3 '%2'").arg(linenum).arg(fields[field_k3]));

      // make sure ev1 and ev2 are unit vectors before storing in cellAxis. It is not the case in the older csv files
      ev1 = ev1 / norm(ev1);
      ev2 = ev2 / norm(ev2);

      cellAxis[label] = SymmetricTensor(ev1, ev2, evals);
    }

    mesh->setCellAxisType(type);

    // Display cell axis according to their type, with the parameters from the GUI

    if(type == "curvature") {
      DisplayTissueCurvature *proc;
      if(!getProcess("Mesh/Cell Axis/Curvature/Display Tissue Curvature", proc))
        throw QString("%1::run() Unable to make DisplayTissueCurvature process").arg(name());
      return proc->run();
    } else if(type == "PDG") {
      DisplayPDGs *proc;
      //if(!getProcess("Mesh/Cell Axis/Deformation Gradient/Display Vertex Gradient", proc))
      if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", proc))
        throw QString("%1::run() Unable to make DisplayPDGs process").arg(name());
			// DisplayPDGs::initialize() changes the display parms in case they commend to draw the PDGs projected
			// on a Bezier and cellAxisBezier doesn't exist yet.
  		// ALR: I don't know how to update the GUI with the new parms.
      proc->initialize(0);
      return proc->run();
    } else if(type == "fibril") {
	    DisplayFibrilOrientations *proc;
      if(!getProcess("Mesh/Cell Axis/Fibril Orientations/Display Fibril Orientations", proc))
        throw QString("%1::run() Unable to make DisplayFibrilOrientations process").arg(name());
      return proc->run();
    } else if(type == "polarization") {
	    DisplaySignalOrientation *proc;
      if(!getProcess("Mesh/Cell Axis/Polarization/Display Signal Orientation", proc))
        throw QString("%1::run() Unable to make DisplaSignalOrientation process").arg(name());
      return proc->run();
    } else if(type == "PCA") {
	    DisplayPCA *proc;
      if(!getProcess("Mesh/Cell Axis/Shape Analysis/Display Shape Axis 2D", proc))
        throw QString("%1::run() Unable to make DisplayPCA process").arg(name());
      return proc->run();
		}

    mesh->setShowAxis("Cell Axis");

    SETSTATUS(QString("Loaded directions for %1 cells").arg(cellAxis.size()));
    return true;
  }
  REGISTER_PROCESS(LoadCellAxis);

  bool ClearCellAxis::run(Mesh* mesh)
  {
    mesh->clearCellAxis();
    mesh->updateTriangles();
    return true;
  }
  REGISTER_PROCESS(ClearCellAxis);


  //Calculate the term for a single triangle of the integrated divergence of the neighborhood
  double triangleDivergence(const Point3d &p1, const Point3d &p2, const Point3d &p3,const Point3d &grad)
  {
        // We compute the summed quantity as
        // (-(e1 . e)/(e1 % e) (e2 . X_j) + (e2 . e)/(e2 % e) (e1 . X_j),
        // where e is the vector connecting the opposite vertices.
        Point3d e1 = p2 - p1, e2 = p3 - p1, e = p3 - p2;

        double cotan1 = (e2 * e) / norm(e2 % e), cotan2 = -(e1 * e) / norm(e1 % e);

        return (cotan1 * (e1 * grad) + cotan2 * (e2 * grad));
  }

float cotWeight(Point3d j,Point3d k,Point3d s,Point3d f){
  float weight = 0;
 Point3d js = j-s;
 Point3d ks = k-s;
//  Vec3f kj = k->pos-j->pos;
 Point3d jf = j-f;
 Point3d kf = k-f;
//  if(S.edge(k,s))//!j->border || !k->border || !s->border)
    weight+= js*ks/norm(js^(ks))/2.0;
//  if(S.edge(k,f))//!j->border || !k->border || !f->border)
    weight+=jf*kf/norm(jf^(kf))/2.0;
  return weight;
}


/*
 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomOrthogHeat::run(Mesh *mesh, bool projectOnSurface, QString normalize, bool display)
  {
  std::cout<<"CreateCustomOrthogHeatV1.1"<<std::endl;
  std::set<int> selectedCells = findAllSelectedLabels(mesh->graph());

  const QString attrMapName = "Custom Directions";
//  Mesh cell_mesh = (*mesh);
  std::cout<<"CreateCustomOrthogHeat: Creating cell graph"<<std::endl;
  vvGraph S;
  vvGraph J;
  cellGraph C;
if(!mesh->tissue().createTissueGraphs(S,J,C))
     throw(QString("Unexpected error converting to 2D cell tissue (MGX2D)"));
//  if(!cell_mesh.tissue().toMgx2d())
//     throw(QString("Unexpected error converting to 2D cell tissue (MGX2D)"));
//  vvGraph& S = mesh->graph();
// create MGXC if not already

 // cellGraph C = cell_mesh.tissue().C;


//  AttrMap<int, SymmetricTensor>& attrDataTest = mesh->attributes().attrMap<int, SymmetricTensor>(attrMapName);
//  if(attrDataTest.size() == 0){
//	CreateCustomDirectionsHeat *procToRun = new CreateCustomDirectionsHeat(*this);
//	procToRun->run(mesh, projectOnSurface,normalize,false);
//  }

//  AttrMap<int, SymmetricTensor> &attrData = mesh->attributes().attrMap<int, SymmetricTensor>(attrMapName);
  std::map<int,int> labToMat; //Label to matrix mapping
  std::map<int,int> matToLab; //Label to matrix mapping

//  AttrMap<int, SymmetricTensor>::iterator it = attrData.begin();
//  int mat_ind = 0;
//  while(it!=attrData.end()){
//  	int label = it->first;
//  	labToMat[label] = mat_ind;
//  	matToLab[mat_ind] = label;
//  	mat_ind++;
//  }

  std::cout<<"CreateCustomOrthogHeat: Computing matrix mappings"<<std::endl;


  AttrMap<int, float>::iterator it = mesh->labelHeat().begin();
  int mat_ind = 0;
  while(it!=mesh->labelHeat().end()){
  	int label = it->first;
  	labToMat[label] = mat_ind;
  	matToLab[mat_ind] = label;
  	mat_ind++;
  	it++;
  }

  //Extract cell-graph neigh (ala 6 coloring code)
//  std::map<IntIntPair, double> neighPair;
//  std::map<int, std::set<int> > neighMap;
//  neighborhood2D(S,neighPair);


  //Construct the matrix to solve (Poisson equation ala geodesics in heat), and solve (ala deformation code)
  //NEED CELL TO ROW/COLUMN MAPPING
  //Initialize
  std::cout<<"CreateCustomOrthogHeat: Initialize matrix"<<std::endl;

  int num_cells = mat_ind;
  std::vector<std::vector<double> > Mat;
  std::vector<double> bvec;
  Mat.resize(num_cells);
  bvec.resize(num_cells);
  for(uint i=0;i<num_cells;i++){
    Mat[i].resize(num_cells);
    bvec[i] = 0;
    for(uint ii=0;ii<num_cells;ii++)
    	Mat[i][ii] = 0;
  }


  //Fill
  std::cout<<"CreateCustomOrthogHeat: Fill matrix and b-vector"<<std::endl;
  forall(const cell &v, C){

  	  std::map<int,int>::iterator v_it = labToMat.find(v->label);
  	  int v_ind;
  	  double diagonal = 0;
  	  double area = 0;
	  	if(v_it == labToMat.end())
  			continue;
  		else
  			v_ind = v_it->second;

	  	forall(const cell &c, C.neighbors(v)){
			const cell &n = C.nextTo(v,c);
			const cell &p = C.prevTo(v,c);
	  		std::map<int,int>::iterator c_it = labToMat.find(c->label);
	  		std::map<int,int>::iterator n_it = labToMat.find(n->label);
	  		std::map<int,int>::iterator p_it = labToMat.find(p->label);

	  		if(c_it == labToMat.end() || n_it == labToMat.end())
	  			continue;
	  		int c_ind = c_it->second;
	  		int n_ind = n_it->second;

	  		//Add divergence to bvec
	  		double s_v,s_c,s_n;
	  		s_v = mesh->labelHeat()[v->label];
	  		s_c = mesh->labelHeat()[c->label];
	  		s_n = mesh->labelHeat()[n->label];
	  		Point3d t_grad = triangleGradient(v->pos,c->pos,n->pos,s_v,s_c,s_n);//s_x,s_c,s_n);
	  		double t_grad_len = norm(t_grad);
	//  		t_grad = t_grad%((c->pos-v->pos)%(n->pos-v->pos));
	  		t_grad /= norm(t_grad);
	  		//t_grad *= t_grad_len/ norm(t_grad);
	  		double t_area = triangleArea(v->pos,c->pos,n->pos);
//	  		std::cout<<"triangle area: "<<t_area<<std::endl;
	  		area+=t_area;

	  		bvec[v_ind]+= triangleDivergence(v->pos,c->pos,n->pos,t_grad);
//	  		std::cout<<"triangle divergence: "<<triangleDivergence(v->pos,c->pos,n->pos,t_grad)<<std::endl;;

	  		if(p_it == labToMat.end())
	  			continue;
	  		int p_ind = p_it->second;

	  		double cot_weight = cotWeight(v->pos,c->pos,n->pos,p->pos);//COT_WEIGHT(pos_v,pos_c,pos_n,pos_p);
//	  		std::cout<<"cot weight: "<<cot_weight<<std::endl;
	  		diagonal += cot_weight;
	  		Mat[v_ind][c_ind] = -cot_weight;
	  	}
		Mat[v_ind][v_ind] = diagonal;
//		std::cout<<bvec[v_ind]<<" pre-mult b-vec for "<<v_ind<<std::endl;
//		std::cout<<area<<std::endl;
		//bvec[v_ind] /= 2.0;
		bvec[v_ind] *= area/2.0;

  }
//  bvec[0] = 1;
//  for(int i=0;i<mat_ind;i++)
//  	Mat[0][i] = 0;
//  Mat[0][0] = 1;
//  for(int i=0;i<mat_ind;i++){
//  	for(int j=0;j<mat_ind;j++)
//  		std::cout<<Mat[i][j]<<" ";
//  	std::cout<<"   b   "<<bvec[i]<<std::endl;
//  }
for(int i=0;i<mat_ind;i++){
	bool non_zero = false;
	for(int j=0;j<mat_ind;j++)
		if(Mat[i][j]!=0)
			non_zero = true;
	int label = matToLab[i];
	if(!non_zero || selectedCells.find(label)!=selectedCells.end()){
		for(int j=0;j<mat_ind;j++)
			Mat[i][j] = 0;
		Mat[i][i] = 1;

		bvec[i] = 0;
	}

}



  std::cout<<"CreateCustomOrthogHeat: Solve matrix"<<std::endl;
  //Solve
  systemSolveGSL(Mat,bvec);

  //Copy values to cells, for heatmap and save the heatmap

  std::cout<<"CreateCustomOrthogHeat: Copy resulting values to heatmap"<<std::endl;

  it = mesh->labelHeat().begin();
  while(it!=mesh->labelHeat().end()){
  	int label = it->first;
  	mesh->labelHeat()[label] = bvec[labToMat[label]];
  	it++;
  }
  std::cout<<"CreateCustomOrthogHeat: Update and display heatmap"<<std::endl;
  mesh->updateAll();

  if(display){

  	//Display the reconstructed orthogonal heat-map

    // Run display, with the parameters from the GUI
//    QStringList parms;
//    DisplayCustomOrientations *proc = getProcessParms<DisplayCustomOrientations>(this, parms);
//    if(!proc)
//      throw(QString("Unable to create display custom orientation process"));
//    proc->run(parms);
  }
  std::cout<<"Returning now"<<std::endl;

  return true;
  }
  REGISTER_PROCESS(CreateCustomOrthogHeat);
*/

 bool CreateIntrinsicHeat::run(Mesh *mesh, bool projectOnSurface, QString normalize, bool display)
  {
  std::cout<<"CreateIntrinsicHeatV1.1"<<std::endl;
  std::set<int> selectedCells = findAllSelectedLabels(mesh->graph());

  const QString attrMapName = "Measure Label Tensor CustomDirections";
//  Mesh cell_mesh = (*mesh);
  std::cout<<"CreateIntrinsicHeat: Creating cell graph"<<std::endl;
  vvGraph S;
  vvGraph J;
  cellGraph C;
if(!mesh->tissue().createTissueGraphs(S,J,C))
     throw(QString("Unexpected error converting to 2D cell tissue (MGX2D)"));
  std::map<int,int> labToMat; //Label to matrix mapping
  std::map<int,int> matToLab; //Label to matrix mapping

  std::cout<<"CreateIntrinsicHeat: Computing matrix mappings"<<std::endl;


  AttrMap<int, float>::iterator it = mesh->labelHeat().begin();
  int mat_ind = 0;
  while(it!=mesh->labelHeat().end()){
  	int label = it->first;
  	labToMat[label] = mat_ind;
  	matToLab[mat_ind] = label;
  	mat_ind++;
  	it++;
  }

  //Extract cell-graph neigh (ala 6 coloring code)
//  std::map<IntIntPair, double> neighPair;
//  std::map<int, std::set<int> > neighMap;
//  neighborhood2D(S,neighPair);


  //Construct the matrix to solve (Poisson equation ala geodesics in heat), and solve (ala deformation code)
  //NEED CELL TO ROW/COLUMN MAPPING
  //Initialize
  std::cout<<"CreateIntrinsicHeat: Initialize matrix"<<std::endl;

  int num_cells = mat_ind;
  std::vector<std::vector<double> > Mat;
  std::vector<double> bvec;
  Mat.resize(num_cells);
  bvec.resize(num_cells);
  for(int i=0;i<num_cells;i++){
    Mat[i].resize(num_cells);
    bvec[i] = 0;
    for(int ii=0;ii<num_cells;ii++)
    	Mat[i][ii] = 0;
  }

        double max_signal = -1e20;
        double min_signal = 1e20;
      	forall(const cell &c, C){
            double s_c = mesh->labelHeat()[c->label];
            if(max_signal<s_c)
                max_signal=s_c;
            if(min_signal>s_c)
                min_signal=s_c;
        }



  //Fill
  std::cout<<"CreateIntrinsicHeat: Fill matrix and b-vector"<<std::endl;
  forall(const cell &v, C){

  	  std::map<int,int>::iterator v_it = labToMat.find(v->label);
  	  int v_ind;
  	  double diagonal = 0;
  	  double area = 0;
	  	if(v_it == labToMat.end())
  			continue;
  		else
  			v_ind = v_it->second;

	  	forall(const cell &c, C.neighbors(v)){
			const cell &n = C.nextTo(v,c);
			const cell &p = C.prevTo(v,c);
	  		std::map<int,int>::iterator c_it = labToMat.find(c->label);
	  		std::map<int,int>::iterator n_it = labToMat.find(n->label);
	  		std::map<int,int>::iterator p_it = labToMat.find(p->label);

	  		if(c_it == labToMat.end() || n_it == labToMat.end())
	  			continue;
	  		int c_ind = c_it->second;
	  		//int n_ind = n_it->second;

	  		//Add divergence to bvec
	  		double s_v,s_c,s_n;
	  		s_v = mesh->labelHeat()[v->label];
	  		s_c = mesh->labelHeat()[c->label];
	  		s_n = mesh->labelHeat()[n->label];
	  		Point3d t_grad = triangleGradient(v->pos,c->pos,n->pos,s_v,s_c,s_n);//s_x,s_c,s_n);
//	  		double t_grad_len = norm(t_grad);
	//  		t_grad = t_grad%((c->pos-v->pos)%(n->pos-v->pos));
	  		t_grad /= norm(t_grad);
	  		//t_grad *= t_grad_len/ norm(t_grad);
	  		double t_area = triangleArea(v->pos,c->pos,n->pos);
//	  		std::cout<<"triangle area: "<<t_area<<std::endl;
	  		area+=t_area;
            bvec[v_ind]=0;
//	  		bvec[v_ind]+= triangleDivergence(v->pos,c->pos,n->pos,t_grad);
//	  		std::cout<<"triangle divergence: "<<triangleDivergence(v->pos,c->pos,n->pos,t_grad)<<std::endl;;

	  		if(p_it == labToMat.end())
	  			continue;
//	  		int p_ind = p_it->second;

	  		double cot_weight = cotWeight(v->pos,c->pos,n->pos,p->pos);//COT_WEIGHT(pos_v,pos_c,pos_n,pos_p);
//	  		std::cout<<"cot weight: "<<cot_weight<<std::endl;
	  		diagonal += cot_weight;
	  		Mat[v_ind][c_ind] = -cot_weight;
	  	}
		Mat[v_ind][v_ind] = diagonal;
//		std::cout<<bvec[v_ind]<<" pre-mult b-vec for "<<v_ind<<std::endl;
//		std::cout<<area<<std::endl;
		//bvec[v_ind] /= 2.0;
//		bvec[v_ind] *= area/2.0;

  }
//  bvec[0] = 1;
//  for(int i=0;i<mat_ind;i++)
//  	Mat[0][i] = 0;
//  Mat[0][0] = 1;
//  for(int i=0;i<mat_ind;i++){
//  	for(int j=0;j<mat_ind;j++)
//  		std::cout<<Mat[i][j]<<" ";
//  	std::cout<<"   b   "<<bvec[i]<<std::endl;
//  }
for(int i=0;i<mat_ind;i++){
	bool non_zero = false;
	for(int j=0;j<mat_ind;j++)
		if(Mat[i][j]!=0)
			non_zero = true;
	int label = matToLab[i];
    double s_c = mesh->labelHeat()[label];
	if(!non_zero || s_c == min_signal){
		for(int j=0;j<mat_ind;j++)
			Mat[i][j] = 0;
		Mat[i][i] = 1;

		bvec[i] = 0;
	}
    if(s_c == max_signal){
		for(int j=0;j<mat_ind;j++)
			Mat[i][j] = 0;
		Mat[i][i] = 1;
		bvec[i] = 1;

    }

}



  std::cout<<"CreateCustomOrthogHeat: Solve matrix"<<std::endl;
  //Solve
  systemSolveGSL(Mat,bvec);

  //Copy values to cells, for heatmap and save the heatmap

  std::cout<<"CreateCustomOrthogHeat: Copy resulting values to heatmap"<<std::endl;

  it = mesh->labelHeat().begin();
  while(it!=mesh->labelHeat().end()){
  	int label = it->first;
  	mesh->labelHeat()[label] = bvec[labToMat[label]];
  	it++;
  }
  std::cout<<"CreateCustomOrthogHeat: Update and display heatmap"<<std::endl;
  mesh->updateAll();

  if(display){

  	//Display the reconstructed orthogonal heat-map

    // Run display, with the parameters from the GUI
//    QStringList parms;
//    DisplayCustomOrientations *proc = getProcessParms<DisplayCustomOrientations>(this, parms);
//    if(!proc)
//      throw(QString("Unable to create display custom orientation process"));
//    proc->run(parms);
  }
  std::cout<<"Returning now"<<std::endl;

  return true;
  }
  REGISTER_PROCESS(CreateIntrinsicHeat);

 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsHeat::run(Mesh *mesh, bool projectOnSurface, QString normalize, bool display)
  {

  vvGraph& S = mesh->graph();

  AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

  // create neighbor graph
  std::map<int, std::set<int> > cellNeighborMap;
  std::map<IntInt, double> neighMap;
  neighborhood2D(S, neighMap, "");

  forall(auto p, neighMap){
    cellNeighborMap[p.first.first].insert(p.first.second);
  }


  // Update the normals and center of the cell for visualisation.
  mesh->updateCentersNormals();

  // cell centers
  IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
  IntPoint3fAttr normals = (mesh->useParents()? mesh->parentNormal():mesh->labelNormal());

  // merge neighbor graph
  std::map<int, std::set<int> > cellNeighborMapParents;
  if(mesh->useParents()){
    std::set<int> meshLabels = findAllLabelsSet(S);
    // go through labels
    forall(int l, meshLabels){
      forall(int neighb, cellNeighborMap[l]){
        if(mesh->parents()[l] != mesh->parents()[neighb])
          cellNeighborMapParents[mesh->parents()[l]].insert(mesh->parents()[neighb]);
      }
    }
    cellNeighborMap = cellNeighborMapParents;
  }

  // for each cell: check heat difference to all of its neighbors.
  forall(auto p, centers){
    Point3d dir(0,0,0);

    int cell = p.first; // cell label
    Point3d cellCenter(centers[cell]);
    Point3d cellNormal(normals[cell]);

    forall(int neighb, cellNeighborMap[p.first]){
      double dif = mesh->labelHeat()[p.first] - mesh->labelHeat()[neighb];
      Point3d dirN = cellCenter-Point3d(centers[neighb]);
      if(normalize == "norm"){
        dirN /= norm(dirN);
      } else if(normalize == "square norm"){
        dirN /= norm(dirN)*norm(dirN);
      }
      dir += dirN * dif;
    }

    dir /= norm(dir);
    attrData[p.first].ev1() = Point3f(dir);

    if(projectOnSurface){
      Point3d projP = projectPointOnPlane(cellCenter+dir,cellCenter,cellNormal);

      attrData[p.first].ev1() = Point3f(projP - cellCenter);
    }
    attrData[p.first].ev2() = attrData[p.first].ev1() % Point3f(cellNormal);
  }

  if(display) {
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations *proc;
    if(!getProcess("Mesh/Cell Axis/Custom/Display Custom Orientations", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }
  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsHeat);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool SmoothCustomDirections::run(Mesh *m, bool weight, bool projectOnSurface, bool display)
  {

  const QString attrMapName = "Measure Label Tensor CustomDirections";
  vvGraph& S = m->graph();

  AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  // create neighbor graph
  std::map<int, std::set<int> > cellNeighborMap;
  std::map<IntInt, double> neighMap;
  neighborhood2D(S, neighMap, "");

  forall(auto p, neighMap){
    cellNeighborMap[p.first.first].insert(p.first.second);
  }

  MeasureArea *ma;
  if(!getProcess("Mesh/Heat Map/Measures/Geometry/Area", ma))
    throw(QString("MeshProcessHeatMap:: Unable to make MeasureArea process"));
  IntFloatAttr cellAreas;
  ma->calculateArea(S, cellAreas);


  // Update the normals and center of the cell for visualisation.
  m->updateCentersNormals();

  // cell centers
  IntPoint3fAttr centers = (m->useParents()? m->parentCenterVis():m->labelCenterVis());
  IntPoint3fAttr normals = (m->useParents()? m->parentNormal():m->labelNormal());

  // merge neighbor graph
  std::map<int, std::set<int> > cellNeighborMapParents;
  IntFloatAttr cellAreasParents;
  if(m->useParents()){
    std::set<int> meshLabels = findAllLabelsSet(S);
    // go through labels
    forall(int l, meshLabels){
      forall(int neighb, cellNeighborMap[l]){
        if(m->parents()[l] != m->parents()[neighb])
          cellNeighborMapParents[m->parents()[l]].insert(m->parents()[neighb]);
      }

      cellAreasParents[m->parents()[l]] += cellAreas[l];
    }
    cellNeighborMap = cellNeighborMapParents;
    cellAreas = cellAreasParents;
  }
  if(!weight){
    forall(auto p, centers){
      cellAreas[p.first] = 1;
    }
  }

  // for each cell: average across neighbors
  std::map<int, SymmetricTensor> averagedDirections;
  std::set<int> labelsToCopy;
  forall(auto p, centers){

    int cellLabel = p.first;

    double totalArea = cellAreas[cellLabel];

    if(attrData.find(cellLabel) == attrData.end()) continue;

    averagedDirections[cellLabel].ev1() = cellAreas[cellLabel] * attrData[cellLabel].ev1();
    averagedDirections[cellLabel].ev2() = cellAreas[cellLabel] * attrData[cellLabel].ev2();
    averagedDirections[cellLabel].evals() = cellAreas[cellLabel] * attrData[cellLabel].evals();

    forall(int nLabel, cellNeighborMap[cellLabel]){
      if(attrData.find(nLabel) == attrData.end()) continue;
      // check for nan
      if(averagedDirections[cellLabel].ev1().x() != averagedDirections[cellLabel].ev1().x()) continue;
      if(averagedDirections[cellLabel].ev2().x() != averagedDirections[cellLabel].ev2().x()) continue;
      if(averagedDirections[cellLabel].evals().x() != averagedDirections[cellLabel].evals().x()) continue;

      totalArea += cellAreas[nLabel];
      averagedDirections[cellLabel].ev1() += cellAreas[nLabel] * attrData[nLabel].ev1();
      averagedDirections[cellLabel].ev2() += cellAreas[nLabel] * attrData[nLabel].ev2();
      averagedDirections[cellLabel].evals() += cellAreas[nLabel] * attrData[nLabel].evals();
    }

    if(norm(totalArea) < 1E-5) continue;

    averagedDirections[cellLabel].ev1() /= totalArea;
    averagedDirections[cellLabel].ev2() /= totalArea;
    averagedDirections[cellLabel].evals() /= totalArea;

    averagedDirections[cellLabel].ev1() /= norm(averagedDirections[cellLabel].ev1());
    averagedDirections[cellLabel].ev2() /= norm(averagedDirections[cellLabel].ev2());
    averagedDirections[cellLabel].evals() /= norm(averagedDirections[cellLabel].evals());

    labelsToCopy.insert(cellLabel);
  }

  if(projectOnSurface){
    forall(auto p, centers){
      int cellLabel = p.first;
      Point3d cellCenter(centers[cellLabel]);
      Point3d cellNormal(normals[cellLabel]);

      Point3d projP = projectPointOnPlane(cellCenter+Point3d(averagedDirections[cellLabel].ev1()),cellCenter,cellNormal);
      averagedDirections[cellLabel].ev1() = Point3f(projP - cellCenter);
      averagedDirections[cellLabel].ev1() /= norm(averagedDirections[cellLabel].ev1());

      projP = projectPointOnPlane(cellCenter+Point3d(averagedDirections[cellLabel].ev2()),cellCenter,cellNormal);
      averagedDirections[cellLabel].ev2() = Point3f(projP - cellCenter);
      averagedDirections[cellLabel].ev2() /= norm(averagedDirections[cellLabel].ev2());

      projP = projectPointOnPlane(cellCenter+Point3d(averagedDirections[cellLabel].evals()),cellCenter,cellNormal);
      averagedDirections[cellLabel].evals() = Point3f(projP - cellCenter);
      averagedDirections[cellLabel].evals() /= norm(averagedDirections[cellLabel].evals());
    }
  }


  forall(int l, labelsToCopy){
    attrData[l] = averagedDirections[l];
  }




  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations *proc;
    if(!getProcess("Mesh/Cell Axis/Custom/Display Custom Orientations", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }
  return true;
  }
  REGISTER_PROCESS(SmoothCustomDirections);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CopyCustomDirections::run(Mesh *m1, Mesh *m2)
  {

  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData1 = m1->attributes().attrMap<int, SymmetricTensor>(attrMapName);  

  AttrMap<int, SymmetricTensor>& attrData2 = m2->attributes().attrMap<int, SymmetricTensor>(attrMapName);
  attrData2.clear();

  forall(auto p, attrData1){
    attrData2[p.first] = p.second;
  }

  return true;
  }
  REGISTER_PROCESS(CopyCustomDirections);


 // TBD
 bool ImportCustomDirections::run(Mesh *m1, Mesh *m2, QString prefix, QString name, QString nameMesh)
  {


  QString attrMapName = prefix + name;
  AttrMap<int, SymmetricTensor>& attrData = m1->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  if(attrData.empty()){
    return setErrorMessage("There is no Tensor Attr Map with this name");
  }

  qDebug() << "Load Tensor Attr Map: " << attrMapName << " of size: " << attrData.size() << "/n";

  



  if(nameMesh == "Active Mesh"){
    AttrMap<int, SymmetricTensor>& customDirs = m1->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    customDirs.clear();

    forall(auto p, attrData){
      customDirs[p.first] = p.second;
    }
  } else {
    AttrMap<int, SymmetricTensor>& customDirs = m2->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    customDirs.clear();

    forall(auto p, attrData){
      int l = p.first;
      if(m2->useParents()) l = m2->parents()[p.first];
      customDirs[l] = p.second;
    }
  }

  return true;
  }
  REGISTER_PROCESS(ImportCustomDirections);

 // TBD
 bool ImportCustomDirectionsVector::run(Mesh *m1, Mesh *m2, QString prefix, QString name, QString dimension, QString nameMesh)
  {


  QString attrMapName = prefix + name;
  AttrMap<int, Point3d>& attrData = m1->attributes().attrMap<int, Point3d>(attrMapName);

  if(attrData.empty()){
    return setErrorMessage("There is no Vector Attr Map with this name");
  }
  qDebug() << "Load Vector Attr Map: " << attrMapName << " of size: " << attrData.size() << "/n";
 
  if(nameMesh == "Active Mesh"){
    AttrMap<int, SymmetricTensor>& customDirs = m1->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    customDirs.clear();

    forall(auto p, attrData){
      if(dimension == "X")
        customDirs[p.first].ev1() = Point3f(p.second);
      else if(dimension == "Y")
        customDirs[p.first].ev2() = Point3f(p.second);
      else
        customDirs[p.first].evals() = Point3f(p.second);
    }
  } else {
    AttrMap<int, SymmetricTensor>& customDirs = m2->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    customDirs.clear();

    forall(auto p, attrData){
      int l = p.first;
      if(m2->useParents()) l = m2->parents()[p.first];
      if(dimension == "X")
        customDirs[l].ev1() = Point3f(p.second);
      else if(dimension == "Y")
        customDirs[l].ev2() = Point3f(p.second);
      else
        customDirs[l].evals() = Point3f(p.second);
    }
  }

  return true;
  }
  REGISTER_PROCESS(ImportCustomDirectionsVector);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsSurface::run(Mesh *m, Mesh *m2, QString customX, QString customY, QString customZ, bool display)
  {

  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  vvGraph& S2 = m2->graph();

  // get cell centroids
  AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

  if(m->useParents()){
    centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
    std::cout << "using the parent labels " << centroids.size() << std::endl;
  }

  if(centroids.empty()){
    return setErrorMessage("No valid cell centroids found!");
  }

  // get nearest surf points
  forall(auto p, centroids){
    vertex sv = nearestVertexOnMesh(p.second, S2);
    Point3d surfDir = sv->pos - p.second;
    surfDir/=norm(surfDir);

    Point3d surfNormal = sv->nrml;

    if(customX == "Surface Direction")
      attrData[p.first].ev1() = Point3f(surfDir);
    else if(customX == "Surface Normal")
      attrData[p.first].ev1() = Point3f(surfNormal);

    if(customY == "Surface Direction")
      attrData[p.first].ev2() = Point3f(surfDir);
    else if(customY == "Surface Normal")
      attrData[p.first].ev2() = Point3f(surfNormal);

    if(customZ == "Surface Direction")
      attrData[p.first].evals() = Point3f(surfDir);
    else if(customZ == "Surface Normal")
      attrData[p.first].evals() = Point3f(surfNormal);

    std::cout << "l " << p.first << "/" << p.second << "/" << surfDir << "/" << attrData[p.first].ev1() << std::endl;
  }

  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Custom/Display Custom Orientations", proc))
      throw(QString("CreateCustomDirectionsSurface:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsSurface);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CustomDirectionsEnforceOrthogonality::run(Mesh *m, QString primX, QString secX, bool changeSecondary, bool display)
  {

  if(primX == secX){
    return setErrorMessage("Primary and Secondary Axis must be different!");
  }

  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  if(attrData.empty()){
    return setErrorMessage("Couldn't find a Custom Directions Attribute Map. Please create Custom Directions first!");
  }

  forall(auto p, attrData){
    Point3d primary;
    if(primX == "X") primary = Point3d(p.second.ev1());
    else if(primX == "Y") primary = Point3d(p.second.ev2());
    else if(primX == "Z") primary = Point3d(p.second.evals());

    if(norm(primary) < 1E-5) continue;

    Point3d secondary;
    if(secX == "X") secondary = Point3d(p.second.ev1());
    else if(secX == "Y") secondary = Point3d(p.second.ev2());
    else if(secX == "Z") secondary = Point3d(p.second.evals());

    if(norm(secondary) < 1E-5) continue;

    Point3d tertiary = primary % secondary;
    tertiary /= norm(tertiary);

    if((primX == "X" and secX == "Y") or (primX == "Y" and secX == "X")) attrData[p.first].evals() = Point3f(tertiary);
    else if((primX == "X" and secX == "Z") or (primX == "Z" and secX == "X")) attrData[p.first].ev2() = Point3f(tertiary);
    else if((primX == "Z" and secX == "Y") or (primX == "Y" and secX == "Z")) attrData[p.first].ev1() = Point3f(tertiary);

    if(!changeSecondary) continue;

    secondary = primary % tertiary;
    if(secX == "X") attrData[p.first].ev1() = Point3f(secondary);
    else if(secX == "Y") attrData[p.first].ev2() = Point3f(secondary);
    else if(secX == "Z") attrData[p.first].evals() = Point3f(secondary);

  }


  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Custom/Display Custom Orientations", proc))
      throw(QString("CustomDirectionsEnforceOrthogonality:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CustomDirectionsEnforceOrthogonality);



 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CustomDirectionsClear::run(Mesh *m, bool clearX, bool clearY, bool clearZ, bool display)
  {

  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  if(clearX and clearY and clearZ){
    attrData.clear();
    return true;
  }

  forall(auto p, attrData){
    Point3f zeroV(0.,0.,0.);

    if(clearX) attrData[p.first].ev1() = zeroV;
    if(clearY) attrData[p.first].ev2() = zeroV;
    if(clearZ) attrData[p.first].evals() = zeroV;
  }

  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Custom/Display Custom Orientations", proc))
      throw(QString("CustomDirectionsClear:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CustomDirectionsClear);


 bool ProjectCustomDirections::run(Mesh *m,QString makeOrthog)
  {
    const QString attrMapName = "Measure Label Tensor CustomDirections";
    if(m->attributes().attrMap<int, SymmetricTensor>(attrMapName).empty()){
        setErrorMessage(QString("No custom directions stored in active mesh!"));
        return false;
    }

    AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

    // Update the normals and center of the cell for visualisation.
    m->updateCentersNormals();

    // cell centers
    IntPoint3fAttr centers = (m->useParents()? m->parentCenterVis():m->labelCenterVis());
    IntPoint3fAttr normals = (m->useParents()? m->parentNormal():m->labelNormal());
    
    forall(auto p, centers){
      int cell = p.first; // cell label
      Point3d cellCenter(centers[cell]);
      Point3d cellNormal(normals[cell]);

      Point3d dir1(attrData[p.first].ev1());
      Point3d dir2(attrData[p.first].ev2());

      Point3d projP = projectPointOnPlane(cellCenter+dir1,cellCenter,cellNormal);
      attrData[p.first].ev1() = Point3f(projP - cellCenter);
      projP = projectPointOnPlane(cellCenter+dir2,cellCenter,cellNormal);
      attrData[p.first].ev2() = Point3f(projP - cellCenter);


      if(makeOrthog == "Y from X")
          attrData[p.first].ev2() = attrData[p.first].ev1() % Point3f(cellNormal);
      if(makeOrthog == "X from Y")
          attrData[p.first].ev1() = attrData[p.first].ev2() % Point3f(cellNormal);
    }


  }
  REGISTER_PROCESS(ProjectCustomDirections);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsBezier::run(Mesh *m, QString customX, QString customY, bool projectOnSurface, bool display)
  {

    const QString attrMapName = "Measure Label Tensor CustomDirections";
    AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

    // get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    int dataPointsBezier = 200; // number of intermediate points
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    Bezier b = cutSurf->bezier();

    std::vector<std::vector<Point3d> > bezGrid;

    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    // Update the normals and center of the cell for visualisation.
    m->updateCentersNormals();

    // cell centers
    IntPoint3fAttr centers = (m->useParents()? m->parentCenterVis():m->labelCenterVis());
    IntPoint3fAttr normals = (m->useParents()? m->parentNormal():m->labelNormal());

    // maps to store cell label -> local bezier vector
    IntP3dMap cellBezVecMapX, cellBezVecMapY;

    // loop through cell centers
    forall(auto p, centers){
      int cell = p.first; // cell label
      Point3d cellCenter(centers[cell]);
      Point3d cellNormal(normals[cell]);

      // find nearest bez grid point and take its derivatives
      Point2i idx;
      Point3d bezInterpPoint = calcNearestPointOnBezierGrid(Point3d(p.second), bezGrid, idx);

      if(idx.x() == dataPointsBezier-1) idx.x()--;
      if(idx.x() == 0) idx.x()++;
      if(idx.y() == dataPointsBezier-1) idx.y()--;
      if(idx.y() == 0) idx.y()++;

      Point3d bezDerivX = bezGrid[idx.x()+1][idx.y()] - bezGrid[idx.x()-1][idx.y()];
      bezDerivX/=norm(bezDerivX);
      Point3d bezDerivY = bezGrid[idx.x()][idx.y()+1] - bezGrid[idx.x()][idx.y()-1];
      bezDerivY/=norm(bezDerivY);
Information::out << "Bez Deriv X:" << bezDerivX << " Bez Deriv Y:" << bezDerivY << endl;

      if(projectOnSurface) {
        Point3d bez, bezPlusDifX, bezPlusDifY;

        bez = projectPointOnPlane(bezInterpPoint,cellCenter,cellNormal);
        bezPlusDifX = projectPointOnPlane(bezInterpPoint + bezDerivX,cellCenter,cellNormal) - bez;
        bezPlusDifY = projectPointOnPlane(bezInterpPoint + bezDerivY,cellCenter,cellNormal) - bez;

        cellBezVecMapX[cell] = bezPlusDifX/norm(bezPlusDifX);
        cellBezVecMapY[cell] = bezPlusDifY/norm(bezPlusDifY);
      }

      if(customX == "Bezier X")
        attrData[cell].ev1() = Point3f(cellBezVecMapX[cell]);
      else if(customX == "Bezier Y")
        attrData[cell].ev1() = Point3f(cellBezVecMapY[cell]);

      if(customY == "Bezier X")
        attrData[cell].ev2() = Point3f(cellBezVecMapX[cell]);
      else if(customY == "Bezier Y")
        attrData[cell].ev2() = Point3f(cellBezVecMapY[cell]);

    }
    if(display) {
      // Run display, with the parameters from the GUI
      DisplayCustomOrientations *proc;
      if(!getProcess("Mesh/Cell Axis/Custom/Display Custom Orientations", proc))
        throw(QString("MeshProcessHeatMap:: Unable to make DisplayCustomOrientations process"));
      proc->run();
    }

  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsBezier);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsBezierLine::run(Mesh *m, QString customX, QString customY, bool projectOnSurface, bool display)
   {

    const QString attrMapName = "Measure Label Tensor CustomDirections";
    AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

    // get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    int dataPointsBezier = 200; // number of intermediate points
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    Bezier b = cutSurf->bezier();

    std::vector<std::vector<Point3d> > bezGrid;
    std::vector<Point3d> discretizedBez, differentialBez;
    double totalLength;

    //b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);
    b.discretizeLineEqualWeight(dataPointsBezier, 0.5, mGLTot, discretizedBez, differentialBez, totalLength);

    // Update the normals and center of the cell for visualisation.
    m->updateCentersNormals();

    // cell centers
    IntPoint3fAttr centers = (m->useParents()? m->parentCenterVis():m->labelCenterVis());
    IntPoint3fAttr normals = (m->useParents()? m->parentNormal():m->labelNormal());

    // maps to store cell label -> local bezier vector
    IntP3dMap cellBezVecMapX, cellBezVecMapY;

    // loop through cell centers
    forall(auto p, centers){
      int cell = p.first; // cell label
      Point3d cellCenter(centers[cell]);
      Point3d cellNormal(normals[cell]);

      // find nearest bez grid point and take its derivatives
      //Point2i idx;
      //Point3d bezInterpPoint = calcNearestPointOnBezierGrid(Point3d(p.second), bezGrid, idx);

      int bezIdx;
      double bezWeight;
      Point3d bezInterpPoint = calcNearestPointOnBezierLine(Point3d(p.second), discretizedBez, bezIdx, bezWeight);

      Point3d bezDeriv = differentialBez[bezIdx];
      bezDeriv/=norm(bezDeriv);

      Point3d radial, longit, circ;

      radial = bezInterpPoint - cellCenter;
      radial/=norm(radial);

      longit = bezDeriv;

      circ = radial % longit;

      if(projectOnSurface){
        Point3d bez, bezPlusDifX, bezPlusDifY;

        Point3d longitP = projectPointOnPlane(cellCenter + longit,cellCenter,cellNormal) - cellCenter;
        longitP /= norm(longitP);

        Point3d radialP = projectPointOnPlane(cellCenter + radial,cellCenter,cellNormal) - cellCenter;
        radialP /= norm(radialP);

        Point3d circP = projectPointOnPlane(cellCenter + circ,cellCenter,cellNormal) - cellCenter;
        circP /= norm(circP);

        circ = circP;
        radial = radialP;
        longit = longitP;
      }

      if(customX == "Radial")
        attrData[cell].ev1() = Point3f(radial);
      else if(customX == "Circumferential")
        attrData[cell].ev1() = Point3f(circ);
      else if(customX == "Longitudinal")
        attrData[cell].ev1() = Point3f(longit);

      if(customY == "Radial")
        attrData[cell].ev2() = Point3f(radial);
      else if(customY == "Circumferential")
        attrData[cell].ev2() = Point3f(circ);
      else if(customY == "Longitudinal")
        attrData[cell].ev2() = Point3f(longit);

    }

  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations *proc;
    if(!getProcess("Mesh/Cell Axis/Custom/Display Custom Orientations", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsBezierLine);


 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsBezierLine3D::run(Mesh *m, QString customX, QString customY, QString customZ, bool display)
   {

    const QString attrMapName = "Measure Label Tensor CustomDirections";
    AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);

    AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    if(m->useParents()){
      centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
      std::cout << "using the parent labels " << centroids.size() << std::endl;
    }

    if(centroids.empty()){
      return setErrorMessage("No valid cell centroids found!");
    }


    // get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    int dataPointsBezier = 200; // number of intermediate points
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    Bezier b = cutSurf->bezier();

    std::vector<std::vector<Point3d> > bezGrid;
    std::vector<Point3d> discretizedBez, differentialBez;
    double totalLength;

    //b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);
    b.discretizeLineEqualWeight(dataPointsBezier, 0.5, mGLTot, discretizedBez, differentialBez, totalLength);

    // maps to store cell label -> local bezier vector
    IntP3dMap cellBezVecMapX, cellBezVecMapY;

    // loop through cell centers
    forall(auto p, centroids){
      int cell = p.first; // cell label
      Point3d cellCenter = p.second;

      // find nearest bez grid point and take its derivatives
      //Point2i idx;
      //Point3d bezInterpPoint = calcNearestPointOnBezierGrid(Point3d(p.second), bezGrid, idx);

      int bezIdx;
      double bezWeight;
      Point3d bezInterpPoint = calcNearestPointOnBezierLine(Point3d(p.second), discretizedBez, bezIdx, bezWeight);

      Point3d bezDeriv = differentialBez[bezIdx];
      bezDeriv/=norm(bezDeriv);

      Point3d radial, longit, circ;

      radial = bezInterpPoint - cellCenter;
      radial/=norm(radial);

      longit = bezDeriv;

      circ = radial % longit;

      if(customX == "Radial")
        attrData[cell].ev1() = Point3f(radial);
      else if(customX == "Circumferential")
        attrData[cell].ev1() = Point3f(circ);
      else if(customX == "Longitudinal")
        attrData[cell].ev1() = Point3f(longit);

      if(customY == "Radial")
        attrData[cell].ev2() = Point3f(radial);
      else if(customY == "Circumferential")
        attrData[cell].ev2() = Point3f(circ);
      else if(customY == "Longitudinal")
        attrData[cell].ev2() = Point3f(longit);

      if(customZ == "Radial")
        attrData[cell].evals() = Point3f(radial);
      else if(customZ == "Circumferential")
        attrData[cell].evals() = Point3f(circ);
      else if(customZ == "Longitudinal")
        attrData[cell].evals() = Point3f(longit);
    }

  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Custom/Display Custom Orientations", proc))
      throw(QString("CreateCustomDirectionsBezierLine3D:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsBezierLine3D);

 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CreateCustomDirectionsBezier3D::run(Mesh *m, QString customX, QString customY, QString customZ, bool display)
  {

    AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

  if(m->useParents()){
    centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
    std::cout << "using the parent labels " << centroids.size() << std::endl;
  }

    if(centroids.empty()){
      return setErrorMessage("No valid cell centroids found!");
    }


  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = m->attributes().attrMap<int, SymmetricTensor>(attrMapName);
  attrData.clear();

// get the bezier grid
    CuttingSurface* cutSurf = cuttingSurface();
    P2iP3dMap bezGridMap, diffBezGridMapX, diffBezGridMapY;
    double bezStart = 0.0;
    double bezEnd = 1.0;
    int dataPointsBezier = 200; // number of intermediate points
    double stepSize = (double)(bezEnd-bezStart)/(double)(dataPointsBezier-1);
    // correct directions, take into account the rotations/translations
    Matrix4d rotMatrixS1, rotMatrixCS;
    m->stack()->getFrame().getMatrix(rotMatrixS1.data());
    //currentMesh()->stack()->getFrame().getMatrix(rotMatrixS1.data());
    cutSurf->frame().getMatrix(rotMatrixCS.data());
    Matrix4d mGLTot = transpose(inverse(rotMatrixS1)) * transpose(rotMatrixCS);

    // create bezier map
    for(int i = 0; i < dataPointsBezier; i++){
      double u = bezStart + i*stepSize;
      for(int j = 0; j < dataPointsBezier; j++){
        double v = bezStart+j*stepSize;
        Point3d p(cutSurf->evalCoord(v,u)); // get bez world coordinate
        p = multMatrix4Point3(mGLTot,p); // correct rotations
        Point2i a(j,i);
        bezGridMap[a] = p;
      }
    }

    // create derivatives in X and Y direction of the bezier
    for(int i=1; i<dataPointsBezier-1; i++){
      for(int j=1; j<dataPointsBezier-1; j++){
        Point2i a (j,i);
        Point2i a1 (j-1,i);
        Point2i a2 (j+1,i);
        Point2i aY1 (j,i-1);
        Point2i aY2 (j,i+1);
        diffBezGridMapX[a] = bezGridMap[a2] - bezGridMap[a1];
        diffBezGridMapY[a] = bezGridMap[aY2] - bezGridMap[aY1];
      }
    }

  forall(auto p, centroids){

    int cell = p.first; // cell label
    Point3d cellCenter(centroids[cell]);
    //Point3d cellNormal(normals[cell]);
    // this code comes from CellNetworks
    double minDis = 1E20;
    Point3d bezInterpPoint (0,0,0);
    Point3d diffBezInterpPointX (0,0,0);
    Point3d diffBezInterpPointY (0,0,0);
    // find nearest bez grid point and take its derivatives
    for(int j=1; j<dataPointsBezier-1; j++){
      for(int k=1; k<dataPointsBezier-1; k++){
        Point2i c(k,j);
        double currentDistance = norm(bezGridMap[c] - cellCenter);
        if(currentDistance < minDis){
          minDis = currentDistance;
          bezInterpPoint = bezGridMap[c];
          diffBezInterpPointX = diffBezGridMapX[c];

          diffBezInterpPointY = diffBezGridMapY[c];

          if(norm(diffBezInterpPointX) < 1E-5){
            std::cout << "x " << diffBezInterpPointX << std::endl;
            diffBezInterpPointX = cellCenter - bezInterpPoint;

          } else if(norm(diffBezInterpPointY) < 1E-5){
            std::cout << "y " << diffBezInterpPointY << std::endl;
            diffBezInterpPointY = cellCenter - bezInterpPoint;
          }

          //coordBezier[cell] = Point3f(c.x(), c.y(), 0);
        }
      }



    }
    // normalize derivative
    Point3f bezX = Point3f(diffBezInterpPointX/norm(diffBezInterpPointX));
    Point3f bezY = Point3f(diffBezInterpPointY/norm(diffBezInterpPointY));
    Point3f bezN = bezX % bezY;
    bezN /= norm(bezN);

    if(customX == "Bezier X")
      attrData[cell].ev1() = bezX;
    else if(customX == "Bezier Y")
      attrData[cell].ev1() = bezY;
    else if(customX == "Bezier Normal")
      attrData[cell].ev1() = bezN;

    if(customY == "Bezier X")
      attrData[cell].ev2() = bezX;
    else if(customY == "Bezier Y")
      attrData[cell].ev2() = bezY;
    else if(customY == "Bezier Normal")
      attrData[cell].ev2() = bezN;

    if(customZ == "Bezier X")
      attrData[cell].evals() = bezX;
    else if(customZ == "Bezier Y")
      attrData[cell].evals() = bezY;
    else if(customZ == "Bezier Normal")
      attrData[cell].evals() = bezN;

    //std::cout << "cust " << cell << "/" << attrData[cell].ev1() << "/" << attrData[cell].ev2() << std::endl;

  }

  if(display){
    // Run display, with the parameters from the GUI
    DisplayCustomOrientations3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Custom/Display Custom Orientations", proc))
      throw(QString("CreateCustomDirectionsBezier3D:: Unable to make DisplayCustomOrientations process"));
    proc->run();
  }

  return true;
  }
  REGISTER_PROCESS(CreateCustomDirectionsBezier3D);

 bool CustomDirectionsAngle::run(Mesh *mesh, QString dirAxis, QString dirCustom)
  {
  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  if(attrData.empty())
    return setErrorMessage("No Custom Directions in the current mesh.");

  // if there is no cell axis stored, display error message
  const IntSymTensorAttr& cellAxis = mesh->cellAxis();
  if(cellAxis.size() == 0){
    setErrorMessage(QString("No cell axis stored in active mesh!"));
    return false;
  }

  // Update the normals and center of the cell for visualisation.
  mesh->updateCentersNormals();

  // cell centers
  IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
  IntPoint3fAttr normals = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());

  std::vector<double> vecValues;


  mesh->labelHeat().clear();
  forall(const IntSymTensorPair& p, cellAxis){
    Point3f axis1, axis2;
    if(dirAxis == "Max"){
      axis1 = p.second.ev1();
    } else {
      axis1 = p.second.ev2();
    }
    if(dirCustom == "X"){
      axis2 = attrData[p.first].ev1();
    } else {
      axis2 = attrData[p.first].ev2();
    }

    double angle = 0;

    if(norm(axis1) > 0 and norm(axis2) > 0){
      double cosAngle = axis1 * axis2 / norm(axis1) / norm(axis2);
      if(cosAngle < 0) cosAngle *= -1;
      angle = acos(cosAngle);

    }
    mesh->labelHeat()[p.first] = angle/M_PI*180.;
    vecValues.push_back(angle);

  }

  mesh->heatMapBounds() = Point2f(0, 90);
  mesh->setShowLabel("Label Heat");
  mesh->updateTriangles();

  return true;
  }
  REGISTER_PROCESS(CustomDirectionsAngle);


 bool CustomDirectionsAngle3D::run(Mesh *mesh, QString dirAxis, QString dirCustom)
  {
  const QString attrMapName = "Measure Label Tensor CustomDirections";
  AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  if(attrData.empty())
    return setErrorMessage("No Custom Directions in the current mesh.");

  // if there is no cell axis stored, display error message
  const IntSymTensorAttr& cellAxis = mesh->cellAxis();
  if(cellAxis.size() == 0){
    setErrorMessage(QString("No cell axis stored in active mesh!"));
    return false;
  }

  // Update the normals and center of the cell for visualisation.
  mesh->updateCentersNormals();

  // cell centers
  IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
  IntPoint3fAttr normals = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());

  std::vector<double> vecValues;


  mesh->labelHeat().clear();
  forall(const IntSymTensorPair& p, cellAxis){
    Point3f axis1, axis2;
    if(dirAxis == "Max"){
      axis1 = p.second.ev1();
    } else if(dirAxis == "Mid"){
      axis1 = p.second.ev2();
    } else {
      axis1 = p.second.evals();
    }
    if(dirCustom == "X"){
      axis2 = attrData[p.first].ev1();
    } else if(dirCustom == "Y"){
      axis2 = attrData[p.first].ev2();
    } else {
      axis2 = attrData[p.first].evals();
    }

    double angle = 0;

    if(norm(axis1) > 0 and norm(axis2) > 0){
      double cosAngle = axis1 * axis2 / norm(axis1) / norm(axis2);
      if(cosAngle < 0) cosAngle *= -1;
      angle = acos(cosAngle);

    }
    mesh->labelHeat()[p.first] = angle/M_PI*180;
    vecValues.push_back(angle);

  }

  mesh->heatMapBounds() = Point2f(0, 90);
  mesh->setShowLabel("Label Heat");
  mesh->updateTriangles();

  return true;
  }
  REGISTER_PROCESS(CustomDirectionsAngle3D);



 // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
 bool CellAxisAttrMapExport::run(Mesh *m, Mesh *m2, QString prefix, QString name, QString key, QString value, QString meshNumber, bool cleanMesh)
  {

  // const QString attrMapName = prefix + " " + name;

  // AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>(attrMapName);

  // forall(const IntSymTensorPair& p, mesh->cellAxis()){
  //   attrData[p.first] = p.second;
  // }


  Mesh* useMesh = m;
  if(meshNumber!="Active Mesh"){
    useMesh = m2;
  }

    AttrMap<int, SymmetricTensor>& attrData = useMesh->attributes().attrMap<int, SymmetricTensor>(prefix + " " + name);
    attrData.clear();

    std::vector<int> labels = findAllLabels(m->graph());

    if(key == "Parent" and value == "Parent Axis"){
      forall(IntSymTensorPair p, m->cellAxis())
          attrData[p.first] = p.second;
    } else if(key == "Parent" and value == "Label Axis"){
      forall(IntSymTensorPair p, m->cellAxis())
          attrData[m->parents()[p.first]] = p.second;
    } else if(key == "Label" and value == "Parent Axis"){
      // go through all labels, take parent value
      forall(const int& l, labels){
        if(l < 1) continue;
        attrData[l] = m->cellAxis()[m->parents()[l]];
      }
    } else { // Label, Label Heat
      std::cout<<"Saving the heat-map"<<std::endl;
      forall(IntSymTensorPair p, m->cellAxis()){
        bool add_value = !cleanMesh;

        if(cleanMesh){
          forall(const int& l, labels)
            if(l == p.first)
              add_value = true;
        }

        if(add_value){
          attrData[p.first] = p.second;
        }

      }
    }


  return true;
  }
  REGISTER_PROCESS(CellAxisAttrMapExport);

  // Display custom directions from the attribute map
  bool DisplayCustomOrientations::run(Mesh *mesh, QString directions, const QColor& qaxisColor,
                                             float axisLineWidth, float axisLineScale, float axisOffset,
                                             float orientationThreshold)
  {
    AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    // if there are no custom directions stored, display error message
    if(attrData.empty()){
      setErrorMessage(QString("No custom directions stored in active mesh!"));
      return false;
    }

    // Check cell axis and prepare them for display. Populate cellAxisScaled in Mesh.
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();
    mesh->setAxisWidth(axisLineWidth);
    mesh->setAxisOffset(axisOffset);
    Colorb axisColor(qaxisColor);

    // Update the normals and center of the cell for visualisation.
    if(mesh->labelCenterVis().empty() or mesh->labelNormalVis().empty())
      mesh->updateCentersNormals();

    IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis() : mesh->labelCenterVis());

    // Check if there is a cell center for each cell axis
    IntPoint3fAttr labelCenterVis = mesh->labelCenterVis();
    int nProblems = 0;
    forall(const IntSymTensorPair &p, attrData) {
      int label = p.first;
      if(labelCenterVis.count(label) == 0)
        nProblems++;
    }
    if(nProblems != 0)
      SETSTATUS("Warning: non-existing cell center found for " << nProblems << " cell axis.");

    // Scale the cell axis for display
    //forall(const IntSymTensorPair &p, cellAxis) {
    forall(const IntSymTensorPair &p, attrData) {
      int cell = p.first;
      const SymmetricTensor& tensor = p.second;
      cellAxisVis[cell][0] = Point3f(0,0,0);
      cellAxisVis[cell][1] = Point3f(0,0,0);
      cellAxisVis[cell][2] = Point3f(0,0,0);
      if(directions == "X" or directions == "XY" or directions == "XYZ") cellAxisVis[cell][0] = tensor.ev1() * axisLineScale;
      if(directions == "Y" or directions == "XY" or directions == "XYZ") cellAxisVis[cell][1] = tensor.ev2() * axisLineScale;
      if(directions == "Z" or directions == "XYZ") cellAxisVis[cell][2] = tensor.evals() * axisLineScale;

      cellAxisColor[cell][0] = axisColor;
      cellAxisColor[cell][1] = axisColor;
      cellAxisColor[cell][2] = axisColor;
    }

    return true;
  }
  REGISTER_PROCESS(DisplayCustomOrientations);


  // // Display custom directions from the attribute map
  // bool DisplayCustomOrientations::run(Mesh *mesh, QString directions, const QColor& qaxisColor,
  //                                            float axisLineWidth, float axisLineScale, float axisOffset)
  // {
    
  //   IntSymTensorAttr &cellAxis = mesh->cellAxis();
  //   AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

  //   cellAxis = attrData;

  //   QString mode = "";
    
  //   if(directions == "X" or directions == "Y" or directions == "Z" or directions == "XY" or directions == "XYZ"){
  //     mode = directions;
  //   } else {
  //     return setErrorMessage("Invalid Directions");
  //   }
  //   bool custom = false;
  //   QColor qaxisColorX = qaxisColor;
  //   QColor qaxisColorY = qaxisColor;
  //   QColor qaxisColorZ = qaxisColor;

  //   double orientationThreshold = -1;

  //   DisplayCellAxisTemplate displayAxisP(*this);
  //   displayAxisP.run(mesh, cellAxis, mode, custom, qaxisColorX, qaxisColorY, qaxisColorZ, axisLineWidth, axisLineScale, axisOffset,
  //     orientationThreshold);


  //   return true;
  // }
  // REGISTER_PROCESS(DisplayCustomOrientations);



  // Display custom directions from the attribute map
  bool DisplayCellAxisFromTensor::run(Mesh *mesh, IntSymTensorAttr &cellAxis,
    QString mode, bool custom, const QColor& qaxisColorX, const QColor& qaxisColorY, const QColor& qaxisColorZ,
    const QColor& qaxisColorShrink,
    double axisLineWidth, double axisLineScale, double axisOffset, double orientationThreshold, bool strain)
  {

  /*
   * modes
   *******
   * X - show only X
   * Y - show only Y
   * Z - show only Z
   * XY - show X & Y
   * XZ - show X & Z
   * YZ - show Y & Z
   * XYZ - show all
   * Difference - X-Y
   * Ratio - X/Y along X dir
   ********
   * if custom then use projection on custom direction instead
   */


    AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    // if there are no custom directions stored, display error message
    if(attrData.empty() and custom){
      std::cout << "No custom directions stored in active mesh!" << std::endl;
      return false;
    }

    // Check cell axis and prepare them for display. Populate cellAxisScaled in Mesh.
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();
    mesh->setAxisWidth(axisLineWidth);
    mesh->setAxisOffset(axisOffset);
    Colorb axisColorX(qaxisColorX);
    Colorb axisColorY(qaxisColorY);
    Colorb axisColorZ(qaxisColorZ);

    // Update the normals and center of the cell for visualisation.
    if(mesh->labelCenterVis().empty() or mesh->labelNormalVis().empty())
      mesh->updateCentersNormals();

    IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis() : mesh->labelCenterVis());

    // Check if there is a cell center for each cell axis
    IntPoint3fAttr labelCenterVis = mesh->labelCenterVis();
    int nProblems = 0;
    forall(const IntSymTensorPair &p, attrData) {
      int label = p.first;
      if(labelCenterVis.count(label) == 0)
        nProblems++;
    }
    if(nProblems != 0)
      SETSTATUS("Warning: non-existing cell center found for " << nProblems << " cell axis.");

    // Scale the cell axis for display
    //forall(const IntSymTensorPair &p, cellAxis) {
    forall(const IntSymTensorPair &p, cellAxis) {
      int cell = p.first;
      SymmetricTensor tensor = p.second;
      double length1st, length2nd, length3rd, anisotropy;
      Point3f dir1st, dir2nd,dir3rd;

      
      // if()
      // if(cellAxisCustom.find(cell) == cellAxisCustom.end())
      //   continue;

      length1st = tensor.evals()[0];
      length2nd = tensor.evals()[1];
      length3rd = tensor.evals()[2];
      if(strain){
        length1st -= 1;
        length2nd -= 1;
        length3rd -= 1;
      }
      dir1st = tensor.ev1();
      dir2nd = tensor.ev2();
      dir3rd = dir1st % dir2nd;
      anisotropy = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);

      Point3b showAxis;
      showAxis[0] = (mode == "X" or mode == "XY" or mode == "XZ" or mode == "XYZ");
      showAxis[1] = (mode == "Y" or mode == "XY" or mode == "YZ" or mode == "XYZ");
      showAxis[2] = (mode == "Z" or mode == "XZ" or mode == "YZ" or mode == "XYZ");

      if(orientationThreshold > 0 and anisotropy < orientationThreshold)
        showAxis = Point3b(false, false, false);

      cellAxisVis[cell][0] = (showAxis[0] ? length1st : 0.f) * dir1st * axisLineScale;
      cellAxisVis[cell][1] = (showAxis[0] ? length2nd : 0.f) * dir2nd * axisLineScale;
      cellAxisVis[cell][2] = (showAxis[0] ? length3rd : 0.f) * dir1st * axisLineScale;

      cellAxisColor[cell][0] = axisColorX;
      cellAxisColor[cell][1] = axisColorY;
      cellAxisColor[cell][2] = axisColorZ;

      for(size_t i = 0; i < 3; ++i) {
        if(tensor.evals()[i] - 1 < 0)
          cellAxisColor[cell][i] = Colorb(qaxisColorShrink);
      }


    }

    mesh->updateTriangles();

    return true;
  }
  //REGISTER_PROCESS(DisplayCellAxisTemplate);


  // Display custom directions from the attribute map
  bool DisplayCustomOrientations3D::run(Mesh *mesh, QString directions, const QColor& qaxisColor,
                                             float axisLineWidth, float axisLineScale,
                                             float orientationThreshold)
  {
    AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    // if there are no custom directions stored, display error message
    if(attrData.empty()){
      setErrorMessage(QString("No custom directions stored in active mesh!"));
      return false;
    }

    AttrMap<int, Point3d>& centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    if(mesh->useParents()){
      centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
      std::cout << "using the parent labels " << centroids.size() << std::endl;
    }

    if(centroids.empty()){
      return setErrorMessage("No valid cell centroids found!");
    }

    // Check cell axis and prepare them for display. Populate cellAxisScaled in Mesh.
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();
    mesh->setAxisWidth(axisLineWidth);
    mesh->setAxisOffset(0.0);
    Colorb axisColor(qaxisColor);

    // Check if there is a cell center for each cell axis
    IntPoint3fAttr& labelCenterVis = mesh->labelCenterVis();
    IntPoint3fAttr& parentCenterVis = mesh->parentCenterVis();
    
    labelCenterVis.clear();
    parentCenterVis.clear();
    int nProblems = 0;
    forall(const IntSymTensorPair &p, attrData) {
      //if(mesh->useParents()){
        parentCenterVis[p.first] = Point3f(centroids[p.first]);
      //} else {
        labelCenterVis[p.first] = Point3f(centroids[p.first]);
      //}
      
      int label = p.first;
      if(centroids.count(label) == 0)
        nProblems++;
    }
    if(nProblems != 0)
      SETSTATUS("Warning: non-existing cell centroid found for " << nProblems << " cells.");

    // Scale the cell axis for display
    forall(const IntSymTensorPair &p, attrData) {
      int cell = p.first;
      const SymmetricTensor& tensor = p.second;
      cellAxisVis[cell][0] = Point3f(0,0,0);
      cellAxisVis[cell][1] = Point3f(0,0,0);
      cellAxisVis[cell][2] = Point3f(0,0,0);
      if(directions == "X" or directions == "XY" or directions == "XZ" or directions == "XYZ") cellAxisVis[cell][0] = tensor.ev1() * axisLineScale;
      if(directions == "Y" or directions == "XY" or directions == "YZ" or directions == "XYZ") cellAxisVis[cell][1] = tensor.ev2() * axisLineScale;
      if(directions == "Z" or directions == "XZ" or directions == "YZ" or directions == "XYZ") cellAxisVis[cell][2] = tensor.evals() * axisLineScale;

      cellAxisColor[cell][0] = axisColor;
      cellAxisColor[cell][1] = axisColor;
      cellAxisColor[cell][2] = axisColor;
    }

    mesh->updateTriangles();
    mesh->setShowAxis("Cell Axis");

    return true;
  }
  REGISTER_PROCESS(DisplayCustomOrientations3D);


  // Project the cell axis on local coordinates defined by Bezier. NB: we update the cell axis themselves, not cellAxisVis!
  bool CellAxisAttrMapImport::run(Mesh *mesh, QString type, QString prefix, QString attrName, bool reproject)
  {
    AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>(prefix + " " + attrName);

    if(attrData.empty()) 
      return setErrorMessage("No Attribute Map with this name!");

    // project onto new cell planes
    if(reproject) {
      mesh->updateCentersNormals();

      // cell centers
      IntPoint3fAttr centers = (mesh->useParents()? mesh->parentCenterVis():mesh->labelCenterVis());
      IntPoint3fAttr normals = (mesh->useParents()? mesh->parentNormal():mesh->labelNormal());

      // for each cell: check heat difference to all of its neighbors.
      forall(auto p, centers){
        Point3d dir(0,0,0);

        int cell = p.first; // cell label
        Point3d cellCenter(centers[cell]);
        Point3d cellNormal(normals[cell]);

        Point3d maxDir = Point3d(attrData[p.first].ev1());
        Point3d minDir = Point3d(attrData[p.first].ev2());

        Point3d projPMax = projectPointOnPlane(cellCenter+maxDir,cellCenter,cellNormal);
        Point3d projPMin = projectPointOnPlane(cellCenter+minDir,cellCenter,cellNormal);

        attrData[p.first].ev1() = Point3f(Point3f(projPMax - cellCenter));
        attrData[p.first].ev2() = Point3f(Point3f(projPMin - cellCenter));
      }
    }

    forall(const IntSymTensorPair& p, attrData)
      mesh->cellAxis()[p.first] = p.second;

    // Display cell axis according to their type, with the parameters from the GUI
    if(type == "curvature") {
      DisplayTissueCurvature *proc;
      if(!getProcess("Mesh/Cell Axis/Curvature/Display Tissue Curvature", proc))
        throw QString("%1:: Unable to make DisplayTissueCurvature process").arg(name());
      return proc->run();
    } else if(type == "PDG") {
      mesh->setCellAxisType("PDG");
    DisplayPDGs *proc;
    if(!getProcess("Mesh/Cell Axis/PDG/Display Growth Directions", proc))
      throw QString("%1:: Unable to make DisplayPDGs process").arg(name());
			// DisplayPDGs::initialize() changes the display parms in case they commend to draw the PDGs projected
			// on a Bezier and cellAxisBezier doesn't exist yet.
  		// ALR: I don't know how to update the GUI with the new parms.
      proc->initialize(0);
      return proc->run();
    } else if(type == "fibril") {
      mesh->setCellAxisType("fibril");
	    DisplayFibrilOrientations *proc;
      if(!getProcess("Mesh/Cell Axis/Fibril Orientations/Display Fibril Orientations", proc))
        throw(QString("%1:: Unable to make DisplayFibrilOrientations process"));
      return proc->run();
    } else if(type == "polarization") {
	    DisplaySignalOrientation *proc;
      if(!getProcess("Mesh/Cell Axis/Polarization/Display Signal Orientation", proc))
        throw QString("%1:: Unable to make DisplaySignalOrientation process").arg(name());
      return proc->run();
    } else if(type == "PCA") {
      mesh->setCellAxisType("PCA");
	    DisplayPCAOld *proc;
      if(!getProcess("Mesh/Cell Axis/Shape Analysis/Display Shape Axis", proc))
        throw QString("%1:: Unable to make DisplayPCA process").arg(name());
      return proc->run();
		}
    return true;
  }
  REGISTER_PROCESS(CellAxisAttrMapImport);

    // bool CellAxisPDGAngle::initialize(QWidget *parent){
    // QString workingDir = QDir::currentPath();
    //   if(parm("T1 mesh").isEmpty()){
    //     setParm("T1 mesh", QFileDialog::getOpenFileName(parent, "Choose mesh for T1", QDir::currentPath(), "Mesh files (*.mgxm);;All files (*.*)"));
    //     QFileInfo fi(parm("T1 mesh"));
    //     workingDir = fi.absoluteFilePath();
    //   }
    //   if(parm("T2 mesh").isEmpty()){
    //     setParm("T2 mesh", QFileDialog::getOpenFileName(parent, "Choose mesh for T2", workingDir, "Mesh files (*.mgxm);;All files (*.*)"));
    //     QFileInfo fi(parm("T2 mesh"));
    //     workingDir = fi.absoluteFilePath();
    //   }
    //   if(parm("T3 mesh").isEmpty()){
    //     setParm("T3 mesh", QFileDialog::getOpenFileName(parent, "Choose mesh for T3", workingDir, "Mesh files (*.mgxm);;All files (*.*)"));
    //     QFileInfo fi(parm("T3 mesh"));
    //     workingDir = fi.absoluteFilePath();
    //   }
    //   if(parm("T1-T2 parents").isEmpty()){
    //    setParm( "T1-T2 parents", QFileDialog::getOpenFileName(parent, "Choose parents T1-T2", workingDir, "CSV files (*.csv)"));
    //     QFileInfo fi(parm("T1-T2 parents"));
    //     workingDir = fi.absoluteFilePath();
    //   }
    //   if(parm("T1-T3 parents").isEmpty()){
    //     setParm("T1-T3 parents", QFileDialog::getOpenFileName(parent, "Choose parents T1-T3", workingDir, "CSV files (*.csv)"));
    //     QFileInfo fi(parm("T1-T3 parents"));
    //     workingDir = fi.absoluteFilePath();
    //   }

    //   if(parm("T1 mesh").isEmpty() or parm("T2 mesh").isEmpty() or parm("T3 mesh").isEmpty() or parm("T1-T2 parents").isEmpty() or parm("T1-T3 parents").isEmpty())
    //       return false;


    //   return true;
    // }


    // bool CellAxisPDGAngle::run()
    // {
    //   return run(parm("T1 mesh"), parm("T2 mesh"), parm("T3 mesh"),parm("T1-T2 parents"),parm("T1-T3 parents"),parm("Isotropy threshold"));
    // }


// bool CellAxisPDGAngle::run(QString meshNameT1, QString meshNameT2, QString meshNameT3,QString T1T2,QString T1T3,QString isoThresh)
//   {

//      MeshLoad *loadMesh = new MeshLoad(*this);
//      LoadParents *loadParents = new LoadParents(*this);
//      CorrespondenceJunctions *checkCorr = new CorrespondenceJunctions(*this);
//      GrowthDirections *growthDir = new GrowthDirections(*this);
//      CellAxisAttrMapExport *exportPDG = new CellAxisAttrMapExport(*this);
//      CopyParentsToLabels *copyParToLab = new CopyParentsToLabels(*this);
//      CellAxisAngle *cellAxisAngle = new CellAxisAngle(*this);

//       Mesh *m1 = mesh(0);
//       Mesh *m2 = mesh(1);
//       m1->setUseParents(false);
//       m2->setUseParents(false);

//       //Load T0 and T1
//       bool m1Loaded = loadMesh->run(m1,meshNameT1,false,false);
//       bool m2Loaded = loadMesh->run(m2,meshNameT2,false,false);
//       if(!m1Loaded || !m2Loaded)
//         return false;

//       QString fileType = T1T2.split(".").last().toUpper();
//       bool p12Loaded = loadParents->run(m2,T1T2,fileType,false);
//       if(!p12Loaded)
//         return false;
//       m2->setUseParents(true);
//       bool checkCor12 = checkCorr->run(m1,m2,false,false);
//       if(!checkCor12)
//         return false;
//       bool growthDir12 = growthDir->run(m1,m2);
//       if(!growthDir12)
//         return false;

//       QString prefix = "Measure Label Tensor";
//       QString attMapName1 = "PDG T1T2 on T2";
//       bool exportedPDG12 = exportPDG->run(m2,prefix,attMapName1);
//       if(!exportedPDG12)
//         return false;
//       attMapName1 = prefix + " " + attMapName1;

//       bool copiedParToLabM2 = copyParToLab->run(m2);
//       m2->setUseParents(false);
//       //Now load T2
//       Mesh *m3 = mesh(0);
//       bool m3Loaded = loadMesh->run(m3,meshNameT3,false,false);
//       if(!m3Loaded)
//         return false;

//       fileType = T1T2.split(".").last().toUpper();
//       bool p13Loaded = loadParents->run(m3,T1T3,fileType,false);
//       if(!p13Loaded)
//         return false;
//       m3->setUseParents(true);
//       bool checkCor13 = checkCorr->run(m2,m3,false,false);
//       if(!checkCor13)
//         return false;
//       bool growthDir23 = growthDir->run(m2,m3);
//       if(!growthDir23)
//         return false;

//       QString attMapName2 = "PDG T2T3 on T2";
//       bool exportedPDG23 = exportPDG->run(m2,prefix,attMapName2);
//       if(!exportedPDG23)
//         return false;
//       attMapName2 = prefix + " " + attMapName2;
//       QString attAngleOutput = "Measure Label Double AxisAngleT2";
//       bool cellAxisAngleWorked = cellAxisAngle->run(m2,attMapName1,attMapName2,attAngleOutput,isoThresh);
//       if(!cellAxisAngleWorked)
//         return false;

//   return true;

//   }
//   REGISTER_PROCESS(CellAxisPDGAngle);



  inline float angle(const Point3f &v1, const Point3f &v2)
  {
    float pi = 4*atan(1);
    float q = (v1*v2) / (sqrt(v1*v1) * sqrt(v2*v2));
    float a = min(acos(q), acos(-q));
    return (180/pi) * a;
  }

  bool CellAxisAngles::run(Mesh* m, bool usePrefix, QString attr1, QString dir1, double isoT1, QString attr2, QString dir2, 
      double isoT2, 
      bool displayAxis,const QColor& color1,const QColor& color2,
      double lineWidth, double lineScale, double lineOffset)
  {
    QString prefix = "";

    if(usePrefix) prefix = "Measure Label Tensor ";

    AttrMap<int, SymmetricTensor>& attrMap1 = m->attributes().attrMap<int, SymmetricTensor>(prefix + attr1);
    AttrMap<int, SymmetricTensor>& attrMap2 = m->attributes().attrMap<int, SymmetricTensor>(prefix + attr2);

    if(attrMap1.empty()) return setErrorMessage("Attr Map 1 empty!");
    if(attrMap2.empty()) return setErrorMessage("Attr Map 2 empty!");
    
    IntVec3ColorbAttr& cellAxisColor = m->cellAxisColor();
    cellAxisColor.clear();

    m->setAxisOffset(lineOffset);
    m->setAxisWidth(lineWidth);
    Colorb c1(color1);
    Colorb c2(color2);

    m->labelHeat().clear();

    IntMatrix3fAttr &cellAxisVis = m->cellAxisVis();
    cellAxisVis.clear();

    forall(auto p, attrMap1){
      int label = p.first;
      if(attrMap2.find(label) == attrMap2.end()) continue;
      double aniso1 = p.second.evals()[1] != 0 ? p.second.evals()[0]/p.second.evals()[1] : 0;
      double aniso2 = attrMap2[p.first].evals()[1] != 0 ? attrMap2[p.first].evals()[0]/attrMap2[p.first].evals()[1] : 0;
      if(aniso1 < isoT1 and isoT1 > 0) continue;
      if(aniso2 < isoT2 and isoT1 > 0) continue;

      Point3f vec1 = p.second.ev1();
      if(dir1 == "Min/Y") vec1 = p.second.ev2();

      Point3f vec2 = attrMap2[label].ev1();
      if(dir2 == "Min/Y") vec2 = attrMap2[label].ev2();

      m->labelHeat()[label] = angle(vec1, vec2);

      if(displayAxis){
        cellAxisVis[label][0] = vec1 * lineScale;
        cellAxisVis[label][1] = vec2 * lineScale;
        cellAxisVis[label][2] = Point3f(0,0,0);
        cellAxisColor[label][0] = c1;
        cellAxisColor[label][1] = c2;
      }

    }


    m->setShowLabel("Label Heat");
    m->heatMapUnit() = "°";
    m->heatMapBounds() = m->calcHeatMapBounds();
    m->updateTriangles();
   
    return true;
  }
  REGISTER_PROCESS(CellAxisAngles);


  bool CellAxisAngles3D::run(Mesh* m, bool usePrefix, QString attr1, QString dir1, double isoT1, 
    QString attr2, QString dir2, double isoT2, bool displayAxis)
  {
    QString prefix = "";

    if(usePrefix) prefix = "Measure Label Tensor ";

    AttrMap<int, SymmetricTensor>& attrMap1 = m->attributes().attrMap<int, SymmetricTensor>(prefix + attr1);
    AttrMap<int, SymmetricTensor>& attrMap2 = m->attributes().attrMap<int, SymmetricTensor>(prefix + attr2);

    if(attrMap1.empty()) return setErrorMessage("Attr Map 1 empty!");
    if(attrMap2.empty()) return setErrorMessage("Attr Map 2 empty!");
    
    m->labelHeat().clear();

    IntMatrix3fAttr &cellAxisVis = m->cellAxisVis();
    cellAxisVis.clear();

    forall(auto p, attrMap1){
      int label = p.first;
      if(attrMap2.find(label) == attrMap2.end()) continue;
      double aniso1 = p.second.evals()[1] != 0 ? p.second.evals()[0]/p.second.evals()[1] : 0;
      double aniso2 = attrMap2[p.first].evals()[1] != 0 ? attrMap2[p.first].evals()[0]/attrMap2[p.first].evals()[1] : 0;
      if(aniso1 < isoT1 and attr1 != "CustomDirections") continue;
      if(aniso2 < isoT2 and attr2 != "CustomDirections") continue;

      Point3f vec1 = p.second.ev1();
      if(dir1 == "Max") vec1 = p.second.ev1();
      else if(dir1 == "Mid") vec1 = p.second.ev2();
      else if(attr1 == "CustomDirections") vec1 = p.second.evals();
      else vec1 = p.second.ev1() ^ p.second.ev2();

      Point3f vec2 = attrMap2[label].ev1();
      if(dir2 == "Max") vec2 = attrMap2[label].ev1();
      else if(dir2 == "Mid") vec2 = attrMap2[label].ev2();
      else if(attr2 == "CustomDirections") vec2 = attrMap2[label].evals();
      else vec2 = attrMap2[label].ev1() ^ attrMap2[label].ev2();

      m->labelHeat()[label] = angle(vec1, vec2);

      if(displayAxis){
        cellAxisVis[label][0] = vec1;
        cellAxisVis[label][1] = vec2;
        cellAxisVis[label][2] = Point3f(0,0,0);
      }

    }


    m->setShowLabel("Label Heat");
    m->heatMapUnit() = "°";
    m->heatMapBounds() = m->calcHeatMapBounds();
    m->updateTriangles();
   
    return true;
  }
  REGISTER_PROCESS(CellAxisAngles3D);



}
