1
module ParasitemiaCore.Morpho
5 open System.Collections.Generic
14 /// Remove M-adjacent pixels. It may be used after thinning.
16 let suppressMAdjacency (img
: Matrix<byte
>) =
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
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
32 let inline findExtremum
(img : Image<Gray, 'TDepth>) (extremumType : ExtremumType) : IEnumerable<Points> when 'TDepth : unmanaged
=
35 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
37 let imgData = img.Data
38 let suppress : bool[,] = Array2D.zeroCreate
h w
40 let result = List<List<Point>>()
42 let flood (start
: Point) : List<List<Point>> =
43 let sameLevelToCheck = Stack<Point>()
44 let betterLevelToCheck = Stack<Point>()
45 betterLevelToCheck.Push(start
)
47 let result' = List<List<Point>>()
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>()
56 let mutable betterExists = false
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
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]
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
75 betterLevelToCheck.Push(Point(nj, ni))
77 if not betterExists then
83 let maxima = flood (Point(j, i))
84 if maxima.Count > 0 then
85 result.AddRange(maxima)
87 result.Select(fun l -> Points(l))
89 let inline findMaxima (img : Image<Gray, 'TDepth>) : IEnumerable<Points> when 'TDepth : unmanaged =
90 findExtremum img ExtremumType.Maxima
92 let inline findMinima (img : Image<Gray, 'TDepth>) : IEnumerable<Points> when 'TDepth : unmanaged =
93 findExtremum img ExtremumType.Minima
95 type PriorityQueue () =
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
101 member this.NextMax () : byte * Point =
103 invalidOp "Queue is empty"
107 l.Remove(next) |> ignore
108 let value = byte highest
111 highest <- highest - 1
112 while highest > lowest && q.[highest].Count = 0 do
113 highest <- highest - 1
114 if highest = lowest then
120 member this.NextMin () : byte * Point =
122 invalidOp "Queue is empty"
124 let l = q.[lowest + 1]
126 l.Remove(next) |> ignore
127 let value = byte (lowest + 1)
131 while lowest < highest && q.[lowest + 1].Count = 0 do
133 if highest = lowest then
145 member this.Add (value : byte) (p : Point) =
153 q.[vi].Add(p) |> ignore
155 member this.Remove (value : byte) (p : Point) =
157 if q.[vi].Remove(p) && q.[vi].Count = 0 then
159 highest <- highest - 1
160 while highest > lowest && q.[highest].Count = 0 do
161 highest <- highest - 1
162 elif vi - 1 = lowest then
164 while lowest < highest && q.[lowest + 1].Count = 0 do
167 if highest = lowest then // The queue is now empty.
171 member this.IsEmpty =
174 member this.Clear () =
175 while highest > lowest do
177 highest <- highest - 1
181 type private AreaState =
186 type private AreaOperation =
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
196 let private areaOperation (img : Image<Gray, byte>) (area : int) (op : AreaOperation) =
199 let imgData = img.Data
200 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
202 let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
204 let pixels : Area[,] = Array2D.create h w null
206 for e in m.Elements do
207 pixels.[e.Y, e.X] <- m
209 let queue = PriorityQueue()
211 let addEdgeToQueue (elements : Points) =
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'
220 // Reverse order is quicker.
221 for i
= areas.Count - 1 downto 0 do
223 if m.Elements.Count <= area
&& m.State <> AreaState.Removed then
225 addEdgeToQueue m.Elements
227 let mutable intensity = if op
= AreaOperation.Opening then queue.Max else queue.Min
228 let nextElements = Points()
230 let mutable stop = false
232 let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
233 let mutable merged = false
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
241 nextElements.Add(p) |> ignore
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
248 if m.Elements.Count = area
then
249 m.State <- AreaState.Validated
250 m.Intensity <- Some (intensity')
253 intensity <- intensity'
255 nextElements.Add(p) |> ignore
258 match pixels.[p.Y, p.X] with
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
275 m.State <- AreaState.Validated
276 m.Intensity <- Some (intensity)
279 if not stop && not merged then
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)
288 elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) then
289 queue.Add (imgData.[ni, nj, 0]) p'
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)
299 if m.State = AreaState.Validated then
300 match m.Intensity with
302 for p in m.Elements do
303 imgData.[p.Y, p.X, 0] <- i
308 /// Area opening on byte image.
310 let areaOpen (img : Image<Gray, byte>) (area : int) =
311 areaOperation img area AreaOperation.Opening
314 /// Area closing on byte image.
316 let areaClose (img : Image<Gray, byte>) (area : int) =
317 areaOperation img area AreaOperation.Closing
319 // A simpler algorithm than 'areaOpen' on byte image but slower.
320 let areaOpen2 (img : Image<Gray, byte>) (area : int) =
323 let imgData = img.Data
324 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
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
332 let flooded : bool[,] = Array2D.zeroCreate h w
334 let pointsChecked = HashSet<Point>()
335 let pointsToCheck = Stack<Point>()
337 for level = 255 downto 0 do
338 let mutable n = histogram.[level]
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))
348 while pointsToCheck.Count > 0 do
349 let next = pointsToCheck.Pop()
350 pointsChecked.Add(next) |> ignore
351 flooded.[next.Y, next.X] <- true
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
363 if int maxNeighborValue < level && pointsChecked.Count <= area then
364 for p in pointsChecked do
365 imgData.[p.Y, p.X, 0] <- maxNeighborValue
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
374 let private areaOperationF (img : Image<Gray, float32>) (areas : (int * 'a
) list
) (f
: ('a -> float32 -> unit) option) (op : AreaOperation) =
378 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
381 if op = AreaOperation.Opening then
382 { new IComparer<float32> with member this.Compare(v1, v2) = v1.CompareTo(v2) }
384 { new IComparer<float32> with member this.Compare(v1, v2) = v2.CompareTo(v1) }
386 let ownership : Island[,] = Array2D.create h w null
388 // Initialize islands with their shore.
389 let islands = List<Island>()
390 let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima
394 Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count)
396 let shorePoints = Points()
398 ownership.[p.Y, p.X] <- island
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
407 for area, obj in areas do
408 for island in islands do
409 let mutable stop = island.Shore.IsEmpty
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)
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 ()
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
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
441 elif comparer.Compare(level, island.Level) > 0 then
444 island.Shore.RemoveNext ()
448 if ni < 0 || ni >= h || nj < 0 || nj >= w then
449 island.Surface <- Int32.MaxValue
452 let neighbor = Point(nj, ni)
453 if not <| ownedOrAdjacent neighbor then
454 island.Shore.Add earth.[ni, nj, 0] neighbor
456 ownership.[next.Y, next.X] <- island
457 island.Level <- level
458 island.Surface <- island.Surface + 1
460 let mutable diff = 0.f
462 for i = 0 to h - 1 do
463 for j = 0 to w - 1 do
464 match ownership.[i, j] with
468 diff <- diff + l - earth.[i, j, 0]
472 | Some f' -> f
' obj diff
477 /// Area opening on float image.
479 let areaOpenF (img : Image<Gray, float32>) (area : int) =
480 areaOperationF img [ area, () ] None AreaOperation.Opening
483 /// Area closing on float image.
485 let areaCloseF (img : Image<Gray, float32>) (area : int) =
486 areaOperationF img [ area, () ] None AreaOperation.Closing
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.
493 let areaOpenFWithFun (img : Image<Gray, float32
>) (areas : (int * 'a) list) (f : 'a
-> float32
-> unit) =
494 areaOperationF
img areas (Some f) AreaOperation.Opening
497 /// Same as 'areaOpenFWithFun' for closing operation.
499 let areaCloseFWithFun (img : Image<Gray, float32
>) (areas : (int * 'a) list) (f : 'a
-> float32
-> unit) =
500 areaOperationF
img areas (Some f) AreaOperation.Closing
503 /// Zhang and Suen thinning algorithm.
504 /// Modify 'mat' in place.
506 let thin (mat
: Matrix<byte
>) =
509 let mutable data1 = mat
.Data
510 let mutable data2 = Array2D.copy
data1
512 let mutable pixelChanged = true
513 let mutable oddIteration = true
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]
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 &&
540 p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
542 p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
549 oddIteration <- not oddIteration
555 /// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
556 /// Modify 'mat' in place.
558 let removeArea (mat
: Matrix<byte
>) (areaSize
: int) =
569 use mat' = new Matrix<byte>(mat.Size)
575 let data' = mat'.Data
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))
585 while neighborsToCheck.Count > 0 do
586 let n = neighborsToCheck.Pop()
588 for (ni, nj) in neighbors do
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
598 let connectedComponents (img : Image<Gray, byte
>) (startPoints
: List<Point>) : Points =
602 let pointChecked = Points()
603 let pointToCheck = Stack<Point>(startPoints
);
607 while pointToCheck.Count > 0 do
608 let next = pointToCheck.Pop()
609 pointChecked.Add(next) |> ignore
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