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