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

#include <SystemProcessSave.hpp>
#include <Triangulate.hpp> // findPolygonPlane
#include <PlyFile.hpp>
#include <Version.hpp>

#include <ui_PlyCellGraphDlg.h>

namespace mgx
{

  // bool savePLYCellGraphMarion(Mesh* mesh, const QString& filename, bool binary)
  // {
  //   QFile file(filename);
  //   if(!file.open(QIODevice::WriteOnly)) {
  //     //setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
  //     std::cout << "MeshExport::Error:Cannot open output file" << std::endl;
  //     return false;
  //   }

  //   const vvGraph& S = mesh->graph();
  //   const CellTissue& T = mesh->tissue();
  //   QTextStream outText(&file);
  //   QDataStream outBinary(&file);
  //   outBinary.setByteOrder(QDataStream::LittleEndian);
  //   outBinary.setFloatingPointPrecision(QDataStream::SinglePrecision);

  //   if(mesh->meshType() != "MGX2D") return false;
  //   progressStart(QString("Saving Cell Graph").arg(mesh->userId()), 0);
  //   // First count the junctions
  //   int cells = 0;
  //   int junctions = 0;
  //   // Cell mesh
  //   forall(const vertex& v, S) {
  //     if(v->type == 'c')
  //       ++cells;
  //     else if(v->type == 'j')
  //       // We are not writing cell centers
  //       v->saveId = junctions++;
  //     else {
  //       //setErrorMessage(QString("MeshExport::Error:Bad vertex type:%1").arg(v->type));
  //       return false;
  //     }
  //   }

  //   //std::vector<int> allLabels = findAllLabels(S);

  //   int cellNr = 0;
  //   int edges = 0;
  //   forall(const cell& c, T.C) {
  //     cellNr++;
  //     forall(const cell& n, T.C.neighbors(c)) {
  //       edges++;
  //     }
  //   }

  //   // Write header
  //   outText << "ply" << endl;
  //   if(binary)
  //     outText << "format binary_little_endian 1.0" << endl;
  //   else
  //     outText << "format ascii 1.0" << endl;
  //   outText << "comment Exported from MorphoGraphX " << VERSION << " " REVISION << endl;
  //   outText << "comment Cell Graph: each vertex is a cell, each edge a wall" << endl;
  //   outText << "element vertex " << cellNr << endl;
  //   outText << "property int label" << endl;
  //   outText << "property float x" << endl;
  //   outText << "property float y" << endl;
  //   outText << "property float z" << endl;
  //   outText << "property int parent" << endl;
  //   outText << "property float cellArea" << endl; // TODO
  //   outText << "property list uchar float cellCurvature" << endl; // TODO
  //   outText << "property float aspectRatio" << endl; // TODO
  //   outText << "property float growthRate" << endl; // TODO
  //   outText << "property list uchar float growthDirections" << endl; // TODO
  //   outText << "element edge " << edges << endl;
  //   outText << "property uint source" << endl;
  //   outText << "property uint target" << endl;
  //   outText << "property float wallLength" << endl; // TODO
  //   outText << "property int newWall" << endl; // TODO
  //   outText << "property float wallElongation" << endl; // TODO
  //   outText << "property float angle" << endl; // TODO
  //   outText << "end_header" << endl;

  //   float placeHolder = 0.f;

  //   AttrMap<int, double>& cellArea                  = mesh->attributes().attrMap<int, double>("Measure Label Double cellArea");
  //   AttrMap<int, SymmetricTensor>& cellCurvature    = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor cellCurvature");
  //   AttrMap<int, double>& aspectRatio               = mesh->attributes().attrMap<int, double>("Measure Label Double aspectRatio");
  //   AttrMap<int, double>& growthRate                = mesh->attributes().attrMap<int, double>("Measure Label Double growthRate");
  //   AttrMap<int, SymmetricTensor>& growthDirections = mesh->attributes().attrMap<int, SymmetricTensor>("Measure Label Tensor growthDirections");

  //   AttrMap<std::pair<int, int>, double>& wallLength    = mesh->attributes().attrMap<std::pair<int, int>, double>("Measure Edge Double wallLength");
  //   AttrMap<std::pair<int, int>, double>& newWall       = mesh->attributes().attrMap<std::pair<int, int>, double>("Measure Edge Double newWall");
  //   AttrMap<std::pair<int, int>, double>& wallElongation= mesh->attributes().attrMap<std::pair<int, int>, double>("Measure Edge Double wallElongation");
  //   AttrMap<std::pair<int, int>, double>& angle         = mesh->attributes().attrMap<std::pair<int, int>, double>("Measure Edge Double angle");

  //   std::cout << "Export File for Marion's R plugin" << std::endl;
  //   if(cellArea.empty()) std::cout << "Attr Map Measure Label Double cellArea not found." << std::endl;
  //   if(cellCurvature.empty()) std::cout << "Attr Map Measure Label Tensor cellCurvature not found." << std::endl;
  //   if(aspectRatio.empty()) std::cout << "Attr Map Measure Label Double aspectRatio not found." << std::endl;
  //   if(growthRate.empty()) std::cout << "Attr Map Measure Label Double growthRate not found." << std::endl;
  //   if(growthDirections.empty()) std::cout << "Attr Map Measure Label Tensor growthDirections not found." << std::endl;

  //   // First write vertices
  //   forall(const cell& c, T.C) {
  //     int lengthTensor = 9;

