# centerline tracing and stroke thinning

example algorithm for converting and simplifying broad brushstroke outlines to thinner lines:
* given that separate brushstrokes are available as outlines
* extract the centerline (medial axis or skeleton) as polylines `(((x y) ...) ...)`
* remove short paths and identify and merge the longest path to eliminate unwanted skeleton branches
* apply curve fitting to smooth and simplify the polylines
* use centroidal derivative scaling to reduce gaps between thinned strokes and bring them closer together uniformly

# coffeescript examples
## extracting the centerline using [skeleton-tracing](https://github.com/LingDong-/skeleton-tracing/tree/master/js)
~~~
centerline_from_canvas = (canvas) ->
  extract_longest_path remove_short_paths(TraceSkeleton.fromCanvas(canvas).polylines, 40)
~~~

## curve fitting using [fit-curve](https://github.com/soswow/fit-curve) and converting to svg path commands
~~~
simplify_to_svg = (polylines) ->
  error = 400.0
  svg_paths = []
  for polyline in polylines
    beziers = fit_curve polyline, error
    path_data = ''
    continue unless beziers.length
    beziers = for a in beziers
      for [x, y] in a
        [x.toFixed(0), y.toFixed(0)]
    [x0, y0] = beziers[0][0]
    path_data += "M#{x0},#{y0}"
    for bezier in beziers
      [p0, p1, p2, p3] = bezier
      [x1, y1] = p1
      [x2, y2] = p2
      [x3, y3] = p3
      path_data += "C#{x1},#{y1},#{x2},#{y2},#{x3},#{y3}"
    svg_paths.push path_data
  svg_paths
~~~

## extracting the longest path
~~~
extract_longest_path = (polylines) ->
  if 2 > polylines.length
    return (if polylines.length then polylines[0] else polylines)
  # Step 1: Build the graph
  nodeMap = new Map() # Map from point string to node ID
  nodeId = 0
  adj = {} # Adjacency list
  pointKey = (pt) ->
    # Round coordinates to 6 decimal places
    x = pt[0].toFixed 6
    y = pt[1].toFixed 6
    "#{x},#{y}"
  getNodeId = (pt) ->
    key = pointKey pt
    unless nodeMap.has key
      nodeMap.set key, nodeId++
    nodeMap.get key
  for polyline in polylines
    for i in [0...polyline.length - 1]
      pt1 = polyline[i]
      pt2 = polyline[i + 1]
      id1 = getNodeId pt1
      id2 = getNodeId pt2
      adj[id1] ?= []
      adj[id2] ?= []
      adj[id1].push id2
      adj[id2].push id1 # Since it's undirected
  # Map from node ID to point coordinates
  idToPoint = {}
  nodeMap.forEach (id, key) ->
    [x, y] = key.split(',').map Number
    idToPoint[id] = [x, y]
  # Step 2: Find the longest path
  bfs = (start) ->
    visited = new Set()
    queue = [[start, 0]]
    farthestNode = start
    maxDistance = 0
    while queue.length > 0
      [node, dist] = queue.shift()
      visited.add node
      if dist > maxDistance
        maxDistance = dist
        farthestNode = node
      for neighbor in adj[node]
        unless visited.has neighbor
          queue.push [neighbor, dist + 1]
          visited.add neighbor
    {node: farthestNode, distance: maxDistance}
  # First BFS to find one end of the longest path
  firstBfs = bfs 0 # Start from any node, say node 0
  # Second BFS from the farthest node found in the first BFS
  secondBfs = bfs firstBfs.node
  # Reconstruct the path from firstBfs.node to secondBfs.node
  bfsWithParents = (start, target) ->
    visited = new Set()
    queue = [start]
    parent = {}
    visited.add start
    while queue.length > 0
      node = queue.shift()
      break if node == target
      for neighbor in adj[node]
        unless visited.has neighbor
          visited.add neighbor
          parent[neighbor] = node
          queue.push neighbor
    # Reconstruct path from target to start
    path = []
    currentNode = target
    while currentNode?
      path.push currentNode
      currentNode = parent[currentNode]
    path.reverse() # From start to target
  longestPathNodeIds = bfsWithParents firstBfs.node, secondBfs.node
  idToPoint[a] for a in longestPathNodeIds
~~~

## removing short paths
~~~
calculate_polyline_length = (polyline) ->
  length = 0
  for i in [0...polyline.length - 1]
    x1 = polyline[i][0]
    y1 = polyline[i][1]
    x2 = polyline[i + 1][0]
    y2 = polyline[i + 1][1]
    dx = x2 - x1
    dy = y2 - y1
    length += Math.sqrt dx * dx + dy * dy
  length

remove_short_paths = (polylines, limit) ->
  for a in polylines
    continue if limit > calculate_polyline_length a
    a
~~~

## centroidal derivative scaling
~~~
polyline_centroid = (polyline) ->
  n = polyline.length
  x_sum = 0
  y_sum = 0
  for [x, y] in polyline
    x_sum += x
    y_sum += y
  [x_sum / n, y_sum / n]

scale_by_centroids = (polylines, factor) ->
  centroids = (polyline_centroid a for a in polylines)
  scaled_centroids = ([x * factor, y * factor] for [x, y] in centroids)
  for polyline, i in polylines
    centroid = centroids[i]
    scaled_centroid = scaled_centroids[i]
    dx = scaled_centroid[0] - centroid[0]
    dy = scaled_centroid[1] - centroid[1]
    [x + dx, y + dy] for [x, y] in polyline
~~~

## drawing an svg path onto canvas using [svg-path-parser](https://github.com/convertSvg/svg-path-parse)
~~~
SVGPathParser = require "svg-path-parser"

canvas_context_draw_svg_commands = (ctx, commands) ->
  ctx.beginPath()
  for a in commands
    switch a.code
      when "M" then ctx.moveTo a.x, a.y
      when "L" then ctx.lineTo a.x, a.y
      when "Q" then ctx.quadraticCurveTo a.x1, a.y1, a.x, a.y
      when "C" then ctx.bezierCurveTo a.x1, a.y1, a.x2, a.y2, a.x, a.y
      when "Z" then ctx.closePath()
  ctx.fillStyle = "#fff"
  ctx.fill()

canvas_context_draw_svg_path = (ctx, path) -> canvas_context_draw_svg_commands ctx, SVGPathParser.parseSVG path

canvas = createCanvas canvas_width, canvas_height
ctx = canvas.getContext "2d"
canvas_context_draw_svg_path ctx, path
~~~