//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//
#include "PCAnalysis.hpp"
#include "Information.hpp"
#include "Progress.hpp"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include "UnorderedMap.hpp"
#include "UnorderedSet.hpp"
#include <algorithm>
#include <GraphUtils.hpp>
#include <MeshProcessPDG.hpp> // custom directions
#include <MeshProcessMeasures.hpp>
#include <MeshProcessCellAxis.hpp> // display cell axis processes

#include <QFile>
#include <QTextStream>

namespace mgx
{
  // Note: empty namespace = content is not exported, and it's better (cleaner and more general) than the static keyword
  namespace
	{
    struct PCAnalysis_result
  	{
      Point3f p1, p2, p3;
      Point3f ev;
      Point3f mean;

      bool valid() const { return normsq(ev) > 0; }

      operator bool() const { return valid(); }
    };

    struct SelectThreshold
  	{
      SelectThreshold(ushort th) : threshold(th) {}

      int operator()(ushort v) const
      {
        if(v > threshold)
          return 0;
        return -1;
      }

      ushort value(int vid) const
      {
        if(vid == 1)
          return threshold;
        return 0;
      }

      int nb_values() const {
        return 1;
      }

      ushort threshold;
    };

    struct SelectLabel
   	{
      SelectLabel(const std::vector<ushort>& ls) : labels(ls)
      {
        for(size_t i = 0; i < labels.size(); ++i)
          inv_labels[labels[i]] = i;
      }

      int operator()(ushort v) const
      {
        std::unordered_map<ushort, int>::const_iterator found = inv_labels.find(v);
        if(found != inv_labels.end())
          return found->second;
        return -1;
      }

      ushort value(int vid) const
      {
        if(vid >= 0 and vid < (int)labels.size())
          return labels[vid];
        return 0;
      }

      int nb_values() const {
        return inv_labels.size();
      }

      std::vector<ushort> labels;
      std::unordered_map<ushort, int> inv_labels;
    };

    template <typename Fct>
    std::vector<PCAnalysis_result> analyzePC(const Stack* stk, const HVecUS& data, const Fct& selection,
                                             Point3f correcting_factor)
    {
      // First, compute mean
      size_t nb_values = selection.nb_values();
      std::vector<Point3f> mean(nb_values, Point3f(0, 0, 0));
      std::vector<unsigned long long> sum(nb_values, 0);
      std::vector<PCAnalysis_result> result(nb_values);

      // Compute the mean for each selected value
      Information::out << "Compute mean" << endl;
      Point3u s = stk->size();
      size_t k = 0;
      for(size_t z = 0; z < s.z(); ++z)
        for(size_t y = 0; y < s.y(); ++y)
          for(size_t x = 0; x < s.x(); ++x, ++k) {
            int vid = selection(data[k]);
            if(vid >= 0) {
              ++sum[vid];
              mean[vid] += stk->imageToWorld(Point3d(x, y, z));
            }
          }
      for(size_t vid = 0; vid < nb_values; ++vid) {
        if(sum[vid] > 0)
          mean[vid] /= double(sum[vid]);
      }
      Information::out << "Compute CC matrix" << endl;
      // Compute the cross-correlation matrix
      std::vector<Matrix3f> corr(nb_values);
      k = 0;
      for(size_t z = 0; z < s.z(); ++z)
        for(size_t y = 0; y < s.y(); ++y)
          for(size_t x = 0; x < s.x(); ++x, ++k) {
            int vid = selection(data[k]);
            if(vid >= 0) {
              Point3f dp = stk->imageToWorld(Point3f(x, y, z)) - mean[vid];
              corr[vid](0, 0) += dp[0] * dp[0];
              corr[vid](1, 1) += dp[1] * dp[1];
              corr[vid](2, 2) += dp[2] * dp[2];
              corr[vid](1, 0) += dp[0] * dp[1];
              corr[vid](2, 0) += dp[0] * dp[2];
              corr[vid](2, 1) += dp[1] * dp[2];
            }
          }
      gsl_matrix* mat = gsl_matrix_alloc(3, 3);
      gsl_vector* eval = gsl_vector_alloc(3);
      gsl_matrix* evec = gsl_matrix_alloc(3, 3);
      gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(3);
      Information::out << "Decompose matrices" << endl;
      // Eigen-decomposition of the matrices
      for(size_t vid = 0; vid < nb_values; ++vid) {
        if(sum[vid] == 0)
          continue;

        corr[vid](0, 1) = corr[vid](1, 0);
        corr[vid](0, 2) = corr[vid](2, 0);
        corr[vid](1, 2) = corr[vid](2, 1);
        corr[vid] /= sum[vid];

        for(int i = 0; i < 3; ++i)
          for(int j = 0; j < 3; ++j)
            gsl_matrix_set(mat, i, j, corr[vid](i, j));
        gsl_eigen_symmv(mat, eval, evec, w);
        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

        Point3f p1(gsl_matrix_get(evec, 0, 0), gsl_matrix_get(evec, 1, 0), gsl_matrix_get(evec, 2, 0));
        Point3f p2(gsl_matrix_get(evec, 0, 1), gsl_matrix_get(evec, 1, 1), gsl_matrix_get(evec, 2, 1));
        Point3f p3(gsl_matrix_get(evec, 0, 2), gsl_matrix_get(evec, 1, 2), gsl_matrix_get(evec, 2, 2));
        Point3f ev(gsl_vector_get(eval, 0), gsl_vector_get(eval, 1), gsl_vector_get(eval, 2));
        if((p1 ^ p2) * p3 < 0)
          p3 = -p3;

        // Correct eigen-vector for shape-factor
        if(correcting_factor.x() > 0)
          ev = multiply(Point3f(correcting_factor), map(sqrtf, ev));

        // Information::out << "Size of the various dimension = " << ev << endl;

        result[vid].p1 = p1;
        result[vid].p2 = p2;
        result[vid].p3 = p3;
        result[vid].ev = ev;
        result[vid].mean = mean[vid];
        // Information::out << "Value: " << selection.value(vid) << " -- Major axis = " << p1 << " eigenvalues = "
        // << ev << endl;
      }

      if(correcting_factor.x() == 0) {
        Information::out << "Compute regions spans" << endl;
        std::vector<Point3f> min_pos(nb_values, Point3f(FLT_MAX));
        std::vector<Point3f> max_pos(nb_values, Point3f(-FLT_MAX));
        // Find the maximum span of each label along each direction
        k = 0;
        for(size_t z = 0; z < s.z(); ++z)
          for(size_t y = 0; y < s.y(); ++y)
            for(size_t x = 0; x < s.x(); ++x, ++k) {
              int vid = selection(data[k]);
              if(vid >= 0) {
                Point3f dp = stk->imageToWorld(Point3d(x, y, z));
                const PCAnalysis_result& res = result[vid];
                double p1 = dp * res.p1;
                double p2 = dp * res.p2;
                double p3 = dp * res.p3;
                Point3f& pmin = min_pos[vid];
                Point3f& pmax = max_pos[vid];
                if(p1 < pmin.x())
                  pmin.x() = p1;
                if(p1 > pmax.x())
                  pmax.x() = p1;
                if(p2 < pmin.y())
                  pmin.y() = p2;
                if(p2 > pmax.y())
                  pmax.y() = p2;
                if(p3 < pmin.z())
                  pmin.z() = p3;
                if(p3 > pmax.z())
                  pmax.z() = p3;
              }
            }
        for(size_t vid = 0; vid < nb_values; ++vid)
          if(sum[vid] > 0) {
            PCAnalysis_result& res = result[vid];
            Point3f pmin = min_pos[vid];
            Point3f pmax = max_pos[vid];
            res.ev = (pmax - pmin) / 2;
            res.mean = (pmax + pmin) / 2;
            res.mean = res.mean.x() * res.p1 + res.mean.y() * res.p2 + res.mean.z() * res.p3;
          }
      } else if(correcting_factor.x() < 0) {
        float percentile = -correcting_factor.x();
        Information::out << QString("Computing the %1% of voxels").arg(percentile * 100) << endl;
        std::vector<std::vector<double> > pos_p1(nb_values);
        std::vector<std::vector<double> > pos_p2(nb_values);
        std::vector<std::vector<double> > pos_p3(nb_values);
        k = 0;
        for(size_t z = 0; z < s.z(); ++z)
          for(size_t y = 0; y < s.y(); ++y)
            for(size_t x = 0; x < s.x(); ++x, ++k) {
              int vid = selection(data[k]);
              if(vid >= 0) {
                Point3f dp = stk->imageToWorld(Point3d(x, y, z));
                const PCAnalysis_result& res = result[vid];
                double p1 = dp * res.p1;
                double p2 = dp * res.p2;
                double p3 = dp * res.p3;
                pos_p1[vid].push_back(p1);
                pos_p2[vid].push_back(p2);
                pos_p3[vid].push_back(p3);
              }
            }
        for(size_t vid = 0; vid < nb_values; ++vid)
          if(sum[vid] > 0) {
            std::vector<double>& ps1 = pos_p1[vid];
            std::vector<double>& ps2 = pos_p2[vid];
            std::vector<double>& ps3 = pos_p3[vid];

            std::sort(ps1.begin(), ps1.end());
            std::sort(ps2.begin(), ps2.end());
            std::sort(ps3.begin(), ps3.end());

            int lower = (1 - percentile) / 2 * (ps1.size() - 1);
            int upper = ps1.size() - 1 - lower;
            // Information::out << "Positions for vid = " << vid << "(" << ps1.size() << ") = [" << lower << ";"
            // << upper << "]\n";
            PCAnalysis_result& res = result[vid];
            Point3f pmin(ps1[lower], ps2[lower], ps3[lower]);
            Point3f pmax(ps1[upper], ps2[upper], ps3[upper]);
            // Information::out << " -> " << pmin << " -- " << pmax << endl;
            res.ev = (pmax - pmin) / 2;
            res.mean = (pmax + pmin) / 2;
            res.mean = res.mean.x() * res.p1 + res.mean.y() * res.p2 + res.mean.z() * res.p3;
          }
      }
      gsl_eigen_symmv_free(w);
      gsl_vector_free(eval);
      gsl_matrix_free(evec);
      gsl_matrix_free(mat);
      return result;
    }

