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

Added two alternative functions for chromatin movement, 

they to the same job as their original version, but do
corrections differently -- no DT transform is used, but
just slow shift towards chromatin centre is used instead,
which is a lot faster (but seems to be a little error prone).

Added:
Cell::ChrMoveWithBrown() and Cell::ChrMoveByFlow() with no dmOfBackground parameter.

M    cell.h
M    cellChr.cpp
parent 7a7ae91a

src/cell.h

+36 −0
Original line number Diff line number Diff line
@@ -546,6 +546,8 @@ template <class MV, class PV> class Cell 
			  * Move chromosome with specified \e chrID (all its molecules) by the vector

			  * \e shift but keep molecules within the \e mask. Uses mask image of the

			  * template type MV.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveByVector(const size_t chrID, 

										 i3d::Vector3d<float> const &shift,
@@ -556,6 +558,8 @@ template <class MV, class PV> class Cell 
			 /**

			  * Move chromosome with specified \e chrID (all its molecules) by the vector

			  * \e shift but keep molecules within the \e mask. Uses binary (bool) mask.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveByVectorBM(const size_t chrID, 

											i3d::Vector3d<float> const &shift,
@@ -567,12 +571,24 @@ template <class MV, class PV> class Cell 
			  * Perform movement of molecules belonging to chromosome specified by \e chrID

			  * while considering the \e mask. The magnitude of performed movement is controlled

			  * with the parameter \e step. Uses mask image of the template type MV.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveWithBrown(const size_t chrID, 

										  const i3d::Vector3d<float> step,

										  i3d::Image3d<MV> const &mask, 

										  i3d::Image3d<float>* &dmOfBackground,

										  const MV value);

 			 /**

			  * Perform movement of molecules belonging to chromosome specified by \e chrID

			  * while considering the \e mask. The magnitude of performed movement is controlled

			  * with the parameter \e step. Uses mask image of the template type MV.

			  * If a chromatin dot gets out of the mask, it is moved towards its centre.

			  */

			 void ChrMoveWithBrown(const size_t chrID, 

				                    const i3d::Vector3d<float> step,

				                    i3d::Image3d<MV> const &mask,

				                    const MV value);



			 /// wrapper for low level function

			 void ChrMoveWithBrown(const size_t chrID,
@@ -592,6 +608,8 @@ template <class MV, class PV> class Cell 
			  * Perform movement of molecules belonging to chromosome specified by \e chrID

			  * while considering the \e mask. The magnitude of performed movement is controlled

			  * with the parameter \e step. Uses binary (bool) mask.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveWithBrownBM(const size_t chrID, 

											 const i3d::Vector3d<float> step,
@@ -617,6 +635,8 @@ template <class MV, class PV> class Cell 
			  * Perform gathering of molecules belonging to chromosome specified by \e chrID

			  * within the given \e mask. The magnitude of performed movement is controlled

			  * with the parameter \e step. Uses mask image of the template type MV.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveWithGravity(const size_t chrID, 

											 const float step,
@@ -628,6 +648,8 @@ template <class MV, class PV> class Cell 
			  * Perform gathering of molecules belonging to chromosome specified by \e chrID

			  * within the given \e mask. The magnitude of performed movement is controlled

			  * with the parameter \e step. Uses binary (bool) mask.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveWithGravityBM(const size_t chrID, 

												const float step,
@@ -638,15 +660,27 @@ template <class MV, class PV> class Cell 
			 /**

			  * Perform movement of all dots using optical flow \e FF but keep the dots

			  * within the \e mask. Uses mask image of the template type MV.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveByFlow(FlowField<float> const &FF,

									  i3d::Image3d<MV> const &mask, 

									  i3d::Image3d<float>* &dmOfBackground,

									  const MV value);

			 /**

			  * Perform movement of all dots using optical flow \e FF but keep the dots

			  * within the \e mask. Uses mask image of the template type MV.

			  * If a chromatin dot gets out of the mask, it is moved towards its centre.

			  */

			 void ChrMoveByFlow(FlowField<float> const &FF,

									  i3d::Image3d<MV> const &mask, 

									  const MV value);

			 /**

			  * Perform movement of dots specified by \e chrID using optical flow \e FF

			  * but keep the dots within the \e mask. Uses mask image of the template 

			  * type MV.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveByFlow(const size_t chrID,

									  FlowField<float> const &FF,
@@ -657,6 +691,8 @@ template <class MV, class PV> class Cell 
			 /**

			  * Perform movement of all dots using optical flow \e FF but keep the dots

			  * within the \e mask. Uses binary (bool) mask.

			  * If a chromatin dot gets out of the \e mask, it is moved in the gradient of the

			  * distance map \e dmOfBackground until it reaches the mask.

			  */

			 void ChrMoveByFlowBM(FlowField<float> const &FF,

										 i3d::Image3d<bool> const &mask,

src/cellChr.cpp

+170 −0
Original line number Diff line number Diff line
@@ -707,6 +707,84 @@ void Cell<MV, PV>::ChrMoveWithBrownBM(const size_t chrID, 
							 value);

}



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

template <class MV, class PV>

void Cell<MV, PV>::ChrMoveWithBrown(const size_t chrID,

												const i3d::Vector3d<float> step,

												i3d::Image3d<MV> const &mask,

												const MV value)

{

	//sanitary checks...

	if (chrDots * chrCount != chrDotList.size()) {

		std::ostringstream err;

		err << "Cell ID=" << ID << " has bad number (" << chrDotList.size()

		    << ") of chrDotList items, should be " << chrDots*chrCount;

		throw ERROR_REPORT(err.str());

	}

	if (chrID >= chrCount) {

		std::ostringstream err;

		err << "Cannot work with the asked chromosome no. " << chrID

		    << ", I have only " << chrCount;

		throw ERROR_REPORT(err.str());

	}



	//time savers:

	const float xRes=mask.GetResolution().GetX();

	const float yRes=mask.GetResolution().GetY();

	const float zRes=mask.GetResolution().GetZ();

	const float xOff=mask.GetOffset().x;

	const float yOff=mask.GetOffset().y;

	const float zOff=mask.GetOffset().z;



	//sweep over given portion of the chrDotList

	for (size_t i=chrDots*chrID; i < chrDots*(chrID+1); ++i) {

		//a working copy of the new coordinate

		//const i3d::Vector3d<float> &D=chrDotList[i].pos;

		const i3d::Vector3d<float> D(chrDotList[i].pos);



		//up to 3x times Brownian motion while we are not in the mask

		i3d::Vector3d<float> V; //testing vector

		char tries=0;		//number of tries

		int x,y,z;		//new pixel position

		do {

			SuggestBrownianVector(V, step);

			++tries;



			//nearest neighbor mask image pixel coordinate

			x=(int)roundf( (D.x+V.x-xOff) *xRes );

			y=(int)roundf( (D.y+V.y-yOff) *yRes );

			z=(int)roundf( (D.z+V.z-zOff) *zRes );

			//REPORT(D << " (tries " << (int)tries << "): + " << V << " = (" << x << "," << y << "," << z << ")");

		} while ( ((!mask.Include(x,y,z)) || (mask.GetVoxel(x,y,z) != value))

			&& (tries < 3) );



		//move the point to its new position

		if (tries < 3) 

			 chrDotList[i].pos=D+V;

		else if (mask.Include(x,y,z) && (mask.GetVoxel(x,y,z) != value))

		{

			//broken chromatin we shall move towards centre:

			//incremental vector from the last suggested position towards the chromatin centre

			const i3d::Vector3d<float> ds=(chrCentres[chrID].pos - (D+V)) / 20.f;



			//iteratively shift the last failing position

			int j=0;

			do {

				++j;

				//nearest neighbor mask image pixel coordinate

				x=(int)roundf( (D.x+V.x+(float)j*ds.x-xOff) *xRes );

				y=(int)roundf( (D.y+V.y+(float)j*ds.y-yOff) *yRes );

				z=(int)roundf( (D.z+V.z+(float)j*ds.z-zOff) *zRes );

			} while ((j < 20) && (mask.GetVoxel(x,y,z) != value));

			//note: it is believed that chrCentre is in the image and the path

			//from current (also in the image) position lies in the image entirely



			chrDotList[i].pos=D+V +(float)j*ds;

			if (j == 20) DEBUG_REPORT("chromatin issue A");

		}

		// else do no motion

	}

}



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

/**
@@ -1041,6 +1119,98 @@ void Cell<MV, PV>::ChrMoveByFlowBM(FlowField<float> const &FF, 
	__ChrMoveByFlow(chrDotList, FF,mask,dmOfBackground, value);

}



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

template <class MV, class PV>

void Cell<MV, PV>::ChrMoveByFlow(FlowField<float> const &FF,

											i3d::Image3d<MV> const &mask,

											const MV value)

{

	//sanitary checks...

	if (!FF.isConsistent())

		throw ERROR_REPORT("Flow field is not consistent.");

	if (!FF.z)

		throw ERROR_REPORT("The input flow field is not 3D.");



	//process the whole chrDotList

	//time savers:

	const float xRes=mask.GetResolution().GetX();

	const float yRes=mask.GetResolution().GetY();

	const float zRes=mask.GetResolution().GetZ();

	const float xOff=mask.GetOffset().x;

	const float yOff=mask.GetOffset().y;

	const float zOff=mask.GetOffset().z;



	const float xResFF=FF.x->GetResolution().GetX();

	const float yResFF=FF.x->GetResolution().GetY();

	const float zResFF=FF.x->GetResolution().GetZ();

	const float xOffFF=FF.x->GetOffset().x;

	const float yOffFF=FF.x->GetOffset().y;

	const float zOffFF=FF.x->GetOffset().z;



	//move dot by dot...

	for (size_t i=0; i < chrDotList.size(); ++i) 

	{

		const float X=chrDotList[i].pos.x;

		const float Y=chrDotList[i].pos.y;

		const float Z=chrDotList[i].pos.z;



		//the position in pixels of the flow field

		const float Xp=(X -xOffFF) *xResFF;

		const float Yp=(Y -yOffFF) *yResFF;

		const float Zp=(Z -zOffFF) *zResFF;



		//new molecule position after moving along the flow field

		//note: GetPixel() returns 0 in case we ask for value outside the image molecule position

		//(both the position and vectors in the flow field are in microns)

		const float nX=X + GetPixel(*FF.x, Xp,Yp,Zp);

		const float nY=Y + GetPixel(*FF.y, Xp,Yp,Zp);

		const float nZ=Z + GetPixel(*FF.z, Xp,Yp,Zp);



		//new position in the pixels of the mask image, using nearest neighbor

		const int nXp=(int)roundf( (nX-xOff) *xRes );

		const int nYp=(int)roundf( (nY-yOff) *yRes );

		const int nZp=(int)roundf( (nZ-zOff) *zRes );



		if (mask.Include(nXp,nYp,nZp))

		{

			 if (mask.GetVoxel(nXp,nYp,nZp)==value) 

			 {

				  //new position fits into the image and even into this cell mask,

				  // update the dot then (with the real unrounded position)

				  chrDotList[i].pos.x=nX;

				  chrDotList[i].pos.y=nY;

				  chrDotList[i].pos.z=nZ;

			 }

			 else

			 {

				//broken chromatin we shall move towards centre:

				//incremental vector from the last suggested position towards the chromatin centre

				const i3d::Vector3d<float> ds=(chrCentres[i/chrDots].pos - i3d::Vector3d<float>(nX,nY,nZ)) / 20.f;



				//iteratively shift the last failing position

				int j=0;

				int x,y,z;

				do {

					++j;

					//nearest neighbor mask image pixel coordinate

					x=(int)roundf( (nX+(float)j*ds.x-xOff) *xRes );

					y=(int)roundf( (nY+(float)j*ds.y-yOff) *yRes );

					z=(int)roundf( (nZ+(float)j*ds.z-zOff) *zRes );

				} while ((j < 20) && (mask.GetVoxel(x,y,z) != value));

				//note: it is believed that chrCentre is in the image and the path

				//from current (also in the image) position lies in the image entirely



				chrDotList[i].pos=i3d::Vector3d<float>(nX,nY,nZ) +(float)j*ds;

				if (j == 20) DEBUG_REPORT("chromatin issue C");

			 }

		}

		else

		{

			 // throw ERROR_REPORT("The movement of dots is too fast.");

		}

	}

}





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

// centrosomes moving functions