-- SPDX-License-Identifier: MPL-2.0

{- | Image measurements use caller-declared intensity ranges and explicit
local windows. All images must share shape/channel count for comparisons.
-}
module HMath.Image.Measure (
    TemplateMetric (..),
    templateMatch,
    rawMoment,
    centralMoment,
    imageCentroid,
    integralImage,
    rectangleSum,
    meanSquaredError,
    peakSignalToNoise,
    structuralSimilarity,
) where

import Data.Array (Array, array, (!))
import HMath.Error
import HMath.Image.Core
import HMath.Image.Filter (scalar2D)

{- | Template matching score: squared difference or normalized correlation. Interpret
whether smaller/larger values are better according to the metric.
-}
data TemplateMetric = SquaredDifference | NormalizedCorrelation deriving (TemplateMetric -> TemplateMetric -> Bool
(TemplateMetric -> TemplateMetric -> Bool)
-> (TemplateMetric -> TemplateMetric -> Bool) -> Eq TemplateMetric
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: TemplateMetric -> TemplateMetric -> Bool
== :: TemplateMetric -> TemplateMetric -> Bool
$c/= :: TemplateMetric -> TemplateMetric -> Bool
/= :: TemplateMetric -> TemplateMetric -> Bool
Eq, Int -> TemplateMetric -> ShowS
[TemplateMetric] -> ShowS
TemplateMetric -> String
(Int -> TemplateMetric -> ShowS)
-> (TemplateMetric -> String)
-> ([TemplateMetric] -> ShowS)
-> Show TemplateMetric
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> TemplateMetric -> ShowS
showsPrec :: Int -> TemplateMetric -> ShowS
$cshow :: TemplateMetric -> String
show :: TemplateMetric -> String
$cshowList :: [TemplateMetric] -> ShowS
showList :: [TemplateMetric] -> ShowS
Show)

{- | Valid 2D, one-channel matching. Each output is one template origin.
NCC centers both patches and divides by their energies; constant templates
or candidate patches have undefined correlation and return MathError.
-}
templateMatch :: TemplateMetric -> Image -> Image -> Either MathError Image
templateMatch :: TemplateMetric -> Image -> Image -> Either MathError Image
templateMatch TemplateMetric
method Image
needle Image
haystack = do
    (th, tw) <- Image -> Either MathError (Int, Int)
scalar2D Image
needle
    (height, width) <- scalar2D haystack
    require
        ( th <= height
            && tw <= width
            && toInteger th * toInteger tw * toInteger (height - th + 1) * toInteger (width - tw + 1) <= 20000000
        )
        "template size/work limit"
    let template = Image -> [Double]
imageValues Image
needle
        origins = [Int] -> [[Int]]
coordinates [Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
th Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
tw Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1]
        evaluate [Int]
origin =
            let patch :: [Double]
patch = [Image -> [Int] -> Int -> Double
at Image
haystack [[Int]
origin [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i, [Int]
origin [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j] Int
0 | Int
i <- [Int
0 .. Int
th Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
tw Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
             in case TemplateMetric
method of
                    TemplateMetric
SquaredDifference -> Double -> Either MathError Double
single ([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 c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
template [Double]
patch))
                    TemplateMetric
NormalizedCorrelation -> do
                        let centered :: [b] -> [b]
centered [b]
values = let mean :: b
mean = [b] -> b
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
forall a. Fractional a => a -> a -> a
/ Int -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([b] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
values)) [b]
values) in (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
forall a. Num a => a -> a -> a
subtract b
mean) [b]
values
                            a :: [Double]
a = [Double] -> [Double]
forall {b}. Fractional b => [b] -> [b]
centered [Double]
template
                            b :: [Double]
b = [Double] -> [Double]
forall {b}. Fractional b => [b] -> [b]
centered [Double]
patch
                            den :: Double
den = Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
a)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
b))
                        Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
den Bool -> Bool -> Bool
&& Double
den Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"NCC requires positive finite patch variance"
                        Double -> Either MathError Double
