/*
 *
 *  This file is part of the Virtual Leaf.
 *
 *  The Virtual Leaf is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  The Virtual Leaf is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with the Virtual Leaf.  If not, see <http://www.gnu.org/licenses/>.
 *
 *  Copyright 2010 Roeland Merks.
 *
 */

#include <QDebug>

#include <string>
#include "cell.h"
#include "node.h"
#include "mesh.h"
#include "tiny.h"
#include "nodeset.h"
#include "cellitem.h"
#include "nodeitem.h"
#include "qcanvasarrow.h"
#include "parameter.h"


static const std::string _module_id("$Id$");

extern Parameter par;

double Cell::factor=1.;
double Cell::offset[3]={0,0,0};

Cell::Cell(void) : CellBase() {

	m=0;

}

Cell::Cell(double x, double y, double z) : CellBase(x,y,z) {

	m=0;
	
}

Cell::Cell(const Cell &src) :  CellBase(src) {
	
	m=src.m;
}

bool Cell::Cmp(Cell *c) const { return this->Index() < c->Index(); }
bool Cell::Eq(Cell *c) const { return this->Index() == c->Index(); }

Cell Cell::operator=(const Cell &src) {
	CellBase::operator=(src);
	m=src.m;
	return *this;
}
//Cell(void) : CellBase() {}

void Cell::DivideOverAxis(Vector axis) {
	
	// Build a wall
	// ->  find the position of the wall
	
	// better look for intersection with a simple line intersection algorithm as below?
	// this leads to some exceptions: e.g. dividing a horizontal rectangle.
	// leaving it like this for the time being
	
	if (dead) return;
	
	Vector centroid=Centroid();
	double prev_cross_z=(axis * (centroid - *(nodes.back()) ) ).z ;
		
	ItList new_node_locations;
	
	for (list<Node *>::iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		// cross product to detect position of division
		Vector cross = axis * (centroid - *(*i));
		
		if (cross.z * prev_cross_z < 0 ) {
			
			new_node_locations.push_back(i);
			
		}		
		prev_cross_z=cross.z;
	}
	
	DivideWalls(new_node_locations, centroid, centroid+axis);
	
}
double Cell::MeanArea(void) {
	return m->MeanArea();
}


void Cell::Apoptose(void) {
	
	// First kill walls
        #ifdef QDEBUG
        qDebug() << "This is cell " << Index() << endl;
	qDebug() << "Number of walls: " << walls.size() << endl;
	#endif
	for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
             #ifdef QDEBUG
	     qDebug() << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = "
		      << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
             #endif
	}
	for (list<Wall *>::iterator w=walls.begin();
		 w!=walls.end();
		 w++) {
		
		bool illegal_flag = false;
		if ((*w)->c1 == (*w)->c2 ) 
			illegal_flag=true;
		if ((*w)->c1 == this) {
			
			// invert wall?
			(*w)->c1 = (*w)->c2;      
			(*w)->c2 = m->boundary_polygon;
			
			Node *n1 = (*w)->n1;
			(*w)->n1 = (*w)->n2;
			(*w)->n2 = n1;
			
		} else {
			(*w)->c2 = m->boundary_polygon;
		}
		
		#ifdef QDEBUG
		if (illegal_flag && (*w)->c1==(*w)->c2) {
		  qDebug() << "I created an illegal wall." << endl;
		}
		#endif

		if ( ((*w)->N1()->DeadP() || (*w)->N2()->DeadP()) ||
			((*w)->C1() == (*w)->C2() ) ){
			// kill wall
		        #ifdef QDEBUG
		        qDebug() << "Killing wall." << endl;
			#endif
			(*w)->Kill();

			#ifdef QDEBUG
			if ((*w)) {
			  qDebug() << "Wall " << (*w)->Index() << " says: c1 = " 
				   << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
			}
			#endif
			(*w)=0;
		} else {
		        #ifdef QDEBUG
		        qDebug() << "Not killing wall." << endl;
			qDebug() << "Wall " << (*w)->Index() << " says: c1 = " 
				 << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;
                        #endif
		}
		
	}
	walls.remove(0);
	
	// Unregister me from my nodes, and delete the node if it no longer belongs to any cells
	list<Node *> superfluous_nodes;
	for (list<Node *>::iterator n=nodes.begin();
		 n!=nodes.end();
		 n++) {
		
		Node &no(*(*n));
		// locate myself in the node's owner list
		list<Neighbor>::iterator cellpos;
		bool cell_found=false;
		for (list<Neighbor>::iterator nb=no.owners.begin();
			 nb!=no.owners.end();
			 nb++) {
			if (nb->cell == this) {
				cellpos = nb;
				cell_found = true;
				break;
			}
		}
		
		if (!cell_found) {
			// I think this cannot happen; but if I am wrong, unpredictable things would happen. So throw an exception.
			throw ("Cell not found in CellBase::Apoptose()\n\rPlease correct the code to handle this situation.");
		}
		
		Neighbor tmp = *cellpos;
		no.owners.erase(cellpos);
		
		// if node has no owners left, or only has a connection to special cell -1 (outside world), mark it as dead.
		
		if (no.owners.size()==0 || (no.owners.size()==1 && no.owners.front().cell->BoundaryPolP()) ) {
			no.MarkDead();
		} else {
			// register node with outside world
			if (find_if( no.owners.begin(), no.owners.end(), 
				     bind2nd ( mem_fun_ref(&Neighbor::CellEquals), m->boundary_polygon->Index() ) ) == no.owners.end() ) {
				
				tmp.cell = m->boundary_polygon;
				no.owners.push_back(tmp);
			}
		}
	}
	
	
	
	// mark cell as dead
	MarkDead();
}

void Cell::ConstructConnections(void) {
	
    // Tie up the nodes of this cell, assuming they are correctly ordered
	
    //cerr << "Constructing connections of cell " << index << endl;
	
    for (list<Node *>::iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		//cerr << "Connecting node " << *i << endl;
		//cerr << "Node " << *i << endl << " = " << *(*i) << endl;
		// 1. Tidy up existing connections (which are part of this cell)
		if ((*i)->owners.size()>0) {
			list<Neighbor>::iterator neighb_with_this_cell=
			// remove myself from the neighbor list of the node
			find_if((*i)->owners.begin(),
					(*i)->owners.end(),
				bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );
			if (neighb_with_this_cell!=(*i)->owners.end()) 
				(*i)->owners.erase(neighb_with_this_cell);
		}
		
		Node *previous;
		if (i!=nodes.begin()) {
			list<Node *>::iterator previous_iterator=i;
			previous_iterator--;
			previous=*previous_iterator;
		} else {
			previous=nodes.back();
		}
		
		Node *next;
		list<Node *>::iterator next_iterator=i;
		next_iterator++;
		if (next_iterator==nodes.end()) {
			next=nodes.front();
		} else {
			next=*next_iterator;
		}
		
		//cerr << "[" << *i << "]";
		//if (*i==10 || *i==11) {
		//cerr << "previous = " << previous << ", next = " << next << endl;
		//}
		//if (*i!=10 && *i!=11)
		//cerr << "Node " << *i << endl << " = " << *(*i) << endl;
		(*i)->owners.push_back( Neighbor( this, previous, next ) );
		// if (*i==50 || *i==51) {
		//cerr << "Node " << *i << ".size() = " << (*i)->owners.size() << endl;
		// }
    }
}


