Merge branch 'ISBI2016' of ssh://xulman@helenos.fi.muni.cz:2200/home/xulman/GIT/meshSurface.git

.gitignore

 BUILD
+sample*

cmath3d/TriangleMesh.cpp

@@ -3,8 +3,16 @@
 #include <map>
 #include <lapacke.h>
 
+#include <i3d/draw.h>
+#include <i3d/morphology.h>
+#include <i3d/DistanceTransform.h>
+#include <i3d/filters.h>
+
 #include "TriangleMesh.h"
 #include "../cmath3d_v/TriangleMesh_v.h"
+#include "../src/rnd_generators.h"
+#undef min
+#undef max
 
 //'multiple' should be ideally 10^desired_decimal_accuracy
 int inline RoundTo(const float val, const float multiple=1000.f)
@@ -174,7 +182,112 @@ int ActiveMesh::ImportSTL(const char *filename)
 }
 
 
-int ActiveMesh::ImportVTK(const char *filename)
+int ActiveMesh::ImportVTK(const char *filename) //surface version
+{
+	Pos.clear();
+	ID.clear();
+	norm.clear();
+
+	//try to open the file
+	std::ifstream file(filename);
+	if (!file.is_open()) return 1;
+
+	//read the "header" line
+	char tmp[1024];
+	file >> tmp >> tmp; //dangerous...
+	//check tmp for "vtk" or complain
+	if (tmp[0] != 'v'
+	 || tmp[1] != 't'
+	 || tmp[2] != 'k') { file.close(); return(2); }
+	//read (and skip) the rest of the header line
+	file.ignore(10240,'\n');
+
+	//ignore "vtk output"
+	file.ignore(10240,'\n');
+
+	//read "ASCII"
+	file >> tmp;
+	if (tmp[0] != 'A'
+	 || tmp[1] != 'S'
+	 || tmp[2] != 'C'
+	 || tmp[3] != 'I'
+	 || tmp[4] != 'I') { file.close(); return(3); }
+	file.ignore(10240,'\n');
+
+	//ignore "DATASET POLYDATA"
+	file.ignore(10240,'\n');
+
+	//read points header
+	int itemCount;
+	file >> tmp >> itemCount;;
+	if (tmp[0] != 'P'
+	 || tmp[1] != 'O'
+	 || tmp[2] != 'I'
+	 || tmp[3] != 'N'
+	 || tmp[4] != 'T'
+	 || tmp[5] != 'S') { file.close(); return(4); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " point coordinates\n";
+	Pos.reserve(itemCount);
+
+	//read all points...
+	float x,y,z;
+	while (itemCount > 0 && file >> x)
+	{
+		file >> y >> z;
+
+		//... and save them
+		Vector3FC v1(x,y,z);
+		v1*=100.f;
+		Pos.push_back(v1);
+
+		--itemCount;
+	}
+	std::cout << "last coordinate was: " << x << "," << y << "," << z << "\n";
+
+	//prepare "information about faces normals"
+	Vector3F fictiveNormal(1.f,0.f,0.f);
+
+	//read polyhedra header
+	file >> tmp >> itemCount;
+	if (tmp[0] != 'P'
+	 || tmp[1] != 'O'
+	 || tmp[2] != 'L'
+	 || tmp[3] != 'Y'
+	 || tmp[4] != 'G'
+	 || tmp[5] != 'O'
+	 || tmp[6] != 'N'
+	 || tmp[7] != 'S') { file.close(); return(5); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " triangles\n";
+	ID.reserve(3*itemCount);
+	norm.reserve(itemCount);
+
+	//read all polyhedra vertices
+	int ignore,v1,v2,v3;
+	while (itemCount > 0 && file >> ignore && ignore == 3)
+	{
+		file >> v1 >> v2 >> v3;
+
+		//save v1,v2,v3 (TODO: if not already saved...)
+		//make triangles use CW winding order
+		ID.push_back(v1);
+		ID.push_back(v3);
+		ID.push_back(v2);
+		norm.push_back(fictiveNormal);
+
+		--itemCount;
+	}
+	std::cout << "last triangle was: " << v1 << "," << v2 << "," << v3 << "\n";
+
+	file.close();
+	return(0);
+}
+
+
+int ActiveMesh::ImportVTK_Volumetric(const char *filename)
 {
 	Pos.clear();
 	ID.clear();
@@ -184,15 +297,368 @@ int ActiveMesh::ImportVTK(const char *filename)
 	std::ifstream file(filename);
 	if (!file.is_open()) return 1;
 
+	//read the "header" line
+	char tmp[1024];
+	file >> tmp >> tmp; //dangerous...
+	//check tmp for "vtk" or complain
+	if (tmp[0] != 'v'
+	 || tmp[1] != 't'
+	 || tmp[2] != 'k') { file.close(); return(2); }
+	//read (and skip) the rest of the header line
+	file.ignore(10240,'\n');
+
+	//ignore "vtk output"
+	file.ignore(10240,'\n');
+
+	//read "ASCII"
+	file >> tmp;
+	if (tmp[0] != 'A'
+	 || tmp[1] != 'S'
+	 || tmp[2] != 'C'
+	 || tmp[3] != 'I'
+	 || tmp[4] != 'I') { file.close(); return(3); }
+	file.ignore(10240,'\n');
+
+	//ignore "DATASET UNSTRUCTURED GRID"
+	file.ignore(10240,'\n');
+
+	//read points header
+	int itemCount;
+	file >> tmp >> itemCount;;
+	if (tmp[0] != 'P'
+	 || tmp[1] != 'O'
+	 || tmp[2] != 'I'
+	 || tmp[3] != 'N'
+	 || tmp[4] != 'T'
+	 || tmp[5] != 'S') { file.close(); return(4); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " point coordinates\n";
+	Pos.reserve(itemCount);
+
+	//read all points...
+	float x,y,z;
+	while (itemCount > 0 && file >> x)
+	{
+		file >> y >> z;
+
+		//... and save them
+		Vector3FC v1(x,y,z);
+		v1*=100.f;
+		Pos.push_back(v1);
+
+		--itemCount;
+	}
+	std::cout << "last coordinate was: " << x << "," << y << "," << z << "\n";
+
+	//prepare "information about faces normals"
+	Vector3F fictiveNormal(1.f,0.f,0.f);
+
+	//read polyhedra header
+	file >> tmp >> itemCount;
+	if (tmp[0] != 'C'
+	 || tmp[1] != 'E'
+	 || tmp[2] != 'L'
+	 || tmp[3] != 'L'
+	 || tmp[4] != 'S') { file.close(); return(5); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " polyhedra\n";
+	ID.reserve(3*itemCount);
+	norm.reserve(itemCount);
+
+	//read all polyhedra vertices
+	int ignore,v1,v2,v3,v4;
+	while (itemCount > 0 && file >> ignore && ignore == 4)
+	{
+		file >> v1 >> v2 >> v3 >> v4;
+
+		//save v1,v2,v3 (TODO: if not already saved...)
+		ID.push_back(v1);
+		ID.push_back(v2);
+		ID.push_back(v3);
+		norm.push_back(fictiveNormal);
+
+		//save v1,v2,v4 (TODO: if not already saved...)
+		ID.push_back(v1);
+		ID.push_back(v2);
+		ID.push_back(v4);
+		norm.push_back(fictiveNormal);
+
+		//save v1,v4,v3 (TODO: if not already saved...)
+		ID.push_back(v1);
+		ID.push_back(v4);
+		ID.push_back(v3);
+		norm.push_back(fictiveNormal);
+
+		//save v4,v2,v3 (TODO: if not already saved...)
+		ID.push_back(v4);
+		ID.push_back(v2);
+		ID.push_back(v3);
+		norm.push_back(fictiveNormal);
+
+		--itemCount;
+	}
+	std::cout << "last polyhedron was: " << v1 << "," << v2 << "," << v3 << "," << v4 << "\n";
+
+	file.close();
+	return(0);
+}
+
+
+int ActiveMesh::ImportVTK_Ftree(const char *filename,bool resetMesh)
+{
+	const int PointsOnRadiusPeriphery=10;
+	const float radiusCorrection=0.1f;
+
+	if (resetMesh)
+	{
+		Pos.clear();
+		ID.clear();
+		norm.clear();
+	}
+	
+	//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();
+
+	//try to open the file
+	std::ifstream file(filename);
+	if (!file.is_open()) return 1;
+
+	//read the "header" line
+	char tmp[1024];
+	file >> tmp >> tmp; //dangerous...
+	//check tmp for "vtk" or complain
+	if (tmp[0] != 'v'
+	 || tmp[1] != 't'
+	 || tmp[2] != 'k') { file.close(); return(2); }
+	//read (and skip) the rest of the header line
+	file.ignore(10240,'\n');
+
+	//ignore "vtk output"
+	file.ignore(10240,'\n');
+
+	//read "ASCII"
+	file >> tmp;
+	if (tmp[0] != 'A'
+	 || tmp[1] != 'S'
+	 || tmp[2] != 'C'
+	 || tmp[3] != 'I'
+	 || tmp[4] != 'I') { file.close(); return(3); }
+	file.ignore(10240,'\n');
+
+	//ignore "DATASET UNSTRUCTURED GRID"
+	file.ignore(10240,'\n');
+
+	//read points header
+	int itemCount;
+	file >> tmp >> itemCount;;
+	if (tmp[0] != 'P'
+	 || tmp[1] != 'O'
+	 || tmp[2] != 'I'
+	 || tmp[3] != 'N'
+	 || tmp[4] != 'T'
+	 || tmp[5] != 'S') { file.close(); return(4); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " point coordinates\n";
+	fPoints.reserve(itemCount);
+
+	//read all points...
+	float x,y,z;
+	while (itemCount > 0 && file >> x)
+	{
+		file >> y >> z;
+
+		//... and save them
+		Vector3FC v1(x,y,z);
+		v1*=100.f;
+		fPoints.push_back(v1);
+
+		--itemCount;
+	}
+	std::cout << "last coordinate was: " << x << "," << y << "," << z << "\n";
+
+	//read segments header
+	file >> tmp >> itemCount;
+	if (tmp[0] != 'C'
+	 || tmp[1] != 'E'
+	 || tmp[2] != 'L'
+	 || tmp[3] != 'L'
+	 || tmp[4] != 'S') { file.close(); return(5); }
+	file.ignore(10240,'\n');
+
+	std::cout << "reading " << itemCount << " segments\n";
+	segFromPoint.reserve(itemCount);
+	  segToPoint.reserve(itemCount);
+	segFromRadius.reserve(itemCount);
+
+	//read all segments vertices
+	int ignore,v1,v2;
+	while (itemCount > 0 && file >> ignore && ignore == 2)
+	{
+		file >> v1 >> v2;
+
+		segFromPoint.push_back(v1);
+		  segToPoint.push_back(v2);
+
+		--itemCount;
+	}
+	std::cout << "last segment was: " << v1 << "," << v2 << "\n";
+
+	//"ignore" cell types header, body 
+	file >> tmp >> itemCount;
+	if (tmp[0] != 'C'
+	 || tmp[1] != 'E'
+	 || tmp[2] != 'L'
+	 || tmp[3] != 'L'
+	 || tmp[4] != '_'
+	 || tmp[5] != 'T'
+	 || tmp[6] != 'Y') { file.close(); return(6); }
+	file.ignore(10240,'\n');
+	while (itemCount > 0)
+	{
+		file.ignore(10240,'\n');
+		--itemCount;
+	}
 
+	//read radii, i.e. CELL_DATA header
+	file >> tmp >> itemCount;
+	if (tmp[0] != 'C'
+	 || tmp[1] != 'E'
+	 || tmp[2] != 'L'
+	 || tmp[3] != 'L'
+	 || tmp[4] != '_'
+	 || tmp[5] != 'D'
+	 || tmp[6] != 'A') { file.close(); return(7); }
+	file.ignore(10240,'\n');
+	file.ignore(10240,'\n'); //ignore SCALARS..
+	file.ignore(10240,'\n'); //ignore LOOKUP_TABLE
 
+	//read the radii
+	v1=0;
+	while (itemCount > 0 && file >> x)
+	{
+		segFromRadius.push_back(x*radiusCorrection);
 
+		--itemCount;
+		++v1;
+	}
 	file.close();
+
+	//now widen segments and save triangles...
+
+	//prepare "information about faces normals"
+	Vector3F fictiveNormal(1.f,0.f,0.f);
+
+	size_t firstRadiusPoint=Pos.size();
+	PointsFirstOffset=firstRadiusPoint;
+
+	//for all bending points except the last one (tip)
+	for (unsigned int i=0; i < segFromPoint.size(); ++i)
+	{
+		//we need to construct rotation matrix that would
+		//rotate points on a circle which lays in the XZ plane
+		//(Y axis is normal to it then) such that the points
+		//would lay in the plane to which this new Y axis
+		//would be normal
+		//
+		//new Y axis 
+		Vector3F nYaxis=fPoints[segToPoint[i]];
+		if (i == 0) nYaxis-=fPoints[segFromPoint[i]];
+		else
+		{
+			nYaxis-=fPoints[segFromPoint[i-1]];
+			nYaxis/=2.f; //unnecessary scaling
+		}
+
+		Vector3F nZaxis(0.f,1.f,0.f); //in fact it is original Yaxis _for now_
+		Vector3F nXaxis;
+		//new X axis is perpendicular to the original and new Y axis
+		Mul(nYaxis,nZaxis,nXaxis);
+
+		//new Z axis is perpendicular to the new X and Y axes
+		Mul(nYaxis,nXaxis,nZaxis);
+
+		//normalize...
+		nXaxis/=nXaxis.Len();
+		nZaxis/=nZaxis.Len();
+
+		//now render the points on the circle and project them into the scene
+		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));
+
+			//rotate
+			Vector3FC V(x*nXaxis);
+			V+=z*nZaxis;
+
+			//shift to the centre and save
+			V+=fPoints[segFromPoint[i]];
+			Pos.push_back(V);
+		}
+	}
+
+	//finally, add the tip point
+	Pos.push_back(fPoints[segToPoint.back()]);
+
+	//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)
+	{
+		int p=0;
+		for (; p < PointsOnRadiusPeriphery-1; ++p)
+		{
+			ID.push_back(firstRadiusPoint+p);
+			ID.push_back(firstRadiusPoint+p+1);
+			ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
+			norm.push_back(fictiveNormal);
+
+			ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery+1);
+			ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
+			ID.push_back(firstRadiusPoint+p+1);
+			norm.push_back(fictiveNormal);
+		}
+
+		ID.push_back(firstRadiusPoint+p);
+		ID.push_back(firstRadiusPoint);
+		ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
+		norm.push_back(fictiveNormal);
+
+		ID.push_back(firstRadiusPoint  +PointsOnRadiusPeriphery);
+		ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
+		ID.push_back(firstRadiusPoint);
+		norm.push_back(fictiveNormal);
+
+		firstRadiusPoint+=PointsOnRadiusPeriphery;
+	}
+
+	//the last segment...
+	int p=0;
+	for (; p < PointsOnRadiusPeriphery-1; ++p)
+	{
+		ID.push_back(firstRadiusPoint+p);
+		ID.push_back(firstRadiusPoint+p+1);
+		ID.push_back(firstRadiusPoint+PointsOnRadiusPeriphery);
+		norm.push_back(fictiveNormal);
+	}
+
+	ID.push_back(firstRadiusPoint+p);
+	ID.push_back(firstRadiusPoint);
+	ID.push_back(firstRadiusPoint+PointsOnRadiusPeriphery);
+	norm.push_back(fictiveNormal);
+
 	return(0);
 }
 
 
-void ActiveMesh::RenderMask(i3d::Image3d<i3d::GRAY16>& mask)
+void ActiveMesh::RenderMask(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles)
 {
 	//time savers: resolution
 	const float xRes=mask.GetResolution().GetX();
@@ -229,9 +695,25 @@ void ActiveMesh::RenderMask(i3d::Image3d<i3d::GRAY16>& mask)
 			const int y=(int)roundf( (v.y-yOff) *yRes);
 			const int z=(int)roundf( (v.z-zOff) *zRes);
 
-			if (mask.Include(x,y,z)) mask.SetVoxel(x,y,z,short(i));
+			short val=(showTriangles)? short(i%5 *30 +100) : 100;
+			if (mask.Include(x,y,z)) mask.SetVoxel(x,y,z,val);
 		}
 	}
