EI/VirtualLeaf Changeset - 8aaffb5565a2 · Centrum Wiskunde & Informatica (CWI)

Changeset - 8aaffb5565a2

Parent rev.

Child rev.

[Not reviewed]

default

0 8 0

Roeland Merks - 15 years ago 2010-10-14 11:45:59
roeland.merks@cwi.nl

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'

8 files changed with 88 insertions and 14 deletions:

src/ChangeLog

src/VirtualLeaf.cpp

src/cellbase.cpp

src/cellbase.h

src/hull.cpp

src/hull.h

src/mesh.cpp

src/mesh.h

0 comments (0 inline, 0 general)

src/ChangeLog

➞

Show inline comments

 -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.
 -10-08    <guravage@caterpie.sen.cwi.nl>
 	* pardialog.h:
 	* pardialog.cpp:
 	* parameter.h:
 	* parameter.cpp: Regenerated to include export_interval and export_fn_prefix.
 	* VirtualLeafpar.tmpl: Appended export_interval and export_fn_prefix.
 	* canvas.h (MainBase): Declared polymorphic exportCellData() functions.
 	* canvas.cpp:
 	(TimeStamp): New private TimeStamp() function.
 	(TimeStepWrap): Added invocation of exportCellData().
 	(exportCellData): Created two polymorphic functions: one with a
 	single QString argument, the other with no argument. The former is
 	called from TimeStepWrap() while the latter is called from the
 	"Export cell areas" item in the file menu.
 -10-07    <guravage@caterpie.sen.cwi.nl>
 	* canvas.cpp (exportCellData): Added a Q3FileDialog to inquire
 	where to write the exportCellData.
 -06-28    <guravage@caterpie.sen.cwi.nl>
 	* VirtualLeaf-install.nsi: Grab gpl3.txt from doc directory.
 	* canvas.cpp (gpl): gpl3.txt can be either in an ancestor doc
 	directory (Linux) or a decedent doc directory (Windows, via the
 	binary installer).
 	* VirtualLeaf-install.nsi: Add VirtualLeaf doc directory.
 -06-25    <guravage@caterpie.sen.cwi.nl>
 	* gpl3.txt: Moved gpl3.txt from doc to src directory.
 	* VirtualLeaf.pro: Added -Wno-write-strings and -Wno-unused-parameter to QMAKE_CXXFLAGS.
 	* libplugin.pro: Ditto.
 	* parameter.cpp: Result of adding datadir changes to make_parameter_source.pl.
 	* parameter.h: Ditto.
 	* output.h: Declared new function (AppendHomeDirIfPathRelative).
 	* output.cpp (AppendHomeDirIfPathRelative): Added new function.
 	* canvas.cpp (gpl): Moving gpl3.txt from doc to src obviates the need to docDir.cd("../doc").
 	* VirtualLeaf-install.nsi: Grab gpl3.txt from src directory.
 	Add missing libiconv/bin, libxml2bin and libz/bin directories.
 	Copy libiconv-2.dll, libxml2.dll and zlib1.dll from relative paths.
 	* VirtualLeaf.pro: copy gpl3.txt as part of QMAKE_POST_LINK.
 -06-24    <guravage@caterpie.sen.cwi.nl>
 	* libplugin.pro: Use correct library path.
 	* VirtualLeaf.pro: Ditto.
 	* VirtualLeaf.cpp (DrawCell): Iterate over NChem to construct info_string.
 -06-23    <guravage@caterpie.sen.cwi.nl>
 	* simitembase.cpp: Removed NULL assignments to unused variables.
 	* VirtualLeaf.cpp: Ditto.
 	* apoplastitem.cpp: Ditto.
 	* canvas.cpp: Ditto.
 	* cell.cpp: Ditto.
 	* cellbase.h: Ditto.
 	* forwardeuler.cpp: Ditto.
 	* mainbase.h: Ditto.
 	* nodeitem.cpp: Ditto.
 	* qcanvasarrow.h: Ditto.
 	* simitembase.cpp: Ditto.
 	* Makefile (clean): Add -f Makefile argument to each make invocation.
 	* VirtualLeaf-install.nsi: New gpl license text.
 	* VirtualLeaf.pro: Disabled console mode.
 	* mesh.cpp (Clear): Added parentheses to qDebug statments.
 	(TestIllegalWalls): Replaced qDebug().
 	* canvas.cpp (mouseReleaseEvent): Replaced qDebug() with cerr since qDebug complains about *node_set.
 	* wall.cpp (CorrectWall): Rplaced gDebug() with cerr in transform call and when printing *this.
 -06-22    <guravage@caterpie.sen.cwi.nl>
 	* Makefile (tutorials): Add tutorials target.
 -06-21    <guravage@caterpie.sen.cwi.nl>
 	* parameter.cpp: Added particular reassignment of datadir.
 	* canvas.cpp (gpl): Added GPL3 License text. Display detail text only if the source text file exists.
 -06-18    <guravage@caterpie.sen.cwi.nl>
 	* canvas.cpp (gpl): Added gpl slot to display GPL license.
 	* VirtualLeaf.pro: Changed default LIBXML2DIR, LIBICONVDIR and LIBZDIR to corresponding distribution lib directories.
 	* libplugin.pro: Ditto.
 	* Makefile (clean): add if stmt not to `touch` on windows.
 -06-17    <guravage@caterpie.sen.cwi.nl>
 	* VirtualLeaf.pro: Removed perl references.
 	* libplugin.pro: Ditto.
 -06-15    <guravage@caterpie.sen.cwi.nl>
 	* VirtualLeaf.pro: Removed xmlwritecode.cpp from SOURCES list.
 	* xmlwrite.cpp (XMLSave): Removed references to XMLWriteLeafSourceCode and XMLWriteReactionsCode.
 	* xmlwrite.h (XMLIO): Ditto!
 	* mesh.cpp (findNextBoundaryNode): Initialize Node *next_boundary_node = NULL;
 	* xmlwrite.cpp (XMLReadSimtime): Removed unused variable cur
 	(XMLReadWalls): viz_flux need not be declared twice; default value of 0.0.
 	(XMLReadCells): Removed unused count variable.
 	(XMLReadSimtime): Removed unused cur variable.
 	(XMLRead): Removed unused v_str variable.
 	* simitembase.cpp (userMove): Use assignment merely to obviate compilation warning.
 	(SimItemBase) Ditto.
 	* qcanvasarrow.h (QGraphicsLineItem): Use assignment merely to obviate compilation warning.
 	* output.cpp (OpenWriteFile): Removed unused par variable.
 	* nodeitem.cpp (paint): Use assignment merely to obviate compilation warning.
 	* forwardeuler.cpp (odeint): Use assignment merely to obviate compilation warning.
 	* cell.cpp (DivideOverGivenLine): Use assignment merely to obviate compilation warning.
 	* canvas.cpp (FigureEditor): Use assignments merely to obviate compilation errors.
 	(mousePressEvent): Removed unused item variable.
 	* apoplastitem.cpp
 	(ApoplastItem): Removed unused par variable.
 	(OnClick): Use NULL assignment merely to obviate compilation warning.
 	* mainbase.h (MainBase): Use assignment merely to obviate compilation warning.
 	* cellbase.h (CellsStaticDatamembers): Use assignment merely to obviate compilation warning.
 	* cell.cpp: Wrapped diagnostic output in QDEBUG blocks.
 	* VirtualLeaf.cpp ditto.
 	* canvas.cpp ditto.
 	* cell.cpp ditto.
 	* data_plot.cpp ditto.
 	* forwardeuler.cpp ditto.
 	* mesh.cpp ditto.
 	* mesh.h
 	* random.cpp ditto.
 	* wall.cpp ditto.
 	* wallbase.cpp ditto.
 	* wallitem.cpp ditto.
 -06-07    <guravage@caterpie.sen.cwi.nl>
 	* VirtualLeaf.pro: Removed explicit perl invocation to regerenerate parameter files.
 	* libplugin.pro: ditto.
 -06-03    <guravage@caterpie.sen.cwi.nl>
 	* pardialog.h: Added default versions of this automatically generated file.
 	* pardialog.cpp: ditto.
 	* parameter.h: ditto.
 	* parameter.cpp: ditto.
 	* VirtualLeaf.pro: delete/generate  parameter.{h,cpp}and pardialog.{h,cpp} only if perl is installed.
  	* libplugin.pro: dito.
 	* Makefile: Added top-level Makefile
 -05-10    <guravage@caterpie.sen.cwi.nl>

