Line Fractal Point Computation

After reading some in Mandelbrot's The Fractal Geometry of Nature I decided to figure out how to create a simple tool to generate line fractals similar to those in the book. I tried several techniques from using a 2d canvas in HTML to using WebGL to render up to 8 million segments.

Here's how I compute the segments from a pattern iterated to a certain depth.

{pow} = Math

# In place concat, modify m1 in place
concat = (m1, m2) ->
  {a, b, c, d, tx, ty} = m1

  m1.a = a * m2.a + c * m2.b
  m1.b = b * m2.a + d * m2.b
  m1.c = a * m2.c + c * m2.d
  m1.d = b * m2.c + d * m2.d
  m1.tx = a * m2.tx + c * m2.ty + tx
  m1.ty = b * m2.tx + d * m2.ty + ty

  return

vecToMatrix = (p1, p2) ->
  a = d = p2.x - p1.x
  c = -(b = p2.y - p1.y)

  {a, b, c, d, tx: p1.x, ty: p1.y}

patternExpand = (pattern, depth=2, output) ->
  segmentCount = pattern.length - 1

  # Output float array of x,y point pairs
  outputLength = (pow(segmentCount, depth) + 1) * 2
  output ?= new Float32Array outputLength

  matrices = pattern.map (p, i) ->
    if i < segmentCount
      vecToMatrix(p, pattern[i+1])
    else
      vecToMatrix(p, pattern[0]) # :P

  m = {a:1,b:0,c:0,d:1,tx:0,ty:0}
  x = y = null

  i = 0
  while i < outputLength
    # set to identity
    m.a = m.d = 1
    m.c = m.b = m.tx = m.ty = 0

    d = depth-1
    k = (i / 2)|0
    # Extract the "address path" for the index
    # index / segmentCount^depth, index / segmentCount^(depth-1), index / segmentCount^(depth-2), ...
    while d > 0
      # TODO: Divisor will be product of pattern lengths when/if we have multiple patterns
      divisor = pow(segmentCount, d)
      n = (k / divisor)|0
      k = k % divisor
      # Math without object allocations
      concat(m, matrices[n])
      d--

    # m.transformPoint without allocations
    x = pattern[k].x
    y = pattern[k].y
    output[i] = x * m.a + y * m.c + m.tx
    output[i+1] = x * m.b + y * m.d + m.ty
    i += 2

  return output

The basic idea is to generate a list of matrix transforms that correspond to the list of points in the pattern. Then to use the index to find out which transforms to multiply in what order to compute the position. The function has been modified to be performant by not allocating any objects in the main loop. It will also be possible to parallelize because aside from setting up the initial n matrices all the other computations are independent.