Parasitemia/ParasitemiaCore/Granulometry.fs

   1 module ParasitemiaCore.Granulometry
   2
   3 open System
   4 open System.Drawing
   5
   6 open Emgu.CV
   7 open Emgu.CV.Structure
   8
   9 open Utils
  10
  11 /// <summary>
  12 /// Granulometry by closing on the image 'img' by testing circle structrual elements with radius range in 'range'.
  13 /// </summary>
  14 /// <param name="img"></param>
  15 /// <param name="range">Minimum radius * maximum radius</param>
  16 /// <param name="scale">le 1.0, to speed up the process.</param>
  17 /// <param name="useOctagon"></param>
  18 let findRadiusByClosing (img : Image<Gray, 'TDepth>) (range : int * int) (scale : float) (useOctagon : bool) : int =
  19     use scaledImg = if scale = 1. then img else img.Resize (scale, CvEnum.Inter.Area)
  20
  21     let r1, r2 = range
  22     let r1', r2' = roundInt (float r1 * scale), roundInt (float r2 * scale)
  23
  24     let patternSpectrum = Array.zeroCreate (r2' - r1')
  25     let intensityImg = scaledImg.GetSum().Intensity
  26
  27     // 's' must be odd.
  28     let octagon (s : int) : Matrix<byte> =
  29         if s % 2 = 0 then failwith "s must be odd"
  30         let m = new Matrix<byte> (Array2D.create s s 1uy)
  31         let r = (float s) / (Math.Sqrt 2. + 2.) |> roundInt
  32         for i = 0 to r - 1 do
  33             for j = 0 to r - 1 do
  34                 if i + j < r then
  35                     m.[i, j] <- 0uy
  36                     m.[s - i - 1, j] <- 0uy
  37                     m.[i, s - j - 1] <- 0uy
  38                     m.[s - i - 1, s - j - 1] <- 0uy
  39         m
  40
  41     let mutable previous_n = Double.NaN
  42     for r = r1' to r2' do
  43         let se =
  44             if useOctagon then
  45                 (octagon (2 * r - 1)).Mat // It doesn't speed up the process.
  46             else
  47                 CvInvoke.GetStructuringElement (CvEnum.ElementShape.Ellipse, Size (2 * r, 2 * r), Point (-1, -1))
  48
  49         use closed = scaledImg.MorphologyEx (CvEnum.MorphOp.Close, se, Point (-1, -1), 1, CvEnum.BorderType.Replicate, MCvScalar (0.0))
  50
  51         let n = closed.GetSum().Intensity
  52
  53         if not (Double.IsNaN previous_n) then
  54             patternSpectrum.[r - r1' - 1] <- abs (n - previous_n)
  55         previous_n <- n
  56
  57     let max, _ = patternSpectrum |> Array.indexed |> Array.fold (fun (iMax, sMax) (i, s) -> if s > sMax then (i, s) else (iMax, sMax)) (0, Double.MinValue)
  58
  59     float (max + r1') / scale |> roundInt
  60
  61 /// <summary>
  62 /// Granulometry by area closing on the image 'img' by testing the circle area into the given radius range.
  63 /// </summary>
  64 let findRadiusByAreaClosing (img : Image<Gray, float32>) (radiusRange : int * int) : int =
  65     let r1, r2 = radiusRange
  66
  67     if r1 > r2 then
  68         failwithf "'radiusRange' invalid : %O" radiusRange
  69
  70     use imgCopy = img.Copy ()
  71
  72     let mutable maxDiff = 0.f
  73     let mutable max_r = r1
  74
  75     Morpho.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (
  76         fun r diff ->
  77             if r <> r1 && diff > maxDiff then
  78                 maxDiff <- diff
  79                 max_r <- r - 1
  80     )
  81
  82     max_r
  83