-- SPDX-License-Identifier: MPL-2.0

{- | Unit-range, per-channel histograms, equalization, thresholds, and CLAHE.
Values must be in [0,1]; exactly one belongs to the last bin. CLAHE clips
tile counts, redistributes excess without exceeding the feasible clip cap,
and bilinearly interpolates mappings at tile centers.
-}
module HMath.Image.Histogram (histogram, equalizeHistogram, clahe, threshold, otsu, adaptiveThreshold) where

import Data.Array (Array, accumArray, elems, listArray, (!))
import HMath.Error
import HMath.Image.Core
import HMath.Image.Filter (boxBlur, gaussianBlur, validateBoundary)

unit :: Int -> Image -> Either MathError ()
unit :: Int -> Image -> Either MathError ()
unit Int
bins Image
img = do
    Bool -> String -> Either MathError ()
require (Int
bins Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2 Bool -> Bool -> Bool
&& Int
bins Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096) String
"histogram bin count outside [2,4096]"
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) (Image -> [Double]
imageValues Image
img)) String
"unit-range image required"
bin :: Int -> Double -> Int
bin :: Int -> Double -> Int
bin Int
bins Double
x = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
bins))
counts :: Int -> [Double] -> [Int]
counts :: Int -> [Double] -> [Int]
counts Int
bins [Double]
values = Array Int Int -> [Int]
forall i e. Array i e -> [e]
elems ((Int -> Int -> Int)
-> Int -> (Int, Int) -> [(Int, Int)] -> Array Int Int
forall i e a.
Ix i =>
(e -> a -> e) -> e -> (i, i) -> [(i, a)] -> Array i e
accumArray Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 (Int
0, Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [(Int -> Double -> Int
bin Int
bins Double
x, Int
1) | Double
x <- [Double]
values] :: Array Int Int)

{- | Given bin count in [2,4096] then image with every sample in [0,1], return one Int count
list per channel. Values at 1 enter the final bin. Invalid count/range fails.
-}
histogram :: Int -> Image -> Either MathError [[Int]]
histogram :: Int -> Image -> Either MathError [[Int]]
histogram Int
bins Image
img = do
    Int -> Image -> Either MathError ()
unit Int
bins Image
img
    [[Int]] -> Either MathError [[Int]]
forall a b. b -> Either a b
Right
        [Int -> [Double] -> [Int]
counts Int
bins [Image -> [Int] -> Int -> Double
at Image
img [Int]
ix Int
c | [Int]
ix <- [Int] -> [[Int]]
coordinates (Image -> [Int]
imageShape Image
img)] | Int
c <- [Int
0 .. Image -> Int
imageChannels Image
img Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

{- | CDF-min equalization with quantized bin mapping. A single occupied bin
is left unchanged, avoiding an arbitrary jump on constant images.
-}
equalizeHistogram :: Int -> Image -> Either MathError Image
equalizeHistogram :: Int -> Image -> Either MathError Image
equalizeHistogram Int
bins Image
img = do
    hist <- Int -> Image -> Either MathError [[Int]]
histogram Int
bins Image
img
    let maps = [(Int, Int) -> [Double] -> Array Int Double
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0, Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) ([Int] -> [Double]
forall {a} {a}. (Fractional a, Integral a, Ord a) => [a] -> [a]
mapping [Int]
h) | [Int]
h <- [[Int]]
hist]
        mapping [a]
h =
            let cdf :: [a]
cdf = (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1) ((a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0 [a]
h); initial :: a
initial = ([a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ((a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0) [a]
cdf); total :: a
total = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
h
             in if a
initial a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
total
                    then [Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) | Int
i <- [Int
0 .. Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                    else [a -> a -> a
forall a. Ord a => a -> a -> a
max a
0 (a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
initial) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
total a -> a -> a
forall a. Num a => a -> a -> a
- a
initial)) | a
n <- [a]
cdf]
        constant Int
c = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0) ([[Int]]
hist [[Int]] -> Int -> [Int]
forall a. HasCallStack => [a] -> Int -> a
!! Int
c)) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
    build
        (imageShape img)
        (imageChannels img)
        (\[Int]
ix Int
c -> let x :: Double
x = Image -> [Int] -> Int -> Double
at Image
img [Int]
ix Int
c in if Int -> Bool
constant Int
c then Double
x else [Array Int Double]
maps [Array Int Double] -> Int -> Array Int Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
c Array Int Double -> Int -> Double
forall i e. Ix i => Array i e -> i -> e
! Int -> Double -> Int
bin Int
bins Double
x)

{- | 2D adaptive equalization. clipFactor>=1 is relative to uniform tile
occupancy; the cap is max(ceil(tileSamples/bins),floor(factor*tileSamples/bins)).
Tiles at image ends may be smaller. Intensities are quantized without
channel coupling. At most 4,096 tiles and 1,048,576 stored histogram bins.
-}
clahe :: Int -> (Int, Int) -> Double -> Image -> Either MathError Image
clahe :: Int -> (Int, Int) -> Double -> Image -> Either MathError Image
clahe Int
bins (Int
tileHeight, Int
tileWidth) Double
clipFactor Image
img = do
    Int -> Image -> Either MathError ()
unit Int
bins Image
img
    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
&& Int
tileHeight Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
tileWidth Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
clipFactor
            Bool -> Bool -> Bool
&& Double
clipFactor Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
1
            Bool -> Bool -> Bool
&& Double
clipFactor Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e6
        )
        String
"invalid CLAHE parameters"
    let height :: Int
height = Image -> [Int]
imageShape Image
img [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
        width :: Int
width = Image -> [Int]
imageShape Image
img [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1
        nh :: Int
nh = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
tileHeight
        nw :: Int
nw = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
tileWidth
    Bool -> String -> Either MathError ()
require
        ( Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nh Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nw Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4096
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nh Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nw Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
bins Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Image -> Int
imageChannels Image
img) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1048576
        )
        String
"CLAHE histogram storage limit"
    let tile :: Int -> Int -> Int -> Array Int e
tile Int
i Int
j Int
c =
            let values :: [Double]
values =
                    [ Image -> [Int] -> Int -> Double
at Image
img [Int
y, Int
x] Int
c
                    | Int
y <- [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileHeight .. Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
height ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileHeight) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    , Int
x <- [Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileWidth .. Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
width ((Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileWidth) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
                n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
values
                cap :: Int
cap = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max ((Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
bins) (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
n (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
clipFactor Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
bins)))
                clipped :: [Int]
clipped = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
cap) (Int -> [Double] -> [Int]
counts Int
bins [Double]
values)
                redistributed :: [Int]
redistributed = Int -> Int -> [Int] -> [Int]
redistribute Int
cap (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Int]
clipped) [Int]
clipped
             in (Int, Int) -> [e] -> Array Int e
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray
                    (Int
0, Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
                    ((Int -> e) -> [Int] -> [e]
forall a b. (a -> b) -> [a] -> [b]
map ((e -> e -> e
forall a. Fractional a => a -> a -> a
/ Int -> e
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (e -> e) -> (Int -> e) -> Int -> e
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> e
forall a b. (Integral a, Num b) => a -> b
fromIntegral) ((Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
1) ((Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) Int
0 [Int]
redistributed)))
        mappings :: Array Int (Array Int Double)
mappings = Int
-> Int
-> Int
-> (Int -> Int -> Int -> Array Int Double)
-> Array Int (Array Int Double)
forall {t} {e}.
(Ix t, Num t, Enum t) =>
t -> t -> t -> (t -> t -> t -> e) -> Array t e
array3 Int
nh Int
nw (Image -> Int
imageChannels Image
img) Int -> Int -> Int -> Array Int Double
forall {e}. Fractional e => Int -> Int -> Int -> Array Int e
tile
        yCenters :: [Double]
yCenters = [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileHeight Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
height ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileHeight) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 | Int
i <- [Int
0 .. Int
nh Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        xCenters :: [Double]
xCenters = [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileWidth Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
width ((Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
tileWidth) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 | Int
j <- [Int
0 .. Int
nw Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        apply :: [Int] -> Int -> Double
apply [Int
y, Int
x] Int
c =
            let (Int
a, Int
b, Double
v) = [Double] -> Double -> (Int, Int, Double)
forall {c}. (Fractional c, Ord c) => [c] -> c -> (Int, Int, c)
bracket [Double]
yCenters (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y)
                (Int
d, Int
e, Double
u) = [Double] -> Double -> (Int, Int, Double)
forall {c}. (Fractional c, Ord c) => [c] -> c -> (Int, Int, c)
bracket [Double]
xCenters (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x)
                k :: Int
k = Int -> Double -> Int
bin Int
bins (Image -> [Int] -> Int -> Double
at Image
img [Int
y, Int
x] Int
c)
                get :: Int -> Int -> Double
get Int
i Int
j = Array Int (Array Int Double)
mappings Array Int (Array Int Double) -> Int -> Array Int Double
forall i e. Ix i => Array i e -> i -> e
! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
nw Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) Int -> Int -> Int
forall a. Num a => a -> a -> a
* Image -> Int
imageChannels Image
img Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
c) Array Int Double -> Int -> Double
forall i e. Ix i => Array i e -> i -> e
! Int
k
             in (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
v) Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
get Int
a Int
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
get Int
a Int
e) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
get Int
b Int
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
get Int
b Int
e)
        apply [Int]
_ Int
_ = Double
0
    [Int] -> Int -> ([Int] -> Int -> Double) -> Either MathError Image
build [Int
height, Int
width] (Image -> Int
imageChannels Image
img) [Int] -> Int -> Double
apply
  where
    array3 :: t -> t -> t -> (t -> t -> t -> e) -> Array t e
array3 t
nh t
nw t
nc t -> t -> t -> e
f =
        (t, t) -> [e] -> Array t e
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray
            (t
0, t
nh t -> t -> t
forall a. Num a => a -> a -> a
* t
nw t -> t -> t
forall a. Num a => a -> a -> a
* t
nc t -> t -> t
forall a. Num a => a -> a -> a
- t
1)
            [t -> t -> t -> e
f t
i t
j t
c | t
i <- [t
0 .. t
nh t -> t -> t
forall a. Num a => a -> a -> a
- t
1], t
j <- [t
0 .. t
nw t -> t -> t
forall a. Num a => a -> a -> a
- t
1], t
c <- [t
0 .. t
nc t -> t -> t
forall a. Num a => a -> a -> a
- t
1]]
    redistribute :: Int -> Int -> [Int] -> [Int]
redistribute Int
_ Int
0 [Int]
hs = [Int]
hs
    redistribute Int
cap Int
remaining [Int]
hs =
        let available :: Int
available = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
cap) [Int]
hs)
            share :: Int
share = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Int
remaining Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
available)
            step :: Int -> [Int] -> (Int, [Int])
step Int
left [] = (Int
left, [])
            step Int
left (Int
h : [Int]
rest) =
                let delta :: Int
delta = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
left (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
share (Int
cap Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
h)); (Int
final, [Int]
tailValues) = Int -> [Int] -> (Int, [Int])
step (Int
left Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
delta) [Int]
rest
                 in (Int
final, (Int
h Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
delta) Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int]
tailValues)
            (Int
unassigned, [Int]
next) = Int -> [Int] -> (Int, [Int])
step Int
remaining [Int]
hs
         in Int -> Int -> [Int] -> [Int]
redistribute Int
cap Int
unassigned [Int]
next
    bracket :: [c] -> c -> (Int, Int, c)
bracket [c]
centers c
value
        | c
value c -> c -> Bool
forall a. Ord a => a -> a -> Bool
<= ([c] -> Int -> c
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [c]
centers = (Int
0, Int
0, c
0)
        | c
value c -> c -> Bool
forall a. Ord a => a -> a -> Bool
>= [c] -> c
forall a. HasCallStack => [a] -> a
last [c]
centers = ([c] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [c]
centers Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, [c] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [c]
centers Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, c
0)
        | Bool
otherwise =
            let b :: Int
b = [c] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((c -> Bool) -> [c] -> [c]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (c -> c -> Bool
forall a. Ord a => a -> a -> Bool
< c
value) [c]
centers); a :: Int
a = Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
             in (Int
a, Int
b, (c
value c -> c -> c
forall a. Num a => a -> a -> a
- [c]
centers [c] -> Int -> c
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) c -> c -> c
forall a. Fractional a => a -> a -> a
/ ([c]
centers [c] -> Int -> c
forall a. HasCallStack => [a] -> Int -> a
!! Int
b c -> c -> c
forall a. Num a => a -> a -> a
- [c]
centers [c] -> Int -> c
forall a. HasCallStack => [a] -> Int -> a
!! Int
a))

-- | Strict comparison x>level; each channel is thresholded independently.
threshold :: Double -> Image -> Either MathError Image
threshold :: Double -> Image -> Either MathError Image
threshold Double
level Image
img = do Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
level) String
"nonfinite threshold"; (Double -> Double) -> Image -> Either MathError Image
mapImage (\Double
x -> if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
level then Double
1 else Double
0) Image
img

{- | Otsu between-class variance maximum. Lowest threshold wins ties. Returns
one bin-edge threshold per channel; a constant channel returns its value.
-}
otsu :: Int -> Image -> Either MathError ([Double], Image)
otsu :: Int -> Image -> Either MathError ([Double], Image)
otsu Int
bins Image
img = do
    hist <- Int -> Image -> Either MathError [[Int]]
histogram Int
bins Image
img
    let choose Int
c [a]
h =
            let values :: [Double]
values = [Image -> [Int] -> Int -> Double
at Image
img [Int]
ix Int
c | [Int]
ix <- [Int] -> [[Int]]
coordinates (Image -> [Int]
imageShape Image
img)]
                n :: Double
n = a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
h) :: Double
                total :: Double
