//
// This file is part of 3DCellAtlas.
// Copyright (C) 2015-2016 George W. Bassel and collaborators.
//
// If you use 3DCellAtlas in your work, please cite:
//   http://dx.doi.org/10.1105/tpc.15.00175
//
// 3DCellAtlas is an AddOn for MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2012 Richard S. Smith and collaborators.
//
// 3DCellAtlas and MorphoGraphX are free software, and are licensed under under the terms of the
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
//
#include <QFileDialog>

#include <CellAtlas.hpp>
#include <Progress.hpp>
#include <Geometry.hpp>

#include <ImageData.hpp>

#include <GraphUtils.hpp>

#include <math.h>
#include <limits>
#include <queue>

#include <Triangle.hpp>

#include <MeshProcessHeatMap.hpp>

#include <MeshProcessPDG.hpp> //jInfo

#include <MeshProcessMeasures3D.hpp>
#include <MeshProcessLineage.hpp>

#include <CellAtlasUtils.hpp>

using namespace std;

namespace mgx
{
  // this process is similar to analyze cells 3d, however a 2D bezier is used
  void writeCellAtlasMeasures(Mesh *m)
  {

  CellAtlasAttr *cellAtlasAttr;
  cellAtlasAttr = &m->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

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

  AttrMap<int, double>& volume = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");
  volume.clear();
  AttrMap<int, double>& wallArea = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");
  wallArea.clear();

  AttrMap<int, double>& lon = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Longitudinal");
  lon.clear();
  AttrMap<int, double>& circ = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Circumferential");
  circ.clear();
  AttrMap<int, double>& rad = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Radial");
  rad.clear();

  AttrMap<int, double>& radAbs = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Radial Abs");
  radAbs.clear();

  AttrMap<int, Point3d>& lonDir = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirLongitudinal");
  lonDir.clear();
  AttrMap<int, Point3d>& circDir = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirCircumferential");
  circDir.clear();
  AttrMap<int, Point3d>& radDir = m->attributes().attrMap<int, Point3d>("Measure Label Vector DirRadial");
  radDir.clear();

  AttrMap<int, double>& lonS = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Longitudinal");
  lonS.clear();
  AttrMap<int, double>& circS = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Circumferential");
  circS.clear();
  AttrMap<int, double>& radS = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Radial");
  radS.clear();

  AttrMap<int, double>& outArea = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area");
  outArea.clear();

  AttrMap<int, double>& lon2 = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Longitudinal Coord Abs");
  lon2.clear();
  AttrMap<int, double>& lon3 = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Longitudinal Coord Min");
  lon3.clear();
  AttrMap<int, double>& lon4 = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Longitudinal Coord Max");
  lon4.clear();

  AttrMap<int, double>& assocCC = m->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Associated Cortical Cell");
  assocCC.clear();

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

  for(uint i = 0; i<(*cellAtlasAttr).size(); i++){
    CellAtlasData data = (*cellAtlasAttr)[i];
    int label = data.cellLabel;
    volume[label] = data.volume;
    wallArea[label] = data.wallArea;

    centroids[label] = data.centroid;

    lon[label] = data.coordLong;
    circ[label] = data.coordCirc;
    rad[label] = data.coordRad;

    lonDir[label] = data.dirLong;
    circDir[label] = data.dirRad;
    radDir[label] = data.dirCirc;

    lonS[label] = data.sizeLong;
    circS[label] = data.sizeCirc;
    radS[label] = data.sizeRad;

    outArea[label] = data.outsideWallArea;

    lon2[label] = data.coordLongUM;
    lon3[label] = data.coordLongMinUM;
    lon4[label] = data.coordLongMaxUM;

    assocCC[label] = data.associatedCorticalCell;
    radAbs[label] = data.coordRadAbs;

    customDirs[label].ev1() = Point3f(data.dirLong);
    customDirs[label].ev2() = Point3f(data.dirRad);
    customDirs[label].evals() = Point3f(data.dirCirc);
  }

  }

// vPos - v position of the bezier grid that is being used
// vRefPos - v position of the refined bezier grid that is being used
void discretizeBezierLineRef(CuttingSurface* cutSurf, double vPos, double vRefPos, Matrix4d rotMat, int bezMapSize,
  std::map<int, Point3d>& bMap, std::map<int, Point3d>& diffbMap, std::map<int, Point3d>& refbMap){

    // create bezier vector and its derivative
    double bezStart = 0.0;
    double bezEnd = 1.0;
    int dataPointsBezier = bezMapSize; // number of intermediate points
    double stepSize = (double)(bezEnd-bezStart)/(double)(dataPointsBezier-1);

    int fineInterpolation = 10000;
    double stepSizeFine = (double)(bezEnd-bezStart)/(double)(fineInterpolation-1);

    std::vector<Point3d> bMapFine, bMapRefFine;

    for(int i=0; i<fineInterpolation; i++){
      double u = bezStart+i*stepSizeFine;
      Point3d p(cutSurf->evalCoord(u,vPos));
      Point3d pRef(cutSurf->evalCoord(u,vRefPos));
      p = multMatrix4Point3(rotMat,p); // correct rotations
      pRef = multMatrix4Point3(rotMat,pRef) - p; // correct rotations
      bMapFine.push_back(p);
      pRef/=norm(pRef);
      bMapRefFine.push_back(pRef);
    }

    double totDis = 0;
    for(int i=1; i<fineInterpolation; i++){
      totDis += norm(bMapFine[i] - bMapFine[i-1]);
    }

    double dis = 0;
    int idxCounter = 1;

    bMap[0] = bMapFine[0];
    bMap[bezMapSize-1] = bMapFine[fineInterpolation-1];

    diffbMap[0] = bMapFine[1] - bMapFine[0];
    diffbMap[dataPointsBezier-1] = bMapFine[fineInterpolation-1] - bMapFine[fineInterpolation-2];

    diffbMap[0]/=norm(diffbMap[0]);
    diffbMap[dataPointsBezier-1]/=norm(diffbMap[dataPointsBezier-1]);

    refbMap[0] = bMapRefFine[0];
    refbMap[bezMapSize-1] = bMapRefFine[fineInterpolation-1];

//cout << "bez " << bMap[idxCounter] << "/" << diffbMap[idxCounter] << endl;

    for(std::vector<Point3d>::size_type i=1; i<bMapFine.size()-1; i++){
      dis += norm(bMapFine[i] - bMapFine[i-1])/totDis;

      if(dis>idxCounter*stepSize){
        bMap[idxCounter] = bMapFine[i];
        diffbMap[idxCounter] = bMapFine[i+1] - bMapFine[i-1];
        diffbMap[idxCounter]/=norm(diffbMap[idxCounter]);
        refbMap[idxCounter] = bMapRefFine[i];

        // not sure what this is doing!
        double sN;
        Point3d radOrth;
        planeLineIntersect(bMap[idxCounter], diffbMap[idxCounter], bMapRefFine[i], bMapRefFine[i]-diffbMap[idxCounter], sN, radOrth);

        refbMap[idxCounter] = radOrth - bMap[idxCounter];
        refbMap[idxCounter]/=norm(refbMap[idxCounter]);

            //    cout << "bez " << idxCounter << "/" << dis << "/" << bMap[idxCounter] << "/" << diffbMap[idxCounter] << endl;
        idxCounter++;

      }
    }
}