+
+	//this floods cell exterior from the image corner
+	i3d::Image3d<i3d::GRAY16> tmpMask(mask);
+	i3d::Dilation(tmpMask,mask,i3d::nb3D_o18);
+	i3d::FloodFill(mask,(i3d::GRAY16)50,0);
+
+	//this removes the flooding while filling
+	//everything else (and closing holes in this was)
+	i3d::GRAY16* p=mask.GetFirstVoxelAddr();
+	i3d::GRAY16* const pL=p+mask.GetImageSize();
+	while (p != pL)
+	{
+		*p=(*p != 50)? std::max(*p,i3d::GRAY16(100)) : 0;
+		++p;
+	}
 }
 
 
@@ -356,6 +838,63 @@ void ActiveMesh::RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask)
 }
 
 
+void ActiveMesh::RenderOneTimeTexture(const i3d::Image3d<i3d::GRAY16>& mask,
+                                      i3d::Image3d<i3d::GRAY16>& texture)
+{
+	i3d::Image3d<float> dt;
+	i3d::GrayToFloat(mask,dt);
+	i3d::EDM(dt,0,100.f,false);
+
+	dt.SaveImage("dt.ics");
+
+	float* p=dt.GetFirstVoxelAddr();
+	float* const pL=p+dt.GetImageSize();
+	int randomCounter=20;
+	while (p != pL)
+	{
+		//are we on a surface shell?
+		if (*p < 0.10f)
+		{
+			//invert the intensity levels (iso-lines)
+			if (*p > 0.f) *p=(0.12f-*p)*100.f;
+
+			//introduce random bumbs?
+			if (randomCounter==0)
+			{
+				*p *= 2.f;
+				randomCounter=30; //something random here...
+			}
+			--randomCounter;
+		}
+		else *p=0.f; //else clear the pixel value
+		++p;
+	}
+
+	dt.SaveImage("dt2.ics");
+	i3d::GaussIIR(dt,3.f,3.f,9.f);
+	p=dt.GetFirstVoxelAddr();
+	for (; p != pL; ++p)
+	{
+		//uncertainty in the number of incoming photons
+		float noiseMean = sqrt(*p), // from statistics: shot noise = sqrt(signal) 
+				noiseVar = noiseMean; // for Poisson distribution E(X) = D(X)
+
+		*p=ceilf(*p + GetRandomPoisson(noiseMean) - noiseVar);
+
+		//constants are parameters of Andor iXon camera provided from vendor:
+		//photon shot noise: dark current
+		*p+=GetRandomPoisson(0.06f);
+
+		//read-out noise:
+		//  variance up to 25.f (old camera on ILBIT)
+		//  variance about 1.f (for the new camera on ILBIT)
+		//*p+=GetRandomGauss(700.f,90.f);
+	}
+
+	i3d::FloatToGrayNoWeight(dt,texture);
+}
+
+
 void ActiveMesh::CenterMesh(const Vector3F& newCentre)
 {
 	//calc geom. centre
@@ -381,6 +920,12 @@ void ActiveMesh::CenterMesh(const Vector3F& newCentre)
 		Pos[i].y-=float(y);
 		Pos[i].z-=float(z);
 	}
