debug du clear des cartes

assert dans les CellMarker et DartMarker si invalides
qq modif dans MarchingCube

Apps/Tuto/tuto_mt.cpp

@@ -70,14 +70,21 @@ AttributeHandler<PFP::VEC3> normal;
 void MyQT::cb_Open()
 {
 	// set some filters
-	std::string filters("all (*.*);; trian (*.trian);; ctm (*.ctm);; off (*.off);; ply (*.ply)");
-
-	std::string filename = selectFile("OpenMesh","",filters);
+//	std::string filters("all (*.*);; trian (*.trian);; ctm (*.ctm);; off (*.off);; ply (*.ply)");
+//
+//	std::string filename = selectFile("OpenMesh","",filters);
+//
+//	std::vector<std::string> attrNames ;
+//	if(!Algo::Import::importMesh<PFP>(myMap, filename.c_str(), attrNames))
+//	{
+//		CGoGNerr << "could not import " << filename << CGoGNendl ;
+//		return;
+//	}
 
 	std::vector<std::string> attrNames ;
-	if(!Algo::Import::importMesh<PFP>(myMap, filename.c_str(), attrNames))
+	if(!Algo::Import::importMesh<PFP>(myMap, "/home/thery/Data/liver.trian", attrNames))
 	{
-		CGoGNerr << "could not import " << filename << CGoGNendl ;
+		CGoGNerr << "could not import xxx" << CGoGNendl ;
 		return;
 	}
 
@@ -133,34 +140,45 @@ void MyQT::cb_initGL()
 	glewInit();
 
 	// init lighting parameters
-	float lightPosition[4]= {0.0f,0.0f,10000.0f,1.0f};
+	float lightPosition[4]= {10.0f,10.0f,10000.0f,1.0f};
 	float lightColor[4]= {0.9f,0.9f,0.9f,1.0f};
 
 	glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
-	glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, 1);
+	glEnable(GL_LIGHT0);
 	glLightfv(GL_LIGHT0, GL_DIFFUSE, lightColor);
 	glLightfv(GL_LIGHT0, GL_POSITION, lightPosition);
-	glEnable(GL_LIGHT0);
-	glEnable(GL_COLOR_MATERIAL);
-	glColorMaterial(GL_FRONT_AND_BACK, GL_DIFFUSE);
+	glEnable(GL_NORMALIZE);
+//	glDisable(GL_CULL_FACE);
+//	glFrontFace(GL_CCW);
+
 }
 
 void MyQT::cb_redraw()
 {
+    GLfloat diff[4]= {0.0f,1.0f,0.1f,1.0f};
+    GLfloat amb[4]= {0.1f,0.0f,0.1f,1.0f};
+    GLfloat spec[4]= {1.0f,1.0f,1.0f,1.0f};
+    float shininess=125.0f;
+
 	// draw the lines
-	glDisable(GL_LIGHTING);
-	glColor3f(0.0f, 0.0f, 0.3f);
-	glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
-	glDisable(GL_LIGHTING);
-	glEnable(GL_SMOOTH);
-	Algo::Render::GL1::renderTriQuadPoly<PFP>(myMap,Algo::Render::GL1::LINE, 1.0f,position, normal);
+//	glDisable(GL_LIGHTING);
+//	glColor3f(0.0f, 0.0f, 0.3f);
+//	glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
+//	glDisable(GL_LIGHTING);
+//
+//	Algo::Render::GL1::renderTriQuadPoly<PFP>(myMap,Algo::Render::GL1::LINE, 1.0f,position, normal);
 
 	// draw the faces
 	glEnable(GL_POLYGON_OFFSET_FILL);
 	glPolygonOffset(1.0f, 1.0f);
-	glColor3f(0.1f, 0.8f, 0.0f);
 	glEnable(GL_LIGHTING);
+	glEnable(GL_SMOOTH);
 	glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+
+    glMaterialfv(GL_FRONT, GL_DIFFUSE, diff);
+    glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
+    glMaterialfv(GL_FRONT,GL_SPECULAR,spec);
+    glMaterialf( GL_FRONT, GL_SHININESS, shininess);
 	Algo::Render::GL1::renderTriQuadPoly<PFP>(myMap,Algo::Render::GL1::SMOOTH, 1.0f,position, normal);
 	glDisable(GL_POLYGON_OFFSET_FILL);
 }
@@ -370,25 +388,27 @@ int main(int argc, char **argv)
  	sqt.statusMsg("Neww to create a sphere or Load for a mesh file");
  	CGoGNStream::allToConsole(&sqt);
 
