diff --git a/DESCRIPTION b/DESCRIPTION
index ba1bfa4..5eadc4b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -24,7 +24,6 @@ Imports:
     wk (>= 0.9),
     sf (>= 1.0),
     geos (>= 0.2.4),
-    sfnetworks (>= 0.6),
     checkmate,
     igraph (>= 2.0.0)
 Suggests:
diff --git a/NEWS.md b/NEWS.md
index 786c787..5ddd9bc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -13,6 +13,12 @@ centerline 0.3.0 (unreleased)
   * `cnt_path()` now reports when start and end points cannot be connected
     through the skeleton graph instead of passing an empty path to GEOS.
 
+### UPDATES
+
+  * Skeleton routing and anchor graphs now use igraph directly; the
+    dependency on sfnetworks has been removed. Public geometry APIs are
+    unchanged.
+
 centerline 0.2.5 (2026-03-12)
 =========================
 
diff --git a/R/cnt_path.R b/R/cnt_path.R
index 13aba71..80c4c3e 100644
--- a/R/cnt_path.R
+++ b/R/cnt_path.R
@@ -7,9 +7,9 @@
 #' \code{start_point} parameters.
 #'
 #' @details
-#' The following function uses the [sfnetworks::st_network_paths()] approach to
-#' connect \code{start_point} with \code{end_point} by using the
-#' \code{skeleton} of a closed polygon as potential routes.
+#' The function connects start and end points with weighted shortest paths
+#' on an undirected skeleton graph using [igraph::shortest_paths()].
+#' The \code{skeleton} of a closed polygon provides the potential routes.
 #'
 #' It is important to note that multiple starting points are permissible,
 #' but there can only be **one ending point**. Should there be two or more
@@ -102,17 +102,11 @@ cnt_path.geos_geometry <- function(skeleton, start_point, end_point) {
   stopifnot(check_points(end_point))
   check_same_class(skeleton, start_point, end_point)
 
-  if (any(get_geom_type(skeleton) == "multilinestring")) {
-    skeleton <- geos::geos_unnest(skeleton, keep_multi = FALSE)
-  }
-
-  # Transform to sf
-  skeleton <- sf::st_as_sf(skeleton)
-  start_point <- sf::st_as_sf(start_point)
-  end_point <- sf::st_as_sf(end_point)
-
-  # Find the paths
-  cnt_path_master(skeleton, start_point, end_point)
+  cnt_path_master(
+    skeleton_geos = skeleton,
+    start_geos = start_point,
+    end_geos = end_point
+  )
 }
 
 #' @export
@@ -123,12 +117,12 @@ cnt_path.sf <- function(skeleton, start_point, end_point) {
   stopifnot(check_points(end_point))
   check_same_class(skeleton, start_point, end_point)
 
-  if (any(get_geom_type(skeleton) == "MULTILINESTRING")) {
-    skeleton <- sf::st_cast(skeleton, "LINESTRING")
-  }
-
   # Find the paths
-  cnt_path_master(skeleton, start_point, end_point) |>
+  cnt_path_master(
+    skeleton_geos = geos::as_geos_geometry(skeleton),
+    start_geos = geos::as_geos_geometry(start_point),
+    end_geos = geos::as_geos_geometry(end_point)
+  ) |>
     sf::st_as_sf() |>
     cbind(sf::st_drop_geometry(start_point))
 }
@@ -140,12 +134,13 @@ cnt_path.sfc <- function(skeleton, start_point, end_point) {
   stopifnot(check_points(start_point))
   stopifnot(check_points(end_point))
 
-  if (any(get_geom_type(skeleton) == "MULTILINESTRING")) {
-    skeleton <- sf::st_cast(skeleton, "LINESTRING")
-  }
-
   # Find the paths
-  cnt_path_master(skeleton, start_point, end_point) |> sf::st_as_sfc()
+  cnt_path_master(
+    skeleton_geos = geos::as_geos_geometry(skeleton),
+    start_geos = geos::as_geos_geometry(start_point),
+    end_geos = geos::as_geos_geometry(end_point)
+  ) |>
+    sf::st_as_sfc()
 }
 
 #' @export
@@ -159,56 +154,34 @@ cnt_path.SpatVector <- function(skeleton, start_point, end_point) {
   # Save CRS
   crs <- terra::crs(skeleton)
 
-  # Transform to sf objects
-  skeleton <- sf::st_as_sf(skeleton)
-  start_point <- sf::st_as_sf(start_point)
-  end_point <- sf::st_as_sf(end_point)
-
-  if (any(get_geom_type(skeleton) == "MULTILINESTRING")) {
-    skeleton <- sf::st_cast(skeleton, "LINESTRING")
-  }
+  start_sf <- sf::st_as_sf(start_point)
 
   # Find the paths
-  cnt_path_master(skeleton, start_point, end_point) |>
+  cnt_path_master(
+    skeleton_geos = terra_to_geos(skeleton),
+    start_geos = terra_to_geos(start_point),
+    end_geos = terra_to_geos(end_point)
+  ) |>
     wk::as_wkt() |>
     as.character() |>
     terra::vect(crs = crs) |>
-    cbind(sf::st_drop_geometry(start_point))
+    cbind(sf::st_drop_geometry(start_sf))
 }
 
 
-cnt_path_master <- function(skeleton_sf, start_point_sf, end_point_sf) {
-  # Convert skeleton sf object to sfnetworks
-  pol_network <- sfnetworks::as_sfnetwork(
-    x = skeleton_sf,
-    directed = FALSE,
-    length_as_weight = TRUE,
-    edges_as_lines = TRUE
-  )
-
-  # Convert sfnetworks to igraph
-  df_graph <- igraph::as_data_frame(pol_network)
-  names(df_graph)[3] <- "geometry"
-  df_graph <- df_graph[, c("weight", "geometry")]
-  df_graph$weight <- as.numeric(df_graph$weight)
+cnt_path_master <- function(skeleton_geos, start_geos, end_geos) {
+  g <- skeleton_graph_build(skeleton_geos, crs = wk::wk_crs(skeleton_geos))
 
-  # Find indices of nearest nodes for start ...
-  start_nodes <- sf::st_nearest_feature(start_point_sf, pol_network)
-  # ... and end points
-  end_nodes <- sf::st_nearest_feature(end_point_sf, pol_network)
+  start_nodes <- skeleton_graph_nearest_nodes(g, start_geos)
+  end_nodes <- skeleton_graph_nearest_nodes(g, end_geos)
 
   # Check if there are several end nodes
   stopifnot("Only one end point is allowed" = length(end_nodes) == 1)
 
-  # Find the shortest path between two points
-  paths <- base::suppressWarnings(sfnetworks::st_network_paths(
-    pol_network,
-    from = end_nodes,
-    to = start_nodes,
-    weights = "weight"
-  ))
+  # Route from the single end ID to all start IDs
+  edge_paths <- skeleton_graph_paths(g, from = end_nodes, to = start_nodes)
 
-  failed_start_indices <- which(lengths(paths$edge_paths) == 0L)
+  failed_start_indices <- which(lengths(edge_paths) == 0L)
   if (length(failed_start_indices) > 0L) {
     if (length(start_nodes) > 1L) {
       stop(
@@ -225,17 +198,13 @@ cnt_path_master <- function(skeleton_sf, start_point_sf, end_point_sf) {
     )
   }
 
-  # Convert to GEOS geometries and create a GEOS collection
-  lines_list_geos <- lapply(paths$edge_paths, function(x) {
-    df_graph[x, "geometry"]
-  }) |>
-    lapply(geos::as_geos_geometry) |>
-    lapply(geos::geos_make_collection) |>
-    lapply(geos::geos_line_merge)
+  lines_list_geos <- lapply(edge_paths, function(epath) {
+    skeleton_graph_path_line(g, epath)
+  })
 
   # Check if we need to reverse the lines
-  rev_lines_list <- reverse_lines_if_needed(lines_list_geos, end_point_sf)
+  rev_lines_list <- reverse_lines_if_needed(lines_list_geos, end_geos)
 
-  # Return pathes binded together as GEOS geometry
+  # Return paths binded together as GEOS geometry
   do.call(c, rev_lines_list)
 }
diff --git a/R/cnt_path_guess.R b/R/cnt_path_guess.R
index 8e48989..9f114f2 100644
--- a/R/cnt_path_guess.R
+++ b/R/cnt_path_guess.R
@@ -244,74 +244,44 @@ cnt_path_guess.SpatVector <- function(
 }
 
 cnt_path_guess_master <- function(skeleton_geos) {
-  if (geos::geos_type(skeleton_geos) == "multilinestring") {
-    skeleton_geos <- geos::geos_unnest(skeleton_geos, keep_multi = FALSE)
-  }
+  g <- skeleton_graph_build(skeleton_geos, crs = wk::wk_crs(skeleton_geos))
 
-  skeleton_sf <- sf::st_as_sf(skeleton_geos)
-  # Convert skeleton to sfnetworks
-  pol_network <- sfnetworks::as_sfnetwork(
-    x = skeleton_sf,
-    directed = FALSE,
-    length_as_weight = TRUE,
-    edges_as_lines = TRUE
-  )
-  # Convert sfnetworks to igraph
-  # pol_graph <- igraph::as.igraph(pol_network)
-  df_graph <- igraph::as_data_frame(pol_network)
-  names(df_graph)[3] <- "geometry"
-  df_graph <- df_graph[, c("weight", "geometry")]
-  df_graph$weight <- as.numeric(df_graph$weight)
+  # find_outer_nodes expects unnested linework
+  outer_lines <- skeleton_geos
+  if (any(geos::geos_type(outer_lines) == "multilinestring")) {
+    outer_lines <- geos::geos_unnest(outer_lines, keep_multi = FALSE)
+  }
 
   # Find border points of skeleton
-  closest_points <- find_closest_nodes(
-    pol_network,
-    find_outer_nodes(skeleton_geos)
-  )
+  terminal_ids <- find_closest_nodes(g, find_outer_nodes(outer_lines))
 
-  # Find the most distant point from center
+  # Find the most distant point from center (lowest closeness).
   # It will serve as the end point
-  closest_end_points <- closest_points[which.min(igraph::closeness(
-    pol_network,
-    vid = closest_points
+  closest_end_points <- terminal_ids[which.min(igraph::closeness(
+    g$graph,
+    vids = terminal_ids,
+    weights = g$edges$weight
   ))]
 
-  # Find paths
-  paths <- base::suppressWarnings(sfnetworks::st_network_paths(
-    pol_network,
-    to = closest_points[closest_points != closest_end_points],
+  remaining <- terminal_ids[terminal_ids != closest_end_points]
+  edge_paths <- skeleton_graph_paths(
+    g,
     from = closest_end_points,
-    weights = "weight"
-  ))
-
-  # Paths lengths in counts
-  paths_length <- lengths(paths$edge_paths)
-
-  # Filter non-zero paths
-  paths_length_flag <- paths_length > 1
-  paths_length_nonzero <- paths_length[paths_length_flag]
-  edge_paths_nonzero <- paths$edge_paths[paths_length_flag]
-
-  # Estimate paths lengths
-  edge_paths_vec <- unlist(edge_paths_nonzero, use.names = FALSE)
-  edge_paths_groups <- rep(
-    seq_along(edge_paths_nonzero),
-    times = paths_length_nonzero
+    to = remaining
   )
-  edge_paths_length <- df_graph[edge_paths_vec, "weight"]
 
-  # Sum paths lengths in meters
-  true_paths_igraph <- tapply(edge_paths_length, edge_paths_groups, FUN = sum)
+  # Filter non-zero multi-edge paths
+  paths_length <- lengths(edge_paths)
+  paths_length_flag <- paths_length > 1L
+  edge_paths_nonzero <- edge_paths[paths_length_flag]
 
-  # Return the longest path
-  longest_path_igraph <- which.max(true_paths_igraph)
-  longest_path_geos <- df_graph[
-    edge_paths_nonzero[[longest_path_igraph]],
-    "geometry"
-  ] |>
-    geos::as_geos_geometry() |>
-    geos::geos_make_collection() |>
-    geos::geos_line_merge()
+  true_paths_igraph <- vapply(
+    edge_paths_nonzero,
+    function(epath) sum(g$edges$weight[epath]),
+    numeric(1),
+    USE.NAMES = FALSE
+  )
 
-  longest_path_geos
+  longest_path_igraph <- which.max(true_paths_igraph)
+  skeleton_graph_path_line(g, edge_paths_nonzero[[longest_path_igraph]])
 }
diff --git a/R/cnt_skeleton_anchors.R b/R/cnt_skeleton_anchors.R
index 459df1d..be5c7c2 100644
--- a/R/cnt_skeleton_anchors.R
+++ b/R/cnt_skeleton_anchors.R
@@ -113,66 +113,48 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
   boundary <- geos::geos_boundary(polygon)
 
   # Raw unrounded graph used for connector endpoints and branch coverage.
-  raw_lines <- ordinary_skeleton |>
-    geos::geos_node() |>
-    geos::geos_unnest(keep_multi = FALSE) |>
-    sf::st_as_sf()
-  raw_network <- sfnetworks::as_sfnetwork(
-    raw_lines,
-    directed = FALSE,
-    edges_as_lines = TRUE
+  # First implementation: use the same raw graph for analysis (single-graph).
+  # If fixture-approved junction/scale assertions require dual-graph behavior,
+  # fall back to a precision-cleaned analysis graph below.
+  raw_graph <- skeleton_graph_build(ordinary_skeleton, crs = crs)
+  raw_degree <- skeleton_graph_degree(raw_graph)
+  raw_xy <- cbind(raw_graph$nodes$x, raw_graph$nodes$y)
+  raw_node_pts <- skeleton_graph_node_points(raw_graph)
+
+  # Precision-cleaned analysis graph for candidate discovery/scoring.
+  # Raw graph remains the connector-target and postcondition authority.
+  # Extreme scales fall back to the unrounded node set for analysis only.
+  analysis_graph <- tryCatch(
+    skeleton_graph_build(
+      ordinary_skeleton,
+      crs = crs,
+      grid_size = grid_size
+    ),
+    error = function(e) raw_graph
   )
-  raw_nodes <- sfnetworks::activate(raw_network, "nodes") |> sf::st_as_sf()
-  raw_degree <- igraph::degree(raw_network)
-  raw_xy <- sf::st_coordinates(raw_nodes)
-
-  # Cleaned analysis graph used only for candidate discovery and scoring.
-  cleaned_lines <- tryCatch(
-    {
-      ordinary_skeleton |>
-        geos::geos_set_precision(
-          grid_size,
-          preserve_topology = TRUE,
-          keep_collapsed = FALSE
-        ) |>
-        geos::geos_node() |>
-        geos::geos_unnest(keep_multi = FALSE) |>
-        sf::st_as_sf() |>
-        sf::st_set_precision(0)
-    },
-    error = function(e) {
-      # Extreme scales can make precision-rounded noding fail to converge.
-      # Fall back to the unrounded node set for analysis only; mapping still
-      # uses grid_size against raw junctions.
-      raw_lines
-    }
-  )
-  cleaned_network <- sfnetworks::as_sfnetwork(
-    cleaned_lines,
-    directed = FALSE,
-    edges_as_lines = TRUE
+
+  analysis_degree <- skeleton_graph_degree(analysis_graph)
+  analysis_xy <- cbind(analysis_graph$nodes$x, analysis_graph$nodes$y)
+  analysis_node_pts <- skeleton_graph_node_points(analysis_graph)
+  analysis_boundary_dist <- as.numeric(
+    geos::geos_distance(analysis_node_pts, boundary)
   )
-  cleaned_nodes <- sfnetworks::activate(cleaned_network, "nodes") |>
-    sf::st_as_sf()
-  cleaned_degree <- igraph::degree(cleaned_network)
-  cleaned_xy <- sf::st_coordinates(cleaned_nodes)
-  cleaned_boundary_dist <- as.numeric(sf::st_distance(
-    cleaned_nodes,
-    sf::st_as_sf(boundary)
-  ))
 
-  cleaned_junction_ids <- which(cleaned_degree >= 3L)
+  analysis_junction_ids <- which(analysis_degree >= 3L)
   raw_junction_ids <- which(raw_degree >= 3L)
 
-  if (length(cleaned_junction_ids) == 0L || length(raw_junction_ids) == 0L) {
+  if (length(analysis_junction_ids) == 0L || length(raw_junction_ids) == 0L) {
     stop(
       "No valid interior junction connector; ",
       "ordinary skeleton has no degree-at-least-three junction"
     )
   }
 
-  map_cleaned_to_raw <- function(cleaned_id) {
-    cxy <- cleaned_xy[cleaned_id, 1:2, drop = FALSE]
+  map_analysis_to_raw <- function(analysis_id) {
+    if (identical(analysis_graph, raw_graph)) {
+      return(as.integer(analysis_id))
+    }
+    cxy <- analysis_xy[analysis_id, 1:2, drop = FALSE]
     d <- sqrt(rowSums(
       (raw_xy[raw_junction_ids, 1:2, drop = FALSE] -
         matrix(cxy, nrow = length(raw_junction_ids), ncol = 2, byrow = TRUE))^2
@@ -189,49 +171,49 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
   }
 
   mapped_raw <- vapply(
-    cleaned_junction_ids,
-    map_cleaned_to_raw,
+    analysis_junction_ids,
+    map_analysis_to_raw,
     integer(1),
     USE.NAMES = FALSE
   )
   keep_map <- !is.na(mapped_raw)
-  cleaned_junction_ids <- cleaned_junction_ids[keep_map]
+  analysis_junction_ids <- analysis_junction_ids[keep_map]
   mapped_raw <- mapped_raw[keep_map]
 
-  if (length(cleaned_junction_ids) == 0L) {
+  if (length(analysis_junction_ids) == 0L) {
     stop(
       "No valid interior junction connector; ",
       "no cleaned junction maps to a raw degree-at-least-three node"
     )
   }
 
-  n_junctions <- length(cleaned_junction_ids)
+  n_junctions <- length(analysis_junction_ids)
   # Junction coordinate matrix and raw targets (numeric hot path).
-  j_xy <- cleaned_xy[cleaned_junction_ids, 1:2, drop = FALSE]
+  j_xy <- analysis_xy[analysis_junction_ids, 1:2, drop = FALSE]
   j_raw_id <- as.integer(mapped_raw)
 
-  # Adjacency once; avoids igraph::neighbors per junction per anchor.
-  cleaned_adj <- vector("list", nrow(cleaned_xy))
-  for (cid in cleaned_junction_ids) {
-    cleaned_adj[[cid]] <- as.integer(igraph::neighbors(cleaned_network, cid))
+  # Adjacency once; avoids neighbors() per junction per anchor.
+  analysis_adj <- vector("list", length(analysis_degree))
+  for (cid in analysis_junction_ids) {
+    analysis_adj[[cid]] <- skeleton_graph_neighbors(analysis_graph, cid)
   }
 
   select_inward_neighbor <- function(cleaned_id, anchor_xy) {
-    neighbours <- cleaned_adj[[cleaned_id]]
+    neighbours <- analysis_adj[[cleaned_id]]
     if (length(neighbours) == 0L) {
       return(NA_integer_)
     }
-    dists <- cleaned_boundary_dist[neighbours]
+    dists <- analysis_boundary_dist[neighbours]
     max_d <- max(dists)
     candidates <- neighbours[dists >= (max_d - grid_size)]
     if (length(candidates) == 1L) {
       return(candidates[[1]])
     }
-    jxy <- cleaned_xy[cleaned_id, 1:2]
+    jxy <- analysis_xy[cleaned_id, 1:2]
     angles <- vapply(
       candidates,
       function(nid) {
-        continuation_angle(anchor_xy, jxy, cleaned_xy[nid, 1:2])
+        continuation_angle(anchor_xy, jxy, analysis_xy[nid, 1:2])
       },
       numeric(1),
       USE.NAMES = FALSE
@@ -277,13 +259,13 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
 
     for (i in seq_len(n)) {
       j <- j_idx[[i]]
-      cid <- cleaned_junction_ids[[j]]
+      cid <- analysis_junction_ids[[j]]
       cxy <- j_xy[j, ]
       nid <- select_inward_neighbor(cid, axy)
       if (is.na(nid)) {
         next
       }
-      ang <- continuation_angle(axy, cxy, cleaned_xy[nid, 1:2])
+      ang <- continuation_angle(axy, cxy, analysis_xy[nid, 1:2])
       if (is.na(ang)) {
         next
       }
@@ -363,16 +345,16 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
   }
 
   already_terminal_on_skeleton <- function(anchor) {
-    d <- as.numeric(sf::st_distance(sf::st_as_sf(anchor), raw_nodes))
-    zero_ids <- which(d == 0)
-    if (length(zero_ids) != 1L) {
+    id <- skeleton_graph_nearest_nodes(raw_graph, anchor)
+    if (
+      as.numeric(geos::geos_distance(anchor, raw_node_pts[id])) != 0
+    ) {
       return(FALSE)
     }
-    id <- zero_ids[[1]]
     if (raw_degree[[id]] != 1L) {
       return(FALSE)
     }
-    neighbours <- as.integer(igraph::neighbors(raw_network, id))
+    neighbours <- skeleton_graph_neighbors(raw_graph, id)
     length(neighbours) == 1L && raw_degree[[neighbours[[1]]]] >= 3L
   }
 
@@ -391,8 +373,10 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
 
     # Anchor on a non-terminal skeleton node cannot become a terminal
     # without pruning branches.
-    d_raw <- as.numeric(sf::st_distance(sf::st_as_sf(anchor), raw_nodes))
-    if (any(d_raw == 0)) {
+    raw_near <- skeleton_graph_nearest_nodes(raw_graph, anchor)
+    if (
+      as.numeric(geos::geos_distance(anchor, raw_node_pts[raw_near])) == 0
+    ) {
       stop(
         "No valid interior junction connector for anchor ",
         ai,
@@ -574,22 +558,19 @@ add_boundary_anchors <- function(skeleton, polygon, anchors) {
 
   # Postconditions: anchors are degree-one terminals; ordinary branches remain.
   aug_edges <- augmented |> geos::geos_unnest(keep_multi = FALSE)
-  aug_lines <- sf::st_as_sf(aug_edges)
-  aug_network <- sfnetworks::as_sfnetwork(
-    aug_lines,
-    directed = FALSE,
-    edges_as_lines = TRUE
-  )
-  aug_nodes <- sfnetworks::activate(aug_network, "nodes") |> sf::st_as_sf()
-  aug_degree <- igraph::degree(aug_network)
+  aug_graph <- skeleton_graph_build(augmented, crs = crs)
+  aug_degree <- skeleton_graph_degree(aug_graph)
+  aug_node_pts <- skeleton_graph_node_points(aug_graph)
 
   for (ai in seq_along(anchors)) {
     if (isTRUE(selected_meta[[ai]]$noop)) {
       next
     }
-    d <- as.numeric(sf::st_distance(sf::st_as_sf(anchors[ai]), aug_nodes))
-    zero_ids <- which(d == 0)
-    if (length(zero_ids) != 1L || aug_degree[[zero_ids[[1]]]] != 1L) {
+    zero_id <- skeleton_graph_nearest_nodes(aug_graph, anchors[ai])
+    on_node <- as.numeric(
+      geos::geos_distance(anchors[ai], aug_node_pts[zero_id])
+    ) == 0
+    if (!on_node || aug_degree[[zero_id]] != 1L) {
       stop(
         "Postcondition failed: anchor ",
         ai,
diff --git a/R/skeleton_graph.R b/R/skeleton_graph.R
new file mode 100644
index 0000000..d1ba90b
--- /dev/null
+++ b/R/skeleton_graph.R
@@ -0,0 +1,195 @@
+# Planar line graph from GEOS-noded skeleton segments --------------------
+# Undirected weighted graph with edge geometry and node strtree. Used by
+# cnt_path(), cnt_path_guess(), and boundary-anchor topology queries.
+
+skeleton_graph_build <- function(lines, crs = wk::wk_crs(lines), grid_size = NULL) {
+  if (!inherits(lines, "geos_geometry")) {
+    stop("lines must be a geos_geometry object")
+  }
+
+  line_types <- geos::geos_type(lines)
+  if (
+    !all(
+      line_types %in%
+        c("linestring", "multilinestring", "linearring")
+    )
+  ) {
+    stop("lines must contain only line geometries")
+  }
+
+  if (any(geos::geos_is_empty(lines))) {
+    lines <- lines[!geos::geos_is_empty(lines)]
+  }
+
+  if (length(lines) == 0L) {
+    stop("no non-empty edges after noding")
+  }
+
+  gs <- NA_real_
+  if (!is.null(grid_size)) {
+    if (
+      length(grid_size) != 1L ||
+        !is.finite(grid_size) ||
+        grid_size <= 0
+    ) {
+      stop("grid_size must be one positive finite number")
+    }
+    gs <- as.numeric(grid_size)
+    lines <- geos::geos_set_precision(
+      lines,
+      gs,
+      preserve_topology = TRUE,
+      keep_collapsed = FALSE
+    )
+  }
+
+  # One collection so intersections across separate input features node together.
+  noded <- geos::geos_node(geos::geos_make_collection(lines))
+  edges <- geos::geos_unnest(noded, keep_multi = FALSE)
+  edges <- edges[!geos::geos_is_empty(edges)]
+  edge_types <- geos::geos_type(edges)
+  edges <- edges[edge_types %in% c("linestring", "linearring")]
+
+  if (length(edges) == 0L) {
+    stop("no non-empty edges after noding")
+  }
+
+  weight <- as.numeric(geos::geos_length(edges))
+  if (any(!is.finite(weight) | weight <= 0)) {
+    stop("edge lengths must be positive and finite")
+  }
+
+  start_coords <- wk::wk_coords(geos::geos_point_start(edges))
+  end_coords <- wk::wk_coords(geos::geos_point_end(edges))
+
+  sx <- start_coords$x
+  sy <- start_coords$y
+  ex <- end_coords$x
+  ey <- end_coords$y
+
+  # Interleave start/end in edge order for deterministic first-occurrence IDs.
+  xs <- as.numeric(rbind(sx, ex))
+  ys <- as.numeric(rbind(sy, ey))
+  # Normalize signed zero so %a keys match.
+  xs[xs == 0] <- 0
+  ys[ys == 0] <- 0
+
+  keys_all <- paste(sprintf("%a", xs), sprintf("%a", ys), sep = "|")
+  key_unique <- unique(keys_all)
+  endpoint_ids <- match(keys_all, key_unique)
+
+  from <- as.integer(endpoint_ids[c(TRUE, FALSE)])
+  to <- as.integer(endpoint_ids[c(FALSE, TRUE)])
+
+  n_nodes <- length(key_unique)
+  # First-occurrence coordinates for each unique key.
+  first_idx <- match(key_unique, keys_all)
+  node_x <- xs[first_idx]
+  node_y <- ys[first_idx]
+
+  graph <- igraph::graph_from_edgelist(
+    cbind(from, to),
+    directed = FALSE
+  )
+  # Ensure isolated-node-free graphs still report full vertex count when
+  # edgelist max id equals n_nodes (always true here).
+  if (igraph::vcount(graph) < n_nodes) {
+    graph <- igraph::add_vertices(graph, n_nodes - igraph::vcount(graph))
+  }
+
+  igraph::E(graph)$weight <- weight
+  igraph::E(graph)$edge_id <- seq_along(weight)
+
+  node_geom <- geos::geos_make_point(node_x, node_y, crs = crs)
+  node_tree <- geos::geos_strtree(node_geom)
+
+  if (!is.null(crs)) {
+    wk::wk_crs(edges) <- crs
+  }
+
+  structure(
+    list(
+      nodes = list(
+        x = node_x,
+        y = node_y,
+        key = key_unique,
+        geom = node_geom,
+        tree = node_tree
+      ),
+      edges = list(
+        from = from,
+        to = to,
+        weight = weight,
+        geom = edges
+      ),
+      graph = graph,
+      crs = crs,
+      meta = list(
+        noded = TRUE,
+        grid_size = gs
+      )
+    ),
+    class = "cnt_skeleton_graph"
+  )
+}
+
+skeleton_graph_degree <- function(g) {
+  as.integer(igraph::degree(g$graph))
+}
+
+skeleton_graph_neighbors <- function(g, id) {
+  as.integer(igraph::neighbors(g$graph, id, mode = "all"))
+}
+
+skeleton_graph_nearest_nodes <- function(g, points) {
+  as.integer(geos::geos_nearest(points, g$nodes$tree))
+}
+
+skeleton_graph_paths <- function(g, from, to) {
+  if (length(from) != 1L) {
+    stop("paths() requires exactly one source node")
+  }
+
+  withCallingHandlers(
+    {
+      sp <- igraph::shortest_paths(
+        g$graph,
+        from = from,
+        to = to,
+        mode = "all",
+        weights = g$edges$weight,
+        output = "epath"
+      )
+    },
+    warning = function(w) {
+      if (grepl("Couldn't reach some vertices", conditionMessage(w), fixed = TRUE)) {
+        invokeRestart("muffleWarning")
+      }
+    }
+  )
+
+  lapply(sp$epath, function(e) {
+    if (length(e) == 0L) {
+      integer()
+    } else {
+      as.integer(igraph::E(g$graph)$edge_id[e])
+    }
+  })
+}
+
+skeleton_graph_path_line <- function(g, epath) {
+  if (length(epath) == 0L) {
+    stop("empty edge path")
+  }
+  line <- g$edges$geom[epath] |>
+    geos::geos_make_collection() |>
+    geos::geos_line_merge()
+  if (!is.null(g$crs)) {
+    wk::wk_crs(line) <- g$crs
+  }
+  line
+}
+
+skeleton_graph_node_points <- function(g) {
+  g$nodes$geom
+}
diff --git a/R/utils.R b/R/utils.R
index e0c02cd..44cc9b6 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -93,13 +93,8 @@ find_outer_nodes <- function(skeleton_geos) {
   c(lonely_end, lonely_start)
 }
 
-find_closest_nodes <- function(sf_graph, nodes_geos) {
-  geos_graph <- sfnetworks::activate(sf_graph, "nodes") |>
-    sf::st_as_sf() |>
-    geos::as_geos_geometry() |>
-    geos::geos_strtree()
-
-  geos::geos_nearest(nodes_geos, geos_graph)
+find_closest_nodes <- function(g, nodes_geos) {
+  skeleton_graph_nearest_nodes(g, nodes_geos)
 }
 
 # Straight skeleton helpers ----------------------------------------------
diff --git a/codemeta.json b/codemeta.json
index 523e514..5d8c7ca 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -189,28 +189,28 @@
     },
     "5": {
       "@type": "SoftwareApplication",
-      "identifier": "sfnetworks",
-      "name": "sfnetworks",
-      "version": ">= 0.6",
+      "identifier": "checkmate",
+      "name": "checkmate",
       "provider": {
         "@id": "https://cran.r-project.org",
         "@type": "Organization",
         "name": "Comprehensive R Archive Network (CRAN)",
         "url": "https://cran.r-project.org"
       },
-      "sameAs": "https://CRAN.R-project.org/package=sfnetworks"
+      "sameAs": "https://CRAN.R-project.org/package=checkmate"
     },
     "6": {
       "@type": "SoftwareApplication",
-      "identifier": "checkmate",
-      "name": "checkmate",
+      "identifier": "igraph",
+      "name": "igraph",
+      "version": ">= 2.0.0",
       "provider": {
         "@id": "https://cran.r-project.org",
         "@type": "Organization",
         "name": "Comprehensive R Archive Network (CRAN)",
         "url": "https://cran.r-project.org"
       },
-      "sameAs": "https://CRAN.R-project.org/package=checkmate"
+      "sameAs": "https://CRAN.R-project.org/package=igraph"
     },
     "SystemRequirements": null
   },
diff --git a/man/cnt_path.Rd b/man/cnt_path.Rd
index 7393a24..55539ee 100644
--- a/man/cnt_path.Rd
+++ b/man/cnt_path.Rd
@@ -23,9 +23,9 @@ or \code{geos_geometry} class objects of a \code{LINESTRING} geometry
 Find the shortest path between start and end points within a polygon
 }
 \details{
-The following function uses the \code{\link[sfnetworks:st_network_paths]{sfnetworks::st_network_paths()}} approach to
-connect \code{start_point} with \code{end_point} by using the
-\code{skeleton} of a closed polygon as potential routes.
+The function connects start and end points with weighted shortest paths
+on an undirected skeleton graph using \code{\link[igraph:shortest_paths]{igraph::shortest_paths()}}.
+The \code{skeleton} of a closed polygon provides the potential routes.
 
 It is important to note that multiple starting points are permissible,
 but there can only be \strong{one ending point}. Should there be two or more
diff --git a/tests/testthat/fixtures/graph-migration-baseline.rds b/tests/testthat/fixtures/graph-migration-baseline.rds
new file mode 100644
index 0000000..861d5f4
Binary files /dev/null and b/tests/testthat/fixtures/graph-migration-baseline.rds differ
diff --git a/tests/testthat/test-cnt_path.R b/tests/testthat/test-cnt_path.R
index facc4a6..3c763f6 100644
--- a/tests/testthat/test-cnt_path.R
+++ b/tests/testthat/test-cnt_path.R
@@ -492,3 +492,62 @@ test_that("cnt_path returns exact ordered endpoints for anchored skeletons", {
   )
   expect_equal(as.data.frame(path_terra), as.data.frame(anchors_terra[1]))
 })
+
+
+
+migration_geom_close <- function(actual, baseline_wkb, baseline_length) {
+  tol <- max(1e-3, 1e-6 * baseline_length)
+  actual_len <- as.numeric(geos::geos_length(actual))
+  expect_equal(actual_len, baseline_length, tolerance = tol)
+  baseline <- geos::as_geos_geometry(baseline_wkb)
+  equal_exact <- tryCatch(
+    isTRUE(geos::geos_equals_exact(actual, baseline, tol)),
+    error = function(e) FALSE
+  )
+  if (!equal_exact) {
+    # Allow reverse orientation of the whole path
+    equal_rev <- tryCatch(
+      isTRUE(geos::geos_equals_exact(actual, geos::geos_reverse(baseline), tol)),
+      error = function(e) FALSE
+    )
+    if (!equal_rev) {
+      hd <- as.numeric(geos::geos_distance_hausdorff(actual, baseline))
+      expect_lte(hd, tol)
+    }
+  }
+  invisible(TRUE)
+}
+
+test_that("cnt_path free and anchored routes match graph-migration baseline", {
+  baseline <- readRDS(test_path("fixtures/graph-migration-baseline.rds"))
+  f <- system.file("extdata/example.gpkg", package = "centerline")
+  p <- geos::as_geos_geometry(sf::st_read(f, layer = "polygon", quiet = TRUE))
+  pts <- geos::as_geos_geometry(
+    sf::st_read(f, layer = "polygon_points", quiet = TRUE)
+  )
+  b <- geos::geos_boundary(p)
+  a <- pts[geos::geos_equals(geos::geos_intersection(b, pts), pts)][1:2]
+
+  ordinary <- cnt_skeleton(p, keep = 1)
+  free1 <- cnt_path(ordinary, pts[2], pts[1])
+  free2 <- cnt_path(ordinary, pts[3], pts[1])
+  migration_geom_close(
+    free1,
+    baseline$free_paths_wkb[1],
+    baseline$lengths[["free_path_1"]]
+  )
+  migration_geom_close(
+    free2,
+    baseline$free_paths_wkb[2],
+    baseline$lengths[["free_path_2"]]
+  )
+
+  anchored <- cnt_skeleton(p, keep = 1, anchors = a)
+  apath <- cnt_path(anchored, a[1], a[2])
+  migration_geom_close(
+    apath,
+    baseline$anchored_path_wkb,
+    baseline$lengths[["anchored_path"]]
+  )
+})
+
diff --git a/tests/testthat/test-cnt_path_guess.R b/tests/testthat/test-cnt_path_guess.R
index b28d7e4..9231519 100644
--- a/tests/testthat/test-cnt_path_guess.R
+++ b/tests/testthat/test-cnt_path_guess.R
@@ -411,3 +411,35 @@ test_that("cnt_path_guess forwards anchors only when building a skeleton", {
   expect_s3_class(result_supplied, "sf")
   expect_contains(get_geom_type(result_supplied), "LINESTRING")
 })
+
+
+test_that("cnt_path_guess shapes[89] matches graph-migration baseline", {
+  baseline <- readRDS(test_path("fixtures/graph-migration-baseline.rds"))
+  shapes <- sf::st_read(
+    system.file("extdata/example.gpkg", package = "centerline"),
+    layer = "shapes",
+    quiet = TRUE
+  )
+  gpath <- cnt_path_guess(shapes[89, ], keep = 1, return_geos = TRUE)
+  tol <- max(1e-3, 1e-6 * baseline$lengths[["guessed_path"]])
+  expect_equal(
+    as.numeric(geos::geos_length(gpath)),
+    baseline$lengths[["guessed_path"]],
+    tolerance = tol
+  )
+  base_geom <- geos::as_geos_geometry(baseline$guessed_path_wkb)
+  equal_exact <- tryCatch(
+    isTRUE(geos::geos_equals_exact(gpath, base_geom, tol)),
+    error = function(e) FALSE
+  )
+  if (!equal_exact) {
+    equal_rev <- tryCatch(
+      isTRUE(geos::geos_equals_exact(gpath, geos::geos_reverse(base_geom), tol)),
+      error = function(e) FALSE
+    )
+    if (!equal_rev) {
+      hd <- as.numeric(geos::geos_distance_hausdorff(gpath, base_geom))
+      expect_lte(hd, tol)
+    }
+  }
+})
diff --git a/tests/testthat/test-cnt_skeleton.R b/tests/testthat/test-cnt_skeleton.R
index c92502a..ca15d29 100644
--- a/tests/testthat/test-cnt_skeleton.R
+++ b/tests/testthat/test-cnt_skeleton.R
@@ -286,22 +286,13 @@ anchor_grid_size <- function(polygon) {
   sqrt((bbox$xmax - bbox$xmin)^2 + (bbox$ymax - bbox$ymin)^2) * 1e-7
 }
 
-skeleton_nodes_network <- function(skeleton) {
-  lines <- skeleton |>
-    geos::geos_node() |>
-    geos::geos_unnest(keep_multi = FALSE) |>
-    sf::st_as_sf()
-  net <- sfnetworks::as_sfnetwork(
-    lines,
-    directed = FALSE,
-    edges_as_lines = TRUE
-  )
-  nodes <- sfnetworks::activate(net, "nodes") |> sf::st_as_sf()
+skeleton_graph_info <- function(skeleton) {
+  g <- skeleton_graph_build(skeleton, crs = wk::wk_crs(skeleton))
   list(
-    network = net,
-    nodes = nodes,
-    degree = igraph::degree(net),
-    xy = sf::st_coordinates(nodes)
+    graph = g,
+    nodes = skeleton_graph_node_points(g),
+    degree = skeleton_graph_degree(g),
+    xy = cbind(g$nodes$x, g$nodes$y)
   )
 }
 
@@ -310,22 +301,19 @@ test_that("boundary anchors attach approved junctions on example.gpkg", {
   polygon <- data$polygon
   anchors <- data$anchors
   ordinary <- cnt_skeleton(polygon, keep = 1)
-  ordinary_net <- skeleton_nodes_network(ordinary)
+  ordinary_net <- skeleton_graph_info(ordinary)
   grid_size <- anchor_grid_size(polygon)
 
   # Pre-augmentation: anchors are not ordinary skeleton nodes.
   for (i in seq_along(anchors)) {
-    d <- as.numeric(sf::st_distance(
-      sf::st_as_sf(anchors[i]),
-      ordinary_net$nodes
-    ))
+    d <- as.numeric(geos::geos_distance(anchors[i], ordinary_net$nodes))
     expect_true(all(d > 0))
   }
 
   anchored <- cnt_skeleton(polygon, keep = 1, anchors = anchors)
   expect_true(geos::geos_covers(anchored, ordinary))
 
-  aug <- skeleton_nodes_network(anchored)
+  aug <- skeleton_graph_info(anchored)
   approved <- list(
     c(1830873.1875, 5453788.018333333),
     c(1830873.699999998, 5453778.95)
@@ -334,12 +322,12 @@ test_that("boundary anchors attach approved junctions on example.gpkg", {
   boundary <- geos::geos_boundary(polygon)
 
   for (i in seq_along(anchors)) {
-    d <- as.numeric(sf::st_distance(sf::st_as_sf(anchors[i]), aug$nodes))
+    d <- as.numeric(geos::geos_distance(anchors[i], aug$nodes))
     zero_ids <- which(d == 0)
     expect_equal(length(zero_ids), 1L)
     expect_equal(unname(aug$degree[[zero_ids[[1]]]]), 1L)
 
-    nbrs <- as.integer(igraph::neighbors(aug$network, zero_ids[[1]]))
+    nbrs <- skeleton_graph_neighbors(aug$graph, zero_ids[[1]])
     expect_equal(length(nbrs), 1L)
     expect_true(aug$degree[[nbrs[[1]]]] >= 3L)
     target_xy <- aug$xy[nbrs[[1]], 1:2]
@@ -425,23 +413,20 @@ test_that("boundary anchors scale with relative precision grid", {
   for (scale in c(1e-6, 1e6)) {
     pair <- scale_pair(scale)
     ordinary <- cnt_skeleton(pair$polygon, keep = 1)
-    ordinary_net <- skeleton_nodes_network(ordinary)
+    ordinary_net <- skeleton_graph_info(ordinary)
     for (i in seq_along(pair$anchors)) {
-      d <- as.numeric(sf::st_distance(
-        sf::st_as_sf(pair$anchors[i]),
-        ordinary_net$nodes
-      ))
+      d <- as.numeric(geos::geos_distance(pair$anchors[i], ordinary_net$nodes))
       expect_true(all(d > 0))
     }
     anchored <- cnt_skeleton(pair$polygon, keep = 1, anchors = pair$anchors)
-    aug <- skeleton_nodes_network(anchored)
+    aug <- skeleton_graph_info(anchored)
     targets <- list()
     for (i in seq_along(pair$anchors)) {
-      d <- as.numeric(sf::st_distance(sf::st_as_sf(pair$anchors[i]), aug$nodes))
+      d <- as.numeric(geos::geos_distance(pair$anchors[i], aug$nodes))
       zero_ids <- which(d == 0)
       expect_equal(length(zero_ids), 1L)
       expect_equal(unname(aug$degree[[zero_ids[[1]]]]), 1L)
-      nbrs <- as.integer(igraph::neighbors(aug$network, zero_ids[[1]]))
+      nbrs <- skeleton_graph_neighbors(aug$graph, zero_ids[[1]])
       targets[[i]] <- aug$xy[nbrs[[1]], 1:2]
     }
     # normalized junction choice: same relative offset order
@@ -577,11 +562,11 @@ test_that("add_boundary_anchors rejects concavity-crossing chords", {
     anchor
   ))
 
-  aug <- skeleton_nodes_network(augmented)
-  d <- as.numeric(sf::st_distance(sf::st_as_sf(anchor), aug$nodes))
+  aug <- skeleton_graph_info(augmented)
+  d <- as.numeric(geos::geos_distance(anchor, aug$nodes))
   zero_ids <- which(d == 0)
   expect_equal(length(zero_ids), 1L)
-  nbrs <- as.integer(igraph::neighbors(aug$network, zero_ids[[1]]))
+  nbrs <- skeleton_graph_neighbors(aug$graph, zero_ids[[1]])
   target <- aug$xy[nbrs[[1]], 1:2]
   expect_equal(target[[1]], 0.5, tolerance = 1e-9)
   expect_equal(target[[2]], 0.5, tolerance = 1e-9)
@@ -603,10 +588,10 @@ test_that("add_boundary_anchors rejects early skeleton crossings", {
   anchor <- geos::as_geos_geometry("POINT (0 5)")
   augmented <- add_boundary_anchors(skeleton, polygon, anchor)
 
-  aug <- skeleton_nodes_network(augmented)
-  d <- as.numeric(sf::st_distance(sf::st_as_sf(anchor), aug$nodes))
+  aug <- skeleton_graph_info(augmented)
+  d <- as.numeric(geos::geos_distance(anchor, aug$nodes))
   zero_ids <- which(d == 0)
-  nbrs <- as.integer(igraph::neighbors(aug$network, zero_ids[[1]]))
+  nbrs <- skeleton_graph_neighbors(aug$graph, zero_ids[[1]])
   target <- aug$xy[nbrs[[1]], 1:2]
   expect_equal(target[[1]], 8, tolerance = 1e-9)
   expect_equal(target[[2]], 8, tolerance = 1e-9)
@@ -639,12 +624,12 @@ test_that("add_boundary_anchors falls back to first skeleton hit", {
   augmented <- add_boundary_anchors(skeleton, polygon, anchor)
 
   expect_true(geos::geos_covers(augmented, skeleton))
-  aug <- skeleton_nodes_network(augmented)
-  d <- as.numeric(sf::st_distance(sf::st_as_sf(anchor), aug$nodes))
+  aug <- skeleton_graph_info(augmented)
+  d <- as.numeric(geos::geos_distance(anchor, aug$nodes))
   zero_ids <- which(d == 0)
   expect_equal(length(zero_ids), 1L)
   expect_equal(unname(aug$degree[[zero_ids[[1]]]]), 1L)
-  nbrs <- as.integer(igraph::neighbors(aug$network, zero_ids[[1]]))
+  nbrs <- skeleton_graph_neighbors(aug$graph, zero_ids[[1]])
   expect_equal(length(nbrs), 1L)
   target <- aug$xy[nbrs[[1]], 1:2]
   # First hit is the blocking segment at x = 1, y = 5.
@@ -701,9 +686,9 @@ test_that("straight skeleton accepts the same anchors contract", {
   path <- expect_no_warning(cnt_path(sk, anchors[1], anchors[2]))
   expect_equal(geos::geos_distance(geos::geos_point_start(path), anchors[1]), 0)
   expect_equal(geos::geos_distance(geos::geos_point_end(path), anchors[2]), 0)
-  aug <- skeleton_nodes_network(sk)
+  aug <- skeleton_graph_info(sk)
   for (i in seq_along(anchors)) {
-    d <- as.numeric(sf::st_distance(sf::st_as_sf(anchors[i]), aug$nodes))
+    d <- as.numeric(geos::geos_distance(anchors[i], aug$nodes))
     expect_equal(unname(aug$degree[[which(d == 0)[[1]]]]), 1L)
   }
 })
@@ -902,3 +887,37 @@ test_that("voronoi anchors at keep < 1 yield zero-distance cnt_path ends", {
   expect_equal(geos::geos_distance(geos::geos_point_start(path), anchors[1]), 0)
   expect_equal(geos::geos_distance(geos::geos_point_end(path), anchors[2]), 0)
 })
+
+
+test_that("ordinary and anchored skeletons match graph-migration baseline", {
+  baseline <- readRDS(test_path("fixtures/graph-migration-baseline.rds"))
+  f <- system.file("extdata/example.gpkg", package = "centerline")
+  p <- geos::as_geos_geometry(sf::st_read(f, layer = "polygon", quiet = TRUE))
+  pts <- geos::as_geos_geometry(
+    sf::st_read(f, layer = "polygon_points", quiet = TRUE)
+  )
+  b <- geos::geos_boundary(p)
+  a <- pts[geos::geos_equals(geos::geos_intersection(b, pts), pts)][1:2]
+
+  ordinary <- cnt_skeleton(p, keep = 1)
+  anchored <- cnt_skeleton(p, keep = 1, anchors = a)
+
+  for (item in list(
+    list(ordinary, baseline$ordinary_skeleton_wkb, baseline$lengths[["ordinary_skeleton"]]),
+    list(anchored, baseline$anchored_skeleton_wkb, baseline$lengths[["anchored_skeleton"]])
+  )) {
+    actual <- item[[1]]
+    base_len <- item[[3]]
+    tol <- max(1e-3, 1e-6 * base_len)
+    expect_equal(as.numeric(geos::geos_length(actual)), base_len, tolerance = tol)
+    base_geom <- geos::as_geos_geometry(item[[2]])
+    equal_exact <- tryCatch(
+      isTRUE(geos::geos_equals_exact(actual, base_geom, tol)),
+      error = function(e) FALSE
+    )
+    if (!equal_exact) {
+      hd <- as.numeric(geos::geos_distance_hausdorff(actual, base_geom))
+      expect_lte(hd, tol)
+    }
+  }
+})
diff --git a/tests/testthat/test-skeleton_graph.R b/tests/testthat/test-skeleton_graph.R
new file mode 100644
index 0000000..2a63e1f
--- /dev/null
+++ b/tests/testthat/test-skeleton_graph.R
@@ -0,0 +1,169 @@
+test_that("triangle of three edges has correct topology and weights", {
+  lines <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 0)",
+    "LINESTRING(1 0, 1 1)",
+    "LINESTRING(1 1, 0 0)"
+  ))
+  g <- skeleton_graph_build(lines)
+
+  expect_s3_class(g, "cnt_skeleton_graph")
+  expect_equal(length(g$nodes$x), 3L)
+  expect_equal(length(g$edges$weight), 3L)
+  expect_equal(as.integer(igraph::vcount(g$graph)), 3L)
+  expect_equal(as.integer(igraph::ecount(g$graph)), 3L)
+  expect_equal(sort(skeleton_graph_degree(g)), c(2L, 2L, 2L))
+  expect_equal(g$edges$weight, as.numeric(geos::geos_length(g$edges$geom)))
+  expect_true(g$meta$noded)
+  expect_true(is.na(g$meta$grid_size))
+})
+
+test_that("two disconnected segments stay disconnected", {
+  lines <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 0)",
+    "LINESTRING(0 1, 1 1)"
+  ))
+  g <- skeleton_graph_build(lines)
+
+  expect_equal(as.integer(igraph::vcount(g$graph)), 4L)
+  expect_equal(as.integer(igraph::ecount(g$graph)), 2L)
+  expect_equal(sort(skeleton_graph_degree(g)), c(1L, 1L, 1L, 1L))
+
+  from <- skeleton_graph_nearest_nodes(
+    g,
+    geos::as_geos_geometry("POINT(0 0)")
+  )
+  to <- skeleton_graph_nearest_nodes(
+    g,
+    geos::as_geos_geometry("POINT(0 1)")
+  )
+  epath <- skeleton_graph_paths(g, from, to)[[1]]
+  expect_equal(length(epath), 0L)
+})
+
+test_that("four separately supplied ring edges node into a closed ring", {
+  lines <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 0)",
+    "LINESTRING(1 0, 1 1)",
+    "LINESTRING(1 1, 0 1)",
+    "LINESTRING(0 1, 0 0)"
+  ))
+  g <- skeleton_graph_build(lines)
+
+  expect_equal(as.integer(igraph::vcount(g$graph)), 4L)
+  expect_equal(as.integer(igraph::ecount(g$graph)), 4L)
+  expect_equal(sort(skeleton_graph_degree(g)), c(2L, 2L, 2L, 2L))
+})
+
+test_that("crossing X yields five nodes and center degree four", {
+  lines <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 1)",
+    "LINESTRING(0 1, 1 0)"
+  ))
+  g <- skeleton_graph_build(lines)
+
+  expect_equal(as.integer(igraph::vcount(g$graph)), 5L)
+  expect_equal(as.integer(igraph::ecount(g$graph)), 4L)
+  expect_equal(sort(skeleton_graph_degree(g)), c(1L, 1L, 1L, 1L, 4L))
+
+  center <- skeleton_graph_nearest_nodes(
+    g,
+    geos::as_geos_geometry("POINT(0.5 0.5)")
+  )
+  expect_equal(skeleton_graph_degree(g)[[center]], 4L)
+  expect_equal(length(skeleton_graph_neighbors(g, center)), 4L)
+})
+
+test_that("exact endpoints share a node; near-touch needs grid_size", {
+  exact <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 0)",
+    "LINESTRING(1 0, 2 0)"
+  ))
+  g_exact <- skeleton_graph_build(exact)
+  expect_equal(as.integer(igraph::vcount(g_exact$graph)), 3L)
+  expect_equal(as.integer(igraph::ecount(g_exact$graph)), 2L)
+
+  eps <- 1e-12
+  near <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 1 0)",
+    sprintf("LINESTRING(%.16f 0, 2 0)", 1 + eps)
+  ))
+  g_near <- skeleton_graph_build(near)
+  # Distinct IEEE doubles keep separate nodes without precision.
+  expect_equal(as.integer(igraph::vcount(g_near$graph)), 4L)
+
+  g_snap <- skeleton_graph_build(near, grid_size = 1e-9)
+  expect_equal(as.integer(igraph::vcount(g_snap$graph)), 3L)
+  expect_equal(g_snap$meta$grid_size, 1e-9)
+})
+
+test_that("edge weights equal GEOS lengths", {
+  lines <- geos::as_geos_geometry(c(
+    "LINESTRING(0 0, 3 0)",
+    "LINESTRING(0 0, 0 4)"
+  ))
+  g <- skeleton_graph_build(lines)
+  expect_equal(sort(g$edges$weight), c(3, 4))
+  expect_equal(g$edges$weight, as.numeric(geos::geos_length(g$edges$geom)))
+})
+
+test_that("parallel edges keep both geometries and select shorter epath", {
+  short <- geos::as_geos_geometry("LINESTRING(0 0, 1 0)")
+  long <- geos::as_geos_geometry("LINESTRING(0 0, 0 1, 1 1, 1 0)")
+  lines <- c(short, long)
+  g <- skeleton_graph_build(lines)
+
+  expect_equal(as.integer(igraph::ecount(g$graph)), 2L)
+  expect_equal(as.integer(igraph::vcount(g$graph)), 2L)
+  expect_equal(sort(g$edges$weight), sort(as.numeric(geos::geos_length(lines))))
+
+  from <- skeleton_graph_nearest_nodes(
+    g,
+    geos::as_geos_geometry("POINT(0 0)")
+  )
+  to <- skeleton_graph_nearest_nodes(
+    g,
+    geos::as_geos_geometry("POINT(1 0)")
+  )
+  epath <- skeleton_graph_paths(g, from, to)[[1]]
+  expect_equal(length(epath), 1L)
+
+  path_line <- skeleton_graph_path_line(g, epath)
+  expect_equal(
+    as.numeric(geos::geos_length(path_line)),
+    min(g$edges$weight),
+    tolerance = 1e-12
+  )
+  # Edge-ID order matches edges$geom: selected edge is the short segment.
+  expect_true(
+    isTRUE(geos::geos_equals_exact(path_line, short, 1e-9)) ||
+      isTRUE(geos::geos_equals_exact(
+        path_line,
+        geos::geos_reverse(short),
+        1e-9
+      ))
+  )
+})
+
+test_that("path_line rejects empty edge paths", {
+  lines <- geos::as_geos_geometry("LINESTRING(0 0, 1 0)")
+  g <- skeleton_graph_build(lines)
+  expect_error(skeleton_graph_path_line(g, integer()), "empty edge path")
+})
+
+test_that("node_points returns stored node geometries", {
+  lines <- geos::as_geos_geometry("LINESTRING(0 0, 1 0)")
+  g <- skeleton_graph_build(lines)
+  pts <- skeleton_graph_node_points(g)
+  expect_equal(length(pts), 2L)
+  expect_identical(pts, g$nodes$geom)
+})
+
+test_that("builder validates grid_size and empty input", {
+  lines <- geos::as_geos_geometry("LINESTRING(0 0, 1 0)")
+  expect_error(skeleton_graph_build(lines, grid_size = 0), "positive finite")
+  expect_error(skeleton_graph_build(lines, grid_size = c(1, 2)), "positive finite")
+  expect_error(
+    skeleton_graph_build(geos::as_geos_geometry("POINT(0 0)")),
+    "line geometries"
+  )
+})