//
// 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 <Mesh.hpp>
#include <ImageData.hpp>
#include <Stack.hpp>
#include <Information.hpp>
#include <GraphUtils.hpp>
#include <Dir.hpp>

#include <algorithm>
#include <iterator>
#include <limits>
#include <QFileInfo>

namespace mgx 
{
  // From SystemProcessLoad
  QList<int> extractVersion(QIODevice &file); 

//  Mesh::Mesh(const Stack* stack) : _id(-1), _stack(stack)
//  {
//    init();
//  }
  
  Mesh::Mesh(int id, const Stack* stack)
    : _id(id), _stack(stack)
  {
    init();
  }
//  
//  void Mesh::setStack(const Stack* s) {
//    _stack = s;
//  }
  
  void Mesh::init()
  {
    changed_surf_function = false;
    changed_heat_function = false;
    _changed = 0;
    _surfView = SURF_VERTEX;
    _surfVertexView = "Signal";
    _surfTriangleView = "None"; 
    _surfLabelView = "Label";
    _culling = true;
    _blending = true;
    _showSurface = true;
    _showMesh = false;
    _opacity = 1.0f;
    _brightness = 1.0f;
    _meshView = "Selected";
    _showMeshLines = false;
    _showMeshPoints = false;
    _showMeshCellMap = false;
    _useParents = false;
    _imgtex = false;
    _scaled = false;
    _transformed = false;
    _showBBox = false;
  }
  
  Mesh::~Mesh() {}
  
  void Mesh::updateAll() 
  {
    _changed = ImgData::RELOAD_VBO;
  }
  
  void Mesh::updateTriangles() 
  {
    _changed |= ImgData::RELOAD_TRIS;
  }
  
  void Mesh::updateLines() 
  {
    _changed |= ImgData::RELOAD_LINES;
  }
  
  void Mesh::updatePositions() 
  {
    _changed |= ImgData::RELOAD_POS;
  }
  
  void Mesh::updateSelection() 
  {
    _changed |= ImgData::UPDATE_SELECTION;
  }

  bool Mesh::setAxisColor(const Colorb &pos, const Colorb &neg, const SymmetricTensor &s, Vec3Colorb &color)
  {
    for(int i = 0; i < 3; ++i)
      if(s.evals()[i] < 0)
        color[i] = neg;
      else
        color[i] = pos;

    return true;
  }

  bool Mesh::setCellAxisColor(const Colorb &pos, const Colorb &neg)
  {
    IntVec3ColorbAttr &color = cellAxisColor();
    color.clear();
    
    forall(const IntSymTensorPair &pr, cellAxis())
      setAxisColor(pos, neg, pr.second, color[pr.first]);

    return true;
  }

  void Mesh::clearCellAxis()
  {
    cellAxis().clear();
    cellAxisVis().clear();
    cellAxisColor().clear();
    cellAxisType().clear();
  }

  bool Mesh::setTriangleAxisColor(const Colorb &pos, const Colorb &neg)
  {
    TriVec3ColorbAttr &color = triangleAxisColor();
    color.clear();
    
    forall(const TriSymTensorPair &pr, triangleAxis())
      setAxisColor(pos, neg, pr.second, color[pr.first]);

    return true;
  }

  bool Mesh::setVertexAxisColor(const Colorb &pos, const Colorb &neg)
  {
    VtxVec3ColorbAttr &color = vertexAxisColor();
    color.clear();
    
    forall(const VtxSymTensorPair &pr, vertexAxis())
      setAxisColor(pos, neg, pr.second, color[pr.first]);

    return true;
  }

  void Mesh::updateWallGeometry(float borderSize)
  {
    // First, mark the borders
    markBorder(borderSize);
  
    _labelNeighbors.clear();
    _wallVId.clear();
  
    vvGraph &S = graph();
    forall(const vertex& v, S) {
      if(v->label != -1)
        continue;
      forall(const vertex& n, S.neighbors(v)) {
        if(n->label != -1 or n > v)
          continue;
        IntSet Sv, Sn;
        // Get list of neighbors for each vertex.
        forall(const vertex& s, S.neighbors(v))
          if(s->label != -1)
            Sv.insert(s->label);
        forall(const vertex& s, S.neighbors(n))
          if(s->label != -1)
            Sn.insert(s->label);
        IntVec labelV;
        forall(int lab, Sv)
          if(Sn.find(lab) != Sn.end())
            labelV.push_back(lab);
        if(labelV.size() <= 1)
          continue;
        // Find wall length, add to neighbors array, and add vertices in wall to set
        for(size_t i = 0; i < labelV.size() - 1; i++) {
          IntIntPair wall;
  
          wall = IntIntPair(labelV[i], labelV[i + 1]);
          if(wall.first > 0) {
            wallGeom()[wall] += (v->pos - n->pos).norm();
            _labelNeighbors[wall.first].insert(wall.second);
            _wallVId[wall].insert(v.id());
            _wallVId[wall].insert(n.id());
          }
  
          wall = IntIntPair(labelV[i + 1], labelV[i]);
          if(wall.first > 0) {
            wallGeom()[wall] += (v->pos - n->pos).norm();
            _labelNeighbors[wall.first].insert(wall.second);
            _wallVId[wall].insert(v.id());
            _wallVId[wall].insert(n.id());
          }
        }
      }
    }
  }

