From 19a7f90ea832b2a4865a23f7359f083a44d92587 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vladim=C3=ADr=20Ulman?= <ulman@mpi-cbg.de>
Date: Sun, 8 Jan 2017 09:31:18 +0100
Subject: [PATCH] Moved filopodia geometry meta-data (fPoints, segFromPoints
 etc.) into the mesh_t structure, so now every filopodium remembers it. This
 is useful for generating filo texture later. Mesh moving and scaling
 functions were updated to reflect this.

Updated ImportVTK_Ftree() and InitDots_filo() so that the former really only
reads the filo geometry and the latter creates a texture for it.

Filopodia skeleton segments are displayed in OpenGL now by default.
---
 cmath3d/TriangleMesh.cpp | 203 +++++++++++++++++++--------------------
 cmath3d/TriangleMesh.h   |  15 ++-
 src/graphics.cpp         |  38 ++++----
 src/main.cpp             |  18 +++-
 4 files changed, 138 insertions(+), 136 deletions(-)

diff --git a/cmath3d/TriangleMesh.cpp b/cmath3d/TriangleMesh.cpp
index 9df2d59..e7193c7 100644
--- a/cmath3d/TriangleMesh.cpp
+++ b/cmath3d/TriangleMesh.cpp
@@ -490,10 +490,10 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	
 	//local/temporary list of vertices of segments that make up the f-tree
 	//later, we convert these to triangles and save these guys instead
-	fPoints.clear();
-	segFromPoint.clear();
-	segToPoint.clear();
-	segFromRadius.clear();
+	filoMesh.fPoints.clear();
+	filoMesh.segFromPoint.clear();
+	filoMesh.segToPoint.clear();
+	filoMesh.segFromRadius.clear();
 
 	//try to open the file
 	std::ifstream file(filename);
@@ -546,7 +546,7 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	file.ignore(10240,'\n');
 
 	//std::cout << "reading " << itemCount << " point coordinates\n";
-	fPoints.reserve(itemCount);
+	filoMesh.fPoints.reserve(itemCount);
 
 	//read all points...
 	float x,y,z;
@@ -557,7 +557,7 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 		//... and save them
 		Vector3FC v1(x,y,z);
 		//v1*=100.f;
-		fPoints.push_back(v1);
+		filoMesh.fPoints.push_back(v1);
 
 		--itemCount;
 	}
@@ -573,9 +573,9 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	file.ignore(10240,'\n');
 
 	//std::cout << "reading " << itemCount << " segments\n";
-	segFromPoint.reserve(itemCount);
-	  segToPoint.reserve(itemCount);
-	segFromRadius.reserve(itemCount);
+	filoMesh.segFromPoint.reserve(itemCount);
+	  filoMesh.segToPoint.reserve(itemCount);
+	filoMesh.segFromRadius.reserve(itemCount);
 
 	//read all segments vertices
 	int ignore,v1,v2;
@@ -583,8 +583,8 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	{
 		file >> v1 >> v2;
 
-		segFromPoint.push_back(v1);
-		  segToPoint.push_back(v2);
+		filoMesh.segFromPoint.push_back(v1);
+		  filoMesh.segToPoint.push_back(v2);
 
 		--itemCount;
 	}
