{-# LANGUAGE BangPatterns , FlexibleContexts , GADTs #-} module Vision.Image.Threshold ( -- * Simple threshold ThresholdType (..), thresholdType , threshold -- * Adaptive threshold , AdaptiveThresholdKernel (..), AdaptiveThreshold , adaptiveThreshold, adaptiveThresholdFilter -- * Other methods , otsu, scw ) where import Data.Int import Foreign.Storable (Storable) import qualified Data.Vector.Storable as V import qualified Data.Vector as VU import Vision.Image.Class ( Image, ImagePixel, FromFunction (..), FunctorImage, (!), shape ) import Vision.Image.Filter.Internal ( Filter (..), BoxFilter, Kernel (..), SeparableFilter, SeparatelyFiltrable , KernelAnchor (KernelAnchorCenter), FilterFold (..) , BorderInterpolate (BorderReplicate) , apply, blur, gaussianBlur, Mean, mean ) import Vision.Image.Type (Manifest, delayed, manifest) import Vision.Histogram ( HistogramShape, PixelValueSpace, ToHistogram, histogram ) import Vision.Primitive (Z (..), (:.) (..), Size, shapeLength) import qualified Vision.Histogram as H import qualified Vision.Image.Class as I -- | Specifies what to do with pixels matching the threshold predicate. -- -- @'BinaryThreshold' a b@ will replace matching pixels by @a@ and non-matchings -- pixels by @b@. -- -- @'Truncate' a@ will replace matching pixels by @a@. -- -- @'TruncateInv' a@ will replace non-matching pixels by @a@. data ThresholdType src res where BinaryThreshold :: res -> res -> ThresholdType src res Truncate :: src -> ThresholdType src src TruncateInv :: src -> ThresholdType src src -- | Given the thresholding method, a boolean indicating if the pixel match the -- thresholding condition and the pixel, returns the new pixel value. thresholdType :: ThresholdType src res -> Bool -> src -> res thresholdType (BinaryThreshold ifTrue ifFalse) match _ | match = ifTrue | otherwise = ifFalse thresholdType (Truncate ifTrue) match pix | match = ifTrue | otherwise = pix thresholdType (TruncateInv ifFalse) match pix | match = pix | otherwise = ifFalse {-# INLINE thresholdType #-} -- ----------------------------------------------------------------------------- -- | Applies the given predicate and threshold policy on the image. threshold :: FunctorImage src res => (ImagePixel src -> Bool) -> ThresholdType (ImagePixel src) (ImagePixel res) -> src -> res threshold !cond !thresType = I.map (\pix -> thresholdType thresType (cond pix) pix) {-# INLINE threshold #-} -- ----------------------------------------------------------------------------- -- | Defines how pixels of the kernel of the adaptive threshold will be -- weighted. -- -- With 'MeanKernel', pixels of the kernel have the same weight. -- -- With @'GaussianKernel' sigma@, pixels are weighted according to their distance -- from the thresholded pixel using a Gaussian function parametred by @sigma@. -- See 'gaussianBlur' for details. data AdaptiveThresholdKernel acc where MeanKernel :: Integral acc => AdaptiveThresholdKernel acc GaussianKernel :: (Floating acc, RealFrac acc) => Maybe acc -> AdaptiveThresholdKernel acc -- | Compares every pixel to its surrounding ones in the kernel of the given -- radius. adaptiveThreshold :: ( Image src, Integral (ImagePixel src) , Ord (ImagePixel src) , FromFunction res, Integral (FromFunctionPixel res) , Storable acc , SeparatelyFiltrable src res acc) => AdaptiveThresholdKernel acc -> Int -- ^ Kernel radius. -> ImagePixel src -- ^ Minimum difference between the pixel and the kernel -- average. The pixel is thresholded if -- @pixel_value - kernel_mean > difference@ where difference -- is this number. Can be negative. -> ThresholdType (ImagePixel src) (FromFunctionPixel res) -> src -> res adaptiveThreshold kernelType radius thres thresType img = adaptiveThresholdFilter kernelType radius thres thresType `apply` img {-# INLINABLE adaptiveThreshold #-} type AdaptiveThreshold src acc res = SeparableFilter src () acc res -- | Creates an adaptive thresholding filter to be used with 'apply'. -- -- Use 'adaptiveThreshold' if you only want to apply the filter on the image. -- -- Compares every pixel to its surrounding ones in the kernel of the given -- radius. adaptiveThresholdFilter :: (Integral src, Ord src, Storable acc) => AdaptiveThresholdKernel acc -> Int -- ^ Kernel radius. -> src -- ^ Minimum difference between the pixel and the kernel -- average. The pixel is thresholded if -- @pixel_value - kernel_mean > difference@ where -- difference is this number. Can be negative. -> ThresholdType src res -> AdaptiveThreshold src acc res adaptiveThresholdFilter !kernelType !radius !thres !thresType = kernelFilter { fPost = post } where !kernelFilter = case kernelType of MeanKernel -> blur radius GaussianKernel sig -> gaussianBlur radius sig post ix pix ini acc = let !acc' = (fPost kernelFilter) ix pix ini acc !cond = (pix - acc') > thres in thresholdType thresType cond pix {-# INLINE adaptiveThresholdFilter #-} -- ----------------------------------------------------------------------------- -- | 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 #-} -- ----------------------------------------------------------------------------- -- | This is a sliding concentric window filter (SCW) that uses the ratio of the -- standard deviations of two sliding windows centered on a same point to detect -- regions of interest (ROI). -- -- > scw sizeWindowA sizeWindowB beta thresType img -- -- Let @σA@ be the standard deviation of a fist window around a pixel and @σB@ -- be the standard deviation of another window around the same pixel. -- Then the pixel will match the threshold if @σB / σA >= beta@, and will be -- thresholded according to the given 'ThresholdType'. -- -- See . scw :: ( Image src, Integral (ImagePixel src), FromFunction dst , Floating stdev, Fractional stdev, Ord stdev, Storable stdev) => Size -> Size -> stdev -> ThresholdType (ImagePixel src) (FromFunctionPixel dst) -> src -> dst scw !sizeA !sizeB !beta !thresType !img = betaThreshold (stdDev sizeA) (stdDev sizeB) where betaThreshold a b = fromFunction (shape img) $ \pt -> let !cond = (b ! pt) / (a ! pt) < beta in thresholdType thresType cond (img ! pt) stdDev size = let filt :: (Integral src, Fractional res) => Mean src Int16 res filt = mean size !meanImg = manifest $ apply filt img !varImg = manifest $ apply (variance size meanImg) img in delayed $ I.map sqrt varImg {-# INLINABLE scw #-} -- | Given a mean image and an original image, computes the variance of the -- kernel of the given size. -- -- @average [ (origPix - mean)^2 | origPix <- kernel pixels on original ]@. variance :: (Integral src, Fractional res, Storable res) => Size -> Manifest res -> BoxFilter src res res res variance !size@(Z :. h :. w) !meanImg = Filter size KernelAnchorCenter (Kernel kernel) (\pt _ -> meanImg ! pt) (FilterFold (const 0)) post BorderReplicate where kernel !kernelMean _ !val !acc = acc + square (fromIntegral val - kernelMean) !nPixsFactor = 1 / (fromIntegral $! h * w) post _ _ _ !acc = acc * nPixsFactor {-# INLINABLE variance #-} -- ----------------------------------------------------------------------------- square :: Num a => a -> a square a = a * a double :: Integral a => a -> Double double = fromIntegral