[2582] | 1 |
|
---|
| 2 | // Written by Vlastimil Havran, December 2007
|
---|
| 3 |
|
---|
| 4 | // This macro allows to use the ray shooting written by
|
---|
| 5 | // Vlastimil Havran, 2007-2008
|
---|
| 6 |
|
---|
| 7 |
|
---|
| 8 | #include "HavranRayCaster.h"
|
---|
| 9 | #include "VssRay.h"
|
---|
| 10 | #include "KdTree.h"
|
---|
| 11 | #include "Preprocessor.h"
|
---|
| 12 | #include "IntersectableWrapper.h"
|
---|
| 13 |
|
---|
| 14 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 15 | #include "timer.h"
|
---|
| 16 | #include "ktball.h"
|
---|
| 17 | #include "ktb.h"
|
---|
| 18 | #include "raypack.h"
|
---|
| 19 | #endif // USE_HAVRAN_RAYCASTER
|
---|
| 20 |
|
---|
| 21 | #define DEBUG_RAYCAST 0
|
---|
| 22 |
|
---|
| 23 | namespace GtpVisibilityPreprocessor {
|
---|
| 24 |
|
---|
| 25 |
|
---|
| 26 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 27 | // static rays
|
---|
| 28 | RayPacket2x2
|
---|
| 29 | HavranRayCaster::raypack;
|
---|
| 30 | #endif // USE_HAVRAN_RAYCASTER
|
---|
| 31 |
|
---|
| 32 | HavranRayCaster::HavranRayCaster(const Preprocessor &preprocessor):
|
---|
| 33 | RayCaster(preprocessor), mKtbtree(0)
|
---|
| 34 | {
|
---|
| 35 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 36 | mKtbtree = new CKTB();
|
---|
| 37 | #endif
|
---|
| 38 | }
|
---|
| 39 |
|
---|
| 40 |
|
---|
| 41 | HavranRayCaster::~HavranRayCaster()
|
---|
| 42 | {
|
---|
| 43 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 44 | delete mKtbtree;
|
---|
| 45 | mKtbtree = 0;
|
---|
| 46 | #endif // USE_HAVRAN_RAYCASTER
|
---|
| 47 | }
|
---|
| 48 |
|
---|
| 49 | void
|
---|
| 50 | HavranRayCaster::Build(ObjectContainer &objlist)
|
---|
| 51 | {
|
---|
| 52 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 53 |
|
---|
| 54 | CTimer timer;
|
---|
| 55 |
|
---|
| 56 | cout << "\nBuilding up kd-tree for Havran ray caster ..."<<endl<<flush;
|
---|
| 57 |
|
---|
| 58 | timer.Start();
|
---|
| 59 | mKtbtree->BuildUp(objlist);
|
---|
| 60 | timer.Stop();
|
---|
| 61 | cout << "\nBuilding up kd-tree is finished, user time = "
|
---|
| 62 | << timer.UserTime() << " real time = " << timer.RealTime() <<
|
---|
| 63 | endl <<flush;
|
---|
| 64 | #endif
|
---|
| 65 | }
|
---|
| 66 |
|
---|
| 67 |
|
---|
| 68 | // Using packet of 4 rays supposing that these are coherent
|
---|
| 69 | // We give a box to which each ray is clipped to before the
|
---|
| 70 | // ray shooting is computed !
|
---|
| 71 | void HavranRayCaster::CastRaysPacket4(Vector3 &boxmax,
|
---|
[2583] | 72 | Vector3 &boxmin,
|
---|
| 73 | Vector3 origin4[],
|
---|
| 74 | Vector3 direction4[],
|
---|
| 75 | int result4[],
|
---|
| 76 | float dist4[])
|
---|
[2582] | 77 | {
|
---|
| 78 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 79 | for (int i = 0; i < 4; i++) {
|
---|
| 80 | result4[i] = -1;
|
---|
| 81 | raypack.SetLoc(i, origin4[i]);
|
---|
| 82 | raypack.SetDir(i, direction4[i]);
|
---|
[2583] | 83 | }
|
---|
| 84 |
|
---|
[2582] | 85 | // The same operations for packets of rays, if implemented by
|
---|
| 86 | // a particular ASDS, otherwise it is emulated by decomposition
|
---|
| 87 | // of a packet to individual rays and traced individually.
|
---|
| 88 | mKtbtree->FindNearestI(raypack, boxmin, boxmax);
|
---|
| 89 |
|
---|
| 90 | for (int i = 0; i < 4; i++) {
|
---|
| 91 | // if ray intersects an object, we set the pointer to
|
---|
| 92 | // this object
|
---|
| 93 | Intersectable* intersectable;
|
---|
| 94 | if ( (intersectable = raypack.GetObject(i)) )
|
---|
[2583] | 95 | // $$JB this is very hacky - > should be replaced with fetching the index of the triangle
|
---|
| 96 | result4[i] = (int)intersectable;
|
---|
[2582] | 97 | }
|
---|
[2583] | 98 |
|
---|
[2582] | 99 | return;
|
---|
| 100 | #endif
|
---|
| 101 | }
|
---|
| 102 |
|
---|
| 103 | int HavranRayCaster::CastRay(const SimpleRay &simpleRay,
|
---|
| 104 | VssRayContainer &vssRays,
|
---|
| 105 | const AxisAlignedBox3 &box,
|
---|
| 106 | const bool castDoubleRay,
|
---|
| 107 | const bool pruneInvalidRays)
|
---|
| 108 | {
|
---|
| 109 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 110 |
|
---|
| 111 | int hits = 0;
|
---|
| 112 | Intersection hitA(simpleRay.mOrigin), hitB(simpleRay.mOrigin);
|
---|
| 113 |
|
---|
| 114 | // inside test for bounding box
|
---|
| 115 | // enlarge box slightly so the view point fits for sure
|
---|
| 116 | // AxisAlignedBox3 sbox = box;
|
---|
| 117 | // sbox.Enlarge(Vector3(-Limits::Small));
|
---|
| 118 | // $$ JB moved here from Validate routine
|
---|
| 119 |
|
---|
| 120 | // if (!box.IsInside(simpleRay.mOrigin)) {
|
---|
| 121 | // cout<<"out of box "<<simpleRay.mOrigin<<" "<<box<<endl;
|
---|
| 122 | // return 0;
|
---|
| 123 | // }
|
---|
| 124 |
|
---|
| 125 | // ray.mFlags &= ~Ray::CULL_BACKFACES;
|
---|
| 126 |
|
---|
| 127 | if (mKtbtree->FindNearestI(simpleRay)) {
|
---|
| 128 | hitA.mObject = SimpleRay::IntersectionRes[0].intersectable;
|
---|
| 129 | float tdist = SimpleRay::IntersectionRes[0].tdist;
|
---|
| 130 | hitA.mPoint = simpleRay.Extrap(tdist);
|
---|
| 131 | hitA.mNormal = SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
|
---|
| 132 | //hitA.mNormal = (dynamic_cast< TriangleIntersectable *>
|
---|
| 133 | // (hitA.mObject))->GetNormal(0);
|
---|
| 134 | }
|
---|
| 135 |
|
---|
| 136 |
|
---|
| 137 | if (castDoubleRay) {
|
---|
| 138 | Vector3 *v = (Vector3*)(&simpleRay.mDirection);
|
---|
| 139 | *v = -(*v);
|
---|
| 140 | // ray.mFlags &= ~Ray::CULL_BACKFACES;
|
---|
| 141 | if (mKtbtree->FindNearestI(simpleRay)) {
|
---|
| 142 | hitB.mObject = SimpleRay::IntersectionRes[0].intersectable;
|
---|
| 143 | float tdist = SimpleRay::IntersectionRes[0].tdist;
|
---|
| 144 | hitB.mPoint = simpleRay.Extrap(tdist);
|
---|
| 145 | hitB.mNormal = SimpleRay::IntersectionRes[0].intersectable->GetNormal(0);
|
---|
| 146 | //hitB.mNormal = (dynamic_cast< TriangleIntersectable *>
|
---|
| 147 | // (hitB.mObject))->GetNormal(0);
|
---|
| 148 | }
|
---|
| 149 | // restore the direction to the original
|
---|
| 150 | *v = -simpleRay.mDirection;
|
---|
| 151 | }
|
---|
| 152 |
|
---|
| 153 | // This code is also in IntelRayCaster.cpp
|
---|
| 154 | return ProcessRay(
|
---|
| 155 | simpleRay,
|
---|
| 156 | hitA,
|
---|
| 157 | hitB,
|
---|
| 158 | vssRays,
|
---|
| 159 | box,
|
---|
| 160 | castDoubleRay,
|
---|
| 161 | pruneInvalidRays
|
---|
| 162 | );
|
---|
| 163 | #else // USE_HAVRAN_RAYCASTER
|
---|
| 164 | return 0;
|
---|
| 165 | #endif
|
---|
| 166 | }
|
---|
| 167 |
|
---|
| 168 | void HavranRayCaster::CastRays16(SimpleRayContainer &rays,
|
---|
| 169 | VssRayContainer &vssRays,
|
---|
| 170 | const AxisAlignedBox3 &sbox,
|
---|
| 171 | const bool castDoubleRay,
|
---|
| 172 | const bool pruneInvalidRays)
|
---|
| 173 | {
|
---|
| 174 | CastRays16(rays, 0, vssRays, sbox, castDoubleRay, pruneInvalidRays);
|
---|
| 175 | }
|
---|
| 176 |
|
---|
| 177 |
|
---|
| 178 | void HavranRayCaster::CastRays16(SimpleRayContainer &rays,
|
---|
| 179 | int offset,
|
---|
| 180 | VssRayContainer &vssRays,
|
---|
| 181 | const AxisAlignedBox3 &sbox,
|
---|
| 182 | const bool castDoubleRays,
|
---|
| 183 | const bool pruneInvalidRays)
|
---|
| 184 | {
|
---|
| 185 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 186 |
|
---|
| 187 | #if DEBUG_RAYCAST
|
---|
| 188 | Debug << "C16 " << flush;
|
---|
| 189 | #endif
|
---|
| 190 |
|
---|
| 191 | #if 0
|
---|
[2588] | 192 | SimpleRayContainer::const_iterator sit = rays.begin() + offset;
|
---|
| 193 | SimpleRayContainer::const_iterator sit_end = sit + 16;
|
---|
| 194 |
|
---|
[2582] | 195 | // no acceleration for ray bundles implemented right now
|
---|
| 196 | for (sit = rays.begin(); sit != sit_end; ++ sit)
|
---|
| 197 | {
|
---|
| 198 | CastRay(*sit, vssRays, sbox, castDoubleRays, pruneInvalidRays);
|
---|
| 199 | }
|
---|
| 200 | #else
|
---|
| 201 | if (castDoubleRays) {
|
---|
| 202 | #if 0
|
---|
| 203 | // Special routine to cast double sided rays
|
---|
| 204 | mKtbtree->SetOffset(0);
|
---|
| 205 | mKtbtree->FindNearestI_16twoDir(rays);
|
---|
| 206 | #else
|
---|
| 207 | // Here we decompose shooting into two phases
|
---|
| 208 | #if 1
|
---|
| 209 | // Here we shoot first forward rays and then backward ones
|
---|
| 210 | mKtbtree->SetOffset(0);
|
---|
| 211 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 212 | SimpleRayContainer::iterator sit = rays.begin() + offset;
|
---|
| 213 | SimpleRayContainer::const_iterator sit_end = rays.begin() + offset + 16;
|
---|
| 214 | for (; sit != sit_end; ++ sit)
|
---|
| 215 | {
|
---|
| 216 | (*sit).mDirection = - (*sit).mDirection;
|
---|
| 217 | }
|
---|
| 218 | // We store the results at different place
|
---|
| 219 | mKtbtree->SetOffset(16);
|
---|
| 220 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 221 | sit = rays.begin() + offset;
|
---|
| 222 | for ( ; sit != sit_end; ++ sit)
|
---|
| 223 | {
|
---|
| 224 | (*sit).mDirection = - (*sit).mDirection;
|
---|
| 225 | }
|
---|
| 226 | #else
|
---|
| 227 | // Here we shoot first backward rays and forward ones
|
---|
| 228 | SimpleRayContainer::iterator sit = rays.begin() + offset;
|
---|
| 229 | SimpleRayContainer::const_iterator sit_end = rays.begin() + offset + 16;
|
---|
| 230 | for ( ; sit != sit_end; ++ sit)
|
---|
| 231 | {
|
---|
| 232 | (*sit).mDirection = - (*sit).mDirection;
|
---|
| 233 | }
|
---|
| 234 | // backward rays to be shot
|
---|
| 235 | mKtbtree->SetOffset(16);
|
---|
| 236 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 237 | sit = rays.begin() + offset;
|
---|
| 238 | for ( ; sit != sit_end; ++ sit)
|
---|
| 239 | {
|
---|
| 240 | (*sit).mDirection = - (*sit).mDirection;
|
---|
| 241 | }
|
---|
| 242 | // forward rays to be shot
|
---|
| 243 | mKtbtree->SetOffset(0);
|
---|
| 244 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 245 | #endif
|
---|
| 246 | #endif
|
---|
| 247 | }
|
---|
| 248 | else {
|
---|
| 249 | // Shoot all 16 rays at the same time using a special algorithm
|
---|
| 250 | mKtbtree->SetOffset(0);
|
---|
| 251 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 252 | }
|
---|
| 253 | #endif
|
---|
| 254 |
|
---|
| 255 |
|
---|
| 256 | for (int i=0, k=offset; i < 16; i++, k++)
|
---|
| 257 | {
|
---|
| 258 | Intersection hitA(rays[k].mOrigin), hitB(rays[k].mOrigin);
|
---|
| 259 |
|
---|
| 260 | #if DEBUG_RAYCAST
|
---|
| 261 | Debug<<"FH\n"<<flush;
|
---|
| 262 | #endif
|
---|
| 263 |
|
---|
| 264 | Intersectable *intersect = SimpleRay::IntersectionRes[i].intersectable;
|
---|
| 265 |
|
---|
| 266 | if (intersect)
|
---|
| 267 | {
|
---|
| 268 | hitA.mObject = intersect;
|
---|
| 269 | // Get the normal of that face
|
---|
| 270 | hitA.mNormal = intersect->GetNormal(0);
|
---|
| 271 |
|
---|
| 272 | //-rays[index+i].mDirection; // $$ temporary
|
---|
| 273 | float tdist = SimpleRay::IntersectionRes[i].tdist;
|
---|
| 274 | hitA.mPoint = rays[k].Extrap(tdist);
|
---|
| 275 | }
|
---|
| 276 |
|
---|
| 277 | #if DEBUG_RAYCAST
|
---|
| 278 | Debug<<"BH\n"<<flush;
|
---|
| 279 | #endif
|
---|
| 280 |
|
---|
| 281 | if (castDoubleRays)
|
---|
| 282 | {
|
---|
| 283 | Intersectable *intersect =
|
---|
[2585] | 284 | SimpleRay::IntersectionRes[i+16].intersectable;
|
---|
[2582] | 285 |
|
---|
| 286 | if (intersect)
|
---|
| 287 | {
|
---|
[2585] | 288 | hitB.mObject = intersect;
|
---|
| 289 | hitB.mNormal = intersect->GetNormal(0);
|
---|
| 290 | // $$ JB : i -> i+16
|
---|
| 291 | float tdist = SimpleRay::IntersectionRes[i+16].tdist;
|
---|
| 292 | hitB.mPoint = rays[k].Extrap(-tdist);
|
---|
[2582] | 293 | }
|
---|
| 294 | }
|
---|
| 295 |
|
---|
| 296 | #if DEBUG_RAYCAST
|
---|
| 297 | Debug<<"PR\n"<<flush;
|
---|
| 298 | #endif
|
---|
| 299 |
|
---|
| 300 | #if 1
|
---|
| 301 | ProcessRay(rays[k],
|
---|
[2585] | 302 | hitA,
|
---|
| 303 | hitB,
|
---|
| 304 | vssRays,
|
---|
| 305 | sbox,
|
---|
| 306 | castDoubleRays,
|
---|
| 307 | pruneInvalidRays
|
---|
| 308 | );
|
---|
[2582] | 309 | #endif
|
---|
| 310 | } // for
|
---|
| 311 |
|
---|
| 312 |
|
---|
| 313 | #if DEBUG_RAYCAST
|
---|
| 314 | Debug<<"C16F\n"<<flush;
|
---|
| 315 | #endif
|
---|
| 316 |
|
---|
| 317 | #endif // USE_HAVRAN_RAYCASTER
|
---|
| 318 | }
|
---|
| 319 |
|
---|
| 320 |
|
---|
| 321 |
|
---|
| 322 | void
|
---|
| 323 | HavranRayCaster::CastSimpleForwardRays(
|
---|
| 324 | SimpleRayContainer &rays,
|
---|
| 325 | const AxisAlignedBox3 &sbox
|
---|
| 326 | )
|
---|
| 327 | {
|
---|
| 328 | int hit_triangles[16];
|
---|
| 329 | float dist[16];
|
---|
| 330 | Vector3 normals[16];
|
---|
| 331 | Vector3 min = sbox.Min();
|
---|
| 332 | Vector3 max = sbox.Max();
|
---|
| 333 |
|
---|
| 334 | int packets = rays.size() / 16;
|
---|
| 335 |
|
---|
| 336 | int i, j, k = 0;
|
---|
| 337 | Vector3 dir;
|
---|
| 338 |
|
---|
| 339 | // By groups of rays
|
---|
| 340 | for (i=0; i < packets; i++) {
|
---|
| 341 | int offset = i * 16;
|
---|
| 342 | mKtbtree->FindNearestI_16oneDir(rays, offset);
|
---|
| 343 | // ??? What to do with the results ? These are
|
---|
| 344 | // not used at the moment also in IntelRayCaster.cpp
|
---|
| 345 | } // for
|
---|
| 346 |
|
---|
| 347 |
|
---|
| 348 | for (; k < rays.size(); k++) {
|
---|
| 349 | double normal[3];
|
---|
| 350 | mKtbtree->FindNearestI(rays[k]);
|
---|
| 351 | // ??? What to do with the results ? These are
|
---|
| 352 | // not used at the moment also in IntelRayCaster.cpp
|
---|
| 353 | }
|
---|
| 354 |
|
---|
| 355 |
|
---|
| 356 | return;
|
---|
| 357 | }
|
---|
| 358 |
|
---|
| 359 | void HavranRayCaster::CastRays(
|
---|
| 360 | SimpleRayContainer &rays,
|
---|
| 361 | VssRayContainer &vssRays,
|
---|
| 362 | const AxisAlignedBox3 &sbox,
|
---|
| 363 | const bool castDoubleRay,
|
---|
| 364 | const bool pruneInvalidRays )
|
---|
| 365 | {
|
---|
| 366 | int buckets = rays.size()/16;
|
---|
| 367 | int offset = 0;
|
---|
| 368 |
|
---|
| 369 | #if 0
|
---|
| 370 | int time = GetTime();
|
---|
| 371 | CastSimpleForwardRays(rays, sbox);
|
---|
| 372 | cout<<1e-3*2*rays.size()/TimeDiff(time, GetTime())<<" Mrays/sec"<<endl;
|
---|
| 373 | #endif
|
---|
| 374 |
|
---|
| 375 | for (int i=0; i < buckets; i++, offset+=16) {
|
---|
| 376 | CastRays16(rays, offset, vssRays, sbox,
|
---|
| 377 | castDoubleRay, pruneInvalidRays);
|
---|
| 378 |
|
---|
| 379 | if ((int)rays.size() > 100000 && i % (100000/16) == 0)
|
---|
| 380 | cout<<"\r"<<offset<<"/"<<(int)rays.size()<<"\r";
|
---|
| 381 | }
|
---|
| 382 |
|
---|
| 383 | for (; offset < (int)rays.size(); offset++)
|
---|
| 384 | CastRay(rays[offset], vssRays, sbox, castDoubleRay, pruneInvalidRays);
|
---|
| 385 | }
|
---|
| 386 |
|
---|
| 387 | #if 0
|
---|
| 388 | // BUG1 41579 196948 1064111
|
---|
| 389 | // BUG2 254 1672 10869
|
---|
| 390 |
|
---|
| 391 | // Just for testing concept
|
---|
| 392 | void
|
---|
| 393 | HavranRayCaster::CastRaysPacket2x2(RayPacket2x2 &raysPack,
|
---|
| 394 | bool castDoubleRay,
|
---|
| 395 | const bool pruneInvalidRays)
|
---|
| 396 | {
|
---|
| 397 | #ifdef USE_HAVRAN_RAYCASTER
|
---|
| 398 |
|
---|
| 399 | if (castDoubleRay) {
|
---|
| 400 | // cast forward rays
|
---|
| 401 | mKtbtree->FindNearestI(raysPack);
|
---|
| 402 | for (int i = 0; i < 4; i++)
|
---|
| 403 | raysPack.SetDir(i, -raysPack.GetDir(i));
|
---|
| 404 | // cast backward rays
|
---|
| 405 | mKtbtree->FindNearestI(raysPack);
|
---|
| 406 | // reverse the rays back
|
---|
| 407 | for (int i = 0; i < 4; i++)
|
---|
| 408 | raysPack.SetDir(i, -raysPack.GetDir(i));
|
---|
| 409 | }
|
---|
| 410 | else {
|
---|
| 411 | // only rays forwards
|
---|
| 412 | mKtbtree->FindNearestI(raysPack);
|
---|
| 413 | #if 0
|
---|
| 414 | // Only verification of correctness by casting single rays
|
---|
| 415 | static int cntBugs = 0;
|
---|
| 416 | SimpleRay ray;
|
---|
| 417 | int cntErrors = 0, cntDistErrors = 0;
|
---|
| 418 | bool newBug = false;
|
---|
| 419 | for (int i = 0; i < 4; i++) {
|
---|
| 420 | ray.mOrigin = raysPack.GetLoc(i);
|
---|
| 421 | ray.mDirection = raysPack.GetDir(i);
|
---|
| 422 | mKtbtree->FindNearestI(ray);
|
---|
| 423 | if (raysPack.GetObject(i) != SimpleRay::IntersectionRes[0].intersectable) {
|
---|
| 424 | float dist = (raysPack.GetT(i) - SimpleRay::IntersectionRes[0].tdist);
|
---|
| 425 | if (fabs(dist) > 0.001f) {
|
---|
| 426 | cntErrors++; newBug = true;
|
---|
| 427 | cntBugs++;
|
---|
| 428 | cout << " BUG1 d= " << dist;
|
---|
| 429 | }
|
---|
| 430 | }
|
---|
| 431 | else {
|
---|
| 432 | float dist = 0.f;
|
---|
| 433 | if (raysPack.GetObject(i) && SimpleRay::IntersectionRes[0].intersectable)
|
---|
| 434 | if (fabs((dist=(fabs (raysPack.GetT(i) - SimpleRay::IntersectionRes[0].tdist)))) > 1.f) {
|
---|
| 435 | cntDistErrors++; newBug = true; cntBugs++;
|
---|
| 436 | cout << " BUG2 distdf= " << dist ;
|
---|
| 437 | }
|
---|
| 438 | }
|
---|
| 439 | } // for
|
---|
| 440 | if (newBug) cout << " CB= " << cntBugs << "\n";
|
---|
| 441 | #endif
|
---|
| 442 | }
|
---|
| 443 |
|
---|
| 444 | return;
|
---|
| 445 | #endif // USE_HAVRAN_RAYCASTER
|
---|
| 446 | }
|
---|
| 447 | #endif
|
---|
| 448 |
|
---|
| 449 |
|
---|
| 450 | } // the namespace
|
---|