    void addPCACuboidShape(int lab, const PCAnalysis_result& res, Mesh* mesh, int slices)
    {
      if(slices < 1)
        return;
      // Information::out << "Draw PCA for label " << lab << endl;
      vvGraph& S = mesh->graph();
      std::vector<vertex> vs(8, vertex(0));
      Point3f center = res.mean;
      for(int i = 0; i < 8; ++i) {
        vertex v;
        v->label = lab;
        vs[i] = v;
      }
      Point3f p1 = res.ev[0] * res.p1;
      Point3f p2 = res.ev[1] * res.p2;
      Point3f p3 = res.ev[2] * res.p3;

      vs[0]->pos = Point3d((center - p1 - p2 - p3));
      vs[1]->pos = Point3d(center + p1 - p2 - p3);
      vs[2]->pos = Point3d(center - p1 + p2 - p3);
      vs[3]->pos = Point3d(center + p1 + p2 - p3);
      vs[4]->pos = Point3d(center - p1 - p2 + p3);
      vs[5]->pos = Point3d(center + p1 - p2 + p3);
      vs[6]->pos = Point3d(center - p1 + p2 + p3);
      vs[7]->pos = Point3d(center + p1 + p2 + p3);

      std::vector<Point3i> triangles(12);
      triangles[0] = Point3i(0, 1, 4);   // 1
      triangles[1] = Point3i(1, 5, 4);

      triangles[2] = Point3i(1, 3, 5);   // 2
      triangles[3] = Point3i(3, 7, 5);

      triangles[4] = Point3i(3, 2, 7);   // 3
      triangles[5] = Point3i(2, 6, 7);

      triangles[6] = Point3i(2, 0, 6);   // 4
      triangles[7] = Point3i(0, 4, 6);

      triangles[8] = Point3i(2, 3, 0);   // 5
      triangles[9] = Point3i(3, 1, 0);

      triangles[10] = Point3i(4, 5, 6);   // 6
      triangles[11] = Point3i(5, 7, 6);

      // Information::out << "Calling meshFromTriangles" << endl;
      Mesh::meshFromTriangles(S, vs, triangles);
    }

    void addPCACylinderShape(int lab, const PCAnalysis_result& res, Mesh* mesh, int slices)
    {
      if(slices < 1)
        return;
      vvGraph& S = mesh->graph();
      std::vector<vertex> vs(5 * slices + 2, vertex(0));
      Point3f p1 = res.ev[0] * res.p1;
      Point3f p2 = res.ev[1] * res.p2;
      Point3f p3 = res.ev[2] * res.p3;
      Point3f center = res.mean;

      vertex c1, c2;
      c1->label = c2->label = lab;
      c1->pos = Point3d(center - p1);
      c2->pos = Point3d(center + p1);
      int nc1 = 5 * slices;
      int nc2 = 5 * slices + 1;
      vs[nc1] = c1;
      vs[nc2] = c2;

      float alpha = 2 * M_PI / slices;
      for(int i = 0; i < slices; ++i) {
        vertex v1, v2, v3, v4, v5;
        v1->label = v2->label = v3->label = v4->label = v5->label = lab;

        Point3f v1p = Point3f(-1, cos(i * alpha), sin(i * alpha));
        v3->pos = v1->pos = Point3d(center + v1p.x() * p1 + v1p.y() * p2 + v1p.z() * p3);
        Point3f v2p = Point3f(1, cos(i * alpha), sin(i * alpha));
        v4->pos = v2->pos = Point3d(center + v2p.x() * p1 + v2p.y() * p2 + v2p.z() * p3);

        v5->pos = Point3d(center + v2p.y() * p2 + v2p.z() * p3);

        vs[i] = v1;
        vs[slices + i] = v2;
        vs[2 * slices + i] = v3;
        vs[3 * slices + i] = v4;
        vs[4 * slices + i] = v5;
      }

      std::vector<Point3i> triangles(6 * slices);
      int i = slices - 1;
      for(int i1 = 0; i1 < slices; ++i1) {
        int j = slices + i;
        int j1 = slices + i1;
        int k = 2 * slices + i;
        int k1 = 2 * slices + i1;
        triangles[i] = Point3i(nc1, i1, i);
        triangles[i + slices] = Point3i(nc2, j, j1);
        triangles[i + 2 * slices] = 2 * slices + Point3i(i, i1, k);
        triangles[i + 3 * slices] = 2 * slices + Point3i(i1, k1, k);
        triangles[i + 4 * slices] = 2 * slices + Point3i(k, k1, j);
        triangles[i + 5 * slices] = 2 * slices + Point3i(k1, j1, j);
        i = i1;
      }

      Mesh::meshFromTriangles(S, vs, triangles);
    }

    struct pair_hash {
      int operator()(const std::pair<int, int>& p) const {
        return p.first ^ p.second;
      }
    };

    void addPCAEllipsoidShape(int lab, const PCAnalysis_result& res, Mesh* mesh, int subdiv)
    {
      static int generated_subdiv = 0;
      static std::vector<Point3i> triangles;
      static std::vector<Point3f> positions;
      if(generated_subdiv != subdiv) {
        generated_subdiv = subdiv;
        if(subdiv < 0) {
          triangles.clear();
          positions.clear();
          return;
        } else {
          positions.resize(12);
          triangles.resize(20);
          const double t = (1 + sqrt(5.)) / 2;
          const double s = sqrt(1 + t * t);
          positions[0] = Point3f(t, 1, 0) / s;
          positions[1] = Point3f(-t, 1, 0) / s;
          positions[2] = Point3f(t, -1, 0) / s;
          positions[3] = Point3f(-t, -1, 0) / s;
          positions[4] = Point3f(1, 0, t) / s;
          positions[5] = Point3f(1, 0, -t) / s;
          positions[6] = Point3f(-1, 0, t) / s;
          positions[7] = Point3f(-1, 0, -t) / s;
          positions[8] = Point3f(0, t, 1) / s;
          positions[9] = Point3f(0, -t, 1) / s;
          positions[10] = Point3f(0, t, -1) / s;
          positions[11] = Point3f(0, -t, -1) / s;

          triangles[0] = Point3i(0, 8, 4);
          triangles[1] = Point3i(1, 10, 7);
          triangles[2] = Point3i(2, 9, 11);
          triangles[3] = Point3i(7, 3, 1);
          triangles[4] = Point3i(0, 5, 10);
          triangles[5] = Point3i(3, 9, 6);
          triangles[6] = Point3i(3, 11, 9);
          triangles[7] = Point3i(8, 6, 4);
          triangles[8] = Point3i(2, 4, 9);
          triangles[9] = Point3i(3, 7, 11);
          triangles[10] = Point3i(4, 2, 0);
          triangles[11] = Point3i(9, 4, 6);
          triangles[12] = Point3i(2, 11, 5);
          triangles[13] = Point3i(0, 10, 8);
          triangles[14] = Point3i(5, 0, 2);
          triangles[15] = Point3i(10, 5, 7);
          triangles[16] = Point3i(1, 6, 8);
          triangles[17] = Point3i(1, 8, 10);
          triangles[18] = Point3i(6, 1, 3);
          triangles[19] = Point3i(11, 7, 5);

    #define MAKE_EDGE(i1, i2) (i1 < i2 ? std::make_pair(i1, i2) : std::make_pair(i2, i1))

          // Subdivide triangles
          for(size_t i = 0; i < (size_t)subdiv; ++i) {
            typedef std::unordered_map<std::pair<int, int>, int, pair_hash> edge_map_t;
            typedef edge_map_t::iterator edge_iterator;
            edge_map_t edges;
            std::vector<Point3i> new_triangles(4 * triangles.size());
            for(size_t j = 0; j < triangles.size(); ++j) {
              const Point3i tr = triangles[j];
              int p1 = tr[0], p2 = tr[1], p3 = tr[2];
              std::pair<int, int> e1 = MAKE_EDGE(p1, p2);
              std::pair<int, int> e2 = MAKE_EDGE(p2, p3);
              std::pair<int, int> e3 = MAKE_EDGE(p1, p3);
              int n1, n2, n3;
              edge_iterator found = edges.find(e1);
              if(found == edges.end()) {
                n1 = positions.size();
                positions.push_back(normalized(positions[p1] + positions[p2]));
                edges[e1] = n1;
              } else
                n1 = found->second;
              found = edges.find(e2);
              if(found == edges.end()) {
                n2 = positions.size();
                positions.push_back(normalized(positions[p2] + positions[p3]));
                edges[e2] = n2;
              } else
                n2 = found->second;
              found = edges.find(e3);
              if(found == edges.end()) {
                n3 = positions.size();
                positions.push_back(normalized(positions[p3] + positions[p1]));
                edges[e3] = n3;
              } else
                n3 = found->second;
              new_triangles[4 * j] = Point3i(p1, n1, n3);
              new_triangles[4 * j + 1] = Point3i(p2, n2, n1);
              new_triangles[4 * j + 2] = Point3i(p3, n3, n2);
              new_triangles[4 * j + 3] = Point3i(n1, n2, n3);
            }
            swap(triangles, new_triangles);
          }

          Information::out << "# triangles = " << triangles.size() << endl;
          Information::out << "# points = " << positions.size() << endl;

    #undef MAKE_EDGE
        }
      }
      if(subdiv < 0)
        return;

      vvGraph& S = mesh->graph();
      std::vector<vertex> vs(positions.size(), vertex(0));
      Point3f p1 = res.ev[0] * res.p1;
      Point3f p2 = res.ev[1] * res.p2;
      Point3f p3 = res.ev[2] * res.p3;
      Point3f center = res.mean;

      for(size_t i = 0; i < vs.size(); ++i) {
        Point3f p = positions[i];
        vertex v;
        v->pos = Point3d(center + p.x() * p1 + p.y() * p2 + p.z() * p3);
        v->label = lab;
        vs[i] = v;
      }

      Mesh::meshFromTriangles(S, vs, triangles);
    }
  }

