> 
with(LinearAlgebra):
with(plots):
with(Statistics):

The principle is always the same:
1/ Let L a straight line which is either defined by its two ending points (xkvd_hline) or taken as the default [0, 0], [1, 0] line.
For xkvd_hline the given line L is firstly rotate to be aligned with the horizontal axis.
2/ Let P1, ..., PN N points on L. Each Pn writes [xn, yn]
3/ A random perturbation rn is added yo the values y1, ..., yN
4/ A stationnary random process RP, with gaussian correlation function is used to build a smooth curve passing through the points
(x1, y1+r1), ..., (xN, yN+rN) (procedure KG where "KG" stands for "Kriging")
5/ The result is drawn or mapped to some predefined shape :
xkcd_hist,
xkcd_polyline,
xkcd_circle
6/ A procedure xkcd_func is also provided to draw functions defined by an explicit relation.
> 
KG := proc(X, Y, psi, sigma)
local NX, DX, K, mu, k, y:
NX := numelems(X);
DX := < seq(Vector[row](NX, X), k=1..NX) >^+:
K := (sigma, psi) > evalf( sigma^2 *~ exp~(((DX  DX^+) /~ psi)^~2) ):
mu := add(Y) / NX;
k := (x, sigma, psi) > evalf( convert(sigma^2 *~ exp~(((x ~ X ) /~ psi)^~2), Vector[row]) ):
y := mu + k(x, sigma, psi) . (K(sigma, psi))^(1) . convert((Y ~ mu), Vector[column]):
return y
end proc:
xkcd_hline := proc(p1::list, p2::list, a::nonnegative, lc::positive, col)
# p1 : first ending point
# p2 : second ending point
# a : amplitude of the random perturbations
# lc : correlation length
# col: color
local roll, NX, LX, X, Z:
roll := rand(1.0 .. 1.0):
NX := 10:
LX := p2[1]p1[1]:
X := [seq(p1[1]..p2[1], LX/(NX1))]:
Z := [p1[2], seq(p1[2]+a*roll(), k=1..NX1)]:
return plot(KG(X, Z, lc*LX, 1), x=min(X)..max(X), color=col, scaling=constrained):
end proc:
xkcd_line := proc(L::list, a::nonnegative, lc::positive, col, {lsty::integer:=1})
# L : list which contains the two ending point
# a : amplitude of the random perturbations
# lc : correlation length
# col: color
local T, roll, NX, DX, DY, LX, A, m, M, X, Z, P:
T := (a, x0, y0, l) >
plottools:transform(
(x,y) > [ x0 + l * (x*cos(a)y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
):
roll := rand(1.0 .. 1.0):
NX := 5:
DX := L[2][1]L[1][1]:
DY := L[2][2]L[1][2]:
LX := sqrt(DX^2+DY^2):
if DX <> 0 then
A := arcsin(DY/LX):
else
A:= Pi/2:
end if:
X := [seq(0..1, 1/(NX1))]:
Z := [ seq(a*roll(), k=1..NX)]:
P := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained, linestyle=lsty):
return T(A, op(L[1]), LX)(P)
end proc:
xkcd_func := proc(f, r::list, NX::posint, a::positive, lc::positive, col)
# f : function to draw
# r : plot range
# NX : number of equidistant "nodes" in the range r (boundaries included)
# a : amplitude of the random perturbations
# lc : correlation length
# col: color
local roll, F, LX, Pf, Xf, Zf:
roll := rand(1.0 .. 1.0):
F := unapply(f, indets(f, name)[1]);
LX := r[2]r[1]:
Pf := [seq(r[1]..r[2], LX/(NX1))]:
Xf := Pf +~ [seq(a*roll(), k=1..numelems(Pf))]:
Zf := F~(Pf) +~ [seq(a*roll(), k=1..numelems(Pf))]:
return plot(KG(Xf, Zf, lc*LX, 1), x=min(Xf)..max(Xf), color=col):
end proc:
xkcd_hist := proc(H, ah, av, ax, ay, lch, lcv, lcx, lcy, colh, colxy)
# H : Histogram
# ah : amplitude of the random perturbations on the horizontal boundaries of the bins
# av : amplitude of the random perturbations on the vertical boundaries of the bins
# ax : amplitude of the random perturbations on the horizontal axis
# ay : amplitude of the random perturbations on the vertical axis
# lch : correlation length on the horizontal boundaries of the bins
# lcv : correlation length on the vertical boundaries of the bins
# lcx : correlation length on the horizontal axis
# lcy : correlation length on the vertical axis
# colh: color of the histogram
# col : color of the axes
local data, horiz, verti, horizontal_lines, vertical_lines, po, rpo, p1, p2:
data := op(1..2, op(1, H)):
verti := sort( [seq(data[n][3..4][], n=1..numelems([data]))] , key=(x>x[1]) );
verti := verti[1],
map(
n > if verti[n][2] > verti[n+1][2] then
verti[n]
else
verti[n+1]
end if,
[seq(2..numelems(verti)2,2)]
)[],
verti[1];
horiz := seq(data[n][[4, 3]], n=1..numelems([data])):
horizontal_lines := NULL:
for po in horiz do
horizontal_lines := horizontal_lines, xkcd_hline(po[1], po[2], ah, lch, colh):
end do:
vertical_lines := NULL:
for po in [verti] do
rpo := po[[2, 1]]:
vertical_lines := vertical_lines, xkcd_hline([0, rpo[2]], rpo, av, lcv, colh):
end do:
p1 := [2*verti[1][1]verti[2][1], 0]:
p2 := [2*verti[1][1]verti[2][1], 0]:
return
display(
horizontal_lines,
T~([vertical_lines]),
xkcd_hline(p1, p2, ax, lcx, colxy),
T(xkcd_hline([0, 0], [1.2*max(op~(2, [verti])), 0], ay, lcy, colxy)),
axes=none,
scaling=unconstrained
);
end proc:
xkcd_polyline := proc(L::list, a::nonnegative, lc::positive, col)
# xkcd_polyline reduces to xkcd_line whebn L has 2 elements
# L : List of points
# a : amplitude of the random perturbations
# lc : correlation length
# col: color
local T, roll, NX, n, DX, DY, LX, A, m, M, X, Z, P:
T := (a, x0, y0, l) >
plottools:transform(
(x,y) > [ x0 + l * (x*cos(a)y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
):
roll := rand(1.0 .. 1.0):
NX := 5:
for n from 1 to numelems(L)1 do
DX := L[n+1][1]L[n][1]:
DY := L[n+1][2]L[n][2]:
LX := sqrt(DX^2+DY^2):
if DX <> 0 then
A := evalf(arcsin(abs(DY)/LX)):
if DX >= 0 and DY <= 0 then A := A end if:
if DX <= 0 and DY > 0 then A := PiA end if:
if DX <= 0 and DY <= 0 then A := Pi+A end if:
else
A:= Pi/2:
if DY < 0 then A := 3*Pi/2 end if:
end if:
if n=1 then
X := [seq(0..1, 1/(NX1))]:
Z := [seq(a*roll(), k=1..NX)]:
else
X := [0 , seq(1/(NX1)..1, 1/(NX1))]:
Z := [Z[NX], seq(a*roll(), k=1..NX1)]:
end if:
P := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained):
Pn := T(A, op(L[n]), LX)(P):
end do;
return seq(Pn, n=1..numelems(L)1)
end proc:
xkcd_circle := proc(a::nonnegative, lc::positive, r::positive, cent::list, col)
# a : amplitude of the random perturbations
# lc : correlation length
# r : redius of the circle
# cent: center of the circle
# col : color
local roll, NX, LX, X, Z, xkg, A:
roll := rand(1.0 .. 1.0):
NX := 10:
X := [seq(0..1, 1/(NX1))]:
Z := [0, seq(a*roll(), k=1..NX1)]:
xkg := KG(X, Z, lc, 1):
A := Pi*roll():
return plot([cent[1]+r*(1+xkg)*cos(2*Pi*x+A), cent[2]+r*(1+xkg)*sin(2*Pi*x+A), x=0..1], color=col)
end proc:
T := plottools:transform((x,y) > [y, x]):

> 
# Axes plot
x_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black):
y_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black):
display(
x_axis,
T(y_axis),
axes=none,
scaling=constrained
)

