//
// 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 <SystemProcessSave.hpp>

#include <QAbstractItemModel>
#include <QComboBox>
#include <QInputDialog>
#include <QItemDelegate>
#include <QRegExp>
#include <QtXml>

#include <CImg.h>

#include <Triangulate.hpp> // findPolygonPlane
#include <Dir.hpp>
#include <Image.hpp>
#include <PlyFile.hpp>
#include <Version.hpp>

#include <SystemProcessLoad.hpp> // Batch compress process

#include <ui_SaveMesh.h>
#include <ui_ExportMesh.h>
#include <ui_ExportStack.h>
#include <ui_PlyCellGraphDlg.h>

// #include "H5Cpp.h" HDF5 now in addon

using namespace cimg_library;
using namespace qglviewer;
using namespace mgx;
// using namespace H5; HDF5 now in addon


#ifdef WIN32
#  define FileDialogOptions QFileDialog::DontUseNativeDialog
#else
#  define FileDialogOptions 0
#endif

namespace mgx 
{
  typedef CImg<ushort> CImgUS;
  
  struct StackHdr
  {
    uint hdrsz;
    uint xsz, ysz, zsz;
    float xum, yum, zum;
    uint datasz;

    StackHdr() : hdrsz(1024) {}
  };
  
  bool StackSave::run()
  {
    int stackId = parm("Stack Number").toInt();
    int compressionLevel = parm("Compression Level").toInt();
    if(!checkState().stack(STACK_ANY, stackId))
      return false;
    Stack* stack = this->stack(stackId);
    Store* store = (stringToWorkStore(parm("Store")) ? stack->work() : stack->main());
    bool res = run(stack, store, parm("Filename"), compressionLevel);
    return res;
  }
  
  static const char *const ELEMENT_SIZE_UM = "element_size_um";