  // directizes the bezier line into bezMapSize points
  // -> first a very fine interpolation is calculated with 10k points
  //    then the final interpolation is derived from it
  // gives also the differential
  void discretizeBezierLine(CuttingSurface* cutSurf, double vPos, Matrix4d rotMat, int bezMapSize,
    std::map<int, Point3d>& bMap, std::map<int, Point3d>& diffbMap){

    std::map<int, Point3d> refbMap;
    discretizeBezierLineRef(cutSurf, vPos, 0.0, rotMat, bezMapSize, bMap, diffbMap, refbMap);
}

// process to analyze the cell properties
 bool AnalyzeCells::run(const Stack *s1, const Stack *s2, Mesh *m1, Mesh *m2, double minVolume, QString mode, bool openmpMode)
  {
    progressStart("Analyze Cells 3D - Preparation", 0);

    if(mode != "Layer" and mode != "Root" and mode != "Cartesian") return setErrorMessage("Invalid mode!");

    CellAtlasAttr *cellAtlasAttr;
    CellAtlasConfigAttr *cellAtlasConfigAttr;

    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    const vvGraph& segmentedMesh = m1->graph();
    const vvGraph& surfaceMesh = m2->graph();

    VVGraphVec cellVertex;
    VtxIntMap vertexCell;

    int selectedLabel1 = 0;

    // get selected cells
    progressStart("Analyze Cells 3D - Find Selected Cell", 0);

    selectedLabel1 = findSelectedLabel(segmentedMesh);

    cout << "Selected Cell: " << selectedLabel1 << endl;
    if(selectedLabel1 == 0){
      setErrorMessage("Warning: no first cell selected");
      return false;
    }

    // prepare data structures
    //////////////////////////

    progressStart("Analyze Cells 3D - Interpolate Bezier", 0);
    // Bezier
    CuttingSurface* cutSurf = cuttingSurface();

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

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

    // create bezier vector and its derivative
    std::map<int, Point3d> bMap,diffbMap;
    int dataPointsBezier = 200;
    //discretizeBezierLine(cutSurf, 0.0, mGLTot, dataPointsBezier, bMap, diffbMap);

    Bezier b = cutSurf->bezier();

    std::vector<Point3d> bMap2,diffbMap2;
    double totLength;
    //std::cout << "bez " << std::endl;
    b.discretizeLineEqualWeight(dataPointsBezier, 0.0, mGLTot, bMap2, diffbMap2, totLength);
    //std::cout << "bez " << bMap.size() << "/" << bMap2.size() << std::endl;
    for(int i=0; i< dataPointsBezier; i++){
    //  std::cout << "bez " << i << "/" << bMap[i] << "/" << bMap2[i] << "//" << diffbMap[i] << "/" << diffbMap2[i] << std::endl;
    //  if(i>0) std::cout << "bez " << i << "/" << norm(bMap[i]-bMap[i-1]) << "/" << norm(bMap2[i]-bMap2[i-1]) << std::endl;
      bMap[i] = bMap2[i];
      diffbMap[i] = diffbMap2[i];
    }

    //// input parameter ////
    /////////////////////////

    // label of first cell -> comes from selection
    // label of last rootcap cell -> comes from 2nd selection

    progressStart("Analyze Cells 3D - Find Cells", 0);
    std::vector<int> uniqueLabels = findAllLabels(segmentedMesh);
    RootCellProcessing rcp;

    rcp.rootData.numCells = uniqueLabels.size();
    rcp.rootData.uniqueLabels = uniqueLabels;

    labelVertexMap lvMap;
    std::vector<vertex> vertexMap;

    // create label - vertex map
    progressStart("Analyze Cells 3D - Create label/vertex map", 0);
    forall(const vertex &v, segmentedMesh) {
      if(!progressAdvance(0)) userCancel();
      lvMap[v->label].push_back(v);
    }

    std::vector<std::vector<Point3d> > bezGrid;
    b.discretizeGrid(dataPointsBezier, mGLTot, bezGrid);

    // analyze all cells
    RootCellAnalyzing rca(cellAtlasAttr,cellAtlasConfigAttr);
    if(!rca.analyzeEmbryo(uniqueLabels, segmentedMesh, surfaceMesh, bezGrid, bMap, diffbMap,
      lvMap, rcp, selectedLabel1, minVolume, mode))
      userCancel();


    //if(!neighborhoodGraphLocal(segmentedMesh, 0.0001, rcp.rootData.cellWallArea, rcp.rootData.wallArea))
    if(!neighborhoodGraphLocal(segmentedMesh, 0.0001, rcp.rootData.cellWallArea, rcp.rootData.wallArea, rcp.rootData.outsideWallArea))
      userCancel();
    progressStart("Analyze Cells 3D - Finalize", 0);


    forall(IntIntDouPair p, rcp.rootData.wallArea){
      rcp.rootData.nrNeighbors[p.first.first]++;
    }

    writeAttrMaps(rcp, cellAtlasAttr, cellAtlasConfigAttr);

    rcp.setDataPointers(cellAtlasAttr, cellAtlasConfigAttr);
    (*cellAtlasConfigAttr)[0].labelFirstCell = selectedLabel1;

    rcp.correctDirections();

    (*cellAtlasConfigAttr)[0].analyzed = true;

    writeCellAtlasMeasures(m1);

    return true;
  }
  REGISTER_PROCESS(AnalyzeCells);

// set a number to a label
void setLabelProperties(QLabel* label, int x, int y, int number, bool reset, QRgb col)
  {
    int xOffset = 10+10;
    int yOffset = 100+7;
    int xSize = 10;
    int ySize = 10;

    int r = qRed(col);
    int g = qGreen(col);
    int b = qBlue(col);

    label->setGeometry(QRect(x+xOffset, y+yOffset, xSize, ySize));
    if(!reset)
    label->setText(QString::number(number+1));
    else
    label->setText("");
    QString s = "QLabel { color : rgb(" + QString::number(r) + ", " + QString::number(g) + ", " + QString::number(b) + "); }";
    label->setStyleSheet(s);

  }

// for returning the right UI label
 QLabel* getLabel(Ui_CellAtlasDialog& ui, int label)
  {
    if(label == 0){
      return ui.labelCluster1;
    } else if(label == 1){
      return ui.labelCluster2;
    } else if(label == 2){
      return ui.labelCluster3;
    } else if(label == 3){
      return ui.labelCluster4;
    } else if(label == 4){
      return ui.labelCluster5;
    } else if(label == 5){
      return ui.labelCluster6;
    } else if(label == 6){
      return ui.labelCluster7;
    } else if(label == 7){
      return ui.labelCluster8;
    } else if(label == 8){
      return ui.labelCluster9;
    } else if(label == 9){
      return ui.labelCluster10;
    } else if(label == 10){
      return ui.labelCluster11;
    } else if(label == 11){
      return ui.labelCluster12;
    } else if(label == 12){
      return ui.labelCluster13;
    } else if(label == 13){
      return ui.labelCluster14;
    } else if(label == 14){
      return ui.labelCluster15;
    } else if(label == 15){
      return ui.labelCluster16;
    } else if(label == 16){
      return ui.labelCluster17;
    } else if(label == 17){
      return ui.labelCluster18;
    } else if(label == 18){
      return ui.labelCluster19;
    }
    return ui.labelCluster20;
  }

// find the label that has to be changed
// requires ui
void drawNumber(Ui_CellAtlasDialog& ui, int x, int y, int label, int number)
  {
  QRgb col = qRgba(120,0,120,255);
  setLabelProperties(getLabel(ui, label), x, y, number,false,col);

  }

// reset the label number
void resetNumber(Ui_CellAtlasDialog& ui, int label)
  {
  QRgb col = qRgba(120,0,120,255);
  setLabelProperties(getLabel(ui, label), 0, 0, 0,true,col);

  }


// create the qimage with the heatmap
QImage DataToQImage(heatMapDataStructure& body, bool drawPoints)
{
    // init variables
    int size = body.gridSize;
    cellMorphLandscape data = body.heatMap;
    double maxValue = body.high;
    std::map<int,Point2d> maxVec = body.maximaHeatMapAll;

    cellMorphLandscape points = body.points;
    IntIntMap labels = body.maximaHeatMapSelectedLabel;

    // init image and colours
    int factor = 3; // scaling factor
    QImage image( size*factor, size*factor, QImage::Format_ARGB32 );
    QRgb value;
    QRgb white = qRgba(255,255,255,255);
    QRgb grey1 = qRgba(170,170,170,255);
    //QRgb grey2 = qRgba(85,85,85,255);
    //QRgb purple = qRgba(120,0,120,255);

    int yOffset = size*factor-1;

    // draw the heatmap and (if option is checked) small points for the cells
    for(int i = 0; i < size; i++) {
      for(int j = 0; j < size; j++) {
        //std::cout << "i/j/data " << i << "/" << j << "/" << data[i][j] << std::endl;
        for(int k = 0; k < factor; k++) {
          for(int l = 0; l < factor; l++) {
            value = calcRGB(maxValue,data[i][j]);
            image.setPixel(i*factor+k, yOffset-(j*factor+l), value);
          }
        }
        if(drawPoints and points[i][j] > 0){
          image.setPixel(i*factor+factor/2, yOffset-j*factor-factor/2, grey1);
        }
      }
    }

    // draw a small cross on all heatmap maxima
    int maxSize = maxVec.size();
    for(int i = 0; i < maxSize; i++) {
      Point2d currentPoint = maxVec[i];
      drawCross(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, 2, white);
    }

    return image;
}


// create the qimage with the cells
QImage DataToQImageNew(heatMapDataStructure& body)
{

// init variables
    int size = body.gridSize;
    cellMorphLandscape data = body.heatMap;
    //double maxValue = body.high;
    std::map<int,Point2d> maxVec = body.maximaHeatMapAll;

    cellMorphLandscape points = body.points;
    std::map<Point2d,int> pointsLabel = body.pointsLabel;
    IntIntMap labels = body.maximaHeatMapSelectedLabel;

    // init image and colours
    int factor = 3; // scaling factor
    QImage image( size*factor, size*factor, QImage::Format_ARGB32 );
    //QRgb value;
    //QRgb white = qRgba(255,255,255,255);
    //QRgb grey1 = qRgba(170,170,170,255);
    QRgb black = qRgba(0,0,0,255);
    //QRgb purple = qRgba(120,0,120,255);

    int yOffset = size*factor-1;

    // draw the heatmap and (if option is checked) small points for the cells
    for(int i = 0; i < size; i++) {
      for(int j = 0; j < size; j++) {
        //std::cout << "i/j/data " << i << "/" << j << "/" << data[i][j] << std::endl;
        for(int k = 0; k < factor; k++) {
          for(int l = 0; l < factor; l++) {
            ///value = calcRGB(maxValue,data[i][j]);
            image.setPixel(i*factor+k, yOffset-(j*factor+l), black);
          }
        }

        QRgb col = getColorFromLabel(pointsLabel[Point2d(i,j)]);

        if(points[i][j] > 0){
          image.setPixel(i*factor+factor/2, yOffset-j*factor-factor/2, col);
          image.setPixel(i*factor+factor/2+1, yOffset-j*factor-factor/2, col);
          image.setPixel(i*factor+factor/2, yOffset-j*factor-factor/2+1, col);
          image.setPixel(i*factor+factor/2+1, yOffset-j*factor-factor/2+1, col);
        }
        if(points[i][j] > 4){
          image.setPixel(i*factor+factor/2+1, yOffset-j*factor-factor/2-1, col);
          image.setPixel(i*factor+factor/2-1, yOffset-j*factor-factor/2+1, col);
          image.setPixel(i*factor+factor/2-1, yOffset-j*factor-factor/2-1, col);
          image.setPixel(i*factor+factor/2, yOffset-j*factor-factor/2-1, col);
          image.setPixel(i*factor+factor/2-1, yOffset-j*factor-factor/2, col);

        }
      }
    }

    // draw a small cross on all heatmap maxima
    //int maxSize = maxVec.size();
    //for(int i = 0; i < maxSize; i++) {
    //  Point2d currentPoint = maxVec[i];
    //  drawCross(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, 2, white);
    //}

    return image;

}

// checks whether there exists a maximum in maxVec close to the pos coordinates
// writes the result in maxIdx
bool chosenMaxExistNearby(std::map<int,Point2d>& maxVec, Point2d pos, int& maxIdx)
  {
    double minDis = 1E20;
    int maxVecSize = maxVec.size();

    for(int i = 0; i < maxVecSize; i++) {
      double dis = norm(pos - maxVec[i]);
      if(dis < minDis){
        minDis = dis;
        maxIdx = i;
      }
    }

    if(minDis < 10){
      return true;
    } else {
      maxIdx = -1;
      return false;
    }
  }

// generate the heatmap image
void CellAtlas::setImage()
  {

    QImage image;
    std::map<int,Point2d> maxVec2;
    int factor = 3;
    QRgb purple = qRgba(120,0,120,255);
    IntIntMap labels;// = body.maximaHeatMapSelectedLabel;

    int size;/// = body.gridSize;

    // check which root part is active
    if(ui.rootPartCmb->currentIndex() == 0){
      // generate image
      ui.spinBoxSigma->setValue((*cellAtlasConfigAttr)[0].sigma1);
      ui.xDimCmb->setCurrentIndex((*cellAtlasConfigAttr)[0].part1ChosenX-1);
      ui.yDimCmb->setCurrentIndex((*cellAtlasConfigAttr)[0].part1ChosenY-1);
      image = DataToQImage(rcp.mainBody, ui.showCellsCheckBox->isChecked() );
      maxVec2 = rcp.mainBody.maximaHeatMapSelected;
      size = rcp.mainBody.gridSize;
      labels = rcp.mainBody.maximaHeatMapSelectedLabel;
    } else {
      // generate image
      ui.spinBoxSigma->setValue((*cellAtlasConfigAttr)[0].sigma2);
      ui.xDimCmb->setCurrentIndex((*cellAtlasConfigAttr)[0].part2ChosenX-1);
      ui.yDimCmb->setCurrentIndex((*cellAtlasConfigAttr)[0].part2ChosenY-1);
      image = DataToQImage(rcp.radicle, ui.showCellsCheckBox->isChecked() );
      maxVec2 = rcp.radicle.maximaHeatMapSelected;
      size = rcp.radicle.gridSize;
      labels = rcp.radicle.maximaHeatMapSelectedLabel;
    }
    int yOffset = size*factor-1;
    // reset the numbers
    for(int i = 0; i < 20; i++) {
      resetNumber(ui, i);
    }

    // draw the user defined cluster center with circle and big cross
    int maxSize = maxVec2.size();
    for(int i = 0; i < maxSize; i++) {
      Point2d currentPoint = maxVec2[i];
      drawCross(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, 4, purple);
      drawCircle(image, currentPoint[0]*factor+factor/2, yOffset-currentPoint[1]*factor-factor/2, size*factor, purple);
      if(i<20)
        drawNumber(ui, currentPoint[0]*factor+factor/2-2, yOffset-currentPoint[1]*factor-factor/2+2, i, labels[i]);
    }

    // set image to image area
    ui.imageHolder->setPixmap(QPixmap::fromImage(image));
    ui.imageHolder->setFocus();
  }

// upon double click: delete if there is a cluster nearby
void CellAtlas::setDeletePosition(const QPoint& p)
  {
    //cout << "Mouse double click" << p.x() << "   " << p.y() << endl;
    Point2d currentPoint;
    currentPoint[0] = p.x()/3.0;
    currentPoint[1] = (160.0*3.0-1-p.y())/3.0;
    int idx;
    if(ui.rootPartCmb->currentIndex() == 0){ // main
      // check whether there already is a maximum at the clicked position
      if(chosenMaxExistNearby(rcp.mainBody.maximaHeatMapSelected, currentPoint, idx)){
        // yes : delete this maximum
        rcp.delIdx(rcp.mainBody, idx);
      }
    } else { // radicle
      if(chosenMaxExistNearby(rcp.radicle.maximaHeatMapSelected, currentPoint, idx)){
        // yes : delete this maximum
        rcp.delIdx(rcp.radicle, idx);
      }
    }

    setImage();
  }

// upon mouse click: if a maximum is nearby: activate the maximum (for drag and drop)
// if no maximum is nearby create a new one
void CellAtlas::setPosition(const QPoint& p)
  {
    Point2d currentPoint;
    currentPoint[0] = p.x()/3.0;
    currentPoint[1] = (160.0*3.0-1-p.y())/3.0;
    if(ui.rootPartCmb->currentIndex() == 0){ // main
      int maxVecSize = rcp.mainBody.maximaHeatMapSelected.size();
      int idx;
      // check whether there already is a maximum at the clicked position
      if(chosenMaxExistNearby(rcp.mainBody.maximaHeatMapSelected, currentPoint, idx)){
      // yes : activate this maximum
        rcp.mainBody.activatedMaximum = idx;
      } else {
      // no : create a new maximum
        rcp.mainBody.maximaHeatMapSelected[maxVecSize] = currentPoint;
        rcp.mainBody.maximaHeatMapSelectedLabel[maxVecSize] = 0;
        rcp.mainBody.activatedMaximum = -1;
        setImage();
      }

    } else { // radicle

      int maxVecSize = rcp.radicle.maximaHeatMapSelected.size();
      int idx;
      // check whether there already is a maximum at the clicked position
      if(chosenMaxExistNearby(rcp.radicle.maximaHeatMapSelected, currentPoint, idx)){
      // yes : activate this maximum
        rcp.radicle.activatedMaximum = idx;
      } else {
      // no : create a new maximum
        rcp.radicle.maximaHeatMapSelected[maxVecSize] = currentPoint;
        rcp.radicle.maximaHeatMapSelectedLabel[maxVecSize] = 0;
        rcp.radicle.activatedMaximum = -1;
        setImage();
      }

    }

  }

// upon mouse release: check whether there is an activated maximum (=nearby maximum) and move it to the new position
void CellAtlas::setReleasePosition(const QPoint& p)
  {
    Point2d currentPoint;
    currentPoint[0] = p.x()/3.0;
    currentPoint[1] = (160.0*3.0-1-p.y())/3.0;


    if(ui.rootPartCmb->currentIndex() == 0){ // main
      if(rcp.mainBody.activatedMaximum != -1){    // check if there is an activated maximum
        //yes : move it to the new position
        rcp.mainBody.maximaHeatMapSelected[rcp.mainBody.activatedMaximum] = currentPoint;
      }
    } else { // rootcap
      if(rcp.radicle.activatedMaximum != -1){    // check if there is an activated maximum
        //yes : move it to the new position
        rcp.radicle.maximaHeatMapSelected[rcp.radicle.activatedMaximum] = currentPoint;
      }
    }
    setImage();

    rcp.mainBody.activatedMaximum = -1;
    rcp.radicle.activatedMaximum = -1;

  }

// automatically assign the highest maxima
void CellAtlas::setAutoCluster()
  {
    if(ui.rootPartCmb->currentIndex() == 0){ // main
      rcp.setAutoCluster(ui.spinBoxAutoClustering->value(), true);
    } else {
      rcp.setAutoCluster(ui.spinBoxAutoClustering->value(), false);
    }
    setImage();
    ui.imageHolder->setFocus();
  }

// upon keypress: check whether there is a cluster nearby and if yes: assign the label to it
void CellAtlas::setClusterLabel(QString label){

    if(label == "1" or label == "2" or label == "3" or label == "4" or label == "5" or label == "6" or label == "7" or label == "8" or label == "9"){
      int labelNr = label.toInt() - 1;
      int idx;
      std::map<int,Point2d> maxVec;
      if(ui.rootPartCmb->currentIndex() == 0){ // main
        //cout << " ggg" << endl;
        maxVec = rcp.mainBody.maximaHeatMapSelected;
        if(chosenMaxExistNearby(rcp.mainBody.maximaHeatMapSelected, mousePos, idx)){
          // assign label to cluster
          rcp.mainBody.maximaHeatMapSelectedLabel[idx] = labelNr;
        }
      } else { // rootcap
        maxVec = rcp.radicle.maximaHeatMapSelected;
        if(chosenMaxExistNearby(rcp.radicle.maximaHeatMapSelected, mousePos, idx)){
          // assign label to cluster
          rcp.radicle.maximaHeatMapSelectedLabel[idx] = labelNr;
        }
      }
    } else {
      cout << "Enter a number between 1 and 9" << endl;
    }

    setImage();

  }

// update the position of the mouse within the heatmap window
void CellAtlas::setMousePosition(const QPoint& p){

    Point2d currentPoint;
    currentPoint[0] = p.x()/3.0;
    currentPoint[1] = (160.0*3.0-1-p.y())/3.0;
    mousePos = currentPoint;

  }

// change the Y coord of the heatmap and redraw with the new data
void CellAtlas::changeHeatmapX(QString option)
  {
    if(ui.rootPartCmb->currentIndex() == 0){ // main
        rcp.setHeatmap(rcp.mainBody, option, (*cellAtlasConfigAttr)[0].part1ChosenX, true);
    } else { // root
      //rcp.setHeatmapX(rcp.radicle, stringX);
        rcp.setHeatmap(rcp.radicle, option, (*cellAtlasConfigAttr)[0].part2ChosenX, true);
    }

    rcp.geometryBasedLabelling();
    setImage();
  }


// change the Y coord of the heatmap and redraw with the new data
void CellAtlas::changeHeatmapY(QString option)
  {
    if(ui.rootPartCmb->currentIndex() == 0){ // main
      rcp.setHeatmap(rcp.mainBody, option, (*cellAtlasConfigAttr)[0].part1ChosenY, false);
    } else { // root
      rcp.setHeatmap(rcp.radicle, option, (*cellAtlasConfigAttr)[0].part2ChosenY, false);
    }

    rcp.geometryBasedLabelling();
    setImage();
  }

void updateParentLabels(Mesh *m1)
  {
    // update parent labels with the generated ones
    m1->parents().clear();

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    int numCells = (*cellAtlasAttr).size();

    for(int i=0; i<numCells; i++){
      int currentLabel = (*cellAtlasAttr)[i].cellLabel;
      m1->parents()[currentLabel] = (*cellAtlasAttr)[i].parentLabel;
    }

    m1->setUseParents(true);
    m1->updateTriangles();
  }

// changes the sigma parameter of the heatmap and redraws the heatmap
void CellAtlas::changeSigma(double sigma)
  {

    if(ui.rootPartCmb->currentIndex() == 0){ // main
      (*cellAtlasConfigAttr)[0].sigma1 = sigma;

    } else { // radicle
      (*cellAtlasConfigAttr)[0].sigma2 = sigma;
    }

    rcp.geometryBasedLabelling();
    setImage();

  }

// switch the heatmap view from main body to radicle or back
void CellAtlas::changeRootPart()
  {

    if(ui.rootPartCmb->currentIndex() == 0){ // main
      ui.spinBoxSigma->setValue((*cellAtlasConfigAttr)[0].sigma1);
    } else { // radicle
      ui.spinBoxSigma->setValue((*cellAtlasConfigAttr)[0].sigma2);
    }

    setImage();
  }

void CellAtlas::resetCluster()
  {
    // reset all preselected clusters

    if(ui.rootPartCmb->currentIndex() == 0){ // main
      rcp.resetHeat(rcp.mainBody);

    } else { // radicle
      rcp.resetHeat(rcp.radicle);
    }
    rcp.geometryBasedLabelling();
    setImage();
  }

void CellAtlas::setPreCluster()
  {
// take all currently selected clusters, look for their nearest maximum, register them as preselected and remove them from the heatmap

    if(ui.rootPartCmb->currentIndex() == 0){ // main
      rcp.preselectCells(rcp.mainBody);

    } else { // radicle
      rcp.preselectCells(rcp.radicle);
    }
    rcp.geometryBasedLabelling();
    setImage();

  }

