Creating Bitmaps from Lisp

Having learned about programming computers later in my education (during university), the internals of computers and binary file formats have been very mysterious to me. Text files (and files with a text format) were usually understandable (or conceivably understandable), but binary files were impenetrable to me. At the time, I had no understanding of bytes, let alone concepts such as little/big-endianness. So, to get a feel for how binary file formats work, I wrote a bitmap writer function.

As reference, I used the Wikipedia page about the BMP file format. The code is released under a BSD license and is available for download here. It can make 24-bit colour bitmaps. It also contains some functions to manipulate colours. Bitmaps are created by calling a write-bmp function.

It was an interesting experience to work with binary streams, and I learned a lot. The code's not that nice, and probably not that efficient either, but it does what I want, which is output pictures.

One of the arguments to write-bmp is used to determine the colour of each pixel, and can be either a list of lists, or a function of two arguments. Here's a call that produces a twelve colour diagonal rainbow:
(write-bmp (path "bitmaps/diagonal-rainbow3.bmp") 80 80
  #'(lambda (row col) (rainbow12 (/ (+ row col) 2))))
The function rainbow12 selects a colour based on the floor of it's input mod 12. The above call produces*:

Here's a fun example using MathP:
(write-bmp (path "bitmaps/picture-perfect.bmp") 80 80
    #M(fn(r c)
       {applR = 2.5
        inRad(x y rad) = (c-x)^2 + (r-y)^2 < rad^2
        if inRad(15 65 6) yellow; inner sun
        if inRad(15 65 8) orange; outer sun
        if inRad(48 55 applR) || inRad(65 45 applR) || inRad(52 40 applR) red
        if (r-50)^2+(c-55)^2<15^2 darkgreen
        if r<20+3*cos(c/10) green; grass
        if 50<c<60 && r<50 brown blue}))
It produces*:

But since computers are good at calculations, let's try something a little more intense (and random, and ... psychedelic):
(write-bmp (path "bitmaps/trigonometric2.bmp")
    360 300
    #'(lambda (row col)
        (rainbow12
            #M{(row^1.1 + 0.6*col + 35*cos(col/15)*sin(row/16)
                + 2 * sqrt((row-250)^2 + (col-200)^2)
                + 20 * cos((col-0.3*row)/70) * (2-sin(row/80))
                + 20 * (1-cos(0.7+col/25)) * cos(0.6-row/50))/15})))
Viola*:

I hope your eyes aren't bleeding too much after that one. ;-)

I've also implemented some anti-aliasing, here's are some anti-aliased pictures:

The last one above is an attempt at creating a 'world' map with things like seas (blue), beaches (yellow), vegetation (green) and mountains (brown and white). The jump from vegetation to beach is sharp, but the anti-aliasing smooths it nicely.

* I've converted all of the pictures displayed here to PNG format using Paint.Net. It feels like cheating, but I was able to get things done. One of these days I'll look into ZPNG and Salza2. The compression used in the PNG format seems very mysterious to me.

I Am With You Always

Last Sunday I presented an exhortation at my church. The title is "I Am With You Always", and was inspired by the last recorded words of Jesus in the gospel of Matthew (28:18-20):

All authority in heaven and on earth has been given to me.

Go therefore and make disciples of all nations,

baptizing them in the name of the Father and of the Son and of the Holy Spirit,

teaching them to observe all that I have commanded you.

And behold, I am with you always, to the end of the age.

You can download my notes here. They're released under a Creative Commons Attribution 3.0 Unported licence. When I presented it, it took about 28 minutes.

I hope you find it an interesting read, and may God be with you.

Cody's Algorithm for ERF

Here's a more substantial example of the use of MathP. Erf(x) is the error function - a function related to the cumulative normal distribution. One way to compute the error function is with Cody's algorithm. I've implemented it in Common Lisp to double precision with MathP.

The implementation here starts by defining two constants (i.e. constants in the following equations), then three helper functions. The helper functions are the three rational functions denoted Rlm in Cody's paper. Here's the code in one big MathP expression:

