-- SPDX-License-Identifier: MPL-2.0

{- | Direct ND reference filtering and 2D edge detection. Kernel centers are
at floor(size/2); odd extents are required. All boundary policies are
explicit. Work is limited before neighborhood expansion or output allocation.
-}
module HMath.Image.Filter (
    convolve,
    correlate,
    gaussianBlur,
    boxBlur,
    medianBlur,
    bilateralBlur,
    sharpen,
    GradientOperator (..),
    gradient,
    laplacian,
    canny,
    Morphology (..),
    morphology,
    morphologicalGradient,
    validateBoundary,
    scalar2D,
    windowOffsets,
) where

import Control.Monad (foldM)
import Data.List (sort)
import qualified Data.Map.Strict as M
import qualified Data.Set as Set
import HMath.Error
import HMath.Image.Core

{- | Validate Constant boundary value as finite; Clamp,Reflect,Wrap are always valid.
Non-finite constant returns Left InvalidArgument.
-}
validateBoundary :: Boundary -> Either MathError ()
validateBoundary :: Boundary -> Either MathError ()
validateBoundary (Constant Double
value) = Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
value) String
"nonfinite boundary constant"
validateBoundary Boundary
_ = () -> Either MathError ()
forall a b. b -> Either a b
Right ()

{- | Require a two-dimensional one-channel image and return (height,width). Other
spatial/channel layouts return Left InvalidArgument.
-}
scalar2D :: Image -> Either MathError (Int, Int)
scalar2D :: Image -> Either MathError (Int, Int)
scalar2D Image
img = do
    Bool -> String -> Either MathError ()
require ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
img) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 Bool -> Bool -> Bool
&& Image -> Int
imageChannels Image
img Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1) String
"one-channel 2D image required"
    (Int, Int) -> Either MathError (Int, Int)
forall a b. b -> Either a b
Right (Image -> [Int]
imageShape Image
img [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0, Image -> [Int]
imageShape Image
img [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)

{- | Correlation uses sample(ix+kernelOffset-center). Convolution reverses all
kernel axes; channels are filtered independently without channel mixing.
-}
correlate :: Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
correlate :: Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
correlate Boundary
boundary [Int]
dims [Double]
weights Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    Bool -> String -> Either MathError ()
require
        ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
dims Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
img) Bool -> Bool -> Bool
&& (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
d -> Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd Int
d Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
65) [Int]
dims)
        String
"kernel rank/odd extent invalid"
    let count :: Integer
count = [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a. Integral a => a -> Integer
toInteger [Int]
dims)
    Bool -> String -> Either MathError ()
require
        ( Integer
count Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4096
            Bool -> Bool -> Bool
&& Integer
count Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
weights)
            Bool -> Bool -> Bool
&& Integer
count Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Double]
imageValues Image
img)) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
20000000
        )
        String
"kernel count/work limit"
    _ <- [Double] -> Either MathError [Double]
checked [Double]
weights
    let offsets = ([Int] -> [Int]) -> [[Int]] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Int]
ix -> (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Int]
ix ((Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) [Int]
dims)) ([Int] -> [[Int]]
coordinates [Int]
dims)
        apply [Int]
ix Int
c = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
img ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
ix [Int]
offset) Int
c | ([Int]
offset, Double
w) <- [[Int]] -> [Double] -> [([Int], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Int]]
offsets [Double]
weights]
    build (imageShape img) (imageChannels img) apply

{- | Given boundary,kernel dimensions,flat weights,then image, reverse all kernel axes and
correlate. Channels are filtered independently; require odd extents\<=65, matching
rank/count, \<=4096 weights and \<=20000000 sample operations. Invalid/non-finite
input/result fails.
-}
convolve :: Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
convolve :: Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
convolve Boundary
boundary [Int]
dims [Double]
weights = Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
correlate Boundary
boundary [Int]
dims ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
weights)

