//
// 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.
//
#ifndef CELL_MESH_HPP
#define CELL_MESH_HPP

#include <CellTissue.hpp>
#include <Mesh.hpp>
#include <Parms.hpp>
#include <GraphUtils.hpp>
#include <DynamXProcessCellTissue.hpp>

#include <Triangulate.hpp>

#include "MGXViewer/qglviewer.h"
using namespace qglviewer;

namespace mgx 
{
  // Changing mesh types should require conversion routines below.
  const QString &CellTissue::setMeshType(const QString &mType)
  {
    _meshType = mType;
    return _meshType;
  }

  // Check graph
  bool CellTissue::chkGraph(const vvGraph &S) const
  {
    int badEdges = 0;
    #pragma omp parallel for
    for(uint i = 0; i < S.size(); i++) {
      vertex v = S[i];
      forall(const vertex &n, S.neighbors(v))
        if(!S.edge(n, v))
          #pragma omp atomic
          badEdges++;
    }
    if(badEdges > 0) {
      Information::out << badEdges << " bad edges found in graph" << endl;
      return true;
    }
    return false;
  }

  // Set area of cells, and length of cell walls
  void CellTissue::updGeometry(int mode) 
  {
    forall(const cell &c, C)
      updGeometry(c, mode);
    if(mode & UPD_NORMAL) {
      forall(const vertex &v, S)
        if(!v->cId)
          setNormal(S, v);
    }
  } 

  void CellTissue::updGeometry(const cell &c, int mode) 
  {
    if(!mode)
      return;

    if(_meshType == "MGX2D") {
      vertex vc = getVtx(c);
      Point3d cpos(0,0,0), cest(0,0,0);
  
      // First get average of points as estimate to center
      forall(const vertex &n, S.neighbors(vc))
        cest += n->pos;
      cest /= S.valence(vc);
  
      // Clear lengths
      if(mode & UPD_LENGTH) {
        forall(const cell &n, C.neighbors(c))
          C.edge(c,n)->length = 0;
      }

      double area = 0.0;
      forall(const vertex &n, S.neighbors(vc)) {
        // Find area
        vertex m = S.nextTo(vc,n);
        double ta = triangleArea(cest, m->pos, n->pos);
        area += ta;
        cpos += ta * (cest + n->pos + m->pos)/3;
    
        if(mode & UPD_LENGTH) {
          // Find length of cell wall
          vertex d = S.prevTo(n,vc,2);
          cell dc = getCell(d);
          if(dc) {
            C.edge(c,dc)->length += norm(m->pos - n->pos); 
            if(C.edge(c,dc)->length > 0) ; else { // Avoid nan
              Information::out << "Bad length " << C.edge(c,dc)->length << endl;
              C.edge(c,dc)->length = 0.1;
            }
          }
        }
      }
      if(area > 0) ; else { // Avoid nan
        Information::out << "Bad area " << area << endl;
        area = 0.1;
      }
      if(mode & UPD_AREA)
        c->area = area;
      if(mode & UPD_CENTER)
        c->pos = vc->pos = cpos / area;
      if(mode & UPD_NORMAL)
        setNormal(S, vc);
  
  
    } else { // if(_meshType == "MGX3D"){ // the "fake" cell division will not work if this is uncommented
  
      Point3d cellCentroid (0,0,0);
      double cellVolume = 0;
      double wallArea = 0;
  
      forall(const vertex &v, c->S) {
        forall(const vertex &n, c->S.neighbors(v)) {
          vertex m = c->S.nextTo(v,n);
          if(c->S.uniqueTri(v,n,m)){
            // Compute x component of center of mass
            float volume = signedTetraVolume(v->pos, n->pos, m->pos);
            cellCentroid += volume * (v->pos + n->pos + m->pos) / 4.0;
            cellVolume += volume;
            wallArea += triangleArea(v->pos, n->pos, m->pos);
          }
        }
      }
      //std::cout << "vol " << cellVolume << "/" << cellCentroid << std::endl;
  
      if(mode & UPD_AREA)
        c->area = wallArea;
      if(mode & UPD_CENTER)
        c->pos = cellCentroid / cellVolume;
      if(mode & UPD_VOLUME)
        c->volume = cellVolume;
    }
  }

  // Return whether Mesh has cell graph or not
  bool CellTissue::hasCellGraph()
  {
    if(_meshType == "MGX2D" or _meshType == "MGX3D")
      return true;
    else
      return false;
  }
 
  // Divide a cell
  bool CellTissue::divideCell(const cell &c, Subdivide *subdiv, DivAlg divAlg)
  {
    std::vector<Point3d> nrmls;
    return divideCell(c, subdiv, divAlg, nrmls);
  }