src/VirtualLeaf.cpp

➞

Show inline comments

 /*
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  VirtualLeaf 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.
+ *
  *  VirtualLeaf 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 <string>
 #include <fstream>
 #include <sstream>
 #include <cstring>
 #include <functional>
 #include <getopt.h>
 #include <cerrno>
 #include "mesh.h"
 #include "parameter.h"
 #include "random.h"
 #include "pi.h"
 #include "cellitem.h"
 #include "canvas.h"
 #include "cell.h"
 #include "output.h"
 #include <qwidget.h>
 #include <q3process.h>
 #include <qapplication.h>
 #include <QDesktopWidget>
 #include <QGraphicsScene>
 #include <QMessageBox>
 //Added by qt3to4:
 #include <QMouseEvent>
 #include <unistd.h>
 #include <q3textstream.h>
 #ifdef HAVE_QWT
 #include "data_plot.h"
 #endif
 #include <QPalette>
 #include <QBrush>
 #include <QToolTip>
 #include "simplugin.h"
 #include <QPluginLoader>
 #include <QDir>
 #include "modelcatalogue.h"
 static const std::string _module_id("$Id$");
 extern Parameter par;
 MainBase *main_window = 0;
 #ifdef XFIGGRAPHICS
 #define TIMESTEP double Graphics::TimeStep(void)
 #endif
 class PrintNode {
 public:
   void operator() (const Node &n) const
+  {
     cerr << n.Index() << ": " << n <<  endl;
+  }
 };
 class EdgeSource {
 public:
   void operator() (Cell &c) {
     if (c.AtBoundaryP()) {
       cerr << "Cell " << c.Index() << " is a source cell.\n";
       c.SetSource(0,par.source);
     } else {
       cerr << "Cell " << c.Index() << " is _not_ a source cell.\n";
+    }
+  }
 };
 class CellInfo {
 public:
   void operator() (Cell &c,std::ostream &os) const {
     os << "Cell " << c.index << " says: " << endl;
     os << "c.nodes.size() = " << c.nodes.size() << endl;
     for (list<Node *>::iterator i=c.nodes.begin(); i!=c.nodes.end(); i++) {
       cerr << (*i)->Index() << " ";
+    }
     cerr << endl;
+  }
 };
 double PINSum(Cell &c) {
   return c.Chemical(1) + c.SumTransporters(1);// + c.ReduceCellAndWalls<double>( complex_PijAj );
+}
 class DrawCell {
 public:
   void operator() (Cell &c,QGraphicsScene &canvas, MainBase &m) const {
     if (m.ShowBorderCellsP() || c.Boundary()==Cell::None) {
       if (!m.ShowBoundaryOnlyP() && !m.HideCellsP()) {
 	if (m.ShowToolTipsP()) {
 	  //QString info_string=QString("Cell %1, chemicals: ( %2, %3, %4, %5, %6)\n %7 of PIN1 at walls.\n Area is %8\n PIN sum is %9\n Circumference is %10\n Boundary type is %11").arg(c.Index()).arg(c.Chemical(0)).arg(c.Chemical(1)).arg(c.Chemical(2)).arg(c.Chemical(3)).arg(c.Chemical(4)).arg(c.SumTransporters(1)).arg(c.Area()).arg(PINSum(c)).arg(c.Circumference()).arg(c.BoundaryStr());
 		QString info_string=QString("Cell %1, chemicals(%2): ").arg(c.Index()).arg(Cell::NChem());
 		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);
 	} else {
 	  c.Draw(&canvas);
+	}
+      }
       if (m.ShowCentersP()){
 	c.DrawCenter(&canvas);
+      }
       if (m.ShowFluxesP()){
 	c.DrawFluxes(&canvas, par.arrowsize);
+      }
+    }
+  }
 };
 Mesh mesh;
 bool batch=false;
 void MainBase::Plot(int resize_stride)
+{
   clear();
   static int count=0;
   if (resize_stride) {
     if ( !((count)%resize_stride) ) {
       FitLeafToCanvas();
+    }
+  }
   mesh.LoopCells(DrawCell(),canvas,*this);
   if (ShowNodeNumbersP())
     mesh.LoopNodes( bind2nd (mem_fun_ref ( &Node::DrawIndex), &canvas ) ) ;
   if (ShowCellNumbersP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawIndex), &canvas ) ) ;
   if (ShowCellAxesP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawAxis), &canvas ) );
   if (ShowCellStrainP())
     mesh.LoopCells( bind2nd (mem_fun_ref ( &Cell::DrawStrain), &canvas ) );
   if (ShowWallsP())
     mesh.LoopWalls( bind2nd( mem_fun_ref( &Wall::Draw ), &canvas ) );
 /*  if (ShowApoplastsP())
     mesh.LoopWalls( bind2nd( mem_fun_ref( &Wall::DrawApoplast ), &canvas ) );
 */
   if (ShowMeshP())
     mesh.DrawNodes(&canvas);
   if (ShowBoundaryOnlyP())
     mesh.DrawBoundary(&canvas);
   if ( ( batch || MovieFramesP() )) {
     static int frame = 0;
     // frame numbers are sequential for the most frequently written file type.
     // for the less frequently written file type they match the other type
     if (!(count%par.storage_stride) )  {
       stringstream fname;
       fname << par.datadir << "/leaf.";
       fname.fill('0');
       fname.width(6);
       fname << frame << ".jpg";
       if (par.storage_stride <= par.xml_storage_stride) {
 	frame++;
+      }
       // Write high-res JPG snapshot every plot step
       Save(fname.str().c_str(), "JPEG",1024,768);
+    }
     if (!(count%par.xml_storage_stride)) {
       stringstream fname;
       fname << par.datadir << "/leaf.";
       fname.fill('0');
       fname.width(6);
       fname << frame << ".xml";
       if (par.xml_storage_stride < par.storage_stride) {
 	frame++;
+      }
       // Write XML file every ten plot steps
       mesh.XMLSave(fname.str().c_str(), XMLSettingsTree());
+    }
+  }
  count++;