#M{e = exp(1d0) irpi = 1/sqrt(pi); irp = inverse root pi
rlmErf(x) = (3.209377589138469472562d3;; The 4th degree polynomial from Table II.
      + x * (3.774852376853020208137d2
      + x * (1.138641541510501556495d2
      + x * (3.161123743870565596947d0
      + x *  1.857777061846031526730d\-1))))
      /     (2.844236833439170622273d3
      + x * (1.282616526077372275645d3
      + x * (2.440246379344441733056d2
      + x * (2.360129095234412093499d1 + x))))
rlmErc(x) = (1.23033935479799725272d3;; The 8th degree polynomial from Table III.
      + x * (2.05107837782607146532d3
      + x * (1.71204761263407058314d3
      + x * (8.81952221241769090411d2
      + x * (2.98635138197400131132d2
      + x * (6.61191906371416294775d1
      + x * (8.88314979438837594118d0
      + x * (5.64188496988670089180d\-1
      + x *  2.15311535474403846343d\-8))))))))
      /     (1.23033935480374942043d3
      + x * (3.43936767414372163696d3
      + x * (4.36261909014324715820d3
      + x * (3.29079923573345962678d3
      + x * (1.62138957456669018874d3
      + x * (5.37181101862009857509d2
      + x * (1.17693950891312499305d2
      + x * (1.57449261107098347253d1 + x))))))))
rlmEr2(x) = (-6.58749161529837803157d\-4;; The 5th degree polynomial from Table IV.
      + x * (-1.60837851487422766278d\-2
      + x * (-1.25781726111229246204d\-1
      + x * (-3.60344899949804439429d\-1
      + x * (-3.05326634961232344035d\-1
      + x *  -1.63153871373020978498d\-2)))))
      /     (2.33520497626869185443d\-3
      + x * (6.05183413124413191178d\-2
      + x * (5.27905102951428412248d\-1
      + x * (1.87295284992346047209d0
      + x * (2.56852019228982242072d0 + x)))))
defun erfc(x) {"Placeholder" x}
defun erf (x) {"Error function computed with Cody's algorithm"
    x = float(x pi); Make sure it's a double
    if x < 0.0 (-erf(-x))
    if x < 0.5 x * rlmErf(x^2)
               1d0 - erfc(x)}
defun erfc(x) {"Complementary error function computed with Cody's algorithm"
    x = float(x pi); Make sure it's a double
    if x < 0.5 1d0+erf(-x); Based on identities.
    if x < 4.0 e^-x^2 * rlmErc(x);; These two avoid subtraction from one.
               (oox2=1d0/x^2  e^-x^2/x * {irpi + oox2*rlmEr2(oox2)})}}

To check that it works, we have here a small selection of values for erf and erfc, sampled from the Wikipedia page, compared to those calculated using the above code:

erf(0.1) ≈ 0.1124629 ≈ 0.112463d0
erfc(15)/2 ≈ 3.6065e-100 ≈ 3.606497d-100
erf(0.95) ≈ 0.8208908 ≈ 0.820891d0
erfc(0.7) ≈ 0.3221988 ≈ 0.322199d0

Going a bit further, here are some values based on the arbitrarily precise values computed by Wolfram Alpha:

erf(1/100)
0.0112834155558496160
0.0112834155558496169...
erf(1/3)
0.3626481117660628000
0.3626481117660629334...
erf(1/2)
0.5204998778130465000
0.52049987781304653768...
erf(1)
0.8427007929497148000
0.8427007929497148693...
erf(4)
0.9999999845827421000
0.9999999845827420997...

So it appears (from this small sample at least) that this implementation is accurate to about 1e-16 (see the value for erf(1/3)), but usually 1e-17, which appears to be the resolution of doubles for Corman Common Lisp.

MathP is in git

Good news! The source code (and other tidbits) for my Common Lisp library MathP is available from my public git repository at bitbucket here. Please, go and fork as desired.

I've given MathP a logo at the same time:

Other logos I came up with and considered were:

Like a specific one? Let me know?

And since I want to show you more MathP, here's an implementation of exponential smoothing on the list of data points 0, 1, 2, 3 and 1:

#M{theta=0.2 d=nil
   s(x) = d := if !d x theta*x+(1-theta)*d
   mapcar(#'s '(0 1 2 3 1))}

It should produce the list:
(0 0.2 0.56 1.048 1.0384)