module ParasitemiaCore.Edges

open System
open System.Drawing
open System.Collections.Generic

open Emgu.CV
open Emgu.CV.Structure

open Const
open Histogram
open Otsu

// Sensibilities of the hysteresis search.
let sensibilityHigh = 0.1f
let sensibilityLow = 0.0f

/// <summary>
/// Find edges of an image by using the Canny approach.
/// The thresholds are automatically defined with otsu on gradient magnitudes.
/// </summary>
/// <param name="img"></param>
let find (img : Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
    let w = img.Width
    let h = img.Height

    use sobelKernel =
        new Matrix<float32> (
            array2D [[ -1.0f; 0.0f; 1.0f ]
                     [ -2.0f; 0.0f; 2.0f ]
                     [ -1.0f; 0.0f; 1.0f ]]
        )

    let xGradient = new Matrix<float32> (img.Size)
    let yGradient = new Matrix<float32> (img.Size)
    CvInvoke.Filter2D (img, xGradient, sobelKernel, Point (1, 1))
    CvInvoke.Filter2D (img, yGradient, sobelKernel.Transpose (), Point (1, 1))

    use magnitudes = new Matrix<float32> (xGradient.Size)
    use angles = new Matrix<float32> (xGradient.Size)
    CvInvoke.CartToPolar (xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. The angles are between 0 and 2 * pi.

    let thresholdHigh, thresholdLow =
        let threshold, _, _ = otsu (histogramMat magnitudes 300)
        threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)

    // Non-maximum suppression.
    use nms = new Matrix<byte> (xGradient.Size)

    let nmsData = nms.Data
    let anglesData = angles.Data
    let magnitudesData = magnitudes.Data
    let xGradientData = xGradient.Data
    let yGradientData = yGradient.Data

    for i = 1 to h - 2 do
        for j = 1 to w - 2 do
            let vx = xGradientData.[i, j]
            let vy = yGradientData.[i, j]
            if vx <> 0.f || vy <> 0.f then
                let angle = anglesData.[i, j]

                let vx', vy' = abs vx, abs vy
                let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
                let ratio1 = 1.f - ratio2

                let mNeigbors (sign : int) : float32 =
                    if angle < PI / 4.f then
                        ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
                    elif angle < PI / 2.f then
                        ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
                    elif angle < 3.f * PI / 4.f then
                        ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
                    elif angle < PI then
                        ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
                    elif angle < 5.f * PI / 4.f then
                        ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
                    elif angle < 3.f * PI / 2.f then
                        ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
                    elif angle < 7.f * PI / 4.f then
                        ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
                    else
                        ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]

                let m = magnitudesData.[i, j]
                if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 then
                    nmsData.[i, j] <- 1uy

    // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).

    let edges = new Matrix<byte> (xGradient.Size)
    let edgesData = edges.Data

    // Hysteresis thresholding.
    let toVisit = Stack<Point> ()
    for i = 0 to h - 1 do
        for j = 0 to w - 1 do
            if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh then
                nmsData.[i, j] <- 0uy
                toVisit.Push (Point (j, i))
                while toVisit.Count > 0 do
                    let p = toVisit.Pop ()
                    edgesData.[p.Y, p.X] <- 1uy
                    for i' = -1 to 1  do
                        for j' = -1 to 1 do
                            if i' <> 0 || j' <> 0 then
                                let ni = p.Y + i'
                                let nj = p.X + j'
                                if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy then
                                    nmsData.[ni, nj] <- 0uy
                                    toVisit.Push (Point (nj, ni))

    edges, xGradient, yGradient