+	for (unsigned int i=0; i < fPoints.size(); ++i)
+	{
+		fPoints[i].x-=float(x);
+		fPoints[i].y-=float(y);
+		fPoints[i].z-=float(z);
+	}
 }
 
 

cmath3d/TriangleMesh.h

@@ -97,14 +97,27 @@ class ActiveMesh
 	///returns 0 on success, otherwise non-zero
 	int ImportSTL(const char* filename);
 	int ImportVTK(const char* filename);
-	void RenderMask(i3d::Image3d<i3d::GRAY16>& mask);
+	int ImportVTK_Volumetric(const char* filename);
+
+	int ImportVTK_Ftree(const char* filename,bool resetMesh=false);
+	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;
+
+	void RenderMask(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles=false);
 	void RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask);
+	void RenderOneTimeTexture(const i3d::Image3d<i3d::GRAY16>& mask,
+	                          i3d::Image3d<i3d::GRAY16>& texture);
 
 	void CenterMesh(const Vector3F& newCentre);
 	void ScaleMesh(const Vector3F& scale);
 
 	//the following functions are all defined in graphics.cpp
 	void displayMesh(void);
+	void displayMeshEdges(void);
 	void displayVertexAndNeigs(void);
 	void displayQuadricSurface(void);
 	void displayTEST(void);

