1
#include <Rcpp.h>
2
#include <deque>
3
using namespace Rcpp;
4

5
// NOTE: This code was copied from the source code of the ggrepl package with 
6
//       the kind permission of its author Kamil Slowikowski.
7
//       ggrepel URL: https://github.com/slowkow/ggrepel
8

9
// Exported convenience functions ---------------------------------------------
10

11
//' Euclidean distance between two points.
12
//' @param a A numeric vector.
13
//' @param b A numeric vector.
14
//' @return The distance between two points.
15
//' @noRd
16
// [[Rcpp::export]]
17 0
double euclid(NumericVector a, NumericVector b) {
18 0
  return sqrt(
19 0
    (a[0] - b[0]) * (a[0] - b[0]) +
20 0
      (a[1] - b[1]) * (a[1] - b[1])
21
  );
22
}
23

24
//' Get the coordinates of the center of a box.
25
//' @param b A box like \code{c(x1, y1, x2, y2)}
26
//' @noRd
27
// [[Rcpp::export]]
28 0
NumericVector centroid(NumericVector b) {
29 0
  return NumericVector::create((b[0] + b[2]) / 2, (b[1] + b[3]) / 2);
30
}
31

32
//' Find the intersections between a line and a rectangle.
33
//' @param p1 A point like \code{c(x, y)}
34
//' @param p2 A point like \code{c(x, y)}
35
//' @param b A rectangle like \code{c(x1, y1, x2, y2)}
36
//' @noRd
37
// [[Rcpp::export]]
38 0
NumericVector intersect_line_rectangle(
39
    NumericVector p1, NumericVector p2, NumericVector b
40
) {
41
  // Sorry for the ugly code :(
42 0
  double dy = p2[1] - p1[1];
43 0
  double dx = p2[0] - p1[0];
44

45 0
  double slope     = dy / dx;
46 0
  double intercept = p2[1] - p2[0] * slope;
47

48 0
  NumericMatrix retval(4, 2);
49 0
  std::fill(retval.begin(), retval.end(), -INFINITY);
50

51
  double x, y;
52

53
  //   +----------+ < b[3]
54
  //   |          |
55
  //   |          | < y
56
  //   |          |
57
  //   +----------+ < b[1]
58
  //   ^    ^     ^
59
  //  b[0]  x    b[2]
60

61 0
  if (dx != 0) {
62
    // Left boundary
63 0
    x = b[0];
64 0
    y = dy == 0 ? p1[1] : slope * x + intercept;
65 0
    if (b[1] <= y && y <= b[3]) {
66 0
      retval(0, _) = NumericVector::create(x, y);
67
    }
68

69
    // Right boundary
70 0
    x = b[2];
71 0
    y = dy == 0 ? p1[1] : slope * x + intercept;
72 0
    if (b[1] <= y && y <= b[3]) {
73 0
      retval(1, _) = NumericVector::create(x, y);
74
    }
75
  }
76

77 0
  if (dy != 0) {
78
      // Bottom boundary
79 0
      y = b[1];
80 0
      x = dx == 0 ? p1[0] : (y - intercept) / slope;
81 0
      if (b[0] <= x && x <= b[2]) {
82 0
        retval(2, _) = NumericVector::create(x, y);
83
      }
84

85
      // Top boundary
86 0
      y = b[3];
87 0
      x = dx == 0 ? p1[0] : (y - intercept) / slope;
88 0
      if (b[0] <= x && x <= b[2]) {
89 0
        retval(3, _) = NumericVector::create(x, y);
90
      }
91
  }
92

93 0
  int i = 0;
94 0
  int imin = 0;
95
  double d;
96 0
  double dmin = INFINITY;
97 0
  for (i = 0; i < 4; i++) {
98 0
    d = euclid(retval(i, _), p1);
99
    // Rcout << i << " euclid = " << d << std::endl;
100 0
    if (d < dmin) {
101 0
      dmin = d;
102 0
      imin = i;
103
    }
104
  }
105

106 0
  return retval(imin, _);
107
}
108

109
// Main code for text label placement -----------------------------------------
110

111
typedef struct {
112
  double x, y;
113
} Point;
114

115 6
Point operator -(const Point& a, const Point& b) {
116 6
  Point p = {a.x - b.x, a.y - b.y};
117 6
  return p;
118
}
119

