TriangleMesh.cpp

#include <iostream>
#include <fstream>
#include <map>

#include <i3d/draw.h>
#include <i3d/morphology.h>
#include <i3d/transform.h>
#include <i3d/DistanceTransform.h>
#include <i3d/filters.h>

#include "TriangleMesh.h"
#include "../cmath3d_v/TriangleMesh_v.h"
#include "../src/rnd_generators.h"
#include "../src/texture/texture.h"

#undef min
#undef max

//some debug-code enabling triggers
//#define SAVE_INTERMEDIATE_IMAGES

//'multiple' should be ideally 10^desired_decimal_accuracy
int inline RoundTo(const float val, const float multiple=1000.f)
{
	return ( int(floorf(val*multiple)) );
}

//puts v1 into Pos, with mPos being a helper structure preventing having
//v1 multiple times inside the Pos
long unsigned int Enlist(
	const Vector3FC& v1,
	std::vector<Vector3FC>& Pos,
	std::map< int,std::map< int,std::map< int,long unsigned int > > >& mPos)
{
	long unsigned int o1; //ret val

	std::map< int,std::map< int,long unsigned int > >& mY=mPos[RoundTo(v1.x)];
	if (mY.empty())
	{
		Pos.push_back(v1);
		o1=Pos.size();

		//add reference to this vertex in the mPos structure
		std::map< int,long unsigned int > mZ;
		mZ[RoundTo(v1.z)]=o1;
		std::map< int,std::map< int,long unsigned int > > my;
		my[RoundTo(v1.y)]=mZ;
		mPos[RoundTo(v1.x)]=my;
	}
	else
	{
		std::map< int,long unsigned int >& mZ=mY[RoundTo(v1.y)];
		if (mZ.empty())
		{
			Pos.push_back(v1);
			o1=Pos.size();

			//add reference to this vertex in the mPos structure
			std::map< int,long unsigned int > mZ;
			mZ[RoundTo(v1.z)]=o1;
			mY[RoundTo(v1.y)]=mZ;
		}
		else
		{
			if (mZ[RoundTo(v1.z)] == 0)
			{
				Pos.push_back(v1);
				o1=Pos.size();

				//add reference to this vertex in the mPos structure
				mZ[RoundTo(v1.z)]=o1;
			}
			else
			{
				o1=mZ[RoundTo(v1.z)];
			}
		}
	}

	return o1;
}


int ActiveMesh::ImportSTL(const char *filename)
{
	Pos.clear();
	ID.clear();
	norm.clear();

	//a helper map to (efficiently) search for already stored vertices inside Pos
	std::map< int,std::map< int,std::map< int,long unsigned int > > > mPos;
	//         x             y             z  offset+1 in Pos

	//try to open the file
	std::ifstream file(filename);
	if (!file.is_open()) return 1;

	//read the "header" line
	char tmp[1024];
	file >> tmp; //dangerous...
	//check tmp for "solid" or complain
	if (tmp[0] != 's'
	 || tmp[1] != 'o'
	 || tmp[2] != 'l'
	 || tmp[3] != 'i'
	 || tmp[4] != 'd') { file.close(); return(2); }
	//read (and skip) the rest of the header line
	file.ignore(10240,'\n');

	//read facet by facet
	while (file >> tmp)
	{
		//check tmp for "facet" or "endsolid" or complain
		if (tmp[0] != 'f'
		 || tmp[1] != 'a'
		 || tmp[2] != 'c'
		 || tmp[3] != 'e'
		 || tmp[4] != 't')
		{
		//no new face starting, end of file then?
			if (tmp[0] != 'e'
			 || tmp[1] != 'n'
			 || tmp[2] != 'd'
			 || tmp[3] != 's'
			 || tmp[4] != 'o') { file.close(); return(3); }
			else break;
		}

		//read normal
		file >> tmp; //"normal" keyword
		float x,y,z;
		file >> x >> y >> z;
		Vector3F normal(x,y,z);

		//read triangle vertices
		file >> tmp;
		//check tmp for "outer" or complain
		if (tmp[0] != 'o'
		 || tmp[1] != 'u'
		 || tmp[2] != 't'
		 || tmp[3] != 'e'
		 || tmp[4] != 'r') { file.close(); return(4); }
		file >> tmp; //"loop" keyword

		file >> tmp; //"vertex" keyword
		file >> x >> y >> z;
		Vector3FC v1(x,y,z);

		file >> tmp;
		file >> x >> y >> z;
		Vector3FC v2(x,y,z);

		file >> tmp;
		file >> x >> y >> z;
		Vector3FC v3(x,y,z);

		file >> tmp; //"endloop" keyword
		file >> tmp; //"endfacet" keyword

		//add this triangle to the ActiveMesh data structures
		//we need to:
		// scale, round and use this for comparison against already
		// discovered vertices to avoid for having the same vertex saved twice
		long unsigned int o1,o2,o3;
		o1=Enlist(v1,Pos,mPos);
		o2=Enlist(v2,Pos,mPos);
		o3=Enlist(v3,Pos,mPos);
		// 
		// three offsets to the Pos array should be output
		// add them to the ID array
		ID.push_back(o1-1);
		ID.push_back(o2-1);
		ID.push_back(o3-1);
		// add normal to the norm array
		norm.push_back(normal);

		/*
		std::cout << "v1: " << v1.x << "," << v1.y << "," << v1.z << " -- o1=" << o1 << "\n";
		std::cout << "v2: " << v2.x << "," << v2.y << "," << v2.z << " -- o2=" << o2 << "\n";
		std::cout << "v3: " << v3.x << "," << v3.y << "," << v3.z << " -- o3=" << o3 << "\n";
		std::cout << "normal: " << normal.x << "," << normal.y << "," << normal.z << "\n\n";
		*/
	}
	
	file.close();
	return(0);
}