  // Divide a cell with a list of possible directions (only required for SHORTEST_DIR)
  bool CellTissue::divideCell(const cell &c, Subdivide *subdiv, DivAlg divAlg, std::vector<Point3d> &nrmls)
  {
    cell cl, cr;
    Point3d pu, pv;
    vertex u1, u2, v1, v2;
    double su = 0, sv = 0;

    // Find dividing wall
    if(!findCellDiv(divAlg, c, pu, pv, u1, u2, v1, v2, su, sv, nrmls))
      return false;

    // Tie the centers to the vertices
    vertex vc = getVtx(c), vcl, vcr;
    cl->vId = vcl.id();
    cr->vId = vcr.id();
    vcl->cId = cl.id();
    vcr->cId = cr.id(); 
    vcl->label = cl->label = nextLabel();
    vcr->label = cr->label = nextLabel();

    // Insert a new vertex on first wall
    vertex nu;
    nu->pos = pu;
    nu->label = -1;
    S.insert(nu); 
    S.insertEdge(nu, u1);
    S.insertEdge(nu, u2);
    S.spliceAfter(u1, u2, nu);
    S.spliceAfter(u2, u1, nu);
    if(subdiv)
      subdiv->updateEdgeData(u1, nu, u2, su);
    S.eraseEdge(u1, u2);
    S.eraseEdge(u2, u1);

    S.spliceBefore(nu,u1,vc);
    S.spliceAfter(vc,u1,nu);
  
    // Insert a new vertex on other side
    vertex nv;
    nv->pos = pv;
    nv->label = -1;
    S.insert(nv); 
    S.insertEdge(nv, v1);
    S.insertEdge(nv, v2);
    S.spliceAfter(v1, v2, nv);
    S.spliceAfter(v2, v1, nv);
    if(subdiv)
      subdiv->updateEdgeData(v1, nv, v2, sv);
    S.eraseEdge(v1, v2);
    S.eraseEdge(v2, v1);

    S.spliceAfter(nv,v2,vc);
    S.spliceAfter(vc,v1,nv);

    // If necessary adjust the connectivity in the adjacent cell.
    vertex uadj = S.prevTo(u1,nu);
    cell cuadj = getCell(uadj);
    if (cuadj and uadj == S.nextTo(u2,nu)) {
      S.spliceAfter(nu,u1,uadj);
      S.spliceAfter(uadj,u2,nu);
    }
  
    // Adjust connectivity in adjacent cell if required
    vertex vadj = S.nextTo(v2,nv);
    cell cvadj = getCell(vadj);
    if (cvadj and vadj == S.prevTo(v1,nv)) {
      S.spliceAfter(nv,v1,vadj);
      S.spliceAfter(vadj,v2,nv);
    }

    S.insert(vcl); S.insert(vcr); 
    C.insert(cl); C.insert(cr);

    // Connect new edge vertices
    S.insertEdge(vcl,nv); S.spliceAfter(vcl,nv, nu);
    S.insertEdge(vcr,nu); S.spliceAfter(vcr,nu, nv);

    S.spliceAfter(nu,vc,vcl);
    S.spliceBefore(nu,vc,vcr);
    S.spliceBefore(nv,vc,vcl);
    S.spliceAfter(nv,vc,vcr);

    S.replace(nu,vc,nv);
    S.replace(nv,vc,nu);

    // Connect new centers
    vertex x = S.nextTo(vc,nu);
    while (x != nv) {
      S.spliceBefore(vcr,nv,x);
      S.replace(x,vc,vcr);
      x = S.nextTo(vc,x);
    }
    vertex y = S.nextTo(vc,nv);
    while (y != nu) {
      S.spliceBefore(vcl,nu,y);
      S.replace(y,vc,vcl);
      y = S.nextTo(vc,y);
    }

    // Connect cells in cell graph, copy topology from S graph
    cell pc(0);
    forall(const vertex &n, S.neighbors(vcl)) {
      vertex d = S.prevTo(n,vcl,2);
      if(d->cId) {
        cell dc = getCell(d);
        if(dc != pc) {
          if(!pc)
            C.insertEdge(cl, dc);
          else
            C.spliceAfter(cl, pc, dc);
          C.spliceAfter(dc, c, cl);
          pc = dc;
        }
      }
    }
    pc = cell(0);
    forall(const vertex &n, S.neighbors(vcr)) {
      vertex d = S.prevTo(n,vcr,2);
      if(d->cId) {
        cell dc = getCell(d);
        if(dc != pc) {
          if(!pc)
            C.insertEdge(cr, dc);
          else
            C.spliceAfter(cr, pc, dc);
          C.spliceAfter(dc, c, cr);
          pc = dc;
        }
      }
    }

    // Update geometry of new cells
    updGeometry(cl); 
    updGeometry(cr);

    // Update geometry of neighbor cells
    CellSet U;
    forall(const cell &n, C.neighbors(cl))
      U.insert(n);
    forall(const cell &n, C.neighbors(cr))
      U.insert(n);
    forall(const cell &n, U)
      updGeometry(n);

    // Call subdivision object
    if(subdiv)
      subdiv->updateCellData(c, cl, cr);

// RSS top and bottom are not enough, maybe we should change the callback footprint?
//    if(!subdiv or subdiv->updateCellData(c, cl, cr, cuadj, cvadj)) {
//      if(cuadj)
//        updGeometry(cuadj); 
//      if(cvadj)
//        updGeometry(cvadj);
//    }

    // Propagate signal and heat to daughter cells
    vcl->signal = vc->signal;
    vcr->signal = vc->signal;

    // Remove the old cell from the graphs
    C.erase(c);
    S.erase(vc);

    return true;
  } 