+}
 INIT {
   //mesh.SetSimPlugin(plugin);
   if (leaffile) {
     xmlNode *settings;
     mesh.XMLRead(leaffile, &settings);
     main_window->XMLReadSettings(settings);
     xmlFree(settings);
     main_window->UserMessage(QString("Ready. Time is %1").arg(mesh.getTimeHours().c_str()));
   } else {
     mesh.StandardInit();
+  }
   Cell::SetMagnification(1);
   Cell::setOffset(0,0);
   FitLeafToCanvas();
   Plot();
+}
 TIMESTEP {
   static int i=0;
   static int t=0;
   static int ncells;
   if (!batch) {
     UserMessage(QString("Time: %1").arg(mesh.getTimeHours().c_str()),0);
+  }
   ncells=mesh.NCells();
   double dh;
   if(DynamicCellsP()) {
     dh = mesh.DisplaceNodes();
     // Only allow for node insertion, cell division and cell growth
     // if the system has equillibrized
     // i.e. cell wall tension equillibrization is much faster
     // than biological processes, including division, cell wall yielding
     // and cell expansion
     mesh.InsertNodes(); // (this amounts to cell wall yielding)
     if ( (-dh) < par.energy_threshold) {
       mesh.IncreaseCellCapacityIfNecessary();
       mesh.DoCellHouseKeeping();
       //mesh.LoopCurrentCells(mem_fun(&plugin->CellHouseKeeping)); // this includes cell division
       // Reaction diffusion
       mesh.ReactDiffuse(par.rd_dt);
       t++;
       Plot(par.resize_stride);
+    }
   } else {
     mesh.ReactDiffuse(par.rd_dt);
     Plot(par.resize_stride);
+  }
   i++;
   return mesh.getTime();
+}
 /* Called if a cell is clicked */
 void Cell::OnClick(QMouseEvent *e){}
 /* Custom message handler - Default appends a newline character to the end of each line. */
 void vlMessageOutput(QtMsgType type, const char *msg)
+{
   switch (type) {
   case QtDebugMsg:
     //fprintf(stderr, "Debug: %s\n", msg);
     cerr << msg << flush;
     break;
   case QtWarningMsg:
     //fprintf(stderr, "Warning: %s\n", msg);
     cerr << "Warning: " << msg << flush;
     break;
   case QtCriticalMsg:
     fprintf(stderr, "Critical: %s\n", msg);
     cerr << "Critical: " << msg << flush;
     break;
   case QtFatalMsg:
     //fprintf(stderr, "Fatal: %s\n", msg);
     cerr << "Fatal: " << msg << flush;
     abort();
+  }
+}

src/cellbase.cpp

➞

Show inline comments

@@ @@ -437,193 +437,216 @@ double CellBase::CalcLength(Vector *long @@
   // see file inertiatensor.nb for explanation of this method
   double my_intgrl_xx=0., my_intgrl_xy=0., my_intgrl_yy=0.;
   double my_intgrl_x=0., my_intgrl_y=0., my_area=0.;
   my_area=0.;
   list<Node *>::const_iterator nb;
   list<Node *>::const_iterator i=nodes.begin();
   for (; i!=(nodes.end()); i++) {
     nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();
     my_area+=(*i)->x*(*nb)->y;
     my_area-=(*nb)->x*(*i)->y;
     my_intgrl_xx+=
       ((*i)->x*(*i)->x+
        (*nb)->x*(*i)->x+
        (*nb)->x*(*nb)->x ) *
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     my_intgrl_xy+=
       ((*nb)->x*(*i)->y-
        (*i)->x*(*nb)->y)*
       ((*i)->x*(2*(*i)->y+(*nb)->y)+
        (*nb)->x*((*i)->y+2*(*nb)->y));
     my_intgrl_yy+=
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y)*
       ((*i)->y*(*i)->y+
        (*nb)->y*(*i)->y+
        (*nb)->y*(*nb)->y );
     my_intgrl_x+=
       ((*nb)->x+(*i)->x)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
     my_intgrl_y+=
       ((*nb)->y+(*i)->y)*
       ((*i)->x*(*nb)->y-
        (*nb)->x*(*i)->y);
+  }
   //my_area/=2.0;
   my_area = fabs(my_area)/2.0;
   double intrx=my_intgrl_x/6.;
   double intry=my_intgrl_y/6.;
   double ixx=(my_intgrl_xx/12.)-(intrx*intrx)/my_area;
   double ixy=(my_intgrl_xy/24.)+(intrx*intry)/my_area;
   double iyy=(my_intgrl_yy/12.)-(intry*intry)/my_area;
   double rhs1=(ixx+iyy)/2., rhs2=sqrt( (ixx-iyy)*(ixx-iyy)+4*ixy*ixy )/2.;
   double lambda_b=rhs1+rhs2;
   // see: http://scienceworld.wolfram.com/physics/MomentofInertiaEllipse.html
   //    cerr << "n = " << n << "\n";
   if (long_axis) {
     *long_axis = Vector(-ixy, lambda_b - ixx, 0);
     //   cerr << "ixx = " << ixx << ", ixy = " << ixy << ", iyy = " << iyy << ", my_area = " << my_area << endl;
+  }
   if (width) {
     *width = 4*sqrt((rhs1-rhs2)/my_area);
+  }
   return 4*sqrt(lambda_b/my_area);
