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

#include <Progress.hpp>
#include <GraphUtils.hpp>
#include <Insert.hpp>

namespace mgx 
{
  bool TransformMesh::run(Mesh* mesh, Point3f translation, Point3f rotation, float angle, float scale)
  {
    Point4f dir = homogeneous(rotation);
    Matrix4f M = Matrix4f::rotation(dir, angle);
    M(0, 3) = translation.x();
    M(1, 3) = translation.y();
    M(2, 3) = translation.z();
    M(3, 3) = 1.f / scale;
    const std::vector<vertex>& vs = mesh->activeVertices();
    if(scale < 0 and vs.size() != mesh->graph().size())
      return setErrorMessage(
        "Error, the scaling can be negative only if all the vertices are transformed at the same time.");
    forall(const vertex& v, vs) {
      v->pos = Point3d(cartesian(M * homogeneous(Point3f(v->pos))));
    }
    if(scale < 0) {
      mesh->graph().reverse();
      mesh->updateTriangles();
    }
    setNormals(mesh->graph());
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(TransformMesh);
  
  bool ReverseMesh::run(Mesh* mesh)
  {
    if(mesh->activeVertices().size() != mesh->graph().size())
      return setErrorMessage("All vertices must be reversed at the same.");
  
    mesh->graph().reverse();
    setNormals(mesh->graph());
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(ReverseMesh);
  
  bool SmoothMesh::run(Mesh* mesh, uint passes, bool wallsOnly)
  {
    const std::vector<vertex>& vs = mesh->activeVertices();
    vvGraph& S = mesh->graph();
    progressStart(QString("Smoothing mesh %1").arg(mesh->userId()), passes);
    for(uint i = 0; i < passes; ++i) {
      if(!progressAdvance(i))
        userCancel();
      // We'll use the normal as a temp position to save space in the vertex struct
      forall(const vertex& v, vs) {
        double wt = 0;
        v->nrml = Point3d(0, 0, 0);
        forall(const vertex &n, S.neighbors(v)) {
          if(wallsOnly) {
            if(n->label < 0) {
              wt += 1;
              v->nrml += n->pos;
            }
          } else if(!v->margin) {
            vertex m = S.nextTo(v, n);
            float a = triangleArea(v->pos, n->pos, m->pos);
            wt += a;
            v->nrml += a * (n->pos + m->pos) / 2.0;
          } else if(n->margin) {
            float a = (n->pos - v->pos).norm();
            v->nrml += a * n->pos;
            wt += a;
          }
        }
        if(wt > 0)
          v->nrml /= wt;
        else
          v->nrml = v->pos;
      }
      forall(const vertex& v, vs)
        v->pos = v->pos * .5 + v->nrml * .5;
    }
  
    setNormals(S);
    mesh->clearImgTex();
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(SmoothMesh);
  
  bool ShrinkMesh::run(Mesh* mesh, float distance)
  {
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    forall(const vertex& v, vs)
      v->pos -= distance * v->nrml;
  
    setNormals(mesh->graph());
    mesh->clearImgTex();
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(ShrinkMesh);
  
  bool SubdivideMesh::run(Mesh* mesh)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGXM)").arg(mesh->meshType()));

    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    VtxPoint2fAttr& tex2d = mesh->texCoord2d();
    bool doTex2d = tex2d.size() > 0;
  
    progressStart(QString("Subdividing mesh %1").arg(mesh->userId()), 0);
    vvGraph SS = S.subgraph(vs);
    vvGraph T, T1;
    T = S;
    Insert<vvGraph> insert;
    std::vector<vertex> new_vertices;
  
    // Add vertex in middle of all existing edges
    forall(const vertex& n, SS) {
      if(!progressAdvance()) {
        S = T;
        userCancel();
      }
      forall(const vertex& m, T.neighbors(n)) {
        if(SS.contains(m) and m < n)
          continue;
        vertex v = insert(n, m, S);
        new_vertices.push_back(v);
        if(n->selected or m->selected)
          v->selected = true;
        v->margin = false;
        v->signal = (m->signal + n->signal) / 2.0;
        v->pos = (m->pos + n->pos) / 2.0;
        if(doTex2d)
          v->*tex2d = (m->*tex2d + n->*tex2d) / 2.0;
        if(m->label == n->label)
          v->label = m->label;
      }
    }
    progressStart("", new_vertices.size() + 1);
    T1 = S;
    // Now add new edges
    int i = 0;
    forall(const vertex& v, new_vertices) {
      ++i;
      if((i % 1000 == 0)and !progressAdvance(i)) {
        S = T;
        userCancel();
      }
      forall(const vertex& n, T1.neighbors(v)) {
        vertex m = T1.nextTo(v, n);
        vertex nv = T1.nextTo(n, v);
        vertex mv = T1.prevTo(m, v);
        vertex nm = T1.nextTo(nv, n);
        vertex mn = T1.prevTo(mv, m);
        if(nm != mn) {
          if(mn == nv)         // it's a border!
          {
            if(norm(mv->pos - v->pos) < norm(m->pos - nv->pos)) {
              S.spliceAfter(v, m, mv);
              S.spliceBefore(mv, m, v);
            } else {
              S.spliceBefore(nv, mv, m);
              S.spliceAfter(m, mv, nv);
            }
          }
          continue;
        }
        S.spliceBefore(nv, nm, mv);
        S.spliceAfter(mv, mn, nv);
      }
    }
  
    SETSTATUS("Mesh " << mesh->userId() << " - Subdivide triangles, total vertices:" << S.size());
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(SubdivideMesh);
  
  bool LoopSubdivisionMesh::run(Mesh* mesh)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGXM)").arg(mesh->meshType()));

    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    vvGraph OldS;
    if(vs.size() != S.size())
      OldS = S.subgraph(S);
    else
      OldS = S;
  
    size_t ss = OldS.size();
    size_t disp = ss / 50;
    VtxPoint2fAttr & tex2d = mesh->texCoord2d();
    bool doTex2d = tex2d.size() > 0;
  
    progressStart(QString("Subdividing mesh %1").arg(mesh->userId()), 2 * ss);
  
    Insert<vvGraph> insert;
    std::vector<vertex> new_vertices;
    std::vector<Point3d> old_pos;
    std::map<vertex, size_t> vertex_index;
    old_pos.reserve(OldS.size());
  
    for(size_t i = 0; i < OldS.size(); ++i) {
      vertex v = OldS[i];
      vertex_index[v] = i;
      old_pos[i] = v->pos;
    }
  
    // Add vertex in middle of all existing edges
    for(size_t i = 0; i < OldS.size(); ++i) {
      if((i % disp == 0)and not progressAdvance(i))
        userCancel();
      bool surface_v = false;     // The vertex is on the border of the surface
      const vertex& v = OldS[i];
      size_t n = OldS.valence(v);
      Point3d q;
  
      float w = 0.375 + 0.25 * std::cos(2.0 * M_PI / float(n));
      w *= w;
      w = 0.625 - w;
  
      forall(const vertex& u, OldS.neighbors(v)) {
        size_t j = vertex_index[u];
        bool boundary = true;
        if((OldS.nextTo(u, v) == OldS.prevTo(v, u))and (OldS.prevTo(u, v) == OldS.nextTo(v, u)))
          boundary = false;
        else {
          if(surface_v)
            q = 0.0;
          surface_v = true;
        }
        if((surface_v and boundary) or not surface_v)
          q += old_pos[j];
  
        // Only sub-divide every undirected edge once
        if(u < v)
          continue;
  
        vertex x = insert(u, v, S);
        new_vertices.push_back(x);
  
        // Set properties of the new vertex
        if(u->selected or v->selected)
          x->selected = true;
        x->margin = boundary;
        if(u->label == v->label)
          x->label = u->label;
        if(boundary) {
          x->pos = (old_pos[i] + old_pos[j]) / 2.0;
          if(doTex2d)
            x->*tex2d = (u->*tex2d + v->*tex2d) / 2.0;
          x->signal = (u->signal + v->signal) / 2.0;
        } else {
          const vertex& nv = OldS.nextTo(v, u);
          const vertex& pv = OldS.prevTo(v, u);
          size_t kn = vertex_index[nv];
          size_t kp = vertex_index[pv];
          x->pos = (old_pos[i] * 3.0 + old_pos[j] * 3.0 + old_pos[kn] + old_pos[kp]) / 8.0;
          if(doTex2d)
            x->*tex2d = (u->*tex2d * 3.0 + v->*tex2d * 3.0 + pv->*tex2d + nv->*tex2d) / 8.0;
          x->signal = (u->signal * 3.0 + v->signal * 3.0 + pv->signal + nv->signal) / 8.0;
        }
      }
      if(surface_v) {
        v->pos = 0.75 * v->pos + 0.125 * q;
      } else {
        v->pos = (1 - w) * v->pos + (w / n) * q;
      }
    }
  
    progressStart("", ss + new_vertices.size());
  
    // Connect the inside triangles
    size_t i = 0;
    disp = new_vertices.size() / 50;
    forall(const vertex& v, new_vertices) {
      if((i % disp == 0)and not progressAdvance(ss + i))
        userCancel();
      ++i;
      const vertex& a = S.anyIn(v);
      const vertex& b = S.nextTo(v, a);
      if(OldS.prevTo(a, b) == OldS.nextTo(b, a)) {
        S.spliceAfter(v, a, S.prevTo(a, v));
        S.spliceBefore(v, b, S.nextTo(b, v));
      }
      if(OldS.nextTo(a, b) == OldS.prevTo(b, a)) {
        S.spliceAfter(v, b, S.prevTo(b, v));
        S.spliceBefore(v, a, S.nextTo(a, v));
      }
    }
  
    SETSTATUS("Mesh " << mesh->userId() << " - Subdivide triangles, total vertices:" << S.size());
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(LoopSubdivisionMesh);
  
  bool AdaptiveSubdivideBorderMesh::run(Mesh* mesh, float cellMaxArea, float borderDist)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGXM)").arg(mesh->meshType()));

    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    // Find triangles on borders
    mesh->markBorder(borderDist);
  
    // Copy to temp graph
    vvGraph T = S;
    forall(vertex v, vs)
      forall(vertex n, T.neighbors(v)) {
        vertex m = T.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
  
        // Subdivide triangles on border and clear label
        if(v->minb == 0 or m->minb == 0 or n->minb == 0)
          continue;
        v->label = n->label = m->label = 0;
        float area = triangleArea(v->pos, n->pos, m->pos);
        if(area < cellMaxArea)
          continue;
  
        if(!subdivideBisect(S, v, n, m)) {
          setErrorMessage("Failed bisecting triangle");
          return false;
        }
      }
    SETSTATUS(S.size() << " total vertices");

    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(AdaptiveSubdivideBorderMesh);
  
  bool AdaptiveSubdivideSignalMesh::run(Mesh* mesh, float cellMaxAreaLow, float cellMaxAreaHigh)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGXM)").arg(mesh->meshType()));

    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    // Find max/min colors
    float maxc = 0, minc = HUGE_VAL;
    forall(const vertex& v, vs) {
      if(v->signal > maxc)
        maxc = v->signal;
      if(v->signal < minc)
        minc = v->signal;
    }
    float rangec = maxc - minc;
    // Return if no real signal difference
    if(fabs(rangec) < 1e-5) {
      setErrorMessage("subdivideTriangles::Error:No signal difference");
      return false;
    }
    // Copy to temp graph
    vvGraph T = S;
    forall(vertex v, vs)
      forall(vertex n, T.neighbors(v)) {
        vertex m = T.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
  
        // Find max signal on triangle
        float maxcol = v->signal;
        if(n->signal > maxcol)
          maxcol = n->signal;
        if(m->signal > maxcol)
          maxcol = m->signal;
  
        float val = ((maxcol - minc) / rangec);
        float area = triangleArea(v->pos, n->pos, m->pos);
        if(area < interpolate(cellMaxAreaLow, cellMaxAreaHigh, val))
          continue;
  
        if(!subdivideBisect(S, v, n, m)) {
          setErrorMessage("Failed bisecting triangle");
          return false;
        }
      }
    SETSTATUS(S.size() << " total vertices");

