1#ifndef RNAPUZZLER_BOUNDING_WEDGE_H
2#define RNAPUZZLER_BOUNDING_WEDGE_H
3
4/*
5 *      RNApuzzler wedges
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 "ViennaRNA/utils/basic.h"
15
16#include "../headers/configtree_struct.h"
17#include "configtree.inc"
18#include "boundingBoxes.inc"
19#include "vector_math.inc"
20
21/**
22 * @brief getBoundingWedge
23 *      calculates the bounding wedge for root's i-th child
24 * @param root
25 *      center node of the bounding wedge
26 * @param childIndex
27 *      index of root's child for which the bounding wedge is calculated
28 * @param minAngle
29 *      pointer used for returning the minimum angle
30 * @param maxAngle
31 *      pointer used for returning the maximum angle
32 */
33PRIVATE void
34getBoundingWedge(const treeNode *root,
35                 int            childIndex,
36                 double         *minAngle,
37                 double         *maxAngle);
38
39
40PRIVATE void
41getBoundingWedgeRec(const treeNode  *root,
42                    const treeNode  *node,
43                    const double    parentAngle,
44                    double          *minAngle,
45                    double          *maxAngle)
46{
47  /*
48   * --- Documentation ---
49   *
50   * How to get the bounding wedge of root's i-th child tree?
51   *
52   * get interesting points of current node
53   * (2) touch points of tangents from root.center to node's loop circle
54   * (n) bulge points of node's stem
55   * (2) corners of node's stem that coincide with root's loop [only for direct child]
56   *
57   * update min-angle
58   * for each interesting point
59   * check if the point is on the left side of the min-angle axis
60   * if so
61   * min-angle -= diff-angle of min-axis and point-axis
62   * update max-angle
63   * for each interesting point
64   * check if the point is on the right side of the min-angle axis
65   * if so
66   * max-angle += diff-angle of max-axis and point-axis
67   */
68
69  const double  distance  = epsilonFix;
70  treeNode      *parent   = getParent(node);
71
72  double        centerRoot[2];
73
74  getLoopCenter(root, centerRoot);
75  double        centerNode[2];
76  getLoopCenter(node, centerNode);
77  double        vRootNode[2];
78  vector(centerRoot, centerNode, vRootNode);
79
80  /*
81   * set appropriate nodeAngle
82   * this could have been done using getChildAngle function
83   * but in terms of performance we can get this for free O(1)
84   * costs of getChildAngle: O( maxDegreeOnPath * ( depth(node) - depth(root) ) ) per call
85   */
86  double nodeAngle;
87
88  if (parent == root) {
89    /*
90     * this happens only for the initial call and not for the recursive calls
91     * initialize min/max with the direct child's angle
92     */
93    nodeAngle = getChildAngle(root, node);
94    *minAngle = nodeAngle;
95    *maxAngle = nodeAngle;
96  } else {
97    /* compare getChildAngle function */
98    double  centerParent[2];
99    getLoopCenter(parent, centerParent);
100    double  vRootParent[2];
101    vector(centerRoot, centerParent, vRootParent);
102    double  diffParent = angleBetweenVectors2D(vRootParent, vRootNode);
103
104    if (!isToTheRightPointVector(centerRoot, vRootParent, centerNode))
105
106      diffParent *= -1;
107
108    nodeAngle = parentAngle + diffParent;
109  }
110
111  /* get all bounding boxes */
112  loopBox *loopNode = node->lBox;
113  stemBox *stemNode = node->sBox;
114
115  /* allocate space for points of interest */
116  int     numPoints = stemNode->bulgeCount;
117
118  if (parent == root)
119    /* for bottom corners of direct child's stem */
120    numPoints += 2;
121
122  double  **points    = (double **)vrna_alloc(numPoints * sizeof(double *));
123  int     pointIndex  = 0;
124
125  /*
126   * points of interest (part 1)
127   * bulge points
128   */
129  for (int i = 0; i < stemNode->bulgeCount; i++) {
130    double  o[2], q[2];    /* o, q are unused but necessary for function call */
131    double  *bulgePoint = (double *)vrna_alloc(2 * sizeof(double));
132    getBulgeCoordinatesExtraDistance(stemNode, i, distance, o, bulgePoint, q);
133
134    points[pointIndex] = bulgePoint;
135    pointIndex++;
136  }
137
138  /*
139   * points of interest (part 2, only for direct child of root)
140   * for the direct child of this computation's root
141   * the corners of that child that coincide with root's loop (-> stem bottom corners)
142   * have a huge effect on the size of the bounding wedge
143   * for all further descendants these points do not matter because
144   * of the greater impact of the loops
145   */
146  if (parent == root) {
147    double *pStemBottomCornerL = (double *)vrna_alloc(2 * sizeof(double));
148    pStemBottomCornerL[0] = stemNode->c[0] - stemNode->e[0] * stemNode->a[0] + stemNode->e[1] *
149                            stemNode->b[0];
150    pStemBottomCornerL[1] = stemNode->c[1] - stemNode->e[0] * stemNode->a[1] + stemNode->e[1] *
151                            stemNode->b[1];
152    points[pointIndex] = pStemBottomCornerL;
153    pointIndex++;
154
155    double *pStemBottomCornerR = (double *)vrna_alloc(2 * sizeof(double));
156    pStemBottomCornerR[0] = stemNode->c[0] - stemNode->e[0] * stemNode->a[0] - stemNode->e[1] *
157                            stemNode->b[0];
158    pStemBottomCornerR[1] = stemNode->c[1] - stemNode->e[0] * stemNode->a[1] - stemNode->e[1] *
159                            stemNode->b[1];
160    points[pointIndex] = pStemBottomCornerR;
161    pointIndex++;
162  }
163
164  /*
165   * we compute the two tangents from root's center to node's circle periphery
166   * using these for our min/max calculation ensures that the whole loop
167   * is contained in the wedge
168   *
169   * we can directly compute the diffAngle for the touching points of the tangents
170   * by using pythagoras' sentence
171   */
172  double  radiusNode        = loopNode->r + distance;
173  double  distanceRootNode  = vectorLength2D(vRootNode);
174
175  /* positive angle and negative angle share their size */
176  double  angle1 = asin(radiusNode / distanceRootNode);
177
178  /* no need for double computation, just flip sign */
179  double  angle2 = (-1) * angle1;
180
181  /* store both angles in an array to avoid code duplication */
182  double  angles[2] = {
183    angle1, angle2
184  };
185
186  /* update min / max for the tangents touching points (angles) */
187  for (int currentAngle = 0; currentAngle < 2; currentAngle++) {
188    double  diffAngle   = angles[currentAngle];
189    double  pointAngle  = nodeAngle + diffAngle;
190
191    /* actual updating for min and max */
192    if (pointAngle < *minAngle)
193
194      *minAngle = pointAngle;
195
196    if (pointAngle > *maxAngle)
197
198      *maxAngle = pointAngle;
199  }
200
201  /* update min / max for bulge points (and stem bottom points in first level of recursion) */
202  for (int currentPoint = 0; currentPoint < numPoints; currentPoint++) {
203    double  *point = points[currentPoint];
204
205    /* vector from root.loop.center to point */
206    double  vCenterPoint[2];
207    vector(centerRoot, point, vCenterPoint);
208
209    /* (positive) offset angle between node.loop.center and point */
210    double  diffAngle = angleBetweenVectors2D(vRootNode, vCenterPoint);
211    double  sign;
212
213    if (isToTheRightPointVector(centerRoot, vRootNode, point))
214      sign = 1;
215    else
216      sign = -1;
217
218    /* now the offset angle has some direction information */
219    diffAngle *= sign;
220
221    /* the current point's angle has to be set with regard to the node's angle */
222    double pointAngle = nodeAngle + diffAngle;
223
224    /* actual updating for min and max */
225    if (pointAngle < *minAngle)
226
227      *minAngle = pointAngle;
228
229    if (pointAngle > *maxAngle)
230
231      *maxAngle = pointAngle;
232  }
233
234  /* free allocated space */
235  for (int currentPoint = 0; currentPoint < numPoints; currentPoint++) {
236    double *point = points[currentPoint];
237    free(point);
238  }
239  free(points);
240
241  /* recursive call */
242  for (int currentChild = 0; currentChild < node->childCount; currentChild++) {
243    treeNode *child = getChild(node, currentChild);
244    getBoundingWedgeRec(root, child, nodeAngle, minAngle, maxAngle);
245  }
246}
247
248
249PRIVATE void
250getBoundingWedge(const treeNode *root,
251                 const int      childIndex,
252                 double         *minAngle,
253                 double         *maxAngle)
254{
255  treeNode *child = getChild(root, childIndex);
256
257  getBoundingWedgeRec(root, child, 0, minAngle, maxAngle);
258}
259
260
261#endif
262