# HG changeset patch
# User Roeland Merks <roeland.merks@cwi.nl>
# Date 2010-10-14 11:45:59
# Node ID 8aaffb5565a2a5d54e8632c7680d5d6dd9c60204
# Parent  fc4e80bf722db80709b7260236233457e91e3fa7

Corrected Compactness algorithm. Added algorithms to calculate circumference of convex hull and of the whole morph/

user: Roeland Merks <roeland.merks@cwi.nl>
branch 'default'

diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,18 @@
+2010-10-14    <merks@cwi.nl>
+
+	* cell.cpp (DivideWalls): accomodated for rename of Circumference -> WallCircumference
+
+	* hull.h: added an operator< to sort Points
+
+	* hull.cpp: added an operator< to sort Points
+
+	* cellbase.cpp (ExactCircumference): I added a new function ExactCircumference, yielding the circumference of the cell along its wall_elements
+
+	* VirtualLeaf.cpp: adjust info_string to accomodate for new name of function CellBase::Circumference -> CellBase::WallCircumference
+	* mesh.cpp: corrected Mesh::Compactness, the boundary coordinates need to be sorted in x,y order for the convex hull algorithm (thanks Margriet!). I updated CSVExportCellData so it exports the circumferences of hull and boundary_polygon.
+	
+
+
 2010-10-08    <guravage@caterpie.sen.cwi.nl>
 
 
diff --git a/src/VirtualLeaf.cpp b/src/VirtualLeaf.cpp
--- a/src/VirtualLeaf.cpp
+++ b/src/VirtualLeaf.cpp
@@ -120,7 +120,7 @@ public:
 		for (int i=0;i<Cell::NChem();i++) {
 			info_string += QString("%1 ").arg(c.Chemical(i));
 		}
-		info_string += QString("\nArea is %1\n Circumference is %2\n Boundary type is %3").arg(c.Area()).arg(c.Circumference()).arg(c.BoundaryStr());
+		info_string += QString("\nArea is %1\n Circumference is %2\n Boundary type is %3").arg(c.Area()).arg(c.WallCircumference()).arg(c.BoundaryStr());
 		
 	  info_string += "\nNodes: " + c.printednodelist();
 	  c.Draw(&canvas, info_string);
diff --git a/src/cellbase.cpp b/src/cellbase.cpp
--- a/src/cellbase.cpp
+++ b/src/cellbase.cpp
@@ -626,4 +626,27 @@ QString CellBase::printednodelist(void)
   return info_string;
 }
 
+double CellBase::ExactCircumference(void) const
+{
+
+  // simply sum length of all edges
+  double circumference=0.;
+    
+  for (list<Node *>::const_iterator i=nodes.begin(); i!=(nodes.end()); i++) {
+
+    list<Node *>::const_iterator i_plus_1=i; i_plus_1++;
+    if (i_plus_1==nodes.end())
+      i_plus_1=nodes.begin();
+
+    double dx=((*i_plus_1)->x-(*i)->x);
+    double dy=((*i_plus_1)->y-(*i)->y);
+    double l=sqrt(dx*dx+dy*dy);
+    //    f << (*i)->x << " " << (*i)->y << " " << (*i_plus_1)->x << " " << (*i_plus_1)->y << " " << l << endl;
+
+    circumference += l;
+  }
+
+  return circumference;
+} 
+
 /* finis*/
diff --git a/src/cellbase.h b/src/cellbase.h
--- a/src/cellbase.h
+++ b/src/cellbase.h
@@ -175,7 +175,7 @@ class CellBase :  public QObject, public
   double Length(Vector *long_axis = 0, double *width = 0) const;
   double CalcLength(Vector *long_axis = 0, double *width = 0) const;
 
-
+  double ExactCircumference(void) const;
   inline int Index(void) const { return index; }
 
 
@@ -203,7 +203,7 @@ class CellBase :  public QObject, public
     flag_for_divide = true;
   }
 
