srcjar/fr/orsay/lri/varna/models/rna/DrawRNATemplate.java

   1 /**
   2  * File written by Raphael Champeimont
   3  * UMR 7238 Genomique des Microorganismes
   4  */
   5 package fr.orsay.lri.varna.models.rna;
   6
   7 import java.awt.geom.AffineTransform;
   8 import java.awt.geom.Line2D;
   9 import java.awt.geom.Point2D;
  10 import java.util.ArrayList;
  11 import java.util.Arrays;
  12 import java.util.HashMap;
  13 import java.util.HashSet;
  14 import java.util.Iterator;
  15 import java.util.Map;
  16 import java.util.Set;
  17
  18 import fr.orsay.lri.varna.exceptions.ExceptionInvalidRNATemplate;
  19 import fr.orsay.lri.varna.models.VARNAConfig;
  20 import fr.orsay.lri.varna.models.geom.ComputeArcCenter;
  21 import fr.orsay.lri.varna.models.geom.ComputeEllipseAxis;
  22 import fr.orsay.lri.varna.models.geom.CubicBezierCurve;
  23 import fr.orsay.lri.varna.models.geom.HalfEllipse;
  24 import fr.orsay.lri.varna.models.geom.LinesIntersect;
  25 import fr.orsay.lri.varna.models.geom.MiscGeom;
  26 import fr.orsay.lri.varna.models.templates.DrawRNATemplateCurveMethod;
  27 import fr.orsay.lri.varna.models.templates.DrawRNATemplateMethod;
  28 import fr.orsay.lri.varna.models.templates.RNANodeValueTemplate;
  29 import fr.orsay.lri.varna.models.templates.RNANodeValueTemplateBasePair;
  30 import fr.orsay.lri.varna.models.templates.RNATemplate;
  31 import fr.orsay.lri.varna.models.templates.RNATemplateAlign;
  32 import fr.orsay.lri.varna.models.templates.RNATemplateDrawingAlgorithmException;
  33 import fr.orsay.lri.varna.models.templates.RNATemplateMapping;
  34 import fr.orsay.lri.varna.models.templates.RNATemplate.EdgeEndPointPosition;
  35 import fr.orsay.lri.varna.models.templates.RNATemplate.In1Is;
  36 import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateElement;
  37 import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateHelix;
  38 import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateUnpairedSequence;
  39 import fr.orsay.lri.varna.models.templates.RNATemplate.RNATemplateElement.EdgeEndPoint;
  40 import fr.orsay.lri.varna.models.treealign.Tree;
  41
  42 public class DrawRNATemplate {
  43         private RNA rna;
  44         private RNATemplateMapping mapping;
  45         private ArrayList<ModeleBase> _listeBases;
  46
  47
  48         public DrawRNATemplate(RNA rna) {
  49                 this.rna = rna;
  50                 this._listeBases = rna.getListeBases();
  51         }
  52
  53         public RNATemplateMapping getMapping() {
  54                 return mapping;
  55         }
  56
  57
  58
  59         /**
  60          * Draw this RNA like the given template.
  61          * The helixLengthAdjustmentMethod argument tells what to do in case
  62          * some helices are of a different length in the template and the
  63          * actual helix. See class DrawRNATemplateMethod above for possible values.
  64          * @param straightBulges
  65          */
  66         public void drawRNATemplate(
  67                         RNATemplate template,
  68                         VARNAConfig conf,
  69                         DrawRNATemplateMethod helixLengthAdjustmentMethod,
  70                         DrawRNATemplateCurveMethod curveMethod, boolean straightBulges)
  71         throws RNATemplateDrawingAlgorithmException {
  72
  73                 // debug
  74 //              try {
  75 //                      RNA perfectMatchingRNA = template.toRNA();
  76 //                      System.out.println("An RNA that would perfectly match this template would be:");
  77 //                      System.out.println(perfectMatchingRNA.getStructDBN());
  78 //              } catch (ExceptionInvalidRNATemplate e) {
  79 //                      e.printStackTrace();
  80 //              }
  81
  82                 mapping = RNATemplateAlign.mapRNAWithTemplate(rna, template);
  83                 //System.out.println(mapping.showCompact(this));
  84
  85                 // debug
  86 //                      RNATemplateAlign.printMapping(mapping, template, getSeq());
  87 //                      try {
  88 //                              TreeGraphviz.treeToGraphvizPostscript(alignment, "alignment_graphviz.ps");
  89 //                      } catch (IOException e) {
  90 //                              e.printStackTrace();
  91 //                      }
  92
  93
  94
  95                 Iterator<RNATemplateElement> iter;
  96                 double globalIncreaseFactor = 1;
  97                 Map<RNATemplateHelix, Point2D.Double> translateVectors = null;
  98                 if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.MAXSCALINGFACTOR) {
  99                         // Compute increase factors for helices.
 100                         Map<RNATemplateHelix, Double> lengthIncreaseFactor = new HashMap<RNATemplateHelix, Double>();
 101                         double maxLengthIncreaseFactor = Double.NEGATIVE_INFINITY;
 102                         //RNATemplateHelix maxIncreaseHelix = null;
 103                         iter = template.rnaIterator();
 104                         while (iter.hasNext()) {
 105                                 RNATemplateElement element = iter.next();
 106                                 if (element instanceof RNATemplateHelix
 107                                                 && mapping.getAncestor(element) != null
 108                                                 && !lengthIncreaseFactor.containsKey(element)) {
 109                                         RNATemplateHelix helix = (RNATemplateHelix) element;
 110                                         int[] basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
 111                                         Arrays.sort(basesInHelixArray);
 112                                         double l = computeLengthIncreaseFactor(basesInHelixArray, helix, straightBulges);
 113                                         lengthIncreaseFactor.put(helix, l);
 114                                         if (l > maxLengthIncreaseFactor) {
 115                                                 maxLengthIncreaseFactor = l;
 116                                                 //maxIncreaseHelix = helix;
 117                                         }
 118                                 }
 119                         }
 120
 121                         // debug
 122                         //System.out.println("Max helix length increase factor = " + maxLengthIncreaseFactor + " reached with helix " + maxIncreaseHelix);;
 123
 124                         globalIncreaseFactor = Math.max(1, maxLengthIncreaseFactor);
 125
 126                 } else if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.HELIXTRANSLATE) {
 127                         try {
 128                                 // Now we need to propagate this helices translations
 129                                 Tree<RNANodeValueTemplate> templateAsTree = template.toTree();
 130                                 translateVectors = computeHelixTranslations(templateAsTree, mapping, straightBulges);
 131
 132                         } catch (ExceptionInvalidRNATemplate e) {
 133                                 throw (new RNATemplateDrawingAlgorithmException("ExceptionInvalidRNATemplate: " + e.getMessage()));
 134                         }
 135                 }
 136
 137                 // Allocate the coords and centers arrays
 138                 // We create Point2D.Double objects in it but the algorithms
 139                 // we use may choose to create new Point2D.Double objects or to
 140                 // modify those created here.
 141                 Point2D.Double[] coords = new Point2D.Double[_listeBases.size()];
 142                 Point2D.Double[] centers = new Point2D.Double[_listeBases.size()];
 143                 double[] angles = new double[_listeBases.size()];
 144                 for (int i = 0; i < _listeBases.size(); i++) {
 145                         coords[i] = new Point2D.Double(0, 0);
 146                         centers[i] = new Point2D.Double(0, 0);
 147                 }
 148
 149                 boolean computeCoords = true;
 150                 while (computeCoords) {
 151                         computeCoords = false;
 152                         // Compute coords and centers
 153                         Set<RNATemplateHelix> alreadyDrawnHelixes = new HashSet<RNATemplateHelix>();
 154                         RNATemplateHelix lastMappedHelix = null;
 155                         EdgeEndPoint howWeGotOutOfLastHelix = null;
 156                         int howWeGotOutOfLastHelixBaseIndex = 0;
 157                         iter = template.rnaIterator();
 158                         RNATemplateElement element = null;
 159                         while (iter.hasNext()) {
 160                                 element = iter.next();
 161                                 if (element instanceof RNATemplateHelix
 162                                                 && mapping.getAncestor(element) != null) {
 163                                         // We have a mapping between an helix in the RNA sequence
 164                                         // and an helix in the template.
 165
 166                                         RNATemplateHelix helix = (RNATemplateHelix) element;
 167                                         boolean firstTimeWeMeetThisHelix;
 168                                         int[] basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
 169                                         Arrays.sort(basesInHelixArray);
 170
 171                                         // Draw this helix if it has not already been done
 172                                         if (!alreadyDrawnHelixes.contains(helix)) {
 173                                                 firstTimeWeMeetThisHelix = true;
 174                                                 drawHelixLikeTemplateHelix(basesInHelixArray, helix, coords, centers,angles, globalIncreaseFactor, translateVectors, straightBulges);
 175                                                 alreadyDrawnHelixes.add(helix);
 176                                         } else {
 177                                                 firstTimeWeMeetThisHelix = false;
 178                                         }
 179
 180                                         EdgeEndPoint howWeGetInCurrentHelix;
 181                                         if (firstTimeWeMeetThisHelix) {
 182                                                 if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
 183                                                         howWeGetInCurrentHelix = helix.getIn1();
 184                                                 } else {
 185                                                         howWeGetInCurrentHelix = helix.getIn2();
 186                                                 }
 187                                         } else {
 188                                                 if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
 189                                                         howWeGetInCurrentHelix = helix.getIn2();
 190                                                 } else {
 191                                                         howWeGetInCurrentHelix = helix.getIn1();
 192                                                 }
 193                                         }
 194
 195                                         Point2D.Double P0 = new Point2D.Double();
 196                                         Point2D.Double P3 = new Point2D.Double();
 197
 198                                         if (lastMappedHelix != null) {
 199                                                 // Now draw the RNA sequence (possibly containing helixes)
 200                                                 // between the last template drawn helix and this one.
 201
 202                                                 if (lastMappedHelix == helix) {
 203                                                         // Last helix is the same as the current one so
 204                                                         // nothing matched (or at best a single
 205                                                         // non-paired sequence) so we will just
 206                                                         // use the Radiate algorithm
 207
 208                                                         Point2D.Double helixVector = new Point2D.Double();
 209                                                         computeHelixEndPointDirections(howWeGotOutOfLastHelix, helixVector, new Point2D.Double());
 210
 211                                                         double angle = MiscGeom.angleFromVector(helixVector);
 212                                                         int b1 = basesInHelixArray[basesInHelixArray.length/2 - 1];
 213                                                         P0.setLocation(coords[b1]);
 214                                                         int b2 = basesInHelixArray[basesInHelixArray.length/2];
 215                                                         P3.setLocation(coords[b2]);
 216                                                         Point2D.Double loopCenter = new Point2D.Double((P0.x + P3.x)/2, (P0.y + P3.y)/2);
 217                                                         rna.drawLoop(b1,
 218                                                                          b2,
 219                                                                          loopCenter.x,
 220                                                                          loopCenter.y,
 221                                                                          angle,
 222                                                                          coords,
 223                                                                          centers,
 224                                                                          angles,
 225                                                                          straightBulges);
 226                                                         // If the helix is flipped, we need to compute the symmetric
 227                                                         // of the whole loop.
 228                                                         if (helix.isFlipped()) {
 229                                                                 symmetric(loopCenter, helixVector, coords, b1, b2);
 230                                                                 symmetric(loopCenter, helixVector, centers, b1, b2);
 231                                                         }
 232                                                 } else {
 233                                                         // No helices matched between the last helix and
 234                                                         // the current one, so we draw what is between
 235                                                         // using the radiate algorithm but on the Bezier curve.
 236
 237                                                         int b1 = howWeGotOutOfLastHelixBaseIndex;
 238                                                         int b2 = firstTimeWeMeetThisHelix ? basesInHelixArray[0] : basesInHelixArray[basesInHelixArray.length/2];
 239                                                         P0.setLocation(coords[b1]);
 240                                                         P3.setLocation(coords[b2]);
 241
 242                                                         Point2D.Double P1, P2;
 243
 244                                                         if (howWeGotOutOfLastHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
 245                                                                         && howWeGetInCurrentHelix.getOtherElement() instanceof RNATemplateUnpairedSequence) {
 246                                                                 // We will draw the bases on a Bezier curve
 247                                                                 P1 = new Point2D.Double();
 248                                                                 computeBezierTangentVectorTarget(howWeGotOutOfLastHelix, P0, P1);
 249
 250                                                                 P2 = new Point2D.Double();
 251                                                                 computeBezierTangentVectorTarget(howWeGetInCurrentHelix, P3, P2);
 252                                                         } else {
 253                                                                 // We will draw the bases on a straight line between P0 and P3
 254                                                                 P1 = null;
 255                                                                 P2 = null;
 256                                                         }
 257
 258                                                         drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, lastMappedHelix.isFlipped(), straightBulges);
 259                                                 }
 260
 261                                         } else if (basesInHelixArray[0] > 0) {
 262                                                 // Here we draw what is before the first mapped helix.
 263
 264                                                 RNATemplateUnpairedSequence templateSequence;
 265                                                 // Try to find our template sequence as the mapped element of base 0
 266                                                 RNATemplateElement templateSequenceCandidate = mapping.getPartner(0);
 267                                                 if (templateSequenceCandidate != null
 268                                                                 && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 269                                                         templateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 270                                                 } else {
 271                                                         // Try other idea: first template element if it is a sequence
 272                                                         templateSequenceCandidate = template.getFirst();
 273                                                         if (templateSequenceCandidate != null
 274                                                                         && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 275                                                                 templateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 276                                                         } else {
 277                                                                 // We don't know where to start
 278                                                                 templateSequence = null;
 279                                                         }
 280                                                 }
 281
 282                                                 int b1 = 0;
 283                                                 int b2 = firstTimeWeMeetThisHelix ? basesInHelixArray[0] : basesInHelixArray[basesInHelixArray.length/2];
 284                                                 P3.setLocation(coords[b2]);
 285
 286                                                 if (templateSequence != null) {
 287                                                         coords[0].setLocation(templateSequence.getVertex5());
 288                                                         coords[0].x *= globalIncreaseFactor;
 289                                                         coords[0].y *= globalIncreaseFactor;
 290                                                 } else {
 291                                                         // Put b1 at an ideal distance from b2, using the "flat exterior loop" method
 292                                                         double idealLength = computeStraightLineIdealLength(b1, b2);
 293                                                         Point2D.Double j = new Point2D.Double();
 294                                                         if (howWeGetInCurrentHelix != null) {
 295                                                                 computeHelixEndPointDirections(howWeGetInCurrentHelix, new Point2D.Double(), j);
 296                                                         } else {
 297                                                                 j.setLocation(1, 0);
 298                                                         }
 299                                                         coords[b1].setLocation(coords[b2].x + j.x*idealLength, coords[b2].y + j.y*idealLength);
 300                                                 }
 301                                                 P0.setLocation(coords[0]);
 302
 303                                                 Point2D.Double P1, P2;
 304
 305                                                 if (howWeGetInCurrentHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
 306                                                                 && templateSequence != null) {
 307                                                         // We will draw the bases on a Bezier curve
 308                                                         P1 = new Point2D.Double();
 309                                                         computeBezierTangentVectorTarget(templateSequence.getIn(), P0, P1);
 310
 311                                                         P2 = new Point2D.Double();
 312                                                         computeBezierTangentVectorTarget(howWeGetInCurrentHelix, P3, P2);
 313                                                 } else {
 314                                                         // We will draw the bases on a straight line between P0 and P3
 315                                                         P1 = null;
 316                                                         P2 = null;
 317                                                 }
 318
 319                                                 drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, false, straightBulges);
 320                                         }
 321
 322                                         lastMappedHelix = helix;
 323                                         howWeGotOutOfLastHelix = howWeGetInCurrentHelix.getNextEndPoint();
 324                                         if (firstTimeWeMeetThisHelix) {
 325                                                 howWeGotOutOfLastHelixBaseIndex = basesInHelixArray[basesInHelixArray.length/2-1];
 326                                         } else {
 327                                                 howWeGotOutOfLastHelixBaseIndex = basesInHelixArray[basesInHelixArray.length-1];
 328                                         }
 329                                 }
 330                         } // end template iteration
 331
 332
 333                         // Now we need to draw what is after the last mapped helix.
 334                         if (howWeGotOutOfLastHelixBaseIndex < coords.length-1
 335                                         && element != null
 336                                         && coords.length > 1) {
 337
 338                                 RNATemplateUnpairedSequence beginTemplateSequence = null;
 339                                 if (lastMappedHelix == null) {
 340                                         // No helix at all matched between the template and RNA!
 341                                         // So the sequence we want to draw is the full RNA.
 342
 343                                         // Try to find our template sequence as the mapped element of base 0
 344                                         RNATemplateElement templateSequenceCandidate = mapping.getPartner(0);
 345                                         if (templateSequenceCandidate != null
 346                                                         && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 347                                                 beginTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 348                                         } else {
 349                                                 // Try other idea: first template element if it is a sequence
 350                                                 templateSequenceCandidate = template.getFirst();
 351                                                 if (templateSequenceCandidate != null
 352                                                                 && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 353                                                         beginTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 354                                                 } else {
 355                                                         // We don't know where to start
 356                                                         beginTemplateSequence = null;
 357                                                 }
 358                                         }
 359
 360                                         if (beginTemplateSequence != null) {
 361                                                 coords[0].setLocation(beginTemplateSequence.getVertex5());
 362                                                 coords[0].x *= globalIncreaseFactor;
 363                                                 coords[0].y *= globalIncreaseFactor;
 364                                         }
 365
 366                                 }
 367
 368                                 RNATemplateUnpairedSequence endTemplateSequence;
 369                                 // Try to find our template sequence as the mapped element of last base
 370                                 RNATemplateElement templateSequenceCandidate = mapping.getPartner(coords.length-1);
 371                                 if (templateSequenceCandidate != null
 372                                                 && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 373                                         endTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 374                                 } else {
 375                                         // Try other idea: last template element if it is a sequence
 376                                         templateSequenceCandidate = element;
 377                                         if (templateSequenceCandidate != null
 378                                                         && templateSequenceCandidate instanceof RNATemplateUnpairedSequence) {
 379                                                 endTemplateSequence = (RNATemplateUnpairedSequence) templateSequenceCandidate;
 380                                         } else {
 381                                                 // We don't know where to end
 382                                                 endTemplateSequence = null;
 383                                         }
 384                                 }
 385
 386                                 int b1 = howWeGotOutOfLastHelixBaseIndex;
 387                                 int b2 = coords.length - 1;
 388
 389                                 if (endTemplateSequence != null) {
 390                                         coords[b2].setLocation(endTemplateSequence.getVertex3());
 391                                         coords[b2].x *= globalIncreaseFactor;
 392                                         coords[b2].y *= globalIncreaseFactor;
 393                                 } else {
 394                                         // Put b2 at an ideal distance from b1, using the "flat exterior loop" method
 395                                         double idealLength = computeStraightLineIdealLength(b1, b2);
 396                                         Point2D.Double j = new Point2D.Double();
 397                                         if (howWeGotOutOfLastHelix != null) {
 398                                                 computeHelixEndPointDirections(howWeGotOutOfLastHelix, new Point2D.Double(), j);
 399                                         } else {
 400                                                 j.setLocation(1, 0);
 401                                         }
 402                                         coords[b2].setLocation(coords[b1].x + j.x*idealLength, coords[b1].y + j.y*idealLength);
 403                                 }
 404
 405
 406                                 Point2D.Double P0 = new Point2D.Double();
 407                                 Point2D.Double P3 = new Point2D.Double();
 408
 409                                 P0.setLocation(coords[b1]);
 410                                 P3.setLocation(coords[b2]);
 411
 412                                 Point2D.Double P1, P2;
 413
 414                                 if (howWeGotOutOfLastHelix != null
 415                                                 && howWeGotOutOfLastHelix.getOtherElement() instanceof RNATemplateUnpairedSequence
 416                                                 && endTemplateSequence != null) {
 417                                         // We will draw the bases on a Bezier curve
 418                                         P1 = new Point2D.Double();
 419                                         computeBezierTangentVectorTarget(howWeGotOutOfLastHelix, P0, P1);
 420
 421                                         P2 = new Point2D.Double();
 422                                         computeBezierTangentVectorTarget(endTemplateSequence.getOut(), P3, P2);
 423                                 } else if (lastMappedHelix == null
 424                                                 && beginTemplateSequence != null
 425                                                 && endTemplateSequence != null) {
 426                                         // We will draw the bases on a Bezier curve
 427                                         P1 = new Point2D.Double();
 428                                         computeBezierTangentVectorTarget(beginTemplateSequence.getIn(), P0, P1);
 429
 430                                         P2 = new Point2D.Double();
 431                                         computeBezierTangentVectorTarget(endTemplateSequence.getOut(), P3, P2);
 432                                 } else {
 433                                         // We will draw the bases on a straight line between P0 and P3
 434                                         P1 = null;
 435                                         P2 = null;
 436                                 }
 437
 438                                 drawAlongCurve(b1, b2, P0, P1, P2, P3, coords, centers, angles, curveMethod, lastMappedHelix != null ? lastMappedHelix.isFlipped() : false, straightBulges);
 439
 440                         }
 441
 442
 443                         if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.NOINTERSECT && coords.length > 3) {
 444                                 // Are we happy with this value of globalIncreaseFactor?
 445                                 Line2D.Double[] lines = new Line2D.Double[coords.length-1];
 446                                 for (int i=0; i<coords.length-1; i++) {
 447                                         lines[i] = new Line2D.Double(coords[i], coords[i+1]);
 448                                 }
 449                                 int intersectLines = 0;
 450                                 for (int i=0; i<lines.length; i++) {
 451                                         for (int j=i+2; j<lines.length; j++) {
 452                                                 if (LinesIntersect.linesIntersect(lines[i], lines[j])) {
 453                                                         intersectLines++;
 454                                                 }
 455                                         }
 456                                 }
 457                                 // If no intersection we keep this globalIncreaseFactor value
 458                                 if (intersectLines > 0) {
 459                                         // Don't increase more than a maximum value
 460                                         if (globalIncreaseFactor < 3) {
 461                                                 globalIncreaseFactor += 0.1;
 462                                                 //System.out.println("globalIncreaseFactor increased to " + globalIncreaseFactor);
 463                                                 // Compute the drawing again
 464                                                 computeCoords = true;
 465                                         }
 466                                 }
 467                         }
 468
 469                 }
 470
 471                 // debug
 472                 if (helixLengthAdjustmentMethod == DrawRNATemplateMethod.MAXSCALINGFACTOR
 473                                 || helixLengthAdjustmentMethod == DrawRNATemplateMethod.NOINTERSECT) {
 474                         //System.out.println("globalIncreaseFactor = " + globalIncreaseFactor);
 475                 }
 476
 477                 // Now we actually move the bases, according to arrays coords and centers
 478                 // and taking in account the space between bases parameter.
 479                 for (int i = 0; i < _listeBases.size(); i++) {
 480                         _listeBases.get(i).setCoords(
 481                                         new Point2D.Double(coords[i].x * conf._spaceBetweenBases,
 482                                                         coords[i].y * conf._spaceBetweenBases));
 483                         _listeBases.get(i).setCenter(
 484                                         new Point2D.Double(centers[i].x * conf._spaceBetweenBases,
 485                                                         centers[i].y * conf._spaceBetweenBases));
 486                 }
 487
 488         }
 489
 490
 491
 492
 493
 494
 495         /**
 496          * IN: Argument helixEndPoint is an IN argument (will be read),
 497          * and must contain an helix edge endpoint.
 498          *
 499          * The other arguments are OUT arguments
 500          * (must be existing objects, content will be overwritten).
 501          *
 502          * OUT: The i argument will contain a vector colinear to the vector
 503          * from the helix startPosition to endPosition or the opposite
 504          * depending on there the endpoint is (the endpoint will be on the
 505          * destination side of the vector). ||i|| = 1
 506          *
 507          * OUT: The j vector will contain an vector that is colinear
 508          * to the last/first base pair connection on the side of this endpoint.
 509          * The vector will be oriented to the side of the given endpoint.
 510          * ||j|| = 1
 511          */
 512         private void computeHelixEndPointDirections(
 513                         EdgeEndPoint helixEndPoint, // IN
 514                         Point2D.Double i, // OUT
 515                         Point2D.Double j // OUT
 516                         ) {
 517                 RNATemplateHelix helix = (RNATemplateHelix) helixEndPoint.getElement();
 518                 Point2D.Double startpos = helix.getStartPosition();
 519                 Point2D.Double endpos = helix.getEndPosition();
 520                 Point2D.Double helixVector = new Point2D.Double();
 521                 switch (helixEndPoint.getPosition()) {
 522                 case IN1:
 523                 case OUT2:
 524                         helixVector.x = startpos.x - endpos.x;
 525                         helixVector.y = startpos.y - endpos.y;
 526                         break;
 527                 case IN2:
 528                 case OUT1:
 529                         helixVector.x = endpos.x - startpos.x;
 530                         helixVector.y = endpos.y - startpos.y;
 531                         break;
 532                 }
 533                 double helixVectorLength = Math.hypot(helixVector.x, helixVector.y);
 534                 // i is the vector which is colinear to helixVector and such that ||i|| = 1
 535                 i.x = helixVector.x / helixVectorLength;
 536                 i.y = helixVector.y / helixVectorLength;
 537                 // Find j such that it is orthogonal to i, ||j|| = 1
 538                 // and j goes to the side where the sequence will be connected
 539                 switch (helixEndPoint.getPosition()) {
 540                 case IN1:
 541                 case IN2:
 542                         // rotation of +pi/2
 543                         j.x = - i.y;
 544                         j.y =   i.x;
 545                         break;
 546                 case OUT1:
 547                 case OUT2:
 548                         // rotation of -pi/2
 549                         j.x =   i.y;
 550                         j.y = - i.x;
 551                         break;
 552                 }
 553                 if (helix.isFlipped()) {
 554                         j.x = - j.x;
 555                         j.y = - j.y;
 556                 }
 557
 558         }
 559
 560         /**
 561          * A cubic Bezier curve can be defined by 4 points,
 562          * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
 563          * For each of the curve end points, there is the last/first point of the
 564          * curve and a point which gives the direction and length of the tangent
 565          * vector on that side. This two points are respectively curveEndPoint
 566          * and curveVectorOtherPoint.
 567          * IN:  Argument helixVector is the vector formed by the helix,
 568          *      in the right direction for our sequence.
 569          * IN:  Argument curveEndPoint is the position of the endpoint on the helix.
 570          * OUT: Argument curveVectorOtherPoint must be allocated
 571          *      and the values will be modified.
 572          */
 573         private void computeBezierTangentVectorTarget(
 574                         EdgeEndPoint endPoint,
 575                         Point2D.Double curveEndPoint,
 576                         Point2D.Double curveVectorOtherPoint)
 577                         throws RNATemplateDrawingAlgorithmException {
 578
 579                 boolean sequenceEndPointIsIn;
 580                 RNATemplateUnpairedSequence sequence;
 581
 582                 if (endPoint.getElement() instanceof RNATemplateHelix) {
 583                         sequence = (RNATemplateUnpairedSequence) endPoint.getOtherElement();
 584                         EdgeEndPointPosition endPointPositionOnHelix = endPoint.getPosition();
 585                         switch (endPointPositionOnHelix) {
 586                         case IN1:
 587                         case IN2:
 588                                 sequenceEndPointIsIn = false;
 589                                 break;
 590                         default:
 591                                 sequenceEndPointIsIn = true;
 592                         }
 593
 594                         EdgeEndPoint endPointOnHelix =
 595                                 sequenceEndPointIsIn ?
 596                                                 sequence.getIn().getOtherEndPoint() :
 597                                                 sequence.getOut().getOtherEndPoint();
 598                         if (endPointOnHelix == null) {
 599                                 throw (new RNATemplateDrawingAlgorithmException("Sequence is not connected to an helix."));
 600                         }
 601                 } else {
 602                         // The endpoint is on an unpaired sequence.
 603                         sequence = (RNATemplateUnpairedSequence) endPoint.getElement();
 604                         if (endPoint == sequence.getIn()) {
 605                                 // endpoint is 5'
 606                                 sequenceEndPointIsIn = true;
 607                         } else {
 608                                 sequenceEndPointIsIn = false;
 609                         }
 610                 }
 611
 612                 double l =
 613                         sequenceEndPointIsIn ?
 614                                 sequence.getInTangentVectorLength() :
 615                                 sequence.getOutTangentVectorLength();
 616
 617                 // Compute the absolute angle our line makes to the helix
 618                 double theta =
 619                         sequenceEndPointIsIn ?
 620                                 sequence.getInTangentVectorAngle() :
 621                                 sequence.getOutTangentVectorAngle();
 622
 623                 // Compute v, the tangent vector of the Bezier curve
 624                 Point2D.Double v = new Point2D.Double();
 625                 v.x = l * Math.cos(theta);
 626                 v.y = l * Math.sin(theta);
 627                 curveVectorOtherPoint.x = curveEndPoint.x + v.x;
 628                 curveVectorOtherPoint.y = curveEndPoint.y + v.y;
 629         }
 630
 631
 632         /**
 633          * Compute (actual helix length / helix length in template).
 634          * @param straightBulges
 635          */
 636         private double computeLengthIncreaseFactor(
 637                         int[] basesInHelixArray,  // IN
 638                         RNATemplateHelix helix,    // IN
 639                         boolean straightBulges
 640                         ) {
 641                 double templateLength = computeHelixTemplateLength(helix);
 642                 double realLength = computeHelixRealLength(basesInHelixArray,straightBulges);
 643                 return realLength / templateLength;
 644         }
 645
 646         /**
 647          * Compute (actual helix vector - helix vector in template).
 648          */
 649         private Point2D.Double computeLengthIncreaseDelta(
 650                         int[] basesInHelixArray,  // IN
 651                         RNATemplateHelix helix,   // IN
 652                         boolean straightBulges
 653                         ) {
 654                 double templateLength = computeHelixTemplateLength(helix);
 655                 double realLength = computeHelixRealLength(basesInHelixArray,straightBulges);
 656                 Point2D.Double i = new Point2D.Double();
 657                 computeTemplateHelixVectors(helix, null, i, null);
 658                 return new Point2D.Double(i.x*(realLength-templateLength), i.y*(realLength-templateLength));
 659         }
 660
 661         /**
 662          * Compute helix interesting vectors from template helix.
 663          * @param helix The template helix you want to compute the vectors from.
 664          * @param o This point coordinates will be set the origin of the helix (or not if null),
 665          *          ie. the point in the middle of the base pair with the two most extreme bases.
 666          * @param i Will be set to the normalized helix vector. (nothing done if null)
 667          * @param j Will be set to the normalized helix base pair vector (5' -> 3'). (nothing done if null)
 668          */
 669         private void computeTemplateHelixVectors(
 670                         RNATemplateHelix helix,  // IN
 671                         Point2D.Double o,        // OUT
 672                         Point2D.Double i,        // OUT
 673                         Point2D.Double j         // OUT
 674                         ) {
 675                 Point2D.Double startpos, endpos;
 676                 if (helix.getIn1Is() == In1Is.IN1_IS_5PRIME) {
 677                         startpos = helix.getStartPosition();
 678                         endpos = helix.getEndPosition();
 679                 } else {
 680                         endpos = helix.getStartPosition();
 681                         startpos = helix.getEndPosition();
 682                 }
 683                 if (o != null) {
 684                         o.x = startpos.x;
 685                         o.y = startpos.y;
 686                 }
 687                 if (i != null || j != null) {
 688                         // (i_x,i_y) is the vector between two consecutive bases of the same side of an helix
 689                         if (i == null)
 690                                 i = new Point2D.Double();
 691                         i.x = (endpos.x - startpos.x);
 692                         i.y = (endpos.y - startpos.y);
 693                         double i_original_norm = Math.hypot(i.x, i.y);
 694                         // change its norm to 1
 695                         i.x = i.x / i_original_norm;
 696                         i.y = i.y / i_original_norm;
 697                         if (j != null) {
 698                                 j.x = - i.y;
 699                                 j.y =   i.x;
 700                                 if (helix.isFlipped()) {
 701                                         j.x = - j.x;
 702                                         j.y = - j.y;
 703                                 }
 704                                 double j_original_norm = Math.hypot(j.x, j.y);
 705                                 // change (j_x,j_y) so that its norm is 1
 706                                 j.x = j.x / j_original_norm;
 707                                 j.y = j.y / j_original_norm;
 708                         }
 709                 }
 710         }
 711
 712
 713         /**
 714          * Estimate bulge arc length.
 715          */
 716         private double estimateBulgeArcLength(int firstBase, int lastBase) {
 717                 if (firstBase + 1 == lastBase)
 718                         return RNA.LOOP_DISTANCE; // there is actually no bulge
 719                 double len = 0.0;
 720                 int k = firstBase;
 721                 while (k < lastBase) {
 722                         int l = _listeBases.get(k).getElementStructure();
 723                         if (k < l && l < lastBase) {
 724                                 len += RNA.BASE_PAIR_DISTANCE;
 725                                 k = l;
 726                         } else {
 727                                 len += RNA.LOOP_DISTANCE;
 728                                 k++;
 729                         }
 730                 }
 731                 return len;
 732         }
 733
 734
 735         /**
 736          * Estimate bulge width, the given first and last bases must be those in the helix.
 737          */
 738         private double estimateBulgeWidth(int firstBase, int lastBase) {
 739                 double len = estimateBulgeArcLength(firstBase, lastBase);
 740                 return 2 * (len / Math.PI);
 741         }
 742
 743
 744         /**
 745          * Get helix length in template.
 746          */
 747         private double computeHelixTemplateLength(RNATemplateHelix helix) {
 748                 return Math.hypot(helix.getStartPosition().x - helix.getEndPosition().x,
 749                                 helix.getStartPosition().y - helix.getEndPosition().y);
 750         }
 751
 752
 753         /**
 754          * Compute helix actual length (as drawHelixLikeTemplateHelix() would draw it).
 755          */
 756         private double computeHelixRealLength(int[] basesInHelixArray, boolean straightBulges) {
 757                 return drawHelixLikeTemplateHelix(basesInHelixArray, null, null, null, null, 0, null,straightBulges);
 758         }
 759
 760
 761         /**
 762          * Result type of countUnpairedLine(), see below.
 763          */
 764         private static class UnpairedLineCounts {
 765                 public int nBP, nLD, total;
 766         }
 767         /**
 768          * If we are drawing an unpaired region that may contains helices,
 769          * and we are drawing it on a line (curve or straight, doesn't matter),
 770          * how many intervals should have a base-pair length (start of an helix)
 771          * and how many should have a consecutive unpaired bases length?
 772          * Returns an array with three elements:
 773          * - answer to first question
 774          * - answer to second question
 775          * - sum of both, ie. total number of intervals
 776          *   (this is NOT lastBase-firstBase because bases deep in helices do not count)
 777          */
 778         private UnpairedLineCounts countUnpairedLine(int firstBase, int lastBase) {
 779                 UnpairedLineCounts counts = new UnpairedLineCounts();
 780                 int nBP = 0;
 781                 int nLD = 0;
 782                 {
 783                         int b = firstBase;
 784                         while (b < lastBase) {
 785                                 int l = _listeBases.get(b).getElementStructure();
 786                                 if (b < l && l < lastBase) {
 787                                         nBP++;
 788                                         b = l;
 789                                 } else {
 790                                         nLD++;
 791                                         b++;
 792                                 }
 793                         }
 794                 }
 795                 counts.nBP = nBP;
 796                 counts.nLD = nLD;
 797                 counts.total = nBP + nLD;
 798                 return counts;
 799         }
 800
 801
 802         /**
 803          * Draw the given helix (given as a *SORTED* array of indexes)
 804          * like defined in the given template helix.
 805          * OUT: The bases positions are not changed in fact,
 806          *      instead the coords and centers arrays are modified.
 807          * IN:  The helix origin position is multiplied by scaleHelixOrigin
 808          *      and translateVectors.get(helix) is added.
 809          * RETURN VALUE:
 810          *      The length of the drawn helix.
 811          * @param straightBulges
 812          *
 813          */
 814         private double drawHelixLikeTemplateHelix(
 815                         int[] basesInHelixArray,  // IN
 816                         RNATemplateHelix helix,   // IN  (optional, ie. may be null)
 817                         Point2D.Double[] coords,  // OUT (optional, ie. may be null)
 818                         Point2D.Double[] centers, // OUT (optional, ie. may be null)
 819                         double[] angles,  // OUT
 820                         double scaleHelixOrigin,  // IN
 821                         Map<RNATemplateHelix, Point2D.Double> translateVectors // IN (optional, ie. may be null)
 822 , boolean straightBulges
 823                         ) {
 824                 int n = basesInHelixArray.length / 2;
 825                 if (n == 0)
 826                         return 0;
 827                  // Default values when not template helix is provided:
 828                 Point2D.Double o = new Point2D.Double(0, 0);
 829                 Point2D.Double i = new Point2D.Double(1, 0);
 830                 Point2D.Double j = new Point2D.Double(0, 1);
 831                 boolean flipped = false;
 832                 if (helix != null) {
 833                         computeTemplateHelixVectors(helix, o, i, j);
 834                         flipped = helix.isFlipped();
 835                 }
 836                 Point2D.Double li = new Point2D.Double(i.x*RNA.LOOP_DISTANCE, i.y*RNA.LOOP_DISTANCE);
 837                 // We want o to be the point where the first base (5' end) is
 838                 o.x = (o.x - j.x * RNA.BASE_PAIR_DISTANCE / 2) * scaleHelixOrigin;
 839                 o.y = (o.y - j.y * RNA.BASE_PAIR_DISTANCE / 2) * scaleHelixOrigin;
 840                 if (translateVectors != null && translateVectors.containsKey(helix)) {
 841                         Point2D.Double v = translateVectors.get(helix);
 842                         o.x = o.x + v.x;
 843                         o.y = o.y + v.y;
 844                 }
 845
 846                 // We need this array so that we can store positions even if coords == null
 847                 Point2D.Double[] helixBasesPositions = new Point2D.Double[basesInHelixArray.length];
 848                 for (int k=0; k<helixBasesPositions.length; k++) {
 849                         helixBasesPositions[k] = new Point2D.Double();
 850                 }
 851                 Point2D.Double accDelta = new Point2D.Double(0, 0);
 852                 for (int k=0; k<n; k++) {
 853                         int kp = 2*n-k-1;
 854                         Point2D.Double p1 = helixBasesPositions[k]; // we assign the point *reference*
 855                         Point2D.Double p2 = helixBasesPositions[kp];
 856                         // Do we have a bulge between previous base pair and this one?
 857                         boolean bulge = k >= 1 && (basesInHelixArray[k] != basesInHelixArray[k-1] + 1
 858                                                            || basesInHelixArray[kp+1] != basesInHelixArray[kp] + 1);
 859                         if (k >= 1) {
 860                                 if (basesInHelixArray[k] < basesInHelixArray[k-1]
 861                                     || basesInHelixArray[kp+1] < basesInHelixArray[kp]) {
 862                                         throw new Error("Internal bug: basesInHelixArray must be sorted");
 863                                 }
 864                                 if (bulge) {
 865                                         // Estimate a good distance (delta) between the previous base pair and this one
 866                                         double delta1 = estimateBulgeWidth(basesInHelixArray[k-1], basesInHelixArray[k]);
 867                                         double delta2 = estimateBulgeWidth(basesInHelixArray[kp], basesInHelixArray[kp+1]);
 868                                         // The biggest bulge defines the width
 869                                         double delta = Math.max(delta1, delta2);
 870
 871                                         if (coords != null) {
 872                                                 // Now, where do we put the bases that are part of the bulge?
 873                                                 for (int side=0; side<2; side++) {
 874                                                         Point2D.Double pstart = new Point2D.Double();
 875                                                         Point2D.Double pend = new Point2D.Double();
 876                                                         Point2D.Double bisectVect = new Point2D.Double();
 877                                                         Point2D.Double is = new Point2D.Double();
 878                                                         int firstBase, lastBase;
 879                                                         double alphasign = flipped ? -1 : 1;
 880                                                         if (side == 0) {
 881                                                                 firstBase = basesInHelixArray[k-1];
 882                                                                 lastBase  = basesInHelixArray[k];
 883                                                                 pstart.setLocation(o.x + accDelta.x,
 884                                                                                    o.y + accDelta.y);
 885                                                                 pend.setLocation(o.x + accDelta.x + i.x*delta,
 886                                                                                  o.y + accDelta.y + i.y*delta);
 887                                                                 bisectVect.setLocation(-j.x, -j.y);
 888                                                                 is.setLocation(i);
 889                                                         } else {
 890                                                                 firstBase = basesInHelixArray[kp];
 891                                                                 lastBase  = basesInHelixArray[kp+1];
 892                                                                 pstart.setLocation(o.x + accDelta.x + i.x*delta + j.x*RNA.BASE_PAIR_DISTANCE,
 893                                                                            o.y + accDelta.y + i.y*delta + j.y*RNA.BASE_PAIR_DISTANCE);
 894                                                                 pend.setLocation(o.x + accDelta.x + j.x*RNA.BASE_PAIR_DISTANCE,
 895                                                                                  o.y + accDelta.y + j.y*RNA.BASE_PAIR_DISTANCE);
 896
 897                                                                 bisectVect.setLocation(j);
 898                                                                 is.setLocation(-i.x, -i.y);
 899                                                         }
 900                                                         double arclen = estimateBulgeArcLength(firstBase, lastBase);
 901                                                         double centerOnBisect = ComputeArcCenter.computeArcCenter(delta, arclen);
 902
 903                                                         // Should we draw the base on an arc or simply use a line?
 904                                                         if (centerOnBisect > -1000) {
 905                                                                 Point2D.Double center = new Point2D.Double(pstart.x + is.x*delta/2 + bisectVect.x*centerOnBisect,
 906                                                                                                            pstart.y + is.y*delta/2 + bisectVect.y*centerOnBisect);
 907                                                                 int b = firstBase;
 908                                                                 double len = 0;
 909                                                                 double r = Math.hypot(pstart.x - center.x, pstart.y - center.y);
 910                                                                 double alpha0 = MiscGeom.angleFromVector(pstart.x - center.x, pstart.y - center.y);
 911                                                                 while (b < lastBase) {
 912                                                                         int l = _listeBases.get(b).getElementStructure();
 913                                                                         if (b < l && l < lastBase) {
 914                                                                                 len += RNA.BASE_PAIR_DISTANCE;
 915                                                                                 b = l;
 916                                                                         } else {
 917                                                                                 len += RNA.LOOP_DISTANCE;
 918                                                                                 b++;
 919                                                                         }
 920                                                                         if (b < lastBase) {
 921                                                                                 coords[b].x = center.x + r*Math.cos(alpha0 + alphasign*len/r);
 922                                                                                 coords[b].y = center.y + r*Math.sin(alpha0 + alphasign*len/r);
 923                                                                         }
 924                                                                 }
 925                                                         } else {
 926                                                                 // Draw on a line
 927
 928                                                                 // Distance between paired bases cannot be changed
 929                                                                 // (imposed by helix width) but distance between other
 930                                                                 // bases can be adjusted.
 931                                                                 UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
 932                                                                 double LD = Math.max((delta - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
 933                                                                 //System.out.println("nBP=" + nBP + " nLD=" + nLD);
 934                                                                 double len = 0;
 935                                                                 {
 936                                                                         int b = firstBase;
 937                                                                         while (b < lastBase) {
 938                                                                                 int l = _listeBases.get(b).getElementStructure();
 939                                                                                 if (b < l && l < lastBase) {
 940                                                                                         len += RNA.BASE_PAIR_DISTANCE;
 941                                                                                         b = l;
 942                                                                                 } else {
 943                                                                                         len += LD;
 944                                                                                         b++;
 945                                                                                 }
 946                                                                                 if (b < lastBase) {
 947                                                                                         coords[b].x = pstart.x + is.x*len;
 948                                                                                         coords[b].y = pstart.y + is.y*len;
 949                                                                                 }
 950                                                                         }
 951                                                                 }
 952                                                                 //System.out.println("len=" + len + " delta=" + delta + " d(pstart,pend)=" + Math.hypot(pend.x-pstart.x, pend.y-pstart.y));
 953                                                         }
 954
 955                                                         // Does the bulge contain an helix?
 956                                                         // If so, use drawLoop() to draw it.
 957                                                         {
 958                                                                 int b = firstBase;
 959                                                                 while (b < lastBase) {
 960                                                                         int l = _listeBases.get(b).getElementStructure();
 961                                                                         if (b < l && l < lastBase) {
 962                                                                                 // Helix present in bulge
 963                                                                                 Point2D.Double b1pos = coords[b];
 964                                                                                 Point2D.Double b2pos = coords[l];
 965                                                                                 double beta = MiscGeom.angleFromVector(b2pos.x - b1pos.x, b2pos.y - b1pos.y) - Math.PI / 2 + (flipped ? Math.PI : 0);
 966                                                                                 Point2D.Double loopCenter = new Point2D.Double((b1pos.x + b2pos.x)/2, (b1pos.y + b2pos.y)/2);
 967                                                                                 rna.drawLoop(b,
 968                                                                                                  l,
 969                                                                                                  loopCenter.x,
 970                                                                                                  loopCenter.y,
 971                                                                                                  beta,
 972                                                                                                  coords,
 973                                                                                                  centers,
 974                                                                                                  angles,
 975                                                                                                  straightBulges);
 976                                                                                 // If the helix is flipped, we need to compute the symmetric
 977                                                                                 // of the whole loop.
 978                                                                                 if (helix.isFlipped()) {
 979                                                                                         Point2D.Double v = new Point2D.Double(Math.cos(beta), Math.sin(beta));
 980                                                                                         symmetric(loopCenter, v, coords, b, l);
 981                                                                                         symmetric(loopCenter, v, centers, b, l);
 982                                                                                 }
 983                                                                                 // Continue
 984                                                                                 b = l;
 985                                                                         } else {
 986                                                                                 b++;
 987                                                                         }
 988                                                                 }
 989                                                         }
 990                                                 }
 991                                         }
 992
 993                                         accDelta.x += i.x*delta;
 994                                         accDelta.y += i.y*delta;
 995                                         p1.x = o.x + accDelta.x;
 996                                         p1.y = o.y + accDelta.y;
 997                                         p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
 998                                         p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
 999
1000                                 } else {
1001                                         accDelta.x += li.x;
1002                                         accDelta.y += li.y;
1003                                         p1.x = o.x + accDelta.x;
1004                                         p1.y = o.y + accDelta.y;
1005                                         p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
1006                                         p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
1007                                 }
1008                         } else {
1009                                 // First base pair
1010                                 p1.x = o.x;
1011                                 p1.y = o.y;
1012                                 p2.x = p1.x + j.x*RNA.BASE_PAIR_DISTANCE;
1013                                 p2.y = p1.y + j.y*RNA.BASE_PAIR_DISTANCE;
1014                         }
1015                 }
1016
1017                 Point2D.Double p1 = helixBasesPositions[0];
1018                 Point2D.Double p2 = helixBasesPositions[n-1];
1019
1020                 if (coords != null) {
1021                         for (int k=0; k<helixBasesPositions.length; k++) {
1022                                 coords[basesInHelixArray[k]] = helixBasesPositions[k];
1023                         }
1024                 }
1025
1026                 return Math.hypot(p2.x-p1.x, p2.y-p1.y);
1027         }
1028
1029
1030
1031
1032         /**
1033          * Return ideal length to draw unpaired region from firstBase to lastBase.
1034          */
1035         private double computeStraightLineIdealLength(int firstBase, int lastBase) {
1036                 UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
1037                 return RNA.BASE_PAIR_DISTANCE*counts.nBP + RNA.LOOP_DISTANCE*counts.nLD;
1038         }
1039
1040         /**
1041          * Like drawOnBezierCurve(), but on a straight line.
1042          */
1043         private boolean drawOnStraightLine(
1044                         int firstBase,
1045                         int lastBase,
1046                         Point2D.Double P0,
1047                         Point2D.Double P3,
1048                         Point2D.Double[] coords,
1049                         Point2D.Double[] centers,
1050                         boolean cancelIfNotGood) {
1051
1052                 Point2D.Double vector = new Point2D.Double(P3.x - P0.x, P3.y - P0.y);
1053                 double vectorNorm = Math.hypot(vector.x, vector.y);
1054
1055                 UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
1056                 double LD = Math.max((vectorNorm - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
1057
1058                 if (cancelIfNotGood && LD < RNA.LOOP_DISTANCE*0.75) {
1059                         // Bases would be drawn "too" near (0.75 is heuristic)
1060                         return false;
1061                 }
1062
1063                 double len = 0;
1064                 {
1065                         int b = firstBase;
1066                         while (b <= lastBase) {
1067                                 coords[b].x = P0.x + vector.x*len/vectorNorm;
1068                                 coords[b].y = P0.y + vector.y*len/vectorNorm;
1069                                 int l = _listeBases.get(b).getElementStructure();
1070                                 if (b < l && l < lastBase) {
1071                                         len += RNA.BASE_PAIR_DISTANCE;
1072                                         b = l;
1073                                 } else {
1074                                         len += LD;
1075                                         b++;
1076                                 }
1077                         }
1078                 }
1079
1080                 return true;
1081         }
1082
1083
1084
1085
1086         /**
1087          * A Bezier curve can be defined by four points,
1088          * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
1089          * Here we give this four points and an array of bases indexes
1090          * (which must be indexes in this RNA sequence) which will be moved
1091          * to be on the Bezier curve.
1092          * The bases positions are not changed in fact, instead the coords and
1093          * centers arrays are modified.
1094          * If cancelIfNotGood, draw only if it would look good,
1095          * otherwise just return false.
1096          */
1097         private boolean drawOnBezierCurve(
1098                         int firstBase,
1099                         int lastBase,
1100                         Point2D.Double P0,
1101                         Point2D.Double P1,
1102                         Point2D.Double P2,
1103                         Point2D.Double P3,
1104                         Point2D.Double[] coords,
1105                         Point2D.Double[] centers,
1106                         boolean cancelIfNotGood) {
1107
1108                 UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
1109
1110                 // Draw the bases of the sequence along a Bezier curve
1111                 int n = counts.total;
1112
1113                 // We choose to approximate the Bezier curve by 10*n straight lines.
1114                 CubicBezierCurve bezier = new CubicBezierCurve(P0, P1, P2, P3, 10*n);
1115                 double curveLength = bezier.getApproxCurveLength();
1116
1117
1118                 double LD = Math.max((curveLength - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
1119
1120                 if (cancelIfNotGood && LD < RNA.LOOP_DISTANCE*0.75) {
1121                         // Bases would be drawn "too" near (0.75 is heuristic)
1122                         return false;
1123                 }
1124
1125                 double[] t = new double[n+1];
1126
1127                 {
1128                         double len = 0;
1129                         int k = 0;
1130                         int b = firstBase;
1131                         while (b <= lastBase) {
1132                                 t[k] = len;
1133                                 k++;
1134
1135                                 int l = _listeBases.get(b).getElementStructure();
1136                                 if (b < l && l < lastBase) {
1137                                         len += RNA.BASE_PAIR_DISTANCE;
1138                                         b = l;
1139                                 } else {
1140                                         len += LD;
1141                                         b++;
1142                                 }
1143                         }
1144                 }
1145
1146                 Point2D.Double[] sequenceBasesCoords = bezier.uniformParam(t);
1147
1148                 {
1149                         int k = 0;
1150                         int b = firstBase;
1151                         while (b <= lastBase) {
1152                                 coords[b] = sequenceBasesCoords[k];
1153                                 int l = _listeBases.get(b).getElementStructure();
1154                                 if (b < l && l < lastBase) {
1155                                         b = l;
1156                                 } else {
1157                                         b++;
1158                                 }
1159                                 k++;
1160                         }
1161                 }
1162
1163                 return true;
1164         }
1165
1166
1167
1168         private void drawOnEllipse(
1169                         int firstBase,
1170                         int lastBase,
1171                         Point2D.Double P0,
1172                         Point2D.Double P3,
1173                         Point2D.Double[] coords,
1174                         Point2D.Double[] centers,
1175                         boolean reverse) {
1176
1177                 UnpairedLineCounts counts = countUnpairedLine(firstBase, lastBase);
1178                 double halfEllipseLength = RNA.BASE_PAIR_DISTANCE*counts.nBP + RNA.LOOP_DISTANCE*counts.nLD;
1179                 double fullEllipseLength = halfEllipseLength*2;
1180
1181                 double axisA = P0.distance(P3)/2;
1182                 double axisB = ComputeEllipseAxis.computeEllipseAxis(axisA, fullEllipseLength);
1183
1184                 if (axisB == 0) {
1185                         // Ellipse is in fact a line, and it is much faster to draw it as such.
1186                         drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, false);
1187                         return;
1188                 }
1189
1190                 //System.out.println("Ellipse a=" + axisA + " b=" + axisB);
1191
1192                 int n = counts.total;
1193
1194                 // We choose to approximate the Bezier curve by 10*n straight lines.
1195                 HalfEllipse curve = new HalfEllipse(axisA, axisB, 10*n);
1196                 double curveLength = curve.getApproxCurveLength();
1197
1198
1199                 double LD = Math.max((curveLength - counts.nBP*RNA.BASE_PAIR_DISTANCE) / (double) counts.nLD, 0);
1200
1201                 double[] t = new double[n+1];
1202
1203                 {
1204                         double len = 0;
1205                         int k = 0;
1206                         int b = firstBase;
1207                         while (b <= lastBase) {
1208                                 t[k] = len;
1209                                 k++;
1210
1211                                 int l = _listeBases.get(b).getElementStructure();
1212                                 if (b < l && l < lastBase) {
1213                                         len += RNA.BASE_PAIR_DISTANCE;
1214                                         b = l;
1215                                 } else {
1216                                         len += LD;
1217                                         b++;
1218                                 }
1219                         }
1220                 }
1221
1222                 // Ellipse with axis = X,Y
1223                 Point2D.Double[] sequenceBasesCoords = curve.uniformParam(t);
1224
1225                 // Compute the symmetric half-ellipse
1226                 if (reverse) {
1227                         AffineTransform tranform1 = new AffineTransform();
1228                         tranform1.scale(1, -1);
1229                         tranform1.transform(sequenceBasesCoords, 0, sequenceBasesCoords, 0, sequenceBasesCoords.length);
1230                 }
1231
1232                 // Translate + rotate such that ellipse is at the right place
1233                 AffineTransform tranform = HalfEllipse.matchAxisA(P0, P3);
1234                 tranform.transform(sequenceBasesCoords, 0, sequenceBasesCoords, 0, sequenceBasesCoords.length);
1235
1236                 {
1237                         int k = 0;
1238                         int b = firstBase;
1239                         while (b <= lastBase) {
1240                                 coords[b] = sequenceBasesCoords[k];
1241                                 //coords[b] = curve.standardParam((double) k/n);
1242                                 //System.out.println(b + " " + coords[b]);
1243                                 int l = _listeBases.get(b).getElementStructure();
1244                                 if (b < l && l < lastBase) {
1245                                         b = l;
1246                                 } else {
1247                                         b++;
1248                                 }
1249                                 k++;
1250                         }
1251                 }
1252
1253         }
1254
1255         /**
1256          * This functions draws the RNA sequence between (including)
1257          * firstBase and lastBase along a curve.
1258          * The sequence may contain helices.
1259          *
1260          * Bezier curve:
1261          * A Bezier curve can be defined by four points,
1262          * see http://en.wikipedia.org/wiki/Bezier_curve#Cubic_B.C3.A9zier_curves
1263          *
1264          * Straight line:
1265          * If P1 and P2 are null, the bases are drawn on a straight line.
1266          *
1267          * OUT: The bases positions are not changed in fact,
1268          * instead the coords and centers arrays are modified.
1269          *
1270          * Argument reverse says if we should reverse the direction of inserted helices.
1271          */
1272         private void drawAlongCurve(
1273                         int firstBase,
1274                         int lastBase,
1275                         Point2D.Double P0,
1276                         Point2D.Double P1,
1277                         Point2D.Double P2,
1278                         Point2D.Double P3,
1279                         Point2D.Double[] coords,
1280                         Point2D.Double[] centers,
1281                         double[] angles,
1282                         DrawRNATemplateCurveMethod curveMethod,
1283                         boolean reverse,
1284                         boolean straightBulges) {
1285
1286                 // First we find the bases which are directly on the Bezier curve
1287                 ArrayList<Integer> alongBezierCurve = new ArrayList<Integer>();
1288                 for (int depth=0, i=firstBase; i<=lastBase; i++) {
1289                         int k = _listeBases.get(i).getElementStructure();
1290                         if (k < 0 || k > lastBase || k < firstBase) {
1291                                 if (depth == 0) {
1292                                         alongBezierCurve.add(i);
1293                                 }
1294                         } else {
1295                                 if (i < k) {
1296                                         if (depth == 0) {
1297                                                 alongBezierCurve.add(i);
1298                                                 alongBezierCurve.add(k);
1299                                         }
1300                                         depth++;
1301                                 } else {
1302                                         depth--;
1303                                 }
1304                         }
1305                 }
1306                 // number of bases along the Bezier curve
1307                 int n = alongBezierCurve.size();
1308                 int[] alongBezierCurveArray = RNATemplateAlign.intArrayFromList(alongBezierCurve);
1309                 if (n > 0) {
1310                         if (curveMethod == DrawRNATemplateCurveMethod.ALWAYS_REPLACE_BY_ELLIPSES) {
1311                                 drawOnEllipse(firstBase, lastBase, P0, P3, coords, centers, reverse);
1312                         } else if (curveMethod == DrawRNATemplateCurveMethod.SMART) {
1313                                 boolean passed;
1314                                 if (P1 != null && P2 != null) {
1315                                         passed = drawOnBezierCurve(firstBase, lastBase, P0, P1, P2, P3, coords, centers, true);
1316                                 } else {
1317                                         passed = drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, true);
1318                                 }
1319                                 if (!passed) {
1320                                         drawOnEllipse(firstBase, lastBase, P0, P3, coords, centers, reverse);
1321                                 }
1322                         } else {
1323                                 if (P1 != null && P2 != null) {
1324                                         drawOnBezierCurve(firstBase, lastBase, P0, P1, P2, P3, coords, centers, false);
1325                                 } else {
1326                                         drawOnStraightLine(firstBase, lastBase, P0, P3, coords, centers, false);
1327                                 }
1328                         }
1329                 }
1330                 // Now use the radiate algorithm to draw the helixes
1331                 for (int k=0; k<n-1; k++) {
1332                         int b1 = alongBezierCurveArray[k];
1333                         int b2 = alongBezierCurveArray[k+1];
1334                         if (_listeBases.get(b1).getElementStructure() == b2) {
1335                                 Point2D.Double b1pos = coords[b1];
1336                                 Point2D.Double b2pos = coords[b2];
1337                                 double alpha = MiscGeom.angleFromVector(b2pos.x - b1pos.x, b2pos.y - b1pos.y);
1338                                 rna.drawLoop(b1,
1339                                                  b2,
1340                                                  (b1pos.x + b2pos.x)/2,
1341                                                  (b1pos.y + b2pos.y)/2,
1342                                                  alpha - Math.PI / 2,
1343                                                  coords,
1344                                                  centers,
1345                                                  angles,
1346                                                  straightBulges);
1347                                 if (reverse) {
1348                                         Point2D.Double symAxisVect = new Point2D.Double(coords[b2].x - coords[b1].x, coords[b2].y - coords[b1].y);
1349                                         symmetric(coords[b1], symAxisVect, coords, b1, b2);
1350                                         symmetric(coords[b1], symAxisVect, centers, b1, b2);
1351                                 }
1352                         }
1353                 }
1354         }
1355
1356
1357         /**
1358          * Compute the symmetric of all the points from first to last in the points array
1359          * relative to the line that goes through p and has director vector v.
1360          * The array is modified in-place.
1361          */
1362         private static void symmetric(
1363                         Point2D.Double p,
1364                         Point2D.Double v,
1365                         Point2D.Double[] points,
1366                         int first,
1367                         int last) {
1368                 // ||v||^2
1369                 double lv = v.x*v.x + v.y*v.y;
1370                 for (int i=first; i<=last; i++) {
1371                         // A is the coordinates of points[i] after moving the origin at p
1372                         Point2D.Double A = new Point2D.Double(points[i].x - p.x, points[i].y - p.y);
1373                         // Symmetric of A
1374                         Point2D.Double B = new Point2D.Double(
1375                                         -(A.x*v.y*v.y - 2*A.y*v.x*v.y - A.x*v.x*v.x) / lv,
1376                                          (A.y*v.y*v.y + 2*A.x*v.x*v.y - A.y*v.x*v.x) / lv);
1377                         // Change the origin back
1378                         points[i].x = B.x + p.x;
1379                         points[i].y = B.y + p.y;
1380                 }
1381         }
1382
1383         private void computeHelixTranslations(
1384                         Tree<RNANodeValueTemplate> tree, // IN
1385                         Map<RNATemplateHelix, Point2D.Double> translateVectors, // OUT (must be initialized)
1386                         RNATemplateMapping mapping,      // IN
1387                         Point2D.Double parentDeltaVector, // IN
1388                         boolean straightBulges
1389         ) {
1390                 RNANodeValueTemplate nvt = tree.getValue();
1391                 Point2D.Double newDeltaVector = parentDeltaVector;
1392                 if (nvt instanceof RNANodeValueTemplateBasePair) {
1393                         RNATemplateHelix helix = ((RNANodeValueTemplateBasePair) nvt).getHelix();
1394                         if (! translateVectors.containsKey(helix)) {
1395                                 translateVectors.put(helix, parentDeltaVector);
1396                                 int[] basesInHelixArray;
1397                                 if (mapping.getAncestor(helix) != null) {
1398                                         basesInHelixArray = RNATemplateAlign.intArrayFromList(mapping.getAncestor(helix));
1399                                         Arrays.sort(basesInHelixArray);
1400                                 } else {
1401                                         basesInHelixArray = new int[0];
1402                                 }
1403                                 Point2D.Double helixDeltaVector = computeLengthIncreaseDelta(basesInHelixArray, helix, straightBulges);
1404                                 newDeltaVector = new Point2D.Double(parentDeltaVector.x+helixDeltaVector.x, parentDeltaVector.y+helixDeltaVector.y);
1405                         }
1406                 }
1407                 for (Tree<RNANodeValueTemplate> subtree: tree.getChildren()) {
1408                         computeHelixTranslations(subtree, translateVectors, mapping, newDeltaVector, straightBulges);
1409                 }
1410         }
1411
1412         private Map<RNATemplateHelix, Point2D.Double> computeHelixTranslations(
1413                         Tree<RNANodeValueTemplate> tree, // IN
1414                         RNATemplateMapping mapping,       // IN
1415                         boolean straightBulges
1416         ) {
1417                 Map<RNATemplateHelix, Point2D.Double> translateVectors = new HashMap<RNATemplateHelix, Point2D.Double>();
1418                 computeHelixTranslations(tree, translateVectors, mapping, new Point2D.Double(0,0), straightBulges);
1419                 return translateVectors;
1420         }
1421 }