  // Delete a cell from the graph
  void CellTissue::deleteCell(const cell &c)
  {
    vertex vc = getVtx(c);
    std::vector<vertex> J;
    forall(const vertex &n, S.neighbors(vc))
      J.push_back(n);
    S.erase(vc);
    C.erase(c);
    for(size_t i = 0; i < J.size(); ++i) {
      bool keep = false;
      forall(const vertex &n, S.neighbors(J[i]))
        if(n->cId) {
          keep = true;
          break;
        }
      if(!keep)
        S.erase(J[i]);
    }
  }

  // Read the parameters for the cell tissue
  void CellTissue::processParmsDivide(const QStringList &parms)
  {
    if(parms.size() < CellDivideMgx2d::pNumParms  )
      throw(QString("CellTissue::processParmsDivide: Error, not enough parameters."));

    CellMaxArea = parms[CellDivideMgx2d::pCellMaxArea].toFloat();
    CellDivAlg = stringToDivAlg(parms[CellDivideMgx2d::pCellDivAlg]);
    CellWallSample = parms[CellDivideMgx2d::pCellWallSample].toDouble();
    CellPinch = parms[CellDivideMgx2d::pCellPinch].toDouble();
    CellMaxPinch = parms[CellDivideMgx2d::pCellMaxPinch].toDouble();
    CellWallMin = parms[CellDivideMgx2d::pCellWallMin].toDouble();
  }

  // Read the parameters for the cell tissue
  void CellTissue::processParms(const QStringList &parms)
  {
    if(parms.size() < CellTissueMgx2d::pNumParms  )
      throw(QString("CellTissue::processParms: Error, not enough parameters."));
    
    CellUpdGeomMode = parms[CellTissueMgx2d::pCellUpdGeomMode].toInt();
    CellWallWidth = parms[CellTissueMgx2d::pCellWallWidth].toDouble();
    CellWallCorner = parms[CellTissueMgx2d::pCellWallCorner].toDouble();
    CellColorBegin = parms[CellTissueMgx2d::pCellColorBegin].toInt();
    CellColorEnd = parms[CellTissueMgx2d::pCellColorEnd].toInt();
    CellColorCenter = parms[CellTissueMgx2d::pCellColorCenter].toDouble();
    CellBorderColor = parms[CellTissueMgx2d::pCellBorderColor].toInt();

    CellPolWidth = parms[CellTissueMgx2d::pCellPolWidth].toDouble();
    CellPolCorner = parms[CellTissueMgx2d::pCellPolCorner].toDouble();
    CellPolColorBegin = parms[CellTissueMgx2d::pCellPolColorBegin].toInt();
    CellPolColorEnd = parms[CellTissueMgx2d::pCellPolColorEnd].toInt();
    
    CellNormalColor = parms[CellTissueMgx2d::pCellNormalColor].toInt();
    CellGraphColor = parms[CellTissueMgx2d::pCellGraphColor].toInt();
    CellLineColor = parms[CellTissueMgx2d::pCellLineColor].toInt();
    CellLineWidth = parms[CellTissueMgx2d::pCellLineWidth].toDouble();
  }

  // Add a cell to the mesh, first position is center
  bool CellTissue::addCell(std::vector<vertex> poly)
  {
    if(poly.size() < 3)
      return false;

    // Center
    vertex vc = poly[0];
    cell c;
    c->vId = vc.id();
    S.insert(vc); C.insert(c);
  
    vertex fv, pv; // First and previous vertex
    for(size_t i = 1; i < poly.size(); ++i) {
      vertex v = poly[i];
      S.insert(v);
  
      if(i == 1) {
        fv = pv = v;
        S.insertEdge(vc,v);
        S.insertEdge(v,vc);
      } else {
        S.spliceBefore(vc,pv,v);
        S.spliceAfter(pv,vc,v);
        S.insertEdge(v,pv); S.spliceAfter(v,pv, vc);
      }
      pv = v;
    }
    S.spliceAfter(pv,vc,fv);
    S.spliceBefore(fv,vc,pv);

    return true;
  }