/*! \brief Divide the cell over the line v1-v2.
 
 \param v1: First vertex of line.
 \param v2: Second vertex of line.
 \param fixed_wall: If true: wall will be set to "fixed" (i.e. not motile)
 \return: true if the cell divided, false if not (i.e. no intersection between v1 and v2, and the cell)
 */
bool Cell::DivideOverGivenLine(const Vector v1, const Vector v2, bool fix_cellwall, NodeSet *node_set ) {
	
	if (dead) return false;
	
	
	
	// check each edge for intersection with the line
	ItList new_node_locations;
	
	#ifdef QDEBUG
	qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl;
	#endif
	for (list<Node *>::iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		Vector v3 = *(*i);
		list<Node *>::iterator nb=i;
		nb++;
		if (nb == nodes.end()) {
			nb = nodes.begin();
		}
		Vector v4 = *(*nb);
		
		double denominator = 
		(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
		
		double ua = 
		((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
		double ub = 
		((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
		
		
		//cerr << "Edge " << *i << " to " << *nb << ": ua = " << ua << ", ub = " << ub << ":  ";
		// this construction with "TINY" should simulate open/closed interval <0,1]
		if ( ( TINY < ua && ua < 1.+TINY ) && ( TINY < ub && ub < 1.+TINY ) ) {
			// yes, intersection detected. Push the location to the list of iterators
			new_node_locations.push_back(nb);
			
		} 
	}
	
        #ifdef QDEBUG
	if (new_node_locations.size()<2) { 
	  qDebug() << "Line does not intersect with two edges of Cell " << Index() << endl;
	  qDebug() << "new_node_locations.size() = " << new_node_locations.size() << endl;
	  return false;
	}
	
	ItList::iterator i = new_node_locations.begin();
	list< Node *>::iterator j;
	qDebug() << "-------------------------------" << endl;
	qDebug() << "Location of new nodes: " << (**i)->Index() << " and ";

	++i;
	j = *i; 
	if (j==nodes.begin()) j=nodes.end(); j--;
	
	qDebug() << (*j)->Index() << endl;
	qDebug() << "-------------------------------" << endl;
    
	if ( **new_node_locations.begin() == *j ) {
	  qDebug() << "Rejecting proposed division (cutting off zero area)." << endl;
		return false;
	}
	#endif

	DivideWalls(new_node_locations, v1, v2, fix_cellwall, node_set);
	
	return true;
	
}

// Core division procedure
void Cell::DivideWalls(ItList new_node_locations, const Vector from, const Vector to, bool fix_cellwall, NodeSet *node_set) {
	
	if (dead) return;
	
	bool boundary_touched_flag=false;
	
	// Step 0: keep some data about the parent before dividing
	
	ParentInfo parent_info;
	parent_info.polarization = ReduceCellAndWalls<Vector>( PINdir );
	parent_info.polarization.Normalise();
	parent_info.PINmembrane = SumTransporters(1);
	parent_info.PINendosome = Chemical(1);
	
	//cerr << "Parent polarization before division: " << parent_info.polarization << endl;
	
	// Step 1: create a daughter cell
	Cell *daughter=m->AddCell(new Cell());
    
	// Step 2: Copy the basics of parent cell to daughter
	for (int i=0;i<NChem();i++) {
		daughter->chem[i]=chem[i];
	}
	
	daughter->cell_type = cell_type;
	
	for (int i=0;i<NChem();i++) {
		daughter->new_chem[i]=new_chem[i];
	}
	
	
	daughter->boundary=boundary;
	daughter->m=m;
	
	daughter->target_area=target_area/2.;
	
	target_area/=2;
	daughter->cellvec=cellvec;
	
	
	// Division currently only works for convex cells: i.e. if the division line
	// intersects the cells at two points only.
	if (new_node_locations.size()!=2) {
		
		// Note: if you would add the possibility of dividing non-convex
		// cells, remember to update the code below. There are some
		// fixed-size arrays over there!
		
		cerr << "Warning in Cell::Division: division of non-convex cells not fully implemented" << endl;
		
		// Reject the daughter cell and decrement the amount of cells
		// again. We can do this here because it is the last cell added.
		// Never, ever try to fully erase a cell elsewhere, because we
		// make heavy use of cell indices in this project; if you erase a
		// Cell somewhere in the middle of Mesh::Cells the indices will
		// get totally messed up...! (e.g. the indices used in Nodes::cells)
		
		#ifdef QDEBUG
		qDebug() << "new_node_locations.size() = " << new_node_locations.size() <<endl;
		qDebug() << "daughter->index = " << daughter->index << endl;
		qDebug() << "cells.size() = " << m->cells.size() << endl;
		#endif

		m->cells.pop_back();
		Cell::NCells()--;
		m->shuffled_cells.pop_back();
		return;
	}
	
	
	// We can be sure we only need two positions here because divisions
	// of non-convex cells are rejected above.
	Vector new_node[2];
	Node *new_node_ind[2];
	
	int new_node_flag[2];
	Edge div_edges[2];
	
	int nnc=0;
	
	Wall *div_wall[4];
	double orig_length[2];
	for (int i=0;i<4;i++) { div_wall[i]=0; orig_length[i/2] = 0.; }
	
	// construct new Nodes at the intersection points
	// unless they coincide with existing points
	for ( ItList::const_iterator i=new_node_locations.begin();
		 i!=new_node_locations.end();
		 i++) {
		
		// intersection between division axis
		// and line from this node to its predecessor
		
		// method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
		Vector v1 = from;
		Vector v2 = to;
		Vector v3 = *(**i);
		
		// get previous node
		list<Node *>::iterator nb=*i;
		if (nb == nodes.begin()) {
			nb = nodes.end();
		} 
		nb--;
		Vector v4=*( *nb ); 
		
		double denominator = 
		(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
		
		double ua = 
		((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
		
		double intersec_x = v1.x + ua*(v2.x-v1.x);
		double intersec_y = v1.y + ua*(v2.y-v1.y);
		
		// construct a new node at intersec
		// we construct a vector temporarily,
		// until we are sure we are going to keep the node...
		// Node count is damaged if we construct superfluous nodes
		Vector *n=new Vector(intersec_x,intersec_y,0);
		
		div_edges[nnc].first=*nb;
		div_edges[nnc].second=**i;
		
		// Insert this new Node if it is far enough (5% of element length)
		// from one of the two existing nodes, else use existing node
		//
		// old, fixed value was: par.collapse_node_threshold = 0.05
		double collapse_node_threshold = 0.05;
#ifdef FLEMING
		collapse_node_threshold = par.collapse_node_threshold;
#endif
		
		double elem_length = ( (*(**i)) - (*(*nb)) ).Norm();
		if ( ( *(**i) - *n ).Norm() < collapse_node_threshold  * elem_length ) {
			new_node_flag[nnc]=1;
			new_node[nnc] = *(**i);
			new_node_ind[nnc] = **i;
			//cerr << **i << endl ;
		} else 
			if ( (*(*nb) - *n).Norm() < collapse_node_threshold * elem_length ) {
				new_node_flag[nnc]=2;
				new_node[nnc] = *(*nb);
				new_node_ind[nnc] = *nb;
			} else {
				new_node_flag[nnc]=0;
				new_node[nnc] = *n;
			}
		
		nnc++;
		delete n;
	}
	
	
	for (int i=0;i<2;i++) {
		
		Cell *neighbor_cell=0; // we need this to split up the "Wall" objects.
		
		// for both divided edges: 
		//      insert its new node into all cells that own the divided edge
		// but only if it really is a new node:
		if (new_node_flag[i]!=0) {
			if (fix_cellwall) {
				(new_node_ind[i])->fixed = true;
				
				// all this we'll do later for the node set :-)
				/* (new_node_ind[i])->boundary = true;
				 (new_node_ind[i])->sam = true;
				 boundary = SAM;
				 daughter->boundary = SAM;
				 boundary_touched_flag = true;
				 */ 
			}
			
		} else {
			
			// (Construct a list of all owners:)
			// really construct the new node (if this is a new node)
			new_node_ind[i] = 
			m->AddNode(new Node (new_node[i]) );
			
			
			
			// if a new node is inserted into a fixed edge (i.e. in the petiole)
			// make the new node fixed as well
			(new_node_ind[i])->fixed = (div_edges[i].first)->fixed &&
			(div_edges[i].second)->fixed;
			
			// Insert Node into NodeSet if the div_edge is part of it.
			if (
				(div_edges[i].first->node_set && div_edges[i].second->node_set) &&
				(div_edges[i].first->node_set == div_edges[i].second->node_set))
			{
				//cerr << "Inserting node into node set\n";
				div_edges[i].first->node_set->AddNode( new_node_ind[i] );
			}
			
			// if the new wall should be fixed (i.e. immobile, or moving as
			// solid body), make it so, and make it part of the boundary. Using
			// this to make a nice initial condition by cutting off part of a
			// growing leaf.
			
			if (fix_cellwall) {
				(new_node_ind[i])->fixed = true;
				
				// All this we'll do later for the node set only
			}
			
			// if new node is inserted into the boundary
			// it will be part of the boundary, too

			new_node_ind[i]->UnsetBoundary();
			if ((div_edges[i].first->BoundaryP() && div_edges[i].second->BoundaryP()) && // Both edge nodes are boundary nodes AND
			     ((m->findNextBoundaryNode(div_edges[i].first))->Index() == div_edges[i].second->Index())){ // The boundary proceeds from first to second.

                                #ifdef QDEBUG
			        qDebug() << "Index of the first node: " << div_edges[i].first->Index() << endl;
			        qDebug() << "Index of the second node: " << div_edges[i].second->Index() << endl;
			        qDebug() << "Boundary proceeds from: " <<  div_edges[i].first->Index() 
				         << "to: " << (m->findNextBoundaryNode(div_edges[i].first))->Index() << endl << endl;
                                #endif
			        new_node_ind[i]->SetBoundary();

				// We will need to repair the boundary polygon later, since we will insert new nodes
				//cerr << "Boundary touched for Node " << new_node_ind[i]->Index() << "\n";
				boundary_touched_flag=true;
				
				// and insert it into the boundary_polygon
				// find the position of the first node in the boundary
				list<Node *>::iterator ins_pos = find
				(m->boundary_polygon->nodes.begin(),
				 m->boundary_polygon->nodes.end(),
				 div_edges[i].first);
				// ... second node comes before or after it ...
				if (*(++ins_pos!=m->boundary_polygon->nodes.end()?
					  ins_pos:m->boundary_polygon->nodes.begin())!=div_edges[i].second) {
					
					m->boundary_polygon->nodes.insert(((ins_pos--)!=m->boundary_polygon->nodes.begin()?ins_pos:(--m->boundary_polygon->nodes.end())), new_node_ind[i]);
					
					// .. set the neighbors of the new node ...
					// in this case e.second and e.first are inverted
				} else {
					// insert before second node, so leave ins_pos as it is,
					// that is: incremented
					m->boundary_polygon->nodes.insert(ins_pos, new_node_ind[i]);	
					// .. set the neighbors of the new node ...
				}
			}
			
			list<Neighbor> owners;
			
			// push all cells owning the two nodes of the divides edges
			// onto a list

			copy((div_edges[i].first)->owners.begin(),
				 (div_edges[i].first)->owners.end(),
				 back_inserter(owners));
			copy((div_edges[i].second)->owners.begin(),
				 (div_edges[i].second)->owners.end(),
				 back_inserter(owners));
			
			
			// find first non-self duplicate in the owners: 
			// cells owning the same two nodes
			// share an edge with me
			owners.sort( mem_fun_ref( &Neighbor::Cmp ) );


                        #ifdef QDEBUG  
			list<Neighbor> unique_owners;
			copy(owners.begin(), owners.end(), back_inserter(unique_owners));
			unique_owners.unique( mem_fun_ref( &Neighbor::Eq ) );
			qDebug() << "The dividing edge nodes: " << div_edges[i].first->Index() << " and " << div_edges[i].second->Index() << " are owned by cells: ";
			// spit out each owners' cell index
			foreach(Neighbor neighbor, unique_owners){
			  qDebug() << neighbor.cell->Index() << "  ";
			}
			qDebug() << endl;
                        #endif

			// Search through the sorted list of edge node owners looking for duplicate pairs. Each pair represents an actual edge owner.
			list<Neighbor> edge_owners;
			list<Neighbor>::iterator it;
			for (it=owners.begin(); it!=owners.end(); it++) {
			  it = adjacent_find(it, owners.end(), neighbor_cell_eq);
			  if (it == owners.end()) break; // bail if reach the end of the list
                          #ifdef QDEBUG
			  qDebug() << "Considering: " << it->cell->Index() << " as a possible edge owner." << endl;
                          #endif
			  if (it->cell->Index() != this->Index()) {
                            #ifdef QDEBUG
			    qDebug() << "Adding: " << it->cell->Index() << " to the list of edge owners." << endl;
                            #endif
			    edge_owners.push_back(*it);
			  }
			} 

			if (edge_owners.size() > 1){
			  // Remove the boundary polygon - if its there
			  list<Neighbor>::iterator it;
			  if ((it = find_if (edge_owners.begin(), edge_owners.end(), bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1))) != edge_owners.end()) {
                            #ifdef QDEBUG
			    qDebug() << "deleating: " << it->cell->Index() << " from the list of edge owners." << endl;
                            #endif
			    edge_owners.erase(it);
			  }
			}

                        #ifdef QDEBUG
			qDebug() << "The edge owners list has: " << edge_owners.size() << " elements" << endl;
			#endif

			// Since the list should always contain exactly one element, pass it on as an iterator
			list<Neighbor>::iterator c = (edge_owners.size() != 0) ? edge_owners.begin() : edge_owners.end();

			// (can we have more than one neighboring cell here??)
			if (c!=owners.end()) { 
				neighbor_cell = c->cell;
				if (!c->cell->BoundaryPolP()) {

					// find correct position in the cells node list
					// to insert the new node
					list<Node *>::iterator ins_pos = find
					(neighbor_cell->nodes.begin(),
					 neighbor_cell->nodes.end(),
					 div_edges[i].first);
					
					neighbor_cell->nodes.insert(ins_pos, new_node_ind[i]);
					neighbor_cell->ConstructConnections();
					
					// give walls to daughter later
			  }
			} else {
				neighbor_cell = 0;
			}
		}
		
		// Split the Wall with the neighboring cell
		
		// if the neighbor cell has not yet been identified above, do it now
		if (neighbor_cell == 0) {
			
			list<Neighbor> owners;
			
			// push all cells owning the two nodes of the divides edges
			// onto a list
			copy((div_edges[i].first)->owners.begin(),
				 (div_edges[i].first)->owners.end(),
				 back_inserter(owners));
			copy((div_edges[i].second)->owners.begin(),
				 (div_edges[i].second)->owners.end(),
				 back_inserter(owners));
			
			
			// find first non-self duplicate in the owners: 
			// cells owning the same two nodes
			// share an edge with me
			owners.sort( mem_fun_ref( &Neighbor::Cmp ) );

			list<Neighbor>::iterator c;
			for (c=owners.begin();
				 c!=owners.end();
				 c++) {
				c=adjacent_find(c,owners.end(),neighbor_cell_eq);
				if (c->cell->Index() != this->Index() || c==owners.end()) break;
			}
			
			if (c!=owners.end())
				neighbor_cell = c->cell;
			else 
				neighbor_cell = 0;
		}
		
		
		if (neighbor_cell /* && !neighbor_cell->BoundaryPolP() */) {
			
			//cerr << "Cell "  << index << " says: neighboring cell is " << neighbor_cell->index << endl;
			
			/*************** 1. Find the correct wall element  ********************/
			
			list<Wall *>::iterator w, start_search;
			w = start_search = walls.begin();
			do {
				// Find wall between this cell and neighbor cell
				w = find_if( start_search, walls.end(), bind2nd (mem_fun( &Wall::is_wall_of_cell_p ), neighbor_cell ) );
				start_search = w; start_search++; // continue searching at next element
			} while ( w!=walls.end() && !(*w)->IntersectsWithDivisionPlaneP( from, to ) ); // go on until we find the right one.
			
			if (w == walls.end()) {
			        #ifdef QDEBUG
			        qDebug() << "Whoops, wall element not found...!" << endl;
			        qDebug() << "Cell ID: " << neighbor_cell->Index() << endl;
				qDebug() << "My cell ID: " << Index() << endl;
                                #endif
			} else {
				
				// 2. Split it up, if we should (sometimes, the new node coincides with an existing node so
				// we should not split up the Wall)
				
				if (new_node_ind[i]!=(*w)->n1 && new_node_ind[i]!=(*w)->n2) {
					
					Wall *new_wall;
					
					// keep the length of the original wall; we need it to equally divide the transporter concentrations
					// over the two daughter walls
					(*w)->SetLength(); // make sure we've got the current length
					orig_length[i] = (*w)->Length();
					//cerr << "Original length is " << orig_length[i] << endl;
					if ((*w)->c1 == this ) {
						
						//  cerr << "Cell " << (*w)->c1->Index() << " splits up wall " << *(*w) << ", into: " << endl;
						new_wall = new Wall( (*w)->n1, new_node_ind[i], this, neighbor_cell);
						(*w)->n1 = new_node_ind[i];
						
						//  cerr << "wall " << *(*w) << ", and new wall " << *new_wall << endl;
						
					} else {
						new_wall = new Wall( (*w)->n1, new_node_ind[i], neighbor_cell, this);
						
						(*w)->n1 = new_node_ind[i];
					}
					
					
					// 3. Give wall elements to appropriate cells
					if (new_wall->n1 != new_wall->n2) {
						
						if (par.copy_wall)
							new_wall->CopyWallContents(**w);
						else {
							// If wall contents are not copied, decide randomly which wall will be the "parent"
							// otherwise we will get biases (to the left), for example in the meristem growth model
							if (RANDOM()<0.5) {
								new_wall->SwapWallContents(*w);
							}
						}
						AddWall(new_wall);
						// cerr << "Building new wall: this=" << Index() << ", neighbor_cell = " << neighbor_cell->Index() << endl;
						
						neighbor_cell->AddWall( new_wall);
						//cerr << "Existing wall: c1 = " << (*w)->c1->Index() << ", neighbor_cell = " << (*w)->c2->Index() << endl;
						
						// Remember the addresses of the new walls
						div_wall[2*i+0] = *w;
						div_wall[2*i+1] = new_wall;
						
						// we will correct the transporter concentrations later in this member function, after division
						// First the new nodes should be inserted into the cells, before we can calculate wall lengths
						// Remember that cell walls can be bent, so have a bigger length than the Euclidean distance n1->n2
						
					} else {
						delete new_wall;
					}
				}
			}
		}
	}  // closing loop over the two divided edges (for (int i=0;i<2;i++) )
	
	// move half of the nodes to the daughter
	{
		//cerr << "Daughter: ";
		list<Node *>::iterator start, stop;
		
		start=new_node_locations.front();
		
		//cerr << "*new_node_locations.front() = " << *new_node_locations.front() << endl;
		if (new_node_flag[0]==1) {
			start++;
			if (start==nodes.end())
				start=nodes.begin();
		}  
		
		stop=new_node_locations.back();
		if (new_node_flag[1]==2) {
			if (stop==nodes.begin())
				stop=nodes.end();
			stop--;
		}
		list<Node *>::iterator i=start;
		while ( i!=stop) {
			
			// give the node to the daughter
			// (find references to parent cell from this node,
			// and remove them)
			list<Neighbor>::iterator neighb_with_this_cell=
			find_if((*i)->owners.begin(),
					(*i)->owners.end(),
				bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );
			if (neighb_with_this_cell==(*i)->owners.end()) {
			        #ifdef QDEBUG
			  qDebug() << "not found" << endl;
                                #endif
				abort();
			}
			
			(*i)->owners.erase(neighb_with_this_cell);
			
			daughter->nodes.push_back( *i );
			
			
			i++;
			if (i==nodes.end())
				i=nodes.begin();
		};
	}
	
	// new node list of parent
	list<Node *> new_nodes_parent;
	
	// half of the nodes stay with the parent
	{
		list<Node *>::iterator start, stop;
		start=new_node_locations.back();
		if (new_node_flag[1]==1) {
			start++;
			if (start==nodes.end())
				start=nodes.begin();
		}
		stop=new_node_locations.front();
		if (new_node_flag[0]==2) {
			if (stop==nodes.begin())
				stop=nodes.end();
			stop--;
		}
		
		list<Node *>::iterator i=start;
		while (i!=stop) {
			new_nodes_parent.push_back( *i );
			
			i++;
			if (i==nodes.end()) 
				i = nodes.begin();
		};
	}

	// insert shared wall
	// insert shared nodes on surface of parent cell
	new_nodes_parent.push_back( new_node_ind[0] );
	daughter->nodes.push_back ( new_node_ind[1] );
	
	// optionally add the new node to the nodeset (passed by pointer)
	// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
	if (node_set) {
		node_set->AddNode( new_node_ind[0] );
	}
	
#define MULTIPLE_NODES
#ifdef MULTIPLE_NODES
	// intermediate, extra nodes
	// Calculate distance between the two new nodes
	double dist=( new_node[1] - new_node[0] ).Norm();
	//bool fixed_wall = (new_node_ind[0])->fixed && (new_node_ind[1])->fixed;
	bool fixed_wall = false;
	
	// Estimate number of extra nodes in wall
	// factor 4 is to keep tension on the walls;
	// this is a hidden parameter and should be made explicit
	// later on.
	int n=(int)((dist/Node::target_length)/4+0.5);
	
	Vector nodevec = ( new_node[1]- new_node[0]).Normalised();
	
	double element_length = dist/(double)(n+1);
	
	// note that wall nodes need to run in inverse order in parent
	list<Node *>::iterator ins_pos = daughter->nodes.end();
	for (int i=1;i<=n;i++) {
		Node *node=
		m->AddNode( new Node( new_node[0] + i*element_length*nodevec ) );
		
		node->fixed=fixed_wall;
		
		if (!fix_cellwall)
			node->boundary = false;
		else {
			node->fixed = true;
		}
		
		ins_pos=daughter->nodes.insert(ins_pos, node );
		new_nodes_parent.push_back( node );
		
		// optionally add the new node to the nodeset (passed by pointer)
		// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
		if (node_set) {
			node_set->AddNode( node );
		}
		
	}
#endif
	daughter->nodes.push_back( new_node_ind[0] );
	new_nodes_parent.push_back( new_node_ind[1] );
	
	// optionally add the new node to the nodeset (passed by pointer)
	// (in this way we can move the NodeSet as a whole; useful for a fixed cutting line)
	if (node_set) {
		node_set->AddNode( new_node_ind[1] );
	}
	
	// move the new nodes to the parent
	nodes.clear();
	copy( new_nodes_parent.begin(), 
		 new_nodes_parent.end(), 
		 back_inserter(nodes) );
	
	
	// Repair cell lists of Nodes, and node connectivities
	ConstructConnections();
	daughter->ConstructConnections();
	
	if (boundary_touched_flag) {
		m->boundary_polygon->ConstructConnections();
	} 
	
	// collecting neighbors of divided cell
	list<CellBase *> broken_neighbors;
	
	// this cell's old neighbors
	copy(neighbors.begin(), neighbors.end(), back_inserter(broken_neighbors) );
	
	// this cell
	broken_neighbors.push_back(this);
	
	// its daughter
	broken_neighbors.push_back(daughter);
	
	
	
	// Recalculate area of parent and daughter
	area = CalcArea();
	daughter->area = daughter->CalcArea();
	
	SetIntegrals();
	daughter->SetIntegrals();
	
    
	
	// Add a "Cell Wall" for diffusion algorithms
	Wall *wall = new Wall( new_node_ind[0], new_node_ind[1], this, daughter );
	
	AddWall( wall );
	
	daughter->AddWall( wall );
	
	//cerr << "Correct walls of cell " << Index() << " and daughter " << daughter->Index() << endl;
	
	// Move Walls to daughter cell
	list <Wall *> copy_walls = walls;
	for (list<Wall *>::iterator w = copy_walls.begin();
		 w!=copy_walls.end();
		 w++) {
		
		//cerr << "Doing wall, before:  " << **w << endl;
		
		//  checks the nodes of the wall and gives it away if appropriate
		(*w)->CorrectWall ( );
		
		//cerr << "and after: " << **w << endl;
		
	}
	
	
	
	// Correct tranporterconcentrations of divided walls
	for (int i=0;i<4;i++) {
		if (div_wall[i]) {
			div_wall[i]->SetLength();
			div_wall[i]->CorrectTransporters(orig_length[i/2]);
		}
	}
	
	//cerr << "Cell " << index << " has been dividing, and gave birth to Cell " << daughter->index << endl;
	
	// now reconstruct neighbor list for all "broken" neighbors
	
	for (list<CellBase *>::iterator i=broken_neighbors.begin();
		 i!=broken_neighbors.end();i++) {
		((Cell *)(*i))->ConstructNeighborList();
	}
	
	
	ConstructNeighborList();
	daughter->ConstructNeighborList();
	
	m->plugin->OnDivide(&parent_info, daughter, this);
	
	daughter->div_counter=(++div_counter);
	
	
}

// Move the whole cell
void Cell::Move(const Vector T) {
	
    for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		*(*i)+=T;
    }
}

double Cell::Displace(double dx, double dy, double dh) {
	
	// Displace whole cell, add resulting energy to dh,
	// and accept displacement if energetically favorable
	// 
	// Method is called if a "fixed" node is displaced
	
	// Warning: length constraint not yet  CORRECTLY implemented for this function
	
	// Attempt to move this cell in a random direction
	//  Vector movement(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
	
	
	dh=0;
	
	Vector movement(dx,dy,0);
	
	vector< pair<Node *, Node *> > length_edges;
	vector<double> cellareas;
	cellareas.reserve(neighbors.size());
	
	// for the length constraint, collect all edges to this cell's nodes,
	// which are not part of the cell
	// the length of these edges will change
	
	double old_length=0.;
	for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		

		for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
			 n!=(*i)->owners.end();
			 n++) {
			
			if (n->getCell()!=this) {
				length_edges.push_back( pair <Node *,Node *> (*i, n->nb1) );
				length_edges.push_back( pair <Node *,Node *> (*i, n->nb2) );
				old_length += 
				DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
				DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
			}
		}
	}
	
	// calculate area energy difference of neighboring cells
	// (this cells' shape remains unchanged)
	double old_area_energy=0., old_length_energy=0.;
	for (list<CellBase *>::const_iterator i=neighbors.begin();
		 i!=neighbors.end();
		 i++) {
		old_area_energy += DSQR((*i)->Area()-(*i)->TargetArea());
		old_length_energy += DSQR((*i)->Length()-(*i)->TargetLength());
	}
	
	Move(movement);
	
	double new_area_energy=0., new_length_energy=0.;
	for (list<CellBase *>::const_iterator i=neighbors.begin();
		 i!=neighbors.end();
		 i++) {
		cellareas.push_back((*i)->CalcArea());
		new_area_energy += DSQR(cellareas.back()-(*i)->TargetArea());
		new_length_energy += DSQR((*i)->CalcLength()-(*i)->TargetLength());
	}
	
	double new_length=0;
	for ( vector< pair< Node *, Node * > >::const_iterator e = length_edges.begin();
		 e != length_edges.end();
		 e++) {
		new_length +=  DSQR(Node::target_length-
							(*(e->first)-*(e->second)).Norm());
	}
	
	
	dh += (new_area_energy - old_area_energy) + (new_length_energy - old_length_energy) * lambda_celllength +
	par.lambda_length * (new_length - old_length);
	
	if (dh<0 || RANDOM()<exp(-dh/par.T)) {
		
		// update areas of cells
		//cerr << "neighbors: ";
		list<CellBase *>::const_iterator nb_it = neighbors.begin();
		for (vector<double>::const_iterator ar_it = cellareas.begin();
			 ar_it!=cellareas.end();
			 ( ar_it++, nb_it++) ) {
			((Cell *)(*nb_it))->area = *ar_it;
			(*nb_it)->SetIntegrals(); 
		}
		
		//cerr << endl;
		
	} else {
		
		Move ( -1*movement);
		
	}
	
	return dh;
}


void Cell::Displace (void) {
	Displace(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);
}

// Get energy level of whole cell (excluding length constraint?)
double Cell::Energy(void) const {
	
	double energy = 0.;
	double length_contribution = 0.;
	
	for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		for (list<Neighbor>::const_iterator n=(*i)->owners.begin();
			 n!=(*i)->owners.end();
			 n++) {
			
			if (n->getCell()==this) {
				
				length_contribution += 
				DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+
				DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());
				
			}
		}
	}
	
	// wall elasticity constraint
	energy += par.lambda_length * length_contribution;
	
	// area constraint
	energy += DSQR(CalcArea() - target_area);
	
	// cell length constraint
	
	
	energy += lambda_celllength * DSQR(Length() - target_length);
	
	
	return energy;
}





bool Cell::SelfIntersect(void) {
	
    // The (obvious) O(N*N) algorithm
	
    // Compare each edge against each other edge
	
    // An O(N log(N)) algorithm by Shamos & Hoey (1976) supposedly exists;
    // it was mentioned on comp.graphics.algorithms
	
    // But I haven't been able to lay my hand on the paper.
    // Let's try whether we need it....
	
    // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
	
    for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		list<Node *>::const_iterator j=i; 
		++j;
		for (;
			 j!=nodes.end();
			 j++) 
		{
			
			Vector v1 = *(*i);
			list<Node *>::const_iterator nb=i;
			nb++;
			if (nb == nodes.end()) {
				nb = nodes.begin();
			} 
			Vector v2 = *(*nb);
			Vector v3 = *(*j);
			nb=j;
			nb++;
			if (nb == nodes.end()) {
				nb = nodes.begin();
			} 
			Vector v4=*( *nb ); 
			
			double denominator = 
			(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
			
			double ua = 
			((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
			double ub = 
			((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
			
			
			if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
				//cerr << "ua = " << ua << ", ub = " << ub << endl;
				return true;
			}
		}
    }
	
    return false;
}


bool Cell::MoveSelfIntersectsP(Node *moving_node_ind, Vector new_pos) {
	
	// Check whether the polygon will self-intersect if moving_node_ind 
	// were displaced to new_pos
	
	// Compare the two new edges against each other edge
	
	// O(2*N)
	
	// method used for segment intersection:
	// http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
	
	Vector neighbor_of_moving_node[2];
	
	//cerr << "list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);\n";
	list<Node *>::const_iterator moving_node_ind_pos = find (nodes.begin(),nodes.end(),moving_node_ind);
	
	list<Node *>::const_iterator nb = moving_node_ind_pos;
	//cerr << "Done\n";
	nb++;
	if (nb == nodes.end()) {
		nb = nodes.begin();
	} 
	
	neighbor_of_moving_node[0]=*(*nb); 
	
	nb=moving_node_ind_pos;
	if (nb == nodes.begin()) {
		nb = nodes.end();
	}
	nb--;
	
	neighbor_of_moving_node[1]=*( *nb ); 
	
	
	for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) {
		
		for (int j=0;j<2;j++) { // loop over the two neighbors of moving node
			list<Node *>::const_iterator nb=i;
			nb++;
			if (nb == nodes.end()) {
				nb = nodes.begin();
			} 
			if (*i == moving_node_ind || *nb == moving_node_ind) {
				// do not compare to self
				continue;
			}
			
			Vector v3 = *(*i);
			Vector v4 = *(*nb);
			
			double denominator = 
			(v4.y - v3.y)*(neighbor_of_moving_node[j].x - new_pos.x) - (v4.x - v3.x)*(neighbor_of_moving_node[j].y - new_pos.y);
			
			double ua = 
			((v4.x - v3.x)*(new_pos.y - v3.y) - (v4.y - v3.y)*(new_pos.x -v3.x))/denominator;
			double ub = 
			((neighbor_of_moving_node[j].x - new_pos.x)*(new_pos.y-v3.y) - (neighbor_of_moving_node[j].y- new_pos.y)*(new_pos.x - v3.x))/denominator;
			
			if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
				//cerr << "ua = " << ua << ", ub = " << ub << endl;
				return true;
			}
		}
	}
	
	return false;
}

/*! \brief Test if this cell intersects with the given line.
 
 */
bool Cell::IntersectsWithLineP(const Vector v1, const Vector v2) {
	
	
	// Compare the line against each edge
	
	// method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
	
	
	
	for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();
		 i++) 
    {
		
		Vector v3 = *(*i);
		list<Node *>::const_iterator nb=i;
		nb++;
		if (nb == nodes.end()) {
			nb = nodes.begin();
		}
		Vector v4 = *(*nb);
		
		double denominator = 
		(v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);
		
		double ua = 
		((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;
		double ub = 
		((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;
		
		if ( ( TINY < ua && ua < 1.-TINY ) && ( TINY < ub && ub < 1.-TINY ) ) {
			return true;
		}
    }
	
	return false;
	
	
}
/*! \brief Constructs Walls, but only one per cell boundary.
 
 Standard method constructs a Wall for each cell wall element,
 making transport algorithms computationally more intensive than needed.
 
 We can remove this? Well, let's leave it in the code in case we need it for something else. E.g. for importing leaf architectures in different formats than our own... :-)
 
 */
void Cell::ConstructWalls(void) {
	
	return;
	if (dead) return;
	
	walls.clear();
	neighbors.clear();
	
	// Get "corner points; i.e. nodes where more than 2 cells are connected
	list<Node *> corner_points;
	
	for (list<Node *>::const_iterator i=nodes.begin();
		 i!=nodes.end();i++) {
		
		// look for nodes belonging to >2 cells
		if ((*i)->owners.size()>2) {
			
			// push onto list
			corner_points.push_back(*i);
		}
		
	}
	
	// Construct Walls between corner points
	
	// previous one in list
	list<Node *>::const_iterator nb = (--corner_points.end());
	
	// loop over list, 
	for (list<Node *>::const_iterator i=corner_points.begin();
		 i!=corner_points.end(); ( i++, nb++) ) {
		
		if (nb==corner_points.end()) nb=corner_points.begin();
		// add owning cells to a list
		list<Cell *> owning_cells;
		Node &n(*(*i));
		
		for (list<Neighbor>::const_iterator j=n.owners.begin();
			 j!=n.owners.end();
			 j++) {
			owning_cells.push_back(j->cell);
		}
		
		Node &n2(*(*nb));
		for (list<Neighbor>::const_iterator j=n2.owners.begin();
			 j!=n2.owners.end();
			 j++) {
			owning_cells.push_back(j->cell);
		}
		
		// sort cell owners
		owning_cells.sort( mem_fun( &Cell::Cmp ));

		// find duplicates
		vector<Cell *> duplicates;
		list<Cell *>::const_iterator prevj = (--owning_cells.end());
		for (list<Cell *>::const_iterator j=owning_cells.begin();
			 j!=owning_cells.end();
			 ( j++, prevj++) ) {
			
			if (prevj==owning_cells.end()) prevj=owning_cells.begin();
			if (*j==*prevj) duplicates.push_back(*j);
			
		}
		
		
		if (duplicates.size()==3) { // ignore cell boundary (this occurs only after the first division, I think)
			vector<Cell *>::iterator dup_it = find_if(duplicates.begin(),duplicates.end(),mem_fun(&Cell::BoundaryPolP) );
			if (dup_it!=duplicates.end()) 
				duplicates.erase(dup_it);
			else {
				return;
			}
			
		}
		
		
		// One Wall for each neighbor, so we should be able to correctly construct neighbor lists here.
		if (duplicates[0]==this) {
			AddWall(  new Wall(*nb,*i,duplicates[0],duplicates[1]) );
			if (!duplicates[1]->BoundaryPolP()) {
				
				neighbors.push_back(duplicates[1]);
			}
		} else {
			AddWall ( new Wall(*nb,*i,duplicates[1],duplicates[0]) );
			if (!duplicates[0]->BoundaryPolP()) {
				neighbors.push_back(duplicates[0]);
				
			}
		}
	}
	
}


void BoundaryPolygon::Draw(QGraphicsScene *c, QString tooltip) {
	
	// Draw the BoundaryPolygon on a QCanvas object
	
	
	CellItem* p = new CellItem(this, c);
	
	QPolygonF pa(nodes.size());
	int cc=0;
	
	for (list<Node *>::const_iterator n=nodes.begin();
		 n!=nodes.end();
		 n++) {
		Node *i=*n;
		
		pa[cc++] = QPoint((int)((Offset().x+i->x)*Factor()),
						  (int)((Offset().y+i->y)*Factor()) );
	}
	
	
	p->setPolygon(pa);
	p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
	p->setBrush( Qt::NoBrush );
	p->setZValue(1);
	
	if (!tooltip.isEmpty())
		p->setToolTip(tooltip);
	
	p->show();
	
}

void Cell::Flux(double *flux, double *D)  {
	
	
	// loop over cell edges
	
	for (int c=0;c<NChem();c++) flux[c]=0.;
	
	for (list<Wall *>::iterator i=walls.begin();
		 i!=walls.end();
		 i++) {
		
		
		// leaf cannot take up chemicals from environment ("no flux boundary")
		if ((*i)->c2->BoundaryPolP()) continue;
		
		
		// flux depends on edge length and concentration difference
		for (int c=0;c<NChem();c++) {
			double phi = (*i)->length * ( D[c] ) * ( ((Cell *)(*i)->c2)->chem[c] - chem[c] );
			
			#ifdef QDEBUG
			if ((*i)->c1!=this) {
			  qDebug() << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl;
			}
			#endif
			
			flux[c] += phi;
		}    
	}
	
}


// graphics stuff, not compiled for batch versions
#ifdef QTGRAPHICS

#include "canvas.h"

void Cell::Draw(QGraphicsScene *c, QString tooltip) {
	
	// Draw the cell on a QCanvas object
	
	if (DeadP()) { 
	        #ifdef QDEBUG
	  qDebug() << "Cell " << index << " not drawn, because dead." << endl;
		#endif
		return;
	}
	
	CellItem* p = new CellItem(this, c);
	
	QPolygonF pa(nodes.size());
	int cc=0;
	
	for (list<Node *>::const_iterator n=nodes.begin();
		 n!=nodes.end();
		 n++) {
		Node *i=*n;
		
		pa[cc++] = QPoint((int)((offset[0]+i->x)*factor),
						  (int)((offset[1]+i->y)*factor) );
	}
	
	
	QColor cell_color;
	
	m->plugin->SetCellColor(this,&cell_color);
	
	p->setPolygon(pa);
	p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));
	p->setBrush( cell_color );
	p->setZValue(1);
	
	if (!tooltip.isEmpty())
		p->setToolTip(tooltip);
	
	p->show();
	
}


void Cell::DrawCenter(QGraphicsScene *c) const {
  // Maginfication derived similarly to that in nodeitem.cpp
  // Why not use Cell::Magnification()?
  const double mag = par.node_mag;
	
	// construct an ellipse
  QGraphicsEllipseItem *disk = new QGraphicsEllipseItem ( -1*mag, -1*mag, 2*mag, 2*mag, 0, c );
	disk->setBrush( QColor("forest green") );
	disk->setZValue(5);
	disk->show();
	Vector centroid=Centroid();
	disk -> setPos((offset[0]+centroid.x)*factor,(offset[1]+centroid.y)*factor);
}

void Cell::DrawNodes(QGraphicsScene *c) const {
	
	for (list<Node *>::const_iterator n=nodes.begin();
		 n!=nodes.end();
		 n++) {
		Node *i=*n;
		

		NodeItem *item = new NodeItem ( &(*i), c );
		item->setColor();
		item->setZValue(5);
		item->show();
		item ->setPos(((offset[0]+i->x)*factor),
					  ((offset[1]+i->y)*factor) );
	}
	
}

void Cell::DrawIndex(QGraphicsScene *c) const {
	
	DrawText( c, QString("%1").arg(index));
}

// Draw any text in the cell's center
void Cell::DrawText(QGraphicsScene *c, const QString &text) const {
    
	Vector centroid = Centroid();
	QGraphicsSimpleTextItem *ctext = new QGraphicsSimpleTextItem ( text, 0, c );
	ctext->setPen( QPen(QColor(par.textcolor)) );
	ctext->setZValue(20);
	ctext->setFont( QFont( "Helvetica", par.cellnumsize, QFont::Bold) );
	ctext->show();
	ctext ->setPos(((offset[0]+centroid.x)*factor),
				   ((offset[1]+centroid.y)*factor) );
	
}


void Cell::DrawAxis(QGraphicsScene *c) const {
	
	Vector long_axis;
	double width;
	Length(&long_axis, &width);
	
	//cerr << "Length is "  << length << endl;
	long_axis.Normalise();
	Vector short_axis=long_axis.Perp2D();
    

	Vector centroid = Centroid();
	Vector from = centroid - 0.5 * width * short_axis;
	Vector to = centroid + 0.5 * width *short_axis;

	
	QGraphicsLineItem *line = new QGraphicsLineItem(0, c);
	line->setPen( QPen(QColor(par.arrowcolor),2) );
	line->setZValue(2);
    
	line->setLine( ( (offset[0]+from.x)*factor ),
				  ( (offset[1]+from.y)*factor ), 
				  ( (offset[0]+to.x)*factor ),
				  ( (offset[1]+to.y)*factor ) );
	line->setZValue(10);
	line->show();
	
}

void Cell::DrawStrain(QGraphicsScene *c) const {
  c = NULL; // assignment merely to obviate compilation warning
	MyWarning::warning("Sorry, Cell::DrawStrain temporarily not implemented.");
}


void Cell::DrawFluxes(QGraphicsScene *c, double arrowsize)  {
	
	// get the mean flux through this cell
	Vector vec_flux = ReduceCellAndWalls<Vector>( PINdir );
	
	vec_flux.Normalise();
	
	vec_flux *= arrowsize;
	
	QGraphicsArrowItem *arrow = new QGraphicsArrowItem(0,c);
	
	Vector centroid = Centroid();
	Vector from = centroid - vec_flux/2.;
	Vector to = centroid + vec_flux/2.;
    
	
	arrow->setPen( QPen(QColor(par.arrowcolor),par.outlinewidth));
	arrow->setZValue(2);
	
	arrow->setLine( ( (offset[0]+from.x)*factor ),
				   ( (offset[1]+from.y)*factor ), 
				   ( (offset[0]+to.x)*factor ),
				   ( (offset[1]+to.y)*factor ) );
	arrow->setZValue(10);
	arrow->show();
	
}


void Cell::DrawWalls(QGraphicsScene *c) const {
	
	for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::Draw ) , c ) );
	
	// to see the cells connected the each wall (for debugging), uncomment the following
	//for_each(walls.begin(), walls.end(), bind2nd ( mem_fun ( &Wall::ShowStructure ), c ) );
}


void Cell::DrawValence(QGraphicsScene *c) const {
	
	DrawText(c, QString("%1").arg(walls.size()) );
	
}

#endif

/*! \brief Recalculate the lengths of the cell's Walls.
 
 Call this function after the Monte Carlo updates, and before doing the reaction-diffusion iterations.
 
 */
void Cell::SetWallLengths(void) {
	
	for (list<Wall *>::iterator de=walls.begin();
		 de!=walls.end();
		 de++) {
		
		// Step 1: find the path of nodes leading along the Wall.
		// A Wall often represents a curved cell wall: we want the total
		// length _along_ the wall here...
		
		
		// Locate first and second nodes of the edge in list of nodes
		list<Node *>::const_iterator first_node_edge = find(nodes.begin(), nodes.end(), (*de)->n1);
		list<Node *>::const_iterator second_node_edge_plus_1 = ++find(nodes.begin(), nodes.end(), (*de)->n2);
		
		double sum_length = 0.;
		
		// Now, walk to the second node of the edge in the list of nodes
		for (list<Node *>::const_iterator n=++first_node_edge;
			 n!=second_node_edge_plus_1;
			 ++n ) {
			
			if (n==nodes.end()) n=nodes.begin(); /* wrap around */ 
			
			
			list<Node *>::const_iterator prev_n = n; 
			if (prev_n==nodes.begin()) prev_n=nodes.end();
			--prev_n;
			
			
			// Note that Node derives from a Vector, so we can do vector calculus as defined in vector.h 
			sum_length += (*(*prev_n) - *(*n)).Norm(); 
			
			//cerr << "Node " << *prev_n << " to " << *n << ", cumulative length = " << sum_length << endl;
		}
		
		// We got the total length of the Wall now, store it:
		(*de)->length = sum_length;
		
		//cerr << endl;
		// goto next de
	}
}


//! Add Wall w to the list of Walls
void Cell::AddWall( Wall *w ) {
	
	// if necessary, we could try later inserting it at the correct position
        #ifdef QDEBUG
	if (w->c1 == w->c2 ){
	  qDebug() << "Wall between identical cells: " << w->c1->Index()<< endl;
	}
	#endif

	// Add Wall to Cell's list
	walls.push_back( w );
	
	// Add wall to Mesh's list if it isn't there yet
	
	if (find (
			  m->walls.begin(), m->walls.end(),
			  w )
		== m->walls.end() ) {
		m->walls.push_back(w);
	}
	
}

//! Remove Wall w from the list of Walls
list<Wall *>::iterator Cell::RemoveWall( Wall *w ) {
	
	// remove wall from Mesh's list
	m->walls.erase(
				   find(
						m->walls.begin(), m->walls.end(),
						w )
				   );
	
	// remove wall from Cell's list
	return 
	walls.erase (
				 find( 
					  walls.begin(), walls.end(),
					  w )
				 );
	
}



void Cell::EmitValues(double t) {
	
	//  cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl;
	emit ChemMonValue(t, chem);
	
}