-  inline double Circumference(void) const {
+  inline double WallCircumference(void) const {
     double sum=0.;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
diff --git a/src/hull.cpp b/src/hull.cpp
--- a/src/hull.cpp
+++ b/src/hull.cpp
@@ -24,6 +24,18 @@ isLeft( Point P0, Point P1, Point P2 )
 {
     return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
 }
+
+// required to sort points (e.g. Qt's qSort())
+bool operator<(const Point & p1, const Point & p2) {
+  if (p1.y<p2.y) return true;
+  else
+    if (p1.y>p2.y) return false;
+    else 
+      if (p1.x<p2.x) return true;
+      else
+	return false;
+}
+
 //===================================================================
  
 
diff --git a/src/hull.h b/src/hull.h
--- a/src/hull.h
+++ b/src/hull.h
@@ -1,3 +1,6 @@
+#ifndef _HULL_H_
+#define _HULL_H_
+
 // Class point needed by 2D convex hull code
 class Point {
 
@@ -11,7 +14,11 @@ public:
   }
   float x,y;
 
+
 };
 
+// required to sort points (e.g. Qt's qSort())
+bool operator<(const Point & p1, const Point & p2);
+
 int chainHull_2D( Point* P, int n, Point* H );
-
+#endif
diff --git a/src/mesh.cpp b/src/mesh.cpp
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1970,29 +1970,42 @@ void Mesh::StandardInit(void) {
 
 #include "hull.h"
 
-double Mesh::Compactness(double *res_compactness, double *res_area, double *res_cell_area) {
+
+double Mesh::Compactness(double *res_compactness, double *res_area, double *res_cell_area, double *res_circumference) {
   
   // Calculate compactness using the convex hull of the cells
   // We use Andrew's Monotone Chain Algorithm (see hull.cpp)
 
   // Step 1. Prepare data for 2D hull code - get boundary polygon
   int pc=0;
-  Point *p=new Point[boundary_polygon->nodes.size()];
+  Point *p=new Point[boundary_polygon->nodes.size()+1];
   for (list<Node *>::const_iterator i = boundary_polygon->nodes.begin(); 
        i!=boundary_polygon->nodes.end(); i++) {
     p[pc++]=Point((*i)->x,(*i)->y);
   }
   
+  // chainHull algorithm requires sorted points
+  qSort( p, p+pc );
+
+ 
   // Step 2: call 2D Hull code
   int np=boundary_polygon->nodes.size();
-  Point *hull=new Point[np];
+  Point *hull=new Point[np+1];
   int nph=chainHull_2D(p,np,hull);
   
   
-  // Step 3: calculate area of convex hull
+  // Step 3: calculate area and circumference of convex hull
   double hull_area=0.;
+  double hull_circumference=0.;
+
   for (int i=0;i<nph-1;i++) {
     hull_area+=hull[i].x * hull[i+1].y - hull[i+1].x * hull[i].y;
+    double s_dx=(hull[i+1].x-hull[i].x);
+    double s_dy=(hull[i+1].y-hull[i].y);
+    double l=sqrt(s_dx*s_dx+s_dy*s_dy);
+    //    f << hull[i].x << " " << hull[i].y << " " << hull[i+1].x << " " << hull[i+1].y << " " << l << endl;
+    hull_circumference+=l;
+      
   }
   hull_area/=2.;
 
@@ -2020,7 +2033,10 @@ double Mesh::Compactness(double *res_com
   if (res_cell_area) {
     *res_cell_area = boundary_pol_area;
   }
-
+  if (res_circumference) {
+    *res_circumference = hull_circumference;
+  }
+  
   // return compactness
   return boundary_pol_area/hull_area;
 
@@ -2054,11 +2070,12 @@ void Mesh::CSVExportCellData(QTextStream
 
 void Mesh::CSVExportMeshData(QTextStream &csv_stream) { 
   
-  csv_stream << "\"Mesh area\",\"Number of cells\",\"Number of nodes\",\"Compactness\",\"Hull area\",\"Cell area\"" << endl;
+  csv_stream << "\"Morph area\",\"Number of cells\",\"Number of nodes\",\"Compactness\",\"Hull area\",\"Morph circumference\",\"Hull circumference\"" << endl;
   
-  double res_compactness, res_area, res_cell_area;
-  Compactness(&res_compactness, &res_area, &res_cell_area);
-  csv_stream << Area() << ", " << NCells() << ", " << NNodes() << ", " << res_compactness << ", " << res_area << ", " << res_cell_area  << endl;
+  double res_compactness, res_area, res_cell_area, hull_circumference;
+  Compactness(&res_compactness, &res_area, &res_cell_area, &hull_circumference);
+  double morph_circumference = boundary_polygon->ExactCircumference();
+  csv_stream << Area() << ", " << NCells() << ", " << NNodes() << ", " << res_compactness << ", " << res_area << ", " << morph_circumference << ", " << hull_circumference << endl;
   
 }
 /* finis */
diff --git a/src/mesh.h b/src/mesh.h
--- a/src/mesh.h
+++ b/src/mesh.h
@@ -385,7 +385,7 @@ class Mesh {
   }
   QString ModelID(void) { return plugin?plugin->ModelID():QString("undefined"); }
   void StandardInit(void);	
-  double Compactness(double *res_compactness=0, double *res_area=0, double *res_cell_area=0);
+  double Compactness(double *res_compactness=0, double *res_area=0, double *res_cell_area=0, double *hull_circumference=0);
   void CSVExportCellData(QTextStream &csv_stream) const;
   void CSVExportMeshData(QTextStream &csv_stream);