+}
 void CellBase::ConstructNeighborList(void)
+{
   neighbors.clear();
   for (//list<Wall *>::const_reverse_iterator wit=walls.rbegin();
        list<Wall *>::const_iterator wit=walls.begin();
        // somehow the reverse_iterator returns by walls needs to be casted to const to let this work.
        // it seems to me it is a bug in the STL implementation...
        wit!=walls.end();
        wit++) {
     if ((*wit)->C1() != this) {
       neighbors.push_back((*wit)->C1());
     } else {
       neighbors.push_back((*wit)->C2());
+    }
+  }
   // remove all boundary_polygons from the list
   list <CellBase *>::iterator e=neighbors.begin();
   at_boundary=false;
   do {
     // Code crashes here after cutting off part of the leaf. I can't find the problem.
     // Leaving the "Illegal" walls in the simulation helps. (c1=-1 && c2=-1)
     // Work-around: define leaf primordium. Save to XML. Restart. Read XML file.
     // Sorry about this; I hope to solve this annoying issue later. RM :-).
     // All cells in neighbors seem to be okay (I might be messing some part of the memory elsewhere
     // during the cutting operation?).
     e = find_if(neighbors.begin(),neighbors.end(),mem_fun(&CellBase::BoundaryPolP));
     if (e!=neighbors.end()) {
       e=neighbors.erase(e);
       at_boundary=true;
     } else {
       break;
+    }
   } while(1);
+}
 // Save the cell to a stream so we can reconstruct its state later
 void CellBase::Dump(ostream &os) const
+{
   os << index << " " << nodes.size() << endl;
   Vector::Dump(os);
   os << endl;
   for (list<Node *>::const_iterator i=nodes.begin();i!=nodes.end();i++) {
     os << *i << " ";
+  }
   os << endl;
   os << index << " " << neighbors.size() << endl;
   for (list<CellBase *>::const_iterator i=neighbors.begin();i!=neighbors.end();i++) {
     os << *i << " ";
+  }
   os << endl << walls.size() << endl << endl;
   os << NChem() << " ";
   for (int i=0;i<NChem();i++) {
     os << chem[i] << " ";
+  }
   os << endl;
   os << NChem() << " ";
   for (int i=0;i<NChem();i++) {
     os << new_chem[i] << " ";
+  }
   os << endl;
   os << boundary << " " << area << " " << target_area << " " << target_length
      << " " << fixed << " " << intgrl_xx << " " << intgrl_xy << " " << intgrl_yy
      << " " << intgrl_x << " " << intgrl_y << " " << source << " ";
   cellvec.Dump(os);
   os << " " << source_conc << " " << source_chem;
   os << endl;
+}
 void CellBase::UnfixNodes(void)
+{
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->Unfix();
+  }
+}
 void CellBase::FixNodes(void)
+{
   for (list<Node *>::const_iterator i=nodes.begin(); i!=nodes.end(); i++) {
     (*i)->Fix();
+  }
+}
 // returns true if cell is at border
 bool CellBase::AtBoundaryP(void) const
+{
   return at_boundary;
+}
 QString CellBase::printednodelist(void)
+{
   QString info_string = "Nodelist = { ";
   for (list<Node *>::const_iterator i =  nodes.begin(); i!=nodes.end(); i++) {
     info_string += QString("%1 ").arg((*i)->Index());
+  }
   info_string += " } ";
   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*/

src/cellbase.h

➞

Show inline comments

 /*
+ *
  *  $Id$
+ *
  *  This file is part of the Virtual Leaf.
+ *
  *  VirtualLeaf 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.
+ *
  *  VirtualLeaf 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.
+ *
  */
 // CellBase derives from Vector, where Vector is simply used as a Vertex
 #ifndef _CELLBASE_H_
 #define _CELLBASE_H_
 #include <list>
 #include <vector>
 #include <iostream>
 #include <QString>
 #include <QDebug>
 #include "vector.h"
 #include "parameter.h"
 #include "wall.h"
 #include "warning.h"
 #include "assert.h"
 extern Parameter par;
 using namespace std;
 class Mesh;
 class Node;
 class CellBase;
 class NodeSet;
 struct ParentInfo {
   Vector polarization;
   double PINmembrane;
   double PINendosome;
 };
 // We need a little trick here, to make sure the plugin and the main application will see the same static datamembers
 // otherwise each have their own instantation.
 // My solution is as follows. I collect all original statics in a class. The main application instantiates it and
 // has a static pointer to it. After loading the plugin I set a static pointer to the same class
 class CellsStaticDatamembers {
  public:
   CellsStaticDatamembers(void) {
     ncells = 0;
     nchem = 0;
     base_area = 0.;
+  }
   ~CellsStaticDatamembers() {
 #ifdef QDEBUG
     qDebug() << "Oops! Desctructor of CellsStaticDatamembers called" << endl;
 #endif
+  }
   int ncells;
   int nchem;
   double base_area;
 };
 class CellBase :  public QObject, public Vector
+{
   Q_OBJECT
   friend class Mesh;
   friend class CellInfo;
   friend class Node;
   friend class WallBase;
   friend class SimPluginInterface;
  public:
   CellBase(QObject *parent=0);
   CellBase(double x,double y,double z=0); // constructor
   virtual ~CellBase() {
     delete[] chem;
     delete[] new_chem;
     if (division_axis) delete division_axis;
     //cerr << "CellBase " << index << " is dying. " << endl;
+  }
   CellBase(const CellBase &src); // copy constructor
   virtual bool BoundaryPolP(void) const { return false; }
   CellBase operator=(const CellBase &src); // assignment operator
   CellBase operator=(const Vector &src);
   void SetChemical(int chem, double conc);
   inline void SetNewChem(int chem, double conc)
+  {
     new_chem[chem] = conc;
+  }
   void SetSource(int chem, double conc)
+  {
     source=true;
     source_chem = chem;
     source_conc = conc;
+  }
   // set chem 1 to conc in all membranes of this cell
   void SetTransporters(int chem, double conc);
   void UnfixNodes(void);
   void FixNodes(void);
   void UnsetSource(void) {
     source = false;
+  }
   inline bool Source(void) { return source; }
   enum boundary_type {None, Noflux, SourceSink, SAM};
   static const char * boundary_type_names[4];
   inline const char *BoundaryStr(void) { return boundary_type_names[boundary]; }
   ostream &print(ostream &os) const;
   inline double Chemical(int c) const { // returns the value of chemical c
 #ifdef _undefined_
     qDebug() << endl << "Entering cellbase::chemical(" << c << "), and nchem is: " << NChem() << "." << endl;
 #endif
     int nchem = NChem();
 #ifdef _undefined_
     if ((c<0) || (c>=nchem))
       MyWarning::warning("CellBase::Chemical says: index c is: %d, but nchem is: %d. Merely return zero", c, nchem);
 #endif
     return ((c<0) || (c>=nchem)) ? 0 : chem[c];
+  }
   //void print_nblist(void) const;
   boundary_type SetBoundary(boundary_type bound)
+  {
     if (bound!=None) {
+    }
     return boundary=bound;
+  }
   boundary_type ResetBoundary(void) { return boundary=None; }
   boundary_type Boundary(void) const { return boundary; }
   static int &NChem(void) { return static_data_members->nchem; }
   double CalcArea(void) const;
   double RecalcArea(void) { return area = CalcArea(); }
   Vector Centroid(void) const;
   void SetIntegrals(void) const;
   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; }
   void SetTargetArea(double tar_ar) { target_area=tar_ar; }
   inline void SetTargetLength(double tar_l) { target_length=tar_l; }
   inline void SetLambdaLength(double lambda_length) { lambda_celllength = lambda_length; }
   inline double TargetArea(void) { return target_area; }
   inline void SetStiffness(double stiff) { stiffness = stiff; }
   inline double Stiffness(void) { return stiffness; }
   inline double EnlargeTargetArea(double da) { return target_area+=da; }
   inline double Area(void) const { return area; }
   inline void Divide(void) { flag_for_divide = true; }
   inline void DivideOverAxis(const Vector &v)
+  {
     division_axis = new Vector(v);
     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();
 	 w++) {
       sum +=  (*w)->Length();
+    }
     return sum;
+  }
   QList<WallBase *> getWalls(void) {
     QList<WallBase *> wall_list;
     for (list<Wall *>::iterator i=walls.begin();
 	 i!=walls.end();
 	 i++) {
       wall_list << *i;
+    }
     return wall_list;
+  }
   void Dump(ostream &os) const;
   QString printednodelist(void);
   inline bool DeadP(void) { return dead; }
   inline void MarkDead(void) { dead  = true; }
   static double &BaseArea(void)
