1#ifndef RNAPUZZLER_COORDINATE_TRANSFORM
2#define RNAPUZZLER_COORDINATE_TRANSFORM
3
4/*
5 *      RNAturtle algorithm
6 *
7 *      c  Daniel Wiegreffe, Daniel Alexander, Dirk Zeckzer
8 *      ViennaRNA package
9 */
10
11#include "../headers/tBaseInformation_struct.h"
12#include "vector_math.inc"
13#include "drawingconfig.inc"
14
15/**
16 * Compute affin angles for all loops
17 */
18PRIVATE void
19computeAffineCoordinates(short const *const       pair_table,
20                         const double             paired,
21                         const double             unpaired,
22                         tBaseInformation *const  baseInformation);
23
24
25/**
26 * Calculate the coordinates for the drawing with the given affin angles
27 */
28PRIVATE void
29affineToCartesianCoordinates(tBaseInformation *const  baseInformation,
30                             unsigned short const     length,
31                             double *const            x,
32                             double *const            y);
33
34
35/* -------------------------------------------------------------------------------------------------------------------------- */
36
37/**
38 * Handle bases of the exterior loop
39 *
40 * @returns first paired base
41 */
42PRIVATE short
43handleExteriorBases(const short *const  pair_table,
44                    short               currentBase,
45                    tBaseInformation    *baseInformation,
46                    const int           direction)
47{
48  short const length = pair_table[0];
49
50  if (currentBase > 1) {
51    baseInformation[currentBase].angle = baseInformation[currentBase].angle + direction *
52                                         MATH_PI_HALF;
53    baseInformation[currentBase].baseType = TYPE_EXTERIOR;
54  }
55
56  while (currentBase < length && pair_table[currentBase] <= 0) {
57    baseInformation[currentBase + 1].angle  = 0.0;
58    baseInformation[currentBase].baseType   = TYPE_EXTERIOR;
59    currentBase++;
60  }
61  /*
62   * For following stem, since it cannot alter the first two bases
63   * which may contain loop exit angles.
64   */
65  if ((currentBase + 1) <= (length)) {
66    baseInformation[currentBase + 1].angle  = direction * MATH_PI_HALF;
67    baseInformation[currentBase].baseType   = TYPE_EXTERIOR;
68  } else {
69    baseInformation[currentBase].baseType = TYPE_EXTERIOR;
70  }
71
72  return currentBase;
73}
74
75
76PRIVATE void
77countLoopPairs(short *const       m,
78               short *const       n,
79               short              i,
80               const short *const pair_table)
81{
82  short const end = pair_table[i++];  /* end is 1st base of 2nd stem half */
83
84  *m = *n = 1;                        /* Base pair of enclosing stem, consecutive pair from last stem base to first loop base */
85  while (i < end) {
86    if (pair_table[i] <= 0 || pair_table[i] < i) {
87      (*n)++;
88      i++;       /* Go to consecutive base */
89    } else {
90      /* base pair ahead */
91      (*m)++;
92      i = pair_table[i];       /* Go to base pair partner */
93    }
94  }
95
96  return;
97}
98
99
100/* -------------------------------------------------------------------------------------------------------------------------- */
101PRIVATE void
102handleStem(const short *const       pair_table,
103           short                    currentBase,
104           const double             paired,
105           const double             unpaired,
106           tBaseInformation *const  baseInformation,
107           int                      direction);
108
109
110/* -------------------------------------------------------------------------------------------------------------------------- */
111
112/*---------------------------------------------------------------------------*/
113
114/**
115 * Detect bulges (internal loops with only one unpaired base)
116 */
117PRIVATE int
118detectBulge(short               start,
119            const short *const  pair_table)
120{
121  int bulge   = 0;
122  int end     = pair_table[start];
123  int iterate = 1;
124  int old     = 0;
125  int i       = start + 1;
126  int side    = 0;
127
128  do {
129    if (pair_table[i] > 0) {
130      /* There was a basepair before */
131      if (iterate > 0) {
132        /* Stem entry */
133        if (pair_table[i] == old) {
134          side++;
135          i++;
136        }        /* Bulge found */
137        else {
138          /* It's only a bulge if the paired base is the start base or the bulge is on the other strand */
139          if (start == pair_table[i] || pair_table[i] == end - 2) {
140            bulge = pair_table[i];
141            break;
142          } else {
143            break;
144          }
145        }
146      }      /* Found first basepair */
147      else {
148        iterate++;
149        old = i;
150        i   = pair_table[i];
151      }
152    }    /* Consecutive Base found */
153    else {
154      /* There was a basepair before */
155      if (iterate > 0) {
156        iterate = 0;
157        i++;
158      }      /* Just the next consecutive base */
159      else {
160        i++;
161      }
162    }
163  } while (i > start);
164
165  return bulge;
166}
167
168
169PRIVATE void
170handleLoop(short              i,
171           const short *const pair_table,
172           const double       paired,
173           const double       unpaired,
174           tBaseInformation   *baseInformation,
175           int                direction)
176{
177  int         start = i;
178  short const end = pair_table[i];  /* end is 1st base of 2nd stem half */
179  short       m, n;                 /* m - no. of base pairs; n - no. of consecutive pairs */
180
181  countLoopPairs(&m, &n, i, pair_table);
182
183
184  int         bulge = detectBulge(i, pair_table);
185
186  if (bulge > 0 && n - m == 1) {
187    /* detected bulge with a single unpaired base */
188    int     length = (unpaired * (n - m + 1)) / 2;
189
190    double  alpha = acos(unpaired / (2 * length));
191
192    if (pair_table[i + 1] == 0) {
193      baseInformation[i + 1].angle            = baseInformation[i + 1].angle + direction * alpha;
194      baseInformation[i].baseType             = TYPE_BULGE;
195      baseInformation[pair_table[i]].baseType = TYPE_BULGE;
196      i++;
197
198      baseInformation[i + 1].angle  = -direction * alpha * 2;
199      baseInformation[i].baseType   = TYPE_BULGE;
200      i++;
201
202      baseInformation[i + 1].angle            = direction * alpha;
203      baseInformation[i].baseType             = TYPE_BULGE;
204      baseInformation[pair_table[i]].baseType = TYPE_BULGE;
205
206      handleStem(pair_table, i, paired, unpaired, baseInformation, direction);
207
208      i = pair_table[i];
209    } else {
210      baseInformation[i + 1].angle  = baseInformation[i + 1].angle + 0.0;
211      baseInformation[i].baseType   = TYPE_BULGE;
212      i++;
213
214      baseInformation[i + 1].angle    = baseInformation[i + 1].angle + 0.0;
215      baseInformation[i + 1].baseType = TYPE_BULGE;
216
217      baseInformation[i + 2].angle    = baseInformation[i + 2].angle + 0.0;
218      baseInformation[i + 1].baseType = TYPE_BULGE;
219
220      handleStem(pair_table, i, paired, unpaired, baseInformation, direction);
221
222      i = pair_table[i];
223
224      baseInformation[i + 1].angle  = baseInformation[i + 1].angle + direction * alpha;
225      baseInformation[i].baseType   = TYPE_BULGE;
226      i++;
227
228      baseInformation[i + 1].angle  = -direction * alpha * 2;
229      baseInformation[i].baseType   = TYPE_BULGE;
230      i++;
231
232      baseInformation[i + 1].angle  = direction * alpha;
233      baseInformation[i].baseType   = TYPE_BULGE;
234    }
235  } else {
236    /* loop is not bulge with a single unpaired base */
237
238    /*
239     * ******************************************************************
240     * ***************** Loop Drawing Algorithm - Start *****************
241     * ******************************************************************
242     */
243
244
245    /* variables for config drawing */
246    double  angle_over_paired;                /* alpha      at calculation */
247    double  current_angle;                    /* phi        at calculation */
248    double  current_bb_angle;                 /* beta       at calculation */
249    double  current_distance;                 /* b          at calculation */
250    double  current_delta_ab;                 /* delta_ab   at calculation */
251    double  current_delta_bb;                 /* delta_bb   at calculation */
252    int     current_stem_count;
253
254    config  *cfg        = baseInformation[start].config;
255    int     currentArc  = 0;
256    double  r           = cfg->radius;
257
258    angle_over_paired = 2 * asin(paired / (2 * r));
259
260    current_angle     = getArcAngle(cfg, currentArc);
261    current_bb_angle  = (current_angle - angle_over_paired) /
262                        (cfg->cfgArcs[currentArc]).numberOfArcSegments;
263    current_distance  = sqrt(2 * r * r * (1 - cos(current_bb_angle)));
264    current_delta_ab  = 0.5 * (MATH_PI + angle_over_paired + current_bb_angle);
265    current_delta_bb  = MATH_PI + current_bb_angle;
266    ++currentArc;
267
268    /* consider the direction that was given before (e.g. because of a bulge) */
269    baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction *
270                                   (MATH_PI - current_delta_ab);
271
272    baseInformation[i].distance = current_distance;
273
274    current_stem_count = 0;
275
276    /* Double Loop Check: Loop is followed directly by a loop */
277    if (baseInformation[i].baseType == TYPE_LOOP1)
278      baseInformation[i].baseType = TYPE_LOOP2;
279    else
280      baseInformation[i].baseType = TYPE_LOOP1;
281
282    i++;
283
284    while (i < end) {
285      if (pair_table[i] <= 0) {
286        /* unpaired */
287
288        /* handle backbone element */
289        baseInformation[i + 1].angle  = -direction * (current_delta_bb - MATH_PI);
290        baseInformation[i].distance   = current_distance;
291
292        baseInformation[i].baseType = TYPE_LOOP1;
293        i++;
294      } else if (pair_table[i] > i) {
295        /* base pair ahead */
296
297        /* i is the beginning of a stem and therefore the end of the current arc */
298        baseInformation[i + 1].angle = direction * (MATH_PI - current_delta_ab);
299        current_stem_count++;
300        baseInformation[i].baseType = TYPE_LOOP1;
301        handleStem(pair_table, i, paired, unpaired, baseInformation, direction);  /* Recurse multiloop */
302        i = pair_table[i];                                                        /* Go to base pair partner */
303      } else {
304        /* base pair behind */
305
306
307        /*
308         * i is the end of a stem that returned to the loop
309         * recalculation if angles and distances is nessecary because the next arc has been reached
310         */
311        if (current_stem_count == 1) {
312          current_stem_count  = 0;
313          current_angle       = getArcAngle(cfg, currentArc);
314          current_bb_angle    = (current_angle - angle_over_paired) /
315                                (cfg->cfgArcs[currentArc]).numberOfArcSegments;
316          current_distance  = sqrt(2 * r * r * (1 - cos(current_bb_angle)));
317          current_delta_ab  = 0.5 * (MATH_PI + angle_over_paired + current_bb_angle);
318          current_delta_bb  = MATH_PI + current_bb_angle;
319          ++currentArc;
320        }
321
322        /* consider the direction that was given before (e.g. because of a bulge) */
323        baseInformation[i + 1].angle = baseInformation[i + 1].angle + direction *
324                                       (MATH_PI - current_delta_ab);
325        baseInformation[i].distance = current_distance;
326        baseInformation[i].baseType = TYPE_LOOP1;
327
328        i++;
329      }
330    }
331
332    if ((i + 1) <= pair_table[0])
333      baseInformation[i + 1].angle = direction * (MATH_PI - current_delta_ab);
334
335    baseInformation[i].baseType = TYPE_LOOP1;
336
337
338    /*
339     * ******************************************************************
340     * ****************** Loop Drawing Algorithm - End ******************
341     * ******************************************************************
342     */
343  }
344
345  return;
346}
347
348
349/* -------------------------------------------------------------------------------------------------------------------------- */
350
351PRIVATE void
352handleStem(const short *const       pair_table,
353           short                    i,
354           const double             paired,
355           const double             unpaired,
356           tBaseInformation *const  baseInformation,
357           int                      direction)
358{
359  short end = pair_table[i] + 1;   /* first base after stem */
360
361  /*
362   * Skip first two bases of stem, angle already written by
363   * either loop (+exit angle) or handle_unpaired
364   * i+=2;
365   */
366  baseInformation[i].baseType = TYPE_STEM;
367  i++;
368
369  while (pair_table[i] > 0 &&                     /* i is base paired */
370         (pair_table[i] == end - 1 ||             /* First position of stem, continue! */
371          pair_table[i] + 1 == pair_table[i - 1]  /* Detect bulges on opposite strand */
372         )
373         ) {
374    baseInformation[i + 1].angle  = 0.0;
375    baseInformation[i].baseType   = TYPE_STEM;
376    i++;
377  }
378  if (pair_table[i] == end - 1) {
379  } else {
380    handleLoop(--i, pair_table, paired, unpaired, baseInformation, direction);     /* i - last base of first stem half */
381  }
382
383  i                           = pair_table[i]; /* set i as base pairing partner of last stem base */
384  baseInformation[i].baseType = TYPE_STEM;
385  i++;
386
387  while (i < end && i < pair_table[0]) {
388    baseInformation[i].baseType = TYPE_STEM;
389    i++;
390  }
391
392  return;
393}
394
395
396/* ------------------------------------------------------------------------------ */
397
398/**
399 * Compute angle angles for all loops
400 */
401PRIVATE void
402computeAffineCoordinates(short const *const       pair_table,
403                         const double             paired,
404                         const double             unpaired,
405                         tBaseInformation *const  baseInformation)
406{
407  /* Initialization */
408  short const length      = pair_table[0];
409  short       currentBase = 1;
410  const int   direction   = -1;
411
412  baseInformation[0].angle = 0.0;
413
414  /*
415   * For first stem, since handle_stem cannot alter the first two bases
416   * which may contain loop exit angles (in other stems).
417   */
418  if (2 <= length) {
419    baseInformation[1].angle  = baseInformation[0].angle;
420    baseInformation[2].angle  = baseInformation[1].angle;
421  }
422
423  /* reset dangle_count */
424  int dangle_count = 0;
425  while (currentBase < length) {
426    if (pair_table[currentBase] <= 0) {
427      if (currentBase > 1)
428
429        baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR;
430
431      currentBase = handleExteriorBases(pair_table, currentBase, baseInformation, direction);
432      /* returns first paired base currentBase */
433      dangle_count++;
434    }
435
436    if (currentBase < length) {
437      /* it was not the dangling end */
438      if (pair_table[currentBase] - pair_table[currentBase - 1] != 1
439          && pair_table[currentBase] != 0
440          && pair_table[currentBase - 1] != 0) {
441        if (currentBase == 1) {
442          if (dangle_count < 1) {
443            baseInformation[0].angle              =
444              baseInformation[1].angle            =
445                baseInformation[2].angle          = -MATH_PI_HALF;
446            baseInformation[currentBase].baseType = TYPE_EXTERIOR;
447          }
448
449          handleStem(pair_table, currentBase, paired, unpaired, baseInformation, direction);
450          currentBase = pair_table[currentBase] + 1;
451
452          /* Check, if there is a minor dangling end at the end */
453          if (currentBase == length) {
454            baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR;
455            baseInformation[currentBase].baseType     = TYPE_EXTERIOR;
456            baseInformation[currentBase].angle        = -MATH_PI_HALF;
457          }
458
459          continue;
460        } else {
461          baseInformation[currentBase].angle = baseInformation[currentBase].angle +
462                                               direction * MATH_PI_HALF;
463          baseInformation[currentBase + 1].distance = unpaired;
464          baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR;
465          baseInformation[currentBase + 1].angle    = baseInformation[currentBase + 1].angle +
466                                                      direction * MATH_PI_HALF;
467          baseInformation[currentBase].baseType = TYPE_EXTERIOR;
468
469          dangle_count++;
470        }
471      }
472
473      handleStem(pair_table, currentBase, paired, unpaired, baseInformation, direction);
474
475      currentBase = pair_table[currentBase] + 1;       /* currentBase is next base after stem */
476
477      if (currentBase == length) {
478        baseInformation[currentBase - 1].baseType = TYPE_EXTERIOR;
479        currentBase                               = handleExteriorBases(pair_table,
480                                                                        currentBase,
481                                                                        baseInformation,
482                                                                        direction);
483      }
484    }
485  }
486  baseInformation[length].baseType = TYPE_EXTERIOR;
487}
488
489
490/* ------------------------------------------------------------------------------ */
491
492/**
493 * Calculate the coordinates for the drawing with the given angle angles
494 */
495PRIVATE void
496affineToCartesianCoordinates(tBaseInformation *const  baseInformation,
497                             unsigned short const     length,
498                             double *const            x,
499                             double *const            y)
500{
501  if (length < 1)
502    return;
503
504  double angle = 0.0;
505  x[0] = y[0] = EXTERIOR_Y;
506  for (int i = 1; i < length; i++) {
507    angle = angle - baseInformation[i + 1].angle;
508
509    x[i]  = x[i - 1] + baseInformation[i].distance * cos(angle);
510    y[i]  = y[i - 1] + baseInformation[i].distance * sin(angle);
511  }
512  return;
513}
514
515
516#endif
517