total = [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
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
count | (Int
i, a
count) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [a]
h]
                choices :: [(Double, Int)]
choices =
                    [ (Double
score, Int
k)
                    | Int
k <- [Int
0 .. Int
bins Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]
                    , let w0 :: Double
w0 = a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
h))
                    , let s0 :: Double
s0 = [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
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
count | (Int
i, a
count) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int .. Int
k] (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
h)]
                    , let score :: Double
score =
                            if Double
w0 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
w0 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
n then -Double
1 else Double
w0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w0) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
s0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
w0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
total Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w0)) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
                    ]
                best :: (Double, Int)
best = ((Double, Int) -> (Double, Int) -> (Double, Int))
-> (Double, Int) -> [(Double, Int)] -> (Double, Int)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\(Double, Int)
a (Double, Int)
b -> if (Double, Int) -> Double
forall a b. (a, b) -> a
fst (Double, Int)
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> (Double, Int) -> Double
forall a b. (a, b) -> a
fst (Double, Int)
a then (Double, Int)
b else (Double, Int)
a) (([(Double, Int)] -> Int -> (Double, Int)
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [(Double, Int)]
choices) ((Int -> [(Double, Int)] -> [(Double, Int)]
forall a. Int -> [a] -> [a]
drop Int
1) [(Double, Int)]
choices)
             in if [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
values Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Double]
values
                    then ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [Double]
values
                    else Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Double, Int) -> Int