  //     Point3f pos = Point3f(0,0,0);//v->pos;//savedPos(v->pos, transform, stack);
  //     if(binary){
  //     	// until cellArea
  //       outBinary << c->label << c->pos.x() << c->pos.y() << c->pos.z() << mesh->parents()[c->label] << cellArea[c->label];
  //       // cellCurvature
  //       outBinary << lengthTensor << cellCurvature[c->label].ev1().x()<< cellCurvature[c->label].ev1().y()<< cellCurvature[c->label].ev1().z();
  //       outBinary << cellCurvature[c->label].ev2().x()<< cellCurvature[c->label].ev2().y()<< cellCurvature[c->label].ev2().z();
  //       outBinary << cellCurvature[c->label].evals().x()<< cellCurvature[c->label].evals().y()<< cellCurvature[c->label].evals().z();
  //       // aspectRatio/growthRate
  //       outBinary<< aspectRatio[c->label] << growthRate[c->label];
  //       // growthDirections
  //       outBinary << lengthTensor << growthDirections[c->label].ev1().x()<< growthDirections[c->label].ev1().y()<< growthDirections[c->label].ev1().z();
  //       outBinary << growthDirections[c->label].ev2().x()<< growthDirections[c->label].ev2().y()<< growthDirections[c->label].ev2().z();
  //       outBinary << growthDirections[c->label].evals().x()<< growthDirections[c->label].evals().y()<< growthDirections[c->label].evals().z() << endl;
  //       //outBinary << c->label << c->pos.x() << c->pos.y() << c->pos.z() << mesh->parents()[c->label];
  //       //outBinary << placeHolder
  //       //          << lengthTensor << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder
  //       //          << placeHolder<< placeHolder
  //       //          << lengthTensor << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << endl;// << v->signal;
  //     } else {
  //       // until cellArea
  //       outText << c->label << " " << c->pos << " " << mesh->parents()[c->label] << " " << cellArea[c->label] << " ";
  //       // cellCurvature
  //       outText << lengthTensor << " " << cellCurvature[c->label].ev1() << " " << cellCurvature[c->label].ev2() << " " << cellCurvature[c->label].evals() << " ";
  //       // aspectRatio/growthRate
  //       outText<< aspectRatio[c->label] << " " << growthRate[c->label] << " ";
  //       // growthDirections
  //       outText << lengthTensor << " " << growthDirections[c->label].ev1() << " " << growthDirections[c->label].ev2() << " " << growthDirections[c->label].evals() << " " << endl;
  //       //if(!progressAdvance(0))
  //       //  userCancel();
  //     }
  //   }
  //   // Now write walls (edges)
  //   forall(const cell& c, T.C) {
  //     forall(const cell& n, T.C.neighbors(c)) {
  //       std::pair<int,int> currentEdge = std::make_pair(c->label, n->label);
  //       if(binary)
  //         outBinary << c->label << n->label << wallLength[currentEdge] << newWall[currentEdge] << wallElongation[currentEdge] << angle[currentEdge] << endl;
  //       else
  //         outText << c->label << " " << n->label << " " << wallLength[currentEdge] << " " << newWall[currentEdge] << " " << wallElongation[currentEdge] << " " << angle[currentEdge] << endl;
  //     }
  //   }

  //   return true;
  // }


  bool savePLYCellGraphGeneral2D(Mesh* mesh, std::vector<PlyCellData>& cellVec, std::vector<PlyWallData>& wallVec, const QString& filename, bool binary)
  {
    //std::cout << "save ply " << std::endl;
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      //setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }
    //std::cout << "go " << std::endl;
    vvGraph& S = mesh->graph();
    const CellTissue& T = mesh->tissue();
    QTextStream outText(&file);
    QDataStream outBinary(&file);
    outBinary.setByteOrder(QDataStream::LittleEndian);
    outBinary.setFloatingPointPrecision(QDataStream::SinglePrecision);

    // find all cells and create neighborhoodgraph

    int cellNr = 0;
    int edges = 0;

    std::map<IntInt, double> neighMap;

    std::set<int> cellLabels;
    std::set<IntInt> wallLabels;

    std::cout << "mesh Type" << std::endl;

    mesh->updateCentersNormals();

    if(mesh->meshType() == "MGX2D"){
      forall(const cell& c, T.C) {
        //cellNr++;
        forall(const cell& n, T.C.neighbors(c)) {
          IntInt cellPair = std::make_pair(c->label, n->label);
          neighMap[cellPair] = 1;
        }
      }
    } else if(mesh->meshType() == "MGXM" || mesh->meshType() == "MGXC"){
      neighborhood2D(S, neighMap);
    } else {
      return false;
    }

    forall(auto p, neighMap){
      cellLabels.insert(p.first.first);
      cellLabels.insert(p.first.second);
      if(p.first.first < p.first.second) wallLabels.insert(p.first);
    }

    cellNr = cellLabels.size();
    edges = wallLabels.size();

    // //progressStart(QString("Saving Cell Graph").arg(mesh->userId()), 0);
    // // First count the junctions
    // int cells = 0;
    // int junctions = 0;
    // // Cell mesh
    // forall(const vertex& v, S) {
    //   if(v->type == 'c')
    //     ++cells;
    //   else if(v->type == 'j')
    //     // We are not writing cell centers
    //     v->saveId = junctions++;
    //   else {
    //     //setErrorMessage(QString("MeshExport::Error:Bad vertex type:%1").arg(v->type));
    //     return false;
    //   }
    // }



    std::cout << "write file " << std::endl;
    // Write header
    outText << "ply" << endl;
    if(binary)
      outText << "format binary_little_endian 1.0" << endl;
    else
      outText << "format ascii 1.0" << endl;
    outText << "comment Exported from MorphoGraphX " << VERSION << "." << REVISION << endl;
    outText << "comment Cell Graph: each vertex is a cell, each edge a wall" << endl;
    outText << "element vertex " << cellNr << endl;
    outText << "property int label" << endl;
    outText << "property float x" << endl;
    outText << "property float y" << endl;
    outText << "property float z" << endl;
    forall(const PlyCellData& data, cellVec){
      if(data.type == "uint" or data.type == "float"){
        outText << "property " << data.type << " " << data.name << endl;
      } else if(data.type == "vector"){
        outText << "property list uchar float " << data.name << endl;
      } else if(data.type == "tensor"){
        outText << "property list uchar float " << data.name << endl;
      }
    }
    outText << "element edge " << edges << endl;
    outText << "property uint source" << endl;
    outText << "property uint target" << endl;
    forall(const PlyWallData& data, wallVec){
      outText << "property " << data.type << " " << data.name << endl;
    }
    outText << "end_header" << endl;

    float placeHolder = 0.f;
    int placeHInt = 0;

    //AttrMap<int, double>& areaMap = mesh->attributes().attrMap<int, double>("Measure Label Double /Geometry/Area");

    std::cout << "write data " << std::endl;
    // First write vertices
    //forall(const cell& c, T.C) {
    forall(int c, cellLabels){
      int lengthTensor = 9;

      Point3f pos = Point3f(0,0,0);//v->pos;//savedPos(v->pos, transform, stack);
      if(binary){
        /*   outBinary << c->label << c->pos.x() << c->pos.y() << c->pos.z() << mesh->parents()[c->label]<< placeHolder
          << lengthTensor << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder
          << placeHolder<< placeHolder
          << lengthTensor << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << placeHolder << endl;// << v->signal;
      */} else {
        outText << c << " " << mesh->parentCenterVis()[c].x() << " " << mesh->parentCenterVis()[c].y() << " " << mesh->parentCenterVis()[c].z() << " ";
        forall(PlyCellData& data, cellVec){
          if(data.type == "uint"){
            outText << data.intData[c] << " ";
          } else if(data.type == "float"){
            outText << data.floatData[c] << " ";
          } else if(data.type == "vector"){
            outText << "3 " << data.pointData[c].x() << " "
                    << data.pointData[c].y() << " "
                    << data.pointData[c].z() << " ";
          } else if(data.type == "tensor"){
            outText << "9 " << data.tensorData[c].ev1() << " "
                    << data.tensorData[c].ev2() << " "
                    << data.tensorData[c].evals() << " ";
          }
        }
        outText << endl;
      }

      //if(!progressAdvance(0))
      //  userCancel();
    }
    // Now write walls (edges)
    forall(IntInt w, wallLabels) {
      if(binary)
        outBinary << w.first << w.second << endl;//<< placeHolder << placeHInt << placeHolder << placeHolder << endl;
      else
        outText << w.first << " " << w.second << endl;//<< " " << placeHolder << " " << placeHInt << " " << placeHolder << " " << placeHolder << endl;
    }

    return true;
  }

