diff --git a/blurRelax.cpp b/blurRelax.cpp index 8077d70..9851aff 100644 --- a/blurRelax.cpp +++ b/blurRelax.cpp @@ -57,6 +57,7 @@ SOFTWARE. #include #include #include +#include #define CHECKSTAT(m) if (!status) {status.perror(m); return status;}; @@ -72,6 +73,9 @@ SOFTWARE. #define GB_PIN 1 #define GB_SLIDE 2 +#define SA_LAPLACIAN 0 +#define SA_TAUBIN 1 + #define DEFORMER_NAME "BlurRelax" // Double vs Float @@ -79,6 +83,7 @@ SOFTWARE. #define point_t MPoint #define pointArray_t MPointArray + void edgeProject( const float_t basePoints[][4], const std::vector &group, @@ -174,8 +179,8 @@ void quickLaplacianSmooth( const std::vector> &neighbors, const std::vector &valence, const std::vector &shiftVal, - const std::vector &shiftComp, - const std::vector &pinPoints + const std::vector &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> &neighbors, std::vector> &hardEdges, std::vector &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge - std::vector &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge std::vector &valence, // as float for vectorizing std::vector &pinPoints, std::vector &creaseCount, @@ -282,7 +280,8 @@ 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 &group, const std::vector &order, @@ -290,7 +289,6 @@ class BlurRelax : public MPxDeformerNode { const std::vector> &neighbors, const std::vector> &hardEdges, const std::vector &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge - const std::vector &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge const std::vector &valence, // as float for vectorizing const std::vector &pinPoints, const std::vector &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> neighbors; std::vector> hardEdges; std::vector shiftVal; // normally 0.5; but it's 0.25 if on a hard edge - std::vector shiftComp; // normally 0.5; but it's 0.75 if on a hard edge std::vector valence; // as float for vectorizing std::vector pinPoints; std::vector 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 weightVals; @@ -548,7 +575,6 @@ void BlurRelax::buildQuickData( std::vector> &neighbors, std::vector> &hardEdges, std::vector &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge - std::vector &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge std::vector &valence, // as float for vectorizing std::vector &pinPoints, std::vector &creaseCount, @@ -665,11 +691,8 @@ void BlurRelax::buildQuickData( // if a vert has a pinned neighbor, remove it std::vector rawShiftVal; - std::vector 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,7 +766,8 @@ 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 &group, const std::vector &order, @@ -754,11 +775,10 @@ void BlurRelax::quickRelax( const std::vector> &neighbors, const std::vector> &hardEdges, const std::vector &shiftVal, // normally 0.5, but it's 0.25 if on a hard edge - const std::vector &shiftComp, // normally 0.5, but it's 0.75 if on a hard edge const std::vector &valence, // as float for vectorizing const std::vector &pinPoints, const std::vector &creaseCount, - float_t(*verts)[4] // already resized + float_t(*verts)[4] ) { bool rpEdges = (borderBehavior == BB_SLIDE) || (hardEdgeBehavior == HB_SLIDE) || (groupEdgeBehavior == GB_SLIDE); std::vector 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; }