-- | Separable Gaussian, truncated at an explicitly supplied integer radius.
gaussianBlur :: Boundary -> Double -> Int -> Image -> Either MathError Image
gaussianBlur :: Boundary -> Double -> Int -> Image -> Either MathError Image
gaussianBlur Boundary
boundary Double
sigma Int
radius Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
sigma Bool -> Bool -> Bool
&& Double
sigma Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Int
radius Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
radius Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
32) String
"invalid Gaussian width/radius"
    let raw :: [Double]
raw = [Double -> Double
forall a. Floating a => a -> a
exp (-Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) | Int
k <- [-Int
radius .. Int
radius]]
        weights :: [Double]
weights = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
raw) [Double]
raw
        rank :: Int
rank = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
img)
    (Image -> Int -> Either MathError Image)
-> Image -> [Int] -> Either MathError Image
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM
        ( \Image
old Int
axis -> Boundary -> [Int] -> [Double] -> Image -> Either MathError Image
correlate Boundary
boundary [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
axis then Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
radius Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 else Int
1 | Int
j <- [Int
0 .. Int
rank Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] [Double]
weights Image
old
        )
        Image
img
        [Int
0 .. Int
rank Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]

{- | Given radius then image, enumerate centered neighborhood offsets. Require radius in
[0,16], \<=4096 neighbors and \<=20000000 sample operations; otherwise fail.
-}
windowOffsets :: Int -> Image -> Either MathError [[Int]]
windowOffsets :: Int -> Image -> Either MathError [[Int]]
windowOffsets Int
radius Image
img = do
    let rank :: Int
rank = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
img); count :: Integer
count = (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
radius Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
rank
    Bool -> String -> Either MathError ()
require
        ( Int
radius Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
radius Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16
            Bool -> Bool -> Bool
&& Integer
count Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4096
            Bool -> Bool -> Bool
&& Integer
count Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Double]
imageValues Image
img)) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
20000000
        )
        String
"neighborhood/work limit"
    [[Int]] -> Either MathError [[Int]]
forall a b. b -> Either a b
Right (([Int] -> [Int]) -> [[Int]] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Num a => a -> a -> a
subtract Int
radius)) ([Int] -> [[Int]]
coordinates (Int -> Int -> [Int]
forall a. Int -> a -> [a]
replicate Int
rank (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
radius Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))))

{- | Given boundary,radius,then image, average a centered full-dimensional box neighborhood
independently per channel. Neighborhood/work bounds and boundary finiteness are
checked; invalid/non-finite output fails.
-}
boxBlur :: Boundary -> Int -> Image -> Either MathError Image
boxBlur :: Boundary -> Int -> Image -> Either MathError Image
boxBlur Boundary
boundary Int
radius Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    offsets <- Int -> Image -> Either MathError [[Int]]
windowOffsets Int
radius Image
img
    build
        (imageShape img)
        (imageChannels img)
        ( \[Int]
ix Int
c ->
            [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                [Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
img ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
ix [Int]
offset) Int
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([[Int]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Int]]
offsets) | [Int]
offset <- [[Int]]
offsets]
        )

{- | Given boundary,radius,then image, choose the median of each centered odd-sized
neighborhood independently per channel. Invalid boundary/neighborhood/work bounds fail.
-}
medianBlur :: Boundary -> Int -> Image -> Either MathError Image
medianBlur :: Boundary -> Int -> Image -> Either MathError Image
medianBlur Boundary
boundary Int
radius Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    offsets <- Int -> Image -> Either MathError [[Int]]
windowOffsets Int
radius Image
img
    build
        (imageShape img)
        (imageChannels img)
        ( \[Int]
ix Int
c ->
            [Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
img ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
ix [Int]
offset) Int
c | [Int]
offset <- [[Int]]
offsets] [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([[Int]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Int]]
offsets Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)
        )

