module ParasitemiaCore.Granulometry

open System
open System.IO
open System.Drawing

open Emgu.CV
open Emgu.CV.Structure

open Utils

/// <summary>
/// Granulometry by closing on the image 'img' by testing circle structrual elements with radius range in 'range'.
/// </summary>
/// <param name="img"></param>
/// <param name="range">Minimum radius * maximum radius</param>
/// <param name="scale">le 1.0, to speed up the process.</param>
let findRadiusByClosing (img : Image<Gray, 'TDepth>) (range : int * int) (scale : float) (useOctagon : bool) : int =
    use scaledImg = if scale = 1. then img else img.Resize (scale, CvEnum.Inter.Area)

    let r1, r2 = range
    let r1', r2' = roundInt (float r1 * scale), roundInt (float r2 * scale)

    let patternSpectrum = Array.zeroCreate (r2' - r1')
    let intensityImg = scaledImg.GetSum().Intensity

    // 's' must be odd.
    let octagon (s : int) : Matrix<byte> =
        if s % 2 = 0 then failwith "s must be odd"
        let m = new Matrix<byte> (Array2D.create s s 1uy)
        let r = (float s) / (Math.Sqrt 2. + 2.) |> roundInt
        for i = 0 to r - 1 do
            for j = 0 to r - 1 do
                if i + j < r then
                    m.[i, j] <- 0uy
                    m.[s - i - 1, j] <- 0uy
                    m.[i, s - j - 1] <- 0uy
                    m.[s - i - 1, s - j - 1] <- 0uy
        m

    let mutable previous_n = Double.NaN
    for r = r1' to r2' do
        let se =
            if useOctagon then
                (octagon (2 * r - 1)).Mat // It doesn't speed up the process.
            else
                CvInvoke.GetStructuringElement (CvEnum.ElementShape.Ellipse, Size (2 * r, 2 * r), Point (-1, -1))

        use closed = scaledImg.MorphologyEx (CvEnum.MorphOp.Close, se, Point (-1, -1), 1, CvEnum.BorderType.Replicate, MCvScalar (0.0))

        let n = closed.GetSum().Intensity

        if not (Double.IsNaN previous_n) then
            patternSpectrum.[r - r1' - 1] <- abs (n - previous_n)
        previous_n <- n

    let max, _ = patternSpectrum |> Array.indexed |> Array.fold (fun (iMax, sMax) (i, s) -> if s > sMax then (i, s) else (iMax, sMax)) (0, Double.MinValue)

    float (max + r1') / scale |> roundInt

/// <summary>
/// Granulometry by area closing on the image 'img' by testing the circle area into the given radius range.
/// </summary>
let findRadiusByAreaClosing (img : Image<Gray, float32>) (radiusRange : int * int) : int =
    let r1, r2 = radiusRange

    if r1 > r2 then
        failwithf "'radiusRange' invalid : %O" radiusRange

    use imgCopy = img.Copy ()

    let mutable maxDiff = 0.f
    let mutable max_r = r1

    Morpho.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (
        fun r diff ->
            if r <> r1 && diff > maxDiff then
                maxDiff <- diff
                max_r <- r - 1
    )

    max_r