Version 1.0.10
[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) 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