1 | /* |
---|---|

2 | --------------------------------------------------------------------------- |

3 | Open Asset Import Library (assimp) |

4 | --------------------------------------------------------------------------- |

5 | |

6 | Copyright (c) 2006-2017, assimp team |

7 | |

8 | |

9 | All rights reserved. |

10 | |

11 | Redistribution and use of this software in source and binary forms, |

12 | with or without modification, are permitted provided that the following |

13 | conditions are met: |

14 | |

15 | * Redistributions of source code must retain the above |

16 | copyright notice, this list of conditions and the |

17 | following disclaimer. |

18 | |

19 | * Redistributions in binary form must reproduce the above |

20 | copyright notice, this list of conditions and the |

21 | following disclaimer in the documentation and/or other |

22 | materials provided with the distribution. |

23 | |

24 | * Neither the name of the assimp team, nor the names of its |

25 | contributors may be used to endorse or promote products |

26 | derived from this software without specific prior |

27 | written permission of the assimp team. |

28 | |

29 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |

30 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |

31 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |

32 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT |

33 | OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, |

34 | SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT |

35 | LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |

36 | DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |

37 | THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |

38 | (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE |

39 | OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |

40 | --------------------------------------------------------------------------- |

41 | */ |

42 | |

43 | /** @file Implementation of the helper class to quickly find vertices close to a given position */ |

44 | |

45 | #include "SpatialSort.h" |

46 | #include <assimp/ai_assert.h> |

47 | |

48 | using namespace Assimp; |

49 | |

50 | // CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8. |

51 | #ifndef CHAR_BIT |

52 | # define CHAR_BIT 8 |

53 | #endif |

54 | |

55 | // ------------------------------------------------------------------------------------------------ |

56 | // Constructs a spatially sorted representation from the given position array. |

57 | SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions, |

58 | unsigned int pElementOffset) |

59 | |

60 | // define the reference plane. We choose some arbitrary vector away from all basic axises |

61 | // in the hope that no model spreads all its vertices along this plane. |

62 | : mPlaneNormal(0.8523f, 0.34321f, 0.5736f) |

63 | { |

64 | mPlaneNormal.Normalize(); |

65 | Fill(pPositions,pNumPositions,pElementOffset); |

66 | } |

67 | |

68 | // ------------------------------------------------------------------------------------------------ |

69 | SpatialSort :: SpatialSort() |

70 | : mPlaneNormal(0.8523f, 0.34321f, 0.5736f) |

71 | { |

72 | mPlaneNormal.Normalize(); |

73 | } |

74 | |

75 | // ------------------------------------------------------------------------------------------------ |

76 | // Destructor |

77 | SpatialSort::~SpatialSort() |

78 | { |

79 | // nothing to do here, everything destructs automatically |

80 | } |

81 | |

82 | // ------------------------------------------------------------------------------------------------ |

83 | void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions, |

84 | unsigned int pElementOffset, |

85 | bool pFinalize /*= true */) |

86 | { |

87 | mPositions.clear(); |

88 | Append(pPositions,pNumPositions,pElementOffset,pFinalize); |

89 | } |

90 | |

91 | // ------------------------------------------------------------------------------------------------ |

92 | void SpatialSort :: Finalize() |

93 | { |

94 | std::sort( mPositions.begin(), mPositions.end()); |

95 | } |

96 | |

97 | // ------------------------------------------------------------------------------------------------ |

98 | void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions, |

99 | unsigned int pElementOffset, |

100 | bool pFinalize /*= true */) |

101 | { |

102 | // store references to all given positions along with their distance to the reference plane |

103 | const size_t initial = mPositions.size(); |

104 | mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2)); |

105 | for( unsigned int a = 0; a < pNumPositions; a++) |

106 | { |

107 | const char* tempPointer = reinterpret_cast<const char*> (pPositions); |

108 | const aiVector3D* vec = reinterpret_cast<const aiVector3D*> (tempPointer + a * pElementOffset); |

109 | |

110 | // store position by index and distance |

111 | ai_real distance = *vec * mPlaneNormal; |

112 | mPositions.push_back( Entry( static_cast<unsigned int>(a+initial), *vec, distance)); |

113 | } |

114 | |

115 | if (pFinalize) { |

116 | // now sort the array ascending by distance. |

117 | Finalize(); |

118 | } |

119 | } |

120 | |

121 | // ------------------------------------------------------------------------------------------------ |

122 | // Returns an iterator for all positions close to the given position. |

123 | void SpatialSort::FindPositions( const aiVector3D& pPosition, |

124 | ai_real pRadius, std::vector<unsigned int>& poResults) const |