src/graphics.cpp

@@ -46,18 +46,26 @@ float zoom = 1.0;
 bool MouseActive = false;
 int last_x = 0, last_y = 0;
 
+bool showAxes=true;
+bool showFrame=false;
+bool showBox=true;
+
 bool showWholeMesh=true;
+bool showWholeMeshFilled=false;
+bool showWholeMeshEdges=false;
 bool showVertexSurface=false;
 bool showQuadricSurface=false;
 bool showTEST=false;
 
+int LoadNewMesh(int fileNo);
+int RenderMesh(int fileNo);
 
 ///manage \e count steps of the simulation and then print and draw
 void step(int count)
 {
 	while ((params.currTime < params.stopTime) && (count > 0))
 	{
-		//do some job
+		LoadNewMesh(int(params.currTime));
 
 		params.currTime += params.incrTime;
 		--count;
@@ -76,7 +84,7 @@ void step(int count)
 }
 
 
-///draws the frame around the playground
+///draws a flat frame around the playground
 void displayFrame(void)
 {
     glColor3f(params.sceneBorderColour.r,
@@ -91,6 +99,45 @@ void displayFrame(void)
     glEnd();
 }
 
+///draws a frame box around the playground
+void displayBox(void)
+{
+    glColor3f(params.sceneBorderColour.r,
+	 			  params.sceneBorderColour.g,
+	 			  params.sceneBorderColour.b);
+    glLineWidth(1);
+	 //front square/plane
+    glBegin(GL_LINE_LOOP);
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y                      , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y                      , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z );
+    glEnd();
+
+	 //rear square/plane
+    glBegin(GL_LINE_LOOP);
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y                      , params.sceneOffset.z + params.sceneSize.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y                      , params.sceneOffset.z + params.sceneSize.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z );
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z );
+    glEnd();
+
+	 //sideways "connectors"
+	 glBegin(GL_LINES);
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y                      , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y                      , params.sceneOffset.z + params.sceneSize.z );
+
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y                      , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y                      , params.sceneOffset.z + params.sceneSize.z );
+
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x + params.sceneSize.x , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z );
+
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z );
+    glVertex3f( params.sceneOffset.x                      , params.sceneOffset.y + params.sceneSize.y , params.sceneOffset.z + params.sceneSize.z );
+    glEnd();
+}
+
 void displayAxes(void)
 {
 	glLineWidth(2);
@@ -107,7 +154,6 @@ void displayAxes(void)
 	glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z);
 	glVertex3f(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z+params.sceneSize.z/4.f);
 	glEnd();