{- | Per-channel bilateral filtering: spatial and range Gaussian weights are
multiplied, normalized locally, and computed relative to the center sample.
-}
bilateralBlur :: Boundary -> Double -> Double -> Int -> Image -> Either MathError Image
bilateralBlur :: Boundary
-> Double -> Double -> Int -> Image -> Either MathError Image
bilateralBlur Boundary
boundary Double
spatialSigma Double
rangeSigma Int
radius Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    Bool -> String -> Either MathError ()
require
        ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double
spatialSigma, Double
rangeSigma])
        String
"positive bilateral scales required"
    offsets <- Int -> Image -> Either MathError [[Int]]
windowOffsets Int
radius Image
img
    let apply [Int]
ix Int
c =
            let center :: Double
center = Image -> [Int] -> Int -> Double
at Image
img [Int]
ix Int
c
                weighted :: [(Double, Double)]
weighted =
                    [ (Double
weight, Double
value)
                    | [Int]
offset <- [[Int]]
offsets
                    , let value :: Double
value = Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
img ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
ix [Int]
offset) Int
c
                    , let weight :: Double
weight =
                            Double -> Double
forall a. Floating a => a -> a
exp
                                ( -Double
0.5
                                    Double -> Double -> Double
forall a. Num a => a -> a -> a
* ( [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
spatialSigma) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Int
d <- [Int]
offset]
                                            Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ((Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
center) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rangeSigma) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
                                      )
                                )
                    ]
                den :: Double
den = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Double)]
weighted)
             in [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v | (Double
w, Double
v) <- [(Double, Double)]
weighted]
    build (imageShape img) (imageChannels img) apply

-- | Unsharp mask x+amount*(x-Gaussian(x)); no implicit range clipping.
sharpen :: Boundary -> Double -> Double -> Int -> Image -> Either MathError Image
sharpen :: Boundary
-> Double -> Double -> Int -> Image -> Either MathError Image
sharpen Boundary
boundary Double
amount Double
sigma Int
radius Image
img = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
amount Bool -> Bool -> Bool
&& Double
amount Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) String
"invalid sharpening amount"
    blurred <- Boundary -> Double -> Int -> Image -> Either MathError Image
gaussianBlur Boundary
boundary Double
sigma Int
radius Image
img
    zipImage (\Double
x Double
y -> Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
amount Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)) img blurred

-- | Image gradient stencil selector, choosing Sobel,Prewitt,or Scharr.
data GradientOperator = Sobel | Prewitt | Scharr deriving (GradientOperator -> GradientOperator -> Bool
(GradientOperator -> GradientOperator -> Bool)
-> (GradientOperator -> GradientOperator -> Bool)
-> Eq GradientOperator
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GradientOperator -> GradientOperator -> Bool
== :: GradientOperator -> GradientOperator -> Bool
$c/= :: GradientOperator -> GradientOperator -> Bool
/= :: GradientOperator -> GradientOperator -> Bool
Eq, Int -> GradientOperator -> ShowS
[GradientOperator] -> ShowS
GradientOperator -> String
(Int -> GradientOperator -> ShowS)
-> (GradientOperator -> String)
-> ([GradientOperator] -> ShowS)
-> Show GradientOperator
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GradientOperator -> ShowS
showsPrec :: Int -> GradientOperator -> ShowS
$cshow :: GradientOperator -> String
show :: GradientOperator -> String
$cshowList :: [GradientOperator] -> ShowS
showList :: [GradientOperator] -> ShowS
Show)

{- | Returns row derivative, column derivative, magnitude, and atan2(row,col)
direction in radians. Smoothing/derivative kernels are normalized so an
interior unit ramp has derivative one.
-}
gradient :: GradientOperator -> Boundary -> Image -> Either MathError (Image, Image, Image, Image)
gradient :: GradientOperator
-> Boundary
-> Image
-> Either MathError (Image, Image, Image, Image)
gradient GradientOperator
operator Boundary
boundary Image
img = do
    _ <- Image -> Either MathError (Int, Int)
