Parasitemia/Parasitemia/ImgTools.fs

   1 module ImgTools
   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 Utils
  11
  12 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
  13     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
  14     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
  15
  16 // Zhang and Suen algorithm.
  17 // Modify 'mat' in place.
  18 let thin (mat: Matrix<byte>) =
  19     let neighbors = [|
  20         (-1,  0) // p2
  21         (-1,  1) // p3
  22         ( 0,  1) // p4
  23         ( 1,  1) // p5
  24         ( 1,  0) // p6
  25         ( 1, -1) // p7
  26         ( 0, -1) // p8
  27         (-1, -1) |] // p9
  28
  29     let w = mat.Width
  30     let h = mat.Height
  31     let mutable data1 = mat.Data
  32     let mutable data2 = Array2D.zeroCreate<byte> h w
  33
  34     // Return the list of neighbor values.
  35     let neighborsValues (p1i, p1j) =
  36         Array.map (fun (ni, nj) ->
  37             let pi = p1i + ni
  38             let pj = p1j + nj
  39             if pi < 0 || pi >= h || pj < 0 || pj >= w then 0uy else data1.[pi, pj]
  40             ) neighbors
  41
  42     // Return the number of 01 pattern in 'values' in a circular way.
  43     let pattern01 (values: byte[]) =
  44         let mutable nb = 0
  45         let mutable lastValue = 255uy
  46         for v in values do
  47             if lastValue = 0uy && v = 1uy
  48             then
  49                 nb <- nb + 1
  50             lastValue <- v
  51         if lastValue = 0uy && values.[0] = 1uy
  52         then
  53             nb <- nb + 1
  54         nb
  55
  56     let mutable pixelChanged = true
  57     let mutable oddIteration = true
  58     while pixelChanged do
  59         pixelChanged <- false
  60         for i in 0..h-1 do
  61             for j in 0..w-1 do
  62                 if data1.[i, j] = 1uy
  63                 then
  64                     let values = neighborsValues (i, j)
  65                     let s = Array.reduce (+) values
  66                     if s >= 2uy && s <= 6uy &&
  67                         pattern01 values = 1 &&
  68                         (not oddIteration || (values.[0] * values.[2] * values.[4] = 0uy && values.[2] * values.[4] * values.[6] = 0uy)) && // Odd iteration.
  69                         (oddIteration || (values.[0] * values.[2] * values.[6] = 0uy && values.[0] * values.[4] * values.[6] = 0uy)) // Even iterations.
  70                     then
  71                         data2.[i, j] <- 0uy
  72                         pixelChanged <- true
  73                     else
  74                         data2.[i, j] <- 1uy
  75                 else
  76                     data2.[i, j] <- 0uy
  77
  78         oddIteration <- not oddIteration
  79         let tmp = data1
  80         data1 <- data2
  81         data2 <- tmp
  82
  83
  84 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
  85 // Modify 'mat' in place.
  86 let removeArea (mat: Matrix<byte>) (areaSize: int) =
  87     let neighbors = [|
  88         (-1,  0) // p2
  89         (-1,  1) // p3
  90         ( 0,  1) // p4
  91         ( 1,  1) // p5
  92         ( 1,  0) // p6
  93         ( 1, -1) // p7
  94         ( 0, -1) // p8
  95         (-1, -1) |] // p9
  96
  97     let mat' = new Matrix<byte>(mat.Size)
  98     let w = mat'.Width
  99     let h = mat'.Height
 100     mat.CopyTo(mat')
 101
 102     let data = mat.Data
 103     let data' = mat'.Data
 104
 105     for i in 0..h-1 do
 106         for j in 0..w-1 do
 107             if data'.[i, j] = 1uy
 108             then
 109                 let neighborhood = List<(int*int)>()
 110                 let neighborsToCheck = List<(int*int)>()
 111                 neighborsToCheck.Add((i, j))
 112                 data'.[i, j] <- 0uy
 113
 114                 let pop (l: List<'a>) : 'a =
 115                     let n = l.[l.Count - 1]
 116                     l.RemoveAt(l.Count - 1)
 117                     n
 118
 119                 while neighborsToCheck.Count > 0 do
 120                     let (ci, cj) = pop neighborsToCheck
 121                     neighborhood.Add((ci, cj))
 122                     for (ni, nj) in neighbors do
 123                         let pi = ci + ni
 124                         let pj = cj + nj
 125                         if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
 126                         then
 127                             neighborsToCheck.Add((pi, pj))
 128                             data'.[pi, pj] <- 0uy
 129                 if neighborhood.Count <= areaSize
 130                 then
 131                     for (ni, nj) in neighborhood do
 132                         data.[ni, nj] <- 0uy
 133
 134
 135 let saveImg (img: Image<'TColor, 'TDepth>) (name: string) =
 136     img.Save("output/" + name)
 137
 138
 139 let saveMat (mat: Matrix<'TDepth>) (name: string) =
 140     use img = new Image<Gray, 'TDeph>(mat.Size)
 141     mat.CopyTo(img)
 142     saveImg img name
 143
 144 (*let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
 145     let e' = Ellipse(PointF(float32 e.cx, float32 e.cy), SizeF(2.0f * float32 e.a, 2.0f * float32 e.b), float32 e.alpha)
 146     img.Draw(e', color)*)
 147
 148 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) =
 149     let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1)
 150
 151     img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1);
 152
 153 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
 154     let cosAlpha = cos e.alpha
 155     let sinAlpha = sin e.alpha
 156
 157     let mutable x0 = 0.0
 158     let mutable y0 = 0.0
 159     let mutable first_iteration = true
 160
 161     let n = 40
 162     let thetaIncrement = 2.0 * Math.PI / (float n)
 163
 164     for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do
 165         let cosTheta = cos theta
 166         let sinTheta = sin theta
 167         let x = e.cx + cosAlpha * e.a * cosTheta - sinAlpha * e.b * sinTheta
 168         let y = e.cy + sinAlpha * e.a * cosTheta + cosAlpha * e.b * sinTheta
 169
 170         if not first_iteration
 171         then
 172             drawLine img color x0 y0 x y
 173         else
 174             first_iteration <- false
 175
 176         x0 <- x
 177         y0 <- y