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