120 6
Point operator +(const Point& a, const Point& b) {
121 6
  Point p = {a.x + b.x, a.y + b.y};
122 6
  return p;
123
}
124

125 6
Point operator /(const Point& a, const double& b) {
126 6
  Point p = {a.x / b, a.y / b};
127 6
  return p;
128
}
129

130 6
Point operator *(const double& b, const Point& a) {
131 6
  Point p = {a.x * b, a.y * b};
132 6
  return p;
133
}
134

135 6
Point operator *(const Point& a, const double& b) {
136 6
  Point p = {a.x * b, a.y * b};
137 6
  return p;
138
}
139

140
typedef struct {
141
  double x1, y1, x2, y2;
142
} Box;
143

144 6
Box operator +(const Box& b, const Point& p) {
145 6
  Box c = {b.x1 + p.x, b.y1 + p.y, b.x2 + p.x, b.y2 + p.y};
146 6
  return c;
147
}
148

149
//' Euclidean distance between two points.
150
//' @param a A point.
151
//' @param b A point.
152
//' @return The distance between two points.
153
//' @noRd
154 0
double euclid(Point a, Point b) {
155 0
  Point dist = a - b;
156 0
  return sqrt(dist.x * dist.x + dist.y * dist.y);
157
}
158

159
//' Squared Euclidean distance between two points.
160
//' @param a A point.
161
//' @param b A point.
162
//' @return The distance between two points.
163
//' @noRd
164 0
double euclid2(Point a, Point b) {
165 0
  Point dist = a - b;
166 0
  return dist.x * dist.x + dist.y * dist.y;
167
}
168

169
//' Move a box into the area specificied by x limits and y limits.
170
//' @param b A box like \code{c(x1, y1, x2, y2)}
171
//' @param xlim A Point with limits on the x axis like \code{c(xmin, xmax)}
172
//' @param ylim A Point with limits on the y axis like \code{c(xmin, xmax)}
173
//' @param force Magnitude of the force (defaults to \code{1e-6})
174
//' @noRd
175 6
Box put_within_bounds(Box b, Point xlim, Point ylim, double force = 1e-5) {
176
  //double d;
177
  //if (b.x1 < xlim.x) {
178
  //  d = std::max(fabs(b.x1 - xlim.x), 0.02);
179
  //  b.x1 += force / pow(d, 2);
180
  //  b.x2 += force / pow(d, 2);
181
  //} else if (b.x2 > xlim.y) {
182
  //  d = std::max(fabs(b.x2 - xlim.y), 0.02);
183
  //  b.x1 -= force / pow(d, 2);
184
  //  b.x2 -= force / pow(d, 2);
185
  //}
186
  //if (b.y1 < ylim.x) {
187
  //  d = std::max(fabs(b.y1 - ylim.x), 0.02);
188
  //  b.y1 += force / pow(d, 2);
189
  //  b.y2 += force / pow(d, 2);
190
  //} else if (b.y2 > ylim.y) {
191
  //  d = std::max(fabs(b.y2 - ylim.y), 0.02);
192
  //  b.y1 -= force / pow(d, 2);
193
  //  b.y2 -= force / pow(d, 2);
194
  //}
195 6
  double width = fabs(b.x1 - b.x2);
196 6
  double height = fabs(b.y1 - b.y2);
197 6
  if (b.x1 < xlim.x) {
198 6
    b.x1 = xlim.x;
199 6
    b.x2 = b.x1 + width;
200 6
  } else if (b.x2 > xlim.y) {
201 6
    b.x2 = xlim.y;
202 6
    b.x1 = b.x2 - width;
203
  }
204 6
  if (b.y1 < ylim.x) {
205 0
    b.y1 = ylim.x;
206 0
    b.y2 = b.y1 + height;
207 6
  } else if (b.y2 > ylim.y) {
208 4
    b.y2 = ylim.y;
209 4
    b.y1 = b.y2 - height;
210
  }
211 6
  return b;
212
}
213

214
//' Get the coordinates of the center of a box.
215
//' @param b A box like \code{c(x1, y1, x2, y2)}
216
//' @noRd
217 6
Point centroid(Box b) {
218 6
  Point p = {(b.x1 + b.x2) / 2, (b.y1 + b.y2) / 2};
219 6
  return p;
220
}
221