- 	int xx = 3;
- 	double yy = 2.5;
- 	bool zz=true;
- 	int kk=32;
- 	int cc=2;
-
- 	{
- 	using namespace CGoGN::Utils::QT;
-
- 	inputValues(	VarInteger(0,20,xx, "Entier",
- 					VarBool(zz, "Bool",
- 					VarDbl(0.314,3.14,yy,"Double",
- 					VarSlider(10,100,kk,"Slider",
- 					VarCombo("Riri;Fifi;Loulou;Donald",cc,"Combo") )))));
- 	}
-
- 	std::cout << "Int:" << xx << "  Bool:"<<zz<< std::endl;
- 	std::cout << "Dbl:" << yy << "  Slider:"<<kk<< std::endl;
- 	std::cout << "Combo:" << cc << std::endl;
+ 	sqt.cb_Open();
+
+// 	int xx = 3;
+// 	double yy = 2.5;
+// 	bool zz=true;
+// 	int kk=32;
+// 	int cc=2;
+//
+// 	{
+// 	using namespace CGoGN::Utils::QT;
+//
+// 	inputValues(	VarInteger(0,20,xx, "Entier",
+// 					VarBool(zz, "Bool",
+// 					VarDbl(0.314,3.14,yy,"Double",
+// 					VarSlider(10,100,kk,"Slider",
+// 					VarCombo("Riri;Fifi;Loulou;Donald",cc,"Combo") )))));
+// 	}
+//
+// 	std::cout << "Int:" << xx << "  Bool:"<<zz<< std::endl;
+// 	std::cout << "Dbl:" << yy << "  Slider:"<<kk<< std::endl;
+// 	std::cout << "Combo:" << cc << std::endl;
 
 	return app.exec();
 }

README.TXT

 Dépendences Linux:
 installer les paquets suivants:
-cmake libXi-dev libXmu-dev freeglut3-dev libdevil-dev libglew-dev libgmp3-dev libxml2-dev libboost-dev libboost-thread-dev libzip-dev libqt4-help qt4-designer qt4-dev-tools uuid-dev
+cmake libXi-dev libXmu-dev freeglut3-dev libdevil-dev libglew-dev libgmp3-dev libxml2-dev libboost-all-dev libzip-dev libqt4-help qt4-designer qt4-dev-tools uuid-dev
 
 Pour compiler CGoGN:
 - aller dans ThirdParty et taper "cmake .", puis make

include/Algo/MC/image.h

@@ -332,6 +332,17 @@ public:
 
 	float computeCurvatureCount(const DataType *ptrVox, const std::vector<int>& sphere, DataType val);
 
+	float computeCurvatureCount3(const DataType *ptrVox, const std::vector<int>& cylX, const std::vector<int>& cylY, const std::vector<int>& cyl2, DataType val);
+
+	void createMaskOffsetCylinders(std::vector<int>& tableX, std::vector<int>& tableY, std::vector<int>& tableZ, int _i32radius);
+
+	Image<DataType>* cropz(unsigned int zmin, unsigned int zmax);
+
+	void createNormalSphere(std::vector<Geom::Vec3f>& table, int _i32radius);
+
+	Geom::Vec3f computeNormal(DataType *ptrVox, const std::vector<Geom::Vec3f>& sphere, DataType val, unsigned int radius);
+
+	bool checkSaddlecomputeNormal(const Geom::Vec3f& P, const Geom::Vec3f& normal, unsigned int radius);
 };
 
 

include/Algo/MC/image.hpp

@@ -479,7 +479,6 @@ Image<DataType>* Image<DataType>::Blur3()
 
 
 	return newImg;
-
 }
 
 template<typename DataType>
@@ -548,6 +547,325 @@ float Image<DataType>::computeCurvatureCount(const DataType *ptrVox, const std::
 }
 
 