+  {
     return static_data_members->base_area;
+  }
   void CheckForDivision(void);
   // write flux from neighboring cells into "flux"
   void Flux(double *flux, double *D);
   inline bool FixedP(void) { return fixed; }
   inline bool Fix(void) {  FixNodes(); return (fixed=true); }
   inline bool Unfix(void) { UnfixNodes(); return (fixed=false);}
   inline void setCellVec(Vector cv) { cellvec = cv; }
   bool AtBoundaryP(void) const;
   static inline int &NCells(void)
+  {
     return static_data_members->ncells;
+  }
   inline void Mark(void)
+  {
     marked=true;
+  }
   inline void Unmark(void)
+  {
     marked=false;
+  }
   inline bool Marked(void) const {
     return marked;
+  }
   //! Returns the sum of chemical "chem" of this CellBase's neighbors
   double SumChemicalsOfNeighbors(int chem)
+  {
     double sum=0.;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum +=  (*w)->Length() * ( (*w)->c1!=this ? (*w)->c1->Chemical(chem) : (*w)->c2->Chemical(chem) );
+    }
     return sum;
+  }
   //! Generalization of the previous member function
   template<class P, class Op> P ReduceNeighbors(Op f) {
     P sum=0;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += (*w)->c1 != this ? f( *((*w)->c1) ) : f ( *((*w)->c2) );
+    }
     return sum;
+  }
   //! The same, but now for the walls
   template<class P, class Op> P ReduceWalls(Op f, P sum) {
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += f( **w );
+    }
     return sum;
+  }
   //! The same, but now for the walls AND neighbors
   template<class P, class Op> P ReduceCellAndWalls(Op f)
+  {
     P sum = 0;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += ((*w)->c1 == this) ?
 	f( ((*w)->c1), ((*w)->c2), *w ) :
 	f( ((*w)->c2), ((*w)->c1), *w );
+    }
     return sum;
+  }
   //! Sum transporters at this CellBase's side of the walls
   double SumTransporters(int ch)
+  {
     double sum=0.;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += (*w)->getTransporter(this, ch);
+    }
     return sum;
+  }
   inline int NumberOfDivisions(void) { return div_counter; }
   //! Sum transporters at this CellBase's side of the walls
   double SumLengthTransporters(int ch)
+  {
     double sum=0.;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += (*w)->getTransporter(this, ch) * (*w)->Length();
+    }
     return sum;
+  }
   double SumLengthTransportersChemical(int trch, int ch)
+  {
     double sum=0.;
     for (list<Wall *>::const_iterator w=walls.begin();
 	 w!=walls.end();
 	 w++) {
       sum += (*w)->getTransporter(this, trch) * ( (*w)->c1!=this ? (*w)->c1->Chemical(ch) : (*w)->c2->Chemical(ch) );
+    }
     return sum;
+  }
   inline int CellType(void) const { return cell_type; }
   inline void SetCellType(int ct) { cell_type = ct; }
   static void SetNChem(int new_nchem)
+  {
     if (NCells()) {
       MyWarning::error("CellBase::SetNChem says: not permitted, call SetNChem after deleting all cells.");
     } else {
       NChem() = new_nchem;
+    }
+  }
   inline double TargetLength() const { return target_length; }
   static inline CellsStaticDatamembers *GetStaticDataMemberPointer(void) { return static_data_members; }
  protected:
   // (define a list of Node* iterators)
   typedef list < list<Node *>::iterator > ItList;
   int index;
   inline void SetChemToNewchem(void)
+  {
     for (int c=0;c<CellBase::NChem();c++) {
       chem[c]=new_chem[c];
+    }
+  }
   inline void SetNewChemToChem(void)
+  {
     for (int c=0;c<CellBase::NChem();c++) {
       new_chem[c]=chem[c];
+    }
+  }
   inline double NewChem(int c) const { return new_chem[c]; }
  protected:
   list<Node *> nodes;
   void ConstructNeighborList(void);
   long wall_list_index (Wall *elem) const;