scalar2D Image
img
    let smooth = case GradientOperator
operator of GradientOperator
Sobel -> [Double
1, Double
2, Double
1]; GradientOperator
Prewitt -> [Double
1, Double
1, Double
1]; GradientOperator
Scharr -> [Double
3, Double
10, Double
3]
        scale = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
smooth
        rowKernel = [Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale | Double
d <- [-Double
1, Double
0, Double
1], Double
s <- [Double]
smooth]
        colKernel = [Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale | Double
s <- [Double]
smooth, Double
d <- [-Double
1, Double
0, Double
1]]
    dr <- correlate boundary [3, 3] rowKernel img
    dc <- correlate boundary [3, 3] colKernel img
    mag <-
        zipImage
            ( \Double
a Double
b ->
                let m :: Double
m = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
a) (Double -> Double
forall a. Num a => a -> a
abs Double
b)
                 in if Double
m Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
            )
            dr
            dc
    angle <- zipImage atan2 dr dc
    Right (dr, dc, mag, angle)

{- | Given boundary then one-channel 2D image, apply the four-neighbor 3-by-3 Laplacian
[0,1,0;1,-4,1;0,1,0]. Invalid image/boundary or non-finite output fails.
-}
laplacian :: Boundary -> Image -> Either MathError Image
laplacian :: Boundary -> Image -> Either MathError Image
laplacian Boundary
boundary Image
img = do _ <- Image -> Either MathError (Int, Int)
scalar2D Image
img; correlate boundary [3, 3] [0, 1, 0, 1, -4, 1, 0, 1, 0] img

{- | Canny: Gaussian smoothing, Sobel gradient, four-direction nonmaximum
suppression, and 8-connected hysteresis. low/high are absolute normalized
gradient magnitudes, not quantiles. Output is exactly zero or one.
-}
canny :: Boundary -> Double -> Int -> Double -> Double -> Image -> Either MathError Image
canny :: Boundary
-> Double
-> Int
-> Double
-> Double
-> Image
-> Either MathError Image
canny Boundary
boundary Double
sigma Int
radius Double
low Double
high Image
img = do
    (height, width) <- Image -> Either MathError (Int, Int)
scalar2D Image
img
    require (all finite [low, high] && low >= 0 && high >= low) "invalid Canny thresholds"
    smoothed <- gaussianBlur boundary sigma radius img
    (_, _, mag, angle) <- gradient Sobel boundary smoothed
    let suppressed Int
i Int
j =
            let a :: Double
a = Image -> [Int] -> Int -> Double
at Image
angle [Int
i, Int
j] Int
0
                direction :: Int
direction = (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4)) :: Int) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
4
                (Int
di, Int
dj) = case Int
direction of Int
0 -> (Int
0, Int
1); Int
1 -> (Int
1, Int
1); Int
2 -> (Int
1, Int
0); Int
_ -> (Int
1, -Int
1)
                v :: Double
v = Image -> [Int] -> Int -> Double
at Image
mag [Int
i, Int
j] Int
0
             in if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Boundary -> Image -> [Int] -> Int -> Double
sample (Double -> Boundary
Constant Double
0) Image
mag [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
di, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
dj] Int
0 Bool -> Bool -> Bool
&& Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Boundary -> Image -> [Int] -> Int -> Double
sample (Double -> Boundary
Constant Double
0) Image
mag [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
di, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
dj] Int
0
                    then Double
v
                    else Double
0
        values = [((Int, Int), Double)] -> Map (Int, Int) Double
