Update coding style.
[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 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