-	glLineWidth(1);
 }
 
 
@@ -214,6 +260,12 @@ void ActiveMesh::displayTEST(void)
 
 void ActiveMesh::displayQuadricSurface(void)
 {
+
+	//zkusit nad dalsimi krivkami, paraboloid?
+	//pripadne udelat jinak dopocitavani na povrch..
+	//bo to nyni dava docela velke vzdalenosti pro vypocitany povrch
+
+
 	//get surface params
 	float surf_coeff[10];
 	//CalcQuadricSurface_Taubin(VertexID,surf_coeff);
@@ -289,25 +341,68 @@ void ActiveMesh::displayQuadricSurface(void)
 
 void ActiveMesh::displayMesh(void)
 {
-	glColor3f(0.6f,0.6f,0.6f);
+	const float gray=0.6f;
+	glColor4f(gray,gray,gray,1.0f);
+	glLineWidth(1);
 
 	long unsigned int* id=ID.data();
 	for (unsigned int i=0; i < ID.size(); i+=3)
 	{
-		//glBegin(GL_TRIANGLES);
-		glBegin(GL_LINE_LOOP);
+		glBegin(GL_TRIANGLES);
 		glNormal3f(norm[i/3].x,norm[i/3].y,norm[i/3].z);
 		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
 		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
 		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
 		glEnd();
 	}
+
+/*
+	//show line segments
+	glLineWidth(1);
+	glBegin(GL_LINES);
+	for (unsigned int i=0; i < segFromPoint.size(); ++i)
+	{
+		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);
+	}
+	glEnd();
+
+	//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();
+*/
+}
+
+void ActiveMesh::displayMeshEdges(void)
+{
+	const float gray=0.8f;
+	glColor4f(gray,gray,gray,1.0f);
+	glLineWidth(3);
+
+	long unsigned int* id=ID.data();
+	for (unsigned int i=0; i < ID.size(); i+=3)
+	{
+		glBegin(GL_LINE_LOOP);
+		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
+		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
+		glVertex3f(Pos[*id].x,Pos[*id].y,Pos[*id].z); ++id;
+		glEnd();
+	}
 }
 
 
 void display(void)
 {
-	glClear(GL_COLOR_BUFFER_BIT);
+	glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
 	glLoadIdentity();
 	glTranslatef(params.sceneCentre.x,params.sceneCentre.y,params.sceneCentre.z);
 	glScalef(zoom,zoom,zoom);
@@ -320,14 +415,19 @@ void display(void)
 	          //0.0, cosf(-uhel_y), sinf(-uhel_y)); // up
 	          0.0, 1.0, 0.0); // up
 
-	displayFrame();
-	displayAxes();
+	glEnable(GL_FOG);
+	if (showBox) displayBox();
+	if (showFrame) displayFrame();
+	glDisable(GL_FOG);
 
 	if (showWholeMesh) mesh.displayMesh();
+	if (showWholeMeshEdges) mesh.displayMeshEdges();
 	if (showVertexSurface) mesh.displayVertexAndNeigs();
 	if (showQuadricSurface) mesh.displayQuadricSurface();
 	if (showTEST) mesh.displayTEST();
 
+	if (showAxes) displayAxes();
+
 	glFlush();
 	glutSwapBuffers();
 }