@@ -623,7 +623,7 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	v1=0;
 	while (itemCount > 0 && file >> x)
 	{
-		segFromRadius.push_back(x*radiusCorrection);
+		filoMesh.segFromRadius.push_back(x*radiusCorrection);
 
 		--itemCount;
 		++v1;
@@ -639,17 +639,16 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	//we need to render the circular-like base (as triangle fan),
 	//and the cylinder-like "cone" towards the tip; this point is for the triangle fan
 	const size_t fanCentrePoint=filoMesh.Pos.size();
-	filoMesh.Pos.push_back(fPoints[segFromPoint[0]]);
+	filoMesh.Pos.push_back(filoMesh.fPoints[filoMesh.segFromPoint[0]]);
 
 	//bookmark the beginning of the cylinder-like surface points
 	size_t firstRadiusPoint=filoMesh.Pos.size();
-	PointsFirstOffset=firstRadiusPoint;
 
 	//will store a former nXaxis
 	Vector3F pXaxis;
 
 	//for all bending points except the last one (tip)
-	for (unsigned int i=0; i < segFromPoint.size(); ++i)
+	for (unsigned int i=0; i < filoMesh.segFromPoint.size(); ++i)
 	{
 		//we need to construct rotation matrix that would
 		//rotate points on a circle which lays in the XZ plane
@@ -658,11 +657,11 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 		//would be normal
 		//
 		//new Y axis 
-		Vector3F nYaxis=fPoints[segToPoint[i]];
-		if (i == 0) nYaxis-=fPoints[segFromPoint[i]];
+		Vector3F nYaxis=filoMesh.fPoints[filoMesh.segToPoint[i]];
+		if (i == 0) nYaxis-=filoMesh.fPoints[filoMesh.segFromPoint[i]];
 		else
 		{
-			nYaxis-=fPoints[segFromPoint[i-1]];
+			nYaxis-=filoMesh.fPoints[filoMesh.segFromPoint[i-1]];
 			nYaxis/=2.f; //unnecessary scaling
 		}
 
@@ -709,21 +708,21 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 		for (int p=0; p < PointsOnRadiusPeriphery; ++p)
 		{
 			//the point in its original position
-			float x=segFromRadius[i]*cosf(6.28f*float(p)/float(PointsOnRadiusPeriphery));
-			float z=segFromRadius[i]*sinf(6.28f*float(p)/float(PointsOnRadiusPeriphery));
+			float x=filoMesh.segFromRadius[i]*cosf(6.28f*float(p)/float(PointsOnRadiusPeriphery));
+			float z=filoMesh.segFromRadius[i]*sinf(6.28f*float(p)/float(PointsOnRadiusPeriphery));
 
 			//rotate
 			Vector3FC V(x*nXaxis);
 			V+=z*nZaxis;
 
 			//shift to the centre and save
-			V+=fPoints[segFromPoint[i]];
+			V+=filoMesh.fPoints[filoMesh.segFromPoint[i]];
 			filoMesh.Pos.push_back(V);
 		}
 	}
 
 	//finally, add the tip point
-	filoMesh.Pos.push_back(fPoints[segToPoint.back()]);
+	filoMesh.Pos.push_back(filoMesh.fPoints[filoMesh.segToPoint.back()]);
 
 	//now, add the triangle fan - aka filopodia base
 	for (size_t p=0; p < PointsOnRadiusPeriphery-1; ++p)
@@ -740,7 +739,7 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 
 	//now create triangles for all segment strips
 	//except the last one (that leads to the tip)
-	for (unsigned int i=1; i < segFromPoint.size(); ++i)
+	for (unsigned int i=1; i < filoMesh.segFromPoint.size(); ++i)
 	{
 		int p=0;
 		for (; p < PointsOnRadiusPeriphery-1; ++p)
@@ -784,63 +783,6 @@ int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
 	filoMesh.ID.push_back(firstRadiusPoint+PointsOnRadiusPeriphery);
 	filoMesh.norm.push_back(fictiveNormal);
 
-	//render the filopodium texture points along its skeleton
-	//the amount of points is determined from the radius, from circle area -> square of radius
-	for (size_t i=0; i < segFromPoint.size(); ++i)
-	{
-		//time-savers...
-		const Vector3F& fP=fPoints[segFromPoint[i]];
-		const Vector3F& tP=fPoints[segToPoint[i]];
-
-		const float length=(tP-fP).Len(); //full segment has length = 0.005
-
-		//stepping along the axis
-		const int S=(int)ceil(5.f*length/0.005); //"fraction factor"
-		//std::cout << "S=" << S << ", length=" << length << "\n";
-		for (int s=0; s < S; ++s)
-		{
-			//ratios "from" and "to"
-			const float fS=float(S-s)/float(S);
-			const float tS=float(s)/float(S);
-
-			//target radius
-			const float target_r=(i+1 < segFromPoint.size())? segFromRadius[i+1] : 0.f;
-
-			//current radius
-			float r=fS*segFromRadius[i] + tS*target_r;
-
-			//current position
-			Vector3F p =fS*fP;
-			         p+=tS*tP;
-
-			//modify coordinates exactly the same as mesh coordinates will be modified
-			//TODO: devise a universal function to shift points to target place
-			//and make all function use this one, so that shift must be changed at only place!
-			p.x*=180.f;
-			p.y*=180.f;
-			p.z*=250.f;
-			p.x+=13.75f;
-			p.y+=13.75f;
-			p.z+=11.00f;
-
-			//insert fl. dots
-			r*=4000.f;
-
-			//std::cout << s << ": r=" << r << ", r*r=" << r*r << "\n";
-
-			//float dotsCount=10.f*r*r;
-			float dotsCount=100.f*r;
-			int qM=(dotsCount > 1300.f)? 1300 : (int)dotsCount;
-
-			for (int q=0; q < qM; ++q)
-			{
-				dots_filo.push_back(p);
-				dots_filo.push_back(p+Vector3F(0.0f,0.0f,+0.1f));
-				dots_filo.push_back(p+Vector3F(0.0f,0.0f,-0.1f));
-			}
-		}
-	}
-
 	return(0);
 }
 
@@ -1057,7 +999,7 @@ void ActiveMesh::RenderMask_filo(i3d::Image3d<i3d::GRAY16>& mask,const int idx,c
 	i3d::FloodFill(tmpI,(i3d::GRAY16)50,0);
 
 	//this removes the flooding while filling
-	//everything else (and closing holes in this was)
+	//everything else (and closing holes in this way)
 	i3d::GRAY16* pI=tmpI.GetFirstVoxelAddr();
 	i3d::GRAY16* pM=mask.GetFirstVoxelAddr();
 	i3d::GRAY16* const pL=pM+mask.GetImageSize();
@@ -1196,17 +1138,30 @@ void ActiveMesh::CenterMesh(const Vector3F& newCentre)
 	z-=newCentre.z;
 
 	//shift the centre to point (0,0,0)
+	//cell body
 	for (unsigned int i=0; i < Pos.size(); ++i)
 	{
 		Pos[i].x-=float(x);
 		Pos[i].y-=float(y);
 		Pos[i].z-=float(z);
 	}
-	for (unsigned int i=0; i < fPoints.size(); ++i)
+
+	//filopodia
+	for (unsigned int f=0; f < Ftree.size(); ++f)
 	{
-		fPoints[i].x-=float(x);
-		fPoints[i].y-=float(y);
-		fPoints[i].z-=float(z);
+		mesh_t& filoMesh=Ftree[f];
+		for (unsigned int i=0; i < filoMesh.Pos.size(); ++i)
+		{
+			filoMesh.Pos[i].x-=float(x);
+			filoMesh.Pos[i].y-=float(y);
+			filoMesh.Pos[i].z-=float(z);
+		}
+		for (unsigned int i=0; i < filoMesh.fPoints.size(); ++i)
+		{
+			filoMesh.fPoints[i].x-=float(x);
+			filoMesh.fPoints[i].y-=float(y);
+			filoMesh.fPoints[i].z-=float(z);
+		}
 	}
 }
 
@@ -1231,6 +1186,12 @@ void ActiveMesh::ScaleMesh(const Vector3F& scale)
 			filoMesh.Pos[i].y*=scale.y;
 			filoMesh.Pos[i].z*=scale.z;
 		}
