Koch’s Snowflake

December 23, 2011

We produce output by writing to a PostScript file. Here’s the main program that writes the file header, calls the koch function to produce the snowflake, then writes the file footer:

(define (draw-koch n)
  (display "%!") (newline)
  (display "newpath") (newline)
  (display "426 316 moveto") (newline)
  (display "90 rotate") (newline)
  (display "2 setlinewidth") (newline)
  (koch n 243)
  (display "stroke") (newline)
  (display "showpage") (newline))

The initial coordinates 426 and 316 are arbitrary, but seem to work well. The koch function computes and writes the snowflake:

(define (koch n len)
  (let koch ((pat (string->list "F++F++F")) (n n) (len len))
    (if (positive? n)
        (let ((new (replace pat #\F (string->list "F-F++F-F"))))
          (koch new (- n 1) (/ len 3)))
        (let draw ((pat pat))
          (when (pair? pat)
            (display "currentpoint translate") (newline)
            (case (car pat)
              ((#\F) (display 0) (display " ") (display len)
                     (display " rlineto") (newline))
              ((#\-) (display "60 rotate") (newline))
              ((#\+) (display "-60 rotate") (newline)))
            (draw (cdr pat)))))))

The outer loop koch computes the rewrite steps to the nth fractal; the inner loop draw writes the PostScript output. The main trick in the PostScript output is the command currentpoint translate, which moves the coordinate axes so each line can start at 0. The replace function modifies a list:

(define (replace xs old new)
  (let loop ((xs xs) (zs (list)))
    (cond ((null? xs) (reverse zs))
          ((equal? (car xs) old)
            (loop (cdr xs) (append (reverse new) zs)))
          (else (loop (cdr xs) (cons (car xs) zs))))))

To produce the 4-fractal PostScript output, say:

> (with-output-to-file "koch.ps"
    (lambda () (draw-koch 4)))

You can run the code and see the PostScript output at http://programmingpraxis.codepad.org/vDdsP7o6; the image produced by the program is shown above right.

Merry Christmas!


Pages: 1 2

6 Responses to “Koch’s Snowflake”

  1. Jussi Piitulainen said

    I noticed there are only six possible directions, and with each axis scaled properly, the arithmetic is very simple. The following computes the coordinates for a given list of moves, F, +, or -, starting at origin and initial angle of 60 degrees. (The x and y coordinates are in separate lists and in the backwards order.)

    (define (kochdinates moves)
      (let loop ((x 0) (y 0) (xs '(0)) (ys '(0))
                 (theta 60) (moves moves))
        (if (null? moves)
            (values xs ys)
            (case (car moves)
              ((F) (call-with-values
                       (lambda ()
                         (case theta
                           (( 0) (values (+ x 2) y))
                           (( 60) (values (+ x 1) (+ y 1)))
                           ((120) (values (- x 1) (+ y 1)))
                           ((180) (values (- x 2) y))
                           ((240) (values (- x 1) (- y 1)))
                           ((300) (values (+ x 1) (- y 1)))))
                     (lambda (x1 y1)
                       (loop x1 y1 (cons x1 xs) (cons y1 ys)
                             theta (cdr moves)))))
              ((+) (loop x y xs ys (modulo (- theta 60) 360) (cdr moves)))
              ((-) (loop x y xs ys (modulo (+ theta 60) 360) (cdr moves)))))))

  2. Jussi Piitulainen said

    Here’s the code that produces the snowflake moves for the coordinate computer above. Nest at will.

    (define (init) '(F + + F + + F))

    (define (rewrite moves)
      (apply append
             (map (lambda (move)
                    (case move
                      ((F) '(F - F + + F - F))
                      ((+) '(+))
                      ((-) '(-))))

  3. Jussi Piitulainen said

    And here’s a compiler that produces an R command to plot a snowflake.

    (define (R xs ys)
      (display "plot(c(") (display (car xs))
      (for-each (lambda (x) (display ",") (display x)) (cdr xs))
      (display "), c(") (display (car ys))
      (for-each (lambda (y) (display ",") (display y)) (cdr ys))
      (display "), type='l', xlab='', ylab='')") (newline))

    For example, this:

    (call-with-values (lambda () (kochdinates (rewrite (init)))) R)

    prints this:

    plot(c(0,2,3,4,6,5,6,4,3,2,0,1,0), c(0,0,-1,0,0,1,2,2,3,2,2,1,0), type='l', xlab='', ylab='')

    which can be pasted in an R session to open the plot window. For scaling on the screen, grab the corner of the plot window and let the computer do the arithmetic. (I ran Guile and R in one Emacs shell session, so cutting and pasting was easy, but I had to close one interpreter to enter the other. And the R command line editor seemed to choke when there were too many coordinates :)

  4. […] lovely programming exercise. I took the task one step further and created an interactive explorer of variants of the Koch […]

  5. Axio said

    Use your favourite LOGO interpreter…

    (defparameter *koch* '(f + + f + + f))
    (defun koch-iterate (koch iterations)
      (loop while (> iterations 0)
         do (setq iterations (1- iterations))
           (setq koch (apply #'append
                                (case v
                                  ((f) '(f - f + + f - f))
                                  (t `(,v))))
    (defun 2logo ()
      (with-output-to-string  (s)
           (case v
             ((f) (format s "fd 5 "))
             ((+) (format s "rt 60 "))
             ((-) (format s "lt 60 "))))
         (koch-iterate *koch* 5))))
  6. Globules said

    Here’s a Haskell version, using version 0.4 of the diagrams library. The heart of the program is the step function. Given an input shape, _, it produces _/\_. This is iterated some number of times (the edge function), starting with a single, straight line segment. Lastly, three copies of the final list of segments are glued together in a triangular shape.

    import Diagrams.Prelude
    import Diagrams.Backend.Cairo.CmdLine
    koch n = let e = edge !! n
                 t = mconcat [e, rotateBy (-1/3) e, rotateBy (-2/3) e]
             in pathFromTrail $ scale (1/3^n) t
      where edge = iterate step $ fromSegments [straight (1,0)]
            step p = mconcat [p, rotateBy (1/6) p, rotateBy (-1/6) p, p]
    main = defaultMain . attrs . stroke $ koch 5
      where attrs = lw 0.002 . pad 1.1

    Here are a few examples of running the program to generate PNG, PDF and SVG output.

    $ ./koch -h 800 -w 800 -o /tmp/koch.png
    $ ./koch -o /tmp/koch.pdf
    $ ./koch -o /tmp/koch.svg

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: