Add documentation.
[master-thesis.git] / Parasitemia / ParasitemiaCore / ImgTools.fs
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