  bool CellAtlas::initialize(QWidget* parent)
  {
    progressStart("Assign Cell Types", 0);

    Mesh *m1 = mesh(0);
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();


    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      std::cout << "out " << std::endl;
      return true;
    }
    //bool hasRootCap = false;
    bool hasRootCap = stringToBool(parm("Has multiple segments"));

    if(!(*cellAtlasConfigAttr)[0].analyzed){
      // do nothing and return error message in operator
    } else {

      const vvGraph& segmentedMesh = m1->graph();

      int selectedLabel1 = 0;
      int selectedLabel2 = 0;

      findTwoSelectedLabels(segmentedMesh,selectedLabel1,selectedLabel2);
      cout << "Selected Cells: " << selectedLabel1 << "  " << selectedLabel2 << endl;
      if(selectedLabel1 == 0){
        // do nothing and return error message in operator
        setParm("Has multiple segments", "error");
      } else {

        dlg = new QDialog(parent);
        ui.setupUi(dlg);
        this->dlg = dlg;

        if(!hasRootCap){
          ui.rootPartCmb->setEnabled(false);
        } else {
          ui.rootPartCmb->setEnabled(true);
        }

        QStringList cellAtlasMeasures;
        cellAtlasMeasures << "Radial Coordinate" << "Circumferential Coordinate" << "Longitudinal Coordinate"
        << "Radial Cell Length" << "Circumferential Cell Length" << "Longitudinal Cell Length"
        << "Cell Volume" << "Cell Wall Area";

        ui.xDimCmb->addItems(cellAtlasMeasures);
        ui.yDimCmb->addItems(cellAtlasMeasures);


        connect(ui.imageHolder, SIGNAL(mousePressed( const QPoint& )), this, SLOT(setPosition(const QPoint&)));
        connect(ui.imageHolder, SIGNAL(mouseReleased( const QPoint& )), this, SLOT(setReleasePosition(const QPoint&)));
        connect(ui.buttonAutoClustering, SIGNAL(released()), this, SLOT(setAutoCluster()));
        connect(ui.preselectButton, SIGNAL(released()), this, SLOT(setPreCluster()));
        connect(ui.imageHolder, SIGNAL(mouseDouble( const QPoint& )), this, SLOT(setDeletePosition(const QPoint&)));
        connect(ui.showCellsCheckBox, SIGNAL(stateChanged(int)), this, SLOT(setImage()));
        connect(ui.rootPartCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeRootPart()));
        connect(ui.imageHolder, SIGNAL(keyPressed( const QString& )), this, SLOT(setClusterLabel(QString)));
        connect(ui.imageHolder, SIGNAL(mouseMove( const QPoint& )), this, SLOT(setMousePosition(const QPoint&)));
        connect(ui.yDimCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeHeatmapY(QString)));
        connect(ui.xDimCmb, SIGNAL(currentIndexChanged(QString)), this, SLOT(changeHeatmapX(QString)));
        connect(ui.spinBoxSigma, SIGNAL(valueChanged(double)), this, SLOT(changeSigma(double)));
        connect(ui.buttonReset, SIGNAL(released()), this, SLOT(resetCluster()));

        rcp.setDataPointers(cellAtlasAttr, cellAtlasConfigAttr);

        rcp.setParameter(hasRootCap, selectedLabel1, selectedLabel2);
        rcp.assignRootRegions();
        rcp.geometryBasedLabelling();
        setImage();

        cout << " main cells  " << rcp.mainBody.cellCounter << " root cells  " << rcp.radicle.cellCounter << " col cells " << rcp.columella.cellCounter << endl;

        if(dlg->exec() == QDialog::Accepted)
        {
          //relabel
          rcp.geometryBasedLabelling();
          // update Labels
          updateParentLabels(m1);// = true;
          (*cellAtlasConfigAttr)[0].cellTypesAssigned = true;
        } else {
          // aborted: do nothing
        }
      }
    }
    return true;
  }

 bool CellAtlas::run(QString rootCap)
  {
    Mesh* m1 = mesh(0);

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or !(*cellAtlasConfigAttr)[0].analyzed or (*cellAtlasConfigAttr).size() == 0){
      return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    if(rootCap == "error"){
      return setErrorMessage("No cells selected. Please select a first cell and if necessary a last root cap cell.");
    }

    writeCellAtlasMeasures(m1);

    return true;
  }
  REGISTER_PROCESS(CellAtlas);

// selects all cells that are labelled as bad (=too small, no centroid or no intersect found)
 bool SelectBadCells::run(Mesh *m1)
  {

    //get cell atlas attribute maps
    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0) return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");

    forall(const vertex& v, m1->graph())
    {
      int idx = (*cellAtlasConfigAttr)[0].labelIdxMap[v->label];
      if((*cellAtlasAttr)[idx].badCell) v->selected = true;
      else v->selected = false;
    }

    m1->correctSelection(true);
    m1->updateSelection();

    return true;
  }
  REGISTER_PROCESS(SelectBadCells);

