Add functions to apply an area opening and to compute granulometry.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index f7e5a4a..e9cbbe9 100644 (file)
@@ -3,11 +3,13 @@
 open System
 open System.Drawing
 open System.Collections.Generic
+open System.Linq
 
 open Emgu.CV
 open Emgu.CV.Structure
 
 open Utils
+open Heap
 
 // Normalize image values between 0uy and 255uy.
 let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
@@ -18,10 +20,266 @@ let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
     img.MinMax(min, max, minLocation, maxLocation)
     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
 
+
 let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> =
     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 
+
+let findMaxima (img: Image<Gray, byte>) : IEnumerable<HashSet<Point>> =
+    use suppress = new Image<Gray, byte>(img.Size)
+    let w = img.Width
+    let h = img.Height
+
+    let imgData = img.Data
+    let suppressData = suppress.Data
+
+    let result = List<List<Point>>()
+
+    let flood (start: Point) : List<List<Point>> =
+        let sameLevelToCheck = Stack<Point>()
+        let betterLevelToCheck = Stack<Point>()
+        betterLevelToCheck.Push(start)
+
+        let result' = List<List<Point>>()
+
+        while betterLevelToCheck.Count > 0 do
+            let p = betterLevelToCheck.Pop()
+            if suppressData.[p.Y, p.X, 0] = 0uy
+            then
+                suppressData.[p.Y, p.X, 0] <- 1uy
+                sameLevelToCheck.Push(p)
+                let current = List<Point>()
+
+                let mutable betterExists = false
+
+                while sameLevelToCheck.Count > 0 do
+                    let p' = sameLevelToCheck.Pop()
+                    let currentLevel = imgData.[p'.Y, p'.X, 0]
+                    current.Add(p') |> ignore
+                    for i in -1 .. 1 do
+                        for j in -1 .. 1 do
+                            if i <> 0 || j <> 0
+                            then
+                                let ni = i + p'.Y
+                                let nj = j + p'.X
+                                if ni >= 0 && ni < h && nj >= 0 && nj < w
+                                then
+                                    let level = imgData.[ni, nj, 0]
+                                    let notSuppressed = suppressData.[ni, nj, 0] = 0uy
+
+                                    if level = currentLevel && notSuppressed
+                                    then
+                                        suppressData.[ni, nj, 0] <- 1uy
+                                        sameLevelToCheck.Push(Point(nj, ni))
+                                    elif level > currentLevel
+                                    then
+                                        betterExists <- true
+                                        if notSuppressed
+                                        then
+                                            betterLevelToCheck.Push(Point(nj, ni))
+
+                if not betterExists
+                then
+                    result'.Add(current)
+        result'
+
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let maxima = flood (Point(j, i))
+            if maxima.Count > 0
+            then result.AddRange(maxima)
+
+    result.Select(fun l -> HashSet<Point>(l))
+
+
+type PriorityQueue () =
+    let q = List<HashSet<Point>>() // TODO: Check performance with an HasSet
+    let mutable highest = -1 // Value of the first elements of 'q'.
+
+    member this.Next () : byte * Point =
+        if this.IsEmpty
+        then
+            invalidOp "Queue is empty"
+        else
+            let l = q.[0]
+            let next = l.First()
+            l.Remove(next) |> ignore
+            let value = byte highest
+            if l.Count = 0
+            then
+                q.RemoveAt(0)
+                highest <- highest - 1
+                while q.Count > 0 && q.[0] = null do
+                    q.RemoveAt(0)
+                    highest <- highest - 1
+            value, next
+
+    member this.Max =
+        highest |> byte
+
+    (*member this.UnionWith (other: PriorityQueue) =
+        while not other.IsEmpty do
+            let p, v = other.Next
+            this.Add p v*)
+
+    member this.Add (value: byte) (p: Point) =
+        let vi = int value
+
+        if this.IsEmpty
+        then
+            highest <- int value
+            q.Insert(0, null)
+        elif vi > highest
+        then
+            for i in highest .. vi - 1  do
+                q.Insert(0, null)
+            highest <- vi
+        elif highest - vi >= q.Count
+        then
+            for i in 0 .. highest - vi - q.Count do
+                q.Add(null)
+
+        let pos = highest - vi
+        if q.[pos] = null
+        then
+            q.[pos] <- HashSet<Point>([p])
+        else
+            q.[pos].Add(p) |> ignore
+
+
+    member this.IsEmpty =
+        q.Count = 0
+
+    member this.Clear () =
+        while highest >= 0 do
+            q.[highest] <- null
+            highest <- highest - 1
+
+
+
+type MaximaState =  Uncertain | Validated | TooBig
+type Maxima = {
+    elements : HashSet<Point>
+    mutable intensity: byte option
+    mutable state: MaximaState }
+
+
+let areaOpen (img: Image<Gray, byte>) (area: int) =
+    let w = img.Width
+    let h = img.Height
+
+    let maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq
+    let toValidated = Stack<Maxima>(maxima)
+
+    while toValidated.Count > 0 do
+        let m = toValidated.Pop()
+        if m.elements.Count <= area
+        then
+            let queue =
+                let q = PriorityQueue()
+                let firstElements = HashSet<Point>()
+                for p in m.elements do
+                    for i in -1 .. 1 do
+                        for j in -1 .. 1 do
+                            if i <> 0 || j <> 0
+                            then
+                                let ni = i + p.Y
+                                let nj = j + p.X
+                                let p' = Point(nj, ni)
+                                if ni >= 0 && ni < h && nj >= 0 && nj < w && not (m.elements.Contains(p')) && not (firstElements.Contains(p'))
+                                then
+                                    firstElements.Add(p') |> ignore
+                                    q.Add (img.Data.[ni, nj, 0]) p'
+                q
+
+            let mutable intensity = queue.Max
+            let nextElements = HashSet<Point>()
+
+            let mutable stop = false
+            while not stop do
+                let intensity', p = queue.Next ()
+
+                if intensity' = intensity // The intensity doesn't change.
+                then
+                    if m.elements.Count + nextElements.Count + 1 > area
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some intensity
+                        stop <- true
+                    else
+                        nextElements.Add(p) |> ignore
+                elif intensity' < intensity
+                then
+                    m.elements.UnionWith(nextElements)
+                    if m.elements.Count = area
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some (intensity')
+                        stop <- true
+                    else
+                        intensity <- intensity'
+                        nextElements.Clear()
+                        nextElements.Add(p) |> ignore
+                else // i' > i
+                    seq {
+                        for m' in maxima do
+                            if m' <> m && m'.elements.Contains(p) then
+                                if m'.elements.Count + m.elements.Count <= area
+                                then
+                                    m'.state <- Uncertain
+                                    m'.elements.UnionWith(m.elements)
+                                    if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list.
+                                    then
+                                        toValidated.Push(m')
+                                    stop <- true
+                                yield false
+                    } |> Seq.forall id |> ignore
+
+                    if not stop
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some (intensity)
+                        stop <- true
+
+                if not stop
+                then
+                    for i in -1 .. 1 do
+                        for j in -1 .. 1 do
+                            if i <> 0 || j <> 0
+                            then
+                                let ni = i + p.Y
+                                let nj = j + p.X
+                                let p' = Point(nj, ni)
+                                if ni < 0 || ni >= h || nj < 0 || nj >= w
+                                then
+                                    m.state <- Validated
+                                    m.intensity <- Some (intensity)
+                                    stop <- true
+                                elif not (m.elements.Contains(p')) && not (nextElements.Contains(p'))
+                                then
+                                    queue.Add (img.Data.[ni, nj, 0]) p'
+
+                if queue.IsEmpty
+                then
+                    if m.elements.Count + nextElements.Count <= area
+                    then
+                        m.state <- Validated
+                        m.intensity <- Some intensity'
+                        m.elements.UnionWith(nextElements)
+                    stop <- true
+
+    for m in maxima do
+        if m.state = Validated
+        then
+            match m.intensity with
+            | Some i ->
+                for p in m.elements do
+                    img.Data.[p.Y, p.X, 0] <- i
+            | _ -> ()
+    ()
+
+
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =
@@ -73,6 +331,7 @@ let thin (mat: Matrix<byte>) =
         data2 <- tmp
 
 
+// FIXME: replace by a queue or stack.
 let pop (l: List<'a>) : 'a =
     let n = l.[l.Count - 1]
     l.RemoveAt(l.Count - 1)
@@ -156,41 +415,37 @@ let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
     mat.CopyTo(img)
     saveImg img filepath
 
-let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) =
-    img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1);
+let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
+    img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
 
-let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) =
-    let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1)
-    drawLine img color x0 y0 x1 y1
+let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
+    img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
 
-let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
-    let cosAlpha = cos e.Alpha
-    let sinAlpha = sin e.Alpha
+let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) =
 