+		for (unsigned int i=0; i < filoMesh.fPoints.size(); ++i)
+		{
+			filoMesh.fPoints[i].x*=scale.x;
+			filoMesh.fPoints[i].y*=scale.y;
+			filoMesh.fPoints[i].z*=scale.z;
+		}
 	}
 }
 
@@ -1255,6 +1216,12 @@ void ActiveMesh::TranslateMesh(const Vector3F& shift)
 			filoMesh.Pos[i].y+=shift.y;
 			filoMesh.Pos[i].z+=shift.z;
 		}
+		for (unsigned int i=0; i < filoMesh.fPoints.size(); ++i)
+		{
+			filoMesh.fPoints[i].x+=shift.x;
+			filoMesh.fPoints[i].y+=shift.y;
+			filoMesh.fPoints[i].z+=shift.z;
+		}
 	}
 }
 
@@ -1564,28 +1531,56 @@ void ActiveMesh::InitDots_body(const i3d::Image3d<i3d::GRAY16>& mask)
 //#endif
 }
 
-void ActiveMesh::InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask)
+void ActiveMesh::InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask,const int idx)
 {
-	//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;
+	//syntactic short-cut to the proper filopodium geometry data
+	mesh_t& filoMesh=Ftree[idx];
 
-	const i3d::GRAY16* p=mask.GetFirstVoxelAddr();
-	for (int z=0; z < (signed)mask.GetSizeZ(); ++z)
-	for (int y=0; y < (signed)mask.GetSizeY(); ++y)
-	for (int x=0; x < (signed)mask.GetSizeX(); ++x, ++p)
-	if (*p > 100)
+	//render the filopodium texture points along its skeleton
+	//the amount of points is determined from the radius, from circle area -> square of radius
+	for (size_t i=0; i < filoMesh.segFromPoint.size(); ++i)
 	{
-		//convert px coords into um
-		const float X=(float)x/xRes + xOff;
-		const float Y=(float)y/yRes + yOff;
-		const float Z=(float)z/zRes + zOff;
+		//time-savers...
+		const Vector3F& fP=filoMesh.fPoints[filoMesh.segFromPoint[i]];
+		const Vector3F& tP=filoMesh.fPoints[filoMesh.segToPoint[i]];
+
+		const float length=(tP-fP).Len(); //full segment has length = 0.005
+
+		//stepping along the axis
+		const int S=(int)ceil(5.f*length/0.005); //"fraction factor"
+		//std::cout << "S=" << S << ", length=" << length << "\n";
+		for (int s=0; s < S; ++s)
+		{
+			//ratios "from" and "to"
+			const float fS=float(S-s)/float(S);
+			const float tS=float(s)/float(S);
 
-		for (int q=0; q < 80; ++q) dots.push_back(Vector3F(X,Y,Z));
+			//target radius
+			const float target_r=(i+1 < filoMesh.segFromPoint.size())? filoMesh.segFromRadius[i+1] : 0.f;
+
+			//current radius
+			float r=fS*filoMesh.segFromRadius[i] + tS*target_r;
+
+			//current position
+			Vector3F p =fS*fP;
+			         p+=tS*tP;
+
+			//insert fl. dots
+			r*=4000.f;
+
+			//std::cout << s << ": r=" << r << ", r*r=" << r*r << "\n";
+
+			//float dotsCount=10.f*r*r;
+			float dotsCount=100.f*r;
+			int qM=(dotsCount > 1300.f)? 1300 : (int)dotsCount;
+
+			for (int q=0; q < qM; ++q)
+			{
+				dots_filo.push_back(p);
+				dots_filo.push_back(p+Vector3F(0.0f,0.0f,+0.1f));
+				dots_filo.push_back(p+Vector3F(0.0f,0.0f,-0.1f));
+			}
+		}
 	}
 }
 