  int Mesh::getLabel(const vertex &v1, const vertex &v2, const vertex &v3, 
                                                         const IntIntAttr &parents)
  {
    int label = mgx::getLabel(v1, v2, v3);

    return getLabel(label, parents);
  }
 
  int Mesh::getLabel(int label, const IntIntAttr &parents)
  {
    if(useParents())
      return parents[label];

    return label;
  }
   
  bool Mesh::isBordTriangle(int label, vertex v, vertex n, vertex m, IntIntPair &wall) const
  {
    if(_wallVId.empty() or v->minb == 0 or n->minb == 0 or m->minb == 0)
      return (false);
  
    // Search through neighbor list and sum signal and area
    IntIntSetMap::const_iterator found = _labelNeighbors.find(label);
    if(found == _labelNeighbors.end())
      return false;
    forall(int lab, found->second) {
      IntIntPair p(label, lab);
      IntIntVIdSetMap::const_iterator found = _wallVId.find(p);
      if(found == _wallVId.end())
        continue;
  
      const VIdSet& VtxSet = found->second;
  
      // If triangle belongs to this side of wall
      int count = 0;
      if(VtxSet.find(v->minb) != VtxSet.end())
        count++;
      if(VtxSet.find(n->minb) != VtxSet.end())
        count++;
      if(VtxSet.find(m->minb) != VtxSet.end())
        count++;
  
      if(count >= 2) {
        wall = p;
        return true;
      }
    }
    return false;
  }
  
  void Mesh::markBorder(float borderSize)
  {
    vvGraph &S = graph();
    if(S.empty())
      return;
  
    // For cell mesh we just need to mark all vertices that are not centers
    if(meshType() == "MGXC") {
      int bcount = 0;
      forall(const vertex& v, S)
        if(v->type == 'j') {
          v->minb = v.id();
          bcount++;
        } else
          v->minb = S.anyIn(v).id();
    } else if(meshType() == "MGXM") {
      // Process non cellular mesh, first all border vertices in R
      std::set<vertex> R;
      forall(const vertex& v, S)
        if(v->label == -1) {
          R.insert(v);
          v->minb = v.id();
        } else
          v->minb = 0;
      // uint bsize = R.size();
  
      // Start with border vertices and propogate
      std::set<vertex> T = R;     // Temp list
      std::set<vertex> D;         // delete list
      std::set<vertex> I;         // insert list
  
      uint rsize;
      do {
        rsize = R.size();
        forall(const vertex& v, T) {
          // If a vertex doesn't add any neighbors, we'll assume it is interior
          bool del = true;
          forall(const vertex& n, S.neighbors(v)) {
            // If not in border front
            if(R.find(n) == R.end()) {
              vertex w = vertex(v->minb);
              if(w.isNull()) {
                Information::out << "markBorder::Error:Null vertex w" << endl;
                continue;
              }
              vertex minb = w;
              double min = (n->pos - w->pos).norm();
              // forall(const vertex &m, S.neighbors(w))
              forall(const vertex& m, S.neighbors(n)) {
                if(m->minb != 0) {
                  vertex mb(m->minb);
                  if(mb.isNull()) {
                    Information::out << "markBorder::Error:Null vertex mb" << endl;
                    continue;
                  }
                  double dis = (n->pos - mb->pos).norm();
                  if(dis < min) {
                    min = dis;
                    minb = mb;
                  }
                }
              }
              if(min <= borderSize) {
                n->minb = minb.id();
                R.insert(n);
                I.insert(n);
              }
              del = false;
            }
          }
          if(del)
            D.insert(v);
        }
        forall(const vertex& v, I)
          T.insert(v);
        forall(const vertex& v, D)
          T.erase(v);
        I.clear();
        D.clear();
      } while(rsize < R.size());     // Do till R stops growing
    }
  }
  
