Merge pull request #6 from tbttfox/taubin

Implement Taubin Smoothing

Author: tbttfox

blurRelax.cpp

@@ -57,6 +57,7 @@ SOFTWARE.
 #include <vector>
 #include <algorithm>
 #include <numeric>
+#include <math.h>
 
 #define CHECKSTAT(m) if (!status) {status.perror(m); return status;};
 
@@ -72,13 +73,17 @@ SOFTWARE.
 #define GB_PIN    1
 #define GB_SLIDE  2
 
+#define SA_LAPLACIAN 0
+#define SA_TAUBIN 1
+
 #define DEFORMER_NAME "BlurRelax"
 
 // Double vs Float
 #define float_t double
 #define point_t MPoint
 #define pointArray_t MPointArray
 
+
 void edgeProject(
 	const float_t basePoints[][4],
 	const std::vector<size_t> &group,
@@ -174,8 +179,8 @@ void quickLaplacianSmooth(
 	const std::vector<std::vector<size_t>> &neighbors,
 	const std::vector<float_t> &valence,
 	const std::vector<float_t> &shiftVal,
-	const std::vector<float_t> &shiftComp,
-	const std::vector<bool> &pinPoints
+	const std::vector<bool> &pinPoints,
+	const float_t taubinBias=1.0
 ) {
 	/*
 	All the crazy hoops I've jumped through are to make auto-vectorization work
@@ -184,23 +189,19 @@ void quickLaplacianSmooth(
 	of a std::vector of neighbors per vert with that vector sorted so the verts
 	with the most neighbors were at the top.
 	neighborOffsets contains the index offsets of vertices with at least [index] of neighbors
-	So a single-subdivided cube would have 8 3-valence, and 18 4-valence for a total of 26 verts, and 96 neighbors
-	So the neighborOffsets would be [0, 26, 52, 78, 96] meaning that
-	verts 0-25 have at least 1 neighbor
-	verts 26-51 have at least 2 neighbors
-	verts 52-77 have at least 3 neighbors
-	verts 78-95 have at least 4 neighbors
 	*/
 
-	// First, get verts as a single pointer to the contiguous memory stored in (verts*)[4]
+	// First, get verts as a single pointer to the contiguous memory stored in (verts2d*)[4]
 	float_t* verts = &(verts2d[0][0]);
 
 	// number of nonzero valence
 	size_t nzv = neighbors[0].size();
+
+	// number of nonzero components
 	size_t nzc = 4 * nzv;
 
 	// The __restrict keyword tells the compiler that *outComp
-	// are not pointed to by any other pointer in this scope
+	// is not pointed to by any other pointer in this scope
 	// This allows for auto-vectorization
 	float_t * __restrict outComp = new float_t[nzc];
 	memset(outComp, 0, nzc*sizeof(float_t));
@@ -216,10 +217,8 @@ void quickLaplacianSmooth(
 		}
 	}
 
-	// Depending on the compiler optimization, it may be faster to break up this line
-	// Gotta test
 	for (size_t i = 0; i < nzc; ++i) {
-		outComp[i] = shiftVal[i] * (outComp[i] / valence[i]) + shiftComp[i] * verts[i];
+		outComp[i] = shiftVal[i] * taubinBias * ((outComp[i] / valence[i]) - verts[i]) + verts[i];
 	}
 
 	memcpy(verts, outComp, nzc*sizeof(float_t));
@@ -244,7 +243,7 @@ class BlurRelax : public MPxDeformerNode {
 		static MObject aHardEdgeBehavior;
 		static MObject aGroupEdgeBehavior;
 		static MObject aReproject;
-
+		static MObject aTaubinBias;
 		static MTypeId id;
 	private:
 
@@ -269,7 +268,6 @@ class BlurRelax : public MPxDeformerNode {
 			std::vector<std::vector<size_t>> &neighbors,
 			std::vector<std::vector<bool>> &hardEdges,
 			std::vector<float_t> &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge
-			std::vector<float_t> &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge
 			std::vector<float_t> &valence, // as float for vectorizing
 			std::vector<bool> &pinPoints,
 			std::vector<UINT> &creaseCount,
@@ -282,15 +280,15 @@ class BlurRelax : public MPxDeformerNode {
 			const short hardEdgeBehavior,
 			const short groupEdgeBehavior,
 			const bool reproject,
-			const UINT iterations,
+			const float taubinBias,
+			const float_t iterations,
 			const UINT numVerts,
 			const std::vector<bool> &group,
 			const std::vector<size_t> &order,
 			const std::vector<size_t> &invOrder,
 			const std::vector<std::vector<size_t>> &neighbors,
 			const std::vector<std::vector<bool>> &hardEdges,
 			const std::vector<float_t> &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge
-			const std::vector<float_t> &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge
 			const std::vector<float_t> &valence, // as float for vectorizing
 			const std::vector<bool> &pinPoints,
 			const std::vector<UINT> &creaseCount,
@@ -306,6 +304,8 @@ MObject BlurRelax::aBorderBehavior;
 MObject BlurRelax::aHardEdgeBehavior;
 MObject BlurRelax::aGroupEdgeBehavior;
 MObject BlurRelax::aReproject;
+MObject BlurRelax::aTaubinBias;
+
 
 BlurRelax::BlurRelax() {}
 BlurRelax::~BlurRelax() {}
@@ -423,9 +423,20 @@ MStatus BlurRelax::initialize() {
 	status = attributeAffects(aReproject, outputGeom);
 	CHECKSTAT("aReproject");
 
-	aIterations = nAttr.create("iterations", "i", MFnNumericData::kInt, 10, &status);
+	// Taubin Bias is divided by 1000 internally
+	aTaubinBias = nAttr.create("preserveVolume", "pv", MFnNumericData::kFloat, 0.0f, &status);
+	CHECKSTAT("aTaubinBias");
+	nAttr.setMin(0.0f);
+	nAttr.setMax(2.0f);
+	nAttr.setChannelBox(true);
+	status = addAttribute(aTaubinBias);
+	CHECKSTAT("aTaubinBias");
+	status = attributeAffects(aTaubinBias, outputGeom);
+	CHECKSTAT("aTaubinBias");
+
+	aIterations = nAttr.create("iterations", "i", MFnNumericData::kFloat, 0, &status);
 	CHECKSTAT("aIterations");
-	nAttr.setMin(0);
+	nAttr.setMin(0.0);
 	nAttr.setChannelBox(true);
 	status = addAttribute(aIterations);
 	CHECKSTAT("aIterations");
@@ -444,7 +455,7 @@ MStatus BlurRelax::deform(MDataBlock& dataBlock, MItGeometry& vertIter, const MM
 	float env = hEnv.asFloat();
 
 	MDataHandle hIter = dataBlock.inputValue(aIterations);
-	int iterations = hIter.asInt();
+	float iterations = hIter.asFloat();
 
 	if (iterations > 0 && env > 0.0f){
 		// Get the data from the node
@@ -457,6 +468,23 @@ MStatus BlurRelax::deform(MDataBlock& dataBlock, MItGeometry& vertIter, const MM
 		MDataHandle hReproject = dataBlock.inputValue(aReproject);
 		bool reproject = hReproject.asBool();
 
+		MDataHandle hTBias = dataBlock.inputValue(aTaubinBias);
+		float tBias = hTBias.asFloat();
+		// volume preservation uses 2 steps per iteration
+		// so half the number of iterations if I'm volume preserving
+		// The iterations interpolating as floats takes care of 99% of the jumping
+		// There *will* be a tiny jump from 0.0 on preserve volume if
+		// reprojection is also turned on. I don't think I can get rid of that one
+		if (tBias > 0.0f) {
+			iterations /= 2.0f;
+		}
+		// So the numbers shown are intuitive to the user
+		// 0 maps to 1.0 and 1 maps to -1.05
+		// -1.05 is used because taubin smoothing needs something
+		// just a little over 1.0 to truly preserve volume,
+		// and that value looked good on my test mesh.
+		tBias = -2.05f * tBias + 1.0f;
+
 		// get the input mesh corresponding to this output
 		MObject thisNode = this->thisMObject();
 		MPlug inPlug(thisNode, input);
@@ -471,7 +499,6 @@ MStatus BlurRelax::deform(MDataBlock& dataBlock, MItGeometry& vertIter, const MM
 		std::vector<std::vector<size_t>> neighbors;
 		std::vector<std::vector<bool>> hardEdges;
 		std::vector<float_t> shiftVal; // normally 0.5; but it's 0.25 if on a hard edge
-		std::vector<float_t> shiftComp; // normally 0.5; but it's 0.75 if on a hard edge
 		std::vector<float_t> valence; // as float for vectorizing
 		std::vector<bool> pinPoints;
 		std::vector<UINT> creaseCount;
@@ -489,7 +516,7 @@ MStatus BlurRelax::deform(MDataBlock& dataBlock, MItGeometry& vertIter, const MM
 		// Populate the variables with *SPECIALLY ORDERED* data
 		// all vertex data is now shuffled by the order vector
 		buildQuickData(mesh, vertIter, bb, hb, gb, reproject,
-			group, order, invOrder, neighbors, hardEdges, shiftVal, shiftComp, valence,
+			group, order, invOrder, neighbors, hardEdges, shiftVal, valence,
 			pinPoints, creaseCount, verts);
 
 		// This can happen if the user is pinning all the points
@@ -500,9 +527,9 @@ MStatus BlurRelax::deform(MDataBlock& dataBlock, MItGeometry& vertIter, const MM
 		}
 
 		// Calculate the relax, and store in verts
-		quickRelax(mesh, bb, hb, gb, reproject,
+		quickRelax(mesh, bb, hb, gb, reproject, tBias,
 			iterations, numVerts, group, order, invOrder, neighbors, hardEdges, shiftVal,
-			shiftComp, valence, pinPoints, creaseCount, verts);
+			valence, pinPoints, creaseCount, verts);
 
 		// Get the painted weight values
 		std::vector<float> weightVals;
@@ -548,7 +575,6 @@ void BlurRelax::buildQuickData(
 	std::vector<std::vector<size_t>> &neighbors,
 	std::vector<std::vector<bool>> &hardEdges,
 	std::vector<float_t> &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge
-	std::vector<float_t> &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge
 	std::vector<float_t> &valence, // as float for vectorizing
 	std::vector<bool> &pinPoints,
 	std::vector<UINT> &creaseCount,
@@ -665,11 +691,8 @@ void BlurRelax::buildQuickData(
 	// if a vert has a pinned neighbor, remove it
 
 	std::vector<float_t> rawShiftVal;
-	std::vector<float_t> rawShiftComp;
 	rawShiftVal.resize(numVertices);
-	rawShiftComp.resize(numVertices);
 	std::fill(rawShiftVal.begin(), rawShiftVal.end(), 0.5);
-	std::fill(rawShiftComp.begin(), rawShiftComp.end(), 0.5);
 
 	for (size_t i = 0; i < rawNeighbors.size(); ++i) {
 		if ((rawCreaseCount[i] != 0) || rawPinPoints[i]) {
@@ -689,7 +712,6 @@ void BlurRelax::buildQuickData(
 			rawNeighbors[i] = newNeigh;
 			rawHardEdges[i] = newHard;
 			rawShiftVal[i] = 0.25;
-			rawShiftComp[i] = 0.75;
 		}
 	}
 
@@ -709,7 +731,6 @@ void BlurRelax::buildQuickData(
 
 	valence.resize(numVertices*4);
 	shiftVal.resize(numVertices*4);
-	shiftComp.resize(numVertices*4);
 
 	invOrder.resize(order.size());
 	for (size_t i = 0; i < order.size(); ++i) {
@@ -734,7 +755,6 @@ void BlurRelax::buildQuickData(
 		for (size_t xx = 0; xx < 4; ++xx) {
 			valence[4 * i + xx] = float_t(vale);
 			shiftVal[4 * i + xx] = rawShiftVal[order[i]];
-			shiftComp[4 * i + xx] = rawShiftComp[order[i]];
 		}
 	}
 	delete [] rawVerts;
@@ -746,19 +766,19 @@ void BlurRelax::quickRelax(
 	const short hardEdgeBehavior,
 	const short groupEdgeBehavior,
 	const bool reproject,
-	const UINT iterations,
+	const float taubinBias,
+	const float_t iterations,
 	const UINT numVerts,
 	const std::vector<bool> &group,
 	const std::vector<size_t> &order,
 	const std::vector<size_t> &invOrder,
 	const std::vector<std::vector<size_t>> &neighbors,
 	const std::vector<std::vector<bool>> &hardEdges,
 	const std::vector<float_t> &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge
-	const std::vector<float_t> &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge
 	const std::vector<float_t> &valence, // as float for vectorizing
 	const std::vector<bool> &pinPoints,
 	const std::vector<UINT> &creaseCount,
-	float_t(*verts)[4] // already resized
+	float_t(*verts)[4]
 ) {
 	bool rpEdges = (borderBehavior == BB_SLIDE) || (hardEdgeBehavior == HB_SLIDE) || (groupEdgeBehavior == GB_SLIDE);
 	std::vector<size_t> groupIdxs;
@@ -773,6 +793,16 @@ void BlurRelax::quickRelax(
 		memcpy(&(baseVerts[0][0]), &(verts[0][0]), 4 * numVerts * sizeof(float_t));
 	}
 
+	float_t(*prevVerts)[4];
+	prevVerts = new float_t[numVerts][4];
+
+	float_t iterT, iterFI;
+	iterT = modf(iterations, &iterFI);
+	UINT iterI = (UINT)iterFI;
+	if (iterT > 0.0) {
+		iterI += 1;
+	}
+	
 	size_t nonzeroValence = neighbors[0].size();
 
 	MStatus status;
@@ -790,11 +820,20 @@ void BlurRelax::quickRelax(
 		octree.create(smoothMesh);
 	}
 
-	for (size_t r = 0; r < iterations; ++r) {
-		quickLaplacianSmooth(verts, numVerts, neighbors, valence, shiftVal, shiftComp, pinPoints);
+	for (size_t r = 0; r < iterI; ++r) {
+		if ((r == iterI - 1) && (iterT > 0.0)){
+			// Store the next-to-last iteration to interpolate with
+			memcpy(&(prevVerts[0][0]), &(verts[0][0]), 4 * numVerts * sizeof(float_t));
+		}
+		quickLaplacianSmooth(verts, numVerts, neighbors, valence, shiftVal, pinPoints);
+		if (taubinBias < 1.0){
+			quickLaplacianSmooth(verts, numVerts, neighbors, valence, shiftVal, pinPoints, taubinBias);
+		}
+
 		if (rpEdges) {
 			edgeProject(baseVerts, groupIdxs, invOrder, neighbors, hardEdges, creaseCount, verts);
 		}
+
 		if (reproject) {
 			#pragma omp parallel for if(numVerts>2000)
 			for (int i = 0; i < nonzeroValence; ++i) {
@@ -812,7 +851,17 @@ void BlurRelax::quickRelax(
 		}
 	}
 
-	if (rpEdges) delete [] baseVerts;
+	// Interpolate between prevVerts and verts based on iterT
+	if (iterT > 0.0) {
+		// This should vectorize
+		float_t * vv = &verts[0][0];
+		float_t * pv = &prevVerts[0][0];
+		for (size_t i = 0; i < numVerts * 4; ++i) {
+			vv[i] = ((vv[i] - pv[i]) * iterT) + pv[i];
+		}
+	}
 
+	if (rpEdges) delete[] baseVerts;
+	delete[] prevVerts;
 }