Parasitemia/ParasitemiaCore/ImgTools/Edges.fs

   1 module ParasitemiaCore.Edges
   2
   3 open System.Drawing
   4 open System.Collections.Generic
   5
   6 open Emgu.CV
   7 open Emgu.CV.Structure
   8
   9 open Const
  10 open Histogram
  11 open Otsu
  12
  13 // Sensibilities of the hysteresis search.
  14 let sensibilityHigh = 0.1f
  15 let sensibilityLow = 0.0f
  16
  17 /// <summary>
  18 /// Find edges of an image by using the Canny approach.
  19 /// The thresholds are automatically defined with otsu on gradient magnitudes.
  20 /// </summary>
  21 /// <param name="img"></param>
  22 let find (img : Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
  23     let w = img.Width
  24     let h = img.Height
  25
  26     use sobelKernel =
  27         new Matrix<float32> (
  28             array2D [[ -1.0f; 0.0f; 1.0f ]
  29                      [ -2.0f; 0.0f; 2.0f ]
  30                      [ -1.0f; 0.0f; 1.0f ]]
  31         )
  32
  33     let xGradient = new Matrix<float32> (img.Size)
  34     let yGradient = new Matrix<float32> (img.Size)
  35     CvInvoke.Filter2D (img, xGradient, sobelKernel, Point (1, 1))
  36     CvInvoke.Filter2D (img, yGradient, sobelKernel.Transpose (), Point (1, 1))
  37
  38     use magnitudes = new Matrix<float32> (xGradient.Size)
  39     use angles = new Matrix<float32> (xGradient.Size)
  40     CvInvoke.CartToPolar (xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. The angles are between 0 and 2 * pi.
  41
  42     let thresholdHigh, thresholdLow =
  43         let threshold, _, _ = otsu (histogramMat magnitudes 300)
  44         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
  45
  46     // Non-maximum suppression.
  47     use nms = new Matrix<byte> (xGradient.Size)
  48
  49     let nmsData = nms.Data
  50     let anglesData = angles.Data
  51     let magnitudesData = magnitudes.Data
  52     let xGradientData = xGradient.Data
  53     let yGradientData = yGradient.Data
  54
  55     for i = 1 to h - 2 do
  56         for j = 1 to w - 2 do
  57             let vx = xGradientData.[i, j]
  58             let vy = yGradientData.[i, j]
  59             if vx <> 0.f || vy <> 0.f then
  60                 let angle = anglesData.[i, j]
  61
  62                 let vx', vy' = abs vx, abs vy
  63                 let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
  64                 let ratio1 = 1.f - ratio2
  65
  66                 let mNeigbors (sign : int) : float32 =
  67                     if angle < PI / 4.f then
  68                         ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
  69                     elif angle < PI / 2.f then
  70                         ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
  71                     elif angle < 3.f * PI / 4.f then
  72                         ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
  73                     elif angle < PI then
  74                         ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
  75                     elif angle < 5.f * PI / 4.f then
  76                         ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
  77                     elif angle < 3.f * PI / 2.f then
  78                         ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
  79                     elif angle < 7.f * PI / 4.f then
  80                         ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
  81                     else
  82                         ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
  83
  84                 let m = magnitudesData.[i, j]
  85                 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 then
  86                     nmsData.[i, j] <- 1uy
  87
  88     // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
  89
  90     let edges = new Matrix<byte> (xGradient.Size)
  91     let edgesData = edges.Data
  92
  93     // Hysteresis thresholding.
  94     let toVisit = Stack<Point> ()
  95     for i = 0 to h - 1 do
  96         for j = 0 to w - 1 do
  97             if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh then
  98                 nmsData.[i, j] <- 0uy
  99                 toVisit.Push (Point (j, i))
 100                 while toVisit.Count > 0 do
 101                     let p = toVisit.Pop ()
 102                     edgesData.[p.Y, p.X] <- 1uy
 103                     for i' = -1 to 1  do
 104                         for j' = -1 to 1 do
 105                             if i' <> 0 || j' <> 0 then
 106                                 let ni = p.Y + i'
 107                                 let nj = p.X + j'
 108                                 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy then
 109                                     nmsData.[ni, nj] <- 0uy
 110                                     toVisit.Push (Point (nj, ni))
 111
 112     edges, xGradient, yGradient