Phase 8b: flow_tree(optimize = TRUE) refinement

NEWS.md

@@ -10,7 +10,7 @@
 * `edge_bundle_hammer()` is now a native C++ implementation of KDE edge bundling (Hurter et al. 2012) and no longer depends on Python/`reticulate`/datashader. `reticulate` was dropped from Imports and `install_bundle_py()` was removed. Bundling output differs from the datashader-based version
 * added `edge_bundle_mingle()`, multilevel agglomerative edge bundling (Gansner et al. 2011); the kNN proximity graph is built with a bundled kd-tree (nanoflann) for O(E log E) scaling, with a `k` parameter for the number of merge candidates per edge
 * `metro_multicriteria()` is deprecated in favour of `graphlayouts::layout_as_metromap()` and will be removed in a future release
-* added `flow_tree()`, a spiral-tree flow map (Verbeek, Buchin & Speckmann 2011): planar, angle-restricted, keeps node positions fixed, and needs no dummy nodes or triangulation. It is the recommended flow map layout; the `tnss_*()` functions remain as an alternative
+* added `flow_tree()`, a spiral-tree flow map (Verbeek, Buchin & Speckmann 2011): planar, angle-restricted, keeps node positions fixed, and needs no dummy nodes or triangulation. It is the recommended flow map layout; the `tnss_*()` functions remain as an alternative. `flow_tree(optimize = TRUE)` additionally refines the tree (smoothness + obstacle avoidance) with join points held fixed
 * `tnss_tree()`: `order = "weight"` now orders leaves by flow magnitude (previously a no-op); fixed a bug in `tnss_dummies()` where one diagonal corner was misplaced; removed dead code. `interp` moved from Imports to Suggests (only `tnss_tree()` needs it)
 
 # edgebundle 0.4.2

R/flow_spiral.R

@@ -11,6 +11,10 @@
 #' @param alpha restricting angle in degrees (0-90). Smaller values bundle more
 #' tightly toward the root; typical values are 20-45.
 #' @param n number of points sampled per tree edge
+#' @param optimize logical. If `TRUE`, refine the tree with an approximate
+#' FLOWTREE optimization (Verbeek et al. 2011, section 5): join points stay
+#' fixed while edge interiors are relaxed for smoothness and to keep clear of
+#' node obstacles. This can slightly relax the strict `alpha` bound.
 #' @return data.frame with columns `x`, `y`, `flow` and `edge` (one smooth arc
 #' per tree edge)
 #' @author David Schoch
@@ -21,7 +25,7 @@
 #' xy <- cbind(state.center$x, state.center$y)[!state.name %in% c("Alaska", "Hawaii"), ]
 #' flow <- flow_tree(cali2010, xy, root = 4, alpha = 40)
 #' @export
