Parasitemia/Parasitemia/ImgTools.fs

   1 module ImgTools
   2
   3 open System
   4 open System.Drawing
   5 open System.Collections.Generic
   6 open System.Linq
   7
   8 open Emgu.CV
   9 open Emgu.CV.Structure
  10
  11 open Utils
  12 open Heap
  13
  14 // Normalize image values between 0uy and 255uy.
  15 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
  16     let min = ref [| 0.0 |]
  17     let minLocation = ref <| [| Point() |]
  18     let max = ref [| 0.0 |]
  19     let maxLocation = ref <| [| Point() |]
  20     img.MinMax(min, max, minLocation, maxLocation)
  21     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
  22
  23
  24 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
  25     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
  26     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
  27
  28
  29 let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
  30     use suppress = new Image<Gray, byte>(img.Size)
  31     let w = img.Width
  32     let h = img.Height
  33
  34     let imgData = img.Data
  35     let suppressData = suppress.Data
  36
  37     let result = List<List<Point>>()
  38
  39     let flood (start: Point) : List<List<Point>> =
  40         let sameLevelToCheck = Stack<Point>()
  41         let betterLevelToCheck = Stack<Point>()
  42         betterLevelToCheck.Push(start)
  43
  44         let result' = List<List<Point>>()
  45
  46         while betterLevelToCheck.Count > 0 do
  47             let p = betterLevelToCheck.Pop()
  48             if suppressData.[p.Y, p.X, 0] = 0uy
  49             then
  50                 suppressData.[p.Y, p.X, 0] <- 1uy
  51                 sameLevelToCheck.Push(p)
  52                 let current = List<Point>()
  53
  54                 let mutable betterExists = false
  55
  56                 while sameLevelToCheck.Count > 0 do
  57                     let p' = sameLevelToCheck.Pop()
  58                     let currentLevel = imgData.[p'.Y, p'.X, 0]
  59                     current.Add(p') |> ignore
  60                     for i in -1 .. 1 do
  61                         for j in -1 .. 1 do
  62                             if i <> 0 || j <> 0
  63                             then
  64                                 let ni = i + p'.Y
  65                                 let nj = j + p'.X
  66                                 if ni >= 0 && ni < h && nj >= 0 && nj < w
  67                                 then
  68                                     let level = imgData.[ni, nj, 0]
  69                                     let notSuppressed = suppressData.[ni, nj, 0] = 0uy
  70
  71                                     if level = currentLevel && notSuppressed
  72                                     then
  73                                         suppressData.[ni, nj, 0] <- 1uy
  74                                         sameLevelToCheck.Push(Point(nj, ni))
  75                                     elif level > currentLevel
  76                                     then
  77                                         betterExists <- true
  78                                         if notSuppressed
  79                                         then
  80                                             betterLevelToCheck.Push(Point(nj, ni))
  81
  82                 if not betterExists
  83                 then
  84                     result'.Add(current)
  85         result'
  86
  87     for i in 0 .. h - 1 do
  88         for j in 0 .. w - 1 do
  89             let maxima = flood (Point(j, i))
  90             if maxima.Count > 0
  91             then result.AddRange(maxima)
  92
  93     result.Select(fun l -> HashSet<Point>(l))
  94
  95
  96 type PriorityQueue () =
  97     let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
  98     let mutable highest = -1 // Value of the first elements of 'q'.
  99
 100     member this.Next () : byte * Point =
 101         if this.IsEmpty
 102         then
 103             invalidOp "Queue is empty"
 104         else
 105             let l = q.[0]
 106             let next = l.First()
 107             l.Remove(next) |> ignore
 108             let value = byte highest
 109             if l.Count = 0
 110             then
 111                 q.RemoveAt(0)
 112                 highest <- highest - 1
 113                 while q.Count > 0 && q.[0] = null do
 114                     q.RemoveAt(0)
 115                     highest <- highest - 1
 116             value, next
 117
 118     member this.Max =
 119         highest |> byte
 120
 121     (*member this.UnionWith (other: PriorityQueue) =
 122         while not other.IsEmpty do
 123             let p, v = other.Next
 124             this.Add p v*)
 125
 126     member this.Add (value: byte) (p: Point) =
 127         let vi = int value
 128
 129         if this.IsEmpty
 130         then
 131             highest <- int value
 132             q.Insert(0, null)
 133         elif vi > highest
 134         then
 135             for i in highest .. vi - 1  do
 136                 q.Insert(0, null)
 137             highest <- vi
 138         elif highest - vi >= q.Count
 139         then
 140             for i in 0 .. highest - vi - q.Count do
 141                 q.Add(null)
 142
 143         let pos = highest - vi
 144         if q.[pos] = null
 145         then
 146             q.[pos] <- HashSet<Point>([p])
 147         else
 148             q.[pos].Add(p) |> ignore
 149
 150
 151     member this.IsEmpty =
 152         q.Count = 0
 153
 154     member this.Clear () =
 155         while highest >= 0 do
 156             q.[highest] <- null
 157             highest <- highest - 1
 158
 159
 160
 161 type MaximaState =  Uncertain | Validated | TooBig
 162 type Maxima = {
 163     elements : HashSet<Point>
 164     mutable intensity: byte option
 165     mutable state: MaximaState }
 166
 167
 168 let areaOpen (img: Image<Gray, byte>) (area: int) =
 169     let w = img.Width
 170     let h = img.Height
 171
 172     let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
 173     let toValidated = Stack<Maxima>(maxima)
 174
 175     while toValidated.Count > 0 do
 176         let m = toValidated.Pop()
 177         if m.elements.Count <= area
 178         then
 179             let queue =
 180                 let q = PriorityQueue()
 181                 let firstElements = HashSet<Point>()
 182                 for p in m.elements do
 183                     for i in -1 .. 1 do
 184                         for j in -1 .. 1 do
 185                             if i <> 0 || j <> 0
 186                             then
 187                                 let ni = i + p.Y
 188                                 let nj = j + p.X
 189                                 let p' = Point(nj, ni)
 190                                 if ni >= 0 && ni < h && nj >= 0 && nj < w && not (m.elements.Contains(p')) && not (firstElements.Contains(p'))
 191                                 then
 192                                     firstElements.Add(p') |> ignore
 193                                     q.Add (img.Data.[ni, nj, 0]) p'
 194                 q
 195
 196             let mutable intensity = queue.Max
 197             let nextElements = HashSet<Point>()
 198
 199             let mutable stop = false
 200             while not stop do
 201                 let intensity', p = queue.Next ()
 202
 203                 if intensity' = intensity // The intensity doesn't change.
 204                 then
 205                     if m.elements.Count + nextElements.Count + 1 > area
 206                     then
 207                         m.state <- Validated
 208                         m.intensity <- Some intensity
 209                         stop <- true
 210                     else
 211                         nextElements.Add(p) |> ignore
 212                 elif intensity' < intensity
 213                 then
 214                     m.elements.UnionWith(nextElements)
 215                     if m.elements.Count = area
 216                     then
 217                         m.state <- Validated
 218                         m.intensity <- Some (intensity')
 219                         stop <- true
 220                     else
 221                         intensity <- intensity'
 222                         nextElements.Clear()
 223                         nextElements.Add(p) |> ignore
 224                 else // i' > i
 225                     seq {
 226                         for m' in maxima do
 227                             if m' <> m && m'.elements.Contains(p) then
 228                                 if m'.elements.Count + m.elements.Count <= area
 229                                 then
 230                                     m'.state <- Uncertain
 231                                     m'.elements.UnionWith(m.elements)
 232                                     if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list.
 233                                     then
 234                                         toValidated.Push(m')
 235                                     stop <- true
 236                                 yield false
 237                     } |> Seq.forall id |> ignore
 238
 239                     if not stop
 240                     then
 241                         m.state <- Validated
 242                         m.intensity <- Some (intensity)
 243                         stop <- true
 244
 245                 if not stop
 246                 then
 247                     for i in -1 .. 1 do
 248                         for j in -1 .. 1 do
 249                             if i <> 0 || j <> 0
 250                             then
 251                                 let ni = i + p.Y
 252                                 let nj = j + p.X
 253                                 let p' = Point(nj, ni)
 254                                 if ni < 0 || ni >= h || nj < 0 || nj >= w
 255                                 then
 256                                     m.state <- Validated
 257                                     m.intensity <- Some (intensity)
 258                                     stop <- true
 259                                 elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
 260                                 then
 261                                     queue.Add (img.Data.[ni, nj, 0]) p'
 262
 263                 if queue.IsEmpty
 264                 then
 265                     if m.elements.Count + nextElements.Count <= area
 266                     then
 267                         m.state <- Validated
 268                         m.intensity <- Some intensity'
 269                         m.elements.UnionWith(nextElements)
 270                     stop <- true
 271
 272     for m in maxima do
 273         if m.state = Validated
 274         then
 275             match m.intensity with
 276             | Some i ->
 277                 for p in m.elements do
 278                     img.Data.[p.Y, p.X, 0] <- i
 279             | _ -> ()
 280     ()
 281
 282
 283 // Zhang and Suen algorithm.
 284 // Modify 'mat' in place.
 285 let thin (mat: Matrix<byte>) =
 286     let w = mat.Width
 287     let h = mat.Height
 288     let mutable data1 = mat.Data
 289     let mutable data2 = Array2D.copy data1
 290
 291     let mutable pixelChanged = true
 292     let mutable oddIteration = true
 293
 294     while pixelChanged do
 295         pixelChanged <- false
 296         for i in 0..h-1 do
 297             for j in 0..w-1 do
 298                 if data1.[i, j] = 1uy
 299                 then
 300                     let p2 = if i = 0 then 0uy else data1.[i-1, j]
 301                     let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
 302                     let p4 = if j = w-1 then 0uy else data1.[i, j+1]
 303                     let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
 304                     let p6 = if i = h-1 then 0uy else data1.[i+1, j]
 305                     let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
 306                     let p8 = if j = 0 then 0uy else data1.[i, j-1]
 307                     let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
 308
 309                     let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
 310                     if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
 311                         (if p2 = 0uy && p3 = 1uy then 1 else 0) +
 312                         (if p3 = 0uy && p4 = 1uy then 1 else 0) +
 313                         (if p4 = 0uy && p5 = 1uy then 1 else 0) +
 314                         (if p5 = 0uy && p6 = 1uy then 1 else 0) +
 315                         (if p6 = 0uy && p7 = 1uy then 1 else 0) +
 316                         (if p7 = 0uy && p8 = 1uy then 1 else 0) +
 317                         (if p8 = 0uy && p9 = 1uy then 1 else 0) +
 318                         (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
 319                         if oddIteration
 320                         then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
 321                         else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
 322                     then
 323                         data2.[i, j] <- 0uy
 324                         pixelChanged <- true
 325                 else
 326                     data2.[i, j] <- 0uy
 327
 328         oddIteration <- not oddIteration
 329         let tmp = data1
 330         data1 <- data2
 331         data2 <- tmp
 332
 333
 334 // FIXME: replace by a queue or stack.
 335 let pop (l: List<'a>) : 'a =
 336     let n = l.[l.Count - 1]
 337     l.RemoveAt(l.Count - 1)
 338     n
 339
 340 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
 341 // Modify 'mat' in place.
 342 let removeArea (mat: Matrix<byte>) (areaSize: int) =
 343     let neighbors = [|
 344         (-1,  0) // p2
 345         (-1,  1) // p3
 346         ( 0,  1) // p4
 347         ( 1,  1) // p5
 348         ( 1,  0) // p6
 349         ( 1, -1) // p7
 350         ( 0, -1) // p8
 351         (-1, -1) |] // p9
 352
 353     let mat' = new Matrix<byte>(mat.Size)
 354     let w = mat'.Width
 355     let h = mat'.Height
 356     mat.CopyTo(mat')
 357
 358     let data = mat.Data
 359     let data' = mat'.Data
 360
 361     for i in 0..h-1 do
 362         for j in 0..w-1 do
 363             if data'.[i, j] = 1uy
 364             then
 365                 let neighborhood = List<(int*int)>()
 366                 let neighborsToCheck = List<(int*int)>()
 367                 neighborsToCheck.Add((i, j))
 368                 data'.[i, j] <- 0uy
 369
 370                 while neighborsToCheck.Count > 0 do
 371                     let (ci, cj) = pop neighborsToCheck
 372                     neighborhood.Add((ci, cj))
 373                     for (ni, nj) in neighbors do
 374                         let pi = ci + ni
 375                         let pj = cj + nj
 376                         if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
 377                         then
 378                             neighborsToCheck.Add((pi, pj))
 379                             data'.[pi, pj] <- 0uy
 380                 if neighborhood.Count <= areaSize
 381                 then
 382                     for (ni, nj) in neighborhood do
 383                         data.[ni, nj] <- 0uy
 384
 385 let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
 386     let w = img.Width
 387     let h = img.Height
 388
 389     let pointChecked = HashSet<Point>()
 390     let pointToCheck = List<Point>(startPoints);
 391
 392     let data = img.Data
 393
 394     while pointToCheck.Count > 0 do
 395         let next = pop pointToCheck
 396         pointChecked.Add(next) |> ignore
 397         for ny in -1 .. 1 do
 398             for nx in -1 .. 1 do
 399                 if ny <> 0 && nx <> 0
 400                 then
 401                     let p = Point(next.X + nx, next.Y + ny)
 402                     if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
 403                     then
 404                         pointToCheck.Add(p)
 405
 406     List<Point>(pointChecked)
 407
 408
 409 let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
 410     img.Save(filepath)
 411
 412
 413 let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
 414     use img = new Image<Gray, 'TDeph>(mat.Size)
 415     mat.CopyTo(img)
 416     saveImg img filepath
 417
 418 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
 419     img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
 420
 421 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
 422     img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
 423
 424 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) =
 425
 426     if alpha >= 1.0
 427     then
 428         img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
 429     else
 430         let windowPosX = e.Cx - e.A - 5.0
 431         let gapX = windowPosX - (float (int windowPosX))
 432
 433         let windowPosY = e.Cy - e.A - 5.0
 434         let gapY = windowPosY - (float (int windowPosY))
 435
 436         let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
 437
 438         img.ROI <- roi
 439         if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
 440         then
 441             use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
 442             i.Draw(Ellipse(PointF(float32 <| (e.A + 5. + gapX) , float32 <| (e.A + 5. + gapY)), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
 443             CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
 444         img.ROI <- Rectangle.Empty
 445
 446
 447 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha: float) =
 448     List.iter (fun e -> drawEllipse img e color alpha) ellipses
 449
 450
 451 let rngCell =  System.Random()
 452 let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
 453     if drawCellContent
 454     then
 455         let colorB = rngCell.Next(20, 70)
 456         let colorG = rngCell.Next(20, 70)
 457         let colorR = rngCell.Next(20, 70)
 458
 459         for y in 0 .. c.elements.Height - 1 do
 460             for x in 0 .. c.elements.Width - 1 do
 461                 if c.elements.[y, x] > 0uy
 462                 then
 463                     let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
 464                     let b = img.Data.[y + dy, x + dx, 0] |> int
 465                     let g = img.Data.[y + dy, x + dx, 1] |> int
 466                     let r = img.Data.[y + dy, x + dx, 2] |> int
 467                     img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
 468                     img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
 469                     img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
 470
 471     let crossColor, crossColor2 =
 472         match c.cellClass with
 473         | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
 474         | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
 475         | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
 476
 477     drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
 478     drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
 479
 480     drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
 481     drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
 482
 483
 484 let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
 485     List.iter (fun c -> drawCell img drawCellContent c) cells