Parasitemia/ParasitemiaCore/ImgTools/Edges.fs

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