// selects all cells that are labelled as bad (=too small, no centroid or no intersect found)
 bool CellAtlasFindCells::run(Mesh *m1, double vol, double sizes, double volratio, double outsideArea, QString logicalCon)
  {

    //get cell atlas attribute maps
    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0) return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");

    vvGraph& S = m1->graph();

    std::set<int> validCells;

    for(int i = 0; i<numCells; i++){
      int currentLabel = (*cellAtlasAttr)[i].cellLabel;

      double volWallRatio = std::pow((*cellAtlasAttr)[i].volume,1./3.)/std::pow((*cellAtlasAttr)[i].wallArea, 0.5);

      if(logicalCon == "OR"){

        if((*cellAtlasAttr)[i].volume < vol) validCells.insert(currentLabel);

        if((*cellAtlasAttr)[i].sizeLong < sizes) validCells.insert(currentLabel);
        if((*cellAtlasAttr)[i].sizeRad < sizes) validCells.insert(currentLabel);
        if((*cellAtlasAttr)[i].sizeCirc < sizes or volWallRatio < volratio) validCells.insert(currentLabel);
        if((*cellAtlasAttr)[i].outsideWallAreaPercent > outsideArea) validCells.insert(currentLabel);
      } else {
        if((*cellAtlasAttr)[i].volume < vol and (*cellAtlasAttr)[i].sizeLong < sizes and
           (*cellAtlasAttr)[i].sizeRad < sizes and (*cellAtlasAttr)[i].sizeCirc < sizes and
           volWallRatio < volratio and (*cellAtlasAttr)[i].outsideWallAreaPercent < outsideArea)
             validCells.insert(currentLabel);
      }

    }

    forall(const vertex& v, S)
    {
      if(validCells.find(v->label) != validCells.end()) v->selected = true;
      else v->selected = false;
    }

    m1->correctSelection(true);
    m1->updateSelection();

    return true;
  }
  REGISTER_PROCESS(CellAtlasFindCells);

  // file dialogue for save data
  bool SaveCellData::initialize(QWidget* parent)
  {

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


  bool saveNeighborhoodFile(QString filename, std::map<IntInt, double>& sharedWallArea, bool parentLabels, Mesh* m)
  {

    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out3(&file);

    out3 << "label, neighbor, wall area" << endl;

    std::map<IntInt, double> sharedWallAreaToFile;

    if(parentLabels){
      forall(const IntIntDouPair &pr, sharedWallArea){
        IntInt p = std::make_pair(m->parents()[pr.first.first],m->parents()[pr.first.second]);
        sharedWallAreaToFile[p] += pr.second;
      }
    } else {
      sharedWallAreaToFile = sharedWallArea;
    }


    forall(const IntIntDouPair &pr, sharedWallAreaToFile){
      out3 << pr.first.first << "," << pr.first.second << "," << pr.second << endl;
    }

    return true;
  }


  bool saveNeighborhoodFile(QString filename, Mesh *m1)
  {

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return false; // setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    QFile file3(filename);
    if(!file3.open(QIODevice::WriteOnly))
    {
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out3(&file3);

    out3 << "label, neighbor, wall area" << endl;

    //forall(const IntIntDouPair &pr, rcp.rootData.wallArea){
    forall(const IntIntDouPair &pr, (*cellAtlasConfigAttr)[0].wallArea){
      out3 << pr.first.first << "," << pr.first.second << "," << pr.second << endl;
    }

    return true;
  }

  bool saveRootFiles(QString filename, bool extended, Mesh *m1)
  {

    cout << "fileroot   " << filename.toLocal8Bit().constData() << endl;

    // CellAtlasAttr *cellAtlasAttr;
    // cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    // CellAtlasConfigAttr *cellAtlasConfigAttr;
    // cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    // int numCells = (*cellAtlasAttr).size();

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

    AttrMap<int, double>& volume = m1->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");
    AttrMap<int, double>& wallArea = m1->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");

    AttrMap<int, double>& lon = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Longitudinal");
    AttrMap<int, double>& circ = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Circumferential");
    AttrMap<int, double>& rad = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Coord Radial");

    AttrMap<int, double>& lonS = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Longitudinal");
    AttrMap<int, double>& circS = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Circumferential");
    AttrMap<int, double>& radS = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Cell Length Radial");

    AttrMap<int, double>& lon2 = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Longitudinal Coord Abs");
    AttrMap<int, double>& radAbs = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Radial Coord Abs");
    
    AttrMap<int, double>& assocCC = m1->attributes().attrMap<int, double>("Measure Label Double CellAtlas/Associated Cortical Cell");

    if(centroids.empty() or  lon.empty()){
      std::cout << "No cell data present. Please run Analyze Cells first." << std::endl;
      return false;
    }

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

    out << "Cell label,Cell type,Associated Cortical Cell,Arclength,Radius,Longitudinal length,Radial length,Circumferential length,Cell Volume,Cell Wall Area";
    if(extended) out << ",Circumferential Angle,Centroid (x),Centroid (y),Centroid (z),Arclength (um),Radial (um)";
    out << endl;

    forall(auto p, centroids){
      int currentLabel = p.first;
      out << currentLabel << "," << m1->parents()[currentLabel] << "," << assocCC[currentLabel] << "," << lon[currentLabel] 
          << "," << rad[currentLabel] << "," << lonS[currentLabel] << "," << radS[currentLabel] << "," << circS[currentLabel]
          << "," << volume[currentLabel] << "," << wallArea[currentLabel];
      if(extended) out << "," << circ[currentLabel] << "," << centroids[currentLabel].x() << "," << centroids[currentLabel].y() << "," << centroids[currentLabel].z() << "," << lon2[currentLabel] << "," << radAbs[currentLabel] ;
      out << endl;
    }
/*
    // file for selected cluster positions and labels
    filename.remove(".csv", Qt::CaseInsensitive);
    filename += "_config";
    filename += ".csv";

    QFile file2(filename);
    if(!file2.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out2(&file2);

    out2 << "x, y" << endl;

    int lengthMain = rcp.mainBody.maximaHeatMapSelected.size();
    int lengthRoot = rcp.radicle.maximaHeatMapSelected.size();

    out2 << -1*rcp.mainBody.chosenX << "," << rcp.mainBody.chosenY << "," << rcp.rootData.labelFirstCell << endl;

    for(int i = 0; i<lengthMain; i++){
      out2 << rcp.mainBody.maximaHeatMapSelected[i][0] << "," << rcp.mainBody.maximaHeatMapSelected[i][1] << "," << rcp.mainBody.maximaHeatMapSelectedLabel[i] << endl;
    }

    out2 << -1*rcp.radicle.chosenX << "," << rcp.radicle.chosenY << "," << rcp.rootData.labelLastRootCapCell << endl;

    for(int i = 0; i<lengthRoot; i++){
      out2 << rcp.radicle.maximaHeatMapSelected[i][0] << "," << rcp.radicle.maximaHeatMapSelected[i][1] << "," << rcp.radicle.maximaHeatMapSelectedLabel[i] << endl;
    }
*/
    // file for neighborhood clusters
    filename.remove("_config.csv", Qt::CaseInsensitive);
    filename += "_neighborhood";
    filename += ".csv";

    saveNeighborhoodFile(filename, m1);

    // QFile file3(filename);
    // if(!file3.open(QIODevice::WriteOnly))
    // {
    //   //setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
    //   std::cout << "File cannot be opened for writing";
    //   return false;
    // }
    // QTextStream out3(&file3);

    // out3 << "label, neighbor, wall area" << endl;

    // //forall(const IntIntDouPair &pr, rcp.rootData.wallArea){
    // forall(const IntIntDouPair &pr, (*cellAtlasConfigAttr)[0].wallArea){
    //   out3 << pr.first.first << "," << pr.first.second << "," << pr.second << endl;
    // }

    return true;
  }

  bool saveMeristemFile(QString filename, bool extended, IntFloatAttr& LabelIntSig, Mesh *m1)
  {

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return false; // setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    // file for root cell data
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
    {
      std::cout << "File cannot be opened for writing";
      return false;
    }
    QTextStream out(&file);

    out << "Cell label,Cell type,Cell wall area,Volume,Interior signal";
    //if(extended) out << ", Circumferential Angle,Centroid (x), Centroid (y), Centroid (z)";
    out << endl;

    for(int i = 0; i<numCells; i++){
      int currentLabel = (*cellAtlasAttr)[i].cellLabel;// rcp.rootData.uniqueLabels[i];
      //out << currentLabel << "," << m1->parents()[currentLabel] << "," << rcp.rootData.cellWallArea[currentLabel] << "," << rcp.rootData.cellVolumes[currentLabel] << "," << LabelIntSig[currentLabel];
      out << currentLabel << "," << m1->parents()[currentLabel] << "," << (*cellAtlasAttr)[i].wallArea << "," <<(*cellAtlasAttr)[i].volume << "," << LabelIntSig[currentLabel];
      out << endl;
    }

    return true;
  }


// saves data in 3 files
  bool SaveCellData::run(QString filename, QString type, bool extended)
  {
    cout << "file   " << filename.toLocal8Bit().constData() << endl;

    Mesh *m1 = mesh(0);

    if(type == "Root") return saveRootFiles(filename, extended, m1);
    if(type == "Meristem"){
      ComputeHeatMap chm(*this);
      ComputeHeatMap::ImageHeatMaps map;
      chm.getLabelMaps(map,m1);
      IntFloatAttr& LabelIntSig = map.LabelIntSig;
      return saveMeristemFile(filename, extended, LabelIntSig, m1);
    }

    return false;
  }
  REGISTER_PROCESS(SaveCellData);

 // file dialogue for save data
 bool SaveNeighborhoodData::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Filename", filename);
    return true;
  }

// saves neighborhood data
  bool SaveNeighborhoodData::run(Mesh* m, QString filename, NhbdGraphInfo& info1, bool parentLabels)
  {
    //cout << "file   " << filename.toLocal8Bit().constData() << endl;

    return saveNeighborhoodFile(filename, info1.sharedWallArea, parentLabels, m);
  }
  REGISTER_PROCESS(SaveNeighborhoodData);

 // file dialogue for save data
 bool SaveNeighborhoodData2D::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Filename", filename);
    return true;
  }