  bool PCAnalysis::run(Stack* stk, Store* store, Mesh* m, QString filename, int threshold, QString shape,
                              Point3f correcting_factor, int slices, bool& draw_result)
  {
    Information::out << "Start PCAnalysis" << endl;

    const HVecUS& data = store->data();
    // Move the stack w.r.t. the PCs

    std::vector<PCAnalysis_result> result;
    std::vector<ushort> labels;
    if(store->labels()) {
      std::unordered_set<ushort> labs;
      forall(const ushort& l, data) {
        if(l > 0)
          labs.insert(l);
      }
      labels = std::vector<ushort>(labs.begin(), labs.end());

      result = analyzePC(stk, data, SelectLabel(labels), correcting_factor);
    } else
      result = analyzePC(stk, data, SelectThreshold(threshold), correcting_factor);

    if(shape == "Cuboids" or shape == "Cylinder" or shape == "Ellipsoids") {
      Information::out << "Drawing " << shape << endl;
      draw_result = true;
      Mesh* mesh = this->mesh(stk->id());
      mesh->reset();
      void (*shape_fct)(int, const PCAnalysis_result&, Mesh*, int) = 0;
      if(shape == "Cuboids")
        shape_fct = addPCACuboidShape;
      else if(shape == "Cylinder")
        shape_fct = addPCACylinderShape;
      else if(shape == "Ellipsoids")
        shape_fct = addPCAEllipsoidShape;
      for(size_t vid = 0; vid < labels.size(); ++vid) {
        if(result[vid]) {
          shape_fct(labels[vid], result[vid], mesh, slices);
        }
      }
      // Clear memory used by addPCAEllipsoidShape
      addPCAEllipsoidShape(0, PCAnalysis_result(), 0, -1);
      setNormals(mesh->graph());
      mesh->updateAll();
    } else {
      if(result.size() == 1) {
        stk->trans().setPositionAndOrientation(qglviewer::Vec(0, 0, 0), qglviewer::Quaternion(0, 0, 0, 1));
        stk->setShowTrans(true);

        CuttingSurface* cf = cuttingSurface();
        qglviewer::Frame& frame = cf->frame();
        // Information::out << "Mean = " << mean << endl;

        PCAnalysis_result res = result[0];
        if(!res) {
          setErrorMessage("Error, not result found, do you have any point?");
          return false;
        }
        double mm[16] = { res.p1.x(), res.p2.x(), res.p3.x(), 0, res.p1.y(),   res.p2.y(),   res.p3.y(),   0,
                          res.p1.z(), res.p2.z(), res.p3.z(), 0, res.mean.x(), res.mean.y(), res.mean.z(), 1 };
        Matrix4d mm_gl(mm);
        Information::out << "Transform matrix = " << mm_gl << endl;
        frame.setFromMatrix(mm_gl.c_data());

        // Now, set the planes
        cf->setMode(CuttingSurface::THREE_AXIS);
        cf->setSize(Point3f(res.ev));
        cf->showGrid();
      }

      draw_result = false;
    }

    if(!filename.isEmpty()) {
      QFile f(filename);
      if(!f.open(QIODevice::WriteOnly)) {
        setErrorMessage(QString("Error, cannot open file '%1' for writing.").arg(filename));
        return false;
      }
      QTextStream ts(&f);
      ts << "Image file," << store->file() << endl;
      if(store->labels()) {
        ts << QString("Label, pos.X (%1), pos.Y (%1), pos.Z (%1), e1.X, e1.Y, e1.Z, e2.X, e2.Y, e2.Z, "
                      "r1 (%1), r2 (%1), r3 (%1)").arg(UM) << endl;
        for(size_t vid = 0; vid < labels.size(); ++vid) {
          const PCAnalysis_result& res = result[vid];
          ts << labels[vid] << ", " << res.mean.x() << ", " << res.mean.y() << ", " << res.mean.z() << ", "
             << res.p1.x() << ", " << res.p1.y() << ", " << res.p1.z() << ", " << res.p2.x() << ", " << res.p2.y()
             << ", " << res.p2.z() << ", " << res.ev[0] << ", " << res.ev[1] << ", " << res.ev[2] << endl;
        }
      } else {
        const PCAnalysis_result& res = result[0];
        ts << "Threshold," << threshold << endl;
        ts << "Radii";
        for(int i = 0; i < 3; ++i)
          ts << "," << res.ev[i];
        ts << endl;
        ts << "Center";
        for(int i = 0; i < 3; ++i)
          ts << "," << res.mean[i];
        ts << "Eigenvectors" << endl;
        ts << QString("X (%1), Y (%1), Z (%1)").arg(UM) << endl;
        ts << res.p1.x() << "," << res.p1.y() << "," << res.p1.z() << endl;
        ts << res.p2.x() << "," << res.p2.y() << "," << res.p2.z() << endl;
        ts << res.p3.x() << "," << res.p3.y() << "," << res.p3.z() << endl;
        ts << endl;
      }
      f.close();
      Information::out << "Written result in file " << filename << endl;
    }

    for(size_t vid = 0; vid < labels.size(); ++vid) {
      const PCAnalysis_result& res = result[vid];

      SymmetricTensor& tensor = m->cellAxis()[labels[vid]];
      tensor.ev1() = res.p1;
      tensor.ev2() = res.p2;
      tensor.evals() = res.ev;
    }

    //if(!runDisplay) return true;
    m->setCellAxisType("PCA");

    // Run PCA display, with the parameters from the GUI
    DisplayPCA3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Shape Analysis/Display Shape Axis 3D", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPCA3D process"));
    proc->run();


    return true;
  }

  REGISTER_PROCESS(PCAnalysis);