> 
# A simple function
f := 1+10*(x/51)^2:
F := xkcd_func(f, [0.5, 9.5], 6, 0.3, 0.4, red):
display(
x_axis,
T(y_axis),
F,
axes=none,
scaling=constrained
)

> 
# An histogram
S := Sample(Normal(0,1),100):
H := Histogram(S, maxbins=6):
xkcd_hist(H, 0, 0.02, 0.001, 0.01, 1, 0.1, 0.01, 1, red, black)

> 
# Axes plus grid with two red straight lines
r := rand(0.1 .. 0.1):
x_axis := xkcd_line([[2, 0], [10, 0]], 0.01, 0.2, black):
y_axis := xkcd_line([[0, 2], [0, 10]], 0.01, 0.2, black):
d1 := xkcd_line([[1, 1], [9, 9]] , 0.01, 0.2, red):
d2 := xkcd_line([[1, 9], [9, 1]], 0.01, 0.2, red):
display(
x_axis, y_axis,
seq( xkcd_line([[2+0.3*r(), u+0.3*r()], [10+0.3*r(), u+0.3*r()]], 0.005, 0.5, gray), u in [seq(1..9, 2)]),
seq( xkcd_line([[u+0.3*r(), 2+0.3*r()], [u+0.3*r(), 10+0.3*r()]], 0.005, 0.5, gray), u in [seq(1..9, 2)]),
d1, d2,
axes=none,
scaling=constrained
)

