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