diff --git a/FFT/src/FFT2d/Main.hs b/FFT/src/FFT2d/Main.hs deleted file mode 100644 --- a/FFT/src/FFT2d/Main.hs +++ /dev/null @@ -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 " - , "" - , " 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 - diff --git a/FFT/src/HighPass/Main.hs b/FFT/src/HighPass/Main.hs deleted file mode 100644 --- a/FFT/src/HighPass/Main.hs +++ /dev/null @@ -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 " - , "" ] - - -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 - diff --git a/Laplace/src/Main.hs b/Laplace/src/Main.hs deleted file mode 100644 --- a/Laplace/src/Main.hs +++ /dev/null @@ -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 :: 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) - - diff --git a/Laplace/src/Solver.hs b/Laplace/src/Solver.hs deleted file mode 100644 --- a/Laplace/src/Solver.hs +++ /dev/null @@ -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 - - diff --git a/MMult/src/Main.hs b/MMult/src/Main.hs deleted file mode 100644 --- a/MMult/src/Main.hs +++ /dev/null @@ -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 Use a random matrix of this size." - , " -file Read a matrix from this file." - , " -out Write resulting matrix to this file." - , "" - , " Format of matrix file:" - , " MATRIX" - , " " - , " " - , "" ] - - --- | 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 - - diff --git a/examples/Blur/src-repa/Main.hs b/examples/Blur/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/Blur/src-repa/Main.hs @@ -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 " ] + +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] +-} + diff --git a/examples/Canny/src-repa/Main.hs b/examples/Canny/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/Canny/src-repa/Main.hs @@ -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 [ ] " + + +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 diff --git a/examples/FFT/FFT2d/src-repa/Main.hs b/examples/FFT/FFT2d/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/FFT/FFT2d/src-repa/Main.hs @@ -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 " + , " " + , " 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 diff --git a/examples/FFT/HighPass2d/src-repa/Main.hs b/examples/FFT/HighPass2d/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/FFT/HighPass2d/src-repa/Main.hs @@ -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 " + , "" + , " 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 + diff --git a/examples/FFT/HighPass3d/src-repa/Main.hs b/examples/FFT/HighPass3d/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/FFT/HighPass3d/src-repa/Main.hs @@ -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 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 + diff --git a/examples/Laplace/src-repa/Main.hs b/examples/Laplace/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/Laplace/src-repa/Main.hs @@ -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 " + , "" + , " 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) + + diff --git a/examples/Laplace/src-repa/SolverGet.hs b/examples/Laplace/src-repa/SolverGet.hs new file mode 100644 --- /dev/null +++ b/examples/Laplace/src-repa/SolverGet.hs @@ -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) + + diff --git a/examples/Laplace/src-repa/SolverStencil.hs b/examples/Laplace/src-repa/SolverStencil.hs new file mode 100644 --- /dev/null +++ b/examples/Laplace/src-repa/SolverStencil.hs @@ -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' + diff --git a/examples/MMult/src-repa/Main.hs b/examples/MMult/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/MMult/src-repa/Main.hs @@ -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 Use a random matrix of this size." + , " -file Read a matrix from this file." + , " -out Write resulting matrix to this file." + , "" + , " Format of matrix file:" + , " MATRIX" + , " " + , " " + , "" ] + + +-- | 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 + diff --git a/examples/Sobel/src-repa/Main.hs b/examples/Sobel/src-repa/Main.hs new file mode 100644 --- /dev/null +++ b/examples/Sobel/src-repa/Main.hs @@ -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 " + + +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 diff --git a/repa-examples.cabal b/repa-examples.cabal --- a/repa-examples.cabal +++ b/repa-examples.cabal @@ -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