@@ -376,7 +476,25 @@ void GetSceneViewSize(const int w,const int h,
 void printKeysHelp(void)
 {
 	std::cout << "\nKeys:\n";
-	std::cout << "ESC or 'q': quit\n";
+	std::cout << "ESC or 'q': quit        ENTER: enters empty line on a console\n";
+	std::cout << "' ','n','z': simulation goes/stop at various steps\n";
+	std::cout << "'r': resets the view\n";
+	std::cout << "'o','O': zoom out,in\n";
+	std::cout << "'h': displays this help\n";
+	std::cout << "'d': displays enabled features,displays status\n";
+	std::cout << "'m': toggle display of the mesh\n";
+	std::cout << "'M': toggle display of the mesh as wireframe or filled\n";
+	std::cout << "'e': toggle display of the mesh wireframe (triangle edges)'\n";
+	std::cout << "'v': toggle display of the vertex supporting region\n";
+	std::cout << "'V': chooses another vertex\n";
+	std::cout << "'s': toggle display of the fitted quadratic surface \n";
+	std::cout << "'t': toggle display of the test/debug function\n";
+	std::cout << "'a': toggle display of the R,G,B <-> x,y,z axes \n";
+	std::cout << "'f': toggle display of the frame/plane az z=0\n";
+	std::cout << "'b': toggle display of the bounding box\n";
+	std::cout << "'c','C': turns off,on z-buffer (off = render in the order of drawing)\n";
+	std::cout << "'x','X': controls for some testing\n";
+	std::cout << "'i': not used\n";
 }
 
 
@@ -407,6 +525,9 @@ void keyboard(unsigned char key, int mx, int my)
         case 'q':  //quit
             exit(EXIT_SUCCESS);
             break;
+			case 13: //Enter key -- just empty lines
+				std::cout << std::endl;
+				break;
 
 		  //control execution of the simulation
         case ' ':  //pause/unpause, refresh every 10 cycles
@@ -430,7 +551,7 @@ void keyboard(unsigned char key, int mx, int my)
 				step(1);
             break;
 
-		  //etc stuff...
+		  //view params
         case 'r':  //set scale 1:1
 				uhel_x = 0.0;
 				uhel_y = 0.0;
@@ -438,12 +559,6 @@ void keyboard(unsigned char key, int mx, int my)
 				glutReshapeWindow(int(params.imgSizeX),int(params.imgSizeY));
 				glutPostRedisplay();
             break;
-			case 'h': //print help
-				printKeysHelp();
-				break;
-			case 'd': //print what is displayed
-				printDispayedFeauters();
-				break;
 			case 'o':
 				zoom -= 0.2f;
 				if (zoom < 0.2f) zoom=0.2f;
@@ -454,27 +569,87 @@ void keyboard(unsigned char key, int mx, int my)
 				glutPostRedisplay();
 				break;
 
-			case 'v':
-				++VertexID;
-				glutPostRedisplay();
+			//control hints/help
+			case 'h': //print help
+				printKeysHelp();
+				break;
+			case 'd': //print what is displayed
+				printDispayedFeauters();
 				break;
 
-			case 'M':
+
+			//display objects
+			case 'm':
 				showWholeMesh^=true;
 				glutPostRedisplay();
 				break;
-			case 'V':
+			case 'M':
+				showWholeMeshFilled^=true;
+				if (showWholeMeshFilled)
+					glPolygonMode(GL_FRONT,GL_FILL);
+				else
+					glPolygonMode(GL_FRONT,GL_LINE);
+				glutPostRedisplay();
+				break;
+			case 'e':
+				showWholeMeshEdges^=true;
+				glutPostRedisplay();
+				break;
+			case 'v':
 				showVertexSurface^=true;
 				glutPostRedisplay();
 				break;
-			case 'Q':
+			case 'V':
+				++VertexID;
+				glutPostRedisplay();
+				break;
+			case 's':
 				showQuadricSurface^=true;
 				glutPostRedisplay();
 				break;
-			case 'T':
+			case 't':
 				showTEST^=true;
 				glutPostRedisplay();
 				break;
+			case 'c':
+				glDisable(GL_DEPTH_TEST);
+				std::cout << "Z-buffer is OFF\n";
+				glutPostRedisplay();
+				break;
+			case 'C':
+				glEnable(GL_DEPTH_TEST);
+				std::cout << "Z-buffer is ON\n";
+				glutPostRedisplay();
+				break;
+
+			//annotation controls
+			case 'a':
+				showAxes^=true;
+				glutPostRedisplay();
+				break;
+			case 'f':
+				showFrame^=true;
+				glutPostRedisplay();
+				break;
+			case 'b':
+				showBox^=true;
+				glutPostRedisplay();
+				break;
+
+			//some testing controls
+			case 'x':
+				glFrontFace(GL_CW);
+				glutPostRedisplay();
+				break;
+			case 'X':
+				glFrontFace(GL_CCW);
+				glutPostRedisplay();
+				/*
+				GLfloat ranges[10];
+				glGetFloatv(GL_DEPTH_RANGE,ranges);
+				std::cout << "near=" << ranges[0] << ", far=" << ranges[1] << "\n";
+				*/
+				break;
 
 			case 'i': //inspect a cell
 				GetSceneViewSize(windowSizeX,windowSizeY, xVisF,xVisT,yVisF,yVisT);
@@ -493,9 +668,6 @@ void keyboard(unsigned char key, int mx, int my)
 					}
 			*/
 				break;