125 | { |

126 | const ai_real dist = pPosition * mPlaneNormal; |

127 | const ai_real minDist = dist - pRadius, maxDist = dist + pRadius; |

128 | |

129 | // clear the array |

130 | poResults.clear(); |

131 | |

132 | // quick check for positions outside the range |

133 | if( mPositions.size() == 0) |

134 | return; |

135 | if( maxDist < mPositions.front().mDistance) |

136 | return; |

137 | if( minDist > mPositions.back().mDistance) |

138 | return; |

139 | |

140 | // do a binary search for the minimal distance to start the iteration there |

141 | unsigned int index = (unsigned int)mPositions.size() / 2; |

142 | unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4; |

143 | while( binaryStepSize > 1) |

144 | { |

145 | if( mPositions[index].mDistance < minDist) |

146 | index += binaryStepSize; |

147 | else |

148 | index -= binaryStepSize; |

149 | |

150 | binaryStepSize /= 2; |

151 | } |

152 | |

153 | // depending on the direction of the last step we need to single step a bit back or forth |

154 | // to find the actual beginning element of the range |

155 | while( index > 0 && mPositions[index].mDistance > minDist) |

156 | index--; |

157 | while( index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist) |

158 | index++; |

159 | |

160 | // Mow start iterating from there until the first position lays outside of the distance range. |

161 | // Add all positions inside the distance range within the given radius to the result aray |

162 | std::vector<Entry>::const_iterator it = mPositions.begin() + index; |

163 | const ai_real pSquared = pRadius*pRadius; |

164 | while( it->mDistance < maxDist) |

165 | { |

166 | if( (it->mPosition - pPosition).SquareLength() < pSquared) |

167 | poResults.push_back( it->mIndex); |

168 | ++it; |

169 | if( it == mPositions.end()) |

170 | break; |

171 | } |

172 | |

173 | // that's it |

174 | } |

175 | |

176 | namespace { |

177 | |

178 | // Binary, signed-integer representation of a single-precision floating-point value. |

179 | // IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are |

180 | // ordered the same way when their bits are reinterpreted as sign-magnitude integers." |

181 | // This allows us to convert all floating-point numbers to signed integers of arbitrary size |

182 | // and then use them to work with ULPs (Units in the Last Place, for high-precision |

183 | // computations) or to compare them (integer comparisons are faster than floating-point |

184 | // comparisons on many platforms). |

185 | typedef ai_int BinFloat; |

186 | |

187 | // -------------------------------------------------------------------------------------------- |

188 | // Converts the bit pattern of a floating-point number to its signed integer representation. |

189 | BinFloat ToBinary( const ai_real & pValue) { |

190 | |

191 | // If this assertion fails, signed int is not big enough to store a float on your platform. |

192 | // Please correct the declaration of BinFloat a few lines above - but do it in a portable, |

193 | // #ifdef'd manner! |

194 | static_assert( sizeof(BinFloat) >= sizeof(ai_real), "sizeof(BinFloat) >= sizeof(ai_real)"); |

195 | |

196 | #if defined( _MSC_VER) |

197 | // If this assertion fails, Visual C++ has finally moved to ILP64. This means that this |

198 | // code has just become legacy code! Find out the current value of _MSC_VER and modify |

199 | // the #if above so it evaluates false on the current and all upcoming VC versions (or |

200 | // on the current platform, if LP64 or LLP64 are still used on other platforms). |

201 | static_assert( sizeof(BinFloat) == sizeof(ai_real), "sizeof(BinFloat) == sizeof(ai_real)"); |

202 | |

203 | // This works best on Visual C++, but other compilers have their problems with it. |

204 | const BinFloat binValue = reinterpret_cast<BinFloat const &>(pValue); |

205 | #else |

206 | // On many compilers, reinterpreting a float address as an integer causes aliasing |

207 | // problems. This is an ugly but more or less safe way of doing it. |

208 | union { |

209 | ai_real asFloat; |

210 | BinFloat asBin; |

211 | } conversion; |

212 | conversion.asBin = 0; // zero empty space in case sizeof(BinFloat) > sizeof(float) |

213 | conversion.asFloat = pValue; |

214 | const BinFloat binValue = conversion.asBin; |

215 | #endif |

216 | |

217 | // floating-point numbers are of sign-magnitude format, so find out what signed number |

218 | // representation we must convert negative values to. |

219 | // See http://en.wikipedia.org/wiki/Signed_number_representations. |

220 | |

221 | // Two's complement? |

222 | if( (-42 == (~42 + 1)) && (binValue & 0x80000000)) |

223 | return BinFloat(1 << (CHAR_BIT * sizeof(BinFloat) - 1)) - binValue; |

224 | // One's complement? |

225 | else if( (-42 == ~42) && (binValue & 0x80000000)) |

226 | return BinFloat(-0) - binValue; |

227 | // Sign-magnitude? |

228 | else if( (-42 == (42 | (-0))) && (binValue & 0x80000000)) // -0 = 1000... binary |

229 | return binValue; |

230 | else |

231 | return binValue; |

232 | } |

233 | |

234 | } // namespace |

