vrshoot

annotate libs/assimp/SpatialSort.cpp @ 0:b2f14e535253

find changesets by author, revision, files, or words in the commit message
initial commit
author John Tsiombikas <nuclear@member.fsf.org>
date Sat, 01 Feb 2014 19:58:19 +0200
parents
children
rev   line source
nuclear@0 1 /*
nuclear@0 2 ---------------------------------------------------------------------------
nuclear@0 3 Open Asset Import Library (assimp)
nuclear@0 4 ---------------------------------------------------------------------------
nuclear@0 5
nuclear@0 6 Copyright (c) 2006-2012, assimp team
nuclear@0 7
nuclear@0 8 All rights reserved.
nuclear@0 9
nuclear@0 10 Redistribution and use of this software in source and binary forms,
nuclear@0 11 with or without modification, are permitted provided that the following
nuclear@0 12 conditions are met:
nuclear@0 13
nuclear@0 14 * Redistributions of source code must retain the above
nuclear@0 15 copyright notice, this list of conditions and the
nuclear@0 16 following disclaimer.
nuclear@0 17
nuclear@0 18 * Redistributions in binary form must reproduce the above
nuclear@0 19 copyright notice, this list of conditions and the
nuclear@0 20 following disclaimer in the documentation and/or other
nuclear@0 21 materials provided with the distribution.
nuclear@0 22
nuclear@0 23 * Neither the name of the assimp team, nor the names of its
nuclear@0 24 contributors may be used to endorse or promote products
nuclear@0 25 derived from this software without specific prior
nuclear@0 26 written permission of the assimp team.
nuclear@0 27
nuclear@0 28 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
nuclear@0 29 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
nuclear@0 30 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
nuclear@0 31 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
nuclear@0 32 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
nuclear@0 33 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
nuclear@0 34 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
nuclear@0 35 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
nuclear@0 36 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
nuclear@0 37 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
nuclear@0 38 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
nuclear@0 39 ---------------------------------------------------------------------------
nuclear@0 40 */
nuclear@0 41
nuclear@0 42 /** @file Implementation of the helper class to quickly find vertices close to a given position */
nuclear@0 43
nuclear@0 44 #include "AssimpPCH.h"
nuclear@0 45 #include "SpatialSort.h"
nuclear@0 46
nuclear@0 47 using namespace Assimp;
nuclear@0 48
nuclear@0 49 // CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8.
nuclear@0 50 #ifndef CHAR_BIT
nuclear@0 51 # define CHAR_BIT 8
nuclear@0 52 #endif
nuclear@0 53
nuclear@0 54 // ------------------------------------------------------------------------------------------------
nuclear@0 55 // Constructs a spatially sorted representation from the given position array.
nuclear@0 56 SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions,
nuclear@0 57 unsigned int pElementOffset)
nuclear@0 58
nuclear@0 59 // define the reference plane. We choose some arbitrary vector away from all basic axises
nuclear@0 60 // in the hope that no model spreads all its vertices along this plane.
nuclear@0 61 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
nuclear@0 62 {
nuclear@0 63 mPlaneNormal.Normalize();
nuclear@0 64 Fill(pPositions,pNumPositions,pElementOffset);
nuclear@0 65 }
nuclear@0 66
nuclear@0 67 // ------------------------------------------------------------------------------------------------
nuclear@0 68 SpatialSort :: SpatialSort()
nuclear@0 69 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
nuclear@0 70 {
nuclear@0 71 mPlaneNormal.Normalize();
nuclear@0 72 }
nuclear@0 73
nuclear@0 74 // ------------------------------------------------------------------------------------------------
nuclear@0 75 // Destructor
nuclear@0 76 SpatialSort::~SpatialSort()
nuclear@0 77 {
nuclear@0 78 // nothing to do here, everything destructs automatically
nuclear@0 79 }
nuclear@0 80
nuclear@0 81 // ------------------------------------------------------------------------------------------------
nuclear@0 82 void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions,
nuclear@0 83 unsigned int pElementOffset,
nuclear@0 84 bool pFinalize /*= true */)
nuclear@0 85 {
nuclear@0 86 mPositions.clear();
nuclear@0 87 Append(pPositions,pNumPositions,pElementOffset,pFinalize);
nuclear@0 88 }
nuclear@0 89
nuclear@0 90 // ------------------------------------------------------------------------------------------------
nuclear@0 91 void SpatialSort :: Finalize()
nuclear@0 92 {
nuclear@0 93 std::sort( mPositions.begin(), mPositions.end());
nuclear@0 94 }
nuclear@0 95
nuclear@0 96 // ------------------------------------------------------------------------------------------------
nuclear@0 97 void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions,
nuclear@0 98 unsigned int pElementOffset,
nuclear@0 99 bool pFinalize /*= true */)
nuclear@0 100 {
nuclear@0 101 // store references to all given positions along with their distance to the reference plane
nuclear@0 102 const size_t initial = mPositions.size();
nuclear@0 103 mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2));
nuclear@0 104 for( unsigned int a = 0; a < pNumPositions; a++)
nuclear@0 105 {
nuclear@0 106 const char* tempPointer = reinterpret_cast<const char*> (pPositions);
nuclear@0 107 const aiVector3D* vec = reinterpret_cast<const aiVector3D*> (tempPointer + a * pElementOffset);
nuclear@0 108
nuclear@0 109 // store position by index and distance
nuclear@0 110 float distance = *vec * mPlaneNormal;
nuclear@0 111 mPositions.push_back( Entry( a+initial, *vec, distance));
nuclear@0 112 }
nuclear@0 113
nuclear@0 114 if (pFinalize) {
nuclear@0 115 // now sort the array ascending by distance.
nuclear@0 116 Finalize();
nuclear@0 117 }
nuclear@0 118 }
nuclear@0 119
nuclear@0 120 // ------------------------------------------------------------------------------------------------
nuclear@0 121 // Returns an iterator for all positions close to the given position.
nuclear@0 122 void SpatialSort::FindPositions( const aiVector3D& pPosition,
nuclear@0 123 float pRadius, std::vector<unsigned int>& poResults) const
nuclear@0 124 {
nuclear@0 125 const float dist = pPosition * mPlaneNormal;
nuclear@0 126 const float minDist = dist - pRadius, maxDist = dist + pRadius;
nuclear@0 127
nuclear@0 128 // clear the array in this strange fashion because a simple clear() would also deallocate
nuclear@0 129 // the array which we want to avoid
nuclear@0 130 poResults.erase( poResults.begin(), poResults.end());
nuclear@0 131
nuclear@0 132 // quick check for positions outside the range
nuclear@0 133 if( mPositions.size() == 0)
nuclear@0 134 return;
nuclear@0 135 if( maxDist < mPositions.front().mDistance)
nuclear@0 136 return;
nuclear@0 137 if( minDist > mPositions.back().mDistance)
nuclear@0 138 return;
nuclear@0 139
nuclear@0 140 // do a binary search for the minimal distance to start the iteration there
nuclear@0 141 unsigned int index = (unsigned int)mPositions.size() / 2;
nuclear@0 142 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
nuclear@0 143 while( binaryStepSize > 1)
nuclear@0 144 {
nuclear@0 145 if( mPositions[index].mDistance < minDist)
nuclear@0 146 index += binaryStepSize;
nuclear@0 147 else
nuclear@0 148 index -= binaryStepSize;
nuclear@0 149
nuclear@0 150 binaryStepSize /= 2;
nuclear@0 151 }
nuclear@0 152
nuclear@0 153 // depending on the direction of the last step we need to single step a bit back or forth
nuclear@0 154 // to find the actual beginning element of the range
nuclear@0 155 while( index > 0 && mPositions[index].mDistance > minDist)
nuclear@0 156 index--;
nuclear@0 157 while( index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist)
nuclear@0 158 index++;
nuclear@0 159
nuclear@0 160 // Mow start iterating from there until the first position lays outside of the distance range.
nuclear@0 161 // Add all positions inside the distance range within the given radius to the result aray
nuclear@0 162 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
nuclear@0 163 const float pSquared = pRadius*pRadius;
nuclear@0 164 while( it->mDistance < maxDist)
nuclear@0 165 {
nuclear@0 166 if( (it->mPosition - pPosition).SquareLength() < pSquared)
nuclear@0 167 poResults.push_back( it->mIndex);
nuclear@0 168 ++it;
nuclear@0 169 if( it == mPositions.end())
nuclear@0 170 break;
nuclear@0 171 }
nuclear@0 172
nuclear@0 173 // that's it
nuclear@0 174 }
nuclear@0 175
nuclear@0 176 namespace {
nuclear@0 177
nuclear@0 178 // Binary, signed-integer representation of a single-precision floating-point value.
nuclear@0 179 // IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are
nuclear@0 180 // ordered the same way when their bits are reinterpreted as sign-magnitude integers."
nuclear@0 181 // This allows us to convert all floating-point numbers to signed integers of arbitrary size
nuclear@0 182 // and then use them to work with ULPs (Units in the Last Place, for high-precision
nuclear@0 183 // computations) or to compare them (integer comparisons are faster than floating-point
nuclear@0 184 // comparisons on many platforms).
nuclear@0 185 typedef signed int BinFloat;
nuclear@0 186
nuclear@0 187 // --------------------------------------------------------------------------------------------
nuclear@0 188 // Converts the bit pattern of a floating-point number to its signed integer representation.
nuclear@0 189 BinFloat ToBinary( const float & pValue) {
nuclear@0 190
nuclear@0 191 // If this assertion fails, signed int is not big enough to store a float on your platform.
nuclear@0 192 // Please correct the declaration of BinFloat a few lines above - but do it in a portable,
nuclear@0 193 // #ifdef'd manner!
nuclear@0 194 BOOST_STATIC_ASSERT( sizeof(BinFloat) >= sizeof(float));
nuclear@0 195
nuclear@0 196 #if defined( _MSC_VER)
nuclear@0 197 // If this assertion fails, Visual C++ has finally moved to ILP64. This means that this
nuclear@0 198 // code has just become legacy code! Find out the current value of _MSC_VER and modify
nuclear@0 199 // the #if above so it evaluates false on the current and all upcoming VC versions (or
nuclear@0 200 // on the current platform, if LP64 or LLP64 are still used on other platforms).
nuclear@0 201 BOOST_STATIC_ASSERT( sizeof(BinFloat) == sizeof(float));
nuclear@0 202
nuclear@0 203 // This works best on Visual C++, but other compilers have their problems with it.
nuclear@0 204 const BinFloat binValue = reinterpret_cast<BinFloat const &>(pValue);
nuclear@0 205 #else
nuclear@0 206 // On many compilers, reinterpreting a float address as an integer causes aliasing
nuclear@0 207 // problems. This is an ugly but more or less safe way of doing it.
nuclear@0 208 union {
nuclear@0 209 float asFloat;
nuclear@0 210 BinFloat asBin;
nuclear@0 211 } conversion;
nuclear@0 212 conversion.asBin = 0; // zero empty space in case sizeof(BinFloat) > sizeof(float)
nuclear@0 213 conversion.asFloat = pValue;
nuclear@0 214 const BinFloat binValue = conversion.asBin;
nuclear@0 215 #endif
nuclear@0 216
nuclear@0 217 // floating-point numbers are of sign-magnitude format, so find out what signed number
nuclear@0 218 // representation we must convert negative values to.
nuclear@0 219 // See http://en.wikipedia.org/wiki/Signed_number_representations.
nuclear@0 220
nuclear@0 221 // Two's complement?
nuclear@0 222 if( (-42 == (~42 + 1)) && (binValue & 0x80000000))
nuclear@0 223 return BinFloat(1 << (CHAR_BIT * sizeof(BinFloat) - 1)) - binValue;
nuclear@0 224 // One's complement?
nuclear@0 225 else if( (-42 == ~42) && (binValue & 0x80000000))
nuclear@0 226 return BinFloat(-0) - binValue;
nuclear@0 227 // Sign-magnitude?
nuclear@0 228 else if( (-42 == (42 | (-0))) && (binValue & 0x80000000)) // -0 = 1000... binary
nuclear@0 229 return binValue;
nuclear@0 230 else
nuclear@0 231 return binValue;
nuclear@0 232 }
nuclear@0 233
nuclear@0 234 } // namespace
nuclear@0 235
nuclear@0 236 // ------------------------------------------------------------------------------------------------
nuclear@0 237 // Fills an array with indices of all positions indentical to the given position. In opposite to
nuclear@0 238 // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units.
nuclear@0 239 void SpatialSort::FindIdenticalPositions( const aiVector3D& pPosition,
nuclear@0 240 std::vector<unsigned int>& poResults) const
nuclear@0 241 {
nuclear@0 242 // Epsilons have a huge disadvantage: they are of constant precision, while floating-point
nuclear@0 243 // values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but
nuclear@0 244 // if you apply it to 0.001, it is enormous.
nuclear@0 245
nuclear@0 246 // The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs
nuclear@0 247 // tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of
nuclear@0 248 // logarithmic precision - around 1, they are 1÷(2^24) and around 10000, they are 0.00125.
nuclear@0 249
nuclear@0 250 // For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The
nuclear@0 251 // incoming vertex positions might have already been transformed, probably using rather
nuclear@0 252 // inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify
nuclear@0 253 // identical vertex positions.
nuclear@0 254 static const int toleranceInULPs = 4;
nuclear@0 255 // An interesting point is that the inaccuracy grows linear with the number of operations:
nuclear@0 256 // multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs
nuclear@0 257 // plus 0.5 ULPs for the multiplication.
nuclear@0 258 // To compute the distance to the plane, a dot product is needed - that is a multiplication and
nuclear@0 259 // an addition on each number.
nuclear@0 260 static const int distanceToleranceInULPs = toleranceInULPs + 1;
nuclear@0 261 // The squared distance between two 3D vectors is computed the same way, but with an additional
nuclear@0 262 // subtraction.
nuclear@0 263 static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1;
nuclear@0 264
nuclear@0 265 // Convert the plane distance to its signed integer representation so the ULPs tolerance can be
nuclear@0 266 // applied. For some reason, VC won't optimize two calls of the bit pattern conversion.
nuclear@0 267 const BinFloat minDistBinary = ToBinary( pPosition * mPlaneNormal) - distanceToleranceInULPs;
nuclear@0 268 const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs;
nuclear@0 269
nuclear@0 270 // clear the array in this strange fashion because a simple clear() would also deallocate
nuclear@0 271 // the array which we want to avoid
nuclear@0 272 poResults.erase( poResults.begin(), poResults.end());
nuclear@0 273
nuclear@0 274 // do a binary search for the minimal distance to start the iteration there
nuclear@0 275 unsigned int index = (unsigned int)mPositions.size() / 2;
nuclear@0 276 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
nuclear@0 277 while( binaryStepSize > 1)
nuclear@0 278 {
nuclear@0 279 // Ugly, but conditional jumps are faster with integers than with floats
nuclear@0 280 if( minDistBinary > ToBinary(mPositions[index].mDistance))
nuclear@0 281 index += binaryStepSize;
nuclear@0 282 else
nuclear@0 283 index -= binaryStepSize;
nuclear@0 284
nuclear@0 285 binaryStepSize /= 2;
nuclear@0 286 }
nuclear@0 287
nuclear@0 288 // depending on the direction of the last step we need to single step a bit back or forth
nuclear@0 289 // to find the actual beginning element of the range
nuclear@0 290 while( index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance) )
nuclear@0 291 index--;
nuclear@0 292 while( index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance))
nuclear@0 293 index++;
nuclear@0 294
nuclear@0 295 // Now start iterating from there until the first position lays outside of the distance range.
nuclear@0 296 // Add all positions inside the distance range within the tolerance to the result aray
nuclear@0 297 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
nuclear@0 298 while( ToBinary(it->mDistance) < maxDistBinary)
nuclear@0 299 {
nuclear@0 300 if( distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength()))
nuclear@0 301 poResults.push_back(it->mIndex);
nuclear@0 302 ++it;
nuclear@0 303 if( it == mPositions.end())
nuclear@0 304 break;
nuclear@0 305 }
nuclear@0 306
nuclear@0 307 // that's it
nuclear@0 308 }
nuclear@0 309
nuclear@0 310 // ------------------------------------------------------------------------------------------------
nuclear@0 311 unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int>& fill,float pRadius) const
nuclear@0 312 {
nuclear@0 313 fill.resize(mPositions.size(),UINT_MAX);
nuclear@0 314 float dist, maxDist;
nuclear@0 315
nuclear@0 316 unsigned int t=0;
nuclear@0 317 const float pSquared = pRadius*pRadius;
nuclear@0 318 for (size_t i = 0; i < mPositions.size();) {
nuclear@0 319 dist = mPositions[i].mPosition * mPlaneNormal;
nuclear@0 320 maxDist = dist + pRadius;
nuclear@0 321
nuclear@0 322 fill[mPositions[i].mIndex] = t;
nuclear@0 323 const aiVector3D& oldpos = mPositions[i].mPosition;
nuclear@0 324 for (++i; i < fill.size() && mPositions[i].mDistance < maxDist
nuclear@0 325 && (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i)
nuclear@0 326 {
nuclear@0 327 fill[mPositions[i].mIndex] = t;
nuclear@0 328 }
nuclear@0 329 ++t;
nuclear@0 330 }
nuclear@0 331
nuclear@0 332 #ifdef _DEBUG
nuclear@0 333
nuclear@0 334 // debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1
nuclear@0 335 for (size_t i = 0; i < fill.size(); ++i) {
nuclear@0 336 ai_assert(fill[i]<mPositions.size());
nuclear@0 337 }
nuclear@0 338
nuclear@0 339 #endif
nuclear@0 340 return t;
nuclear@0 341 }
nuclear@0 342