@@ -1593,8 +1588,10 @@ void ActiveMesh::InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask)
 void ActiveMesh::BrownDots(const i3d::Image3d<i3d::GRAY16>& mask)
 {
 	//TODO REMOVE ME
+	/*
 	for (size_t i=0; i < dots.size(); ++i)
 		dots[i].x+=1.0f;
+	*/
 }
 
 
diff --git a/cmath3d/TriangleMesh.h b/cmath3d/TriangleMesh.h
index a466e4b..1741ae0 100644
--- a/cmath3d/TriangleMesh.h
+++ b/cmath3d/TriangleMesh.h
@@ -10,6 +10,12 @@ struct mesh_t {
 	std::vector<long unsigned int> ID; //list of triangles organized in tripplets of indices into the Pos array of vertices;
                                       //the array length should be 3x the number of triangles in the mesh
 	std::vector<Vector3F> norm;        //list of vectors that are normal/orthogonal to the planes given by the individual triangles
+
+	///metadata to keep the skeleton information
+	std::vector<Vector3FC> fPoints;    //skel. vertex coordinates
+	std::vector<int> segFromPoint;     //indices of starting points of all segments
+	std::vector<int> segToPoint;       //...of ending points
+	std::vector<float> segFromRadius;  //radii at starting points for all segments
 };
 
 /**
@@ -125,13 +131,6 @@ class ActiveMesh
 	int ImportVTK_Volumetric(const char* filename,bool saveAlsoTetrahedra=false);
 
 	int ImportVTK_Ftree(const char* filename,const int idx,const float stretch=1.f,bool resetMesh=true);
-	std::vector<Vector3FC> fPoints;
-	std::vector<int> segFromPoint; //segment description, quite ugly solution for now
-	std::vector<int> segToPoint;
-	std::vector<float> segFromRadius;
-	std::vector<float> segToRadius;
-	size_t PointsFirstOffset;
-
 	///naive, but it works
 	//renders first aside, then overlays the given mask
 	void RenderMask_body(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles=false);
@@ -146,7 +145,7 @@ class ActiveMesh
 	///time-lapse (coherent) manipulation with fluorescent dots,
 	///always with the given mask
 	void InitDots_body(const i3d::Image3d<i3d::GRAY16>& mask);
-	void InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask); //TODO: if it would work, keep skeleton of every filopodium and introduce points accordingly
+	void InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask,const int idx);
 
 	void BrownDots(const i3d::Image3d<i3d::GRAY16>& mask);
 	void FFDots(const i3d::Image3d<i3d::GRAY16>& mask,
diff --git a/src/graphics.cpp b/src/graphics.cpp
index 684fdf6..6ba1df6 100644
--- a/src/graphics.cpp
+++ b/src/graphics.cpp
@@ -356,32 +356,28 @@ void ActiveMesh::displayMesh(void)
 		glEnd();
 	}
 
-/*
-	//show line segments
-	glLineWidth(1);
+/**/
+	//show filopodia line segments
+	glLineWidth(4);
 	glBegin(GL_LINES);