  // Convert to MGXM
  bool CellTissue::toMgxm()
  {
    if(_meshType == "" or _meshType == "MGXC")
      _meshType = "MGXM"; // Nothing to do but change the type
    else if(_meshType == "MGX2D") {
      if(!toMgxc())
        return false;
      _meshType = "MGXM";
    } else if(_meshType == "MGX3D") {
      S.clear();
      // First insert vertices, we need to duplicate shared ones.
      // We'll use the saveId to mark them
      forall(const cell &c, C)
        forall(const vertex &v, c->S)
          v->saveId = 0;

      VtxVtxMap vMap;
      VtxCellMap cMap;
      typedef std::pair<cell, vertex> CellVtxPair;
      std::unordered_map<CellVtxPair, vertex> cvMap;
      forall(const cell &c, C)
        forall(const vertex &v, c->S)
          if(v->saveId == 0) {
            vMap[v] = v;
            cMap[v] = c;
            cvMap[CellVtxPair(c, v)] = v;
            S.insert(v);
            v->saveId++;
            v->label = c->label;
          } else {
            vertex t;
            // Copy the data
            // RSS Can we do anything with attributes here?
            *t = *v;
            vMap[t] = v;
            cMap[t] = c;
            cvMap[CellVtxPair(c, v)] = t;
            S.insert(t);
            t->label = c->label;
          }

      // Now insert neighborhoods
      forall(const vertex &v, S) {
        vertex u = vMap[v];
        cell c = cMap[v];
        vertex pn;
        forall(const vertex &m, c->S.neighbors(u)) {
          vertex n = cvMap[CellVtxPair(c, m)];
          if(S.valence(v) == 0)
            S.insertEdge(v, n);
          else
            S.spliceAfter(v, pn, n);
          pn = n;
        }
      }
      C.clear();
      _meshType = "MGXM";
    } else
      throw(QString("Cannot convert mesh type (%1) to MGXM").arg(_meshType));

		return true;
  }

  // Convert to MGXC
  bool CellTissue::toMgxc(double wallMax)
  {
    if(_meshType == "MGXM") {
      if(!mgxmToMgxc(S, S, wallMax))
				return false;
			_meshType = "MGXC";
	  } else if(_meshType == "MGX2D") {
      // Just delete the cell graph and change the type
      C.clear();
      _meshType = "MGXC";
		} else
      throw(QString("Cannot convert mesh type (%1) to MGXC").arg(_meshType));

    return true;
  }

  // Convert to MGX2D (from MGXC)
  bool CellTissue::toMgx2d(double wallMax)
  {
    if(_meshType == "MGXM") {
			if(!toMgxc(wallMax))
        throw(QString("Cannot convert mesh to cell mesh (MGXC)"));
		} else if(_meshType != "MGXC")
      throw(QString("Cannot convert mesh type (%1) to MGX2D").arg(_meshType));

    std::unordered_map<int, vertex> labelV;
    std::unordered_map<int, cell> labelC;

    // Create cell centers and tie to vertices
    C.clear();
    forall(const vertex &v, S) {
      v->cId = 0;
      if(v->label > 0) {
        labelV[v->label] = v;
        cell c;
        C.insert(c);
        c->label = v->label;
        labelC[c->label] = c;
        v->cId = c.id();
        c->vId = v.id();
      }
    }

    // Now connect C graph based on S
    forall(const cell &c, C) {
      cell pc(0);
      vertex v = getVtx(c);
      forall(const vertex &n, S.neighbors(v)) {
        vertex d = S.prevTo(n,v,2);
        if(d->cId) {
          cell dc = getCell(d);
          if(dc != pc) {
            if(!pc)
              C.insertEdge(c, dc);
            else
              C.spliceAfter(c, pc, dc);
            pc = dc;
          }
        }
      }
    }      

    // Change the mesh type
    _meshType = "MGX2D";

    // Update geometry
    updGeometry();

    // Clear modified on S graph
    S.setModified(false);

    return true;
  }

bool CellTissue::createTissueGraphs(vvGraph &cS,vvGraph &cJ,cellGraph & cC, double wallMax)
  {
    if(_meshType == "MGXM") {
			if(!toMgxc(wallMax))
        throw(QString("Cannot convert mesh to cell mesh (MGXC)"));
		} else if(_meshType != "MGXC")
      throw(QString("Cannot convert mesh type (%1) to MGX2D").arg(_meshType));

    cS = S;
    cJ = J;
    cC = C;

    std::unordered_map<int, vertex> labelV;
    std::unordered_map<int, cell> labelC;

    // Create cell centers and tie to vertices
    cC.clear();
    forall(const vertex &v, cS) {
      v->cId = 0;
      if(v->label > 0) {
        labelV[v->label] = v;
        cell c;
        cC.insert(c);
        c->label = v->label;
        labelC[c->label] = c;
        v->cId = c.id();
        c->vId = v.id();
      }
    }

    // Now connect C graph based on S
    forall(const cell &c, cC) {
      cell pc(0);
      vertex v = getVtx(c);
      forall(const vertex &n, cS.neighbors(v)) {
        vertex d = cS.prevTo(n,v,2);
        if(d->cId) {
          cell dc = getCell(d);
          if(dc != pc) {
            if(!pc)
              cC.insertEdge(c, dc);
            else
              cC.spliceAfter(c, pc, dc);
            pc = dc;
          }
        }
      }
    }      


   forall(const cell &c, cC) {  
      vertex vc = getVtx(c);
      Point3d cpos(0,0,0), cest(0,0,0);
  
      // First get average of points as estimate to center
      forall(const vertex &n, cS.neighbors(vc))
        cest += n->pos;
      cest /= cS.valence(vc);
  
      // Clear lengths
      //if(mode & UPD_LENGTH) {
        forall(const cell &n, cC.neighbors(c))
          cC.edge(c,n)->length = 0;
      //}

      double area = 0.0;
      forall(const vertex &n, cS.neighbors(vc)) {
        // Find area
        vertex m = cS.nextTo(vc,n);
        double ta = triangleArea(cest, m->pos, n->pos);
        area += ta;
        cpos += ta * (cest + n->pos + m->pos)/3;
    
//        if(mode & UPD_LENGTH) {
          // Find length of cell wall
          vertex d = cS.prevTo(n,vc,2);
          cell dc = getCell(d);
         if(dc) {
            cC.edge(c,dc)->length += norm(m->pos - n->pos); 
            if(cC.edge(c,dc)->length > 0) ; else { // Avoid nan
              Information::out << "Bad length " << C.edge(c,dc)->length << endl;
              cC.edge(c,dc)->length = 0.1;
            }
          }
        //}
      }

      if(area > 0) ; else { // Avoid nan
        Information::out << "Bad area " << area << endl;
        area = 0.1;
      }
//      if(mode & UPD_AREA)
        c->area = area;
//      if(mode & UPD_CENTER)
        c->pos = vc->pos = cpos / area;
//      if(mode & UPD_NORMAL)
        setNormal(cS, vc);
  
  
    }

    return true;
  }