forall k a. Ord k => [(k, a)] -> Map k a
M.fromList [((Int
i, Int
j), Int -> Int -> Double
suppressed Int
i Int
j) | Int
i <- [Int
0 .. Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        weak = [(Int, Int)] -> Set (Int, Int)
forall a. Ord a => [a] -> Set a
Set.fromList [(Int, Int)
ij | ((Int, Int)
ij, Double
v) <- Map (Int, Int) Double -> [((Int, Int), Double)]
forall k a. Map k a -> [(k, a)]
M.toList Map (Int, Int) Double
values, Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
low]
        strong = [(Int, Int)] -> Set (Int, Int)
forall a. Ord a => [a] -> Set a
Set.fromList [(Int, Int)
ij | ((Int, Int)
ij, Double
v) <- Map (Int, Int) Double -> [((Int, Int), Double)]
forall k a. Map k a -> [(k, a)]
M.toList Map (Int, Int) Double
values, Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
high]
        flood [] Set (Int, Int)
visited = Set (Int, Int)
visited
        flood ((Int
i, Int
j) : [(Int, Int)]
queue) Set (Int, Int)
visited =
            let next :: [(Int, Int)]
next =
                    [ (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
di, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
dj)
                    | Int
di <- [-Int
1 .. Int
1]
                    , Int
dj <- [-Int
1 .. Int
1]
                    , (Int
di, Int
dj) (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
/= (Int
0, Int
0)
                    , (Int, Int) -> Set (Int, Int) -> Bool
forall a. Ord a => a -> Set a -> Bool
Set.member (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
di, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
dj) Set (Int, Int)
weak
                    , (Int, Int) -> Set (Int, Int) -> Bool
forall a. Ord a => a -> Set a -> Bool
Set.notMember (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
di, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
dj) Set (Int, Int)
visited
                    ]
             in [(Int, Int)] -> Set (Int, Int) -> Set (Int, Int)
flood ([(Int, Int)]
next [(Int, Int)] -> [(Int, Int)] -> [(Int, Int)]
forall a. [a] -> [a] -> [a]
++ [(Int, Int)]
queue) (((Int, Int) -> Set (Int, Int) -> Set (Int, Int))
-> Set (Int, Int) -> [(Int, Int)] -> Set (Int, Int)
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Int, Int) -> Set (Int, Int) -> Set (Int, Int)
forall a. Ord a => a -> Set a -> Set a
Set.insert Set (Int, Int)
visited [(Int, Int)]
next)
        edges = [(Int, Int)] -> Set (Int, Int) -> Set (Int, Int)
flood (Set (Int, Int) -> [(Int, Int)]
forall a. Set a -> [a]
Set.toList Set (Int, Int)
strong) Set (Int, Int)
strong
    build [height, width] 1 (\[Int]
ix Int
_ -> if (Int, Int) -> Set (Int, Int) -> Bool
forall a. Ord a => a -> Set a -> Bool
Set.member ([Int]
ix [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0, [Int]
ix [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) Set (Int, Int)
edges then Double
1 else Double
0)

-- | Grayscale footprint operation selector: erosion,dilation,opening,or closing.
data Morphology = Erosion | Dilation | Opening | Closing deriving (Morphology -> Morphology -> Bool
(Morphology -> Morphology -> Bool)
-> (Morphology -> Morphology -> Bool) -> Eq Morphology
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Morphology -> Morphology -> Bool
== :: Morphology -> Morphology -> Bool
$c/= :: Morphology -> Morphology -> Bool
/= :: Morphology -> Morphology -> Bool
Eq, Int -> Morphology -> ShowS
[Morphology] -> ShowS
Morphology -> String
(Int -> Morphology -> ShowS)
-> (Morphology -> String)
-> ([Morphology] -> ShowS)
-> Show Morphology
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Morphology -> ShowS
showsPrec :: Int -> Morphology -> ShowS
$cshow :: Morphology -> String
show :: Morphology -> String
$cshowList :: [Morphology] -> ShowS
showList :: [Morphology] -> ShowS
Show)

{- | Flat grayscale morphology with an explicit odd-rank Boolean footprint.
Erosion samples x+offset; dilation samples x-offset (reflected footprint),
so asymmetric footprints obey the usual adjunction. Opening/closing compose
these operators under the same caller-selected boundary extension.
-}
morphology :: Morphology -> Boundary -> [Int] -> [Bool] -> Image -> Either MathError Image
morphology :: Morphology
-> Boundary -> [Int] -> [Bool] -> Image -> Either MathError Image
morphology Morphology
operation Boundary
boundary [Int]
dims [Bool]
footprint Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    Bool -> String -> Either MathError ()
require
        ( [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
dims Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
img)
            Bool -> Bool -> Bool
&& (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Int
d -> Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd Int
d Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
33) [Int]
dims
            Bool -> Bool -> Bool
&& [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a. Integral a => a -> Integer
toInteger [Int]
dims) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Bool] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Bool]
footprint)
            Bool -> Bool -> Bool
&& (Bool -> Bool) -> [Bool] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any Bool -> Bool
forall a. a -> a
id [Bool]
footprint
            Bool -> Bool -> Bool
&& [Bool] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Bool]
footprint Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096
        )
        String