222
//' Test if a box overlaps another box.
223
//' @param a A box like \code{c(x1, y1, x2, y2)}
224
//' @param b A box like \code{c(x1, y1, x2, y2)}
225
//' @noRd
226 6
bool overlaps(Box a, Box b) {
227
  return
228 6
    b.x1 <= a.x2 &&
229 6
    b.y1 <= a.y2 &&
230 6
    b.x2 >= a.x1 &&
231 6
    b.y2 >= a.y1;
232
}
233

234

235

236 6
Point repel_force_both(
237
    Point a, Point b, double force = 0.000001
238
) {
239 6
  double dx = fabs(a.x - b.x);
240 6
  double dy = fabs(a.y - b.y);
241
  // Constrain the minimum distance, so it is never 0.
242 6
  double d2 = std::max(dx * dx + dy * dy, 0.0004);
243
  // Compute a unit vector in the direction of the force.
244 6
  Point v = (a - b) / sqrt(d2);
245
  // Divide the force by the squared distance.
246 6
  Point f = force * v / d2;
247 6
  if (dx > dy) {
248
    // f.y = f.y * dx / dy;
249 6
    f.y = f.y * 2;
250
  } else {
251
    // f.x = f.x * dy / dx;
252 0
    f.x = f.x * 2;
253
  }
254 6
  return f;
255
}
256

257

258 0
Point repel_force_y(
259
    Point a, Point b, double force = 0.000001
260
) {
261 0
  double dx = fabs(a.x - b.x);
262 0
  double dy = fabs(a.y - b.y);
263
  // Constrain the minimum distance, so it is never 0.
264 0
  double d2 = std::max(dx * dx + dy * dy, 0.0004);
265
  // Compute a unit vector in the direction of the force.
266 0
  Point v = {0, (a.y - b.y) / sqrt(d2)};
267
  // Divide the force by the squared distance.
268 0
  Point f = force * v / d2 * 2;
269 0
  return f;
270
}
271

272 0
Point repel_force_x(
273
    Point a, Point b, double force = 0.000001
274
) {
275 0
  double dx = fabs(a.x - b.x);
276 0
  double dy = fabs(a.y - b.y);
277
  // Constrain the minimum distance, so it is never 0.
278 0
  double d2 = std::max(dx * dx + dy * dy, 0.0004);
279
  // Compute a unit vector in the direction of the force.
280 0
  Point v = {(a.x - b.x) / sqrt(d2), 0};
281
  // Divide the force by the squared distance.
282 0
  Point f = force * v / d2 * 2;
283 0
  return f;
284
}
285

286
//' Compute the repulsion force upon point \code{a} from point \code{b}.
287
//'
288
//' The force decays with the squared distance between the points, similar
289
//' to the force of repulsion between magnets.
290
//'
291
//' @param a A point like \code{c(x, y)}
292
//' @param b A point like \code{c(x, y)}
293
//' @param force Magnitude of the force (defaults to \code{1e-6})
294
//' @param direction direction in which to exert force, either "both", "x", or "y"
295
//' @noRd
296 6
Point repel_force(
297
    Point a, Point b, double force = 0.000001, std::string direction = "both"
298
) {
299
  Point out;
300 6
  if (direction == "x"){
301 0
    out = repel_force_x(a, b, force);
302 6
  } else if (direction == "y"){
303 0
    out = repel_force_y(a, b, force);
304
  } else{
305 6
    out = repel_force_both(a, b, force);
306
  }
307 6
  return out;
308
}
309

310

311

312 6
Point spring_force_both(
313
    Point a, Point b, double force = 0.000001
314
) {
315 6
  double dx = fabs(a.x - b.x);
316 6
  double dy = fabs(a.y - b.y);
317 6
  double d = sqrt(dx * dx + dy * dy);
318 6
  Point f = {0, 0};
319 6
  if (d > 0.02) {
320
    // Compute a unit vector in the direction of the force.
321 6
    Point v = (a - b) / d;
322 6
    f = force * v * d;
323 6
    if (dx < dy) {
324 0
      f.y = f.y * 1.5;
325 0
      f.x = f.x * 0.5;
326
    } else {
327 6
      f.y = f.y * 0.5;
328 6
      f.x = f.x * 1.5;
329
    }
330
  }
331 6
  return f;
332
}
333

334 0
Point spring_force_y(
335
    Point a, Point b, double force = 0.000001
336
) {
337 0
  double dx = fabs(a.x - b.x);
338 0
  double dy = fabs(a.y - b.y);
339 0
  double d = sqrt(dx * dx + dy * dy);
340 0
  Point f = {0, 0};
341 0
  if (d > 0.02) {
342
    // Compute a unit vector in the direction of the force.
343 0
    Point v = {0, (a.y - b.y) / d};
344 0
    f = force * v * d;
345 0
    f.y = f.y * 1.5;
346
  }
347 0
  return f;
348
}
349

350 0
Point spring_force_x(
351
    Point a, Point b, double force = 0.000001
352
) {
353 0
  double dx = fabs(a.x - b.x);
354 0
  double dy = fabs(a.y - b.y);
355 0
  double d = sqrt(dx * dx + dy * dy);
356 0
  Point f = {0, 0};
357 0
  if (d > 0.02) {
358
    // Compute a unit vector in the direction of the force.
359 0
    Point v = {(a.x - b.x) / d, 0};
360 0
    f = force * v * d;
361 0
    f.x = f.x * 1.5;
362
  }
363 0
  return f;
364
}
365

366
//' Compute the spring force upon point \code{a} from point \code{b}.
367
//'
368
//' The force increases with the distance between the points, similar
369
//' to Hooke's law for springs.
370
//'
371
//' @param a A point like \code{c(x, y)}
372
//' @param b A point like \code{c(x, y)}
373
//' @param force Magnitude of the force (defaults to \code{1e-6})
374
//' @param direction direction in which to exert force, either "both", "x", or "y"
375
//' @noRd
376 6
Point spring_force(
377
    Point a, Point b, double force = 0.000001, std::string direction = "both"
378
) {
379
  Point out;
380 6
  if (direction == "x"){
381 0
    out = spring_force_x(a, b, force);
382 6
  } else if (direction == "y"){
383 0
    out = spring_force_y(a, b, force);
384
  } else{
385 6
    out = spring_force_both(a, b, force);
386
  }
387 6
  return out;
388
}
389