typedef Vector<3, vertex> triangle;

  bool PCAnalysis2D::run(Mesh* m1, QString filename, double axisLineScale, double axisLineWidth, QString outp, Point3f correcting_factor, IntFloatAttr &outputHeatMap, bool runDisplay)
  {

    vvGraph& S = m1->graph();

    m1->updateCentersNormals();

    // get all vertices of a cell
    std::map<int, std::set<Point3d> > labelPointMap;
    std::map<int, Point3d> centroids;
    std::map<int, double> labelAreaMap;

    std::unordered_map<int, std::set<triangle> > TriangleMap;
    std::map<int, int> labelEntryMap;
    typedef std::pair<int,int> IntIntP;

    // List all the internal triangles (non border) first.
    forall(const vertex &v, S) {
      forall(const vertex &n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
        int label = getLabel(v, n, m);
        if(m1->useParents()) label = m1->parents()[label];
        if(label <= 0)
          continue;

        float area = triangleArea(v->pos, n->pos, m->pos);
        Point3d triCenter = (v->pos + n->pos + m->pos)/3.;

        //Areas[label] += area;
        labelAreaMap[label] += area;
        centroids[label] += area * triCenter;
        if(n->minb == 0 and v->minb == 0 and m->minb == 0) {
          //InsideAreas[label] += area;
          triangle t(v,n,m);
          TriangleMap[label].insert(t);

        }
      }
    }
    // normalize cell centers
    for(std::map<int, Point3d>::iterator it = centroids.begin(); it!=centroids.end(); it++){
      it->second /= labelAreaMap[it->first];
    }

    std::unordered_map<int, std::set<triangle> >::iterator it;
    for(it = TriangleMap.begin() ; it != TriangleMap.end() ; ++it) {
    }

    int nb_values = TriangleMap.size();
    std::vector<Matrix3f> corr(nb_values);
    std::vector<PCAnalysis_result> result(nb_values);

    int counter = 0;

    // assemble correlation matrix
    for(it = TriangleMap.begin() ; it != TriangleMap.end() ; ++it) {
      std::set<triangle> Triangles = it->second;
      std::set<triangle>::iterator t;
      int vid = counter;
      int label;
      for(t = Triangles.begin() ; t != Triangles.end() ; ++t) {
        triangle tri = *t;
        vertex v = tri[0];
        vertex n = tri[1];
        vertex m = tri[2];
        Point3d triCenter = (v->pos + n->pos + m->pos)/3.;
        label = getLabel(v, n, m);
        if(m1->useParents()) label = m1->parents()[label];
        float area = triangleArea(v->pos, n->pos, m->pos);

        Point3d dp = triCenter - centroids[label];
        labelEntryMap[label] = vid;
        corr[vid](0, 0) += dp[0] * dp[0] * area;
        corr[vid](1, 1) += dp[1] * dp[1] * area;
        corr[vid](2, 2) += dp[2] * dp[2] * area;
        corr[vid](1, 0) += dp[0] * dp[1] * area;
        corr[vid](2, 0) += dp[0] * dp[2] * area;
        corr[vid](2, 1) += dp[1] * dp[2] * area;
      }
        corr[vid](0, 0) /= labelAreaMap[label];
        corr[vid](1, 1) /= labelAreaMap[label];
        corr[vid](2, 2) /= labelAreaMap[label];
        corr[vid](1, 0) /= labelAreaMap[label];
        corr[vid](2, 0) /= labelAreaMap[label];
        corr[vid](2, 1) /= labelAreaMap[label];

      counter++;

    }

    gsl_matrix* mat = gsl_matrix_alloc(3, 3);
    gsl_vector* eval = gsl_vector_alloc(3);
    gsl_matrix* evec = gsl_matrix_alloc(3, 3);
    gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(3);
    Information::out << "Decompose matrices" << endl;
      // Eigen-decomposition of the matrices
    forall(IntIntP p, labelEntryMap) {
        int vid = p.second;

        corr[vid](0, 1) = corr[vid](1, 0);
        corr[vid](0, 2) = corr[vid](2, 0);
        corr[vid](1, 2) = corr[vid](2, 1);
        corr[vid] /= TriangleMap[p.first].size();

        for(int i = 0; i < 3; ++i)
          for(int j = 0; j < 3; ++j)
            gsl_matrix_set(mat, i, j, corr[vid](i, j));
        gsl_eigen_symmv(mat, eval, evec, w);
        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

        Point3f p1(gsl_matrix_get(evec, 0, 0), gsl_matrix_get(evec, 1, 0), gsl_matrix_get(evec, 2, 0));
        Point3f p2(gsl_matrix_get(evec, 0, 1), gsl_matrix_get(evec, 1, 1), gsl_matrix_get(evec, 2, 1));
        Point3f p3(gsl_matrix_get(evec, 0, 2), gsl_matrix_get(evec, 1, 2), gsl_matrix_get(evec, 2, 2));
        Point3f ev(gsl_vector_get(eval, 0), gsl_vector_get(eval, 1), gsl_vector_get(eval, 2));
        if((p1 ^ p2) * p3 < 0)
          p3 = -p3;

        // Correct eigen-vector for shape-factor
        //if(correcting_factor.x() > 0)
        double (*sqrtd)(double) = &std::sqrt; // FIXME: this declaration allows the compiler to find right template of map()
        ev = Point3f(multiply(Point3d(correcting_factor), map(sqrtd, Point3d(ev))));
        // Information::out << "Size of the various dimension = " << ev << endl;

        result[vid].p1 = p1;
        result[vid].p2 = p2;
        result[vid].p3 = p3;
        result[vid].ev = ev * /*sqrt(labelAreaMap[p.first])*/ sqrt(TriangleMap[p.first].size());
        result[vid].mean = Point3f(centroids[p.first]);
        // Information::out << "Value: " << selection.value(vid) << " -- Major axis = " << p1 << " eigenvalues = "
        // << ev << endl;
      }

      //if(correcting_factor.x() == 0) {
        Information::out << "Compute regions spans" << endl;
        std::vector<Point3f> min_pos(nb_values, Point3f(FLT_MAX));
        std::vector<Point3f> max_pos(nb_values, Point3f(-FLT_MAX));
        // Find the maximum span of each label along each direction


      gsl_eigen_symmv_free(w);
      gsl_vector_free(eval);
      gsl_matrix_free(evec);
      gsl_matrix_free(mat);


    if(!filename.isEmpty()) {
      QFile f(filename);
      if(!f.open(QIODevice::WriteOnly)) {
        setErrorMessage(QString("Error, cannot open file '%1' for writing.").arg(filename));
        return false;
      }
      QTextStream ts(&f);
      ts << "Image file," << endl;

        ts << QString("Label, pos.X (%1), pos.Y (%1), pos.Z (%1), e1.X, e1.Y, e1.Z, e2.X, e2.Y, e2.Z, "
                      "r1 (%1), r2 (%1), r3 (%1)").arg(UM) << endl;
        forall(IntIntP p, labelEntryMap) {
          int vid = p.second;
          int label = p.first;
          const PCAnalysis_result& res = result[vid];
          ts << label << ", " << res.mean.x() << ", " << res.mean.y() << ", " << res.mean.z() << ", "
             << res.p1.x() << ", " << res.p1.y() << ", " << res.p1.z() << ", " << res.p2.x() << ", " << res.p2.y()
             << ", " << res.p2.z() << ", " << res.ev[0] << ", " << res.ev[1] << ", " << res.ev[2] << endl;
        }

      f.close();
      Information::out << "Written result in file " << filename << endl;
    }

    AttrMap<int, SymmetricTensor>& shapeAttr = m1->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor Shape2D");
    shapeAttr.clear();

    outputHeatMap.clear();
    forall(IntIntP p, labelEntryMap) {
      int vid = p.second;
      const PCAnalysis_result& res = result[vid];
      SymmetricTensor& tensor = m1->cellAxis()[p.first];
      tensor.ev1() = res.p1;//pca.p1;
      tensor.ev2() = res.p2;
      tensor.evals() = res.ev;
      shapeAttr[p.first] = tensor;
      // write to map
      if(outp == "ratio"){
        outputHeatMap[p.first] = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
      } else if(outp == "max"){
        outputHeatMap[p.first] = max(tensor.evals()[0],tensor.evals()[1]);
      } else if(outp == "min"){
        outputHeatMap[p.first] = min(tensor.evals()[0], tensor.evals()[1]);
      }
    }

		if(!runDisplay) return true;
		m1->setCellAxisType("PCA");

		// Run PCA display, with the parameters from the GUI

    DisplayPCAOld *proc;
    if(!getProcess("Mesh/Cell Axis/Shape Analysis/Display Shape Axis", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPCA process"));
  	proc->run();

    return true;

  }
	REGISTER_PROCESS(PCAnalysis2D);

  void DisplayPCAOld::processParms()
  {
    DisplayHeatMap = parm("Heatmap");
    scaleHeatMap = parm("ScaleHeat");
    RangeHeat = Point2d(parm("Heat min").toDouble(), parm("Heat max").toDouble());
    DisplayAxis = parm("Show Axis");
    AxisColor = parm("Axis Color");
    AxisWidth = parm("Line Width").toFloat();
    AxisScale = parm("Line Scale").toFloat();
    AxisOffset = parm("Line Offset").toFloat();
    AspectRatioThreshold = parm("Threshold").toFloat();
    CreateAttrs = stringToBool(parm("Create Attr Maps"));
  }

	bool DisplayPCAOld::initialize(QWidget* parent)
  {
	  if(!checkState().mesh(MESH_NON_EMPTY))
		  return false;

	  // Check values of parameters and change them in the GUI if needed.
		// const IntSymTensorAttr& cellAxisCustom = currentMesh()->attributes().attrMap<int, SymmetricTensor>("Custom Directions");
		// // If PCA not projected on Bezier yet, make sure we won't try to display the projection:
		// if(cellAxisCustom.size() == 0) {
    // QString heatMap = parm("Heatmap");
		// if(heatMap == "ShearPCABezier" or heatMap == "PCABezierX" or heatMap == "PCABezierY")
		//   setParm( "Heatmap", "PCAMax");
    //
    // QString axis = parm("Show Axis");
		// if(axis == "PCABezierX" or axis == "PCABezierY")
		//   setParm("Show Axis", "PCAMax");
		// }

    return true;
  }

  // Display PCA (separatly from each other) and heatmaps of aspectRatio etc.
	// The stretch vectors are stored in cellAxis, so that the norm of the vector give us the stretch (always positive)
	// The strain (= stretch -1) vectors are visualized in cellAxisVis. No deformation ==> strain = 0 ==> cell axis null
  bool DisplayPCAOld::run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
	  const QString displayAxis, const QColor& axisColor, float axisWidth, float axisScale,
		float axisOffset, float aspectRatioThreshold, bool customShear, bool createAttrMaps)
  {
    // if there is no PCA stored, display error message
    const IntSymTensorAttr& cellAxis = mesh->cellAxis();
    IntSymTensorAttr& cellAxisCustom = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    AttrMap<int, double>& shapeProduct = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Product");
    AttrMap<int, double>& shapeMax = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Max");
    AttrMap<int, double>& shapeMin = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Min");
    AttrMap<int, double>& shapeAniso = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/ShapeAnisotropy");
    AttrMap<int, double>& shapeCustomX = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomX");
    AttrMap<int, double>& shapeCustomY = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomY");
    AttrMap<int, double>& shapeCustomZ = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomZ");


    if(cellAxis.size() == 0)
      throw(QString("No cell axis stored in active mesh"));
    if(mesh->cellAxisType() != QString("PCA"))
      throw(QString("No PCA stored in active mesh"));

    // Scale PCA for visualization
    mesh->setAxisWidth(axisWidth);
    mesh->setAxisOffset(axisOffset);
    Colorb colorExpansion(axisColor);

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

    // Scale the cell axis for display
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();

    std::vector<float> vecValues;
    IntFloatAttr& labelHeatMap = mesh->labelHeat();
    if(displayHeatMap != "None") {
      mesh->labelHeat().clear();
    }

    forall(const IntSymTensorPair &p, cellAxis) {
      int cell = p.first;
      SymmetricTensor tensor = p.second;
      float pca1st, pca2nd, aspectRatio;
			Point3f dir1st, dir2nd;

      // if any "custom" option is selected, calculate the custom directions
      if(displayAxis == "CustomX" or displayAxis == "CustomY" or displayAxis == "BothCustom" or
        displayHeatMap == "CustomX" or displayHeatMap == "CustomY" or displayHeatMap == "Shear"
        or displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio" or createAttrMaps){

        if(cellAxisCustom.empty())
          throw(QString("No custom cell axis. Create a custom cell axis first."));
        if(cellAxisCustom.find(cell) == cellAxisCustom.end())
          continue;

        Point3f customX = cellAxisCustom[cell].ev1();
        Point3f customY = cellAxisCustom[cell].ev2();
        cellAxisCustom[cell] = customDirectionTensor(tensor, customX, customY, customShear);

      }

      double axisLengthMax = tensor.evals()[0];
      double axisLengthMin = tensor.evals()[1];
      double axisLengthCustomX = cellAxisCustom[cell].evals()[0];
      double axisLengthCustomY = cellAxisCustom[cell].evals()[1];
      double axisLengthCustomZ = cellAxisCustom[cell].evals()[2];

      // if custom cell axis required, grab it
      if(displayAxis == "CustomX" or displayAxis == "CustomY" or displayAxis == "BothCustom"){
        tensor = cellAxisCustom[cell];
      }

      pca1st = tensor.evals()[0];
      pca2nd = tensor.evals()[1];
			dir1st = tensor.ev1();
			dir2nd = tensor.ev2();
      aspectRatio = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);

      Point3b showAxis;
      showAxis[0] = (displayAxis == "Max" or displayAxis == "CustomX" or displayAxis == "Both" or displayAxis == "BothCustom");
      showAxis[1] = (displayAxis == "Min" or displayAxis == "CustomY" or displayAxis == "Both" or displayAxis == "BothCustom");
      showAxis[2] = false;

      if(aspectRatioThreshold != 0 and aspectRatio < aspectRatioThreshold)
        showAxis = Point3b(false, false, false);

			cellAxisVis[cell][0] = (showAxis[0] ? pca1st : 0.f) * dir1st * axisScale;
      cellAxisVis[cell][1] = (showAxis[1] ? pca2nd : 0.f) * dir2nd * axisScale;
      cellAxisVis[cell][2] = Point3f(0, 0, 0);

      cellAxisColor[cell][0] = Colorb(axisColor);
      cellAxisColor[cell][1] = Colorb(axisColor);
      cellAxisColor[cell][2] = Colorb(axisColor);


	  // Fill in heat map table
    
		if(displayHeatMap != "None") {
      //mesh->labelHeat().clear();
      //IntFloatAttr& labelHeatMap = mesh->labelHeat();
      //forall(const IntSymTensorPair &p, cellAxis) {
        //int cell = p.first;
        //SymmetricTensor tensor = p.second;
        float value = 0;

        if(displayHeatMap == "Max")
          value = max(tensor.evals()[0], tensor.evals()[1]) *2;
        else if(displayHeatMap == "Min")
          value = min(tensor.evals()[0], tensor.evals()[1]) *2;
        else if(displayHeatMap == "Ratio Max/Min")
          value = max(tensor.evals()[0], tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
        else if(displayHeatMap == "Product Max*Min")
          value = tensor.evals()[0] * tensor.evals()[1] * 4; // times 4 to make it equivalent to the area
        else if(displayHeatMap == "Anisotropy Max/(Max+Min)" and tensor.evals()[1] != 0)
          value = tensor.evals()[0] / (tensor.evals()[0]+tensor.evals()[1]); // times 4 to make it equivalent to the area

        // if custom heat map is required, grab it
        if(displayHeatMap == "CustomX" or displayHeatMap == "CustomY" or displayHeatMap == "Shear" or
          displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio" or displayHeatMap == "RatioCustomYX"){
          if(cellAxisCustom.size() == 0)
            throw(QString("No custom cell axis. Run 'Cell Axis Custom Directions' first."));
          if(cellAxisCustom.find(cell) == cellAxisCustom.end())
            continue;
          tensor = cellAxisCustom.find(cell)->second;
          if(displayHeatMap == "CustomX"){
            value = tensor.evals()[0]*2;;
          } else if(displayHeatMap == "CustomY"){
            value = tensor.evals()[1]*2;;
          } else if(displayHeatMap == "Shear"){
            value = fabs(tensor.evals()[2]);
          } else if(displayHeatMap == "AnisoCustom"){
            double minT = min(tensor.evals()[0], tensor.evals()[1]);
            if(minT != 0) value = max(tensor.evals()[0], tensor.evals()[1])/minT;
          } else if(displayHeatMap == "RatioCustomYX"){
            value = tensor.evals()[0] != 0 ? tensor.evals()[1]/tensor.evals()[0] : 0;
          } else if(displayHeatMap == "AnisoRatio"){
            double minT = min(p.second.evals()[0], p.second.evals()[1]);
            double minC = min(tensor.evals()[0], tensor.evals()[1]);
            if(minT != 0 and minC != 0){
              double aniso = max(p.second.evals()[0], p.second.evals()[1])/minT;
              double anisoCus = max(tensor.evals()[0], tensor.evals()[1])/minC;
              value = max(aniso, anisoCus)/min(aniso, anisoCus);
            }

          }
        }

        labelHeatMap[cell] = value;
		   	vecValues.push_back(value);
		}

      if(createAttrMaps){
        shapeProduct[cell] = axisLengthMax * axisLengthMin *4;
        shapeMax[cell] = axisLengthMax *2;
        shapeMin[cell] = axisLengthMin *2;
        shapeAniso[cell] = axisLengthMax / axisLengthMin;
        shapeCustomX[cell] = axisLengthCustomX *2;
        shapeCustomY[cell] = axisLengthCustomY *2;
        shapeCustomZ[cell] = axisLengthCustomZ *2;
      }
    }

    // Find bounds for heatmap scale.
		float lowBound = 0, highBound = 1;
		// Manual scaling by user input
		if(scaleHeatMap == "Manual"){
		  lowBound  = rangeHeat.x();
		  highBound = rangeHeat.y();
		}
		// Automatic scaling
	  if(vecValues.size() != 0){
		  std::sort(vecValues.begin(), vecValues.end());
	  	// Default scaling: take min and max values
	  	if(scaleHeatMap == "None"){
	  	  lowBound = vecValues.front();
	  	  highBound = vecValues.back();
	  	}
	  	// Auto scale: take low 5 and high 5 percentile
	  	else if(scaleHeatMap == "Auto"){
	  		lowBound = *(vecValues.begin()+int(vecValues.size()*.05));
	  		highBound = *(vecValues.begin()+int(vecValues.size()*.95));
	  	}
		}
    // Round the bounds values of heat map, to make sure we have the same bound in case
  	// of loading PCA on both meshes.
    lowBound = floor(lowBound * 1000) / 1000;
    highBound = ceil(highBound * 1000) / 1000;

    mesh->heatMapBounds() = Point2f(lowBound, highBound);
    mesh->heatMapUnit() = UM;
    mesh->setShowLabel("Label Heat");
    mesh->updateTriangles();

    if(displayAxis != "None")
      mesh->setShowAxis("Cell Axis");

    return true;
  }
  REGISTER_PROCESS(DisplayPCAOld);


  // Display PCA (separatly from each other) and heatmaps of aspectRatio etc.
  // The stretch vectors are stored in cellAxis, so that the norm of the vector give us the stretch (always positive)
  // The strain (= stretch -1) vectors are visualized in cellAxisVis. No deformation ==> strain = 0 ==> cell axis null
  bool DisplayPCA::run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
    const QString displayAxis, const QColor& axisColor, float axisWidth, float axisScale,
    float axisOffset, float aspectRatioThreshold, bool customShear, bool createAttrMaps)
  {
    // if there is no PCA stored, display error message
    // const IntSymTensorAttr& cellAxis = mesh->cellAxis();
    // IntSymTensorAttr& cellAxisCustom = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");

    AttrMap<int, double>& shapeProduct = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Product");
    AttrMap<int, double>& shapeMax = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Max");
    AttrMap<int, double>& shapeMin = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/Min");
    AttrMap<int, double>& shapeAniso = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/ShapeAnisotropy");
    AttrMap<int, double>& shapeCustomX = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomX");
    AttrMap<int, double>& shapeCustomY = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomY");
    AttrMap<int, double>& shapeCustomZ = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis/CustomZ");

    IntSymTensorAttr &cellAxis = mesh->cellAxis();

    QString mode = "";
    bool custom = false;

    if(displayAxis == "Both"){
      mode = "XY";
    } else if(displayAxis == "Max"){
      mode = "X";
    } else if(displayAxis == "Min"){
      mode = "Y";
    } else if(displayAxis == "CustomX"){
      mode = "X";
      custom = true;
    } else if(displayAxis == "CustomY"){
      mode = "Y";
      custom = true;
    } else if(displayAxis == "BothCustom"){
      mode = "XY";
      custom = true;
    } else if(displayAxis == "None"){
      mode = "None";
    } else {
      return setErrorMessage("Invalid Directions");
    }

    QColor qaxisColorX = axisColor;
    QColor qaxisColorY = axisColor;
    QColor qaxisColorZ = axisColor;
    bool strain = false;


    DisplayCellAxisFromTensor displayAxisP(*this);
    displayAxisP.run(mesh, cellAxis, mode, custom, qaxisColorX, qaxisColorY, qaxisColorZ, qaxisColorX, axisWidth, axisScale, axisOffset,
      aspectRatioThreshold, strain);

    // Fill in heat map table
    
    // if(displayHeatMap != "None") {
    //   forall(const IntSymTensorPair &p, cellAxis) {
    //   //mesh->labelHeat().clear();
    //   //IntFloatAttr& labelHeatMap = mesh->labelHeat();
    //   //forall(const IntSymTensorPair &p, cellAxis) {
    //     //int cell = p.first;
    //     //SymmetricTensor tensor = p.second;
    //     float value = 0;

    //     if(displayHeatMap == "Max")
    //       value = max(tensor.evals()[0], tensor.evals()[1]) *2;
    //     else if(displayHeatMap == "Min")
    //       value = min(tensor.evals()[0], tensor.evals()[1]) *2;
    //     else if(displayHeatMap == "Ratio Max/Min")
    //       value = max(tensor.evals()[0], tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
    //     else if(displayHeatMap == "Product Max*Min")
    //       value = tensor.evals()[0] * tensor.evals()[1] * 4; // times 4 to make it equivalent to the area
    //     else if(displayHeatMap == "Anisotropy Max/(Max+Min)" and tensor.evals()[1] != 0)
    //       value = tensor.evals()[0] / tensor.evals()[1]; // times 4 to make it equivalent to the area

    //     // if custom heat map is required, grab it
    //     if(displayHeatMap == "CustomX" or displayHeatMap == "CustomY" or displayHeatMap == "Shear" or
    //       displayHeatMap == "AnisoCustom" or displayHeatMap == "AnisoRatio" or displayHeatMap == "RatioCustomYX"){
    //       if(cellAxisCustom.size() == 0)
    //         throw(QString("No custom cell axis. Run 'Cell Axis Custom Directions' first."));
    //       if(cellAxisCustom.find(cell) == cellAxisCustom.end())
    //         continue;
    //       tensor = cellAxisCustom.find(cell)->second;
    //       if(displayHeatMap == "CustomX"){
    //         value = tensor.evals()[0];
    //       } else if(displayHeatMap == "CustomY"){
    //         value = tensor.evals()[1];
    //       } else if(displayHeatMap == "Shear"){
    //         value = fabs(tensor.evals()[2]);
    //       } else if(displayHeatMap == "AnisoCustom"){
    //         double minT = min(tensor.evals()[0], tensor.evals()[1]);
    //         if(minT != 0) value = max(tensor.evals()[0], tensor.evals()[1])/minT;
    //       } else if(displayHeatMap == "RatioCustomYX"){
    //         value = tensor.evals()[0] != 0 ? tensor.evals()[1]/tensor.evals()[0] : 0;
    //       } else if(displayHeatMap == "AnisoRatio"){
    //         double minT = min(p.second.evals()[0], p.second.evals()[1]);
    //         double minC = min(tensor.evals()[0], tensor.evals()[1]);
    //         if(minT != 0 and minC != 0){
    //           double aniso = max(p.second.evals()[0], p.second.evals()[1])/minT;
    //           double anisoCus = max(tensor.evals()[0], tensor.evals()[1])/minC;
    //           value = max(aniso, anisoCus)/min(aniso, anisoCus);
    //         }

    //       }
    //     }

    //     labelHeatMap[cell] = value;
    //     vecValues.push_back(value);
    // }

    //   if(createAttrMaps){
    //     shapeProduct[cell] = axisLengthMax * axisLengthMin *4;
    //     shapeMax[cell] = axisLengthMax *2;
    //     shapeMin[cell] = axisLengthMin *2;
    //     shapeAniso[cell] = axisLengthMax / axisLengthMin;
    //     shapeCustomX[cell] = axisLengthCustomX *2;
    //     shapeCustomY[cell] = axisLengthCustomY *2;
    //     shapeCustomZ[cell] = axisLengthCustomZ *2;
    //   }
    // }

    // // Find bounds for heatmap scale.
    // float lowBound = 0, highBound = 1;
    // // Manual scaling by user input
    // if(scaleHeatMap == "Manual"){
    //   lowBound  = rangeHeat.x();
    //   highBound = rangeHeat.y();
    // }
    // // Automatic scaling
    // if(vecValues.size() != 0){
    //   std::sort(vecValues.begin(), vecValues.end());
    //   // Default scaling: take min and max values
    //   if(scaleHeatMap == "None"){
    //     lowBound = vecValues.front();
    //     highBound = vecValues.back();
    //   }
    //   // Auto scale: take low 5 and high 5 percentile
    //   else if(scaleHeatMap == "Auto"){
    //     lowBound = *(vecValues.begin()+int(vecValues.size()*.05));
    //     highBound = *(vecValues.begin()+int(vecValues.size()*.95));
    //   }
    // }
    // // Round the bounds values of heat map, to make sure we have the same bound in case
    // // of loading PCA on both meshes.
    // lowBound = floor(lowBound * 1000) / 1000;
    // highBound = ceil(highBound * 1000) / 1000;

    // mesh->heatMapBounds() = Point2f(lowBound, highBound);
    // mesh->heatMapUnit() = QString("a.u.");
    // mesh->setShowLabel("Label Heat");
    // mesh->updateTriangles();

    if(displayAxis != "None")
      mesh->setShowAxis("Cell Axis");

    return true;
  }
  REGISTER_PROCESS(DisplayPCA);

  std::unordered_map<int, std::set<triangle> > createLabelTriangleMap(vvGraph& S)
  {
    std::unordered_map<int, std::set<triangle> > TriangleMap;

    // create a label->triangle map
    forall(const vertex &v, S) {
      forall(const vertex &n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
        int label = getLabel(v, n, m);
        if(label <= 0)
          continue;

        triangle t(v,n,m);
        TriangleMap[label].insert(t);

      }
    }
    return TriangleMap;
  }

  void calcCentroidVolume(std::set<triangle>& triSet, Point3d& cellCentroid, double& cellVolume)
  {

    cellCentroid = Point3d(0,0,0);
    cellVolume = 0;

    for(auto t = triSet.begin() ; t != triSet.end() ; ++t) {
      triangle tri = *t;
      vertex v = tri[0];
      vertex n = tri[1];
      vertex m = tri[2];

      double volume = signedTetraVolume(v->pos, n->pos, m->pos);
      cellCentroid += volume * (v->pos + n->pos + m->pos) / 4.0;
      cellVolume += volume;
    }


    cellCentroid /= cellVolume;
  }

  bool PCAnalysis3D::run(Mesh* m, Point3d correcting_factor, QString inputMode)
  {
    Information::out << "Start PCAnalysis" << endl;

    vvGraph& S = m->graph();

    // get all vertices of a cell
    std::map<int, std::set<Point3d> > labelPointMap;
    std::map<int, double> labelAreaMap;
    AttrMap<int, Point3d> &centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");

    std::unordered_map<int, std::set<triangle> > TriangleMap = createLabelTriangleMap(S);
    std::map<int, int> labelEntryMap;
    typedef std::pair<int,int> IntIntP;

    int nb_values = TriangleMap.size();
    std::vector<Matrix3f> corr(nb_values);
    std::vector<PCAnalysis_result> result(nb_values);

    std::map<int,int> sum;

    int counter = 0;

    // assemble correlation matrix
    if(inputMode == "Stack Voxels"){

      Stack* stk = currentStack();
      Store* store = stk->currentStore();
      const HVecUS& data = store->data();

      // First, compute mean
      
      std::set<int> nb_labels;

      // Compute the mean for each selected value
      Information::out << "Compute mean" << endl;
      Point3u s = stk->size();
      size_t k = 0;
      for(size_t z = 0; z < s.z(); ++z)
        for(size_t y = 0; y < s.y(); ++y)
          for(size_t x = 0; x < s.x(); ++x, ++k) {
            int vid = data[k];
            if(vid > 0) {
              nb_labels.insert(vid);
              if(labelEntryMap.find(vid)==labelEntryMap.end()){
                labelEntryMap[vid]=counter;
                counter++;
              }
              sum[vid]++;
              centroids[vid] += Point3d(stk->imageToWorld(Point3d(x, y, z)));
            }
          }
      forall(auto p, sum) {
        if(p.second > 0)
          centroids[p.first] /= double(p.second);
      }
      Information::out << "Compute CC matrix " << nb_labels.size() << endl;

      for(auto p : centroids){
        Information::out << "Label/Pos " << p.first << "/" << p.second << endl;
      }

      nb_values = nb_labels.size();
      corr = std::vector<Matrix3f>(nb_values);
      result = std::vector<PCAnalysis_result>(nb_values);

Information::out << "corr " << corr.size() << endl;

      // Compute the cross-correlation matrix
      k = 0;
      for(size_t z = 0; z < s.z(); ++z)
        for(size_t y = 0; y < s.y(); ++y)
          for(size_t x = 0; x < s.x(); ++x, ++k) {
            int label = data[k];
            if(label > 0) {
              int vid = labelEntryMap[label];
              Point3f dp = Point3f(Point3d(stk->imageToWorld(Point3f(x, y, z))) - centroids[label]);
              corr[vid](0, 0) += dp[0] * dp[0];
              corr[vid](1, 1) += dp[1] * dp[1];
              corr[vid](2, 2) += dp[2] * dp[2];
              corr[vid](1, 0) += dp[0] * dp[1];
              corr[vid](2, 0) += dp[0] * dp[2];
              corr[vid](2, 1) += dp[1] * dp[2];
            }
          }
Information::out << "Matrix done" << endl;
    } else { // inputMode == Mesh Triangles

      for(auto it = TriangleMap.begin() ; it != TriangleMap.end() ; ++it) {
        std::set<triangle> Triangles = it->second;
        int vid = counter;
        int label = it->first;

        double cellVol = 0;
        Point3d cellCentroid(0,0,0);

        calcCentroidVolume(Triangles, cellCentroid, cellVol);

        centroids[label] = cellCentroid;

        double totalWeight = 0.0;

        for(auto t = Triangles.begin() ; t != Triangles.end() ; ++t) {
          triangle tri = *t;
          vertex v = tri[0];
          vertex n = tri[1];
          vertex m = tri[2];
          Point3d triCenter = (v->pos + n->pos + m->pos)/3.;

          double tetraVol = ((v->pos-cellCentroid) * ((n->pos-cellCentroid) ^ (m->pos-cellCentroid))) / 6.;

          label = getLabel(v, n, m);
          float area = triangleArea(v->pos, n->pos, m->pos);
          float signal = (v->signal + n->signal + m->signal)/3.;

          double weight = area; //tetraVol;//1.0;
          // if(parm("Weight") == "Area") weight = area;
          // else if(parm("Weight") == "Signal") weight = signal;
          // else if(parm("Weight") == "Area*Signal") weight = area*signal;

          totalWeight += weight;

          Point3d dp = triCenter - centroids[label];

          //if(parm("Normalize Directions") == "Yes") dp /= norm(dp);

          labelEntryMap[label] = vid;
          corr[vid](0, 0) += dp[0] * dp[0] * weight;
          corr[vid](1, 1) += dp[1] * dp[1] * weight;
          corr[vid](2, 2) += dp[2] * dp[2] * weight;
          corr[vid](1, 0) += dp[0] * dp[1] * weight;
          corr[vid](2, 0) += dp[0] * dp[2] * weight;
          corr[vid](2, 1) += dp[1] * dp[2] * weight;
        }
          corr[vid](0, 0) /= totalWeight;
          corr[vid](1, 1) /= totalWeight;
          corr[vid](2, 2) /= totalWeight;
          corr[vid](1, 0) /= totalWeight;
          corr[vid](2, 0) /= totalWeight;
          corr[vid](2, 1) /= totalWeight;

        counter++;

      }

    }

    gsl_matrix* mat = gsl_matrix_alloc(3, 3);
    gsl_vector* eval = gsl_vector_alloc(3);
    gsl_matrix* evec = gsl_matrix_alloc(3, 3);
    gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc(3);
    Information::out << "Decompose matrices" << endl;
      // Eigen-decomposition of the matrices
    forall(IntIntP p, labelEntryMap) {

        int label = p.first;
        int vid = p.second;
        double corFac = 1.;

        corr[vid](0, 1) = corr[vid](1, 0);
        corr[vid](0, 2) = corr[vid](2, 0);
        corr[vid](1, 2) = corr[vid](2, 1);
        if(inputMode == "Stack Voxels"){
          if(sum[label] == 0)
            continue;
          corr[vid] /= sum[label];
        } else { // Mesh Triangles
          corr[vid] /= TriangleMap[label].size();
          corFac = sqrt(TriangleMap[p.first].size());
        }
        
        for(int i = 0; i < 3; ++i)
          for(int j = 0; j < 3; ++j)
            gsl_matrix_set(mat, i, j, corr[vid](i, j));
        gsl_eigen_symmv(mat, eval, evec, w);
        gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_DESC);

        Point3f p1(gsl_matrix_get(evec, 0, 0), gsl_matrix_get(evec, 1, 0), gsl_matrix_get(evec, 2, 0));
        Point3f p2(gsl_matrix_get(evec, 0, 1), gsl_matrix_get(evec, 1, 1), gsl_matrix_get(evec, 2, 1));
        Point3f p3(gsl_matrix_get(evec, 0, 2), gsl_matrix_get(evec, 1, 2), gsl_matrix_get(evec, 2, 2));
        Point3f ev(gsl_vector_get(eval, 0), gsl_vector_get(eval, 1), gsl_vector_get(eval, 2));
        if((p1 ^ p2) * p3 < 0)
          p3 = -p3;

        // Correct eigen-vector for shape-factor
        double (*sqrtd)(double) = &std::sqrt; // FIXME: this declaration allows the compiler to find right template of map()
        ev = Point3f(multiply(correcting_factor, map(sqrtd, Point3d(ev))));
        // Information::out << "Size of the various dimension = " << ev << endl;

        result[vid].p1 = p1;
        result[vid].p2 = p2;
        result[vid].p3 = p3;
        result[vid].ev = ev * corFac;//sqrt(TriangleMap[p.first].size());
        result[vid].mean = Point3f(centroids[p.first]);
        // Information::out << "Value: " << selection.value(vid) << " -- Major axis = " << p1 << " eigenvalues = "
        // << ev << endl;
      }

      //if(correcting_factor.x() == 0) {
        Information::out << "Compute regions spans" << endl;
        std::vector<Point3f> min_pos(nb_values, Point3f(FLT_MAX));
        std::vector<Point3f> max_pos(nb_values, Point3f(-FLT_MAX));
        // Find the maximum span of each label along each direction


      gsl_eigen_symmv_free(w);
      gsl_vector_free(eval);
      gsl_matrix_free(evec);
      gsl_matrix_free(mat);

    //outputHeatMap.clear();
    forall(IntIntP p, labelEntryMap) {
      int vid = p.second;
      const PCAnalysis_result& res = result[vid];
      SymmetricTensor& tensor = m->cellAxis()[p.first];
      tensor.ev1() = res.p1;
      tensor.ev2() = res.p2;
      tensor.evals() = res.ev; // times two to get the full length of the cell
      // write to map
      // if(outp == "ratio"){
      //   outputHeatMap[p.first] = max(tensor.evals()[0],tensor.evals()[1])/min(tensor.evals()[0], tensor.evals()[1]);
      // } else if(outp == "max"){
      //   outputHeatMap[p.first] = max(tensor.evals()[0],tensor.evals()[1]);
      // } else if(outp == "min"){
      //   outputHeatMap[p.first] = min(tensor.evals()[0], tensor.evals()[1]);
      // }
    }


    //if(!runDisplay) return true;
    m->setCellAxisType("PCA");

    // Run PCA display, with the parameters from the GUI
    DisplayPCA3D *proc;
    if(!getProcess("Mesh/Cell Axis 3D/Shape Analysis/Display Shape Axis 3D", proc))
      throw(QString("MeshProcessHeatMap:: Unable to make DisplayPCA3D process"));
    proc->run();

    return true;
  }

  REGISTER_PROCESS(PCAnalysis3D);


  void DisplayPCA3D::processParms()
  {
    DisplayHeatMap = parm("Heatmap");
    ScaleHeatMap = parm("ScaleHeat");
    RangeHeat = Point2d(parm("Heat min").toDouble(), parm("Heat max").toDouble());
    DisplayAxis = parm("Show Axis");
    AxisColor = parm("Axis Color");
    AxisWidth = parm("Line Width").toFloat();
    AxisScale = parm("Line Scale").toFloat();
    AxisOffset = parm("Line Offset").toFloat();
    AspectRatioThreshold = parm("Threshold").toFloat();
    attrMaps = stringToBool(parm("Create Attr Maps"));
  }

  // Display PCA (separatly from each other) and heatmaps of aspectRatio etc.
  // The stretch vectors are stored in cellAxis, so that the norm of the vector give us the stretch (always positive)
  // The strain (= stretch -1) vectors are visualized in cellAxisVis. No deformation ==> strain = 0 ==> cell axis null
  bool DisplayPCA3D::run(Mesh* mesh, const QString displayHeatMap, const QString scaleHeatMap, const Point2d &rangeHeat,
    const QString displayAxis, const QColor& axisColor, float axisWidth, float axisScale,
    float axisOffset, float aspectRatioThreshold, bool attrMaps)
  {
    IntPoint3fAttr &centers = mesh->labelCenterVis();//(mesh->useParents() ? mesh->parentCenterVis() : mesh->labelCenterVis());
    AttrMap<int, SymmetricTensor>& customDirections = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor CustomDirections");
    AttrMap<int, Point3d> &centroids = mesh->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    AttrMap<int, double> &shape3dProduct = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/Product");
    AttrMap<int, double> &shape3dMax = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/Max");
    AttrMap<int, double> &shape3dMid = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/Mid");
    AttrMap<int, double> &shape3dMin = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/Min");
    AttrMap<int, double> &shape3dMaxMid = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/MaxMid");
    AttrMap<int, double> &shape3dMaxMin = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/MaxMin");
    AttrMap<int, double> &shape3dMidMin = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/MidMin");
    AttrMap<int, double> &shape3dAniso = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/ShapeAnisotropy");
    AttrMap<int, double> &shape3dCustomX = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/CustomX");
    AttrMap<int, double> &shape3dCustomY = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/CustomY");
    AttrMap<int, double> &shape3dCustomZ = mesh->attributes().attrMap<int, double>("Measure Label Double ShapeAnalysis3D/CustomZ");

    if(centroids.empty())
      throw QString("%1::run No cell centroids found").arg(name());

    centers.clear();
    for(auto p : centroids)
      centers[p.first] = Point3f(p.second);

    IntFloatAttr& labelHeatMap = mesh->labelHeat();
    if(displayHeatMap != "None")
      labelHeatMap.clear();

    std::vector<float> vecValues;

    // if there is no PCA stored, display error message
    const IntSymTensorAttr& cellAxis = mesh->cellAxis();

    if(cellAxis.size() == 0)
      throw(QString("No cell axis stored in active mesh"));
    if(mesh->cellAxisType() != QString("PCA"))
      throw(QString("No PCA stored in active mesh"));

    // Scale PCA for visualization
    mesh->setAxisWidth(axisWidth);
    mesh->setAxisOffset(axisOffset);
    Colorb colorExpansion(axisColor);

    // Scale the cell axis for display
    IntMatrix3fAttr &cellAxisVis = mesh->cellAxisVis();
    cellAxisVis.clear();
    IntVec3ColorbAttr &cellAxisColor = mesh->cellAxisColor();
    cellAxisColor.clear();

    forall(const IntSymTensorPair &p, cellAxis) {

      int cell = p.first;
      SymmetricTensor tensor = p.second;
      float pca1st, pca2nd, pca3rd, aspectRatio;


      Point3f axisDirX, axisDirY, axisDirZ;
      float axisLengthX, axisLengthY, axisLengthZ;
      Point3f axisDirCustomX, axisDirCustomY, axisDirCustomZ;
      float axisLengthCustomX = -1, axisLengthCustomY = -1, axisLengthCustomZ = -1;

       // if any "custom" option is selected, calculate the custom directions
      if(displayAxis == "Custom X" or displayAxis == "Custom Y" or displayAxis == "Custom Z" or displayAxis == "Custom All" or
        displayHeatMap == "Custom X" or displayHeatMap == "Custom Y" or displayHeatMap == "Custom Z" or (attrMaps and !customDirections.empty())){

        if(customDirections.empty())
          throw(QString("No custom cell axis found. Run 'Cell Axis Custom Directions' first."));
        if(customDirections.find(cell) == customDirections.end())
          continue;

        Point3f customX = customDirections[cell].ev1();
        Point3f customY = customDirections[cell].ev2();
        Point3f customZ = customDirections[cell].evals();

        // length of vectors
        axisLengthCustomX = lengthCustomDirTensor3D(tensor, customX);
        axisLengthCustomY = lengthCustomDirTensor3D(tensor, customY);
        axisLengthCustomZ = lengthCustomDirTensor3D(tensor, customZ);

        // direction of vectors
        axisDirCustomX = customX;
        axisDirCustomX /= (norm(axisDirCustomX));
        axisDirCustomY = customY;
        axisDirCustomY /= (norm(axisDirCustomY));
        axisDirCustomZ = customZ;
        axisDirCustomZ /= (norm(axisDirCustomZ));
      }

      // now the non-custom cell axis
      // length of vectors
      axisLengthX = tensor.evals()[0];
      axisLengthY = tensor.evals()[1];
      axisLengthZ = tensor.evals()[2];

      // direction of vectors
      axisDirX = tensor.ev1();
      axisDirY = tensor.ev2();
      axisDirZ = axisDirX ^ axisDirY;
      axisDirZ /= (norm(axisDirZ));

      // set the axes
      Point3b showAxis;
      showAxis[0] = (displayAxis == "Max" or displayAxis == "Custom X" or displayAxis == "All" or displayAxis == "Custom All");
      showAxis[1] = (displayAxis == "Mid" or displayAxis == "Custom Y" or displayAxis == "All" or displayAxis == "Custom All");
      showAxis[2] = (displayAxis == "Min" or displayAxis == "Custom Z" or displayAxis == "All" or displayAxis == "Custom All");

      // if(aspectRatioThreshold != 0 and aspectRatio < aspectRatioThreshold)
      //   showAxis = Point3b(false, false, false);

      cellAxisVis.erase(cell);
      if(displayAxis == "Custom X" or displayAxis == "Custom Y" or displayAxis == "Custom Z"  or displayAxis == "Custom All"){
        cellAxisVis[cell][0] = (showAxis[0] ? axisLengthCustomX : 0.f) * axisDirCustomX * axisScale;
        cellAxisVis[cell][1] = (showAxis[1] ? axisLengthCustomY : 0.f) * axisDirCustomY * axisScale;
        cellAxisVis[cell][2] = (showAxis[2] ? axisLengthCustomZ : 0.f) * axisDirCustomZ * axisScale;
      } else {
        cellAxisVis[cell][0] = (showAxis[0] ? axisLengthX : 0.f) * axisDirX * axisScale;
        cellAxisVis[cell][1] = (showAxis[1] ? axisLengthY : 0.f) * axisDirY * axisScale;
        cellAxisVis[cell][2] = (showAxis[2] ? axisLengthZ : 0.f) * axisDirZ * axisScale;
      }

      cellAxisColor[cell][0] = Colorb(axisColor);
      cellAxisColor[cell][1] = Colorb(axisColor);
      cellAxisColor[cell][2] = Colorb(axisColor);

    // set the heat map

      if(displayHeatMap != "None") {

        int cell = p.first;
        float value = 0;

        if(displayHeatMap == "Max")
          value = axisLengthX *2; // times 2 to capture the full length of the cell (not 1/2 or the radius of it)
        else if(displayHeatMap == "Mid")
          value = axisLengthY *2;
        else if(displayHeatMap == "Min")
          value = axisLengthZ *2;
        else if(displayHeatMap == "Elongation (Max/Mid)")
          value = axisLengthX/axisLengthY;
        else if(displayHeatMap == "Max/Min")
          value = axisLengthX/axisLengthZ;
        else if(displayHeatMap == "Flatness (Mid/Min)")
          value = axisLengthY/axisLengthZ;
        else if(displayHeatMap == "Product")
          value = axisLengthX * axisLengthY * axisLengthZ *8;
        else if(displayHeatMap == "Shape Anisotropy")
          value = (axisLengthX-0.5*axisLengthY-0.5*axisLengthZ)/(axisLengthX+axisLengthY+axisLengthZ);
        else if(displayHeatMap == "Custom X")
          value = axisLengthCustomX *2;
        else if(displayHeatMap == "Custom Y")
          value = axisLengthCustomY *2;
        else if(displayHeatMap == "Custom Z")
          value = axisLengthCustomZ *2;

        labelHeatMap[cell] = value;
        vecValues.push_back(value);
      }

      if(attrMaps){

        shape3dProduct[cell] = axisLengthX * axisLengthY * axisLengthZ *8;
        shape3dMax[cell] = axisLengthX *2;
        shape3dMid[cell] = axisLengthY * 2;
        shape3dMin[cell] = axisLengthZ *2;
        shape3dMaxMid[cell] = axisLengthX/axisLengthY;
        shape3dMaxMin[cell] = axisLengthX/axisLengthZ;
        shape3dMidMin[cell] = axisLengthY/axisLengthZ;
        shape3dAniso[cell] = (axisLengthX-0.5*axisLengthY-0.5*axisLengthZ)/(axisLengthX+axisLengthY+axisLengthZ);
        shape3dCustomX[cell] = axisLengthCustomX *2;
        shape3dCustomY[cell] = axisLengthCustomY *2;
        shape3dCustomZ[cell] = axisLengthCustomZ *2;

      }
    }

    // Find bounds for heatmap scale.
    float lowBound = 0, highBound = 1;
    // Manual scaling by user input
    if(scaleHeatMap == "Manual"){
      lowBound  = rangeHeat.x();
      highBound = rangeHeat.y();
    }
    // Automatic scaling
    if(vecValues.size() != 0){
      sort(vecValues.begin(), vecValues.end());
      // Default scaling: take min and max values
      if(scaleHeatMap == "None"){
        lowBound = vecValues.front();
        highBound = vecValues.back();
      }
      // Auto scale: take low 5 and high 5 percentile
      else if(scaleHeatMap == "Auto"){
        lowBound = *(vecValues.begin()+int(vecValues.size()*.05));
        highBound = *(vecValues.begin()+int(vecValues.size()*.95));
      }
    }
    // Round the bounds values of heat map, to make sure we have the same bound in case
    // of loading PCA on both meshes.
    lowBound = floor(lowBound * 1000) / 1000;
    highBound = ceil(highBound * 1000) / 1000;

    mesh->heatMapBounds() = Point2f(lowBound, highBound);
    mesh->heatMapUnit() = QString("a.u.");
    mesh->setShowLabel("Label Heat");
    mesh->updateTriangles();

    if(displayAxis != "None")
      mesh->setShowAxis("Cell Axis");

    return true;
  }
  REGISTER_PROCESS(DisplayPCA3D);

}