// saves neighborhood data
  bool SaveNeighborhoodData2D::run(Mesh* m, Mesh* m2, QString filename, bool parentLabels)
  {
    //cout << "file   " << filename.toLocal8Bit().constData() << endl;

    std::map<IntInt, double> neighMap, neighMap2;

    neighborhood2D(m->graph(), neighMap);

    neighborhood2D(m2->graph(), neighMap2);

    std::map<IntInt, double> neighMapP, neighMapP2;

    if(m->useParents()){
      forall(const IntIntDouPair &pr, neighMap){
        IntInt p = std::make_pair(m->parents()[pr.first.first],m->parents()[pr.first.second]);
        neighMapP[p] += pr.second;
      }
    } else {
      neighMapP = neighMap;
    }

    if(m2->useParents()){
      forall(const IntIntDouPair &pr, neighMap2){
        IntInt p = std::make_pair(m2->parents()[pr.first.first],m2->parents()[pr.first.second]);
        neighMapP2[p] += pr.second;
      }
    } else {
      neighMapP2 = neighMap2;
    }

    std::map<IntInt, double> changeMap;

    forall(const IntIntDouPair &pr, neighMapP){
      if(neighMapP2.find(pr.first) == neighMapP2.end()) continue;
      if(pr.second < 0.0001) changeMap[pr.first] = -1;
      changeMap[pr.first] = neighMapP2[pr.first] / pr.second;
    }

    if(parm("Save") == "Active Mesh")
      saveNeighborhoodFile(filename, neighMap, parentLabels, m);
    else if(parm("Save") == "Other Mesh")
      saveNeighborhoodFile(filename, neighMap2, parentLabels, m2);
    else{
      QFile file(filename);
      if(!file.open(QIODevice::WriteOnly))
      {
        std::cout << "File cannot be opened for writing";
        return false;
      }
      QTextStream out3(&file);

      out3 << "label,neighbor,wall area,change" << endl;

      forall(const IntIntDouPair &pr, changeMap){
        out3 << pr.first.first << "," << pr.first.second << "," << neighMapP[pr.first] << "," << pr.second << endl;
      }
    }

    return true;
  }
  REGISTER_PROCESS(SaveNeighborhoodData2D);

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

 // TODO doesnt work for Config and Neighborhood!
 bool LoadCellData::run(QString filename, QString type)
  {
  std::cout << "load" << endl;
    Mesh* m1 = mesh(0);

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    //int numCells = (*cellAtlasAttr).size();

   // if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
   //   return false; // setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
   // }

    RootCellProcessing rcp;
    cout << "file   " << filename.toLocal8Bit().constData() << endl;

    bool ok;
    int counter = 0;
    if(type == "All" or type == "Data"){

    // open file
    QFile file(filename);
    if(!file.open(QIODevice::ReadOnly))
    {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }
    QTextStream ss(&file);
    QString line = ss.readLine();
    QStringList fields = line.split(",");

    // variables to be filled

    std::vector<int> uniqueLabels;
    std::map<int,double> parents, scaledRadialDis, arclengths;
    std::map<int,double> lengthRad, lengthLong, lengthCirc, vol, wall;
    IntIntMap assocCC;


    while(ss.status() == QTextStream::Ok)
    {
      fields = ss.readLine().split(",");
      if(fields[0].isEmpty()) break;

      int label = fields[0].toInt(&ok);
      int parent = fields[1].toInt(&ok);

      //cout << "lll  " << label << " parent" << parent << " long " << cellLengthLong << " rad " << cellLengthRad << " circ " << cellLengthCirc << endl;
      uniqueLabels.push_back(label);
      assocCC[label] = fields[2].toInt(&ok);
      arclengths[label] = fields[3].toDouble(&ok);
      scaledRadialDis[label] = fields[4].toDouble(&ok);
      lengthLong[label] = fields[5].toDouble(&ok);
      lengthRad[label] = fields[6].toDouble(&ok);
      lengthCirc[label] = fields[7].toDouble(&ok);
      vol[label] = fields[8].toDouble(&ok);
      wall[label] = fields[9].toDouble(&ok);
      parents[label] = parent;
      (*cellAtlasAttr)[counter].parentLabel = parent;
      counter++;
    }

    // fill rootData structure
    rcp.rootData.associatedCorticalCell = assocCC;
    rcp.rootData.scaledRadialDis = scaledRadialDis;
    rcp.rootData.arclengths = arclengths;
    rcp.rootData.lengthLong = lengthLong;
    rcp.rootData.lengthRad = lengthRad;
    rcp.rootData.lengthCirc = lengthCirc;
    rcp.rootData.cellWallArea = wall;
    rcp.rootData.uniqueLabels = uniqueLabels;
    rcp.rootData.numCells = counter;
    rcp.rootData.cellVolumes = vol;

    //rcp.cellLabel = parents;

    // file for selected clusters
    filename.remove(".csv", Qt::CaseInsensitive);
    filename += "_config";
    filename += ".csv";

    }

    if(type == "All" or type == "Config"){

cout << "file   " << filename.toLocal8Bit().constData() << endl;
    QFile file2(filename);
    if(!file2.open(QIODevice::ReadOnly))
    {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }

    QTextStream ss2(&file2);
    QString line = ss2.readLine();
    QStringList fields = line.split(",");

    counter = 0;
    std::map<int,Point2d> maxVecRootCap, maxVecMainBody;
    IntIntMap labelRadicle, labelMain;
    int xDimRootCap = 0, xDimMainBody = 0;
    int yDimRootCap = 0, yDimMainBody = 0;
    bool root = false;
    int labelFirst=0, labelRoot=0;

    while(ss2.status() == QTextStream::Ok)
    {
      fields = ss2.readLine().split(",");
      if(fields[0].isEmpty()) break;
      double xCoord = fields[0].toDouble(&ok);
      double yCoord = fields[1].toDouble(&ok);
      int number = fields[2].toInt(&ok);
      if(!root){ // main
        if(xCoord < 0){ // break point
          if(yDimMainBody>0){
            root = true;
            xDimRootCap = -1*xCoord;
            yDimRootCap = yCoord;
            labelRoot = fields[2].toDouble(&ok);
            counter = 0;
          } else {
            xDimMainBody = -1*xCoord;
            yDimMainBody = yCoord;
            labelFirst = fields[2].toDouble(&ok);
          }
        } else {
          maxVecMainBody[counter][0] = xCoord;
          maxVecMainBody[counter][1] = yCoord;
          labelMain[counter] = number;
          counter++;
        }
      } else { // root cap
        maxVecRootCap[counter][0] = xCoord;
        maxVecRootCap[counter][1] = yCoord;
        labelRadicle[counter] = number;
        counter++;
      }
    }

    rcp.mainBody.maximaHeatMapSelected = maxVecMainBody;
    rcp.mainBody.maximaHeatMapSelectedLabel = labelMain;
    rcp.mainBody.chosenX = xDimMainBody;
    rcp.mainBody.chosenY = yDimMainBody;
    rcp.radicle.maximaHeatMapSelected = maxVecRootCap;
    rcp.radicle.maximaHeatMapSelectedLabel = labelRadicle;
    rcp.radicle.chosenX = xDimRootCap;
    rcp.radicle.chosenY = yDimRootCap;

    // file for neighborhood
    filename.remove("_config.csv", Qt::CaseInsensitive);
    filename += "_neighborhood";
    filename += ".csv";


    // selected first and last radicle cell
    forall(const vertex& v, m1->graph())
    {
      if(v->label == labelRoot or v->label == labelFirst) v->selected = true;
      else v->selected = false;
    }


    }

    if(type == "All" or type == "Neighborhood"){

cout << "file   " << filename.toLocal8Bit().constData() << endl;
    QFile file3(filename);
    if(!file3.open(QIODevice::ReadOnly))
    {
      setErrorMessage(QString("File '%1' cannot be opened for reading").arg(filename));
      return false;
    }

    QTextStream ss3(&file3);
    QString line = ss3.readLine();
    QStringList fields = line.split(",");

    std::map<IntInt, double> wallArea;
    while(ss3.status() == QTextStream::Ok){
      fields = ss3.readLine().split(",");
      if(fields[0].isEmpty()) break;
      int neigh1 = fields[0].toInt(&ok);
      int neigh2 = fields[1].toInt(&ok);
      double area = fields[2].toDouble(&ok);

      wallArea[make_pair(neigh1, neigh2)] = area;
      wallArea[make_pair(neigh2, neigh1)] = area;

    }

    rcp.rootData.wallArea = wallArea;
    }


    m1->correctSelection(true);
    m1->updateSelection();

    // set parent labels
    updateParentLabels(m1);



    writeAttrMaps(rcp, cellAtlasAttr, cellAtlasConfigAttr);

    //loaded = true;
    return true;
  }
  REGISTER_PROCESS(LoadCellData);

 bool DisplayCellData::run(Mesh* m, QString choice)
  {

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0) return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");

    double maxValue = -1E20;
    double minValue = 1E20;
    //rcp.calcMinMax(heatMap, minValue, maxValue);

    //float delta = maxValue - minValue;
    m->wallHeat().clear();
    m->labelHeat().clear();
    for(int i = 0; i<numCells; i++){
      int currentLabel = (*cellAtlasAttr)[i].cellLabel;
      double value;
      if(choice == "Longitudinal (Arclengths)"){
        value = (*cellAtlasAttr)[i].coordLong;
      } else if(choice == "Longitudinal (Absolute)"){
        value = (*cellAtlasAttr)[i].coordLongUM;
      } else if(choice == "Radial"){
        value = (*cellAtlasAttr)[i].coordRad;
      } else if(choice == "Radial (Absolute)"){
        value = (*cellAtlasAttr)[i].coordRadAbs;
      } else if(choice == "Longitudinal Cell Length"){
        value = (*cellAtlasAttr)[i].sizeLong;
      } else if(choice == "Radial Cell Length"){
        value = (*cellAtlasAttr)[i].sizeRad;
      } else if(choice == "Circumferential Cell Length"){
        value = (*cellAtlasAttr)[i].sizeCirc;
      } else if(choice == "Volume"){
        value = (*cellAtlasAttr)[i].volume;
      } else if(choice == "Circumferential Angle"){
        value = (*cellAtlasAttr)[i].coordCirc;
      } else if(choice == "Cell Wall Area"){
        value = (*cellAtlasAttr)[i].wallArea;
      } else if(choice == "Nr of Neighbors"){
        value = (*cellAtlasAttr)[i].nrNeighbors;
      } else if(choice == "Outside Wall Area"){
        value = (*cellAtlasAttr)[i].outsideWallArea;
      } else if(choice == "Outside Wall Area percent"){
        value = (*cellAtlasAttr)[i].outsideWallAreaPercent;
      } else if(choice == "Associated Cortical Cell"){
        value = (*cellAtlasAttr)[i].associatedCorticalCell;
      } else if(choice == "Long Min"){
        value = (*cellAtlasAttr)[i].coordLongMinUM;
      } else if(choice == "Long Max"){
        value = (*cellAtlasAttr)[i].coordLongMaxUM;
      } else {
        return setErrorMessage("Wrong Heatmap option.");
      }
      m->labelHeat()[currentLabel] = value;
      if(maxValue < value) maxValue = value;
      if(minValue > value) minValue = value;
    }
    m->heatMapBounds() = Point2f(minValue, maxValue);
    m->setShowLabel("Label Heat");
    m->updateTriangles();
    SETSTATUS(QString("Loaded heat map with values ranging from %1 to %2").arg(minValue).arg(maxValue));

    return true;
  }
  REGISTER_PROCESS(DisplayCellData);


 bool TopologicalCheck::run(QString selection, double threshVol, double threshWallArea, int rootCapLabel, int airBubbleLabel, int vascLabel, int errors)
  {
  bool selectedArea = stringToBool(selection);
  Mesh *m1 = mesh(0);

  CellAtlasAttr *cellAtlasAttr;
  cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

  CellAtlasConfigAttr *cellAtlasConfigAttr;
  cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

  int numCells = (*cellAtlasAttr).size();

  if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
    return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
  }

  // 11 for parent labels 0 to 10
  typedef Vector<11, double> Vec11;
  typedef Vector<11, bool> Vec11Bool;
  typedef Vector<11, Vec11> Mat11;
  typedef Vector<11, Vec11Bool> Mat11Bool;

  Vec11 nrOfParentLabels;
  Vec11 avgVol;
  Mat11Bool valid;
  Mat11 nrOfParentConnections;

  std::set<int> mislabelled;

  intToVec cellConnections;
  intToVec cellWallAreas;

  // create a connection list
  forall(const IntIntDouPair &pr, (*cellAtlasConfigAttr)[0].wallArea/*rcp.rootData.wallArea*/){
    int idx1 = (*cellAtlasConfigAttr)[0].labelIdxMap[pr.first.first];
    int idx2 = (*cellAtlasConfigAttr)[0].labelIdxMap[pr.first.second];
    //if(pr.second > threshWallArea and rcp.rootData.cellVolumes[pr.first.first] > threshVol and rcp.rootData.cellVolumes[pr.first.second] > threshVol){
    if(pr.second > threshWallArea and (*cellAtlasAttr)[idx1].volume > threshVol and (*cellAtlasAttr)[idx2].volume > threshVol){
      cellConnections[pr.first.first].push_back(pr.first.second);
      cellWallAreas[pr.first.first].push_back(pr.second);
    }
  }

  std::map<int, bool> cellSelected;

  // check for selection
  const vvGraph& segmentedMesh = m1->graph();
  forall(const vertex& v, segmentedMesh) {
    if(v->selected or !selectedArea)
      cellSelected[v->label] = true;
  }


  double avgNeigh=0, zeroNeigh=0;
  // check whether a cell has neighbors of a wrong label
  for(int i = 0; i<numCells; i++){
    int currentLabel = (*cellAtlasAttr)[i].cellLabel;//rcp.rootData.uniqueLabels[i];
    if(cellSelected[currentLabel]){
      int currentNeighbors = cellConnections[currentLabel].size();
      avgNeigh+=currentNeighbors;
      if(currentNeighbors==0) zeroNeigh++;
      //mislabelled.insert(currentLabel);
      int currentParent = m1->parents()[currentLabel];
      nrOfParentLabels[currentParent]++;
      avgVol[currentParent]+=(*cellAtlasAttr)[i].volume;//rcp.rootData.cellVolumes[currentLabel];

      for(int j = 0; j<currentNeighbors; j++){
        // fill frequency matrix
        int currentNeighborParent = m1->parents()[cellConnections[currentLabel][j]];
        nrOfParentConnections[currentParent][currentNeighborParent]++;
      }
    }
  }

  // find the two highest connection values and mark them as labelled
  for(int i = 0; i<11; i++){
    double high1=0, high2=0;
    double high1pos=0, high2pos=0;
    if(i!=airBubbleLabel and i!=10 and i!=0){ // ignore airbubbles, columella and unassigned
      for(int j = 0; j<11; j++){
        valid[i][j] = false;
        if(nrOfParentLabels[i]>0){
          double nr = (double)nrOfParentConnections[i][j]/(double)nrOfParentLabels[i];
          if(i!=j and j!=airBubbleLabel and j!=10 and j!=0){ // ignore same layer, airbubbles, columella and unassigned
            if(nr>high1) {
              high2 = high1;
              high2pos = high1pos;
              high1 = nr;
              high1pos = j;
            } else if(nr>high2) {
              high2 = nr;
              high2pos = j;
            }
          }
        }
      }
      valid[i][high1pos] = true;
      if(i!=rootCapLabel and i!=vascLabel) // only one valid connection for rootcap and vasculature
        valid[i][high2pos] = true;
    }

  }

  for(int i = 0; i<numCells; i++){
    int currentLabel = (*cellAtlasAttr)[i].cellLabel;//rcp.rootData.uniqueLabels[i];
    if(cellSelected[currentLabel]){
      int currentNeighbors = cellConnections[currentLabel].size();
      int currentParent = m1->parents()[currentLabel];
      int errorCount = 0;
      for(int j = 0; j<currentNeighbors; j++){
        int currentNeighborParent = m1->parents()[cellConnections[currentLabel][j]];
        if(valid[currentParent][currentNeighborParent] or currentParent==10 or currentNeighborParent==10 or currentParent==airBubbleLabel or
           currentNeighborParent==airBubbleLabel or currentParent==currentNeighborParent){
        // perfect
        } else { // error detected
          errorCount++;
        }
      }
      //std::cout << "cell " << currentLabel << "/" << currentNeighbors << "//" << errorCount << "///" << errors << std::endl;
      if(errorCount>=errors){
        mislabelled.insert(currentLabel);
      }
    }
  }

  cout << "Allowed Connections:" << endl;
  for(int i = 0; i<11; i++){
    if(nrOfParentLabels[i]>0 and i!=airBubbleLabel and i!=10 and i!=0){
      cout << "Parent Label " << i << " with";
      for(int j = 0; j<11; j++){
      //cout << (double)nrOfParentConnections[i][j]/(double)nrOfParentLabels[i] << "  ";
        if(valid[i][j])
          cout << " " << j << " ";
      }
      cout << " " << endl;
    }

  }

  int counterMis = 0;

  for(int i = 0; i<numCells; i++){
    int currentLabel = (*cellAtlasAttr)[i].cellLabel;//rcp.rootData.uniqueLabels[i];
    if(mislabelled.find(currentLabel) != mislabelled.end()){
      counterMis++;
      (*cellAtlasAttr)[i].mislabelled = true;
    }

  }

  cout << "There are  " << counterMis << " mislabelled cells." << endl;

  // select all wrong cells
  forall(const vertex& v, m1->graph()){
    if(mislabelled.find(v->label)!= mislabelled.end()){
      v->selected = true;
    }
    else v->selected = false;
  }

  m1->correctSelection(true);
  m1->updateSelection();

    return true;
  }
  REGISTER_PROCESS(TopologicalCheck);

  REGISTER_PROCESS(CollapseBezierCellAtlas);

// new december 2014/ january 2015
// Data Analysis Tool
 bool DataAnalysis::initialize(QWidget* parent)
  {
    QString filename = parm("Input folder");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getExistingDirectory(0, "Choose folder for control", QDir::currentPath(), QFileDialog::ShowDirsOnly);
    setParm("Input folder", filename);

    QString filename2 = parm("Reporter file string");
    if(filename2.isEmpty() and parent)
      filename2 = QFileDialog::getExistingDirectory(0, "Choose folder for treatment", QDir::currentPath(), QFileDialog::ShowDirsOnly);
    setParm("Reporter file string", filename2);

    QString filename3 = parm("Merge With File");
    if(filename3.isEmpty() and parent)
      filename3 = QFileDialog::getExistingDirectory(0, "Choose folder for output", QDir::currentPath(), QFileDialog::ShowDirsOnly);
    setParm("Merge With File", filename3);

    return true;
  }

  bool DataAnalysis::run(QString folderControl, QString folderTreatment, QString folderOutput, QString outputFileType, QString avg, double window, bool perCell)
  {

    DataAnalyzer da;
    bool avgBool = stringToBool(avg);

    bool treatment = true;
    if(outputFileType == "Control")
      treatment = false;

    // everything happens in the DataAnalyzer class
    da.analyze(folderControl, folderTreatment, folderOutput, treatment, avgBool, window, perCell);

    return true;
  }

  REGISTER_PROCESS(DataAnalysis);

  bool DataAnalysisGUS::initialize(QWidget* parent)
  {
    QString filename = parm("Input folder");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getExistingDirectory(0, "Choose folder for input", QDir::currentPath(), QFileDialog::ShowDirsOnly);
    setParm("Input folder", filename);
    bool mergeBool = stringToBool(parm("Merge With File"));
    if(mergeBool){
      QString filename4 = parm("Output File");
      if(filename4.isEmpty() and parent)
        filename4 = QFileDialog::getOpenFileName(0, "Choose output file", QDir::currentPath(), "CSV files (*.csv)");
      setParm("Output File", filename4);
    }

    return true;
  }

bool DataAnalysisGUS::run(Mesh *m1, QString folderInput, QString gusFileExt,/* QString folderOutput, */QString mergeWithFile, QString fileToMerge, QString avg, double window, QString upperFilter, double upperFilterValue, QString lowerFilter, double lowerFilterValue)
  {

    DataAnalyzer da;
    bool mergeBool = stringToBool(mergeWithFile);
    bool avgBool = stringToBool(avg);

    IntDoubleMap GUSInfo;
    // everything important happens in the DataAnalyzer class
    da.analyzeGUS(folderInput, /*folderOutput,*/ gusFileExt, mergeBool, fileToMerge, avgBool, window, upperFilter, upperFilterValue, lowerFilter, lowerFilterValue, GUSInfo);

    // now update the heatmap
    m1->labelHeat().clear();
    //IntDoubleMap::iterator it;

    double minVal = 1E20, maxVal = -1E20;

    for (IntDoubleMap::iterator it = GUSInfo.begin(); it != GUSInfo.end(); ++it){
      if(it->second > maxVal) maxVal = it->second;
      if(it->second < minVal) minVal = it->second;
    }

    for (IntDoubleMap::iterator it = GUSInfo.begin(); it != GUSInfo.end(); ++it){
      m1->labelHeat()[it->first] = (it->second - minVal)/(maxVal - minVal);
    }

    m1->heatMapBounds() = Point2f(minVal, maxVal);
    m1->heatMapUnit() = QString("GUS");
    m1->showHeat();
    m1->updateTriangles();


    return true;
  }

  REGISTER_PROCESS(DataAnalysisGUS);



// Assign the Columella cells to rootcap and (inner) columella
bool AssignColumella::run(int labelRootCap, int labelCol, double rat_val, double sca_val)
  {
    Mesh *m1 = mesh(0);

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    RootCellProcessing rcp;

    rcp.setDataPointers(cellAtlasAttr, cellAtlasConfigAttr);
    rcp.assignColumella(labelRootCap, labelCol, rat_val, sca_val);
    updateParentLabels(m1);

    return true;
  }

  REGISTER_PROCESS(AssignColumella);