 // create a map of label,label and their shared wall Area
 void CellTissue::getNeighborhoodMap3D(std::map<IntInt, double>& neighborMap)
 {

  typedef std::pair<vertex, vertex> VtxVtx;
  typedef std::pair<VtxVtx, vertex> Tri;

  std::map<IntInt, double> result;

  std::map<Tri, std::vector<int> > triLabelsMap;
  typedef std::pair<Tri, std::vector<int> > TriLabelPair;

  // go through all cells
  forall(const cell& c, C){
    // write triangles of cell to a map
    forall(const vertex &v, c->S){
      forall(const vertex &n, c->S.neighbors(v)){
        vertex m = c->S.nextTo(v,n);
        if(S.uniqueTri(v,n,m)){
          VtxVtx vv = std::make_pair(v,n);
          Tri tt = std::make_pair(vv,m);
          triLabelsMap[tt].push_back(c->label);
        }
        vertex m2 = c->S.prevTo(v,n); // needed as triangles are oriented the other way in neighboring cells
        if(S.uniqueTri(v,n,m2)){
          VtxVtx vv = std::make_pair(v,n);
          Tri tt = std::make_pair(vv,m2);
          triLabelsMap[tt].push_back(c->label);
        }
      }
    }
  }

  forall(const TriLabelPair &p, triLabelsMap){
    //if(p.second.size()>2) std::cout << "sth weird happened" << std::endl;
    if(p.second.size()==2){
      std::pair<int, int> neighborPair1 = std::make_pair(p.second[0], p.second[1]);
      std::pair<int, int> neighborPair2 = std::make_pair(p.second[1], p.second[0]);
      double triArea = triangleArea(p.first.first.first->pos, p.first.first.second->pos, p.first.second->pos)/2.; // divide by 2 as everything was counted twice
      result[neighborPair1] += triArea;
      result[neighborPair2] += triArea;
    }
  }

  neighborMap = result;

  }