src/hull.cpp

➞

Show inline comments

 // Copyright 2001, softSurfer (www.softsurfer.com)
 // This code may be freely used and modified for any purpose
 // providing that this copyright notice is included with it.
 // SoftSurfer makes no warranty for this code, and cannot be held
 // liable for any real or imagined damage resulting from its use.
 // Users of this code must verify correctness for their application.
 // Assume that a class is already given for the object:
 //    Point with coordinates {float x, y;}
 //===================================================================
 // isLeft(): tests if a point is Left|On|Right of an infinite line.
 //    Input:  three points P0, P1, and P2
 //    Return: >0 for P2 left of the line through P0 and P1
 //            =0 for P2 on the line
 //            <0 for P2 right of the line
 //    See: the January 2001 Algorithm on Area of Triangles
 #include "hull.h"
 inline float
 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;
+}
 //===================================================================
 // chainHull_2D(): Andrew's monotone chain 2D convex hull algorithm
 //     Input:  P[] = an array of 2D points
 //                   presorted by increasing x- and y-coordinates
 //             n = the number of points in P[]
 //     Output: H[] = an array of the convex hull vertices (max is n)
 //     Return: the number of points in H[]
 int
 chainHull_2D( Point* P, int n, Point* H )
+{
     // the output array H[] will be used as the stack
     int    bot=0, top=(-1);  // indices for bottom and top of the stack
     int    i;                // array scan index
     // Get the indices of points with min x-coord and min|max y-coord
     int minmin = 0, minmax;
     float xmin = P[0].x;
     for (i=1; i<n; i++)
         if (P[i].x != xmin) break;
     minmax = i-1;
     if (minmax == n-1) {       // degenerate case: all x-coords == xmin
         H[++top] = P[minmin];
         if (P[minmax].y != P[minmin].y) // a nontrivial segment
             H[++top] = P[minmax];
         H[++top] = P[minmin];           // add polygon endpoint
         return top+1;
+    }
     // Get the indices of points with max x-coord and min|max y-coord
     int maxmin, maxmax = n-1;
     float xmax = P[n-1].x;
     for (i=n-2; i>=0; i--)
         if (P[i].x != xmax) break;
     maxmin = i+1;
     // Compute the lower hull on the stack H
     H[++top] = P[minmin];      // push minmin point onto stack
     i = minmax;
     while (++i <= maxmin)
+    {
         // the lower line joins P[minmin] with P[maxmin]
         if (isLeft( P[minmin], P[maxmin], P[i]) >= 0 && i < maxmin)
             continue;          // ignore P[i] above or on the lower line
         while (top > 0)        // there are at least 2 points on the stack
+        {
             // test if P[i] is left of the line at the stack top
             if (isLeft( H[top-1], H[top], P[i]) > 0)
                 break;         // P[i] is a new hull vertex
             else
                 top--;         // pop top point off stack
+        }
         H[++top] = P[i];       // push P[i] onto stack
+    }
     // Next, compute the upper hull on the stack H above the bottom hull
     if (maxmax != maxmin)      // if distinct xmax points
         H[++top] = P[maxmax];  // push maxmax point onto stack
     bot = top;                 // the bottom point of the upper hull stack
     i = maxmin;
     while (--i >= minmax)
+    {
         // the upper line joins P[maxmax] with P[minmax]
         if (isLeft( P[maxmax], P[minmax], P[i]) >= 0 && i > minmax)
             continue;          // ignore P[i] below or on the upper line
         while (top > bot)    // at least 2 points on the upper stack
+        {
             // test if P[i] is left of the line at the stack top
             if (isLeft( H[top-1], H[top], P[i]) > 0)
                 break;         // P[i] is a new hull vertex
             else
                 top--;         // pop top point off stack
+        }
         H[++top] = P[i];       // push P[i] onto stack
+    }
     if (minmax != minmin)
         H[++top] = P[minmin];  // push joining endpoint onto stack
     return top+1;
+}

src/hull.h

➞

Show inline comments

 #ifndef _HULL_H_
 #define _HULL_H_
 // Class point needed by 2D convex hull code
 class Point {
 public:
   Point(float xx, float yy) {
     x=xx;
     y=yy;
+  }
   Point(void) {
     x=0; y=0;
+  }
   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

src/mesh.cpp

➞

Show inline comments

@@ @@ -1781,284 +1781,301 @@ void Mesh::setValues(double x, double *y @@
   emit_count++;
+}
 double *Mesh::getValues(int *neqs) {
   int nwalls = walls.size();
   int ncells = cells.size();
   int nchems = Cell::NChem();
   // two eqs per chemical for each wall, and one eq per chemical for each cell
   // This is for generality. For a specific model you may optimize
   // this by removing superfluous (empty) equations.
   (*neqs) = 2 * nwalls * nchems + ncells * nchems;
   // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
   static double *values = 0;
   if (values!=0) { delete[] values; }
   values = new double[*neqs];
   int i=0;
   for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*c)->Chemical(ch);
+    }
     i+=nchems;
+  }
   for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*w)->Transporters1(ch);
+    }
     i+=nchems;
     for (int ch=0;ch<nchems;ch++) {
       values[i+ch]=(*w)->Transporters2(ch);
+    }
     i+=nchems;
+  }
   return values;
+}
 void Mesh::DrawNodes(QGraphicsScene *c) const {
   for (vector<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(((Cell::offset[0]+i->x)*Cell::factor),
 		  ((Cell::offset[1]+i->y)*Cell::factor) );
+  }
+}
 /*! Returns the sum of protein "ch" of a cycling protein in cells and walls */
 double Mesh::CalcProtCellsWalls(int ch) const {
   double sum_prot=0.;
   // At membranes
   for (list<Wall *>::const_iterator w=walls.begin(); w!=walls.end(); w++) {
     sum_prot += (*w)->Transporters1(ch);
     sum_prot += (*w)->Transporters2(ch);
+  }
   // At cells
   for (vector<Cell *>::const_iterator c=cells.begin(); c!=cells.end(); c++) {
     sum_prot += (*c)->Chemical(ch);
+  }
   return sum_prot;
+}
 void Mesh::SettoInitVals(void) {
   vector<double> clean_chem(Cell::NChem());
   vector<double> clean_transporters(Cell::NChem());
   for (int i=0;i<Cell::NChem();i++) {
     clean_transporters[i]=0.;
     clean_chem[i]=par.initval[i];
+  }
   CleanChemicals(clean_chem);
   CleanTransporters(clean_transporters);
+}
 string Mesh::getTimeHours(void) const {
   int hours = (int)(time / 3600);
   int mins = (int)((time - hours * 3600)/60);
   int secs = (int)((time - hours * 3600 - mins * 60));
   ostringstream tstr;
   tstr << hours << " h " << mins << " m " << secs << " s";
   return tstr.str();
+}
 QVector<qreal> Mesh::VertexAngles(void) {
   QVector<qreal> angles;
   for (vector<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
       angles+=(*n)->NeighbourAngles();
+    }
+  }
   return angles;