+// TESTING NEW CURVATURE METHODS
+template<typename DataType>
+void Image<DataType>::createMaskOffsetCylinders(std::vector<int>& tableX, std::vector<int>& tableY, std::vector<int>& tableZ, int _i32radius)
+{
+	// compute the width of the sphere for memory allocation
+	int i32Width = 2*_i32radius + 1;
+	// squared radius
+    float fRad2 = (float)(_i32radius*_i32radius);
+
+	// memory allocation
+	// difficult to know how many voxels before computing,
+	// so the reserve for the BB
+	tableX.reserve(i32Width*i32Width*7);
+	tableX.clear();
+	tableY.reserve(i32Width*i32Width*7);
+	tableY.clear();
+	tableZ.reserve(i32Width*i32Width*7);
+	tableZ.clear();
+
+
+	// scan all the BB of the sphere
+	for (int z = -_i32radius;  z<=_i32radius; z++)
+	{
+		for (int y = -_i32radius;  y<=_i32radius; y++)
+		{
+			for (int x = -_i32radius;  x<=_i32radius; x++)
+			{
+				Geom::Vec3f v((float)x,(float)y,(float)z);
+				float fLength =  v.norm2();
+				// if inside the sphere
+				if (fLength<=fRad2)
+				{
+					// the the index of the voxel
+					int index = z * m_WXY + y * m_WX + x;
+
+					if ((x<=3) && (x>=-3))
+						tableX.push_back(index);
+					if ((y<=3) && (y>=-3))
+						tableY.push_back(index);
+					if ((z<=3) && (z>=-3))
+						tableZ.push_back(index);
+
+				}
+			}
+		}
+	}
+}
+
+template<typename DataType>
+float Image<DataType>::computeCurvatureCount3(const DataType *ptrVox, const std::vector<int>& cylX, const std::vector<int>& cylY, const std::vector<int>& cylZ, DataType val)
+{
+	int noir = 0;
+	int blanc = 0;
+
+	float vals[3];
+
+	for (std::vector<int>::const_iterator it=cylX.begin(); it!=cylX.end();it++)
+	{
+		const DataType *data = ptrVox + *it;
+		if (*data != val)
+		{
+			blanc++;
+		}
+		else
+		{
+			noir++;
+		}
+	}
+
+
+	if (blanc >= noir)
+	{
+		vals[0] = 1.0f - ((float)noir) / ((float)blanc);
+	}
+	else
+	{
+		vals[0] =  -1.0f + ((float)blanc)/((float)noir);
+	}
+
+	noir = 0;
+	blanc = 0;
+
+	for (std::vector<int>::const_iterator it=cylY.begin(); it!=cylY.end();it++)
+	{
+		const DataType *data = ptrVox + *it;
+		if (*data != val)
+		{
+			blanc++;
+		}
+		else
+		{
+			noir++;
+		}
+	}
+
+	if (blanc >= noir)
+	{
+		vals[1] = 1.0f - ((float)noir) / ((float)blanc);
+	}
+	else
+	{
+		vals[1] =  -1.0f + ((float)blanc)/((float)noir);
+	}
+
+
+	noir = 0;
+	blanc = 0;
+
+	for (std::vector<int>::const_iterator it=cylZ.begin(); it!=cylZ.end();it++)
+	{
+		const DataType *data = ptrVox + *it;
+		if (*data != val)
+		{
+			blanc++;
+		}
+		else
+		{
+			noir++;
+		}
+	}
+
+	if (blanc >= noir)
+	{
+		vals[2] = 1.0f - ((float)noir) / ((float)blanc);
+	}
+	else
+	{
+		vals[2] =  -1.0f + ((float)blanc)/((float)noir);
+	}
+
+//	if ((valZ>valX) && (valZ>valY))
+//		return (valX + valY)/2.0f;
+//
+//	if ((valY>valX) && (valY>valZ))
+//		return (valX + valZ)/2.0f;
+//
+//	return (valY + valZ)/2.0f;
+
+	unsigned short m1,m2;
+	if ((abs(vals[0]) < abs(vals[1])) && (abs(vals[0]) < abs(vals[2])))
+	{
+		m1 =1;
+		m2 =2;
+	}
+	else
+	{
+		m1=0;
+		if (abs(vals[1]) < abs(vals[2]))
+			m2 = 2;
+		else
+			m2 = 1;
+	}
+
+	if (vals[m1]>0.0f)
+		if (vals[m2]<0.0f)
+			if ((vals[m1] - vals[m2])>0.8f)
+				return 1.0f;
+
+	if (vals[m1]<0.0f)
+		if (vals[m2]>0.0f)
+			if ((vals[m2] - vals[m1])>0.8f)
+			return 1.0f;
+
+	return std::max(std::max(abs(vals[0]),abs(vals[1])),abs(vals[2]));
+
+}
+
+
+
+template< typename  DataType >
+Image<DataType>* Image<DataType>::cropz(unsigned int zmin, unsigned int nb)
+{
+
+	unsigned int zmax = zmin+nb;
+	if (zmax> m_WZ)
+	{
+		zmax = m_WZ;
+		nb = zmax - zmin;
+	}
+
+	DataType* data2 = new DataType[150*m_WY*m_WZ];
+	Image<DataType>* newImg = new Image<DataType>(data2,150,m_WY,nb,getVoxSizeX(),getVoxSizeY(),getVoxSizeZ());
+	newImg->m_Alloc=true;
+	// set origin of real data in image ??
+
+	for(unsigned int z=zmin; z< zmax; ++z)
+	{
+		for(int y=0; y<m_WY; ++y)
+		{
+			for(int x=0; x<150; ++x)
+			{
+				DataType* ori = getVoxelPtr(x,y,z);
+				DataType* dest = newImg->getVoxelPtr(x,y,z-zmin);
+				*dest = *ori;
+			}
+		}
+	}
+
+	return newImg;
+}
+
+
+
+template<typename DataType>
+void Image<DataType>::createNormalSphere(std::vector<Geom::Vec3f>& table, int _i32radius)
+{
+	// compute the width of the sphere for memory allocation
+	int i32Width = 2*_i32radius + 1;
+
+	table.reserve(i32Width*i32Width*i32Width);
+	table.clear();
+
+	// scan all the BB of the sphere
+	for (int z = -_i32radius;  z<=_i32radius; z++)
+	{
+		for (int y = -_i32radius;  y<=_i32radius; y++)
+		{
+			for (int x = -_i32radius;  x<=_i32radius; x++)
+			{
+				Geom::Vec3f v((float)x,(float)y,(float)z);
+				float fLength =  v.normalize();
+				// if inside the sphere
+				if (fLength<=_i32radius)
+					table.push_back(v);
+				else
+					table.push_back(Geom::Vec3f(0.0f,0.0f,0.0f));
+			}
+		}
+	}
+}
+
+
+template<typename DataType>
+Geom::Vec3f Image<DataType>::computeNormal(DataType *ptrVox, const std::vector<Geom::Vec3f>& sphere, DataType val, unsigned int radius)
+{
+
+	DataType *ucDat1 = ptrVox - radius*(m_WX + m_WX*m_WY);
+	Geom::Vec3f normal(0.0f,0.0f,0.0f);
+	unsigned int width = 2*radius+1;
+
+	unsigned int ind = 0;
+	for (unsigned int i=0; i<width;++i)
+	{
+		DataType *ucDat2 = ucDat1;
+		for (unsigned int j=0; j<width;++j)
+		{
+			DataType *ucDat3 = ucDat2;
+			for (unsigned int k=0; k<width;++k)
+			{
+				if (*ucDat3 == val)
+				{
+					normal += sphere[ind];
+				}
+				++ucDat3;
+				++ind;
+			}
+
+			ucDat2 += m_WX;
+		}
+		ucDat1 += m_WXY;
+	}
+
+	normal.normalize();
+	return normal;
+}
+
+template<typename DataType>
+bool Image<DataType>::checkSaddlecomputeNormal(const Geom::Vec3f& P, const Geom::Vec3f& normal, unsigned int radius)
+{
+	Geom::Vec3f V1;
+	if ( (fabs(normal[0]) <fabs(normal[1])) && (fabs(normal[0]) <fabs(normal[2])) )
+	{
+		V1 = normal ^ Geom::Vec3f(1.0f,0.0f,0.0);
+	}
+	else if (fabs(normal[1]) <fabs(normal[2]))
+	{
+		V1 = normal ^ Geom::Vec3f(0.0f,1.0f,0.0);
+	}
+	else
+	{
+		V1 = normal ^ Geom::Vec3f(0.0f,0.0f,1.0);
+	}
+
+	Geom::Vec3f V2 = normal ^ V1;
+
+	Geom::Vec3f Q = P + (0.5f * radius)*normal;
+
+	float le = 0.866f * radius; // (cos30)
+	Geom::Vec3f Q1 = Q + le*V1;
+	Geom::Vec3f Q2 = Q - le*V1;
+	DataType X1 = getVoxel(int(floor(Q1[0])), int(floor(Q1[1])), int(floor(Q1[2])));
+	DataType X2 = getVoxel(int(floor(Q2[0])), int(floor(Q2[1])), int(floor(Q2[2])));
+	Q1 = Q + le*V2;
+	Q2 = Q - le*V2;
+	DataType X3 = getVoxel(int(floor(Q1[0])), int(floor(Q1[1])), int(floor(Q1[2])));
+	DataType X4 = getVoxel(int(floor(Q2[0])), int(floor(Q2[1])), int(floor(Q2[2])));
+
+	le *= 0.707f; // (sqrt(2)/2)
+	Q1 = Q + le*(V1+V2);
+	Q2 = Q - le*(V1+V2);
+	DataType X5 = getVoxel(int(floor(Q1[0])), int(floor(Q1[1])), int(floor(Q1[2])));
+	DataType X6 = getVoxel(int(floor(Q2[0])), int(floor(Q2[1])), int(floor(Q2[2])));
+
+	Q1 = Q + le*(V1-V2);
+	Q2 = Q - le*(V1-V2);
+	DataType X7 = getVoxel(int(floor(Q1[0])), int(floor(Q1[1])), int(floor(Q1[2])));
+	DataType X8 = getVoxel(int(floor(Q2[0])), int(floor(Q2[1])), int(floor(Q2[2])));
+
+	if  ((X1 == X2) && (X3==X4) && (X1!=X3))
+		return true;
+
+	if  ((X5 == X6) && (X7==X8) && (X5!=X7))
+		return true;
+
+	return false;
+}
+
+
+
 } // end namespace
 } // end namespace
 } // end namespace

