{-# LANGUAGE BangPatterns
           , FlexibleContexts
           , FlexibleInstances
           , ParallelListComp
           , TypeFamilies
           , TypeOperators #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

-- | Contains functions to compute and manipulate histograms as well as some
-- images transformations which are histogram-based.
--
-- Every polymorphic function is specialised for histograms of 'Int32', 'Double'
-- and 'Float'. Other types can be specialized as every polymorphic function is
-- declared @INLINABLE@.
module Vision.Histogram (
    -- * Types & helpers
      Histogram (..), HistogramShape (..), ToHistogram (..)
    , index, (!), linearIndex, map, assocs, pixToBin
    -- * Histogram computations
    , histogram,  histogram2D, reduce, resize, cumulative, normalize
    -- * Images processing
    , equalizeImage
    -- * Histogram comparisons
    , compareCorrel, compareChi, compareIntersect, compareEMD
    ) where

import Data.Int
import Data.Vector.Storable (Vector)
import Foreign.Storable (Storable)
import Prelude hiding (map)

import qualified Data.Vector.Storable as V

import Vision.Image.Class (Pixel, MaskedImage, Image, ImagePixel, FunctorImage)
import Vision.Image.Grey.Type (GreyPixel (..))
import Vision.Image.HSV.Type  (HSVPixel (..))
import Vision.Image.RGBA.Type (RGBAPixel (..))
import Vision.Image.RGB.Type  (RGBPixel (..))
import Vision.Primitive (
      Z (..), (:.) (..), Shape (..), DIM1, DIM3, DIM4, DIM5, ix1, ix3, ix4
    )

import qualified Vision.Image.Class as I

-- There is no rule to simplify the conversion from Int32 to Double and Float
-- when using realToFrac. Both conversions are using a temporary yet useless
-- Rational value.

{-# RULES
"realToFrac/Int32->Double" realToFrac = fromIntegral :: Int32 -> Double
"realToFrac/Int32->Float"  realToFrac = fromIntegral :: Int32 -> Float
  #-}

-- Types -----------------------------------------------------------------------

data Histogram sh a = Histogram {
      forall sh a. Histogram sh a -> sh
shape  :: !sh
    , forall sh a. Histogram sh a -> Vector a
vector :: !(Vector a) -- ^ Values of the histogram in row-major order.
    } deriving (Histogram sh a -> Histogram sh a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
forall sh a.
(Storable a, Eq sh, Eq a) =>
Histogram sh a -> Histogram sh a -> Bool
/= :: Histogram sh a -> Histogram sh a -> Bool
$c/= :: forall sh a.
(Storable a, Eq sh, Eq a) =>
Histogram sh a -> Histogram sh a -> Bool
== :: Histogram sh a -> Histogram sh a -> Bool
$c== :: forall sh a.
(Storable a, Eq sh, Eq a) =>
Histogram sh a -> Histogram sh a -> Bool
Eq, Histogram sh a -> Histogram sh a -> Bool
Histogram sh a -> Histogram sh a -> Ordering
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall {sh} {a}. (Storable a, Ord sh, Ord a) => Eq (Histogram sh a)
forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Bool
forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Ordering
forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Histogram sh a
min :: Histogram sh a -> Histogram sh a -> Histogram sh a
$cmin :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Histogram sh a
max :: Histogram sh a -> Histogram sh a -> Histogram sh a
$cmax :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Histogram sh a
>= :: Histogram sh a -> Histogram sh a -> Bool
$c>= :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Bool
> :: Histogram sh a -> Histogram sh a -> Bool
$c> :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Bool
<= :: Histogram sh a -> Histogram sh a -> Bool
$c<= :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Bool
< :: Histogram sh a -> Histogram sh a -> Bool
$c< :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Bool
compare :: Histogram sh a -> Histogram sh a -> Ordering
$ccompare :: forall sh a.
(Storable a, Ord sh, Ord a) =>
Histogram sh a -> Histogram sh a -> Ordering
Ord, Int -> Histogram sh a -> ShowS
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
forall sh a.
(Show sh, Show a, Storable a) =>
Int -> Histogram sh a -> ShowS
forall sh a.
(Show sh, Show a, Storable a) =>
[Histogram sh a] -> ShowS
forall sh a.
(Show sh, Show a, Storable a) =>
Histogram sh a -> String
showList :: [Histogram sh a] -> ShowS
$cshowList :: forall sh a.
(Show sh, Show a, Storable a) =>
[Histogram sh a] -> ShowS
show :: Histogram sh a -> String
$cshow :: forall sh a.
(Show sh, Show a, Storable a) =>
Histogram sh a -> String
showsPrec :: Int -> Histogram sh a -> ShowS
$cshowsPrec :: forall sh a.
(Show sh, Show a, Storable a) =>
Int -> Histogram sh a -> ShowS
Show)

-- | Subclass of 'Shape' which defines how to resize a shape so it will fit
-- inside a resized histogram.
class Shape sh => HistogramShape sh where
    -- | Given a number of bins of an histogram, reduces an index so it will be
    -- mapped to a bin.
    toBin :: sh -- ^ The number of bins we are mapping to.
          -> sh -- ^ The number of possible values of the original index.
          -> sh -- ^ The original index.
          -> sh -- ^ The index of the bin in the histogram.

instance HistogramShape Z where
    toBin :: Z -> Z -> Z -> Z
toBin Z
_ Z
_ Z
_ = Z
Z
    {-# INLINE toBin #-}

instance HistogramShape sh => HistogramShape (sh :. Int) where
    toBin :: (sh :. Int) -> (sh :. Int) -> (sh :. Int) -> sh :. Int
toBin !(sh
shBins :. Int
bins) !(sh
shMaxBins :. Int
maxBins) !(sh
shIx :. Int
ix)
        | Int
bins forall a. Eq a => a -> a -> Bool
== Int
maxBins = sh
inner forall tail head. tail -> head -> tail :. head
:. Int
ix
        | Bool
otherwise       = sh
inner forall tail head. tail -> head -> tail :. head
:. (Int
ix forall a. Num a => a -> a -> a
* Int
bins forall a. Integral a => a -> a -> a
`quot` Int
maxBins)
      where
        inner :: sh
inner = forall sh. HistogramShape sh => sh -> sh -> sh -> sh
toBin sh
shBins sh
shMaxBins sh
shIx
    {-# INLINE toBin #-}

-- | This class defines how many dimensions a histogram will have and what will
-- be the default number of bins.
class (Pixel p, Shape (PixelValueSpace p)) => ToHistogram p where
    -- | Gives the value space of a pixel. Single-channel pixels will be 'DIM1'
    -- whereas three-channels pixels will be 'DIM3'.
    -- This is used to determine the rank of the generated histogram.
    type PixelValueSpace p

    -- | Converts a pixel to an index.
    pixToIndex :: p -> PixelValueSpace p

    -- | Returns the maximum number of different values an index can take for
    -- each dimension of the histogram (aka. the maximum index returned by
    -- 'pixToIndex' plus one).
    domainSize :: p -> PixelValueSpace p

instance ToHistogram GreyPixel where
    type PixelValueSpace GreyPixel = DIM1

    pixToIndex :: GreyPixel -> PixelValueSpace GreyPixel
pixToIndex !(GreyPixel Word8
val) = Int -> DIM1
ix1 forall a b. (a -> b) -> a -> b
$ forall a. Integral a => a -> Int
int Word8
val
    {-# INLINE pixToIndex #-}

    domainSize :: GreyPixel -> PixelValueSpace GreyPixel
domainSize GreyPixel
_ = Int -> DIM1
ix1 Int
256

instance ToHistogram RGBAPixel where
    type PixelValueSpace RGBAPixel = DIM4

    pixToIndex :: RGBAPixel -> PixelValueSpace RGBAPixel
pixToIndex !(RGBAPixel Word8
r Word8
g Word8
b Word8
a) = Int -> Int -> Int -> Int -> DIM4
ix4 (forall a. Integral a => a -> Int
int Word8
r) (forall a. Integral a => a -> Int
int Word8
g) (forall a. Integral a => a -> Int
int Word8
b) (forall a. Integral a => a -> Int
int Word8
a)
    {-# INLINE pixToIndex #-}

    domainSize :: RGBAPixel -> PixelValueSpace RGBAPixel
domainSize RGBAPixel
_ = Int -> Int -> Int -> Int -> DIM4
ix4 Int
256 Int
256 Int
256 Int
256

instance ToHistogram RGBPixel where
    type PixelValueSpace RGBPixel = DIM3

    pixToIndex :: RGBPixel -> PixelValueSpace RGBPixel
pixToIndex !(RGBPixel Word8
r Word8
g Word8
b) = Int -> Int -> Int -> DIM3
ix3 (forall a. Integral a => a -> Int
int Word8
r) (forall a. Integral a => a -> Int
int Word8
g) (forall a. Integral a => a -> Int
int Word8
b)
    {-# INLINE pixToIndex #-}

    domainSize :: RGBPixel -> PixelValueSpace RGBPixel
domainSize RGBPixel
_ = Int -> Int -> Int -> DIM3
ix3 Int
256 Int
256 Int
256

instance ToHistogram HSVPixel where
    type PixelValueSpace HSVPixel = DIM3

    pixToIndex :: HSVPixel -> PixelValueSpace HSVPixel
pixToIndex !(HSVPixel Word8
h Word8
s Word8
v) = Int -> Int -> Int -> DIM3
ix3 (forall a. Integral a => a -> Int
int Word8
h) (forall a. Integral a => a -> Int
int Word8
s) (forall a. Integral a => a -> Int
int Word8
v)
    {-# INLINE pixToIndex #-}

    domainSize :: HSVPixel -> PixelValueSpace HSVPixel
domainSize HSVPixel
_ = Int -> Int -> Int -> DIM3
ix3 Int
180 Int
256 Int
256

-- Functions -------------------------------------------------------------------

index :: (Shape sh, Storable a) => Histogram sh a -> sh -> a
index :: forall sh a. (Shape sh, Storable a) => Histogram sh a -> sh -> a
index !Histogram sh a
hist = forall sh a. (Shape sh, Storable a) => Histogram sh a -> Int -> a
linearIndex Histogram sh a
hist forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh. Shape sh => sh -> sh -> Int
toLinearIndex (forall sh a. Histogram sh a -> sh
shape Histogram sh a
hist)
{-# INLINE index #-}

-- | Alias of 'index'.
(!) :: (Shape sh, Storable a) => Histogram sh a -> sh -> a
! :: forall sh a. (Shape sh, Storable a) => Histogram sh a -> sh -> a
(!) = forall sh a. (Shape sh, Storable a) => Histogram sh a -> sh -> a
index
{-# INLINE (!) #-}

-- | Returns the value at the index as if the histogram was a single dimension
-- vector (row-major representation).
linearIndex :: (Shape sh, Storable a) => Histogram sh a -> Int -> a
linearIndex :: forall sh a. (Shape sh, Storable a) => Histogram sh a -> Int -> a
linearIndex !Histogram sh a
hist = forall a. Storable a => Vector a -> Int -> a
(V.!) (forall sh a. Histogram sh a -> Vector a
vector Histogram sh a
hist)
{-# INLINE linearIndex #-}

map :: (Storable a, Storable b) => (a -> b) -> Histogram sh a -> Histogram sh b
map :: forall a b sh.
(Storable a, Storable b) =>
(a -> b) -> Histogram sh a -> Histogram sh b
map a -> b
f !(Histogram sh
sh Vector a
vec) = forall sh a. sh -> Vector a -> Histogram sh a
Histogram sh
sh (forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map a -> b
f Vector a
vec)
{-# INLINE map #-}

-- | Returns all index/value pairs from the histogram.
assocs :: (Shape sh, Storable a) => Histogram sh a -> [(sh, a)]
assocs :: forall sh a. (Shape sh, Storable a) => Histogram sh a -> [(sh, a)]
assocs !(Histogram sh
sh Vector a
vec) = [ (sh
ix, a
v) | sh
ix <- forall sh. Shape sh => sh -> [sh]
shapeList sh
sh
                                       | a
v <- forall a. Storable a => Vector a -> [a]
V.toList Vector a
vec ]
{-# INLINE assocs #-}

-- | Given the number of bins of an histogram and a given pixel, returns the
-- corresponding bin.
pixToBin :: (HistogramShape (PixelValueSpace p), ToHistogram p)
         => PixelValueSpace p -> p -> PixelValueSpace p
pixToBin :: forall p.
(HistogramShape (PixelValueSpace p), ToHistogram p) =>
PixelValueSpace p -> p -> PixelValueSpace p
pixToBin PixelValueSpace p
size p
p =
    let !domain :: PixelValueSpace p
domain = forall p. ToHistogram p => p -> PixelValueSpace p
domainSize p
p
    in forall sh. HistogramShape sh => sh -> sh -> sh -> sh
toBin PixelValueSpace p
size PixelValueSpace p
domain forall a b. (a -> b) -> a -> b
$! forall p. ToHistogram p => p -> PixelValueSpace p
pixToIndex p
p
{-# INLINE pixToBin #-}

-- | Computes an histogram from a (possibly) multi-channel image.
--
-- If the size of the histogram is not given, there will be as many bins as the
-- range of values of pixels of the original image (see 'domainSize').
--
-- If the size of the histogram is specified, every bin of a given dimension
-- will be of the same size (uniform histogram).
histogram :: ( MaskedImage i, ToHistogram (ImagePixel i), Storable a, Num a
             , HistogramShape (PixelValueSpace (ImagePixel i)))
         => Maybe (PixelValueSpace (ImagePixel i)) -> i
         -> Histogram (PixelValueSpace (ImagePixel i)) a
histogram :: forall i a.
(MaskedImage i, ToHistogram (ImagePixel i), Storable a, Num a,
 HistogramShape (PixelValueSpace (ImagePixel i))) =>
Maybe (PixelValueSpace (ImagePixel i))
-> i -> Histogram (PixelValueSpace (ImagePixel i)) a
histogram Maybe (PixelValueSpace (ImagePixel i))
mSize i
img =
    let initial :: Vector a
initial = forall a. Storable a => Int -> a -> Vector a
V.replicate Int
nBins a
0
        ones :: Vector a
ones    = forall a. Storable a => Int -> a -> Vector a
V.replicate Int
nPixs a
1
        ixs :: Vector Int
ixs     = forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map ImagePixel i -> Int
toIndex (forall i. MaskedImage i => i -> Vector (ImagePixel i)
I.values i
img)
    in forall sh a. sh -> Vector a -> Histogram sh a
Histogram PixelValueSpace (ImagePixel i)
size (forall a b.
(Storable a, Storable b) =>
(a -> b -> a) -> Vector a -> Vector Int -> Vector b -> Vector a
V.accumulate_ forall a. Num a => a -> a -> a
(+) Vector a
initial Vector Int
ixs Vector a
ones)
  where
    !size :: PixelValueSpace (ImagePixel i)
size = case Maybe (PixelValueSpace (ImagePixel i))
mSize of Just PixelValueSpace (ImagePixel i)
s  -> PixelValueSpace (ImagePixel i)
s
                          Maybe (PixelValueSpace (ImagePixel i))
Nothing -> forall p. ToHistogram p => p -> PixelValueSpace p
domainSize (forall i. MaskedImage i => i -> ImagePixel i
I.pixel i
img)
    !nChans :: Int
nChans = forall i. (Pixel (ImagePixel i), MaskedImage i) => i -> Int
I.nChannels i
img
    !nPixs :: Int
nPixs = forall sh. Shape sh => sh -> Int
shapeLength (forall i. MaskedImage i => i -> DIM1 :. Int
I.shape i
img) forall a. Num a => a -> a -> a
* Int
nChans
    !nBins :: Int
nBins = forall sh. Shape sh => sh -> Int
shapeLength PixelValueSpace (ImagePixel i)
size

    toIndex :: ImagePixel i -> Int
toIndex !ImagePixel i
p = forall sh. Shape sh => sh -> sh -> Int
toLinearIndex PixelValueSpace (ImagePixel i)
size forall a b. (a -> b) -> a -> b
$!
        case Maybe (PixelValueSpace (ImagePixel i))
mSize of Just PixelValueSpace (ImagePixel i)
_  -> forall p.
(HistogramShape (PixelValueSpace p), ToHistogram p) =>
PixelValueSpace p -> p -> PixelValueSpace p
pixToBin   PixelValueSpace (ImagePixel i)
size ImagePixel i
p
                      Maybe (PixelValueSpace (ImagePixel i))
Nothing -> forall p. ToHistogram p => p -> PixelValueSpace p
pixToIndex ImagePixel i
p
    {-# INLINE toIndex #-}
{-# INLINABLE histogram #-}

-- | Similar to 'histogram' but adds two dimensions for the y and x-coordinates
-- of the sampled points. This way, the histogram will map different regions of
-- the original image.
--
-- For example, an 'RGB' image will be mapped as
-- @'Z' ':.' red channel ':.' green channel ':.' blue channel ':.' y region
-- ':.' x region@.
--
-- As there is no reason to create an histogram as large as the number of pixels
-- of the image, a size is always needed.
histogram2D :: ( Image i, ToHistogram (ImagePixel i), Storable a, Num a
               , HistogramShape (PixelValueSpace (ImagePixel i)))
            => (PixelValueSpace (ImagePixel i)) :. Int :. Int -> i
            -> Histogram ((PixelValueSpace (ImagePixel i)) :. Int :. Int) a
histogram2D :: forall i a.
(Image i, ToHistogram (ImagePixel i), Storable a, Num a,
 HistogramShape (PixelValueSpace (ImagePixel i))) =>
((PixelValueSpace (ImagePixel i) :. Int) :. Int)
-> i
-> Histogram ((PixelValueSpace (ImagePixel i) :. Int) :. Int) a
histogram2D (PixelValueSpace (ImagePixel i) :. Int) :. Int
size i
img =
    let initial :: Vector a
initial = forall a. Storable a => Int -> a -> Vector a
V.replicate Int
nBins a
0
        ones :: Vector a
ones    = forall a. Storable a => Int -> a -> Vector a
V.replicate Int
nPixs a
1
        imgIxs :: Vector (DIM1 :. Int)
imgIxs  = forall a. Storable a => Int -> (a -> a) -> a -> Vector a
V.iterateN Int
nPixs (forall sh. Shape sh => sh -> sh -> sh
shapeSucc DIM1 :. Int
imgSize) forall sh. Shape sh => sh
shapeZero
        ixs :: Vector Int
ixs     = forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (DIM1 :. Int) -> ImagePixel i -> Int
toIndex Vector (DIM1 :. Int)
imgIxs (forall i. Image i => i -> Vector (ImagePixel i)
I.vector i
img)
    in forall sh a. sh -> Vector a -> Histogram sh a
Histogram (PixelValueSpace (ImagePixel i) :. Int) :. Int
size (forall a b.
(Storable a, Storable b) =>
(a -> b -> a) -> Vector a -> Vector Int -> Vector b -> Vector a
V.accumulate_ forall a. Num a => a -> a -> a
(+) Vector a
initial Vector Int
ixs Vector a
ones)
  where
    !imgSize :: DIM1 :. Int
imgSize@(Z
Z :. Int
h :. Int
w) = forall i. MaskedImage i => i -> DIM1 :. Int
I.shape i
img
    !maxSize :: (PixelValueSpace (ImagePixel i) :. Int) :. Int
maxSize = forall p. ToHistogram p => p -> PixelValueSpace p
domainSize (forall i. MaskedImage i => i -> ImagePixel i
I.pixel i
img) forall tail head. tail -> head -> tail :. head
:. Int
h forall tail head. tail -> head -> tail :. head
:. Int
w
    !nChans :: Int
nChans = forall i. (Pixel (ImagePixel i), MaskedImage i) => i -> Int
I.nChannels i
img
    !nPixs :: Int
nPixs = forall sh. Shape sh => sh -> Int
shapeLength (forall i. MaskedImage i => i -> DIM1 :. Int
I.shape i
img) forall a. Num a => a -> a -> a
* Int
nChans
    !nBins :: Int
nBins = forall sh. Shape sh => sh -> Int
shapeLength (PixelValueSpace (ImagePixel i) :. Int) :. Int
size

    toIndex :: (DIM1 :. Int) -> ImagePixel i -> Int
toIndex !(Z
Z :. Int
y :. Int
x) !ImagePixel i
p =
        let !ix :: (PixelValueSpace (ImagePixel i) :. Int) :. Int
ix = (forall p. ToHistogram p => p -> PixelValueSpace p
pixToIndex ImagePixel i
p) forall tail head. tail -> head -> tail :. head
:. Int
y forall tail head. tail -> head -> tail :. head
:. Int
x
        in forall sh. Shape sh => sh -> sh -> Int
toLinearIndex (PixelValueSpace (ImagePixel i) :. Int) :. Int
size forall a b. (a -> b) -> a -> b
$! forall sh. HistogramShape sh => sh -> sh -> sh -> sh
toBin (PixelValueSpace (ImagePixel i) :. Int) :. Int
size (PixelValueSpace (ImagePixel i) :. Int) :. Int
maxSize (PixelValueSpace (ImagePixel i) :. Int) :. Int
ix
    {-# INLINE toIndex #-}
{-# INLINABLE histogram2D #-}

-- Reshaping -------------------------------------------------------------------

-- | Reduces a 2D histogram to its linear representation. See 'resize' for a
-- reduction of the number of bins of an histogram.
--
-- @'histogram' == 'reduce' . 'histogram2D'@
reduce :: (HistogramShape sh, Storable a, Num a)
       => Histogram (sh :. Int :. Int) a -> Histogram sh a
reduce :: forall sh a.
(HistogramShape sh, Storable a, Num a) =>
Histogram ((sh :. Int) :. Int) a -> Histogram sh a
reduce !(Histogram (sh :. Int) :. Int
sh Vector a
vec) =
    let !(sh
sh' :. Int
h :. Int
w) = (sh :. Int) :. Int
sh
        !len2D :: Int
len2D = Int
h forall a. Num a => a -> a -> a
* Int
w
        !vec' :: Vector a
vec' = forall a b.
Storable a =>
Int -> (b -> Maybe (a, b)) -> b -> Vector a
V.unfoldrN (forall sh. Shape sh => sh -> Int
shapeLength sh
sh') Vector a -> Maybe (a, Vector a)
step Vector a
vec
        step :: Vector a -> Maybe (a, Vector a)
step !Vector a
rest = let (!Vector a
channels, !Vector a
rest') = forall a. Storable a => Int -> Vector a -> (Vector a, Vector a)
V.splitAt Int
len2D Vector a
rest
                     in forall a. a -> Maybe a
Just (forall a. (Storable a, Num a) => Vector a -> a
V.sum Vector a
channels, Vector a
rest')
    in forall sh a. sh -> Vector a -> Histogram sh a
Histogram sh
sh' Vector a
vec'
{-# SPECIALIZE reduce :: Histogram DIM5 Int32  -> Histogram DIM3 Int32
                      ,  Histogram DIM5 Double -> Histogram DIM3 Double
                      ,  Histogram DIM5 Float  -> Histogram DIM3 Float
                      ,  Histogram DIM3 Int32  -> Histogram DIM1 Int32
                      ,  Histogram DIM3 Double -> Histogram DIM1 Double
                      ,  Histogram DIM3 Float  -> Histogram DIM1 Float #-}
{-# INLINABLE reduce #-}

-- | Resizes an histogram to another index shape. See 'reduce' for a reduction
-- of the number of dimensions of an histogram.
resize :: (HistogramShape sh, Storable a, Num a)
       => sh -> Histogram sh a -> Histogram sh a
resize :: forall sh a.
(HistogramShape sh, Storable a, Num a) =>
sh -> Histogram sh a -> Histogram sh a
resize !sh
sh' (Histogram sh
sh Vector a
vec) =
    let initial :: Vector a
initial = forall a. Storable a => Int -> a -> Vector a
V.replicate (forall sh. Shape sh => sh -> Int
shapeLength sh
sh') a
0
        -- TODO: In this scheme, indexes are computed for each bin of the
        -- original histogram. It's sub-optimal as some parts of those indexes
        -- (lower dimensions) don't change at each bin.
        reIndex :: Int -> Int
reIndex = forall sh. Shape sh => sh -> sh -> Int
toLinearIndex sh
sh' forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh. HistogramShape sh => sh -> sh -> sh -> sh
toBin sh
sh' sh
sh forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall sh. Shape sh => sh -> Int -> sh
fromLinearIndex sh
sh
        ixs :: Vector Int
ixs = forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map Int -> Int
reIndex forall a b. (a -> b) -> a -> b
$ forall a. (Storable a, Num a) => a -> Int -> Vector a
V.enumFromN Int
0 (forall sh. Shape sh => sh -> Int
shapeLength sh
sh)
    in forall sh a. sh -> Vector a -> Histogram sh a
Histogram sh
sh' (forall a b.
(Storable a, Storable b) =>
(a -> b -> a) -> Vector a -> Vector Int -> Vector b -> Vector a
V.accumulate_ forall a. Num a => a -> a -> a
(+) Vector a
initial Vector Int
ixs Vector a
vec)

-- Normalisation ---------------------------------------------------------------

-- | Computes the cumulative histogram of another single dimension histogram.
--
-- @C(i) = SUM H(j)@ for each @j@ in @[0..i]@ where @C@ is the cumulative
-- histogram, and @H@ the original histogram.
cumulative :: (Storable a, Num a) => Histogram DIM1 a -> Histogram DIM1 a
cumulative :: forall a.
(Storable a, Num a) =>
Histogram DIM1 a -> Histogram DIM1 a
cumulative (Histogram DIM1
sh Vector a
vec) = forall sh a. sh -> Vector a -> Histogram sh a
Histogram DIM1
sh (forall a. Storable a => (a -> a -> a) -> Vector a -> Vector a
V.scanl1' forall a. Num a => a -> a -> a
(+) Vector a
vec)
{-# SPECIALIZE cumulative :: Histogram DIM1 Int32  -> Histogram DIM1 Int32
                          ,  Histogram DIM1 Double -> Histogram DIM1 Double
                          ,  Histogram DIM1 Float  -> Histogram DIM1 Float #-}
{-# INLINABLE cumulative #-}

-- | Normalizes the histogram so that the sum of the histogram bins is equal to
-- the given value (normalisation by the @L1@ norm).
--
-- This is useful to compare two histograms which have been computed from images
-- with a different number of pixels.
normalize :: (Storable a, Real a, Storable b, Fractional b)
          => b -> Histogram sh a -> Histogram sh b
normalize :: forall a b sh.
(Storable a, Real a, Storable b, Fractional b) =>
b -> Histogram sh a -> Histogram sh b
normalize b
norm !hist :: Histogram sh a
hist@(Histogram sh
_ Vector a
vec) =
    let !ratio :: b
ratio = b
norm forall a. Fractional a => a -> a -> a
/ forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. (Storable a, Num a) => Vector a -> a
V.sum Vector a
vec)
        equalizeVal :: a -> b
equalizeVal !a
val = forall a b. (Real a, Fractional b) => a -> b
realToFrac a
val forall a. Num a => a -> a -> a
* b
ratio
        {-# INLINE equalizeVal #-}
    in forall a b sh.
(Storable a, Storable b) =>
(a -> b) -> Histogram sh a -> Histogram sh b
map a -> b
equalizeVal Histogram sh a
hist
{-# SPECIALIZE normalize :: Double -> Histogram sh Int32  -> Histogram sh Double
                         ,  Float  -> Histogram sh Int32  -> Histogram sh Float
                         ,  Double -> Histogram sh Double -> Histogram sh Double
                         ,  Float  -> Histogram sh Double -> Histogram sh Float
                         ,  Double -> Histogram sh Float  -> Histogram sh Double
                         ,  Float  -> Histogram sh Float  -> Histogram sh Float
                         #-}
{-# INLINABLE normalize #-}

-- | Equalizes a single channel image by equalising its histogram.
--
-- The algorithm equalizes the brightness and increases the contrast of the
-- image by mapping each pixel values to the value at the index of the
-- cumulative @L1@-normalized histogram :
--
-- @N(x, y) = H(I(x, y))@ where @N@ is the equalized image, @I@ is the image and
-- @H@ the cumulative of the histogram normalized over an @L1@ norm.
--
-- See <https://en.wikipedia.org/wiki/Histogram_equalization>.
equalizeImage :: ( FunctorImage i i, Integral (ImagePixel i)
                 , ToHistogram (ImagePixel i)
                 , PixelValueSpace (ImagePixel i) ~ DIM1)
              => i -> i
equalizeImage :: forall i.
(FunctorImage i i, Integral (ImagePixel i),
 ToHistogram (ImagePixel i),
 PixelValueSpace (ImagePixel i) ~ DIM1) =>
i -> i
equalizeImage i
img =
    forall src res.
FunctorImage src res =>
(ImagePixel src -> ImagePixel res) -> src -> res
I.map ImagePixel i -> ImagePixel i
equalizePixel i
img
  where
    hist :: Histogram DIM1 Int32
hist            = forall i a.
(MaskedImage i, ToHistogram (ImagePixel i), Storable a, Num a,
 HistogramShape (PixelValueSpace (ImagePixel i))) =>
Maybe (PixelValueSpace (ImagePixel i))
-> i -> Histogram (PixelValueSpace (ImagePixel i)) a
histogram forall a. Maybe a
Nothing i
img             :: Histogram DIM1 Int32
    Z
Z :. Int
nBins      = forall sh a. Histogram sh a -> sh
shape Histogram DIM1 Int32
hist
    cumNormalized :: Histogram DIM1 Double
cumNormalized   = forall a.
(Storable a, Num a) =>
Histogram DIM1 a -> Histogram DIM1 a
cumulative forall a b. (a -> b) -> a -> b
$ forall a b sh.
(Storable a, Real a, Storable b, Fractional b) =>
b -> Histogram sh a -> Histogram sh b
normalize (forall a. Integral a => a -> Double
double Int
nBins) Histogram DIM1 Int32
hist
    !cumNormalized' :: Histogram DIM1 Int32
cumNormalized' = forall a b sh.
(Storable a, Storable b) =>
(a -> b) -> Histogram sh a -> Histogram sh b
map forall a b. (RealFrac a, Integral b) => a -> b
round Histogram DIM1 Double
cumNormalized           :: Histogram DIM1 Int32
    equalizePixel :: ImagePixel i -> ImagePixel i
equalizePixel !ImagePixel i
val = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ Histogram DIM1 Int32
cumNormalized' forall sh a. (Shape sh, Storable a) => Histogram sh a -> sh -> a
! Int -> DIM1
ix1 (forall a. Integral a => a -> Int
int ImagePixel i
val)
    {-# INLINE equalizePixel #-}
{-# INLINABLE equalizeImage #-}

-- Comparisons -----------------------------------------------------------------

-- | Computes the /Pearson\'s correlation coefficient/ between each
-- corresponding bins of the two histograms.
--
-- A value of 1 implies a perfect correlation, a value of -1 a perfect
-- opposition and a value of 0 no correlation at all.
--
-- @'compareCorrel' = SUM  [ (H1(i) - µ(H1)) (H1(2) - µ(H2)) ]
--                  / (   SQRT [ SUM [ (H1(i) - µ(H1))^2 ] ]
--                      * SQRT [ SUM [ (H2(i) - µ(H2))^2 ] ] )@
--
-- Where @µ(H)@ is the average value of the histogram @H@.
--
-- See <http://en.wikipedia.org/wiki/Pearson_correlation_coefficient>.
compareCorrel :: (Shape sh, Storable a, Real a, Storable b, Eq b, Floating b)
              => Histogram sh a -> Histogram sh a -> b
compareCorrel :: forall sh a b.
(Shape sh, Storable a, Real a, Storable b, Eq b, Floating b) =>
Histogram sh a -> Histogram sh a -> b
compareCorrel (Histogram sh
sh1 Vector a
vec1) (Histogram sh
sh2 Vector a
vec2)
    | sh
sh1 forall a. Eq a => a -> a -> Bool
/= sh
sh2     = forall a. HasCallStack => String -> a
error String
"Histograms are not of equal size."
    | b
denominat forall a. Eq a => a -> a -> Bool
== b
0 = b
1
    | Bool
otherwise      = b
numerat forall a. Fractional a => a -> a -> a
/ b
denominat
  where
    numerat :: b
numerat   = forall a. (Storable a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(*) Vector b
diff1 Vector b
diff2
    denominat :: b
denominat =   forall a. Floating a => a -> a
sqrt (forall a. (Storable a, Num a) => Vector a -> a
V.sum (forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map forall a. Num a => a -> a
square Vector b
diff1))
                forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (forall a. (Storable a, Num a) => Vector a -> a
V.sum (forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map forall a. Num a => a -> a
square Vector b
diff2))

    diff1 :: Vector b
diff1 = forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map (\a
v1 -> forall a b. (Real a, Fractional b) => a -> b
realToFrac a
v1 forall a. Num a => a -> a -> a
- b
avg1) Vector a
vec1
    diff2 :: Vector b
diff2 = forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
V.map (\a
v2 -> forall a b. (Real a, Fractional b) => a -> b
realToFrac a
v2 forall a. Num a => a -> a -> a
- b
avg2) Vector a
vec2

    (b
avg1, b
avg2) = (forall {a} {a}. (Fractional a, Real a, Storable a) => Vector a -> a
avg Vector a
vec1, forall {a} {a}. (Fractional a, Real a, Storable a) => Vector a -> a
avg Vector a
vec2)
    avg :: Vector a -> a
avg !Vector a
vec = forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. (Storable a, Num a) => Vector a -> a
V.sum Vector a
vec) forall a. Fractional a => a -> a -> a
/ forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Storable a => Vector a -> Int
V.length Vector a
vec)
{-# SPECIALIZE compareCorrel
    :: Shape sh => Histogram sh Int32  -> Histogram sh Int32  -> Double
    ,  Shape sh => Histogram sh Int32  -> Histogram sh Int32  -> Float
    ,  Shape sh => Histogram sh Double -> Histogram sh Double -> Double
    ,  Shape sh => Histogram sh Double -> Histogram sh Double -> Float
    ,  Shape sh => Histogram sh Float  -> Histogram sh Float  -> Double
    ,  Shape sh => Histogram sh Float  -> Histogram sh Float  -> Float  #-}
{-# INLINABLE compareCorrel #-}

-- | Computes the Chi-squared distance between two histograms.
--
-- A value of 0 indicates a perfect match.
--
-- @'compareChi' = SUM (d(i))@ for each indice @i@ of the histograms where
-- @d(i) = 2 * ((H1(i) - H2(i))^2 / (H1(i) + H2(i)))@.
compareChi :: (Shape sh, Storable a, Real a, Storable b, Fractional b)
           => Histogram sh a -> Histogram sh a -> b
compareChi :: forall sh a b.
(Shape sh, Storable a, Real a, Storable b, Fractional b) =>
Histogram sh a -> Histogram sh a -> b
compareChi (Histogram sh
sh1 Vector a
vec1) (Histogram sh
sh2 Vector a
vec2)
    | sh
sh1 forall a. Eq a => a -> a -> Bool
/= sh
sh2 = forall a. HasCallStack => String -> a
error String
"Histograms are not of equal size."
    | Bool
otherwise  = (forall a. (Storable a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall {a} {a}. (Fractional a, Real a) => a -> a -> a
step Vector a
vec1 Vector a
vec2) forall a. Num a => a -> a -> a
* b
2
  where
    step :: a -> a -> a
step !a
v1 !a
v2 = let !denom :: a
denom = a
v1 forall a. Num a => a -> a -> a
+ a
v2
                   in if a
denom forall a. Eq a => a -> a -> Bool
== a
0
                        then a
0
                        else forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Num a => a -> a
square (a
v1 forall a. Num a => a -> a -> a
- a
v2)) forall a. Fractional a => a -> a -> a
/ forall a b. (Real a, Fractional b) => a -> b
realToFrac a
denom
    {-# INLINE step #-}
{-# SPECIALIZE compareChi
    :: Shape sh => Histogram sh Int32  -> Histogram sh Int32  -> Double
    ,  Shape sh => Histogram sh Int32  -> Histogram sh Int32  -> Float
    ,  Shape sh => Histogram sh Double -> Histogram sh Double -> Double
    ,  Shape sh => Histogram sh Double -> Histogram sh Double -> Float
    ,  Shape sh => Histogram sh Float  -> Histogram sh Float  -> Double
    ,  Shape sh => Histogram sh Float  -> Histogram sh Float  -> Float  #-}
{-# INLINABLE compareChi #-}

-- | Computes the intersection of the two histograms.
--
-- The higher the score is, the best the match is.
--
-- @'compareIntersect' = SUM (min(H1(i), H2(i))@ for each indice @i@ of the
-- histograms.
compareIntersect :: (Shape sh, Storable a, Num a, Ord a)
                 => Histogram sh a -> Histogram sh a -> a
compareIntersect :: forall sh a.
(Shape sh, Storable a, Num a, Ord a) =>
Histogram sh a -> Histogram sh a -> a
compareIntersect (Histogram sh
sh1 Vector a
vec1) (Histogram sh
sh2 Vector a
vec2)
    | sh
sh1 forall a. Eq a => a -> a -> Bool
/= sh
sh2 = forall a. HasCallStack => String -> a
error String
"Histograms are not of equal size."
    | Bool
otherwise  = forall a. (Storable a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Ord a => a -> a -> a
min Vector a
vec1 Vector a
vec2
{-# SPECIALIZE compareIntersect
    :: Shape sh => Histogram sh Int32  -> Histogram sh Int32  -> Int32
    ,  Shape sh => Histogram sh Double -> Histogram sh Double -> Double
    ,  Shape sh => Histogram sh Float  -> Histogram sh Float  -> Float #-}
{-# INLINABLE compareIntersect #-}

-- | Computed the /Earth mover's distance/ between two histograms.
--
-- Current algorithm only supports histograms of one dimension.
--
-- See <https://en.wikipedia.org/wiki/Earth_mover's_distance>.
compareEMD :: (Num a, Storable a)
           => Histogram DIM1 a -> Histogram DIM1 a -> a
compareEMD :: forall a.
(Num a, Storable a) =>
Histogram DIM1 a -> Histogram DIM1 a -> a
compareEMD hist1 :: Histogram DIM1 a
hist1@(Histogram DIM1
sh1 Vector a
_) hist2 :: Histogram DIM1 a
hist2@(Histogram DIM1
sh2 Vector a
_)
    | DIM1
sh1 forall a. Eq a => a -> a -> Bool
/= DIM1
sh2 = forall a. HasCallStack => String -> a
error String
"Histograms are not of equal size."
    | Bool
otherwise  = let Histogram DIM1
_ Vector a
vec1 = forall a.
(Storable a, Num a) =>
Histogram DIM1 a -> Histogram DIM1 a
cumulative Histogram DIM1 a
hist1
                       Histogram DIM1
_ Vector a
vec2 = forall a.
(Storable a, Num a) =>
Histogram DIM1 a -> Histogram DIM1 a
cumulative Histogram DIM1 a
hist2
                   in forall a. (Storable a, Num a) => Vector a -> a
V.sum forall a b. (a -> b) -> a -> b
$ forall a b c.
(Storable a, Storable b, Storable c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (\a
v1 a
v2 -> forall a. Num a => a -> a
abs (a
v1 forall a. Num a => a -> a -> a
- a
v2)) Vector a
vec1 Vector a
vec2
{-# SPECIALIZE compareEMD
    :: Histogram DIM1 Int32  -> Histogram DIM1 Int32  -> Int32
    ,  Histogram DIM1 Double -> Histogram DIM1 Double -> Double
    ,  Histogram DIM1 Float  -> Histogram DIM1 Float  -> Float #-}
{-# INLINABLE compareEMD #-}

square :: Num a => a -> a
square :: forall a. Num a => a -> a
square a
a = a
a forall a. Num a => a -> a -> a
* a
a

double :: Integral a => a -> Double
double :: forall a. Integral a => a -> Double
double= forall a b. (Integral a, Num b) => a -> b
fromIntegral

int :: Integral a => a -> Int
int :: forall a. Integral a => a -> Int
int = forall a b. (Integral a, Num b) => a -> b
fromIntegral