-			case 13: //Enter key -- just empty lines
-				std::cout << std::endl;
-				break;
     }
 }
 
@@ -600,7 +772,7 @@ void initializeGL(void)
 
     glutInit(&Argc, Argv);
 
-    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
     glutInitWindowSize(800, 800); 
     glutInitWindowPosition(600, 100);
     glutCreateWindow(Argv[0]);
@@ -614,10 +786,48 @@ void initializeGL(void)
     glutMotionFunc(MouseMotion);
 
     glClearColor(0.0, 0.0, 0.0, 0.0);
+
+	 //enables alfa-blending for every object/drawing
     glEnable(GL_BLEND);
     glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+
+	 //filter rendered points (should look nicer)
     glEnable(GL_POINT_SMOOTH);
+
+	 //default displaying mode for triangles
+	 //frontfaces are just lines, backfaces are opaque
+	 glPolygonMode(GL_BACK, GL_FILL);
+	 if (showWholeMeshFilled)
+		glPolygonMode(GL_FRONT,GL_FILL);
+	 else
+		glPolygonMode(GL_FRONT,GL_LINE);
+
+	 //but, well, do not draw backfaces at all
     glEnable(GL_CULL_FACE);
+	 glCullFace(GL_BACK);
+
+	 //enable z-buffer;
+	 //otherwise it would depend on the order of primitives drawing
+	 glEnable(GL_DEPTH_TEST);
+
+	 //set up the fog (is enabled/disabled before/after drawing frames):
+	 GLfloat fogColor[4]={0.4f*params.sceneBorderColour.r,
+	                      0.4f*params.sceneBorderColour.g,
+	                      0.4f*params.sceneBorderColour.b,
+								 1.0f}; //has no effect
+	 glFogfv(GL_FOG_COLOR,fogColor);
+	 //
+	 //density is not required in linear mode
+	 //glFogf(GL_FOG_DENSITY,0.3f);
+	 glFogi(GL_FOG_MODE,GL_LINEAR);
+	 //
+	 glFogf(GL_FOG_START,params.sceneOffset.z);
+	 glFogf(GL_FOG_END,params.sceneOffset.z+params.sceneSize.z);
+
+	 //some OpenGL info, for fun...
+	 std::cout << "Vendor  : " << glGetString(GL_VENDOR) << "\n";
+	 std::cout << "Renderer: " << glGetString(GL_RENDERER) << "\n";
+	 std::cout << "Version : " << glGetString(GL_VERSION) << "\n";
 }
 
 void loopGL(void)

src/main.cpp

@@ -10,12 +10,16 @@ void ParamsSetup(void);
 //
 ActiveMesh mesh;
 