  // Create a 3D cell mesh from a normal 2D one. Vertices with identical positions will
  // be merged, and a cell graph will be created.
  bool CellTissue::toMgx3d(double tolerance, double neighborMinArea)
  {
    if(_meshType != "MGXM")
      throw(QString("Cannot convert mesh type (%1) to cells").arg(_meshType));

    std::unordered_map<int, cell> labelC;
    std::unordered_map<vertex, vertex> vertexV;
    std::map<Point3d, vertex> posV;  // vertices with same pos are the same, must be sorted
    std::unordered_map<int, VtxSet> labelV;
    std::set<IntIntPair> wallP;

    // Create unique vertex list and cell list
    C.clear();
    progressStart("Creating cells", S.size());
    int count = 0;
    forall(const vertex &u, S) {
      if(u->label <= 0)
        throw(QString("Cannot convert to cells, unlabeled or border vertices in mesh"));

      labelV[u->label].insert(u);
      vertex v(0);

      if(posV.count(u->pos) > 0)
        v = posV[u->pos];
      else if(posV.empty() or tolerance <= 0) {
        posV[u->pos] = u;
        v = u;
      } else {
        Point3d p = u->pos;
        p.x() -= tolerance;
        std::map<Point3d, vertex>::iterator it = posV.lower_bound(p);
        if(it == posV.end())
          --it;
        bool found = false;
        while(it != posV.end() and it->first.x() - u->pos.x() < tolerance) {
          if(norm(u->pos - it->first) < tolerance) {
            v = it->second;
            found = true;
            break;
          }
          ++it;
        }
        if(!found) {
          posV[u->pos] = u;
          v = u;
        }
      }
      // store vertex mapping
      vertexV[u] = v;

      // store label pairs for walls (edges in cell graph)
      if(v->label != u->label) {
        wallP.insert(IntIntPair(v->label, u->label));
        wallP.insert(IntIntPair(u->label, v->label));
      }

      // Check if cell exists
      cell c;
      if(labelC.count(u->label) > 0)
        c = labelC[u->label];
      else {
        c->label = u->label;
        C.insert(c);
        labelC[c->label] = c;
      }
      // Insert vertex
      c->S.insert(v);

      progressAdvance(++count);
    }

    // Now create neighborhoods
    progressStart("Creating neighborhoods", C.size());
    count = 0;
    forall(const cell c, C) {
      forall(const vertex &u, labelV[c->label]) {
        vertex v = vertexV[u];
        v->label = 0;
        v->margin = false; // for now we assume they are closed
        v->minb = 0;
        v->cId = 0;
        vertex pn;
        forall(const vertex &m, S.neighbors(u)) {
          vertex n = vertexV[m];
          if(c->S.valence(v) == 0)
            c->S.insertEdge(v, n);
          else
            c->S.spliceAfter(v, pn, n);
          pn = n;
        }
      }
      progressAdvance(++count);
    }

    // old method
    //forall(const IntIntPair &p, wallP) 
    //  C.insertEdge(labelC[p.first], labelC[p.second]);

    // Change the mesh type
    _meshType = "MGX3D";

    // Update the all graph
    updateAllGraph();

    // Now create the cell graph, in 3D it can't normally be oriented
    std::map<IntInt, double> neighborMap;
    getNeighborhoodMap3D(neighborMap);
 
    forall(const IntIntDouble &p, neighborMap){
      //std::cout << "l " << p.first.first << "/" << p.first.second << "//" << p.second << std::endl;
      if(p.second > neighborMinArea)
        C.insertEdge(labelC[p.first.first], labelC[p.first.second]);
    }

    // Update geometry
    // RSS need to fix this, not everything implemented yet for 3D
    updGeometry();
 
    // Clear modified flag on S graph
    S.setModified(false);

    return true;
  }

  // Create/update the all graph from the 3D cell mesh (stored in S)
  void CellTissue::updateAllGraph()
  {
    if(_meshType != "MGX3D")
      return;

    // Clear S graph and replace with all graph
    S.clear();
    
    forall(const cell &c, C)
      forall(const vertex &v, c->S) {
        S.insert(v);
        // For a 3D mesh, the neighborhoods cannot be oriented properly 
        forall(const vertex &n, c->S.neighbors(v))
          if(!S.edge(v, n))
            S.insertEdge(v, n);
      }
  }

//  // Read cell mesh from a text file.
//  bool CellTissue::loadCellTissueFromText(const QString &fileName)
//  {
//    QFile file(fileName);
//    if (!file.open(QIODevice::ReadOnly | QIODevice::Text))
//      throw(QString("loadCellTissueFromText:Error, Cannot open input file: %1").arg(fileName));
//
//    S.clear(); C.clear();
//    IntVtxMap vmap;
//
//    QTextStream in(&file);
//    unsigned long vCnt;
//    in >> vCnt;
//
//    for(unsigned long i = 0; i < vCnt; i++) {
//      if(in.atEnd())
//        throw QString("Premature end of file reading vertices");
//      vertex v;
//      S.insert(v);
//      in >> v->saveId >> v->pos >> v->label;
//      vmap[v->saveId] = v;
//      // Bump up label if required
//      if(v->label > viewLabel())
//        setLabel(v->label);
//    }
//    for(size_t i = 0; i < vCnt; i++) {
//      if(in.atEnd()) 
//        throw QString("Premature end of file reading neighborhoods");
//      int vId, nCnt;
//      in >> vId >> nCnt;
//      typename 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(0);
//      for(int j = 0; j < nCnt; j++) {
//        int nId;
//        in >> nId;
//        typename IntVtxMap::const_iterator nit = vmap.find(nId);
//        if(nit == 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);
//        pn = n;
//      }
//    }
//    file.close();
//
//    // Convert to cell tissue
//    toMgx2d();
//
//    // Update the geometry
//    updGeometry();
//
//    Information::out << QString("Loaded mesh, file: %1, vertices: %2, cells: %3.")
//                                           .arg(fileName).arg(S.size()).arg(C.size()) << endl;
//    return true;
//  }

