Add documentation.
authorGreg Burri <greg.burri@gmail.com>
Wed, 20 Jan 2016 08:51:17 +0000 (09:51 +0100)
committerGreg Burri <greg.burri@gmail.com>
Wed, 20 Jan 2016 08:51:17 +0000 (09:51 +0100)
Parasitemia/ParasitemiaCore/Classifier.fs
Parasitemia/ParasitemiaCore/Granulometry.fs
Parasitemia/ParasitemiaCore/Heap.fs
Parasitemia/ParasitemiaCore/ImgTools.fs
Parasitemia/ParasitemiaCore/MatchingEllipses.fs
Parasitemia/ParasitemiaCore/Types.fs
Parasitemia/ParasitemiaCore/Utils.fs

index 479ba59..c523790 100644 (file)
@@ -57,7 +57,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
                     let c' = PointF(e'.Cx, e'.Cy)
                     let v = pointFromTwoLines d1 (lineFromTwoPoints c c')
                     let case1 = sign (v.X - c.X) <> sign (v.X - c'.X) || Utils.squaredDistanceTwoPoints v c > Utils.squaredDistanceTwoPoints v c'
-                    if d2.Valid
+                    if not (Single.IsInfinity d2.A)
                     then
                         let p' = Utils.pointFromTwoLines d1 d2
                         let delta, delta' =
index ed8bda6..0afe5cc 100644 (file)
@@ -8,9 +8,13 @@ open Emgu.CV.Structure
 
 open Utils
 
-// 'range': a minimum and maximum radius.
-// 'scale': <= 1.0, to speed up the process.
-let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: float) : int =
+/// <summary>
+/// Granulometry by closing on the image 'img' by testing circle structrual elements with radius range in 'range'.
+/// </summary>
+/// <param name="img"></param>
+/// <param name="range">Minimum radius * maximum radius</param>
+/// <param name="scale">le 1.0, to speed up the process.</param>
+let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: float) (useOctagon: bool) : int =
     use scaledImg = if scale = 1. then img else img.Resize(scale, CvEnum.Inter.Area)
 
     let r1, r2 = range
@@ -20,7 +24,7 @@ let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: f
     let intensityImg = scaledImg.GetSum().Intensity
 
     // 's' must be odd.
-    let octagon (s: int) : Matrix<byte> =
+    let octagon (s: int) : Mat =
         if s % 2 = 0 then failwith "s must be odd"
         let m = new Matrix<byte>(Array2D.create s s 1uy)
         let r = (float s) / (Math.Sqrt 2. + 2.) |> roundInt
@@ -32,12 +36,13 @@ let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: f
                     m.[s - i - 1, j] <- 0uy
                     m.[i, s - j - 1] <- 0uy
                     m.[s - i - 1, s - j - 1] <- 0uy
-        m
+        m.Mat
 
     let mutable previous_n = Double.NaN
     for r in r1' .. r2' do
-        let se = CvInvoke.GetStructuringElement(CvEnum.ElementShape.Ellipse, Size(2 * r, 2 * r), Point(-1, -1))
-        //let se = octagon (2 * r - 1)
+        let se = if useOctagon
+                 then octagon (2 * r - 1) // It doesnd't speed up the process.
+                 else CvInvoke.GetStructuringElement(CvEnum.ElementShape.Ellipse, Size(2 * r, 2 * r), Point(-1, -1))
 
         use closed = scaledImg.MorphologyEx(CvEnum.MorphOp.Close, se, Point(-1, -1), 1, CvEnum.BorderType.Replicate, MCvScalar(0.0))
 
@@ -52,8 +57,14 @@ let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: f
 
     float (max + r1') / scale |> roundInt
 
-let findRadiusByAreaClosing (img: Image<Gray, float32>) (range: int * int) : int =
-    let r1, r2 = range
+/// <summary>
+/// Granulometry by area closing on the image 'img' by testing the circle area into the given radius range.
+/// </summary>
+let findRadiusByAreaClosing (img: Image<Gray, float32>) (radiusRange: int * int) : int =
+    let r1, r2 = radiusRange
+
+    if r1 > r2
+    then failwithf "'radiusRange' invalid : %A" radiusRange
 
     use imgCopy = img.Copy()
 
@@ -64,6 +75,6 @@ let findRadiusByAreaClosing (img: Image<Gray, float32>) (range: int * int) : int
         if r <> r1 && diff > maxDiff
         then
             maxDiff <- diff
-            max_r <- r - 1 )
+            max_r <- r - 1)
     max_r
 