bool AssignCorticalCells::run(int labelCort)
  {

    Mesh *m1 = mesh(0);

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    if(!(*cellAtlasConfigAttr)[0].cellTypesAssigned){
      return setErrorMessage("Warning: Assign Cell Types first.");
    }

    // check whether there are some cortical cells
    bool existCortCells = false;
    for(int i = 0; i<numCells; i++){
      //int currentLabel = (*cellAtlasAttr)[i].cellLabel; //rcp.rootData.uniqueLabels[i];
      if((*cellAtlasAttr)[i].parentLabel == labelCort){
        existCortCells = true;
      }
    }

    if(!existCortCells){
      return setErrorMessage("Warning: no cortical cells found. Please make sure you have clustered the root and specified the correct label of the cortical cells");
    }

    // get bezier in high resolution and minmax of the arclengths
    double minS, maxS;

    RootCellProcessing rcp;

    rcp.setDataPointers(cellAtlasAttr, cellAtlasConfigAttr);
    rcp.assignCortCellsGetMinMax(labelCort, minS, maxS);

    std::vector<double> sSmooth;
    int lengthS = 1000;
    double stepS = (maxS-minS)/((double)lengthS-1.0);

    for(int i = 0; i<lengthS; i++){
      sSmooth.push_back(minS+stepS*i);
    }

    std::map<int, Point3d> bMapS,diffbMapS;
    double bezStart = 0.0;
    double bezEnd = 1.0;
    double stepSize = (double)(bezEnd-bezStart)/(double)(lengthS-1);

    CuttingSurface* cutSurf = cuttingSurface();

    for(int i=0; i<lengthS; i++){
      double u = bezStart+i*stepSize;
      bMapS[i] = Point3d(cutSurf->evalCoord(u,0));
    }
    diffbMapS[0] = bMapS[1] - bMapS[0];
    diffbMapS[lengthS-1] = bMapS[lengthS-1] - bMapS[lengthS-2];

    for(int i=1; i<lengthS; i++){
      diffbMapS[i] = (bMapS[i+1] - bMapS[i-1])/(double)(2);
    }

    // run the main function
    rcp.assignCortCells(labelCort, sSmooth, bMapS, lengthS);

    writeCellAtlasMeasures(m1);

    return true;
  }

  REGISTER_PROCESS(AssignCorticalCells);

bool ExamineBadVasculature::run(Mesh *m1, int labelVasc, int labelAirBubble, int labelEndo, int methodAvg, double sigma)
  {
    std::map<int,double> landscape1D;
    int sizeLandsc = 160;
    double minC=1E20, maxC=-1E20;

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }


    // go through all mislabelled cells and find vasc
    for(int i = 0; i<numCells; i++){
      /*int currentLabel = rcp.rootData.uniqueLabels[i];
      if(rcp.rootData.mislabelled[currentLabel] and rcp.cellLabel[currentLabel] == labelVasc){
        if(minC > rcp.rootData.lengthCirc[currentLabel]) minC = rcp.rootData.lengthCirc[currentLabel];
        if(maxC < rcp.rootData.lengthCirc[currentLabel]) maxC = rcp.rootData.lengthCirc[currentLabel];
      }*/
      if((*cellAtlasAttr)[i].mislabelled and (*cellAtlasAttr)[i].parentLabel == labelVasc){
        if(minC > (*cellAtlasAttr)[i].sizeCirc) minC = (*cellAtlasAttr)[i].sizeCirc;
        if(maxC < (*cellAtlasAttr)[i].sizeCirc) maxC = (*cellAtlasAttr)[i].sizeCirc;
      }
    }

    //double sigma = 1.5*3; // fix parameter
    double stepX = 1.0/(double)(sizeLandsc-1.0);
    for(int j = 0; j<numCells; j++){
      //int currentLabel = (*cellAtlasAttr)[j].cellLabel;
      if((*cellAtlasAttr)[j].mislabelled and (*cellAtlasAttr)[j].parentLabel == labelVasc){

        // generate the gaussian function

        // calc center of gaussian
        int binC = interpolateArrayIndex((*cellAtlasAttr)[j].sizeCirc,minC,maxC,sizeLandsc-1);

        // cut off radius around the center
        int radius = std::max(10.0,sigma*2.0);
        for(int i = -radius; i<=radius; i++){
          // calc gauss
          if(binC+i>=0 and binC+i<sizeLandsc){
            double xValue = interpolateArrayIndex((*cellAtlasAttr)[j].sizeCirc+i*stepX,minC,maxC,sizeLandsc-1); // i or j ???
            landscape1D[binC+i] += gauss1D(xValue, binC, sigma);
          }
        }
      }
    }

    // check for minima in that function
    std::vector<int> minVec;

    // check beginning and end
    if(landscape1D[0] < landscape1D[1])  minVec.push_back(0);
    if(landscape1D[sizeLandsc-1] < landscape1D[sizeLandsc-2])  minVec.push_back(sizeLandsc-1);

    // check everything else
    for(int i = 1; i<sizeLandsc-1; i++){
      if(landscape1D[i] < landscape1D[i+1] and landscape1D[i] > landscape1D[i-1]) minVec.push_back(i);
    }

    // more than 2 minima: take the lowest as a separation point:
    // big cells are now endo, small cells are now bubbles (=matlab procedure)
    if(minVec.size() > 2){
      double min = 1E20;
      int minIdx = 0;
      for(uint j = 0; j<minVec.size(); j++){
        if(min > landscape1D[minVec[j]]){
          min = landscape1D[minVec[j]];
          minIdx = j;
        }
      }

      double separationPoint = minC + minVec[minIdx]*(maxC-minC)/(double)(sizeLandsc-1.0);

      double avgCirc = 0;
      int counter = 0;

      for(int j = 0; j<numCells; j++){
        //int currentLabel = rcp.rootData.uniqueLabels[j];
        if((*cellAtlasAttr)[j].mislabelled and (*cellAtlasAttr)[j].parentLabel == labelEndo){
        //if(!rcp.rootData.mislabelled[currentLabel] and rcp.cellLabel[currentLabel] == labelEndo){
          counter++;
          avgCirc+=(*cellAtlasAttr)[j].sizeCirc;
        }
      }
      avgCirc/=counter;

      // matlab procedure doesnt work very well -> new simple method
      if(methodAvg==1){
        separationPoint = avgCirc/2.0;
      }
      for(int j = 0; j<numCells; j++){
        int currentLabel = (*cellAtlasAttr)[j].cellLabel;
        if((*cellAtlasAttr)[j].mislabelled and (*cellAtlasAttr)[j].parentLabel == labelVasc){
        //if(rcp.rootData.mislabelled[currentLabel] && rcp.cellLabel[currentLabel] == labelVasc){
          if((*cellAtlasAttr)[j].sizeCirc>separationPoint){
            //rcp.cellLabel[currentLabel] = labelEndo;// some rcp left!
            (*cellAtlasAttr)[j].parentLabel = labelEndo;
            cout << "relabelled cell " << currentLabel << " to Endo" << endl;
          } else {
            //rcp.cellLabel[currentLabel] = labelAirBubble;// some rcp left!
            (*cellAtlasAttr)[j].parentLabel = labelAirBubble;
            cout << "relabelled cell " << currentLabel << " to Air Bubble" << endl;
          }
        }
      }
    }

    updateParentLabels(m1);
    return true;
  }

  REGISTER_PROCESS(ExamineBadVasculature);

// some new (experimental) processes (post-publication)
/////////////////////////////////

  // displays heatmap for how far cells are away (in terms of neighborhood graph) compared to selected cells
bool DisplayShortestPath::run(Mesh *m1, double wallThreshold, QString weight)
  {

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    int numCells = (*cellAtlasAttr).size();

    if(numCells == 0 or (*cellAtlasConfigAttr).size() == 0){
      return setErrorMessage("The Cell Atlas Data Attribute Map is empty. Please run Analyze Cells first.");
    }

    vvGraph& S = m1->graph();
    // find two selected cells

    int label1, label2;

    findTwoSelectedLabels(S, label1, label2);

    bool equalWeights = weight == "Euclidian" ? false : true;

    // create maps for aStar function

    std::vector<int> uniqueLabels;
    std::map<int,Point3d> cellCentroids;
    for(int i = 0; i<numCells; i++){
      int currentLabel  = (*cellAtlasAttr)[i].cellLabel;
      uniqueLabels.push_back((*cellAtlasAttr)[i].cellLabel);
      cellCentroids[currentLabel] = (*cellAtlasAttr)[i].centroid;
    }

    // A*
    std::set<int> travelledLabels = aStar(uniqueLabels, cellCentroids, (*cellAtlasConfigAttr)[0].wallArea, label1, label2, wallThreshold, equalWeights);

    // set heat map
    //float delta = maxValue - minValue;
    m1->wallHeat().clear();
    m1->labelHeat().clear();
    for(int i = 0; i<numCells; i++){
      //int currentLabel = rcp.rootData.uniqueLabels[i];
      int currentLabel = (*cellAtlasAttr)[i].cellLabel;
      int value = 0;
      if(travelledLabels.find(currentLabel) != travelledLabels.end()) value = 1;
      m1->labelHeat()[currentLabel] = value;
    }
    m1->heatMapBounds() = Point2f(0, 1);
    m1->updateTriangles();
    m1->showHeat();
    SETSTATUS(QString("Loaded heat map with values ranging from %1 to %2").arg(0).arg(1));

    return true;
  }

  REGISTER_PROCESS(DisplayShortestPath);

  void getWallArea2D(RootCellProcessing& rcp, Mesh* m)
  {

    const vvGraph& S = m->graph();

    forall(const vertex& v, S){
      if(v->label < 1) continue;
      forall(const vertex& n, S.neighbors(v)){
        forall(const vertex& m, S.neighbors(n)){
          if(m->label > 0){
            rcp.rootData.wallArea[std::make_pair(v->label,m->label)] = 1;
            rcp.rootData.wallArea[std::make_pair(m->label,v->label)] = 1;
          }
        }
      }
    }

  }

   // file dialogue for save data
   bool WriteGrowthMapData::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Filename", filename);
    return true;
  }

   bool WriteGrowthMapData::run(Mesh *mesh, QString filename)
  {

	  AttrMap<int, double>& lon3 = mesh->attributes().attrMap<int, double>("Measure Label Double /Cell Atlas/Longitudinal Coord Min");
    AttrMap<int, double>& lon4 = mesh->attributes().attrMap<int, double>("Measure Label Double /Cell Atlas/Longitudinal Coord Max");
    CellTissue &T = mesh->tissue();
    //double Xpos, Ypos, Zpos = 0;
    mesh->updateAll();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);
    out << "Label,maxX, minX,LabelHeat" << endl;
    forall(const cell &c, T.C) {

    double maxX = -HUGE_VAL;
    double minX = HUGE_VAL;
    vertex center = T.getVtx(c);
    forall(const vertex& m, T.S.neighbors(center))  {

      Point3d vCord= m->pos;
    if (vCord[0] > maxX)
      maxX = vCord[0];
    if (vCord[0] < minX)
      minX = vCord[0];
      }
    SymmetricTensor& cellAxisValues = mesh->cellAxis()[c->label]; // TODO check for parent
    Point3f& eig = cellAxisValues.evals();
    double cellGrowth = eig.x();
    QTextStream out(&file);

    if (mesh->useParents())
      out << c->label << "," << lon4[c->label] << " , " << lon3[c->label] <<" , "<<  mesh->labelHeat()[mesh->parents()[c->label]] << endl;
    else
      out << c->label << "," << lon4[c->label] << " , " << lon3[c->label] <<" , "<<  mesh->labelHeat()[c->label] << endl;
  }
    return true;
  }
  REGISTER_PROCESS(WriteGrowthMapData);

