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 let findMaxima (img
: Image<Gray, byte
>) : IEnumerable<HashSet<Point>> =
30 use suppress = new Image<Gray, byte
>(img
.Size)
34 let imgData = img
.Data
35 let suppressData = suppress.Data
37 let result = List<List<Point>>()
39 let flood (start
: Point) : List<List<Point>> =
40 let sameLevelToCheck = Stack<Point>()
41 let betterLevelToCheck = Stack<Point>()
42 betterLevelToCheck.Push(start
)
44 let result' = List<List<Point>>()
46 while betterLevelToCheck.Count > 0 do
47 let p = betterLevelToCheck.Pop()
48 if suppressData.[p.Y, p.X, 0] = 0uy
50 suppressData.[p.Y, p.X, 0] <- 1uy
51 sameLevelToCheck.Push(p)
52 let current = List<Point>()
54 let mutable betterExists = false
56 while sameLevelToCheck.Count > 0 do
57 let p' = sameLevelToCheck.Pop()
58 let currentLevel = imgData.[p'.Y, p'.X, 0]
59 current.Add(p') |> ignore
66 if ni >= 0 && ni < h && nj >= 0 && nj < w
68 let level = imgData.[ni, nj, 0]
69 let notSuppressed = suppressData.[ni, nj, 0] = 0uy
71 if level = currentLevel && notSuppressed
73 suppressData.[ni, nj, 0] <- 1uy
74 sameLevelToCheck.Push(Point(nj, ni))
75 elif level > currentLevel
80 betterLevelToCheck.Push(Point(nj, ni))
87 for i in 0 .. h - 1 do
88 for j in 0 .. w - 1 do
89 let maxima = flood (Point(j, i))
91 then result.AddRange(maxima)
93 result.Select(fun l -> HashSet<Point>(l))
96 type PriorityQueue () =
97 let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
98 let mutable highest = -1 // Value of the first elements of 'q'.
100 member this.Next () : byte * Point =
103 invalidOp "Queue is empty"
107 l.Remove(next) |> ignore
108 let value = byte highest
112 highest <- highest - 1
113 while q.Count > 0 && q.[0] = null do
115 highest <- highest - 1
121 (*member this.UnionWith (other: PriorityQueue) =
122 while not other.IsEmpty do
123 let p, v = other.Next
126 member this.Add (value: byte) (p: Point) =
135 for i in highest .. vi - 1 do
138 elif highest - vi >= q.Count
140 for i in 0 .. highest - vi - q.Count do
143 let pos = highest - vi
146 q.[pos] <- HashSet<Point>([p])
148 q.[pos].Add(p) |> ignore
151 member this.IsEmpty =
154 member this.Clear () =
155 while highest >= 0 do
157 highest <- highest - 1
161 type MaximaState = Uncertain | Validated | TooBig
163 elements : HashSet<Point>
164 mutable intensity: byte option
165 mutable state: MaximaState }
168 let areaOpen (img: Image<Gray, byte>) (area: int) =
172 let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
173 let toValidated = Stack<Maxima>(maxima)
175 while toValidated.Count > 0 do
176 let m = toValidated.Pop()
177 if m.elements.Count <= area
180 let q = PriorityQueue()
181 let firstElements = HashSet<Point>()
182 for p in m.elements do
189 let p' = Point(nj, ni)
190 if ni >= 0 && ni < h && nj >= 0 && nj < w && not (m.elements
.Contains(p')) && not (firstElements.Contains(p'))
192 firstElements.Add(p') |> ignore
193 q.Add (img.Data.[ni, nj, 0]) p'
196 let mutable intensity = queue.Max
197 let nextElements = HashSet<Point>()
199 let mutable stop = false
201 let intensity', p = queue.Next ()
203 if intensity' = intensity // The intensity doesn't change.
205 if m.elements
.Count + nextElements.Count + 1 > area
208 m.intensity <- Some intensity
211 nextElements.Add(p) |> ignore
212 elif
intensity' < intensity
214 m.elements.UnionWith(nextElements)
215 if m.elements.Count = area
218 m.intensity <- Some (intensity')
221 intensity <- intensity'
223 nextElements.Add(p) |> ignore
227 if m' <> m && m'.elements.Contains(p) then
228 if m'.elements
.Count + m.elements
.Count <= area
230 m'.state <- Uncertain
231 m'.elements
.UnionWith(m.elements
)
232 if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list.
237 } |> Seq.forall
id |> ignore
242 m.intensity <- Some (intensity)
253 let p' = Point(nj, ni)
254 if ni < 0 || ni >= h || nj < 0 || nj >= w
257 m.intensity <- Some (intensity)
259 elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
261 queue.Add (img.Data.[ni, nj, 0]) p'
265 if m.elements
.Count + nextElements.Count <= area
268 m.intensity <- Some intensity'
269 m.elements.UnionWith(nextElements)
273 if m.state = Validated
275 match m.intensity with
277 for p in m.elements do
278 img.Data.[p.Y, p.X, 0] <- i
283 // Zhang and Suen algorithm.
284 // Modify 'mat
' in place.
285 let thin (mat: Matrix<byte>) =
288 let mutable data1 = mat.Data
289 let mutable data2 = Array2D.copy data1
291 let mutable pixelChanged = true
292 let mutable oddIteration = true
294 while pixelChanged do
295 pixelChanged <- false
298 if data1.[i, j] = 1uy
300 let p2 = if i = 0 then 0uy else data1.[i-1, j]
301 let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
302 let p4 = if j = w-1 then 0uy else data1.[i, j+1]
303 let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
304 let p6 = if i = h-1 then 0uy else data1.[i+1, j]
305 let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
306 let p8 = if j = 0 then 0uy else data1.[i, j-1]
307 let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
309 let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
310 if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
311 (if p2 = 0uy && p3 = 1uy then 1 else 0) +
312 (if p3 = 0uy && p4 = 1uy then 1 else 0) +
313 (if p4 = 0uy && p5 = 1uy then 1 else 0) +
314 (if p5 = 0uy && p6 = 1uy then 1 else 0) +
315 (if p6 = 0uy && p7 = 1uy then 1 else 0) +
316 (if p7 = 0uy && p8 = 1uy then 1 else 0) +
317 (if p8 = 0uy && p9 = 1uy then 1 else 0) +
318 (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
320 then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
321 else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
328 oddIteration <- not oddIteration
334 // FIXME: replace by a queue or stack.
335 let pop (l: List<'a
>) : 'a =
336 let n = l.[l.Count - 1]
337 l.RemoveAt(l.Count - 1)
340 // Remove all 8-connected pixels with an area equal or greater than 'areaSize
'.
341 // Modify 'mat
' in place.
342 let removeArea (mat: Matrix<byte>) (areaSize: int) =
353 let mat' = new Matrix<byte
>(mat.Size)
359 let data' = mat'.Data
363 if data'.[i
, j] = 1uy
365 let neighborhood = List<(int*int)>()
366 let neighborsToCheck = List<(int*int)>()
367 neighborsToCheck.Add((i
, j))
370 while neighborsToCheck.Count > 0 do
371 let (ci, cj) = pop neighborsToCheck
372 neighborhood.Add((ci, cj))
373 for (ni, nj) in neighbors do
376 if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
378 neighborsToCheck.Add((pi, pj))
379 data'.[pi, pj] <- 0uy
380 if neighborhood.Count <= areaSize
382 for (ni, nj) in neighborhood do
385 let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
389 let pointChecked = HashSet<Point>()
390 let pointToCheck = List<Point>(startPoints);
394 while pointToCheck.Count > 0 do
395 let next = pop pointToCheck
396 pointChecked.Add(next) |> ignore
399 if ny <> 0 && nx <> 0
401 let p = Point(next.X + nx, next.Y + ny)
402 if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
406 List<Point>(pointChecked)
409 let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
413 let saveMat (mat: Matrix<'TDepth>) (filepath
: string) =
414 use img = new Image<Gray, 'TDeph>(mat.Size)
418 let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0
: int) (y0
: int) (x1
: int) (y1
: int) (thickness
: int) =
419 img.Draw(LineSegment2D(Point(x0
, y0
), Point(x1
, y1
)), color
, thickness
);
421 let drawLineF (img: Image<'TColor, 'TDepth>) (color
: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
422 img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
424 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha
: float) =
428 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)
430 let windowPosX = e.Cx - e.A - 5.0
431 let gapX = windowPosX - (float (int windowPosX))
433 let windowPosY = e.Cy - e.A - 5.0
434 let gapY = windowPosY - (float (int windowPosY))
436 let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
439 if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
441 use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
442 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)
443 CvInvoke.AddWeighted(img, 1.0, i, alpha
, 0.0, img)
444 img.ROI <- Rectangle.Empty
447 let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses
: Types.Ellipse list) (color
: 'TColor) (alpha: float) =
448 List.iter (fun e -> drawEllipse img e color alpha) ellipses
451 let rngCell = System.Random()
452 let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
455 let colorB = rngCell.Next(20, 70)
456 let colorG = rngCell.Next(20, 70)
457 let colorR = rngCell.Next(20, 70)
459 for y in 0 .. c.elements.Height - 1 do
460 for x in 0 .. c.elements.Width - 1 do
461 if c.elements.[y, x] > 0uy
463 let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
464 let b = img.Data.[y + dy, x + dx, 0] |> int
465 let g = img.Data.[y + dy, x + dx, 1] |> int
466 let r = img.Data.[y + dy, x + dx, 2] |> int
467 img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
468 img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
469 img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
471 let crossColor, crossColor2 =
472 match c.cellClass with
473 | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
474 | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
475 | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
477 drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
478 drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
480 drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
481 drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
484 let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
485 List.iter (fun c -> drawCell img drawCellContent c) cells