Generator tracks reliably (branching) filopodia tips and reports them.

Tips are addressed as filoID+branchID (latter unique within only within its filopodium).

src/main_Generator.cpp

@@ -19,7 +19,27 @@ void ParamsSetup(void);
 ActiveMesh mesh;
 
 //for filopodia tip tracking:
-std::vector< std::vector<int> > tipOffsets;
+// GM addresses filopodia branches with pairs: filoTreeID + branchID;
+// we linearize them with tipIndex() (e.g., as filoTreeID*100+branchID);
+// this is the first index (primary key) in tipOffsets;
+// the second index (secondary key) is time point and
+// the value is actually a pixel offset (encoding a tip position)
+//
+//time-course (vector index) of pixel offsets (vector value)
+class timeVector: public std::vector<int>
+{
+  public:
+	void Init(const size_t length=200)
+	{ this->insert(this->end(),length,-1); }
+};
+//
+//this one holds the entire time-lapse of pixel offsets of all filopodia tips
+//use tipIndex() below to address filopodia branches
+std::map<int,timeVector> tipOffsets;
+//
+//use exclusively this one to address filopodia branches in tipOffsets above
+int tipIndex(const int filoID,const int branchID)
+{ return (filoID*100 +branchID); }
 
 bool LoadNewMesh(const char* path,const int ID,const int fileNo,const bool keepTipTrajectories=false);
 int  RenderMesh(const char* path,const int ID,int fileNo);
@@ -53,13 +73,6 @@ int main(int argc,char **argv)
 	//I'm called as generator/simulator, init the situation...
 	ParamsSetup();
 
-	//for filopodia tip tracking:
-	//container for 200 (all possible) filopodia tips,
-	//initiated for 200 timeframes with position undefined (yet)
-	tipOffsets.reserve(200);
-	for (unsigned int qq=0; qq < 200; ++qq)
-		tipOffsets.push_back( std::vector<int>(200,-1) );
-
 	std::cout << "FIRST RUN: reading all meshes to figure out minimal bounding box.\n";
 	//doubling iteration variables not to loose their values
 	int ii(i);                           //starting index
@@ -111,7 +124,7 @@ int main(int argc,char **argv)
 	          << maxBB.x << "," << maxBB.y << "," << maxBB.z << "] in microns\n";
 
 	std::cout << "SECOND RUN: reading all meshes to conduct the texture simulation and rendering.\n";
-	while (LoadNewMesh(argv[1],ID,i) && (framesToGoHome > 0))
+	while (LoadNewMesh(argv[1],ID,i,true) && (framesToGoHome > 0))
 	{
 		RenderMesh(argv[2],ID,i);
 		++i;
@@ -120,12 +133,14 @@ int main(int argc,char **argv)
 
 	//report the tip trajectories
 	std::ofstream tFile("tip_trajectories.txt");
-	for (unsigned int qq=0; qq < mesh.GetFilopodiaNumber(); ++qq)
+	std::map<int,timeVector>::const_iterator tipIter=tipOffsets.begin();
+	while (tipIter != tipOffsets.end())
 	{
-		tFile << qq << ":";
+		tFile << tipIter->first << ":";
 		for (int t=0; t < i; ++t)
-			tFile << tipOffsets[qq][(unsigned)t] << ",";
+			tFile << tipIter->second.at((unsigned)t) << ",";
 		tFile << "\n";
+		++tipIter;
 	}
 	tFile.close();
 
@@ -211,6 +226,28 @@ bool LoadNewMesh(const char* path,const int ID,const int fileNo,const bool keepT
 
 		retval=mesh.ImportVTK_Ftree(fn,999999);
 		if (retval != 1) std::cout << "reading filopodium: " << fn << "\n";
+		if (retval == 0 && keepTipTrajectories == true)
+		{
+			//managed to import filopodium (possibly with its offspring branches)
+			//remember at what indices in mesh.Ftree are these stored
+			//(as we cannot yet remember directly the pixel offset, we at least
+			//remember index within mesh.Ftree and translate to pixel offset later)
+			//
+			//iterate over the internal mapping (just_branchID -> Ftree indices)
+			std::map<int,unsigned int>::const_iterator iterTips=mesh.tipOffsets.begin();
+			while (iterTips != mesh.tipOffsets.end())
+			{
+				const int idx=tipIndex(filoNo,iterTips->first);
+				if (tipOffsets[idx].size() == 0)
+				{
+					//we have found a new pair filoID+branchID
+					tipOffsets[idx].Init();
+				}
+				tipOffsets[idx].at((unsigned)fileNo)=(signed)iterTips->second;
+
+				++iterTips;
+			}
+		}
 	} while (retval < 2 && filoNo < 100);
 	//NB: retval == 0 when read successfully; == 1 when file open error
 
@@ -401,14 +438,23 @@ int RenderMesh(const char* path,const int ID,int fileNo)
 	const float xOff=mask.GetOffset().x;
 	const float yOff=mask.GetOffset().y;
 	const float zOff=mask.GetOffset().z;
-	for (unsigned int i=0; i < mesh.GetFilopodiaNumber(); ++i)
+
+	//scan over all branches we recognize at this time-point
+	std::map<int,timeVector>::iterator tipIter=tipOffsets.begin();
+	while (tipIter != tipOffsets.end())
 	{
-		const Vector3FC& realPos=mesh.Ftree[i].fPoints[mesh.Ftree[i].segToPoint.back()];
-		//get pixel coordinate
-		const size_t x=size_t((realPos.x-xOff) *xRes);
-		const size_t y=size_t((realPos.y-yOff) *yRes);
-		const size_t z=size_t((realPos.z-zOff) *zRes);
-		tipOffsets[i][(unsigned)fileNo]=(int)mask.GetIndex(x,y,z);
+		//if there is a record for the present time-point, translate it
+		if (tipIter->second.at((unsigned)fileNo) > -1)
+		{
+			const unsigned int i=(unsigned)tipIter->second.at((unsigned)fileNo);
+			const Vector3FC& realPos=mesh.Ftree[i].fPoints[mesh.Ftree[i].segToPoint.back()];
+			//get pixel coordinate
+			const size_t x=size_t((realPos.x-xOff) *xRes);
+			const size_t y=size_t((realPos.y-yOff) *yRes);
+			const size_t z=size_t((realPos.z-zOff) *zRes);
+			tipIter->second.at((unsigned)fileNo)=(int)mask.GetIndex(x,y,z);
+		}
+		++tipIter;
 	}
 
 	return(0);