1#ifndef RNAPUZZLER_INTERSECTLEVELLINES_H
2#define RNAPUZZLER_INTERSECTLEVELLINES_H
3
4/*
5 *      RNApuzzler intersect lines
6 *
7 *      c  Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer
8 *      ViennaRNA package
9 */
10#include <math.h>
11#include <stdlib.h>
12#include <stdio.h>
13
14#include "ViennaRNA/utils/basic.h"
15
16#include "definitions.inc"
17#include "vector_math.inc"
18
19/**
20 * @brief intersectLineArc
21 *      checks for intersection of a circle arc and straight line
22 * @param point_1 first point of line
23 * @param point_2 second point of line
24 * @param arc arc
25 * @return 1 if intersecting, 0 otherwise
26 */
27PRIVATE short
28intersectLineArc(const double point_1[2],
29                 const double point_2[2],
30                 const double arc[6]);
31
32
33PRIVATE short
34intersectArcArc(const double  arc1[6],
35                const double  arc2[6]);
36
37
38PRIVATE short
39intersectCircleCircle(const double  c1[2],
40                      const double  c1r,
41                      const double  c2[2],
42                      const double  c2r);
43
44
45PRIVATE short
46matchLinePoint(const double pLine[2],
47               const double dirLine[2],
48               const double p[2]);
49
50
51/**
52 * @brief intersectLineSegments
53 *      Determines if two line segments defined by the points A,B and X,Y each intersect.
54 *      If they do while both lines are not parallel the coordinates of their cut point are returned in P.
55 * @param A start point of first line segment
56 * @param B end point of first line segment
57 * @param X start point of second line segment
58 * @param Y end point of second line segment
59 * @param P return value for cut point of both line segments
60 * @return 1 if intersecting , 0 otherwise
61 */
62PRIVATE short
63intersectLineSegments(const double  A[2],
64                      const double  B[2],
65                      const double  X[2],
66                      const double  Y[2],
67                      double        P[2]);
68
69
70/**
71 * @brief matchPointArc
72 *      checks if the given point is situated on the given wedge
73 * @param point double[2] array with point coordinates
74 * @param arc double[6] array with arc coordinates [circleX, circleY, circleR, arcFrom, arcTo, clockwise]
75 * @return 1 if point is situated on the wedge, 0 otherwise
76 */
77PRIVATE short
78matchPointArc(const double  point[2],
79              const double  arc[6])
80{
81  double  center[2] = {
82    arc[0], arc[1]
83  };
84  double  from      = toRad(arc[3]);
85  double  to        = toRad(arc[4]);
86  short   clockwise = (arc[5] > 0.5);
87
88  double  v_center_point[2];
89
90  vector(center, point, v_center_point);
91
92  double  zero_degree[2] = {
93    1, 0
94  };
95  double  angle = angleBetweenVectors2D(v_center_point, zero_degree);
96
97  if (point[1] < center[1])
98
99    angle = MATH_TWO_PI - angle;
100
101  short ret = 0;
102  if (clockwise) {
103    if (from > to)
104      /*
105       * normal case
106       * check: from < angle < to
107       */
108      ret = (from >= angle && angle >= to);
109    else
110      /*
111       * interval surpasses 360° border
112       * check: from < angle < 360 or 0 < angle < to
113       */
114      ret = (from >= angle && angle >= 0) || (MATH_TWO_PI >= angle && angle >= to);
115  } else {
116    if (from < to)
117      /*
118       * normal case
119       * check: to < angle < from
120       */
121      ret = (from <= angle && angle <= to);
122    else
123      /*
124       * interval surpasses 360° border
125       * check: from > angle > 0 || MATH_TWO_PI > angle > to
126       */
127      ret = (from <= angle && angle <= MATH_TWO_PI) || (0 <= angle && angle <= to);
128  }
129
130  return ret;
131}
132
133
134PRIVATE short
135matchLinePoint(const double pLine[2],
136               const double dirLine[2],
137               const double p[2])
138{
139  double t = -1.0;                /* init as not matching */
140
141  if (fabs(dirLine[0]) > 0.0001)  /* != 0.0 */
142
143    t = (p[0] - pLine[0]) / dirLine[0];
144
145  else if (fabs(dirLine[1]) > 0.0001)        /* != 0.0 */
146
147    t = (p[1] - pLine[1]) / dirLine[1];
148
149  else
150
151    return 0;     /* this is not even a line since dirLine == (0.0, 0.0) */
152
153  return 0.0 <= t && t <= 1.0;
154}
155
156
157PRIVATE short
158intersectCircleCircle(const double  c1[2],
159                      const double  c1r,
160                      const double  c2[2],
161                      const double  c2r)
162{
163  double  v_c1_c2[2];
164
165  vector(c1, c2, v_c1_c2);
166  double  distance = vectorLength2D(v_c1_c2);
167
168  short   intersect = (distance < (c1r + c2r));
169
170  return intersect;
171}
172
173
174PRIVATE short
175intersectLineSegments(const double  A[2],
176                      const double  B[2],
177                      const double  X[2],
178                      const double  Y[2],
179                      double        P[2])
180{
181  if ((X[0] < A[0] - EPSILON_7 && X[0] < B[0] - EPSILON_7
182       && Y[0] < A[0] - EPSILON_7 && Y[0] < B[0] - EPSILON_7)
183      ||
184      (X[0] > A[0] + EPSILON_7 && X[0] > B[0] + EPSILON_7
185       && Y[0] > A[0] + EPSILON_7 && Y[0] > B[0] + EPSILON_7)
186      )
187    /*
188     * Check if the x-coordinates of X and Y are smaller than
189     * the x-coordinates of A and B -> lines can not intersect
190     */
191    return 0;
192
193  if ((X[1] < A[1] - EPSILON_7 && X[1] < B[1] - EPSILON_7
194       && Y[1] < A[1] - EPSILON_7 && Y[1] < B[1] - EPSILON_7)
195      ||
196      (X[1] > A[1] + EPSILON_7 && X[1] > B[1] + EPSILON_7
197       && Y[1] > A[1] + EPSILON_7 && Y[1] > B[1] + EPSILON_7)
198      )
199    /*
200     * Check if the y-coordinates of X and Y are smaller than
201     * the y-coordinates of A and B -> lines can not intersect
202     */
203    return 0;
204
205  double denominator = (B[0] - A[0]) * (X[1] - Y[1]) - (B[1] - A[1]) * (X[0] - Y[0]);
206
207  if (fabs(denominator) < EPSILON_7) {
208    /*
209     * lines are parallel
210     * check if X is situated on AB line
211     */
212    double  sX, sY;
213
214    double  dx  = B[0] - A[0];
215    double  dy  = B[1] - A[1];
216
217    if (fabs(dx) > EPSILON_7) {
218      sX = (X[0] - A[0]) / (dx);
219      double refXy = A[1] + sX * (dy);
220      if (fabs(refXy - X[1]) > EPSILON_7)
221        /* AB and XY are not part of the same line */
222        return 0;
223
224      sY = (Y[0] - A[0]) / (dx);
225    } else {
226      sX = (X[1] - A[1]) / (dy);
227      double refXx = A[0] + sX * (dx);
228      if (fabs(refXx - X[0]) > EPSILON_7)
229        /* AB and XY are not part of the same line */
230        return 0;
231
232      sY = (Y[1] - A[1]) / (dy);
233    }
234
235    /* check if X or Y are situated directly on AB */
236    if ((0.0 <= sX && sX <= 1.0)
237        || (0.0 <= sY && sY <= 1.0))
238
239      return 1;
240
241    /* check if XY encloses AB */
242    if ((sX < 0.0 && 1.0 < sY)
243        || (sY < 0.0 && 1.0 < sX))
244
245      return 1;
246  } else {
247    /*
248     * lines are not parallel and might intersect
249     * (default case)
250     */
251    double  nominatorS  = (X[0] - Y[0]) * (A[1] - X[1]) - (X[1] - Y[1]) * (A[0] - X[0]);
252    double  nominatorT  = (A[0] - X[0]) * (B[1] - A[1]) - (A[1] - X[1]) * (B[0] - A[0]);
253    double  s           = nominatorS / denominator;
254    double  t           = nominatorT / denominator;
255
256    if (0.0 <= s && s <= 1.0 && 0.0 <= t && t <= 1.0) {
257      double  Ps[2];
258      Ps[0] = A[0] + s * (B[0] - A[0]);
259      Ps[1] = A[1] + s * (B[1] - A[1]);
260
261      double  Pt[2];
262      Pt[0] = X[0] + t * (Y[0] - X[0]);
263      Pt[1] = X[1] + t * (Y[1] - X[1]);
264
265      if (fabs(Ps[0] - Pt[0]) < EPSILON_7 && fabs(Ps[1] - Pt[1]) < EPSILON_7) {
266        if (P != NULL) {
267          P[0]  = Ps[0];
268          P[1]  = Ps[1];
269        }
270
271        return 1;
272      }
273    }
274  }
275
276  return 0;
277}
278
279
280PRIVATE short
281intersectLineArc(const double point_1[2],
282                 const double point_2[2],
283                 const double arc[6])
284{
285  char    *fnName = "intersectLineArc";
286  double  cut[2][2];
287  double  center[2] = {
288    arc[0], arc[1]
289  };
290  double  radius    = arc[2];
291  double  anchor[2] = {
292    point_1[0], point_1[1]
293  };
294  double  direction[2];
295
296  vector(point_1, point_2, direction);
297  short   num_points = getCutPointsOfCircleAndLine(center,
298                                                   radius,
299                                                   anchor,
300                                                   direction,
301                                                   cut[0],
302                                                   cut[1]);
303  short ret = 0;
304
305  for (int i = 0; i < num_points; i++) {
306    /* check if the computed intersection point is situated on the line segment */
307    double  A[2] = {
308      point_1[0], point_1[1]
309    };
310    double  B[2] = {
311      point_2[0], point_2[1]
312    };
313    double  line[2];
314    vector(A, B, line);
315    double  length = vectorLength2D(line);
316    double  v_A_cut[2];
317    double  v_B_cut[2];
318    vector(A, cut[i], v_A_cut);
319    vector(B, cut[i], v_B_cut);
320
321    if (fabs(length - vectorLength2D(v_A_cut) - vectorLength2D(v_B_cut)) > 0.01)
322
323      continue;
324
325    /* check if the computed intersection point is situated between angle_from and angle_to */
326    ret = ret || matchPointArc(cut[i], arc);
327
328    if (ret)
329
330      break;
331  }
332
333  return ret;
334}
335
336
337PRIVATE short
338intersectArcArc(const double  arc1[6],
339                const double  arc2[6])
340{
341  char    *fnName = "intersectArcArc";
342
343  double  c1[2] = {
344    arc1[0], arc1[1]
345  };
346  double  r1    = arc1[2];
347  double  c2[2] = {
348    arc2[0], arc2[1]
349  };
350  double  r2 = arc2[2];
351
352  if (!intersectCircleCircle(c1, r1, c2, r2))
353
354    return 0;
355
356  double  cut[2][2];
357  short   num_points  = getCutPointsOfCircles(c1, r1, c2, r2, cut[0], cut[1]);
358  short   ret         = 0;
359
360  for (int i = 0; i < num_points; i++) {
361    short hit1  = matchPointArc(cut[i], arc1);
362    short hit2  = matchPointArc(cut[i], arc2);
363
364    ret = ret || (hit1 && hit2);
365  }
366
367  return ret;
368}
369
370
371#endif
372