  void Mesh::reset()
  {
    T.S.clear();
    T.C.clear();
    _attributes.clear();

    _labelNeighbors.clear();
    _wallVId.clear();
    _VVCorrespondence.clear();
  
    setFile();
    clearImgTex();
    setScaled(true);
    setTransformed(false);
    setMeshType("MGXM");
  }
  
  // Set the mesh filename
  void Mesh::setFile(const QString& file)
  {
    if(file.isEmpty()) {
      _file = file;
    } else {
      _file = absoluteFilePath(file);
    }
  }

  // Get set of vertices that are connected
  void Mesh::getConnectedVertices(const vvGraph &S, vvGraph &V) const
  {
    vvGraph vNew(V);
    do {
      vvGraph N;
      forall(const vertex& v, vNew)
        forall(const vertex& n, S.neighbors(v))
          if(not V.contains(n)) {
            V.insert(n);
            N.insert(n);
          }
      vNew = N;
      N.clear();
    } while(!vNew.empty());
  }
  
  // Get list of connected regions. Most often used to find 3D cells.
  // cellV is a vector of vvGraphs representing all the regions (cells)
  // vCell maps vertices to region (cell) index numbers
  int Mesh::getConnectedRegions(const vvGraph& S, VVGraphVec& cellV, VtxIntMap& vCell) const
  {
    int count = 0;
    vvGraph T(S.begin(), S.end(), false);
    while(!T.empty()) {
      vvGraph V;
      V.insert(*T.begin());
      getConnectedVertices(S, V);
      forall(const vertex& v, V) {
        vCell[v] = count;
        T.erase(v);
      }
      cellV.push_back(V);
      count++;
    }
    // Return number of cells
    return count;
  }
  
  // Get the neighborhood of a vertex, Euclidean distance, must be connected
  void Mesh::getNhbd(const vvGraph& S, const vertex& v, float radius, VtxSet& nhbd) // static
  {
    VtxVec VLast, VNew;
    nhbd.insert(v);
    VLast.push_back(v);
    Point3d pos = v->pos;
    do {
      forall(const vertex& w, VLast)
        forall(const vertex& n, S.neighbors(w))
          if(nhbd.find(n) == nhbd.end() and norm(n->pos - pos) < radius) {
            VNew.push_back(n);
            nhbd.insert(n);
          }
  
      std::swap(VLast, VNew);
      VNew.clear();
    } while(!VLast.empty());
  }
  
  std::vector<vertex> Mesh::selectedVertices() const
  {
    const vvGraph &S = graph();
    std::vector<vertex> vs;
    if(showMesh() and (showMeshLines() or showMeshPoints())) {
      forall(const vertex &v, S)
        if(v->selected)
          vs.push_back(v);
    }

    return vs;
  }

  int Mesh::selectedCount() const
  {
    int count = 0;
    const vvGraph &S = graph();
    forall(const vertex &v, S)
      if(v->selected)
        count++;

    return count;
  }
  
  std::vector<vertex> Mesh::activeVertices() const
  {
    const vvGraph &S = graph();
    std::vector<vertex> vs;
    vs.reserve(S.size());
    if(showMesh() and (showMeshLines() or showMeshPoints())) {
      forall(const vertex& v, S) {
        if(v->selected)
          vs.push_back(v);
      }
    }
    if(vs.empty())
      std::copy(S.begin(), S.end(), std::back_inserter(vs));
    return vs;
  }

  // Clear the selection
  bool Mesh::clearSelection()
  {
    vvGraph &S = graph();

    if(S.empty())
      return false;

    bool changed = false;

    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      if(v->selected) {
        v->selected = false;
        changed = true;
      }
    }