index c23cdbb..32887a2 100644 (file)
@@ -14,8 +14,8 @@ type private Node<'k, 'v> =
     override this.ToString () = sprintf "%A -> %A" this.key this.value
 
 /// <summary>
-/// An heap min or max depending of the comparer.
-/// The goal is to have a set of data and be able to get the value associated with the min (or max) key.
+/// An heap min or max depending of the given key comparer.
+/// The goal is to have a set of data and be able to get the value associated with the min (or max) key value.
 /// </summary>
 type Heap<'k, 'v> (kComparer : IComparer<'k>) =
     let a = List<Node<'k, 'v>>()
@@ -98,5 +98,3 @@ type Heap<'k, 'v> (kComparer : IComparer<'k>) =
         max.key, max.value
 
     member this.Clear () = a.Clear()
-
-
index abc8714..cc15af7 100644 (file)
@@ -10,6 +10,7 @@ open Emgu.CV.Structure
 
 open Heap
 open Const
+open Types
 open Utils
 
 // Normalize image values between 0uy and 255uy.
@@ -145,7 +146,10 @@ let otsu (hist: Histogram) : float32 * float32 * float32 =
 
     toFloat level, toFloat mean1, toFloat mean2
 
-let suppressMConnections (img: Matrix<byte>) =
+/// <summary>
+/// Remove M-adjacent pixels. It may be used after thinning.
+/// </summary>
+let suppressMAdjacency (img: Matrix<byte>) =
     let w = img.Width
     let h = img.Height
     for i in 1 .. h - 2 do
@@ -159,6 +163,11 @@ let suppressMConnections (img: Matrix<byte>) =
             then
                 img.[i, j] <- 0uy
 
+/// <summary>
+/// Find edges of an image by using the Canny approach.
+/// The thresholds are automatically defined with otsu on gradient magnitudes.
+/// </summary>
+/// <param name="img"></param>
 let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32> * Image<Gray, float32> =
     let w = img.Width
     let h = img.Height
@@ -187,7 +196,7 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
 
     use magnitudes = new Matrix<float32>(xGradient.Size)
     use angles = new Matrix<float32>(xGradient.Size)
-    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles).
+    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles.
 
     let thresholdHigh, thresholdLow =
         let sensibilityHigh = 0.1f
@@ -205,8 +214,6 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
     let xGradientData = xGradient.Data
     let yGradientData = yGradient.Data
 
-    let PI = float32 Math.PI
-
     for i in 0 .. h - 1 do
         nmsData.[i, 0] <- 0uy
         nmsData.[i, w - 1] <- 0uy
@@ -282,8 +289,6 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) :
     let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
     img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
 
-type Points = HashSet<Point>
-
 let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
     for p in points do
         img.Data.[p.Y, p.X, 0] <- intensity
@@ -601,12 +606,73 @@ let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperatio
             | _ -> ()
     ()
 