  bool loadCellGraphPLYGeneral(Mesh* mesh, const QString& filename)
  {
    // Define the reader and read header
    PlyFile ply;
    if(!ply.parseHeader(filename))
      throw(QString("Cannot parse PLY file header"));

    if(!ply.parseContent())
      throw(QString("Unable to parse contents of PLY file"));

    // go through all cells (vertex)
    PlyFile::Element* vtx = ply.element("vertex");

    // get label property and fill idxLabelMap
    PlyFile::Property* cLabel = vtx->property("label");
    std::map<int,int> idxLabelMap;

    std::vector<int>& labels = *cLabel->value<int>();

    for(size_t i = 0; i<vtx->size(); i++){
      idxLabelMap[i] = labels[i];
    }

    // go through all properties and create data
    size_t nrProp = vtx->nbProperties();

    for(size_t i = 0; i<nrProp; i++){
      PlyFile::Property* p = vtx->property(i);
      if(p->name() == "label") continue;
      if(p->kind() == PlyFile::Property::VALUE){
        if(p->memType() == PlyFile::INT){
          std::vector<int>& values = *p->value<int>();
          QString attrName = "Measure Label Int " + p->name();
          AttrMap<int, int>& attrData = mesh->attributes().attrMap<int, int>(attrName);
          attrData.clear();
          for(size_t j = 0; j<vtx->size(); j++){
            attrData[labels[j]] = values[j];
          }

        } else if(p->memType() == PlyFile::FLOAT){

          //PlyFile::Property* t = vtx->property(p->name());
          const std::vector<float>& values = *p->value<float>();
          QString attrName = "Measure Label Double " + p->name();
          AttrMap<int, double>& attrData = mesh->attributes().attrMap<int, double>(attrName);
          attrData.clear();
          for(size_t j = 0; j<vtx->size(); j++){
            attrData[labels[j]] = values[j];
          }
        }
      } else { // list
        if(p->memType() == PlyFile::FLOAT){
          std::vector<std::vector<float> >& values = *p->list<float>();
          if(values[0].size() > 0 and values[0].size() <=3){ // assume vector
            QString attrName = "Measure Label Vector " + p->name();
            AttrMap<int, Point3d>& attrData = mesh->attributes().attrMap<int, Point3d>(attrName);
            attrData.clear();
            for(size_t j = 0; j<vtx->size(); j++){
              if(values[0].size()>0) attrData[labels[j]].x() = values[j][0];
              if(values[0].size()>1) attrData[labels[j]].y() = values[j][1];
              if(values[0].size()>2) attrData[labels[j]].z() = values[j][2];
            }

          } else if(values[0].size() > 3){ // assume tensor
            QString attrName = "Measure Label Tensor " + p->name();
            AttrMap<int, SymmetricTensor>& attrData = mesh->attributes().attrMap<int, SymmetricTensor>(attrName);
            attrData.clear();
            for(size_t j = 0; j<vtx->size(); j++){
              Point3f ev1, ev2, evals;
              ev1.x() = values[j][0];
              ev1.y() = values[j][1];
              ev1.z() = values[j][2];
              ev2.x() = values[j][3];
              if(values[0].size()>4) ev2.y() = values[j][4];
              if(values[0].size()>5) ev2.z() = values[j][5];
              if(values[0].size()>6) evals.x() = values[j][6];
              if(values[0].size()>7) evals.y() = values[j][7];
              if(values[0].size()>8) evals.z() = values[j][8];
              attrData[labels[j]] = SymmetricTensor(ev1,ev2,evals);
            }

          }

        }
      }
    }

    // now cell walls (edge)
    PlyFile::Element* edges = ply.element("edge");

    nrProp = edges->nbProperties();

    // TODO
    for(size_t i = 0; i<nrProp; i++){

    }
    return true;
  }

  // fill a QTreeWidget with all existing measures
  void fillTreeWidgetWithAttrMaps(Mesh* m, QTreeWidget* tree, QStringList *selected)
  {
    Attributes *attributes;
    attributes = &m->attributes();
    QStringList attrList = attributes->getAttrList();
    QStringList measureAttr;

    attrList.sort(Qt::CaseSensitive);

    for(const QString &text : attrList) {
      QStringList list = text.split(" ");

      if(list.size() < 4) 
        continue;
      if(list[0] != "Measure") 
        continue;
      QString name, attrName;
      for(int i = 1; i < list.size(); i++){
        if(i == 1) 
          name = list[i];
        else
          name += " " + list[i];
        if(i == 3)
          attrName = list[i];
        else if(i > 3)
          attrName += " " + list[i];
      }

//qDebug() << "Name:" << name << endl;
//qDebug() << "AttrName:" << attrName << endl;

      QTreeWidgetItem *item = new QTreeWidgetItem(QStringList() << "" << name);
      item->setFlags(Qt::ItemIsEnabled | Qt::ItemIsSelectable | Qt::ItemIsUserCheckable);
      if(selected and selected->contains(attrName))
        item->setCheckState(0, Qt::Checked);
      else
        item->setCheckState(0, Qt::Unchecked);
      tree->addTopLevelItem(item);
    }

  }

  void SavePlyFile::selectAll()
  {
    for(int itemCount = 0; itemCount < ui->attrMapTree->topLevelItemCount(); ++itemCount)
      ui->attrMapTree->topLevelItem(itemCount)->setCheckState(0, Qt::Checked);
  }

