-- The University of New Mexico's Haskell Image Processing Library
-- Copyright (C) 2013 Joseph Collard
--
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program. If not, see .
{-# LANGUAGE TypeFamilies, ViewPatterns, FlexibleContexts, FlexibleInstances #-}
{-# OPTIONS_GHC -O2 #-}
module Data.Image.Boxed(
BoxedImage,
-- * Gray Images
GrayImage, Gray, readImage,
grayToComplex, makeHotImage,
ref',
-- * Color Images
ColorImage, Color(..), readColorImage,
colorImageRed, colorImageGreen, colorImageBlue,
colorImageToRGB, rgbToColorImage,
colorImageHue, colorImageSaturation, colorImageIntensity,
colorImageToHSI, hsiToColorImage,
-- * Complex Images
ComplexImage, Complex,
CI.makeFilter,
fft, ifft,
realPart, imagPart,
magnitude, angle,
complex, complexImageToRectangular,
complexImageToPolar,
shrink,
-- * Binary Images
distanceTransform, label,
-- * Additional Modules
-- | Contains functionality for performing arithmetic operations on images with scalar values.
module Data.Image.Arithmetic,
-- | Contains functionality related to Binary Images
module Data.Image.Binary,
-- | Contains functionality for convolution of images
module Data.Image.Convolution,
-- | Contains basic functionality for Images
module Data.Image.Internal,
-- | Contains functionality for writing images and displaying with an external program
module Data.Image.IO) where
import Data.Image.Arithmetic
import Data.Image.Binary hiding (distanceTransform, label)
import qualified Data.Image.Binary as Bin
import qualified Data.Image.Complex as CI
import Data.Image.Convolution
import Data.Image.Internal
import Data.Image.IO
--base>=4
import Control.Applicative
import qualified Data.Complex as C
import Data.Maybe(fromJust)
import Data.Monoid
--bytestring-0.10.0.2
import qualified Data.ByteString.Char8 as B
--vector>=0.10.0.1
import qualified Data.Vector as V
type Vector = V.Vector
-- Error Messages
differentDimensionsError = error "The images must have the same dimensions."
-- BoxedImage
-- | BoxedImage is a concrete implementation of Image using a boxed internal structure. This allows for it to be installed nicely in Functor and Applicative.
data BoxedImage a = Image { rs :: Int,
cs :: Int,
pixels :: Vector a}
instance Image (BoxedImage a) where
type Pixel (BoxedImage a) = a
rows = rs
cols = cs
ref i r c = (pixels i) V.! (r * (cols i) + c)
makeImage rows cols f = Image rows cols (V.fromList px) where
px = [ f r c | r <- [0..rows-1], c <- [0..cols-1]]
pixelList = V.toList . pixels
imageOp = liftA2
instance Functor BoxedImage where
fmap f (Image rows cols pixels) = Image rows cols (fmap f pixels)
instance Applicative BoxedImage where
pure a = Image 1 1 (V.fromList [a])
(<*>) (Image rows cols partial) (Image rows' cols' toApply)
| rows /= rows' && cols /= cols' = error "Cannot apply images of unequal dimensions."
| otherwise = Image rows cols (V.fromList applied) where
indices = [ r*cols + c | r <- [0..rows-1], c <- [0..cols-1]]
applied = map func indices
func i = (partial V.! i) (toApply V.! i)
instance Show (BoxedImage a) where
show (Image rows cols _) = "< Image " ++ (show rows) ++ "x" ++ (show cols) ++ " >"
instance Num a => Num (BoxedImage a) where
(+) = liftA2 (+)
(-) = liftA2 (-)
(*) = liftA2 (*)
abs = fmap abs
signum = fmap signum
fromInteger i = pure $ fromInteger i
instance Fractional a => Fractional (BoxedImage a) where
(/) = liftA2 (/)
recip = fmap recip
fromRational i = pure $ fromRational i
instance Ord a => Ord (BoxedImage a) where
(<=) img0 img1
| (rows img0) /= (rows img1) = differentDimensionsError
| (cols img0) /= (cols img1) = differentDimensionsError
| otherwise = and . zipWith (<=) (pixelList img0) . pixelList $ img1
instance Eq a => Eq (BoxedImage a) where
(==) img0 img1
| (rows img0) /= (rows img1) = False
| (cols img0) /= (cols img1) = False
| otherwise = and . zipWith (==) (pixelList img0) $ pixelList img1
-- GrayImage
{-| A concrete instance of Image representing a gray scale image.
This instance is installed in DisplayFormat as a gray PGM.
-}
type GrayImage = BoxedImage Gray
type Gray = Double
instance DisplayFormat GrayImage where
format = toPGM
instance GrayPixel Gray where
type GrayVal Gray = Gray
toGray = id
instance RGBPixel Gray where
type ColorVal Gray = Gray
toRGB px = (px, px, px)
instance HSIPixel Gray where
toHSI = toHSI . RGB . toRGB
instance BinaryPixel Gray where
toBinary 0.0 = False
toBinary _ = True
on = 1.0
off = 0.0
instance CI.ComplexPixel Gray where
type Value Gray = Double
toComplex i = i C.:+ 0.0
instance Monoid Gray where
mempty = 0.0
mappend = (+)
instance MaxMin Gray where
maximal = maximum
minimal = minimum
-- ColorImage
{-| A concrete instance of Image that represents images with color values.
This instance is installed in DisplayFormat and can be written to
a color PPM -}
type ColorImage = BoxedImage Color
class HSIPixel px where
toHSI :: px -> (Double, Double, Double)
instance DisplayFormat ColorImage where
format = toPPM
-- | A color encoding scheme
data Color =
-- | Red, Green, Blue encoding
RGB (Double, Double, Double)
-- | Hue, Saturation, Intensity encoding
| HSI (Double, Double, Double) deriving (Show, Eq)
instance GrayPixel Color where
type GrayVal Color = Double
toGray (RGB (r, g, b)) = (r + g + b) / 3.0
toGray (toRGB -> (r, g, b)) = (r + g + b) / 3.0
instance RGBPixel Color where
type ColorVal Color = Double
toRGB (RGB px) = px
toRGB (HSI (h, s, i)) = (r, g, b) where
r = i + v1
g = i - (v1/2) + v2
b = i - (v1/2) - v2
v1 = const*s*(cos h)/3
v2 = const*s*(sin h)/2
const = 2.44948974278318
instance HSIPixel Color where
toHSI (RGB (r, g, b)) = (h, s, i) where
h = if (v1 /= 0.0) then atan2 v2 v1 else 0
s = sqrt( (v1*v1) + (v2*v2) )
i = (r+g+b)/3
v1 = (2.0*r-g-b) / const
v2 = (g - b) / const
const = 2.44948974278318
toHSI (HSI px) = px
--Requires the image to be scaled
--instance ComplexPixel Color where
-- toComplex = undefined
instance BinaryPixel Color where
toBinary (toRGB -> (r, g, b))
| r == 0 && g == 0 && b == 0 = False
| otherwise = True
on = RGB (1.0, 1.0, 1.0)
off = RGB (0.0, 0.0, 0.0)
instance Monoid Color where
mempty = RGB (0.0, 0.0, 0.0)
mappend (toRGB -> (a,b,c)) (toRGB -> (d,e,f)) = RGB (a+d,b+e,c+f)
instance MaxMin Color where
maximal = helper max mempty . map toRGB
minimal = helper min (RGB (10e10, 10e10, 10e10)) . map toRGB
helper :: (Double -> Double -> Double) -> Color -> [(Double, Double, Double)] -> Color
helper compare (RGB (r,g,b)) [] = let i = foldr1 compare [r, g, b] in RGB (i,i,i)
helper compare (RGB (r, g, b)) ((r', g', b'):xs) = helper compare acc' xs where
acc' = (RGB (compare r r', compare g g', compare b b'))
instance Num Color where
(+) = colorOp (+)
(-) = colorOp (-)
(*) = colorOp (*)
abs (toRGB -> (r, g, b)) = RGB (abs r, abs g, abs b)
signum (toRGB -> (r, g, b)) = RGB (signum r, signum g, signum b)
fromInteger (fromIntegral -> i) = RGB (i,i,i)
instance Fractional Color where
(/) = colorOp (/)
recip (toRGB -> (r,g,b)) = RGB (recip r, recip g, recip b)
fromRational _ = error "Could not create Color from Rational."
colorOp :: (Double -> Double -> Double) -> Color -> Color -> Color
colorOp op (toRGB -> (a, b, c)) (toRGB -> (d, e, f)) = RGB (op a d, op b e, op c f)
-- ComplexImage
{-| A concrete instance of Image representing pixels as complex values.
This instance can be written to file as a color PPM.
-}
type ComplexImage = BoxedImage Complex
type Complex = C.Complex Double
instance BinaryPixel Complex where
toBinary (0.0 C.:+ 0.0) = False
toBinary _ = True
on = (1.0 C.:+ 0.0)
off = (0.0 C.:+ 0.0)
instance DisplayFormat ComplexImage where
format (complexImageToColorImage -> rgb) = toPPM rgb
instance CI.ComplexPixel Complex where
type Value Complex = Double
toComplex = id
complexImageToColorImage :: ComplexImage -> ColorImage
complexImageToColorImage img = fmap rgb img where
scale = complexScale img
rgb comp = if radius < 1 then RGB (red', grn', blu') else RGB (red, grn, blu) where
[red, grn, blu] = map (+d') [r',g',b']
[red', grn', blu'] = map (flip (-) d') [r',g',b']
[x, y] = map (*scale) [C.realPart comp, C.imagPart comp]
radius = sqrt((x*x) + (y*y))
a = onedivsqrtsix*x
b = sqrttwodivtwo*y
d = 1.0/(1.0 + (radius*radius))
d' = 0.5 - radius*d
r' = 0.5 + (twodivsqrtsix * x * d)
b' = 0.5 - (d * (a - b))
g' = 0.5 - (d * (a + b))
complexScale :: ComplexImage -> Double
complexScale (CI.complexImageToRectangular -> (real, imag)) = 2.0/(maxv - minv) where
maxr = maximum . pixelList $ (real :: GrayImage)
maxi = maximum . pixelList $ imag
minr = minimum . pixelList $ real
mini = minimum . pixelList $ imag
maxv = max maxr maxi
minv = min minr mini
twodivsqrtsix = 0.81649658092772603273 --2.0/sqrt(6)
onedivsqrtsix = 0.40824829046386301636 --1.0/sqrt(6)
sqrttwodivtwo = 0.70710678118654752440 --sqrt(2)/2.0
getComponent to component img = fmap (component . to) img
getRGB = getComponent toRGB
{-| Given a ColorImage, returns a GrayImage representing the Red color component
>>>let red = colorImageRed cactii
-}
colorImageRed :: ColorImage -> GrayImage
colorImageRed = getRGB (\ (r, _, _) -> r)
{-| Given a ColorImage, returns a GrayImage representing the Green color component
>>>let green = colorImageGreen cactii
-}
colorImageGreen :: ColorImage -> GrayImage
colorImageGreen = getRGB (\ (_,g,_) -> g)
{-| Given a ColorImage, returns a GrayImage representing the Blue color component
>>>let blue = colorImageBlue cactii
-}
colorImageBlue :: ColorImage -> GrayImage
colorImageBlue = getRGB (\ (_,_,b) -> b)
{-| Given a ColorImage, returns a triple containing three GrayImages each
containing one of the color components (red, green, blue)
>>>leftToRight' . colorImageToRGB $ cactii
-}
colorImageToRGB :: ColorImage -> (GrayImage, GrayImage, GrayImage)
colorImageToRGB img = (colorImageRed img, colorImageGreen img, colorImageBlue img)
{-| Given a triple containing three GrayImages each containing one of the
color components (red, green, blue), returns a ColorImage
>>>rgbToColorImage (red,green,blue)
-}
rgbToColorImage :: (GrayImage, GrayImage, GrayImage) -> ColorImage
rgbToColorImage (red, green, blue) = createRGB <$> red <*> green <*> blue where
createRGB r g b = RGB (r, g, b)
getHSI = getComponent toHSI
{-| Given a ColorImage, returns a GrayImage representing the Hue component
>>>let h = colorImageHue cactii
-}
colorImageHue :: ColorImage -> GrayImage
colorImageHue = getHSI (\ (h, _, _) -> h)
{-| Given a ColorImage, returns a GrayImage representing the Saturation component
>>>let s = colorImageSaturation cactii
-}
colorImageSaturation :: ColorImage -> GrayImage
colorImageSaturation = getHSI (\ (_,s,_) -> s)
{-| Given a ColorImage, returns a GrayImage representing the Intensity component
>>>let i = colorImageIntensity cactii
-}
colorImageIntensity :: ColorImage -> GrayImage
colorImageIntensity = getHSI (\ (_,_,i) -> i)
{-| Given a triple containing three GrayImages each containing one of the
color components (hue, saturation, ), returns a ColorImage
>>> hsiToColorImage (h, s, i)
-}
hsiToColorImage :: (GrayImage, GrayImage, GrayImage) -> ColorImage
hsiToColorImage (h, s, i) = toHSI <$> h <*> s <*> i where
toHSI h s i = HSI (h, s, i)
{-| Given a ColorImage, returns a triple containing three GrayImages each
containing one of the components (hue, saturation, intensity)
>>>let (h, s, i) = colorImageToHSI $ cactii
-}
colorImageToHSI :: ColorImage -> (GrayImage, GrayImage, GrayImage)
colorImageToHSI img = (colorImageHue img, colorImageSaturation img, colorImageIntensity img)
{-| Reads in an ASCI PPM file as a ColorImage
>>>cactii <- readColorImage "images/cactii.ppm"
-}
readColorImage :: FilePath -> IO ColorImage
readColorImage fileName =
do
y <- B.readFile fileName
return $ parseRGBPixelImage . B.intercalate (B.pack " ") . stripComments . B.lines $ y
parseRGBPixelImage :: B.ByteString -> ColorImage
parseRGBPixelImage string = Image rows cols (V.fromList rgbs)
where ws = B.words string
getInt = fst. fromJust . B.readInt
px = map (fromIntegral . getInt) $ drop 4 ws
cols = getInt $ ws !! 1
rows = getInt $ ws !! 2
maxi = fromIntegral . getInt $ ws !! 3
[r, g, b] = colors px
rgbs = map rgb3 . zip3 r g $ b
rgb3 (r, g, b) = RGB (r, g, b)
colors :: [Int] -> [[Gray]]
colors xs = helper xs [] [] []
where helper [] red green blue = map (map fromIntegral) $ map reverse [red, green, blue]
helper (r:g:b:cs) red green blue = helper cs (r:red) (g:green) (b:blue)
{-| Coerces a GrayImage to a ComplexImage where the imaginary
part for all pixels is 0.
>>>grayToComplex frog
-}
grayToComplex :: GrayImage -> ComplexImage
grayToComplex img = fmap (C.:+ 0.0) img
{-| Given a GrayImage, makeHotImage returns a ColorImage with the same
dimensions. The R, G, B values of the result image at (i, j) are
determined by using the value of the ColorImage at (i, j) to index
three lookup tables. These lookup tables implement a false coloring
scheme which maps small values to black, large values to white, and
intermediate values to shades of red, orange, and yellow (in that order).
>>>makeHotImage frog
-}
makeHotImage :: GrayImage -> ColorImage
makeHotImage img = fmap (toHot max min) img where
max = maxIntensity img
min = minIntensity img
toHot max min pixel = RGB (r, g, b) where
px = (pixel - min)/(max-min)
r = if px < 0.333333333 then (px*3.0) else 1.0
g = if px < 0.333333333 then 0.0 else
if px < 0.666666667 then (px - 0.333333333)*3 else 1.0
b = if px < 0.666666667 then 0.0 else (px - 0.666666667)*3
{-| Performs bilinear interpolation of a GrayImage at the coordinates provided. -}
ref' :: GrayImage -> Double -> Double -> Double
ref' im x y = if inside im x y then interpolate im x y else 0
inside im x y = x >= 0 && y >= 0 && y < r - 1 && x < c - 1
where r = fromIntegral $ rows im
c = fromIntegral $ cols im
interpolate :: GrayImage -> Double -> Double -> Double
interpolate im x y = (f01 - f00)*x' + (f10 - f00)*y' + (f11 + f00 - f10 - f01)*x'*y' + f00
where x' = x - (fromIntegral i0);
y' = y - (fromIntegral j0);
f00 = ref im i0 j0
f01 = ref im i0 j1
f10 = ref im i1 j0
f11 = ref im i1 j1
i1 = i0 + 1
j1 = j0 + 1
i0 = floor x
j0 = floor y
{-| Given a complex image, returns a real image representing
the real part of the image.
@
harmonicSignal :: Double -> Double -> Int -> Int -> C.Complex Double
harmonicSignal u v m n = exp (-pii*2.0 * var) where
pii = 0.0 C.:+ pi
var = (u*m' + v*n') C.:+ 0.0
[m',n'] = map fromIntegral [m, n]
@
>>> let signal = makeImage 128 128 (harmonicSignal (3/128) (2/128)) :: ComplexImage
>>>let cosine = realPart signal
>>>realPart realPart . ifft $ (fft frogpart) * (fft d2g)
>>>realPart realPart . ifft $ (fft frogpart) * (fft g)
-}
realPart :: ComplexImage -> GrayImage
realPart = CI.realPart
{-| Given a complex image, returns a real image representing
the imaginary part of the image
>>>let sine = imagPart signal
-}
imagPart :: ComplexImage -> GrayImage
imagPart = CI.imagPart
{-| Given a complex image, returns a real image representing
the magnitude of the image.
>>>magnitude signal
-}
magnitude :: ComplexImage -> GrayImage
magnitude = CI.magnitude
{-| Given a complex image, returns a real image representing
the angle of the image
>>>angle signal
-}
angle :: ComplexImage -> GrayImage
angle = CI.angle
{-| Given a complex image, returns a pair of real images each
representing the component (magnitude, phase) of the image
>>>leftToRight' . complexImageToPolar $ signal
-}
complexImageToPolar :: ComplexImage -> (GrayImage, GrayImage)
complexImageToPolar = CI.complexImageToPolar
{-| Given an image representing the real part of a complex image, and
an image representing the imaginary part of a complex image, returns
a complex image.
>>>complex cosine sine
-}
complex :: GrayImage -> GrayImage -> ComplexImage
complex = CI.complex
{-| Given a complex image, return a pair of real images each representing
a component of the complex image (real, imaginary).
>>>leftToRight' . complexImageToRectangular $ signal
-}
complexImageToRectangular :: ComplexImage -> (GrayImage, GrayImage)
complexImageToRectangular = CI.complexImageToRectangular
{-| Given a complex image and a real positive number x, shrink returns
a complex image with the same dimensions. Let z be the value of the
image at location (i, j). The value of the complex result image at
location (i, j) is zero if |z| < x, otherwise the result has the
same phase as z but the amplitude is decreased by x.
-}
shrink :: (Num a,
Image img,
CI.ComplexPixel (Pixel img),
CI.Value (Pixel img) ~ Double) => a -> img -> ComplexImage
shrink = CI.shrink
{-| Given an image whose pixels can be converted to a complex value,
fft returns an image with complex pixels representing its 2D discrete
Fourier transform (DFT). Because the DFT is computed using the Fast Fourier
Transform (FFT) algorithm, the number of rows and columns of the image
must both be powers of two, i.e., 2K where K is an integer.
>>>frog <- readImage "images/frog.pgm"
>>>let frogpart = crop 64 64 128 128 frog
>>>imageMap log . fft $ frogpart :: ComplexImage
>>>fft d2g
>>>fft g
-}
fft :: (Image img,
CI.ComplexPixel (Pixel img),
CI.Value (Pixel img) ~ Double) => img -> ComplexImage
fft = CI.fft
{-| Given an image, ifft returns a complex image representing its 2D
inverse discrete Fourier transform (DFT). Because the inverse DFT is
computed using the Fast Fourier Transform (FFT) algorithm, the number
of rows and columns of must both be powers of two, i.e., 2K
where K is an integer.
>>>ifft ((fft frogpart) * (fft d2g))
>>>ifft ((fft frogpart) * (fft g))
-}
ifft :: (Image img,
CI.ComplexPixel (Pixel img),
CI.Value (Pixel img) ~ Double) => img -> ComplexImage
ifft = CI.ifft
-- Binary Images
{-| Given a binary image, distanceTransform returns an image
representing the 2D distance transform of the image.
The distance transform is accurate to within a 2% error for euclidean
distance.
>>>distanceTransform binaryStop :: GrayImage
< Image 86x159 >
-}
distanceTransform :: (Image img,
BinaryPixel (Pixel img)) => img -> GrayImage
distanceTransform = Bin.distanceTransform
{-| Given a binary image, label returns an image where pixels in
distinct connected components (based on 4-neighbor connectivity)
have distinct integer values. These values range from 1 to n where
n is the number of connected components in image.
>>> label binaryStop
< Image 86x159 >
-}
label :: (Image img,
BinaryPixel (Pixel img)) => img -> GrayImage
label = Bin.label
{-| Reads in a ASCII PGM image located at fileName as a GrayImage
>>>frog <- readImage "images/frog.pgm"
-}
readImage :: FilePath -> IO GrayImage
readImage fileName =
do
y <- B.readFile fileName
return $ parseImage . B.intercalate (B.pack " ") . stripComments . B.lines $ y
parseImage :: B.ByteString -> GrayImage
parseImage string = img
where ws = B.words string
getInt = fst . fromJust . B.readInt
px = map (fromIntegral . getInt) $ drop 4 ws
cols = getInt $ ws !! 1
rows = getInt $ ws !! 2
maxi = fromIntegral . getInt $ ws !! 3
img = Image rows cols (V.fromList px)
stripComments :: [B.ByteString] -> [B.ByteString]
stripComments xs = filter pred xs
where pred x
| B.null x = False
| B.head x == '#' = False
| otherwise = True