Changeset - 2accd83924da
[Not reviewed]
default
0 20 0
Roeland Merks - 15 years ago 2010-11-25 14:57:33
roeland.merks@cwi.nl
Some last fine-tuning
- File dialogs start from appropriate data directories
- XML files contain right values for datadir
- Reordered parameter dialog
- Added new "make perl" heading to the makefile
- Corrected segfault for batch mode

user: Roeland Merks <roeland.merks@cwi.nl>
branch 'default'
changed data/leaves/auxin_growth.xml
changed data/leaves/meinhardt_init.xml
changed data/leaves/tutorial1_init.xml
changed data/leaves/tutorial2_init.xml
changed data/leaves/tutorial3_init.xml
changed data/leaves/tutorial4_init.xml
changed data/leaves/tutorial5_init.xml
changed src/Makefile
changed src/VirtualLeaf.cpp
changed src/VirtualLeafpar.tmpl
changed src/canvas.cpp
changed src/canvas.h
changed src/mesh.cpp
changed src/mesh.h
changed src/modelcatalogue.cpp
changed src/modelcatalogue.h
changed src/parameter.cpp
changed src/parameter.h
changed src/pardialog.cpp
changed src/pardialog.h
13 files changed:
0 comments (0 inline, 0 general)
data/leaves/auxin_growth.xml
Show inline comments
 
binary diff not shown
data/leaves/meinhardt_init.xml
Show inline comments
 
binary diff not shown
data/leaves/tutorial1_init.xml
Show inline comments
 
binary diff not shown
data/leaves/tutorial2_init.xml
Show inline comments
 
binary diff not shown
data/leaves/tutorial3_init.xml
Show inline comments
 
binary diff not shown
data/leaves/tutorial4_init.xml
Show inline comments
 
binary diff not shown
data/leaves/tutorial5_init.xml
Show inline comments
 
binary diff not shown
src/Makefile
Show inline comments
 
# $Id$
 

	
 
# Set an environment varible $(MAKE) to the name of the make command
 
# on your machine, e.g. 'make' on UNIX and 'mingw32-make' on windows.
 
# Fallback default is UNIX make.
 

	
 
ifeq ($(MAKE),)
 
	MAKE=make
 
endif
 

	
 
# Ditto for QMAKE
 
ifeq ($(QMAKE),)
 
	QMAKE=qmake
 
endif
 

	
 
all: VirtualLeaf libplugin plugins tutorials
 

	
 
VirtualLeaf: Makefile.VirtualLeaf
 
	$(MAKE) -f Makefile.VirtualLeaf
 

	
 
Makefile.VirtualLeaf: VirtualLeaf.pro
 
	$(QMAKE) -o $@ $< 
 

	
 
libplugin: Makefile.libplugin
 
	$(MAKE) -f Makefile.libplugin
 

	
 
Makefile.libplugin: libplugin.pro
 
	$(QMAKE) -o $@ $< 
 

	
 
plugins:
 
	$(MAKE) -C build_models -f Makefile
 

	
 
tutorials:
 
	$(MAKE) -C build_models -f Makefile
 
	$(MAKE) -C TutorialCode -f Makefile
 

	
 
perl: Makefile.perl
 
	$(MAKE) -f Makefile.perl
 

	
 
clean:
 
	$(MAKE) -f Makefile.libplugin clean
 
	$(MAKE) -f Makefile.VirtualLeaf clean
 
	$(MAKE) -C build_models -f Makefile clean
 
	$(MAKE) -C TutorialCode -f Makefile clean
 
ifeq ($(MAKE),make)
 
	touch VirtualLeaf.pro
 
	touch libplugin.pro
 
else
 
	copy /b VirtualLeaf.pro +,,
 
	copy /b libplugin.pro +,,
 
endif
 

	
 
#finis
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.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();
 

	
 
  int count=(int)mesh.getTime();
 

	
 
  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() )) {
 
    if (!(count%par.storage_stride) ) {
 
      stringstream fname;
 
      fname << par.datadir << "/leaf.";
 
      fname.fill('0');
 
      fname.width(6);
 
      fname << count << ".png";
 
      // Write high-res PNG snapshot every plot step
 
      Save(fname.str().c_str(), "PNG", 1024, 768);
 
    }
 

	
 
    if (!(count%par.xml_storage_stride)) {
 
      stringstream fname;
 
      fname << par.datadir << "/leaf.";
 
      fname.fill('0');
 
      fname.width(6);
 
      fname << count << ".xml";
 
      // Write XML file every ten plot steps
 
      mesh.XMLSave(fname.str().c_str(), XMLSettingsTree());
 
    }
 
  }
 
}
 

	
 

	
 
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()));
 

	
 
    if (qApp->type()==QApplication::Tty) {
 
      
 
      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 {
 
      ((Main *)main_window)->readStateXML(leaffile);
 
    }
 
    
 
  } 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();
 
  }
 
}
 

	
 

	
 
Parameter par;
 

	
 
int main(int argc,char **argv) {
 

	
 
  try {
 
    int c;
 
    char *leaffile=0;
 
    char *modelfile=0;
 

	
 
    while (1) {
 

	
 
      //int this_option_optind = optind ? optind : 1;
 
      int option_index = 0;
 
      static struct option long_options[] = {
 
	{"batch", no_argument, NULL, 'b'},
 
	{"leaffile", required_argument, NULL, 'l'},
 
	{"model", required_argument, NULL, 'm'} 
 
      };
 

	
 
      // short option 'p' creates trouble for non-commandline usage on MacOSX. Option -p changed to -P (capital)
 
      static char *short_options = "blm";
 
      c = getopt_long (argc, argv, "bl:m:",
 
		       long_options, &option_index);
 
      if (c == -1)
 
	break;
 

	
 

	
 
      if (c==0) {
 
	printf ("option %s", long_options[option_index].name);
 
	if (optarg)
 
	  printf (" with arg %s", optarg);
 
	printf ("\n");
 

	
 
	c = short_options[option_index];
 
      }
 

	
 
      switch (c) {
 
      case 'b':
 
	cerr << "Running in batch mode\n";
 
	batch=true;
 
	break;
 

	
 
      case 'l':
 
	leaffile=strdup(optarg);
 
	if (!leaffile) {
 
	  throw("Out of memory");
 
	}
 
	printf("Reading leaf state file '%s'\n", leaffile);
 
	break;
 

	
 
      case 'm':
 
	modelfile=strdup(optarg);
 
	if (!modelfile) {
 
	  throw("Out of memory");
 
	}
 
	break;
 

	
 
      case '?':
 
	break;
 

	
 
      default:
 
	printf ("?? getopt returned character code 0%o ??\n", c);
 
      }
 
    }
 

	
 

	
 
    if (optind < argc) {
 
      printf ("non-option ARGV-elements: ");
 
      while (optind < argc)
 
	printf ("%s ", argv[optind++]);
 
      printf ("\n");
 
    }
 

	
 
    bool useGUI = !batch;
 
    qInstallMsgHandler(vlMessageOutput); // custom message handler
 
    QApplication app(argc,argv,useGUI);
 

	
 

	
 

	
 
    QPalette tooltippalette = QToolTip::palette();
 
    QColor transparentcolor = QColor(tooltippalette.brush(QPalette::Window).color());
 

	
 
    tooltippalette.setBrush (QPalette::Window, QBrush (transparentcolor) );
 
    QToolTip::setPalette( tooltippalette );
 

	
 
    QGraphicsScene canvas(0,0,8000,6000);
 

	
 

	
 
    if (useGUI) {
 
      main_window=new Main(canvas, mesh);
 
      if ( QApplication::desktop()->width() > ((Main *)main_window)->width() + 10
 
	   && QApplication::desktop()->height() > ((Main *)main_window)->height() +30 ) {
 

	
 
	((Main *)main_window)->show();
 
	((Main *)main_window)->resize( ((Main *)main_window)->sizeHint());
 
      } else {
 
        ((Main *)main_window)->showMaximized();
 
      }
 
      
 
      // show "About" window at start up
 
      ((Main *)main_window)->about();
 
    } else {
 
      main_window=new MainBase(canvas, mesh);
 
    }
 

	
 
    canvas.setSceneRect(QRectF());
 
    if (!batch) {
 
      QObject::connect( qApp, SIGNAL(lastWindowClosed()), qApp, SLOT(quit()) );
 
    }
 

	
 
    //    main_window->Init(leaffile);
 

	
 
    // Install model or read catalogue of models
 
    ModelCatalogue model_catalogue(&mesh, useGUI?(Main *)main_window:0,modelfile);
 
    ModelCatalogue model_catalogue(&mesh, main_window,modelfile);
 

	
 

	
 
    if (useGUI)
 
      model_catalogue.PopulateModelMenu();
 
    model_catalogue.InstallFirstModel();
 
    
 
    
 

	
 
    /*    Cell::SetMagnification(1);
 
    Cell::setOffset(0,0);
 

	
 
    main_window->FitLeafToCanvas();
 

	
 
    main_window->Plot();
 
    */
 
    if (batch) {
 
      double t=0.;
 
      do {
 
	t = main_window->TimeStep();
 
      } while (t < par.maxt);
 
    } else
 
      return app.exec();
 

	
 
  } catch (const char *message) {
 
    if (batch) { 
 
      cerr << "Exception caught:" << endl;
 
      cerr << message << endl;
 
      abort();
 
    } else {
 
      QString qmess=QString("Exception caught: %1").arg(message);
 
      QMessageBox::critical(0, "Critical Error", qmess);
 
      abort();
 
    }
 
  } catch (ios_base::failure) {
 
    stringstream error_message;
 
    error_message << "I/O failure: " << strerror(errno);
 
    if (batch) {
 
      cerr << error_message.str() <<endl;
 
      abort();
 
    } else {
 
      QString qmess(error_message.str().c_str());
 
      QMessageBox::critical(0, "I/O Error", qmess );
 
      abort();
 
    }
 
  }
 
}
 

	
 
/* finis */
src/VirtualLeafpar.tmpl
Show inline comments
 
title = Parameter values for The Virtual Leaf / title
 
label = <b>Visualization</b> / label
 
arrowcolor = white / string
 
arrowsize = 100 / double
 
textcolor = red / string
 
cellnumsize = 1 / int
 
nodenumsize = 1 / int
 
node_mag = 1.0 / double 
 
outlinewidth = 1.0 / double
 
cell_outline_color = forestgreen / string
 
resize_stride = 0 / int
 
label = / label
 
label = <b>Data Export</b> / label
 
export_interval = 0 / int
 
export_fn_prefix = cell. / string
 
storage_stride = 10 / int
 
xml_storage_stride = 500 / int
 
datadir = . / directory 
 
label = / label
 
label = <b>Cell mechanics</b> / label
 
T = 1.0 / double
 
lambda_length = 100. / double
 
yielding_threshold = 4. / double
 
lambda_celllength = 0. / double
 
target_length = 60. / double 
 
cell_expansion_rate = 1. / double
 
cell_div_expansion_rate = 0. / double
 
auxin_dependent_growth = true / bool
 
ode_accuracy = 1e-4 / double
 
mc_stepsize = 0.4 / double
 
mc_cell_stepsize = 0.2 / double
 
energy_threshold = 1000. / double
 
bend_lambda = 0. / double
 
alignment_lambda = 0. / double
 
rel_cell_div_threshold = 2. / double
 
rel_perimeter_stiffness = 2 / double
 
collapse_node_threshold = 0.05 / double
 
morphogen_div_threshold = 0.2 / double
 
morphogen_expansion_threshold = 0.01 / double
 
copy_wall = true / bool
 
label = / label
 
label = <b>Auxin transport and PIN1 dynamics</b> / label
 
source = 0. / double
 
D = 0., 0.0, 0.0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 
initval = 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 
k1 = 1. / double
 
k2 = 0.3 / double
 
r = 1. / double
 
kr = 1. / double
 
#k3 = 1 / double
 
#k4 = 0.3 / double
 
km = 1. / double
 
Pi_tot = 1. / double
 
transport = 0.036 / double
 
ka = 1 / double
 
# d = 1e-06 / double
 
# e = 0.01 / double
 
# f = 10 / double
 
pin_prod = 0.001 / double
 
pin_prod_in_epidermis = 0.1 / double
 
pin_breakdown = 0.001 / double
 
pin_breakdown_internal = 0.001 / double
 