"invalid morphology footprint"
    Bool -> String -> Either MathError ()
require
        (Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Bool] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Bool]
footprint) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Double]
imageValues Image
img)) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
20000000)
        String
"morphology work limit"
    let center :: [Int]
center = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) [Int]
dims
        offsets :: [[Int]]
offsets = [(Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Int]
ix [Int]
center | ([Int]
ix, Bool
True) <- [[Int]] -> [Bool] -> [([Int], Bool)]
forall a b. [a] -> [b] -> [(a, b)]
zip ([Int] -> [[Int]]
coordinates [Int]
dims) [Bool]
footprint]
        erode :: Image -> Either MathError Image
erode Image
old =
            [Int] -> Int -> ([Int] -> Int -> Double) -> Either MathError Image
build
                (Image -> [Int]
imageShape Image
old)
                (Image -> Int
imageChannels Image
old)
                (\[Int]
ix Int
c -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
old ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
ix [Int]
o) Int
c | [Int]
o <- [[Int]]
offsets])
        dilate :: Image -> Either MathError Image
dilate Image
old =
            [Int] -> Int -> ([Int] -> Int -> Double) -> Either MathError Image
build
                (Image -> [Int]
imageShape Image
old)
                (Image -> Int
imageChannels Image
old)
                (\[Int]
ix Int
c -> [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Boundary -> Image -> [Int] -> Int -> Double
sample Boundary
boundary Image
old ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Int]
ix [Int]
o) Int
c | [Int]
o <- [[Int]]
offsets])
    case Morphology
operation of
        Morphology
Erosion -> Image -> Either MathError Image
erode Image
img
        Morphology
Dilation -> Image -> Either MathError Image
dilate Image
img
        Morphology
Opening -> Image -> Either MathError Image
erode Image
img Either MathError Image
-> (Image -> Either MathError Image) -> Either MathError Image
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Image -> Either MathError Image
dilate
        Morphology
Closing -> Image -> Either MathError Image
dilate Image
img Either MathError Image
-> (Image -> Either MathError Image) -> Either MathError Image
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Image -> Either MathError Image
erode

{- | Given boundary,footprint dimensions,Bool footprint,then image, return dilation minus
erosion. Shape/footprint/work/boundary validation and finite output checks propagate
from morphology.
-}
morphologicalGradient :: Boundary -> [Int] -> [Bool] -> Image -> Either MathError Image
morphologicalGradient :: Boundary -> [Int] -> [Bool] -> Image -> Either MathError Image
morphologicalGradient Boundary
boundary [Int]
dims [Bool]
footprint Image
img = do
    dilated <- Morphology
-> Boundary -> [Int] -> [Bool] -> Image -> Either MathError Image
morphology Morphology
Dilation Boundary
boundary [Int]
dims [Bool]
footprint Image
img
    eroded <- morphology Erosion boundary dims footprint img
    zipImage (-) dilated eroded