alekrutkowski · July 2, 2026 14:56
diff --git a/examples.md b/examples.md
diff --git a/functionToSTL.R b/functionToSTL.R
 #' Convert a bivariate function surface to ASCII STL
 #'
 #' Create an ASCII STL representation of a triangulated surface defined by a
 #' function `f(x, y)`. The output can optionally be wrapped in a GitHub Markdown
 #' `stl` code block for interactive 3D rendering.
 #'
 #' The surface is sampled over rectangular `x` and `y` ranges. The resulting
 #' `z` values are clipped to `zlim`, if supplied, and triangulated. Optional
 #' axes, tick marks, and numeric tick labels are added as STL geometry because
 #' STL has no native support for text or annotations.
 #'
 #' @param f A function taking numeric `x` and `y` values and returning a numeric
 #'   `z` value. By default, `f` is called point-by-point. If `vectorized = TRUE`,
 #'   `f` must accept equal-length numeric vectors `x` and `y` and return one
 #'   numeric `z` value per pair.
 #' @param xlim Numeric vector of length 2 giving the lower and upper limits for
 #'   the `x` axis.
 #' @param ylim Numeric vector of length 2 giving the lower and upper limits for
 #'   the `y` axis.
 #' @param zlim Optional numeric vector of length 2 giving the lower and upper
 #'   limits for the `z` axis. If `NULL`, the range of finite sampled `z` values
 #'   is used.
 #' @param resolution Numeric value between 0.1 and 1 controlling mesh detail.
 #'   Larger values create a finer triangulation and a larger STL string.
 #' @param name Character string used as the STL solid name. Characters other
 #'   than letters, numbers, underscores, dots, and hyphens are replaced with
 #'   underscores.
 #' @param n_ticks Integer giving the approximate number of tick marks to draw on
 #'   each axis. Tick positions are chosen with [pretty()].
 #' @param solid Logical. If `TRUE`, add side walls and a gridded bottom face at
 #'   `zlim[1]`. If `FALSE`, only the function surface is drawn.
 #' @param axes Logical. If `TRUE`, add 3D axes, tick marks, axis labels, and
 #'   numeric tick labels as STL geometry.
 #' @param markdown Logical. If `TRUE`, wrap the ASCII STL in a Markdown fenced
 #'   code block using the `stl` language identifier.
 #' @param vectorized Logical. If `TRUE`, call `f` once with vector inputs created
 #'   from all sampled `x, y` pairs. If `FALSE`, call `f` point-by-point.
 #' @param digits Integer giving the number of significant digits used when
 #'   formatting STL vertex coordinates and facet normals.
 #' @param label_digits Integer giving the number of significant digits used when
 #'   formatting numeric axis tick labels.
 #'
 #' @return A character scalar containing ASCII STL code. If `markdown = TRUE`,
 #'   the returned string is wrapped in a Markdown fenced `stl` block.
 #'
 #' @details
 #' This function produces a visual STL mesh rather than a mathematically exact
 #' surface. Axis text is approximated using stroke geometry, so numbers and
 #' letters are part of the mesh itself. Large values of `resolution` can create
 #' very large strings, especially when `solid = TRUE`.
 functionToSTL <- function(f, xlim, ylim, zlim = NULL,
                          resolution = 0.5,
                          name = "function_surface",
                          n_ticks = 5L,
                          solid = TRUE,
                          axes = TRUE,
                          markdown = FALSE,
                          vectorized = FALSE,
                          digits = 7L,
                          label_digits = 4L) {
  stopifnot(is.function(f), length(xlim) == 2L, length(ylim) == 2L)

  resolution <- as.numeric(resolution)[1L]
  if (!is.finite(resolution) || resolution < 0.1 || resolution > 1) {
    stop("resolution must be a number between 0.1 and 1.")
  }

  xlim <- sort(as.numeric(xlim))
  ylim <- sort(as.numeric(ylim))
  if (diff(xlim) <= 0 || diff(ylim) <= 0) {
    stop("xlim and ylim must each contain two distinct values.")
  }

  n <- as.integer(round(10 + 90 * resolution))
  nx <- ny <- n

  x <- seq(xlim[1L], xlim[2L], length.out = nx)
  y <- seq(ylim[1L], ylim[2L], length.out = ny)

  if (isTRUE(vectorized)) {
    xy <- expand.grid(x = x, y = y)
    zz <- f(xy$x, xy$y)
    if (length(zz) != nrow(xy)) {
      stop("When vectorized = TRUE, f(x, y) must return one z value per x/y pair.")
    }
    z <- matrix(as.numeric(zz), nrow = nx, ncol = ny)
  } else {
    z <- outer(x, y, Vectorize(function(a, b) as.numeric(f(a, b))[1L]))
  }

  z[!is.finite(z)] <- NA_real_

  if (is.null(zlim)) {
    if (!any(is.finite(z))) stop("f returned no finite z values.")
    zlim <- range(z, na.rm = TRUE)
  } else {
    zlim <- sort(as.numeric(zlim))
  }

  if (length(zlim) != 2L || any(!is.finite(zlim)) || diff(zlim) <= 0) {
    stop("zlim must contain two distinct finite values.")
  }

  z <- pmin(pmax(z, zlim[1L]), zlim[2L])

  cross <- function(a, b) {
    c(
      a[2L] * b[3L] - a[3L] * b[2L],
      a[3L] * b[1L] - a[1L] * b[3L],
      a[1L] * b[2L] - a[2L] * b[1L]
    )
  }

  unit <- function(a) {
    s <- sqrt(sum(a * a))
    if (!is.finite(s) || s == 0) c(0, 0, 0) else a / s
  }

  fmt <- function(v) {
    formatC(v, format = "fg", digits = digits, flag = "#")
  }

  vertex_line <- function(p) {
    sprintf("      vertex %s %s %s", fmt(p[1L]), fmt(p[2L]), fmt(p[3L]))
  }

  facet <- function(a, b, c) {
    if (!all(is.finite(c(a, b, c)))) return(character())

    nr <- cross(b - a, c - a)
    ns <- sqrt(sum(nr * nr))
    if (!is.finite(ns) || ns == 0) return(character())

    n <- nr / ns

    c(
      sprintf("  facet normal %s %s %s", fmt(n[1L]), fmt(n[2L]), fmt(n[3L])),
      "    outer loop",
      vertex_line(a),
      vertex_line(b),
      vertex_line(c),
      "    endloop",
      "  endfacet"
    )
  }

  line_prism <- function(a, b, width) {
    if (!is.finite(width) || width <= 0) return(character())

    d <- b - a
    len <- sqrt(sum(d * d))
    if (!is.finite(len) || len == 0) return(character())

    e <- d / len
    ref <- if (abs(e[3L]) < 0.85) c(0, 0, 1) else c(0, 1, 0)

    n1 <- unit(cross(e, ref))
    if (all(n1 == 0)) n1 <- c(1, 0, 0)
    n2 <- unit(cross(e, n1))

    r <- width / 2

    v <- list(
      a + n1 * r + n2 * r,
      a - n1 * r + n2 * r,
      a - n1 * r - n2 * r,
      a + n1 * r - n2 * r,
      b + n1 * r + n2 * r,
      b - n1 * r + n2 * r,
      b - n1 * r - n2 * r,
      b + n1 * r - n2 * r
    )

    faces <- list(
      c(1, 2, 6), c(1, 6, 5),
      c(2, 3, 7), c(2, 7, 6),
      c(3, 4, 8), c(3, 8, 7),
      c(4, 1, 5), c(4, 5, 8),
      c(1, 4, 3), c(1, 3, 2),
      c(5, 6, 7), c(5, 7, 8)
    )

    unlist(
      lapply(faces, function(ii) {
        facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
      }),
      use.names = FALSE
    )
  }

  glyph <- local({
    seg <- function(x1, y1, x2, y2) {
      list(c(x1, y1, x2, y2))
    }

    arc <- function(cx, cy, rx, ry, a0, a1, n = 10L) {
      th <- seq(a0, a1, length.out = n + 1L)
      p <- cbind(cx + rx * cos(th), cy + ry * sin(th))

      lapply(seq_len(nrow(p) - 1L), function(i) {
        c(p[i, 1L], p[i, 2L], p[i + 1L, 1L], p[i + 1L, 2L])
      })
    }

    join <- function(...) {
      unlist(list(...), recursive = FALSE, use.names = FALSE)
    }

    function(ch) {
      switch(toupper(ch),
        "0" = arc(0.5, 0.5, 0.42, 0.50, 0, 2 * pi, 18L),

        "1" = join(
          seg(0.35, 0.82, 0.55, 1.00),
          seg(0.55, 1.00, 0.55, 0.00),
          seg(0.35, 0.00, 0.75, 0.00)
        ),

        "2" = join(
          arc(0.50, 0.72, 0.40, 0.28, pi, 0, 10L),
          seg(0.90, 0.72, 0.12, 0.00),
          seg(0.12, 0.00, 0.90, 0.00)
        ),

        "3" = join(
          seg(0.22, 1.00, 0.42, 1.00),
          arc(0.42, 0.75, 0.42, 0.25, pi / 2, -pi / 2, 10L),
          arc(0.42, 0.25, 0.42, 0.25, pi / 2, -pi / 2, 10L),
          seg(0.22, 0.00, 0.42, 0.00)
        ),

        "4" = join(
          seg(0.78, 1.00, 0.78, 0.00),
          seg(0.18, 0.55, 0.90, 0.55),
          seg(0.18, 0.55, 0.72, 1.00)
        ),

        "5" = join(
          seg(0.88, 1.00, 0.20, 1.00),
          seg(0.20, 1.00, 0.20, 0.56),
          seg(0.20, 0.56, 0.50, 0.56),
          arc(0.50, 0.28, 0.38, 0.28, pi / 2, -pi / 2, 10L),
          seg(0.50, 0.00, 0.22, 0.00)
        ),

        "6" = join(
          seg(0.76, 0.96, 0.28, 0.58),
          arc(0.50, 0.32, 0.36, 0.32, 0, 2 * pi, 16L)
        ),

        "7" = join(
          seg(0.12, 1.00, 0.90, 1.00),
          seg(0.90, 1.00, 0.35, 0.00)
        ),

        "8" = join(
          arc(0.50, 0.73, 0.35, 0.27, 0, 2 * pi, 14L),
          arc(0.50, 0.27, 0.38, 0.27, 0, 2 * pi, 14L)
        ),

        "9" = join(
          arc(0.50, 0.68, 0.36, 0.32, 0, 2 * pi, 16L),
          seg(0.78, 0.45, 0.34, 0.04)
        ),

        "-" = seg(0.16, 0.50, 0.84, 0.50),

        "+" = join(
          seg(0.16, 0.50, 0.84, 0.50),
          seg(0.50, 0.16, 0.50, 0.84)
        ),

        "." = arc(0.50, 0.06, 0.08, 0.08, 0, 2 * pi, 8L),

        "X" = join(
          seg(0.12, 0.00, 0.88, 1.00),
          seg(0.12, 1.00, 0.88, 0.00)
        ),

        "Y" = join(
          seg(0.12, 1.00, 0.50, 0.52),
          seg(0.88, 1.00, 0.50, 0.52),
          seg(0.50, 0.52, 0.50, 0.00)
        ),

        "Z" = join(
          seg(0.12, 1.00, 0.88, 1.00),
          seg(0.88, 1.00, 0.12, 0.00),
          seg(0.12, 0.00, 0.88, 0.00)
        ),

        list()
      )
    }
  })

  flat_stroke <- function(a, b, width) {
    d <- b - a
    len <- sqrt(sum(d * d))
    if (!is.finite(len) || len == 0) return(character())

    e <- d / len
    n <- unit(cross(e, c(0, 0, 1)))

    if (all(n == 0)) {
      n <- unit(cross(e, c(0, 1, 0)))
    }

    r <- width / 2

    p1 <- a + n * r
    p2 <- a - n * r
    p3 <- b - n * r
    p4 <- b + n * r

    c(
      facet(p1, p2, p3),
      facet(p1, p3, p4)
    )
  }

  draw_text <- function(txt, center, u, v, size, width,
                        flat_text = TRUE) {
    chars <- strsplit(as.character(txt), "", fixed = TRUE)[[1L]]
    if (!length(chars)) return(character())

    u <- unit(u)
    v <- unit(v)

    char_w <- 1
    gap <- 0.35
    total_w <- length(chars) * char_w + max(0, length(chars) - 1L) * gap

    x0 <- -total_w / 2
    y0 <- -0.5

    out <- vector("list", 0L)
    stroke_fun <- if (isTRUE(flat_text)) flat_stroke else line_prism

    for (ii in seq_along(chars)) {
      segs <- glyph(chars[ii])
      if (!length(segs)) next

      shift <- x0 + (ii - 1L) * (char_w + gap)

      for (sg in segs) {
        p1 <- center + u * ((shift + sg[1L]) * size) + v * ((y0 + sg[2L]) * size)
        p2 <- center + u * ((shift + sg[3L]) * size) + v * ((y0 + sg[4L]) * size)
        out[[length(out) + 1L]] <- stroke_fun(p1, p2, width)
      }
    }

    unlist(out, use.names = FALSE)
  }

  tick_values <- function(lim) {
    vals <- pretty(lim, n = n_ticks)
    eps <- sqrt(.Machine$double.eps) * max(1, diff(lim))
    vals[vals >= lim[1L] - eps & vals <= lim[2L] + eps]
  }

  tick_labels <- function(vals) {
    vals[abs(vals) < sqrt(.Machine$double.eps)] <- 0
    out <- format(signif(vals, label_digits), trim = TRUE,
                  scientific = FALSE, drop0trailing = TRUE)
    out <- gsub(" ", "", out, fixed = TRUE)
    out[out == "-0"] <- "0"
    out
  }

  pieces <- list()

  add <- function(lines) {
    if (length(lines)) {
      pieces[[length(pieces) + 1L]] <<- lines
    }
  }

  P <- function(i, j) c(x[i], y[j], z[i, j])
  Q <- function(i, j) c(x[i], y[j], zlim[1L])

  for (i in seq_len(nx - 1L)) {
    for (j in seq_len(ny - 1L)) {
      p00 <- P(i, j)
      p10 <- P(i + 1L, j)
      p01 <- P(i, j + 1L)
      p11 <- P(i + 1L, j + 1L)

      add(facet(p00, p10, p11))
      add(facet(p00, p11, p01))
    }
  }

  if (isTRUE(solid)) {
    for (i in seq_len(nx - 1L)) {
      p0 <- P(i, 1L)
      p1 <- P(i + 1L, 1L)
      q0 <- Q(i, 1L)
      q1 <- Q(i + 1L, 1L)

      add(facet(p0, q0, q1))
      add(facet(p0, q1, p1))

      p0 <- P(i, ny)
      p1 <- P(i + 1L, ny)
      q0 <- Q(i, ny)
      q1 <- Q(i + 1L, ny)

      add(facet(p0, q1, q0))
      add(facet(p0, p1, q1))
    }

    for (j in seq_len(ny - 1L)) {
      p0 <- P(1L, j)
      p1 <- P(1L, j + 1L)
      q0 <- Q(1L, j)
      q1 <- Q(1L, j + 1L)

      add(facet(p0, q1, q0))
      add(facet(p0, p1, q1))

      p0 <- P(nx, j)
      p1 <- P(nx, j + 1L)
      q0 <- Q(nx, j)
      q1 <- Q(nx, j + 1L)

      add(facet(p0, q0, q1))
      add(facet(p0, q1, p1))
    }

    for (i in seq_len(nx - 1L)) {
      for (j in seq_len(ny - 1L)) {
        q00 <- c(x[i], y[j], zlim[1L])
        q10 <- c(x[i + 1L], y[j], zlim[1L])
        q01 <- c(x[i], y[j + 1L], zlim[1L])
        q11 <- c(x[i + 1L], y[j + 1L], zlim[1L])

        add(facet(q00, q11, q10))
        add(facet(q00, q01, q11))
      }
    }
  }

  if (isTRUE(axes)) {
    span <- max(diff(xlim), diff(ylim), diff(zlim))
    pad <- 0.08 * span
    axis_w <- 0.006 * span
    tick_len <- 0.035 * span
    label_size <- 0.045 * span
    label_w <- 0.45 * axis_w

    anchor <- c(xlim[1L] - pad, ylim[1L] - pad, zlim[1L])

    x_axis_end <- c(xlim[2L] + pad / 2, anchor[2L], zlim[1L])
    y_axis_end <- c(anchor[1L], ylim[2L] + pad / 2, zlim[1L])
    z_axis_end <- c(anchor[1L], anchor[2L], zlim[2L] + pad / 2)

    add(line_prism(anchor, x_axis_end, axis_w))
    add(line_prism(anchor, y_axis_end, axis_w))
    add(line_prism(anchor, z_axis_end, axis_w))

    xt <- tick_values(xlim)
    yt <- tick_values(ylim)
    zt <- tick_values(zlim)

    xl <- tick_labels(xt)
    yl <- tick_labels(yt)
    zl <- tick_labels(zt)

    for (k in seq_along(xt)) {
      xx <- xt[k]

      add(line_prism(
        c(xx, anchor[2L] - tick_len / 2, zlim[1L]),
        c(xx, anchor[2L] + tick_len / 2, zlim[1L]),
        axis_w * 0.75
      ))

      add(draw_text(
        xl[k],
        c(xx, anchor[2L] - 1.7 * label_size, zlim[1L]),
        u = c(1, 0, 0),
        v = c(0, 1, 0),
        size = label_size,
        width = label_w
      ))
    }

    for (k in seq_along(yt)) {
      yy <- yt[k]

      add(line_prism(
        c(anchor[1L] - tick_len / 2, yy, zlim[1L]),
        c(anchor[1L] + tick_len / 2, yy, zlim[1L]),
        axis_w * 0.75
      ))

      add(draw_text(
        yl[k],
        c(anchor[1L] - 1.7 * label_size, yy, zlim[1L]),
        u = c(0, 1, 0),
        v = c(-1, 0, 0),
        size = label_size,
        width = label_w
      ))
    }

    z_label_offset <- (max(nchar(zl), 1L) * 0.45 + 1.9) * label_size

    for (k in seq_along(zt)) {
      zz <- zt[k]

      add(line_prism(
        c(anchor[1L] - tick_len / 2, anchor[2L], zz),
        c(anchor[1L] + tick_len / 2, anchor[2L], zz),
        axis_w * 0.75
      ))

      add(draw_text(
        zl[k],
        c(anchor[1L] - z_label_offset, anchor[2L], zz),
        u = c(1, 0, 0),
        v = c(0, 0, 1),
        size = label_size,
        width = label_w
      ))
    }

    add(draw_text(
      "X",
      c(xlim[2L] + 1.8 * label_size, anchor[2L], zlim[1L]),
      u = c(1, 0, 0),
      v = c(0, 1, 0),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))

    add(draw_text(
      "Y",
      c(anchor[1L], ylim[2L] + 1.8 * label_size, zlim[1L]),
      u = c(0, 1, 0),
      v = c(-1, 0, 0),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))

    add(draw_text(
      "Z",
      c(anchor[1L], anchor[2L], zlim[2L] + 1.8 * label_size),
      u = c(1, 0, 0),
      v = c(0, 0, 1),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))
  }

  body <- unlist(pieces, use.names = FALSE)
  name <- gsub("[^A-Za-z0-9_.-]", "_", name)

  out <- paste(
    c(sprintf("solid %s", name), body, sprintf("endsolid %s", name)),
    collapse = "\n"
  )

  if (isTRUE(markdown)) {
    out <- paste0("```stl\n", out, "\n```")
  }

  out
 }