-    let mutable x0 = 0.0
-    let mutable y0 = 0.0
-    let mutable first_iteration = true
+    if alpha >= 1.0
+    then
+        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)
+    else
+        let windowPosX = e.Cx - e.A - 5.0
+        let gapX = windowPosX - (float (int windowPosX))
 
-    let n = 40
-    let thetaIncrement = 2.0 * Math.PI / (float n)
+        let windowPosY = e.Cy - e.A - 5.0
+        let gapY = windowPosY - (float (int windowPosY))
 
-    for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do
-        let cosTheta = cos theta
-        let sinTheta = sin theta
-        let x = e.Cx + cosAlpha * e.A * cosTheta - sinAlpha * e.B * sinTheta
-        let y = e.Cy + sinAlpha * e.A * cosTheta + cosAlpha * e.B * sinTheta
+        let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
 
-        if not first_iteration
+        img.ROI <- roi
+        if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
         then
-            drawLineF img color x0 y0 x y
-        else
-            first_iteration <- false
+            use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
+            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)
+            CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
+        img.ROI <- Rectangle.Empty
 
-        x0 <- x
-        y0 <- y
 
-let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) =
-    List.iter (fun e -> drawEllipse img e color) ellipses
+let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha: float) =
+    List.iter (fun e -> drawEllipse img e color alpha) ellipses
 
 
 let rngCell =  System.Random()
@@ -213,13 +468,18 @@ let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
                     img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
                     img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
 
-    let crossColor = match c.cellClass with
-                     | Types.HealthyRBC -> Bgr(255.0, 0.0, 0.0)
-                     | Types.InfectedRBC -> Bgr(0.0, 0.0, 255.0)
-                     | Types.Peculiar -> Bgr(0.0, 0.0, 0.0)
+    let crossColor, crossColor2 =
+        match c.cellClass with
+        | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
+        | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
+        | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
+
+    drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
+    drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
+
+    drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
+    drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
 
-    drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y
-    drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3)
 
 let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
     List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file