    mesh->updateAll();

    return true;
  }
  REGISTER_PROCESS(AdaptiveSubdivideSignalMesh);
  
  bool SubdivideBisectMesh::run(Mesh* mesh, float cellMaxArea)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGXM)").arg(mesh->meshType()));

    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    // Copy to temp graph
    vvGraph T = S;
    forall(vertex v, vs)
      forall(vertex n, T.neighbors(v)) {
        vertex m = T.nextTo(v, n);
        if(!S.uniqueTri(v, n, m))
          continue;
  
        float area = triangleArea(v->pos, n->pos, m->pos);
        if(area < cellMaxArea)
          continue;
  
        if(!subdivideBisect(S, v, n, m)) {
          setErrorMessage("Failed bisecting triangle");
          return false;
        }
      }
    SETSTATUS(S.size() << " total vertices");

    mesh->updateAll();

    return true;
  }
  REGISTER_PROCESS(SubdivideBisectMesh);

  typedef std::set<vertex> tri;


  void longestSide(const tri& t, vertex &a, vertex &b, vertex &c)
  {
    std::set<vertex>::iterator it = t.begin();
    vertex v1 = *it; it++;
    vertex v2 = *it; it++;
    vertex v3 = *it;
    a = v1; b = v2; c = v3;
    if(norm(v2->pos - v3->pos) > norm(a->pos - b->pos)) {
      a = v2; b = v3; c = v1;
    }
    if(norm(v3->pos - v1->pos) > norm(a->pos - b->pos)) {
      a = v3; b = v1; c = v2;
    }
  }

  void subdivideInsert(CellTissue& T, const tri& t, std::map<std::pair<vertex,vertex>, vertex>& splits)
  {

    // Find longest edge and insert vertex in all graphs
    vertex a(0), b(0), c(0);
    longestSide(t, a, b, c);

    // Create a vertex and set position
    vertex d;
    d->pos = (a->pos + b->pos)/2.0;
    d->selected = true;

    // Put in list of divided edges, leave if already processed
    std::pair<vertex,vertex> newP = std::make_pair(a,b);
    if(b<a) newP = std::make_pair(b,a);
    if(splits.find(newP) == splits.end())
      splits[newP] = d;
    else
      return;

    // Put in all graphs that contain it
    forall(const cell &ce, T.C)
      if(ce->S.contains(a) and ce->S.contains(b) and ce->S.edge(a, b)) {
        ce->S.insert(d);
        ce->S.insertEdge(d, a);
        ce->S.insertEdge(d, b);
        ce->S.replace(a, b, d);
        ce->S.replace(b, a, d);
      }


  }


  // Do triangle subdivsion from 1-4 triangles depending on how many edges are split
  void subdivideTri(CellTissue& T, const tri& t, std::map<std::pair<vertex,vertex>, vertex>& splits)
  {
    // Vertex which was inserted in longest wall (always split)
    vertex d(0);

    // Find longest edge and divide
    vertex a(0), a0(0), b(0), b0(0), c(0);
    longestSide(t, a0, b0, c);

    forall(const cell &ce, T.C) {
      a = a0; b = b0;
      if(!ce->S.contains(a) or !ce->S.contains(b) or !ce->S.contains(c) /*or !ce->S.edge(a, b) or !ce->S.edge(b, c) or !ce->S.edge(a, c)*/){
        //std::cout << "cell doesnt contain tri " << !ce->S.contains(a) << !ce->S.contains(b) << !ce->S.contains(c) << "/" << !ce->S.edge(a, b) << !ce->S.edge(b, c) << !ce->S.edge(a, c)<< std::endl;
        continue;
      }

      if(ce->S.edge(a,b)) {
        std::cout << "long wall not divided" << std::endl;
        continue;
      }
      std::pair<vertex,vertex> newP = std::make_pair(a,b);
      if(b<a) newP = std::make_pair(b,a);
      if(splits.find(newP) == splits.end())
        std::cout << "Error d not found" << std::endl;

      // Get vertex (d) in long wall
      d = splits[newP];
      // Order a-b-c correctly
      vertex e = c;
      if(!ce->S.edge(a,e)) {
        std::pair<vertex,vertex> newP2 = std::make_pair(a,e);
        if(e<a) newP2 = std::make_pair(e,a);
        if(splits.find(newP2) == splits.end())
          std::cout << "Error a-t not found" << std::endl;
        e = splits[newP2];
      }
      if(e == ce->S.prevTo(a, d)) {
        a = b0; b = a0;
      }
      // Check if d OK
      if(!ce->S.edge(a, d) or !ce->S.edge(b, d)) {
        std::cout << "Missing edge, d not between a and b, d:" << d->pos << " a:" << a->pos << " b:" << b->pos << std::endl;
        continue;
      }
      // Check pos of d
      if(norm((a->pos + b->pos)/2.0 - d->pos) > .00001) {
        std::cout << "d in wrong position" << std::endl;
        continue;
      }
      // This triangle is already divided
      if(ce->S.edge(c, d)) {
        std::cout << "Triangle already divided (c-d) edge exists" << std::endl;
        continue;
      }
      // If all three vertices inserted, then divide into 4
      if(!ce->S.edge(a, c) and !ce->S.edge(b, c)) {
        vertex ac = ce->S.nextTo(a, d);
        vertex bc = ce->S.prevTo(b, d);
        ce->S.spliceBefore(d, a, ac);
        ce->S.spliceAfter(d, b, bc);
        ce->S.spliceAfter(ac, a, d);
        ce->S.spliceBefore(ac, c, bc);
        ce->S.spliceBefore(bc, b, d);
        ce->S.spliceAfter(bc, c, ac);
      } else if(ce->S.edge(a, c) and ce->S.edge(b, c)) {
        ce->S.spliceAfter(c, a, d);
        ce->S.spliceBefore(d, a, c);
      } else if(ce->S.edge(a, c)) {
        ce->S.spliceAfter(c, a, d);
        ce->S.spliceBefore(d, a, c);
        vertex v = ce->S.nextTo(c, d);
        ce->S.spliceAfter(d, b, v);
        ce->S.spliceBefore(v, b, d);
      } else if(ce->S.edge(b, c)) {
        ce->S.spliceBefore(c, b, d);
        ce->S.spliceAfter(d, b, c);
        vertex v = ce->S.prevTo(c, d);
        ce->S.spliceBefore(d, a, v);
        ce->S.spliceAfter(v, a, d);
      } else
        std::cout << "Triangle division case not handled" << std::endl;
        
    }
    
  }


  // Subdivide triangles by bisection propogating split edge to all meshes
  void subdivideTris(CellTissue& T, std::set<tri> &tris)
  {
    bool done = false;
    while(!done) {
      done = true;
      forall(tri t, tris){
        vertex v1(0), v2(0), v3(0);
        longestSide(t,v1,v2,v3);
        forall(const cell& c, T.C){
          if(!c->S.contains(v1) || !c->S.contains(v2) || !c->S.contains(v3) || !c->S.edge(v1, v2) || !c->S.edge(v2, v3) || !c->S.edge(v1, v3))
            continue;
          forall(const cell &c2, T.C) {
            if(!c2->S.contains(v1) || !c2->S.contains(v2) || !c2->S.edge(v1, v2))
              continue;
            // First look on one side
            vertex d = c2->S.nextTo(v1, v2);
            if(d == c2->S.prevTo(v2, v1) and d != v3) {
              tri t2;
              t2.insert(v1);
              t2.insert(v2);
              t2.insert(d);
              t = t2;
              if(tris.find(t) == tris.end()) {
                tris.insert(t);
                done = false;
              }
            }
            // Now look on the other
            d = c2->S.prevTo(v1, v2);
            if(d == c2->S.nextTo(v2, v1) and d != v3) {
              tri t2;
              t2.insert(v1);
              t2.insert(v2);
              t2.insert(d);
              t = t2;
              if(tris.find(t) == tris.end()) {
                tris.insert(t);
                done = false;
              }
            }
          }
      }
    }
  }

    // Insert vertices between edges in all graphs
    std::map<std::pair<vertex,vertex>, vertex> splits;   // List of edges to split
    forall(const tri &t, tris) {
      subdivideInsert(T, t, splits);
    }

    // Connect new vertices
    forall(const tri &t, tris) 
      subdivideTri(T, t, splits);

  }

  bool SubdivideBisectMesh3D::run(Mesh* mesh, float cellMaxArea)
  {
    if(mesh->meshType() != "MGX3D")
      throw(QString("Cannot subdivide mesh type (%1), must be (MGX3D)").arg(mesh->meshType()));

    CellTissue& T = mesh->tissue();
    vvGraph& S = mesh->graph();

    const std::vector<vertex>& vs = mesh->activeVertices();
    bool keepSel = vs.size() != S.size();

    forall(const vertex& v, vs){
      v->selected = true;
    }

    std::set<tri> tris; // set of all its (unique) tris

    // go through all unique triangles with area greate the maxArea
    forall(const cell& c, T.C){
      forall(const vertex& v, c->S){
        if(!v->selected) continue;
        forall(const vertex& n, c->S.neighbors(v)){
          if(!n->selected) continue;
          vertex m = c->S.nextTo(v,n);
          if(!c->S.uniqueTri(v,n,m)) continue;
          if(!m->selected) continue;
          // size too big?
          if(triangleArea(v->pos, n->pos, m->pos) < cellMaxArea) continue;
          tri t;
          t.insert(v);
          t.insert(n);
          t.insert(m);
          tris.insert(t);
        }
      }
    }

    subdivideTris(T, tris);

    T.updateAllGraph();
    T.updGeometry();

    if(!keepSel) {
      forall(const vertex& v, S)
        v->selected = false;
    }

    mesh->updateAll();

    return true;
  }
  REGISTER_PROCESS(SubdivideBisectMesh3D);
  
  bool ScaleMesh::run(Mesh* mesh, float scaleX, float scaleY, float scaleZ)
  {
    vvGraph& S = mesh->graph();
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    if(scaleX * scaleY * scaleZ < 0) {
      if(vs.size() != S.size())
        return setErrorMessage(
          "Error, symmetries have to be applied globally. Make sure the product of the scalings are positive.");
      S.reverse();
    }
  
    // Copy to temp graph
    forall(vertex v, vs) {
      if(scaleX != 1.0)
        v->pos.x() *= scaleX;
      if(scaleY != 1.0)
        v->pos.y() *= scaleY;
      if(scaleZ != 1.0)
        v->pos.z() *= scaleZ;
    }
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(ScaleMesh);
  
  bool MeshDeleteValence::run(Mesh* mesh, int startValence, int endValence)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot delete vertices in mesh type (%1), must be (MGXM)")
                                                                     .arg(mesh->meshType()));

    if(startValence < 0 or endValence < 0)
      return setErrorMessage("Start and end valence must be >= 0.");
    vvGraph& S = mesh->graph();
    std::vector<vertex> D;
    forall(const vertex& v, S)
      if(S.valence(v) >= (size_t)startValence and S.valence(v) <= (size_t)endValence)
        D.push_back(v);
    forall(const vertex& v, D)
      S.erase(v);
  
    SETSTATUS(D.size() << " vertex(es) deleted, total vertices:" << S.size());

    mesh->updateAll();

    return true;
  }
  REGISTER_PROCESS(MeshDeleteValence);

  bool MergeVertices::run(Mesh* mesh)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot merge vertices in mesh type (%1), must be (MGXM)")
                                                                     .arg(mesh->meshType()));
  
    vvGraph& S = mesh->graph();
    const VtxVec &vs = mesh->activeVertices();
    if(vs.size() < 2)
      throw QString("Error, Must select at least two connected vertices.");
    if(mesh->meshType() != "MGXM")
      throw(QString("Do not know how to delete vertices in mesh type %1").arg(mesh->meshType()));

    progressStart("Merging vertices", vs.size() - 1); // Correct if only one area
    int count = mergeVertices(S, vs);
    if(count > 0)
      mesh->updateAll();
  
    Information::out << count << " vertices deleted from mesh" << endl;
    return true;
  }
  REGISTER_PROCESS(MergeVertices);
 
  bool SplitEdges::run(Mesh* mesh, double maxDist, Subdivide *sDiv)
  {
    vvGraph& S = mesh->graph();

    // Does not yet work on 3D meshes
    if(mesh->meshType() == "MGX3D")
      throw(QString("Cannot split edges in 3D cell tissue (MGX3D)"));

    // Split the edges
    std::vector<VtxVtxPair> edges;
    do {
      edges.clear();
      const VtxVec &vs = mesh->activeVertices();
      bool isCellMesh = (mesh->meshType() == "MGXC" or mesh->meshType() == "MGX2D");
      forall(const vertex &v, vs)
        forall(const vertex &n, S.neighbors(v)) {
          if(!S.uniqueLine(v, n))
            continue;
          // Only do edges on cell mesh
          if(isCellMesh and (v->label > 0 or n->label > 0))
            continue;
        if(maxDist == 0 or norm(v->pos - n->pos) > maxDist)
          edges.push_back(std::make_pair(v, n));
      }

      // Split the edges
      forall(const VtxVtxPair &pr, edges)
        splitEdge(S, pr.first, pr.second, sDiv);

    } while(progressAdvance() and maxDist > 0 and edges.size() > 0);

    // Update the GUI
    mesh->updateAll();
 
    return true;
  }
  REGISTER_PROCESS(SplitEdges);
  
  bool DeleteEdge::run(Mesh* mesh)
  {
    if(mesh->meshType() != "MGXM")
      throw(QString("Cannot delete edges in mesh type (%1), must be (MGXM)")
                                                                     .arg(mesh->meshType()));
    vvGraph& S = mesh->graph();
    const VtxVec &vs = mesh->activeVertices();
    if(vs.size() != 2)
      throw QString("Error, Must have exactly two vertices selected.");
  
    vertex v1 = vs[0], v2 = vs[1];
    if(!S.edge(v1, v2))
      throw QString("Error, Vertices have no edge between them.");
  
    S.eraseEdge(v1, v2);
    S.eraseEdge(v2, v1);
  
    mesh->updateAll();
    Information::out << "Edge deleted from mesh" << endl;
    return true;
  }
  REGISTER_PROCESS(DeleteEdge);
  
  bool MeshDeleteSelection::run(Mesh* m)
  {
    m->correctSelection(false);
    if(m->meshType() == "MGXM" or m->meshType() == "MGXC") {
      vvGraph& S = m->graph();
      size_t n = S.size();
      size_t count = 0;
      for(size_t i = 0; i < n; ++i) {
        size_t j = n - i - 1;
        if(S[j]->selected) {
          S.erase(S.begin() + j);
          count++;
        }
      }
      Information::out << count << " vertices deleted, " << S.size() 
                                << " vertices remaining in mesh" << endl;
    } else if(m->meshType() == "MGX2D") {
      CellTissue &T = m->tissue();
      CellSet D;
      forall(const vertex &v, T.S)
        if(v->selected) {
          cell c = T.getCell(v);
          if(c)
            D.insert(c);
        }
      forall(const cell &c, D)
        T.deleteCell(c);
      Information::out << D.size() << " cells deleted, " << T.C.size() 
                                  << " cells remaining in mesh" << endl;
    } else
      throw(QString("Cannot delete vertices in mesh type (%1), must be (MGXM, MGXC, MGX2D)")
                                                                     .arg(m->meshType()));
    m->updateAll();
  
    return true;
  }
  REGISTER_PROCESS(MeshDeleteSelection);
  
  bool MeshKeepVertices::run(Mesh* mesh, bool keep)
  {
    const std::vector<vertex>& vs = mesh->activeVertices();
  
    forall(const vertex& v, vs)
      if(keep)
        v->type = 'l';
      else if(v->type == 'l')
        v->type = 'j';
  
    mesh->updateAll();
    return true;
  }
  REGISTER_PROCESS(MeshKeepVertices);
}