aux1prod = 0.001 / double
 
aux1prodmeso = 0. / double
 
aux1decay = 0.001 / double
 
aux1decaymeso = 0.1 / double
 
aux1transport = 0.036 / double
 
aux_cons = 0. / double
 
aux_breakdown = 0. / double
 
kaux1 = 1 / double
 
kap = 1 / double
 
leaf_tip_source = 0.001 / double
 
sam_efflux = 0.0001 / double
 
sam_auxin = 10. / double
 
# gf_prod = 1e-3 / double
 
# gf_decay = 1e-3 / double
 
sam_auxin_breakdown = 0 / double
 
#label = / label
 
#label = <b>miscellaneous</b> / label
 
van3prod = 0.002 / double
 
van3autokat = 0.1 / double
 
van3sat = 10 / double
 
k2van3 = 0.3 / double
 
#glvprod = 0.0 / double
 
#glvbreakdown = 0. / double
 
label = / label
 
label = <b>Integration parameters</b> / label
 
dt = 0.1 / double 
 
rd_dt = 1.0 / double
 
datadir = . / directory 
 
movie = false / bool
 
nit = 100000 / int
 
maxt = 1000. / double
 
storage_stride = 10 / int
 
xml_storage_stride = 500 / int
 
rseed = -1 / int
 
#label = / label
 
#label = <b>Parameters for new chemical</b> / label
 
#glvproduction = 0. / double
 
#glvdecay = 0. / double
 
#glvdiffthreshold = 0. / double
 
label = / label
 
label = <b>Meinhardt leaf venation model</b> / label
 
constituous_expansion_limit = 16 / int
 
vessel_inh_level = 1 / double
 
vessel_expansion_rate = 0.25 / double
 
d = 0. / double
 
e = 0. / double
 
f = 0. / double
 
c = 0. / double
 
mu = 0. / double
 
nu = 0. / double
 
rho0 = 0. / double
 
rho1 = 0. / double
 
c0 = 0. / double
 
gamma = 0. / double
 
eps = 0. / double
 
label = / label
 
label = <b>User-defined parameters</b> / label
 
k = 0., 0.0, 0.0, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. / doublelist
 
i1 = 0 / int
 
i2 = 0 / int
 
i3 = 0 / int
 
i4 = 0 / int
 
i5 = 0 / int
 
s1 =  / string
 
s2 =  / string
 
s3 =  / string
 
b1 = false / bool
 
b2 = false / bool
 
b3 = false / bool
 
b4 = false / bool
 
dir1 = . / directory
 
dir2 = . / directory
 
export_interval = 0 / int
 
export_fn_prefix = cell. / string
src/canvas.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 <time.h>
 
#include <string>
 
#include <fstream>
 
#include <streambuf>
 
#include <QGraphicsScene>
 
#include <QGraphicsView>
 
#include <qdatetime.h>
 
#include <q3mainwindow.h>
 
#include <qstatusbar.h>
 
#include <qmessagebox.h>
 
#include <qmenubar.h>
 
#include <qapplication.h>
 
#include <qpainter.h>
 
#include <qprinter.h>
 
#include <qlabel.h>
 
#include <qimage.h>
 
#include <q3progressdialog.h>
 
#include <qtimer.h>
 
#include <qslider.h>
 
#include <qpixmap.h>
 
#include <qfile.h>
 
#include <qdir.h>
 
#include <q3filedialog.h>
 
#include <QGraphicsItem>
 
#include <QList>
 
#include <QDir>
 
#include <QFileInfo>
 
#include <QDebug>
 

	
 
#include <set>
 

	
 
//Added by qt3to4:
 
#include <Q3ValueList>
 
#include <Q3PopupMenu>
 
#include <QMouseEvent>
 
#include <typeinfo>
 
#include <cstring>
 
#include <q3process.h>
 
#include <qlayout.h>
 
#include <qspinbox.h>
 
#include <fstream>
 
#include <sstream>
 
#include "pardialog.h"
 
#include "parameter.h"
 
#include "canvas.h"
 
#include "node.h"
 
#include "nodeset.h"
 
#include "nodeitem.h"
 
#include "cellitem.h"
 
#include "wallitem.h"
 
#include "mesh.h"
 
#include "xmlwrite.h"
 
#include "OptionFileDialog.h"
 
#include <cstdlib>
 
#include <cstdio>
 
#include "modelcatalogue.h"
 

	
 
#include <algorithm>
 

	
 
// Include VIB and PSB logos
 
#include "psb.xpm"
 
#include "cwi.xpm"
 

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

	
 
using namespace std;
 

	
 
// We use a global variable to save memory - all the brushes and pens in
 
// the mesh are shared.
 

	
 
#define FNAMESIZE 100
 
#define QUOTE_ME(s) QUOTE_ME_2NDLEV(s)
 
#define QUOTE_ME_2NDLEV(s) #s
 

	
 
static QColor dark_red("darkRed");
 

	
 

	
 
static const int imageRTTI = 984376;
 
extern Parameter par;
 
const QString Main::caption("Virtual leaf");
 
const QString Main::caption_with_file("Virtual leaf: %1");
 

	
 
FigureEditor::FigureEditor(
 
			   QGraphicsScene& c, Mesh &m, QWidget* parent,
 
			   const char* name, Qt::WindowFlags f) :
 
  QGraphicsView(&c,parent), mesh(m)
 
{
 
  intersection_line = 0;
 
  //angle_line = 0;
 
  setInteractive(true);
 
  moving = 0;  
 
  rotation_mode = false;
 
}
 

	
 

	
 
void FigureEditor::clear()
 
{
 
  QList<QGraphicsItem *> list = scene()->items();
 
  QList<QGraphicsItem *>::Iterator it = list.begin();
 
  for (; it != list.end(); ++it) {
 
    delete *it;
 
  }
 
}
 

	
 
void FigureEditor::wheelEvent(QWheelEvent *event)
 
{
 
  scaleView(pow((double)2, -event->delta() / 240.0));
 
}
 

	
 

	
 
void FigureEditor::scaleView (qreal scaleFactor)
 
{
 
  qreal factor = matrix().scale(scaleFactor, scaleFactor).mapRect(QRectF(0, 0, 1, 1)). width();
 
  if (factor < 0.07 || factor > 100) return;
 
  scale (scaleFactor, scaleFactor);
 
}
 

	
 
void FigureEditor::Save(const char *fname, const char *format, int sizex, int sizey)
 
{
 

	
 
  QImage *image = new QImage(sizex, sizey, QImage::Format_RGB32);
 
  image->fill(QColor(Qt::white).rgb());
 
  QPainter *painter=new QPainter(image);
 

	
 
  render(painter);
 

	
 
  image->save(QString(fname),format);
 
  delete painter;
 
  delete image;
 
}
 

	
 
void FigureEditor::mousePressEvent(QMouseEvent* e)
 
{
 
  static QList<Node*> selected;
 
  emit MousePressed();
 

	
 
  //QPointF p = matrix().inverted().map(e->pos());
 
  QPointF p = mapToScene(e->pos());
 

	
 
#ifdef QDEBUG  
 
  qDebug() << endl << "MousePressEvent location: (" << p.x() << "," << p.y() << ")." << endl;
 
  qDebug() << "Magnification:  " << Cell::Magnification() << endl;
 
  qDebug() << "Offsets:  " << Cell::Offset() << endl;
 
#endif
 

	
 

	
 
  QList<QGraphicsItem *> l=scene()->items(p);
 

	
 
#ifdef QDEBUG  
 
  qDebug() << "MousePressEvents, items: " << l.size() << endl;
 
  qDebug() << "Mouse button modifier: " << e->modifiers() << endl;
 
#endif
 

	
 
  if (e->button()==Qt::RightButton || l.size()==0) {
 

	
 
    //cerr << "Drawing an intersection line from " << p.x() << ", " << p.y() << endl;
 
    intersection_line = new QGraphicsLineItem( 0, scene() );
 
    intersection_line->setPen( QPen( QColor("red"), 3, Qt::DashLine ) );
 
    intersection_line->setLine( QLineF(p,p) );
 
    intersection_line->setZValue( 100 );
 
    intersection_line->show();
 
  }
 

	
 
  for (QList<QGraphicsItem *>::Iterator it=l.begin(); it!=l.end(); ++it) {
 
#ifdef QDEBUG
 
    qDebug() << typeid(**it).name() << endl;
 
#endif
 

	
 
    if ( !strcmp(typeid(**it).name(),"8NodeItem")) {
 

	
 
      stringstream data_strstream;
 
      data_strstream << (dynamic_cast<NodeItem*>(*it))->getNode();
 
      dynamic_cast<Main *>(parent())->UserMessage(QString(data_strstream.str().c_str()));
 

	
 
      (dynamic_cast<NodeItem*>(*it))->OnClick(e->button());
 
    }
 
    else 
 
      if ( !strcmp(typeid(**it).name(),"8CellItem") ) {
 

	
 
	Cell &c=((dynamic_cast<CellItem *>(*it))->getCell());      
 
	// OnClick to be defined in end-user code
 
	c.OnClick(e);
 
      } else {
 
	if ( !strcmp(typeid(**it).name(),"8WallItem") ) {
 
	  (dynamic_cast<WallItem *>(*it))->OnClick(e);
 
	} 
 
      }
 
  }
 

	
 
  FullRedraw();
 
  moving = 0;
 
}
 

	
 
void FigureEditor::mouseMoveEvent(QMouseEvent* e)
 
{
 

	
 
  // User chooses "rotation mode" and we can rotate the object around its center of mass
 
  if (dynamic_cast<Main *>(parent())->RotationModeP()) {
 

	
 
    QPointF p = mapToScene(e->pos());
 
    p.setX(p.x() / Cell::Magnification());
 
    p.setY(p.y() / Cell::Magnification());
 

	
 

	
 
    // get object's center of mass
 
    QPointF rotation_midpoint = mesh.Centroid()*Cell::Factor() - Cell::Offset();
 

	
 

	
 
    // calculate rotation angle
 
    double dy = (rotation_midpoint.y() - p.y());
 
    double dx = (rotation_midpoint.x() - p.x());
 
    double new_rot_angle = atan2(dx, dy);
 
    double d_alpha = new_rot_angle - rot_angle;
 
    rot_angle = new_rot_angle;
 

	
 
    mesh.Rotate(d_alpha, ( Vector(rotation_midpoint) + Cell::Offset() ) / Cell::Factor() );
 

	
 
    dynamic_cast<Main *>(parent())->Plot(0);
 
    FullRedraw();
 
    return;
 
  }
 
 
 
  
 
  if ( moving ) {
 

	
 
    QPointF p = mapToScene(e->pos());
 
    moving->userMove(p.x() - moving_start.x(),
 
		     p.y() - moving_start.y());
 
    moving_start = p;
 
    scene()->update();
 

	
 
  }
 

	
 
  //cerr << "event";
 

	
 
  // keep track of intersection line to interactively cut a growing leaf
 

	
 
  if ( intersection_line ) {
 

	
 
    QPointF sp = intersection_line -> line().p1(); // startpoint
 
    QPointF ep = mapToScene(e->pos()); // endpoint
 
    intersection_line -> setLine( QLineF(sp, ep) ); 
 
    scene()->update();
 
    // Need this for Mac
 
    FullRedraw();
 
  }
 
}
 

	
 
//void FigureEditor::contentsMouseReleaseEvent(QMouseEvent* e)
 
void FigureEditor::mouseReleaseEvent(QMouseEvent* e)
 
