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                             for l, p in other.Shore do
 433                                 let mutable currentY = p.Y + 1
 434                                 while currentY < h && ownership.[currentY, p.X] = other do
 435                                     ownership.[currentY, p.X] <- island
 436                                     currentY <- currentY + 1
 437                                 island.Shore.Add l p
 438                             other.Shore.Clear ()
 439
 440                     elif comparer.Compare (level, island.Level) > 0 then
 441                         stop <- true
 442                     else
 443                         island.Shore.RemoveNext ()
 444                         for i, j in se do
 445                             let ni = i + next.Y
 446                             let nj = j + next.X
 447                             if ni < 0 || ni >= h || nj < 0 || nj >= w then
 448                                 island.Surface <- Int32.MaxValue
 449                                 stop <- true
 450                             else
 451                                 let neighbor = Point (nj, ni)
 452                                 if not <| ownedOrAdjacent neighbor then
 453                                     island.Shore.Add earth.[ni, nj, 0] neighbor
 454                         if not stop then
 455                             ownership.[next.Y, next.X] <- island
 456                             island.Level <- level
 457                             island.Surface <- island.Surface + 1
 458
 459         let mutable diff = 0.f
 460
 461         for i = 0 to h - 1 do
 462             for j = 0 to w - 1 do
 463                 match ownership.[i, j] with
 464                 | null -> ()
 465                 | island ->
 466                     let l = island.Level
 467                     diff <- diff + l - earth.[i, j, 0]
 468                     earth.[i, j, 0] <- l
 469
 470         match f with
 471         | Some f' -> f' obj diff
 472         | _ -> ()
 473     ()
 474
 475 /// <summary>
 476 /// Area opening on float image.
 477 /// </summary>
 478 let areaOpenF (img : Image<Gray, float32>) (area : int) =
 479     areaOperationF img [ area, () ] None AreaOperation.Opening
 480
 481 /// <summary>
 482 /// Area closing on float image.
 483 /// </summary>
 484 let areaCloseF (img : Image<Gray, float32>) (area : int) =
 485     areaOperationF img [ area, () ] None AreaOperation.Closing
 486
 487 /// <summary>
 488 /// Area closing on float image with different areas. Given areas must be sorted increasingly.
 489 /// For each area the function 'f' is called with the associated area value of type 'a and the volume difference
 490 /// Between the previous and the current closing.
 491 /// </summary>
 492 let areaOpenFWithFun (img : Image<Gray, float32>) (areas : (int * 'a) list) (f : 'a -> float32 -> unit) =
 493     areaOperationF img areas (Some f) AreaOperation.Opening
 494
 495 /// <summary>
 496 /// Same as 'areaOpenFWithFun' for closing operation.
 497 /// </summary>
 498 let areaCloseFWithFun (img : Image<Gray, float32>) (areas : (int * 'a) list) (f : 'a -> float32 -> unit) =
 499     areaOperationF img areas (Some f) AreaOperation.Closing
 500
 501 /// <summary>
 502 /// Zhang and Suen thinning algorithm.
 503 /// Modify 'mat' in place.
 504 /// </summary>
 505 let thin (mat : Matrix<byte>) =
 506     let w = mat.Width
 507     let h = mat.Height
 508     let mutable data1 = mat.Data
 509     let mutable data2 = Array2D.copy data1
 510
 511     let mutable pixelChanged = true
 512     let mutable oddIteration = true
 513
 514     while pixelChanged do
 515         pixelChanged <- false
 516         for i = 0 to h - 1 do
 517             for j = 0 to w - 1 do
 518                 if data1.[i, j] = 1uy then
 519                     let p2 = if i = 0 then 0uy else data1.[i-1, j]
 520                     let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
 521                     let p4 = if j = w-1 then 0uy else data1.[i, j+1]
 522                     let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
 523                     let p6 = if i = h-1 then 0uy else data1.[i+1, j]
 524                     let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
 525                     let p8 = if j = 0 then 0uy else data1.[i, j-1]
 526                     let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
 527
 528                     let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
 529                     if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
 530                         (if p2 = 0uy && p3 = 1uy then 1 else 0) +
 531                         (if p3 = 0uy && p4 = 1uy then 1 else 0) +
 532                         (if p4 = 0uy && p5 = 1uy then 1 else 0) +
 533                         (if p5 = 0uy && p6 = 1uy then 1 else 0) +
 534                         (if p6 = 0uy && p7 = 1uy then 1 else 0) +
 535                         (if p7 = 0uy && p8 = 1uy then 1 else 0) +
 536                         (if p8 = 0uy && p9 = 1uy then 1 else 0) +
 537                         (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
 538                         if oddIteration then
 539                             p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
 540                         else
 541                             p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
 542                     then
 543                         data2.[i, j] <- 0uy
 544                         pixelChanged <- true
 545                 else
 546                     data2.[i, j] <- 0uy
 547
 548         oddIteration <- not oddIteration
 549         let tmp = data1
 550         data1 <- data2
 551         data2 <- tmp
 552
 553 /// <summary>
 554 /// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
 555 /// Modify 'mat' in place.
 556 /// </summary>
 557 let removeArea (mat : Matrix<byte>) (areaSize : int) =
 558     let neighbors = [|
 559         (-1,  0) // p2
 560         (-1,  1) // p3
 561         ( 0,  1) // p4
 562         ( 1,  1) // p5
 563         ( 1,  0) // p6
 564         ( 1, -1) // p7
 565         ( 0, -1) // p8
 566         (-1, -1) |] // p9
 567
 568     use mat' = new Matrix<byte> (mat.Size)
 569     let w = mat'.Width
 570     let h = mat'.Height
 571     mat.CopyTo mat'
 572
 573     let data = mat.Data
 574     let data' = mat'.Data
 575
 576     for i = 0 to h - 1 do
 577         for j = 0 to w - 1 do
 578             if data'.[i, j] = 1uy then
 579                 let neighborhood = List<Point> ()
 580                 let neighborsToCheck = Stack<Point> ()
 581                 neighborsToCheck.Push (Point (j, i))
 582                 data'.[i, j] <- 0uy
 583
 584                 while neighborsToCheck.Count > 0 do
 585                     let n = neighborsToCheck.Pop ()
 586                     neighborhood.Add n
 587                     for (ni, nj) in neighbors do
 588                         let pi = n.Y + ni
 589                         let pj = n.X + nj
 590                         if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy then
 591                             neighborsToCheck.Push (Point (pj, pi))
 592                             data'.[pi, pj] <- 0uy
 593                 if neighborhood.Count <= areaSize then
 594                     for n in neighborhood do
 595                         data.[n.Y, n.X] <- 0uy
 596
 597 let connectedComponents (img : Image<Gray, byte>) (startPoints : List<Point>) : Points =
 598     let w = img.Width
 599     let h = img.Height
 600
 601     let pointChecked = Points ()
 602     let pointToCheck = Stack<Point> startPoints;
 603
 604     let data = img.Data
 605
 606     while pointToCheck.Count > 0 do
 607         let next = pointToCheck.Pop ()
 608         pointChecked.Add next |> ignore
 609         for ny = -1 to 1 do
 610             for nx = -1 to 1 do
 611                 if ny <> 0 && nx <> 0 then
 612                     let p = Point (next.X + nx, next.Y + ny)
 613                     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
 614                         pointToCheck.Push p
 615
 616     pointChecked