From cf7681c97c735b6f6ffcce1b95b647aa1558cb36 Mon Sep 17 00:00:00 2001
From: Sylvain Thery <sylvain@thery.pipo>
Date: Fri, 14 Sep 2012 09:23:49 +0200
Subject: [PATCH] enhancement of curvature computation

---
 include/Algo/MC/image.h   | 16 ++++++++++++++--
 include/Algo/MC/image.hpp | 31 ++++++++++++++++++-------------
 2 files changed, 32 insertions(+), 15 deletions(-)

diff --git a/include/Algo/MC/image.h b/include/Algo/MC/image.h
index 9c0e8f276..d7da4f843 100644
--- a/include/Algo/MC/image.h
+++ b/include/Algo/MC/image.h
@@ -329,10 +329,22 @@ public:
 	 */
 	Image<DataType>* Blur3();
 
-
+	/**
+	 * create a virtual sphere for computing curvature
+	 * @param table the vector a ptr offset of voxels of sphere
+	 * @param _i32radius radius of sphere
+	 */
 	void createMaskOffsetSphere(std::vector<int>& table, int _i32radius);
 
-	float computeCurvatureCount(const DataType *ptrVox, const std::vector<int>& sphere, DataType val);
+	/**
+	 * compute the curvature with method that count voxels in a virtual sphere
+	 * @param x position in x
+	 * @param y position in y
+	 * @param z position in z
+	 * @param sphere the precomputed sphere
+	 * @param outsideVal value of outside object
+	 */
+	float computeCurvatureCount(float x, float y, float z, const std::vector<int>& sphere, DataType outsideVal);
 
 	float computeCurvatureCount3(const DataType *ptrVox, const std::vector<int>& cylX, const std::vector<int>& cylY, const std::vector<int>& cyl2, DataType val);
 
diff --git a/include/Algo/MC/image.hpp b/include/Algo/MC/image.hpp
index 200531f5e..eb38ce95e 100644
--- a/include/Algo/MC/image.hpp
+++ b/include/Algo/MC/image.hpp
@@ -517,29 +517,29 @@ void Image<DataType>::createMaskOffsetSphere(std::vector<int>& table, int _i32ra
 	float ys = m_SY/smin;
 	float zs = m_SZ/smin;
 
+	int radX = ceil(float(_i32radius)/xs);
+	int radY = ceil(float(_i32radius)/ys);
+	int radZ = ceil(float(_i32radius)/zs);
 
-	// compute the width of the sphere for memory allocation
-	int i32Width = 2*_i32radius + 1;
-	// squared radius
-    float fRad2 = (float)(_i32radius*_i32radius);
+	float sRadius = sqrt( double(_i32radius)/xs*double(_i32radius)/xs +  double(_i32radius)/ys*double(_i32radius)/ys +  double(_i32radius)/zs*double(_i32radius)/zs);
 
 	// memory allocation
 	// difficult to know how many voxels before computing,
-	// so the reserve for the BB
-	table.reserve(i32Width*i32Width*i32Width);
+	// so the reserve for the BB (not really a pb)
+	table.reserve(radX*radY*radZ*8);
 	table.clear();
 
 	// scan all the BB of the sphere
-	for (int z = -_i32radius;  z<=_i32radius; z++)
+	for (int z = -radZ; z<radZ; z++)
 	{
-		for (int y = -_i32radius;  y<=_i32radius; y++)
+		for (int y = -radY; y<radY; y++)
 		{
-			for (int x = -_i32radius;  x<=_i32radius; x++)
+			for (int x = -radX; x<radX; x++)
 			{
 				Geom::Vec3f v(float(x)*xs,float(y)*ys,float(z)*zs);
-				float fLength =  v.norm2();
+				float fLength =  v.norm();
 				// if inside the sphere
-				if (fLength<=fRad2)
+				if (fLength<=sRadius)
 				{
 					// the the index of the voxel
 					int index = z * m_WXY + y * m_WX + x;
@@ -551,12 +551,17 @@ void Image<DataType>::createMaskOffsetSphere(std::vector<int>& table, int _i32ra
 }
 
 
+//template<typename DataType>
+//float Image<DataType>::computeCurvatureCount(const DataType *ptrVox, const std::vector<int>& sphere, DataType val)
+
 template<typename DataType>
-float Image<DataType>::computeCurvatureCount(const DataType *ptrVox, const std::vector<int>& sphere, DataType val)
+float Image<DataType>::computeCurvatureCount(float x, float y, float z, const std::vector<int>& sphere, DataType val)
 {
+//	std::cout << "XYZ"<< x << " , "<<y<<" , "<<z <<"   => "<< int(round(x/m_SX)) << " , "<<int(round(y/m_SY))<<" , "<<int(round(z/m_SZ)) << std::endl;
+
 	int noir = 0;
 	int blanc = 0;
-
+	const DataType *ptrVox = this->getVoxelPtr(int(ceil(x)), int(ceil(y)), int(ceil(z)));
 
 	for (std::vector<int>::const_iterator it=sphere.begin(); it!=sphere.end();it++)
 	{
-- 
GitLab