> 
# Axes and a couple of polygonal lines
d1 := xkcd_polyline([[0, 0], [1, 3], [3, 5], [7, 1], [9, 7]], 0.01, 1, red):
d2 := xkcd_polyline([[0, 9], [2, 8], [5, 2], [8, 3], [10, 1]], 0.01, 1, blue):
display(
x_axis, y_axis,
d1, d2,
axes=none,
scaling=constrained
)

> 
# A few polygonal shapes
display(
xkcd_polyline([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]], 0.01, 1, red),
xkcd_polyline([[1/3, 1/3], [2/3, 1/3], [2/3, 4/3], [1, 4/3], [1/3, 1/3]], 0.01, 1, blue),
xkcd_polyline([[1/3, 1/3], [4/3, 1/2], [1/2, 1/2], [1/2,1], [1/3, 1/3]], 0.01, 1, green),
axes=none,
scaling=constrained
)

> 
# A few circles
cols := [red, green, blue, gold, black]: # colors
cents := convert( Statistics:Sample(Uniform(1, 3), [5, 2]), listlist): # centers
radii := Statistics:Sample(Uniform(1/2, 2), 5): # radii
lcs := Statistics:Sample(Uniform(0.2, 0.7), 5): # correlations lengths
display(
seq(
xkcd_circle(0.02, lcs[n], radii[n], cents[n], cols[n]),
n=1..5
),
axes=none,
scaling=constrained
)



> 
# A 3D drawing
x_axis := xkcd_line([[0, 0], [5, 0]], 0.01, 0.2, black):
y_axis := xkcd_line([[0, 0], [4, 2]], 0.01, 0.2, black):
z_axis := xkcd_line([[0, 0], [0, 5]], 0.01, 0.2, black):
f1 := 4*cos(x/6)1:
F1 := xkcd_func(f1, [0.5, 5], 6, 0.001, 0.8, red):
F2 := xkcd_line([[0.5, eval(f1, x=0.5)], [3, 4]], 0.01, 0.2, red):
f3 := 4*cos((x2)/6):
F3 := xkcd_func(f3, [3, 7], 6, 0.001, 0.8, red):
F4 := xkcd_line([[5, eval(f1, x=5)], [7, eval(f3, x=7)]], 0.01, 0.2, red):
dx := xkcd_line([[2, 1], [4, 1]], 0.01, 0.2, gray, lsty=3):
dy := xkcd_line([[2, 0], [4, 1]], 0.01, 0.2, gray, lsty=3):
dz := xkcd_line([[4, 1], [4, 3]], 0.01, 0.2, gray, lsty=3):
po := xkcd_circle(0.02, 0.3, 0.1, [4, 3], blue):
# Numerical value come from "probe info + copy/paste"
nvect := xkcd_polyline([[4, 3], [4.57, 4.26], [4.35, 4.1], [4.57, 4.26], [4.58, 4.02]], 0.01, 1, blue):
tg1vect := xkcd_polyline([[4, 3], [4.78, 2.59], [4.49, 2.87], [4.78, 2.59], [4.46, 2.57]], 0.01, 1, blue):
tg2vect := xkcd_polyline([[4, 3], [4.79, 3.35], [4.70, 3.13], [4.79, 3.35], [4.46, 3.35]], 0.01, 1, blue):
rec1 := xkcd_polyline([[4.118, 3.286], [4.365, 3.396], [4.257, 3.108]], 0.01, 1, blue):
rec2 := xkcd_polyline([[4.257, 3.108], [4.476, 2.985], [4.259, 2.876]], 0.01, 1, blue):
display(
x_axis, y_axis, z_axis,
F1, F2, F3, F4,
dx, dy, dz,
po,
nvect, tg1vect, tg2vect, rec1, rec2,
axes=none,
scaling=constrained
)



> 
# Arrow
d1 := xkcd_polyline([[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, 0.05]], 0.01, 1, red):
T := (a, x0, y0, l) >
plottools:transform(
(x,y) > [ x0 + l * (x*cos(a)y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
):
display(
seq( T(2*Pi*n/10, 0.5, 0, 1/2)(
display(
xkcd_polyline(
[[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, 0.05]],
0.01,
1,
ColorTools:Color([rand()/10^12, rand()/10^12, rand()/10^12])
)
)
),
n=1..10
),
axes=none,
scaling=constrained
)

