588622e2c0fd216d4ffec7a6f65d687784f87a9d
[master-thesis.git] / 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