{
 

	
 
  emit MouseReleased();
 
  // intersection line for leaf was finished now.
 

	
 
  if (e->button()==Qt::LeftButton) { 
 
    if (intersection_line ) {
 
#ifdef QDEBUG
 
      qDebug() << "Trying to cut leaf" << endl;
 
#endif
 
      QPointF sp = intersection_line -> line().p1(); // startpoint
 
      QPointF ep = mapToScene(e->pos());
 

	
 
      intersection_line -> setLine( QLineF(sp, ep) ); 
 
      intersection_line -> show();
 

	
 
      vector <CellItem *> intersected_cells = getIntersectedCells();
 

	
 
      // no cells selected, do nothing
 
      if (intersected_cells.size()==0) {
 
#ifdef QDEBUG
 
	qDebug() << "No cells detected :-(" << endl;
 
#endif
 
	return;
 
      }
 

	
 

	
 
      Vector startpoint = Vector(sp.x(), sp.y()) / Cell::Factor() - Cell::Offset();
 
      Vector endpoint = Vector(ep.x(), ep.y()) / Cell::Factor() - Cell::Offset();
 

	
 
      NodeSet *node_set = new NodeSet;
 

	
 
      for (vector<CellItem *>::iterator it = intersected_cells.begin(); it != intersected_cells.end(); it++) {
 
	(*it)->setBrush(QBrush("purple"));
 
	
 
	Cell &c=(*it)->getCell();
 

	
 
	// sometimes the cell hasn't properly divided yet before the
 
	// next division is called?  so check for it?  let's find a way
 
	// to do this later. Note that this functionality currently
 
	// might result in a segmentation fault for users who are
 
	// quickly dragging and releasing division lines...
 
	scene()->update();
 

	
 
#ifdef QDEBUG
 
	qDebug() << "Dividing Cell " << c.Index() << endl;
 
#endif
 

	
 
	c.DivideOverGivenLine( startpoint, endpoint, true, node_set);
 
      }
 

	
 
      node_set->CleanUp();
 
      mesh.AddNodeSet(node_set);
 

	
 
#ifdef QDEBUG
 
      qDebug() << "Done DivideOverGivenLine" << endl;
 
#endif
 

	
 
      mesh.TestIllegalWalls();
 
      // Do the actual cutting and removing
 
      if (intersected_cells.size()) {
 
	mesh.CutAwayBelowLine( startpoint, endpoint ); 
 

	
 
	// Correct flags of nodeset
 
	for (NodeSet::iterator i = node_set->begin(); i != node_set->end(); i++) {
 
	  (*i)->SetSAM();
 
	  (*i)->SetBoundary();
 
	}
 

	
 
	// Make cells attached to nodeset part of the boundary
 
	// This is a temporary solution for the following:
 
	//   If a cell attached to a NodeSet divides (with a division plane intersecting the node set), 
 
	//   we must insert a new node into the node set.
 
	// For now, we add a layer of "virtual cells" inbetween. 
 
	list<Cell *> cells_attached_to_nodeset = node_set->getCells();
 
	for ( list<Cell *>::iterator c = cells_attached_to_nodeset.begin(); c != cells_attached_to_nodeset.end(); c++) {
 
	  (*c)->SetBoundary(Cell::SAM);
 
	}
 

	
 

	
 

	
 
#ifdef QDEBUG
 
	qDebug() << "Done CutAwayBelowLine" << endl;
 
#endif
 
	mesh.TestIllegalWalls();
 
	mesh.RepairBoundaryPolygon();
 
#ifdef QDEBUG
 
	qDebug() << "Done RepairBoundaryPolygon" << endl;
 
#endif
 
	mesh.TestIllegalWalls();
 
	mesh.CleanUpWalls();
 
#ifdef QDEBUG
 
	qDebug() << "Done CleanUpWalls" << endl;
 
#endif
 
	mesh.TestIllegalWalls();
 
      }
 

	
 
      dynamic_cast<Main *>(parent())->Plot();
 

	
 
#ifdef QDEBUG
 
      cerr << "NodeSet of cutting line: " << *node_set << endl;
 
#endif
 
    }
 
  } else 
 
    if (e->button()==Qt::RightButton) {
 

	
 
      if (intersection_line /* && !angle_line*/) {
 

	
 
	QPointF p = mapToScene(e->pos());
 
	QPointF sp = intersection_line->line().p1();
 

	
 
	viewport()->setMouseTracking( TRUE );
 
      } 
 
    }
 
}
 

	
 

	
 
// returns a vector of pointer to cells colliding with intersection line
 
vector <CellItem *> FigureEditor::getIntersectedCells(void)
 
{ 
 
  vector <CellItem *> colliding_cells;
 

	
 
  QList<QGraphicsItem *> l = intersection_line->collidingItems( );
 

	
 
#ifdef QDEBUG
 
  qDebug() <<  "l.size() = " << l.size() << endl;
 
#endif
 

	
 
  for (QList<QGraphicsItem *>::Iterator it=l.begin(); it!=l.end(); ++it) {
 

	
 
#ifdef QDEBUG
 
    qDebug() << typeid(**it).name() << endl;
 
#endif
 

	
 
    if ( !strcmp(typeid(**it).name(),"8CellItem") ) {
 
      colliding_cells.push_back(dynamic_cast<CellItem *>(*it));
 
    }
 
  }
 

	
 
  delete intersection_line;
 
  intersection_line = 0;
 
  return colliding_cells;
 
}
 

	
 
void FigureEditor::FullRedraw(void)
 
{
 
  QList<QRectF> rl;
 
  rl << sceneRect();
 
  updateScene(rl);   
 
}
 

	
 

	
 
NodeItem *FigureEditor::selectedNodeItem(QList<QGraphicsItem *> graphicItems) const
 
{
 
  NodeItem *nodeItem = 0;
 
  // graphicItems is a list of the QgraphicItems under the mouse click event 
 
  QList<QGraphicsItem *>::Iterator it = graphicItems.begin();
 
  for (; it != graphicItems.end(); ++it) {
 
    if ( !strcmp(typeid(**it).name(),"8NodeItem")) {
 
      // the object under the mouse click is a Nodeitem
 
      nodeItem = dynamic_cast<NodeItem *>(*it);
 
      // indicate we've selected it
 
      nodeItem->setBrush(dark_red);
 
      break;
 
    }
 
  }
 
  return nodeItem;
 
}
 

	
 

	
 
void FigureEditor::insertNode(QPointF p)
 
{
 
  Node *node = new Node(p.x(), p.y(), 0);
 
  mesh.AddNode(node);
 
  scene()->clearSelection();
 
  dynamic_cast<Main *>(parent())->Plot();
 
  FullRedraw();
 
#ifdef QDEBUG  
 
  qDebug() << "Node: " << p << endl;
 
#endif
 
}
 

	
 
static uint mainCount = 0;
 

	
 
Main::Main(QGraphicsScene& c, Mesh &m, QWidget* parent, const char* name, Qt::WindowFlags f) :
 
  Q3MainWindow(parent,name,f),
 
  MainBase(c,m),
 
  mesh(m)
 
{
 
  editor = new FigureEditor(canvas,mesh, this);
 

	
 
#ifdef QDEBUG
 
  qDebug() << "Interactive = " << editor->isEnabled();
 
#endif
 

	
 
  working_dir = 0;
 
  QObject::connect( editor, SIGNAL(MousePressed()), this, SLOT(PauseIfRunning()));
 
  QObject::connect( editor, SIGNAL(MouseReleased()), this, SLOT(ContIfRunning()));
 
  QMenuBar* menu = menuBar();
 

	
 
  Q3PopupMenu* file = new Q3PopupMenu( menu );
 

	
 
  file->insertItem("&Read leaf", this, SLOT(readStateXML()));
 
  file->insertItem("&Save leaf", this, SLOT(saveStateXML()));
 
  file->insertItem("Snapshot", this, SLOT(snapshot()), Qt::CTRL+Qt::SHIFT+Qt::Key_S);
 

	
 
  file->insertSeparator();
 
  file->insertItem("Read next leaf", this, SLOT(readNextStateXML()), Qt::Key_PageDown);
 
  file->insertItem("Read previous leaf", this, SLOT(readPrevStateXML()), Qt::Key_PageUp);
 
  file->insertItem("Read last leaf", this, SLOT(readLastStateXML()), Qt::Key_End);
 
  file->insertItem("Read first leaf", this, SLOT(readFirstStateXML()), Qt::Key_Home);
 
  file->insertItem("Export cell data", this, SLOT(exportCellData()));
 

	
 
  file->insertSeparator();
 
  file->insertItem("&Print...", this, SLOT(print()), Qt::CTRL+Qt::Key_P);
 
  file->insertSeparator();
 
  file->insertItem("E&xit", qApp, SLOT(quit()), Qt::CTRL+Qt::Key_Q);
 
  menu->insertItem("&File", file);
 

	
 
  Q3PopupMenu* edit = new Q3PopupMenu( menu );
 
  edit->insertItem("Reset Chemicals and Transporters", this, SLOT( CleanMesh()), Qt::CTRL+Qt::Key_R );
 
  edit->insertItem("Reset Chemicals", this, SLOT( CleanMeshChemicals()) );
 
  edit->insertItem("Reset Transporters", this, SLOT( CleanMeshTransporters()) );
 
  edit->insertItem("Randomize PIN1 Transporters", this, SLOT( RandomizeMesh()) );
 
  edit->insertItem("Cut away SAM", this, SLOT( CutSAM() ));
 
  menu->insertItem("&Edit", edit);
 

	
 
  run = new Q3PopupMenu( menu );
 
  running = false;
 
  paused_id = run->insertItem("&Simulation paused", this, SLOT(togglePaused()), Qt::Key_S);
 
  run->setItemChecked(paused_id, FALSE);
 

	
 
  menu->insertItem("&Run", run);
 

	
 
  view = new Q3PopupMenu( menu );
 
  view->insertItem("&Zoom in", this, SLOT(zoomIn()), Qt::CTRL+Qt::Key_Equal);
 
  view->insertItem("Zoom &out", this, SLOT(zoomOut()), Qt::CTRL+Qt::Key_Minus);
 
  view->insertSeparator();
 
  com_id = view->insertItem("Show cell &centers", this, SLOT(toggleShowCellCenters()));
 
  view->setItemChecked(com_id, FALSE);
 

	
 
  mesh_id = view->insertItem("Show &nodes", this, SLOT(toggleShowNodes()), Qt::CTRL+Qt::SHIFT+Qt::Key_N);
 
  view->setItemChecked(mesh_id, TRUE);
 
  node_number_id = view->insertItem("Show node numbers", this, SLOT(toggleNodeNumbers()), Qt::CTRL+Qt::SHIFT+Qt::Key_M);
 
  view->setItemChecked(node_number_id, FALSE);
 
  cell_number_id = view->insertItem("Show cell numbers", this, SLOT(toggleCellNumbers()));
 
  view->setItemChecked(cell_number_id, FALSE);
 
  hide_cells_id = view->insertItem("Hide cells", this, SLOT(toggleHideCells()));
 
  view->setItemChecked(hide_cells_id, FALSE);
 
  border_id = view->insertItem("Show &border cells", this, SLOT(toggleShowBorderCells()));
 
  view->setItemChecked(border_id, FALSE);
 
  cell_axes_id = view->insertItem("Show cell &axes", this, SLOT(toggleCellAxes()));
 
  cell_strain_id = view->insertItem("Show cell &strain", this, SLOT(toggleCellStrain()));
 
  view->setItemChecked(cell_axes_id, FALSE);
 
  fluxes_id = view->insertItem("Show &fluxes", this, SLOT(toggleShowFluxes()));
 
  view->setItemChecked(fluxes_id, FALSE);
 
  cell_walls_id = view->insertItem("Show transporters", this, SLOT(toggleShowWalls()));
 
  view->setItemChecked(cell_walls_id, FALSE);
 
 // apoplasts_id = view->insertItem("Show apoplasts", this, SLOT(toggleShowApoplasts()));
 
 // view->setItemChecked(apoplasts_id, FALSE);
 
  view->insertSeparator();
 
  only_boundary_id = view->insertItem("Show only leaf &boundary", this, SLOT(toggleLeafBoundary()));
 
  view->insertSeparator();
 
  movie_frames_id = view->insertItem("Start saving movie &frames", this, SLOT(toggleMovieFrames()));
 
  view->setItemChecked(movie_frames_id, par.movie);
 

	
 
  view->setItemChecked(only_boundary_id, FALSE);
 
  menu->insertItem("&View", view);
 

	
 

	
 
  options = new Q3PopupMenu( menu );
 
  dyn_cells_id = options->insertItem("Cell growth", this, SLOT(toggleDynCells()));
 
  options->setItemChecked(dyn_cells_id, true);
 

	
 
  options->insertItem("Edit &parameters", this, SLOT(EditParameters()), Qt::CTRL+Qt::Key_E);
 

	
 
  rotation_mode_id = options->insertItem("Rotate leaf", this, SLOT(EnterRotationMode()), Qt::CTRL + Qt::SHIFT + Qt::Key_R);
 
  options->setItemChecked(rotation_mode_id, false);
 

	
 
  menu->insertItem("&Options",options);
 

	
 
  // Menu of models
 
  modelmenu = new QMenu( menu );
 
  menu->insertItem("&Models", modelmenu);
 

	
 

	
 
  menu->insertSeparator();
 

	
 
  helpmenu = new Q3PopupMenu( menu );
 
  tooltips_id = helpmenu->insertItem("Show Cell&Info", this, SLOT(Refresh()));
 
  helpmenu->setItemChecked(tooltips_id, true);
 
  helpmenu->insertSeparator();
 
    //helpmenu->insertSeparator();
 
  helpmenu->insertItem("&LICENSE", this, SLOT(gpl()) );
 
  helpmenu->insertItem("About", this, SLOT(about()) ); //, Key_F1);
 
	
 
	menu->insertItem("&Help",helpmenu);
 
  statusBar();
 
  setCentralWidget(editor);
 
  printer = 0;
 
  init();
 

	
 
  // Start timer which repetitively invokes
 
  // a simulation time step
 
  timer = new QTimer( this );
 
  connect( timer, SIGNAL(timeout()), SLOT(TimeStepWrap()) );
 

	
 
  stopSimulation();
 
  statusBar()->addWidget(new QLabel("Ready."));
 
  setCaption(caption);
 
  gifanim = 0;
 

	
 
  infobar = new InfoBar();
 
  addDockWindow(infobar);
 

	
 
}
 

	
 