    return changed;
  }

  // Correct the selection for 2.5D labelled meshes. This makes sure cell boundaries
  // are not deleted and that the boundaries between cells are labeled correctly.
  void Mesh::correctSelection(bool inclusive)
  {
    vvGraph &S = graph();
    if(S.empty())
      return;

    std::vector<uchar> result(S.size()); // std::vector<bool> not thread safe
    #pragma omp parallel for
    for(size_t i = 0; i < S.size(); ++i) {
      vertex v = S[i];
      if(v->label < 0 and v->type != 'l') {
        bool has_non_border = false;
        bool has_nb_selected = false;
        bool has_nb_non_selected = false;
        bool has_border_selected = false;
				bool has_nb_labeled = false;
        forall(const vertex &n, S.neighbors(v)) {
          if(n->label >= 0) {
            has_non_border = true;
            if(n->label > 0)
              has_nb_labeled = true;
            if(n->selected)
              has_nb_selected = true;
            else
              has_nb_non_selected = true;
          } else if(n->selected)
            has_border_selected = true;
        }
        if(inclusive and has_nb_selected)
          result[i] = true;
        else if(has_non_border and not has_nb_non_selected)
          result[i] = true;
        else if(!has_non_border and has_border_selected)
          result[i] = true;
        else if(v->selected and !has_nb_labeled)
          result[i] = true;
        else
          result[i] = false;
      } else
        result[i] = v->selected;
    }

    #pragma omp parallel for
    for(size_t i = 0; i < S.size(); ++i)
      S[i]->selected = result[i];
  }
  
  void Mesh::resetModified()
  {
    changed_surf_function = false;
    changed_heat_function = false;
    _changed = 0;
  }
  
  void Mesh::updateBBox()
  {
    _bbox = BoundingBox3f(Point3f(FLT_MAX, FLT_MAX, FLT_MAX),
                                             Point3f(-FLT_MAX, -FLT_MAX, -FLT_MAX));
    vvGraph &S = graph();
    forall(const vertex& v, S)
      _bbox |= Point3f(v->pos);
  }

  static bool addTriangle(vvGraph& S, const std::vector<vertex> &vertices, const Point3i &tri)
  {
    vertex x = vertices[tri.x()];
    vertex y = vertices[tri.y()];
    vertex z = vertices[tri.z()];
    if(S.valence(x) == 0) {
      S.insertEdge(x, y);
      S.spliceAfter(x, y, z);
      return true;
    } else {
      bool yinx = S.edge(x, y);
      bool zinx = S.edge(x, z);
  
      if(yinx and zinx) {     // Triangle already there
        if(S.nextTo(x, y) != z) {
          Information::out << "addTriangle::Error in nhbd" << endl;
          return false;
        }
        return true;
      } else if(yinx) {
        S.spliceAfter(x, y, z);
        return true;
      } else if(zinx) {
        S.spliceBefore(x, z, y);
        return true;
      }
    }
    return false;   // Did not insert triangle
  }
  
  bool Mesh::meshFromTriangles(vvGraph& S, const std::vector<vertex> &vertices, 
                                const std::vector<Point3i> &triangles, bool checkUnique)
  {
    // Insert the vertices in the mesh
    S.insert(vertices.begin(), vertices.end(), checkUnique);

    // Make triangle list with all orientations.
    std::vector<Point3i> tris;
    tris.resize(3 * triangles.size());
    #pragma omp parallel for
    for(uint i = 0; i < triangles.size(); i++) {
      int idx = i * 3;
      const Point3i &tri = triangles[i];
      tris[idx++] = tri;
      tris[idx++] = Point3i(tri.y(), tri.z(), tri.x());
      tris[idx] = Point3i(tri.z(), tri.x(), tri.y());
    }

    // Extract neighborhood from triangles, make multiple passes until list empty
    while(!tris.empty()) {
      std::vector<Point3i> new_tris;
      new_tris.reserve(tris.size() / 4);
      bool success = false;
      for(size_t i = 0; i < tris.size(); i++) {
        if(!addTriangle(S, vertices, tris[i]))
          new_tris.push_back(tris[i]);
        else
          success = true;
      }
      if(!success)
        return false;       // Failed to add any triangle on this turn ... won't help again
      swap(tris, new_tris);
    }
    return true;
  }
  
  
  // Update the mesh labelCenter and labelNormal, to display the cell axis on mesh surface.
  // if parent labels are selected, the centers and normals are computed based on parents
  // if mesh is cellular, simply take the average positions and normals
  void Mesh::updateCentersNormals()
  {
		CellTissue &T = tissue();

    IntPoint3fAttr &centers = useParents() ? parentCenter() : labelCenter();
    IntPoint3fAttr &normals = useParents() ? parentNormal() : labelNormal();
    IntPoint3fAttr &centersVis = useParents() ? parentCenterVis() : labelCenterVis();
    IntPoint3fAttr &normalsVis = useParents() ? parentNormalVis() : labelNormalVis();

    AttrMap<int, Point3d>& cellCentersAttr = _attributes.attrMap<int, Point3d>("Measure Label Vector CellCenters");
    cellCentersAttr.clear();

    centers.clear();
    normals.clear();
    centersVis.clear();
    normalsVis.clear();
  
    std::map<int, std::set<vertex> > vertices;
  
    // Compute label centers and normals
    if(T.hasCellGraph()) { // MGX2D and MGX3D
      forall(const cell& c, T.C) {
        int label = c->label;
        centers[label] = Point3f(c->pos);
        normals[label] = Point3f(c->nrml);
        centersVis[label] = Point3f(c->pos);
        normalsVis[label] = Point3f(c->nrml);
      }
	  } else if(meshType() == "MGXC" and !useParents()) {
		  // If the mesh is simplified, then the centers have label > 0 
      forall(const vertex& v, T.S) {
        if(v->label > 0) {
          int label = v->label;
          centers[label] = Point3f(v->pos);
          normals[label] = Point3f(v->nrml);
          centersVis[label] = Point3f(v->pos);
          normalsVis[label] = Point3f(v->nrml);
          cellCentersAttr[label] = v->pos;
        }
      }
    } else if(meshType() == "MGXM" or (meshType() == "MGXC" and useParents())) { // For mgxm or mgxc with parents we need to compute it
      IntFloatMap areas;
      // Find cell centers, areas and normals.
      forall(const vertex& v, T.S) {
        forall(const vertex& n, T.S.neighbors(v)) {
          const vertex& m = T.S.nextTo(v, n);
          if(!T.S.uniqueTri(v, n, m))
            continue;
          int label = getLabel(v, n, m, parents());
          if(label == 0)
            continue;
          vertices[label].insert(v);
          vertices[label].insert(n);
          vertices[label].insert(m);
          float area = triangleArea(v->pos, n->pos, m->pos);
          areas[label] += area;
          centers[label] += Point3f((v->pos + n->pos + m->pos) / 3.0 * area);
          normals[label] += Point3f((v->nrml + n->nrml + m->nrml) / 3.0 * area);
        }
      }
  
      forall(const IntFloatPair& p, areas) {
        int label = p.first;
        float area = p.second;
        if(area > 0) {
          centers[label] /= area;
          normals[label] /= area;
          Point3f center = centers[label];
          cellCentersAttr[label] = Point3d(centers[label]);
  		    // Bump the centers to the surface for visualization
          vertex minv = *(vertices[label].begin());
          float mindis = norm(Point3f(minv->pos) - center);
          forall(const vertex& v, vertices[label]) {
            float dis = norm(Point3f(v->pos) - center);
            if(dis < mindis) {
              mindis = dis;
              minv = v;
            }
          }
          centersVis[label] = Point3f(minv->pos);
          normalsVis[label] = Point3f(minv->nrml);
        }
      }
    }
  }

  // Scale the heat map bounds to the existing range.
  Point2f Mesh::calcHeatMapBounds() 
  {
    IntFloatAttr &heat = labelHeat();
		Point2f bounds;

    bounds[0] = std::numeric_limits<float>::max();
    bounds[1] = -bounds[0];
    forall(const IntFloatPair &pr, heat) {
      if(bounds[0] > pr.second)
        bounds[0] = pr.second;
      if(bounds[1] < pr.second)
        bounds[1] = pr.second;
    }
		return bounds;
  } 

  // Read vertices into a vector
  bool Mesh::readVertices(QIODevice &file, VtxVec &V, bool transform)
  {
    // Get the size
    int vCnt;
    readFile(file, vCnt);
    progressStart(QString("Loading Mesh %1 Vertices").arg(userId()), vCnt);

    // Reserve space
    V.reserve(V.size() + vCnt);

    // Read in the vertices
    for(int i = 0; i < vCnt; ++i) {
      vertex v;
      V.push_back(v);

      // VertexData
      readFile(file, v->pos);
      if(transform)       // Transform the point if required
        v->pos = Point3d(stack()->frame().coordinatesOf(qglviewer::Vec(v->pos)));
      int cId = 0;
      readFile(file, cId);
      v->cId = cId;
      readFile(file, v->signal);
      readFile(file, v->label);
      readFile(file, v->color);
      readFile(file, v->saveId);
      readFile(file, v->type);
      readFile(file, v->selected);

      _vMap[v->saveId] = v;
      if(v->label > viewLabel())
        setLabel(v->label);
      if(!progressAdvance(i))
        throw(QString("Mesh load canceled"));
    }
    return true;
  }

  bool Mesh::readVertexList(QIODevice &file, vvGraph &S)
  {
    // Read in the  vertex list
    IntVec vIds;
    readFile(file, vIds);

    // Now convert saveIds to vertices and insert, they should be unique
    for(uint i = 0; i < vIds.size(); ++i) {
      IntVtxMap::const_iterator vIt = _vMap.find(vIds[i]);
      if(vIt == _vMap.end())
        throw QString("Invalid vertex id in cell list: %1").arg(vIds[i]);
      S.insert(vIt->second, false);
    }
    return true;
  }

  bool Mesh::readNhbds(QIODevice &file, vvGraph &S) 
  {
    int vCnt;
    readFile(file, vCnt);
    progressStart(QString("Loading Mesh %1 Neighborhoods").arg(userId()), vCnt);

    for(int i = 0; i < vCnt; i++) {
      if(file.atEnd())
        throw(QString("Premature end of file reading neighborhoods"));

      int vId, nCnt;
      readFile(file, vId);
      readFile(file, nCnt);
      
      IntVtxMap::const_iterator vIt = _vMap.find(vId);
      if(vIt == _vMap.end())
        throw QString("Invalid vertex id: %1").arg(vId);
      vertex v = vIt->second;
      vertex pn;
      for(int j = 0; j < nCnt; j++) {
        int nId;
        readFile(file, nId);
        IntVtxMap::const_iterator nIt = _vMap.find(nId);
        if(vIt == _vMap.end())
          throw QString("Invalid vertex id: %1").arg(nId);
        vertex n = nIt->second;

        if(S.valence(v) == 0)
          S.insertEdge(v, n);
        else
          S.spliceAfter(v, pn, n);

        // Read in EdgeData
        readFile(file, S.edge(v, n)->color);

        pn = n;
      }
      if(!progressAdvance(i))
        throw(QString("Mesh load canceled"));
    }
    return true;
  }

  bool Mesh::readCellMesh(QIODevice &file, bool transform)
  {
    // Get the graphs
    vvGraph &S = tissue().S;
    vvGraph &J = tissue().J;
    cellGraph &C = tissue().C;

    // Clear the vertex and cell maps. These are used to map saveIds to vertexes and cells
    // RSS If a mesh is added, and there are unloaded vertex/cell attributes, this will be problem
    _vMap.clear();
    _cMap.clear();

    // Read in the vertices
    VtxVec T;
    readVertices(file, T, transform);
    // Insert vertices into graph without checking for duplicates, this is much faster
    S.insert(T.begin(), T.end(), false);

    // Read in the neighborhoods for the main graph
    readNhbds(file, S);

    // Read in the vertex list and neighborhoods for the junction graph
    readVertexList(file, J);
    readNhbds(file, J); 

    // Read in the cells
    int cCnt;
    readFile(file, cCnt);
    Information::out << QString("Loading %1 cells for Mesh %2, Type: %3")
                         .arg(cCnt).arg(userId()).arg(meshType()) << endl;
    progressStart(QString("Reading Mesh %1 Cells").arg(userId()), 0, true);
    C.reserve(cCnt);
    for(int i = 0; i < cCnt; i++) {
      if(file.atEnd())
        throw(QString("Premature end of file reading cells"));
      cell c;
      C.insert(c);

      // CellData
      readFile(file, c->pos);
      if(transform)       // Transform the point if required
        c->pos = Point3d(stack()->frame().coordinatesOf(qglviewer::Vec(c->pos)));
      readFile(file, c->area);
      readFile(file, c->volume);
      int vId = 0;
      readFile(file, vId);
      if(vId >= 0)
        c->vId = _vMap[vId].id();
      else
        c->vId = 0;

      readFile(file, c->label);
      readFile(file, c->color);
      readFile(file, c->saveId);

      // Read in the vertex list and neighborhoods for the S graph for this cell
      readVertexList(file, c->S);
      readNhbds(file, c->S);
  
      // Add the cell to the map, and bump the label if required
      _cMap[c->saveId] = c;
      if(c->label > viewLabel())
        setLabel(c->label);
      if(!progressAdvance())
        throw(UserCancelException());
    }
    // Read the neighborhoods of the cell graph
    for(int i = 0; i < cCnt; i++) {
      if(file.atEnd())
        throw(QString("Premature end of file reading cell neighborhoods"));
  
      int cId, nCnt;
      readFile(file, cId);
      readFile(file, nCnt);
        
      IntCellMap::const_iterator cIt = _cMap.find(cId);
      if(cIt == _cMap.end())
        throw QString("Invalid cell id: %1").arg(cId);
      cell c = cIt->second;
      cell pn(0);
      for(int j = 0; j < nCnt; j++) {
        int nId;
        readFile(file, nId);
        IntCellMap::const_iterator nIt = _cMap.find(nId);
        if(cIt == _cMap.end())
          throw QString("Invalid cell id: %1").arg(nId);
        cell n = nIt->second;

        if(C.valence(c) == 0)
          C.insertEdge(c, n);
        else
          C.spliceAfter(c, pn, n);
        pn = n;

        // WallData
        readFile(file, C.edge(c, n)->length);
        readFile(file, C.edge(c, n)->area);
        readFile(file, C.edge(c, n)->color);
      }
      if(!progressAdvance())
        throw(UserCancelException());
    }

    // Remap the cell Ids in the main graph, cell graphs should all be 0
//    forall(const vertex &v, S) {
//      if(v->label > 0 and v->cId >= 0)
//        v->cId = _cMap[v->cId].id();
//      else
//        v->cId = 0;
//    }

    // Remap the cell Ids in all graphs
    forall(const IntVtxPair &pr, _vMap) {
      vertex v = pr.second; 
      if(v->label > 0 and v->cId >= 0)
        v->cId = _cMap[v->cId].id();
      else
        v->cId = 0;
    }

    return true;
  }

  bool Mesh::read(QIODevice &file, bool transform)
  {
    // Read the header
    char magic[5];
    char magicCheck[5] = "MGXM";
    file.read(magic, 5);
    if(strncmp(magic, magicCheck, 4))
      throw(QString("readCellMesh: Error, bad header"));

    // Get the revision numbers, these will change when the structure changes
    QList<int> version = extractVersion(file);
    int major = version[0];
    int minor = version[1];
    if(major < 2)
      throw(QString("readCellMesh: Error, bad major version number: %1").arg(major));
  
    // Start progress bar, we'll allow cancel so to catch errors and clean up
    progressStart(QString("Loading Mesh %1").arg(userId()), 0, true);
    Information::out << QString("Loading Mesh %1, version %2.%3")
                                       .arg(userId()).arg(major).arg(minor) << endl;
    try {
      // Read in the mesh type
      QString mType;
      readFile(file, mType);
      setMeshType(mType);

      readCellMesh(file, transform);

      // This version had cell complexes, skip them
      if(major == 2 and minor == 1) {
        size_t t;
        readFile(file, t);
      }

      attributes().read(file, &_vMap, &_cMap);

      // Read any attributes that are not live
      if(hasImgTex()) {
        QByteArray ba =
          attributes().attrMap<QString, QByteArray>("MeshQByteArray")["ImageTex"];
        if(ba.size() == 0)
          clearImgTex();
        else {
          QImage image;
          image.loadFromData(ba);
          Information::out << QString("Loaded image texture size %1 x %2.")
                                                .arg(image.width()).arg(image.height()) << endl;
          if(image.width() * image.height() > 0)
            setImgTex(image);
          else
            clearImgTex();
        } 
      }

    } catch(const QString &s) {
      reset();
      throw(s);
    } catch(UserCancelException &e) { 
      reset();
      throw(e);
    } catch(...) {
      reset();
      throw("Unknown exception loading mesh");
    }
    return true;
  }

  bool Mesh::writeVertices(QIODevice &file, const vvGraph &S, bool transform)
  {
    // Write the number of vertices
    int vCnt = S.size();
    writeFile(file, vCnt);
  
    // Write the vertices
    forall(const vertex &v, S) {
      Point3d pos = v->pos;
      if(transform)
        pos = Point3d(stack()->getFrame().inverseCoordinatesOf(qglviewer::Vec(pos)));

      // VertexData
      writeFile(file, pos);
      cell c = T.getCell(v);
      int cId = 0;
      if(c)
        cId = c->saveId;
      else
        cId = -1;
      writeFile(file, cId);
      writeFile(file, v->signal);
      writeFile(file, v->label); 
      writeFile(file, v->color); 
      writeFile(file, v->saveId);
      writeFile(file, v->type);
      writeFile(file, v->selected);

      if(!progressAdvance())
        throw(UserCancelException());
    }
    return true;
  }

  bool Mesh::writeVertexList(QIODevice &file, const vvGraph &S) 
  {
    IntVec vIds;
    vIds.resize(S.size());
    for(uint i = 0; i < S.size(); i++)
      vIds[i] = S[i]->saveId;
    writeFile(file, vIds);
    return true;
  }

  bool Mesh::writeNhbds(QIODevice &file, const vvGraph &S) 
  {
    // Write the number of neighborhoods
    int vCnt = S.size();
    writeFile(file, vCnt);

    // Write the neighborhoods
    forall(const vertex &v, S) {
      writeFile(file, v->saveId);
      int val = S.valence(v);
      writeFile(file, val);

      forall(const vertex& n, S.neighbors(v)) {
        writeFile(file, n->saveId);

        // Write EdgeData
        writeFile(file, S.edge(v, n)->color);
      }

      if(!progressAdvance())
        throw(UserCancelException());
    }
    return true;
  }

  bool Mesh::writeCellMesh(QIODevice &file, bool transform)
  {
    // Get the graphs
    CellTissue &T = tissue();
    vvGraph &S = T.S;
    vvGraph &J = T.J;
    cellGraph &C = T.C;

    // Make sure all vertices are in the all graph
    T.updateAllGraph();

    // Set the saveId on the vertices in S graph
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++)
      S[i]->saveId = i; 

    // Set the saveId on the cells
    for(uint i = 0; i < C.size(); i++)
      C[i]->saveId = i; 

    // Write the vertices from the S graph
    writeVertices(file, S, transform);

    // Write the neighborhoods of the main graph
    writeNhbds(file, S);

    // Write the vertex list and neighborhoods for the junction graph
    writeVertexList(file, J);
    writeNhbds(file, J);

    // Write the cell count
    int cCnt = C.size();
    writeFile(file, cCnt);
    
    // Write the cells
    forall(const cell &c, C) {
      Point3d pos = c->pos;
      if(transform)
        pos = Point3d(stack()->getFrame().inverseCoordinatesOf(qglviewer::Vec(pos)));
      writeFile(file, pos);
      writeFile(file, c->area); 
      writeFile(file, c->volume); 
      int vId = 0;
      vertex v = T.getVtx(c);
      if(v)
        vId = v->saveId;
      else {
        // Only a warning for 2D cells
        // Information::out << "writeCellMesh: Warning, cell has no vertex id" << endl;
        vId = -1;
      }
      writeFile(file, vId); 
      writeFile(file, c->label); 
      writeFile(file, c->color); 
      writeFile(file, c->saveId); 

      // Write the S graph for this cell
      writeVertexList(file, c->S);
      writeNhbds(file, c->S);
    
      if(!progressAdvance())
        throw(UserCancelException());
    }

    // Write the cell neighborhoods
    forall(const cell &c, C) {
      writeFile(file, c->saveId);
      int val = C.valence(c);
      writeFile(file, val);
      forall(const cell& n, C.neighbors(c)) {
        writeFile(file, n->saveId);
        // WallData
        writeFile(file, C.edge(c, n)->length);
        writeFile(file, C.edge(c, n)->area);
        writeFile(file, C.edge(c, n)->color);
      }

      if(!progressAdvance())
        throw(UserCancelException());
    }

    return true;
  }

  bool Mesh::write(QIODevice &file, bool transform)
  {
    // Start progress bar, we'll allow cancel so to catch errors and clean up
    progressStart(QString("Saving Mesh %1").arg(userId()), 0, true);
    Information::out << QString("Saving Mesh %1") .arg(userId()) << endl;

    try {
      // Write the header
      // RSS If this is changes, be sure to skip version 2.1 (see read())
      file.write("MGXM 2.0 ", 9);

      // Write the mesh type
      writeFile(file, meshType());

      // Write out the cellular mesh
      writeCellMesh(file, transform);

      // Add any attributes which are not live
      if(hasImgTex()) {
        QByteArray ba;
        QBuffer buffer(&ba);
        buffer.open(QIODevice::WriteOnly);
        imgTex().save(&buffer, "TIFF");  
        attributes().attrMap<QString, QByteArray>("MeshQByteArray")["ImageTex"] = ba; 
      }

      // Now write the attributes
      attributes().write(file);

    } catch(const QString &s) {
      // try to remove half-written file
      file.close();
      file.open(QIODevice::Truncate);
      // also need to restore All graph for 3D cells
      file.close();
      throw(s);
    } catch(UserCancelException &e) {
      file.close();
      file.open(QIODevice::Truncate);
      throw(e);
    } catch(...) {
      file.close();
      file.open(QIODevice::Truncate);
      throw("Unknown exception writing mesh");
    }
    return true;
  }
}
