A quick naïve Haskell implementation

turing_morphogenesis.hs

import Control.Monad (forM_)
import Data.Array
import Data.Ix
import System.Random (randomIO)

-- A 'Size' is just a pair of Ints
type Pair = (Int, Int)
-- An 'Image' is a packed, pair-indexed array
type Image a = Array Pair a

-- Keep a value between two min/max values
clamp :: Ord a => a -> a -> a -> a
clamp min max x
  | x < min   = min
  | x > max   = max
  | otherwise = x

-- Create a grid of random floats
mkInitial :: Pair -> IO (Image Float)
mkInitial (w,h) = do
  rs <- sequence [ sequence ((x, y), randomIO)
                 | x <- [1..w]
                 , y <- [1..h]
                 ]
  return (array ((1,1), (w,h)) rs)

-- Create a new image from an old image by choosing a new value
-- for each grid square based on the averages of two neighborhoods
-- around the old value
step :: Image Float -> Image Float
step img = img // [ (i,
                     clamp 0.0 1.0 $
                       if getAround i 5 img >
                          getAround i 10 img
                       then v + 0.05
                       else v - 0.05)
                  | i <- indices img
                  , let v = img ! i
                  ]

-- Get the average of the Floats within a given radius
-- of a point, calculated via Manhattan distance
getAround :: Pair -> Int -> Image Float -> Float
getAround (x, y) n img =
  let bs = bounds img
  in average [ img ! idx
             | j <- [-n..n]
             , k <- [-n..n]
             , let idx = (x + j, y + k)
             , inRange bs idx
             ]

-- Take the average of a list of Floats
average :: [Float] -> Float
average xs = sum xs / fromIntegral (length xs)

-- Convert a grid of Floats into a grid of Ints
discretize :: Int -> Image Float -> Image Int
discretize n = fmap (floor . (* (fromIntegral n)))

-- Print an 'Image' as a PBM file
pPBM :: Int -> Image Int -> IO ()
pPBM max arr = do
  putStrLn "P2"
  let ((wb, hb), (wm, hm)) = bounds arr
      (w, h) = (wm - wb + 1, hm - hb + 1)
  putStrLn (show w ++ " " ++ show h)
  putStrLn (show max)
  forM_ [wb..wm] $ \x -> do
    forM_ [hb..hm] $ \y -> do
      putStr (show (arr ! (x,y)))
      putStr " "
    putStrLn ""

-- Repeatedly apply a function
iter :: Int -> (a -> a) -> a -> a
iter 0 _ x = x
iter n f x = f (iter (n-1) f x)

main :: IO ()
main = do
  i <- mkInitial (256,256)
  let j = iter 25 step i
  pPBM 128 (discretize 128 j)
  return ()