diff --git a/pointsToSTL.R b/pointsToSTL.R
 #' Convert 3D points to ASCII STL
 #'
 #' Create an ASCII STL representation of a 3D point cloud. Each point is drawn as
 #' a small cube or sphere, and optional axes, tick marks, and numeric tick labels
 #' are added as STL geometry.
 #'
 #' The output can optionally be wrapped in a GitHub Markdown `stl` code block for
 #' interactive 3D rendering.
 #'
 #' @param x Numeric vector of `x` coordinates.
 #' @param y Numeric vector of `y` coordinates. Must have the same length as `x`.
 #' @param z Numeric vector of `z` coordinates. Must have the same length as `x`.
 #' @param shape Character string specifying how points are drawn. One of
 #'   `"cube"` or `"sphere"`.
 #' @param xlim Optional numeric vector of length 2 giving the lower and upper
 #'   limits for the `x` axis. If `NULL`, the finite range of `x` is used.
 #' @param ylim Optional numeric vector of length 2 giving the lower and upper
 #'   limits for the `y` axis. If `NULL`, the finite range of `y` is used.
 #' @param zlim Optional numeric vector of length 2 giving the lower and upper
 #'   limits for the `z` axis. If `NULL`, the finite range of `z` is used.
 #' @param resolution Numeric value between 0.1 and 1 controlling marker detail.
 #'   For spheres, larger values create smoother spheres. If `point_size = NULL`,
 #'   larger values may also be used to create smaller default markers, depending
 #'   on the implementation.
 #' @param point_size Optional positive numeric value giving the marker diameter
 #'   in data units. If `NULL`, a size is chosen automatically from the largest
 #'   axis span and `resolution`.
 #' @param name Character string used as the STL solid name. Characters other
 #'   than letters, numbers, underscores, dots, and hyphens are replaced with
 #'   underscores.
 #' @param n_ticks Integer giving the approximate number of tick marks to draw on
 #'   each axis. Tick positions are chosen with [pretty()]. Use `0` to suppress
 #'   tick marks and tick labels while keeping axes.
 #' @param axes Logical. If `TRUE`, add 3D axes, tick marks, axis labels, and
 #'   numeric tick labels as STL geometry.
 #' @param markdown Logical. If `TRUE`, wrap the ASCII STL in a Markdown fenced
 #'   code block using the `stl` language identifier.
 #' @param digits Integer giving the number of significant digits used when
 #'   formatting STL vertex coordinates and facet normals.
 #' @param label_digits Integer giving the number of significant digits used when
 #'   formatting numeric axis tick labels.
 #' @param flat_text Logical. If `TRUE`, draw axis labels and tick labels as flat
 #'   ribbon strokes. If `FALSE`, draw text strokes as small rectangular prisms.
 #'   Flat text usually looks cleaner in STL wireframe renderers.
 #' @param clip Logical. If `TRUE`, drop points outside `xlim`, `ylim`, and
 #'   `zlim`. If `FALSE`, all finite points are rendered, even if they lie outside
 #'   the displayed axis limits.
 #'
 #' @return A character scalar containing ASCII STL code. If `markdown = TRUE`,
 #'   the returned string is wrapped in a Markdown fenced `stl` block.
 #'
 #' @details
 #' STL has no native point, text, axis, or tick-mark primitives, so all elements
 #' are represented as triangular mesh facets. Cubes are much lighter than
 #' spheres. Spheres can produce large STL strings, especially with many points
 #' and high `resolution`.
 pointsToSTL <- function(x, y, z,
                        shape = c("cube", "sphere"),
                        xlim = NULL, ylim = NULL, zlim = NULL,
                        resolution = 0.5,
                        point_size = NULL,
                        name = "point_cloud",
                        n_ticks = 5L,
                        axes = TRUE,
                        markdown = FALSE,
                        digits = 7L,
                        label_digits = 4L,
                        flat_text = TRUE,
                        clip = TRUE) {
  shape <- match.arg(shape)

  if (length(x) != length(y) || length(x) != length(z)) {
    stop("x, y, and z must have equal lengths.")
  }

  if (!length(x)) stop("x, y, and z must not be empty.")

  resolution <- as.numeric(resolution)[1L]
  if (!is.finite(resolution) || resolution < 0.1 || resolution > 1) {
    stop("resolution must be a number between 0.1 and 1.")
  }

  x <- as.numeric(x)
  y <- as.numeric(y)
  z <- as.numeric(z)

  ok <- is.finite(x) & is.finite(y) & is.finite(z)
  if (!any(ok)) stop("x, y, and z contain no finite complete points.")

  if (!all(ok)) {
    warning("Dropping non-finite points.")
    x <- x[ok]
    y <- y[ok]
    z <- z[ok]
  }

  make_lim <- function(v, lim, nm) {
    if (is.null(lim)) {
      lim <- range(v)
    } else {
      if (length(lim) != 2L) stop(nm, " must contain two values.")
      lim <- sort(as.numeric(lim))
    }

    if (any(!is.finite(lim))) stop(nm, " must contain finite values.")

    if (diff(lim) == 0) {
      d <- max(1, abs(lim[1L])) * 0.05
      lim <- lim + c(-d, d)
    }

    lim
  }

  xlim <- make_lim(x, xlim, "xlim")
  ylim <- make_lim(y, ylim, "ylim")
  zlim <- make_lim(z, zlim, "zlim")

  if (isTRUE(clip)) {
    keep <- x >= xlim[1L] & x <= xlim[2L] &
      y >= ylim[1L] & y <= ylim[2L] &
      z >= zlim[1L] & z <= zlim[2L]

    if (!any(keep)) stop("No points remain after clipping to xlim, ylim, and zlim.")

    x <- x[keep]
    y <- y[keep]
    z <- z[keep]
  }

  span <- max(diff(xlim), diff(ylim), diff(zlim))
  if (!is.finite(span) || span <= 0) span <- 1

  if (is.null(point_size)) {
    # Higher resolution means smaller markers.
    # resolution = 0.1 -> diameter about 5.5% of the largest axis span
    # resolution = 1.0 -> diameter about 1.2% of the largest axis span
    point_size <- span * (0.012 + (1 - resolution) / 0.9 * (0.055 - 0.012))
  } else {
    point_size <- as.numeric(point_size)[1L]
  }

  if (!is.finite(point_size) || point_size <= 0) {
    stop("point_size must be a positive finite number.")
  }

  n_ticks <- as.integer(n_ticks)[1L]
  if (!is.finite(n_ticks) || n_ticks < 0L) {
    stop("n_ticks must be a non-negative integer.")
  }

  cross <- function(a, b) {
    c(
      a[2L] * b[3L] - a[3L] * b[2L],
      a[3L] * b[1L] - a[1L] * b[3L],
      a[1L] * b[2L] - a[2L] * b[1L]
    )
  }

  unit <- function(a) {
    s <- sqrt(sum(a * a))
    if (!is.finite(s) || s == 0) c(0, 0, 0) else a / s
  }

  fmt <- function(v) {
    formatC(v, format = "fg", digits = digits, flag = "#")
  }

  vertex_line <- function(p) {
    sprintf("      vertex %s %s %s", fmt(p[1L]), fmt(p[2L]), fmt(p[3L]))
  }

  facet <- function(a, b, d) {
    if (!all(is.finite(c(a, b, d)))) return(character())

    nr <- cross(b - a, d - a)
    ns <- sqrt(sum(nr * nr))
    if (!is.finite(ns) || ns == 0) return(character())

    n <- nr / ns

    c(
      sprintf("  facet normal %s %s %s", fmt(n[1L]), fmt(n[2L]), fmt(n[3L])),
      "    outer loop",
      vertex_line(a),
      vertex_line(b),
      vertex_line(d),
      "    endloop",
      "  endfacet"
    )
  }

  line_prism <- function(a, b, width) {
    if (!is.finite(width) || width <= 0) return(character())

    d <- b - a
    len <- sqrt(sum(d * d))
    if (!is.finite(len) || len == 0) return(character())

    e <- d / len
    ref <- if (abs(e[3L]) < 0.85) c(0, 0, 1) else c(0, 1, 0)

    n1 <- unit(cross(e, ref))
    if (all(n1 == 0)) n1 <- c(1, 0, 0)
    n2 <- unit(cross(e, n1))

    r <- width / 2

    v <- list(
      a + n1 * r + n2 * r,
      a - n1 * r + n2 * r,
      a - n1 * r - n2 * r,
      a + n1 * r - n2 * r,
      b + n1 * r + n2 * r,
      b - n1 * r + n2 * r,
      b - n1 * r - n2 * r,
      b + n1 * r - n2 * r
    )

    faces <- list(
      c(1, 2, 6), c(1, 6, 5),
      c(2, 3, 7), c(2, 7, 6),
      c(3, 4, 8), c(3, 8, 7),
      c(4, 1, 5), c(4, 5, 8),
      c(1, 4, 3), c(1, 3, 2),
      c(5, 6, 7), c(5, 7, 8)
    )

    unlist(
      lapply(faces, function(ii) {
        facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
      }),
      use.names = FALSE
    )
  }

  flat_stroke <- function(a, b, width, plane_normal) {
    d <- b - a
    len <- sqrt(sum(d * d))
    if (!is.finite(len) || len == 0) return(character())

    e <- d / len
    plane_normal <- unit(plane_normal)
    if (all(plane_normal == 0)) plane_normal <- c(0, 0, 1)

    n <- unit(cross(plane_normal, e))
    if (all(n == 0)) n <- unit(cross(c(0, 0, 1), e))
    if (all(n == 0)) n <- c(1, 0, 0)

    r <- width / 2

    p1 <- a + n * r
    p2 <- a - n * r
    p3 <- b - n * r
    p4 <- b + n * r

    c(
      facet(p1, p2, p3),
      facet(p1, p3, p4)
    )
  }

  glyph <- local({
    arc_steps <- max(6L, as.integer(round(5 + 15 * resolution)))

    seg <- function(x1, y1, x2, y2) {
      list(c(x1, y1, x2, y2))
    }

    arc <- function(cx, cy, rx, ry, a0, a1, n = arc_steps) {
      n <- max(1L, as.integer(n))
      th <- seq(a0, a1, length.out = n + 1L)
      p <- cbind(cx + rx * cos(th), cy + ry * sin(th))

      lapply(seq_len(nrow(p) - 1L), function(i) {
        c(p[i, 1L], p[i, 2L], p[i + 1L, 1L], p[i + 1L, 2L])
      })
    }

    join <- function(...) {
      unlist(list(...), recursive = FALSE, use.names = FALSE)
    }

    function(ch) {
      switch(toupper(ch),
        "0" = arc(0.5, 0.5, 0.42, 0.50, 0, 2 * pi, arc_steps),

        "1" = join(
          seg(0.35, 0.82, 0.55, 1.00),
          seg(0.55, 1.00, 0.55, 0.00),
          seg(0.35, 0.00, 0.75, 0.00)
        ),

        "2" = join(
          arc(0.50, 0.72, 0.40, 0.28, pi, 0, arc_steps),
          seg(0.90, 0.72, 0.12, 0.00),
          seg(0.12, 0.00, 0.90, 0.00)
        ),

        "3" = join(
          seg(0.22, 1.00, 0.42, 1.00),
          arc(0.42, 0.75, 0.42, 0.25, pi / 2, -pi / 2, arc_steps),
          arc(0.42, 0.25, 0.42, 0.25, pi / 2, -pi / 2, arc_steps),
          seg(0.22, 0.00, 0.42, 0.00)
        ),

        "4" = join(
          seg(0.78, 1.00, 0.78, 0.00),
          seg(0.18, 0.55, 0.90, 0.55),
          seg(0.18, 0.55, 0.72, 1.00)
        ),

        "5" = join(
          seg(0.88, 1.00, 0.20, 1.00),
          seg(0.20, 1.00, 0.20, 0.56),
          seg(0.20, 0.56, 0.50, 0.56),
          arc(0.50, 0.28, 0.38, 0.28, pi / 2, -pi / 2, arc_steps),
          seg(0.50, 0.00, 0.22, 0.00)
        ),

        "6" = join(
          seg(0.76, 0.96, 0.28, 0.58),
          arc(0.50, 0.32, 0.36, 0.32, 0, 2 * pi, arc_steps)
        ),

        "7" = join(
          seg(0.12, 1.00, 0.90, 1.00),
          seg(0.90, 1.00, 0.35, 0.00)
        ),

        "8" = join(
          arc(0.50, 0.73, 0.35, 0.27, 0, 2 * pi, arc_steps),
          arc(0.50, 0.27, 0.38, 0.27, 0, 2 * pi, arc_steps)
        ),

        "9" = join(
          arc(0.50, 0.68, 0.36, 0.32, 0, 2 * pi, arc_steps),
          seg(0.78, 0.45, 0.34, 0.04)
        ),

        "-" = seg(0.16, 0.50, 0.84, 0.50),

        "+" = join(
          seg(0.16, 0.50, 0.84, 0.50),
          seg(0.50, 0.16, 0.50, 0.84)
        ),

        "." = arc(0.50, 0.06, 0.08, 0.08, 0, 2 * pi, max(6L, arc_steps %/% 2L)),

        "X" = join(
          seg(0.12, 0.00, 0.88, 1.00),
          seg(0.12, 1.00, 0.88, 0.00)
        ),

        "Y" = join(
          seg(0.12, 1.00, 0.50, 0.52),
          seg(0.88, 1.00, 0.50, 0.52),
          seg(0.50, 0.52, 0.50, 0.00)
        ),

        "Z" = join(
          seg(0.12, 1.00, 0.88, 1.00),
          seg(0.88, 1.00, 0.12, 0.00),
          seg(0.12, 0.00, 0.88, 0.00)
        ),

        list()
      )
    }
  })

  draw_text <- function(txt, center, u, v, size, width) {
    chars <- strsplit(as.character(txt), "", useBytes = TRUE)[[1L]]
    if (!length(chars)) return(character())

    u <- unit(u)
    v <- unit(v)
    plane_normal <- unit(cross(u, v))

    char_w <- 1
    gap <- 0.35
    total_w <- length(chars) * char_w + max(0, length(chars) - 1L) * gap

    x0 <- -total_w / 2
    y0 <- -0.5

    out <- vector("list", 0L)

    for (ii in seq_along(chars)) {
      segs <- glyph(chars[ii])
      if (!length(segs)) next

      shift <- x0 + (ii - 1L) * (char_w + gap)

      for (sg in segs) {
        p1 <- center + u * ((shift + sg[1L]) * size) + v * ((y0 + sg[2L]) * size)
        p2 <- center + u * ((shift + sg[3L]) * size) + v * ((y0 + sg[4L]) * size)

        out[[length(out) + 1L]] <- if (isTRUE(flat_text)) {
          flat_stroke(p1, p2, width, plane_normal)
        } else {
          line_prism(p1, p2, width)
        }
      }
    }

    unlist(out, use.names = FALSE)
  }

  cube_facets <- function(center, size) {
    r <- size / 2

    v <- list(
      center + c(-r, -r, -r),
      center + c( r, -r, -r),
      center + c( r,  r, -r),
      center + c(-r,  r, -r),
      center + c(-r, -r,  r),
      center + c( r, -r,  r),
      center + c( r,  r,  r),
      center + c(-r,  r,  r)
    )

    faces <- list(
      c(1, 3, 2), c(1, 4, 3),
      c(5, 6, 7), c(5, 7, 8),
      c(1, 2, 6), c(1, 6, 5),
      c(2, 3, 7), c(2, 7, 6),
      c(3, 4, 8), c(3, 8, 7),
      c(4, 1, 5), c(4, 5, 8)
    )

    unlist(
      lapply(faces, function(ii) {
        facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
      }),
      use.names = FALSE
    )
  }

  sphere_facets <- function(center, diameter) {
    radius <- diameter / 2
    segments <- max(6L, as.integer(round(5 + 19 * resolution)))
    rings <- max(4L, as.integer(round(3 + 9 * resolution)))

    theta <- seq(0, pi, length.out = rings + 1L)
    phi <- seq(0, 2 * pi, length.out = segments + 1L)

    point <- function(i, j) {
      center + radius * c(
        sin(theta[i]) * cos(phi[j]),
        sin(theta[i]) * sin(phi[j]),
        cos(theta[i])
      )
    }

    out <- vector("list", 0L)

    for (i in seq_len(rings)) {
      for (j in seq_len(segments)) {
        p00 <- point(i, j)
        p01 <- point(i, j + 1L)
        p10 <- point(i + 1L, j)
        p11 <- point(i + 1L, j + 1L)

        if (i == 1L) {
          out[[length(out) + 1L]] <- facet(p00, p10, p11)
        } else if (i == rings) {
          out[[length(out) + 1L]] <- facet(p00, p10, p01)
        } else {
          out[[length(out) + 1L]] <- c(
            facet(p00, p10, p11),
            facet(p00, p11, p01)
          )
        }
      }
    }

    unlist(out, use.names = FALSE)
  }

  tick_values <- function(lim) {
    if (n_ticks == 0L) return(numeric())
    vals <- pretty(lim, n = n_ticks)
    eps <- sqrt(.Machine$double.eps) * max(1, diff(lim))
    vals[vals >= lim[1L] - eps & vals <= lim[2L] + eps]
  }

  tick_labels <- function(vals) {
    vals[abs(vals) < sqrt(.Machine$double.eps)] <- 0
    out <- format(signif(vals, label_digits), trim = TRUE,
                  scientific = FALSE, drop0trailing = TRUE)
    out <- gsub(" ", "", out, fixed = TRUE)
    out[out == "-0"] <- "0"
    out
  }

  pieces <- list()

  add <- function(lines) {
    if (length(lines)) {
      pieces[[length(pieces) + 1L]] <<- lines
    }
  }

  point_fun <- switch(
    shape,
    cube = function(center) cube_facets(center, point_size),
    sphere = function(center) sphere_facets(center, point_size)
  )

  for (i in seq_along(x)) {
    add(point_fun(c(x[i], y[i], z[i])))
  }

  if (isTRUE(axes)) {
    pad <- 0.08 * span
    axis_w <- 0.006 * span
    tick_len <- 0.035 * span
    label_size <- 0.045 * span
    label_w <- 0.45 * axis_w

    anchor <- c(xlim[1L] - pad, ylim[1L] - pad, zlim[1L])

    x_axis_end <- c(xlim[2L] + pad / 2, anchor[2L], anchor[3L])
    y_axis_end <- c(anchor[1L], ylim[2L] + pad / 2, anchor[3L])
    z_axis_end <- c(anchor[1L], anchor[2L], zlim[2L] + pad / 2)

    add(line_prism(anchor, x_axis_end, axis_w))
    add(line_prism(anchor, y_axis_end, axis_w))
    add(line_prism(anchor, z_axis_end, axis_w))

    xt <- tick_values(xlim)
    yt <- tick_values(ylim)
    zt <- tick_values(zlim)

    xl <- tick_labels(xt)
    yl <- tick_labels(yt)
    zl <- tick_labels(zt)

    for (k in seq_along(xt)) {
      xx <- xt[k]

      add(line_prism(
        c(xx, anchor[2L] - tick_len / 2, anchor[3L]),
        c(xx, anchor[2L] + tick_len / 2, anchor[3L]),
        axis_w * 0.75
      ))

      add(draw_text(
        xl[k],
        c(xx, anchor[2L] - 1.7 * label_size, anchor[3L]),
        u = c(1, 0, 0),
        v = c(0, 1, 0),
        size = label_size,
        width = label_w
      ))
    }

    for (k in seq_along(yt)) {
      yy <- yt[k]

      add(line_prism(
        c(anchor[1L] - tick_len / 2, yy, anchor[3L]),
        c(anchor[1L] + tick_len / 2, yy, anchor[3L]),
        axis_w * 0.75
      ))

      add(draw_text(
        yl[k],
        c(anchor[1L] - 1.7 * label_size, yy, anchor[3L]),
        u = c(0, 1, 0),
        v = c(-1, 0, 0),
        size = label_size,
        width = label_w
      ))
    }

    z_label_offset <- (max(nchar(zl), 1L) * 0.45 + 1.9) * label_size

    for (k in seq_along(zt)) {
      zz <- zt[k]

      add(line_prism(
        c(anchor[1L] - tick_len / 2, anchor[2L], zz),
        c(anchor[1L] + tick_len / 2, anchor[2L], zz),
        axis_w * 0.75
      ))

      add(draw_text(
        zl[k],
        c(anchor[1L] - z_label_offset, anchor[2L], zz),
        u = c(1, 0, 0),
        v = c(0, 0, 1),
        size = label_size,
        width = label_w
      ))
    }

    add(draw_text(
      "X",
      c(xlim[2L] + 1.8 * label_size, anchor[2L], anchor[3L]),
      u = c(1, 0, 0),
      v = c(0, 1, 0),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))

    add(draw_text(
      "Y",
      c(anchor[1L], ylim[2L] + 1.8 * label_size, anchor[3L]),
      u = c(0, 1, 0),
      v = c(-1, 0, 0),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))

    add(draw_text(
      "Z",
      c(anchor[1L], anchor[2L], zlim[2L] + 1.8 * label_size),
      u = c(1, 0, 0),
      v = c(0, 0, 1),
      size = 1.25 * label_size,
      width = 1.2 * label_w
    ))
  }

  body <- unlist(pieces, use.names = FALSE)
  name <- gsub("[^A-Za-z0-9_.-]", "_", name)

  out <- paste(
    c(sprintf("solid %s", name), body, sprintf("endsolid %s", name)),
    collapse = "\n"
  )

  if (isTRUE(markdown)) {
    out <- paste0("```stl\n", out, "\n```")
  }

  out
 }
	#' Convert a bivariate function surface to ASCII STL
	#'
	#' Create an ASCII STL representation of a triangulated surface defined by a
	#' function `f(x, y)`. The output can optionally be wrapped in a GitHub Markdown
	#' `stl` code block for interactive 3D rendering.
	#'
	#' The surface is sampled over rectangular `x` and `y` ranges. The resulting
	#' `z` values are clipped to `zlim`, if supplied, and triangulated. Optional
	#' axes, tick marks, and numeric tick labels are added as STL geometry because
	#' STL has no native support for text or annotations.
	#'
	#' @param f A function taking numeric `x` and `y` values and returning a numeric
	#' `z` value. By default, `f` is called point-by-point. If `vectorized = TRUE`,
	#' `f` must accept equal-length numeric vectors `x` and `y` and return one
	#' numeric `z` value per pair.
	#' @param xlim Numeric vector of length 2 giving the lower and upper limits for
	#' the `x` axis.
	#' @param ylim Numeric vector of length 2 giving the lower and upper limits for
	#' the `y` axis.
	#' @param zlim Optional numeric vector of length 2 giving the lower and upper
	#' limits for the `z` axis. If `NULL`, the range of finite sampled `z` values
	#' is used.
	#' @param resolution Numeric value between 0.1 and 1 controlling mesh detail.
	#' Larger values create a finer triangulation and a larger STL string.
	#' @param name Character string used as the STL solid name. Characters other
	#' than letters, numbers, underscores, dots, and hyphens are replaced with
	#' underscores.
	#' @param n_ticks Integer giving the approximate number of tick marks to draw on
	#' each axis. Tick positions are chosen with [pretty()].
	#' @param solid Logical. If `TRUE`, add side walls and a gridded bottom face at
	#' `zlim[1]`. If `FALSE`, only the function surface is drawn.
	#' @param axes Logical. If `TRUE`, add 3D axes, tick marks, axis labels, and
	#' numeric tick labels as STL geometry.
	#' @param markdown Logical. If `TRUE`, wrap the ASCII STL in a Markdown fenced
	#' code block using the `stl` language identifier.
	#' @param vectorized Logical. If `TRUE`, call `f` once with vector inputs created
	#' from all sampled `x, y` pairs. If `FALSE`, call `f` point-by-point.
	#' @param digits Integer giving the number of significant digits used when
	#' formatting STL vertex coordinates and facet normals.
	#' @param label_digits Integer giving the number of significant digits used when
	#' formatting numeric axis tick labels.
	#'
	#' @return A character scalar containing ASCII STL code. If `markdown = TRUE`,
	#' the returned string is wrapped in a Markdown fenced `stl` block.
	#'
	#' @details
	#' This function produces a visual STL mesh rather than a mathematically exact
	#' surface. Axis text is approximated using stroke geometry, so numbers and
	#' letters are part of the mesh itself. Large values of `resolution` can create
	#' very large strings, especially when `solid = TRUE`.
	functionToSTL <- function(f, xlim, ylim, zlim = NULL,
	resolution = 0.5,
	name = "function_surface",
	n_ticks = 5L,
	solid = TRUE,
	axes = TRUE,
	markdown = FALSE,
	vectorized = FALSE,
	digits = 7L,
	label_digits = 4L) {
	stopifnot(is.function(f), length(xlim) == 2L, length(ylim) == 2L)

	resolution <- as.numeric(resolution)[1L]
	if (!is.finite(resolution) \|\| resolution < 0.1 \|\| resolution > 1) {
	stop("resolution must be a number between 0.1 and 1.")
	}

	xlim <- sort(as.numeric(xlim))
	ylim <- sort(as.numeric(ylim))
	if (diff(xlim) <= 0 \|\| diff(ylim) <= 0) {
	stop("xlim and ylim must each contain two distinct values.")
	}

	n <- as.integer(round(10 + 90 * resolution))
	nx <- ny <- n

	x <- seq(xlim[1L], xlim[2L], length.out = nx)
	y <- seq(ylim[1L], ylim[2L], length.out = ny)

	if (isTRUE(vectorized)) {
	xy <- expand.grid(x = x, y = y)
	zz <- f(xy$x, xy$y)
	if (length(zz) != nrow(xy)) {
	stop("When vectorized = TRUE, f(x, y) must return one z value per x/y pair.")
	}
	z <- matrix(as.numeric(zz), nrow = nx, ncol = ny)
	} else {
	z <- outer(x, y, Vectorize(function(a, b) as.numeric(f(a, b))[1L]))
	}

	z[!is.finite(z)] <- NA_real_

	if (is.null(zlim)) {
	if (!any(is.finite(z))) stop("f returned no finite z values.")
	zlim <- range(z, na.rm = TRUE)
	} else {
	zlim <- sort(as.numeric(zlim))
	}

	if (length(zlim) != 2L \|\| any(!is.finite(zlim)) \|\| diff(zlim) <= 0) {
	stop("zlim must contain two distinct finite values.")
	}

	z <- pmin(pmax(z, zlim[1L]), zlim[2L])

	cross <- function(a, b) {
	c(
	a[2L] * b[3L] - a[3L] * b[2L],
	a[3L] * b[1L] - a[1L] * b[3L],
	a[1L] * b[2L] - a[2L] * b[1L]
	)
	}

	unit <- function(a) {
	s <- sqrt(sum(a * a))
	if (!is.finite(s) \|\| s == 0) c(0, 0, 0) else a / s
	}

	fmt <- function(v) {
	formatC(v, format = "fg", digits = digits, flag = "#")
	}

	vertex_line <- function(p) {
	sprintf(" vertex %s %s %s", fmt(p[1L]), fmt(p[2L]), fmt(p[3L]))
	}

	facet <- function(a, b, c) {
	if (!all(is.finite(c(a, b, c)))) return(character())

	nr <- cross(b - a, c - a)
	ns <- sqrt(sum(nr * nr))
	if (!is.finite(ns) \|\| ns == 0) return(character())

	n <- nr / ns

	c(
	sprintf(" facet normal %s %s %s", fmt(n[1L]), fmt(n[2L]), fmt(n[3L])),
	" outer loop",
	vertex_line(a),
	vertex_line(b),
	vertex_line(c),
	" endloop",
	" endfacet"
	)
	}

	line_prism <- function(a, b, width) {
	if (!is.finite(width) \|\| width <= 0) return(character())

	d <- b - a
	len <- sqrt(sum(d * d))
	if (!is.finite(len) \|\| len == 0) return(character())

	e <- d / len
	ref <- if (abs(e[3L]) < 0.85) c(0, 0, 1) else c(0, 1, 0)

	n1 <- unit(cross(e, ref))
	if (all(n1 == 0)) n1 <- c(1, 0, 0)
	n2 <- unit(cross(e, n1))

	r <- width / 2

	v <- list(
	a + n1 * r + n2 * r,
	a - n1 * r + n2 * r,
	a - n1 * r - n2 * r,
	a + n1 * r - n2 * r,
	b + n1 * r + n2 * r,
	b - n1 * r + n2 * r,
	b - n1 * r - n2 * r,
	b + n1 * r - n2 * r
	)

	faces <- list(
	c(1, 2, 6), c(1, 6, 5),
	c(2, 3, 7), c(2, 7, 6),
	c(3, 4, 8), c(3, 8, 7),
	c(4, 1, 5), c(4, 5, 8),
	c(1, 4, 3), c(1, 3, 2),
	c(5, 6, 7), c(5, 7, 8)
	)

	unlist(
	lapply(faces, function(ii) {
	facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
	}),
	use.names = FALSE
	)
	}

	glyph <- local({
	seg <- function(x1, y1, x2, y2) {
	list(c(x1, y1, x2, y2))
	}

	arc <- function(cx, cy, rx, ry, a0, a1, n = 10L) {
	th <- seq(a0, a1, length.out = n + 1L)
	p <- cbind(cx + rx * cos(th), cy + ry * sin(th))

	lapply(seq_len(nrow(p) - 1L), function(i) {
	c(p[i, 1L], p[i, 2L], p[i + 1L, 1L], p[i + 1L, 2L])
	})
	}

	join <- function(...) {
	unlist(list(...), recursive = FALSE, use.names = FALSE)
	}

	function(ch) {
	switch(toupper(ch),
	"0" = arc(0.5, 0.5, 0.42, 0.50, 0, 2 * pi, 18L),

	"1" = join(
	seg(0.35, 0.82, 0.55, 1.00),
	seg(0.55, 1.00, 0.55, 0.00),
	seg(0.35, 0.00, 0.75, 0.00)
	),

	"2" = join(
	arc(0.50, 0.72, 0.40, 0.28, pi, 0, 10L),
	seg(0.90, 0.72, 0.12, 0.00),
	seg(0.12, 0.00, 0.90, 0.00)
	),

	"3" = join(
	seg(0.22, 1.00, 0.42, 1.00),
	arc(0.42, 0.75, 0.42, 0.25, pi / 2, -pi / 2, 10L),
	arc(0.42, 0.25, 0.42, 0.25, pi / 2, -pi / 2, 10L),
	seg(0.22, 0.00, 0.42, 0.00)
	),

	"4" = join(
	seg(0.78, 1.00, 0.78, 0.00),
	seg(0.18, 0.55, 0.90, 0.55),
	seg(0.18, 0.55, 0.72, 1.00)
	),

	"5" = join(
	seg(0.88, 1.00, 0.20, 1.00),
	seg(0.20, 1.00, 0.20, 0.56),
	seg(0.20, 0.56, 0.50, 0.56),
	arc(0.50, 0.28, 0.38, 0.28, pi / 2, -pi / 2, 10L),
	seg(0.50, 0.00, 0.22, 0.00)
	),

	"6" = join(
	seg(0.76, 0.96, 0.28, 0.58),
	arc(0.50, 0.32, 0.36, 0.32, 0, 2 * pi, 16L)
	),

	"7" = join(
	seg(0.12, 1.00, 0.90, 1.00),
	seg(0.90, 1.00, 0.35, 0.00)
	),

	"8" = join(
	arc(0.50, 0.73, 0.35, 0.27, 0, 2 * pi, 14L),
	arc(0.50, 0.27, 0.38, 0.27, 0, 2 * pi, 14L)
	),

	"9" = join(
	arc(0.50, 0.68, 0.36, 0.32, 0, 2 * pi, 16L),
	seg(0.78, 0.45, 0.34, 0.04)
	),

	"-" = seg(0.16, 0.50, 0.84, 0.50),

	"+" = join(
	seg(0.16, 0.50, 0.84, 0.50),
	seg(0.50, 0.16, 0.50, 0.84)
	),

	"." = arc(0.50, 0.06, 0.08, 0.08, 0, 2 * pi, 8L),

	"X" = join(
	seg(0.12, 0.00, 0.88, 1.00),
	seg(0.12, 1.00, 0.88, 0.00)
	),

	"Y" = join(
	seg(0.12, 1.00, 0.50, 0.52),
	seg(0.88, 1.00, 0.50, 0.52),
	seg(0.50, 0.52, 0.50, 0.00)
	),

	"Z" = join(
	seg(0.12, 1.00, 0.88, 1.00),
	seg(0.88, 1.00, 0.12, 0.00),
	seg(0.12, 0.00, 0.88, 0.00)
	),

	list()
	)
	}
	})

	flat_stroke <- function(a, b, width) {
	d <- b - a
	len <- sqrt(sum(d * d))
	if (!is.finite(len) \|\| len == 0) return(character())

	e <- d / len
	n <- unit(cross(e, c(0, 0, 1)))

	if (all(n == 0)) {
	n <- unit(cross(e, c(0, 1, 0)))
	}

	r <- width / 2

	p1 <- a + n * r
	p2 <- a - n * r
	p3 <- b - n * r
	p4 <- b + n * r

	c(
	facet(p1, p2, p3),
	facet(p1, p3, p4)
	)
	}

	draw_text <- function(txt, center, u, v, size, width,
	flat_text = TRUE) {
	chars <- strsplit(as.character(txt), "", fixed = TRUE)[[1L]]
	if (!length(chars)) return(character())

	u <- unit(u)
	v <- unit(v)

	char_w <- 1
	gap <- 0.35
	total_w <- length(chars) * char_w + max(0, length(chars) - 1L) * gap

	x0 <- -total_w / 2
	y0 <- -0.5

	out <- vector("list", 0L)
	stroke_fun <- if (isTRUE(flat_text)) flat_stroke else line_prism

	for (ii in seq_along(chars)) {
	segs <- glyph(chars[ii])
	if (!length(segs)) next

	shift <- x0 + (ii - 1L) * (char_w + gap)

	for (sg in segs) {
	p1 <- center + u * ((shift + sg[1L]) * size) + v * ((y0 + sg[2L]) * size)
	p2 <- center + u * ((shift + sg[3L]) * size) + v * ((y0 + sg[4L]) * size)
	out[[length(out) + 1L]] <- stroke_fun(p1, p2, width)
	}
	}

	unlist(out, use.names = FALSE)
	}

	tick_values <- function(lim) {
	vals <- pretty(lim, n = n_ticks)
	eps <- sqrt(.Machine$double.eps) * max(1, diff(lim))
	vals[vals >= lim[1L] - eps & vals <= lim[2L] + eps]
	}

	tick_labels <- function(vals) {
	vals[abs(vals) < sqrt(.Machine$double.eps)] <- 0
	out <- format(signif(vals, label_digits), trim = TRUE,
	scientific = FALSE, drop0trailing = TRUE)
	out <- gsub(" ", "", out, fixed = TRUE)
	out[out == "-0"] <- "0"
	out
	}

	pieces <- list()

	add <- function(lines) {
	if (length(lines)) {
	pieces[[length(pieces) + 1L]] <<- lines
	}
	}

	P <- function(i, j) c(x[i], y[j], z[i, j])
	Q <- function(i, j) c(x[i], y[j], zlim[1L])

	for (i in seq_len(nx - 1L)) {
	for (j in seq_len(ny - 1L)) {
	p00 <- P(i, j)
	p10 <- P(i + 1L, j)
	p01 <- P(i, j + 1L)
	p11 <- P(i + 1L, j + 1L)

	add(facet(p00, p10, p11))
	add(facet(p00, p11, p01))
	}
	}

	if (isTRUE(solid)) {
	for (i in seq_len(nx - 1L)) {
	p0 <- P(i, 1L)
	p1 <- P(i + 1L, 1L)
	q0 <- Q(i, 1L)
	q1 <- Q(i + 1L, 1L)

	add(facet(p0, q0, q1))
	add(facet(p0, q1, p1))

	p0 <- P(i, ny)
	p1 <- P(i + 1L, ny)
	q0 <- Q(i, ny)
	q1 <- Q(i + 1L, ny)

	add(facet(p0, q1, q0))
	add(facet(p0, p1, q1))
	}

	for (j in seq_len(ny - 1L)) {
	p0 <- P(1L, j)
	p1 <- P(1L, j + 1L)
	q0 <- Q(1L, j)
	q1 <- Q(1L, j + 1L)

	add(facet(p0, q1, q0))
	add(facet(p0, p1, q1))

	p0 <- P(nx, j)
	p1 <- P(nx, j + 1L)
	q0 <- Q(nx, j)
	q1 <- Q(nx, j + 1L)

	add(facet(p0, q0, q1))
	add(facet(p0, q1, p1))
	}

	for (i in seq_len(nx - 1L)) {
	for (j in seq_len(ny - 1L)) {
	q00 <- c(x[i], y[j], zlim[1L])
	q10 <- c(x[i + 1L], y[j], zlim[1L])
	q01 <- c(x[i], y[j + 1L], zlim[1L])
	q11 <- c(x[i + 1L], y[j + 1L], zlim[1L])

	add(facet(q00, q11, q10))
	add(facet(q00, q01, q11))
	}
	}
	}

	if (isTRUE(axes)) {
	span <- max(diff(xlim), diff(ylim), diff(zlim))
	pad <- 0.08 * span
	axis_w <- 0.006 * span
	tick_len <- 0.035 * span
	label_size <- 0.045 * span
	label_w <- 0.45 * axis_w

	anchor <- c(xlim[1L] - pad, ylim[1L] - pad, zlim[1L])

	x_axis_end <- c(xlim[2L] + pad / 2, anchor[2L], zlim[1L])
	y_axis_end <- c(anchor[1L], ylim[2L] + pad / 2, zlim[1L])
	z_axis_end <- c(anchor[1L], anchor[2L], zlim[2L] + pad / 2)

	add(line_prism(anchor, x_axis_end, axis_w))
	add(line_prism(anchor, y_axis_end, axis_w))
	add(line_prism(anchor, z_axis_end, axis_w))

	xt <- tick_values(xlim)
	yt <- tick_values(ylim)
	zt <- tick_values(zlim)

	xl <- tick_labels(xt)
	yl <- tick_labels(yt)
	zl <- tick_labels(zt)

	for (k in seq_along(xt)) {
	xx <- xt[k]

	add(line_prism(
	c(xx, anchor[2L] - tick_len / 2, zlim[1L]),
	c(xx, anchor[2L] + tick_len / 2, zlim[1L]),
	axis_w * 0.75
	))

	add(draw_text(
	xl[k],
	c(xx, anchor[2L] - 1.7 * label_size, zlim[1L]),
	u = c(1, 0, 0),
	v = c(0, 1, 0),
	size = label_size,
	width = label_w
	))
	}

	for (k in seq_along(yt)) {
	yy <- yt[k]

	add(line_prism(
	c(anchor[1L] - tick_len / 2, yy, zlim[1L]),
	c(anchor[1L] + tick_len / 2, yy, zlim[1L]),
	axis_w * 0.75
	))

	add(draw_text(
	yl[k],
	c(anchor[1L] - 1.7 * label_size, yy, zlim[1L]),
	u = c(0, 1, 0),
	v = c(-1, 0, 0),
	size = label_size,
	width = label_w
	))
	}

	z_label_offset <- (max(nchar(zl), 1L) * 0.45 + 1.9) * label_size

	for (k in seq_along(zt)) {
	zz <- zt[k]

	add(line_prism(
	c(anchor[1L] - tick_len / 2, anchor[2L], zz),
	c(anchor[1L] + tick_len / 2, anchor[2L], zz),
	axis_w * 0.75
	))

	add(draw_text(
	zl[k],
	c(anchor[1L] - z_label_offset, anchor[2L], zz),
	u = c(1, 0, 0),
	v = c(0, 0, 1),
	size = label_size,
	width = label_w
	))
	}

	add(draw_text(
	"X",
	c(xlim[2L] + 1.8 * label_size, anchor[2L], zlim[1L]),
	u = c(1, 0, 0),
	v = c(0, 1, 0),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))

	add(draw_text(
	"Y",
	c(anchor[1L], ylim[2L] + 1.8 * label_size, zlim[1L]),
	u = c(0, 1, 0),
	v = c(-1, 0, 0),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))

	add(draw_text(
	"Z",
	c(anchor[1L], anchor[2L], zlim[2L] + 1.8 * label_size),
	u = c(1, 0, 0),
	v = c(0, 0, 1),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))
	}

	body <- unlist(pieces, use.names = FALSE)
	name <- gsub("[^A-Za-z0-9_.-]", "_", name)

	out <- paste(
	c(sprintf("solid %s", name), body, sprintf("endsolid %s", name)),
	collapse = "\n"
	)

	if (isTRUE(markdown)) {
	out <- paste0("```stl\n", out, "\n```")
	}

	out
	}
	#' Convert 3D points to ASCII STL
	#'
	#' Create an ASCII STL representation of a 3D point cloud. Each point is drawn as
	#' a small cube or sphere, and optional axes, tick marks, and numeric tick labels
	#' are added as STL geometry.
	#'
	#' The output can optionally be wrapped in a GitHub Markdown `stl` code block for
	#' interactive 3D rendering.
	#'
	#' @param x Numeric vector of `x` coordinates.
	#' @param y Numeric vector of `y` coordinates. Must have the same length as `x`.
	#' @param z Numeric vector of `z` coordinates. Must have the same length as `x`.
	#' @param shape Character string specifying how points are drawn. One of
	#' `"cube"` or `"sphere"`.
	#' @param xlim Optional numeric vector of length 2 giving the lower and upper
	#' limits for the `x` axis. If `NULL`, the finite range of `x` is used.
	#' @param ylim Optional numeric vector of length 2 giving the lower and upper
	#' limits for the `y` axis. If `NULL`, the finite range of `y` is used.
	#' @param zlim Optional numeric vector of length 2 giving the lower and upper
	#' limits for the `z` axis. If `NULL`, the finite range of `z` is used.
	#' @param resolution Numeric value between 0.1 and 1 controlling marker detail.
	#' For spheres, larger values create smoother spheres. If `point_size = NULL`,
	#' larger values may also be used to create smaller default markers, depending
	#' on the implementation.
	#' @param point_size Optional positive numeric value giving the marker diameter
	#' in data units. If `NULL`, a size is chosen automatically from the largest
	#' axis span and `resolution`.
	#' @param name Character string used as the STL solid name. Characters other
	#' than letters, numbers, underscores, dots, and hyphens are replaced with
	#' underscores.
	#' @param n_ticks Integer giving the approximate number of tick marks to draw on
	#' each axis. Tick positions are chosen with [pretty()]. Use `0` to suppress
	#' tick marks and tick labels while keeping axes.
	#' @param axes Logical. If `TRUE`, add 3D axes, tick marks, axis labels, and
	#' numeric tick labels as STL geometry.
	#' @param markdown Logical. If `TRUE`, wrap the ASCII STL in a Markdown fenced
	#' code block using the `stl` language identifier.
	#' @param digits Integer giving the number of significant digits used when
	#' formatting STL vertex coordinates and facet normals.
	#' @param label_digits Integer giving the number of significant digits used when
	#' formatting numeric axis tick labels.
	#' @param flat_text Logical. If `TRUE`, draw axis labels and tick labels as flat
	#' ribbon strokes. If `FALSE`, draw text strokes as small rectangular prisms.
	#' Flat text usually looks cleaner in STL wireframe renderers.
	#' @param clip Logical. If `TRUE`, drop points outside `xlim`, `ylim`, and
	#' `zlim`. If `FALSE`, all finite points are rendered, even if they lie outside
	#' the displayed axis limits.
	#'
	#' @return A character scalar containing ASCII STL code. If `markdown = TRUE`,
	#' the returned string is wrapped in a Markdown fenced `stl` block.
	#'
	#' @details
	#' STL has no native point, text, axis, or tick-mark primitives, so all elements
	#' are represented as triangular mesh facets. Cubes are much lighter than
	#' spheres. Spheres can produce large STL strings, especially with many points
	#' and high `resolution`.
	pointsToSTL <- function(x, y, z,
	shape = c("cube", "sphere"),
	xlim = NULL, ylim = NULL, zlim = NULL,
	resolution = 0.5,
	point_size = NULL,
	name = "point_cloud",
	n_ticks = 5L,
	axes = TRUE,
	markdown = FALSE,
	digits = 7L,
	label_digits = 4L,
	flat_text = TRUE,
	clip = TRUE) {
	shape <- match.arg(shape)

	if (length(x) != length(y) \|\| length(x) != length(z)) {
	stop("x, y, and z must have equal lengths.")
	}

	if (!length(x)) stop("x, y, and z must not be empty.")

	resolution <- as.numeric(resolution)[1L]
	if (!is.finite(resolution) \|\| resolution < 0.1 \|\| resolution > 1) {
	stop("resolution must be a number between 0.1 and 1.")
	}

	x <- as.numeric(x)
	y <- as.numeric(y)
	z <- as.numeric(z)

	ok <- is.finite(x) & is.finite(y) & is.finite(z)
	if (!any(ok)) stop("x, y, and z contain no finite complete points.")

	if (!all(ok)) {
	warning("Dropping non-finite points.")
	x <- x[ok]
	y <- y[ok]
	z <- z[ok]
	}

	make_lim <- function(v, lim, nm) {
	if (is.null(lim)) {
	lim <- range(v)
	} else {
	if (length(lim) != 2L) stop(nm, " must contain two values.")
	lim <- sort(as.numeric(lim))
	}

	if (any(!is.finite(lim))) stop(nm, " must contain finite values.")

	if (diff(lim) == 0) {
	d <- max(1, abs(lim[1L])) * 0.05
	lim <- lim + c(-d, d)
	}

	lim
	}

	xlim <- make_lim(x, xlim, "xlim")
	ylim <- make_lim(y, ylim, "ylim")
	zlim <- make_lim(z, zlim, "zlim")

	if (isTRUE(clip)) {
	keep <- x >= xlim[1L] & x <= xlim[2L] &
	y >= ylim[1L] & y <= ylim[2L] &
	z >= zlim[1L] & z <= zlim[2L]

	if (!any(keep)) stop("No points remain after clipping to xlim, ylim, and zlim.")

	x <- x[keep]
	y <- y[keep]
	z <- z[keep]
	}

	span <- max(diff(xlim), diff(ylim), diff(zlim))
	if (!is.finite(span) \|\| span <= 0) span <- 1

	if (is.null(point_size)) {
	# Higher resolution means smaller markers.
	# resolution = 0.1 -> diameter about 5.5% of the largest axis span
	# resolution = 1.0 -> diameter about 1.2% of the largest axis span
	point_size <- span * (0.012 + (1 - resolution) / 0.9 * (0.055 - 0.012))
	} else {
	point_size <- as.numeric(point_size)[1L]
	}

	if (!is.finite(point_size) \|\| point_size <= 0) {
	stop("point_size must be a positive finite number.")
	}

	n_ticks <- as.integer(n_ticks)[1L]
	if (!is.finite(n_ticks) \|\| n_ticks < 0L) {
	stop("n_ticks must be a non-negative integer.")
	}

	cross <- function(a, b) {
	c(
	a[2L] * b[3L] - a[3L] * b[2L],
	a[3L] * b[1L] - a[1L] * b[3L],
	a[1L] * b[2L] - a[2L] * b[1L]
	)
	}

	unit <- function(a) {
	s <- sqrt(sum(a * a))
	if (!is.finite(s) \|\| s == 0) c(0, 0, 0) else a / s
	}

	fmt <- function(v) {
	formatC(v, format = "fg", digits = digits, flag = "#")
	}

	vertex_line <- function(p) {
	sprintf(" vertex %s %s %s", fmt(p[1L]), fmt(p[2L]), fmt(p[3L]))
	}

	facet <- function(a, b, d) {
	if (!all(is.finite(c(a, b, d)))) return(character())

	nr <- cross(b - a, d - a)
	ns <- sqrt(sum(nr * nr))
	if (!is.finite(ns) \|\| ns == 0) return(character())

	n <- nr / ns

	c(
	sprintf(" facet normal %s %s %s", fmt(n[1L]), fmt(n[2L]), fmt(n[3L])),
	" outer loop",
	vertex_line(a),
	vertex_line(b),
	vertex_line(d),
	" endloop",
	" endfacet"
	)
	}

	line_prism <- function(a, b, width) {
	if (!is.finite(width) \|\| width <= 0) return(character())

	d <- b - a
	len <- sqrt(sum(d * d))
	if (!is.finite(len) \|\| len == 0) return(character())

	e <- d / len
	ref <- if (abs(e[3L]) < 0.85) c(0, 0, 1) else c(0, 1, 0)

	n1 <- unit(cross(e, ref))
	if (all(n1 == 0)) n1 <- c(1, 0, 0)
	n2 <- unit(cross(e, n1))

	r <- width / 2

	v <- list(
	a + n1 * r + n2 * r,
	a - n1 * r + n2 * r,
	a - n1 * r - n2 * r,
	a + n1 * r - n2 * r,
	b + n1 * r + n2 * r,
	b - n1 * r + n2 * r,
	b - n1 * r - n2 * r,
	b + n1 * r - n2 * r
	)

	faces <- list(
	c(1, 2, 6), c(1, 6, 5),
	c(2, 3, 7), c(2, 7, 6),
	c(3, 4, 8), c(3, 8, 7),
	c(4, 1, 5), c(4, 5, 8),
	c(1, 4, 3), c(1, 3, 2),
	c(5, 6, 7), c(5, 7, 8)
	)

	unlist(
	lapply(faces, function(ii) {
	facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
	}),
	use.names = FALSE
	)
	}

	flat_stroke <- function(a, b, width, plane_normal) {
	d <- b - a
	len <- sqrt(sum(d * d))
	if (!is.finite(len) \|\| len == 0) return(character())

	e <- d / len
	plane_normal <- unit(plane_normal)
	if (all(plane_normal == 0)) plane_normal <- c(0, 0, 1)

	n <- unit(cross(plane_normal, e))
	if (all(n == 0)) n <- unit(cross(c(0, 0, 1), e))
	if (all(n == 0)) n <- c(1, 0, 0)

	r <- width / 2

	p1 <- a + n * r
	p2 <- a - n * r
	p3 <- b - n * r
	p4 <- b + n * r

	c(
	facet(p1, p2, p3),
	facet(p1, p3, p4)
	)
	}

	glyph <- local({
	arc_steps <- max(6L, as.integer(round(5 + 15 * resolution)))

	seg <- function(x1, y1, x2, y2) {
	list(c(x1, y1, x2, y2))
	}

	arc <- function(cx, cy, rx, ry, a0, a1, n = arc_steps) {
	n <- max(1L, as.integer(n))
	th <- seq(a0, a1, length.out = n + 1L)
	p <- cbind(cx + rx * cos(th), cy + ry * sin(th))

	lapply(seq_len(nrow(p) - 1L), function(i) {
	c(p[i, 1L], p[i, 2L], p[i + 1L, 1L], p[i + 1L, 2L])
	})
	}

	join <- function(...) {
	unlist(list(...), recursive = FALSE, use.names = FALSE)
	}

	function(ch) {
	switch(toupper(ch),
	"0" = arc(0.5, 0.5, 0.42, 0.50, 0, 2 * pi, arc_steps),

	"1" = join(
	seg(0.35, 0.82, 0.55, 1.00),
	seg(0.55, 1.00, 0.55, 0.00),
	seg(0.35, 0.00, 0.75, 0.00)
	),

	"2" = join(
	arc(0.50, 0.72, 0.40, 0.28, pi, 0, arc_steps),
	seg(0.90, 0.72, 0.12, 0.00),
	seg(0.12, 0.00, 0.90, 0.00)
	),

	"3" = join(
	seg(0.22, 1.00, 0.42, 1.00),
	arc(0.42, 0.75, 0.42, 0.25, pi / 2, -pi / 2, arc_steps),
	arc(0.42, 0.25, 0.42, 0.25, pi / 2, -pi / 2, arc_steps),
	seg(0.22, 0.00, 0.42, 0.00)
	),

	"4" = join(
	seg(0.78, 1.00, 0.78, 0.00),
	seg(0.18, 0.55, 0.90, 0.55),
	seg(0.18, 0.55, 0.72, 1.00)
	),

	"5" = join(
	seg(0.88, 1.00, 0.20, 1.00),
	seg(0.20, 1.00, 0.20, 0.56),
	seg(0.20, 0.56, 0.50, 0.56),
	arc(0.50, 0.28, 0.38, 0.28, pi / 2, -pi / 2, arc_steps),
	seg(0.50, 0.00, 0.22, 0.00)
	),

	"6" = join(
	seg(0.76, 0.96, 0.28, 0.58),
	arc(0.50, 0.32, 0.36, 0.32, 0, 2 * pi, arc_steps)
	),

	"7" = join(
	seg(0.12, 1.00, 0.90, 1.00),
	seg(0.90, 1.00, 0.35, 0.00)
	),

	"8" = join(
	arc(0.50, 0.73, 0.35, 0.27, 0, 2 * pi, arc_steps),
	arc(0.50, 0.27, 0.38, 0.27, 0, 2 * pi, arc_steps)
	),

	"9" = join(
	arc(0.50, 0.68, 0.36, 0.32, 0, 2 * pi, arc_steps),
	seg(0.78, 0.45, 0.34, 0.04)
	),

	"-" = seg(0.16, 0.50, 0.84, 0.50),

	"+" = join(
	seg(0.16, 0.50, 0.84, 0.50),
	seg(0.50, 0.16, 0.50, 0.84)
	),

	"." = arc(0.50, 0.06, 0.08, 0.08, 0, 2 * pi, max(6L, arc_steps %/% 2L)),

	"X" = join(
	seg(0.12, 0.00, 0.88, 1.00),
	seg(0.12, 1.00, 0.88, 0.00)
	),

	"Y" = join(
	seg(0.12, 1.00, 0.50, 0.52),
	seg(0.88, 1.00, 0.50, 0.52),
	seg(0.50, 0.52, 0.50, 0.00)
	),

	"Z" = join(
	seg(0.12, 1.00, 0.88, 1.00),
	seg(0.88, 1.00, 0.12, 0.00),
	seg(0.12, 0.00, 0.88, 0.00)
	),

	list()
	)
	}
	})

	draw_text <- function(txt, center, u, v, size, width) {
	chars <- strsplit(as.character(txt), "", useBytes = TRUE)[[1L]]
	if (!length(chars)) return(character())

	u <- unit(u)
	v <- unit(v)
	plane_normal <- unit(cross(u, v))

	char_w <- 1
	gap <- 0.35
	total_w <- length(chars) * char_w + max(0, length(chars) - 1L) * gap

	x0 <- -total_w / 2
	y0 <- -0.5

	out <- vector("list", 0L)

	for (ii in seq_along(chars)) {
	segs <- glyph(chars[ii])
	if (!length(segs)) next

	shift <- x0 + (ii - 1L) * (char_w + gap)

	for (sg in segs) {
	p1 <- center + u * ((shift + sg[1L]) * size) + v * ((y0 + sg[2L]) * size)
	p2 <- center + u * ((shift + sg[3L]) * size) + v * ((y0 + sg[4L]) * size)

	out[[length(out) + 1L]] <- if (isTRUE(flat_text)) {
	flat_stroke(p1, p2, width, plane_normal)
	} else {
	line_prism(p1, p2, width)
	}
	}
	}

	unlist(out, use.names = FALSE)
	}

	cube_facets <- function(center, size) {
	r <- size / 2

	v <- list(
	center + c(-r, -r, -r),
	center + c( r, -r, -r),
	center + c( r, r, -r),
	center + c(-r, r, -r),
	center + c(-r, -r, r),
	center + c( r, -r, r),
	center + c( r, r, r),
	center + c(-r, r, r)
	)

	faces <- list(
	c(1, 3, 2), c(1, 4, 3),
	c(5, 6, 7), c(5, 7, 8),
	c(1, 2, 6), c(1, 6, 5),
	c(2, 3, 7), c(2, 7, 6),
	c(3, 4, 8), c(3, 8, 7),
	c(4, 1, 5), c(4, 5, 8)
	)

	unlist(
	lapply(faces, function(ii) {
	facet(v[[ii[1L]]], v[[ii[2L]]], v[[ii[3L]]])
	}),
	use.names = FALSE
	)
	}

	sphere_facets <- function(center, diameter) {
	radius <- diameter / 2
	segments <- max(6L, as.integer(round(5 + 19 * resolution)))
	rings <- max(4L, as.integer(round(3 + 9 * resolution)))

	theta <- seq(0, pi, length.out = rings + 1L)
	phi <- seq(0, 2 * pi, length.out = segments + 1L)

	point <- function(i, j) {
	center + radius * c(
	sin(theta[i]) * cos(phi[j]),
	sin(theta[i]) * sin(phi[j]),
	cos(theta[i])
	)
	}

	out <- vector("list", 0L)

	for (i in seq_len(rings)) {
	for (j in seq_len(segments)) {
	p00 <- point(i, j)
	p01 <- point(i, j + 1L)
	p10 <- point(i + 1L, j)
	p11 <- point(i + 1L, j + 1L)

	if (i == 1L) {
	out[[length(out) + 1L]] <- facet(p00, p10, p11)
	} else if (i == rings) {
	out[[length(out) + 1L]] <- facet(p00, p10, p01)
	} else {
	out[[length(out) + 1L]] <- c(
	facet(p00, p10, p11),
	facet(p00, p11, p01)
	)
	}
	}
	}

	unlist(out, use.names = FALSE)
	}

	tick_values <- function(lim) {
	if (n_ticks == 0L) return(numeric())
	vals <- pretty(lim, n = n_ticks)
	eps <- sqrt(.Machine$double.eps) * max(1, diff(lim))
	vals[vals >= lim[1L] - eps & vals <= lim[2L] + eps]
	}

	tick_labels <- function(vals) {
	vals[abs(vals) < sqrt(.Machine$double.eps)] <- 0
	out <- format(signif(vals, label_digits), trim = TRUE,
	scientific = FALSE, drop0trailing = TRUE)
	out <- gsub(" ", "", out, fixed = TRUE)
	out[out == "-0"] <- "0"
	out
	}

	pieces <- list()

	add <- function(lines) {
	if (length(lines)) {
	pieces[[length(pieces) + 1L]] <<- lines
	}
	}

	point_fun <- switch(
	shape,
	cube = function(center) cube_facets(center, point_size),
	sphere = function(center) sphere_facets(center, point_size)
	)

	for (i in seq_along(x)) {
	add(point_fun(c(x[i], y[i], z[i])))
	}

	if (isTRUE(axes)) {
	pad <- 0.08 * span
	axis_w <- 0.006 * span
	tick_len <- 0.035 * span
	label_size <- 0.045 * span
	label_w <- 0.45 * axis_w

	anchor <- c(xlim[1L] - pad, ylim[1L] - pad, zlim[1L])

	x_axis_end <- c(xlim[2L] + pad / 2, anchor[2L], anchor[3L])
	y_axis_end <- c(anchor[1L], ylim[2L] + pad / 2, anchor[3L])
	z_axis_end <- c(anchor[1L], anchor[2L], zlim[2L] + pad / 2)

	add(line_prism(anchor, x_axis_end, axis_w))
	add(line_prism(anchor, y_axis_end, axis_w))
	add(line_prism(anchor, z_axis_end, axis_w))

	xt <- tick_values(xlim)
	yt <- tick_values(ylim)
	zt <- tick_values(zlim)

	xl <- tick_labels(xt)
	yl <- tick_labels(yt)
	zl <- tick_labels(zt)

	for (k in seq_along(xt)) {
	xx <- xt[k]

	add(line_prism(
	c(xx, anchor[2L] - tick_len / 2, anchor[3L]),
	c(xx, anchor[2L] + tick_len / 2, anchor[3L]),
	axis_w * 0.75
	))

	add(draw_text(
	xl[k],
	c(xx, anchor[2L] - 1.7 * label_size, anchor[3L]),
	u = c(1, 0, 0),
	v = c(0, 1, 0),
	size = label_size,
	width = label_w
	))
	}

	for (k in seq_along(yt)) {
	yy <- yt[k]

	add(line_prism(
	c(anchor[1L] - tick_len / 2, yy, anchor[3L]),
	c(anchor[1L] + tick_len / 2, yy, anchor[3L]),
	axis_w * 0.75
	))

	add(draw_text(
	yl[k],
	c(anchor[1L] - 1.7 * label_size, yy, anchor[3L]),
	u = c(0, 1, 0),
	v = c(-1, 0, 0),
	size = label_size,
	width = label_w
	))
	}

	z_label_offset <- (max(nchar(zl), 1L) * 0.45 + 1.9) * label_size

	for (k in seq_along(zt)) {
	zz <- zt[k]

	add(line_prism(
	c(anchor[1L] - tick_len / 2, anchor[2L], zz),
	c(anchor[1L] + tick_len / 2, anchor[2L], zz),
	axis_w * 0.75
	))

	add(draw_text(
	zl[k],
	c(anchor[1L] - z_label_offset, anchor[2L], zz),
	u = c(1, 0, 0),
	v = c(0, 0, 1),
	size = label_size,
	width = label_w
	))
	}

	add(draw_text(
	"X",
	c(xlim[2L] + 1.8 * label_size, anchor[2L], anchor[3L]),
	u = c(1, 0, 0),
	v = c(0, 1, 0),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))

	add(draw_text(
	"Y",
	c(anchor[1L], ylim[2L] + 1.8 * label_size, anchor[3L]),
	u = c(0, 1, 0),
	v = c(-1, 0, 0),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))

	add(draw_text(
	"Z",
	c(anchor[1L], anchor[2L], zlim[2L] + 1.8 * label_size),
	u = c(1, 0, 0),
	v = c(0, 0, 1),
	size = 1.25 * label_size,
	width = 1.2 * label_w
	))
	}

	body <- unlist(pieces, use.names = FALSE)
	name <- gsub("[^A-Za-z0-9_.-]", "_", name)

	out <- paste(
	c(sprintf("solid %s", name), body, sprintf("endsolid %s", name)),
	collapse = "\n"
	)

	if (isTRUE(markdown)) {
	out <- paste0("```stl\n", out, "\n```")
	}

	out
	}