module ParasitemiaCore.Otsu

open Histogram

let otsu (hist: Histogram) : float32 * float32 * float32 =
    let mutable sumB = 0
    let mutable wB = 0
    let mutable maximum = 0.0
    let mutable level = 0
    let sum = hist.data |> Array.mapi (fun i v -> i * v |> float) |> Array.sum

    for i in 0 .. hist.data.Length - 1 do
        wB <- wB + hist.data.[i]
        if wB <> 0
        then
            let wF = hist.total - wB
            if wF <> 0
            then
                sumB <- sumB + i * hist.data.[i]
                let mB = (float sumB) / (float wB)
                let mF = (sum - float sumB) / (float wF)
                let between = (float wB) * (float wF) * (mB - mF) ** 2.;
                if between >= maximum
                then
                    level <- i
                    maximum <- between

    let mean1 =
        let mutable sum = 0
        let mutable nb = 0
        for i in 0 .. level - 1 do
            sum <- sum + i * hist.data.[i]
            nb <- nb + hist.data.[i]
        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)

    let mean2 =
        let mutable sum = 0
        let mutable nb = 0
        for i in level + 1 .. hist.data.Length - 1 do
            sum <- sum + i * hist.data.[i]
            nb <- nb + hist.data.[i]
        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)

    let toFloat l =
        float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min

    toFloat level, toFloat mean1, toFloat mean2