single ([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 c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
a [Double]
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den)
    values <- mapM evaluate origins
    image [height - th + 1, width - tw + 1] 1 values

-- | Raw M_pq=sum(row^p*column^q*intensity), orders 0..8 each.
rawMoment :: Int -> Int -> Image -> Either MathError Double
rawMoment :: Int -> Int -> Image -> Either MathError Double
rawMoment Int
p Int
q Image
img = do
    (height, width) <- Image -> Either MathError (Int, Int)
scalar2D Image
img
    require (p >= 0 && q >= 0 && p <= 8 && q <= 8) "moment orders outside [0,8]"
    single
        ( sum
            [ fromIntegral i ^ p * fromIntegral j ^ q * at img [i, j] 0
            | i <- [0 .. height - 1]
            , j <- [0 .. width - 1]
            ]
        )

{- | Intensity-weighted (row,column) centroid. Signed intensities are allowed;
a zero total mass is undefined and reported explicitly.
-}
imageCentroid :: Image -> Either MathError (Double, Double)
imageCentroid :: Image -> Either MathError (Double, Double)
imageCentroid Image
img = do
    mass <- Int -> Int -> Image -> Either MathError Double
rawMoment Int
0 Int
0 Image
img
    require (mass /= 0) "zero image mass"
    row <- rawMoment 1 0 img >>= single . (/ mass)
    col <- rawMoment 0 1 img >>= single . (/ mass)
    Right (row, col)

{- | Given row exponent p,column exponent q,then scalar 2D image, compute the unnormalized
intensity-weighted central moment around its centroid. Require orders in [0,8]; invalid
image/zero mass or non-finite result fails.
-}
centralMoment :: Int -> Int -> Image -> Either MathError Double
centralMoment :: Int -> Int -> Image -> Either MathError Double
centralMoment Int
p Int
q Image
img = do
    (height, width) <- Image -> Either MathError (Int, Int)
scalar2D Image
img
    require (p >= 0 && q >= 0 && p <= 8 && q <= 8) "moment orders outside [0,8]"
    (cy, cx) <- imageCentroid img
    single
        ( sum
            [ (fromIntegral i - cy) ^ p * (fromIntegral j - cx) ^ q * at img [i, j] 0
            | i <- [0 .. height - 1]
            , j <- [0 .. width - 1]
            ]
        )

{- | Padded 2D summed-area image, per channel: output shape (height+1,width+1)
has a leading zero row/column. Rectangle queries use half-open bounds.
-}
integralImage :: Image -> Either MathError Image
integralImage :: Image -> Either MathError Image
integralImage 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) String
"2D integral image required"
    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; nc :: Int
nc = Image -> Int
imageChannels Image
img
    Bool -> String -> Either MathError ()
require
        (Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nc Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
262144)
        String
"integral image output storage limit"
    let values :: Array (Int, Int, Int) Double
        values :: Array (Int, Int, Int) Double
values =
            ((Int, Int, Int), (Int, Int, Int))
