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