Commit d515ee65 authored by Vladimír Ulman's avatar Vladimír Ulman
Browse files

Added two functions Cell::SOFA_InitialExport() and Cell::SOFA_NextImport(), 

which are implemented in the new file src/cellSOFA.cpp. This is related to
the SOFA framework interconnection.

The functions provide easy way to export various data from the simulator
and to import subset of the data back to the simulator. Most of it is
controlled via the configuration file into which a new section "[sofa]"
has been introduced. This section is only mandatory iff SOFA feature is
turned on in the CMake.

M    CMakeLists.txt		- acknowledges the new cellSOFA.cpp file
A    src/cellSOFA.cpp	- SOFA export and import functions
M    src/cell.h			- added their declarations (with poor doc for now)
								- initialization of Cell::chrCentres is now in
								  the cell's constructor
M    src/cellChr.cpp		- the init of chrCentres has been removed from
								  the Cell::ChrInitialDistribution()

M    ISBI2015/config-smallVOI-lowSNR.ini	- a sample of the [sofa] section
parent cd4a9377

CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -134,6 +134,7 @@ SET (SRCS 
		  src/cell_08_G2Phase.cpp

		  src/cellChr.cpp

		  src/cellScm.cpp

		  src/cellSOFA.cpp

		  src/toolbox/rnd_generators.cpp 

		  src/toolbox/flowfields.cpp 

		  src/toolbox/collisions.cpp

ISBI2015/config-smallVOI-lowSNR.ini

+90 −5
Original line number Diff line number Diff line
@@ -114,6 +114,13 @@ cell nonrigid strong deformation = 0.3 
; short sequences:

cell sparse nonrigid deformations = false



;

; Refer to the [sofa] section for parameters related to performing

; non-rigid deformations with the SOFA framework.

;

; Note that the SOFA option must be also enabled in the CMake.

;



;

; How far a cell can "see" when it is trying to move

; to region with lower density of cells (in microns).
@@ -149,6 +156,84 @@ number of frames per cell cycle = 50 
intensity of nonspecific staining = 2





[sofa]

;

; Note that the SOFA option must be also enabled in the CMake in order

; to make the simulator consider the content of this section.

;



;

; The simulation with the SOFA works in two modes:

; - Export mode:	the simulator starts as usual, creates its internal structures etc.,

; 						but then, instead of actually starting the simulation itself, it

; 						rather exports relevant data for the SOFA including the initial scene*

; 						images, and quits

; - Play mode:		the simulator starts as usual except that some structures/data is

;						not created but rather imported from the SOFA data, the simulation is

;						then started as usual (except for some new SOFA stuff in the atomic

;						functions... ;-)

;

; The play mode is default unless this option equals to 1 (one).

; Set this one to 0 (zero) to work in the play mode.

export mode = 0



;

; For creating new data for a new frame of the simulator, the SOFA usually requires to

; calculate several (however, fixed) internal iterations (which are, however, exported

; as well). This parameter specifies multiplication factor to match simulator frame no.

; with corresponding frame no. of the SOFA. In other words, this parameter minus one tells

; how many internal iterations were used in the SOFA.

frame skip step = 5



;

; Provide template file names used when exporting data to the SOFA:

; Include one '%u' for the cell ID. One can disable certain export

; by disabling/commenting respective parameter.

;

; Note: BP stands for "Boundary Points" -- yes, it has relevance

; to the Cell::BPLists.

;

; Note: "mask voxels" make up an initial cell mask. The file will

; contain the same number of lines as there are voxels in the mask.

;

; The syntax follows that of the standard C library printf() function.

export name nucleus mask img    = ID%02u_init_nucleusMask.mhd

;export name nucleus mask voxels = ID%02u_init_nucleusMask_coordinates.txt

;export name nucleus mask grid   = ID%02u_init_gridPoints.txt

export name nucleus BP          = ID%02u_init_nucleusOuterBoundaryPoints.txt

export name nucleoli1 BP        = ID%02u_init_nucleoli1OuterBoundaryPoints.txt

export name nucleoli2 BP        = ID%02u_init_nucleoli2OuterBoundaryPoints.txt

export name chromatin molecules = ID%02u_init_chromatinPoints.txt

;export name nucleus deform img  = ID%02u_init_nucleusOuterBoundaryPoints_DeformationMagnitudes.ics



;

; Similar stuff but for deformation hints for particular cells.

; Use somewhere first '%u' for the cell ID, later '%lu' for frame

; number.

export name nucleus deformation = ID%.2u_T%03lu_nucleusOuterBoundaryPoints_DeformationMagnitudes.txt



;

; This number minus one gives number of voxels between two

; nearby grid lines.

export mask grid stepping = 7



;

; Provide template file names used when importing data from the SOFA:

; Include first '%u' for the cell ID and second '%lu' for frame number.

; One can disable certain import by disabling/commenting respective parameter.

;

; The syntax follows that of the standard C library printf() function.

import name nucleus BP          = ../deformationsID01-05/all/ID%02u_T%03lu_nucleusOuterBoundaryPoints.txt

import name nucleoli1 BP        = ../deformationsID01-05/all/ID%02u_T%03lu_nucleoli1OuterBoundaryPoints.txt

import name nucleoli2 BP        = ../deformationsID01-05/all/ID%02u_T%03lu_nucleoli2OuterBoundaryPoints.txt

import name chromatin molecules = ../deformationsID01-05/all/ID%02u_T%03lu_chromatinPoints.txt



; import name nucleus BP          = ID%02u_T%03lu_nucleusOuterBoundaryPoints.txt

; import name nucleoli1 BP        = ID%02u_T%03lu_nucleoli1OuterBoundaryPoints.txt

; import name nucleoli2 BP        = ID%02u_T%03lu_nucleoli2OuterBoundaryPoints.txt

; import name chromatin molecules = ID%02u_T%03lu_gridPoints.txt





;

;-------------- characteristics of scheduler -------------------

;
@@ -250,8 +335,8 @@ level = 1 
; receive their own IDs, independently of the mask image.

;

; Used in the Scheduler constructor. Should be positive natural numbers.

init mask min label = 1

init mask max label = 1

;init mask min label = 1

;init mask max label = 1



;

; For non-rigid deformations, the simulator uses a coherent deformation pool
@@ -273,8 +358,8 @@ pnArray z length = 52 
; 'strong deformations' is normally inferred the duration of the G2 cell

; cycle phase. This parameter may override this average delay [in frames].

;

; Used in the Cell constructor and Cell::ScmSuggestNonrigidDeformation().

; Should be positive natural number.

; Used in the Cell constructor, Cell::ScmSuggestNonrigidDeformation(),

; and Cell::SOFA_InitialExport_Pt2(). Should be positive natural number.

strong deformation delay = 10



;
@@ -292,5 +377,5 @@ single run length = 1 
; value.

;

; Used in both Scheduler constructors. Should be positive natural number.

cell initial G2 life = 25

cell initial G2 life = 50

src/cell.h

+19 −0
Original line number Diff line number Diff line
@@ -75,6 +75,15 @@ template <class MV, class PV> class Cell 
						cellCycleLength = 

								  (int) configIni["cell"]["number of frames per cell cycle"];

						chrCount = (int) configIni["cell"]["chromosome count"];



						//resize the vectors in advance to avoid their incremental increasing, which is slow

						//(and copies the early dots over and over again... which breaks fluency of IDs :-)

						chrCentres.reserve(chrCount);

						//

						// this way we create continuous ascending sequence of unique IDs

						for (size_t i=0; i<chrCount; i++)

							 chrCentres.push_back(Dot());



						chrSpreads = std::vector<float>(chrCount, 0.0f);

						chrDensities = std::vector<float>(chrCount, 0.0f);
@@ -286,6 +295,16 @@ template <class MV, class PV> class Cell 
			 /// returns comet type

			 StateOfComet GetCometState(void) const { return scmCometState; };



			 /**

			  * SOFA: initial export

			  */

			 void SOFA_InitialExport(void);



			 /**

			  * SOFA: next import (for regular imports during simulation)

			  */

			 void SOFA_NextImport(const size_t T);



  private:

			 /// structure loaded external ini file

			 const IniHandler &configIni;

src/cellChr.cpp

+0 −12
Original line number Diff line number Diff line
@@ -1486,25 +1486,13 @@ void Cell<MV,PV>::ChrInitialDistribution(const i3d::Image3d<bool> &mask, 
	chrDots = numberOfDots/chrCount;

	DEBUG_REPORT("going to use " << chrDots << " molecules for a chromosome, as " << numberOfDots << " is required total");



	//resize the vector in advance to avoid its incremental increasing,

	//which copies the early dots over and over again...

	chrCentres.reserve(chrCount);



	// this way we create continuous ascending sequence of unique IDs

	for (size_t i=0; i<chrCount; i++)

	{

		 chrCentres.push_back(Dot());

	}



	//resize the vector in advance to avoid its incremental increasing,

	//which copies the early dots over and over again...

	chrDotList.reserve(numberOfDots);



	// this way we create continuous ascending sequence of unique IDs

	for (size_t i=0; i<numberOfDots; i++)

	{

		 chrDotList.push_back(Molecule());

	}



	//time savers: resolution [px/um]

	const i3d::Vector3d<float> res = mask.GetResolution().GetRes();

src/cellSOFA.cpp

0 → 100644
+225 −0
Original line number Diff line number Diff line 
/**

 * This file contains implementation of Cell::SOFA_* methods which handle

 * the interconnection between the simulator and the SOFA framework.

 * The framework is used to provide a nice non-rigid deformations.

 */



#ifdef GTGEN_DEBUG

  #include <sys/time.h>

#endif

#include <fstream>



#include "cell.h"

#include "toolbox/bp_lists.h"

#include "toolbox/importexport.h"



template <class MV, class PV>

void Cell<MV, PV>::SOFA_InitialExport(void)

{

	DEBUG_REPORT("doing initial export for cell ID=" << this->ID);



	//for keeping file names

	char fn[100];



	//should we export the mask image, list of its voxels, or list of grid points?

	if (  (configIni["sofa"].present("export name nucleus mask img"))

		|| (configIni["sofa"].present("export name nucleus mask voxels"))

		|| (configIni["sofa"].present("export name nucleus mask grid")) )

	{

		//first, create the local small mask image

		i3d::Image3d<MV> maskNucleus;

		maskNucleus.SetResolution(scheduler.sceneRes);

		maskNucleus.SetOffset(scheduler.sceneMasks[timePoint]->GetOffset());

		AllocateAndZeroNewMask(maskNucleus, scmCellBPList, 2, cellLookDistance);



		//now render the cell into it

		RenderBPListIntoMask(maskNucleus,(MV)1,

		  scmNucleusBPList,scmNucleusBPCentre,scmNucleusOuterBPNumber);

#ifdef ENABLE_NUCLEOLI

		RenderBPListIntoMask(maskNucleus,(MV)2,

		  scmNucleoli1BPList,scmNucleoli1BPCentre,scmNucleoli1OuterBPNumber);

		RenderBPListIntoMask(maskNucleus,(MV)3,

		  scmNucleoli2BPList,scmNucleoli2BPCentre,scmNucleoli2OuterBPNumber);

		//now in the maskNucleus image: 0 means outside,

		//1 nucleus, 2 and 3 each nucleolus

#else

		throw ERROR_REPORT("There are no nucleoli in the nucleus.");

#endif



		//save/export the mask image with the 3 labels ?

		if (configIni["sofa"].present("export name nucleus mask img"))

		{

			sprintf(fn,

				(const char*)configIni["sofa"]["export name nucleus mask img"], this->ID);

			maskNucleus.SaveImage(fn);

		}



		//export the list of mask coordinates ?

		if (configIni["sofa"].present("export name nucleus mask voxels"))

		{

			sprintf(fn,

				(const char*)configIni["sofa"]["export name nucleus mask voxels"], this->ID);

			ExportAllMaskVoxelCoordinates(maskNucleus,(MV)0,fn);

		}



		//export the list of grid coordinates within the mask ?

		if (configIni["sofa"].present("export name nucleus mask grid"))

		{

			sprintf(fn,

				(const char*)configIni["sofa"]["export name nucleus mask grid"], this->ID);

			ExportGridMaskVoxelCoordinates(maskNucleus,(MV)1,

								(int)configIni["sofa"]["export mask grid stepping"],fn);

		}

	}



	//export the "outer" list of nucleus boundary/countour points ?

	if (configIni["sofa"].present("export name nucleus BP"))

	{

		sprintf(fn,

			(const char*)configIni["sofa"]["export name nucleus BP"], this->ID);

		ExportPointList(scmNucleusBPList,scmNucleusOuterBPNumber,fn);

	}



	//export the "outer" list of the first nucleoli boundary/countour points ?

	if (configIni["sofa"].present("export name nucleoli1 BP"))

	{

		sprintf(fn, 

			(const char*)configIni["sofa"]["export name nucleoli1 BP"], this->ID);

		ExportPointList(scmNucleoli1BPList,scmNucleoli1OuterBPNumber,fn);

	}



	//export the "outer" list of the first nucleoli boundary/countour points ?

	if (configIni["sofa"].present("export name nucleoli2 BP"))

	{

		sprintf(fn, 

			(const char*)configIni["sofa"]["export name nucleoli2 BP"], this->ID);

		ExportPointList(scmNucleoli2BPList,scmNucleoli2OuterBPNumber,fn);

	}



	//export the list chromatin points ?

	if (configIni["sofa"].present("export name chromatin molecules"))

	{

		sprintf(fn, 

			(const char*)configIni["sofa"]["export name chromatin molecules"], this->ID);

		ExportPointList(chrDotList,chrDotList.size(),fn);

	}





	//prepare non-rigid deformation plan for export

	//but only if really needed..

	if (  (configIni["sofa"].present("export name nucleus deform img"))

		|| (configIni["sofa"].present("export name nucleus deformation")) )

					scheduler.PnGetMeNonrigidHints(scmNucleusBPList,

											scmNucleusOuterBPNumber,scmPnHints);

	//note that the previous function makes sure:

	//		scmPnHints.GetSizeX() == scmNucleusOuterBPNumber



	// save the deformation hints image ?

	if (configIni["sofa"].present("export name nucleus deform img"))

	{

		sprintf(fn, 

			(const char*)configIni["sofa"]["export name nucleus deform img"], this->ID);

		scmPnHints.SaveImage(fn);

	}



	// save the sequence of deformation recipies ?

	if (configIni["sofa"].present("export name nucleus deformation"))

	{

		const float *myHintPtr=scmPnHints.GetFirstVoxelAddr();

		for (size_t curTime=0; curTime < scmPnHints.GetSizeY(); ++curTime)

		{

			sprintf(fn,

				(const char*)configIni["sofa"]["export name nucleus deformation"],

				this->ID,curTime);

			std::ofstream Plist(fn);



			if (!Plist.is_open()) {

				std::ostringstream err;

				err << "Can't open the output file " << fn << ".";

				throw ERROR_REPORT(err.str());

			}



			float mult=1.f;

			if (curTime == scmPnNextStrongUsage)

			{

				//the code is a copy from Cell::ScmSuggestNonrigidDeformation()

				float PnStrongInitDelay=(float)cellCycleLength/4.f;



				if (configIni["testing"].present("strong deformation delay"))

				{

					REPORT("TESTING MODE: enforcing average strong deformation delay to "

						<< (int)configIni["testing"]["strong deformation delay"] << " frames.");

					PnStrongInitDelay=(float)configIni["testing"]["strong deformation delay"];

				}



				scmPnNextStrongUsage+=(size_t)std::max(GetRandomGauss(PnStrongInitDelay,PnStrongInitDelay/4.f),1.f);



				DEBUG_REPORT("using stronger hints this time; next such event will be at " << scmPnNextStrongUsage);



				mult=cellNonRigidStrongDeformation/cellNonRigidDeformation;

			}



			for ( size_t BPindex=0;

					(BPindex < scmPnHints.GetSizeX()) && (Plist.good());

					++BPindex, ++myHintPtr )

				Plist << *myHintPtr *mult << "\n";



			if (!Plist.good())

			{

				Plist.close();

				std::ostringstream err;

				err << "Error while writing file " << fn << ".";

				throw ERROR_REPORT(err.str());

			}

			Plist.close();

		}

	}

}



//-----------------------------------------------------------------------------

template <class MV, class PV>

void Cell<MV, PV>::SOFA_NextImport(const size_t T)

{

	DEBUG_REPORT("doing next import for cell ID=" << this->ID

						<< " at time " << T);



	//for keeping file names

	char fn[100];



	//"cache" the skipStep

	const size_t skipStep = (int)configIni["sofa"]["frame skip step"];



	//nucleus import

	sprintf(fn,

			(const char*)configIni["sofa"]["import name nucleus BP"],

			this->ID,skipStep*T);

	ImportPointList(scmNucleusBPList,&scmNucleusBPCentre,scmNucleusOuterBPNumber,fn);



#ifdef ENABLE_NUCLEOLI

	//nucleoli1 import

	sprintf(fn,

			(const char*)configIni["sofa"]["import name nucleoli1 BP"],

			this->ID,skipStep*T);

	ImportPointList(scmNucleoli1BPList,&scmNucleoli1BPCentre,scmNucleoli1OuterBPNumber,fn);



	//nucleoli2 import

	sprintf(fn,

			(const char*)configIni["sofa"]["import name nucleoli2 BP"],

			this->ID,skipStep*T);

	ImportPointList(scmNucleoli2BPList,&scmNucleoli2BPCentre,scmNucleoli2OuterBPNumber,fn);

#endif



	//just movement of dots:

	//new positions are given by the SOFA, import the list

	sprintf(fn,

			(const char*)configIni["sofa"]["import name chromatin molecules"],

			this->ID,skipStep*T);

	ImportPointList(chrDotList,NULL,chrDotList.size(),fn);

}







// explicit instantiations to force explicit build of the functions in this file

// because otherwise the linker couldn't find them

template class Cell<i3d::GRAY8,  i3d::GRAY16>;

template class Cell<i3d::GRAY16, i3d::GRAY16>;