nuclear@0: /* nuclear@0: --------------------------------------------------------------------------- nuclear@0: Open Asset Import Library (assimp) nuclear@0: --------------------------------------------------------------------------- nuclear@0: nuclear@0: Copyright (c) 2006-2012, assimp team nuclear@0: nuclear@0: All rights reserved. nuclear@0: nuclear@0: Redistribution and use of this software in source and binary forms, nuclear@0: with or without modification, are permitted provided that the following nuclear@0: conditions are met: nuclear@0: nuclear@0: * Redistributions of source code must retain the above nuclear@0: copyright notice, this list of conditions and the nuclear@0: following disclaimer. nuclear@0: nuclear@0: * Redistributions in binary form must reproduce the above nuclear@0: copyright notice, this list of conditions and the nuclear@0: following disclaimer in the documentation and/or other nuclear@0: materials provided with the distribution. nuclear@0: nuclear@0: * Neither the name of the assimp team, nor the names of its nuclear@0: contributors may be used to endorse or promote products nuclear@0: derived from this software without specific prior nuclear@0: written permission of the assimp team. nuclear@0: nuclear@0: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS nuclear@0: "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT nuclear@0: LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR nuclear@0: A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT nuclear@0: OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, nuclear@0: SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT nuclear@0: LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, nuclear@0: DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY nuclear@0: THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT nuclear@0: (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE nuclear@0: OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. nuclear@0: --------------------------------------------------------------------------- nuclear@0: */ nuclear@0: nuclear@0: /** @file Implementation of the helper class to quickly find vertices close to a given position */ nuclear@0: nuclear@0: #include "AssimpPCH.h" nuclear@0: #include "SpatialSort.h" nuclear@0: nuclear@0: using namespace Assimp; nuclear@0: nuclear@0: // CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8. nuclear@0: #ifndef CHAR_BIT nuclear@0: # define CHAR_BIT 8 nuclear@0: #endif nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: // Constructs a spatially sorted representation from the given position array. nuclear@0: SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions, nuclear@0: unsigned int pElementOffset) nuclear@0: nuclear@0: // define the reference plane. We choose some arbitrary vector away from all basic axises nuclear@0: // in the hope that no model spreads all its vertices along this plane. nuclear@0: : mPlaneNormal(0.8523f, 0.34321f, 0.5736f) nuclear@0: { nuclear@0: mPlaneNormal.Normalize(); nuclear@0: Fill(pPositions,pNumPositions,pElementOffset); nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: SpatialSort :: SpatialSort() nuclear@0: : mPlaneNormal(0.8523f, 0.34321f, 0.5736f) nuclear@0: { nuclear@0: mPlaneNormal.Normalize(); nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: // Destructor nuclear@0: SpatialSort::~SpatialSort() nuclear@0: { nuclear@0: // nothing to do here, everything destructs automatically nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions, nuclear@0: unsigned int pElementOffset, nuclear@0: bool pFinalize /*= true */) nuclear@0: { nuclear@0: mPositions.clear(); nuclear@0: Append(pPositions,pNumPositions,pElementOffset,pFinalize); nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: void SpatialSort :: Finalize() nuclear@0: { nuclear@0: std::sort( mPositions.begin(), mPositions.end()); nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions, nuclear@0: unsigned int pElementOffset, nuclear@0: bool pFinalize /*= true */) nuclear@0: { nuclear@0: // store references to all given positions along with their distance to the reference plane nuclear@0: const size_t initial = mPositions.size(); nuclear@0: mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2)); nuclear@0: for( unsigned int a = 0; a < pNumPositions; a++) nuclear@0: { nuclear@0: const char* tempPointer = reinterpret_cast (pPositions); nuclear@0: const aiVector3D* vec = reinterpret_cast (tempPointer + a * pElementOffset); nuclear@0: nuclear@0: // store position by index and distance nuclear@0: float distance = *vec * mPlaneNormal; nuclear@0: mPositions.push_back( Entry( a+initial, *vec, distance)); nuclear@0: } nuclear@0: nuclear@0: if (pFinalize) { nuclear@0: // now sort the array ascending by distance. nuclear@0: Finalize(); nuclear@0: } nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: // Returns an iterator for all positions close to the given position. nuclear@0: void SpatialSort::FindPositions( const aiVector3D& pPosition, nuclear@0: float pRadius, std::vector& poResults) const nuclear@0: { nuclear@0: const float dist = pPosition * mPlaneNormal; nuclear@0: const float minDist = dist - pRadius, maxDist = dist + pRadius; nuclear@0: nuclear@0: // clear the array in this strange fashion because a simple clear() would also deallocate nuclear@0: // the array which we want to avoid nuclear@0: poResults.erase( poResults.begin(), poResults.end()); nuclear@0: nuclear@0: // quick check for positions outside the range nuclear@0: if( mPositions.size() == 0) nuclear@0: return; nuclear@0: if( maxDist < mPositions.front().mDistance) nuclear@0: return; nuclear@0: if( minDist > mPositions.back().mDistance) nuclear@0: return; nuclear@0: nuclear@0: // do a binary search for the minimal distance to start the iteration there nuclear@0: unsigned int index = (unsigned int)mPositions.size() / 2; nuclear@0: unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4; nuclear@0: while( binaryStepSize > 1) nuclear@0: { nuclear@0: if( mPositions[index].mDistance < minDist) nuclear@0: index += binaryStepSize; nuclear@0: else nuclear@0: index -= binaryStepSize; nuclear@0: nuclear@0: binaryStepSize /= 2; nuclear@0: } nuclear@0: nuclear@0: // depending on the direction of the last step we need to single step a bit back or forth nuclear@0: // to find the actual beginning element of the range nuclear@0: while( index > 0 && mPositions[index].mDistance > minDist) nuclear@0: index--; nuclear@0: while( index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist) nuclear@0: index++; nuclear@0: nuclear@0: // Mow start iterating from there until the first position lays outside of the distance range. nuclear@0: // Add all positions inside the distance range within the given radius to the result aray nuclear@0: std::vector::const_iterator it = mPositions.begin() + index; nuclear@0: const float pSquared = pRadius*pRadius; nuclear@0: while( it->mDistance < maxDist) nuclear@0: { nuclear@0: if( (it->mPosition - pPosition).SquareLength() < pSquared) nuclear@0: poResults.push_back( it->mIndex); nuclear@0: ++it; nuclear@0: if( it == mPositions.end()) nuclear@0: break; nuclear@0: } nuclear@0: nuclear@0: // that's it nuclear@0: } nuclear@0: nuclear@0: namespace { nuclear@0: nuclear@0: // Binary, signed-integer representation of a single-precision floating-point value. nuclear@0: // IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are nuclear@0: // ordered the same way when their bits are reinterpreted as sign-magnitude integers." nuclear@0: // This allows us to convert all floating-point numbers to signed integers of arbitrary size nuclear@0: // and then use them to work with ULPs (Units in the Last Place, for high-precision nuclear@0: // computations) or to compare them (integer comparisons are faster than floating-point nuclear@0: // comparisons on many platforms). nuclear@0: typedef signed int BinFloat; nuclear@0: nuclear@0: // -------------------------------------------------------------------------------------------- nuclear@0: // Converts the bit pattern of a floating-point number to its signed integer representation. nuclear@0: BinFloat ToBinary( const float & pValue) { nuclear@0: nuclear@0: // If this assertion fails, signed int is not big enough to store a float on your platform. nuclear@0: // Please correct the declaration of BinFloat a few lines above - but do it in a portable, nuclear@0: // #ifdef'd manner! nuclear@0: BOOST_STATIC_ASSERT( sizeof(BinFloat) >= sizeof(float)); nuclear@0: nuclear@0: #if defined( _MSC_VER) nuclear@0: // If this assertion fails, Visual C++ has finally moved to ILP64. This means that this nuclear@0: // code has just become legacy code! Find out the current value of _MSC_VER and modify nuclear@0: // the #if above so it evaluates false on the current and all upcoming VC versions (or nuclear@0: // on the current platform, if LP64 or LLP64 are still used on other platforms). nuclear@0: BOOST_STATIC_ASSERT( sizeof(BinFloat) == sizeof(float)); nuclear@0: nuclear@0: // This works best on Visual C++, but other compilers have their problems with it. nuclear@0: const BinFloat binValue = reinterpret_cast(pValue); nuclear@0: #else nuclear@0: // On many compilers, reinterpreting a float address as an integer causes aliasing nuclear@0: // problems. This is an ugly but more or less safe way of doing it. nuclear@0: union { nuclear@0: float asFloat; nuclear@0: BinFloat asBin; nuclear@0: } conversion; nuclear@0: conversion.asBin = 0; // zero empty space in case sizeof(BinFloat) > sizeof(float) nuclear@0: conversion.asFloat = pValue; nuclear@0: const BinFloat binValue = conversion.asBin; nuclear@0: #endif nuclear@0: nuclear@0: // floating-point numbers are of sign-magnitude format, so find out what signed number nuclear@0: // representation we must convert negative values to. nuclear@0: // See http://en.wikipedia.org/wiki/Signed_number_representations. nuclear@0: nuclear@0: // Two's complement? nuclear@0: if( (-42 == (~42 + 1)) && (binValue & 0x80000000)) nuclear@0: return BinFloat(1 << (CHAR_BIT * sizeof(BinFloat) - 1)) - binValue; nuclear@0: // One's complement? nuclear@0: else if( (-42 == ~42) && (binValue & 0x80000000)) nuclear@0: return BinFloat(-0) - binValue; nuclear@0: // Sign-magnitude? nuclear@0: else if( (-42 == (42 | (-0))) && (binValue & 0x80000000)) // -0 = 1000... binary nuclear@0: return binValue; nuclear@0: else nuclear@0: return binValue; nuclear@0: } nuclear@0: nuclear@0: } // namespace nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: // Fills an array with indices of all positions indentical to the given position. In opposite to nuclear@0: // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units. nuclear@0: void SpatialSort::FindIdenticalPositions( const aiVector3D& pPosition, nuclear@0: std::vector& poResults) const nuclear@0: { nuclear@0: // Epsilons have a huge disadvantage: they are of constant precision, while floating-point nuclear@0: // values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but nuclear@0: // if you apply it to 0.001, it is enormous. nuclear@0: nuclear@0: // The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs nuclear@0: // tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of nuclear@0: // logarithmic precision - around 1, they are 1÷(2^24) and around 10000, they are 0.00125. nuclear@0: nuclear@0: // For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The nuclear@0: // incoming vertex positions might have already been transformed, probably using rather nuclear@0: // inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify nuclear@0: // identical vertex positions. nuclear@0: static const int toleranceInULPs = 4; nuclear@0: // An interesting point is that the inaccuracy grows linear with the number of operations: nuclear@0: // multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs nuclear@0: // plus 0.5 ULPs for the multiplication. nuclear@0: // To compute the distance to the plane, a dot product is needed - that is a multiplication and nuclear@0: // an addition on each number. nuclear@0: static const int distanceToleranceInULPs = toleranceInULPs + 1; nuclear@0: // The squared distance between two 3D vectors is computed the same way, but with an additional nuclear@0: // subtraction. nuclear@0: static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1; nuclear@0: nuclear@0: // Convert the plane distance to its signed integer representation so the ULPs tolerance can be nuclear@0: // applied. For some reason, VC won't optimize two calls of the bit pattern conversion. nuclear@0: const BinFloat minDistBinary = ToBinary( pPosition * mPlaneNormal) - distanceToleranceInULPs; nuclear@0: const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs; nuclear@0: nuclear@0: // clear the array in this strange fashion because a simple clear() would also deallocate nuclear@0: // the array which we want to avoid nuclear@0: poResults.erase( poResults.begin(), poResults.end()); nuclear@0: nuclear@0: // do a binary search for the minimal distance to start the iteration there nuclear@0: unsigned int index = (unsigned int)mPositions.size() / 2; nuclear@0: unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4; nuclear@0: while( binaryStepSize > 1) nuclear@0: { nuclear@0: // Ugly, but conditional jumps are faster with integers than with floats nuclear@0: if( minDistBinary > ToBinary(mPositions[index].mDistance)) nuclear@0: index += binaryStepSize; nuclear@0: else nuclear@0: index -= binaryStepSize; nuclear@0: nuclear@0: binaryStepSize /= 2; nuclear@0: } nuclear@0: nuclear@0: // depending on the direction of the last step we need to single step a bit back or forth nuclear@0: // to find the actual beginning element of the range nuclear@0: while( index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance) ) nuclear@0: index--; nuclear@0: while( index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance)) nuclear@0: index++; nuclear@0: nuclear@0: // Now start iterating from there until the first position lays outside of the distance range. nuclear@0: // Add all positions inside the distance range within the tolerance to the result aray nuclear@0: std::vector::const_iterator it = mPositions.begin() + index; nuclear@0: while( ToBinary(it->mDistance) < maxDistBinary) nuclear@0: { nuclear@0: if( distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength())) nuclear@0: poResults.push_back(it->mIndex); nuclear@0: ++it; nuclear@0: if( it == mPositions.end()) nuclear@0: break; nuclear@0: } nuclear@0: nuclear@0: // that's it nuclear@0: } nuclear@0: nuclear@0: // ------------------------------------------------------------------------------------------------ nuclear@0: unsigned int SpatialSort::GenerateMappingTable(std::vector& fill,float pRadius) const nuclear@0: { nuclear@0: fill.resize(mPositions.size(),UINT_MAX); nuclear@0: float dist, maxDist; nuclear@0: nuclear@0: unsigned int t=0; nuclear@0: const float pSquared = pRadius*pRadius; nuclear@0: for (size_t i = 0; i < mPositions.size();) { nuclear@0: dist = mPositions[i].mPosition * mPlaneNormal; nuclear@0: maxDist = dist + pRadius; nuclear@0: nuclear@0: fill[mPositions[i].mIndex] = t; nuclear@0: const aiVector3D& oldpos = mPositions[i].mPosition; nuclear@0: for (++i; i < fill.size() && mPositions[i].mDistance < maxDist nuclear@0: && (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i) nuclear@0: { nuclear@0: fill[mPositions[i].mIndex] = t; nuclear@0: } nuclear@0: ++t; nuclear@0: } nuclear@0: nuclear@0: #ifdef _DEBUG nuclear@0: nuclear@0: // debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1 nuclear@0: for (size_t i = 0; i < fill.size(); ++i) { nuclear@0: ai_assert(fill[i]