235 | |

236 | // ------------------------------------------------------------------------------------------------ |

237 | // Fills an array with indices of all positions identical to the given position. In opposite to |

238 | // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units. |

239 | void SpatialSort::FindIdenticalPositions( const aiVector3D& pPosition, |

240 | std::vector<unsigned int>& poResults) const |

241 | { |

242 | // Epsilons have a huge disadvantage: they are of constant precision, while floating-point |

243 | // values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but |

244 | // if you apply it to 0.001, it is enormous. |

245 | |

246 | // The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs |

247 | // tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of |

248 | // logarithmic precision - around 1, they are 1*(2^24) and around 10000, they are 0.00125. |

249 | |

250 | // For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The |

251 | // incoming vertex positions might have already been transformed, probably using rather |

252 | // inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify |

253 | // identical vertex positions. |

254 | static const int toleranceInULPs = 4; |

255 | // An interesting point is that the inaccuracy grows linear with the number of operations: |

256 | // multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs |

257 | // plus 0.5 ULPs for the multiplication. |

258 | // To compute the distance to the plane, a dot product is needed - that is a multiplication and |

259 | // an addition on each number. |

260 | static const int distanceToleranceInULPs = toleranceInULPs + 1; |

261 | // The squared distance between two 3D vectors is computed the same way, but with an additional |

262 | // subtraction. |

263 | static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1; |

264 | |

265 | // Convert the plane distance to its signed integer representation so the ULPs tolerance can be |

266 | // applied. For some reason, VC won't optimize two calls of the bit pattern conversion. |

267 | const BinFloat minDistBinary = ToBinary( pPosition * mPlaneNormal) - distanceToleranceInULPs; |

268 | const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs; |

269 | |

270 | // clear the array in this strange fashion because a simple clear() would also deallocate |

271 | // the array which we want to avoid |

272 | poResults.resize( 0 ); |

273 | |

274 | // do a binary search for the minimal distance to start the iteration there |

275 | unsigned int index = (unsigned int)mPositions.size() / 2; |

276 | unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4; |

277 | while( binaryStepSize > 1) |

278 | { |

279 | // Ugly, but conditional jumps are faster with integers than with floats |

280 | if( minDistBinary > ToBinary(mPositions[index].mDistance)) |

281 | index += binaryStepSize; |

282 | else |

283 | index -= binaryStepSize; |

284 | |

285 | binaryStepSize /= 2; |

286 | } |

287 | |

288 | // depending on the direction of the last step we need to single step a bit back or forth |

289 | // to find the actual beginning element of the range |

290 | while( index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance) ) |

291 | index--; |

292 | while( index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance)) |

293 | index++; |

294 | |

295 | // Now start iterating from there until the first position lays outside of the distance range. |

296 | // Add all positions inside the distance range within the tolerance to the result aray |

297 | std::vector<Entry>::const_iterator it = mPositions.begin() + index; |

298 | while( ToBinary(it->mDistance) < maxDistBinary) |

299 | { |

300 | if( distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength())) |

301 | poResults.push_back(it->mIndex); |

302 | ++it; |

303 | if( it == mPositions.end()) |

304 | break; |

305 | } |

306 | |

307 | // that's it |

308 | } |

309 | |

310 | // ------------------------------------------------------------------------------------------------ |

311 | unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int>& fill, ai_real pRadius) const |

312 | { |

313 | fill.resize(mPositions.size(),UINT_MAX); |

314 | ai_real dist, maxDist; |

315 | |

316 | unsigned int t=0; |

317 | const ai_real pSquared = pRadius*pRadius; |

318 | for (size_t i = 0; i < mPositions.size();) { |

319 | dist = mPositions[i].mPosition * mPlaneNormal; |

320 | maxDist = dist + pRadius; |

321 | |

322 | fill[mPositions[i].mIndex] = t; |

323 | const aiVector3D& oldpos = mPositions[i].mPosition; |

324 | for (++i; i < fill.size() && mPositions[i].mDistance < maxDist |

325 | && (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i) |

326 | { |

327 | fill[mPositions[i].mIndex] = t; |

328 | } |

329 | ++t; |

330 | } |

331 | |

332 | #ifdef ASSIMP_BUILD_DEBUG |

333 | |

334 | // debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1 |

335 | for (size_t i = 0; i < fill.size(); ++i) { |

336 | ai_assert(fill[i]<mPositions.size()); |

337 | } |

338 | |

339 | #endif |

340 | return t; |

341 | } |

342 |