void Main::RefreshInfoBar(void)
 
{
 
  infobar->SetText(mesh.ModelID());
 
}
 

	
 

	
 
void Main::UserMessage(QString message, int timeout)
 
{
 
  statusBar()->showMessage(message, timeout);
 
}
 

	
 

	
 
void Main::init()
 
{
 
  clear();
 

	
 
  static int r=24;
 
  srand(++r);
 

	
 
  mainCount++;
 
}
 

	
 
Main::~Main()
 
{
 
  delete printer;
 
  if ( !--mainCount ) {
 
  }
 
  //EndGifAnim();
 
}
 

	
 
void Main::newView()
 
{
 
  // Open a new view... have it delete when closed.
 
  Main *m = new Main(canvas, mesh, 0, 0, Qt::WDestructiveClose);
 
  qApp->setMainWidget(m);
 
  m->show();
 
  qApp->setMainWidget(0);
 
}
 

	
 

	
 
void Main::EditParameters()
 
{
 

	
 
  ParameterDialog *pardial = new ParameterDialog(this, "stridediag");
 

	
 
  // Make sure the values in the parameter dialog are updated after a file is read 
 
  // each method changing the parameters (reading XML or PAR files) should
 
  // emit this signal
 
  QObject::connect(   this, SIGNAL( ParsChanged() ), pardial, SLOT( Reset() ) );
 
}
 

	
 
void Main::savePars()
 
{
 

	
 
  stopSimulation();
 

	
 
  Q3FileDialog *fd = new Q3FileDialog( this, "file dialog", TRUE );
 
  fd->setMode( Q3FileDialog::AnyFile );
 
  fd->setFilter( "Parameter files (*.par)");
 
  fd->setDir(par.datadir);
 

	
 
  QString fileName;
 
  if ( fd->exec() == QDialog::Accepted ) {
 
    fileName = fd->selectedFile();
 
    ofstream parfile((const char *)fileName);
 
    par.Write(parfile);
 
  }
 

	
 
  startSimulation();
 
}
 

	
 
void Main::readPars()
 
{
 

	
 
  stopSimulation();
 

	
 
  Q3FileDialog *fd = new Q3FileDialog( this, "file dialog", TRUE );
 
  fd->setMode( Q3FileDialog::ExistingFile );
 
  fd->setFilter( "Parameter files (*.par)");
 
  fd->setDir(par.datadir);
 

	
 
  QString fileName;
 
  if ( fd->exec() == QDialog::Accepted ) {
 
    fileName = fd->selectedFile();
 
    par.Read((const char *)fileName);
 
  }
 

	
 
  emit ParsChanged();
 
}
 

	
 

	
 
void Main::saveStateXML()
 
{
 

	
 
  stopSimulation();
 
  Q3FileDialog *fd = new Q3FileDialog( this, "file dialog", TRUE );
 
  fd->setMode( Q3FileDialog::AnyFile );
 
  fd->setFilter( "LeafML files (*.xml)");
 
  fd->setDir(par.datadir);
 
  QString fileName;
 

	
 
  if ( fd->exec() == QDialog::Accepted ) {
 
    fileName = fd->selectedFile();
 
    if ( QFile::exists( fileName ) &&
 
	 QMessageBox::question(
 
			       this,
 
			       tr("Overwrite File? -- Leaf Growth"),
 
			       tr("A file called %1 already exists."
 
				  " Do you want to overwrite it?")
 
			       .arg( fileName ),
 
			       tr("&Yes"), tr("&No"),
 
			       QString::null, 1, 1 ) ) {
 
      return saveStateXML();
 

	
 
    } else {
 

	
 
      mesh.XMLSave((const char *)fileName, XMLSettingsTree());
 

	
 
    }
 
  }
 
}
 

	
 

	
 

	
 
void Main::snapshot()
 
{
 

	
 

	
 
  stopSimulation();
 
  Q3FileDialog *fd = new Q3FileDialog( this, "Save snapshot", TRUE );
 
  fd->setDir(par.datadir);
 
  fd->setMode( Q3FileDialog::AnyFile );
 
  fd->setFilter( "Image files (*.pdf *.png *.jpg *.tif *.bpm)");
 
  QString fileName;
 

	
 
  if ( fd->exec() == QDialog::Accepted ) {
 
    fileName = fd->selectedFile();
 
    
 
    // extract extension from filename
 
    QFileInfo fi(fileName);
 
    QString extension = fi.suffix();
 
    
 
    if (extension.isEmpty()) {
 
      extension = "jpg";
 
      fileName += ".";
 
      fileName += extension;
 
    }
 
    
 
    
 
    if ( QFile::exists( fileName ) &&
 
	 QMessageBox::question(
 
			       this,
 
			       tr("Overwrite file? -- Cell data export"),
 
			       tr("A file called %1 already exists."
 
				  " Do you want to overwrite it?").arg( fileName ),
 
			       QMessageBox::Yes, QMessageBox::No 
 
			       ) == QMessageBox::No
 
	 ) {
 
      
 
      return snapshot();
 
      
 
    } else {
 

	
 
      // Save bitmaps at 1024x768
 
      Save((const char *)fileName, extension, 1024, 768);
 
      QString status_message = QString("Wrote snapshot to %1").arg(fileName);
 
      cerr << status_message.toStdString() << endl;
 
      statusBar()->showMessage(status_message);
 

	
 
    }
 
  }
 
}
 

	
 

	
 

	
 
void Main::readPrevStateXML()
 
{
 

	
 
  // if we have already read a file, read the next file
 
  if (!currentFile.isEmpty() && working_dir) {
 
    QString next_file;
 

	
 
    QStringList xml_files = working_dir->entryList("*.xml");
 
    QString currentFile_nopath = currentFile.section( '/', -1 );
 
    QString currentFile_path = currentFile.section( '/', 0, -2 );
 

	
 
    QList<QString>::iterator f = xml_files.find( currentFile_nopath );
 

	
 
    if (f == xml_files.end()) {
 
      return;
 
    }
 

	
 
    if (f==xml_files.begin()) {
 
      QMessageBox mb( "Read previous leaf",
 
		      "No more files",
 
		      QMessageBox::Information,
 
		      QMessageBox::Ok | QMessageBox::Default,
 
		      QMessageBox::NoButton,
 
		      QMessageBox::NoButton);
 
      mb.exec();
 
      return;
 
    }
 
    next_file = *(--f);
 
    next_file = currentFile_path+"/"+next_file;
 

	
 
    readStateXML((const char *)next_file);
 

	
 
  }
 
}
 

	
 
int Main::readStateXML(const char *filename, bool geometry, bool pars, bool simtime)
 
{
 

	
 
  try {
 
    xmlNode *settings;
 
    mesh.XMLRead((const char *)filename, &settings, geometry, pars, simtime);
 
#ifdef QDEBUG
 
    qDebug() << "Reading done."<< endl;
 
#endif
 
    XMLReadSettings(settings);
 
    xmlFree(settings);
 
    Cell::SetMagnification(1);
 
    Cell::setOffset(0,0);
 

	
 
    FitLeafToCanvas();
 
    
 
    currentFile =  QString(filename);
 

	
 
    Plot();
 
    QString status_message = QString("Successfully read leaf from file %1. Time is %2 h.").arg(currentFile).arg(mesh.getTimeHours().c_str());
 
    cerr << status_message.toStdString() << endl;
 
    setCaption(caption_with_file.arg(filename));
 
    statusBar()->message(status_message);
 
    emit ParsChanged();
 
#ifdef QDEBUG
 
    qDebug() << "Done. Returning 0." << endl;
 
#endif
 
    return 0;
 
  } catch (const char *error_message) {
 
    QMessageBox mb( "Read leaf from XML file",
 
		    QString(error_message),
 
		    QMessageBox::Critical,
 
		    QMessageBox::Ok | QMessageBox::Default,
 
		    Qt::NoButton,
 
		    Qt::NoButton);
 
    mb.exec();
 
    return 1;
 
  }
 
}
 

	
 

	
 
void Main::readNextStateXML()
 
{
 

	
 
  // if we have already read a file, read the next file
 
  if (!currentFile.isEmpty() && working_dir) {
 
    QString next_file;
 

	
 
    QStringList xml_files = working_dir->entryList("*.xml");
 
    QString currentFile_nopath = currentFile.section( '/', -1 );
 
    QString currentFile_path = currentFile.section( '/', 0, -2 );
 

	
 

	
 
    QList<QString>::iterator f = xml_files.find( currentFile_nopath );
 

	
 
    if (f == xml_files.end()) {
 
      return;
 
    }
 

	
 
    ++f;
 
    if (f==xml_files.end()) {
 
      QMessageBox mb( "Read next leaf",
 
		      "No more files",
 
		      QMessageBox::Information,
 
		      QMessageBox::Ok | QMessageBox::Default,
 
		      QMessageBox::NoButton,
 
		      QMessageBox::NoButton);
 
      mb.exec();
 
      return;
 
    }
 
    next_file = *f;
 
    next_file = currentFile_path+"/"+next_file;
 

	
 
    readStateXML((const char*)next_file);
 
  }
 
}
 

	
 
void Main::readLastStateXML()
 
{
 

	
 
  // if we have already read a file, read the next file
 
  if (!currentFile.isEmpty() && working_dir) {
 
    QString next_file;
 

	
 
    QStringList xml_files = working_dir->entryList("*.xml");
 
    QString currentFile_nopath = currentFile.section( '/', -1 );
 
    QString currentFile_path = currentFile.section( '/', 0, -2 );
 

	
 

	
 
    next_file = xml_files.back();
 

	
 
    next_file = currentFile_path+"/"+next_file;
 

	
 
    readStateXML((const char*)next_file);
 
  }
 
}
 

	
 

	
 
void Main::readFirstStateXML()
 
{
 

	
 
  // if we have already read a file, read the next file
 
  if (!currentFile.isEmpty() && working_dir) {
 
    QString next_file;
 

	
 
    QStringList xml_files = working_dir->entryList("*.xml");
 
    QString currentFile_nopath = currentFile.section( '/', -1 );
 
    QString currentFile_path = currentFile.section( '/', 0, -2 );
 

	
 

	
 
    next_file = xml_files.front();
 

	
 
    next_file = currentFile_path+"/"+next_file;
 

	
 
    readStateXML((const char*)next_file);
 
  }
 
}
 

	
 
QDir Main::GetLeafDir(void) {
 
   
 
  QDir LeafDir(QApplication::applicationDirPath()); 
 
  QStringList plugin_filters; // filter for plugins, i.e "*.dll", "*.dylib"
 
  
 
  
 
#if defined(Q_OS_WIN) 
 
  if (LeafDir.dirName().toLower() =="debug" 
 
      ||LeafDir.dirName().toLower() =="release") 
 
    LeafDir.cdUp(); 
 
  //plugin_filters << "*.dll";
 
#elif defined(Q_OS_MAC) 
 
  if (LeafDir.dirName() =="MacOS"){ 
 
    LeafDir.cdUp(); 
 
    LeafDir.cdUp(); 
 
    LeafDir.cdUp(); 
 
  }
 
  
 
#endif
 
  // for all OS-es. Move from "bin" directory to root application folder.
 
  if (LeafDir.dirName() == "bin") {
 
    LeafDir.cdUp();
 
  }
 

	
 
  LeafDir.cd("data/leaves");
 
  /* if (!LeafDir.cd("data/leaves")) {
 
    QString status_message = QString("No directory data/leaves");
 
    statusBar()->showMessage(status_message);
 
    
 
    return LeafDir;
 
  } 
 
  */
 
  return LeafDir;
 
      
 
}
 
 
 

	
 
