Pen Plotting with Pluto

By Paul Butler (@paulgb) for PlutoCon2021

A plotter is a robot that draws based on instructions you provide it.

There are basically three varieties of plotters:

Vintage machines (like the HP 7470A) lovingly restored by hobbiests.
DIY machines, often of the polargraph variety.
Modern pre-built desktop pen plotters, like the AxiDraw.

At a low level, plotters are controlled by a series of commands, the most important being “move to the location (x, y)”, “raise the pen”, and “lower the pen”. In addition to these low-level commands, there are tools out there for converting various vector graphics files into the raw commands.

In this talk, I will use a tiny library I wrote called PenPlots.jl which provides basic data structures for representing points and paths, as well as generating an SVG (scalable vector graphics) representation of one or more paths. Pleasingly, SVG is supported by my plotter driver of choice, Saxi, and is also supported directly in the browser for easy previews in Pluto.

I'll start by generating a few basic shapes to show you how PenPlots.jl works, and then I'll show how Pluto makes it easy to explore the parameter space of more intricate plots.

xxxxxxxxxx

28.6 μs

Some Basics

Drawing a line

PenPlots.jl uses a Point data structure to represent an (x, y) coordinate on the plotting surface. These are screen coordinates rather than Cartesian coordinates, so the origin is in the upper-left hand corner and the y coordinate increases as we go down the plotting surface.

The Path data structure represents a path through two or more points. Its constructor takes a list of points.

Many types in PenPlots.jl, including Path, implement Base.show for text/html, so that the Pluto notebook will automatically preview them when they are the return value of a cell.

xxxxxxxxxx

13.5 μs

xxxxxxxxxx
 
Path([
    Point(1, 0),
    Point(0, 1),
])

4.2 μs

Drawing Two Lines

Multiple separate paths can be combined by putting them in a vector. In pen plotter terms, separate paths mean that the pen plotter will draw one path, then lift the pen, move to the beginning of the next path, and then draw it.

I'm avoiding using the terms “first path” and “second path” here. That's because even though we are using an ordered data structure (a vector), the paths will usually be reordered to reduce drawing time. This used to be a manual step (I wrote about it here), but recently it's become a built-in feature of the driver, as is the case with Saxi.

xxxxxxxxxx

8.7 μs

xxxxxxxxxx
 
[
    Path([
        Point(1, 0),
        Point(0, 1),
    ]),
    Path([
        Point(-0.5, 0.),
        Point(1, 0.5),
    ]),
]

6.8 μs

Drawing Many Lines

This technique plays really nicely with the map function to draw a bunch of lines.

xxxxxxxxxx

5.6 μs

xxxxxxxxxx
 
map(0:40) do i
    Path([
            Point(i / 2, 0),
            Point(i, 100 - i * 2)
            ])
end

214 ms

Drawing A Circle

map can also be used within a Path to create a vector of points to visit. For example, we can create a circle by rotating the unit vector around the origin in tiny increments.

frac_rotation is a helper function to generate a rotation matrix from a fraction, where 1.0 represents a full rotation. PenPlots.jl provides it along with degree_rotation and radian_rotation.

Rotation matrices can be multiplied by Point, Path, and Vector{Path}. Rotations are always about the origin (0, 0).

xxxxxxxxxx

6.9 μs

xxxxxxxxxx
 
Path(map(0:0.01:1) do i
    frac_rotation(i) * unitvec
end)

13.8 μs

Spiral

xxxxxxxxxx

3.4 μs

xxxxxxxxxx
 
Path(map(0:0.01:30) do i
    frac_rotation(i) * unitvec * i
end)

154 μs

Step Size:

xxxxxxxxxx
 
md"Step Size: $(@bind stepsize PlutoUI.Slider(0.01:0.001:1, default=0.401))"

34.7 ms

xxxxxxxxxx
 
Path(map(0:stepsize:200) do i
    frac_rotation(i) * unitvec * i
end)

35.0 μs

Step Size1:

Step Size2:

Rotation:

xxxxxxxxxx

15.4 ms

xxxxxxxxxx
 
PenPlot(
    [Path(map(0:stepsize1:200) do i
        frac_rotation(i) * unitvec * i
    end)],
    frac_rotation(rotation) * [Path(map(0:stepsize2:200) do i
        frac_rotation(i) * unitvec * i
    end)]
)

159 ms

Representing an Image with Amplitude Modulation

In addition to abstract generative art, people often use real images as an input to plotters. Since we can only draw lines, we have to be creative with how to represent the image.

Fortunately, our brains are adept at recognizing images (especially faces) in patterns of light and dark, so we have a huge creative space of potential techniques to explore. As long as we make sure to put more ink on the page for the darker regions of the image, the image will appear.

The technique I'm going to use is an amplitude-modulating spiral. I'll draw a spiral just as before, but with a sine wave added to it. Then I'll overlay the spiral onto the image and vary the amplitude of the wave by the darkness of the pixel that each point in the spiral falls on to. As a result, the pen will cover more distance around the dark regions of the plot, and the source image will emerge.

xxxxxxxxxx

6.1 μs

xxxxxxxxxx
 
# The size of the wave.
@bind amplitude Slider(0.1:0.1:4, default=1)

14.9 ms

xxxxxxxxxx
 
# Control for the number of rotations.
# This is actually proportional to the _square_ of the number of rotations,
# because we take the square root later to keep the frequency consistent along
# the spiral.
@bind spirals Slider(1:1000, default=500)

9.2 ms

xxxxxxxxxx
 
# The frequency of the wave along the spiral.
@bind frequency Slider(1:100, default=40)