  // Save cell mesh to a text file.
  bool CellTissue::saveCellTissueToText(const QString &fileName)
  {
    QFile file(fileName);
    if (!file.open(QIODevice::WriteOnly)) {
      Information::out << "saveCellTissue::Error:Cannot open output file:" << fileName << endl;
      return false;
    }

    QTextStream out(&file);
    unsigned long saveId = 0;

    // Write number of vertices
    out << S.size() << endl;

    // Write vertices and update saveId
    forall(const vertex &v, S) {
      v->saveId = saveId++;
      out << v->saveId << " " << v->pos << " " << v->label << endl;
    }
    forall(const vertex &v, S) {
      out << v->saveId << " " << S.valence(v);
      forall(const vertex &n, S.neighbors(v))
        out << " " << n->saveId;
      out << endl;
    }
    file.close();
    return true;
  }

  // Apply Min distance from corner and return difference
  double CellTissue::findWallMin(Point3d &v, Point3d v1, Point3d v2, double mw)
  {
    Point3d v1v2 = v2 - v1;
    double v1v2l = norm(v1v2);
    v1v2 /= v1v2l;
  
    double len = (v - v1) * v1v2;
    double dis;
    if(v1v2l <= mw * 2.0)
      dis = (v1v2l/2.0) - len;
    else if(len < mw)
      dis = mw - len;
    else if(len > v1v2l - mw)
      dis = (v1v2l - mw) - len;
    else
      dis = 0.0;
    v = v1 + v1v2 * (len + dis);
    return(fabs(dis));
  }

