Data.Array.Accelerate And Katy Perry, Part II: Convection and Stencils

This time we'll take a really simple convection model (See e.g. Step 5 here) and use it to let our subject flow. We're using grayscale blackness as density and finding derivatives using a Stencil object and the simplest central difference approximation (naturally we can generalise this to higher orders without much difficulty by using a larger stencil and the new coefficients given in the link).

Starting from the same original image as  the previous post we obtain the following:

(The line used to produce the animation was $ convert -loop 0 lol1*1.bmp convect2.gif, with the 1*1 just matching every tenth frame so that the file isn't too big.)

Watch out for the filamentation on either side of Katy's head and the funny thing happening in the top left at the end. This is a pretty simple model, and I haven't been careful with making it make sense, but the dynamics are already starting to look interesting and with any luck we'll get something more realistic going in the next post.

Here's the code (takes a couple of seconds to run on my GeForce GTX 650 Ti, using much larger timesteps causes numerical errors, but it definitely could be made faster):

	import Prelude hiding (map)
	import Data.Array.Accelerate hiding ((++),round,fromIntegral,tail,iterate)
	import Data.Array.Accelerate.IO
	import Data.Array.Accelerate.CUDA
	import Control.Monad
	type RGBA = Word32
	type Image a = Array DIM2 a
	type Local = Stencil3x3 Double -> Exp Double
	totalFrames = 1000
	main = do
	initial_condition <- liftM makeInit $ readImageFromBMP "katy.bmp"
	process initial_condition totalFrames
	print "All done!"
	makeInit (Right x) = x
	process init n = do
	writeFrame (totalFrames-n) result
	print $ "poop: " ++ show n
	if (n==0) then return result else process result (n-1)
	where
	result = (run.again frameskip advance.use) init
	again 1 f = f
	again n f = f. again (n-1) f
	frameskip = 1
	writeFrame n = writeImageToBMP ("lol"++show (10000+n)++".bmp")
	advance :: Acc (Image RGBA) -> Acc (Image RGBA)
	advance = bw (stencil convect Clamp)
	bw :: (Acc (Array DIM2 Double) -> Acc (Array DIM2 Double))
	-> Acc (Image RGBA) -> Acc (Image RGBA)
	bw f a = map rgba32OfLuminance .f $ map luminanceOfRGBA32 a
	convect :: Local
	convect s@((_,t,_),(l,c,r),(_,b,_)) =
	c+c0*(dx s +dy s)
	where c0 = 1e-2
	dx :: Local
	dx ((_,t,_),(l,c,r),(_,b,_)) = 0.5*(l-r)
	dy :: Local
	dy ((_,t,_),(l,c,r),(_,b,_)) = 0.5*(b-t)
	toGreyscale :: Acc (Image RGBA) -> Acc (Image Double)
	toGreyscale = map luminanceOfRGBA32

view raw convectKaty.hs hosted with ❤ by GitHub