Parasitemia/Parasitemia/Granulometry.fs

   1 module 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 // 'range': a minimum and maximum radius.
  12 // 'scale': <= 1.0, to speed up the process.
  13 let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: float) : int =
  14     use scaledImg = if scale = 1. then img else img.Resize(scale, CvEnum.Inter.Area)
  15
  16     let r1, r2 = range
  17     let r1', r2' = roundInt (float r1 * scale), roundInt (float r2 * scale)
  18
  19     let patternSpectrum = Array.zeroCreate (r2' - r1')
  20     let intensityImg = scaledImg.GetSum().Intensity
  21
  22     // 's' must be odd.
  23     let octagon (s: int) : Matrix<byte> =
  24         if s % 2 = 0 then failwith "s must be odd"
  25         let m = new Matrix<byte>(Array2D.create s s 1uy)
  26         let r = (float s) / (Math.Sqrt 2. + 2.) |> roundInt
  27         for i in 0 .. r - 1 do
  28             for j in 0 .. r - 1 do
  29                 if i + j < r
  30                 then
  31                     m.[i, j] <- 0uy
  32                     m.[s - i - 1, j] <- 0uy
  33                     m.[i, s - j - 1] <- 0uy
  34                     m.[s - i - 1, s - j - 1] <- 0uy
  35         m
  36
  37     let mutable previous_n = Double.NaN
  38     for r in r1' .. r2' do
  39         let se = CvInvoke.GetStructuringElement(CvEnum.ElementShape.Ellipse, Size(2 * r, 2 * r), Point(-1, -1))
  40         //let se = octagon (2 * r - 1)
  41
  42         use closed = scaledImg.MorphologyEx(CvEnum.MorphOp.Close, se, Point(-1, -1), 1, CvEnum.BorderType.Replicate, MCvScalar(0.0))
  43
  44         let n = closed.GetSum().Intensity
  45
  46         if not (Double.IsNaN previous_n)
  47         then
  48             patternSpectrum.[r - r1' - 1] <- abs (n - previous_n)
  49         previous_n <- n
  50
  51     let max, _ = patternSpectrum |> Array.indexed |> Array.fold (fun (iMax, sMax) (i, s) -> if s > sMax then (i, s) else (iMax, sMax)) (0, Double.MinValue)
  52
  53     float (max + r1') / scale |> roundInt
  54
  55
  56 let findRadiusByAreaClosing (img: Image<Gray, float32>) (range: int * int) : int =
  57     let r1, r2 = range
  58
  59     use imgCopy = img.Copy()
  60
  61     let mutable maxDiff = 0.f
  62     let mutable max_r = r1
  63
  64     ImgTools.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (fun r diff ->
  65         if r <> r1 && diff > maxDiff
  66         then
  67             maxDiff <- diff
  68             max_r <- r - 1 )
  69     max_r
  70