+/// <summary>
+/// Area opening on byte image.
+/// </summary>
 let areaOpen (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Opening
 
+/// <summary>
+/// Area closing on byte image.
+/// </summary>
 let areaClose (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Closing
 
+// A simpler algorithm than 'areaOpen' on byte image but slower.
+let areaOpen2 (img: Image<Gray, byte>) (area: int) =
+    let w = img.Width
+    let h = img.Height
+    let imgData = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let histogram = Array.zeroCreate 256
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let v = imgData.[i, j, 0] |> int
+            histogram.[v] <- histogram.[v] + 1
+
+    let flooded : bool[,] = Array2D.zeroCreate h w
+
+    let pointsChecked = HashSet<Point>()
+    let pointsToCheck = Stack<Point>()
+
+    for level in 255 .. -1 .. 0 do
+        let mutable n = histogram.[level]
+        if n > 0
+        then
+            for i in 0 .. h - 1 do
+                for j in 0 .. w - 1 do
+                    if not flooded.[i, j] && imgData.[i, j, 0] = byte level
+                    then
+                        let mutable maxNeighborValue = 0uy
+                        pointsChecked.Clear()
+                        pointsToCheck.Clear()
+                        pointsToCheck.Push(Point(j, i))
+
+                        while pointsToCheck.Count > 0 do
+                            let next = pointsToCheck.Pop()
+                            pointsChecked.Add(next) |> ignore
+                            flooded.[next.Y, next.X] <- true
+
+                            for nx, ny in se do
+                                let p = Point(next.X + nx, next.Y + ny)
+                                if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
+                                then
+                                    let v = imgData.[p.Y, p.X, 0]
+                                    if v = byte level
+                                    then
+                                        if not (pointsChecked.Contains(p))
+                                        then
+                                            pointsToCheck.Push(p)
+                                    elif v > maxNeighborValue
+                                    then
+                                        maxNeighborValue <- v
+
+                        if int maxNeighborValue < level && pointsChecked.Count <= area
+                        then
+                            for p in pointsChecked do
+                                imgData.[p.Y, p.X, 0] <- maxNeighborValue
+
 [<AllowNullLiteral>]
 type Island (cmp: IComparer<float32>) =
     member val Shore = Heap.Heap<float32, Point>(cmp) with get
@@ -720,75 +786,36 @@ let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list)
         | _ -> ()
     ()
 
+/// <summary>
+/// Area opening on float image.
+/// </summary>
 let areaOpenF (img: Image<Gray, float32>) (area: int) =
     areaOperationF img [ area, () ] None AreaOperation.Opening
 
+/// <summary>
+/// Area closing on float image.
+/// </summary>
 let areaCloseF (img: Image<Gray, float32>) (area: int) =
     areaOperationF img [ area, () ] None AreaOperation.Closing
 
+/// <summary>
+/// Area closing on float image with different areas. Given areas must be sorted increasingly.
+/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference
+/// Between the previous and the current closing.
+/// </summary>
 let areaOpenFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
     areaOperationF img areas (Some f) AreaOperation.Opening
 
+/// <summary>
+/// Same as 'areaOpenFWithFun' for closing operation.
+/// </summary>
 let areaCloseFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
     areaOperationF img areas (Some f) AreaOperation.Closing
 
-// A simpler algorithm than 'areaOpen' but slower.
-let areaOpen2 (img: Image<Gray, byte>) (area: int) =
-    let w = img.Width
-    let h = img.Height
-    let imgData = img.Data
-    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
-
-    let histogram = Array.zeroCreate 256
-    for i in 0 .. h - 1 do
-        for j in 0 .. w - 1 do
-            let v = imgData.[i, j, 0] |> int
-            histogram.[v] <- histogram.[v] + 1
-
-    let flooded : bool[,] = Array2D.zeroCreate h w
-
-    let pointsChecked = HashSet<Point>()
-    let pointsToCheck = Stack<Point>()
-
-    for level in 255 .. -1 .. 0 do
-        let mutable n = histogram.[level]
-        if n > 0
-        then
-            for i in 0 .. h - 1 do
-                for j in 0 .. w - 1 do
-                    if not flooded.[i, j] && imgData.[i, j, 0] = byte level
-                    then
-                        let mutable maxNeighborValue = 0uy
-                        pointsChecked.Clear()
-                        pointsToCheck.Clear()
-                        pointsToCheck.Push(Point(j, i))
-
-                        while pointsToCheck.Count > 0 do
-                            let next = pointsToCheck.Pop()
-                            pointsChecked.Add(next) |> ignore
-                            flooded.[next.Y, next.X] <- true
-
-                            for nx, ny in se do
-                                let p = Point(next.X + nx, next.Y + ny)
-                                if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
-                                then
-                                    let v = imgData.[p.Y, p.X, 0]
-                                    if v = byte level
-                                    then
-                                        if not (pointsChecked.Contains(p))
-                                        then
-                                            pointsToCheck.Push(p)
-                                    elif v > maxNeighborValue
-                                    then
-                                        maxNeighborValue <- v
-
-                        if int maxNeighborValue < level && pointsChecked.Count <= area
-                        then
-                            for p in pointsChecked do
-                                imgData.[p.Y, p.X, 0] <- maxNeighborValue
-
-// Zhang and Suen algorithm.
-// Modify 'mat' in place.
+/// <summary>
+/// Zhang and Suen thinning algorithm.
+/// Modify 'mat' in place.
+/// </summary>
 let thin (mat: Matrix<byte>) =
     let w = mat.Width
     let h = mat.Height
@@ -837,8 +864,10 @@ let thin (mat: Matrix<byte>) =
         data1 <- data2
         data2 <- tmp
 
-// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
-// Modify 'mat' in place.
+/// <summary>
+/// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
+/// Modify 'mat' in place.
+/// </summary>
 let removeArea (mat: Matrix<byte>) (areaSize: int) =
     let neighbors = [|
         (-1,  0) // p2
@@ -882,7 +911,7 @@ let removeArea (mat: Matrix<byte>) (areaSize: int) =
                     for n in neighborhood do
                         data.[n.Y, n.X] <- 0uy
 
-let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
+let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Points =
     let w = img.Width
     let h = img.Height
 
@@ -903,7 +932,7 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
                     then
                         pointToCheck.Push(p)
 
-    List<Point>(pointChecked)
+    pointChecked
 
 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);
@@ -911,10 +940,10 @@ let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int)
 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) (alpha: float) =