include/Algo/MC/marchingcube.h

@@ -260,6 +260,8 @@ public:
 	 */
 	Geom::Vec3f boundMax() const {return m_Image->boundMax();}
 
+	void removeFacesOfBoundary(AttributeHandler<unsigned char>& boundVertices, unsigned int frameWidth);
+
 };
 
 

include/Algo/MC/marchingcube.hpp

@@ -22,7 +22,8 @@
 *                                                                              *
 *******************************************************************************/
 
-#include "windowing.h"
+#include "Algo/MC/windowing.h"
+#include "Topology/generic/dartmarker.h"
 
 namespace CGoGN
 {
@@ -1192,6 +1193,64 @@ void MarchingCube<DataType, Windowing, PFP>::createLocalFaces(const unsigned cha
 
 
 
+template< typename  DataType, template < typename D2 > class Windowing, typename PFP >
+void MarchingCube<DataType, Windowing, PFP>::removeFacesOfBoundary(AttributeHandler<unsigned char>& boundVertices, unsigned int frameWidth)
+{
+
+	float xmin = frameWidth;
+	float xmax = m_Image->getWidthX() - frameWidth;
+	float ymin = frameWidth;
+	float ymax = m_Image->getWidthY() - frameWidth;
+	float zmin = frameWidth;
+	float zmax = m_Image->getWidthZ() - frameWidth;
+
+
+
+	// traverse position and create bound attrib
+	for(unsigned int it = m_positions.begin(); it != m_positions.end(); m_positions.next(it))
+	{
+		bool bound = (m_positions[it][0]<=xmin) || (m_positions[it][0]>=xmax) || \
+					 (m_positions[it][1]<=ymin) || (m_positions[it][1]>=ymax) || \
+					 (m_positions[it][2]<=zmin) || (m_positions[it][2]>=zmax);
+
+		if (bound)
+		{
+			boundVertices[it] = 1;
+		}
+		else
+			boundVertices[it] = 0;
+	}
+
+	// traverse face and check if all vertices are bound
+	DartMarker mf(*m_map);
+	for (Dart d = m_map->begin(); d != m_map->end();)
+	{
+		if (!mf.isMarked(d))
+		{
+			Dart dd = d;
+			Dart e = m_map->phi1(d);
+			Dart f = m_map->phi1(e);
+			m_map->next(d);
+			while ((d==e) || (d==f))
+			{
+				m_map->next(d);
+			}
+			if ((boundVertices[dd]!=0) && (boundVertices[e]!=0) && (boundVertices[f]!=0))
+				m_map->deleteFace(dd);
+			else
+				mf.markOrbit(FACE,dd);
+		}
+	}
+
+
+
+}
+
+
+
+
+
+
 } // end namespace
 } // end namespace
 } // end namespace

include/Topology/generic/cellmarker.h

@@ -56,12 +56,12 @@ public:
 	{
 		if(!map.isOrbitEmbedded(cell))
 			map.addEmbedding(cell) ;
-		m_mark = m_map.m_marksets[m_cell][m_thread].getNewMark() ;
+		m_mark = m_map.getMarkerSet(m_cell,m_thread).getNewMark() ;
 	}
 
 	virtual ~CellMarkerGen()
 	{
-		m_map.m_marksets[m_cell][m_thread].releaseMark(m_mark) ;
+		m_map.getMarkerSet(m_cell,m_thread).releaseMark(m_mark) ;
 	}
 
 protected:
@@ -75,6 +75,9 @@ public:
 	 */
 	virtual void mark(Dart d)
 	{
+		assert(m_map.getMarkerSet(m_cell,m_thread).testMark(m_mark));
+		assert(m_map.getMarkVector(m_cell, m_thread) != NULL);
+
 		unsigned int a = m_map.getEmbedding(m_cell, d) ;
 		if (a == EMBNULL)
 			a = m_map.embedNewCell(m_cell, d) ;
@@ -87,6 +90,9 @@ public:
 	 */
 	virtual void unmark(Dart d)