b9a31f8a5da2e9e1dd66adaaa30e59c9c741258c
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
1 module ImgTools
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 Utils
12 open Heap
13
14 // Normalize image values between 0uy and 255uy.
15 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
16 let min = ref [| 0.0 |]
17 let minLocation = ref <| [| Point() |]
18 let max = ref [| 0.0 |]
19 let maxLocation = ref <| [| Point() |]
20 img.MinMax(min, max, minLocation, maxLocation)
21 ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
22
23
24 let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
25 img.Save(filepath)
26
27
28 let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
29 use img = new Image<Gray, 'TDeph>(mat.Size)
30 mat.CopyTo(img)
31 saveImg img filepath
32
33
34 let suppressMConnections (img: Matrix<byte>) =
35 let w = img.Width
36 let h = img.Height
37 for i in 1 .. h - 2 do
38 for j in 1 .. w - 2 do
39 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)
40 then
41 img.[i, j] <- 0uy
42 for i in 1 .. h - 2 do
43 for j in 1 .. w - 2 do
44 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)
45 then
46 img.[i, j] <- 0uy
47
48
49 let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
50 let w = img.Width
51 let h = img.Height
52
53 use sobelKernel =
54 new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
55 [ 2.0f; 0.0f; -2.0f ]
56 [ 1.0f; 0.0f; -1.0f ]], Point(1, 1))
57
58 let xGradient = img.Convolution(sobelKernel).Convert<Gray, float>()
59 let yGradient = img.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
60
61 let xGradientData = xGradient.Data
62 let yGradientData = yGradient.Data
63 for r in 0 .. h - 1 do
64 xGradientData.[r, 0, 0] <- 0.0
65 xGradientData.[r, w - 1, 0] <- 0.0
66 yGradientData.[r, 0, 0] <- 0.0
67 yGradientData.[r, w - 1, 0] <- 0.0
68
69 for c in 0 .. w - 1 do
70 xGradientData.[0, c, 0] <- 0.0
71 xGradientData.[h - 1, c, 0] <- 0.0
72 yGradientData.[0, c, 0] <- 0.0
73 yGradientData.[h - 1, c, 0] <- 0.0
74
75 use magnitudes = new Matrix<float>(xGradient.Size)
76 use angles = new Matrix<float>(xGradient.Size)
77 CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles).
78
79 let thresholdHigh, thresholdLow =
80 let sensibilityHigh = 0.1
81 let sensibilityLow = 0.1
82 use magnitudesByte = magnitudes.Convert<byte>()
83 let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
84 threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
85
86 // Non-maximum suppression.
87 use nms = new Matrix<byte>(xGradient.Size)
88
89 for i in 0 .. h - 1 do
90 nms.Data.[i, 0] <- 0uy
91 nms.Data.[i, w - 1] <- 0uy
92
93 for j in 0 .. w - 1 do
94 nms.Data.[0, j] <- 0uy
95 nms.Data.[h - 1, j] <- 0uy
96
97 for i in 1 .. h - 2 do
98 for j in 1 .. w - 2 do
99 let vx = xGradient.Data.[i, j, 0]
100 let vy = yGradient.Data.[i, j, 0]
101 if vx <> 0. || vy <> 0.
102 then
103 let angle = angles.[i, j]
104
105 let vx', vy' = abs vx, abs vy
106 let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
107 let ratio1 = 1. - ratio2
108
109 let mNeigbors (sign: int) : float =
110 if angle < Math.PI / 4.
111 then ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign]
112 elif angle < Math.PI / 2.
113 then ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j]
114 elif angle < 3.0 * Math.PI / 4.
115 then ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign]
116 elif angle < Math.PI
117 then ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign]
118 elif angle < 5. * Math.PI / 4.
119 then ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign]
120 elif angle < 3. * Math.PI / 2.
121 then ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j]
122 elif angle < 7. * Math.PI / 4.
123 then ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign]
124 else ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign]
125
126 let m = magnitudes.Data.[i, j]
127 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
128 then
129 nms.Data.[i, j] <- 1uy
130
131 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
132
133 let edges = new Matrix<byte>(xGradient.Size)
134
135 // Histeresis thresholding.
136 let toVisit = Stack<Point>()
137 for i in 0 .. h - 1 do
138 for j in 0 .. w - 1 do
139 if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh
140 then
141 nms.Data.[i, j] <- 0uy
142 toVisit.Push(Point(j, i))
143 while toVisit.Count > 0 do
144 let p = toVisit.Pop()
145 edges.Data.[p.Y, p.X] <- 1uy
146 for i' in -1 .. 1 do
147 for j' in -1 .. 1 do
148 if i' <> 0 || j' <> 0
149 then
150 let ni = p.Y + i'
151 let nj = p.X + j'
152 if ni >= 0 && ni < h && nj >= 0 && nj < w && nms.Data.[ni, nj] = 1uy
153 then
154 nms.Data.[ni, nj] <- 0uy
155 toVisit.Push(Point(nj, ni))
156
157 edges, xGradient, yGradient
158
159
160 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
161 let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
162 img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
163
164
165 type Points = HashSet<Point>
166
167 let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
168 for p in points do
169 img.Data.[p.Y, p.X, 0] <- intensity
170
171 type ExtremumType =
172 | Maxima = 1
173 | Minima = 2
174
175 let findExtremum (img: Image<Gray, 'TDepth>) (extremumType: ExtremumType) : IEnumerable<Points> =
176 let w = img.Width
177 let h = img.Height
178 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
179
180 let imgData = img.Data
181 let suppress: bool[,] = Array2D.zeroCreate h w
182
183 let result = List<List<Point>>()
184
185 let flood (start: Point) : List<List<Point>> =
186 let sameLevelToCheck = Stack<Point>()
187 let betterLevelToCheck = Stack<Point>()
188 betterLevelToCheck.Push(start)
189
190 let result' = List<List<Point>>()
191
192 while betterLevelToCheck.Count > 0 do
193 let p = betterLevelToCheck.Pop()
194 if not suppress.[p.Y, p.X]
195 then
196 suppress.[p.Y, p.X] <- true
197 sameLevelToCheck.Push(p)
198 let current = List<Point>()
199
200 let mutable betterExists = false
201
202 while sameLevelToCheck.Count > 0 do
203 let p' = sameLevelToCheck.Pop()
204 let currentLevel = imgData.[p'.Y, p'.X, 0]
205 current.Add(p') |> ignore
206 for i, j in se do
207 let ni = i + p'.Y
208 let nj = j + p'.X
209 if ni >= 0 && ni < h && nj >= 0 && nj < w
210 then
211 let level = imgData.[ni, nj, 0]
212 let notSuppressed = not suppress.[ni, nj]
213
214 if level = currentLevel && notSuppressed
215 then
216 suppress.[ni, nj] <- true
217 sameLevelToCheck.Push(Point(nj, ni))
218 elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel
219 then
220 betterExists <- true
221 if notSuppressed
222 then
223 betterLevelToCheck.Push(Point(nj, ni))
224
225 if not betterExists
226 then
227 result'.Add(current)
228 result'
229
230 for i in 0 .. h - 1 do
231 for j in 0 .. w - 1 do
232 let maxima = flood (Point(j, i))
233 if maxima.Count > 0
234 then
235 result.AddRange(maxima)
236
237 result.Select(fun l -> Points(l))
238
239
240 let findMaxima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
241 findExtremum img ExtremumType.Maxima
242
243
244 let findMinima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
245 findExtremum img ExtremumType.Minima
246
247
248 type PriorityQueue () =
249 let size = 256
250 let q: Points[] = Array.init size (fun i -> Points())
251 let mutable highest = -1 // Value of the first elements of 'q'.
252 let mutable lowest = size
253
254 member this.NextMax () : byte * Point =
255 if this.IsEmpty
256 then
257 invalidOp "Queue is empty"
258 else
259 let l = q.[highest]
260 let next = l.First()
261 l.Remove(next) |> ignore
262 let value = byte highest
263
264 if l.Count = 0
265 then
266 highest <- highest - 1
267 while highest > lowest && q.[highest].Count = 0 do
268 highest <- highest - 1
269 if highest = lowest
270 then
271 highest <- -1
272 lowest <- size
273
274 value, next
275
276 member this.NextMin () : byte * Point =
277 if this.IsEmpty
278 then
279 invalidOp "Queue is empty"
280 else
281 let l = q.[lowest + 1]
282 let next = l.First()
283 l.Remove(next) |> ignore
284 let value = byte (lowest + 1)
285
286 if l.Count = 0
287 then
288 lowest <- lowest + 1
289 while lowest < highest && q.[lowest + 1].Count = 0 do
290 lowest <- lowest + 1
291 if highest = lowest
292 then
293 highest <- -1
294 lowest <- size
295
296 value, next
297
298 member this.Max =
299 highest |> byte
300
301 member this.Min =
302 lowest + 1 |> byte
303
304 member this.Add (value: byte) (p: Point) =
305 let vi = int value
306
307 if vi > highest
308 then
309 highest <- vi
310 if vi <= lowest
311 then
312 lowest <- vi - 1
313
314 q.[vi].Add(p) |> ignore
315
316 member this.Remove (value: byte) (p: Point) =
317 let vi = int value
318 if q.[vi].Remove(p) && q.[vi].Count = 0
319 then
320 if vi = highest
321 then
322 highest <- highest - 1
323 while highest > lowest && q.[highest].Count = 0 do
324 highest <- highest - 1
325 elif vi - 1 = lowest
326 then
327 lowest <- lowest + 1
328 while lowest < highest && q.[lowest + 1].Count = 0 do
329 lowest <- lowest + 1
330
331 if highest = lowest // The queue is now empty.
332 then
333 highest <- -1
334 lowest <- size
335
336 member this.IsEmpty =
337 highest = -1
338
339 member this.Clear () =
340 while highest > lowest do
341 q.[highest].Clear()
342 highest <- highest - 1
343 highest <- -1
344 lowest <- size
345
346
347 type private AreaState =
348 | Removed = 1
349 | Unprocessed = 2
350 | Validated = 3
351
352 type private AreaOperation =
353 | Opening = 1
354 | Closing = 2
355
356 [<AllowNullLiteral>]
357 type private Area (elements: Points) =
358 member this.Elements = elements
359 member val Intensity = None with get, set
360 member val State = AreaState.Unprocessed with get, set
361
362 let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperation) =
363 let w = img.Width
364 let h = img.Height
365 let imgData = img.Data
366 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
367
368 let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
369
370 let pixels: Area[,] = Array2D.create h w null
371 for m in areas do
372 for e in m.Elements do
373 pixels.[e.Y, e.X] <- m
374
375 let queue = PriorityQueue()
376
377 let addEdgeToQueue (elements: Points) =
378 for p in elements do
379 for i, j in se do
380 let ni = i + p.Y
381 let nj = j + p.X
382 let p' = Point(nj, ni)
383 if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p'))
384 then
385 queue.Add (imgData.[ni, nj, 0]) p'
386
387 // Reverse order is quicker.
388 for i in areas.Count - 1 .. -1 .. 0 do
389 let m = areas.[i]
390 if m.Elements.Count <= area && m.State <> AreaState.Removed
391 then
392 queue.Clear()
393 addEdgeToQueue m.Elements
394
395 let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min
396 let nextElements = Points()
397
398 let mutable stop = false
399 while not stop do
400 let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
401 let mutable merged = false
402
403 if intensity' = intensity // The intensity doesn't change.
404 then
405 if m.Elements.Count + nextElements.Count + 1 > area
406 then
407 m.State <- AreaState.Validated
408 m.Intensity <- Some intensity
409 stop <- true
410 else
411 nextElements.Add(p) |> ignore
412
413 elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity
414 then
415 m.Elements.UnionWith(nextElements)
416 for e in nextElements do
417 pixels.[e.Y, e.X] <- m
418
419 if m.Elements.Count = area
420 then
421 m.State <- AreaState.Validated
422 m.Intensity <- Some (intensity')
423 stop <- true
424 else
425 intensity <- intensity'
426 nextElements.Clear()
427 nextElements.Add(p) |> ignore
428
429 else
430 let m' = pixels.[p.Y, p.X]
431 if m' <> null
432 then
433 if m'.Elements.Count + m.Elements.Count <= area
434 then
435 m'.State <- AreaState.Removed
436 for e in m'.Elements do
437 pixels.[e.Y, e.X] <- m
438 queue.Remove imgData.[e.Y, e.X, 0] e
439 addEdgeToQueue m'.Elements
440 m.Elements.UnionWith(m'.Elements)
441 let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min
442 if intensityMax <> intensity
443 then
444 intensity <- intensityMax
445 nextElements.Clear()
446 merged <- true
447
448 if not merged
449 then
450 m.State <- AreaState.Validated
451 m.Intensity <- Some (intensity)
452 stop <- true
453
454 if not stop && not merged
455 then
456 for i, j in se do
457 let ni = i + p.Y
458 let nj = j + p.X
459 let p' = Point(nj, ni)
460 if ni < 0 || ni >= h || nj < 0 || nj >= w
461 then
462 m.State <- AreaState.Validated
463 m.Intensity <- Some (intensity)
464 stop <- true
465 elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
466 then
467 queue.Add (imgData.[ni, nj, 0]) p'
468
469 if queue.IsEmpty
470 then
471 if m.Elements.Count + nextElements.Count <= area
472 then
473 m.State <- AreaState.Validated
474 m.Intensity <- Some intensity'
475 m.Elements.UnionWith(nextElements)
476 stop <- true
477
478 for m in areas do
479 if m.State = AreaState.Validated
480 then
481 match m.Intensity with
482 | Some i ->
483 for p in m.Elements do
484 imgData.[p.Y, p.X, 0] <- i
485 | _ -> ()
486 ()
487
488
489 let areaOpen (img: Image<Gray, byte>) (area: int) =
490 areaOperation img area AreaOperation.Opening
491
492 let areaClose (img: Image<Gray, byte>) (area: int) =
493 areaOperation img area AreaOperation.Closing
494
495 [<AllowNullLiteral>]
496 type Island (cmp: IComparer<float32>) =
497 member val Shore = Heap.Heap<float32, Point>(cmp) with get
498 member val Level = 0.f with get, set
499 member val Surface = 0 with get, set
500
501
502 let private areaOperationF (img: Image<Gray, float32>) (area: int) (op: AreaOperation) =
503 let w = img.Width
504 let h = img.Height
505 let earth = img.Data
506 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
507
508 let comparer = if op = AreaOperation.Opening
509 then { new IComparer<float32> with member this.Compare(v1, v2) = v1.CompareTo(v2) }
510 else { new IComparer<float32> with member this.Compare(v1, v2) = v2.CompareTo(v1) }
511
512 let ownership: Island[,] = Array2D.create h w null
513
514 // Initialize islands with their shore.
515 let islands = List<Island>()
516 let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima
517 for e in extremum do
518 let island =
519 let p = e.First()
520 Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count)
521 islands.Add(island)
522 let shorePoints = Points()
523 for p in e do
524 ownership.[p.Y, p.X] <- island
525 for i, j in se do
526 let ni = i + p.Y
527 let nj = j + p.X
528 let neighbor = Point(nj, ni)
529 if ni >= 0 && ni < h && nj >= 0 && nj < w && ownership.[ni, nj] = null && not (shorePoints.Contains(neighbor))
530 then
531 shorePoints.Add(neighbor) |> ignore
532 island.Shore.Add earth.[ni, nj, 0] neighbor
533
534 for island in islands do
535 let mutable stop = island.Shore.IsEmpty
536
537 // 'true' if 'p' is owned or adjacent to 'island'.
538 let ownedOrAdjacent (p: Point) : bool =
539 ownership.[p.Y, p.X] = island ||
540 (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
541 (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
542 (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
543 (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
544
545 while not stop && island.Surface < area do
546 let level, next = island.Shore.Max
547 let other = ownership.[next.Y, next.X]
548 if other = island // During merging, some points on the shore may be owned by the island itself -> ignored.
549 then
550 island.Shore.RemoveNext ()
551 else
552 if other <> null
553 then // We touching another island.
554 if island.Surface + other.Surface >= area
555 then
556 stop <- true
557 else // We can merge 'other' into 'surface'.
558 island.Surface <- island.Surface + other.Surface
559 island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level
560 for l, p in other.Shore do
561 let mutable currentY = p.Y + 1
562 while currentY < h && ownership.[currentY, p.X] = other do
563 ownership.[currentY, p.X] <- island
564 currentY <- currentY + 1
565 island.Shore.Add l p
566 other.Shore.Clear()
567
568 elif comparer.Compare(level, island.Level) > 0
569 then
570 stop <- true
571 else
572 island.Shore.RemoveNext ()
573 for i, j in se do
574 let ni = i + next.Y
575 let nj = j + next.X
576 if ni < 0 || ni >= h || nj < 0 || nj >= w
577 then
578 island.Surface <- Int32.MaxValue
579 stop <- true
580 else
581 let neighbor = Point(nj, ni)
582 if not <| ownedOrAdjacent neighbor
583 then
584 island.Shore.Add earth.[ni, nj, 0] neighbor
585 if not stop
586 then
587 ownership.[next.Y, next.X] <- island
588 island.Level <- level
589 island.Surface <- island.Surface + 1
590
591 for i in 0 .. h - 1 do
592 for j in 0 .. w - 1 do
593 let island = ownership.[i, j]
594 if island <> null
595 then
596 earth.[i, j, 0] <- island.Level
597 ()
598
599
600 let areaOpenF (img: Image<Gray, float32>) (area: int) =
601 areaOperationF img area AreaOperation.Opening
602
603 let areaCloseF (img: Image<Gray, float32>) (area: int) =
604 areaOperationF img area AreaOperation.Closing
605
606 // A simpler algorithm than 'areaOpen' but slower.
607 let areaOpen2 (img: Image<Gray, byte>) (area: int) =
608 let w = img.Width
609 let h = img.Height
610 let imgData = img.Data
611 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
612
613 let histogram = Array.zeroCreate 256
614 for i in 0 .. h - 1 do
615 for j in 0 .. w - 1 do
616 let v = imgData.[i, j, 0] |> int
617 histogram.[v] <- histogram.[v] + 1
618
619 let flooded : bool[,] = Array2D.zeroCreate h w
620
621 let pointsChecked = HashSet<Point>()
622 let pointsToCheck = Stack<Point>()
623
624 for level in 255 .. -1 .. 0 do
625 let mutable n = histogram.[level]
626 if n > 0
627 then
628 for i in 0 .. h - 1 do
629 for j in 0 .. w - 1 do
630 if not flooded.[i, j] && imgData.[i, j, 0] = byte level
631 then
632 let mutable maxNeighborValue = 0uy
633 pointsChecked.Clear()
634 pointsToCheck.Clear()
635 pointsToCheck.Push(Point(j, i))
636
637 while pointsToCheck.Count > 0 do
638 let next = pointsToCheck.Pop()
639 pointsChecked.Add(next) |> ignore
640 flooded.[next.Y, next.X] <- true
641
642 for nx, ny in se do
643 let p = Point(next.X + nx, next.Y + ny)
644 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
645 then
646 let v = imgData.[p.Y, p.X, 0]
647 if v = byte level
648 then
649 if not (pointsChecked.Contains(p))
650 then
651 pointsToCheck.Push(p)
652 elif v > maxNeighborValue
653 then
654 maxNeighborValue <- v
655
656 if int maxNeighborValue < level && pointsChecked.Count <= area
657 then
658 for p in pointsChecked do
659 imgData.[p.Y, p.X, 0] <- maxNeighborValue
660
661
662 // Zhang and Suen algorithm.
663 // Modify 'mat' in place.
664 let thin (mat: Matrix<byte>) =
665 let w = mat.Width
666 let h = mat.Height
667 let mutable data1 = mat.Data
668 let mutable data2 = Array2D.copy data1
669
670 let mutable pixelChanged = true
671 let mutable oddIteration = true
672
673 while pixelChanged do
674 pixelChanged <- false
675 for i in 0..h-1 do
676 for j in 0..w-1 do
677 if data1.[i, j] = 1uy
678 then
679 let p2 = if i = 0 then 0uy else data1.[i-1, j]
680 let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
681 let p4 = if j = w-1 then 0uy else data1.[i, j+1]
682 let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
683 let p6 = if i = h-1 then 0uy else data1.[i+1, j]
684 let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
685 let p8 = if j = 0 then 0uy else data1.[i, j-1]
686 let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
687
688 let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
689 if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
690 (if p2 = 0uy && p3 = 1uy then 1 else 0) +
691 (if p3 = 0uy && p4 = 1uy then 1 else 0) +
692 (if p4 = 0uy && p5 = 1uy then 1 else 0) +
693 (if p5 = 0uy && p6 = 1uy then 1 else 0) +
694 (if p6 = 0uy && p7 = 1uy then 1 else 0) +
695 (if p7 = 0uy && p8 = 1uy then 1 else 0) +
696 (if p8 = 0uy && p9 = 1uy then 1 else 0) +
697 (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
698 if oddIteration
699 then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
700 else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
701 then
702 data2.[i, j] <- 0uy
703 pixelChanged <- true
704 else
705 data2.[i, j] <- 0uy
706
707 oddIteration <- not oddIteration
708 let tmp = data1
709 data1 <- data2
710 data2 <- tmp
711
712
713 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
714 // Modify 'mat' in place.
715 let removeArea (mat: Matrix<byte>) (areaSize: int) =
716 let neighbors = [|
717 (-1, 0) // p2
718 (-1, 1) // p3
719 ( 0, 1) // p4
720 ( 1, 1) // p5
721 ( 1, 0) // p6
722 ( 1, -1) // p7
723 ( 0, -1) // p8
724 (-1, -1) |] // p9
725
726 let mat' = new Matrix<byte>(mat.Size)
727 let w = mat'.Width
728 let h = mat'.Height
729 mat.CopyTo(mat')
730
731 let data = mat.Data
732 let data' = mat'.Data
733
734 for i in 0..h-1 do
735 for j in 0..w-1 do
736 if data'.[i, j] = 1uy
737 then
738 let neighborhood = List<Point>()
739 let neighborsToCheck = Stack<Point>()
740 neighborsToCheck.Push(Point(j, i))
741 data'.[i, j] <- 0uy
742
743 while neighborsToCheck.Count > 0 do
744 let n = neighborsToCheck.Pop()
745 neighborhood.Add(n)
746 for (ni, nj) in neighbors do
747 let pi = n.Y + ni
748 let pj = n.X + nj
749 if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
750 then
751 neighborsToCheck.Push(Point(pj, pi))
752 data'.[pi, pj] <- 0uy
753 if neighborhood.Count <= areaSize
754 then
755 for n in neighborhood do
756 data.[n.Y, n.X] <- 0uy
757
758 let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
759 let w = img.Width
760 let h = img.Height
761
762 let pointChecked = Points()
763 let pointToCheck = Stack<Point>(startPoints);
764
765 let data = img.Data
766
767 while pointToCheck.Count > 0 do
768 let next = pointToCheck.Pop()
769 pointChecked.Add(next) |> ignore
770 for ny in -1 .. 1 do
771 for nx in -1 .. 1 do
772 if ny <> 0 && nx <> 0
773 then
774 let p = Point(next.X + nx, next.Y + ny)
775 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
776 then
777 pointToCheck.Push(p)
778
779 List<Point>(pointChecked)
780
781
782 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
783 img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
784
785
786 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
787 img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
788
789
790 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) =
791
792 if alpha >= 1.0
793 then
794 img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
795 else
796 let windowPosX = e.Cx - e.A - 5.0
797 let gapX = windowPosX - (float (int windowPosX))
798
799 let windowPosY = e.Cy - e.A - 5.0
800 let gapY = windowPosY - (float (int windowPosY))
801
802 let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
803
804 img.ROI <- roi
805 if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
806 then
807 use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
808 i.Draw(Ellipse(PointF(float32 <| (e.A + 5. + gapX) , float32 <| (e.A + 5. + gapY)), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
809 CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
810 img.ROI <- Rectangle.Empty
811
812
813 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha: float) =
814 List.iter (fun e -> drawEllipse img e color alpha) ellipses
815
816
817 let rngCell = System.Random()
818 let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
819 if drawCellContent
820 then
821 let colorB = rngCell.Next(20, 70)
822 let colorG = rngCell.Next(20, 70)
823 let colorR = rngCell.Next(20, 70)
824
825 for y in 0 .. c.elements.Height - 1 do
826 for x in 0 .. c.elements.Width - 1 do
827 if c.elements.[y, x] > 0uy
828 then
829 let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
830 let b = img.Data.[y + dy, x + dx, 0] |> int
831 let g = img.Data.[y + dy, x + dx, 1] |> int
832 let r = img.Data.[y + dy, x + dx, 2] |> int
833 img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
834 img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
835 img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
836
837 let crossColor, crossColor2 =
838 match c.cellClass with
839 | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
840 | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
841 | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
842
843 drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
844 drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
845
846 drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
847 drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
848
849
850 let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
851 List.iter (fun c -> drawCell img drawCellContent c) cells