42.6 μs

xxxxxxxxxx
 
Path(map(0:0.01:spirals) do i
        distance = sqrt(i)
        
        rot = frac_rotation(distance)
        pt = rot * unitvec * distance
            
        pt + (rot * unitvec * amplitude * sin(i * frequency))
end)

4.4 ms

Now, we need a base image to draw. My wife Sarah is more photogenic than I am, so I'm using her image.

xxxxxxxxxx

3.8 μs

rawimage

xxxxxxxxxx
 
rawimage = load(download("https://user-images.githubusercontent.com/6933510/113714134-61a95b00-96e8-11eb-9c67-6170996bc6c6.png"))

1.7 s

For simplicity, I'll stick to a single-pen plot, so I'm only interested in the lightness of each pixel. I'll extract that by converting to grayscale.

xxxxxxxxxx

3.2 μs

grayimage

xxxxxxxxxx
 
grayimage = Gray.(rawimage)

117 ms

I'm also applying a blur to the image. This way, every time we probe a pixel in the image we are probing a sample of the pixels in its neighborhood. It mitigates the problem of noise.

xxxxxxxxxx

4.9 μs

xxxxxxxxxx
 
@bind blur Slider(1:10, default=4)

60.1 μs

image

xxxxxxxxxx
 
image = imfilter(grayimage, Kernel.gaussian(blur))

2.4 s

Now all we need is to map from points in the spiral to a pixel value in the image. The probe_val helper function does just that.

xxxxxxxxxx

4.0 μs

probe_val (generic function with 1 method)

xxxxxxxxxx
 
function probe_val(point)
    h, w = size(image)
    image[
        clamp(Int(floor(point[2] + h/2)), 1, h),
        clamp(Int(floor(point[1] + w/2)), 1, w)
    ]
end

36.6 μs

I'll also use a slider for the scale of the image, which essentially acts as a way to ''crop'' the image inside the circle.

xxxxxxxxxx

4.4 μs

xxxxxxxxxx
 
@bind scale Slider(1:20, default=15)

52.7 μs

xxxxxxxxxx
 
Path(map(0:0.01:spirals) do i
        distance = sqrt(i)
        
        rot = frac_rotation(distance)
        pt = rot * unitvec * distance
        
        value = 1 - probe_val(pt * scale)
        
        pt + (rot * unitvec * amplitude * value * sin(i * frequency))
end)

144 ms

The following are some topics that I did not cover in my talk, but have provided for your exploration and experimentation.

Recursion

Drawing a Koch Curve

xxxxxxxxxx

9.5 μs

koch (generic function with 2 methods)

xxxxxxxxxx
 
function koch(i=6)
    scale = Point(1/3, 1/3)
    
    if i == 1
        [Path([Point(0, 0), Point(1, 0)])]
    else
        c = koch(i-1)
        vcat(
            scale * c,
            Point(1/3, 0) + degree_rotation(60) * (scale * c),
            Point(1/2, sqrt(3)/6) + degree_rotation(-60) * (scale * c),
            Point(2/3, 0) + scale * c,
        )
    end
end

56.1 μs

xxxxxxxxxx
 
koch()

126 ms

Drawing a Tree

xxxxxxxxxx

3.3 μs

Angle 1:

Angle 2:

xxxxxxxxxx

78.2 μs

tree (generic function with 2 methods)

xxxxxxxxxx
 
function tree(angle1, angle2, i=8)
    if i == 1
        [Path([Point(0, 0), Point(0, -1)])]
    else
        c = tree(angle1, angle2, i-1)
        vcat(
            [Path([Point(0, 0), Point(0, -1)])],
            Point(0, -1/2) + degree_rotation(angle1) * (Point(0.8, 0.8) * c),
            Point(0, -1) + degree_rotation(angle2) * (Point(0.6, 0.6) * c),
        )
    end
end

50.5 μs

xxxxxxxxxx
 
tree(angle1, angle2)

38.8 ms

Noise Spiral

Seed:

Big Period:

Small Period:

Outer Radius:

xxxxxxxxxx

22.9 ms

noise_spiral (generic function with 1 method)

xxxxxxxxxx
 
function noise_spiral(seed, small_period, big_period, radius)
    noise = random_vector_matrix(MersenneTwister(seed), big_period, small_period)
    
    map(0:0.004:1-eps()) do j
        r = 1 + perlin_noise(noise, Point(j*big_period, 0))
        center = radius * frac_rotation(j) * unitvec * r
​
        Path(map(0:0.01:1) do i
            r = 1 + perlin_noise(noise, Point(i*big_period, j*small_period))
            center + frac_rotation(i) * unitvec * r
      end)
    end
end

97.2 μs

xxxxxxxxxx
 
noise_spiral(seed, big_period, small_period, radius)

277 ms

xxxxxxxxxx
 
PenPlot(
    noise_spiral(6, 2, 5, 5),
    noise_spiral(4, 6, 6, 5)
)

45.5 ms

xxxxxxxxxx
 
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
            Pkg.PackageSpec(url="https://github.com/paulgb/PenPlots.jl"),
            Pkg.PackageSpec(name="PlutoUI", version="0.7"),
            Pkg.PackageSpec(name="Images", version="0.23"),
            Pkg.PackageSpec(name="ImageMagick", version="1"),
            Pkg.PackageSpec(name="ImageFiltering", version="0.6"),
            Pkg.PackageSpec(name="FileIO", version="1"),
            ])
​
    using PenPlots
    using Random
    using Base.Iterators
    using PlutoUI
    using Images, FileIO
    using ImageFiltering
end

55.7 s