int ActiveMesh::ExportSTL(const char *filename)
{
	//try to open the file
	std::ofstream file(filename);
	if (!file.is_open()) return(1);

	file << "solid Vladimir Ulman - meshSurface testing app\n";

	for (unsigned int i=0; i < ID.size(); i+=3)
	{
		file << "facet normal " << norm[i/3].x << " " << norm[i/3].y << " " << norm[i/3].z << "\n";
		file << "outer loop\n";
		file << "vertex " << Pos[ID[i+0]].x << " " << Pos[ID[i+0]].y << " " << Pos[ID[i+0]].z << "\n";
		file << "vertex " << Pos[ID[i+1]].x << " " << Pos[ID[i+1]].y << " " << Pos[ID[i+1]].z << "\n";
		file << "vertex " << Pos[ID[i+2]].x << " " << Pos[ID[i+2]].y << " " << Pos[ID[i+2]].z << "\n";
		file << "endloop\nendfacet\n";
	}

	file.close();
	return(0);
}


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');

	//search until DATASET lines is found
	int counter=0;
	file >> tmp;
	file.ignore(10240,'\n');
	while (tmp[0] != 'D' || tmp[1] != 'A' || tmp[2] != 'T'
	    || tmp[3] != 'A' || tmp[4] != 'S' || tmp[5] != 'E')
	{
		file >> tmp;
		file.ignore(10240,'\n');
		++counter;
		if (counter == 10) { file.close(); return(35); }
	}

	//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')
	         &&
	    (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 << " 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,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;

	//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');

	//search until DATASET lines is found
	int counter=0;
	file >> tmp;
	file.ignore(10240,'\n');
	while (tmp[0] != 'D' || tmp[1] != 'A' || tmp[2] != 'T'
	    || tmp[3] != 'A' || tmp[4] != 'S' || tmp[5] != 'E')
	{
		file >> tmp;
		file.ignore(10240,'\n');
		++counter;
		if (counter == 10) { file.close(); return(35); }
	}

	//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);

		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";

	file.close();
	return(0);
}


