Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs

   1 module ParasitemiaCore.Morpho
   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 Types
  12
  13 /// <summary>
  14 /// Remove M-adjacent pixels. It may be used after thinning.
  15 /// </summary>
  16 let suppressMAdjacency (img : Matrix<byte>) =
  17     let w = img.Width
  18     let h = img.Height
  19     for i = 1 to h - 2 do
  20         for j = 1 to w - 2 do
  21             if img.[i, j] > 0uy && img.Data.[i + 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i - 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i - 1, j - 1] = 0uy) then
  22                 img.[i, j] <- 0uy
  23     for i = 1 to h - 2 do
  24         for j = 1 to w - 2 do
  25             if img.[i, j] > 0uy && img.Data.[i - 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i + 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i + 1, j - 1] = 0uy) then
  26                 img.[i, j] <- 0uy
  27
  28 type ExtremumType =
  29     | Maxima = 1
  30     | Minima = 2
  31
  32 let inline findExtremum (img : Image<Gray, 'TDepth>) (extremumType : ExtremumType) : IEnumerable<Points>  when 'TDepth : unmanaged =
  33     let w = img.Width
  34     let h = img.Height
  35     let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
  36
  37     let imgData = img.Data
  38     let suppress : bool[,] = Array2D.zeroCreate h w
  39
  40     let result = List<List<Point>>()
  41
  42     let flood (start : Point) : List<List<Point>> =
  43         let sameLevelToCheck = Stack<Point>()
  44         let betterLevelToCheck = Stack<Point>()
  45         betterLevelToCheck.Push(start)
  46
  47         let result' = List<List<Point>>()
  48
  49         while betterLevelToCheck.Count > 0 do
  50             let p = betterLevelToCheck.Pop()
  51             if not suppress.[p.Y, p.X] then
  52                 suppress.[p.Y, p.X] <- true
  53                 sameLevelToCheck.Push(p)
  54                 let current = List<Point>()
  55
  56                 let mutable betterExists = false
  57
  58                 while sameLevelToCheck.Count > 0 do
  59                     let p' = sameLevelToCheck.Pop()
  60                     let currentLevel = imgData.[p'.Y, p'.X, 0]
  61                     current.Add(p') |> ignore
  62                     for i, j in se do
  63                         let ni = i + p'.Y
  64                         let nj = j + p'.X
  65                         if ni >= 0 && ni < h && nj >= 0 && nj < w then
  66                             let level = imgData.[ni, nj, 0]
  67                             let notSuppressed = not suppress.[ni, nj]
  68
  69                             if notSuppressed && level = currentLevel then
  70                                 suppress.[ni, nj] <- true
  71                                 sameLevelToCheck.Push(Point(nj, ni))
  72                             elif (if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel) then
  73                                 betterExists <- true
  74                                 if notSuppressed then
  75                                     betterLevelToCheck.Push(Point(nj, ni))
  76
  77                 if not betterExists then
  78                     result'.Add(current)
  79         result'
  80
  81     for i = 0 to h - 1 do
  82         for j = 0 to w - 1 do
  83             let maxima = flood (Point(j, i))
  84             if maxima.Count > 0 then
  85                 result.AddRange(maxima)
  86
  87     result.Select(fun l -> Points(l))
  88
  89 let inline findMaxima (img : Image<Gray, 'TDepth>) : IEnumerable<Points> when 'TDepth : unmanaged =
  90     findExtremum img ExtremumType.Maxima
  91
  92 let inline findMinima (img : Image<Gray, 'TDepth>) : IEnumerable<Points> when 'TDepth : unmanaged =
  93     findExtremum img ExtremumType.Minima
  94
  95 type PriorityQueue () =
  96     let size = 256
  97     let q : Points[] = Array.init size (fun i -> Points())
  98     let mutable highest = -1 // Value of the first elements of 'q'.
  99     let mutable lowest = size
 100
 101     member this.NextMax () : byte * Point =
 102         if this.IsEmpty then
 103             invalidOp "Queue is empty"
 104         else
 105             let l = q.[highest]
 106             let next = l.First()
 107             l.Remove(next) |> ignore
 108             let value = byte highest
 109
 110             if l.Count = 0 then
 111                 highest <- highest - 1
 112                 while highest > lowest && q.[highest].Count = 0 do
 113                     highest <- highest - 1
 114                 if highest = lowest then
 115                     highest <- -1
 116                     lowest <- size
 117
 118             value, next
 119
 120     member this.NextMin () : byte * Point =
 121         if this.IsEmpty then
 122             invalidOp "Queue is empty"
 123         else
 124             let l = q.[lowest + 1]
 125             let next = l.First()
 126             l.Remove(next) |> ignore
 127             let value = byte (lowest + 1)
 128
 129             if l.Count = 0 then
 130                 lowest <- lowest + 1
 131                 while lowest < highest && q.[lowest + 1].Count = 0 do
 132                     lowest <- lowest + 1
 133                 if highest = lowest then
 134                     highest <- -1
 135                     lowest <- size
 136
 137             value, next
 138
 139     member this.Max =
 140         highest |> byte
 141
 142     member this.Min =
 143         lowest + 1 |> byte
 144
 145     member this.Add (value : byte) (p : Point) =
 146         let vi = int value
 147
 148         if vi > highest then
 149             highest <- vi
 150         if vi <= lowest then
 151             lowest <- vi - 1
 152
 153         q.[vi].Add(p) |> ignore
 154
 155     member this.Remove (value : byte) (p : Point) =
 156         let vi = int value
 157         if q.[vi].Remove(p) && q.[vi].Count = 0 then
 158             if vi = highest then
 159                 highest <- highest - 1
 160                 while highest > lowest && q.[highest].Count = 0 do
 161                     highest <- highest - 1
 162             elif vi - 1 = lowest then
 163                 lowest <- lowest + 1
 164                 while lowest < highest && q.[lowest + 1].Count = 0 do
 165                     lowest <- lowest + 1
 166
 167             if highest = lowest then // The queue is now empty.
 168                 highest <- -1
 169                 lowest <- size
 170
 171     member this.IsEmpty =
 172         highest = -1
 173
 174     member this.Clear () =
 175         while highest > lowest  do
 176             q.[highest].Clear()
 177             highest <- highest - 1
 178         highest <- -1
 179         lowest <- size
 180
 181 type private AreaState =
 182     | Removed = 1
 183     | Unprocessed = 2
 184     | Validated = 3
 185
 186 type private AreaOperation =
 187     | Opening = 1
 188     | Closing = 2
 189
 190 [<AllowNullLiteral>]
 191 type private Area (elements : Points) =
 192     member this.Elements = elements
 193     member val Intensity = None with get, set
 194     member val State = AreaState.Unprocessed with get, set
 195
 196 let private areaOperation (img : Image<Gray, byte>) (area : int) (op : AreaOperation) =
 197     let w = img.Width
 198     let h = img.Height
 199     let imgData = img.Data
 200     let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
 201
 202     let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
 203
 204     let pixels : Area[,] = Array2D.create h w null
 205     for m in areas do
 206         for e in m.Elements do
 207             pixels.[e.Y, e.X] <- m
 208
 209     let queue = PriorityQueue()
 210
 211     let addEdgeToQueue (elements : Points) =
 212         for p in elements do
 213             for i, j in se do
 214                 let ni = i + p.Y
 215                 let nj = j + p.X
 216                 let p' = Point(nj, ni)
 217                 if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p')) then
 218                     queue.Add (imgData.[ni, nj, 0]) p'
 219
 220     // Reverse order is quicker.
 221     for i = areas.Count - 1 downto 0 do
 222         let m = areas.[i]
 223         if m.Elements.Count <= area && m.State <> AreaState.Removed then
 224             queue.Clear()
 225             addEdgeToQueue m.Elements
 226
 227             let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min
 228             let nextElements = Points()
 229
 230             let mutable stop = false
 231             while not stop do
 232                 let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
 233                 let mutable merged = false
 234
 235                 if intensity' = intensity then // The intensity doesn't change.
 236                     if m.Elements.Count + nextElements.Count + 1 > area then
 237                         m.State <- AreaState.Validated
 238                         m.Intensity <- Some intensity
 239                         stop <- true
 240                     else
 241                         nextElements.Add(p) |> ignore
 242
 243                 elif (if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity) then
 244                     m.Elements.UnionWith(nextElements)
 245                     for e in nextElements do
 246                         pixels.[e.Y, e.X] <- m
 247
 248                     if m.Elements.Count = area then
 249                         m.State <- AreaState.Validated
 250                         m.Intensity <- Some (intensity')
 251                         stop <- true
 252                     else
 253                         intensity <- intensity'
 254                         nextElements.Clear()
 255                         nextElements.Add(p) |> ignore
 256
 257                 else
 258                     match pixels.[p.Y, p.X] with
 259                     | null -> ()
 260                     | m' ->
 261                         if m'.Elements.Count + m.Elements.Count <= area 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 = if op = AreaOperation.Opening then queue.Max else queue.Min
 269                             if intensityMax <> intensity then
 270                                 intensity <- intensityMax
 271                                 nextElements.Clear()
 272                             merged <- true
 273
 274                     if not merged then
 275                         m.State <- AreaState.Validated
 276                         m.Intensity <- Some (intensity)
 277                         stop <- true
 278
 279                 if not stop && not merged then
 280                     for i, j in se do
 281                         let ni = i + p.Y
 282                         let nj = j + p.X
 283                         let p' = Point(nj, ni)
 284                         if ni < 0 || ni >= h || nj < 0 || nj >= w then
 285                             m.State <- AreaState.Validated
 286                             m.Intensity <- Some (intensity)
 287                             stop <- true
 288                         elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) then
 289                             queue.Add (imgData.[ni, nj, 0]) p'
 290
 291                 if queue.IsEmpty then
 292                     if m.Elements.Count + nextElements.Count <= area then
 293                         m.State <- AreaState.Validated
 294                         m.Intensity <- Some intensity'
 295                         m.Elements.UnionWith(nextElements)
 296                     stop <- true
 297
 298     for m in areas do
 299         if m.State = AreaState.Validated then
 300             match m.Intensity with
 301             | Some i ->
 302                 for p in m.Elements do
 303                     imgData.[p.Y, p.X, 0] <- i
 304             | _ -> ()
 305     ()
 306
 307 /// <summary>
 308 /// Area opening on byte image.
 309 /// </summary>
 310 let areaOpen (img : Image<Gray, byte>) (area : int) =
 311     areaOperation img area AreaOperation.Opening
 312
 313 /// <summary>
 314 /// Area closing on byte image.
 315 /// </summary>
 316 let areaClose (img : Image<Gray, byte>) (area : int) =
 317     areaOperation img area AreaOperation.Closing
 318
 319 // A simpler algorithm than 'areaOpen' on byte image but slower.
 320 let areaOpen2 (img : Image<Gray, byte>) (area : int) =
 321     let w = img.Width
 322     let h = img.Height
 323     let imgData = img.Data
 324     let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
 325
 326     let histogram = Array.zeroCreate 256
 327     for i = 0 to h - 1 do
 328         for j = 0 to w - 1 do
 329             let v = imgData.[i, j, 0] |> int
 330             histogram.[v] <- histogram.[v] + 1
 331
 332     let flooded : bool[,] = Array2D.zeroCreate h w
 333
 334     let pointsChecked = HashSet<Point>()
 335     let pointsToCheck = Stack<Point>()
 336
 337     for level = 255 downto 0 do
 338         let mutable n = histogram.[level]
 339         if n > 0 then
 340             for i = 0 to h - 1 do
 341                 for j = 0 to w - 1 do
 342                     if not flooded.[i, j] && imgData.[i, j, 0] = byte level then
 343                         let mutable maxNeighborValue = 0uy
 344                         pointsChecked.Clear()
 345                         pointsToCheck.Clear()
 346                         pointsToCheck.Push(Point(j, i))
 347
 348                         while pointsToCheck.Count > 0 do
 349                             let next = pointsToCheck.Pop()
 350                             pointsChecked.Add(next) |> ignore
 351                             flooded.[next.Y, next.X] <- true
 352
 353                             for nx, ny in se do
 354                                 let p = Point(next.X + nx, next.Y + ny)
 355                                 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h then
 356                                     let v = imgData.[p.Y, p.X, 0]
 357                                     if v = byte level then
 358                                         if not (pointsChecked.Contains(p)) then
 359                                             pointsToCheck.Push(p)
 360                                     elif v > maxNeighborValue then
 361                                         maxNeighborValue <- v
 362
 363                         if int maxNeighborValue < level && pointsChecked.Count <= area then
 364                             for p in pointsChecked do
 365                                 imgData.[p.Y, p.X, 0] <- maxNeighborValue
 366
 367 [<AllowNullLiteral>]
 368 type Island (cmp : IComparer<float32>) =
 369     member val Shore = Heap.Heap<float32, Point>(cmp) with get
 370     member val Level = 0.f with get, set
 371     member val Surface = 0 with get, set
 372     member this.IsInfinite = this.Surface = Int32.MaxValue
 373
 374 let private areaOperationF (img : Image<Gray, float32>) (areas : (int * 'a) list) (f : ('a -> float32 -> unit) option) (op : AreaOperation) =
 375     let w = img.Width
 376     let h = img.Height
 377     let earth = img.Data
 378     let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
 379
 380     let comparer =
 381         if op = AreaOperation.Opening then
 382             { new IComparer<float32> with member this.Compare(v1, v2) = v1.CompareTo(v2) }
 383         else
 384             { new IComparer<float32> with member this.Compare(v1, v2) = v2.CompareTo(v1) }
 385
 386     let ownership : Island[,] = Array2D.create h w null
 387
 388     // Initialize islands with their shore.
 389     let islands = List<Island>()
 390     let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima
 391     for e in extremum do
 392         let island =
 393             let p = e.First()
 394             Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count)
 395         islands.Add(island)
 396         let shorePoints = Points()
 397         for p in e do
 398             ownership.[p.Y, p.X] <- island
 399             for i, j in se do
 400                 let ni = i + p.Y
 401                 let nj = j + p.X
 402                 let neighbor = Point(nj, ni)
 403                 if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor)) then
 404                     shorePoints.Add(neighbor) |> ignore
 405                     island.Shore.Add earth.[ni, nj, 0] neighbor
 406
 407     for area, obj in areas do
 408         for island in islands do
 409             let mutable stop = island.Shore.IsEmpty
 410
 411             // 'true' if 'p' is owned or adjacent to 'island'.
 412             let inline ownedOrAdjacent (p : Point) : bool =
 413                 ownership.[p.Y, p.X] = island ||
 414                 (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
 415                 (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
 416                 (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
 417                 (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
 418
 419             while not stop && island.Surface < area do
 420                 let level, next = island.Shore.Max
 421                 let other = ownership.[next.Y, next.X]
 422                 if other = island then // During merging, some points on the shore may be owned by the island itself -> ignored.
 423                     island.Shore.RemoveNext ()
 424                 else
 425                     if not <| Object.ReferenceEquals(other, null) then
 426                         // We touching another island.
 427                         if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0 then
 428                             stop <- true
 429                         else // We can merge 'other' into 'surface'.
 430                             island.Surface <- island.Surface + other.Surface
 431                             island.Level <- other.Level
 432                             // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.Level
 433                             for l, p in other.Shore do
 434                                 let mutable currentY = p.Y + 1
 435                                 while currentY < h && ownership.[currentY, p.X] = other do
 436                                     ownership.[currentY, p.X] <- island
 437                                     currentY <- currentY + 1
 438                                 island.Shore.Add l p
 439                             other.Shore.Clear()
 440
 441                     elif comparer.Compare(level, island.Level) > 0 then
 442                         stop <- true
 443                     else
 444                         island.Shore.RemoveNext ()
 445                         for i, j in se do
 446                             let ni = i + next.Y
 447                             let nj = j + next.X
 448                             if ni < 0 || ni >= h || nj < 0 || nj >= w then
 449                                 island.Surface <- Int32.MaxValue
 450                                 stop <- true
 451                             else
 452                                 let neighbor = Point(nj, ni)
 453                                 if not <| ownedOrAdjacent neighbor then
 454                                     island.Shore.Add earth.[ni, nj, 0] neighbor
 455                         if not stop then
 456                             ownership.[next.Y, next.X] <- island
 457                             island.Level <- level
 458                             island.Surface <- island.Surface + 1
 459
 460         let mutable diff = 0.f
 461
 462         for i = 0 to h - 1 do
 463             for j = 0 to w - 1 do
 464                 match ownership.[i, j] with
 465                 | null -> ()
 466                 | island ->
 467                     let l = island.Level
 468                     diff <- diff + l - earth.[i, j, 0]
 469                     earth.[i, j, 0] <- l
 470
 471         match f with
 472         | Some f' -> f' obj diff
 473         | _ -> ()
 474     ()
 475
 476 /// <summary>
 477 /// Area opening on float image.
 478 /// </summary>
 479 let areaOpenF (img : Image<Gray, float32>) (area : int) =
 480     areaOperationF img [ area, () ] None AreaOperation.Opening
 481
 482 /// <summary>
 483 /// Area closing on float image.
 484 /// </summary>
 485 let areaCloseF (img : Image<Gray, float32>) (area : int) =
 486     areaOperationF img [ area, () ] None AreaOperation.Closing
 487
 488 /// <summary>
 489 /// Area closing on float image with different areas. Given areas must be sorted increasingly.
 490 /// For each area the function 'f' is called with the associated area value of type 'a and the volume difference
 491 /// Between the previous and the current closing.
 492 /// </summary>
 493 let areaOpenFWithFun (img : Image<Gray, float32>) (areas : (int * 'a) list) (f : 'a -> float32 -> unit) =
 494     areaOperationF img areas (Some f) AreaOperation.Opening
 495
 496 /// <summary>
 497 /// Same as 'areaOpenFWithFun' for closing operation.
 498 /// </summary>
 499 let areaCloseFWithFun (img : Image<Gray, float32>) (areas : (int * 'a) list) (f : 'a -> float32 -> unit) =
 500     areaOperationF img areas (Some f) AreaOperation.Closing
 501
 502 /// <summary>
 503 /// Zhang and Suen thinning algorithm.
 504 /// Modify 'mat' in place.
 505 /// </summary>
 506 let thin (mat : Matrix<byte>) =
 507     let w = mat.Width
 508     let h = mat.Height
 509     let mutable data1 = mat.Data
 510     let mutable data2 = Array2D.copy data1
 511
 512     let mutable pixelChanged = true
 513     let mutable oddIteration = true
 514
 515     while pixelChanged do
 516         pixelChanged <- false
 517         for i = 0 to h - 1 do
 518             for j = 0 to w - 1 do
 519                 if data1.[i, j] = 1uy then
 520                     let p2 = if i = 0 then 0uy else data1.[i-1, j]
 521                     let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
 522                     let p4 = if j = w-1 then 0uy else data1.[i, j+1]
 523                     let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
 524                     let p6 = if i = h-1 then 0uy else data1.[i+1, j]
 525                     let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
 526                     let p8 = if j = 0 then 0uy else data1.[i, j-1]
 527                     let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
 528
 529                     let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
 530                     if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
 531                         (if p2 = 0uy && p3 = 1uy then 1 else 0) +
 532                         (if p3 = 0uy && p4 = 1uy then 1 else 0) +
 533                         (if p4 = 0uy && p5 = 1uy then 1 else 0) +
 534                         (if p5 = 0uy && p6 = 1uy then 1 else 0) +
 535                         (if p6 = 0uy && p7 = 1uy then 1 else 0) +
 536                         (if p7 = 0uy && p8 = 1uy then 1 else 0) +
 537                         (if p8 = 0uy && p9 = 1uy then 1 else 0) +
 538                         (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
 539                         if oddIteration then
 540                             p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
 541                         else
 542                             p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
 543                     then
 544                         data2.[i, j] <- 0uy
 545                         pixelChanged <- true
 546                 else
 547                     data2.[i, j] <- 0uy
 548
 549         oddIteration <- not oddIteration
 550         let tmp = data1
 551         data1 <- data2
 552         data2 <- tmp
 553
 554 /// <summary>
 555 /// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
 556 /// Modify 'mat' in place.
 557 /// </summary>
 558 let removeArea (mat : Matrix<byte>) (areaSize : int) =
 559     let neighbors = [|
 560         (-1,  0) // p2
 561         (-1,  1) // p3
 562         ( 0,  1) // p4
 563         ( 1,  1) // p5
 564         ( 1,  0) // p6
 565         ( 1, -1) // p7
 566         ( 0, -1) // p8
 567         (-1, -1) |] // p9
 568
 569     use mat' = new Matrix<byte>(mat.Size)
 570     let w = mat'.Width
 571     let h = mat'.Height
 572     mat.CopyTo(mat')
 573
 574     let data = mat.Data
 575     let data' = mat'.Data
 576
 577     for i = 0 to h - 1 do
 578         for j = 0 to w - 1 do
 579             if data'.[i, j] = 1uy then
 580                 let neighborhood = List<Point>()
 581                 let neighborsToCheck = Stack<Point>()
 582                 neighborsToCheck.Push(Point(j, i))
 583                 data'.[i, j] <- 0uy
 584
 585                 while neighborsToCheck.Count > 0 do
 586                     let n = neighborsToCheck.Pop()
 587                     neighborhood.Add(n)
 588                     for (ni, nj) in neighbors do
 589                         let pi = n.Y + ni
 590                         let pj = n.X + nj
 591                         if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy then
 592                             neighborsToCheck.Push(Point(pj, pi))
 593                             data'.[pi, pj] <- 0uy
 594                 if neighborhood.Count <= areaSize then
 595                     for n in neighborhood do
 596                         data.[n.Y, n.X] <- 0uy
 597
 598 let connectedComponents (img : Image<Gray, byte>) (startPoints : List<Point>) : Points =
 599     let w = img.Width
 600     let h = img.Height
 601
 602     let pointChecked = Points()
 603     let pointToCheck = Stack<Point>(startPoints);
 604
 605     let data = img.Data
 606
 607     while pointToCheck.Count > 0 do
 608         let next = pointToCheck.Pop()
 609         pointChecked.Add(next) |> ignore
 610         for ny = -1 to 1 do
 611             for nx = -1 to 1 do
 612                 if ny <> 0 && nx <> 0 then
 613                     let p = Point(next.X + nx, next.Y + ny)
 614                     if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p) then
 615                         pointToCheck.Push(p)
 616
 617     pointChecked