-> [((Int, Int, Int), Double)] -> Array (Int, Int, Int) Double
forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array
                ((Int
0, Int
0, Int
0), (Int
height, Int
width, Int
nc Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
                [ ( (Int
i, Int
j, Int
c)
                  , if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                        then Double
0
                        else
                            Image -> [Int] -> Int -> Double
at Image
img [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] Int
c
                                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Array (Int, Int, Int) Double
values Array (Int, Int, Int) Double -> (Int, Int, 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
1, Int
j, Int
c)
                                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Array (Int, Int, Int) Double
values Array (Int, Int, Int) Double -> (Int, Int, Int) -> Double
forall i e. Ix i => Array i e -> i -> e
! (Int
i, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, Int
c)
                                Double -> Double -> Double
forall a. Num a => a -> a -> a
- Array (Int, Int, Int) Double
values Array (Int, Int, Int) Double -> (Int, Int, 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
1, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1, Int
c)
                  )
                | Int
i <- [Int
0 .. Int
height]
                , Int
j <- [Int
0 .. Int
width]
                , Int
c <- [Int
0 .. Int
nc Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                ]
    [Int] -> Int -> ([Int] -> Int -> Double) -> Either MathError Image
build [Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1, Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1] Int
nc (\[Int]
ix Int
c -> Array (Int, Int, Int) Double
values Array (Int, Int, Int) Double -> (Int, Int, Int) -> Double
forall i e. Ix i => Array i e -> i -> e
! ([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, Int
c))

{- | Given (row,column) origin,(height,width) extent,then padded integralImage output,
return per-channel sum of the half-open rectangle. Reject invalid/out-of-range bounds;
zero extents are allowed.
-}
rectangleSum :: (Int, Int) -> (Int, Int) -> Image -> Either MathError [Double]
rectangleSum :: (Int, Int) -> (Int, Int) -> Image -> Either MathError [Double]
rectangleSum (Int
row, Int
col) (Int
height, Int
width) Image
summed = 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
summed) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
            Bool -> Bool -> Bool
&& Int
row Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
col Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
height Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
row Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
height Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Image -> [Int]
imageShape Image
summed [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
col Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
width Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Image -> [Int]
imageShape Image
summed [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)
        )
        String
"summed-area rectangle out of bounds"
    [Double] -> Either MathError [Double]
checked
        [ Image -> [Int] -> Int -> Double
at Image
summed [Int
row Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
height, Int
col Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
width] Int
c
            Double -> Double -> Double
forall a. Num a => a -> a -> a
- Image -> [Int] -> Int -> Double
at Image
summed [Int
row, Int
col Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
width] Int
c
            Double -> Double -> Double
forall a. Num a => a -> a -> a
- Image -> [Int] -> Int -> Double
at Image
summed [Int
row Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
height, Int
col] Int
c
            Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Image -> [Int] -> Int -> Double
at Image
summed [Int
row, Int
col] Int
c
        | Int
c <- [Int
0 .. Image -> Int
imageChannels Image
summed Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        ]

same :: Image -> Image -> Either MathError ()
same :: Image -> Image -> Either MathError ()
same Image
a Image
b =
    Bool -> String -> Either MathError ()
require
        (Image -> [Int]
imageShape Image
a [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== Image -> [Int]
imageShape Image
b Bool -> Bool -> Bool
&& Image -> Int
imageChannels Image
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Image -> Int
imageChannels Image
b)
        String
"comparison shape mismatch"

{- | Given two images of equal spatial shape and channels, return average squared sample
difference. Shape mismatch or non-finite arithmetic returns Left MathError.
-}
meanSquaredError :: Image -> Image -> Either MathError Double
meanSquaredError :: Image -> Image -> Either MathError Double
meanSquaredError Image
a Image
b = do
    Image -> Image -> Either MathError ()
same Image
a Image
b
    let n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Double]
imageValues Image
a))
    Double -> Either MathError Double
single ([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 c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> ((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) (Image -> [Double]
imageValues Image
a) (Image -> [Double]
imageValues Image
b)))

{- | PSNR in decibels. Exactly identical images return positive infinity by
definition; this is the sole intentional nonfinite measurement result.
-}
peakSignalToNoise :: Double -> Image -> Image -> Either MathError Double
peakSignalToNoise :: Double -> Image -> Image -> Either MathError Double
peakSignalToNoise Double
dataRange Image
a Image
b = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
dataRange Bool -> Bool -> Bool
&& Double
dataRange Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"positive PSNR range required"
    mse <- Image -> Image -> Either MathError Double
meanSquaredError Image
a Image
b
    if mse == 0 then Right (1 / 0) else single (20 * logBase 10 dataRange - 10 * logBase 10 mse)

{- | Mean local SSIM with uniform valid windows and population covariance.
C1=(0.01*range)^2, C2=(0.03*range)^2; all channels are measured separately.
Windows do not include padding and both spatial extents must fit the window.
Wang et al., Image Quality Assessment: From Error Visibility to Structural
Similarity (2004), https://ece.uwaterloo.ca/~z70wang/publications/ssim.html .
-}
structuralSimilarity :: Double -> Int -> Image -> Image -> Either MathError Double
structuralSimilarity :: Double -> Int -> Image -> Image -> Either MathError Double
structuralSimilarity Double
dataRange Int
window Image
a Image
b = do
    Image -> Image -> Either MathError ()
same Image
a Image
b
    Bool -> String -> Either MathError ()
require
        ([Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Image -> [Int]
imageShape Image
a) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
dataRange Bool -> Bool -> Bool
&& Double
dataRange Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Int
window Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int
window Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
63)
        String
"SSIM range/window/rank"
    let height :: Int
height = Image -> [Int]
imageShape Image
a [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0; width :: Int
width = Image -> [Int]
imageShape Image
a [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
1; nc :: Int
nc = Image -> Int
imageChannels Image
a
    Bool -> String -> Either MathError ()
require
        ( Int
window Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
height
            Bool -> Bool -> Bool
&& Int
window Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
width
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
window
                Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
window
                Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
window Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
window Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nc
                Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
20000000
        )
        String
"SSIM window/work limit"
    let c1 :: Double
c1 = (Double
0.01 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dataRange) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
        c2 :: Double
c2 = (Double
0.03 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dataRange) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
        score :: Int -> Int -> Int -> Double
score Int
i Int
j Int
c =
            let patch :: Image -> [Double]
patch Image
img = [Image -> [Int] -> Int -> Double
at Image
img [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
c | Int
di <- [Int
0 .. Int
window Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
dj <- [Int
0 .. Int
window Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                x :: [Double]
x = Image -> [Double]
patch Image
a
                y :: [Double]
y = Image -> [Double]
patch Image
b
                n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x)
                mx :: Double
mx = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n) [Double]
x)
                my :: Double
my = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n) [Double]
y)
                vx :: Double
vx = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [((Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mx) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
z <- [Double]
x]
                vy :: Double
vy = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [((Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
my) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
z <- [Double]
y]
                covariance :: Double
covariance = [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 c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
u Double
v -> ((Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mx) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n) Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
my) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
n)) [Double]
x [Double]
y)
             in ((Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
mx Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
my Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
covariance Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c2)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ((Double
mx Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
mx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
my Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
my Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
vx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
vy Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c2))
        scores :: [Double]
scores = [Int -> Int -> Int -> Double
score Int
i Int
j Int
c | Int
i <- [Int
0 .. Int
height Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
window], Int
j <- [Int
0 .. Int
width Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
window], Int
c <- [Int
0 .. Int
nc Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    _ <- [Double] -> Either MathError [Double]
checked (Double
c1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Double
c2 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
scores)
    single (sum (map (/ fromIntegral (length scores)) scores))

single :: Double -> Either MathError Double
single :: Double -> Either MathError Double
single Double
value = do _ <- [Double] -> Either MathError [Double]
checked [Double
value]; Right value