/*
 // not needed anymore, function build into image data
 bool SeedParents::run(Mesh *m)
  {

  vvGraph& S = m->graph();

  // Get cells by connected regions
  VVGraphVec cellVertex;
  VtxIntMap vertexCell;
  m->getConnectedRegions(S, cellVertex, vertexCell);

  // Relabel vertices
  forall(const vvGraph& S, cellVertex) {
    bool relabel = false;
    std::unordered_map<int,int> labelCount;
    int parentLabel = 0;
    forall(const vertex& v, S){
      labelCount[v->label]++;
      if(m->parents()[v->label] == v->label){
        relabel = true;
        parentLabel = v->label;
      }
    }
    if(relabel){
      // find highest p.second
      int maxValue = -1;
      int label = 0;
      forall(IntInt p, labelCount){
        if(p.first > 0 and p.second > maxValue){
          maxValue = p.second;
          label = p.first;
        }
      }
      forall(const vertex& v, S){
        m->parents()[v->label] = parentLabel;
        v->label = label;
      }
    }
  }

  // Tell the system the mesh color has changed
  m->updateTriangles();
  m->updateSelection();

  }
  REGISTER_PROCESS(SeedParents);
*/

   // file dialogue for save data
   bool LabelledStackComparison::initialize(QWidget* parent)
  {
    QString filename = parm("Filename");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "CSV files (*.csv)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".csv", Qt::CaseInsensitive))
      filename += ".csv";
    setParm("Filename", filename);
    return true;
  }

   bool LabelledStackComparison::run(Stack* s1, Stack* s2, QString filename, bool overwrite)
  {
    Store* store1 = s1->currentStore();
    Store* store2 = s2->currentStore();

    // convert world coord stack1 into prim coords celldata1
    const HVecUS& data1 = store1->data();
    HVecUS& data2 = store2->data();
    Point3i imgSize(s1->size());
    Point3i imgSize2(s2->size());

    std::map<Point3i, Point3i> pointMap;

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

    std::map<int, Point3d> labelPosMap;
    std::map<int, int> voxelCounterMap;

    #pragma omp parallel for schedule(guided)
    for(int x=0; x<=xDim; x++){
      for(int y=0; y<yDim; y++){
        for(int z=0; z<zDim; z++){
          Point3i pOut (x,y,z);
          int label = data1[s1->offset(pOut)];
          labelPosMap[label] += Point3d(pOut);
          voxelCounterMap[label]++;
        }
      }
    }

    typedef std::pair<int, Point3d> IntP3dPair;
    forall(IntP3dPair pa, labelPosMap){
      if(voxelCounterMap[pa.first] == 0) continue;
      Point3d pos = pa.second/voxelCounterMap[pa.first];
    }

    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("File '%1' cannot be opened for writing").arg(filename));
      return false;
    }
    QTextStream out(&file);
    out << "Label Stack 1,Associated Label Stack 2,Nr of Voxels,Pos X,Pos Y,Pos Z" << endl;

    std::map<int,int> labelS2toS1;

    forall(IntP3dPair pa, labelPosMap){
      if(voxelCounterMap[pa.first] == 0) continue;
      Point3d pos = pa.second/voxelCounterMap[pa.first];
      Point3i pos2 (round(pos.x()), round(pos.y()), round(pos.z()));
      int labelStack2 = 0;
      if(imgSize2.x() > pos2.x() or imgSize2.y() > pos2.y() or imgSize2.z() > pos2.z()){
        labelStack2 = data1[s1->offset(pos2)];
        labelS2toS1[labelStack2] = pa.first;
      }

      out << pa.first << "," << labelStack2 << " , " << voxelCounterMap[pa.first] << " , " << pos2.x() <<" , "<<  pos2.y() << " , " << pos2.z() << endl;
    }

    if(overwrite){
      int xDim2 = imgSize2.x();
      int yDim2 = imgSize2.y();
      int zDim2 = imgSize2.z();

      #pragma omp parallel for schedule(guided)
      for(int x=0; x<=xDim2; x++){
        for(int y=0; y<yDim2; y++){
          for(int z=0; z<zDim2; z++){
            Point3i pOut (x,y,z);
            int label = data2[s2->offset(pOut)];
            if(label != 0 and labelS2toS1[label] != 0)
              data2[s2->offset(pOut)] = labelS2toS1[label];
          }
        }
      }

    }

    store2->changed();
    return true;
  }
  REGISTER_PROCESS(LabelledStackComparison);

  Point3d computeNearestSurfacePoint(const vvGraph& surfaceMesh, Point3d centroid, Point3d minNrml)
  {

    Point3d nearestSurf, nearestNrml;
    double minDis = HUGE_VAL;


    forall(const vertex &v, surfaceMesh){

        double dis = norm(centroid-v->pos);
        if(dis < minDis){
          minDis = dis;
          nearestSurf = v->pos;
          nearestNrml = v->nrml;
        }
    }

    return nearestSurf;
  }


  bool CellAnalysis3D::run(Mesh* m, Mesh* m2, bool useSurface)
  {

    // make sure rcp is empty
    RootCellProcessing rcpNew;

    this->rcp = rcpNew;

    vvGraph& S = m->graph();
    const vvGraph S2 = m2->graph();

    CellAtlasAttr *cellAtlasAttr;
    cellAtlasAttr = &m->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

    CellAtlasConfigAttr *cellAtlasConfigAttr;
    cellAtlasConfigAttr = &m->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    RootCellAnalyzing rca(cellAtlasAttr,cellAtlasConfigAttr);

    uniqueLabels = findAllLabels(S);

    // create label - vertex map
    forall(const vertex &v, S) {
      lvMap[v->label].push_back(v);
    }
    std::cout << "Generate Label/Tri Map" << std::endl;
    generateLabelTriangleMap(S, cellTriangles);

    //rca.generateLabelTriangleMap(uniqueLabels.size(), uniqueLabels, S, lvMap, cellTriangles);

    double minVolume = 0;

    rcp.rootData.uniqueLabels = uniqueLabels;

    std::cout << "Compute Volumes" << std::endl;
    rca.analyzeCellCentroidsVolumes(cellTriangles, rcp, minVolume);

    double tolerance = 0.0001;

    std::cout << "Compute Cell Neighborhood" << std::endl;
    NhbdGraphInfo nbhInfo;
    neighborhoodGraph(S, tolerance, nbhInfo, false);

    rcp.rootData.wallArea = nbhInfo.sharedWallArea;
    rcp.rootData.cellWallArea = nbhInfo.cellWallArea;
    rcp.rootData.outsideWallArea = nbhInfo.outsideWallArea;

    std::cout << "Compute Cell Sizes" << std::endl;
    rca.createCartesianCoordSystem(rcp);
    rca.estimateAllCellLengths(rcp, cellTriangles);

    AttrMap<IntInt, double>& wallArea = m->attributes().attrMap<IntInt, double>("Shared Wall Areas");
    wallArea.clear();

    std::cout << "Write Attr Maps" << std::endl;

    // generate measure attr maps
    AttrMap<int, double>& attrDataVol = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Volume");
    attrDataVol.clear();
    AttrMap<int, double>& attrDataWallArea = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Wall Area");
    attrDataWallArea.clear();
    AttrMap<int, double>& attrDataOutWallArea = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area");
    attrDataOutWallArea.clear();
    AttrMap<int, double>& attrDataOutWallAreaRatio = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Outside Wall Area Ratio");
    attrDataOutWallAreaRatio.clear();
    AttrMap<int, double>& attrDataNeigh = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Neighbors");
    attrDataNeigh.clear();
    AttrMap<int, Point3d>& attrDataCentr = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
    attrDataCentr.clear();
    AttrMap<int, Point3d>& attrDataCentrP = m->attributes().attrMap<int, Point3d>("Measure Label Vector ParentCentroids");
    attrDataCentrP.clear();

    AttrMap<int, Point3d>& attrDataSurfP = m->attributes().attrMap<int, Point3d>("Measure Label Vector SurfacePoint");
    attrDataSurfP.clear();
    AttrMap<int, Point3d>& attrDataSurfN = m->attributes().attrMap<int, Point3d>("Measure Label Vector SurfaceNormal");
    attrDataSurfN.clear();

    AttrMap<int, Point3d>& attrDataSurfPP = m->attributes().attrMap<int, Point3d>("Measure Label Vector SurfacePointParents");
    attrDataSurfPP.clear();
    AttrMap<int, Point3d>& attrDataSurfNP = m->attributes().attrMap<int, Point3d>("Measure Label Vector SurfaceNormalParents");
    attrDataSurfNP.clear();

    AttrMap<int, double>& attrDataSizeX = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length X");
    attrDataSizeX.clear();
    AttrMap<int, double>& attrDataSizeY = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length Y");
    attrDataSizeY.clear();
    AttrMap<int, double>& attrDataSizeZ = m->attributes().attrMap<int, double>("Measure Label Double Geometry/Cell Length Z");
    attrDataSizeZ.clear();

    std::map<int, std::set<int> > countNeighbors;

    forall(IntIntDouPair p, nbhInfo.sharedWallArea){
      wallArea[p.first] = p.second;

      if(p.second > 0.001){
        countNeighbors[p.first.first].insert(p.first.second);
        countNeighbors[p.first.second].insert(p.first.first);
      }
    }

    std::map<int, double> parentVolMap;

    typedef std::pair<int, double> IntDoubleP;
    forall(IntDoubleP p, rcp.rootData.cellVolumes){
      int cellLabel = p.first;
      attrDataVol[cellLabel] = p.second;

      attrDataCentr[cellLabel] = rcp.rootData.cellCentroids[cellLabel];

      

      attrDataOutWallArea[cellLabel] = rcp.rootData.outsideWallArea[cellLabel];
      attrDataWallArea[cellLabel] = rcp.rootData.cellWallArea[cellLabel];
      attrDataOutWallAreaRatio[cellLabel] = rcp.rootData.outsideWallArea[cellLabel]/rcp.rootData.cellWallArea[cellLabel];

      attrDataSizeX[cellLabel] = rcp.rootData.lengthLong[cellLabel];
      attrDataSizeY[cellLabel] = rcp.rootData.lengthCirc[cellLabel];
      attrDataSizeZ[cellLabel] = rcp.rootData.lengthRad[cellLabel];

      attrDataNeigh[cellLabel] = countNeighbors[cellLabel].size();

      if(useSurface){
        Point3d sNrml;
        attrDataSurfP[cellLabel] = computeNearestSurfacePoint(S2, attrDataCentr[cellLabel], sNrml);
        attrDataSurfN[cellLabel] = sNrml;
      }
      // parent label maps
      int cellParent = m->parents()[cellLabel];
      if(cellParent == 0) continue;
      parentVolMap[cellParent] += attrDataVol[cellLabel];

      attrDataCentrP[cellParent] += attrDataCentr[cellLabel] * attrDataVol[cellLabel];
    }

    std::set<int> allParents;
    forall(auto p, m->parents()){
      int cellParent = p.second;
      int cellLabel = p.first;

      if(cellParent == 0) continue;

      allParents.insert(cellParent);
    }
    forall(int cellParent, allParents){

      if(parentVolMap[cellParent] < 1E-5) attrDataCentrP[cellParent] = Point3d(0,0,0);
//std::cout << "parent centroid1 " << cellParent << "/" << attrDataCentrP[cellParent] << std::endl;
      attrDataCentrP[cellParent] /= parentVolMap[cellParent];

      //std::cout << "parent centroid2 " << parentVolMap[cellParent] << "/" << attrDataCentrP[cellParent] << std::endl;

      if(useSurface){
        Point3d sNrml;
        attrDataSurfPP[cellParent] = computeNearestSurfacePoint(S2, attrDataCentrP[cellParent], sNrml);
        attrDataSurfNP[cellParent] = sNrml;
      }
    }




    return true;
  }
  REGISTER_PROCESS(CellAnalysis3D);

     bool GrowthAnalysis3D::run(Mesh* m, Mesh* m2, QString name, QString t1, QString t2)
    {

      // run Cell analysis 3D on both meshes
      std::cout<<"Run Cell Analysis 3D on Mesh1"<<std::endl;
      CellAnalysis3D *cellanalysis3D_1;
      if(!getProcess("Mesh/Heat Map/Analysis/Cell Analysis 3D", cellanalysis3D_1))
        throw(QString("MeshProcessCellAnalysis3D:: Unable to make CellAnalysis3D process"));
      cellanalysis3D_1->run();
      std::cout<<"... done!"<<std::endl;

      Measure3DVolume *Vol1;
      if(!getProcess("Mesh/Heat Map/Measures 3D/Geometry/Volume", Vol1))
        throw(QString("MeshProcessHeatMap:: Unable to make MeasureVolume process"));
      Vol1->run();

      //Export to attr map
      //default + save to other mesh (= the first mesh)
      CreateAttrMap *save_attr;
      if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
        throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
      save_attr-> setParm("Name", name + "CellVolume" + t1);
      save_attr->setParm("Key", "Label");
      save_attr->setParm("Heat Map Type", "Label Heat");
      save_attr->setParm("Save to", "Active Mesh");
      save_attr->run();

      std::cout<<"Run Cell Analysis 3D on Mesh2"<<std::endl;
      CellAnalysis3D *cellanalysis3D_2;
      if(!getProcess("Mesh/Heat Map/Analysis/Cell Analysis 3D", cellanalysis3D_2))
        throw(QString("MeshProcessCellAnalysis3D:: Unable to make CellAnalysis3D process"));

      // switch to the other mesh
      if(m == mesh(0)){
        cellanalysis3D_2->setCurrentMeshId(1);
      } else {
        cellanalysis3D_2->setCurrentMeshId(0);
      }
      cellanalysis3D_2->run();
      std::cout<<"... done!"<<std::endl;

      //Save T2 area
      std::cout<<"Create clean and save AreaT2 on the other mesh"<<std::endl;
      HeatAndChangeMap *Vol2;
      if(!getProcess("Mesh/Heat Map/Heat Map", Vol2))
        throw(QString("MeshProcessHeatMap:: Unable to make HeatAndChangeMap process"));
      //Area2->setParm("Parent based", "Yes"); // parent based heat map
      Vol2->setParm("Measure", "Geometry/Volume"); // parent based heat map
      Vol2->setParm("Use Labels Active Mesh", "No");
        if(m == mesh(0)){
          Vol2->setCurrentMeshId(1);
        } else {
          Vol2->setCurrentMeshId(0);
        }
        m2->setUseParents(true);
        m2->updateTriangles();
      Vol2->run();

      //Export to attr map
      //default + save to other mesh (= the first mesh)
      if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
        throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
      save_attr->setParm("Name", name + "CellVolume" + t2);
      save_attr->setParm("Key", "Parent");
      save_attr->setParm("Heat Map Type", "Parent Heat");
      save_attr->setParm("Save to", "Other Mesh");
      if(m == mesh(0))
        save_attr->setCurrentMeshId(1);
      else
        save_attr->setCurrentMeshId(0);
      save_attr->run();

      // switch back to the first mesh
      if(m == mesh(0))
        save_attr->setCurrentMeshId(0);
      else
        save_attr->setCurrentMeshId(1);

      std::cout<<"Create and save volume increase"<<std::endl;
      CombineAttrMaps *VolRatio;
      if(!getProcess("Mesh/Heat Map/Operators/Heat Map Transform", VolRatio))
        throw(QString("MeshProcessHeatMap:: Unable to make CombineAttrMaps process"));

      VolRatio->setParm("Attr Map 1", name + "CellVolume" + t2);
      VolRatio->setParm("Attr Map 2", name + "CellVolume" + t1);
      VolRatio->setParm("Combination", "Yes");
      VolRatio->setParm("Combination Type", "Ratio (A/B)");
      VolRatio->setParm("Transformation", "No");
      VolRatio->setParm("New Attr Map", name + "GrowthVolume" + t1 + t2);
      VolRatio->run(); //Should save growth

      HeatMapProliferation *Prolif;
      if(!getProcess("Mesh/Lineage Tracking/Heat Map Proliferation", Prolif))
        throw(QString("MeshProcessHeatMap:: Unable to make HeatMapProliferation process"));

      Prolif->setParm("Use Parents On Other Mesh", "Yes");
      Prolif->run();

      //Export to attr map
      //default + save to other mesh (= the first mesh)
      //CreateAttrMap *save_attr;
      if(!getProcess("Mesh/Heat Map/Operators/Export Heat to Attr Map", save_attr))
        throw(QString("MeshProcessHeatMap:: Unable to make CreateAttrMap process"));
      save_attr-> setParm("Name", name + "Proliferation" + t1 + t2);
      save_attr->setParm("Key", "Label");
      save_attr->setParm("Heat Map Type", "Label Heat");
      save_attr->setParm("Save to", "Active Mesh");
      save_attr->run();

      return true;
    }
    REGISTER_PROCESS(GrowthAnalysis3D);