-	for (unsigned int i=0; i < segFromPoint.size(); ++i)
+	for (size_t f=0; f < Ftree.size(); ++f)
 	{
-		glVertex3f(fPoints[segFromPoint[i]].x,
-		           fPoints[segFromPoint[i]].y,
-		           fPoints[segFromPoint[i]].z);
-		glVertex3f(fPoints[segToPoint[i]].x,
-		           fPoints[segToPoint[i]].y,
-		           fPoints[segToPoint[i]].z);
+		mesh_t& filoMesh=Ftree[f];
+		for (unsigned int i=0; i < filoMesh.segFromPoint.size(); ++i)
+		{
+			glVertex3f(filoMesh.fPoints[filoMesh.segFromPoint[i]].x,
+						  filoMesh.fPoints[filoMesh.segFromPoint[i]].y,
+						  filoMesh.fPoints[filoMesh.segFromPoint[i]].z);
+			glVertex3f(filoMesh.fPoints[filoMesh.segToPoint[i]].x,
+						  filoMesh.fPoints[filoMesh.segToPoint[i]].y,
+						  filoMesh.fPoints[filoMesh.segToPoint[i]].z);
+		}
 	}
 	glEnd();
+	glLineWidth(1);
+/**/
 
-	//show points after rotation
-	glPointSize(3);
-	glBegin(GL_POINTS);
-	for (unsigned int i=0; i < 10*segFromPoint.size()+1; ++i)
-		glVertex3f(Pos[PointsFirstOffset+i].x,
-		           Pos[PointsFirstOffset+i].y,
-		           Pos[PointsFirstOffset+i].z);
-	glEnd();
-*/
-
-	//display filopodias
+	//display filopodia
 	for (size_t f=0; f < Ftree.size(); ++f)
 	{
 		int color=f%2;
diff --git a/src/main.cpp b/src/main.cpp
index dca2504..adab43e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -192,8 +192,6 @@ int LoadNewMesh(int fileNo)
 		return(1);
 	}
 
-	mesh.dots_filo.clear();
-	mesh.dots_filo.reserve(300000);
 	retval=mesh.ImportVTK_Ftree(VTK_fTreeA,0);
 	if (retval) 
 	{
@@ -244,7 +242,12 @@ int RenderMesh(int fileNo)
 	{
 		std::cout << "rendering texture first run\n";
 		mesh.InitDots_body(mask);
-		//mesh.InitDots_filo(mask); -- already re-created with ImportVTK_Ftree()
+
+		//second, create texture for filopodia
+		mesh.dots_filo.clear();
+		mesh.dots_filo.reserve(300000);
+		mesh.InitDots_filo(mask,0); //filopodium one by one
+		mesh.InitDots_filo(mask,1);
 		std::cout << "filo dots size=" << mesh.dots_filo.size() << "\n";
 	}
 	else
@@ -290,7 +293,14 @@ int RenderMesh(int fileNo)
 		//apply the FF on the fl. dots
 		mesh.FFDots(mask,FF);
 
-		//FF destructor will delete dyn. allocated images...
+		//FF destructor will delete dyn. allocated images in the FF...
+
+		//also update filopodia (by re-creating them from the scratch)
+		mesh.dots_filo.clear();
+		mesh.dots_filo.reserve(300000);
+		mesh.InitDots_filo(mask,0); //filopodium one by one
+		mesh.InitDots_filo(mask,1);
+		std::cout << "filo dots size=" << mesh.dots_filo.size() << "\n";
 	}
 
 	//renders the isotropic texture
-- 
GitLab