void Main::readStateXML()
 
{
 

	
 
  //  extern Mesh mesh;
 

	
 
  stopSimulation();
 
#ifdef QDEBUG
 
  qDebug() << "Trying to open an OptionFileDialog" << endl;
 
#endif
 
  OptionFileDialog *fd = new OptionFileDialog( this, "read dialog", TRUE );
 
  fd->setMode( OptionFileDialog::ExistingFile );
 
  fd->setFilter( "LeafML files (*.xml)");
 
  if (working_dir) {
 
    fd->setDir(*working_dir);
 
  } else {
 
    fd->setDir(par.datadir);
 
  }
 
  QString fileName;
 
  if ( fd->exec() == QDialog::Accepted ) {
 

	
 
    fileName = fd->selectedFile();
 
    if (working_dir) {
 
      delete working_dir;
 
    }
 
    working_dir = fd->dir();
 

	
 
    if (readStateXML((const char *)fileName,fd->readGeometryP(), fd->readParametersP()) )
 
      return readStateXML(); // user can try again
 
  }
 
}
 

	
 

	
 
void Main::clear()
 
{
 
  editor->clear();
 
}
 

	
 
void Main::about()
 
{
 
  static QMessageBox* about = new QMessageBox
 
    ( "VirtualLeaf V1.0",
 
      "<h3>VirtualLeaf V1.0</h3>"
 
      "<p>"
 
      "An Open Source framework for cell-based modeling of plant tissue growth and development <br>"
 
      "(c) 2005-2008, Roeland Merks <i>et al.</i><br>"
 
      "    <a href=\"http://www.psb.vib-ugent.be\">VIB Department Plant Systems Biology</a>, "
 
      "Ghent, Belgium <br>"
 
      "(c) 2008-2010, <a href=\"http://www.cwi.nl/~merks\">Roeland Merks <i>et al.</i></a> <br>"
 
      "    <a href=\"http://www.cwi.nl\">Centrum Wiskunde & Informatica</a> and <a href=\"http://www.ncsb.nl\">Netherlands Consortium for Systems Biology</a>, Amsterdam, Netherlands <br>"
 
"<br>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.<br>"
 
      "<br>If you use this code for your projects, please cite our paper:"
 
"<br>Merks, Guravage, Inze, and Beemster. An Open Source framework for cell-based modeling of plant tissue growth and development. <i>Plant Physiology,</i> submitted.<br> . <br>"
 
      "<br>Please share your model plugins and extensions at <a href=\"http://virtualleaf.googlecode.com\">http://virtualleaf.googlecode.com</a>.", 
 
      QMessageBox::Information, 1, 0, 0, this, 0, FALSE );
 
  about->setButtonText( 1, "Dismiss" );
 
  about->show();
 
}
 

	
 

	
 
void Main::gpl()
 
{
 
  static QMessageBox* gpl = new QMessageBox ( "GPL License", "", 
 
      QMessageBox::Information, 1, 0, 0, this, 0, FALSE );
 

	
 
  QDir docDir(QApplication::applicationDirPath());
 
  docDir.cd("../doc"); // Where Linux expects gpl3.txt
 
  QString path = docDir.filePath("gpl3.txt");
 
  if (!docDir.exists("gpl3.txt")){
 
    docDir = QApplication::applicationDirPath();
 
    docDir.cd("doc"); // Where Windows expects gpl3.txt
 
    path = docDir.filePath("gpl3.txt");
 
  }
 

	
 
  // At this point path points either to the linux variant, which
 
  // exists, or the windows variant, which may exist. Testing the
 
  // ifstream object will determine whether we've found gpl3.txt.
 

	
 
  std::ifstream file(path.toStdString().c_str());
 
  std::string str;
 

	
 
  if (file) {
 
    file.seekg(0, std::ios::end);   
 
    str.reserve(file.tellg());
 
    file.seekg(0, std::ios::beg);
 

	
 
    str.assign((std::istreambuf_iterator<char>(file)),
 
	       std::istreambuf_iterator<char>());
 

	
 
    gpl->setDetailedText(QString(str.c_str()));
 
  }
 

	
 
  gpl->setText(QString( "<h3>GNU GENERAL PUBLIC LICENSE</h3>"
 
			"<p>Version 3, 29 June 2007</p>"
 
			"<p>Copyright &copy; 2007 Free Software Foundation, Inc. "
 
			"&lt;<a href=\"http://fsf.org/\">http://fsf.org/</a>&gt;</p><p>"
 
			"Everyone is permitted to copy and distribute verbatim copies "
 
			"of this license document, but changing it is not allowed.</p>"
 
			"<h2>GNU GENERAL PUBLIC LICENSE</h2>"));
 

	
 
  gpl->setButtonText( 1, "Dismiss" );
 
  gpl->show();
 
}
 

	
 
void Main::aboutQt(){
 
  QMessageBox::aboutQt( this, "Virtual Leaf" );
 
}
 

	
 
void Main::toggleShowCellCenters()
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleShowWalls()
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleShowApoplasts()
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleShowNodes()
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleNodeNumbers(void)
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleCellNumbers(void)
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleCellAxes(void)
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleCellStrain(void)
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleShowFluxes(void)
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleShowBorderCells()
 
{
 
  Plot();
 
}
 

	
 
void Main::toggleHideCells(void)
 
{
 
  Plot();
 
  editor->FullRedraw();
 
}
 

	
 
void Main::toggleMovieFrames(){}
 

	
 
void Main::toggleLeafBoundary(){}
 

	
 
void Main::toggleDynCells() {}
 

	
 
void Main::startSimulation(void)
 
{
 
  timer->start( 0 );
 
  statusBar()->message("Simulation started");
 
  running = true;
 
}
 

	
 
void Main::stopSimulation(void)
 
{
 
  timer->stop();
 
  cerr << "Stopping simulation" << endl;
 
  statusBar()->message("Simulation paused");
 
  running = false;
 
}
 

	
 
void Main::togglePaused()
 
{
 
  bool s = run->isItemChecked(paused_id);
 
  if (s) {
 
    cerr << "Calling start simulation" << endl;
 
    startSimulation();
 
  } else {
 
    cerr << "Calling stop simulation" << endl;
 
    stopSimulation();
 
  }
 
}
 

	
 
void Main::setFluxArrowSize(int size)
 
{
 
  flux_arrow_size = size/100.;
 
}
 

	
 

	
 
void Main::enlarge()
 
{
 
  canvas.setSceneRect( QRectF( 0,0, canvas.width()*4./3., canvas.height()*4./3.) );
 
}
 

	
 
void Main::shrink()
 
{
 
  canvas.setSceneRect( QRectF( 0,0, canvas.width()*3/4, canvas.height()*3/4) );
 
}
 

	
 

	
 
void Main::scale(double factor)
 
{
 
  QMatrix m = editor->matrix();
 
  m.scale(factor, factor);
 
  editor->setMatrix( m );
 
}
 

	
 
void Main::zoomIn()
 
{
 
  QMatrix m = editor->matrix();
 
  m.scale( 1.1, 1.1 );
 
  editor->setMatrix( m );
 
}
 

	
 
void Main::zoomOut()
 
{
 
  QMatrix m = editor->matrix();
 
  m.scale( 0.9, 0.9 );
 
  editor->setMatrix( m );
 
}
 

	
 

	
 
void Main::print()
 
{
 
  if ( !printer ) printer = new QPrinter;
 

	
 
  if ( printer->setup(this) ) {
 

	
 
    //    extern Mesh mesh;
 
    Vector bbll,bbur;
 
    mesh.BoundingBox(bbll,bbur);
 

	
 
#ifdef QDEBUG
 
    qDebug() << "bbll = " << bbll << endl;
 
    qDebug() << "bbur = " << bbur << endl;
 
#endif
 
    double cw = (bbur.x - bbll.x);
 
    double ch = (bbur.y - bbll.y);
 
    QPainter pp(printer);
 
    QRect vp=pp.viewport();
 
#ifdef QDEBUG
 
    qDebug() << "Paper width = " << vp.width() << " x " << vp.height() << endl;
 
#endif
 

	
 
    // Note that Cell is also translated...
 
    pp.translate(-bbur.x,-bbur.y);
 
    if (cw>ch) {
 
      pp.scale(vp.width()/(2*cw*Cell::Magnification()), vp.width()/(2*cw*Cell::Magnification()));
 
    } else {
 
      pp.scale(vp.height()/(2*ch*Cell::Magnification()), vp.height()/(2*ch*Cell::Magnification()));
 
    }
 
    canvas.render(&pp, QRectF(), QRectF(0.,0.,canvas.width(),canvas.height()));
 
  }
 
}
 

	
 

	
 
void Main::TimeStepWrap(void)
 
{
 
  static int t;
 
  stringstream fname;
 

	
 
  TimeStep();
 

	
 
  t = (int)mesh.getTime();
 

	
 
  if ((par.export_interval > 0) && !(t%par.export_interval)){
 
    fname << par.datadir << "/" << par.export_fn_prefix;
 
    fname.fill('0');
 
    fname.width(6);
 
    fname << t << ".csv";
 
    this->exportCellData(QString(fname.str().c_str()));
 
  }
 

	
 
  // check number of timesteps
 
  if (t == par.nit) {
 
    emit SimulationDone();
 
  }
 
}
 

	
 

	
 
void Main::RestartSim(void)
 
{
 

	
 
  stopSimulation();
 
  if ( QMessageBox::question(
 
			     this,
 
			     tr("Restart simulation?"),
 
			     tr("Restart simulation.\n"
 
				"Are you sure?"),
 
			     QMessageBox::Yes | QMessageBox::No, QMessageBox::NoButton ) == QMessageBox::Yes ) {
 

	
 
    cerr << "Restarting simulation" << endl;
 
    //    extern Mesh mesh;
 
    mesh.Clear();
 
    Init();
 
    Plot();
 
    editor->FullRedraw();
 
  } 
 
  //startSimulation();
 
}
 

	
 

	
 
void Main::FitCanvasToWindow(void)
 
{
 

	
 
  double scale_factor_x = (double)editor->width()/(double)canvas.width();
 
  double scale_factor_y = (double)editor->height()/(double)canvas.height();
 
  double scale_factor = scale_factor_x > scale_factor_y ? scale_factor_x : scale_factor_y;
 
  QMatrix m = editor->matrix();
 

	
 
#ifdef QDEBUG  
 
  qDebug() << "editor->width() = " << editor->width() << endl;
 
  qDebug() << "editor->height() = " << editor->height() << endl;
 

	
 
  qDebug() << "scale_factor = " << scale_factor << endl;
 
  qDebug() << "scale_factor_x = " << scale_factor_x << endl;
 
  qDebug() << "scale_factor_y = " << scale_factor_y << endl;
 
#endif
 
  m.scale( scale_factor, scale_factor );
 
  editor->setMatrix( m );
 
  editor->show();
 
}
 

	
 
void Main::PauseIfRunning(void)
 
{
 
  if (running) {
 
    timer->stop();
 
  }
 
}
 

	
 
void Main::ContIfRunning(void)
 
{
 
  if (running) {
 
    timer->start( 0 );
 
  }
 
}
 

	
 
void Main::FitLeafToCanvas(void) 
 
{
 

	
 
  Vector ll,ur;
 
  mesh.BoundingBox(ll, ur);
 

	
 
  ll*=Cell::Magnification(); ur*=Cell::Magnification();
 

	
 
  // give the leaf some space
 
  Vector border = ((ur-ll)/5.);
 

	
 
  QRectF bb( ll.x - border.x, ll.y - border.y, ur.x-ll.x + 2*border.x, ur.y-ll.y + 2*border.y );
 

	
 

	
 
  // cerr << ur << ", " << ll << endl;
 
  // editor->fitInView(bb, Qt::KeepAspectRatio);
 
  editor->ensureVisible(bb);
 
  //editor->setTransform(viewport);
 
}
 

	
 