+}
 QVector< QPair<qreal,int> > Mesh::VertexAnglesValues(void) {
   QVector< QPair<qreal,int> > anglesvalues;
   for (vector<Node *>::const_iterator n=nodes.begin(); n!=nodes.end(); n++) {
     if ((*n)->Value()>2 && !(*n)->BoundaryP() ) {
       QVector<qreal> angles = (*n)->NeighbourAngles();
       int value_vertex = angles.size();
       for (QVector<qreal>::ConstIterator i=angles.begin(); i!=angles.end(); i++) {
 	anglesvalues += QPair< qreal, int > (*i, value_vertex);
+      }
+    }
+  }
   return anglesvalues;
+}
 void Mesh::Clean(void) {
 #ifdef QDEBUG
   qDebug() << "Freeing nodes" << endl;
 #endif
   for (vector<Node *>::iterator i=nodes.begin(); i!=nodes.end(); i++) {
     delete *i;
+  }
   nodes.clear();
   Node::nnodes=0;
 #ifdef QDEBUG
   qDebug() << "Freeing node sets" << endl;
 #endif
   for (vector<NodeSet *>::iterator i=node_sets.begin(); i!=node_sets.end(); i++) {
     delete *i;
+  }
   node_sets.clear();
 #ifdef QDEBUG
   qDebug() << "Freeing cells" << endl;
 #endif
   //CellsStaticDatamembers *old_static_data_mem = Cell::GetStaticDataMemberPointer();
   for (vector<Cell *>::iterator i=cells.begin(); i!=cells.end(); i++) {
     delete *i;
+  }
   //Cell::static_data_members = old_static_data_mem;
   cells.clear();
   Cell::NCells()=0;
   if (boundary_polygon) {
     delete boundary_polygon; // (already deleted during cleaning of cells?)
     boundary_polygon=0;
+  }
 #ifdef QDEBUG
   qDebug() << "Freeing walls" << endl;
 #endif
   for (list<Wall *>::iterator i=walls.begin(); i!=walls.end(); i++) {
     delete *i;
+  }
   walls.clear();
   Wall::nwalls=0;
   node_insertion_queue.clear();
   shuffled_nodes.clear();
   shuffled_cells.clear();
   time = 0.0;
+}
 void Mesh::StandardInit(void) {
   boundary_polygon = new BoundaryPolygon();
   Cell &circle=CircularCell(0,0,10,10);
   circle.SetTargetArea(circle.CalcArea());
   circle.SetTargetLength(par.target_length);
   circle.SetLambdaLength(par.lambda_celllength);
   SetBaseArea();
   // clean up chemicals
   for (int c=0; c<Cell::NChem(); c++) {
     circle.SetChemical(c, 0.);
+  }
+}
 #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.;
   // Step 4: get area of bounary polygon
   double boundary_pol_area = boundary_polygon->CalcArea();
   /*  ofstream datastr("hull.dat");
   for (int i=0;i<nph<i++) {
     datastr << hull.x << " " << hull.y << endl;
+  }
   ofstream polstr("pol.dat");
   for (int i=0;i<np;h*/
   delete[] p;
   delete[] hull;
   // put intermediate results into optional pointers
   if (res_compactness) {
     *res_compactness = boundary_pol_area/hull_area;
+  }
   if (res_area) {
     *res_area = hull_area;
+  }
   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;
+}
 // DataExport
 void Mesh::CSVExportCellData(QTextStream &csv_stream) const {
   csv_stream << "\"Cell Index\",\"Center of mass (x)\",\"Center of mass (y)\",\"Cell area\",\"Cell length\"";
   for (int c=0;c<Cell::NChem(); c++) {
     csv_stream << ",\"Chemical " << c << "\"";
+  }
   csv_stream << endl;
   for (vector<Cell *>::const_iterator i=cells.begin();
        i!=cells.end();
        i++) {
     Vector centroid = (*i)->Centroid();
     csv_stream << (*i)->Index() << ", "
 	       << centroid.x << ", "
 	       << centroid.y << ", "
 	       <<  (*i)->Area() << ", "
 	       <<(*i)->Length();
     for (int c=0;c<Cell::NChem(); c++) {
       csv_stream << ", " << (*i)->Chemical(c);
+    }
     csv_stream << endl;
+  }
+}
 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 */

src/mesh.h

➞

Show inline comments