+int LoadNewMesh(int fileNo);
+int RenderMesh(int fileNo);
 
 int main(void)
 {
 	ParamsSetup();
 
-	//load mesh
+/*
+	//SURFACE FITTING:
+
 	//char filename[]="../sample_input/samplecell.stl";
 	//char filename[]="../sample_input/torus_100_30.stl";
 	char filename[]="../sample_input/sphere.stl";
@@ -35,27 +39,24 @@ int main(void)
 	mesh.CenterMesh(params.sceneCentre);
 
 	//work with mesh
-/*
 	float surf_coeff[10];
 	mesh.CalcQuadricSurface_Taubin(10,surf_coeff);
 */
 
-/*
-	//render mesh
-	i3d::Image3d<i3d::GRAY16> mask;
-	mask.SetOffset(i3d::Offset(0.0,0.0,0.0));
-	//mask.SetResolution(i3d::Resolution(20.0,20.0,20.0)); //for torus
-	//mask.MakeRoom(300,300,300);
-	mask.SetResolution(i3d::Resolution(2.0,2.0,2.0));    //for sample cell
-	mask.MakeRoom(600,600,600);
-	mask.GetVoxelData()=0;
-	mesh.RenderMask(mask);
-	mask.SaveImage("mesh.ics");
-*/
 
+	//ISBI:
+	//for (int i=0; i < 50; ++i)
+	for (int i=10; i < 11; ++i)
+	{
+		LoadNewMesh(i);
+		RenderMesh(i);
+	}
+
+/*
 	initializeGL();
 	loopGL();
 	closeGL();
+*/
 
 	return(0);
 }
@@ -73,9 +74,9 @@ void ParamsSetup(void)
 	params.sceneCentre+=params.sceneOffset;
 
 	params.sceneOuterBorder=Vector3d<float>(2.f);
-	params.sceneBorderColour.r=0.5f;
-	params.sceneBorderColour.g=0.5f;
-	params.sceneBorderColour.b=0.5f;
+	params.sceneBorderColour.r=0.7f;
+	params.sceneBorderColour.g=0.7f;
+	params.sceneBorderColour.b=0.0f;
 
 	params.inputCellsFilename="cells/cell%d.txt";
 	params.numberOfAgents=2;
@@ -86,9 +87,9 @@ void ParamsSetup(void)
 	//adapted for cell cycle length of 24 hours
 	params.cellCycleLength=14.5*60; //[min]
 
-	params.currTime=0.f; //all three are in units of minutes
-	params.incrTime=0.1f;
-	params.stopTime=5000.f;
+	params.currTime=0.0f; //all three are in units of minutes
+	params.incrTime=1.0f;
+	params.stopTime=50.0f;
 
 	params.imgSizeX=500; //pixels
 	params.imgSizeY=500;
@@ -102,3 +103,68 @@ void ParamsSetup(void)
 	params.imgPhCFilename="PhCImg%05d.tif";
 	params.tracksFilename="tracks.txt";
 }
+
+
+int LoadNewMesh(int fileNo)
+{
+	char VTK_nucleusSurf[1024];
+	char VTK_nucleus[1024];
+	char VTK_fTreeA[1024];
+	char VTK_fTreeB[1024];
+
+	sprintf(VTK_nucleusSurf,"../sample_sequence/ID05_t%03d_nucleusSurf.vtk",fileNo);
+	sprintf(VTK_nucleus,    "../sample_sequence/ID05_t%03d_nucleus.vtk",fileNo);
+	sprintf(VTK_fTreeA,     "../sample_sequence/ID05_t%03d_fTree01.vtk",fileNo);
+	sprintf(VTK_fTreeB,     "../sample_sequence/ID05_t%03d_fTree02.vtk",fileNo);
+
+	int retval=mesh.ImportVTK(VTK_nucleusSurf);
+	//int retval=mesh.ImportVTK_Volumetric(VTK_nucleus);
+	if (retval) 
+	{
+		std::cout << "some error reading nucleus: " << retval << "\n";
+		return(1);
+	}
+
+	retval=mesh.ImportVTK_Ftree(VTK_fTreeA);
+	if (retval) 
+	{
+		std::cout << "some error reading f-tree: " << retval << "\n";
+		return(2);
+	}
+
+	retval=mesh.ImportVTK_Ftree(VTK_fTreeB);
+	if (retval) 
+	{
+		std::cout << "some error reading f-tree: " << retval << "\n";
+		return(2);
+	}
+
+	mesh.CenterMesh(params.sceneCentre);
+	std::cout << "loaded mesh #" << fileNo << "\n";
+
+	return(0);
+}
+
+
+int RenderMesh(int fileNo)
+{
+	i3d::Image3d<i3d::GRAY16> mask;
+	mask.SetOffset(i3d::Offset(0.0,0.0,0.0));
+	mask.SetResolution(i3d::Resolution(20.0,20.0,20.0));    //for sample cell
+	mask.MakeRoom(300,300,300);
+	mask.GetVoxelData()=0;
+	mesh.RenderMask(mask);
+
+	char fileName[1024];
+	sprintf(fileName,"mesh_t%03d.tif",fileNo);
+	mask.SaveImage(fileName);
+	std::cout << "rendered mesh #" << fileNo << "\n";
+
+	i3d::Image3d<i3d::GRAY16> texture;
+	mesh.RenderOneTimeTexture(mask,texture);
+	sprintf(fileName,"text_t%03d.tif",fileNo);
+	texture.SaveImage(fileName);
+	std::cout << "textured mesh #" << fileNo << "\n";
+
+	return(0);
+}