-flow_tree <- function(object, xy, root, alpha = 40, n = 20) {
+flow_tree <- function(object, xy, root, alpha = 40, n = 20, optimize = FALSE) {
     g <- .as_igraph(object)
     if (is.null(g)) {
         stop("flow_tree requires an `igraph` or `tbl_graph` object")
@@ -47,7 +51,11 @@ flow_tree <- function(object, xy, root, alpha = 40, n = 20) {
     }
 
     st <- .spiral_build(xy, root, leaves, wvec[leaves], alpha)
-    .spiral_render(st, n)
+    if (optimize) {
+        .spiral_optimize(st, n)
+    } else {
+        .spiral_render(st, n)
+    }
 }
 
 # Greedy angle-restricted spiral tree. Returns a mutable-env tree description.
@@ -126,28 +134,121 @@ flow_tree <- function(object, xy, root, alpha = 40, n = 20) {
     list(env = e, edges = edges, rootxy = rootxy)
 }
 
-# Render each tree edge as a logarithmic-spiral arc (straight to the root).
+# Sample one tree edge (child -> parent) as n points along its log-spiral arc
+# (straight to the root). Returns an n x 2 matrix.
+.spiral_arc <- function(st, ed, n) {
+    e <- st$env
+    child <- ed[1]
+    parent <- ed[2]
+    u <- seq(0, 1, length.out = n)
+    if (parent == -1L) {
+        xs <- e$x[child] + u * (st$rootxy[1] - e$x[child])
+        ys <- e$y[child] + u * (st$rootxy[2] - e$y[child])
+    } else {
+        dphi <- ((e$phi[parent] - e$phi[child] + pi) %% (2 * pi)) - pi
+        lnt <- log(e$t[child]) + u * (log(e$t[parent]) - log(e$t[child]))
+        ph <- e$phi[child] + u * dphi
+        xs <- st$rootxy[1] + exp(lnt) * cos(ph)
+        ys <- st$rootxy[2] + exp(lnt) * sin(ph)
+    }
+    cbind(xs, ys)
+}
+
+# Render each tree edge as a logarithmic-spiral arc.
 .spiral_render <- function(st, n) {
     e <- st$env
     res <- vector("list", length(st$edges))
+    for (k in seq_along(st$edges)) {
+        arc <- .spiral_arc(st, st$edges[[k]], n)
+        res[[k]] <- data.frame(x = arc[, 1], y = arc[, 2], flow = e$w[st$edges[[k]][1]], edge = k)
+    }
+    do.call(rbind, res)
+}
+
+# Approximate FLOWTREE refinement (Verbeek et al. 2011, section 5): tree nodes
+# (leaves, root, join points) are held fixed so the planar skeleton is
+# preserved; only the interior points of each edge are relaxed to (a) be smooth
+# (FS) and (b) keep clear of node obstacles (Fobs), anchored to the spiral arc
+# so the angle restriction is not badly violated.
+.spiral_optimize <- function(st, n, iterations = 60, buffer = NULL,
+                             w_smooth = 1, w_obs = 0.6, w_anchor = 0.15, lr = 0.25) {
+    e <- st$env
+    ntn <- length(e$t)
+    rootid <- ntn + 1L
+    X <- c(e$x, st$rootxy[1])
+    Y <- c(e$y, st$rootxy[2])
+    fixed <- rep(TRUE, ntn + 1L)
+
+    obs_ids <- c(which(e$isleaf), rootid)
+    ox <- X[obs_ids]
+    oy <- Y[obs_ids]
+
+    chains <- vector("list", length(st$edges))
+    own_ends <- vector("list", 0)
     for (k in seq_along(st$edges)) {
         ed <- st$edges[[k]]
-        child <- ed[1]
-        parent <- ed[2]
-        u <- seq(0, 1, length.out = n)
-        if (parent == -1L) {
-            xs <- e$x[child] + u * (st$rootxy[1] - e$x[child])
-            ys <- e$y[child] + u * (st$rootxy[2] - e$y[child])
+        parent <- if (ed[2] == -1L) rootid else ed[2]
+        arc <- .spiral_arc(st, ed, n)
+        if (n > 2) {
+            ids <- length(X) + seq_len(n - 2)
+            X <- c(X, arc[2:(n - 1), 1])
+            Y <- c(Y, arc[2:(n - 1), 2])
+            fixed <- c(fixed, rep(FALSE, n - 2))
+            for (id in ids) own_ends[[id]] <- c(ed[1], parent)
+            chains[[k]] <- c(ed[1], ids, parent)
         } else {
-            tc <- e$t[child]
-            tp <- e$t[parent]
-            dphi <- ((e$phi[parent] - e$phi[child] + pi) %% (2 * pi)) - pi
-            lnt <- log(tc) + u * (log(tp) - log(tc))
-            ph <- e$phi[child] + u * dphi
-            xs <- st$rootxy[1] + exp(lnt) * cos(ph)
-            ys <- st$rootxy[2] + exp(lnt) * sin(ph)
+            chains[[k]] <- c(ed[1], parent)
+        }
+    }
+    anchorX <- X
+    anchorY <- Y
+    Np <- length(X)
+    nbr <- vector("list", Np)
+    for (ch in chains) {
+        for (i in seq_along(ch)) {
+            if (i > 1) nbr[[ch[i]]] <- c(nbr[[ch[i]]], ch[i - 1])
+            if (i < length(ch)) nbr[[ch[i]]] <- c(nbr[[ch[i]]], ch[i + 1])
         }
-        res[[k]] <- data.frame(x = xs, y = ys, flow = e$w[child], edge = k)
+    }
+
+    if (is.null(buffer)) {
+        lx <- e$x[e$isleaf]
+        ly <- e$y[e$isleaf]
+        nn <- vapply(seq_along(lx), function(i) {
+            d <- sqrt((lx[i] - lx[-i])^2 + (ly[i] - ly[-i])^2)
+            if (length(d)) min(d) else 0
+        }, numeric(1))
+        buffer <- 0.5 * stats::median(nn)
+    }
+
+    movable <- which(!fixed)
+    for (it in seq_len(iterations)) {
+        nX <- X
+        nY <- Y
+        for (p in movable) {
+            nb <- nbr[[p]]
+            mvx <- w_smooth * (mean(X[nb]) - X[p]) + w_anchor * (anchorX[p] - X[p])
+            mvy <- w_smooth * (mean(Y[nb]) - Y[p]) + w_anchor * (anchorY[p] - Y[p])
+            dx <- X[p] - ox
+            dy <- Y[p] - oy
+            d <- sqrt(dx * dx + dy * dy)
+            keep <- d < buffer & d > 1e-9 & !(obs_ids %in% own_ends[[p]])
+            if (any(keep)) {
+                f <- w_obs * (buffer - d[keep]) / buffer
+                mvx <- mvx + sum(f * dx[keep] / d[keep])
+                mvy <- mvy + sum(f * dy[keep] / d[keep])
+            }
+            nX[p] <- X[p] + lr * mvx
+            nY[p] <- Y[p] + lr * mvy
+        }
+        X <- nX
+        Y <- nY
+    }
+
+    res <- vector("list", length(chains))
+    for (k in seq_along(chains)) {
+        ch <- chains[[k]]
+        res[[k]] <- data.frame(x = X[ch], y = Y[ch], flow = e$w[st$edges[[k]][1]], edge = k)
     }
     do.call(rbind, res)
 }

tests/testthat/test-flow_spiral.R

@@ -86,6 +86,24 @@ test_that("flow_tree keeps leaf endpoints fixed at their coordinates", {
     expect_true(all(hit))
 })
 
+test_that("flow_tree(optimize = TRUE) refines while keeping the skeleton", {
+    xy <- cali_xy()
+    base <- flow_tree(cali2010, xy, root = 4, alpha = 40, n = 20)
+    opt <- flow_tree(cali2010, xy, root = 4, alpha = 40, n = 20, optimize = TRUE)
+    expect_named(opt, c("x", "y", "flow", "edge"))
+    expect_equal(sort(unique(opt$edge)), sort(unique(base$edge)))
+    expect_true(all(is.finite(opt$x)) && all(is.finite(opt$y)))
+    # join points and leaves are fixed: edge endpoints must match the skeleton
+    for (k in unique(base$edge)) {
+        b <- base[base$edge == k, ]
+        o <- opt[opt$edge == k, ]
+        expect_equal(unname(as.matrix(b[c(1, nrow(b)), c("x", "y")])),
+                     unname(as.matrix(o[c(1, nrow(o)), c("x", "y")])))
+    }
+    # refinement must not introduce crossings
+    expect_equal(count_crossings(opt), 0)
+})
+
 test_that("flow_tree requires a graph object", {
     expect_error(flow_tree(list(), cbind(0, 0), root = 1), "requires an .igraph")
 })