390
//' Adjust the layout of a list of potentially overlapping boxes.
391
//' @param data_points A numeric matrix with rows representing points like
392
//'   \code{rbind(c(x, y), c(x, y), ...)}
393
//' @param point_padding_x Padding around each data point on the x axis.
394
//' @param point_padding_y Padding around each data point on the y axis.
395
//' @param boxes A numeric matrix with rows representing boxes like
396
//'   \code{rbind(c(x1, y1, x2, y2), c(x1, y1, x2, y2), ...)}
397
//' @param xlim A numeric vector representing the limits on the x axis like
398
//'   \code{c(xmin, xmax)}
399
//' @param ylim A numeric vector representing the limits on the y axis like
400
//'   \code{c(ymin, ymax)}
401
//' @param force Magnitude of the force (defaults to \code{1e-6})
402
//' @param maxiter Maximum number of iterations to try to resolve overlaps
403
//'   (defaults to 2000)
404
//' @noRd
405
// [[Rcpp::export]]
406 6
DataFrame repel_boxes(
407
    NumericMatrix data_points,
408
    double point_padding_x, double point_padding_y,
409
    NumericMatrix boxes,
410
    NumericVector xlim, NumericVector ylim,
411
    double force = 1e-6, int maxiter = 2000,
412
    std::string direction = "both"
413
) {
414 6
  int n_points = data_points.nrow();
415 6
  int n_texts = boxes.nrow();
416
  // assert(n_points >= n_texts);
417 6
  int iter = 0;
418 6
  bool any_overlaps = true;
419

420 6
  if (NumericVector::is_na(force)) {
421 0
    force = 1e-6;
422
  }
423

424
  Point xbounds, ybounds;
425 6
  xbounds.x = xlim[0];
426 6
  xbounds.y = xlim[1];
427 6
  ybounds.x = ylim[0];
428 6
  ybounds.y = ylim[1];
429

430
  // Each data point gets a bounding box.
431 6
  std::vector<Box> DataBoxes(n_points);
432 6
  for (int i = 0; i < n_points; i++) {
433 6
    DataBoxes[i].x1 = data_points(i, 0) - point_padding_x;
434 6
    DataBoxes[i].y1 = data_points(i, 1) - point_padding_y;
435 6
    DataBoxes[i].x2 = data_points(i, 0) + point_padding_x;
436 6
    DataBoxes[i].y2 = data_points(i, 1) + point_padding_y;
437
  }
438

439 6
  std::vector<Point> Points(n_points);
440 6
  for (int i = 0; i < n_points; i++) {
441 6
    Points[i].x = data_points(i, 0);
442 6
    Points[i].y = data_points(i, 1);
443
  }
444

445
  // Add a tiny bit of jitter to each text box at the start.
446 6
  NumericVector r = rnorm(n_texts, 0, force);
447 6
  std::vector<Box> TextBoxes(n_texts);
448 6
  std::vector<double> ratios(n_texts);
449 6
  std::vector<Point> original_centroids(n_texts);
450 6
  for (int i = 0; i < n_texts; i++) {
451 6
    TextBoxes[i].x1 = boxes(i, 0);
452 6
    TextBoxes[i].x2 = boxes(i, 2);
453 6
    TextBoxes[i].y1 = boxes(i, 1);
454 6
    TextBoxes[i].y2 = boxes(i, 3);
455
    // Don't add jitter if the user wants to repel in just one direction.
456 6
    if (direction != "y") {
457 6
      TextBoxes[i].x1 += r[i];
458 6
      TextBoxes[i].x2 += r[i];
459
    }
460 6
    if (direction != "x") {
461 6
      TextBoxes[i].y1 += r[i];
462 6
      TextBoxes[i].y2 += r[i];
463
    }
464
    // height over width
465 6
    ratios[i] = (TextBoxes[i].y2 - TextBoxes[i].y1)
466 6
      / (TextBoxes[i].x2 - TextBoxes[i].x1);
467 6
    original_centroids[i] = centroid(TextBoxes[i]);
468
  }
469

470
  Point f, ci, cj;
471

472 6
  while (any_overlaps && iter < maxiter) {
473 6
    iter += 1;
474 6
    any_overlaps = false;
475

476 6
    for (int i = 0; i < n_texts; i++) {
477 6
      f.x = 0;
478 6
      f.y = 0;
479

480 6
      ci = centroid(TextBoxes[i]);
481

482 6
      for (int j = 0; j < n_points; j++) {
483

484 6
        if (i == j) {
485
          // Skip the data points if the padding is 0.
486 6
          if (point_padding_x == 0 && point_padding_y == 0) {
487 6
            continue;
488
          }
489
          // Repel the box from its data point.
490 0
          if (overlaps(DataBoxes[i], TextBoxes[i])) {
491 0
            any_overlaps = true;
492 0
            f = f + repel_force(ci, Points[i], force, direction);
493
          }
494
        } else {
495
          // Repel the box from overlapping boxes.
496 6
          if (j < n_texts && overlaps(TextBoxes[i], TextBoxes[j])) {
497 6
            any_overlaps = true;
498 6
            cj = centroid(TextBoxes[j]);
499 6
            f = f + repel_force(ci, cj, force * 3, direction);
500
          }
501
          // Skip the data points if the padding is 0.
502 6
          if (point_padding_x == 0 && point_padding_y == 0) {
503 6
            continue;
504
          }
505
          // Repel the box from other data points.
506 0
          if (overlaps(DataBoxes[j], TextBoxes[i])) {
507 0
            any_overlaps = true;
508 0
            f = f + repel_force(ci, Points[j], force, direction);
509
          }
510
        }
511
      }
512

513
      // Pull the box toward its original position.
514 6
      if (!any_overlaps) {
515 6
        f = f + spring_force(original_centroids[i], ci, force * 2e3, direction);
516
      }
517

518
      // Dampen the forces.
519 6
      f = f * (1 - 1e-3);
520

521 6
      TextBoxes[i] = TextBoxes[i] + f;
522 6
      TextBoxes[i] = put_within_bounds(TextBoxes[i], xbounds, ybounds);
523
    }
524
  }
525

526 6
  NumericVector xs(n_texts);
527 6
  NumericVector ys(n_texts);
528

529 6
  for (int i = 0; i < n_texts; i++) {
530 6
    xs[i] = (TextBoxes[i].x1 + TextBoxes[i].x2) / 2;
531 6
    ys[i] = (TextBoxes[i].y1 + TextBoxes[i].y2) / 2;
532
  }
533

534
  return Rcpp::DataFrame::create(
535 6
    Rcpp::Named("x") = xs,
536 6
    Rcpp::Named("y") = ys
537
  );
538
}
539

540

Read our documentation on viewing source code .

Loading