int ActiveMesh::ImportVTK_Ftree(const char *filename,const int idx,
                                const float stretch,bool resetMesh)
{
	const int PointsOnRadiusPeriphery=10;
	const float radiusCorrection=stretch;

	//make sure the F-tree array is long enough to host the requested index
	while (Ftree.size() < idx+1) Ftree.push_back(mesh_t());
	mesh_t& filoMesh=Ftree[idx];

	if (resetMesh)
	{
		filoMesh.Pos.clear();
		filoMesh.ID.clear();
		filoMesh.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');

	//search until DATASET lines is found
	int counter=0;
	file >> tmp;
	file.ignore(10240,'\n');
	while (tmp[0] != 'D' || tmp[1] != 'A' || tmp[2] != 'T'
	    || tmp[3] != 'A' || tmp[4] != 'S' || tmp[5] != 'E')
	{
		file >> tmp;
		file.ignore(10240,'\n');
		++counter;
		if (counter == 10) { file.close(); return(35); }
	}

	//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=filoMesh.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]];
			filoMesh.Pos.push_back(V);
		}
	}

	//finally, add the tip point
	filoMesh.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)
		{
			filoMesh.ID.push_back(firstRadiusPoint+p);
			filoMesh.ID.push_back(firstRadiusPoint+p+1);
			filoMesh.ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
			filoMesh.norm.push_back(fictiveNormal);

			filoMesh.ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery+1);
			filoMesh.ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
			filoMesh.ID.push_back(firstRadiusPoint+p+1);
			filoMesh.norm.push_back(fictiveNormal);
		}

		filoMesh.ID.push_back(firstRadiusPoint+p);
		filoMesh.ID.push_back(firstRadiusPoint);
		filoMesh.ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
		filoMesh.norm.push_back(fictiveNormal);

		filoMesh.ID.push_back(firstRadiusPoint  +PointsOnRadiusPeriphery);
		filoMesh.ID.push_back(firstRadiusPoint+p+PointsOnRadiusPeriphery);
		filoMesh.ID.push_back(firstRadiusPoint);
		filoMesh.norm.push_back(fictiveNormal);

		firstRadiusPoint+=PointsOnRadiusPeriphery;
	}

	//the last segment...
	int p=0;
	for (; p < PointsOnRadiusPeriphery-1; ++p)
	{
		filoMesh.ID.push_back(firstRadiusPoint+p);
		filoMesh.ID.push_back(firstRadiusPoint+p+1);
		filoMesh.ID.push_back(firstRadiusPoint+PointsOnRadiusPeriphery);
		filoMesh.norm.push_back(fictiveNormal);
	}

	filoMesh.ID.push_back(firstRadiusPoint+p);
	filoMesh.ID.push_back(firstRadiusPoint);
	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
			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=200.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);
}


void ActiveMesh::RenderMask_body(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles)
{
	//time savers: resolution
	const float xRes=mask.GetResolution().GetX();
	const float yRes=mask.GetResolution().GetY();
	const float zRes=mask.GetResolution().GetZ();

	//time savers: offset
	const float xOff=mask.GetOffset().x;
	const float yOff=mask.GetOffset().y;
	const float zOff=mask.GetOffset().z;

	//create "aside" empty image similar to the input-output one
	i3d::Image3d<i3d::GRAY16> tmpI;
	tmpI.CopyMetaData(mask);
	tmpI.GetVoxelData()=0;

	//over all triangles
	for (unsigned int i=0; i < ID.size()/3; ++i)
	//for (unsigned int i=0; i < 15; ++i)
	{
		const Vector3F& v1=Pos[ID[3*i+0]];
		const Vector3F& v2=Pos[ID[3*i+1]];
		const Vector3F& v3=Pos[ID[3*i+2]];

		//sweeping (and rendering) the triangle
		for (float c=0.f; c <= 1.0f; c += 0.1f)
		for (float b=0.f; b <= (1.0f-c); b += 0.1f)
		{
			float a=1.0f -b -c;

			Vector3F v=a*v1;
			v+=b*v2;
			v+=c*v3;
			
			//std::cout << "ID #" << i << ": v=(" << v.x << "," << v.y << "," << v.z << ")\n";

			//nearest neighbor pixel coordinate
			const int x=(int)roundf( (v.x-xOff) *xRes);
			const int y=(int)roundf( (v.y-yOff) *yRes);
			const int z=(int)roundf( (v.z-zOff) *zRes);

			short val=(showTriangles)? short(i%5 *30 +100) : 100;
			if (tmpI.Include(x,y,z)) tmpI.SetVoxel(x,y,z,val);
		}
	}

	//this floods cell exterior from the image corner
	i3d::Image3d<i3d::GRAY16> tmpMask(tmpI);
	i3d::Dilation(tmpMask,tmpI,i3d::nb3D_o18);
	i3d::FloodFill(tmpI,(i3d::GRAY16)50,0);

	//this removes the flooding while filling
	//everything else (and closing holes in this was)
	i3d::GRAY16* pI=tmpI.GetFirstVoxelAddr();
	i3d::GRAY16* pM=mask.GetFirstVoxelAddr();
	i3d::GRAY16* const pL=pM+mask.GetImageSize();
	while (pM != pL)
	{
		*pM=(*pI != 50)? std::max(*pI,i3d::GRAY16(100)) : *pM;
		++pI; ++pM;
	}
}