void Main::CleanMesh(void) 
 
{
 
  vector<double> clean_chem(Cell::NChem());
 
  vector<double> clean_transporters(Cell::NChem());
 

	
 
  for (int i=0;i<Cell::NChem();i++) {
 
    clean_chem[i]=par.initval[i];
 
    clean_transporters[i]=0.;		
 
  }
 

	
 
  mesh.CleanChemicals(clean_chem);
 
  mesh.CleanTransporters(clean_transporters);
 

	
 
  mesh.setTime(0);
 
  Plot();
 

	
 
  editor->FullRedraw();
 
}
 

	
 
void Main::CleanMeshChemicals(void) 
 
{
 

	
 
  vector<double> clean_chem(Cell::NChem());
 

	
 
  for (int i=0;i<Cell::NChem();i++) {
 
    clean_chem[i]=par.initval[i];
 
  }
 

	
 
  mesh.CleanChemicals(clean_chem);
 
  mesh.setTime(0);
 
  Plot();
 

	
 
  editor->FullRedraw();
 
}
 

	
 
void Main::CleanMeshTransporters(void) 
 
{
 
  vector<double> clean_transporters(Cell::NChem());
 
  for (int i=0;i<Cell::NChem();i++) {
 
    clean_transporters[i]=0.;
 
  }
 

	
 
  mesh.CleanTransporters(clean_transporters);
 

	
 
  mesh.setTime(0);
 
  Plot();
 

	
 
  editor->FullRedraw();
 
}
 

	
 
void Main::RandomizeMesh(void) 
 
{
 

	
 
  vector<double> max_chem(Cell::NChem());
 
  vector<double> max_transporters(Cell::NChem());
 

	
 
  for (int i=0;i<Cell::NChem();i++) {
 
    max_transporters[i]=0.;
 
    max_chem[i]=par.initval[i];
 
  }
 

	
 
  // Amount of PIN1 at walls
 
  max_transporters[1] = 0.01;
 

	
 
  mesh.RandomizeChemicals(max_chem, max_transporters);
 

	
 
  Plot();
 
}
 

	
 
void Main::XMLReadSettings(xmlNode *settings) 
 
{
 

	
 
  MainBase::XMLReadSettings(settings);
 
  
 
  view->setItemChecked(com_id, showcentersp);
 
  view->setItemChecked(mesh_id, showmeshp);
 
  view->setItemChecked(border_id, showbordercellp);
 
  view->setItemChecked(node_number_id, shownodenumbersp);
 
  view->setItemChecked(cell_number_id, showcellnumbersp);
 
  view->setItemChecked(cell_axes_id, showcellsaxesp);
 
  view->setItemChecked(cell_strain_id, showcellstrainp);
 
  view->setItemChecked(movie_frames_id, movieframesp);
 
  view->setItemChecked(only_boundary_id, showboundaryonlyp);
 
  view->setItemChecked(fluxes_id, showfluxesp);
 
  view->setItemChecked(hide_cells_id, hidecellsp);
 
  options->setItemChecked(dyn_cells_id, dynamicscellsp);
 
  view->setItemChecked( cell_walls_id, showwallsp);
 
 // view->setItemChecked( apoplasts_id, showapoplastsp);
 
  
 
  editor->setTransform(viewport);
 
}
 

	
 
xmlNode *Main::XMLSettingsTree(void) 
 
{
 

	
 
  showcentersp = view->isItemChecked(com_id);
 
  showmeshp = view->isItemChecked(mesh_id);
 
  showbordercellp =  view->isItemChecked(border_id);
 
  shownodenumbersp =  view->isItemChecked(node_number_id);
 
  showcellnumbersp =  view->isItemChecked(cell_number_id);
 
  showcellsaxesp = view->isItemChecked( cell_axes_id );
 
  showcellstrainp = view->isItemChecked( cell_strain_id );
 
  movieframesp = view->isItemChecked(movie_frames_id);;
 
  showboundaryonlyp =  view->isItemChecked(only_boundary_id);
 
  showfluxesp = view->isItemChecked(fluxes_id);
 
  dynamicscellsp = options->isItemChecked(dyn_cells_id);
 
  showwallsp = view->isItemChecked( cell_walls_id);
 
  //showapoplastsp = view->isItemChecked( apoplasts_id);
 
  hidecellsp = view->isItemChecked( hide_cells_id);
 

	
 
  xmlNode *settings = MainBase::XMLSettingsTree();
 
  QTransform viewport(editor->transform());
 
  xmlAddChild(settings, XMLViewportTree(viewport));
 
  return settings;
 
}
 

	
 

	
 
#define QDEBUG
 

	
 
void Main::exportCellData(QString fileName) {
 
  
 
#ifdef QDEBUG  
 
  qDebug() << "exportCellData fileName: " << fileName << endl;
 
#endif
 

	
 
  QFile file(fileName);
 
  if ( file.open( IO_WriteOnly ) ) {
 
    QTextStream stream( &file );
 
    mesh.CSVExportCellData(stream);
 
    mesh.CSVExportWallData(stream);
 
    mesh.CSVExportMeshData(stream);
 
    file.close();
 
  }
 
}
 

	
 

	
 
void Main::exportCellData() {
 
  QString fileName;
 
  Q3FileDialog *fd = new Q3FileDialog( this, "file dialog", TRUE );
 

	
 
  fd->setDir(par.datadir); 
 
  stopSimulation();
 
  fd->setMode( Q3FileDialog::AnyFile );
 
  if ( fd->exec() == QDialog::Accepted ) {
 
    fileName = fd->selectedFile();
 

	
 
    // extract extension from filename
 
    QFileInfo fi(fileName);
 
    QString extension = fi.suffix();
 
    
 
    if (extension.isEmpty()) {
 
      extension = "csv";
 
      fileName += ".";
 
      fileName += extension;
 
    }
 

	
 
    if (extension!="csv" && extension!="CSV") {
 
      
 
      if (
 
	  QMessageBox::question(
 
				 this,
 
				 tr("Change extension? -- Cell data export"),
 
				 tr("VirtualLeaf can only export data to CSV format (Excel-compatible)."
 
				    "Do you want to change the file extension to \".csv\"?"),
 
				 QMessageBox::Yes, QMessageBox::No ) == QMessageBox::No 
 
	  ) {
 

	
 
	return exportCellData();
 
      } else {
 
	fileName.replace(extension,"csv");
 
	extension="csv";
 
      }
 
    }
 
    if ( QFile::exists( fileName ) &&
 
	 QMessageBox::question(
 
			       this,
 
			       tr("Overwrite file? -- Cell data export"),
 
			       tr("A file called %1 already exists."
 
				  " Do you want to overwrite it?").arg( fileName ),
 
			       QMessageBox::Yes, QMessageBox::No 
 
			       ) == QMessageBox::No
 
	 ) {
 
      return exportCellData();
 
    } else {
 
      exportCellData(fileName);
 
      QString status_message = QString("Wrote data file to %1").arg(fileName);
 
      statusBar()->showMessage(status_message);
 
    }
 
  }
 
}
 

	
 

	
 
  /* finis */
src/canvas.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.
 
 *
 
 */
 

	
 
#ifndef _CANVAS_H_
 
#define _CANVAS_H_
 

	
 
#include <q3popupmenu.h>
 
#include <q3mainwindow.h>
 
#include <q3intdict.h>
 
#include <QGraphicsScene>
 
#include <QGraphicsView>
 
#include <QList>
 
#include <QDir>
 

	
 
#include <string>
 
#include <sstream>
 

	
 
//Added by qt3to4:
 
#include <QMouseEvent>
 
#include <vector>
 
#include "simitembase.h"
 
#include "mainbase.h"
 
#include "cellitem.h"
 
#include "infobar.h"
 

	
 
#include "nodeitem.h"
 
#include "cell.h"
 

	
 
#ifdef HAVE_QWT
 
#include "data_plot.h"
 
#endif
 

	
 
#include <libxml/xpath.h>
 
#include <libxml/xmlreader.h>
 

	
 
#if defined(Q_OS_MAC)
 
#define PREFIX "cmd"
 
#else
 
#define PREFIX "crtl"
 
#endif
 

	
 

	
 
class QFile;
 
class QDir;
 
class ModelCatalogue;
 
class InfoBar;
 

	
 
class FigureEditor : public QGraphicsView {
 
  Q_OBJECT
 

	
 
    friend class Main;
 
 public:
 
  FigureEditor(QGraphicsScene&, Mesh&, QWidget* parent=0, const char* name=0, Qt::WFlags f=0);
 
  void clear();
 
  void Save(const char *fname, const char *format, int sizex=640, int sizey=480);
 
  void FullRedraw(void);
 

	
 
 protected:
 
  void mousePressEvent(QMouseEvent*);
 
  void mouseMoveEvent(QMouseEvent*);
 
  void mouseReleaseEvent(QMouseEvent*);
 
  void wheelEvent(QWheelEvent *event);
 
  void scaleView(qreal scaleFactor);
 

	
 
  vector <CellItem *> getIntersectedCells(void);
 
  void insertNode(QPointF p);
 
  NodeItem *selectedNodeItem(QList<QGraphicsItem *> l) const;  
 

	
 
 signals:
 
  void status(const QString&);
 
  void MousePressed(void);
 
  void MouseReleased(void);
 
 protected:
 
  Mesh &mesh;
 

	
 
 private:
 
  //NodeItem* moving;
 
  SimItemBase *moving;
 
  QPointF moving_start;
 

	
 
  QGraphicsLineItem *intersection_line;
 
  bool rotation_mode;
 
  bool insert_mode;
 
  double rot_angle;
 
};
 

	
 
class Main : public Q3MainWindow, public MainBase {
 
  Q_OBJECT
 
    friend class ModelCatalogue; // needs to populate menu item models
 
 public:
 
  Main(QGraphicsScene&, Mesh&, QWidget* parent=0, const char* name=0, Qt::WFlags f=0);
 
  ~Main();
 
  virtual bool ShowCentersP(void) {return view->isItemChecked(com_id);}
 
  virtual bool ShowMeshP(void) {return view->isItemChecked(mesh_id);}
 
  virtual bool ShowBorderCellsP(void) {return view->isItemChecked(border_id);}
 
  virtual bool PausedP(void) {return run->isItemChecked(paused_id);}
 
  virtual bool ShowNodeNumbersP(void) {return view->isItemChecked(node_number_id);}
 
  virtual bool ShowCellNumbersP(void) {return view->isItemChecked(cell_number_id);}
 
  virtual bool ShowCellAxesP(void) {return view->isItemChecked(cell_axes_id);}
 
  virtual bool ShowCellStrainP(void) {return view->isItemChecked(cell_strain_id);}
 
  virtual bool MovieFramesP(void) {return view->isItemChecked(movie_frames_id);}
 
  virtual bool ShowBoundaryOnlyP(void) {return view->isItemChecked(only_boundary_id);}
 
  virtual bool ShowWallsP(void) {return view->isItemChecked(cell_walls_id);}
 
 // virtual bool ShowApoplastsP(void) { return view->isItemChecked(apoplasts_id);}
 
  virtual bool ShowFluxesP(void) { return view->isItemChecked(fluxes_id); }
 
  virtual bool DynamicCellsP(void) { return options->isItemChecked(dyn_cells_id); }
 
  virtual bool RotationModeP(void) { return options->isItemChecked(rotation_mode_id); }
 
  virtual bool InsertModeP(void) { return options->isItemChecked(insert_mode_id); }
 
  virtual bool ShowToolTipsP(void) { return helpmenu->isItemChecked(tooltips_id); }
 
  virtual bool HideCellsP(void) { return view->isItemChecked(hide_cells_id); }
 
  void scale(double factor); 
 
  virtual double getFluxArrowsize(void)
 
  {
 
    return flux_arrow_size;
 
  }
 

	
 
  void FitCanvasToWindow();
 
  void FitLeafToCanvas(void);
 

	
 
  public slots:
 

	
 
  void about();
 
  void gpl();
 
  void TimeStepWrap();
 
  void togglePaused();
 
  void setFluxArrowSize(int size);
 
  void RestartSim(void);
 
  void toggleShowCellCenters(void);
 
  void toggleShowNodes(void);
 
  void toggleShowBorderCells(void);
 
  void toggleShowFluxes(void);
 
