module Data.Image.Binary(
BinaryPixel(..),
toBinaryImage,
(<.), (.<),
(>.), (.>),
(==.), (.==),
(/=.), (./=),
compareImage,
(.<.), (.>.), (.==.), (./=.),
applyMask,
erode, erode',
dilate, dilate',
open, open',
close, close',
label,
areas,
perimeters,
boundingBoxes,
centersOfMass,
distanceTransform,
outline, outline') where
--base>=4
import Control.Monad
import Control.Monad.ST
import Data.STRef
import Data.Monoid
--containers>=0.5.0.0
import qualified Data.Map as M
--vector>=0.10.0.1
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector.Unboxed.Mutable as VM
import qualified Data.Vector as BV
import qualified Data.Vector.Mutable as BVM
import Data.Image.Internal
import Data.Image.Convolution
class BinaryPixel px where
toBinary :: px -> Bool
off :: px
on :: px
toBinaryImage :: (Image img,
BinaryPixel (Pixel img)) => (Pixel img -> Bool) -> img -> img
toBinaryImage pred img@(dimensions -> (rows, cols)) = makeImage rows cols bin where
bin r c = if pred (ref img r c) then on else off
(.<) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img),
Pixel img ~ a) => img -> a -> img
(.<) img num = toBinaryImage pred img where
pred p = p < num
(<.) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img),
Pixel img ~ a) => a -> img -> img
(<.) = flip (.>)
(.>) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img),
Pixel img ~ a) => img -> a -> img
(.>) img num = toBinaryImage pred img where
pred p = p > num
(>.) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img),
Pixel img ~ a) => a -> img -> img
(>.) = flip (.<)
(.==) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img),
Pixel img ~ a) => img -> a -> img
(.==) img num = toBinaryImage pred img where
pred p = p == num
(==.) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img),
Pixel img ~ a) => a -> img -> img
(==.) = flip (.==)
(./=) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img),
Pixel img ~ a) => img -> a -> img
(./=) img num = toBinaryImage pred img where
pred p = p /= num
(/=.) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img),
Pixel img ~ a) => a -> img -> img
(/=.) = flip (./=)
compareImage :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img)) => ((Pixel img) -> (Pixel img) -> Bool) -> img -> img -> img
compareImage comp img0@(dimensions -> (rows, cols)) img1 = makeImage rows cols bin where
bin r c = if p0 `comp` p1 then on else off where
p0 = ref img0 r c
p1 = ref img1 r c
(.>.) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img)) => img -> img -> img
(.>.) = compareImage (>)
(.<.) :: (Image img,
BinaryPixel (Pixel img),
Ord (Pixel img)) => img -> img -> img
(.<.) = compareImage (<)
(.==.) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img)) => img -> img -> img
(.==.) img0@(dimensions -> (rows, cols)) img1 = makeImage rows cols img where
img r c = if p0 == p1 then on else off where
p0 = ref img0 r c
p1 = ref img1 r c
(./=.) :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img)) => img -> img -> img
(./=.) img0@(dimensions -> (rows, cols)) img1 = makeImage rows cols img where
img r c = if p0 /= p1 then on else off where
p0 = ref img0 r c
p1 = ref img1 r c
applyMask :: (Image img,
Monoid (Pixel img),
Image mask,
BinaryPixel (Pixel mask)) => mask -> img -> img
applyMask mask img@(dimensions -> (rows, cols)) = makeImage rows cols mask' where
mask' r c = if (toBinary (ref mask r c)) then (ref img r c) else mempty
erode :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => [[Pixel img]] -> img -> img
erode ls img = (convolve ls img) .== (sum . concat $ ls)
erode' :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => img -> img
erode' = erode [[1,1],[1,1]]
dilate :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => [[Pixel img]] -> img -> img
dilate ls img = toBinaryImage (0 /=) (convolve ls img)
dilate' :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => img -> img
dilate' = dilate [[1,1],[1,1]]
open :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => [[Pixel img]] -> img -> img
open ls = dilate ls . erode ls
open' :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => img -> img
open' = dilate' . erode'
close :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => [[Pixel img]] -> img -> img
close ls = erode ls . dilate ls
close' :: (Image img,
BinaryPixel (Pixel img),
Num (Pixel img),
Eq (Pixel img)) => img -> img
close' = erode' . dilate'
label :: (Image img,
BinaryPixel (Pixel img),
Image img',
Pixel img' ~ Double) => img -> img'
label img@(dimensions -> (rows, cols)) = makeImage rows cols labels where
labels r c = arr V.! (r*cols+c)
arr = getLabels img :: V.Vector Double
getLabels img@(dimensions -> (rows, cols)) = runST $ do
currentLabel <- newSTRef 1 :: ST s (STRef s Int)
equivalences <- newSTRef M.empty :: ST s (STRef s (M.Map Double Double))
labels <- VM.replicate (rows*cols) 0 :: ST s (VM.STVector s Double)
forM_ [ (r, c) | r <- [0..rows1], c <- [0..cols1] ] (\ (r, c) -> do
if toBinary . ref img r $ c then do
smn <- neighbors labels rows cols r c
if (smn == []) then newLabel labels currentLabel cols r c equivalences
else writeLabel labels cols r c smn equivalences
else return ())
eq <- readSTRef equivalences
let eq' = reduceEquivalence . buildEquivalence $ eq
forM_ [ (r, c) | r <- [0..rows1], c <- [0..cols1] ] (\ (r, c) -> do
if toBinary . ref img r $ c then do
currLabel <- VM.read labels (r*cols + c)
let newLabel = (eq' M.! currLabel)
VM.write labels (r*cols + c) newLabel
else return ())
V.freeze labels
reduceEquivalence :: M.Map Double Double -> M.Map Double Double
reduceEquivalence map = reduceEquivalence' (M.keys map) map reduced M.empty where
reduced = reductionMap (M.elems map)
reductionMap :: [Double] -> M.Map Double Double
reductionMap = reductionMap' M.empty 1
reductionMap' :: M.Map Double Double -> Double -> [Double] -> M.Map Double Double
reductionMap' acc _ [] = acc
reductionMap' acc nextLabel (l:ls) = reductionMap' acc' nextLabel' ls where
acc' = if containsL then acc else M.insert l nextLabel acc
nextLabel' = if containsL then nextLabel else (nextLabel+1)
containsL = M.member l acc
reduceEquivalence' :: [Double] -> M.Map Double Double -> M.Map Double Double -> M.Map Double Double -> M.Map Double Double
reduceEquivalence' [] _ _ acc = acc
reduceEquivalence' (k:ks) orig lookup acc = reduceEquivalence' ks orig lookup acc' where
acc' = M.insert k val' acc
val = orig M.! k
val' = lookup M.! val
buildEquivalence :: M.Map Double Double -> M.Map Double Double
buildEquivalence map = buildEquivalence' (M.keys map) map M.empty
buildEquivalence' :: [Double] -> M.Map Double Double -> M.Map Double Double -> M.Map Double Double
buildEquivalence' [] _ acc = acc
buildEquivalence' (k:ks) lookup acc = buildEquivalence' ks lookup (M.insert k equiv acc) where
equiv = findEquivalence lookup k
findEquivalence :: M.Map Double Double -> Double -> Double
findEquivalence eqLookup key = val' where
val = eqLookup M.! key
val' = if val == key then key else findEquivalence eqLookup val
writeLabel labels cols r c smn equiv = do
oldMap <- readSTRef equiv
let min = minimum smn
insert k m = M.insert k min m
newMap = foldr insert oldMap smn
VM.write labels (r*cols + c) min
writeSTRef equiv newMap
newLabel labels currentLabel cols r c equivalences = do
l <- readSTRef currentLabel
modifySTRef currentLabel (+1)
VM.write labels (cols*r + c) (fromIntegral l)
eq <- readSTRef equivalences
let newMap = M.insert (fromIntegral l) (fromIntegral l) eq
writeSTRef equivalences newMap
neighbors :: VM.STVector s Double -> Int -> Int -> Int -> Int -> ST s [Double]
neighbors labels rows cols r c = do
let n' = neighbor labels rows cols
center <- n' r c
north <- n' (r1) c
ne <- n' (r1) (c+1)
east <- n' r (c+1)
se <- n' (r+1) (c+1)
south <- n' (r+1) c
sw <- n' (r+1) (c1)
west <- n' r (c1)
nw <- n' (r1) (c1)
return $ filter (> 0) [center, north, ne, east, se, south, sw, west, nw]
neighbor :: VM.STVector s Double -> Int -> Int -> Int -> Int -> ST s Double
neighbor labels rows cols r c
| r >= 0 && r < rows && c >= 0 && c < cols = VM.read labels (r*cols + c)
| otherwise = return 0.0
areas :: (Image img,
MaxMin (Pixel img),
RealFrac (Pixel img)) => img -> V.Vector Double
areas img@(dimensions -> (rows, cols)) = runST $ do
histogram <- VM.replicate ((floor $ maxIntensity img)+1) 0 :: ST s (VM.STVector s Double)
forM_ (pixelList img) (\ (floor -> p) -> do
currVal <- VM.read histogram p
VM.write histogram p (currVal + 1))
V.freeze histogram
perimeters :: (Image img,
MaxMin (Pixel img),
Pixel img ~ Double) => img -> V.Vector Double
perimeters img@(dimensions -> (rows, cols)) = runST $ do
vector <- VM.replicate ((floor $ maxIntensity img)+1) 0 :: ST s (VM.STVector s Double)
forM_ [ (r, c) | r <- [0..rows1], c <- [0..cols1]] (\ (r, c) -> do
let ns = filter (== 0) . neighborList img r $ c
if null ns then return ()
else do
let (floor -> group) = ref img r c
currPerimeter <- VM.read vector group
VM.write vector group (currPerimeter+1)
)
VM.write vector 0 0
vals <- V.freeze vector
VM.write vector 0 (V.sum vals)
V.freeze vector
neighborList :: (Image img) => img -> Int -> Int -> [Pixel img]
neighborList img@(dimensions -> (rows, cols)) r c =
[ ref img r' c' | r' <- [r1..r+1], c' <- [c1..c+1],
r' /= r && c' /= c, r' >= 0, r' < rows, c' >= 0, c' < cols]
boundingBoxes :: (Image img,
MaxMin (Pixel img),
Pixel img ~ Double) => img -> [(Int, Int, Int, Int)]
boundingBoxes img@(dimensions -> (rows, cols)) = runST $ do
let n = floor . maxIntensity $ img
boxes <- VM.new (n*4) :: ST s (VM.STVector s Int)
forM_ [0..n1] (\ n -> do
VM.write boxes (n*4) rows
VM.write boxes (n*4 + 1) cols
VM.write boxes (n*4 + 2) 0
VM.write boxes (n*4 + 3) 0)
forM_ [ (r, c) | r <- [0..rows1], c <- [0..cols1]] (\ (r, c) ->
let (floor -> px) = ref img r c in updateAt px r c boxes)
boxes' <- V.freeze boxes
return . toQuads . V.toList $ boxes'
updateAt :: Int -> Int -> Int -> VM.STVector s Int -> ST s ()
updateAt ((flip () 1) -> group) r c boxes
| group < 0 = return ()
| otherwise = do
let read = VM.read boxes
write = VM.write boxes
minR <- read (group*4)
minC <- read (group*4 + 1)
maxR <- read (group*4 + 2)
maxC <- read (group*4 + 3)
let minR' = min minR r
minC' = min minC c
maxR' = max maxR r
maxC' = max maxC c
write (group*4) minR'
write (group*4 + 1) minC'
write (group*4 + 2) maxR'
write (group*4 + 3) maxC'
toQuads :: [a] -> [(a,a,a,a)]
toQuads = toQuads' [] where
toQuads' acc [] = reverse acc
toQuads' acc (a:b:c:d:xs) = toQuads' ((a,b,c,d):acc) xs
centersOfMass :: (Image img,
MaxMin (Pixel img),
Pixel img ~ Double) => img -> [(Double, Double)]
centersOfMass img@(dimensions -> (rows, cols)) = runST $ do
let n = floor . maxIntensity $ img
centers <- VM.replicate n (0,0,0) :: ST s (VM.STVector s (Int, Int, Int))
forM_ [ (r, c) | r <- [0..rows1], c <- [0..cols1]] (\ (r, c) ->
let (floor -> px) = ref img r c in updateMass px r c centers)
centers' <- V.freeze centers
return . map averageMass . V.toList $ centers'
updateMass :: Int -> Int -> Int -> VM.STVector s (Int, Int, Int) -> ST s ()
updateMass ((flip () 1) -> group) r c centers
| group < 0 = return ()
| otherwise = do
(pixels, totalRow, totalCol) <- VM.read centers group
VM.write centers group (pixels+1, totalRow+r, totalCol+c)
averageMass :: (Int, Int, Int) -> (Double, Double)
averageMass ((fromIntegral -> total), (fromIntegral -> rows), (fromIntegral -> cols)) = (rows/total, cols/total)
distanceTransform :: (Image img,
BinaryPixel (Pixel img),
Image img',
Pixel img' ~ Double) => img -> img'
distanceTransform img@(dimensions -> (rows, cols)) = makeImage rows cols func
where arr = getDistanceTransformArray img
func r c = arr V.! ((cols*r) + c)
getDistanceTransformArray :: (Image img,
BinaryPixel (Pixel img)) => img -> V.Vector Double
getDistanceTransformArray img@(dimensions -> (rows, cols)) = runST $ do
let size = rows * cols
imgdata = V.fromList . map (boolToDouble . toBinary) . pixelList $ img :: V.Vector Double
on = 10000
let maskRight = [on, 11, on, 11, on, 11, 7, 5, 7, 11, on, 5, 0, on, on]
let maskLeft = reverse maskRight
dtImg <- VM.replicate size 0 :: ST s (VM.STVector s Double)
forM_ [0..size1] $ \ i -> do
let val = if ((imgdata V.! i) == 0) then 0 else on
VM.write dtImg i val
let pass rs cs tr br mask = do
forM_ rs $ \ r -> do
forM_ cs $ \ c -> do
let imgPixels = [ (i, j) | i <- [r+tr..r+br], j <- [c2..c+2]]
pxVals <- forM imgPixels $ \ (i, j) -> do
if (i < 0 || i > (rows1) || j < 0 || j > (cols1))
then return on
else do
val <- VM.read dtImg ((i*cols)+j)
return val
let sums = map (\ (x, y) -> x+y) $ zip pxVals mask
let min = minimum sums
VM.write dtImg ((r*cols) + c) min
pass [0..rows1] [0..cols1] (2) 0 maskRight
pass [rows1,rows2..0] [cols1,cols2..0] 0 2 maskLeft
V.freeze dtImg
boolToDouble :: Bool -> Double
boolToDouble True = 1.0
boolToDouble _ = 0.0
outline :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img)) => img -> img
outline img = outline' off on img
outline' :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img)) => Pixel img -> Pixel img -> img -> img
outline' nonEdge edge img@(dimensions -> (rows, cols)) = makeImage rows cols func
where arr = getOutlineArray img edge nonEdge
func r c = arr BV.! ((cols*r)+c)
getOutlineArray :: (Image img,
BinaryPixel (Pixel img),
Eq (Pixel img)) => img -> Pixel img -> Pixel img -> BV.Vector (Pixel img)
getOutlineArray img@(dimensions -> (rows, cols)) edge nonEdge = runST $ do
let data1 = BV.fromList . map (boolToDouble . toBinary) . pixelList $ img :: BV.Vector Double
data4 <- BVM.replicate (rows*cols) off
forM [0..rows1] $ \ i -> do
let index0 = i*cols
forM [0..cols1] $ \ j -> do
BVM.write data4 (index0+j) nonEdge
forM [0..rows2] $ \ i -> do
let index0 = i*cols
let index1 = (i+1)*cols
forM [0..cols1] $ \ j -> do
let val = (data1 BV.! (index0+j)) + (data1 BV.! (index1+j))
if (val == 1)
then
BVM.write data4 (index0+j) edge
else return ()
let index0 = (rows1)*cols
forM [0..cols1] $ \ j -> do
let val = (data1 BV.! (index0+j)) + (data1 BV.! j)
if (val == 1)
then BVM.write data4 (index0+j) edge
else return ()
forM [0..rows1] $ \ i -> do
let index0 = i*cols
forM [1..cols2] $ \ j -> do
let val = (data1 BV.! (index0+j)) + (data1 BV.! (index0 + j + 1))
if (val == 1)
then BVM.write data4 (index0+j) edge
else return ()
forM [0..rows1] $ \ i -> do
let index0 = i*cols
let val = (data1 BV.! (index0+cols1)) + (data1 BV.! index0)
if (val == 1)
then BVM.write data4 (index0+cols1) edge
else return ()
BV.freeze data4