1#ifndef RNAPUZZLER_INTERSECTLEVELBOUNDINGBOXES_H 2#define RNAPUZZLER_INTERSECTLEVELBOUNDINGBOXES_H 3 4/* 5 * RNApuzzler intersect bounding boxes 6 * 7 * c Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer 8 * ViennaRNA package 9 */ 10 11#include <stdlib.h> 12#include <math.h> 13 14#include "intersectLevelLines.inc" 15#include "boundingBoxes.inc" 16#include "configtree.inc" 17#include "vector_math.inc" 18#include "definitions.inc" 19 20/** 21 * @brief intersectStemStem 22 * Checks for intersection of two given stems. 23 * @param stem1 24 * @param stem2 25 * @return 1 if intersecting, 0 otherwise 26 */ 27PRIVATE short 28intersectStemStem(const stemBox *stem1, 29 const stemBox *stem2); 30 31 32/** 33 * @brief intersectLoopLoop 34 * Checks for intersection of two given loops. 35 * @param loop1 36 * @param loop2 37 * @return 1 if intersecting, 0 otherwise 38 */ 39PRIVATE short 40intersectLoopLoop(const loopBox *loop1, 41 const loopBox *loop2); 42 43 44/** 45 * @brief intersectStemLoop 46 * Check for intersection of given stem and given loop. 47 * @param stem 48 * @param loop 49 * @return 1 if intersecting, 0 otherwise 50 */ 51PRIVATE short 52intersectStemLoop(const stemBox *stem, 53 const loopBox *loop); 54 55 56/** 57 * @brief intersectLoopBulges 58 * Check for intersection of given loop and given stem's bulges. 59 * @param loop 60 * @param stem 61 * @param bulge 62 * @return 1 if intersecting, 0 otherwise 63 */ 64PRIVATE short 65intersectLoopBulges(const loopBox *loop, 66 const stemBox *stem, 67 int *bulge); 68 69 70/** 71 * @brief intersectBulgesBulges 72 * @return 73 */ 74PRIVATE short 75intersectBulgesBulges(const stemBox *stem1, 76 const stemBox *stem2, 77 int *bulge1, 78 int *bulge2); 79 80 81/** 82 * @brief intersectStemBulges 83 * @param stem1 84 * @param stem2 85 * @param bulge2 86 * @return 87 */ 88PRIVATE short 89intersectStemBulges(const stemBox *stem1, 90 const stemBox *stem2, 91 int *bulge2); 92 93 94PRIVATE short 95intersectLoopLoop(const loopBox *loop1, 96 const loopBox *loop2) 97{ 98 double c1[2] = { 99 loop1->c[0], loop1->c[1] 100 }; 101 double r1 = loop1->r + 0.5 * epsilonRecognize; 102 double c2[2] = { 103 loop2->c[0], loop2->c[1] 104 }; 105 double r2 = loop2->r + 0.5 * epsilonRecognize; 106 107 return intersectCircleCircle(c1, r1, c2, r2); 108} 109 110 111PRIVATE void 112projectPointOntoLine(const double a[2], 113 const double b[2], 114 const double p[2], 115 double ret_p[2]) 116{ 117 double u[2]; 118 119 vector(a, p, u); 120 double v[2]; 121 vector(a, b, v); 122 double w[2] = { 123 -v[1], v[0] 124 }; 125 126 /* compute r which grants all needed information */ 127 double r = (u[1] - u[0] * w[1] / w[0]) / (v[1] - v[0] * w[1] / w[0]); 128 129 if (r < 0.0) { 130 ret_p[0] = a[0]; 131 ret_p[1] = a[1]; 132 } else if (r > 1.0) { 133 ret_p[0] = b[0]; 134 ret_p[1] = b[1]; 135 } else { 136 ret_p[0] = a[0] + r * v[0]; 137 ret_p[1] = a[1] + r * v[1]; 138 } 139 140 return; 141} 142 143 144PRIVATE void 145ClosestPtPointBulge(const double p[2], 146 const double a[2], 147 const double b[2], 148 const double c[2], 149 double ret_p[2]) 150{ 151 char *fnName = "CLOSEST ON BULGE"; 152 153 /* 154 * Note: 155 * 156 * In contrast to ClosestPtPointTriangle (taken from book Real-Time Collision Detection) 157 * this function does not work for general triangles. 158 * We implicitely make use of the fact that our bulges are equiliteral triangles 159 * or to be more precise there are no angles greater than 90° in our triangles. 160 * This allows for only checking for the first occurance of a side where point p 161 * is on the outer side of the triangle and not checking any further. 162 * 163 * In fact applying this function to a triangle with some angle being greater than 90° 164 * may lead to wrong results as follows as described in the mentioned book 165 * (at ClosestPtPointTriangle) as well. 166 */ 167 168 169 /* check if p on outer side of AB */ 170 short orientABC = isToTheRightPointPoint(a, b, c); 171 short orientABP = isToTheRightPointPoint(a, b, p); 172 173 if (orientABC != orientABP) { 174 /* p is on outer side of AB */ 175 projectPointOntoLine(a, b, p, ret_p); 176 return; 177 } 178 179 /* check if p on outer side of BC */ 180 short orientBCA = isToTheRightPointPoint(b, c, a); 181 short orientBCP = isToTheRightPointPoint(b, c, p); 182 183 if (orientBCA != orientBCP) { 184 /* p is on outer side of BC */ 185 projectPointOntoLine(b, c, p, ret_p); 186 return; 187 } 188 189 /* check if p on outer side of CA */ 190 short orientCAB = isToTheRightPointPoint(c, a, b); 191 short orientCAP = isToTheRightPointPoint(c, a, p); 192 193 if (orientCAB != orientCAP) { 194 /* p is on outer side of CA */ 195 projectPointOntoLine(c, a, p, ret_p); 196 return; 197 } 198 199 /* p is inside ABC */ 200 ret_p[0] = p[0]; 201 ret_p[1] = p[1]; 202 return; 203} 204 205 206PRIVATE void 207ClosestPtPointTriangle(const double p[2], 208 const double a[2], 209 const double b[2], 210 const double c[2], 211 double ret_p[2]) 212{ 213 /* Check if P in vertex region outside A */ 214 double ab[2]; 215 double ac[2]; 216 double ap[2]; 217 218 vector(a, b, ab); 219 vector(a, c, ac); 220 vector(a, p, ap); 221 double d1 = scalarProduct2D(ab, ap); 222 double d2 = scalarProduct2D(ac, ap); 223 224 if (d1 <= 0.0 && d2 <= 0.0) { 225 /* barycentric coordinates (1,0,0) */ 226 ret_p[0] = a[0]; 227 ret_p[1] = a[1]; 228 return; 229 } 230 231 /* Check if P in vertex region outside B */ 232 double bp[2]; 233 vector(b, p, bp); 234 double d3 = scalarProduct2D(ab, bp); 235 double d4 = scalarProduct2D(ac, bp); 236 237 if (d3 >= 0.0 && d4 <= 0.0) { 238 /* barycentric coordinates (0,1,0) */ 239 ret_p[0] = b[0]; 240 ret_p[1] = b[1]; 241 return; 242 } 243 244 /* Check if P in edge region of AB, if so return projection of P onto AB */ 245 double vc = d1 * d4 - d3 * d2; 246 247 if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) { 248 /* barycentric coordinates (1-v,v,0) */ 249 double v = d1 / (d1 - d3); 250 ret_p[0] = a[0] + v * ab[0]; 251 ret_p[1] = a[1] + v * ab[1]; 252 return; 253 } 254 255 /* Check if P in vertex region outside C */ 256 double cp[2]; 257 vector(c, p, cp); 258 double d5 = scalarProduct2D(ab, cp); 259 double d6 = scalarProduct2D(ac, cp); 260 261 if (d6 >= 0.0 && d5 <= d6) { 262 /* barycentric coordinates (0,0,1) */ 263 ret_p[0] = c[0]; 264 ret_p[1] = c[1]; 265 return; 266 } 267 268 /* Check if P in edge region of AC, if so return projection of P onto AC */ 269 double vb = d5 * d2 - d1 * d6; 270 if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) { 271 /* barycentric coordinates (1-w,0,w) */ 272 double w = d2 / (d2 - d6); 273 ret_p[0] = a[0] + w * ac[0]; 274 ret_p[1] = a[1] + w * ac[1]; 275 return; 276 } 277 278 /* Check if P in edge region of BC, if so return projection of P onto BC */ 279 double va = d3 * d6 - d5 * d4; 280 281 if (va <= 0.0) { 282 if ((d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) { 283 /* barycentric coordinates (0,1-w,w) */ 284 double w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); 285 ret_p[0] = b[0] + w * (c[0] - b[0]); 286 ret_p[1] = b[1] + w * (c[1] - b[1]); 287 return; 288 } else if ((d4 - d3) <= 0.0 && (d5 - d6) >= 0.0) { 289 /* barycentric coordinates (0,1-w,w) */ 290 double w = (d3 - d4) / ((d3 - d4) + (d5 - d6)); 291 ret_p[0] = b[0] + w * (c[0] - b[0]); 292 ret_p[1] = b[1] + w * (c[1] - b[1]); 293 return; 294 } 295 } 296 297 /* P inside face region. Compute Q through its barycentric coordinates (u,v,w) */ 298 double denom = 1.0 / (va + vb + vc); 299 double v = vb * denom; 300 double w = vc * denom; 301 ret_p[0] = a[0] + ab[0] * v + ac[0] * w; 302 ret_p[1] = a[1] + ab[1] * v + ac[1] * w; 303 return; 304} 305 306 307/** 308 * @brief ClosestPtPointOBB 309 * Implementation of ClostestPtPointOBB from the book "Real Time Collision Detection". 310 * Calculates the point on a rectangle (stem) that is closest to a given point p. 311 * @param stem - rectangle / stem bounding box 312 * @param p_x - x value of point p 313 * @param p_y - y value of point p 314 * @param ret_p_x - pointer to closest point's x value (used as return) 315 * @param ret_p_y - pointer to closest point's y value (used as return) 316 */ 317PRIVATE void 318ClosestPtPointOBB(const stemBox stem, 319 const double p[2], 320 double ret_p[2]) 321{ 322 double u0[2] = { 323 stem.a[0], stem.a[1] 324 }; 325 double u1[2] = { 326 stem.b[0], stem.b[1] 327 }; 328 double dv[2]; 329 330 vector(stem.c, p, dv); 331 332 double dist_0 = scalarProduct2D(dv, u0); 333 double dist_1 = scalarProduct2D(dv, u1); 334 335 short sign_d0 = dist_0 < 0 ? -1 : 1; 336 short sign_d1 = dist_1 < 0 ? -1 : 1; 337 short sign_e0 = stem.e[0] < 0 ? -1 : 1; 338 short sign_e1 = stem.e[1] < 0 ? -1 : 1; 339 double abs_d0 = sign_d0 * dist_0; 340 double abs_d1 = sign_d1 * dist_1; 341 double abs_e0 = sign_e0 * stem.e[0]; 342 double abs_e1 = sign_e1 * stem.e[1]; 343 344 /* clamp dist_0/1 to the extents of the OBB */ 345 double clamped_0 = abs_d0 > abs_e0 ? sign_d0 * abs_e0 : sign_d0 * abs_d0; 346 double clamped_1 = abs_d1 > abs_e1 ? sign_d1 * abs_e1 : sign_d1 * abs_d1; 347 348 ret_p[0] = stem.c[0] + clamped_0 * stem.a[0] + clamped_1 * stem.b[0]; 349 ret_p[1] = stem.c[1] + clamped_0 * stem.a[1] + clamped_1 * stem.b[1]; 350} 351 352 353PRIVATE short 354intersectStemLoop(const stemBox *stem, 355 const loopBox *loop) 356{ 357 short intersect = 0; 358 359 double p[2]; 360 361 ClosestPtPointOBB(*stem, loop->c, p); 362 363 double v_c_to_p[2]; 364 vector(loop->c, p, v_c_to_p); 365 366 double distanceSquared = vectorLength2DSquared(v_c_to_p); 367 intersect = (distanceSquared < 368 ((loop->r + epsilonRecognize) * 369 (loop->r + epsilonRecognize))); 370 /* add epsilon to classify nearby objects as intersecting */ 371 return intersect; 372} 373 374 375PRIVATE short 376intersectStemStem(const stemBox *stem1, 377 const stemBox *stem2) 378{ 379 /* brute force approach for intersecting two rectangles */ 380 double stem1_ea[2] = { 381 stem1->e[0] * stem1->a[0], 382 stem1->e[0] * stem1->a[1] 383 }; 384 double stem1_eb[2] = { 385 stem1->e[1] * stem1->b[0], 386 stem1->e[1] * stem1->b[1] 387 }; 388 double B1[2] = { 389 stem1->c[0] + stem1_ea[0] + stem1_eb[0], 390 stem1->c[1] + stem1_ea[1] + stem1_eb[1] 391 }; 392 double C1[2] = { 393 stem1->c[0] + stem1_ea[0] - stem1_eb[0], 394 stem1->c[1] + stem1_ea[1] - stem1_eb[1] 395 }; 396 double D1[2] = { 397 stem1->c[0] - stem1_ea[0] - stem1_eb[0], 398 stem1->c[1] - stem1_ea[1] - stem1_eb[1] 399 }; 400 double A1[2] = { 401 stem1->c[0] - stem1_ea[0] + stem1_eb[0], 402 stem1->c[1] - stem1_ea[1] + stem1_eb[1] 403 }; 404 405 double stem2_ea[2] = { 406 stem2->e[0] * stem2->a[0], 407 stem2->e[0] * stem2->a[1] 408 }; 409 double stem2_eb[2] = { 410 stem2->e[1] * stem2->b[0], 411 stem2->e[1] * stem2->b[1] 412 }; 413 double B2[2] = { 414 stem2->c[0] + stem2_ea[0] + stem2_eb[0], 415 stem2->c[1] + stem2_ea[1] + stem2_eb[1] 416 }; 417 double C2[2] = { 418 stem2->c[0] + stem2_ea[0] - stem2_eb[0], 419 stem2->c[1] + stem2_ea[1] - stem2_eb[1] 420 }; 421 double D2[2] = { 422 stem2->c[0] - stem2_ea[0] - stem2_eb[0], 423 stem2->c[1] - stem2_ea[1] - stem2_eb[1] 424 }; 425 double A2[2] = { 426 stem2->c[0] - stem2_ea[0] + stem2_eb[0], 427 stem2->c[1] - stem2_ea[1] + stem2_eb[1] 428 }; 429 430 /* 431 * Only the sides of the stems (AB, CD) need to be intersected against 432 * each other 433 */ 434 if (intersectLineSegments(A1, B1, A2, B2, NULL) 435 || intersectLineSegments(A1, B1, C2, D2, NULL) 436 || intersectLineSegments(C1, D1, A2, B2, NULL) 437 || intersectLineSegments(C1, D1, C2, D2, NULL) 438 ) 439 return 1; 440 else 441 return 0; 442} 443 444 445/* Returns true if circle circ intersects triangle ABC, false otherwise. */ 446PRIVATE short 447TestCircleTriangle(const double circ_c[2], 448 const double circ_r, 449 const double A[2], 450 const double B[2], 451 const double C[2], 452 double p[2]) 453{ 454 /* Find point P on triangle ABC closest to circle center */ 455 ClosestPtPointBulge(circ_c, A, B, C, p); 456 457 /* 458 * circle and triangle intersect if the (squared) distance from circle 459 * center to point is less than the (squared) circle radius 460 */ 461 double v[2]; 462 vector(p, circ_c, v); 463 464 short ret = scalarProduct2D(v, v) <= (circ_r * circ_r); 465 466 return ret; 467} 468 469 470PRIVATE short 471TestLoopBulge(const double c[2], 472 const double r, 473 const double pPrev[2], 474 const double pThis[2], 475 const double pNext[2]) 476{ 477 /* check if bulge point is inside loop */ 478 double vCenterBulge[2]; 479 480 vector(c, pThis, vCenterBulge); 481 /* compare r and length of vCenterBulge; using squared distances is less expensive that sqrt */ 482 if (r * r > vectorLength2DSquared(vCenterBulge)) 483 484 return 1; 485 486 double vPrevThis[2]; 487 vector(pPrev, pThis, vPrevThis); 488 double vThisNext[2]; 489 vector(pThis, pNext, vThisNext); 490 491 double cut1[2], cut2[2]; 492 short numCut; 493 494 /* evaluate line pPrev->pThis */ 495 numCut = getCutPointsOfCircleAndLine(c, r, pPrev, vPrevThis, cut1, cut2); 496 497 if (numCut > 0) { 498 if (matchLinePoint(pPrev, vPrevThis, cut1)) 499 500 return 1; 501 } 502 503 if (numCut > 1) 504 if (matchLinePoint(pPrev, vPrevThis, cut2)) 505 return 1; 506 507 /* evaluate line pThis->pNext */ 508 numCut = getCutPointsOfCircleAndLine(c, r, pThis, vThisNext, cut1, cut2); 509 if (numCut > 0) 510 if (matchLinePoint(pThis, vThisNext, cut1)) 511 return 1; 512 513 if (numCut > 1) { 514 if (matchLinePoint(pThis, vThisNext, cut2)) 515 516 return 1; 517 } 518 519 return 0; 520} 521 522 523PRIVATE short 524intersectLoopBulges(const loopBox *loop, 525 const stemBox *stem, 526 int *bulge) 527{ 528 *bulge = -1; 529 530 double c[2] = { 531 loop->c[0], loop->c[1] 532 }; 533 double r = loop->r + epsilonRecognize; 534 535 for (int currentBulge = 0; currentBulge < stem->bulgeCount; currentBulge++) { 536 double A[2], B[2], C[2]; 537 getBulgeCoordinates(stem, currentBulge, A, B, C); 538 539 double p[2]; 540 541 if (TestCircleTriangle(c, r, A, B, C, p)) { 542 *bulge = currentBulge; 543 return 1; 544 } 545 } 546 547 return 0; 548} 549 550 551PRIVATE short 552intersectBulgesBulges(const stemBox *stem1, 553 const stemBox *stem2, 554 int *bulge1, 555 int *bulge2) 556{ 557 *bulge1 = -1; 558 *bulge2 = -1; 559 560 double distance = 0.5 * epsilonRecognize; 561 562 for (int currentBulge1 = 0; currentBulge1 < stem1->bulgeCount; currentBulge1++) { 563 double piPrev[2], piThis[2], piNext[2]; 564 getBulgeCoordinatesExtraDistance(stem1, currentBulge1, distance, piPrev, piThis, piNext); 565 566 for (int currentBulge2 = 0; currentBulge2 < stem2->bulgeCount; currentBulge2++) { 567 double pjPrev[2], pjThis[2], pjNext[2]; 568 getBulgeCoordinatesExtraDistance(stem2, currentBulge2, distance, pjPrev, pjThis, pjNext); 569 570 if (intersectLineSegments(piPrev, piThis, pjPrev, pjThis, NULL) 571 || intersectLineSegments(piPrev, piThis, pjThis, pjNext, NULL) 572 || intersectLineSegments(piThis, piNext, pjPrev, pjThis, NULL) 573 || intersectLineSegments(piThis, piNext, pjThis, pjNext, NULL) 574 ) { 575 *bulge1 = currentBulge1; 576 *bulge2 = currentBulge2; 577 return 1; 578 } 579 } 580 } 581 582 return 0; 583} 584 585 586PRIVATE short 587intersectStemBulges(const stemBox *stem1, 588 const stemBox *stem2, 589 int *bulge2) 590{ 591 *bulge2 = -1; 592 593 if (stem2->bulgeCount == 0) 594 595 return 0; 596 597 /* 598 * simplify to only check bulge lines against left and right stem lines 599 * 600 * if the bulge is surrounded by the stem then there is a Stem vs. Stem intersection 601 * if the bulge intersects the stem's bottom or top line then there is an intersection with the adjacent loop 602 * -> no need to checks those cases 603 */ 604 605 /* 606 * N - North, E - East, S - South, W - West 607 * north is direction to loop 608 */ 609 double pNW[2]; 610 pNW[0] = stem1->c[0] + stem1->e[0] * stem1->a[0] - stem1->e[1] * stem1->b[0]; 611 pNW[1] = stem1->c[1] + stem1->e[0] * stem1->a[1] - stem1->e[1] * stem1->b[1]; 612 613 double pSW[2]; 614 pSW[0] = stem1->c[0] - stem1->e[0] * stem1->a[0] - stem1->e[1] * stem1->b[0]; 615 pSW[1] = stem1->c[1] - stem1->e[0] * stem1->a[1] - stem1->e[1] * stem1->b[1]; 616 617 double pNE[2]; 618 pNE[0] = stem1->c[0] + stem1->e[0] * stem1->a[0] + stem1->e[1] * stem1->b[0]; 619 pNE[1] = stem1->c[1] + stem1->e[0] * stem1->a[1] + stem1->e[1] * stem1->b[1]; 620 621 double pSE[2]; 622 pSE[0] = stem1->c[0] - stem1->e[0] * stem1->a[0] + stem1->e[1] * stem1->b[0]; 623 pSE[1] = stem1->c[1] - stem1->e[0] * stem1->a[1] + stem1->e[1] * stem1->b[1]; 624 625 double distance = epsilonRecognize; 626 627 for (int currentBulge2 = 0; currentBulge2 < stem2->bulgeCount; currentBulge2++) { 628 double pPrev[2], pThis[2], pNext[2]; 629 getBulgeCoordinatesExtraDistance(stem2, currentBulge2, distance, pPrev, pThis, pNext); 630 631 if (intersectLineSegments(pNW, pSW, pPrev, pThis, NULL) 632 || intersectLineSegments(pNW, pSW, pThis, pNext, NULL) 633 || intersectLineSegments(pNE, pSE, pPrev, pThis, NULL) 634 || intersectLineSegments(pNE, pSE, pThis, pNext, NULL) 635 ) { 636 *bulge2 = currentBulge2; 637 return 1; 638 } 639 } 640 641 return 0; 642} 643 644 645#endif 646