Added ActiveMesh::VolID -- to calculate more accurate FF later.

Added one new param to ActiveMesh::ImportVTK_Volumetric()

Gaussian-smoothing-based FF creation finished -- ActiveMesh::ConstructFF()

cmath3d/TriangleMesh.cpp

@@ -315,12 +315,14 @@ int ActiveMesh::ImportVTK(const char *filename) //surface version
 }
 
 
-int ActiveMesh::ImportVTK_Volumetric(const char *filename)
+int ActiveMesh::ImportVTK_Volumetric(const char *filename,bool saveAlsoTetrahedra)
 {
 	Pos.clear();
 	ID.clear();
 	norm.clear();
 
+	if (saveAlsoTetrahedra) VolID.clear();
+
 	//try to open the file
 	std::ifstream file(filename);
 	if (!file.is_open()) return 1;
@@ -425,6 +427,14 @@ int ActiveMesh::ImportVTK_Volumetric(const char *filename)
 		ID.push_back(v3);
 		norm.push_back(fictiveNormal);
 
+		if (saveAlsoTetrahedra)
+		{
+			VolID.push_back(v1);
+			VolID.push_back(v2);
+			VolID.push_back(v3);
+			VolID.push_back(v4);
+		}
+
 		--itemCount;
 	}
 	std::cout << "last polyhedron was: " << v1 << "," << v2 << "," << v3 << "," << v4 << "\n";
@@ -1571,43 +1581,33 @@ void ActiveMesh::ConstructFF(FlowField<float> &FF)
 	float* const ffy=FF.y->GetFirstVoxelAddr();
 	float* const ffz=FF.z->GetFirstVoxelAddr();
 
-	//iterations:
-	for (int iters=0; iters < 15; ++iters)
+	//inject displacements
+	for (size_t i=0; i < oldVolPos.size(); ++i)
 	{
-		std::cout << "ConstructFF(): iteration=" << iters << "\n";
+		const int x=(int)roundf( (oldVolPos[i].x-xOff) *xRes);
+		const int y=(int)roundf( (oldVolPos[i].y-yOff) *yRes);
+		const int z=(int)roundf( (oldVolPos[i].z-zOff) *zRes);
 
-		//inject displacements
-		for (size_t i=0; i < oldVolPos.size(); ++i)
-		{
-			const int x=(int)roundf( (oldVolPos[i].x-xOff) *xRes);
-			const int y=(int)roundf( (oldVolPos[i].y-yOff) *yRes);
-			const int z=(int)roundf( (oldVolPos[i].z-zOff) *zRes);
-
-			const float dx=newVolPos[i].x - oldVolPos[i].x;
-			const float dy=newVolPos[i].y - oldVolPos[i].y;
-			const float dz=newVolPos[i].z - oldVolPos[i].z;
+		const float dx=newVolPos[i].x - oldVolPos[i].x;
+		const float dy=newVolPos[i].y - oldVolPos[i].y;
+		const float dz=newVolPos[i].z - oldVolPos[i].z;
 
-			if ((x > 0) && (y > 0) && (z > 0)
-				&& (x < maxX) && (y < maxY) && (z < maxZ))
-			{
-				ffx[z*Slice +y*xLine +x]=dx;
-				ffy[z*Slice +y*xLine +x]=dy;
-				ffz[z*Slice +y*xLine +x]=dz;
-			}
+		if ((x > 0) && (y > 0) && (z > 0)
+			&& (x < maxX) && (y < maxY) && (z < maxZ))
+		{
+			ffx[z*Slice +y*xLine +x]=dx;
+			ffy[z*Slice +y*xLine +x]=dy;
+			ffz[z*Slice +y*xLine +x]=dz;
 		}
-
-		//DEBUG before smoothing
-		char n[1024];
-		sprintf(n,"Vx_ConstructFF_it%03d.ics",iters);
-		FF.x->SaveImage(n);
-		sprintf(n,"Vy_ConstructFF_it%03d.ics",iters);
-		FF.y->SaveImage(n);
-		sprintf(n,"Vz_ConstructFF_it%03d.ics",iters);
-		FF.z->SaveImage(n);
-
-		//smooth
-		i3d::GaussFIR(*FF.x,1.0f);
-		i3d::GaussFIR(*FF.y,1.0f);
-		i3d::GaussFIR(*FF.z,1.0f);
 	}
+
+	//smooth
+	i3d::GaussIIR(*FF.x,15.0f);
+	i3d::GaussIIR(*FF.y,15.0f);
+	i3d::GaussIIR(*FF.z,15.0f);
+
+	//multiply (to "correct" after normalized smoothing)
+	FF.x->GetVoxelData()*=1240.f;
+	FF.y->GetVoxelData()*=1240.f;
+	FF.z->GetVoxelData()*=1240.f;
 }

cmath3d/TriangleMesh.h

@@ -27,12 +27,15 @@ class ActiveMesh
 
 	std::vector<Vector3FC> oldVolPos;  //list of volumetric mesh vertices, previous situation
 	std::vector<Vector3FC> newVolPos;  //list of volumetric mesh vertices, new situation
+	std::vector<long unsigned int> VolID; //list of tetrahedra, filled possibly with this->ImportVTK_Volumetric()
 	void RotateVolPos(void)
 	{
 		oldVolPos=newVolPos;
 		newVolPos=Pos;
 	}
 	void ConstructFF(FlowField<float> &FF);
+	size_t GetVolIDsize(void)
+	{ return (VolID.size()); }
 
 	/**
 	 * Determine coeffs that represent the local surface around vertex \e vertexID.
@@ -106,7 +109,7 @@ class ActiveMesh
 	int ImportSTL(const char* filename);
 	int ExportSTL(const char* filename);
 	int ImportVTK(const char* filename);
-	int ImportVTK_Volumetric(const char* filename);
+	int ImportVTK_Volumetric(const char* filename,bool saveAlsoTetrahedra=false);
 
 	int ImportVTK_Ftree(const char* filename,bool resetMesh=false);
 	std::vector<Vector3FC> fPoints;

src/main.cpp

@@ -106,8 +106,12 @@ int LoadNewMesh(int fileNo)
 	sprintf(VTK_fTreeA,     "../sample_sequence/ID05_t%03d_fTree01.vtk",fileNo);
 	sprintf(VTK_fTreeB,     "../sample_sequence/ID05_t%03d_fTree02.vtk",fileNo);
 
-	//first, obtain a fresh volumetric mesh (before it gets overwritten in the following code)
-	int retval=mesh.ImportVTK_Volumetric(VTK_nucleus);
+	//first, obtain a fresh volumetric mesh of the cell body
+	//(before it gets overwritten in the following code)
+	//
+	//keep the volumetric configuration of the first load
+	bool saveAlsoTetrahedra=(mesh.GetVolIDsize() == 0);
+	int retval=mesh.ImportVTK_Volumetric(VTK_nucleus,saveAlsoTetrahedra);
 	if (retval) 
 	{
 		std::cout << "some error reading nucleus: " << retval << "\n";