@@ @@ -196,241 +196,241 @@ class Mesh { @@
     for (vector<Cell *>::iterator i=cells.begin();
 	 i!=cells.end();
 	 i++) {
       f(**i,g,h);
+    }
+  }
   void DoCellHouseKeeping(void) {
     vector<Cell *> current_cells = cells;
     for (vector<Cell *>::iterator i = current_cells.begin();
 	 i != current_cells.end();
 	 i ++) {
       plugin->CellHouseKeeping(*i);
       // Call functions of Cell that cannot be called from CellBase, including Division
       if ((*i)->flag_for_divide) {
 	if ((*i)->division_axis) {
 	  (*i)->DivideOverAxis(*(*i)->division_axis);
 	  delete (*i)->division_axis;
 	  (*i)->division_axis = 0;
 	} else {
 	  (*i)->Divide();
+	}
 	(*i)->flag_for_divide=false;
+      }
+    }
+  }
   // Apply "f" to cell i
   // i.e. this is an adapter which allows you to call a function
   // operating on Cell on its numeric index index
   template<class Op> void cell_index_adapter(Op f,int i) {
     f(cells[i]);
+  }
   double DisplaceNodes(void);
   void BoundingBox(Vector &LowerLeft, Vector &UpperRight);
   int NEqs(void) {     int nwalls = walls.size();
     int ncells =cells.size();
     int nchems = Cell::NChem();
     // two eqs per chemical for each walls, and one eq per chemical for each cell
     // This is for generality. For a specific model you may optimize
     // this by removing superfluous (empty) equations.
     int neqs = 2 * nwalls * nchems + ncells * nchems;
     return neqs;
+  }
   void IncreaseCellCapacityIfNecessary(void) {
     return;
     // cerr << "Entering Mesh::IncreaseCellCapacityIfNecessary \n";
     // make sure we always have enough space
     // to have each cell divide at least once
     //
     // Note that we must do this, because Cell::Divide pushes a new Cell
     // onto Mesh::cells. As a result, Mesh::cells might be relocated
     // if we are _within_ a Cell object: i.e. pointer "this" will be changed!!
     //
     // An alternative solution could be to make "Mesh::cells" a list,
     // but this won't work because we need random access for
     // the Monte Carlo algorithm.
     if (2*cells.size()>cells.capacity()) {
       cerr << "Increasing capacity to "  << 2*cells.capacity() << endl;
       cerr << "Current capacity is " << cells.capacity() << endl;
       cells.reserve(cells.capacity()*2);
+    }
+  }
   void ReserveMoreCells(int n) {
     if (nodes.size()+n>nodes.capacity()) {
       nodes.reserve(size()+n);
+    }
+  }
   double Area(void);
   double MeanArea(void) {
     double sum=0.;
     for (vector<Cell *>::const_iterator i=cells.begin();
 	 i!=cells.end();
 	 i++) {
       sum+=(*i)->Area();
+    }
     return sum/(double)NCells();
+  }
   void SetBaseArea(void);
   int NCells(void) const {
     return cells.size();
+  }
   inline int NNodes(void) const {
     return nodes.size();
+  }
   void PrintQueue(ostream &os) {
     while (!node_insertion_queue.empty()) {
       os << node_insertion_queue.front() << endl;
       node_insertion_queue.pop();
+    }
+  }
   void InsertNodes(void) {
     // insert the nodes in the insertion queue
     while (!node_insertion_queue.empty()) {
       //cerr << node_insertion_queue.front() << endl;
       InsertNode(node_insertion_queue.front());
       node_insertion_queue.pop();
+    }
+  }
   void Clear();
   void ReactDiffuse( double delta_t = 1 );
   double SumChemical(int ch);
   void SetChemical(int ch, double value) {
     for (vector<Cell *>::iterator c=cells.begin();
 	 c!=cells.end();
 	 c++) {
       (*c)->chem[ch]=value;
+    }
+  }
   // used for interacing with ODE-solvers (e.g. NRCRungeKutta)
   void setValues(double x, double *y);
   double *getValues(int *neqs);
   void Derivatives(double *derivs);
 #ifdef QTGRAPHICS
   inline void DrawBoundary(QGraphicsScene *c) {
     boundary_polygon->Draw(c);
+  }
   void DrawNodes(QGraphicsScene *c) const;
 #endif
   double max_chem;
   void XMLSave(const char *docname, xmlNode *settings=0) const;
   void XMLRead(const char *docname, xmlNode **settings=0, bool geometry = true, bool pars = true, bool simtime = true);
   void XMLReadPars(const xmlNode * root_node);
   void XMLReadGeometry(const xmlNode *root_node);
   void XMLReadSimtime(const xmlNode *root_node);
   void XMLReadNodes(xmlNode *cur);
   void XMLReadCells(xmlNode *cur);
   void XMLParseTree(const xmlNode * root_node);
   void XMLReadWalls(xmlNode *cur, vector<Wall *> *tmp_cells);
   void XMLReadWallsToCells(xmlNode *root, vector<Wall *> *tmp_walls);
   void XMLReadNodeSets(xmlNode *root);
   void XMLReadNodeSetsToNodes(xmlNode *root);
   void PerturbChem(int chemnum, double range);
   void CleanUpCellNodeLists(void);
   void CleanUpWalls(void);
   void CutAwayBelowLine( Vector startpoint, Vector endpoint );
   void CutAwaySAM(void);
   void RepairBoundaryPolygon(void);
   void Rotate(double angle, Vector center);
   void PrintWallList( void );
   void TestIllegalWalls(void);
   Vector FirstConcMoment(int chem);
   inline Vector Centroid(void) {
     return boundary_polygon->Centroid();
+  }
   inline Vector Offset(void) {
     return boundary_polygon->Offset();
+  }
   inline double Factor(void) {
     return boundary_polygon->Factor();
+  }
   void DeleteLooseWalls(void);
   void FitLeafToCanvas(double width, double height);
   void AddNodeSet(NodeSet *node_set) {
     node_sets.push_back(node_set);
+  }
   void CleanChemicals(const vector<double> &clean_chem);
   void CleanTransporters(const vector<double> &clean_transporters);
   void RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters);
   inline double getTime(void) const { return time; }
   string getTimeHours(void) const;
   inline void setTime(double t) { time = t; }
   double CalcProtCellsWalls(int ch) const;
   void SettoInitVals(void);
   QVector<qreal> VertexAngles(void);
   QVector< QPair<qreal,int> > VertexAnglesValues(void);
   void SetSimPlugin(SimPluginInterface *new_plugin) {
     plugin=new_plugin;
+  }
   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);
   Node* findNextBoundaryNode(Node*);
  private:
   // Data members
   vector<Cell *> cells;
   vector<Node *> nodes;
   list<Wall *> walls; // we need to erase elements from this container frequently, hence a list.
  public:
   vector<NodeSet *> node_sets;
  private:
   vector<Node *> shuffled_nodes;
   vector<Cell *> shuffled_cells;
   unique_queue<Edge> node_insertion_queue;
   BoundaryPolygon *boundary_polygon;
   double time;
   SimPluginInterface *plugin;
   // Private member functions
   void AddNodeToCell(Cell *c, Node *n, Node *nb1 , Node *nb2);
   void AddNodeToCellAtIndex(Cell *c, Node *n, Node *nb1 , Node *nb2, list<Node *>::iterator ins_pos);
   void InsertNode(Edge &e);
   inline Node *AddNode(Node *n) {
     nodes.push_back(n);
     shuffled_nodes.push_back(n);
     n->m=this;
     return n;
+  }
   inline Cell *AddCell(Cell *c) {
     cells.push_back(c);
     shuffled_cells.push_back(c);
     //cerr << "Shuffled cell indices:  ";
     /*copy(shuffled_cells.begin(),shuffled_cells.end(),ostream_iterator<int>(cerr," "));
       cerr << endl;*/
     c->m=this;
     return c;
+  }
   void CircumCircle(double x1,double y1,double x2,double y2,double x3,double y3,
 		    double *xc,double *yc,double *r);
 };
 #endif
 /* finis */

0 comments (0 inline, 0 general)