void ActiveMesh::RenderMask_filo(i3d::Image3d<i3d::GRAY16>& mask,const int idx,const bool showTriangles)
{
	//if index out of range, do nothing
	if (idx >= Ftree.size()) return;

	//time savers: resolution
	const float xRes=mask.GetResolution().GetX();
	const float yRes=mask.GetResolution().GetY();
	const float zRes=mask.GetResolution().GetZ();

	//time savers: offset
	const float xOff=mask.GetOffset().x;
	const float yOff=mask.GetOffset().y;
	const float zOff=mask.GetOffset().z;

	//create "aside" empty image similar to the input-output one
	i3d::Image3d<i3d::GRAY16> tmpI;
	tmpI.CopyMetaData(mask);
	tmpI.GetVoxelData()=0;

	//over all triangles of the given filopodia
	const mesh_t& filoMesh=Ftree[idx];
	for (unsigned int i=0; i < filoMesh.ID.size()/3; ++i)
	{
		const Vector3F& v1=filoMesh.Pos[filoMesh.ID[3*i+0]];
		const Vector3F& v2=filoMesh.Pos[filoMesh.ID[3*i+1]];
		const Vector3F& v3=filoMesh.Pos[filoMesh.ID[3*i+2]];

		//sweeping (and rendering) the triangle
		for (float c=0.f; c <= 1.0f; c += 0.1f)
		for (float b=0.f; b <= (1.0f-c); b += 0.1f)
		{
			float a=1.0f -b -c;

			Vector3F v=a*v1;
			v+=b*v2;
			v+=c*v3;
			
			//nearest neighbor pixel coordinate
			const int x=(int)roundf( (v.x-xOff) *xRes);
			const int y=(int)roundf( (v.y-yOff) *yRes);
			const int z=(int)roundf( (v.z-zOff) *zRes);

			short val=(showTriangles)? short(i%5 *30 +100) : 101+idx;
			if (tmpI.Include(x,y,z)) tmpI.SetVoxel(x,y,z,val);
		}
	}

	//this floods cell exterior from the image corner
	i3d::Image3d<i3d::GRAY16> tmpMask(tmpI);
	i3d::Dilation(tmpMask,tmpI,i3d::nb3D_o18);
	i3d::FloodFill(tmpI,(i3d::GRAY16)50,0);

	//this removes the flooding while filling
	//everything else (and closing holes in this was)
	i3d::GRAY16* pI=tmpI.GetFirstVoxelAddr();
	i3d::GRAY16* pM=mask.GetFirstVoxelAddr();
	i3d::GRAY16* const pL=pM+mask.GetImageSize();
	while (pM != pL)
	{
		*pM=(*pI != 50)? std::max(*pI,i3d::GRAY16(101+idx)) : *pM;
		++pI; ++pM;
	}
}


void ActiveMesh::RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask)
{
	//time savers: resolution
	const float xRes=mask.GetResolution().GetX();
	const float yRes=mask.GetResolution().GetY();
	const float zRes=mask.GetResolution().GetZ();

	//time savers: offset
	const float xOff=mask.GetOffset().x;
	const float yOff=mask.GetOffset().y;
	const float zOff=mask.GetOffset().z;

	//over all triangles
	//for (unsigned int i=0; i < ID.size()/3; ++i)
	for (unsigned int i=0; i < 15; ++i)
	{
		const Vector3F& v1=Pos[ID[3*i+0]];
		const Vector3F& v2=Pos[ID[3*i+1]];
		const Vector3F& v3=Pos[ID[3*i+2]];

		//Vector3F e12=v2-v1; //ID=0
		//Vector3F e13=v3-v1; //ID=1
		//Vector3F e23=v3-v2; //ID=2
		Vector3F edges[3]={v2-v1,v3-v1,v3-v2};

		short order[3]={0,1,2};
		if (edges[1].LenQ() < edges[0].LenQ()) { order[0]=1; order[1]=0; }
		if (edges[2].LenQ() < edges[order[1]].LenQ())
		{
			order[2]=order[1];
			order[1]=2;

			if (edges[2].LenQ() < edges[order[0]].LenQ())
			{
				order[1]=order[0];
				order[0]=2;
			}
		}

		//point within a triangle will of the form: vv + b*vb + c*vc
		//vb is the shortest edge, vc is the longest edge
		//vv is their common vertex
		//vb,vc are oriented outward from this vertex
		Vector3F vv,vb,vc;
		if (order[0]==0 && order[2]==1)
		{
			//order[0] "points" at shortest edge
			vv=v1;
			vb=edges[0];
			vc=edges[1];
		}
		else
		if (order[0]==0 && order[2]==2)
		{
			vv=v2;
			vb=-edges[0];
			vc=edges[2];
		}
		else
		if (order[0]==1 && order[2]==0)
		{
			vv=v1;
			vb=edges[1];
			vc=edges[0];
		}
		else
		if (order[0]==1 && order[2]==2)
		{
			vv=v3;
			vb=-edges[1];
			vc=-edges[2];
		}
		else
		if (order[0]==2 && order[2]==0)
		{
			vv=v2;
			vb=edges[2];