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