/*
  bool SmoothCellBorder::run(Mesh* m, double pointDis)
  {

  // get cell border segments from junction to junction

  std:vector<Point3d> cellBorder, cellBorderBez;

  // consider those points as Bezier supporting points
  Bezier b();
  b.setBezPoints(Point2u(bezPointsX, 1));
  b.initBez();
  b.setBezierV(cellBorder);

  // resolve points along this Bezier curve with the given resolution
  for(uint i=0; i < cellBorder.size(); i++){

    Point3d bezP = b.evalCoord(double(i)/double(cellBorder.size()-1),1);
    cellBorderBez.push_back(bezP);
  }

  // problem: more than one border to a cell or outside? -> must treat them separately
  // problem: create new vertices along the border (for starters: do it only on the cell mesh or do not add points)

  }
  return true;

  }
  REGISTER_PROCESS(SmoothCellBorder);
*/

  // TODO Test process
  bool BezierFromCellFile::run(Mesh* m, int maxP)
  {

  if(maxP < 2) return setErrorMessage("More Points Needed");
  if(maxP > 20) return setErrorMessage("Enter Max Point smaller 20");

  AttrMap<int, Point3d>& centroids = m->attributes().attrMap<int, Point3d>("Measure Label Vector CellCentroids");
  if(centroids.empty()) return setErrorMessage("No centroids found.");

  vvGraph& S = m->graph();

  // get neighborhood graph of 3D cells
  double tolerance = 0.0001;
  NhbdGraphInfo nbhInfo;
  neighborhoodGraph(S, tolerance, nbhInfo, false);

  // get selected cells
  std::set<int> selCells = findAllSelectedLabels(S);

  std::map<int, std::set<int> > selectedLabelNeighborMap;
  // find start cell (has only 1 neighbor), or pick random cell (if a ring)
  forall(auto p, nbhInfo.sharedWallArea){
    if(selCells.find(p.first.first) == selCells.end()) continue;
    if(selCells.find(p.first.second) == selCells.end()) continue;

    // here both of the neighboring cells are selected
    selectedLabelNeighborMap[p.first.first].insert(p.first.second);
    selectedLabelNeighborMap[p.first.second].insert(p.first.first);
  }

  bool isRing = true;

  int startLabel = -1;
  int lastLabel = -1;
  int countLabelsOneNeighbor = 0;
  forall(auto p, selectedLabelNeighborMap){
    if(p.second.size() > 2) return setErrorMessage("There are cells with more than 2 neighbors!");
    if(p.second.size() == 1){
      if(startLabel == -1) startLabel = p.first;
      countLabelsOneNeighbor++;
      isRing = false;
    }
    lastLabel = p.first;
  }
  if(countLabelsOneNeighbor > 0 and countLabelsOneNeighbor != 2) return setErrorMessage("Not valid input!");

  if(startLabel == -1) startLabel = lastLabel;

  std::set<int> doneLabels;

  std::vector<int> directedLabels;

  // move through selected cells to end cell
  while(startLabel > 0){
    std::cout << "s " << startLabel << std::endl;
    doneLabels.insert(startLabel);
    directedLabels.push_back(startLabel);

    // get an unprocessed neighbor
    int newLabel = -1;
    forall(int l, selectedLabelNeighborMap[startLabel]){
      if(doneLabels.find(l) != doneLabels.end()) continue;
      newLabel = l;
    }
    if(newLabel == -1){
      // end the loop
      startLabel = -1;
    } else
      startLabel = newLabel;
  }

  while(directedLabels.size() > maxP){
    double minDis = HUGE_VAL;
    int idx1 = -1;
    for(int i = 1; i < directedLabels.size()-2; i++){
      double dis = norm(centroids[directedLabels[i]] - centroids[directedLabels[i+1]]);
      if(dis < minDis){
        minDis = dis;
        idx1 = i;
      }
    }
    std::vector<int> directedLabelsNew;

    for(int i = 0; i < directedLabels.size(); i++){
      if(i == idx1) continue;
      directedLabelsNew.push_back(directedLabels[i]);
    }
    directedLabels = directedLabelsNew;
  }

  // if ring, add the first also to the end to close the ring
  if(isRing) directedLabels.push_back(lastLabel);

  // create a Bezier with the number of points and with centroids as supporting points
  NewBezier nb(*this);
  nb.run(directedLabels.size(), 1, 100, 1, 10, 10);

  CuttingSurface* cutSurf = cuttingSurface();
  std::vector<Point3d>& bezVec = cutSurf->bezier().bezierV();

  uint counter = 0;
  // for(int i = 0; i < directedLabels.size(); i++){
  //   bezVec[i] = centroids[directedLabels[i]];
  //   bezVec[i+directedLabels.size()] = centroids[directedLabels[i]];
  //   bezVec[i+2*directedLabels.size()] = centroids[directedLabels[i]];
  //   bezVec[i+3*directedLabels.size()] = centroids[directedLabels[i]];
  // }
  forall(Point3d &bezPoint, bezVec){
    bezPoint = centroids[directedLabels[counter]];
    counter++;
  }


  return true;

  }
  REGISTER_PROCESS(BezierFromCellFile);


  // TODO Test process
  bool SurfaceAnalyzeBezierLine::run(const Stack* s1, Mesh* m1, bool flipBezier)
  {


  CellAtlasAttr *cellAtlasAttr;
  cellAtlasAttr = &m1->attributes().attrMap<int, CellAtlasData>("Cell Atlas 3D Data");

  CellAtlasConfigAttr *cellAtlasConfigAttr;
  cellAtlasConfigAttr = &m1->attributes().attrMap<int, CellAtlasConfig>("Cell Atlas 3D Config");

    RootCellProcessing rcp;

    // assign longitudinal, circumferential and radial to a leaf prim

    // meshes
    const vvGraph& segmentedMesh = m1->graph();
    //vvGraph& S2 = m2->graph();

    VVGraphVec cellVertex;
    VtxIntMap vertexCell;

    rcp.rootData.labelFirstCell = findSelectedLabel(segmentedMesh);

    progressStart("Analyze Cells 3D - Find Cells", 0);
    std::vector<int> uniqueLabels = findAllLabels(segmentedMesh);

    rcp.rootData.numCells = uniqueLabels.size();
    rcp.rootData.uniqueLabels = uniqueLabels;

    std::map<int,std::vector<vertex> > lvMap;
    //vertMap vertexMap;

    // create label - vertex map
    progressStart("Analyze Cells 3D - Create label/vertex map", 0);
    forall(const vertex &v, segmentedMesh) {
      lvMap[v->label].push_back(v);
    }

    progressStart("Analyze Cells 3D - Interpolate Bezier", 0);
    // Bezier
    CuttingSurface* cutSurf = cuttingSurface();

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

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

    // create bezier vector and its derivative
    std::map<int, Point3d> bMap,diffbMap,bMapRadRef;
    std::vector<std::vector<Point3d> > bezGridMap;
    int dataPointsBezier = 200; // number of intermediate points
    int midP = dataPointsBezier/2-1;

    discretizeBezierLineRef(cutSurf, 0.5, 0.0, mGLTot, dataPointsBezier, bMap, diffbMap, bMapRadRef);

    // analyze all cells
    std::map<int, triVector> cellTriangles;
    RootCellAnalyzing rca(cellAtlasAttr, cellAtlasConfigAttr);

    if(m1->meshType() == "MGXM"){
      if(!rca.generateLabelTriangleMap(rcp.rootData.numCells, uniqueLabels, segmentedMesh, lvMap, cellTriangles)) return false;
      if(!rca.analyzeCellCentroids2D(cellTriangles, rcp)) return false;
    } else if(m1->meshType() == "MGXC"){
      forall(const vertex& v, segmentedMesh){
        if(v->label < 1) continue;
        rcp.rootData.cellCentroids[v->label] = v->pos;
      }
    }

std::cout << "bez analysis" << std::endl;
    rca.analyzeBezierLine(uniqueLabels, rcp.rootData.cellCentroids, bMap, diffbMap, flipBezier);
    rca.minMaxToBezierPerCell(m1->graph(), uniqueLabels, bMap, flipBezier);

    rca.firstLong = diffbMap[midP];
    rca.firstRad = bMapRadRef[midP];
    rca.firstLong/=norm(rca.firstLong);
    rca.firstRad/=norm(rca.firstRad);

    Point3d planeNormal = rca.firstLong ^ rca.firstRad;
    rca.calcCircumferentialEqual(uniqueLabels, rcp.rootData.cellCentroids);
    rca.writeDataFields(rcp);
    rca.calcSimpleRadial(uniqueLabels, rcp, rcp.rootData.cellCentroids);

  int numCells = rcp.rootData.uniqueLabels.size();

  for(int i=0; i<numCells; i++){
    int currentLabel = rcp.rootData.uniqueLabels[i];
    int index;
    double weight;
    rcp.rootData.dirRad[currentLabel] = rcp.rootData.cellCentroids[currentLabel] - calcNearestPointOnBezierLine(rcp.rootData.cellCentroids[currentLabel], bMap, index, weight);
    rcp.rootData.dirRad[currentLabel] /= norm(rcp.rootData.dirRad[currentLabel]);
  }

std::cout << "generate data" << std::endl;
    rca.calcCoordSystem(uniqueLabels, rca.diffBezInterp, rcp.rootData.dirRad, rcp.rootData.dirLong, rcp.rootData.dirCirc);

  
    rca.writeDataFields(rcp);

    if(m1->meshType() == "MGXC")
      getWallArea2D(rcp, m1);

std::cout << "write attr" << std::endl;
    writeAttrMaps(rcp, cellAtlasAttr, cellAtlasConfigAttr);

    //rcp.setDataPointers(cellAtlasAttr, cellAtlasConfigAttr);

    writeCellAtlasMeasures(m1);
std::cout << "done" << std::endl;

  return true;

  }
  REGISTER_PROCESS(SurfaceAnalyzeBezierLine);
  //  bool DetectSurfaceWalls::run(Mesh* m, Mesh* m2, double maxDis, QString mode)
  // {

  //   // m1 segmented
  //   // m2 surface

  //   // build map of vtxs

  //   vvGraph& S = m->graph();
  //   vvGraph& S2 = m2->graph();

  //   typedef std::pair<double, vertex> DoubVtxP;
  //   std::map<double, vertex> xPosVtxMapSeg, xPosVtxMapSurf;

  //   // forall(const vertex& v, S){
  //   //   while(xPosVtxMapSeg.find(v->pos.x()) != xPosVtxMapSeg.end()){
  //   //     v->pos.x()+=1E-10;
  //   //   }
  //   //   xPosVtxMapSeg[v->pos.x()] = v;
  //   // }

  //   forall(const vertex& v, S2){
  //     while(xPosVtxMapSurf.find(v->pos.x()) != xPosVtxMapSurf.end()){
  //       v->pos.x()+=1E-10;
  //     }
  //     xPosVtxMapSurf[v->pos.x()] = v;
  //   }

  //   std::set<int> surfaceCells;
  //   std::set<vertex> surfaceVtxs;

  //   forall(const vertex& v, S){

  //     if(mode == "Label Cells" and surfaceCells.find(v->label) != surfaceCells.end()) continue;

  //     std::map<double, vertex>::iterator it = xPosVtxMapSurf.lower_bound(v->pos.x());
  //     if(it == xPosVtxMapSurf.end())
  //       --it;
  //     while(it != xPosVtxMapSurf.end() and it->first - v->pos.x() < maxDis) {
  //       if(norm(v->pos - it->second->pos) < maxDis){
  //         surfaceVtxs.insert(v);
  //         surfaceCells.insert(v->label);
  //         break;
  //       }
  //     }

  //     // forall(const vertex& n, S2){
  //     //   if(norm(v->pos - n->pos) < maxDis){

  //     //     break;
  //     //   }
  //     // }
  //   }

  //   if(mode == "Label Walls"){
  //     forall(const vertex& v, S){
  //       if(surfaceVtxs.find(v) != surfaceVtxs.end())
  //         v->label = 1;
  //       else
  //         v->label = 2;
  //     }
  //   } else {
  //     m->parents().clear();
  //     forall(int l, surfaceCells){
  //       m->parents()[l] = 1;
  //     }
  //   }

  //   m->updateTriangles();

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


}
