24c6998ed97cd505204d4b253849d3f35dc56456
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
35 let findMaxima (img
: Image<Gray, byte
>) : IEnumerable<Points> =
36 use suppress = new Image<Gray, byte
>(img
.Size)
40 let imgData = img
.Data
41 let suppressData = suppress.Data
43 let result = List<List<Point>>()
45 let flood (start
: Point) : List<List<Point>> =
46 let sameLevelToCheck = Stack<Point>()
47 let betterLevelToCheck = Stack<Point>()
48 betterLevelToCheck.Push(start
)
50 let result' = List<List<Point>>()
52 while betterLevelToCheck.Count > 0 do
53 let p = betterLevelToCheck.Pop()
54 if suppressData.[p.Y, p.X, 0] = 0uy
56 suppressData.[p.Y, p.X, 0] <- 1uy
57 sameLevelToCheck.Push(p)
58 let current = List<Point>()
60 let mutable betterExists = false
62 while sameLevelToCheck.Count > 0 do
63 let p' = sameLevelToCheck.Pop()
64 let currentLevel = imgData.[p'.Y, p'.X, 0]
65 current.Add(p') |> ignore
72 if ni >= 0 && ni < h && nj >= 0 && nj < w
74 let level = imgData.[ni, nj, 0]
75 let notSuppressed = suppressData.[ni, nj, 0] = 0uy
77 if level = currentLevel && notSuppressed
79 suppressData.[ni, nj, 0] <- 1uy
80 sameLevelToCheck.Push(Point(nj, ni))
81 elif level > currentLevel
86 betterLevelToCheck.Push(Point(nj, ni))
93 for i in 0 .. h - 1 do
94 for j in 0 .. w - 1 do
95 let maxima = flood (Point(j, i))
97 then result.AddRange(maxima)
99 result.Select(fun l -> Points(l))
102 type PriorityQueue () =
104 let q: Points[] = Array.init size (fun i -> Points()) // TODO: Check performance with an HasSet
105 let mutable highest = -1 // Value of the first elements of 'q'.
106 let mutable lowest = size
108 member this.Next () : byte * Point =
111 invalidOp "Queue is empty"
115 l.Remove(next) |> ignore
116 let value = byte highest
120 highest <- highest - 1
121 while highest > lowest && q.[highest].Count = 0 do
122 highest <- highest - 1
133 member this.Add (value: byte) (p: Point) =
143 q.[vi].Add(p) |> ignore
145 member this.Remove (value: byte) (p: Point) =
147 if q.[vi].Remove(p) && q.[vi].Count = 0
151 highest <- highest - 1
152 while highest > lowest && q.[highest].Count = 0 do
153 highest <- highest - 1
157 while lowest < highest && q.[lowest + 1].Count = 0 do
160 if highest = lowest // The queue is now empty.
165 member this.IsEmpty =
168 member this.Clear () =
169 while highest > lowest do
171 highest <- highest - 1
182 type Area (elements: Points) =
183 member this.Elements = elements
184 member val Intensity = None with get, set
185 member val State = AreaState.Unprocessed with get, set
187 let areaOpen (img: Image<Gray, byte>) (area: int) =
190 let imgData = img.Data
192 let areas = List<Area>(findMaxima img |> Seq.map Area)
194 let pixels: Area[,] = Array2D.create h w null
196 for e in m.Elements do
197 pixels.[e.Y, e.X] <- m
199 let queue = PriorityQueue()
201 let addEdgeToQueue (elements: Points) =
209 let p' = Point(nj, ni)
210 if ni >= 0 && ni < h && nj >= 0 && nj < w && not
(elements
.Contains(p'))
212 queue.Add (imgData.[ni, nj, 0]) p'
214 // Reverse order is quicker.
215 for i
in areas.Count - 1 .. -1 .. 0 do
217 if m.Elements.Count <= area
&& m.State <> AreaState.Removed
220 addEdgeToQueue m.Elements
222 let mutable intensity = queue.Max
223 let nextElements = Points()
225 let mutable stop = false
227 let intensity', p = queue.Next ()
228 let mutable merged = false
230 if intensity' = intensity // The intensity doesn't change.
232 if m.Elements.Count + nextElements.Count + 1 > area
234 m.State <- AreaState.Validated
235 m.Intensity <- Some intensity
238 nextElements.Add(p) |> ignore
240 elif
intensity' < intensity
242 m.Elements.UnionWith(nextElements)
243 for e in nextElements do
244 pixels.[e.Y, e.X] <- m
246 if m.Elements.Count = area
248 m.State <- AreaState.Validated
249 m.Intensity <- Some (intensity')
252 intensity <- intensity'
254 nextElements.Add(p) |> ignore
257 let m' = pixels.[p.Y, p.X]
260 if m'.Elements.Count + m.Elements.Count <= area
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 = queue.Max
269 if intensityMax <> intensity
271 intensity <- intensityMax
277 m.State <- AreaState.Validated
278 m.Intensity <- Some (intensity)
281 if not stop && not merged
289 let p' = Point(nj, ni)
290 if ni < 0 || ni >= h || nj < 0 || nj >= w
292 m.State <- AreaState.Validated
293 m.Intensity <- Some (intensity)
295 elif
not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
297 queue.Add (imgData.[ni, nj, 0]) p'
301 if m.Elements.Count + nextElements.Count <= area
303 m.State <- AreaState.Validated
304 m.Intensity <- Some intensity'
305 m.Elements.UnionWith(nextElements)
309 if m.State = AreaState.Validated
311 match m.Intensity with
313 for p in m.Elements do
314 imgData.[p.Y, p.X, 0] <- i
319 // Zhang and Suen algorithm.
320 // Modify 'mat' in place.
321 let thin (mat
: Matrix<byte
>) =
324 let mutable data1 = mat
.Data
325 let mutable data2 = Array2D.copy
data1
327 let mutable pixelChanged = true
328 let mutable oddIteration = true
330 while pixelChanged do
331 pixelChanged <- false
334 if data1.[i
, j] = 1uy
336 let p2 = if i
= 0 then 0uy else data1.[i
-1, j]
337 let p3 = if i
= 0 || j = w-1 then 0uy else data1.[i
-1, j+1]
338 let p4 = if j = w-1 then 0uy else data1.[i
, j+1]
339 let p5 = if i
= h-1 || j = w-1 then 0uy else data1.[i
+1, j+1]
340 let p6 = if i
= h-1 then 0uy else data1.[i
+1, j]
341 let p7 = if i
= h-1 || j = 0 then 0uy else data1.[i
+1, j-1]
342 let p8 = if j = 0 then 0uy else data1.[i
, j-1]
343 let p9 = if i
= 0 || j = 0 then 0uy else data1.[i
-1, j-1]
345 let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
346 if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
347 (if p2 = 0uy && p3 = 1uy then 1 else 0) +
348 (if p3 = 0uy && p4 = 1uy then 1 else 0) +
349 (if p4 = 0uy && p5 = 1uy then 1 else 0) +
350 (if p5 = 0uy && p6 = 1uy then 1 else 0) +
351 (if p6 = 0uy && p7 = 1uy then 1 else 0) +
352 (if p7 = 0uy && p8 = 1uy then 1 else 0) +
353 (if p8 = 0uy && p9 = 1uy then 1 else 0) +
354 (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
356 then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
357 else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
364 oddIteration <- not oddIteration
370 // FIXME: replace by a queue or stack.
371 let pop (l: List<'a>) : 'a
=
372 let n = l.[l.Count - 1]
373 l.RemoveAt(l.Count - 1)
376 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
377 // Modify 'mat' in place.
378 let removeArea (mat
: Matrix<byte
>) (areaSize
: int) =
389 let mat' = new Matrix<byte>(mat.Size)
395 let data' = mat'.Data
399 if data'.[i, j] = 1uy
401 let neighborhood = List<(int*int)>()
402 let neighborsToCheck = List<(int*int)>()
403 neighborsToCheck.Add((i, j))
406 while neighborsToCheck.Count > 0 do
407 let (ci
, cj
) = pop neighborsToCheck
408 neighborhood.Add((ci
, cj
))
409 for (ni, nj) in neighbors do
412 if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
414 neighborsToCheck.Add((pi, pj))
415 data'.[pi, pj] <- 0uy
416 if neighborhood.Count <= areaSize
418 for (ni, nj) in neighborhood do
421 let connectedComponents (img
: Image<Gray, byte
>) (startPoints
: List<Point>) : List<Point> =
425 let pointChecked = Points()
426 let pointToCheck = List<Point>(startPoints
);
430 while pointToCheck.Count > 0 do
431 let next = pop pointToCheck
432 pointChecked.Add(next) |> ignore
435 if ny
<> 0 && nx
<> 0
437 let p = Point(next.X + nx
, next.Y + ny
)
438 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
442 List<Point>(pointChecked)
445 let saveImg (img
: Image<'TColor, 'TDepth>) (filepath
: string) =
449 let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
450 use img = new Image<Gray, 'TDeph>(mat.Size)
454 let drawLine (img: Image<'TColor, 'TDepth>) (color
: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
455 img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
457 let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0
: float) (y0
: float) (x1
: float) (y1
: float) (thickness
: int) =
458 img.Draw(LineSegment2DF(PointF(float32 x0
, float32 y0
), PointF(float32 x1
, float32 y1
)), color
, thickness
, CvEnum.LineType.AntiAlias);
460 let drawEllipse (img: Image<'TColor, 'TDepth>) (e
: Types.Ellipse) (color
: 'TColor) (alpha: float) =
464 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)
466 let windowPosX = e.Cx - e.A - 5.0
467 let gapX = windowPosX - (float (int windowPosX))
469 let windowPosY = e.Cy - e.A - 5.0
470 let gapY = windowPosY - (float (int windowPosY))
472 let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
475 if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
477 use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
478 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)
479 CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
480 img.ROI <- Rectangle.Empty
483 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha
: float) =
484 List.iter
(fun e -> drawEllipse img e color alpha
) ellipses
487 let rngCell = System.Random()
488 let drawCell (img: Image<Bgr, byte
>) (drawCellContent
: bool) (c
: Types.Cell) =
491 let colorB = rngCell.Next(20, 70)
492 let colorG = rngCell.Next(20, 70)
493 let colorR = rngCell.Next(20, 70)
495 for y
in 0 .. c
.elements
.Height - 1 do
496 for x
in 0 .. c
.elements
.Width - 1 do
497 if c
.elements
.[y
, x
] > 0uy
499 let dx, dy
= c
.center
.X - c
.elements
.Width / 2, c
.center
.Y - c
.elements
.Height / 2
500 let b = img.Data.[y
+ dy
, x
+ dx, 0] |> int
501 let g = img.Data.[y
+ dy
, x
+ dx, 1] |> int
502 let r = img.Data.[y
+ dy
, x
+ dx, 2] |> int
503 img.Data.[y
+ dy
, x
+ dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
504 img.Data.[y
+ dy
, x
+ dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
505 img.Data.[y
+ dy
, x
+ dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
507 let crossColor, crossColor2
=
508 match c
.cellClass
with
509 | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
510 | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
511 | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
513 drawLine img crossColor2
(c
.center
.X - 3) c
.center
.Y (c
.center
.X + 3) c
.center
.Y 2
514 drawLine img crossColor2 c
.center
.X (c
.center
.Y - 3) c
.center
.X (c
.center
.Y + 3) 2
516 drawLine img crossColor (c
.center
.X - 3) c
.center
.Y (c
.center
.X + 3) c
.center
.Y 1
517 drawLine img crossColor c
.center
.X (c
.center
.Y - 3) c
.center
.X (c
.center
.Y + 3) 1
520 let drawCells (img: Image<Bgr, byte>) (drawCellContent
: bool) (cells
: Types.Cell list) =
521 List.iter
(fun c
-> drawCell img drawCellContent
c) cells