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