  void SavePlyFile::unselectAll()
  {
    for(int itemCount = 0; itemCount < ui->attrMapTree->topLevelItemCount(); ++itemCount)
      ui->attrMapTree->topLevelItem(itemCount)->setCheckState(0, Qt::Unchecked);
  }

  // file dialogue for save data
  bool SavePlyFile::initialize(QWidget* parent)
  {
    Ui_PlyCellGraphDlg ui;
    Mesh *m = currentMesh();
    //if(m->meshType()!="MGX2D" and m->meshType()!="MGX3D") return true;
    QString filename = parm("Filename");
    if(filename.isEmpty() and parent)
      filename = QFileDialog::getSaveFileName(0, "Choose spreadsheet file to save", QDir::currentPath(), "PLY files (*.ply)");
    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".ply", Qt::CaseInsensitive))
      filename += ".ply";
    setParm("Filename", filename);

    if(!parent) return true;

    // now the new GUI
    QDialog dlg(parent);
    ui.setupUi(&dlg);
    this->ui = &ui;

    connect(ui.selectAllButton, SIGNAL(clicked()), this, SLOT(selectAll()));
    connect(ui.unselectAllButton, SIGNAL(clicked()), this, SLOT(unselectAll()));


    fillTreeWidgetWithAttrMaps(m,ui.attrMapTree);

    if(stringToBool(parm("Mesh Type"))) return true;

    if(dlg.exec() == QDialog::Accepted){
      QStringList newList;
      for(int itemCount = 0; itemCount < ui.attrMapTree->topLevelItemCount(); ++itemCount){
        QString currentAttr = ui.attrMapTree->topLevelItem(itemCount)->text(1);
        if(Qt::Checked == ui.attrMapTree->topLevelItem(itemCount)->checkState(0)){
          newList << currentAttr;
        }
      }
      attrMapsToBeSaved = newList;

    }

    return true;
  }




  void createPlyDataFromAttr(Mesh* m, QStringList attrMapsToBeSaved, std::vector<PlyCellData>& cellVec, std::vector<PlyWallData>& wallVec)
  {

    forall(QString attr, attrMapsToBeSaved){
      PlyCellData cellData;
      PlyWallData wallData;
      QStringList parts = attr.split(" ");

      //qDebug() << "attr " << attr << "\n";

      QString measure = "Measure ";

      //if(parts[0] != "Measure") continue; //return setErrorMessage("not a measure");
      if(parts[0] == "Label"){

        cellData.name = parts[2];
        for(int i=4; i<parts.size(); i++) cellData.name += " " + parts[i];

        if(parts[1] == "Double"){
          AttrMap<int, double>& currentMap = m->attributes().attrMap<int, double>(measure + attr);
          qDebug() << attr;
          //std::cout << " attr " << currentMap.size() << std::endl;
          forall(IntDoublePair p, currentMap){
            //std::cout << "p " << p.first << "/" << p.second << std::endl;
            cellData.floatData[p.first] = p.second;
            cellData.type = "float";
          }
        } else if(parts[1] == "Int"){
          AttrMap<int, int>& currentMap = m->attributes().attrMap<int, int>(measure + attr);
          qDebug() << attr;
          //std::cout << " attr " << currentMap.size() << std::endl;
          forall(IntIntPair p, currentMap){
            cellData.intData[p.first] = p.second;
            cellData.type = "float";
          }
        } else if(parts[1] == "Vector"){
          AttrMap<int, Point3d>& currentMap = m->attributes().attrMap<int, Point3d>(measure + attr);
          qDebug() << attr;
          //std::cout << " attr " << currentMap.size() << std::endl;
          forall(IntPoint3dPair p, currentMap){
            cellData.pointData[p.first] = Point3f(p.second);
            cellData.type = "vector";
          }
        } else if(parts[1] == "Tensor"){
          AttrMap<int, SymmetricTensor>& currentMap = m->attributes().attrMap<int, SymmetricTensor>(measure + attr);
          qDebug() << attr;
          //std::cout << " attr " << currentMap.size() << std::endl;
          typedef std::pair<int, SymmetricTensor> IntSymmetricTensorPair;
          forall(IntSymmetricTensorPair p, currentMap){
            cellData.tensorData[p.first] = p.second;
            cellData.type = "tensor";
          }

        } else {
          continue; //return setErrorMessage("wrong type");
        }
        cellVec.push_back(cellData);
      } else if(parts[0] == "Wall"){
      } else {
        continue; //return setErrorMessage("not a label or wall map");
      }
    }


  }


  bool SavePlyFile::run(Mesh* m, QString filename, QString meshType)
  {
    //if(m->meshType()!="MGX2D" and m->meshType()!="MGX3D") return setErrorMessage("Mesh Type has to be MGX2D or MGX3D");

    // save mesh information to two ply files

    //std::cout << "lets save the ply" << std::endl;

    // 1st file standard mesh exported by mgx
    MeshExport me(*this);
    me.savePLY(m, filename, "PLY Ascii", false); // TODO make these parameters?

    // 2nd file contains information about cells (cell graph) and walls

    filename.remove(".ply", Qt::CaseInsensitive);
    QString fileCell = filename + "_cellGraph.ply";

    //std::cout << "create cell graph data" << std::endl;

    // if(marionsFile){
    //   savePLYCellGraphMarion(m, fileCell, false);
    // } else {
    std::vector<PlyCellData> cellVec;
    std::vector<PlyWallData> wallVec;

    createPlyDataFromAttr(m, attrMapsToBeSaved, cellVec, wallVec);

    //std::cout << "now cell graph ply" << std::endl;

      if(!savePLYCellGraphGeneral2D(m,cellVec,wallVec,fileCell,false)) return setErrorMessage("Error writing the file!");
    // }

    return true;
  }
  REGISTER_PROCESS(SavePlyFile);

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


  bool LoadPlyFileCellGraph::run(Mesh* m, const QString& filename)
  {

    // save mesh information to two ply files

    // 1st file standard mesh exported by mgx
    /*   MeshExport me(*this);
    me.savePLY(m, filename, false, false); // TODO make these parameters?

    // 2nd file contains information about cells (cell graph) and walls

    QString fileCell = filename + "_cellGraph.ply";

    savePLYCellGraph(m,fileCell, false);
*/
    loadCellGraphPLYGeneral(m, filename);

    return true;
  }
  REGISTER_PROCESS(LoadPlyFileCellGraph);

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


  bool SaveCellContourFile::run(Mesh* m, const QString& filename)
  {

    if(m->meshType() != "MGXC" and m->meshType() != "MGXM")
      return setErrorMessage("MGXM or MGXC Mesh required.");

    // create cell mesh (only 1 vtx per cell)
    vvGraph T;
    vvGraph &S = m->graph();

    // find selected cells
    std::set<int> selectedCells = findAllSelectedLabels(S);

    if(selectedCells.empty())
      return setErrorMessage("No cells selected.");

    if(m->meshType() == "MGXM")
      mgxmToMgxc(S, T, -1);
    else
      T = S;

    std::map<int, std::vector<Point2d> > labelContourMap;


    // cell normals
    std::map<int, Point3d> cellNormals;
    vvGraph T2;
    // extract border vtxs of the selected cells - in order!
    forall(const vertex& v, T){
      if(v->label < 1) continue;
      if(selectedCells.find(v->label) == selectedCells.end()) continue;

      cellNormals[v->label] = v->nrml;

      Point3d planeNrml = v->nrml;
      double counter = 1;
      std::vector<Point3d> pointsPolygon;
      forall(const vertex &n, T.neighbors(v)) {
        counter++;
        pointsPolygon.push_back(n->pos);
      }
      //pointsPolygon = labelPolygonMap[v->label];

      //if(pcaNrmls){ // generate the cell polygon plane with a pca
      Point3d planePos, planeNrml2;
      findPolygonPlane(pointsPolygon, planePos, planeNrml2);
      if(planeNrml*planeNrml2 < 0)
        planeNrml = planeNrml2;
      else
        planeNrml = -planeNrml2;
      //}

      // now project the cell vtxs onto the cell plane
      projectPointsOnPlane(pointsPolygon, planePos, planeNrml);

      // rotate the cell plane into XY
      Matrix3d rotMat;
      rotatePointsIntoXY(planeNrml, pointsPolygon, rotMat);

      Point2d avgPos(0,0);
      forall(Point3d& p, pointsPolygon){
        //if(norm(p.z()) > 0.00001) std::cout << "sth wrong here" << std::endl;
        Point2d p2d(p.x(), p.y());
        avgPos+=p2d;
      }
      avgPos /= pointsPolygon.size();
      forall(Point3d& p, pointsPolygon){
        //if(norm(p.z()) > 0.00001) std::cout << "sth wrong here" << std::endl;
        Point2d p2d(p.x(), p.y());
        labelContourMap[v->label].push_back(p2d-avgPos);
        //vertex v2;
        //v2->pos = p;
        //T2.insert(v2);
      }
    }

    //S=T2;
    //m->updateAll();
    // save coords to file

    std::cout << "save txt " << labelContourMap.size() << std::endl;
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("CellContour::Error:Cannot open output file: %1").arg(filename));
      return false;
    }
    QTextStream outText(&file);

    typedef std::pair<int, std::vector<Point2d> > IntVecP2DPair;
    forall(IntVecP2DPair p, labelContourMap){
      outText << p.first << "," << m->parents()[p.first];
      forall(Point2d c, p.second){
        outText << "," << c.x() << "," << c.y();
      }
      outText << endl;
    }



    return true;
  }
  REGISTER_PROCESS(SaveCellContourFile);

  std::vector<Point3d> mergeClosePoints(std::vector<Point3d>& vIn, double threshold)
  {
    std::vector<Point3d> vOut;

    if(vIn.empty()) return vOut;

    Point3d last = vIn[vIn.size()-1];

    for(uint i=0;i<vIn.size();i++){
      double dis = norm(last - vIn[i]);
      if(dis > threshold){
        vOut.push_back(last);
        last = vIn[i];
      } else {
        Point3d avg = (last + vIn[i])/2.;
        vOut.push_back(avg);
      }
    }

    return vOut;
  }



   // file dialogue for save data
   bool ExportInternalAngles::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;
    }


    // TODO Test process
    bool ExportInternalAngles::run(Mesh* m, QString filename, double threshold)
    {

      if(m->meshType() != "MGXC")
        return setErrorMessage("MGXC Mesh required.");


    vvGraph& S = m->graph();

    std::map<int, std::vector<double> > labelAnglesMap;

    forall(const vertex& v, S){
      if(v->label < 1) continue;

      std::vector<Point3d> cellContour, cellContourMerged;
      forall(const vertex& n, S.neighbors(v)){
        cellContour.push_back(n->pos);
      }
      cellContourMerged = mergeClosePoints(cellContour,threshold);

      if(cellContourMerged.size() < 3) continue;

      Point3d last, current;
      last = cellContourMerged[cellContourMerged.size()-2];
      current = cellContourMerged[cellContourMerged.size()-1];
      for(uint i=0;i<cellContourMerged.size();i++){
        Point3d next = cellContourMerged[i];
        Point3d vec1, vec2;
        vec1 = next - current;
        vec2 = last - current;
        double a = angleVectors(vec1, vec2, true);
        labelAnglesMap[v->label].push_back(a);
        last = current;
        current = next;
      }



    }

    // save result to file

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

    out << "label, corners, angles" << endl;

    forall(auto p, labelAnglesMap){
      out << p.first << "," << p.second.size() << ",";
      forall(auto values, p.second)
        out << values << ",";

      out << endl;
    }


    return true;

    }

    REGISTER_PROCESS(ExportInternalAngles);



  double angleCellMeshJunction(vertex j, std::vector<vertex>& neighborVtxs, int vectorPosCell, std::map<vertex, std::vector<vertex> >& neighborVtxsNeighbors, double thresholdDis, QString method)
  {

    Point3d vec1, vec2;

    int upper = vectorPosCell+1;
    if(upper >= neighborVtxs.size()) upper = 0;
    vertex cellJ1 = neighborVtxs[upper];

    int lower = vectorPosCell-1;
    if(lower < 0) lower = neighborVtxs.size()-1;
    vertex cellJ2 = neighborVtxs[lower];

    Point3d cell1Dir, cell2Dir;
    Point3d cell1DirFirst, cell2DirFirst;
    Point3d cell1DirLast, cell2DirLast;

    cell1DirFirst = neighborVtxsNeighbors[cellJ1][0]->pos;
    cell1DirLast = neighborVtxsNeighbors[cellJ1][neighborVtxsNeighbors[cellJ1].size()-1]->pos;

    cell2DirFirst = neighborVtxsNeighbors[cellJ2][0]->pos;
    cell2DirLast = neighborVtxsNeighbors[cellJ2][neighborVtxsNeighbors[cellJ2].size()-1]->pos;

    // check neighbors
    int c1 = 0;
    for(uint i=0; i<neighborVtxsNeighbors[cellJ1].size();i++){
      if(norm(neighborVtxsNeighbors[cellJ1][i]->pos - j->pos) < thresholdDis){
        cell1Dir+= neighborVtxsNeighbors[cellJ1][i]->pos;
        c1++;
        //neighborVtxsNeighbors[cellJ1][i]->selected = true;
      }
    }
    cell1Dir+=cellJ1->pos;
    cell1Dir/=(c1+1);

    int c2 = 0;
    for(uint i=0; i<neighborVtxsNeighbors[cellJ2].size();i++){
      if(norm(neighborVtxsNeighbors[cellJ2][i]->pos - j->pos) < thresholdDis){
        cell2Dir+= neighborVtxsNeighbors[cellJ2][i]->pos;
        c2++;
        //neighborVtxsNeighbors[cellJ2][i]->selected = true;
      }
    }
    cell2Dir+=cellJ2->pos;
    cell2Dir/=(c2+1);

    if(method == "Average"){
      vec1 = cell1Dir - j->pos;
      vec2 = cell2Dir - j->pos;
    } else if(method == "Closest"){
      vec1 = cell1DirFirst - j->pos;
      vec2 = cell2DirFirst - j->pos;
    } else { // Furthest
      vec1 = cell1DirLast - j->pos;
      vec2 = cell2DirLast - j->pos;
    }

    return angleVectors(vec1, vec2, true);
  }

  // helper function to calculate internal cell angles at junctions
  std::vector<std::pair<int, double> > cellAngles(vertex j, std::vector<vertex>& neighborVtxs, std::map<vertex, std::vector<vertex> >& neighborVtxsNeighbors, double thresholdDis, QString method)
  {

    std::vector<std::pair<int, double> > result;

    if(neighborVtxs.size() == 6){ // 3 way junction
      vertex cellV1 = neighborVtxs[0];
      int start = 0;
      if(cellV1->label == -1){
        cellV1 = neighborVtxs[1];
        start = 1;
      }
      vertex cellV2 = neighborVtxs[2+start];
      vertex cellV3 = neighborVtxs[4+start];

      double angle1 = angleCellMeshJunction(j, neighborVtxs, start, neighborVtxsNeighbors, thresholdDis, method);
      double angle2 = angleCellMeshJunction(j, neighborVtxs, start+2, neighborVtxsNeighbors, thresholdDis, method);
      double angle3 = angleCellMeshJunction(j, neighborVtxs, start+4, neighborVtxsNeighbors, thresholdDis, method);

      result.push_back(std::make_pair(cellV1->label,angle1));
      result.push_back(std::make_pair(cellV2->label,angle2));
      result.push_back(std::make_pair(cellV3->label,angle3));


    } else if(neighborVtxs.size() == 8){ // 4 way junction
      vertex cellV1 = neighborVtxs[0];
      int start = 0;
      if(cellV1->label == -1){
        cellV1 = neighborVtxs[1];
        start = 1;
      }
      vertex cellV2 = neighborVtxs[2+start];
      vertex cellV3 = neighborVtxs[4+start];
      vertex cellV4 = neighborVtxs[6+start];

      double angle1 = angleCellMeshJunction(j, neighborVtxs, start, neighborVtxsNeighbors, thresholdDis, method);
      double angle2 = angleCellMeshJunction(j, neighborVtxs, start+2, neighborVtxsNeighbors, thresholdDis, method);
      double angle3 = angleCellMeshJunction(j, neighborVtxs, start+4, neighborVtxsNeighbors, thresholdDis, method);
      double angle4 = angleCellMeshJunction(j, neighborVtxs, start+6, neighborVtxsNeighbors, thresholdDis, method);

      result.push_back(std::make_pair(cellV1->label,angle1));
      result.push_back(std::make_pair(cellV2->label,angle2));
      result.push_back(std::make_pair(cellV3->label,angle3));
      result.push_back(std::make_pair(cellV4->label,angle4));
    } else { // do nothing

    }


    return result;
  }



   // file dialogue for save data
   bool ExportJunctionAngles::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 ExportJunctionAngles::run(Mesh* m, Mesh* m2, QString filename, double thresholdDis, QString method)
    {

      if(m->meshType() != "MGXC")
        return setErrorMessage("MGXC Mesh required.");


    vvGraph& S = m->graph();

    JunctionInformation jInfo;
    findCellJunctions(m, S, jInfo);

    JunctionCellsMap jCells = jInfo.junctionLabels1;

    std::map<std::pair<vertex, int>, double> junctionCellPairAnglesMap;


    forall(auto p, jCells){
      vertex j = p.first;
      //std::cout << "j " << j->saveId << "/" << p.second.size() << std::endl;
      // ignore junctions with fewer than 3 or more than 4 cells
      if(p.second.size() < 3 or p.second.size() > 4) continue;

      j->selected = true;

      std::vector<vertex> neighborVtxs;

      std::map<vertex, std::vector<vertex> > neighborVtxsNeighbors;

      // find neighbor vertices of junctions and follow the cell border for up to 4 vertices
      forall(vertex n, S.neighbors(j)){
        //if(n->label > -1) continue;
        neighborVtxs.push_back(n);
        //n->selected=true;
        forall(vertex n2, S.neighbors(n)){
          if(n2->label > -1) continue;
          if(n2 == j) continue;
          neighborVtxsNeighbors[n].push_back(n2);
          //n2->selected=true;
          forall(vertex n3, S.neighbors(n2)){
            if(n3->label > -1) continue;
            if(n3 == n) continue;
            if(jCells.find(n3) != jCells.end()) continue;
            neighborVtxsNeighbors[n].push_back(n3);
            forall(vertex n4, S.neighbors(n3)){
              if(n4->label > -1) continue;
              if(n4 == n2) continue;
              if(jCells.find(n4) != jCells.end()) continue;
              neighborVtxsNeighbors[n].push_back(n4);
            }
          }
        }
      }

      std::vector<std::pair<int, double> > r;
      r = cellAngles(j, neighborVtxs, neighborVtxsNeighbors, thresholdDis, method);


      std::cout << "j " << j->saveId << "/" << r.size() << std::endl;


      forall(auto v, r){
        std::pair<vertex, int> p = std::make_pair(j, v.first);
        junctionCellPairAnglesMap[p] = v.second;
        std::cout << "r " << v.first << "/" << v.second << std::endl;
      }
    }

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

    out << "junction, cell, angle" << endl;

    forall(auto p, junctionCellPairAnglesMap){
      out << p.first.first->saveId << "," << p.first.second << "," << p.second;
      out << endl;
    }


    m->updateAll();

    return true;

    }

    REGISTER_PROCESS(ExportJunctionAngles);


  std::set<int> commonElements(const std::set<int>& s1, const std::set<int>& s2)
  {
    std::set<int> res;

    forall(int el, s1){
      if(s2.find(el) == s2.end()) continue;
      res.insert(el);
    }
    
    return res;
  }

  // helper function to calculate internal cell angles at junctions
  std::vector<std::pair<int, double> > cellAnglesAtJunction(Mesh* m, vvGraph& S, std::vector<std::vector<vertex> >& allBordersJunction, 
    double thresholdDis, QString method, bool useParentLabels)
  {
    std::vector<std::pair<int, double> > cellAngles;

    if(allBordersJunction.empty()) return cellAngles;
    if(allBordersJunction[0].empty()) return cellAngles;

    // cells of border segments
    std::map<int, std::set<int> > borderSegmentCells;

    // endpoints of the segments
    std::map<int, Point3d> borderSegmentEndpoint;

    vertex centerJunctionVtx = allBordersJunction[0][0];
    centerJunctionVtx->selected = true;

    // find neighbor cells of center junction
    std::set<int> centerJunctionNeighbors;
    forall(vertex n, S.neighbors(centerJunctionVtx)){
      if(n->label == -1) continue;
      if(useParentLabels) centerJunctionNeighbors.insert(m->parents()[n->label]);
      else centerJunctionNeighbors.insert(n->label);
    }

    for(int i = 0; i < allBordersJunction.size(); i++){
      std::vector<vertex> vec = allBordersJunction[i];

      // take the 2nd vertex and find its neighbors
      vertex vertex2 = allBordersJunction[i][1];
      std::set<int> vertex2Neighbors;
      forall(vertex n, S.neighbors(vertex2)){
        if(n->label == -1) continue;
        if(useParentLabels) vertex2Neighbors.insert(m->parents()[n->label]);
        else vertex2Neighbors.insert(n->label);
      }

      // find common neighbors of center junction and 2nd junction -> the wall is between the resulting two cells
      std::set<int> cellsSegment = commonElements(centerJunctionNeighbors, vertex2Neighbors);

      if(cellsSegment.size() != 2) continue;
      borderSegmentCells[i] = cellsSegment;

      double totalDis = 0;
      Point3d pointAvg(0,0,0);
      int totalCount = 0;
      vertex last = centerJunctionVtx;

      for(int j = 1; j < allBordersJunction[i].size(); j++){
        if(totalDis > thresholdDis) continue;
        vertex lastNew = allBordersJunction[i][j];
        totalDis += norm(last->pos - lastNew->pos);
        pointAvg += lastNew->pos;
        totalCount++;
        last = lastNew;
      }

      Point3d finalPoint(0,0,0);

      if(method == "Closest"){
        finalPoint = allBordersJunction[i][1]->pos;
      } else if(method == "Furthest"){
        finalPoint = last->pos;
      } else { // Average
        finalPoint = pointAvg / totalCount;
      }

      borderSegmentEndpoint[i] = finalPoint;

    }

    // we already have foreach bordersegment the endpoint and the neighboring cells
    // next we compute the angles between neighboring segments

    for(int i = 0; i < allBordersJunction.size(); i++){
      for(int j = 0; j < allBordersJunction.size(); j++){
        if(i >= j) continue;

        std::set<int> sharedCell = commonElements(borderSegmentCells[i],borderSegmentCells[j]);
        if(sharedCell.empty()) continue;
        if(sharedCell.size() > 1) continue;

        // the two border segments share exactly 1 cell
        int cellLabel = *(sharedCell.begin());

        // compute the angle
        Point3d vec1 = borderSegmentEndpoint[i] - centerJunctionVtx->pos;
        Point3d vec2 = borderSegmentEndpoint[j] - centerJunctionVtx->pos;

        //std::cout << "endPs " << borderSegmentEndpoint[i] << "/" << borderSegmentEndpoint[j] << std::endl;

        std::pair<int, double> cellAngleP = std::make_pair(cellLabel,angleVectors(vec1, vec2, true));

        cellAngles.push_back(cellAngleP);
        //std::cout << "cell angles " << cellAngleP.first << "/" << cellAngleP.second << std::endl;
      }
    }

    return cellAngles;
  }
  

   QString createJunctionString(std::vector<int>& centerJunctionNeighbors)
   {
    QString junctionLabel;

      std::sort(centerJunctionNeighbors.begin(),centerJunctionNeighbors.end());

      // create the junction label
      junctionLabel = QString::number(centerJunctionNeighbors[0]);
      if(centerJunctionNeighbors.size() > 1){
        junctionLabel = junctionLabel + "x";
        junctionLabel = junctionLabel + QString::number(centerJunctionNeighbors[1]);
      } 
      if(centerJunctionNeighbors.size() > 2){
        junctionLabel = junctionLabel + "x";
        junctionLabel = junctionLabel + QString::number(centerJunctionNeighbors[2]);
      } 
      if(centerJunctionNeighbors.size() > 3){
        junctionLabel = junctionLabel + "x";
        junctionLabel = junctionLabel + QString::number(centerJunctionNeighbors[3]);
      } 

    return junctionLabel;
   }

   QString convertJunctionStringLabelToParent(QString junctionString, Mesh* m)
   {
    QStringList pieces = junctionString.split("x");
    std::vector<int> parentLabels;
    forall(QString p, pieces){
      int label = p.toInt();
      parentLabels.push_back(m->parents()[label]);
    }
    return createJunctionString(parentLabels);
   }


   // file dialogue for save data
   bool ExportJunctionAnglesNew::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;
    }

  
  std::map<std::pair<QString, int>, double> createJunctionAngleMap(Mesh* m, vvGraph& S, JunctionCellsMap jCells, 
    double thresholdDis, QString method, bool useParentLabels)
    {
      std::map<std::pair<QString, int>, double> junctionCellPairAnglesMap;

    forall(auto p, jCells){
      vertex j = p.first;
      //std::cout << "j " << j->saveId << "/" << p.second.size() << std::endl;
      // ignore junctions with fewer than 3 or more than 4 cells
      if(p.second.size() < 3 or p.second.size() > 4) continue;

      //j->selected = true;

      std::vector<std::vector<vertex> > allBordersJunction;

      forall(vertex n, S.neighbors(j)){
        if(n->label != -1) continue;
        std::vector<vertex> borderSegment;
        borderSegment.push_back(j);
        borderSegment.push_back(n);
        // now follow the border until the next junction is hit
        vertex currentBorderVtx = n;
        vertex lastBorderVtx = j;
        while(jCells.find(currentBorderVtx) == jCells.end()){

          vertex k;
          k->label = -2;

          forall(vertex m, S.neighbors(currentBorderVtx)){
            if(m->label != -1) continue;
            if(m == lastBorderVtx) continue;
            if(k->label!=-2){
              std::cout << "sth weird happened 1" << std::endl;
            }
            k=m;
          }

          if(k==0){
            std::cout << "sth weird happened 2" << std::endl;
            break;
          } else {
            borderSegment.push_back(k);
            lastBorderVtx = currentBorderVtx;
            currentBorderVtx = k;
          }
        } // while
        allBordersJunction.push_back(borderSegment);
      }

      //std::cout << "j " << j->saveId << "/" << allBordersJunction.size() << std::endl;
      std::vector<std::pair<int, double> > r;
      r = cellAnglesAtJunction(m, S, allBordersJunction, thresholdDis, method, useParentLabels);

      // create unique label for junction based on cell labels

      // find neighbor cells of center junction
      std::vector<int> centerJunctionNeighbors;
      forall(vertex n, S.neighbors(j)){
        if(n->label == -1) continue;
        if(useParentLabels) centerJunctionNeighbors.push_back(m->parents()[n->label]);
        else centerJunctionNeighbors.push_back(n->label);
        
      }

      // create the junction label
      QString junctionLabel = createJunctionString(centerJunctionNeighbors);//QString::number(centerJunctionNeighbors[0]);

      // write output map
      forall(auto v, r){
        std::pair<QString, int> p = std::make_pair(junctionLabel, v.first);
        junctionCellPairAnglesMap[p] = v.second;
        //std::cout << "r " << v.first << "/" << v.second << std::endl;
      }


    }


      return junctionCellPairAnglesMap;
    }



    bool ExportJunctionAnglesNew::run(Mesh* m, Mesh* m2, QString filename, double thresholdDis, QString method, bool changeMap)
    {

      if(m->meshType() != "MGXC")
        return setErrorMessage("MGXC Mesh required in mesh1.");


    vvGraph& S = m->graph();

    JunctionInformation jInfo;
    findCellJunctions(m, S, jInfo);

    JunctionCellsMap jCells = jInfo.junctionLabels1;

    std::map<std::pair<QString, int>, double> junctionCellPairAnglesMap = createJunctionAngleMap(m, S, jCells, thresholdDis, method, m->useParents());
    std::map<std::pair<QString, int>, double> junctionCellPairAnglesMap2, junctionCellPairAnglesMap2_Labels;
    if(changeMap){

      if(m2->meshType() != "MGXC")
        return setErrorMessage("MGXC Mesh required in mesh2.");
      vvGraph& S2 = m2->graph();

      JunctionInformation jInfo2;
      findCellJunctions(m2, S2, jInfo2);

      JunctionCellsMap jCells2 = jInfo2.junctionLabels1;

      //junctionCellPairAnglesMap2 = createJunctionAngleMap(m2, S2, jCells2, thresholdDis, method, m2->useParents());
      //if(m2->useParents()) // also get the label map, there might be additional junctions from divided cells
        junctionCellPairAnglesMap2_Labels = createJunctionAngleMap(m2, S2, jCells2, thresholdDis, method, false);
    }

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

    out << "junction,cell,angle";
    if(changeMap) out << ",junction2,cell2,angle2,angleRatio";
    out << endl;

    std::set<std::pair<QString, int> > finishedJunctions;

    forall(auto p, junctionCellPairAnglesMap){
      finishedJunctions.insert(p.first);
      out << p.first.first << "," << p.first.second << "," << p.second;
      if(changeMap){
        bool m2EntryFound = false;

        if(m2->useParents()){ // search for a fitting entry
          std::pair<QString, int> map2_Pair;
          forall(auto p2, junctionCellPairAnglesMap2_Labels){
            map2_Pair = p2.first;
            QString convertedString = convertJunctionStringLabelToParent(p2.first.first, m2);
            std::pair<QString, int> newP = std::make_pair(convertedString,m2->parents()[p2.first.second]);
            if(p.first == newP){
              m2EntryFound = true;
              break;
            }
          }
          if(m2EntryFound){
            if(p.second == 0){
              out << "," << map2_Pair.first << "," << map2_Pair.second << "," << junctionCellPairAnglesMap2_Labels[map2_Pair] << ",NA";
            } else {
              out << "," << map2_Pair.first << "," << map2_Pair.second << "," << junctionCellPairAnglesMap2_Labels[map2_Pair] << "," << junctionCellPairAnglesMap2_Labels[map2_Pair]/p.second;
            }
          } else {
            out << ",NA,NA,NA,NA";
          }

        } else {// look if entry exists
          if(junctionCellPairAnglesMap2_Labels.find(p.first) == junctionCellPairAnglesMap2_Labels.end()){
            out << ",NA,NA,NA,NA";
          } else {
            if(p.second == 0){
              out << ",NA,NA," << junctionCellPairAnglesMap2_Labels[p.first] << ",NA";
            } else {
              out << ",NA,NA," << junctionCellPairAnglesMap2_Labels[p.first] << "," << junctionCellPairAnglesMap2_Labels[p.first]/p.second;
            }
          }
        }
        
      } 
      out << endl;
    }

    if(changeMap){

      forall(auto p, junctionCellPairAnglesMap2_Labels){
        std::pair<QString, int> newP;
        if(m2->useParents()){
          QString convertedString = convertJunctionStringLabelToParent(p.first.first, m2);
          newP = std::make_pair(convertedString,m2->parents()[p.first.second]);
        } else {
          newP = p.first;
        }
        if(finishedJunctions.find(newP) != finishedJunctions.end()) continue;
        out <<  "NA,NA,NA," << p.first.first << "," << p.first.second << "," << p.second << ",NA";
        out << endl;
      }

      
    }

    m->updateAll();
    return true;
  }

  REGISTER_PROCESS(ExportJunctionAnglesNew);
}
