repa-examples 1.1.1.0 → 2.0.0.1
raw patch · 16 files changed
+1349/−519 lines, 16 filesdep +template-haskelldep +vectordep −dph-basedep −dph-prim-pardep ~repadep ~repa-algorithmsdep ~repa-ionew-component:exe:repa-blurnew-component:exe:repa-cannynew-component:exe:repa-fft3d-highpassnew-component:exe:repa-sobel
Dependencies added: template-haskell, vector
Dependencies removed: dph-base, dph-prim-par
Dependency ranges changed: repa, repa-algorithms, repa-io
Files
- FFT/src/FFT2d/Main.hs +0/−52
- FFT/src/HighPass/Main.hs +0/−73
- Laplace/src/Main.hs +0/−152
- Laplace/src/Solver.hs +0/−82
- MMult/src/Main.hs +0/−124
- examples/Blur/src-repa/Main.hs +98/−0
- examples/Canny/src-repa/Main.hs +384/−0
- examples/FFT/FFT2d/src-repa/Main.hs +55/−0
- examples/FFT/HighPass2d/src-repa/Main.hs +91/−0
- examples/FFT/HighPass3d/src-repa/Main.hs +105/−0
- examples/Laplace/src-repa/Main.hs +181/−0
- examples/Laplace/src-repa/SolverGet.hs +74/−0
- examples/Laplace/src-repa/SolverStencil.hs +51/−0
- examples/MMult/src-repa/Main.hs +108/−0
- examples/Sobel/src-repa/Main.hs +85/−0
- repa-examples.cabal +117/−36
− FFT/src/FFT2d/Main.hs
@@ -1,52 +0,0 @@---- | Perform the 2D FFT on a BMP image.-import Data.Array.Repa.Algorithms.FFT-import Data.Array.Repa.Algorithms.DFT.Center-import Data.Array.Repa.Algorithms.Complex-import Data.Array.Repa.IO.BMP-import Data.Array.Repa as A-import System.Environment-import Control.Monad--main :: IO ()-main - = do args <- getArgs- case args of- [fileIn, clipMag, fileMag, filePhase]- -> mainWithArgs fileIn (read clipMag) fileMag filePhase-- _ -> putStr $ unlines- [ "Usage: repa-fft2d <fileIn.bmp> <clip mag :: Int> <fileOutMag.bmp> <fileOutPhase.bmp>"- , ""- , " The output magnitude has a high dynamic range. We need to clip it otherwise"- , " most of the pixels in the output BMP will be black. Start with a value equal"- , " to about the width of the image (eg 512)"- , "" ]- - -mainWithArgs fileIn clipMag fileMag filePhase- = do - -- Load in the matrix.- arrReal <- liftM (either (\e -> error $ show e) id)- $ readMatrixFromGreyscaleBMP fileIn- let arrComplex = force $ A.map (\r -> r :*: 0) arrReal- - -- Apply the centering transform so that the output has the zero- -- frequency in the middle of the image.- let arrCentered = centerMatrix arrComplex- - -- Do the 2d transform.- let arrFreq = fft2d arrCentered- - -- Write out the magnitude of the transformed array, - -- clipping it at the given value.- let clip m = if m >= clipMag then clipMag else m- let arrMag = A.map (clip . mag) arrFreq- writeMatrixToGreyscaleBMP fileMag arrMag-- -- Write out the phase of the transofmed array, - -- scaling it to make full use of the 8 bit greyscale.- let scaledArg x = (arg x + pi) * (255 / (2 * pi))- let arrPhase = A.map scaledArg arrFreq- writeMatrixToGreyscaleBMP filePhase arrPhase-
− FFT/src/HighPass/Main.hs
@@ -1,73 +0,0 @@-{-# LANGUAGE PatternGuards #-}---- | Perform high pass filtering on a BMP image.-import Data.Array.Repa.Algorithms.FFT-import Data.Array.Repa.Algorithms.DFT.Center-import Data.Array.Repa.Algorithms.Complex-import Data.Array.Repa.IO.BMP-import Data.Array.Repa as A-import System.Environment-import Control.Monad--main :: IO ()-main - = do args <- getArgs- case args of- [cutoff, fileIn, fileOut]- -> mainWithArgs (read cutoff) fileIn fileOut-- _ -> putStr $ unlines- [ "Usage: repa-fft-highpass <cutoff::Int> <fileIn.bmp> <fileOut.bmp>"- , "" ]- - -mainWithArgs cutoff fileIn fileOut- = do - -- Load in the matrix.- (arrRed, arrGreen, arrBlue)- <- liftM (either (\e -> error $ show e) id)- $ readComponentsFromBMP fileIn- - -- Do the transform on each component individually- let arrRed' = transform cutoff arrRed- let arrGreen' = transform cutoff arrGreen- let arrBlue' = transform cutoff arrBlue- - -- Write it back to file.- writeComponentsToBMP fileOut- arrRed' arrGreen' arrBlue' - --transform cutoff arrReal- = let arrComplex = force $ A.map (\r -> (fromIntegral r) :*: 0) arrReal- - -- Do the 2d transform.- arrCentered = centerMatrix arrComplex- arrFreq = fft2d arrCentered-- -- Zap out the low frequency components.- _ :. height :. width = extent arrFreq- centerX = width `div` 2- centerY = height `div` 2- - {-# INLINE highpass #-}- highpass get ix@(_ :. y :. x)- | x > centerX + cutoff- || x < centerX - cutoff- || y > centerY + cutoff- || y < centerY - cutoff- = get ix- - | otherwise- = 0- - arrFilt = traverse arrFreq id highpass-- -- Do the inverse transform to get back to image space.- arrInv = ifft2d arrFilt- - -- Write out the magnitude of the transformed array, - arrMag = A.map (truncate . mag) arrInv-- in arrMag-
− Laplace/src/Main.hs
@@ -1,152 +0,0 @@-{-# LANGUAGE BangPatterns #-}---- | Solver for the Laplace equation--- You supply a BMP files specifying the boundary conditions.--- The output is written back to another BMP file.----import Solver-import Data.Array.Repa as A-import Data.Array.Repa.IO.BMP -import Data.Array.Repa.IO.ColorRamp-import System.Environment-import Data.Word-import Control.Monad--main :: IO ()-main - = do args <- getArgs- case args of- [steps, fileInput, fileOutput] - -> laplace (read steps) fileInput fileOutput-- _ -> do- putStr usage- return ()----- | Command line usage information.-usage :: String-usage = unlines- [ "Usage: laplace <iterations> <input.bmp> <output.bmp>"- , ""- , " iterations :: Int Number of iterations to use in the solver."- , " input.bmp :: FileName Uncompressed RGB24 or RGBA32 BMP file for initial and boundary values."- , " output.bmp :: FileName BMP file to write output to."- , "" - , " Format of input file:"- , " Boundary values are indicated in greyscale,"- , " ie from the list [(x, x, x) | x <- [0 .. 255]]"- , " Non-boundary values are indicated in blue,"- , " ie (0, 0, 255)"- , " Any other pixel value is an error." - , ""- ]- ---- | Solve it.-laplace :: Int -- ^ Number of iterations to use.- -> FilePath -- ^ Input file.- -> FilePath -- ^ Output file- -> IO ()--laplace steps fileInput fileOutput- = do- -- Load up the file containing boundary conditions.- arrImage <- liftM (either (error . show) id)- $ readImageFromBMP fileInput-- let arrBoundValue = force $ slurpDoublesFromImage slurpBoundValue arrImage- let arrBoundMask = force $ slurpDoublesFromImage slurpBoundMask arrImage- - -- Use the boundary condition values as the initial matrix.- let arrInitial = arrBoundValue-- -- Run the solver.- let arrFinal = solveLaplace- steps- arrBoundMask- arrBoundValue- arrInitial-- arrFinal `deepSeqArray` return ()-- -- Make the result image- let arrImageOut = makeImageFromDoubles (rampColorHotToCold 0.0 1.0) arrFinal-- -- Write out the image to file. - writeImageToBMP- fileOutput- arrImageOut----slurpDoublesFromImage- :: (Word8 -> Word8 -> Word8 -> Double)- -> Array DIM3 Word8- -> Array DIM2 Double- -{-# INLINE slurpDoublesFromImage #-}-slurpDoublesFromImage mkDouble arrBound- = traverse arrBound- (\(Z :. height :. width :. _) - -> Z :. height :. width)-- (\get (Z :. y :. x)- -> mkDouble- (get (Z :. y :. x :. 0))- (get (Z :. y :. x :. 1))- (get (Z :. y :. x :. 2)))---makeImageFromDoubles- :: (Double -> (Double, Double, Double))- -> Array DIM2 Double- -> Array DIM3 Word8- -{-# INLINE makeImageFromDoubles #-}-makeImageFromDoubles fnColor arrDoubles- = traverse arrDoubles- (\(Z :. height :. width)- -> Z :. height :. width :. 4)- - (\get (Z :. y :. x :. c)- -> let (r, g, b) = fnColor (get (Z :. y :. x))- in case c of- 0 -> truncate (r * 255)- 1 -> truncate (g * 255)- 2 -> truncate (b * 255)- 3 -> 0)----- | Extract the boundary value from a RGB triple.-slurpBoundValue :: Word8 -> Word8 -> Word8 -> Double-{-# INLINE slurpBoundValue #-}-slurpBoundValue r g b- -- A non-boundary value.- | r == 0 && g == 0 && b == 255 - = 0-- -- A boundary value.- | (r == g) && (r == b) - = fromIntegral r / 255- - | otherwise- = error $ "Unhandled pixel value in input " ++ show (r, g, b)----- | Extract boundary mask from a RGB triple.-slurpBoundMask :: Word8 -> Word8 -> Word8 -> Double-{-# INLINE slurpBoundMask #-}-slurpBoundMask r g b- -- A non-boundary value.- | r == 0 && g == 0 && b == 255 - = 1-- -- A boundary value.- | (r == g) && (r == b) - = 0- - | otherwise- = error $ "Unhandled pixel value in input " ++ show (r, g, b)- -
− Laplace/src/Solver.hs
@@ -1,82 +0,0 @@-{-# LANGUAGE BangPatterns #-}-module Solver - (solveLaplace)-where -import Data.Array.Repa as A-import qualified Data.Array.Repa.Shape as S-import Data.Array.Parallel.Unlifted (Elt)-import qualified Data.Array.Parallel.Unlifted as U----- | Solver for the Laplace equation.-solveLaplace- :: Int -- ^ Number of iterations to use.- -> Array DIM2 Double -- ^ Boundary value mask.- -> Array DIM2 Double -- ^ Boundary values.- -> Array DIM2 Double -- ^ Initial state.- -> Array DIM2 Double--{-# INLINE solveLaplace #-}-solveLaplace steps arrBoundMask@Manifest{} arrBoundValue@Manifest{} arrInit- = go steps arrInit- where go i arr@Manifest{}- | i == 0 = arr- | otherwise = go (i - 1) - $ force- $ applyBoundary arrBoundMask arrBoundValue- $ relaxLaplace arr----- | Perform matrix relaxation for the Laplace equation,--- using a stencil function.------ Computation fn is--- u'(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1)) / 4----relaxLaplace- :: Array DIM2 Double- -> Array DIM2 Double--{-# INLINE relaxLaplace #-}-relaxLaplace !arr- = traverse arr id elemFn- where- _ :. height :. width - = extent arr-- {-# INLINE elemFn #-}- elemFn get d@(sh :. i :. j)- = if isBorder i j- then get d- else (get (sh :. (i-1) :. j)- + get (sh :. i :. (j-1))- + get (sh :. (i+1) :. j)- + get (sh :. i :. (j+1))) / 4-- -- Check if this element is on the border of the matrix.- -- If so we can't apply the stencil because we don't have all the neighbours.- {-# INLINE isBorder #-}- isBorder i j- = (i == 0) || (i >= width - 1) - || (j == 0) || (j >= height - 1) ----- | Apply the boundary conditions to this matrix.--- The mask matrix has 0 in places where boundary conditions hold--- and 1 otherwise.------ The value matrix has the boundary condition value in places where it holds,--- and 0 otherwise.--- -applyBoundary- :: Array DIM2 Double -- ^ Boundary condition mask.- -> Array DIM2 Double -- ^ Boundary condition values.- -> Array DIM2 Double -- ^ Initial matrix.- -> Array DIM2 Double -- ^ Matrix with boundary conditions applied.--{-# INLINE applyBoundary #-}-applyBoundary arrBoundMask arrBoundValue arr- = A.zipWith (+) arrBoundValue- $ A.zipWith (*) arrBoundMask arr--
− MMult/src/Main.hs
@@ -1,124 +0,0 @@-{-# LANGUAGE PatternGuards #-}--import Data.Array.Repa as A-import Data.Array.Repa.IO.Matrix-import Data.Array.Repa.Algorithms.Matrix-import Data.Maybe-import System.Environment-import Control.Monad-import System.Random-import qualified Data.Array.Parallel.Unlifted as U---- Arg Parsing -------------------------------------------------------------------------------------data Arg- = ArgSolver String- | ArgMatrixRandom Int Int- | ArgMatrixFile FilePath- | ArgOutFile FilePath- deriving Show--isArgMatrix arg- = case arg of- ArgMatrixRandom{} -> True- ArgMatrixFile{} -> True- _ -> False--parseArgs [] = []-parseArgs (flag:xx)- | "-file" <- flag- , f:rest <- xx- = ArgMatrixFile f : parseArgs rest-- | "-out" <- flag- , f:rest <- xx- = ArgOutFile f : parseArgs rest- - | "-random" <- flag- , x:y:rest <- xx- = ArgMatrixRandom (read x) (read y) : parseArgs rest- - | otherwise - = error $ "bad arg " ++ flag ++ "\n"--printHelp- = putStr - $ unlines- [ "Usage: mmult [args..]"- , ""- , " -random <height> <width> Use a random matrix of this size."- , " -file <filename> Read a matrix from this file."- , " -out <filename> Write resulting matrix to this file."- , ""- , " Format of matrix file:"- , " MATRIX"- , " <width> <height>"- , " <whitespace separated values..>"- , "" ]----- | Get a matrix from a file, or generate a random one.-getMatrix :: Arg -> IO (Array DIM2 Double)-getMatrix arg- = case arg of- ArgMatrixFile fileName - -> readMatrixFromTextFile fileName-- ArgMatrixRandom height width - -> genRandomMatrix (Z :. height :. width) ----- Random -------------------------------------------------------------------------------------------- | Generate a random(ish) matrix.-genRandomMatrix - :: DIM2 - -> IO (Array DIM2 Double)--genRandomMatrix sh- = do uarr <- genRandomUArray (A.size sh)- return $ fromUArray sh uarr---- | Generate a random(ish) UArray of doubles.--- The std random function is too slow to generate really big vectors--- with. Instead, we generate a short random vector and repeat that.-genRandomUArray :: Int -> IO (U.Array Double)-genRandomUArray n - = do let k = 1000- rg <- newStdGen- let randvec = U.randomRs k (-100, 100) rg- let vec = U.map (\i -> randvec U.!: (i `mod` k)) (U.enumFromTo 0 (n-1))- return vec-- --- Main --------------------------------------------------------------------------------------------main :: IO ()-main - = do args <- liftM parseArgs $ getArgs- main' args--main' args- | [argMat1, argMat2] <- filter isArgMatrix args- , mArgOut <- listToMaybe [s | ArgOutFile s <- args]- = do - -- Get matrices from files, - -- or generate random ones we were asked to.- mat1 <- getMatrix argMat1- mat2 <- getMatrix argMat2-- mat1- `deepSeqArray` mat2- `deepSeqArray` return ()- - -- Run the solver.- let matResult = multiplyMM mat1 mat2-- matResult `deepSeqArray` return ()-- -- Write the output to file if requested.- case mArgOut of - Nothing -> return ()- Just fileOut -> writeMatrixToTextFile fileOut matResult- - | otherwise- = printHelp--
+ examples/Blur/src-repa/Main.hs view
@@ -0,0 +1,98 @@+{-# LANGUAGE PackageImports, BangPatterns, TemplateHaskell, QuasiQuotes #-}+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}++import Data.List+import Control.Monad+import System.Environment+import Data.Word+import Data.Array.Repa.IO.BMP+import Data.Array.Repa.IO.Timing+import Data.Array.Repa.Algorithms.Iterate+import Data.Array.Repa as A+import Data.Array.Repa.Stencil as A+import Prelude as P++main + = do args <- getArgs+ case args of+ [iterations, fileIn, fileOut] -> run (read iterations) fileIn fileOut+ _ -> usage++usage = putStr $ unlines+ [ "repa-blur <iterations::Int> <fileIn.bmp> <fileOut.bmp>" ]+ +run iterations fileIn fileOut+ = do comps <- liftM (either (error . show) id) + $ readComponentsListFromBMP fileIn++ comps `deepSeqArrays` return ()+ + (comps', tElapsed)+ <- time $ let comps' = P.map (process iterations) comps+ in comps' `deepSeqArrays` return comps'+ + putStr $ prettyTime tElapsed+ + writeComponentsListToBMP fileOut comps'++{-# NOINLINE process #-}+process :: Int -> Array DIM2 Word8 -> Array DIM2 Word8+process iterations = demote . blur iterations . promote++ +{-# NOINLINE promote #-}+promote :: Array DIM2 Word8 -> Array DIM2 Double+promote arr@(Array _ [Region RangeAll (GenManifest _)])+ = arr `deepSeqArray` force+ $ A.map ffs arr++ where {-# INLINE ffs #-}+ ffs :: Word8 -> Double+ ffs x = fromIntegral (fromIntegral x :: Int)+++{-# NOINLINE demote #-}+demote :: Array DIM2 Double -> Array DIM2 Word8+demote arr@(Array _ [Region RangeAll (GenManifest _)])+ = arr `deepSeqArray` force+ $ A.map ffs arr++ where {-# INLINE ffs #-}+ ffs :: Double -> Word8+ ffs x = fromIntegral (truncate x :: Int)+++{-# NOINLINE blur #-}+blur :: Int -> Array DIM2 Double -> Array DIM2 Double+blur !iterations arr@(Array _ [Region RangeAll (GenManifest _)])+ = arr `deepSeqArray` force2+ $ iterateBlockwise' iterations arr+ $ A.map (/ 159)+ . mapStencil2 BoundClamp+ [stencil2| 2 4 5 4 2+ 4 9 12 9 4+ 5 12 15 12 5+ 4 9 12 9 4+ 2 4 5 4 2 |]+ ++{- version using convolveOut+-- | Run several iterations of blurring.+blur :: Int -> Array DIM2 Double -> Array DIM2 Double+blur 0 arr = arr+blur n arr = blurs (n - 1) (force $ blur arr)++-- | Run a single iteration of blurring.+blur :: Array DIM2 Double -> Array DIM2 Double+blur input@Manifest{}+ = convolveOut outClamp kernel input+ where kernel = force + $ A.fromList (Z :. 5 :. 5) + $ Data.List.map (\x -> x / 159) + [2.0, 4.0, 5.0, 4.0, 2.0,+ 4.0, 9.0, 12.0, 9.0, 4.0,+ 5.0, 12.0, 15.0, 12.0, 5.0,+ 4.0, 9.0, 12.0, 9.0, 4.0,+ 2.0, 4.0, 5.0, 4.0, 2.0]+-} +
+ examples/Canny/src-repa/Main.hs view
@@ -0,0 +1,384 @@+{-# LANGUAGE PackageImports, BangPatterns, QuasiQuotes, PatternGuards #-}+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}++-- | Canny edge detector.+--+-- NOTE: for best performance this needs to be compiled with the following GHC options:+-- -fllvm -optlo-O3 -Odph -fno-liberate-case+-- -funfolding-use-threshold100 -funfolding-keeness-factor100+--+import Data.List+import Data.Word+import Data.Int+import Control.Monad+import System.Environment+import Data.Array.Repa as R+import Data.Array.Repa.Stencil+import Data.Array.Repa.Specialised.Dim2+import Data.Array.Repa.IO.BMP+import Data.Array.Repa.IO.Timing+import Prelude hiding (compare)++import System.IO.Unsafe+import qualified Data.Vector.Unboxed.Mutable as VM+import qualified Data.Vector.Unboxed as V+import Prelude as P++type Image a = Array DIM2 a++-- Constants ------------------------------------------------------------------+-- TODO: It would be better to use Word8 to represent the edge orientations,+-- but doing so currently triggers a bug in the LLVM mangler.+orientUndef = 0 :: Float+orientPosDiag = 64 :: Float+orientVert = 128 :: Float+orientNegDiag = 192 :: Float+orientHoriz = 255 :: Float++data Edge = None | Weak | Strong+edge None = 0 :: Word8+edge Weak = 128 :: Word8+edge Strong = 255 :: Word8++-- Main routine ---------------------------------------------------------------+main + = do args <- getArgs+ case args of+ [fileIn, fileOut] + -> run 0 50 100 fileIn fileOut++ [loops, threshLow, threshHigh, fileIn, fileOut]+ -> run (read loops) (read threshLow) (read threshHigh) fileIn fileOut++ _ -> putStrLn "repa-canny [<loops::Int> <threshLow::Int> <threshHigh::Int>] <fileIn.bmp> <fileOut.bmp>"+++run loops threshLow threshHigh fileIn fileOut+ = do arrInput <- liftM (force . either (error . show) id) + $ readImageFromBMP fileIn++ arrInput `deepSeqArray` return ()+ (arrResult, tTotal)+ <- time+ $ do {-let stage str arr + = timeStage loops str $ return arr-}+ + arrGrey <- timeStage loops "toGreyScale" $ return $ toGreyScale arrInput+ arrBluredX <- timeStage loops "blurX" $ return $ blurSepX arrGrey+ arrBlured <- timeStage loops "blurY" $ return $ blurSepY arrBluredX+ arrDX <- timeStage loops "diffX" $ return $ gradientX arrBlured+ arrDY <- timeStage loops "diffY" $ return $ gradientY arrBlured+ + arrMagOrient <- timeStage loops "magOrient" + $ return $ gradientMagOrient threshLow arrDX arrDY++ arrSuppress <- timeStage loops "suppress" + $ return $ suppress threshLow threshHigh arrMagOrient++ arrStrong <- timeStage loops "select" $ return $ selectStrong arrSuppress + arrEdges <- timeStage loops "wildfire" $ return $ wildfire arrSuppress arrStrong ++ return arrEdges++ when (loops >= 1)+ $ putStrLn $ "\nTOTAL\n"+ + putStr $ prettyTime tTotal+ + writeComponentsToBMP fileOut arrResult arrResult arrResult+++-- | Wrapper to time each stage of the algorithm.+timeStage+ :: (Shape sh, Elt a)+ => Int+ -> String + -> (IO (Array sh a))+ -> (IO (Array sh a))++{-# NOINLINE timeStage #-}+timeStage loops name fn+ = do let burn !n+ = do arr <- fn+ arr `deepSeqArray` return ()+ if n <= 1 then return arr+ else burn (n - 1)+ + (arrResult, t)+ <- time $ do arrResult' <- burn loops+ arrResult' `deepSeqArray` return arrResult'++ when (loops >= 1) + $ putStr $ name P.++ "\n"+ P.++ unlines [ " " P.++ l | l <- lines $ prettyTime t ]++ return arrResult+++-------------------------------------------------------------------------------+-- | The default conversions supplied via the prelude go via a GMP function+-- call instead of just using the appropriate primop.+{-# INLINE floatOfWord8 #-}+floatOfWord8 :: Word8 -> Float+floatOfWord8 w8+ = fromIntegral (fromIntegral w8 :: Int)+++{-# INLINE word8OfFloat #-}+word8OfFloat :: Float -> Word8+word8OfFloat f+ = fromIntegral (truncate f :: Int)+++-------------------------------------------------------------------------------+-- | RGB to greyscale conversion.+{-# NOINLINE toGreyScale #-}+toGreyScale :: Array DIM3 Word8 -> Array DIM2 Word8+toGreyScale + arr@(Array _ [Region RangeAll (GenManifest _)])+ = arr `deepSeqArray` force2+ $ unsafeTraverse arr+ (\(sh :. _) -> sh)+ (\get ix -> rgbToLuminance + (get (ix :. 0))+ (get (ix :. 1))+ (get (ix :. 2)))++ where {-# INLINE rgbToLuminance #-}+ rgbToLuminance :: Word8 -> Word8 -> Word8 -> Word8+ rgbToLuminance r g b + = word8OfFloat + ( floatOfWord8 r * 0.3+ + floatOfWord8 g * 0.59+ + floatOfWord8 b * 0.11)+++-- | Separable Gaussian blur in the X direction.+{-# NOINLINE blurSepX #-}+blurSepX :: Array DIM2 Word8 -> Array DIM2 Float+blurSepX arr@(Array _ [Region RangeAll (GenManifest _)]) + = arr `deepSeqArray` force2+ $ forStencilFrom2 BoundClamp arr floatOfWord8+ [stencil2| 1 4 6 4 1 |] +++-- | Separable Gaussian blur in the Y direction.+{-# NOINLINE blurSepY #-}+blurSepY :: Array DIM2 Float -> Array DIM2 Float+blurSepY arr@(Array _ [Region RangeAll (GenManifest _)]) + = arr `deepSeqArray` force2+ $ R.map (/ 256)+ $ forStencil2 BoundClamp arr+ [stencil2| 1+ 4+ 6+ 4+ 1 |]+++-- | Compute gradient in the X direction.+{-# NOINLINE gradientX #-}+gradientX :: Array DIM2 Float -> Array DIM2 Float+gradientX img@(Array _ [Region RangeAll (GenManifest _)])+ = img `deepSeqArray` force2+ $ forStencil2 (BoundConst 0) img+ [stencil2| -1 0 1+ -2 0 2+ -1 0 1 |]+++-- | Compute gradient in the Y direction.+{-# NOINLINE gradientY #-}+gradientY :: Array DIM2 Float -> Array DIM2 Float+gradientY img@(Array _ [Region RangeAll (GenManifest _)])+ = img `deepSeqArray` force2+ $ forStencil2 (BoundConst 0) img+ [stencil2| 1 2 1+ 0 0 0+ -1 -2 -1 |] +++-- | Classify the magnitude and orientation of the vector gradient.+{-# NOINLINE gradientMagOrient #-}+gradientMagOrient :: Float -> Image Float -> Image Float -> Image (Float, Float)+gradientMagOrient !threshLow + dX@(Array _ [Region RangeAll (GenManifest _)])+ dY@(Array _ [Region RangeAll (GenManifest _)])++ = [dX, dY] `deepSeqArrays` force2+ $ R.zipWith magOrient dX dY++ where {-# INLINE magOrient #-}+ magOrient :: Float -> Float -> (Float, Float)+ magOrient !x !y+ = (magnitude x y, orientation x y)+ + {-# INLINE magnitude #-}+ magnitude :: Float -> Float -> Float+ magnitude !x !y+ = sqrt (x * x + y * y)+ + {-# INLINE orientation #-}+ orientation :: Float -> Float -> Float+ orientation !x !y++ -- Don't bother computing orientation if vector is below threshold.+ | x >= negate threshLow, x < threshLow+ , y >= negate threshLow, y < threshLow+ = orientUndef++ | otherwise+ = let -- Determine the angle of the vector and rotate it around a bit+ -- to make the segments easier to classify.+ !d = atan2 y x + !dRot = (d - (pi/8)) * (4/pi)+ + -- Normalise angle to beween 0..8+ !dNorm = if dRot < 0 then dRot + 8 else dRot++ -- Doing explicit tests seems to be faster than using the FP floor function.+ in if dNorm >= 4+ then if dNorm >= 6+ then if dNorm >= 7+ then orientHoriz -- 7+ else orientNegDiag -- 6++ else if dNorm >= 5+ then orientVert -- 5+ else orientPosDiag -- 4++ else if dNorm >= 2+ then if dNorm >= 3+ then orientHoriz -- 3+ else orientNegDiag -- 2++ else if dNorm >= 1+ then orientVert -- 1+ else orientPosDiag -- 0+++-- | Suppress pixels that are not local maxima, and use the magnitude to classify maxima+-- into strong and weak (potential) edges.+{-# NOINLINE suppress #-}+suppress :: Float -> Float -> Image (Float, Float) -> Image Word8+suppress threshLow threshHigh+ dMagOrient@(Array shSrc [Region RangeAll (GenManifest _)]) + = dMagOrient `deepSeqArray` force2 + $ makeBordered2 shSrc 1 + (GenCursor id addDim (const 0))+ (GenCursor id addDim comparePts)++ where {-# INLINE comparePts #-}+ comparePts d@(sh :. i :. j)+ | o == orientUndef = edge None+ | o == orientHoriz = isMax (getMag (sh :. i :. j-1)) (getMag (sh :. i :. j+1)) + | o == orientVert = isMax (getMag (sh :. i-1 :. j)) (getMag (sh :. i+1 :. j)) + | o == orientNegDiag = isMax (getMag (sh :. i-1 :. j-1)) (getMag (sh :. i+1 :. j+1)) + | o == orientPosDiag = isMax (getMag (sh :. i-1 :. j+1)) (getMag (sh :. i+1 :. j-1)) + | otherwise = edge None+ + where+ !o = getOrient d + !m = getMag (Z :. i :. j)++ getMag = fst . (R.unsafeIndex dMagOrient)+ getOrient = snd . (R.unsafeIndex dMagOrient)++ {-# INLINE isMax #-}+ isMax intensity1 intensity2+ | m < threshLow = edge None+ | m < intensity1 = edge None+ | m < intensity2 = edge None+ | m < threshHigh = edge Weak+ | otherwise = edge Strong+++-- | Select indices of strong edges.+-- TODO: If would better if we could medge this into the above stage, and record the strong edge+-- during non-maximum suppression, but Repa doesn't provide a fused mapFilter primitive yet.+{-# NOINLINE selectStrong #-}+selectStrong :: Image Word8 -> Array DIM1 Int+selectStrong img@(Array _ [Region RangeAll (GenManifest vec)])+ = img `deepSeqArray` + let {-# INLINE match #-}+ match ix = vec `V.unsafeIndex` ix == edge Strong++ {-# INLINE process #-}+ process ix = ix+ + in select match process (size $ extent img)+++-- | Trace out strong edges in the final image. +-- Also trace out weak edges that are connected to strong edges.+{-# NOINLINE wildfire #-}+wildfire + :: Image Word8 -- ^ Image with strong and weak edges set.+ -> Array DIM1 Int -- ^ Array containing flat indices of strong edges.+ -> Image Word8++wildfire img@(Array _ [Region RangeAll (GenManifest _)]) arrStrong+ = unsafePerformIO + $ do (sh, vec) <- wildfireIO + return $ sh `seq` vec `seq` + Array sh [Region RangeAll (GenManifest vec)]++ where lenImg = R.size $ R.extent img+ lenStrong = R.size $ R.extent arrStrong+ vStrong = toVector arrStrong+ + wildfireIO+ = do -- Stack of image indices we still need to consider.+ vStrong' <- V.thaw vStrong+ vStack <- VM.grow vStrong' (lenImg - lenStrong)+ + -- Burn in new edges.+ vImg <- VM.unsafeNew lenImg+ VM.set vImg 0+ burn vImg vStack lenStrong+ vImg' <- V.unsafeFreeze vImg+ return (extent img, vImg')++ + burn :: VM.IOVector Word8 -> VM.IOVector Int -> Int -> IO ()+ burn !vImg !vStack !top+ | top == 0+ = return ()+ + | otherwise+ = do let !top' = top - 1+ n <- VM.unsafeRead vStack top'+ let (Z :. y :. x) = fromIndex (R.extent img) n++ let {-# INLINE push #-}+ push t = pushWeak vImg vStack t+ + VM.write vImg n (edge Strong)+ >> push (Z :. y - 1 :. x - 1) top'+ >>= push (Z :. y - 1 :. x )+ >>= push (Z :. y - 1 :. x + 1)++ >>= push (Z :. y :. x - 1)+ >>= push (Z :. y :. x + 1)++ >>= push (Z :. y + 1 :. x - 1)+ >>= push (Z :. y + 1 :. x )+ >>= push (Z :. y + 1 :. x + 1)++ >>= burn vImg vStack++ -- If this ix is weak in the source then set it to strong in the+ -- result and push the ix onto the stack.+ {-# INLINE pushWeak #-}+ pushWeak vImg vStack ix top+ = do let n = toIndex (extent img) ix+ xDst <- VM.unsafeRead vImg n+ let xSrc = img `R.unsafeIndex` ix++ if xDst == edge None + && xSrc == edge Weak+ then do+ VM.unsafeWrite vStack top (toIndex (extent img) ix)+ return (top + 1)+ + else return top
+ examples/FFT/FFT2d/src-repa/Main.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE ScopedTypeVariables #-}++-- | Perform the 2D FFT on a BMP image.+import Data.Array.Repa.Algorithms.FFT+import Data.Array.Repa.Algorithms.DFT.Center+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa.IO.BMP+import Data.Array.Repa as A+import System.Environment+import Control.Monad++main :: IO ()+main + = do args <- getArgs+ case args of+ [fileIn, clipMag, fileMag, filePhase]+ -> mainWithArgs fileIn (read clipMag) fileMag filePhase++ _ -> putStr $ unlines+ [ "Usage: repa-fft2d <fileIn.bmp> <clip mag :: Int> <fileOutMag.bmp> <fileOutPhase.bmp>"+ , " "+ , " Image dimensions must be powers of two, eg 128x512 or 64x256"+ , ""+ , " The output magnitude has a high dynamic range. We need to clip it otherwise"+ , " most of the pixels in the output BMP will be black. Start with a value equal"+ , " to about the width of the image (eg 512)"+ , "" ]+ + +mainWithArgs fileIn (clipMag :: Double) fileMag filePhase+ = do + -- Load in the matrix.+ arrReal <- liftM (either (\e -> error $ show e) force)+ $ readMatrixFromGreyscaleBMP fileIn++ let arrComplex = force $ A.map (\r -> (r, 0 :: Double)) arrReal++ -- Apply the centering transform so that the output has the zero+ -- frequency in the middle of the image.+ let arrCentered = center2d arrComplex+ + -- Do the 2d transform.+ let arrFreq = fft2d Forward arrCentered+ + -- Write out the magnitude of the transformed array, + -- clipping it at the given value.+ let clip m = if m >= clipMag then clipMag else m+ let arrMag = force $ A.map (clip . mag) arrFreq+ writeMatrixToGreyscaleBMP fileMag arrMag++ -- Write out the phase of the transofmed array, + -- scaling it to make full use of the 8 bit greyscale.+ let scaledArg x = (arg x + pi) * (255 / (2 * pi))+ let arrPhase = force $ A.map scaledArg arrFreq+ writeMatrixToGreyscaleBMP filePhase arrPhase
+ examples/FFT/HighPass2d/src-repa/Main.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE PatternGuards #-}++-- | Perform high pass filtering on a BMP image.+import Data.Array.Repa.Algorithms.FFT+import Data.Array.Repa.Algorithms.DFT.Center+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa.IO.BMP+import Data.Array.Repa.IO.Timing+import Data.Array.Repa as A+import System.Environment+import Control.Monad+++main :: IO ()+main + = do args <- getArgs+ case args of+ [cutoff, fileIn, fileOut]+ -> mainWithArgs (read cutoff) fileIn fileOut++ _ -> putStr $ unlines+ [ "Usage: repa-fft-highpass <cutoff::Int> <fileIn.bmp> <fileOut.bmp>"+ , ""+ , " Image dimensions must be powers of two, eg 128x512 or 64x256"+ , "" ]+ + +mainWithArgs cutoff fileIn fileOut+ = do + -- Load in the matrix.+ (arrRed, arrGreen, arrBlue)+ <- liftM (either (\e -> error $ show e) id)+ $ readComponentsFromBMP fileIn+ + -- The deepSeqs are to make sure we're just measuring the transform time.+ arrRed + `deepSeqArray` arrGreen+ `deepSeqArray` arrBlue+ `deepSeqArray` return ()+ + -- Do the transform on each component individually+ ((arrRed', arrGreen', arrBlue'), t)+ <- time+ $ let arrRed' = transform cutoff arrRed+ arrGreen' = transform cutoff arrGreen+ arrBlue' = transform cutoff arrBlue+ in arrRed' + `deepSeqArray` arrGreen'+ `deepSeqArray` arrBlue'+ `deepSeqArray` return (arrRed', arrGreen', arrBlue')+ + putStr (prettyTime t)+ + -- Write it back to file.+ writeComponentsToBMP fileOut+ arrRed' arrGreen' arrBlue' + ++transform cutoff arrReal+ = let arrComplex = force $ A.map (\r -> (fromIntegral r, 0)) arrReal+ + -- Do the 2d transform.+ arrCentered = center2d arrComplex+ arrFreq = fft2d Forward arrCentered++ -- Zap out the low frequency components.+ _ :. height :. width = extent arrFreq+ centerX = width `div` 2+ centerY = height `div` 2+ + {-# INLINE highpass #-}+ highpass get ix@(_ :. y :. x)+ | x > centerX + cutoff+ || x < centerX - cutoff+ || y > centerY + cutoff+ || y < centerY - cutoff+ = get ix+ + | otherwise+ = 0+ + arrFilt = traverse arrFreq id highpass++ -- Do the inverse transform to get back to image space.+ arrInv = fft2d Inverse arrFilt+ + -- Get the magnitude of the transformed array, + arrMag = A.map (truncate . mag) arrInv++ in arrMag+
+ examples/FFT/HighPass3d/src-repa/Main.hs view
@@ -0,0 +1,105 @@+{-# LANGUAGE PatternGuards #-}++import Data.Array.Repa.Algorithms.FFT+import Data.Array.Repa.Algorithms.DFT.Center+import Data.Array.Repa.Algorithms.Complex+import Data.Array.Repa.IO.BMP+import Data.Array.Repa.IO.ColorRamp+import Data.Array.Repa.IO.Timing+import Data.Array.Repa as A+import Data.Word+import System.Environment+import Control.Monad+import Prelude as P++main :: IO ()+main + = do args <- getArgs+ case args of+ [size, prefix] -> mainWithArgs (read size) prefix++ _ -> putStr $ unlines+ [ "Usage: repa-fft3d-highpass <size::Int> <prefix>"+ , ""+ , " Size must be a power of two."+ , " You get a stack of prefix###.bmp files resulting from high-pass filtering a cube."+ , "" ]+ + +mainWithArgs size prefixOut+ = do+ -- Generate a cube for initial data.+ let shape = Z :. size :. size :. size+ let cubeSize = size `div` 4+ let center = size `div` 2+ let cutoff = 4++ let arrInit = force + $ fromFunction shape + (\ix -> if isInCenteredCube center cubeSize ix + then (1, 0) else (0, 0))++ arrInit `deepSeqArray` return ()++ (arrFinal, t)+ <- time + $ let arrFinal' = transform arrInit center cutoff+ in arrFinal' `deepSeqArray` return arrFinal'++ putStr (prettyTime t)++ mapM_ (dumpSlice prefixOut arrFinal) [0..size - 1]+++-- | To the high pass transform.+transform+ :: Array DIM3 Complex+ -> Int+ -> Int+ -> Array DIM3 Complex+transform arrInit center cutoff+ = let -- Transform to frequency space.+ arrCentered = center3d arrInit+ arrFreq = fft3d Forward arrCentered+ + -- Zap out the high frequency components+ arrFilt = traverse arrFreq id (highpass center cutoff)+ + -- Do the inverse transform to get back to image space.+ arrInv = fft3d Inverse arrFilt+ in arrInv+++-- | Dump a numbered slice of this array to a BMP file.+dumpSlice + :: FilePath+ -> Array DIM3 Complex+ -> Int+ -> IO ()++dumpSlice prefix arr sliceNum+ = do let arrSlice = slice arr (Any :. sliceNum :. All)+ let arrGrey = A.map (truncate . (* 255) . mag) arrSlice+ let fileName = prefix P.++ (pad0 3 (show sliceNum)) P.++ ".bmp"++ writeComponentsToBMP fileName+ arrGrey arrGrey arrGrey++pad0 n str+ = P.replicate (n - length str) '0' P.++ str+++{-# INLINE isInCenteredCube #-}+isInCenteredCube center cutoff ix@(_ :. z :. y :. x)+ = let high = center + cutoff+ low = center - cutoff+ in x >= low && x <= high+ && y >= low && y <= high+ && z >= low && z <= high+++{-# INLINE highpass #-}+highpass center cutoff get ix+ | isInCenteredCube center cutoff ix = 0+ | otherwise = get ix+
+ examples/Laplace/src-repa/Main.hs view
@@ -0,0 +1,181 @@+{-# LANGUAGE BangPatterns #-}++-- | Solver for the Laplace equation+-- You supply a BMP files specifying the boundary conditions.+-- The output is written back to another BMP file.+--+import SolverGet as SG+import SolverStencil as SS+import Data.Array.Repa as A+import Data.Array.Repa.IO.BMP +import Data.Array.Repa.IO.ColorRamp+import Data.Array.Repa.IO.Timing+import System.Environment+import Data.Word+import Control.Monad+import Prelude as P++type Solver + = Int -- ^ Number of iterations to use.+ -> Array DIM2 Double -- ^ Boundary value mask.+ -> Array DIM2 Double -- ^ Boundary values.+ -> Array DIM2 Double -- ^ Initial state.+ -> Array DIM2 Double++solvers + = [ ("get", SG.solveLaplace)+ , ("stencil", SS.solveLaplace) ]+++main :: IO ()+main + = do args <- getArgs+ case args of+ [strSolver, steps, fileInput, fileOutput] + -> let Just solver = lookup strSolver solvers+ in laplace solver (read steps) fileInput fileOutput++ _ -> do+ putStr usage+ return ()+++-- | Command line usage information.+usage :: String+usage = unlines+ [ "Usage: laplace <solver> <iterations> <input.bmp> <output.bmp>"+ , ""+ , " iterations :: Int Number of iterations to use in the solver."+ , " input.bmp :: FileName Uncompressed RGB24 or RGBA32 BMP file for initial and boundary values."+ , " output.bmp :: FileName BMP file to write output to."+ , "" + , " solver = one of " P.++ show (P.map fst solvers)+ , ""+ , " Format of input file:"+ , " Boundary values are indicated in greyscale,"+ , " ie from the list [(x, x, x) | x <- [0 .. 255]]"+ , " Non-boundary values are indicated in blue,"+ , " ie (0, 0, 255)"+ , " Any other pixel value is an error." + , ""+ , " NOTE: For GHC 7.0.3, this runs better when you turn off the parallel"+ , " garbage collector. Run with +RTS -qg"+ , "" ]+ ++-- | Solve it.+laplace :: Solver+ -> Int -- ^ Number of iterations to use.+ -> FilePath -- ^ Input file.+ -> FilePath -- ^ Output file+ -> IO ()++laplace solve steps fileInput fileOutput+ = do+ -- Load up the file containing boundary conditions.+ arrImage <- liftM (either (error . show) force)+ $ readImageFromBMP fileInput+ + arrImage `deepSeqArray` return ()+ let arrBoundValue = force $ slurpDoublesFromImage slurpBoundValue arrImage+ let arrBoundMask = force $ slurpDoublesFromImage slurpBoundMask arrImage+ arrBoundValue `deepSeqArray` return ()+ arrBoundMask `deepSeqArray` return ()++ let arrInitial = arrBoundValue + + -- Run the solver.+ (arrFinal, t)+ <- time+ $ let arrFinal = solve steps arrBoundMask arrBoundValue arrInitial+ in arrFinal `deepSeqArray` return arrFinal++ -- Print how long it took+ putStr (prettyTime t)++ -- Make the result image+ let arrImageOut + = makeImageFromDoubles (rampColorHotToCold 0.0 1.0) arrFinal++ -- Write out the image to file. + writeImageToBMP+ fileOutput+ arrImageOut++++slurpDoublesFromImage+ :: (Word8 -> Word8 -> Word8 -> Double)+ -> Array DIM3 Word8+ -> Array DIM2 Double+ +{-# INLINE slurpDoublesFromImage #-}+slurpDoublesFromImage mkDouble + arrBound@(Array _ [Region _ (GenManifest _)])+ = arrBound `deepSeqArray` force+ $ unsafeTraverse arrBound+ (\(Z :. height :. width :. _) + -> Z :. height :. width)++ (\get (Z :. y :. x)+ -> mkDouble+ (get (Z :. y :. x :. 0))+ (get (Z :. y :. x :. 1))+ (get (Z :. y :. x :. 2)))+++makeImageFromDoubles+ :: (Double -> (Double, Double, Double))+ -> Array DIM2 Double+ -> Array DIM3 Word8+ +{-# INLINE makeImageFromDoubles #-}+makeImageFromDoubles fnColor + arrDoubles@(Array _ [Region _ (GenManifest _)])+ = arrDoubles `deepSeqArray` force+ $ unsafeTraverse arrDoubles+ (\(Z :. height :. width)+ -> Z :. height :. width :. 4)+ + (\get (Z :. y :. x :. c)+ -> let (r, g, b) = fnColor (get (Z :. y :. x))+ in case c of+ 0 -> fromIntegral $ (truncate (r * 255) :: Int)+ 1 -> fromIntegral $ (truncate (g * 255) :: Int)+ 2 -> fromIntegral $ (truncate (b * 255) :: Int)+ 3 -> 0)++++-- | Extract the boundary value from a RGB triple.+slurpBoundValue :: Word8 -> Word8 -> Word8 -> Double+{-# INLINE slurpBoundValue #-}+slurpBoundValue r g b+ -- A non-boundary value.+ | r == 0 && g == 0 && b == 255 + = 0++ -- A boundary value.+ | (r == g) && (r == b) + = fromIntegral (fromIntegral r :: Int) / 255+ + | otherwise+ = error $ "Unhandled pixel value in input " P.++ show (r, g, b)+++-- | Extract boundary mask from a RGB triple.+slurpBoundMask :: Word8 -> Word8 -> Word8 -> Double+{-# INLINE slurpBoundMask #-}+slurpBoundMask r g b+ -- A non-boundary value.+ | r == 0 && g == 0 && b == 255 + = 1++ -- A boundary value.+ | (r == g) && (r == b) + = 0+ + | otherwise+ = error $ "Unhandled pixel value in input " P.++ show (r, g, b)+ +
+ examples/Laplace/src-repa/SolverGet.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE BangPatterns #-}+module SolverGet+ (solveLaplace)+where +import Data.Array.Repa as A+import qualified Data.Array.Repa.Shape as S++-- | Solver for the Laplace equation.+solveLaplace+ :: Int -- ^ Number of iterations to use.+ -> Array DIM2 Double -- ^ Boundary value mask.+ -> Array DIM2 Double -- ^ Boundary values.+ -> Array DIM2 Double -- ^ Initial state.+ -> Array DIM2 Double++{-# NOINLINE solveLaplace #-}+solveLaplace steps arrBoundMask arrBoundValue arrInit+ = go steps arrInit+ where go !i !arr+ | i == 0 = arr+ | otherwise + = let arr' = relaxLaplace arrBoundMask arrBoundValue arr+ in arr' `deepSeqArray` go (i - 1) arr'+++-- | Perform matrix relaxation for the Laplace equation,+-- using a stencil function.+--+-- Computation fn is+-- u'(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1)) / 4+--+-- Apply the boundary conditions to this matrix.+-- The mask matrix has 0 in places where boundary conditions hold+-- and 1 otherwise.+--+-- The value matrix has the boundary condition value in places where it holds,+-- and 0 otherwise.+-- +relaxLaplace+ :: Array DIM2 Double -- ^ Boundary condition mask+ -> Array DIM2 Double -- ^ Boundary condition values+ -> Array DIM2 Double -- ^ Initial matrix+ -> Array DIM2 Double ++{-# INLINE relaxLaplace #-}+relaxLaplace + arrBoundMask@(Array _ [Region RangeAll (GenManifest _)])+ arrBoundValue@(Array _ [Region RangeAll (GenManifest _)])+ arr@(Array _ [Region RangeAll (GenManifest _)])+ = [arrBoundMask, arrBoundValue, arr] `deepSeqArrays` force2+ $ A.zipWith (+) arrBoundValue+ $ A.zipWith (*) arrBoundMask+ $ unsafeTraverse arr id elemFn+ where+ _ :. height :. width + = extent arr++ {-# INLINE elemFn #-}+ elemFn !get !d@(sh :. i :. j)+ = if isBorder i j+ then get d+ else (get (sh :. (i-1) :. j)+ + get (sh :. i :. (j-1))+ + get (sh :. (i+1) :. j)+ + get (sh :. i :. (j+1))) / 4++ -- Check if this element is on the border of the matrix.+ -- If so we can't apply the stencil because we don't have all the neighbours.+ {-# INLINE isBorder #-}+ isBorder !i !j+ = (i == 0) || (i >= width - 1) + || (j == 0) || (j >= height - 1) ++
+ examples/Laplace/src-repa/SolverStencil.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE BangPatterns, TemplateHaskell, QuasiQuotes #-}+module SolverStencil+ (solveLaplace)+where +import Data.Array.Repa as A+import Data.Array.Repa.Stencil as A+import Data.Array.Repa.Algorithms.Iterate as A+import qualified Data.Array.Repa.Shape as S+import Language.Haskell.TH+import Language.Haskell.TH.Quote++-- | Solver for the Laplace equation.+solveLaplace+ :: Int -- ^ Number of iterations to use.+ -> Array DIM2 Double -- ^ Boundary value mask.+ -> Array DIM2 Double -- ^ Boundary values.+ -> Array DIM2 Double -- ^ Initial state.+ -> Array DIM2 Double++{-# NOINLINE solveLaplace #-}+solveLaplace !steps !arrBoundMask !arrBoundValue !arrInit+ = go steps arrInit+ where go 0 !arr = arr+ go n !arr = go (n - 1) (relaxLaplace arrBoundMask arrBoundValue arr)+ ++{-# INLINE relaxLaplace #-}+relaxLaplace+ :: Array DIM2 Double -- ^ Boundary value mask.+ -> Array DIM2 Double -- ^ Boundary values.+ -> Array DIM2 Double -- ^ Initial state.+ -> Array DIM2 Double++relaxLaplace + arrBoundMask@(Array _ [Region RangeAll (GenManifest _)])+ arrBoundValue@(Array _ [Region RangeAll (GenManifest _)])+ arr@(Array _ [Region RangeAll (GenManifest _)])+ = [arrBoundMask, arrBoundValue, arr] `deepSeqArrays` + let ex = extent arr+ arrBoundMask' = reshape ex arrBoundMask+ arrBoundValue' = reshape ex arrBoundValue+ arr' = reshape ex arr+ in force+ $ A.zipWith (+) arrBoundValue'+ $ A.zipWith (*) arrBoundMask'+ $ A.map (/ 4)+ $ mapStencil2 BoundClamp+ [stencil2| 0 1 0+ 1 0 1 + 0 1 0 |] arr'+
+ examples/MMult/src-repa/Main.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE PatternGuards, PackageImports #-}++import Data.Array.Repa as A+import Data.Array.Repa.IO.Matrix+import Data.Array.Repa.IO.Timing+import Data.Array.Repa.Algorithms.Matrix+import Data.Array.Repa.Algorithms.Randomish+import Data.Maybe+import System.Environment+import Control.Monad+import System.Random+import Prelude as P++-- Arg Parsing ------------------------------------------------------------------------------------+data Arg+ = ArgSolver String+ | ArgMatrixRandom Int Int+ | ArgMatrixFile FilePath+ | ArgOutFile FilePath+ deriving Show++isArgMatrix arg+ = case arg of+ ArgMatrixRandom{} -> True+ ArgMatrixFile{} -> True+ _ -> False++parseArgs [] = []+parseArgs (flag:xx)+ | "-file" <- flag+ , f:rest <- xx+ = ArgMatrixFile f : parseArgs rest++ | "-out" <- flag+ , f:rest <- xx+ = ArgOutFile f : parseArgs rest+ + | "-random" <- flag+ , x:y:rest <- xx+ = ArgMatrixRandom (read x) (read y) : parseArgs rest+ + | otherwise + = error $ "bad arg " P.++ flag P.++ "\n"++printHelp+ = putStr + $ unlines+ [ "Usage: mmult [args..]"+ , ""+ , " -random <height> <width> Use a random matrix of this size."+ , " -file <filename> Read a matrix from this file."+ , " -out <filename> Write resulting matrix to this file."+ , ""+ , " Format of matrix file:"+ , " MATRIX"+ , " <width> <height>"+ , " <whitespace separated values..>"+ , "" ]+++-- | Get a matrix from a file, or generate a random one.+getMatrix :: Arg -> IO (Array DIM2 Double)+getMatrix arg+ = case arg of+ ArgMatrixFile fileName + -> readMatrixFromTextFile fileName++ ArgMatrixRandom height width + -> return $ randomishDoubleArray (Z :. height :. width) (-100) 100 12345+++-- Main -------------------------------------------------------------------------------------------+main :: IO ()+main + = do args <- liftM parseArgs $ getArgs+ main' args++main' args+ | [argMat1, argMat2] <- filter isArgMatrix args+ , mArgOut <- listToMaybe [s | ArgOutFile s <- args]+ = do + -- Get matrices from files, + -- or generate random ones we were asked to.+ mat1 <- getMatrix argMat1+ mat2 <- getMatrix argMat2++ mat1 `deepSeqArray` mat2 `deepSeqArray` return ()+ + -- Run the solver.+ (matResult, t) + <- time + $ let matResult = multiplyMM mat1 mat2+ in matResult `deepSeqArray` return matResult++ -- Print how long it took.+ putStr (prettyTime t)++ -- Print a checksum of all the elements+ putStrLn $ "checkSum = " P.++ show (A.sumAll matResult)++ -- Write the output to file if requested.+ case mArgOut of + Nothing -> return ()+ Just fileOut -> writeMatrixToTextFile fileOut matResult+ + | otherwise+ = printHelp+
+ examples/Sobel/src-repa/Main.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE PackageImports, BangPatterns, QuasiQuotes, PatternGuards #-}+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}++-- | Apply Sobel operators to an image.+import Data.Word+import Control.Monad+import System.Environment+import Data.Array.Repa as Repa+import Data.Array.Repa.IO.BMP+import Data.Array.Repa.IO.Timing+import Prelude hiding (compare)++import Solver++-- Main routine ---------------------------------------------------------------+main + = do args <- getArgs+ case args of+ [iterations, fileIn, fileOut] + -> run (read iterations) fileIn fileOut+ _ -> putStrLn "Usage: sobel <iterations::Int> <fileIn.bmp> <fileOut.bmp>"+++run iterations fileIn fileOut+ = do inputImage <- liftM (force . either (error . show) id) + $ readImageFromBMP fileIn+ + let greyImage = toGreyScale inputImage+ greyImage `deepSeqArray` return ()+ + (result, tElapsed)+ <- time $ let (gX, gY) = loop iterations greyImage+ in gX `deepSeqArray` gY `deepSeqArray` return (gX, gY)++ putStr $ prettyTime tElapsed+ + let (gX, gY) = result+ let outImage = force2 $ Repa.zipWith magnitude gX gY ++ outImage `seq` return ()++ -- TODO: The image normalization in this write fn eats up most of the runtime.+ writeMatrixToGreyscaleBMP fileOut outImage+++loop :: Int -> Image -> (Image, Image)+loop n + = withManifest $ \img ->+ if n == 0+ then (img, img)+ else do + let gX = gradientX img+ let gY = gradientY img + if (n == 1) + then gX `deepSeqArray` gY `deepSeqArray` (gX, gY)+ else gX `deepSeqArray` gY `deepSeqArray` loop (n - 1) img+++-- | Determine the squared magnitude of a vector.+magnitude :: Float -> Float -> Double+{-# INLINE magnitude #-}+magnitude x y+ = fromRational $ toRational $ sqrt (x * x + y * y)+++-- | RGB to greyscale conversion.+toGreyScale :: Array DIM3 Word8 -> Image+{-# NOINLINE toGreyScale #-}+toGreyScale + = withManifest $ \arr ->+ arr `seq` force2 $ traverse arr+ (\(sh :. _) -> sh)+ (\get ix -> rgbToLuminance + (get (ix :. 0))+ (get (ix :. 1))+ (get (ix :. 2)))+++-- | Convert a RGB value to a luminance.+rgbToLuminance :: Word8 -> Word8 -> Word8 -> Float+{-# INLINE rgbToLuminance #-}+rgbToLuminance r g b + = fromIntegral r * 0.3+ + fromIntegral g * 0.59+ + fromIntegral b * 0.11
repa-examples.cabal view
@@ -1,5 +1,5 @@ Name: repa-examples-Version: 1.1.1.0+Version: 2.0.0.1 License: BSD3 License-file: LICENSE Author: The DPH Team@@ -16,59 +16,140 @@ Synopsis: Examples using the Repa array library. ++Executable repa-mmult+ Build-depends: + base == 4.*,+ repa == 2.0.*,+ repa-io == 2.0.*,+ repa-algorithms == 2.0.*,+ random == 1.0.*++ Main-is: examples/MMult/src-repa/Main.hs+ hs-source-dirs: examples/MMult/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30+++Executable repa-laplace+ Build-depends: + base == 4.*,+ repa == 2.0.*,+ repa-io == 2.0.*++ Main-is: examples/Laplace/src-repa/Main.hs+ other-modules: SolverGet SolverStencil+ hs-source-dirs: examples/Laplace/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30++ Executable repa-fft2d Build-depends: base == 4.*,- dph-prim-par == 0.4.*,- dph-base == 0.4.*,- repa == 1.1.*,- repa-algorithms == 1.1.*,- repa-io == 1.1.*+ repa == 2.0.*,+ repa-algorithms == 2.0.*,+ repa-io == 2.0.* - Main-is: FFT/src/FFT2d/Main.hs- hs-source-dirs: FFT/src .- ghc-options: -Odph -threaded + Main-is: examples/FFT/FFT2d/src-repa/Main.hs+ hs-source-dirs: examples/FFT/FFT2d/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30 Executable repa-fft2d-highpass Build-depends: base == 4.*,- dph-prim-par == 0.4.*,- dph-base == 0.4.*,- repa == 1.1.*,- repa-algorithms == 1.1.*,- repa-io == 1.1.*+ repa == 2.0.*,+ repa-algorithms == 2.0.*,+ repa-io == 2.0.* - Main-is: FFT/src/HighPass/Main.hs- hs-source-dirs: FFT/src .- ghc-options: -Odph -threaded + Main-is: examples/FFT/HighPass2d/src-repa/Main.hs+ hs-source-dirs: examples/FFT/HighPass2d/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30 -Executable repa-laplace+Executable repa-fft3d-highpass Build-depends: base == 4.*,- dph-prim-par == 0.4.*,- dph-base == 0.4.*,- repa == 1.1.*,- repa-io == 1.1.*+ repa == 2.0.*,+ repa-algorithms == 2.0.* - Main-is: Laplace/src/Main.hs- other-modules: Solver- hs-source-dirs: Laplace/src .- ghc-options: -Odph -threaded + Main-is: examples/FFT/HighPass3d/src-repa/Main.hs+ hs-source-dirs: examples/FFT/HighPass3d/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30 -Executable repa-mmult+Executable repa-blur Build-depends: base == 4.*,- dph-prim-par == 0.4.*,- dph-base == 0.4.*,- repa == 1.1.*,- repa-io == 1.1.*,- repa-algorithms == 1.1.*,- random == 1.0.*+ repa == 2.0.*,+ repa-algorithms == 2.0.*,+ vector == 0.7.* - Main-is: MMult/src/Main.hs- hs-source-dirs: MMult/src .- ghc-options: -Odph -threaded + Main-is: examples/Blur/src-repa/Main.hs+ hs-source-dirs: examples/Blur/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30+++Executable repa-sobel+ Build-depends: + base == 4.*,+ repa == 2.0.*,+ repa-algorithms == 2.0.*,+ template-haskell >= 2.5 && < 2.6++ Main-is: examples/Sobel/src-repa/Main.hs+ hs-source-dirs: examples/Sobel/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold30+++Executable repa-canny+ Build-depends: + base == 4.*,+ repa == 2.0.*,+ repa-algorithms == 2.0.*,+ template-haskell >= 2.5 && < 2.6++ Main-is: examples/Canny/src-repa/Main.hs+ hs-source-dirs: examples/Canny/src-repa .+ ghc-options: + -threaded + -rtsopts + -Odph -fllvm -optlo-O3+ -fno-liberate-case+ -funfolding-use-threshold100+ -funfolding-keeness-factor100