+let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Ellipse) (color: 'TColor) (alpha: float) =
     if alpha >= 1.0
     then
-        img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+        img.Draw(Emgu.CV.Structure.Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
     else
         let windowPosX = e.Cx - e.A - 5.f
         let gapX = windowPosX - (float32 (int windowPosX))
@@ -928,15 +957,15 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo
         if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
         then
             use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
-            i.Draw(Ellipse(PointF(float32 <| (e.A + 5.f + gapX) , float32 <| (e.A + 5.f + gapY)), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+            i.Draw(Emgu.CV.Structure.Ellipse(PointF(float32 <| (e.A + 5.f + gapX) , float32 <| (e.A + 5.f + gapY)), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
             CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
         img.ROI <- Rectangle.Empty
 
-let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha: float) =
+let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Ellipse list) (color: 'TColor) (alpha: float) =
     List.iter (fun e -> drawEllipse img e color alpha) ellipses
 
 let rngCell =  System.Random()
-let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
+let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Cell) =
     if drawCellContent
     then
         let colorB = rngCell.Next(20, 70)
@@ -957,9 +986,9 @@ let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
 
     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.)
+        | HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
+        | InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
+        | 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
@@ -968,5 +997,5 @@ let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
     drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
 
 
-let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
+let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Cell list) =
     List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file
index 65ae3e9..4cc3d2b 100644 (file)
@@ -1,6 +1,7 @@
 module ParasitemiaCore.MatchingEllipses
 
 open System
+open System.Drawing
 open System.Linq
 open System.Collections
 open System.Collections.Generic
index 8a3bd4f..917e7d8 100644 (file)
@@ -2,12 +2,15 @@
 
 open System
 open System.Drawing
+open System.Collections.Generic
 
 open Emgu.CV
 open Emgu.CV.Structure
 
 open Const
 
+type Points = HashSet<Point>
+
 type Ellipse (cx: float32, cy: float32, a: float32, b: float32, alpha: float32) =
     member this.Cx = cx
     member this.Cy = cy
@@ -16,7 +19,7 @@ type Ellipse (cx: float32, cy: float32, a: float32, b: float32, alpha: float32)
     member this.Alpha = alpha
     member this.Area = a * b * PI
 
-    // Does the ellipse contain the point (x, y)?.
+    // Does the ellipse contain the point (x, y)?
     member this.Contains x y =
         ((x - cx) * cos alpha + (y - cy) * sin alpha) ** 2.f / a ** 2.f + ((x - cx) * sin alpha - (y - cy) * cos alpha) ** 2.f / b ** 2.f <= 1.f
 
@@ -55,10 +58,5 @@ type Cell = {
 type Line (a: float32, b: float32) =
     member this.A = a
     member this.B = b
-    member this.Valid = not (Single.IsInfinity this.A)
 
-[<Struct>]
-type PointF (x: float32, y: float32) =
-    member this.X = x
-    member this.Y = y
 
index 2af692d..d83bf8d 100644 (file)
@@ -1,6 +1,7 @@
 module ParasitemiaCore.Utils
 
 open System
+open System.Drawing
 
 open Types