Merge pull request #4130 from MalcolmTyrrell/MalcolmTyrrell/spatialSortImprovements

SpatialSort improvements

code/Common/SpatialSort.cpp

@@ -58,14 +58,16 @@ const aiVector3D PlaneInit(0.8523f, 0.34321f, 0.5736f);
 // define the reference plane. We choose some arbitrary vector away from all basic axes
 // in the hope that no model spreads all its vertices along this plane.
 SpatialSort::SpatialSort(const aiVector3D *pPositions, unsigned int pNumPositions, unsigned int pElementOffset) :
-        mPlaneNormal(PlaneInit) {
+        mPlaneNormal(PlaneInit),
+        mFinalized(false) {
     mPlaneNormal.Normalize();
     Fill(pPositions, pNumPositions, pElementOffset);
 }
 
 // ------------------------------------------------------------------------------------------------
 SpatialSort::SpatialSort() :
-        mPlaneNormal(PlaneInit) {
+        mPlaneNormal(PlaneInit),
+        mFinalized(false) {
     mPlaneNormal.Normalize();
 }
 
@@ -80,28 +82,41 @@ void SpatialSort::Fill(const aiVector3D *pPositions, unsigned int pNumPositions,
         unsigned int pElementOffset,
         bool pFinalize /*= true */) {
     mPositions.clear();
+    mFinalized = false;
     Append(pPositions, pNumPositions, pElementOffset, pFinalize);
+    mFinalized = pFinalize;
+}
+
+// ------------------------------------------------------------------------------------------------
+ai_real SpatialSort::CalculateDistance(const aiVector3D &pPosition) const {
+    return (pPosition - mCentroid) * mPlaneNormal;
 }
 
 // ------------------------------------------------------------------------------------------------
 void SpatialSort::Finalize() {
+    const ai_real scale = 1.0f / mPositions.size();
+    for (unsigned int i = 0; i < mPositions.size(); i++) {
+        mCentroid += scale * mPositions[i].mPosition; 
+    }
+    for (unsigned int i = 0; i < mPositions.size(); i++) {
+        mPositions[i].mDistance = CalculateDistance(mPositions[i].mPosition);
+    }
     std::sort(mPositions.begin(), mPositions.end());
+    mFinalized = true;
 }
 
 // ------------------------------------------------------------------------------------------------
 void SpatialSort::Append(const aiVector3D *pPositions, unsigned int pNumPositions,
         unsigned int pElementOffset,
         bool pFinalize /*= true */) {
+    ai_assert(!mFinalized && "You cannot add positions to the SpatialSort object after it has been finalized.");
     // store references to all given positions along with their distance to the reference plane
     const size_t initial = mPositions.size();
-    mPositions.reserve(initial + (pFinalize ? pNumPositions : pNumPositions * 2));
+    mPositions.reserve(initial + pNumPositions);
     for (unsigned int a = 0; a < pNumPositions; a++) {
         const char *tempPointer = reinterpret_cast<const char *>(pPositions);
         const aiVector3D *vec = reinterpret_cast<const aiVector3D *>(tempPointer + a * pElementOffset);
-
-        // store position by index and distance
-        ai_real distance = *vec * mPlaneNormal;
-        mPositions.push_back(Entry(static_cast<unsigned int>(a + initial), *vec, distance));
+        mPositions.push_back(Entry(static_cast<unsigned int>(a + initial), *vec));
     }
 
     if (pFinalize) {
@@ -114,7 +129,8 @@ void SpatialSort::Append(const aiVector3D *pPositions, unsigned int pNumPosition
 // Returns an iterator for all positions close to the given position.
 void SpatialSort::FindPositions(const aiVector3D &pPosition,
         ai_real pRadius, std::vector<unsigned int> &poResults) const {
-    const ai_real dist = pPosition * mPlaneNormal;
+    ai_assert(mFinalized && "The SpatialSort object must be finalized before FindPositions can be called.");
+    const ai_real dist = CalculateDistance(pPosition);
     const ai_real minDist = dist - pRadius, maxDist = dist + pRadius;
 
     // clear the array
@@ -229,6 +245,7 @@ BinFloat ToBinary(const ai_real &pValue) {
 // Fills an array with indices of all positions identical to the given position. In opposite to
 // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units.
 void SpatialSort::FindIdenticalPositions(const aiVector3D &pPosition, std::vector<unsigned int> &poResults) const {
+    ai_assert(mFinalized && "The SpatialSort object must be finalized before FindIdenticalPositions can be called.");
     // Epsilons have a huge disadvantage: they are of constant precision, while floating-point
     //  values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but
     //  if you apply it to 0.001, it is enormous.
@@ -254,7 +271,7 @@ void SpatialSort::FindIdenticalPositions(const aiVector3D &pPosition, std::vecto
 
     // Convert the plane distance to its signed integer representation so the ULPs tolerance can be
     //  applied. For some reason, VC won't optimize two calls of the bit pattern conversion.
-    const BinFloat minDistBinary = ToBinary(pPosition * mPlaneNormal) - distanceToleranceInULPs;
+    const BinFloat minDistBinary = ToBinary(CalculateDistance(pPosition)) - distanceToleranceInULPs;
     const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs;
 
     // clear the array in this strange fashion because a simple clear() would also deallocate
@@ -297,13 +314,14 @@ void SpatialSort::FindIdenticalPositions(const aiVector3D &pPosition, std::vecto
 
 // ------------------------------------------------------------------------------------------------
 unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int> &fill, ai_real pRadius) const {
+    ai_assert(mFinalized && "The SpatialSort object must be finalized before GenerateMappingTable can be called.");
     fill.resize(mPositions.size(), UINT_MAX);
     ai_real dist, maxDist;
 
     unsigned int t = 0;
     const ai_real pSquared = pRadius * pRadius;
     for (size_t i = 0; i < mPositions.size();) {
-        dist = mPositions[i].mPosition * mPlaneNormal;
+        dist = (mPositions[i].mPosition - mCentroid) * mPlaneNormal;
         maxDist = dist + pRadius;
 
         fill[mPositions[i].mIndex] = t;

include/assimp/SpatialSort.h

@@ -51,6 +51,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include <assimp/types.h>
 #include <vector>
+#include <limits>
 
 namespace Assimp {
 
@@ -142,24 +143,35 @@ public:
             ai_real pRadius) const;
 
 protected:
-    /** Normal of the sorting plane, normalized. The center is always at (0, 0, 0) */
+    /** Return the distance to the sorting plane. */
+    ai_real CalculateDistance(const aiVector3D &pPosition) const;
+
+protected:
+    /** Normal of the sorting plane, normalized.
+     */
     aiVector3D mPlaneNormal;
 
+    /** The centroid of the positions, which is used as a point on the sorting plane
+     * when calculating distance. This value is calculated in Finalize.
+    */
+    aiVector3D mCentroid;
+
     /** An entry in a spatially sorted position array. Consists of a vertex index,
      * its position and its pre-calculated distance from the reference plane */
     struct Entry {
         unsigned int mIndex; ///< The vertex referred by this entry
         aiVector3D mPosition; ///< Position
-        ai_real mDistance; ///< Distance of this vertex to the sorting plane
+        /// Distance of this vertex to the sorting plane. This is set by Finalize.
+        ai_real mDistance; 
 
         Entry() AI_NO_EXCEPT
-                : mIndex(999999999),
+                : mIndex(std::numeric_limits<unsigned int>::max()),
                   mPosition(),
-                  mDistance(99999.) {
+                  mDistance(std::numeric_limits<ai_real>::max()) {
             // empty
         }
-        Entry(unsigned int pIndex, const aiVector3D &pPosition, ai_real pDistance) :
-                mIndex(pIndex), mPosition(pPosition), mDistance(pDistance) {
+        Entry(unsigned int pIndex, const aiVector3D &pPosition) :
+                mIndex(pIndex), mPosition(pPosition), mDistance(std::numeric_limits<ai_real>::max()) {
             // empty
         }
 
@@ -168,6 +180,9 @@ protected:
 
     // all positions, sorted by distance to the sorting plane
     std::vector<Entry> mPositions;
+
+    /// false until the Finalize method is called.
+    bool mFinalized;
 };
 
 } // end of namespace Assimp

test/unit/Common/utSpatialSort.cpp

@@ -82,3 +82,39 @@ TEST_F(utSpatialSort, findPositionsTest) {
     sSort.FindPositions(vecs[0], 0.01f, indices);
     EXPECT_EQ(1u, indices.size());
 }
+
+TEST_F(utSpatialSort, highlyDisplacedPositionsTest) {
+    // Make a cube of positions, and then query it using the SpatialSort object.
+    constexpr unsigned int verticesPerAxis = 10;
+    constexpr ai_real step = 0.001f;
+    // Note the large constant offset here.
+    constexpr ai_real offset = 5000.0f - (0.5f * verticesPerAxis * step);
+    constexpr unsigned int totalNumPositions = verticesPerAxis * verticesPerAxis * verticesPerAxis;
+    aiVector3D* positions = new aiVector3D[totalNumPositions];
+    for (unsigned int x = 0; x < verticesPerAxis; ++x) {
+        for (unsigned int y = 0; y < verticesPerAxis; ++y) {
+            for (unsigned int z = 0; z < verticesPerAxis; ++z) {
+                const unsigned int index = (x * verticesPerAxis * verticesPerAxis) + (y * verticesPerAxis) + z;
+                positions[index] = aiVector3D(offset + (x * step), offset + (y * step), offset + (z * step));
+            }
+        }
+    }
+
+    SpatialSort sSort;
+    sSort.Fill(positions, totalNumPositions, sizeof(aiVector3D));
+
+    // Enough to find a point and its 6 immediate neighbors, but not any other point.
+    const ai_real epsilon = 1.1f * step;
+    std::vector<unsigned int> indices;
+    // Iterate through the _interior_ points of the cube.
+    for (unsigned int x = 1; x < verticesPerAxis - 1; ++x) {
+        for (unsigned int y = 1; y < verticesPerAxis - 1; ++y) {
+            for (unsigned int z = 1; z < verticesPerAxis - 1; ++z) {
+                const unsigned int index = (x * verticesPerAxis * verticesPerAxis) + (y * verticesPerAxis) + z;
+                sSort.FindPositions(positions[index], epsilon, indices);
+                ASSERT_EQ(7u, indices.size());
+            }
+        }
+    }
+    delete[] positions;
+}