  // Find Cell divide point
  bool CellTissue::findCellDiv(DivAlg divAlg, cell c, Point3d &pu, Point3d &pv, vertex &p_u1, 
    vertex &p_u2, vertex &p_v1, vertex &p_v2, double &su, double &sv, std::vector<Point3d> &nrmls)
  {
    // u1-u2, first edge, v1-v2 opposite edge
    Point3d minu, minv;
    vertex vc = getVtx(c);
    Point3d cpos = vc->pos;
  
    if(divAlg == CLOSEST_MID) {
      // Find closest midpoint on wall
      double mindis = HUGE_VAL;
      forall(const vertex &tu1, S.neighbors(vc)) { 
        // other end of wall
        vertex tu2 = S.nextTo(vc,tu1);
        Point3d u1 = tu1->pos, u2 = tu2->pos;
        Point3d u;
        u = (u1 + u2)/2.0;
        double dis = norm(u - cpos);
        if(dis < mindis) {
          mindis = dis;
          minu = u;
          p_u1 = tu1;
          p_u2 = tu2;
        }
      }
    } else if(divAlg == SHORTEST_WALL) {
      // Find shortest line through center
      double mindis = HUGE_VAL;
      Point3d n(0.0, 0.0, 0.0);
      forall(const vertex &n1, S.neighbors(vc)) {
        vertex n2 = S.nextTo(vc,n1);
        Point3d cn1 = n1->pos - cpos;
        Point3d cn2 = n2->pos - cpos;
        cn1 /= norm(cn1);
        cn2 /= norm(cn2);
        n += cn1 ^ cn2;
      }
      n /= S.valence(vc);
      forall(const vertex &tu1, S.neighbors(vc)) {
        vertex tu2 = S.nextTo(vc,tu1);
        Point3d u1 = tu1->pos, u2 = tu2->pos;
        Point3d u1u2 = u2 - u1;
        double u1u2l = norm(u1u2);
        u1u2 /= u1u2l;
        vertex tv1 = tu2;
        while(true) {
          vertex tv2 = S.nextTo(vc,tv1);
          if(tv2 == tu2)
            break;
          Point3d v1 = tv1->pos, v2 = tv2->pos;
            
          double ule = u1u2l - CellWallMin;
          for(double ul = CellWallMin; ul < ule; ul += CellWallSample) { 
            Point3d u = u1 + ul * u1u2;
            Point3d uc = cpos - u;
            Point3d nu = n ^ uc;
            Point3d v;
            double s;
            if(planeLineIntersect(u, nu, v1, v2, s, v)) {
              if(s >= 0  && s <= 1.0) {
                double dis = findWallMin(v, v1, v2, CellWallMin);
                dis += norm(u - v);
                if(dis < mindis) {
                  mindis = dis;
                  minu = u;
                  minv = v;
                  p_u1 = tu1;
                  p_u2 = tu2;
                  p_v1 = tv1;
                  p_v2 = tv2;
                }
              }
            }
          } // Samples
          tv1 = tv2;
        } // v1v2
      } // u1u2
      if(mindis == HUGE_VAL) {
        Information::out << "FindCellDiv::Error:Failed to find divide line" << endl;
        return(false);
      }
    } else if(divAlg == CLOSEST_WALL) {
      // Find closest point on wall
      double mindis = HUGE_VAL;
      forall(const vertex &tu1, S.neighbors(vc)) { 
        vertex tu2 = S.nextTo(vc,tu1);
    
        Point3d u1 = tu1->pos, u2 = tu2->pos;
        Point3d u1u2 = u2 - u1;
        u1u2 /= norm(u1u2);
        Point3d u = u1 + ((cpos - u1) * u1u2) * u1u2;
        findWallMin(u, u1, u2, CellWallMin);
        double dis = norm(u - cpos);
        if(dis < mindis) {
          mindis = dis;
          minu = u;
          p_u1 = tu1;
          p_u2 = tu2;
        }
      }
    } else if(divAlg == SHORTEST_DIR) {
      // Divide on a specified direction through the center
      double mindis = HUGE_VAL;
      vertex ttu1, ttu2;
      Point3d uu;
      for(size_t i = 0; i < nrmls.size(); ++i) {
        int count = 0;
        double dis = 0;
        forall(const vertex &tu1, S.neighbors(vc)) { 
          vertex tu2 = S.nextTo(vc,tu1);
          Point3d u, u1 = tu1->pos, u2 = tu2->pos;
          double s;
          if(planeLineIntersect(cpos, nrmls[i], u1, u2, s, u))
            if(s >= 0 and s <= 1.0) {
              ttu1 = tu1;
              ttu2 = tu2;
              uu = u;
              count++;
              dis += norm(u - cpos);
            }
        }
        if(count == 2 and dis < mindis) {
          mindis = dis;
          minu = uu;
          p_u1 = ttu1; 
          p_u2 = ttu2;
        }
      }
      if(mindis == HUGE_VAL) {
        Information::out << "FindCellDiv::Error:Failed to find divide line (SHORTEST_DIR)" << endl;
        return(false);
      }
    } else {
      Information::out << "FindCellDiv::Error:Bad Division Algorithm" << endl;
      return(false);
    }
    
    // Algorithms that only find pu in first step
    if(divAlg == CLOSEST_MID || divAlg == CLOSEST_WALL || divAlg == SHORTEST_DIR) {
      double mindis = HUGE_VAL;
      Point3d cnml = vc->nrml;
      vertex tv1 = p_u2;
      while(true) {
        vertex tv2 = S.nextTo(vc,tv1);
        Point3d v1 = tv1->pos, v2 = tv2->pos;
        if(tv2 == p_u2)
          break;
        // Keep uc close to same plane as v1v2
        Point3d uc = cpos - minu;
        uc -= (uc * cnml) * cnml;
        Point3d n = cnml ^ uc;
        uc += minu;
        double s;
        Point3d v;
        bool found = planeLineIntersect(uc, n, v1, v2, s, v);
        if(found) {
          double dis = findWallMin(v, v1, v2, CellWallMin);
          if(dis < mindis) {
            mindis = dis;
            minv = v;
            p_v1 = tv1;
            p_v2 = tv2;
          }
          if(dis == 0)
            break;
        }
        tv1 = tv2;
      }
    }

    // Find distance along walls values
    su = norm(minu - p_u1->pos)/norm(p_u2->pos - p_u1->pos);
    sv = norm(minv - p_v1->pos)/norm(p_v2->pos - p_v1->pos);
  
    // Calculate cell pinch
    double mind = 0.0;
    double pinch = 0.0;
  
    // If on edge, don't pinch
    vertex t1 = S.nextTo(p_u1,vc,2);
    vertex t2 = S.nextTo(p_u2,vc,2);
    if(t1->cId and t2->cId) {
      Point3d u1 = p_u1->pos;
      Point3d u2 = p_u2->pos;
      Point3d uc = (cpos - minu);
      double ucl = norm(uc);
      uc /= ucl;
      double uu1 = norm(minu - u1);
      double uu2 = norm(minu - u2);
      if(uu1 > uu2)
        mind = uu2;
      else
        mind = uu1;
      pinch = std::min(mind, ucl) * CellPinch;
      if(pinch > CellMaxPinch)
        pinch = CellMaxPinch;
      minu += uc * pinch;
    }
    t1 = S.nextTo(p_v1,vc,2);
    t2 = S.nextTo(p_v2,vc,2);
    if(t1->cId and t2->cId) {
      Point3d v1 = p_v1->pos;
      Point3d v2 = p_v2->pos;
      Point3d vc = (cpos - minv);
      double vcl = norm(vc);
      vc /= vcl;
      double vv1 = norm(minv - v1);
      double vv2 = norm(minv - v2);
      if(vv1 > vv2)
        mind = vv2;
      else
        mind = vv1;
      pinch = std::min(mind, vcl) * CellPinch;
      if(pinch > CellMaxPinch)
        pinch = CellMaxPinch;
      minv += vc * pinch;
    }
  
    // Return points
    pu = minu;
    pv = minv;

    return(true);
  }

  // Find Cell divide point
  bool CellTissue::findCellDiv3D(QString divAlg, cell c, Point3d &pu, Point3d &pv, vertex &p_u1, 
    vertex &p_u2, vertex &p_v1, vertex &p_v2, double &su, double &sv, std::vector<Point3d> &nrmls)
  {
    if(divAlg == "SHORTEST_DIR") {


    } else { // no other div alg implemented yet
      return false;
    }
    return true;
  }
}

#endif
