diff --git a/bench/Benchmark.hs b/bench/Benchmark.hs --- a/bench/Benchmark.hs +++ b/bench/Benchmark.hs @@ -124,6 +124,7 @@ , bgroup "threshold" [ bench "simple threshold" $ whnf threshold' grey , bench "adaptive threshold" $ whnf adaptiveThreshold' grey + , bench "Otsu's method" $ whnf otsu' grey ] @@ -175,6 +176,9 @@ filt = I.adaptiveThreshold (I.GaussianKernel Nothing) 1 0 (I.BinaryThreshold 0 255) in filt `I.apply` img + + otsu' :: Grey -> Grey + otsu' !img = I.otsu (I.BinaryThreshold 0 255) img miniature !rgb = let Z :. h :. w = I.shape rgb diff --git a/example/Threshold.hs b/example/Threshold.hs new file mode 100644 --- /dev/null +++ b/example/Threshold.hs @@ -0,0 +1,18 @@ +import System.Environment (getArgs) + +import Vision.Image + +-- Thresholds an image by applying the Otsu's method. +-- +-- usage: ./threshold .png output.png +main :: IO () +main = do + [input, output] <- getArgs + io <- either (\x -> error $ "Load failed: " ++ show x) return + =<< load Nothing input + + let grey = convert io :: Grey + thresholded = otsu (BinaryThreshold 0 255) grey :: Grey + + _ <- save output thresholded + return () diff --git a/friday.cabal b/friday.cabal --- a/friday.cabal +++ b/friday.cabal @@ -1,5 +1,5 @@ name: friday -version: 0.1.4 +version: 0.1.5 synopsis: A functionnal image processing library for Haskell. homepage: https://github.com/RaphaelJ/friday license: LGPL-3 @@ -50,18 +50,22 @@ Vision.Image Vision.Image.Grey Vision.Image.Grey.Conversion + Vision.Image.Grey.Specialize Vision.Image.Grey.Type Vision.Image.Filter Vision.Image.HSV Vision.Image.HSV.Conversion + Vision.Image.HSV.Specialize Vision.Image.HSV.Type Vision.Image.Interpolate Vision.Image.Mutable Vision.Image.RGBA Vision.Image.RGBA.Conversion + Vision.Image.RGBA.Specialize Vision.Image.RGBA.Type Vision.Image.RGB Vision.Image.RGB.Conversion + Vision.Image.RGB.Specialize Vision.Image.RGB.Type Vision.Image.Storage Vision.Image.Threshold @@ -146,6 +150,17 @@ build-depends: base >= 4 && < 5 , friday +executable threshold + if !flag(examples) + Buildable: False + + main-is: Threshold.hs + ghc-options: -Wall -O2 -rtsopts + hs-source-dirs: example/ + default-language: Haskell2010 + + build-depends: base >= 4 && < 5 + , friday Benchmark benchmark type: exitcode-stdio-1.0 diff --git a/src/Vision/Histogram.hs b/src/Vision/Histogram.hs --- a/src/Vision/Histogram.hs +++ b/src/Vision/Histogram.hs @@ -26,15 +26,14 @@ import Foreign.Storable (Storable) import Prelude hiding (map) -import Vision.Image ( - Pixel, MaskedImage, Image, ImagePixel, FunctorImage - , Grey, GreyPixel (..), HSV, HSVPixel (..), RGBA, RGBAPixel (..) - , RGB, RGBPixel (..) - ) -import qualified Vision.Image as I +import Vision.Image.Grey.Type (GreyPixel (..)) +import Vision.Image.HSV.Type (HSVPixel (..)) +import Vision.Image.RGBA.Type (RGBAPixel (..)) +import Vision.Image.RGB.Type (RGBPixel (..)) +import Vision.Image.Type (Pixel, MaskedImage, Image, ImagePixel, FunctorImage) +import qualified Vision.Image.Type as I import Vision.Primitive ( - Z (..), (:.) (..), Shape (..), DIM1, DIM3, DIM4, DIM5, DIM6 - , ix1, ix3, ix4 + Z (..), (:.) (..), Shape (..), DIM1, DIM3, DIM4, DIM5, ix1, ix3, ix4 ) -- There is no rule to simplify the conversion from Int32 to Double and Float @@ -161,8 +160,8 @@ -- -- If the size of the histogram is specified, every bin of a given dimension -- will be of the same size (uniform histogram). -histogram :: (MaskedImage i, ToHistogram (ImagePixel i), Storable a, Num a - , HistogramShape (PixelValueSpace (ImagePixel i))) +histogram :: ( MaskedImage i, ToHistogram (ImagePixel i), Storable a, Num a + , HistogramShape (PixelValueSpace (ImagePixel i))) => Maybe (PixelValueSpace (ImagePixel i)) -> i -> Histogram (PixelValueSpace (ImagePixel i)) a histogram mSize img = @@ -181,18 +180,6 @@ case mSize of Just _ -> pixToBin size p Nothing -> pixToIndex p {-# INLINE toIndex #-} -{-# SPECIALIZE histogram :: Maybe DIM1 -> Grey -> Histogram DIM1 Int32 - , Maybe DIM1 -> Grey -> Histogram DIM1 Double - , Maybe DIM1 -> Grey -> Histogram DIM1 Float - , Maybe DIM3 -> HSV -> Histogram DIM3 Int32 - , Maybe DIM3 -> HSV -> Histogram DIM3 Double - , Maybe DIM3 -> HSV -> Histogram DIM3 Float - , Maybe DIM4 -> RGBA -> Histogram DIM4 Int32 - , Maybe DIM4 -> RGBA -> Histogram DIM4 Double - , Maybe DIM4 -> RGBA -> Histogram DIM4 Float - , Maybe DIM3 -> RGB -> Histogram DIM3 Int32 - , Maybe DIM3 -> RGB -> Histogram DIM3 Double - , Maybe DIM3 -> RGB -> Histogram DIM3 Float #-} {-# INLINABLE histogram #-} -- | Similar to 'histogram' but adds two dimensions for the y and x-coordinates @@ -205,8 +192,8 @@ -- -- As there is no reason to create an histogram as large as the number of pixels -- of the image, a size is always needed. -histogram2D :: (Image i, ToHistogram (ImagePixel i), Storable a, Num a - , HistogramShape (PixelValueSpace (ImagePixel i))) +histogram2D :: ( Image i, ToHistogram (ImagePixel i), Storable a, Num a + , HistogramShape (PixelValueSpace (ImagePixel i))) => (PixelValueSpace (ImagePixel i)) :. Int :. Int -> i -> Histogram ((PixelValueSpace (ImagePixel i)) :. Int :. Int) a histogram2D size img = @@ -226,18 +213,6 @@ let !ix = (pixToIndex p) :. y :. x in toLinearIndex size $! toBin size maxSize ix {-# INLINE toIndex #-} -{-# SPECIALIZE histogram2D :: DIM3 -> Grey -> Histogram DIM3 Int32 - , DIM3 -> Grey -> Histogram DIM3 Double - , DIM3 -> Grey -> Histogram DIM3 Float - , DIM5 -> HSV -> Histogram DIM5 Int32 - , DIM5 -> HSV -> Histogram DIM5 Double - , DIM5 -> HSV -> Histogram DIM5 Float - , DIM6 -> RGBA -> Histogram DIM6 Int32 - , DIM6 -> RGBA -> Histogram DIM6 Double - , DIM6 -> RGBA -> Histogram DIM6 Float - , DIM5 -> RGB -> Histogram DIM5 Int32 - , DIM5 -> RGB -> Histogram DIM5 Double - , DIM5 -> RGB -> Histogram DIM5 Float #-} {-# INLINABLE histogram2D #-} -- Reshaping ------------------------------------------------------------------- @@ -320,9 +295,10 @@ -- @H@ the cumulative of the histogram normalized over an @L1@ norm. -- -- See . -equalizeImage :: (FunctorImage i i, Integral (ImagePixel i) +equalizeImage :: ( FunctorImage i i, Integral (ImagePixel i) , ToHistogram (ImagePixel i) - , PixelValueSpace (ImagePixel i) ~ DIM1) => i -> i + , PixelValueSpace (ImagePixel i) ~ DIM1) + => i -> i equalizeImage img = I.map equalizePixel img where @@ -332,8 +308,6 @@ !cumNormalized' = map round cumNormalized :: Histogram DIM1 Int32 equalizePixel !val = fromIntegral $ cumNormalized' `index` ix1 (int val) {-# INLINE equalizePixel #-} --- FIXME: GHC 7.8.2 fails to specialize -{-# SPECIALIZE equalizeImage :: Grey -> Grey #-} {-# INLINABLE equalizeImage #-} -- Comparisons ----------------------------------------------------------------- diff --git a/src/Vision/Image/Grey.hs b/src/Vision/Image/Grey.hs --- a/src/Vision/Image/Grey.hs +++ b/src/Vision/Image/Grey.hs @@ -2,5 +2,6 @@ module Vision.Image.Grey.Type ) where -import Vision.Image.Grey.Type import Vision.Image.Grey.Conversion () +import Vision.Image.Grey.Specialize () +import Vision.Image.Grey.Type diff --git a/src/Vision/Image/Grey/Specialize.hs b/src/Vision/Image/Grey/Specialize.hs new file mode 100644 --- /dev/null +++ b/src/Vision/Image/Grey/Specialize.hs @@ -0,0 +1,33 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} + +-- | @SPECIALIZE@ pragma declarations for grey-scale images. +module Vision.Image.Grey.Specialize () where + +import Data.Int + +import Vision.Histogram (Histogram, histogram, histogram2D, equalizeImage) +import Vision.Image.Grey.Type (Grey, GreyPixel) +import Vision.Image.Threshold (ThresholdType, otsu) +import Vision.Image.Transform ( + InterpolMethod, crop, resize, horizontalFlip, verticalFlip + ) +import Vision.Primitive (DIM1, DIM3, Rect, Size) + +{-# SPECIALIZE histogram :: Maybe DIM1 -> Grey -> Histogram DIM1 Int32 + , Maybe DIM1 -> Grey -> Histogram DIM1 Double + , Maybe DIM1 -> Grey -> Histogram DIM1 Float #-} + +{-# SPECIALIZE histogram2D :: DIM3 -> Grey -> Histogram DIM3 Int32 + , DIM3 -> Grey -> Histogram DIM3 Double + , DIM3 -> Grey -> Histogram DIM3 Float #-} + +-- FIXME: GHC 7.8.2 fails to specialize +{-# SPECIALIZE equalizeImage :: Grey -> Grey #-} + +{-# SPECIALIZE crop :: Rect -> Grey -> Grey #-} +{-# SPECIALIZE resize :: InterpolMethod -> Size -> Grey -> Grey #-} +{-# SPECIALIZE horizontalFlip :: Grey -> Grey #-} +{-# SPECIALIZE verticalFlip :: Grey -> Grey #-} + +{-# SPECIALIZE otsu :: ThresholdType GreyPixel GreyPixel -> Grey + -> Grey #-} diff --git a/src/Vision/Image/Grey/Type.hs b/src/Vision/Image/Grey/Type.hs --- a/src/Vision/Image/Grey/Type.hs +++ b/src/Vision/Image/Grey/Type.hs @@ -9,11 +9,7 @@ import Foreign.Storable (Storable) import Vision.Image.Interpolate (Interpolable (..)) -import Vision.Image.Transform ( - InterpolMethod, crop, resize, horizontalFlip, verticalFlip - ) import Vision.Image.Type (Pixel (..), Manifest, Delayed) -import Vision.Primitive (Rect, Size) newtype GreyPixel = GreyPixel Word8 deriving (Bits, Bounded, Enum, Eq, Integral, Num, Ord, Real, Read, Show @@ -35,8 +31,3 @@ instance Interpolable GreyPixel where interpol f (GreyPixel a) (GreyPixel b) = GreyPixel $ f a b {-# INLINE interpol #-} - -{-# SPECIALIZE crop :: Rect -> Grey -> Grey #-} -{-# SPECIALIZE resize :: InterpolMethod -> Size -> Grey -> Grey #-} -{-# SPECIALIZE horizontalFlip :: Grey -> Grey #-} -{-# SPECIALIZE verticalFlip :: Grey -> Grey #-} diff --git a/src/Vision/Image/HSV.hs b/src/Vision/Image/HSV.hs --- a/src/Vision/Image/HSV.hs +++ b/src/Vision/Image/HSV.hs @@ -2,5 +2,6 @@ module Vision.Image.HSV.Type ) where -import Vision.Image.HSV.Type import Vision.Image.HSV.Conversion () +import Vision.Image.HSV.Specialize () +import Vision.Image.HSV.Type diff --git a/src/Vision/Image/HSV/Specialize.hs b/src/Vision/Image/HSV/Specialize.hs new file mode 100644 --- /dev/null +++ b/src/Vision/Image/HSV/Specialize.hs @@ -0,0 +1,26 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} + +-- | @SPECIALIZE@ pragma declarations for HSV images. +module Vision.Image.HSV.Specialize () where + +import Data.Int + +import Vision.Histogram (Histogram, histogram, histogram2D) +import Vision.Image.HSV.Type (HSV) +import Vision.Image.Transform ( + InterpolMethod, crop, resize, horizontalFlip, verticalFlip + ) +import Vision.Primitive (DIM3, DIM5, Rect, Size) + +{-# SPECIALIZE histogram :: Maybe DIM3 -> HSV -> Histogram DIM3 Int32 + , Maybe DIM3 -> HSV -> Histogram DIM3 Double + , Maybe DIM3 -> HSV -> Histogram DIM3 Float #-} + +{-# SPECIALIZE histogram2D :: DIM5 -> HSV -> Histogram DIM5 Int32 + , DIM5 -> HSV -> Histogram DIM5 Double + , DIM5 -> HSV -> Histogram DIM5 Float #-} + +{-# SPECIALIZE crop :: Rect -> HSV -> HSV #-} +{-# SPECIALIZE resize :: InterpolMethod -> Size -> HSV -> HSV #-} +{-# SPECIALIZE horizontalFlip :: HSV -> HSV #-} +{-# SPECIALIZE verticalFlip :: HSV -> HSV #-} diff --git a/src/Vision/Image/HSV/Type.hs b/src/Vision/Image/HSV/Type.hs --- a/src/Vision/Image/HSV/Type.hs +++ b/src/Vision/Image/HSV/Type.hs @@ -10,11 +10,7 @@ import Foreign.Ptr (castPtr, plusPtr) import Vision.Image.Interpolate (Interpolable (..)) -import Vision.Image.Transform ( - InterpolMethod, crop, resize, horizontalFlip, verticalFlip - ) import Vision.Image.Type (Pixel (..), Manifest, Delayed) -import Vision.Primitive (Rect, Size) data HSVPixel = HSVPixel { hsvHue :: {-# UNPACK #-} !Word8, hsvSat :: {-# UNPACK #-} !Word8 @@ -74,8 +70,3 @@ , hsvValue = f aVal bVal } {-# INLINE interpol #-} - -{-# SPECIALIZE crop :: Rect -> HSV -> HSV #-} -{-# SPECIALIZE resize :: InterpolMethod -> Size -> HSV -> HSV #-} -{-# SPECIALIZE horizontalFlip :: HSV -> HSV #-} -{-# SPECIALIZE verticalFlip :: HSV -> HSV #-} diff --git a/src/Vision/Image/RGB.hs b/src/Vision/Image/RGB.hs --- a/src/Vision/Image/RGB.hs +++ b/src/Vision/Image/RGB.hs @@ -2,5 +2,6 @@ module Vision.Image.RGB.Type ) where -import Vision.Image.RGB.Type import Vision.Image.RGB.Conversion () +import Vision.Image.RGB.Specialize () +import Vision.Image.RGB.Type diff --git a/src/Vision/Image/RGB/Specialize.hs b/src/Vision/Image/RGB/Specialize.hs new file mode 100644 --- /dev/null +++ b/src/Vision/Image/RGB/Specialize.hs @@ -0,0 +1,26 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} + +-- | @SPECIALIZE@ pragma declarations for RGB images. +module Vision.Image.RGB.Specialize () where + +import Data.Int + +import Vision.Histogram (Histogram, histogram, histogram2D) +import Vision.Image.RGB.Type (RGB) +import Vision.Image.Transform ( + InterpolMethod, crop, resize, horizontalFlip, verticalFlip + ) +import Vision.Primitive (DIM3, DIM5, Rect, Size) + +{-# SPECIALIZE histogram :: Maybe DIM3 -> RGB -> Histogram DIM3 Int32 + , Maybe DIM3 -> RGB -> Histogram DIM3 Double + , Maybe DIM3 -> RGB -> Histogram DIM3 Float #-} + +{-# SPECIALIZE histogram2D :: DIM5 -> RGB -> Histogram DIM5 Int32 + , DIM5 -> RGB -> Histogram DIM5 Double + , DIM5 -> RGB -> Histogram DIM5 Float #-} + +{-# SPECIALIZE crop :: Rect -> RGB -> RGB #-} +{-# SPECIALIZE resize :: InterpolMethod -> Size -> RGB -> RGB #-} +{-# SPECIALIZE horizontalFlip :: RGB -> RGB #-} +{-# SPECIALIZE verticalFlip :: RGB -> RGB #-} diff --git a/src/Vision/Image/RGB/Type.hs b/src/Vision/Image/RGB/Type.hs --- a/src/Vision/Image/RGB/Type.hs +++ b/src/Vision/Image/RGB/Type.hs @@ -10,11 +10,7 @@ import Foreign.Ptr (castPtr, plusPtr) import Vision.Image.Interpolate (Interpolable (..)) -import Vision.Image.Transform ( - InterpolMethod, crop, resize, horizontalFlip, verticalFlip - ) import Vision.Image.Type (Pixel (..), Manifest, Delayed) -import Vision.Primitive (Rect, Size) data RGBPixel = RGBPixel { rgbRed :: {-# UNPACK #-} !Word8, rgbGreen :: {-# UNPACK #-} !Word8 @@ -65,8 +61,3 @@ , rgbBlue = f aBlue bBlue } {-# INLINE interpol #-} - -{-# SPECIALIZE crop :: Rect -> RGB -> RGB #-} -{-# SPECIALIZE resize :: InterpolMethod -> Size -> RGB -> RGB #-} -{-# SPECIALIZE horizontalFlip :: RGB -> RGB #-} -{-# SPECIALIZE verticalFlip :: RGB -> RGB #-} diff --git a/src/Vision/Image/RGBA.hs b/src/Vision/Image/RGBA.hs --- a/src/Vision/Image/RGBA.hs +++ b/src/Vision/Image/RGBA.hs @@ -2,5 +2,6 @@ module Vision.Image.RGBA.Type ) where -import Vision.Image.RGBA.Type import Vision.Image.RGBA.Conversion () +import Vision.Image.RGBA.Specialize () +import Vision.Image.RGBA.Type diff --git a/src/Vision/Image/RGBA/Specialize.hs b/src/Vision/Image/RGBA/Specialize.hs new file mode 100644 --- /dev/null +++ b/src/Vision/Image/RGBA/Specialize.hs @@ -0,0 +1,26 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} + +-- | @SPECIALIZE@ pragma declarations for RGBA images. +module Vision.Image.RGBA.Specialize () where + +import Data.Int + +import Vision.Histogram (Histogram, histogram, histogram2D) +import Vision.Image.RGBA.Type (RGBA) +import Vision.Image.Transform ( + InterpolMethod, crop, resize, horizontalFlip, verticalFlip + ) +import Vision.Primitive (DIM4, DIM6, Rect, Size) + +{-# SPECIALIZE histogram :: Maybe DIM4 -> RGBA -> Histogram DIM4 Int32 + , Maybe DIM4 -> RGBA -> Histogram DIM4 Double + , Maybe DIM4 -> RGBA -> Histogram DIM4 Float #-} + +{-# SPECIALIZE histogram2D :: DIM6 -> RGBA -> Histogram DIM6 Int32 + , DIM6 -> RGBA -> Histogram DIM6 Double + , DIM6 -> RGBA -> Histogram DIM6 Float #-} + +{-# SPECIALIZE crop :: Rect -> RGBA -> RGBA #-} +{-# SPECIALIZE resize :: InterpolMethod -> Size -> RGBA -> RGBA #-} +{-# SPECIALIZE horizontalFlip :: RGBA -> RGBA #-} +{-# SPECIALIZE verticalFlip :: RGBA -> RGBA #-} diff --git a/src/Vision/Image/RGBA/Type.hs b/src/Vision/Image/RGBA/Type.hs --- a/src/Vision/Image/RGBA/Type.hs +++ b/src/Vision/Image/RGBA/Type.hs @@ -10,11 +10,7 @@ import Foreign.Ptr (castPtr, plusPtr) import Vision.Image.Interpolate (Interpolable (..)) -import Vision.Image.Transform ( - InterpolMethod, crop, resize, horizontalFlip, verticalFlip - ) import Vision.Image.Type (Pixel (..), Manifest, Delayed) -import Vision.Primitive (Rect, Size) data RGBAPixel = RGBAPixel { rgbaRed :: {-# UNPACK #-} !Word8, rgbaGreen :: {-# UNPACK #-} !Word8 @@ -67,8 +63,3 @@ , rgbaBlue = f aBlue bBlue, rgbaAlpha = f aAlpha bAlpha } {-# INLINE interpol #-} - -{-# SPECIALIZE crop :: Rect -> RGBA -> RGBA #-} -{-# SPECIALIZE resize :: InterpolMethod -> Size -> RGBA -> RGBA #-} -{-# SPECIALIZE horizontalFlip :: RGBA -> RGBA #-} -{-# SPECIALIZE verticalFlip :: RGBA -> RGBA #-} diff --git a/src/Vision/Image/Threshold.hs b/src/Vision/Image/Threshold.hs --- a/src/Vision/Image/Threshold.hs +++ b/src/Vision/Image/Threshold.hs @@ -1,17 +1,24 @@ -{-# LANGUAGE BangPatterns, GADTs #-} +{-# LANGUAGE BangPatterns, FlexibleContexts, GADTs #-} module Vision.Image.Threshold ( ThresholdType (..) , threshold , AdaptiveThresholdKernel (..), adaptiveThreshold + , otsu ) where import Foreign.Storable (Storable) import Vision.Image.Filter (Filter (..), SeparableFilter, blur, gaussianBlur) import Vision.Image.Type (ImagePixel, FunctorImage) +import Vision.Histogram +import Vision.Histogram as H +import Vision.Primitive.Shape (shapeLength) import qualified Vision.Image.Type as I +import qualified Data.Vector.Storable as V +import qualified Data.Vector as VU + -- | Specifies what to do with pixels matching the threshold predicate. -- -- @'BinaryThreshold' a b@ will replace matching pixels by @a@ and non-matchings @@ -46,6 +53,7 @@ => Maybe acc -> AdaptiveThresholdKernel acc -- | Applies a thresholding adaptively. +-- -- Compares every pixel to its surrounding ones in the kernel of the given -- radius. adaptiveThreshold :: (Integral src, Num src, Ord src, Storable acc) @@ -71,3 +79,41 @@ else ifFalse Truncate ifTrue -> if cond then ifTrue else pix {-# INLINE adaptiveThreshold #-} + +-- | Applies a clustering-based image thresholding using the Otsu's method. +-- +-- See . +otsu :: ( HistogramShape (PixelValueSpace (ImagePixel src)) + , ToHistogram (ImagePixel src), FunctorImage src res + , Ord (ImagePixel src), Num (ImagePixel src), Enum (ImagePixel src)) + => ThresholdType (ImagePixel src) (ImagePixel res) -> src -> res +otsu !thresType !img = + threshold (<= thresh) thresType img + where + !thresh = + let hist = histogram Nothing img + histV = H.vector hist + tot = shapeLength (I.shape img) + runningMul = V.zipWith (\v i -> v * i) histV (V.fromList [0..255]) + sm = double (V.sum $ V.drop 1 runningMul) + wB = V.scanl1 (+) histV + wF = V.map (\x -> tot - x) wB + sumB = V.scanl1 (+) runningMul + mB = V.zipWith (\n d -> if d == 0 then 1 + else double n / double d) + sumB wB + mF = V.zipWith (\b f -> if f == 0 then 1 + else (sm - double b) + / double f) + sumB wF + between = V.zipWith4 (\x y b f -> double x * double y + * (b - f)^two) + wB wF mB mF + in snd $ VU.maximum (VU.zip (VU.fromList $ V.toList between) + (VU.fromList [0..255])) + + !two = 2 :: Int +{-# INLINABLE otsu #-} + +double :: Integral a => a -> Double +double = fromIntegral