forall a b. (a, b) -> b
snd (Double, Int)
best Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
bins
        levels = (Int -> [Int] -> Double) -> [Int] -> [[Int]] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> [Int] -> Double
forall {a}. Integral a => Int -> [a] -> Double
choose [Int
0 ..] [[Int]]
hist
    result <-
        build (imageShape img) (imageChannels img) (\[Int]
ix Int
c -> if Image -> [Int] -> Int -> Double
at Image
img [Int]
ix Int
c Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double]
levels [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
c then Double
1 else Double
0)
    Right (levels, result)

{- | Local mean or Gaussian threshold: x > localMean-offset. Gaussian uses
sigma=max(0.5,radius/2); boundary extension is explicit.
-}
adaptiveThreshold :: Bool -> Boundary -> Int -> Double -> Image -> Either MathError Image
adaptiveThreshold :: Bool
-> Boundary -> Int -> Double -> Image -> Either MathError Image
adaptiveThreshold Bool
gaussian Boundary
boundary Int
radius Double
offset Image
img = do
    Boundary -> Either MathError ()
validateBoundary Boundary
boundary
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
offset) String
"nonfinite adaptive offset"
    smoothed <-
        if Bool
gaussian
            then Boundary -> Double -> Int -> Image -> Either MathError Image
gaussianBlur Boundary
boundary (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0.5 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
radius Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) Int
radius Image
img
            else Boundary -> Int -> Image -> Either MathError Image
boxBlur Boundary
boundary Int
radius Image
img
    zipImage (\Double
x Double
mean -> if Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
mean Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
offset then Double
1 else Double
0) img smoothed