  void toggleNodeNumbers(void);
 
  void toggleCellNumbers(void);
 
  void toggleCellAxes(void);
 
  void toggleCellStrain(void);
 
  void toggleShowWalls(void);
 
  void toggleShowApoplasts(void);
 
  void toggleDynCells(void);
 
  void toggleMovieFrames(void);
 
  void toggleLeafBoundary(void);
 
  void toggleHideCells(void);
 
  void print();
 
  void startSimulation(void);
 
  void stopSimulation(void);
 
  void RefreshInfoBar(void);
 

	
 
  int readStateXML(const char *filename, bool geometry = true, bool pars=true, bool simtime = true);
 
  
 
  void EnterRotationMode(void)
 
  {
 

	
 
    UserMessage("Rotation mode. Click mouse to exit.");
 
    if (editor) {
 
      editor->rot_angle = 0. ; 
 

	
 
      // Exit rotation mode if mouse is clicked
 
      connect(editor, SIGNAL(MousePressed()), this, SLOT(ExitRotationMode()));
 
      
 
      editor->setMouseTracking(true);
 

	
 
    }
 

	
 
  }
 
  void ExitRotationMode(void)
 
  { 
 
    UserMessage("Exited rotation mode.",2000);
 

	
 
    options->setItemChecked(rotation_mode_id, false); 
 
    if (editor)
 
      disconnect(editor, SIGNAL(MousePressed()), this, SLOT(ExitRotationMode()));
 
    editor->setMouseTracking(false);
 
 
 
  }
 

	
 
  virtual void UserMessage(QString message, int timeout=0);
 
  void Refresh(void) { Plot(); }
 
  void PauseIfRunning(void);
 
  void ContIfRunning(void);
 
  virtual void XMLReadSettings(xmlNode *settings);
 

	
 
  private slots:
 
  void aboutQt();
 
  void newView();
 
  void EditParameters();
 
  QDir GetLeafDir(void);
 
  void readStateXML();
 
  int readStateXML(const char *filename, bool geometry = true, bool pars=true, bool simtime = true);
 
  void readNextStateXML();
 
  void readPrevStateXML();
 
  void readFirstStateXML();
 
  void readLastStateXML();
 
  void exportCellData();
 
  void exportCellData(QString);
 
  void saveStateXML();
 
  void snapshot();
 
  void savePars();
 
  void readPars();
 
  void clear();
 
  void init();
 
  virtual void CutSAM() { MainBase::CutSAM(); Refresh();}
 

	
 
  void enlarge();
 
  void shrink();
 
  void zoomIn();
 
  void zoomOut();
 

	
 
  void CleanMesh();
 
  void CleanMeshChemicals(void);
 
  void CleanMeshTransporters(void);
 

	
 
  void RandomizeMesh();
 

	
 
 signals:
 
  void SimulationDone(void);
 
  void ParsChanged(void);
 

	
 
 protected:
 
  Mesh &mesh;
 

	
 
 private:
 
  NodeSet *node_set;
 
  FigureEditor *editor;
 
  Q3PopupMenu* options;
 
  Q3PopupMenu *view;
 
  Q3PopupMenu *run;
 
  QMenu *modelmenu;
 
  Q3PopupMenu *helpmenu;
 

	
 
  QPrinter* printer;
 
  const QDir *working_dir;
 
  QString currentFile;
 
  //  toggle item states 
 
  int dbf_id; // options->Double Buffer
 
  int com_id; // view->Show centers
 
  int mesh_id; // view->Show mesh
 
  int node_number_id; // view->Show Node numbers
 
  int cell_number_id; // view->Show Cell numbers
 
  int border_id; // view->Show border cells
 
  int paused_id; // run->Simulation paused
 
  int cell_axes_id; // view->Show cell axes
 
  int cell_strain_id; // view->Show cell strain
 
  int only_boundary_id; // view ->Show only leaf boundary
 
  int cell_walls_id; // view -> Show transporters
 
  //int apoplasts_id; // view -> Show apoplasts
 
  int tooltips_id; // help -> Show Cell Info
 
  int hide_cells_id; // view->Hide Cells
 
  double flux_arrow_size;
 
  int movie_frames_id;
 
  int fluxes_id;
 
  int dyn_cells_id;
 
  int rotation_mode_id;
 
  int insert_mode_id;
 
  QTimer *timer;
 
  QFile *gifanim;
 
  bool running;
 
  virtual xmlNode *XMLSettingsTree(void);
 
  static const QString caption;
 
  static const QString caption_with_file;
 
  InfoBar *infobar;
 
};
 

	
 
#endif
 

	
 
/* finis*/
src/mesh.cpp
Show inline comments
 
@@ -1300,782 +1300,818 @@ void Mesh::RepairBoundaryPolygon(void) {
 
  // So, this function just puts boundary nodes into the boundary
 
  // polygon in the right order; it cannot detect boundaries from
 
  // scratch.
 

	
 
  Node *boundary_node=0, *next_boundary_node=0, *internal_node;
 
  set<int> original_boundary_nodes, repaired_boundary_nodes;
 
  vector<int> difference; // set difference result
 

	
 
  // Step 0: print out current boundary polygon
 
#ifdef QDEBUG
 
  qDebug() << endl << "Original Boundary Polygon node indices: ";
 
  foreach (Node* node, boundary_polygon->nodes) {
 
    qDebug() << node->Index() << " " ;
 
  }    
 
  qDebug() << endl << endl;
 
#endif
 

	
 
  // Step 1a: Create a set containing the current boundary polygon nodes' Indices.
 
  foreach (Node* node, boundary_polygon->nodes) {
 
    original_boundary_nodes.insert(node->Index());
 
  }
 

	
 
  // Step 1b: remove all nodes from boundary polygon
 
  boundary_polygon->nodes.clear();
 

	
 
  // Step 2: Remove all references to the boundary polygon from the Mesh's current list of nodes
 
  foreach (Node* node, nodes) {
 
    node->Unmark(); // remove marks, we need them to determine if we have closed the circle
 
    list<Neighbor>::iterator boundary_ref_pos;
 
    if ((boundary_ref_pos = find_if (node->owners.begin(), node->owners.end(), 
 
				     bind2nd(mem_fun_ref(&Neighbor::CellEquals), -1))) != node->owners.end()) {
 
      // i.e. if one of the node's owners is the boundary polygon 
 
      node->owners.erase(boundary_ref_pos); // remove the reference
 
    }
 
  }
 

	
 
  // Step 3: Search for the first boundary node.  We reconstruct the
 
  // boundary polygon by moving along the boundary nodes until we've
 
  // encircled the polygon. Since manually adding nodes may have
 
  // turned nodes previously along the boundary into internal nodes,
 
  // we search through all the node until we find first boundary node
 
  // and proceed from there. If findNextBoundaryNode() returns a node
 
  // other than the one passed to it, the original node is the first
 
  // boundary node.
 
  foreach (Node* node, nodes) {
 
    if ((findNextBoundaryNode(node))->index != node->index){
 
      next_boundary_node = node;
 
      break;
 
    }
 
  }
 

	
 
  // We have a problem if we arrive here without having found a boundary node.
 
  if (!next_boundary_node) throw("Cannot find a boundary node!.");
 

	
 
  // Reconstruct the list of boundary polygon nodes.
 
  do {
 
    boundary_node = next_boundary_node;
 
    boundary_node->Mark();
 
    boundary_polygon->nodes.push_back(boundary_node);
 
    next_boundary_node = findNextBoundaryNode(boundary_node);
 
  } while ( !next_boundary_node->Marked() );
 

	
 

	
 
  // Create a set containing the reconstructed boundary polygon nodes' Indices.
 
  for (list<Node *>::iterator it = boundary_polygon->nodes.begin(); it!=boundary_polygon->nodes.end(); ++it) {
 
    repaired_boundary_nodes.insert((*it)->Index());
 
  }
 

	
 
  // Calculate the difference between the original and repaired sets of boundary nodes
 
  // yielding the set of nodes that are no longer part of the boundary polygon.
 
  set_difference(original_boundary_nodes.begin(), original_boundary_nodes.end(),
 
                 repaired_boundary_nodes.begin(), repaired_boundary_nodes.end(), back_inserter(difference));
 

	
 
  // Tell each node in the difference that it's no longer part of the boundary polygon
 
  vector<Node *>::iterator internal_node_it;
 
  foreach (int i, difference){
 
    internal_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), i));
 
    internal_node = *internal_node_it; // dereference the itterator to get to the node pointer
 
    if (!internal_node) throw("Found a null Node pointer.");
 
    internal_node->UnsetBoundary();
 
  }
 

	
 
  boundary_polygon->ConstructConnections();
 
  for (list<Wall *>::iterator w=boundary_polygon->walls.begin(); w!=boundary_polygon->walls.end(); w++) {
 
    if ((*w)->DeadP()) {
 
      (*w)=0;
 
    }
 
  }
 
  boundary_polygon->walls.remove(0);
 
  boundary_polygon->ConstructNeighborList();
 

	
 
#ifdef QDEBUG
 
  qDebug() << "Repaired Boundary Polygon node indices: ";
 
  foreach (Node* node, boundary_polygon->nodes){
 
    qDebug() << node->Index() << " " ;
 
  }
 
  qDebug() << endl ;
 

	
 
#ifdef _undefined_
 
  qDebug() << "NODES:" << endl;
 
  foreach(Node* node, nodes) {
 
    qDebug() << *node;
 
  }
 
  qDebug() << endl;
 

	
 
  qDebug() << "WALLS:" << endl;
 
  foreach(Wall* wall, walls) {
 
    qDebug() << *wall;
 
  }
 
  qDebug() << endl;
 

	
 
  qDebug() << "CELLS:" << endl;
 
  foreach(Cell* cell, cells) {
 
    qDebug() << *cell;
 
  }
 
  qDebug() << endl;
 
#endif
 
#endif
 
}
 

	
 

	
 
Node* Mesh::findNextBoundaryNode(Node* boundary_node) {
 
  bool found_next_boundary_node = false;
 
  Node *next_boundary_node = 0;
 
  set<int> boundary_node_owners; // This is a list of the current boundary node's owners' Ids
 
  vector<int> neighborIds; // A list of the current boundary node's owners' 2nd neighbor Ids
 
  vector<set<int> *>  nodeOwners; // A vector of set pointers where each set contains the owner Ids of the nodes in the neighborIds list.
 
  vector<int> intersection; // set intersection result
 

	
 
  // The next boundary node is that which has only one owner in common with the current boundary node
 
  for (list<Neighbor>::iterator it=boundary_node->owners.begin(); it!=boundary_node->owners.end(); ++it) {
 
    if (it->cell->Index() != -1) boundary_node_owners.insert(it->cell->Index()); // Save each of the current boundary node's owners' Ids - except the boundary polygon 
 
    set<int> *owners = new set<int>; // create a set to hold a 2nd neighbor's owners' Ids
 
    nodeOwners.push_back(owners);
 
    neighborIds.push_back(it->nb2->Index());
 
    foreach(Neighbor neighbor, it->nb2->owners){
 
      if (neighbor.cell->Index() != -1) owners->insert(neighbor.cell->Index()); // Save second neighbors' owners' Ids - except the boundary polygon 
 
    }
 
  }
 
  vector<int>::iterator itt = neighborIds.begin();
 
  vector<set<int> *>::iterator it = nodeOwners.begin();
 

	
 
#ifdef QDEBUG  
 
  qDebug() << "Boundary node: " <<  boundary_node->Index() << " is owned by the following cells: ";
 
  foreach (int i, boundary_node_owners){
 
    qDebug() << i << "  ";
 
  }
 
  qDebug() << endl;
 
#endif
 

	
 
  for (; it < nodeOwners.end(); it++, itt++) {
 
    intersection.clear();
 
    set_intersection(boundary_node_owners.begin(), boundary_node_owners.end(), (*it)->begin(), (*it)->end(), back_inserter(intersection));
 

	
 
#ifdef QDEBUG  
 
    qDebug() << "The intersection of the boundary node(" << boundary_node->Index() << ") owners and its 2nd neighbor(" <<  *itt << ") owners is: ";
 
    foreach (int i, intersection){
 
      qDebug() << i << "  ";
 
    }
 
    qDebug() << endl;
 
#endif
 

	
 
    if (intersection.size() == 1){
 
      found_next_boundary_node = true;
 
      vector<Node *>::iterator next_boundary_node_it = find_if (nodes.begin(), nodes.end(), bind2nd(mem_fun(&Node::IndexEquals), *itt));
 
      next_boundary_node = *next_boundary_node_it; // defeference the itterator to get to the node pointer
 

	
 
#ifdef QDEBUG  
 
      qDebug() << "The Current boundary node is: " << boundary_node->Index()
 
	       << ". The Next boundary node is: " << *itt << ((next_boundary_node->Marked()) ? " Marked" : " Unmarked") << endl << endl;
 
#endif
 

	
 
      break;
 
    }
 
  }
 

	
 
#ifdef QDEBUG  
 
  if (!found_next_boundary_node) {
 
    qDebug() << "OOPS! Didn't find the next boundrary node!" << endl;
 
  }
 
#endif
 

	
 
  return next_boundary_node;
 
}
 

	
 

	
 
