module Vision.Image.Filter (
Filterable (..), Filter (..)
, BoxFilter, BoxFilter1, SeparableFilter, SeparableFilter1
, KernelAnchor (..)
, Kernel (..)
, SeparableKernel (..), SeparatelyFiltrable (..)
, FilterFold (..), FilterFold1 (..)
, BorderInterpolate (..)
, kernelAnchor, borderInterpolate
, dilate, erode
, blur, gaussianBlur
, Derivative (..), scharr, sobel
) where
import Data.List
import Data.Ratio
import qualified Data.Vector.Storable as V
import Data.Word
import Foreign.Storable (Storable)
import Vision.Image.Type (
MaskedImage (..), Image (..), FromFunction (..)
, Manifest, Delayed
)
import Vision.Primitive (Z (..), (:.) (..), DIM1, DIM2, Size, ix1, ix2)
class Filterable src res f where
apply :: f -> src -> res
data Filter src kernel init acc res = Filter {
fKernelSize :: !Size
, fKernelCenter :: !KernelAnchor
, fKernel :: !kernel
, fInit :: !init
, fPost :: !(src -> acc -> res)
, fInterpol :: !(BorderInterpolate src)
}
type BoxFilter src acc res = Filter src (Kernel src acc) (FilterFold acc)
acc res
type BoxFilter1 src res = Filter src (Kernel src src) FilterFold1 src
res
type SeparableFilter src acc res = Filter src (SeparableKernel src acc)
(FilterFold acc) acc res
type SeparableFilter1 src res = Filter src (SeparableKernel src src)
FilterFold1 src res
data KernelAnchor = KernelAnchor !DIM2 | KernelAnchorCenter
newtype Kernel src acc = Kernel (DIM2 -> src -> acc -> acc)
data SeparableKernel src acc = SeparableKernel {
skVertical :: !(DIM1 -> src -> acc -> acc)
, skHorizontal :: !(DIM1 -> acc -> acc -> acc)
}
class SeparatelyFiltrable src res acc where
type SeparableFilterAccumulator src res acc
instance SeparatelyFiltrable src (Manifest p) acc where
type SeparableFilterAccumulator src (Manifest p) acc = Manifest acc
instance SeparatelyFiltrable src (Delayed p) acc where
type SeparableFilterAccumulator src (Delayed p) acc = Delayed acc
data FilterFold acc = FilterFold acc
data FilterFold1 = FilterFold1
data BorderInterpolate a =
BorderReplicate
| BorderReflect
| BorderWrap
| BorderConstant !a
instance (Image src, FromFunction res, src_p ~ ImagePixel src
, res_p ~ FromFunctionPixel res)
=> Filterable src res (BoxFilter src_p acc res_p) where
apply !(Filter ksize anchor (Kernel kernel) ini post interpol) !img =
let !(FilterFold acc) = ini
in fromFunction size $ \(!pt@(Z :. iy :. ix)) ->
let !iy0 = iy kcy
!ix0 = ix kcx
!safe = iy0 >= 0 && iy0 + kh <= ih
&& ix0 >= 0 && ix0 + kw <= iw
!pix = img `index` pt
in post pix $! if safe then goColumnSafe (iy0 * iw) ix0 0 acc
else goColumn iy0 ix0 0 acc
where
!size@(Z :. ih :. iw) = shape img
!(Z :. kh :. kw) = ksize
!(Z :. kcy :. kcx) = kernelAnchor anchor ksize
goColumn !iy !ix !ky !acc
| ky < kh = case borderInterpolate interpol ih iy of
Left iy' -> goLine iy (iy' * iw) ix ix ky 0 acc
Right val -> goLineConst iy ix ky 0 val acc
| otherwise = acc
goColumnSafe !linearIY !ix !ky !acc
| ky < kh = goLineSafe linearIY ix ix ky 0 acc
| otherwise = acc
goLine !iy !linearIY !ix0 !ix !ky !kx !acc
| kx < kw =
let !val = case borderInterpolate interpol iw ix of
Left ix' -> img `linearIndex` (linearIY + ix')
Right val' -> val'
!acc' = kernel (ix2 ky kx) val acc
in goLine iy linearIY ix0 (ix + 1) ky (kx + 1) acc'
| otherwise = goColumn (iy + 1) ix0 (ky + 1) acc
goLineSafe !linearIY !ix0 !ix !ky !kx !acc
| kx < kw =
let !val = img `linearIndex` (linearIY + ix)
!acc' = kernel (ix2 ky kx) val acc
in goLineSafe linearIY ix0 (ix + 1) ky (kx + 1) acc'
| otherwise = goColumnSafe (linearIY + iw) ix0 (ky + 1) acc
goLineConst !iy !ix !ky !kx !val !acc
| kx < kw = let !acc' = kernel (ix2 ky kx) val acc
in goLineConst iy ix ky (kx + 1) val acc'
| otherwise = goColumn (iy + 1) ix (ky + 1) acc
instance (Image src, FromFunction res, src_p ~ ImagePixel src
, res_p ~ FromFunctionPixel res)
=> Filterable src res (BoxFilter1 src_p res_p) where
apply !(Filter ksize anchor (Kernel kernel) _ post interpol) !img
| kh == 0 || kw == 0 =
error "Using FilterFold1 with an empty kernel."
| otherwise =
fromFunction size $ \(!pt@(Z :. iy :. ix)) ->
let !iy0 = iy kcy
!ix0 = ix kcx
!safe = iy0 >= 0 && iy0 + kh <= ih
&& ix0 >= 0 && ix0 + kw <= iw
!pix = img `index` pt
in post pix $! if safe then goColumn1Safe iy0 ix0
else goColumn1 iy0 ix0
where
!size@(Z :. ih :. iw) = shape img
!(Z :. kh :. kw) = ksize
!(Z :. kcy :. kcx) = kernelAnchor anchor ksize
goColumn1 !iy !ix =
case borderInterpolate interpol ih iy of
Left iy' ->
let !linearIY = iy' * iw
!acc = safeIndex linearIY ix
in goLine iy linearIY ix (ix + 1) 0 1 acc
Right val -> goLineConst iy ix 0 1 val val
goColumn1Safe !iy !ix =
let !linearIY = iy * iw
!acc = img `linearIndex` (linearIY + ix)
in goLineSafe linearIY ix (ix + 1) 0 1 acc
goColumn !iy !ix !ky !acc
| ky < kh = case borderInterpolate interpol ih iy of
Left iy' -> goLine iy (iy' * iw) ix ix ky 0 acc
Right val -> goLineConst iy ix ky 0 val acc
| otherwise = acc
goColumnSafe !linearIY !ix !ky !acc
| ky < kh = goLineSafe linearIY ix ix ky 0 acc
| otherwise = acc
goLine !iy !linearIY !ix0 !ix !ky !kx !acc
| kx < kw =
let !val = safeIndex linearIY ix
!acc' = kernel (ix2 ky kx) val acc
in goLine iy linearIY ix0 (ix + 1) ky (kx + 1) acc'
| otherwise = goColumn (iy + 1) ix0 (ky + 1) acc
goLineSafe !linearIY !ix0 !ix !ky !kx !acc
| kx < kw =
let !val = img `linearIndex` (linearIY + ix)
!acc' = kernel (ix2 ky kx) val acc
in goLineSafe linearIY ix0 (ix + 1) ky (kx + 1) acc'
| otherwise = goColumnSafe (linearIY + iw) ix0 (ky + 1) acc
goLineConst !iy !ix !ky !kx !val !acc
| kx < kw = let !acc' = kernel (ix2 ky kx) val acc
in goLineConst iy ix ky (kx + 1) val acc'
| otherwise = goColumn (iy + 1) ix (ky + 1) acc
safeIndex !linearIY !ix =
case borderInterpolate interpol iw ix of
Left ix' -> img `linearIndex` (linearIY + ix')
Right val -> val
instance (Image src, FromFunction res, SeparatelyFiltrable src res acc
, src_p ~ ImagePixel src, res_p ~ FromFunctionPixel res
, FromFunction (SeparableFilterAccumulator src res acc)
, FromFunctionPixel (SeparableFilterAccumulator src res acc) ~ acc
, Image (SeparableFilterAccumulator src res acc)
, ImagePixel (SeparableFilterAccumulator src res acc) ~ acc)
=> Filterable src res (SeparableFilter src_p acc res_p)
where
apply !f !img =
fst $! wrapper img f
where
wrapper :: (Image src, FromFunction res
, FromFunction (SeparableFilterAccumulator src res acc)
, FromFunctionPixel (SeparableFilterAccumulator src res acc) ~ acc
, Image (SeparableFilterAccumulator src res acc)
, ImagePixel (SeparableFilterAccumulator src res acc) ~ acc)
=> src
-> SeparableFilter (ImagePixel src) acc (FromFunctionPixel res)
-> (res, SeparableFilterAccumulator src res acc)
wrapper !src !(Filter ksize anchor kernel ini post interpol) =
(res, tmp)
where
!size@(Z :. ih :. iw) = shape src
!(Z :. kh :. kw) = ksize
!(Z :. kcy :. kcx) = kernelAnchor anchor ksize
!(SeparableKernel vert horiz) = kernel
!(FilterFold acc0) = ini
!tmp = fromFunction size $ \(!(Z :. iy :. ix)) ->
let !iy0 = iy kcy
in if iy0 >= 0 && iy0 + kh <= ih
then goColumnSafe iy0 ix 0 acc0
else goColumn iy0 ix 0 acc0
!res = fromFunction size $ \(!pt@(Z :. iy :. ix)) ->
let !ix0 = ix kcx
!pix = src `index` pt
in post pix $! if ix0 >= 0 && ix0 + kw <= iw
then goLineSafe (iy * iw) ix0 0 acc0
else goLine (iy * iw) ix0 0 acc0
goColumn !iy !ix !ky !acc
| ky < kh =
let !val = case borderInterpolate interpol ih iy of
Left iy' -> src `index` ix2 iy' ix
Right val' -> val'
!acc' = vert (ix1 ky) val acc
in goColumn (iy + 1) ix (ky + 1) acc'
| otherwise = acc
goColumnSafe !iy !ix !ky !acc
| ky < kh =
let !val = src `index` ix2 iy ix
!acc' = vert (ix1 ky) val acc
in goColumnSafe (iy + 1) ix (ky + 1) acc'
| otherwise = acc
goLine !linearIY !ix !kx !acc
| kx < kw =
let !val =
case borderInterpolate interpol iw ix of
Left ix'-> tmp `linearIndex` (linearIY + ix')
Right _ -> constLine
!acc' = horiz (ix1 kx) val acc
in goLine linearIY (ix + 1) (kx + 1) acc'
| otherwise = acc
goLineSafe !linearIY !ix !kx !acc
| kx < kw =
let !val = tmp `linearIndex` (linearIY + ix)
!acc' = horiz (ix1 kx) val acc
in goLineSafe linearIY (ix + 1) (kx + 1) acc'
| otherwise = acc
constLine | BorderConstant val <- interpol =
foldl' (\acc ky -> vert (ix1 ky) val acc) acc0 [0..kh1]
| otherwise = undefined
instance (Image src, FromFunction res, SeparatelyFiltrable src res src_p
, src_p ~ ImagePixel src, res_p ~ FromFunctionPixel res
, FromFunction (SeparableFilterAccumulator src res src_p)
, FromFunctionPixel (SeparableFilterAccumulator src res src_p) ~ src_p
, Image (SeparableFilterAccumulator src res src_p)
, ImagePixel (SeparableFilterAccumulator src res src_p) ~ src_p)
=> Filterable src res (SeparableFilter1 src_p res_p)
where
apply !f !img =
fst $! wrapper img f
where
wrapper :: (Image src, FromFunction res, acc ~ ImagePixel src
, FromFunction (SeparableFilterAccumulator src res acc)
, FromFunctionPixel (SeparableFilterAccumulator src res acc) ~ acc
, Image (SeparableFilterAccumulator src res acc)
, ImagePixel (SeparableFilterAccumulator src res acc) ~ acc)
=> src
-> SeparableFilter1 (ImagePixel src) (FromFunctionPixel res)
-> (res, SeparableFilterAccumulator src res acc)
wrapper !src !(Filter ksize anchor kernel _ post interpol)
| kh == 0 || kw == 0 =
error "Using FilterFold1 with an empty kernel."
| otherwise =
(res, tmp)
where
!size@(Z :. ih :. iw) = shape src
!(Z :. kh :. kw) = ksize
!(Z :. kcy :. kcx) = kernelAnchor anchor ksize
!(SeparableKernel vert horiz) = kernel
!tmp = fromFunction size $ \(!(Z :. iy :. ix)) ->
let !iy0 = iy kcy
in if iy0 >= 0 && iy0 + kh <= ih
then goColumn1Safe iy0 ix
else goColumn1 iy0 ix
!res = fromFunction size $ \(!pt@(Z :. iy :. ix)) ->
let !ix0 = ix kcx
!pix = src `index` pt
in post pix $! if ix0 >= 0 && ix0 + kw <= iw
then goLine1Safe (iy * iw) ix0
else goLine1 (iy * iw) ix0
goColumn1 !iy !ix =
case borderInterpolate interpol ih iy of
Left iy' ->
let !acc = src `index` ix2 iy' ix
in goColumn (iy + 1) ix 1 acc
Right val ->
goColumn (iy + 1) ix 1 val
goColumn1Safe !iy !ix =
let !linearIY = iy * iw
!acc = src `linearIndex` (linearIY + ix)
in goColumnSafe (linearIY + iw) ix 1 acc
goColumn !iy !ix !ky !acc
| ky < kh =
let !val = case borderInterpolate interpol ih iy of
Left iy' -> src `index` ix2 iy' ix
Right val' -> val'
!acc' = vert (ix1 ky) val acc
in goColumn (iy + 1) ix (ky + 1) acc'
| otherwise = acc
goColumnSafe !linearIY !ix !ky !acc
| ky < kh =
let !val = src `linearIndex` (linearIY + ix)
!acc' = vert (ix1 ky) val acc
in goColumnSafe (linearIY + iw) ix (ky + 1) acc'
| otherwise = acc
goLine1 !linearIY !ix =
let !acc =
case borderInterpolate interpol iw ix of
Left ix' -> tmp `linearIndex` (linearIY + ix')
Right _ -> columnConst
in goLine linearIY (ix + 1) 1 acc
goLine1Safe !linearIY !ix =
let !linearIX = linearIY + ix
!acc = tmp `linearIndex` linearIX
in goLineSafe (linearIX + 1) 1 acc
goLine !linearIY !ix !kx !acc
| kx < kw =
let !val =
case borderInterpolate interpol iw ix of
Left ix'-> tmp `linearIndex` (linearIY + ix')
Right _ -> columnConst
!acc' = horiz (ix1 kx) val acc
in goLine linearIY (ix + 1) (kx + 1) acc'
| otherwise = acc
goLineSafe !linearIX !kx !acc
| kx < kw =
let !val = tmp `linearIndex` linearIX
!acc' = horiz (ix1 kx) val acc
in goLineSafe (linearIX + 1) (kx + 1) acc'
| otherwise = acc
columnConst
| BorderConstant val <- interpol = goColumnConst 1 val val
| otherwise = undefined
goColumnConst !ky !val !acc
| ky < kh = goColumnConst (ky + 1) val (vert (ix1 ky) acc val)
| otherwise = acc
kernelAnchor :: KernelAnchor -> Size -> DIM2
kernelAnchor (KernelAnchor ix) _ = ix
kernelAnchor (KernelAnchorCenter) (Z :. kh :. kw) = ix2 (round $ (kh 1) % 2)
(round $ (kw 1) % 2)
borderInterpolate :: BorderInterpolate a
-> Int
-> Int
-> Either Int a
borderInterpolate !interpol !len !ix
| word ix < word len = Left ix
| otherwise =
case interpol of
BorderReplicate | ix < 0 -> Left 0
| otherwise -> Left $! len 1
BorderReflect -> Left $! goReflect ix
BorderWrap -> Left $! ix `mod` len
BorderConstant i -> Right i
where
goReflect !ix' | ix' < 0 = goReflect (ix' 1)
| ix' >= len = goReflect ((len 1) (ix' len))
| otherwise = ix'
dilate :: Ord src => Int -> SeparableFilter1 src src
dilate radius =
Filter (ix2 size size) KernelAnchorCenter (SeparableKernel kernel kernel)
FilterFold1 (\_ acc -> acc) BorderReplicate
where
!size = radius * 2 + 1
kernel _ = max
erode :: Ord src => Int -> SeparableFilter1 src src
erode radius =
Filter (ix2 size size) KernelAnchorCenter (SeparableKernel kernel kernel)
FilterFold1 (\_ acc -> acc) BorderReplicate
where
!size = radius * 2 + 1
kernel _ = min
blur :: (Integral src, Integral acc, Num res)
=> Int
-> SeparableFilter src acc res
blur radius =
Filter (ix2 size size) KernelAnchorCenter (SeparableKernel vert horiz)
(FilterFold 0) post BorderReplicate
where
!size = radius * 2 + 1
!nPixs = fromIntegral $ square size
vert _ !val !acc = acc + fromIntegral val
horiz _ !acc' !acc = acc + acc'
post _ acc = fromIntegral $ acc `div` nPixs
gaussianBlur :: (Integral src, Floating acc, RealFrac acc, Storable acc
, Integral res)
=> Int
-> Maybe acc
-> SeparableFilter src acc res
gaussianBlur !radius !mSig =
Filter (ix2 size size) KernelAnchorCenter (SeparableKernel vert horiz)
(FilterFold 0) (\_ !acc -> round acc) BorderReplicate
where
!size = radius * 2 + 1
!sig = case mSig of Just s -> s
Nothing -> (0.5 + fromIntegral radius) / 3
vert !(Z :. y) !val !acc = let !coeff = kernelVec V.! y
in acc + fromIntegral val * coeff
horiz !(Z :. x) !val !acc = let !coeff = kernelVec V.! x
in acc + val * coeff
!kernelVec =
let !unormalized = V.generate size $ \x ->
gaussian $! fromIntegral $! abs $! x radius
!kernelSum = V.sum unormalized
in V.map (/ kernelSum) unormalized
gaussian !x = invSigSqrt2Pi * exp (inv2xSig2 * square x)
!invSigSqrt2Pi = 1 / (sig * sqrt (2 * pi))
!inv2xSig2 = 1 / (2 * square sig)
data Derivative = DerivativeX | DerivativeY
scharr :: (Integral src, Integral res)
=> Derivative -> SeparableFilter src res res
scharr der =
let !kernel =
case der of
DerivativeX -> SeparableKernel kernel1 kernel2
DerivativeY -> SeparableKernel kernel2 kernel1
in Filter (ix2 3 3) KernelAnchorCenter kernel (FilterFold 0) (\_ acc -> acc)
BorderReplicate
where
kernel1 !(Z :. 1) !val !acc = acc + 10 * fromIntegral val
kernel1 !(Z :. _) !val !acc = acc + 3 * fromIntegral val
kernel2 !(Z :. 0) !val !acc = acc fromIntegral val
kernel2 !(Z :. 1) !_ !acc = acc
kernel2 !(Z :. ~2) !val !acc = acc + fromIntegral val
sobel :: (Integral src, Integral res, Storable res)
=> Int
-> Derivative
-> SeparableFilter src res res
sobel radius der =
Filter (ix2 size size) KernelAnchorCenter (SeparableKernel vert horiz)
(FilterFold 0) (\_ acc -> acc) BorderReplicate
where
!size = radius * 2 + 1
vert !(Z :. x) !val !acc = let !coeff = vec1 V.! x
in acc + fromIntegral val * coeff
horiz !(Z :. x) !val !acc = let !coeff = vec2 V.! x
in acc + fromIntegral val * coeff
!radius' = fromIntegral radius
(!vec1, !vec2) = case der of DerivativeX -> (vec1', vec2')
DerivativeY -> (vec2', vec1')
!vec1' = let pows = [ 2^i | i <- [0..radius'] ]
in V.fromList $ pows ++ (tail (reverse pows))
!vec2' = V.fromList $ map negate [1..radius'] ++ [0] ++ [1..radius']
square :: Num a => a -> a
square a = a * a
word :: Integral a => a -> Word
word = fromIntegral