5 open System.Collections.Generic
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
>()
24 let gaussianFilter (img
: Image<'TColor, 'TDepth>) (standardDeviation
: float) : Image<'TColor, 'TDepth> =
25 let size = 2 * int (ceil
(4.0 * standardDeviation
)) + 1
26 img
.SmoothGaussian(size, size, standardDeviation
, standardDeviation
)
29 type Points = HashSet<Point>
31 let drawPoints (img
: Image<Gray, byte
>) (points
: Points) (intensity
: byte
) =
33 img
.Data.[p
.Y, p
.X, 0] <- intensity
40 let findExtremum (img
: Image<Gray, byte
>) (extremumType
: ExtremumType) : IEnumerable<Points> =
43 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
45 let imgData = img
.Data
46 let suppress: bool[,] = Array2D.zeroCreate
h w
48 let result = List<List<Point>>()
50 let flood (start
: Point) : List<List<Point>> =
51 let sameLevelToCheck = Stack<Point>()
52 let betterLevelToCheck = Stack<Point>()
53 betterLevelToCheck.Push(start
)
55 let result' = List<List<Point>>()
57 while betterLevelToCheck.Count > 0 do
58 let p = betterLevelToCheck.Pop()
59 if not suppress.[p.Y, p.X]
61 suppress.[p.Y, p.X] <- true
62 sameLevelToCheck.Push(p)
63 let current = List<Point>()
65 let mutable betterExists = false
67 while sameLevelToCheck.Count > 0 do
68 let p' = sameLevelToCheck.Pop()
69 let currentLevel = imgData.[p'.Y, p'.X, 0]
70 current.Add(p') |> ignore
74 if ni >= 0 && ni < h && nj >= 0 && nj < w
76 let level = imgData.[ni, nj, 0]
77 let notSuppressed = not suppress.[ni, nj]
79 if level = currentLevel && notSuppressed
81 suppress.[ni, nj] <- true
82 sameLevelToCheck.Push(Point(nj, ni))
83 elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel
88 betterLevelToCheck.Push(Point(nj, ni))
95 for i in 0 .. h - 1 do
96 for j in 0 .. w - 1 do
97 let maxima = flood (Point(j, i))
100 result.AddRange(maxima)
102 result.Select(fun l -> Points(l))
105 let findMaxima (img: Image<Gray, byte>) : IEnumerable<Points> =
106 findExtremum img ExtremumType.Maxima
108 let findMinima (img: Image<Gray, byte>) : IEnumerable<Points> =
109 findExtremum img ExtremumType.Minima
112 type PriorityQueue () =
114 let q: Points[] = Array.init size (fun i -> Points())
115 let mutable highest = -1 // Value of the first elements of 'q'.
116 let mutable lowest = size
118 member this.NextMax () : byte * Point =
121 invalidOp "Queue is empty"
125 l.Remove(next) |> ignore
126 let value = byte highest
130 highest <- highest - 1
131 while highest > lowest && q.[highest].Count = 0 do
132 highest <- highest - 1
140 member this.NextMin () : byte * Point =
143 invalidOp "Queue is empty"
145 let l = q.[lowest + 1]
147 l.Remove(next) |> ignore
148 let value = byte (lowest + 1)
153 while lowest < highest && q.[lowest + 1].Count = 0 do
168 member this.Add (value: byte) (p: Point) =
178 q.[vi].Add(p) |> ignore
180 member this.Remove (value: byte) (p: Point) =
182 if q.[vi].Remove(p) && q.[vi].Count = 0
186 highest <- highest - 1
187 while highest > lowest && q.[highest].Count = 0 do
188 highest <- highest - 1
192 while lowest < highest && q.[lowest + 1].Count = 0 do
195 if highest = lowest // The queue is now empty.
200 member this.IsEmpty =
203 member this.Clear () =
204 while highest > lowest do
206 highest <- highest - 1
211 type private AreaState =
216 type private AreaOperation =
221 type private Area (elements: Points) =
222 member this.Elements = elements
223 member val Intensity = None with get, set
224 member val State = AreaState.Unprocessed with get, set
226 let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperation) =
229 let imgData = img.Data
230 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
232 let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
234 let pixels: Area[,] = Array2D.create h w null
236 for e in m.Elements do
237 pixels.[e.Y, e.X] <- m
239 let queue = PriorityQueue()
241 let addEdgeToQueue (elements: Points) =
246 let p' = Point(nj, ni)
247 if ni >= 0 && ni < h && nj >= 0 && nj < w && not
(elements
.Contains(p'))
249 queue.Add (imgData.[ni, nj, 0]) p'
251 // Reverse order is quicker.
252 for i
in areas.Count - 1 .. -1 .. 0 do
254 if m.Elements.Count <= area
&& m.State <> AreaState.Removed
257 addEdgeToQueue m.Elements
259 let mutable intensity = if op
= AreaOperation.Opening then queue.Max else queue.Min
260 let nextElements = Points()
262 let mutable stop = false
264 let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
265 let mutable merged = false
267 if intensity' = intensity // The intensity doesn't change.
269 if m.Elements.Count + nextElements.Count + 1 > area
271 m.State <- AreaState.Validated
272 m.Intensity <- Some intensity
275 nextElements.Add(p) |> ignore
277 elif
if op
= AreaOperation.Opening then intensity' < intensity else intensity' > intensity
279 m.Elements.UnionWith(nextElements)
280 for e
in nextElements do
281 pixels.[e
.Y, e
.X] <- m
283 if m.Elements.Count = area
285 m.State <- AreaState.Validated
286 m.Intensity <- Some (intensity')
289 intensity <- intensity'
291 nextElements.Add(p) |> ignore
294 let m' = pixels.[p.Y, p.X]
297 if m'.Elements.Count + m.Elements.Count <= area
299 m'.State <- AreaState.Removed
300 for e
in m'.Elements do
301 pixels.[e.Y, e.X] <- m
302 queue.Remove imgData.[e.Y, e.X, 0] e
303 addEdgeToQueue m'.Elements
304 m.Elements.UnionWith(m'.Elements)
305 let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min
306 if intensityMax <> intensity
308 intensity <- intensityMax
314 m.State <- AreaState.Validated
315 m.Intensity <- Some (intensity)
318 if not stop && not merged
323 let p' = Point(nj, ni)
324 if ni < 0 || ni >= h || nj < 0 || nj >= w
326 m.State <- AreaState.Validated
327 m.Intensity <- Some (intensity)
329 elif
not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
331 queue.Add (imgData.[ni, nj, 0]) p'
335 if m.Elements.Count + nextElements.Count <= area
337 m.State <- AreaState.Validated
338 m.Intensity <- Some intensity'
339 m.Elements.UnionWith(nextElements)
343 if m.State = AreaState.Validated
345 match m.Intensity with
347 for p in m.Elements do
348 imgData.[p.Y, p.X, 0] <- i
353 let areaOpen (img
: Image<Gray, byte
>) (area
: int) =
354 areaOperation img area
AreaOperation.Opening
356 let areaClose (img
: Image<Gray, byte
>) (area
: int) =
357 areaOperation img area
AreaOperation.Closing
359 let areaOpen2 (img
: Image<Gray, byte
>) (area
: int) =
362 let imgData = img
.Data
363 let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
365 let histogram = Array.zeroCreate
256
366 for i
in 0 .. h - 1 do
367 for j
in 0 .. w - 1 do
368 let v = imgData.[i
, j
, 0] |> int
369 histogram.[v] <- histogram.[v] + 1
371 let flooded : bool[,] = Array2D.zeroCreate
h w
373 let pointsChecked = HashSet<Point>()
374 let pointsToCheck = Stack<Point>()
376 for level in 255 .. -1 .. 0 do
377 let mutable n = histogram.[level]
380 for i
in 0 .. h - 1 do
381 for j
in 0 .. w - 1 do
382 if not flooded.[i
, j
] && imgData.[i
, j
, 0] = byte
level
384 let mutable maxNeighborValue = 0uy
385 pointsChecked.Clear()
386 pointsToCheck.Clear()
387 pointsToCheck.Push(Point(j
, i
))
389 while pointsToCheck.Count > 0 do
390 let next = pointsToCheck.Pop()
391 pointsChecked.Add(next) |> ignore
392 flooded.[next.Y, next.X] <- true
395 let p = Point(next.X + nx
, next.Y + ny
)
396 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
398 let v = imgData.[p.Y, p.X, 0]
401 if not (pointsChecked.Contains(p))
403 pointsToCheck.Push(p)
404 elif
v > maxNeighborValue
406 maxNeighborValue <- v
408 if int maxNeighborValue < level && pointsChecked.Count <= area
410 for p in pointsChecked do
411 imgData.[p.Y, p.X, 0] <- maxNeighborValue
414 // Zhang and Suen algorithm.
415 // Modify 'mat' in place.
416 let thin (mat
: Matrix<byte
>) =
419 let mutable data1 = mat
.Data
420 let mutable data2 = Array2D.copy
data1
422 let mutable pixelChanged = true
423 let mutable oddIteration = true
425 while pixelChanged do
426 pixelChanged <- false
429 if data1.[i
, j
] = 1uy
431 let p2 = if i
= 0 then 0uy else data1.[i
-1, j
]
432 let p3 = if i
= 0 || j = w-1 then 0uy else data1.[i
-1, j+1]
433 let p4 = if j = w-1 then 0uy else data1.[i
, j+1]
434 let p5 = if i
= h-1 || j = w-1 then 0uy else data1.[i
+1, j+1]
435 let p6 = if i
= h-1 then 0uy else data1.[i
+1, j]
436 let p7 = if i
= h-1 || j = 0 then 0uy else data1.[i
+1, j-1]
437 let p8 = if j = 0 then 0uy else data1.[i
, j-1]
438 let p9 = if i
= 0 || j = 0 then 0uy else data1.[i
-1, j-1]
440 let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
441 if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
442 (if p2 = 0uy && p3 = 1uy then 1 else 0) +
443 (if p3 = 0uy && p4 = 1uy then 1 else 0) +
444 (if p4 = 0uy && p5 = 1uy then 1 else 0) +
445 (if p5 = 0uy && p6 = 1uy then 1 else 0) +
446 (if p6 = 0uy && p7 = 1uy then 1 else 0) +
447 (if p7 = 0uy && p8 = 1uy then 1 else 0) +
448 (if p8 = 0uy && p9 = 1uy then 1 else 0) +
449 (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
451 then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
452 else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
459 oddIteration <- not oddIteration
465 // FIXME: replace by a queue or stack.
466 let pop (l: List<'a>) : 'a
=
467 let n = l.[l.Count - 1]
468 l.RemoveAt(l.Count - 1)
471 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
472 // Modify 'mat' in place.
473 let removeArea (mat
: Matrix<byte
>) (areaSize
: int) =
484 let mat' = new Matrix<byte>(mat.Size)
490 let data' = mat'.Data
494 if data'.[i, j] = 1uy
496 let neighborhood = List<(int*int)>()
497 let neighborsToCheck = List<(int*int)>()
498 neighborsToCheck.Add((i, j))
501 while neighborsToCheck.Count > 0 do
502 let (ci
, cj
) = pop neighborsToCheck
503 neighborhood.Add((ci
, cj
))
504 for (ni, nj) in neighbors do
507 if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
509 neighborsToCheck.Add((pi, pj))
510 data'.[pi, pj] <- 0uy
511 if neighborhood.Count <= areaSize
513 for (ni, nj) in neighborhood do
516 let connectedComponents (img
: Image<Gray, byte
>) (startPoints
: List<Point>) : List<Point> =
520 let pointChecked = Points()
521 let pointToCheck = List<Point>(startPoints
);
525 while pointToCheck.Count > 0 do
526 let next = pop pointToCheck
527 pointChecked.Add(next) |> ignore
530 if ny
<> 0 && nx
<> 0
532 let p = Point(next.X + nx
, next.Y + ny
)
533 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
537 List<Point>(pointChecked)
540 let saveImg (img
: Image<'TColor, 'TDepth>) (filepath
: string) =
544 let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
545 use img = new Image<Gray, 'TDeph>(mat.Size)
549 let drawLine (img: Image<'TColor, 'TDepth>) (color
: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
550 img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
552 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0
: float) (y0
: float) (x1
: float) (y1
: float) (thickness
: int) =
553 img.Draw(LineSegment2DF(PointF(float32 x0
, float32 y0
), PointF(float32 x1
, float32 y1
)), color
, thickness
, CvEnum.LineType.AntiAlias);
555 let drawEllipse (img: Image<'TColor, 'TDepth>) (e
: Types.Ellipse) (color
: 'TColor) (alpha: float) =
559 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)
561 let windowPosX = e.Cx - e.A - 5.0
562 let gapX = windowPosX - (float (int windowPosX))
564 let windowPosY = e.Cy - e.A - 5.0
565 let gapY = windowPosY - (float (int windowPosY))
567 let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
570 if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
572 use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
573 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)
574 CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
575 img.ROI <- Rectangle.Empty
578 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha
: float) =
579 List.iter
(fun e -> drawEllipse img e color alpha
) ellipses
582 let rngCell = System.Random()
583 let drawCell (img: Image<Bgr, byte
>) (drawCellContent
: bool) (c
: Types.Cell) =
586 let colorB = rngCell.Next(20, 70)
587 let colorG = rngCell.Next(20, 70)
588 let colorR = rngCell.Next(20, 70)
590 for y
in 0 .. c
.elements
.Height - 1 do
591 for x
in 0 .. c
.elements
.Width - 1 do
592 if c
.elements
.[y
, x
] > 0uy
594 let dx, dy
= c
.center
.X - c
.elements
.Width / 2, c
.center
.Y - c
.elements
.Height / 2
595 let b = img.Data.[y
+ dy
, x
+ dx, 0] |> int
596 let g = img.Data.[y
+ dy
, x
+ dx, 1] |> int
597 let r = img.Data.[y
+ dy
, x
+ dx, 2] |> int
598 img.Data.[y
+ dy
, x
+ dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
599 img.Data.[y
+ dy
, x
+ dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
600 img.Data.[y
+ dy
, x
+ dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
602 let crossColor, crossColor2
=
603 match c
.cellClass
with
604 | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
605 | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
606 | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
608 drawLine img crossColor2
(c
.center
.X - 3) c
.center
.Y (c
.center
.X + 3) c
.center
.Y 2
609 drawLine img crossColor2 c
.center
.X (c
.center
.Y - 3) c
.center
.X (c
.center
.Y + 3) 2
611 drawLine img crossColor (c
.center
.X - 3) c
.center
.Y (c
.center
.X + 3) c
.center
.Y 1
612 drawLine img crossColor c
.center
.X (c
.center
.Y - 3) c
.center
.X (c
.center
.Y + 3) 1
615 let drawCells (img: Image<Bgr, byte>) (drawCellContent
: bool) (cells
: Types.Cell list) =
616 List.iter
(fun c
-> drawCell img drawCellContent
c) cells