  bool StackSave::run(Stack* stack, Store* store, const QString& filename, int compressionLevel)
  {
    if(filename.isEmpty()) {
      setErrorMessage("Error, trying to save a stack with an empty filename.");
      return false;
    }

    actingFile(filename);

    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly))
      throw QString("Cannot open output file - %1").arg(filename);

    store->setFile(stripCurrentDir(filename));

    HVecUS& data = store->data();

    progressStart(QString("Saving %1 Stack %2").arg(store == stack->main() ? "main" : "work").arg(stack->userId()),
                  0);
    qint64 qty;

    // Process the various types
    if(filename.endsWith(".mgxs", Qt::CaseInsensitive)) {
      if(compressionLevel > 9)
        compressionLevel = 9;
      if(compressionLevel < -1)
        compressionLevel = -1;
      Point3u size = stack->size();
      Point3f step = stack->step();
      Point3f origin = stack->origin();
      quint32 xsz = size.x();
      quint32 ysz = size.y();
      quint32 zsz = size.z();
      quint8 cl = (quint8)compressionLevel;
      float xum = step.x();
      float yum = step.y();
      float zum = step.z();
      float sxum = origin.x();
      float syum = origin.y();
      float szum = origin.z();
      bool isLabel = store->labels();
      quint64 datasz = quint64(data.size()) * 2ul;

      static const char version[] = "MGXS 1.3 ";
      qty = file.write(version, sizeof(version) - 1);
      if(qty != sizeof(version) - 1)
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(version)
                                                                                                 - 1);
      qty = file.write((const char*)&isLabel, sizeof(bool));
      if(qty != sizeof(bool))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(bool));
      qty = file.write((const char*)&sxum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&syum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&szum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&xsz, sizeof(quint32));
      if(qty != sizeof(quint32))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(quint32));
      qty = file.write((const char*)&ysz, sizeof(quint32));
      if(qty != sizeof(quint32))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(quint32));
      qty = file.write((const char*)&zsz, sizeof(quint32));
      if(qty != sizeof(quint32))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(quint32));
      qty = file.write((const char*)&xum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&yum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&zum, sizeof(float));
      if(qty != sizeof(float))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(float));
      qty = file.write((const char*)&datasz, sizeof(quint64));
      if(qty != sizeof(quint64))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(quint64));
      qty = file.write((const char*)&cl, sizeof(quint8));
      if(qty != sizeof(quint8))
        throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(sizeof(quint8));

      if(compressionLevel > 0) {
        // QByteArray stkdata = qCompress((const uchar *)data.data(), data.size() * 2ul, 1);
        // file.write(stkdata);
        // Make slices of 64MB
        quint64 nb_slices = datasz >> 26;
        quint64 slice_size = 1ul << 26;
        progressStart("Writing Slices", nb_slices);
        for(quint64 i = 0; i < nb_slices; ++i) {
          QByteArray stkdata
              = qCompress((const uchar*)data.data() + i * slice_size, slice_size, compressionLevel);
          quint32 sz = stkdata.size();
          qty = file.write((const char*)&sz, sizeof(quint32));
          if(qty != sizeof(quint32)) {
            throw QString("Could not write enough data to file (written = %1, expected = %2)").arg(qty).arg(
                  sizeof(quint32));
          }
          qty = file.write(stkdata);
          if(qty != stkdata.size()) {
            throw QString("Could not write enough data to file (written = %1, expected = %2)").arg(qty).arg(
                  stkdata.size());
          }
          if(!progressAdvance(i))
            userCancel();
        }
        quint64 left_size = datasz - (nb_slices << 26);
        QByteArray stkdata = qCompress((const uchar*)data.data() + (nb_slices << 26), left_size, compressionLevel);
        quint32 sz = stkdata.size();
        qint64 qty = file.write((const char*)&sz, sizeof(quint32));
        if(qty != sizeof(quint32)) {
          throw QString("Could not write enough data to file (written = %1, expected = %2)").arg(qty).arg(
                sizeof(quint32));
        }
        qty = file.write(stkdata);
        if(qty != stkdata.size()) {
          throw QString("Could not write enough data from MGXS file (written = %1, expected = %2)").arg(qty).arg(
                stkdata.size());
        }
        if(!progressAdvance(nb_slices))
          userCancel();
      } else {
        // Ok, if greater than 64MB, write in bits
        quint64 nb_slices = datasz >> 26;
        quint64 slice_size = 1ul << 26;
        for(quint64 i = 0; i < nb_slices; ++i) {
          qty = file.write((char*)data.data() + i * slice_size, slice_size);
          if(quint64(qty) != slice_size) {
            throw QString("Could not write enough data to file (written = %1, expected = %2)").arg(qty).arg(
                  slice_size);
          }
        }
        quint64 left_size = datasz - (nb_slices << 26);
        qty = file.write((char*)data.data() + (nb_slices << 26), left_size);
        if(quint64(qty) != left_size) {
          throw QString("Could not write enough data from MGXS file (written = %1, expected = %2)").arg(qty).arg(
                left_size);
        }
      }
      progressAdvance(0);

      file.close();
    } else if(filename.endsWith(".inr", Qt::CaseInsensitive)) {
      file.close();
      Point3u size = stack->size();
      Point3f step = stack->step();
      CImgUS image(&data[0], size.x(), size.y(), size.z(), 1, true);
      image.save_inr(filename.toStdString().data(), step.data());
    } else if(filename.endsWith(".tif", Qt::CaseInsensitive) or filename.endsWith(".tiff", Qt::CaseInsensitive)) {
      file.close();
      Point3u size = stack->size();
      Point3f step = stack->step();
      Image3D image(data, size, step, store->labels());
      saveTIFFImage(filename, image);
    // } else if(filename.endsWith(".h5", Qt::CaseInsensitive) or filename.endsWith(".hdf5", Qt::CaseInsensitive) or filename.endsWith(".hdf", Qt::CaseInsensitive)) { // h5 stack save
    //   file.close();
    //   Point3u size = stack->size();
    //   Point3f step = stack->step();
    //   Image3D image(data, size, step, store->labels());

    //   std::string filenameString = filename.toStdString();
    //   const char* fileNameH5 = filenameString.c_str();

    //   QString qDataSetName = parm("HDF5 DataSetName");
    //   std::string datasetNameString = qDataSetName.toStdString();
    //   const char* datasetName = datasetNameString.c_str();

    //   // qDebug() << qDataSetName;
    //   // std::cout << __LINE__ << "/" << datasetName << std::endl;

    //   if (qDataSetName.isEmpty()) {
    //       throw QString("StackSave::Error: H5 dataset name cannot be empty");
    //   }

    //   // Create HDF5 file and dataset
    //   H5File file(fileNameH5, H5F_ACC_TRUNC);
    //   // hardcode dataset dim
    //   const int rank = 3;
    //   hsize_t* dims = new hsize_t[rank];
    //   hsize_t dims_1d = 1;
    //   for(int i=0; i < rank; i++) {
    //       dims[rank-i-1] = size[i];
    //       dims_1d *= size[i];
    //   }
    //   const hsize_t* rankHdf5 = new hsize_t(rank);

    //   DSetCreatPropList  *plist = new DSetCreatPropList;
    //   hsize_t chunk_dims[3] = { 20, 20, 20 };
    //   plist->setChunk(3, chunk_dims);
    //   plist->setDeflate(compressionLevel);

    //   DataSpace dataspace(rank, dims);
    //   DataSet dataset = file.createDataSet(datasetName, PredType::NATIVE_UINT16, dataspace, *plist);

    //   // Create a dataset attribute. 
    //   DataSpace attr_dataspace = DataSpace (1, rankHdf5);
    //   Attribute attribute = dataset.createAttribute(ELEMENT_SIZE_UM, PredType::NATIVE_FLOAT, attr_dataspace);

    //   // Write the attribute data. 
    //   float* element_size_um = new float[3];
    //   element_size_um[2] = stack->step().x();
    //   element_size_um[1] = stack->step().y();
    //   element_size_um[0] = stack->step().z();
    //   attribute.write( PredType::NATIVE_FLOAT, element_size_um);

    //   // TODO: copy image data into 1d array
    //   std::vector<short> buf(dims_1d);

    //   HVecUS &data = store->data();
    //   Point3i imgSize(stack->size());
    //   int xDim = imgSize.x();
    //   int yDim = imgSize.y();
    //   int zDim = imgSize.z();
    //   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);
    //         Point3i pIn(z, y, x);
    //         if (stack->boundsOK(x,y,z))
    //           buf[stack->offset(pOut)] = data[stack->offset(pOut)];
    //       }
    //     }
    //   }

      
    //   dataset.write(buf.data(), PredType::NATIVE_UINT16, dataspace, dataspace);

    //   delete plist;
    //   file.close();
    } else
      throw QString("Invalid file type %1").arg(filename);

    store->setFile(stripCurrentDir(filename));

    SETSTATUS("Saved Stack " << store->file() << ", Size: " << stack->size() << " voxel size:" << stack->step());
    return true;
  }
  
  bool StackSave::initialize(QWidget* parent)
  {
    int stackId = parm("Stack Number").toInt();
    bool is_work_store = stringToWorkStore(parm("Store"));
    if(!checkState().store((is_work_store ? STORE_WORK : STORE_MAIN), stackId))
      return false;
    QString mgxsFilter("MorphoGraphX Stack files (*.mgxs)"), inrFilter("Inria Stack files (*.inr)")/*, h5Filter("HDF5 (*.h5 *.hdf5 *.hdf)")*/;
    QString tiffFilter("TIFF Stack files (*.tif *.tiff)"), filter;
    Stack* s = stack(stackId);
    Store* store = (stringToWorkStore(parm("Store")) ? s->work() : s->main());
    QString filename;
    // If calling from GUI, use the active file
    if(parent)
      filename = store->file();
    if(filename.isEmpty())
      filename = parm("Filename");
    // If we loaded a file before, strip of the ending
    if(filename.right(4) == ".txt") {
      filename = filename.left(filename.length() - 4);
      filename += ".mgxs";
    }

    if(filename.endsWith(".inr", Qt::CaseInsensitive))
      filter = inrFilter;
    else if(filename.endsWith(".tif", Qt::CaseInsensitive) or filename.endsWith(".tiff", Qt::CaseInsensitive))
      filter = tiffFilter;
    // else if(filename.endsWith(".h5", Qt::CaseInsensitive) or filename.endsWith(".hdf", Qt::CaseInsensitive) or filename.endsWith(".hdf5", Qt::CaseInsensitive))
    //   filter = h5Filter;
    else
      filter = mgxsFilter;

    // Get the file name
    if(parent)
      filename = QFileDialog::getSaveFileName(parent, QString("Select stack file"), filename,
                                              QString(mgxsFilter + ";;" + inrFilter + ";;" + tiffFilter/* + ";;" + h5Filter*/), &filter,
                                              FileDialogOptions);
    if(filename.isEmpty())
      return false;

    // Check ending
    if(filter == inrFilter and !filename.endsWith(".inr", Qt::CaseInsensitive))
      filename += ".inr";
    else if(filter == tiffFilter
            and !(filename.endsWith(".tif", Qt::CaseInsensitive) or filename.endsWith(".tiff", Qt::CaseInsensitive)))
      filename += ".tif";
    else if(filter == mgxsFilter and !filename.endsWith(".mgxs", Qt::CaseInsensitive))
      filename += ".mgxs";
     // else if(filter == h5Filter and !(filename.endsWith(".h5", Qt::CaseInsensitive) or filename.endsWith(".hdf", Qt::CaseInsensitive) or filename.endsWith(".hdf5", Qt::CaseInsensitive)))
     //  filename += ".h5";

    setParm("Filename",filename);
    return true;
  }
  
  bool StackExport::initialize(QWidget* parent)
  {
    int stackId = parm("Stack Number").toInt();
    bool is_work_store = stringToWorkStore(parm("Store"));
    if(!checkState().store((is_work_store ? STORE_WORK : STORE_MAIN), stackId))
      return false;
    Stack* s = stack(stackId);
    Store* store = (stringToWorkStore(parm("Store")) ? s->work() : s->main());
    QString filename = store->file();

    if(!parent) return true;

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

    connect(ui.SelectImageFile, SIGNAL(clicked()), this, SLOT(selectImageFile()));
    connect(ui.ImageType, SIGNAL(currentIndexChanged(const QString &)), this, SLOT(setImageType(const QString &)));

    QStringList suppTypes = supportedImageWriteFormats();
    QStringList types;
    foreach(QString f, suppTypes) {
      types << f;
    }

    ui.ImageType->addItems(types);
    ui.generateVoxelSpacing->setChecked(stringToBool(parm("Generate Voxel Spacing")));

    if(dlg.exec() == QDialog::Accepted) {
      setParm("Format",suppTypes[ui.ImageType->currentIndex()]);
      setParm("Filename", ui.ImageFile->text());
      setParm("Generate Voxel Spacing", (ui.generateVoxelSpacing->isChecked()) ? "Yes" : "No");
      setParm("Num Digits", QString::number(ui.nbDigits->value()));
      return true;
    }
    return false;
  }
  
  void StackExport::selectImageFile()
  {
    QStringList filter;
    foreach(QString f, supportedImageWriteFormats()) {
      if(f != "CImg Auto")
        filter << QString("%1 images (*.%2 *.%1)").arg(f.toUpper()).arg(f.toLower());
      else
        filter << "CImg Auto (*.*)";
    }
    QString filename = ui->ImageFile->text();
    QString selected;
    filename = QFileDialog::getSaveFileName(dlg, "Filename to export to", filename, filter.join(";;"), &selected);
    if(!filename.isEmpty()) {
      int idx = filter.indexOf(selected);
      ui->ImageType->setCurrentIndex(idx);
      if(!selected.startsWith("CImg Auto")) {
        QString type = selected.left(selected.indexOf(" ")).toLower();
        if(!filename.endsWith("." + type, Qt::CaseInsensitive)) {
          filename += "." + type;
        }
      }
      ui->ImageFile->setText(filename);
    }
  }
  
  void StackExport::setImageType(const QString& type)
  {
    if(!type.startsWith("CImg Auto")) {
      QString t = type.toLower();
      QString fn = ui->ImageFile->text();
      int idx = fn.lastIndexOf(".");
      if(idx != -1)
        fn = fn.left(idx + 1) + t;
      else
        fn += "." + t;
      ui->ImageFile->setText(fn);
    }
  }
  
  bool StackExport::run()
  {
    int stackId = parm("Stack Number").toInt();
    bool is_work_store = stringToWorkStore(parm("Store"));
    if(!checkState().store((is_work_store ? STORE_WORK : STORE_MAIN), stackId))
      return false;
    Stack* s = stack(stackId);
    Store* store = (stringToWorkStore(parm("Store")) ? s->work() : s->main());
    QString filename = parm("Filename");
    QString type = parm("Format");
    uint nb_digits = parm("Num Digits").toUInt();
    if(filename.isEmpty()) {
      setErrorMessage("Error, filename provided is empty");
      return false;
    }
    if(type.isEmpty()) {
      setErrorMessage("Error, type provided is empty");
      return false;
    }
    if(nb_digits > 500) {
      setErrorMessage("Error, the number of digits asked for is greater than 500. It seems unreasonnably high.");
      return false;
    }
    return run(s, store, filename, type, nb_digits, stringToBool(parm("Generate Voxel Spacing")));
  }
  
  bool StackExport::run(Stack* stack, Store* store, const QString& filename, const QString& type,
                        unsigned int nb_digits, bool generate_voxel_spacing)
  {
    actingFile(filename);
    Point3u size = stack->size();
    Point3f step = stack->step();
    Image3D image(store->data(), size, step);
    saveImage(filename, image, type, nb_digits);
    if(generate_voxel_spacing) {
      QFileInfo fi(filename);
      QDir d = fi.dir();
      QFile file(d.filePath("voxelspacing.txt"));
      if(!file.open(QIODevice::WriteOnly)) {
        setErrorMessage(QString("Error, cannot open file '%1' for writing.").arg(file.fileName()));
        return false;
      }
      QTextStream ss(&file);
      ss << "x " << stack->step().x() << endl << "y " << stack->step().y() << endl << "z " << stack->step().z()
         << endl;
      file.close();
    }
    return true;
  }
  
  QString MeshExport::properFile(QString filename, const QString& type) const
  {
    QFileInfo fi(filename);
    QString suf = fi.suffix();
    if(!suf.isEmpty())
      filename = filename.left(filename.size() - suf.size() - 1);
    if(type == "Text" or type == "Cells")
      filename += ".txt";
    else if(type == "MeshEdit")
      filename += ".mesh";
    else if(type == "STL")
      filename += ".stl";
    else if(type.startsWith("VTK Mesh"))
      filename += ".vtu";
    else if(type == "OBJ")
      filename += ".obj";
    else if(type.startsWith("PLY"))
      filename += ".ply";
    return filename;
  }
  
  bool MeshExport::initialize(QWidget* parent)
  {

    int stackId = parm("Mesh Number").toInt();
    if(!checkState().mesh(MESH_ANY, stackId))
      return false;

    if(!parent) return true;

    QDialog dlg(parent);
    Ui_ExportMeshDialog ui;
    ui.setupUi(&dlg);

    this->ui = &ui;
    this->dlg = &dlg;

    connect(ui.SelectMeshFile, SIGNAL(clicked()), this, SLOT(selectMeshFile()));
    connect(ui.MeshType, SIGNAL(currentIndexChanged(const QString &)), this, SLOT(selectMeshType(const QString &)));

    Mesh* m = mesh(stackId);
    QString filename = m->file();
    if(filename.isEmpty())
      filename = parm("Filename");

    setMeshFile(filename);
    ui.Transform->setChecked(stringToBool(parm("Transform")));

    bool res = false;
    if(dlg.exec() == QDialog::Accepted) {
      setParm("Kind",ui.MeshType->currentText());
      setParm("Filename", ui.MeshFile->text());
      setParm("Transform", boolToString(ui.Transform->isChecked()));
      res = true;
    }
    this->ui = 0;
    this->dlg = 0;
    return res;
  }
  
  void MeshExport::selectMeshType(const QString& type)
  {
    QString filename = ui->MeshFile->text();
    if(filename.isEmpty())
      return;
    ui->MeshFile->setText(properFile(filename, type));
  }
  
  void MeshExport::selectMeshFile()
  {
    QString filename = ui->MeshFile->text();
    QStringList filters;
    filters << "Stanford Polygon file (*.ply)"
            << "VTK Mesh File (*.vtu)"
            << "Text or cell files (*.txt)"
            << "MeshEdit files (*.mesh)"
            << "STL CAD files (*.stl)"
            << "Wavefront Object files (*.obj)"
            << "All Mesh Files (*.ply *.vtu *.txt *.mesh *.stl *.obj)";
    QString filter;
    {
      QFileInfo fi(filename);
      QString suf = fi.suffix();
      if(suf == "ply")
        filter = filters[0];
      else if(suf == "vtu")
        filter = filters[1];
      else if(suf == "txt")
        filter = filters[2];
      else if(suf == "mesh")
        filter = filters[3];
      else if(suf == "stl")
        filter = filters[4];
      else if(suf == "obj")
        filter = filters[5];
      else
        filter = filters[6];
    }
    filename = QFileDialog::getSaveFileName(dlg, QString("Select mesh file"), filename, filters.join(";;"), &filter,
                                            FileDialogOptions);
    if(!filename.isEmpty())
      setMeshFile(filename);
  }
  
  void MeshExport::setMeshFile(const QString& filename)
  {
    QFileInfo fi(filename);
    QString suf = fi.suffix();
    // This is nasty, if you add/change anything here you'll have to change the exportmesh.ui file
    if(suf == "ply")
      ui->MeshType->setCurrentIndex(0);
    else if(suf == "vtu")
      ui->MeshType->setCurrentIndex(2);
    else if(suf == "txt")
      ui->MeshType->setCurrentIndex(4);
    else if(suf == "mesh")
      ui->MeshType->setCurrentIndex(6);
    else if(suf == "stl")
      ui->MeshType->setCurrentIndex(7);
    else if(suf == "obj")
      ui->MeshType->setCurrentIndex(8);
    else {
      ui->MeshType->setCurrentIndex(0);
      ui->MeshFile->setText(properFile(filename, "PLY Binary"));
      return;
    }
    ui->MeshFile->setText(filename);
  }
  
  bool MeshExport::run()
  {
    int meshId = parm("Mesh Number").toInt();
    if(!checkState().mesh(STACK_ANY, meshId))
      return false;
    Mesh* mesh = this->mesh(meshId);
    QString type = parm("Kind");
    if(parm("Kind") != "PLY Binary" and parm("Kind") != "PLY Ascii" and parm("Kind") != "VTK Mesh Binary" and
       parm("Kind") != "VTK Mesh Ascii" and parm("Kind") != "Text" and parm("Kind") != "Cells" and
       parm("Kind") != "MeshEdit" and parm("Kind") != "STL" and parm("Kind") != "OBJ") {
      setErrorMessage("Error, type must be one of 'PLY Binary', 'PLY Ascii',"
                      " 'VTK Mesh Binary', 'VTK Mesh Ascii', 'Text', 'Cells', 'MeshEdit', 'STL', 'OBJ'.");
      return false;
    }
    bool transform = stringToBool(parm("Transform"));
    bool res = run(mesh, parm("Filename"), type, transform);
    return res;
  }
  
  bool MeshExport::run(Mesh* mesh, const QString& filename, const QString& type, bool transform)
  {
    if(filename.isEmpty()) {
      setErrorMessage("Filename is empty, cannot save mesh");
      return false;
    }
    actingFile(filename);
    progressStart(QString("Saving Mesh %1").arg(mesh->userId()), 0, false);
    QFileInfo fi(filename);
    QString suf = fi.suffix();
    bool success = false;
    extendedPLY = stringToBool(parm("Extended PLY"));
    if(type == "PLY Binary")
      success = savePLY(mesh, filename, transform, true);
    else if(type == "PLY Ascii")
      success = savePLY(mesh, filename, transform, false);
    if(type == "VTK Mesh Binary")
      success = saveVTU(mesh, filename, transform, true);
    else if(type == "VTK Mesh Ascii")
      success = saveVTU(mesh, filename, transform, false);
    else if(type == "Text")
      success = saveText(mesh, filename, transform);
    else if(type == "Cells")
      success = saveCells(mesh, filename, transform);
    else if(type == "MeshEdit")
      success = saveMeshEdit(mesh, filename, transform);
    else if(type == "STL")
      success = saveMeshSTL(mesh, filename, transform);
    else if(type == "OBJ")
      success = saveOBJ(mesh, filename, transform);
    if(success) {
      mesh->setTransformed(transform);
      mesh->updateAll();
    }
    return success;
  }
  
  Point3d MeshExport::savedPos(Point3d pos, bool transform, const Stack* stack)
  {
    if(transform)
      pos = Point3d(stack->getFrame().inverseCoordinatesOf(Vec(pos)));
    return pos;
  }

  static QString vtkBinaryEncoding(const char* data, int size)
  {
    QByteArray ba = (QByteArray::fromRawData((const char*)&size, 4)
                     + QByteArray::fromRawData(data, size)).toBase64();
    return QString::fromLocal8Bit(ba);
  }
  
  template <typename T>
  static QDomText vtkBinaryEncoding(const std::vector<T>& data, QDomDocument& doc)
  {
    QDomText text = doc.createTextNode(vtkBinaryEncoding((const char*)&data[0], data.size()
                                       * sizeof(T)));
    return text;
  }
  
  QString vtkAsciiEncoding(const int& i)
  {
    return QString::number(i);
  }
  
  QString vtkAsciiEncoding(const float& f)
  {
    return QString::number(f, 'g', 10);
  }
  
  QString vtkAsciiEncoding(const Point3f& p)
  {
    return QString("%1 %2 %3").arg(p.x(), 0, 'g', 10).arg(p.y(), 0, 'g', 10).arg(p.z(), 0, 'g', 10);
  }
  
  QString vtkAsciiEncoding(const Point2i& p)
  {
    return QString("%1 %2").arg(p.x()).arg(p.y());
  }
  
  QString vtkAsciiEncoding(const Point3i& p)
  {
    return QString("%1 %2 %3").arg(p.x()).arg(p.y()).arg(p.z());
  }
  
  template <typename T>
  static QDomText vtkAsciiEncoding(const std::vector<T>& values, QDomDocument& doc)
  {
    QStringList strings;
#if QT_VERSION >= 0x040700
    strings.reserve(values.size());
#endif
    forall(const T& v, values)
        strings << vtkAsciiEncoding(v);
    return doc.createTextNode(strings.join("\n"));
  }
  
  template <typename T>
  static QDomText vtkEncoding(const std::vector<T>& values, QDomDocument& doc, bool binary)
  {
    if(binary)
      return vtkBinaryEncoding(values, doc);
    else
      return vtkAsciiEncoding(values, doc);
  }
  
  bool MeshExport::saveOBJ(Mesh* mesh, const QString& filename, bool transform)
  {
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

    QTextStream ts(&file);

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

    std::unordered_map<vertex, int> inv_vs;

    ts << "# Triangular mesh created by MorphoGraphX" << endl;
    ts << QString::fromWCharArray(L"# Length unit: \xb5m") << endl << endl;

    ts << "# Vertices" << endl;

    // write vertices
    int i = 0;
    forall(const vertex& v, S) {
      inv_vs[v] = i + 1;
      ts << "v " << savedPos(v->pos, transform, mesh->stack()) << endl;
      ts << "vn " << normalized(savedPos(v->nrml, transform, mesh->stack())) << endl;
      ++i;
    }

    ts << endl << "# Triangles" << endl;

    // write triangles
    forall(const vertex& v, S) {
      forall(const vertex& n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v, n);
        if(S.uniqueTri(v, n, m)) {
          ts << "f " << inv_vs[v] << " " << inv_vs[n] << " " << inv_vs[m] << endl;
        }
      }
    }

    file.close();
    return true;
  }
  
  bool MeshExport::saveVTU(Mesh* mesh, const QString& filename, bool transform, bool binary)
  {
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

    QDomDocument doc("VTKFile");
    QDomElement root = doc.createElement("VTKFile");
    root.setAttribute("type", "UnstructuredGrid");
    root.setAttribute("byte_order", "LittleEndian");
    root.setAttribute("version", "0.1");
    doc.appendChild(root);
    // Add the grid
    QDomElement grid = doc.createElement("UnstructuredGrid");
    root.appendChild(grid);
    // Add the pieces
    QDomElement piece = doc.createElement("Piece");
    grid.appendChild(piece);
    const vvGraph& S = mesh->graph();
    piece.setAttribute("NumberOfPoints", (uint)S.size());
    QString format;
    if(binary)
      format = "binary";
    else
      format = "ascii";
    // Points
    QDomElement points = doc.createElement("Points");
    piece.appendChild(points);
    QDomElement points_pos = doc.createElement("DataArray");
    points_pos.setAttribute("NumberOfComponents", "3");
    points_pos.setAttribute("type", "Float32");
    points_pos.setAttribute("format", format);
    std::vector<Point3f> points_pos_data(S.size());
    points.appendChild(points_pos);
    // Points attributes
    QDomElement points_data = doc.createElement("PointData");
    piece.appendChild(points_data);
    QString signal_name;
    if(mesh->signalUnit().isEmpty())
      signal_name = "Signal";
    else
      signal_name = QString("Signal (%1)").arg(mesh->signalUnit());
    points_data.setAttribute("Scalars", signal_name);
    points_data.setAttribute("Normals", "Normals");
    QDomElement points_signal = doc.createElement("DataArray");
    points_signal.setAttribute("Name", signal_name);
    points_signal.setAttribute("type", "Float32");
    points_signal.setAttribute("NumberOfComponents", "1");
    points_signal.setAttribute("format", format);
    std::vector<float> points_signal_data(S.size());
    QDomElement points_normal = doc.createElement("DataArray");
    points_normal.setAttribute("Name", "Normals");
    points_normal.setAttribute("type", "Float32");
    points_normal.setAttribute("NumberOfComponents", "3");
    points_normal.setAttribute("format", format);
    std::vector<Point3f> points_normal_data(S.size());
    QDomElement points_label = doc.createElement("DataArray");
    points_label.setAttribute("Name", "Label");
    points_label.setAttribute("type", "Int32");
    points_label.setAttribute("NumberOfComponents", "1");
    points_label.setAttribute("format", format);
    std::vector<int> points_label_data(S.size());
    points_data.appendChild(points_signal);
    points_data.appendChild(points_normal);
    points_data.appendChild(points_label);
    // Initialize vertex number and counts cells
    int saveId = 0;
    size_t count_cells = 0;
    forall(const vertex& v, S) {
      v->saveId = saveId++;
      forall(const vertex& n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v, n);
        if(S.uniqueTri(v, n, m))
          ++count_cells;
      }
    }
    // Cells
    piece.setAttribute("NumberOfCells", (uint)count_cells);
    QDomElement cells = doc.createElement("Cells");
    piece.appendChild(cells);
    QDomElement cells_connectivity = doc.createElement("DataArray");
    cells_connectivity.setAttribute("type", "Int32");
    cells_connectivity.setAttribute("Name", "connectivity");
    cells_connectivity.setAttribute("format", format);
    std::vector<Point3i> cells_connectivity_data(count_cells);
    QDomElement cells_offsets = doc.createElement("DataArray");
    cells_offsets.setAttribute("type", "Int32");
    cells_offsets.setAttribute("Name", "offsets");
    cells_offsets.setAttribute("format", format);
    std::vector<int> cells_offsets_data(count_cells);
    QDomElement cells_types = doc.createElement("DataArray");
    cells_types.setAttribute("type", "UInt8");
    cells_types.setAttribute("Name", "types");
    cells_types.setAttribute("format", format);
    std::vector<quint8> cells_types_data(count_cells, quint8(5));
    cells.appendChild(cells_connectivity);
    cells.appendChild(cells_offsets);
    cells.appendChild(cells_types);
    // Cells attributes
    QDomElement cells_data = doc.createElement("CellData");
    piece.appendChild(cells_data);
    QDomElement cells_label = doc.createElement("DataArray");
    cells_label.setAttribute("Name", "Label");
    cells_label.setAttribute("type", "Int32");
    cells_label.setAttribute("NumberOfComponents", "1");
    cells_label.setAttribute("format", format);
    cells_data.appendChild(cells_label);
    std::vector<int> cells_label_data(count_cells);
    QDomElement cells_heat = doc.createElement("DataArray");
    const IntFloatAttr& mesh_cell_heat = mesh->labelHeat();
    const IntIntFloatAttr& mesh_wall_heat = mesh->wallHeat();
    bool has_cell_heat = !mesh_cell_heat.empty();
    bool has_wall_heat = !mesh_wall_heat.empty();
    QString heat_name = "Heat";
    if(has_cell_heat)
      heat_name = "Cell " + heat_name;
    else if(has_wall_heat)
      heat_name = "Wall " + heat_name;
    if(!mesh->heatMapUnit().isEmpty())
      heat_name += QString(" (%1)").arg(mesh->heatMapUnit());

    float minHeat = mesh->heatMapBounds()[0];
    float maxHeat = mesh->heatMapBounds()[1];

    cells_heat.setAttribute("Name", heat_name);
    cells_heat.setAttribute("type", "Float32");
    cells_heat.setAttribute("NumberOfComponents", "1");
    cells_heat.setAttribute("format", format);
    cells_heat.setAttribute("RangeMin", minHeat);
    cells_heat.setAttribute("RangeMax", maxHeat);
    std::vector<float> cells_heat_data;
    if(has_cell_heat or has_wall_heat) {
      cells_heat_data.resize(count_cells);
      cells_data.appendChild(cells_heat);
      cells_data.setAttribute("Scalars", heat_name);
    } else
      cells_data.setAttribute("Scalars", "Label");

    std::unordered_set<int> labels;

    // Fills in everything
    int cell_id = 0;
    forall(const vertex& v, S) {
      int id = v->saveId;
      Point3f pos(savedPos(v->pos, transform, mesh->stack()));
      Point3f nrml(normalized(savedPos(v->nrml, transform, mesh->stack())));
      points_pos_data[id] = pos;
      points_label_data[id] = v->label;
      points_signal_data[id] = v->signal;
      points_normal_data[id] = nrml;
      forall(const vertex& n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v, n);
        if(S.uniqueTri(v, n, m)) {
          cells_connectivity_data[cell_id] = Point3i(id, n->saveId, m->saveId);
          cells_offsets_data[cell_id] = 3 * (cell_id + 1);
          int label = getLabel(v, n, m);
          labels.insert(label);
          cells_label_data[cell_id] = label;
          if(has_cell_heat) {
            IntFloatAttr::const_iterator found = mesh_cell_heat.find(label);
            if(found != mesh_cell_heat.end())
              cells_heat_data[cell_id] = found->second;
            else
              cells_heat_data[cell_id] = std::numeric_limits<float>::quiet_NaN();
          } else if(has_wall_heat) {
            IntIntPair wall;
            if(mesh->isBordTriangle(label, v, n, m, wall)) {
              IntIntFloatAttr::const_iterator found = mesh_wall_heat.find(wall);
              if(found != mesh_wall_heat.end())
                cells_heat_data[cell_id] = found->second;
              else
                cells_heat_data[cell_id] = std::numeric_limits<float>::quiet_NaN();
            } else
              cells_heat_data[cell_id] = std::numeric_limits<float>::quiet_NaN();
          }
          cell_id++;
        }
      }
    }

    // Global fields
    QDomElement global_fields = doc.createElement("FieldData");
    grid.appendChild(global_fields);

    QDomElement scale = doc.createElement("DataArray");
    scale.setAttribute("Name", "Scale");
    scale.setAttribute("type", "Float32");
    scale.setAttribute("NumberOfTuples", 1);
    scale.setAttribute("format", "ascii");
    scale.appendChild(doc.createTextNode("1e-6"));
    global_fields.appendChild(scale);

    QDomElement labels_field = doc.createElement("DataArray");
    labels_field.setAttribute("Name", "Labels");
    labels_field.setAttribute("type", "Int32");
    labels_field.setAttribute("NumberOfTuples", (int)labels.size());
    labels_field.setAttribute("format", format);
    std::vector<int> labels_fields_data(labels.begin(), labels.end());
    std::sort(labels_fields_data.begin(), labels_fields_data.end());
    labels_field.appendChild(vtkEncoding(labels_fields_data, doc, binary));
    global_fields.appendChild(labels_field);

    if(has_cell_heat) {
      QDomElement cell_heat = doc.createElement("DataArray");
      QString cell_heat_name = "Cell Heat";
      if(!mesh->heatMapUnit().isEmpty())
        cell_heat_name += QString(" (%1)").arg(mesh->heatMapUnit());
      cell_heat.setAttribute("Name", cell_heat_name);
      cell_heat.setAttribute("type", "Float32");
      cell_heat.setAttribute("NumberOfTuples", (int)labels.size());
      cell_heat.setAttribute("format", format);
      cell_heat.setAttribute("RangeMin", minHeat);
      cell_heat.setAttribute("RangeMax", maxHeat);
      std::vector<float> cell_heat_data(labels.size());
      for(size_t i = 0; i < labels.size(); ++i) {
        int label = labels_fields_data[i];
        IntFloatAttr::const_iterator found = mesh_cell_heat.find(label);
        if(found == mesh_cell_heat.end())
          cell_heat_data[i] = std::numeric_limits<float>::quiet_NaN();
        else
          cell_heat_data[i] = found->second;
      }
      cell_heat.appendChild(vtkEncoding(cell_heat_data, doc, binary));
      global_fields.appendChild(cell_heat);
    } else if(has_wall_heat) {
      QDomElement walls_fields = doc.createElement("DataArray");
      walls_fields.setAttribute("Name", "Walls");
      walls_fields.setAttribute("type", "Int32");
      walls_fields.setAttribute("NumberOfTuples", (int)mesh_wall_heat.size());
      walls_fields.setAttribute("format", format);
      std::vector<Point2i> walls_fields_data(mesh_wall_heat.size());

      QDomElement wall_heat = doc.createElement("DataArray");
      QString wall_heat_name = "Wall Heat";
      if(!mesh->heatMapUnit().isEmpty())
        wall_heat_name += QString(" (%1)").arg(mesh->heatMapUnit());
      wall_heat.setAttribute("Name", wall_heat_name);
      wall_heat.setAttribute("type", "Float32");
      wall_heat.setAttribute("NumberOfTuples", (int)mesh_wall_heat.size());
      wall_heat.setAttribute("format", format);
      wall_heat.setAttribute("RangeMin", minHeat);
      wall_heat.setAttribute("RangeMax", maxHeat);
      std::vector<float> wall_heat_data(labels.size());
      int i = 0;
      forall(const IntIntFloatPair& vs, mesh_wall_heat) {
        const IntIntPair& w = vs.first;
        float value = vs.second;
        walls_fields_data[i] = Point2i(w.first, w.second);
        wall_heat_data[i] = value;
        ++i;
      }
      walls_fields.appendChild(vtkEncoding(walls_fields_data, doc, binary));
      wall_heat.appendChild(vtkEncoding(wall_heat_data, doc, binary));
      global_fields.appendChild(walls_fields);
      global_fields.appendChild(wall_heat);
    }

    points_pos.appendChild(vtkEncoding(points_pos_data, doc, binary));
    points_label.appendChild(vtkEncoding(points_label_data, doc, binary));
    points_signal.appendChild(vtkEncoding(points_signal_data, doc, binary));
    points_normal.appendChild(vtkEncoding(points_normal_data, doc, binary));

    cells_connectivity.appendChild(vtkEncoding(cells_connectivity_data, doc, binary));
    cells_offsets.appendChild(vtkEncoding(cells_offsets_data, doc, binary));
    cells_types.appendChild(vtkEncoding(cells_types_data, doc, binary));

    if(has_cell_heat or has_wall_heat)
      cells_heat.appendChild(vtkEncoding(cells_heat_data, doc, binary));

    cells_label.appendChild(vtkEncoding(cells_label_data, doc, binary));

    QString xmlString = doc.toString();
    QByteArray xmlArray = xmlString.toUtf8();
    qint64 qty = file.write(xmlArray);
    if(qty != xmlArray.size())
      throw QString("Could not write data to file (written = %1, expected = %2)").arg(qty).arg(xmlArray.size());

    file.close();

    return true;
  }
  
  bool MeshExport::saveText(Mesh* mesh, const QString& filename, bool transform)
  {
    const Stack* stack = mesh->stack();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

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

    QTextStream out(&file);
    int saveId = 0;
    int i = 0;

    out << S.size() << endl;

    progressStart(QString("Saving Text Mesh %1").arg(mesh->userId()), S.size() * 2);
    forall(const vertex& v, S) {
      v->saveId = saveId++;
      Point3f pos(savedPos(v->pos, transform, stack));
      out << v->saveId << " " << pos << " " << v->label << endl;
      if(!progressAdvance(i))
        userCancel();
      ++i;
    }
    forall(const vertex& v, S) {
      out << v->saveId << " " << S.valence(v);
      forall(const vertex& n, S.neighbors(v))
          out << " " << n->saveId;
      out << endl;
      if(!progressAdvance(i))
        userCancel();
      ++i;
    }
    if(!progressAdvance(S.size() * 2))
      userCancel();
    file.close();
    SETSTATUS("Saved mesh, file:" << mesh->file() << ", vertices:" << S.size());
    return true;
  }
  
  bool MeshExport::saveCells(Mesh* mesh, const QString& filename, bool transform)
  {
    const Stack* stack = mesh->stack();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

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

    QTextStream out(&file);
    int saveId = 0;
    progressStart(QString("Saving Cells %1").arg(mesh->userId()), S.size() * 2);
    out << S.size() << endl;
    int i = 0;
    forall(const vertex& v, S) {
      v->saveId = saveId++;
      Point3f pos(savedPos(v->pos, transform, stack));
      out << v->saveId << " " << pos << " " << v->label << " " << v->type << endl;
      if(!progressAdvance(i))
        userCancel();
      ++i;
    }
    forall(const vertex& v, S) {
      out << v->saveId << " " << S.valence(v);
      forall(const vertex& n, S.neighbors(v))
          out << " " << n->saveId;
      out << endl;
      if(!progressAdvance(i))
        userCancel();
      ++i;
    }
    if(!progressAdvance(S.size() * 2))
      userCancel();
    SETSTATUS("Saved cellular mesh, file:" << mesh->file() << ", vertices:" << S.size());
    return true;
  }
  
  bool MeshExport::saveMeshEdit(Mesh* mesh, const QString& filename, bool transform)
  {
    const Stack* stack = mesh->stack();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

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

    QTextStream out(&file);
    int saveId = 1;
    out << "MeshVersionFormatted 1" << endl;
    out << "Dimension 3" << endl;
    out << "Vertices" << endl;
    out << S.size() << endl;
    progressStart(QString("Saving Mesh Edit %1").arg(mesh->userId()), S.size());
    int i = 0;
    forall(const vertex& v, S) {
      v->saveId = saveId++;
      Point3f pos(savedPos(v->pos, transform, stack));
      out << pos << " " << v->label << endl;
    }
    // Count the triangles
    int count = 0;
    forall(const vertex& v, S)
        forall(const vertex& n, S.neighbors(v)) {
      vertex m = S.nextTo(v, n);
      if(!S.uniqueTri(v, n, m))
        continue;
      count++;
    }
    out << "Triangles" << endl;
    out << count << endl;
    forall(const vertex& v, S) {
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
        out << v->saveId << " " << n->saveId << " " << m->saveId << " ";
        out << getLabel(v, n, m) << endl;
      }
      if(!progressAdvance(i))
        userCancel();
    }
    if(!progressAdvance(S.size()))
      userCancel();
    out << "End" << endl;
    SETSTATUS("Saving triangle mesh, file:" << mesh->file() << ", vertices:" << S.size());
    return true;
  }
  
  bool MeshExport::saveMeshSTL(Mesh* mesh, const QString& filename, bool transform)
  {
    const Stack* stack = mesh->stack();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

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

    QTextStream out(&file);
    progressStart(QString("Saving STL Mesh %1").arg(mesh->userId()), S.size());
    int i = 0;

    out << "solid mgx" << endl;
    const std::vector<vertex>& av = mesh->activeVertices();
    forall(const vertex& v, av) {
      forall(const vertex& n, S.neighbors(v)) {
        vertex m = S.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
        Point3f vpos(savedPos(v->pos, transform, stack));
        Point3f npos(savedPos(n->pos, transform, stack));
        Point3f mpos(savedPos(m->pos, transform, stack));
        Point3f nrml(((npos - vpos) ^ (mpos - vpos)).normalize());
        out << "facet normal " << nrml.x() << " " << nrml.y() << " " << nrml.z() << endl;
        out << "  outer loop" << endl;
        out << "    vertex " << vpos.x() << " " << vpos.y() << " " << vpos.z() << endl;
        out << "    vertex " << npos.x() << " " << npos.y() << " " << npos.z() << endl;
        out << "    vertex " << mpos.x() << " " << mpos.y() << " " << mpos.z() << endl;
        out << "  endloop" << endl;
        out << "endfacet" << endl;
      }
      if(!progressAdvance(i))
        userCancel();
    }
    out << "endsolid mgx" << endl;
    if(!progressAdvance(S.size()))
      userCancel();
    SETSTATUS("Saving STL mesh, file:" << mesh->file() << ", vertices:" << S.size());
    return true;
  }
  
  bool MeshExport::savePLY(Mesh* mesh, const QString& filename, bool transform, bool binary)
  {
    const Stack* stack = mesh->stack();
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshExport::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

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

    if(mesh->meshType() == "MGXC") {
      progressStart(QString("Saving Cell Mesh").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;
        }
      }

      // 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 Cellular mesh" << endl;
      outText << "element vertex " << junctions << endl;
      outText << "property float x" << endl;
      outText << "property float y" << endl;
      outText << "property float z" << endl;
      outText << "property float signal" << endl;    
      outText << "element face " << cells << endl;
      outText << "property list int int vertex_index" << endl;
      outText << "property int label" << endl;
      outText << "end_header" << endl;

      // First write vertices
      forall(const vertex& v, S) {
        if(v->type == 'c')
          continue;

        Point3f pos(savedPos(v->pos, transform, stack));
        if(binary)
          outBinary << pos.x() << pos.y() << pos.z() << v->signal;
        else
          outText << pos << " " << v->signal << endl;
        if(!progressAdvance(0))
          userCancel();
      }

      // Now write cells (faces)
      forall(const vertex& v, S) {
        if(v->type == 'j')
          continue;

        if(binary)
          outBinary << qint32(S.valence(v));
        else
          outText << S.valence(v);

        // Write cell polygon
        forall(const vertex& n, S.neighbors(v))
            if(binary)
            outBinary << qint32(n->saveId);
        else
        outText << " " << n->saveId;
        
        if(binary)
          outBinary << qint32(v->label);
        else
          outText << " " << v->label << endl;

        if(!progressAdvance(0))
          userCancel();
      }
      if(!progressAdvance(0))
        userCancel();

      SETSTATUS("Saved cellular mesh to file:"
                << mesh->file() << ", cells:" << cells << ", vertices:" << junctions);

    } else {
      // Save triangle mesh
      progressStart(
            QString("Saving Triangle Mesh %1").arg(mesh->userId()), 0);

      // First count the junctions
      int vertices = 0;
      int triangles = 0;
      forall(const vertex& v, S) {
        v->saveId = vertices++;
        forall(const vertex& n, S.neighbors(v)) {
          const vertex& m = S.nextTo(v, n);
          if(!S.uniqueTri(v, n, m))
            continue;
          triangles++;
        }
      }

      // 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 Triangle mesh" << endl;
      outText << "element vertex " << vertices << endl;
      outText << "property float x" << endl;
      outText << "property float y" << endl;
      outText << "property float z" << endl;
      outText << "property int label" << endl;
      outText << "property float signal" << endl;
      if(extendedPLY){
        outText << "property int parent" << endl; // parent label
        outText << "property int cellNr" << endl; // marion's "vertex belong to so many cells"
        outText << "property float xNrml" << endl;
        outText << "property float yNrml" << endl;
        outText << "property float zNrml" << endl;
        outText << "property uchar red" << endl;
        outText << "property uchar green" << endl;
        outText << "property uchar blue" << endl;
      }
      outText << "element face " << triangles << endl;
      outText << "property list uchar int vertex_index" << endl;
      outText << "property int label" << endl;
      outText << "end_header" << endl;

      
      // TODO this is just a temporary hack for marion to save an attribute map
      typedef AttrMap<vertex, int> VtxCounterAttr;
      VtxCounterAttr *vtxC;
      vtxC = &mesh->attributes().attrMap<vertex, int>("Vertex Cell Counter");


      // First write vertices
      forall(const vertex& v, S) {
        Point3f pos(savedPos(v->pos, transform, stack));
        int r = ImgData::LabelColors[v->label % 16].r() * 255;
        int b = ImgData::LabelColors[v->label % 16].b() * 255;
        int g = ImgData::LabelColors[v->label % 16].g() * 255;
        if(binary){
          outBinary << pos.x() << pos.y() << pos.z() << qint32(v->label) << v->signal;
        
          if(extendedPLY){
            outBinary << qint32(mesh->parents()[v->label]) 
                      << qint32(v->*vtxC) 
                      << v->nrml.x() << v->nrml.y() << v->nrml.z()
                      << quint8(r) << quint8(g) << quint8(b);
          }
        } else {
          outText << pos << " " << v->label << " " << v->signal;
          if(extendedPLY){
            outText << " " << mesh->parents()[v->label] << " " << v->*vtxC << " " 
                    << v->nrml << " "
                    << r << " " << g << " " << b;
          }
          outText << endl;
        }
        if(!progressAdvance(0))
          userCancel();
      }
      // Now write cells or triangles (faces)
      forall(const vertex& v, S) {
        // Write fan of triangles
        forall(const vertex& n, S.neighbors(v)) {
          const vertex& m = S.nextTo(v, n);
          if(!S.uniqueTri(v, n, m))
            continue;
          if(binary)
            outBinary << quint8(3) << qint32(v->saveId) << qint32(n->saveId) << qint32(m->saveId)
                      << qint32(getLabel(v, n, m));
          else
            outText << "3 " << v->saveId << " " << n->saveId << " " << m->saveId << " "
                    << getLabel(v, n, m) << endl;
        }
        if(!progressAdvance(0))
          userCancel();
      }
      if(!progressAdvance(0))
        userCancel();

      SETSTATUS("Saved triangle mesh to file:" << mesh->file() << ", triangles:" << triangles
                << ", vertices:" << vertices);
    }
    return true;
  }
  
  bool MeshSave::initialize(QWidget* parent)
  {
    int meshId = parm("Mesh Number").toInt();
    if(!checkState().mesh(MESH_ANY, meshId))
      return false;

    Mesh* m = mesh(meshId);
    QString filename = m->file();
    if(filename.isEmpty())
      filename = parm("Filename");

    if(!parent) return true;

    QDialog dlg(parent);
    Ui_SaveMeshDialog ui;
    ui.setupUi(&dlg);

    this->ui = &ui;
    this->dlg = &dlg;

    connect(ui.SelectMeshFile, SIGNAL(clicked()), this, SLOT(selectMeshFile()));

    setMeshFile(filename);
    ui.Transform->setChecked(stringToBool(parm("Transform")));

    bool res = false;
    if(dlg.exec() == QDialog::Accepted) {
      setParm("Filename", properFile(ui.MeshFile->text()));
      setParm("Transform", (ui.Transform->isChecked()) ? "yes" : "no");
      res = true;
    }
    this->ui = 0;
    this->dlg = 0;
    return res;
  }
  
  QString MeshSave::properFile(QString filename) const
  {
    QFileInfo fi(filename);
    QString suf = fi.suffix();
    if(!suf.isEmpty())
      filename = filename.left(filename.size() - suf.size() - 1);
    return filename + ".mgxm";
  }
  
  void MeshSave::selectMeshFile()
  {
    QString filename = ui->MeshFile->text();
    filename = QFileDialog::getSaveFileName(dlg, QString("Save mesh as ..."), filename, "Mesh files (*.mgxm)", 0,
                                            FileDialogOptions);
    if(!filename.isEmpty())
      setMeshFile(filename);
  }
  
  void MeshSave::setMeshFile(const QString& filename)
  {
    if(!filename.endsWith(".mgxm", Qt::CaseInsensitive))
      ui->MeshFile->setText(filename + ".mgxm");
    else
      ui->MeshFile->setText(filename);
  }
  
  bool MeshSave::run()
  {
    int meshId = parm("Mesh Number").toInt();
    if(!checkState().mesh(STACK_ANY, meshId))
      return false;
    Mesh* mesh = this->mesh(meshId);
    bool transform = stringToBool(parm("Transform"));
    bool res = run(mesh, parm("Filename"), transform);
    return res;
  }
  
  Point3d MeshSave::savedPos(Point3d pos, bool transform, const Stack* stack)
  {
    if(transform)
      pos = Point3d(stack->getFrame().inverseCoordinatesOf(Vec(pos)));
    return pos;
  }
  
  bool MeshSave::run(Mesh* mesh, const QString& filename, bool transform)
  {
    actingFile(filename);
    QFile file(filename);
    if(!file.open(QIODevice::WriteOnly)) {
      setErrorMessage(QString("MeshSave::Error:Cannot open output file: %1").arg(filename));
      return false;
    }

    // Write out the mesh
    mesh->write(file, transform);

    file.close();
    mesh->setFile(stripCurrentDir(filename));

    SETSTATUS("Saved mesh, file:" << mesh->file());

    return true;
  }
  
  bool SaveViewFile::run() {
    return run(parm("Filename"));
  }
  
  bool SaveViewFile::run(const QString& filename)
  {qDebug() << filename;
    if(filename.isEmpty()) {
      setErrorMessage("Filename is empty. Cannot save view file.");
      return false;
    }
    actingFile(filename, true);
    bool res = saveView(filename);
    SETSTATUS("saving view to file " << filename);
    if(!res) {
      setErrorMessage("Error while saving view file.");
      return false;
    }
    return true;
  }
  
  bool SaveViewFile::initialize(QWidget* parent)
  {
    QString filename(parm("Filename"));

    // Get the file name
    if(parent)
      filename = QFileDialog::getSaveFileName(parent, QString("Select view file"), filename,
                                              "MorphoGraphX View files (*.mgxv);;All files (*.*)", 0, FileDialogOptions);

    if(filename.isEmpty())
      return false;
    if(!filename.endsWith(".mgxv", Qt::CaseInsensitive))
      filename += ".mgxv";

    setParm("Filename",filename);
    return true;
  }
 
  QString QuickSaveFile::addIncSuffix(QString filename, QString fileExtension)
  {
    QString newFile;

    QString suffix = "_MGX_";

    // remove the extension
    filename.remove(fileExtension, Qt::CaseInsensitive);

    // check if the ending already has _MGX_xxxx, if yes remove it

    QString leftString = filename.left(filename.length() - 9);
    QString suffixString = filename.right(filename.length() - leftString.length());

    QString suffixLeft = suffixString.left(suffixString.length() - 4);
    QString suffixRight = suffixString.right(suffixString.length() - 5);

    int counter = 0;

    counter += suffixRight.toInt();
    
    counter = extCounter;
    QString counterString;
    QString d0 = QString::number(counter % 10);
    counter/=10;
    QString d1 = QString::number(counter % 10);
    counter/=10;
    QString d2 = QString::number(counter % 10);
    counter/=10;
    QString d3 = QString::number(counter % 10);
    counter/=10;
    if(suffixLeft == suffix){
      newFile = leftString + suffix + d3 + d2 + d1 + d0 + fileExtension;
    } else {
      newFile = filename + suffix + d3 + d2 + d1 + d0 + fileExtension;
    }

    //qDebug() << "newFile " << newFile;
    return newFile;
  }


  void QuickSaveFile::saveMesh(int meshNr)
  {

    QString filename = mesh(meshNr)->file();
    
    if(filename != ""){
      MeshSave ms(*this);
      if(saveWithExt)
        ms.run(mesh(meshNr), addIncSuffix(filename, ".mgxm"), false);
      else
        ms.run(mesh(meshNr), filename, false);
    }

  }

  void QuickSaveFile::saveStore(int stackNr, Store* store, int compressionLevel)
  {

    QString filename = store->file();
    
    if(filename != ""){
      StackSave ss(*this);
      if(saveWithExt)
        ss.run(stack(stackNr), store, addIncSuffix(filename, ".mgxs"), compressionLevel);
      else
        ss.run(stack(stackNr), store, filename, compressionLevel);
      
    }

  }

  bool QuickSaveFile::initialize(QWidget* parent)
  {

    setParm("Extension Counter", QString::number(extCounter+1));
    return true;
  }

  bool QuickSaveFile::run()
  {
    if(stringToBool(parm("Save with File Extensions")))
      saveWithExt = true;
    else
      saveWithExt = false;

    extCounter = parm("Extension Counter").toInt();

    std::cout << "Quick Save!" << std::endl;

    // save mgxv file
    SaveViewFile saveViewProc(*this);

    if(saveWithExt) {
      QString newFile = addIncSuffix(saveViewProc.parm("Filename"), ".mgxv");
      saveView(newFile);
    } else
      saveView(saveViewProc.parm("Filename"));

    saveMesh(0);
    saveMesh(1);

    StackSave stackSave(*this);

    int compressionLevel = stackSave.parm("Compression Level").toInt();

    saveStore(0, stack(0)->main(), compressionLevel);
    saveStore(0, stack(0)->work(), compressionLevel);
    saveStore(1, stack(1)->main(), compressionLevel);
    saveStore(1, stack(1)->work(), compressionLevel);

    std::cout << "Quick Save finished" << std::endl;

    return true;
  }
  REGISTER_PROCESS(QuickSaveFile);

  class MeshListModel : public QAbstractTableModel {
  public:
    MeshListModel(Process* p, QObject* parent = 0)
      : QAbstractTableModel(parent)
      , proc(p)
    {
      for(int i = 0; i < proc->meshCount(); ++i) {
        names << QString("Mesh%1").arg(i + 1);
        Mesh* m = proc->mesh(i);
        QString shortFile = stripCurrentDir(m->file());
        selected << not (shortFile.isEmpty());
        files << shortFile;
        transforms << false;
      }
    }

    int rowCount(const QModelIndex& parent = QModelIndex()) const
    {
      if(parent.isValid())
        return 0;
      return names.size();
    }

    int columnCount(const QModelIndex& /*parent*/ = QModelIndex()) const {
      return 4;
    }

    QVariant headerData(int section, Qt::Orientation orientation, int role = Qt::DisplayRole) const
    {
      if(role != Qt::DisplayRole)
        return QVariant();
      if(orientation == Qt::Vertical)
        return QVariant();
      switch(section) {
      case 0:
        return "Mesh";
      case 1:
        return "Filename";
      case 2:
        return "Transformed";
      default:
        return QVariant();
      }
    }

    Qt::ItemFlags flags(const QModelIndex& index) const
    {
      if(!index.isValid())
        return Qt::ItemIsEnabled;
      int i = index.row();
      if(i < 0 or i >= names.size())
        return Qt::ItemIsEnabled;
      switch(index.column()) {
      case 0:
        return Qt::ItemIsEnabled | Qt::ItemIsUserCheckable | Qt::ItemIsSelectable;
      case 1:
        if(selected[i])
          return Qt::ItemIsEnabled | Qt::ItemIsEditable | Qt::ItemIsSelectable;
        else
          return 0;
      case 2:
        if(selected[i])
          return Qt::ItemIsEnabled | Qt::ItemIsUserCheckable | Qt::ItemIsSelectable;
        else
          return 0;
      default:
        return Qt::ItemIsEnabled;
      }
    }

    QString properFile(const QString& filename)
    {
      if(!filename.endsWith(".mgxm", Qt::CaseInsensitive))
        return filename + ".mgxm";
      return filename;
    }

    bool setData(const QModelIndex& index, const QVariant& value, int role = Qt::EditRole)
    {
      if(!index.isValid() or index.column() > 3)
        return false;

      int i = index.row();

      if(i < 0 or i >= names.size())
        return false;

      switch(index.column()) {
      case 0:
        if(role != Qt::CheckStateRole)
          return false;
        selected[i] = value.toBool();
        emit dataChanged(this->index(i, 0), this->index(i, 4));
        break;
      case 1:
        if(role != Qt::EditRole)
          return false;
        files[i] = properFile(value.toString());
        emit dataChanged(index, index);
        break;
      case 2:
        if(role != Qt::CheckStateRole)
          return false;
        transforms[i] = value.toBool();
        emit dataChanged(index, index);
        break;
      default:
        return false;
      }
      return true;
    }

    QVariant data(const QModelIndex& index, int role) const
    {
      if(!index.isValid() or index.row() >= names.size() or index.column() > 3)
        return QVariant();
      int i = index.row();
      switch(role) {
      case Qt::CheckStateRole:
        switch(index.column()) {
        case 0:
          return selected[i] ? Qt::Checked : Qt::Unchecked;
        case 2:
          return transforms[i] ? Qt::Checked : Qt::Unchecked;
        case 1:
        default:
          return QVariant();
        }
      case Qt::EditRole:
        if(index.column() == 1)
          return files[i];
        else
          return QVariant();
      case Qt::DisplayRole:
        if(index.column() == 0)
          return names[i];
        else if(index.column() == 1)
          return files[i];
        else
          return QVariant();
      default:
        return QVariant();
      }
    }

    void changeMesh(const QString& name, const QString& file, bool trans)
    {
      int idx = names.indexOf(name);
      if(idx >= 0) {
        files[idx] = file;
        selected[idx] = true;
        transforms[idx] = trans;
        emit dataChanged(index(idx, 0), index(idx, 4));
      }
    }

    QStringList typeList;
    QStringList names;
    QStringList files;
    QList<bool> selected, transforms;

  protected:
    Process* proc;
  };
  
  class StackListModel : public QAbstractTableModel {
  public:
    StackListModel(Process* p, QObject* parent = 0)
      : QAbstractTableModel(parent)
      , proc(p)
    {
      for(int i = 0; i < proc->stackCount(); ++i) {
        Stack* s = proc->stack(i);
        if(!s)
          continue;
        Store* main = s->main();
        if(main) {
          names << QString("MainStack%1").arg(i + 1);
          QString shortFile = stripCurrentDir(main->file());
          files << properFile(shortFile);
          selected << not shortFile.isEmpty();
        }
        Store* work = s->work();
        if(work) {
          names << QString("WorkStack%1").arg(i + 1);
          QString shortFile = stripCurrentDir(work->file());
          files << properFile(shortFile);
          selected << not shortFile.isEmpty();
        }
      }
    }

    QString properFile(const QString& filename)
    {
      QString result = filename;
      QFileInfo fi(filename);
      QString suf = fi.suffix();
      if(suf == "txt")
        result = filename.left(result.size() - 3) + "mgxs";
      else if(suf != "inr" and suf != "mgxs")
        result = filename + ".mgxs";
      return result;
    }

    int rowCount(const QModelIndex& parent = QModelIndex()) const
    {
      if(parent.isValid())
        return 0;
      return names.size();
    }

    int columnCount(const QModelIndex& /*parent*/ = QModelIndex()) const {
      return 2;
    }

    QVariant data(const QModelIndex& index, int role) const
    {
      if(!index.isValid() or index.row() >= names.size() or index.column() > 2)
        return QVariant();
      int i = index.row();
      if(role == Qt::CheckStateRole and index.column() == 0)
        return selected[i] ? Qt::Checked : Qt::Unchecked;
      if(role == Qt::EditRole and index.column() == 1)
        return files[i];
      if(role == Qt::DisplayRole) {
        if(index.column() == 0)
          return names[i];
        if(index.column() == 1)
          return files[i];
      }
      return QVariant();
    }

    QVariant headerData(int section, Qt::Orientation orientation, int role = Qt::DisplayRole) const
    {
      if(role != Qt::DisplayRole)
        return QVariant();
      if(orientation == Qt::Vertical)
        return QVariant();
      switch(section) {
      case 0:
        return "Stack";
      case 1:
        return "Filename";
      default:
        return QVariant();
      }
    }

    Qt::ItemFlags flags(const QModelIndex& index) const
    {
      if(!index.isValid())
        return Qt::ItemIsEnabled;
      int i = index.row();
      if(i < 0 or i >= names.size())
        return Qt::ItemIsEnabled;
      switch(index.column()) {
      case 0:
        return Qt::ItemIsEnabled | Qt::ItemIsUserCheckable | Qt::ItemIsSelectable;
      case 1:
        if(selected[i])
          return Qt::ItemIsEnabled | Qt::ItemIsEditable | Qt::ItemIsSelectable;
        else
          return 0;
      default:
        return Qt::ItemIsEnabled;
      }
    }

    bool setData(const QModelIndex& index, const QVariant& value, int role = Qt::EditRole)
    {
      if(!index.isValid() or index.column() > 1)
        return false;

      int i = index.row();

      if(i < 0 or i >= names.size())
        return false;

      if(index.column() == 0 and role == Qt::CheckStateRole) {
        selected[i] = value.toBool();
        emit dataChanged(this->index(i, 0), this->index(i, 1));
      } else if(index.column() == 1 and role == Qt::EditRole) {
        files[i] = properFile(value.toString());
        emit dataChanged(index, index);
      } else
        return false;

      return true;
    }

    void changeStack(const QString& name, const QString& file)
    {
      int idx = names.indexOf(name);
      if(idx >= 0) {
        files[idx] = file;
        selected[idx] = true;
        emit dataChanged(index(idx, 0), index(idx, 1));
      }
    }

    QStringList names;
    QStringList files;
    QList<bool> selected;

  protected:
    Process* proc;
  };
  
  bool TakeSnapshot::run(QString filename, float overSampling, int width, int height, int quality,
                         bool expand_frustum)
  {
    return takeSnapshot(filename, overSampling, width, height, quality, expand_frustum);
  }



  bool StackOpenSaveFolder::initialize(QWidget* parent)
  {
    QDir directory = QFileDialog::getExistingDirectory(parent, tr("select directory"));
    images = directory.entryList(QStringList() << "*.tiff" << "*.TIFF"<< "*.tif" << "*.TIF",QDir::Files);

    setParm("Folder", directory.absolutePath());
    
    return true;
  }
  
  bool StackOpenSaveFolder::run()
  {

    QString path = parm("Folder");
    QString prefix = parm("File Prefix");
    int compression = parm("Compression Level").toInt();

    if(compression < -1) compression = -1;
    if(compression > 9) compression = 9;

    Stack* stack = currentStack();
    Store* store = currentStack()->currentStore();

    foreach(QString filename, images) {
      //do whatever you need to do
      qDebug() << path << "/" << filename << "\n";

      QString fullPath = path+"/";

      // open
      StackOpen so(*this);
      so.run(stack, store, fullPath+filename);

      // add prefix to filename
      QString filenameNew = prefix + filename;

      // save
      StackSave ss(*this);
      ss.run(stack, store, fullPath+filenameNew, 9);

    }
    return true;
  }
  REGISTER_PROCESS(StackOpenSaveFolder);

  typedef Vector<3, vertex> Point3v;
  typedef std::map<Point2i, double> Point2iDoubleMap;
  typedef std::pair<Point2i, double> Point2iDoublePair;
  typedef std::map<Point3v, std::set<int> > Point3vIntSetMap;
  typedef std::pair<Point3v, std::set<int> > Point3vIntSetPair;
  typedef std::map<Point3d, vertex> Point3dVtxMap;
  typedef std::pair<Point3d, vertex> Point3dVtxPair;
  typedef std::map<Point2i, VtxSet> Point2iVtxSetMap;
  typedef std::pair<Point2i, VtxSet> Point2iVtxSetPair;
  
  // Put triangle in standard form
  Point3v triIndex(Point3v t)
  {
    if(t.x() > t.y())
      std::swap(t.x(), t.y());
    if(t.y() > t.z())
      std::swap(t.y(), t.z());
    if(t.x() > t.y())
      std::swap(t.x(), t.y());
    return t;
  }
  
  // Get vertex index, combine positions closer that the tolerance
  const vertex &vIndex(const vertex &v, Point3dVtxMap &vMap, double tolerance)
  {
    Point3d pos = v->pos;
    if(vMap.count(pos) > 0)
      return vMap[pos];
    else if(vMap.empty() or tolerance <= 0)
      return (vMap[pos] = v);
    else {
      Point3d p = pos;
      p.x() -= tolerance;
      Point3dVtxMap::iterator it = vMap.lower_bound(p);
      if(it == vMap.end())
        --it;
      while(it != vMap.end() and it->first.x() - pos.x() <= tolerance) {
        if(norm(pos - it->first) <= tolerance)
          return it->second;
        ++it;
      }
      return (vMap[pos] = v);
    }
  }
  
  // Compute the cell graph and 3D position/volume of the cells.
  bool CellGraph3D::run(Mesh *mesh, const QString &fileName, double minArea, double tolerance)
  {
    vvGraph &S = mesh->graph();
  
    // Get cells by connected regions
    VVGraphVec cellVertex;
    VtxIntMap vertexCell;
    mesh->getConnectedRegions(S, cellVertex, vertexCell);
  
    // Build cell to label map and check labeling
    IntIntMap cellLabel;
    for(uint cell = 0; cell < cellVertex.size(); ++cell) {
      int label = 0;
      for(const vertex &v : cellVertex[cell])
        if(label == 0)
          label = v->label;
        else if(label != v->label)
          throw QString("%1: The mesh is not labeled correctly.").arg(name());
  
      if(label == 0)
        throw QString("%1: the mesh has unlabeled cell(s).").arg(name());
  
      cellLabel[cell] = label;
    }
  
    // Find the volume and centers of the cells and create triangle map for cell graph
    IntDoubleMap cellVolume;
    IntDoubleMap cellArea;
    IntPoint3dMap cellCenter;
    Point3dVtxMap vMap;
    Point3vIntSetMap triCell; // This is set of cells for each triangle
    for(const vertex &v : S) {
      if(!progressAdvance())
        userCancel();
      forall(const vertex &n, S.neighbors(v)) {
        const vertex& m = S.nextTo(v,n);
        if(!S.uniqueTri(v,n,m)) // Only once per triangle in the mesh
          continue;
        // If the labels are not the same, this means we have more than one label on the same cell
        if(v->label != n->label or v->label != m->label)
          throw QString("%1: The mesh has more than one label on a cell.").arg(name());
        int cell = vertexCell[v];
  
        double vol = signedTetraVolume(v->pos, n->pos, m->pos);
        cellVolume[cell] += vol;
        cellArea[cell] += triangleArea(v->pos, n->pos, m->pos);
        cellCenter[cell] += vol * (v->pos + n->pos + m->pos)/ 4.0;
  
  	    // Add the cell to the triangle map, points with similar positions will get the same index
  	    triCell[triIndex(Point3v(vIndex(v, vMap, tolerance), vIndex(n, vMap, tolerance), 
                               vIndex(m, vMap, tolerance)))].insert(cell);
      }
    }
  
    // Now calculate wall areas from triangle list
    Point2iDoubleMap wallArea;
    Point2iDoubleMap wallSignal;
    for(const Point3vIntSetPair &pr : triCell) {
      if(!progressAdvance())
        userCancel();
      if(pr.second.size() > 2) {
        Information::out << "Error, triangle belongs to more than 2 cells:";
        for(int cell : pr.second)
          Information::out << " " << cell;
        Information::out << " Area:" << triangleArea(pr.first.x()->pos, pr.first.y()->pos, pr.first.z()->pos) << endl;
        continue;
          
        //throw QString("Error, triangle belongs to more that 2 cells.");
      } else if(pr.second.size() != 2)
        continue;
  
      IntSet::iterator cell1 = pr.second.begin();
      IntSet::iterator cell2 = cell1;
      cell2++;
      double area = triangleArea(pr.first.x()->pos, pr.first.y()->pos, pr.first.z()->pos);
      wallArea[Point2i(*cell1, *cell2)] += area;
      wallArea[Point2i(*cell2, *cell1)] += area;
      double signal = (pr.first.x()->signal + pr.first.y() + pr.first.z()->signal)/3.0 * area;
      wallSignal[Point2i(*cell1, *cell2)] += signal;
      wallSignal[Point2i(*cell2, *cell1)] += signal;
    }
    std::cout << "wallArea size:" << wallArea.size() << std::endl;

    // Calculate standard deviation
    Point2iDoubleMap wallStdDev;
    for(const Point3vIntSetPair &pr : triCell) {
      if(!progressAdvance())
        userCancel();
      if(pr.second.size() > 2) {
        Information::out << "Error, triangle belongs to more than 2 cells:";
        for(int cell : pr.second)
          Information::out << " " << cell;
        Information::out << " Area:" << triangleArea(pr.first.x()->pos, pr.first.y()->pos, pr.first.z()->pos) << endl;
        continue;
          
        //throw QString("Error, triangle belongs to more that 2 cells.");
      } else if(pr.second.size() != 2)  
        continue; // Outside triangle

      IntSet::iterator cell1 = pr.second.begin();
      IntSet::iterator cell2 = cell1;
      cell2++;
      Point2i wall1(*cell1, *cell2);
      Point2i wall2(*cell2, *cell1);
      double area = triangleArea(pr.first.x()->pos, pr.first.y()->pos, pr.first.z()->pos);
      double signal = (pr.first.x()->signal + pr.first.y() + pr.first.z()->signal)/3.0;
      wallStdDev[wall1] += fabs(signal - wallSignal[wall1]/wallArea[wall1]) * area;
      wallStdDev[wall2] += fabs(signal - wallSignal[wall2]/wallArea[wall2]) * area;
    }

    // Calculate Moran's I
    Point2iDoubleMap wallMorans;
    Point2iVtxSetMap wallVs;
    for(const Point3vIntSetPair &pr : triCell) {
      if(pr.second.size() != 2)
        continue;
      IntSet::iterator cell1 = pr.second.begin();
      IntSet::iterator cell2 = cell1;
      cell2++;
      Point2i wall1(*cell1, *cell2);
      Point2i wall2(*cell2, *cell1);

      for(int i = 0; i < 3; i++) {
        wallVs[wall1].insert(pr.first[i]);
        wallVs[wall2].insert(pr.first[i]);
      }
    }
    for(const Point2iVtxSetPair &pr : wallVs) {
      double N = 0, W = 0, numer = 0, denom = 0;
      double mean = wallSignal[pr.first]/wallArea[pr.first];
      for(const vertex &v : pr.second) {
        N += 1;
        forall(const vertex &n, S.neighbors(v)) {
          if(pr.second.count(n) == 0)
            continue;
          W += 1;
          numer += (v->signal - mean) * (n->signal - mean);
        }
        denom += (v->signal - mean) * (v->signal - mean);
      }
      wallMorans[pr.first] = (N * numer) / (W * denom);
    }
 
    // Calculate centers, update the heat map min/max and mormalize
    double vol_min = std::numeric_limits<double>::max();
    double vol_max = 0;
  
    // Go over all the volumes to find the min and max, and calculate the centers
    for(const IntDoublePair& vl : cellVolume) {
      if(vl.second != 0)
        cellCenter[vl.first] /= vl.second;
      if(vl.second < vol_min) 
        vol_min = vl.second;
      if(vl.second > vol_max) 
        vol_max = vl.second;
    }
  
    // Write data to .csv file
    QFile file(fileName);
    if (!file.open(QIODevice::WriteOnly | QIODevice::Text))
      return setErrorMessage("Error:Cannot open output file:" + fileName);
    QTextStream csvout(&file);
  
    // Write out cell data 
    csvout << "Cell Data:" << endl ;
    csvout << "Label, Area (" << UM2 << "), Volume (" << UM3 << "), CenterX (" 
        << UM << "), CenterY (" << UM << "), CenterZ (" << UM << ")" << endl;
  
    for(const IntDoublePair &pr : cellVolume) 
      csvout << cellLabel[pr.first] << ", " << cellArea[pr.first] << ", " 
             << pr.second << ", " << cellCenter[pr.first].x() << ", " << cellCenter[pr.first].y() << ", " 
             << cellCenter[pr.first].z() << endl;
    csvout << endl << endl;
  
    // Write out wall data 
    csvout << "Cell Wall Data:" << endl ;
    csvout << "Cell, Cell, Area (" << UM2 << "), Total Signal, Avg Signal, Std Dev, Coeff Variance, Distance(" << UM << "), Moran's I" << endl;
  
    for(const Point2iDoublePair &pr : wallArea) 
      if(minArea == 0 or pr.second >= minArea)
        csvout << cellLabel[pr.first.x()] << ", " << cellLabel[pr.first.y()] << ", " << pr.second 
               << ", " << wallSignal[pr.first] << ", " << wallSignal[pr.first]/pr.second 
               << ", " << wallStdDev[pr.first]/pr.second << ", " << wallStdDev[pr.first]/wallSignal[pr.first]
               << ", " << norm(cellCenter[pr.first.x()] - cellCenter[pr.first.y()])
               << ", " << wallMorans[pr.first] << endl;
  
    // Generate heat-map, normalize, and set the bounds in the mesh
    auto &labelHeat = mesh->labelHeat();
    labelHeat.clear();
    for(const Point2iDoublePair &pr : wallArea) 
      labelHeat[cellLabel[pr.first.x()]] += 1.0;
  
    double minNbs = 1e20, maxNbs = 0;
    for(const IntDoublePair &pr : labelHeat) {
      if(maxNbs < pr.second)
        maxNbs = pr.second;
      if(minNbs > pr.second)
        minNbs = pr.second;
    }
    // Set the map bounds. A heat of 0 will be interpreted has vol_min, and a heat of 1 as vol_max.
    mesh->heatMapBounds() = Point2f(minNbs, maxNbs);
    // Unit of the heat map
    mesh->heatMapUnit() = "Neighbors";
    // Tell the system the triangles have changed (here, the heat map)
    mesh->updateTriangles();
  
    return true;
  }
  
  REGISTER_PROCESS(CellGraph3D);
}
