diff --git a/FFT/src/DFT.hs b/FFT/src/DFT.hs deleted file mode 100644 --- a/FFT/src/DFT.hs +++ /dev/null @@ -1,61 +0,0 @@ -{-# LANGUAGE TypeOperators #-} - -module DFT - ( dft - , idft - , dftWithRoots) -where -import Data.Array.Repa as A -import Data.Ratio -import StrictComplex -import Roots - --- | Compute the (non-fast) Discrete Fourier Transform of a vector. -dft :: Shape sh - => Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -dft v - = let rofu = calcRofu (extent v) - in force $ dftWithRoots rofu v - - --- | Compute the inverse (non-fast) Discrete Fourier Transform of a vector. -idft :: Shape sh - => Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -idft v - = let _ :. len = extent v - scale = fromIntegral len :*: 0 - rofu = calcInverseRofu (extent v) - in force $ A.map (/ scale) $ dftWithRoots rofu v - - --- | Generic function for computation of forward or inverse Discrete Fourier Transforms. -dftWithRoots - :: Shape sh - => Array (sh :. Int) Complex -- ^ Roots of unity for this vector length. - -> Array (sh :. Int) Complex -- ^ Input vector. - -> Array (sh :. Int) Complex - -dftWithRoots rofu arr - = traverse arr id (\_ k -> dftK rofu arr k) - - --- | Compute one value of the DFT. -dftK :: Shape sh - => Array (sh :. Int) Complex -- ^ Roots of unity for this vector length. - -> Array (sh :. Int) Complex -- ^ Input vector. - -> (sh :. Int) -- ^ Index of the value we want. - -> Complex - -dftK rofu arrX (_ :. k) - = A.sumAll $ A.zipWith (*) arrX wroots - where sh@(_ :. len) = extent arrX - - -- All the roots we need to multiply with. - wroots = fromFunction sh elemFn - elemFn (sh :. n) - = rofu !: (sh :. (k * n) `mod` len) - diff --git a/FFT/src/FFT.hs b/FFT/src/FFT.hs deleted file mode 100644 --- a/FFT/src/FFT.hs +++ /dev/null @@ -1,136 +0,0 @@ -{-# LANGUAGE TypeOperators, PatternGuards #-} - --- | Computation of Fast Fourier Transforms using the Cooley-Tuckey algorithm. -module FFT - ( fft - , ifft - , fft2d - , fftWithRoots ) -where -import Data.Array.Repa as A -import Data.Ratio -import StrictComplex -import Roots - --- Vector Transform ------------------------------------------------------------------------------- --- | Compute the (fast) Discrete Fourier Transform of a vector. -fft :: Shape sh - => Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -fft v - = let rofu = calcRofu (extent v) - in force $ fftWithRoots rofu v - - --- | Compute the (fast) Inverse Discrete Fourier Transform of a vector. -ifft :: Shape sh - => Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -ifft v - = let _ :. len = extent v - scale = fromIntegral len :*: 0 - rofu = calcInverseRofu (extent v) - in force $ A.map (/ scale) $ fftWithRoots rofu v - - --- Matrix Transform ------------------------------------------------------------------------------- --- | Compute the (fast) Discrete Fourier Transform of a square matrix. -fft2d :: Array DIM2 Complex - -> Array DIM2 Complex - -fft2d arr - = let rofu = calcRofu (extent arr) - fftTrans = transpose . fftWithRoots rofu - in fftTrans $ fftTrans arr - - --- Cube Transform --------------------------------------------------------------------------------- --- | Compute the (fast) Discrete Fourier Transform of a 3d array. -fft3d :: Array DIM3 Complex - -> Array DIM3 Complex - -fft3d arr - = let rofu = calcRofu (extent arr) - - transpose3 arr - = traverse arr - (\(Z :. k :. l :. m) -> (Z :. l :. m :. k)) - (\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) - - fftTrans = transpose3 . fftWithRoots rofu - - in fftTrans $ fftTrans $ fftTrans arr - - --- Worker ----------------------------------------------------------------------------------------- --- | Generic function for computation of forward or inverse Discrete Fourier Transforms. --- The length of the roots vector must be the same as the values vector. --- The length of these vectors must be a power of two. -fftWithRoots - :: Shape sh - => Array (sh :. Int) Complex -- ^ Roots of unity. - -> Array (sh :. Int) Complex -- ^ Input values. - -> Array (sh :. Int) Complex - -fftWithRoots rofu v - | not $ (denominator $ toRational (logBase 2 $ fromIntegral vLen)) == 1 - = error $ "fft: vector length of " ++ show vLen ++ " is not a power of 2" - - | rLen /= vLen - = error $ "fft: length of vector is not the length of the roots" - - | otherwise - = fftWithRoots' rofu v - - where _ :. rLen = extent rofu - _ :. vLen = extent v - - -fftWithRoots' - :: Shape sh - => Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -> Array (sh :. Int) Complex - -{-# INLINE fftWithRoots' #-} -fftWithRoots' rofu v - = case extent v of - _ :. 2 -> fft_two rofu v - dim -> fft_split rofu v dim - -{-# INLINE fft_two #-} -fft_two rofu v - = let vFn' vFn (sh :. 0) = vFn (sh :. 0) + vFn (sh :. 1) - vFn' vFn (sh :. 1) = vFn (sh :. 0) - vFn (sh :. 1) - in traverse v id vFn' - -{-# INLINE fft_split #-} -fft_split rofu v vLen - = let fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v) - - fft_l = traverse2 fft_lr rofu - (\(sh :. 2 :. n) _ -> sh :. n) - (\f r (sh :. i) -> f (sh :. 0 :. i) + r (sh :. i) * f (sh :. 1 :. i)) - - fft_r = traverse2 fft_lr rofu - (\(sh :. 2 :. n) _ -> sh :. n) - (\f r (sh :. i) -> f (sh :. 0 :. i) - r (sh :. i) * f (sh :. 1 :. i)) - - in fft_l +:+ fft_r - -{-# INLINE splitRofu #-} -splitRofu rofu - = traverse rofu - (\(rSh :. rLen) -> rSh :. (2::Int) :. (rLen `div` 2)) - (\rFn (sh :. _ :. i) -> rFn (sh :. 2*i)) - -{-# INLINE splitVector #-} -splitVector v - = let vFn' vFn (sh :. 0 :. i) = vFn (sh :. 2*i) - vFn' vFn (sh :. 1 :. i) = vFn (sh :. 2*i+1) - in traverse v - (\(vSh :. vLen) -> vSh :. 2 :. (vLen `div` 2)) - vFn' - diff --git a/FFT/src/FFT2d/Main.hs b/FFT/src/FFT2d/Main.hs new file mode 100644 --- /dev/null +++ b/FFT/src/FFT2d/Main.hs @@ -0,0 +1,52 @@ + +-- | 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 new file mode 100644 --- /dev/null +++ b/FFT/src/HighPass/Main.hs @@ -0,0 +1,73 @@ +{-# 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/FFT/src/Main.hs b/FFT/src/Main.hs deleted file mode 100644 --- a/FFT/src/Main.hs +++ /dev/null @@ -1,194 +0,0 @@ -{-# LANGUAGE ParallelListComp, PatternGuards, ScopedTypeVariables #-} - -module Main where - -import FFT -import DFT -import Roots -import StrictComplex -import Vector -import Matrix -import PPM -import ColorRamp - -import Data.Array.Repa as A -import Data.List as L -import Data.Maybe -import Prelude as P -import Control.Monad -import System.Environment - - --- Arg Parsing ------------------------------------------------------------------------------------ -data Arg - -- | Use DFT instead of FFT - = ArgDFT - - -- | Compute the inverse transform. - | ArgInverse - - -- | Use a vector step function for input. - | ArgVectorRealStep - Int -- length of 'on' part - Int -- total length of vector - - -- | Read the a real valued vector from this file. - | ArgVectorReal - FilePath - - -- | Read a PPM file as input. - | ArgPPM - FilePath - - -- | Write the result to this file. - | ArgOutMagnitude - FilePath - - -- | Write the magnitude of the transformed matrix as a PPM file. - | ArgOutPPMMagnitude - FilePath - - -- | Clip tranformed values to this level. - | ArgOutPPMClip - Double - - deriving (Show, Eq) - - -parseArgs [] = [] -parseArgs (flag:xx) - | "-dft" <- flag - = ArgDFT : parseArgs xx - - | "-inverse" <- flag - = ArgInverse : parseArgs xx - - | "-vector-real-step" <- flag - , onlen:len:rest <- xx - = ArgVectorRealStep (read onlen) (read len) : parseArgs rest - - | "-vector-real" <- flag - , fileName:rest <- xx - = ArgVectorReal fileName : parseArgs rest - - | "-ppm" <- flag - , fileName:rest <- xx - = ArgPPM fileName : parseArgs rest - - | "-ppm-clip" <- flag - , num:rest <- xx - = ArgOutPPMClip (read num) : parseArgs rest - - | "-out-magnitude" <- flag - , fileName:rest <- xx - = ArgOutMagnitude fileName : parseArgs rest - - | "-out-ppm-magnitude" <- flag - , fileName:rest <- xx - = ArgOutPPMMagnitude fileName : parseArgs rest - - | "-out-ppm-clip" <- flag - , level:rest <- xx - = ArgOutPPMClip (read level) : parseArgs rest - - | otherwise - = error $ "bad arg " ++ flag ++ "\n" - - -help = unlines - [ "Usage: fft [args..]" - , "" - , " -dft Use (slow) DFT instead of (fast) FFT." - , " -inverse Compute the inverse transform." - , "" - , "INPUT:" - , " -vector-real Read a real valued input vector from this file." - , " -vector-real-step Use a real valued step function for input." - , " -ppm Use a PPM file for input." - , "" - , "OUTPUT:" - , " -out-vector-magnitude Write the magnitute of the transformed vector to file." - , " -out-ppm-magnitude Write transformed matrix to a ppm file." - , " -out-ppm-clip ... while clipping transformed values to this level." - , "" - , "NOTE:" - , " - For the fast algorithm, the length of the input vector/dimensions of the array" - , " must be powers of two." - , "" - , " - When using PPM input files, pixels are converted to grey-scale and then used as" - , " the real values of the input matrix." ] - --- Main ------------------------------------------------------------------------------------------- -main :: IO () -main - = do args <- liftM parseArgs $ getArgs - - -- Decide which algorithm to use - let alg | elem ArgDFT args - , elem ArgInverse args = idft - - | elem ArgDFT args = dft - - | elem ArgInverse args = ifft - - | otherwise = fft - - main' args alg - - -main' args alg - - -- | Transform a real-valued step function - | [(onLength, vecLength)] <- [(ol, vl) | ArgVectorRealStep ol vl <- args] - = let offLength = vecLength - onLength - step_real = P.replicate onLength 1 ++ P.replicate offLength 0 - step = P.map (:*: 0) step_real - - arr = A.fromList (Z :. vecLength) step - arrT = alg arr - in outVector args arrT - - -- | Transform some vector from a file - | [fileName] <- [f | ArgVectorReal f <- args] - = do arr_real :: Array DIM1 Double <- readVectorFromTextFile fileName - let arr = A.map (\r -> r :*: 0) arr_real - let arrT = alg arr - outVector args arrT - - -- | Transform a PPM file. - | [fileName] <- [f | ArgPPM f <- args] - = do let loadPixel r g b = sqrt (fromIntegral r^2 + fromIntegral g^2 + fromIntegral b^2) - arr_double <- readPPMAsMatrix loadPixel fileName - let arr_real = (A.map (\r -> r :*: 0) arr_double) :: Array DIM2 Complex - let arrT = fft2d arr_real - - outPPM args arrT - - -- Not sure what you mean... - | otherwise - = putStr help - - -outVector args vec - | [fileName] <- [f | ArgOutMagnitude f <- args ] - = writeVectorAsTextFile (A.map mag vec) fileName - - | otherwise - = return () - -outPPM :: [Arg] -> Array DIM2 Complex -> IO () -outPPM args arr - | [fileName] <- [f | ArgOutPPMMagnitude f <- args ] - , mClipLevel <- listToMaybe [l | ArgOutPPMClip l <- args] - = do let arr_mag = A.map mag arr - let arr_clipped = maybe arr_mag - (\level -> A.map (\x -> if x > level then level else x) arr_mag) - mClipLevel - - writeMatrixAsNormalisedPPM fileName - pixelGrey arr_clipped - - -pixelColor x = (x, x, x) -pixelGrey x = (x, x, x) - diff --git a/FFT/src/Roots.hs b/FFT/src/Roots.hs deleted file mode 100644 --- a/FFT/src/Roots.hs +++ /dev/null @@ -1,43 +0,0 @@ -{-# LANGUAGE TypeOperators #-} - -module Roots - ( calcRofu - , calcInverseRofu) -where -import Data.Array.Repa -import StrictComplex - - --- Roots of Unity --------------------------------------------------------------------------------- - --- | Fill a vector with roots of unity (Rofu) -calcRofu - :: Shape sh - => (sh :. Int) -- ^ Length of resulting vector. - -> Array (sh :. Int) Complex - -calcRofu sh@(_ :. n) - = force $ fromFunction sh f - where - f :: Shape sh => (sh :. Int) -> Complex - f (_ :. i) = (cos (2 * pi * (fromIntegral i) / len)) - :*: (- sin (2 * pi * (fromIntegral i) / len)) - - len = fromIntegral n - - --- | Fill a vector with roots of unity (Rofu) --- for the inverse transform. -calcInverseRofu - :: Shape sh - => (sh :. Int) -- ^ Length of resulting vector. - -> Array (sh :. Int) Complex - -calcInverseRofu sh@(_ :. n) - = force $ fromFunction sh f - where - f :: Shape sh => (sh :. Int) -> Complex - f (_ :. i) = (cos (2 * pi * (fromIntegral i) / len)) - :*: (sin (2 * pi * (fromIntegral i) / len)) - - len = fromIntegral n diff --git a/FFT/src/StrictComplex.hs b/FFT/src/StrictComplex.hs deleted file mode 100644 --- a/FFT/src/StrictComplex.hs +++ /dev/null @@ -1,32 +0,0 @@ -{-# LANGUAGE TypeOperators, TypeSynonymInstances #-} - -module StrictComplex - ( Complex(..) - , mag - , (:*:)(..)) -where -import Data.Array.Parallel.Base ((:*:)(..)) - --- | Strict complex doubles. -type Complex - = Double :*: Double - -instance Num Complex where - (r :*: i) + (r' :*: i') = r+r' :*: i+i' - (r :*: i) - (r' :*: i') = r-r' :*: i-i' - (r :*: i) * (r' :*: i') = r*r' - i*i' :*: r*i' + r'*i - fromInteger n = fromInteger n :*: 0.0 - -instance Fractional Complex where - (a :*: b) / (c :*: d) - = let den = c^2 + d^2 - re = (a * c + b * d) / den - im = (b * c - a * d) / den - in re :*: im - - --- | Take the magnitude of a complex number. -mag :: Complex -> Double -mag (r :*: i) = sqrt (r * r + i * i) - - diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,13 +1,24 @@ -Copyright (c) 2010 The DPH Team +Copyright (c) 2010, University of New South Wales. +All rights reserved. - Permission is hereby granted, free of charge, to any person - obtaining a copy of this software and associated documentation - files (the "Software"), to deal in the Software without - restriction, including without limitation the rights to use, - copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the - Software is furnished to do so, subject to the following - condition: +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the University of New South Wales nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Laplace/src/Main.hs b/Laplace/src/Main.hs --- a/Laplace/src/Main.hs +++ b/Laplace/src/Main.hs @@ -1,16 +1,16 @@ {-# LANGUAGE BangPatterns #-} -- | Solver for the Laplace equation --- Writes a PPM file of the completed solution. --- --- You can use the ImageMagick convert program to make a png --- with "convert out.ppm out.png" +-- 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 ColorRamp -import PPM +import Data.Array.Repa.IO.BMP +import Data.Array.Repa.IO.ColorRamp import System.Environment +import Data.Word +import Control.Monad main :: IO () main @@ -27,11 +27,11 @@ -- | Command line usage information. usage :: String usage = unlines - [ "Usage: laplace " + [ "Usage: laplace " , "" , " iterations :: Int Number of iterations to use in the solver." - , " input.ppm :: FileName ASCII 8 bit RGB PPM file for initial and boundary values." - , " output.ppm :: FileName PPM file to write output to." + , " 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," @@ -52,38 +52,95 @@ laplace steps fileInput fileOutput = do -- Load up the file containing boundary conditions. - (matBoundMask, matBoundValue) - <- readPPMAsMatrix2 loadPixel fileInput + 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 matInitial = matBoundValue + let arrInitial = arrBoundValue -- Run the solver. - let matFinal = solveLaplace + let arrFinal = solveLaplace steps - matBoundMask - matBoundValue - matInitial + arrBoundMask + arrBoundValue + arrInitial - matFinal `deepSeqArray` return () + arrFinal `deepSeqArray` return () - -- Write out the matrix as a colorised PPM image - writeMatrixAsNormalisedPPM + -- Make the result image + let arrImageOut = makeImageFromDoubles (rampColorHotToCold 0.0 1.0) arrFinal + + -- Write out the image to file. + writeImageToBMP fileOutput - (rampColorHotToCold 0.0 1.0) - matFinal + arrImageOut --- | Extract boundary mask and value from a pixel value. -loadPixel :: Int -> Int -> Int -> (Double, Double) -loadPixel r g b + +slurpDoublesFromImage + :: (Word8 -> Word8 -> Word8 -> Double) + -> Array DIM3 Word8 + -> Array DIM2 Double + +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 + +makeImageFromDoubles fnColor arrDoubles + = traverse arrDoubles + (\(Z :. height :. width) + -> Z :. height :. width :. 4) + + (\get (Z :. y :. x :. c) + -> let (r, g, b) = rampColorHotToCold 0 1 (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 +slurpBoundValue r g b -- A non-boundary value. | r == 0 && g == 0 && b == 255 - = (1, 0) + = 0 -- A boundary value. | (r == g) && (r == b) - = (0, (fromIntegral r) / 255) + = 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 +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/MMult/src/Main.hs b/MMult/src/Main.hs --- a/MMult/src/Main.hs +++ b/MMult/src/Main.hs @@ -1,11 +1,13 @@ {-# LANGUAGE PatternGuards #-} -import Solver -import Data.Array.Repa +import Data.Array.Repa as A +import Data.Array.Repa.IO.Matrix +import Data.Array.Repa.Algorithms.Matrix import Data.Maybe -import Matrix import System.Environment import Control.Monad +import System.Random +import qualified Data.Array.Parallel.Unlifted as U -- Arg Parsing ------------------------------------------------------------------------------------ data Arg @@ -64,6 +66,28 @@ 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 () @@ -85,18 +109,16 @@ `deepSeqArray` return () -- Run the solver. - let matResult = mmMult mat1 mat2 + let matResult = multiplyMM mat1 mat2 matResult `deepSeqArray` return () -- Write the output to file if requested. case mArgOut of Nothing -> return () - Just fileOut -> writeMatrixAsTextFile matResult fileOut + Just fileOut -> writeMatrixToTextFile fileOut matResult | otherwise = printHelp - - diff --git a/MMult/src/Solver.hs b/MMult/src/Solver.hs deleted file mode 100644 --- a/MMult/src/Solver.hs +++ /dev/null @@ -1,14 +0,0 @@ -{-# LANGUAGE TypeOperators #-} - -module Solver where -import Data.Array.Repa - -mmMult :: Array DIM2 Double - -> Array DIM2 Double - -> Array DIM2 Double - -mmMult arr1@Manifest{} arr2@Manifest{} - = fold (+) 0 - $ traverse2 arr1 (force $ transpose arr2) - (\(sh :. m1 :. n1) -> \(_ :. n2 :. m2) -> (sh :. m1 :. n2 :. n1)) - (\f1 -> \f2 -> \(sh :. i :. j :. k) -> f1 (sh :. i :. k) * f2 (sh :. j :. k)) diff --git a/lib/ColorRamp.hs b/lib/ColorRamp.hs deleted file mode 100644 --- a/lib/ColorRamp.hs +++ /dev/null @@ -1,47 +0,0 @@ - -module ColorRamp - (rampColorHotToCold) -where - - --- Color Ramps ----------------------------------------------------------------------------------- --- | Standard Hot -> Cold hypsometric color ramp. --- Sequence is red, yellow, green, cyan, blue. -rampColorHotToCold - :: (Ord a, Floating a) - => a - -> a - -> a - -> (a, a, a) - -{-# INLINE rampColorHotToCold #-} -rampColorHotToCold vmin vmax vNotNorm - = let - v | vNotNorm < vmin = vmin - | vNotNorm > vmax = vmax - | otherwise = vNotNorm - - dv = vmax - vmin - - result | v < vmin + 0.25 * dv - = ( 0 - , 4 * (v - vmin) / dv - , 1.0) - - | v < vmin + 0.5 * dv - = ( 0 - , 1.0 - , 1 + 4 * (vmin + 0.25 * dv - v) / dv) - - | v < vmin + 0.75 * dv - = ( 4 * (v - vmin - 0.5 * dv) / dv - , 1.0 - , 0.0) - - | otherwise - = ( 1.0 - , 1 + 4 * (vmin + 0.75 * dv - v) / dv - , 0) - - in result - diff --git a/lib/Matrix.hs b/lib/Matrix.hs deleted file mode 100644 --- a/lib/Matrix.hs +++ /dev/null @@ -1,86 +0,0 @@ - --- | Reading and writing matricies as ASCII files. --- We use ASCII so we can generate and check simple test data by hand, --- and we don't want to fool around with byte order issues. --- --- Matrix file format is like: --- --- MATRIX -- header --- 100 100 -- width and height --- 1.23 1.56 1.23 ... -- data, separated by whitespace --- .... --- --- TODO: Merge this with PPM.hs --- TODO: Merge this with Vector.hs. --- We should really have fns that read and write arrays of arbitrary dimension. --- -module Matrix - ( readMatrixFromTextFile - , writeMatrixAsTextFile - , genRandomMatrix) -where -import Data.List as L -import Data.Array.Repa as A -import qualified Data.Array.Parallel.Unlifted as U -import Prelude as P -import System.IO -import Control.Monad -import Data.Char -import System.Random -import Vector - - --- Reading ---------------------------------------------------------------------------------------- --- | Read a matrix from a text file. -readMatrixFromTextFile - :: (U.Elt a, Num a, Read a) - => FilePath -- ^ File name of matrix file. - -> IO (Array DIM2 a) - -readMatrixFromTextFile fileName - = do handle <- openFile fileName ReadMode - - "MATRIX" <- hGetLine handle - [width, height] <- liftM (P.map read . words) $ hGetLine handle - str <- hGetContents handle - let vals = readValues str - - let dims = Z :. width :. height - let mat = fromList dims vals - - return mat - - - --- | Write a matrix as a text file. -writeMatrixAsTextFile - :: (U.Elt a, Show a) - => Array DIM2 a -- ^ Matrix to write. - -> FilePath -- ^ File name of output file. - -> IO () - -writeMatrixAsTextFile arr fileName - = do file <- openFile fileName WriteMode - - hPutStrLn file "MATRIX" - - let Z :. width :. height - = extent arr - - hPutStrLn file $ show width ++ " " ++ show height - - hWriteValues file $ toList arr - hClose file - - --- | Generate a random(ish) matrix. -genRandomMatrix - :: DIM2 - -> IO (Array DIM2 Double) - -genRandomMatrix sh - = do uarr <- genRandomUArray (A.size sh) - return $ fromUArray sh uarr - - - diff --git a/lib/PPM.hs b/lib/PPM.hs deleted file mode 100644 --- a/lib/PPM.hs +++ /dev/null @@ -1,171 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables, BangPatterns #-} - --- | Writing out matricies as PPM image files. -module PPM - ( writeMatrixAsNormalisedPPM - , writeMatrixAsPPM - , readPPMAsMatrix - , readPPMAsMatrix2) -where -import qualified Data.Array.Parallel.Unlifted as U -import Data.List as L -import Data.Array.Repa as A -import Prelude as P -import System.IO -import Control.Monad -import Data.Char - - --- Write ------------------------------------------------------------------------------------------ --- | Convert a matrix to a PPM image, --- while normalising it to the maximum value present in the matrix. -writeMatrixAsNormalisedPPM - :: FilePath -- ^ Filename of output file. - -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values, - -> Array DIM2 Double -- ^ Matrix of values (need not be normalised). - -> IO () - -writeMatrixAsNormalisedPPM fileName colorFn arr - = let -- Use the maximum elem in the array as the white value. - vals = U.toList $ toUArray arr - maxVal = maximum vals - - -- Normalise the array to the range [0..1] for display. - arrNorm = A.map (/ maxVal) arr - - in writeMatrixAsPPM fileName colorFn arrNorm - - --- | Convert a matrix to a PPM image. --- Matrix elements should be normalised to [0..1] -writeMatrixAsPPM - :: FilePath - -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values. - -> Array DIM2 Double -- ^ Matrix of values, normalised to [0..1] - -> IO () - -writeMatrixAsPPM fileName colorFn arr - = let - -- Break flat array data into individual rows - Z :. width :. height - = extent arr - - in do - file <- openFile fileName WriteMode - hPutStrLn file $ "P3" - hPutStrLn file $ show width ++ " " ++ show height - hPutStrLn file $ "255" - - hWritePixels file colorFn $ toList arr - hClose file - - --- | Write out pixel values to a file. -hWritePixels - :: Handle - -> (Double -> (Double, Double, Double)) -- ^ Function for producing colors from data values. - -> [Double] -- ^ Data values. - -> IO () - -hWritePixels h colorFn !xx - = go xx - where - go [] = return () - go (x:xs) - = do let (r, g, b) = colorFn x - hPutStr h $ showInt $ truncate (r * 255) - hPutStr h $ " " - hPutStr h $ showInt $ truncate (g * 255) - hPutStr h $ " " - hPutStr h $ showInt $ truncate (b * 255) - hPutStr h $ "\n" - go xs - -showInt :: Int -> String -showInt i = show i - - --- Read ------------------------------------------------------------------------------------------- -readPPMAsMatrix - :: (Int -> Int -> Int -> Double) -- ^ Function for producing array values from RGB pixel values - -> FilePath -- ^ File name of ppm file. - -> IO (Array DIM2 Double) -- ^ Loaded matrix. - -readPPMAsMatrix pointFn fileName - = do file <- openFile fileName ReadMode - - "P3" <- hGetLine file - [width, height] <- liftM (P.map read . words) $ hGetLine file - _maxVal <- hGetLine file - vals <- loadPixels pointFn file - - let dims = Z :. width :. height - let mat = fromList dims vals - - return mat - - --- | Read the values of two separate matricies that are encoded in a single PPM file. -readPPMAsMatrix2 - :: (Int -> Int -> Int -> (Double, Double)) - -- ^ Function for producing array values from RGB pixel values - -> FilePath -- ^ File name of ppm file. - -> IO ( Array DIM2 Double -- ^ Loaded matrix. - , Array DIM2 Double) - -readPPMAsMatrix2 pointFn fileName - = do file <- openFile fileName ReadMode - - "P3" <- hGetLine file - [width, height] <- liftM (P.map read . words) $ hGetLine file - _maxVal <- hGetLine file - vals <- loadPixels pointFn file - - let dims = Z :. width :. height - let mat1 = fromList dims $ P.map fst vals - let mat2 = fromList dims $ P.map snd vals - - return (mat1, mat2) - - --- | Load list of array data from the data part of a PPM file. -loadPixels - :: (Int -> Int -> Int -> a) -- ^ Function for producing array values from RGB pixel values. - -> Handle -- ^ Handle of file to load from. - -> IO [a] - -loadPixels pointFn handle - = do str <- hGetContents handle - let ints = readInts str - let vals = convertLine pointFn ints - return vals - - --- | Convert groups of RGB components to values. -convertLine - :: (Int -> Int -> Int -> a) -- ^ Function for producing array values from RGB pixel values. - -> [Int] -- ^ Ints read from the PPM file. - -> [a] -- ^ Output values. - -convertLine fn vs - = case vs of - [] -> [] - r : g : b : rest -> fn r g b : convertLine fn rest - _ -> error "PPM.convertLine: bad pixel" - --- | Read a string containing ints separated by whitespace. -readInts :: String -> [Int] -readInts cs = readInts' [] cs - where readInts' _ [] = [] - readInts' acc (c : rest) - | isSpace c - = if null acc - then readInts' [] rest - else read (reverse acc) : readInts' [] rest - - | isDigit c - = readInts' (c : acc) rest - - | otherwise - = error $ "unexpected char in PPM file " ++ show (ord c) - diff --git a/lib/Vector.hs b/lib/Vector.hs deleted file mode 100644 --- a/lib/Vector.hs +++ /dev/null @@ -1,137 +0,0 @@ - --- | Reading and writing vectors as ASCII files. --- We use ASCII so we can generate and check simple test data by hand, --- and we don't want to fool around with byte order issues. --- --- Vector file format is like: --- --- VECTOR -- header --- 100 -- length of vector --- 1.23 1.56 1.23 ... -- data, separated by whitespace --- .... --- --- -module Vector - ( readVectorFromTextFile - , writeVectorAsTextFile - , genRandomVector - - , readValues - , hWriteValues - , genRandomUArray) -where -import Data.List as L -import Data.Array.Repa as A -import qualified Data.Array.Parallel.Unlifted as U -import Prelude as P -import System.IO -import Control.Monad -import Data.Char -import System.Random - - --- | Read a vector from a text file. -readVectorFromTextFile - :: (U.Elt a, Num a, Read a) - => FilePath -- ^ File name of vector file. - -> IO (Array DIM1 a) - -readVectorFromTextFile fileName - = do handle <- openFile fileName ReadMode - - "VECTOR" <- hGetLine handle - [len] <- liftM (P.map readInt . words) $ hGetLine handle - str <- hGetContents handle - let vals = readValues str - - let dims = Z :. len - let vec = fromList dims vals - - return vec - -readInt :: String -> Int -readInt str - | and $ P.map isDigit str - = read str - - --- | Write a vector as a text file. -writeVectorAsTextFile - :: (U.Elt a, Show a) - => Array DIM1 a -- ^ Vector to write. - -> FilePath -- ^ File name of output file. - -> IO () - -writeVectorAsTextFile arr fileName - = do file <- openFile fileName WriteMode - - hPutStrLn file "VECTOR" - - let Z :. len - = extent arr - - hPutStrLn file $ show len - hWriteValues file $ toList arr - hClose file - - --- | Generate a random(ish) vector. -genRandomVector - :: DIM2 - -> IO (Array DIM2 Double) - -genRandomVector sh - = do uarr <- genRandomUArray (A.size sh) - return $ fromUArray sh uarr - - - --- Stuff shared with Matrix module ------------------------------------------------------------- --- | Write out values to a file. -hWriteValues - :: Show a - => Handle - -> [a] -- ^ Data values. - -> IO () - -hWriteValues handle xx - = go xx - where - go [] = return () - go (x:xs) - = do hPutStr handle $ show x - hPutStr handle $ "\n" - go xs - - --- | Read a string containing ints separated by whitespace. -readValues :: (Num a, Read a) => String -> [a] -readValues cs = readValues' [] cs - where readValues' _ [] = [] - readValues' acc (c : rest) - | isSpace c - = if null acc - then readValues' [] rest - else read (reverse acc) : readValues' [] rest - - | isDigit c || c == '.' || c == 'e' || c == '-' - = readValues' (c : acc) rest - - | otherwise - = error $ "unexpected char in Matrix file " ++ show (ord c) - - --- | 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 - - - - diff --git a/repa-examples.cabal b/repa-examples.cabal --- a/repa-examples.cabal +++ b/repa-examples.cabal @@ -1,6 +1,6 @@ Name: repa-examples -Version: 1.0.0.0 -License: MIT +Version: 1.1.0.0 +License: BSD3 License-file: LICENSE Author: The DPH Team Maintainer: Ben Lippmeier @@ -16,44 +16,59 @@ Synopsis: Examples using the Repa array library. - -Executable repa-laplace +Executable repa-fft2d Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1, - dph-prim-par >= 0.4 && < 0.5, - repa >= 1.0 && < 1.1 + base == 4.*, + dph-prim-par == 0.4.*, + dph-base == 0.4.*, + repa == 1.1.*, + repa-algorithms == 1.1.*, + repa-io == 1.1.* - Main-is: Laplace/src/Main.hs - other-modules: Solver Matrix Vector ColorRamp PPM - hs-source-dirs: Laplace/src . lib + Main-is: FFT/src/FFT2d/Main.hs + hs-source-dirs: FFT/src . ghc-options: -Odph -threaded -Executable repa-mmult +Executable repa-fft2d-highpass Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1, - dph-prim-par >= 0.4 && < 0.5, - repa >= 1.0 && < 1.1 + base == 4.*, + dph-prim-par == 0.4.*, + dph-base == 0.4.*, + repa == 1.1.*, + repa-algorithms == 1.1.*, + repa-io == 1.1.* - Main-is: MMult/src/Main.hs - other-modules: Solver Matrix Vector ColorRamp PPM - hs-source-dirs: MMult/src . lib + Main-is: FFT/src/HighPass/Main.hs + hs-source-dirs: FFT/src . ghc-options: -Odph -threaded -Executable repa-fft +Executable repa-laplace Build-depends: - base >= 3 && < 5, - random >= 1.0 && < 1.1, - dph-prim-par >= 0.4 && < 0.5, - dph-base >= 0.4 && < 0.5, - repa >= 1.0 && < 1.1 + base == 4.*, + dph-prim-par == 0.4.*, + dph-base == 0.4.*, + repa == 1.1.*, + repa-io == 1.1.* - Main-is: FFT/src/Main.hs - other-modules: FFT DFT Roots StrictComplex Matrix Vector ColorRamp PPM - hs-source-dirs: FFT/src . lib + Main-is: Laplace/src/Main.hs + other-modules: Solver + hs-source-dirs: Laplace/src . ghc-options: -Odph -threaded + +Executable repa-mmult + 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.* + + Main-is: MMult/src/Main.hs + hs-source-dirs: MMult/src . + ghc-options: -Odph -threaded