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 // Normalize image values between 0uy and 255uy.
  13 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
  14     let min = ref [| 0.0 |]
  15     let minLocation = ref <| [| Point() |]
  16     let max = ref [| 0.0 |]
  17     let maxLocation = ref <| [| Point() |]
  18     img.MinMax(min, max, minLocation, maxLocation)
  19     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
  20
  21 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
  22     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
  23     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
  24
  25 // Zhang and Suen algorithm.
  26 // Modify 'mat' in place.
  27 let thin (mat: Matrix<byte>) =
  28     let neighbors = [|
  29         (-1,  0) // p2
  30         (-1,  1) // p3
  31         ( 0,  1) // p4
  32         ( 1,  1) // p5
  33         ( 1,  0) // p6
  34         ( 1, -1) // p7
  35         ( 0, -1) // p8
  36         (-1, -1) |] // p9
  37
  38     let w = mat.Width
  39     let h = mat.Height
  40     let mutable data1 = mat.Data
  41     let mutable data2 = Array2D.zeroCreate<byte> h w
  42
  43     // Return the list of neighbor values.
  44     let neighborsValues (p1i, p1j) =
  45         Array.map (fun (ni, nj) ->
  46             let pi = p1i + ni
  47             let pj = p1j + nj
  48             if pi < 0 || pi >= h || pj < 0 || pj >= w then 0uy else data1.[pi, pj]
  49             ) neighbors
  50
  51     // Return the number of 01 pattern in 'values' in a circular way.
  52     let pattern01 (values: byte[]) =
  53         let mutable nb = 0
  54         let mutable lastValue = 255uy
  55         for v in values do
  56             if lastValue = 0uy && v = 1uy
  57             then
  58                 nb <- nb + 1
  59             lastValue <- v
  60         if lastValue = 0uy && values.[0] = 1uy
  61         then
  62             nb <- nb + 1
  63         nb
  64
  65     let mutable pixelChanged = true
  66     let mutable oddIteration = true
  67     while pixelChanged do
  68         pixelChanged <- false
  69         for i in 0..h-1 do
  70             for j in 0..w-1 do
  71                 if data1.[i, j] = 1uy
  72                 then
  73                     let values = neighborsValues (i, j)
  74                     let s = Array.reduce (+) values
  75                     if s >= 2uy && s <= 6uy &&
  76                         pattern01 values = 1 &&
  77                         (not oddIteration || (values.[0] * values.[2] * values.[4] = 0uy && values.[2] * values.[4] * values.[6] = 0uy)) && // Odd iteration.
  78                         (oddIteration || (values.[0] * values.[2] * values.[6] = 0uy && values.[0] * values.[4] * values.[6] = 0uy)) // Even iterations.
  79                     then
  80                         data2.[i, j] <- 0uy
  81                         pixelChanged <- true
  82                     else
  83                         data2.[i, j] <- 1uy
  84                 else
  85                     data2.[i, j] <- 0uy
  86
  87         oddIteration <- not oddIteration
  88         let tmp = data1
  89         data1 <- data2
  90         data2 <- tmp
  91
  92
  93
  94 let pop (l: List<'a>) : 'a =
  95     let n = l.[l.Count - 1]
  96     l.RemoveAt(l.Count - 1)
  97     n
  98
  99 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
 100 // Modify 'mat' in place.
 101 let removeArea (mat: Matrix<byte>) (areaSize: int) =
 102     let neighbors = [|
 103         (-1,  0) // p2
 104         (-1,  1) // p3
 105         ( 0,  1) // p4
 106         ( 1,  1) // p5
 107         ( 1,  0) // p6
 108         ( 1, -1) // p7
 109         ( 0, -1) // p8
 110         (-1, -1) |] // p9
 111
 112     let mat' = new Matrix<byte>(mat.Size)
 113     let w = mat'.Width
 114     let h = mat'.Height
 115     mat.CopyTo(mat')
 116
 117     let data = mat.Data
 118     let data' = mat'.Data
 119
 120     for i in 0..h-1 do
 121         for j in 0..w-1 do
 122             if data'.[i, j] = 1uy
 123             then
 124                 let neighborhood = List<(int*int)>()
 125                 let neighborsToCheck = List<(int*int)>()
 126                 neighborsToCheck.Add((i, j))
 127                 data'.[i, j] <- 0uy
 128
 129                 while neighborsToCheck.Count > 0 do
 130                     let (ci, cj) = pop neighborsToCheck
 131                     neighborhood.Add((ci, cj))
 132                     for (ni, nj) in neighbors do
 133                         let pi = ci + ni
 134                         let pj = cj + nj
 135                         if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
 136                         then
 137                             neighborsToCheck.Add((pi, pj))
 138                             data'.[pi, pj] <- 0uy
 139                 if neighborhood.Count <= areaSize
 140                 then
 141                     for (ni, nj) in neighborhood do
 142                         data.[ni, nj] <- 0uy
 143
 144 let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
 145     let w = img.Width
 146     let h = img.Height
 147
 148     let pointChecked = HashSet<Point>()
 149     let pointToCheck = List<Point>(startPoints);
 150
 151     let data = img.Data
 152
 153     while pointToCheck.Count > 0 do
 154         let next = pop pointToCheck
 155         pointChecked.Add(next) |> ignore
 156         for ny in -1 .. 1 do
 157             for nx in -1 .. 1 do
 158                 if ny <> 0 && nx <> 0
 159                 then
 160                     let p = Point(next.X + nx, next.Y + ny)
 161                     if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
 162                     then
 163                         pointToCheck.Add(p)
 164
 165     List<Point>(pointChecked)
 166
 167
 168 let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
 169     img.Save(filepath)
 170
 171
 172 let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
 173     use img = new Image<Gray, 'TDeph>(mat.Size)
 174     mat.CopyTo(img)
 175     saveImg img filepath
 176
 177 (*let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
 178     let e' = Ellipse(PointF(float32 e.cx, float32 e.cy), SizeF(2.0f * float32 e.a, 2.0f * float32 e.b), float32 e.alpha)
 179     img.Draw(e', color)*)
 180
 181 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) =
 182     img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1);
 183
 184 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) =
 185     let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1)
 186     drawLine img color x0 y0 x1 y1
 187
 188 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
 189     let cosAlpha = cos e.Alpha
 190     let sinAlpha = sin e.Alpha
 191
 192     let mutable x0 = 0.0
 193     let mutable y0 = 0.0
 194     let mutable first_iteration = true
 195
 196     let n = 40
 197     let thetaIncrement = 2.0 * Math.PI / (float n)
 198
 199     for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do
 200         let cosTheta = cos theta
 201         let sinTheta = sin theta
 202         let x = e.Cx + cosAlpha * e.A * cosTheta - sinAlpha * e.B * sinTheta
 203         let y = e.Cy + sinAlpha * e.A * cosTheta + cosAlpha * e.B * sinTheta
 204
 205         if not first_iteration
 206         then
 207             drawLineF img color x0 y0 x y
 208         else
 209             first_iteration <- false
 210
 211         x0 <- x
 212         y0 <- y
 213
 214 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) =
 215     List.iter (fun e -> drawEllipse img e color) ellipses
 216
 217
 218 let rngCell =  System.Random()
 219 let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
 220     if drawCellContent
 221     then
 222         let colorB = rngCell.Next(20, 70)
 223         let colorG = rngCell.Next(20, 70)
 224         let colorR = rngCell.Next(20, 70)
 225
 226         for y in 0 .. c.elements.Height - 1 do
 227             for x in 0 .. c.elements.Width - 1 do
 228                 if c.elements.[y, x] > 0uy
 229                 then
 230                     let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
 231                     let b = img.Data.[y + dy, x + dx, 0] |> int
 232                     let g = img.Data.[y + dy, x + dx, 1] |> int
 233                     let r = img.Data.[y + dy, x + dx, 2] |> int
 234                     img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
 235                     img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
 236                     img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
 237
 238     let crossColor = match c.cellClass with
 239                      | Types.HealthyRBC -> Bgr(255.0, 0.0, 0.0)
 240                      | Types.InfectedRBC -> Bgr(0.0, 0.0, 255.0)
 241                      | Types.Peculiar -> Bgr(0.0, 0.0, 0.0)
 242
 243     drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y
 244     drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3)
 245
 246 let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
 247     List.iter (fun c -> drawCell img drawCellContent c) cells