void Mesh::CleanUpWalls(void) {
 
  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 
    if ((*w)->DeadP()) {
 
      delete *w;
 
      (*w)=0;      
 
    }
 
  }
 
  walls.remove(0);
 
}
 

	
 
void Mesh::Rotate(double angle, Vector center) {
 

	
 
  // Rotate the mesh over the angle "angle", relative to center point "center".
 

	
 
  Matrix rotmat;
 

	
 
  rotmat.Rot2D(angle);
 

	
 
  for (vector<Node *>::iterator n=nodes.begin(); n!=nodes.end(); n++) {
 
    (*n)->setPos ( rotmat * ( *(*n) - center ) + center );  
 
  }
 
}
 

	
 

	
 
void Mesh::PrintWallList( void ) {
 

	
 
  transform ( walls.begin(), walls.end(), ostream_iterator<Wall>(cerr, "\n"), deref_ptr<Wall> );
 
}
 

	
 
#include <QString>
 
//#include "forwardeuler.h"
 
#include "rungekutta.h"
 

	
 
class SolveMesh : public RungeKutta {
 

	
 
private:
 
  SolveMesh(void);
 

	
 
public:
 
  SolveMesh(Mesh *m_) {
 

	
 
    m = m_;
 

	
 
    kmax=0;
 
    kount=0;
 
    xp=0; yp=0; dxsav=0;
 

	
 

	
 
  }
 

	
 
protected:
 
  virtual void derivs(double x, double *y, double *dydx) {
 

	
 
    // set mesh with new values given by ODESolver
 
    // (we must do this, because only mesh knows the connections
 
    // between the variables)
 

	
 
    m->setValues(x,y);
 
    m->Derivatives(dydx);
 

	
 
    //cerr << "Calculated derivatives at " << x << "\n";    
 
  }
 

	
 
private:
 
  Mesh *m;
 
  int kmax,kount;
 
  double *xp,**yp,dxsav;
 
  bool monitor_window;
 
};
 

	
 

	
 

	
 
void Mesh::ReactDiffuse(double delta_t) {
 

	
 
  // Set Lengths of Walls
 
  for_each ( walls.begin(), walls.end(), 
 
	     mem_fun( &Wall::SetLength ) );
 

	
 
  static SolveMesh *solver = new SolveMesh(this);
 

	
 
  int nok, nbad, nvar;
 
  double *ystart = getValues(&nvar);
 

	
 
  solver->odeint(ystart, nvar, getTime(), getTime() + delta_t, 
 
		 par.ode_accuracy, par.dt, 1e-10, &nok, &nbad);
 

	
 
  setTime(getTime()+delta_t);
 
  setValues(getTime(),ystart);
 
}
 

	
 

	
 
Vector Mesh::FirstConcMoment(int chem) {
 

	
 
  Vector moment;
 
  for (vector<Cell *>::const_iterator c=cells.begin(); c!=cells.end(); c++) {
 
    moment += (*c)->Chemical(chem) * (*c)->Centroid();
 
  }
 
  return moment / (double)cells.size();
 
}
 

	
 
/*! This member function deletes all walls connected to two dead cells from the mesh.
 
  It should be called before the Cells are actually removed.
 
  If the cell is connect to one dead cell only, that reference is substituted for a reference 
 
  to the boundary polygon.
 
*/
 
void Mesh::DeleteLooseWalls(void) {
 

	
 
  list<Wall *>::iterator w=walls.begin();
 

	
 
  while (w!=walls.end()) {
 

	
 
    // if both cells of the wall are dead, remove the wall
 
    if ((*w)->C1()->DeadP() || (*w)->C2()->DeadP()) {
 
      if ((*w)->C1()->DeadP() && (*w)->C2()->DeadP()) {
 
	delete *w;
 
	w=walls.erase(w);
 
      } else {
 
	if ((*w)->C1()->DeadP())
 
	  (*w)->c1 = boundary_polygon;
 
	else
 
	  (*w)->c2 = boundary_polygon;
 
	w++;
 
      }
 
    } else {
 
      w++;
 
    }
 

	
 
  }
 
}
 

	
 
/*void Mesh::FitLeafToCanvas(double width, double height) {
 

	
 
  Vector bbll,bbur;
 
  BoundingBox(bbll,bbur);
 

	
 
  double scale_x = width/(bbur.x-bbll.x);
 
  double scale_y = height/(bbur.y-bbll.y);
 

	
 
  double factor = scale_x<scale_y ? scale_x:scale_y;
 

	
 
  Cell::SetMagnification(factor); // smallest of scale_x and scale_y
 

	
 
  double offset_x = (width/Cell::Magnification()-(bbur.x-bbll.x))/2.;  
 
  double offset_y = (height/Cell::Magnification()-(bbur.y-bbll.y))/2.;
 

	
 
  Cell::setOffset(offset_x, offset_y);
 

	
 
  }*/
 

	
 

	
 

	
 
void Mesh::CleanChemicals(const vector<double> &clean_chem) {
 

	
 
  if (clean_chem.size()!=(unsigned)Cell::NChem()) {
 
    throw "Run time error in Mesh::CleanChemicals: size of clean_chem should be equal to Cell::NChem()";
 
  }
 
  for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*c)->SetChemical(i,clean_chem[i]);
 
    }
 
    (*c)->SetNewChemToChem();
 
  }
 
}
 

	
 

	
 
void Mesh::CleanTransporters(const vector<double> &clean_transporters) {
 

	
 
  if (clean_transporters.size()!=(unsigned)Cell::NChem()) {
 
    throw "Run time error in Mesh::CleanTransporters: size ofclean_transporters should be equal to Cell::NChem()";
 
  }
 

	
 

	
 
  // clean transporters
 
  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*w)->setTransporters1(i,clean_transporters[i]); (*w)->setNewTransporters1(i,clean_transporters[i]);
 
      (*w)->setTransporters2(i,clean_transporters[i]); (*w)->setNewTransporters2(i,clean_transporters[i]);
 
    }
 
  }
 
}
 

	
 

	
 
void Mesh::RandomizeChemicals(const vector<double> &max_chem, const vector<double> &max_transporters) {
 

	
 
  if (max_chem.size()!=(unsigned)Cell::NChem() || max_transporters.size()!=(unsigned)Cell::NChem()) {
 
    throw "Run time error in Mesh::CleanChemicals: size of max_chem and max_transporters should be equal to Cell::NChem()";
 
  }
 

	
 
  for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*c)->SetChemical(i,max_chem[i]*RANDOM());
 
    }
 
    (*c)->SetNewChemToChem();
 
  }
 

	
 
  // randomize transporters
 
  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 
    for (int i=0;i<Cell::NChem();i++) {
 
      (*w)->setTransporters1(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters1(i, (*w)->Transporters1(i) );
 
      (*w)->setTransporters2(i,max_transporters[i] * RANDOM()); (*w)->setNewTransporters2(i, (*w)->Transporters1(i) );
 
    }
 
  }
 
}
 

	
 
//!\brief Calculates a vector with derivatives of all variables, which
 
// we can pass to an ODESolver. 
 
void Mesh::Derivatives(double *derivs) {
 

	
 
  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;
 

	
 
  //static double *derivs = 0; 
 
  // derivs is allocated by RungeKutta class.
 

	
 
  for (int i=0;i<neqs;i++) {
 
    derivs[i]=0.;
 
  }
 

	
 
  // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
 

	
 
  int i=0;
 

	
 
  for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
 
    plugin->CellDynamics(*c, &(derivs[i]));
 
    i+=nchems;
 
  }
 

	
 
  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 
    // (*wr)(*w, &(derivs[i]), &(derivs[i+nchems]));
 
    plugin->WallDynamics(*w,  &(derivs[i]), &(derivs[i+nchems]));
 
    // Transport function adds to derivatives of cell chemicals
 
    double *dchem_c1 = &(derivs[(*w)->c1->Index() * nchems]);
 
    double *dchem_c2 = &(derivs[(*w)->c2->Index() * nchems]);
 
    //plugin->CelltoCellTransport(*w, &(derivs[(*w)->c1->Index() * nchems]),
 
    //	  &(derivs[(*w)->c2->Index() * nchems]));
 
    // quick fix: dummy values to prevent end user from writing into outer space and causing a crash :-)
 
    // start here if you want to implement chemical input/output into environment over boundaries
 
    double dummy1, dummy2;
 
    if ((*w)->c1->Index()<0) { // tests if c1 is the boundary pol
 
      dchem_c1 = &dummy1;
 
    }
 
    if ((*w)->c2->Index()<0) {
 
      dchem_c2 = &dummy2;
 
    }
 
    plugin->CelltoCellTransport(*w, dchem_c1, dchem_c2); 
 

	
 
    //(*tf)(*w, &(derivs[(*w)->c1->Index() * nchems]),
 
    //&(derivs[(*w)->c2->Index() * nchems] ) );
 
    i+=2*nchems;
 
  }
 
}
 

	
 
void Mesh::setValues(double x, double *y) {
 

	
 
  //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;
 

	
 
  // Layout of derivatives: cells [ chem1 ... chem n]  walls [ [ w1(chem 1) ... w1(chem n) ] [ w2(chem 1) ... w2(chem n) ] ]
 

	
 
  int i=0;
 
  static int emit_count=0;
 
  const int stride = 100;
 
  for (vector<Cell *>::iterator c=cells.begin(); c!=cells.end(); c++) {
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*c)->SetChemical(ch, y[i+ch]);
 
    }
 
    if ( !(emit_count%stride)) {
 
      (*c)->EmitValues(x);
 
    }
 
    i+=nchems;
 
  }
 

	
 
  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*w)->setTransporters1(ch,y[i+ch]);
 
    }
 
    i+=nchems;
 

	
 
    for (int ch=0;ch<nchems;ch++) {
 
      (*w)->setTransporters2(ch,y[i+ch]);
 
    }
 
    i+=nchems;
 
  }
 
  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 *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()+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+1];
 
  int nph=chainHull_2D(p,np,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;
 
  }
 

	
 
  
 
}
 

	
 

	
 
// DataExport
 
void Mesh::CSVExportWallData(QTextStream &csv_stream) const {
 

	
 
  csv_stream << "\"Wall Index\",\"Cell A\",\"Cell B\",\"Length\"";
 
  
 
  for (int c=0;c<Cell::NChem(); c++) {
 
    csv_stream << ",\"Transporter A:" << c << "\"";
 
  }
 
  for (int c=0;c<Cell::NChem(); c++) {
 
    csv_stream << ",\"Transporter B:" << c << "\"";
 
  }
 
  csv_stream << endl;
 
  for (list<Wall *>::const_iterator i=walls.begin();
 
       i!=walls.end();
 
       i++) {
 
    csv_stream << (*i)->Index() << ","
 
	       << (*i)->C1()->Index() << ","
 
	       << (*i)->C2()->Index() << ","
 
	       << (*i)->Length();
 
	      
 
    for (int c=0;c<Cell::NChem(); c++) {
 
      csv_stream << "," << (*i)->Transporters1(c);
 
    }
 
    for (int c=0;c<Cell::NChem(); c++) {
 
      csv_stream << "," << (*i)->Transporters2(c);
 
    }
 

	
 
    csv_stream << endl;
 
  }
 

	
 
  
 
}
 

	
 

	
 
void Mesh::CSVExportMeshData(QTextStream &csv_stream) { 
 
  
 
  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, 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 */

Changeset was too big and was cut off... Show full diff anyway

0 comments (0 inline, 0 general)