Added initial (just contour yet) rendering of meshes into mask images.

CMakeLists.txt

@@ -55,6 +55,12 @@ else (WIN32)
 		  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
 endif (WIN32)
 
+#
+# due to Lapack+i3dalgo
+#
+
+SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Xlinker -defsym -Xlinker MAIN__=main")
+
 #
 # user defined options
 #
@@ -102,6 +108,13 @@ include_directories(${INC_GSL})
 find_library(LIB_LAPACK lapack)
 set(LIBS ${LIBS} ${LIB_LAPACK})
 
+# i3dlibs libraries
+find_path(INC_I3D "i3d/image3d.h")
+include_directories(${INC_I3D})
+find_library(LIB_I3DCORE i3dcore)
+find_library(LIB_I3DALGO i3dalgo)
+set(LIBS ${LIBS} ${LIB_I3DCORE} ${LIB_I3DALGO})
+
 #--------------
 # SOURCE FILES
 #--------------

cmath3d/TriangleMesh.cpp

@@ -190,6 +190,170 @@ int ActiveMesh::ImportVTK(const char *filename)
 }
 
 
+void ActiveMesh::RenderMask(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]];
+
+		//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);
+
+			if (mask.Include(x,y,z)) mask.SetVoxel(x,y,z,i);
+		}
+	}
+}
+
+
+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];
+			vc=-edges[0];
+		}
+		else
+		if (order[0]==2 && order[2]==1)
+		{
+			vv=v3;
+			vb=-edges[2];
+			vc=-edges[1];
+		}
+
+		//optimal increment for the short and long edge, respectively
+		float db=(vb.x*xRes > vb.y*yRes)? vb.x*xRes : vb.y*yRes;
+		db=(vb.z*zRes > db)? vb.z*zRes : db;
+		db=1.f/db;
+
+		float dc=(vc.x*xRes > vc.y*yRes)? vc.x*xRes : vc.y*yRes;
+		dc=(vc.z*zRes > dc)? vc.z*zRes : dc;
+		dc=1.f/dc;
+
+		std::cout << "\nID #" << i << ": v1=(" << v1.x << "," << v1.y << "," << v1.z << ")\n";
+		std::cout << "ID #" << i << ": v2=(" << v2.x << "," << v2.y << "," << v2.z << ")\n";
+		std::cout << "ID #" << i << ": v3=(" << v3.x << "," << v3.y << "," << v3.z << ")\n";
+		
+		//sweeping (and rendering) the triangle
+		for (float c=0.f; c <= 1.0f; c += dc)
+		for (float b=0.f; b <= (1.0f-c); b += db)
+		{
+			Vector3F v=vv;
+			v+=b*vb;
+			v+=c*vc;
+			
+			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);
+
+			if (mask.Include(x,y,z)) mask.SetVoxel(x,y,z,i);
+		}
+	}
+}
+
+
 void ActiveMesh::CenterMesh(const Vector3F& newCentre)
 {
 	//calc geom. centre

cmath3d/TriangleMesh.h

@@ -2,6 +2,7 @@
 #define __TRIANGLEMESH__FOR_MT__
 
 #include <vector>
+#include <i3d/image3d.h>
 #include "../src/params.h"
 
 /**
@@ -28,6 +29,8 @@ 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);
+	void RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask);
 
 	void CenterMesh(const Vector3F& newCentre);
 	void ScaleMesh(const Vector3F& scale);

src/main.cpp

-#include <list>
-#include <string.h>
-#include <stdlib.h>		//for atoi()
+#include <i3d/image3d.h>
 
 #include "params.h"
 #include "graphics.h"
@@ -36,6 +34,15 @@ int main(void)
 	mesh.CenterMesh(params.sceneCentre);
 
 	//work with 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");
 
 	initializeGL();
 	loopGL();

src/params.h

@@ -71,6 +71,9 @@ class Vector3d {
 		return( *this );
 	}
 
+	/// unary minus = negate
+	Vector3d<T> operator-() const { return Vector3d<T>(-x,-y,-z); }
+
 	Vector3d<T>& operator=(const T scalar)
 	{
 		this->x=scalar;