* Add an exact method to compute an ellipse from three points and two tangents.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 3cfdc89..096fd94 100644 (file)
@@ -31,7 +31,6 @@ let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
     mat.CopyTo(img)
     saveImg img filepath
 
-
 type Histogram = { data: int[]; total: int; sum: int; min: float32; max: float32 }
 
 let histogramImg (img: Image<Gray, float32>) (nbSamples: int) : Histogram =
@@ -148,7 +147,6 @@ let otsu (hist: Histogram) : float32 * float32 * float32 =
 
     toFloat level, toFloat mean1, toFloat mean2
 
-
 let suppressMConnections (img: Matrix<byte>) =
     let w = img.Width
     let h = img.Height
@@ -163,7 +161,6 @@ let suppressMConnections (img: Matrix<byte>) =
             then
                 img.[i, j] <- 0uy
 
-
 let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32> * Image<Gray, float32> =
     let w = img.Width
     let h = img.Height
@@ -198,9 +195,7 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
         let sensibilityHigh = 0.1f
         let sensibilityLow = 0.0f
         use magnitudesByte = magnitudes.Convert<byte>()
-        let threshold = float32 <| CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
         let threshold, _, _ = otsu (histogramMat magnitudes 300)
-
         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
 
     // Non-maximum suppression.
@@ -285,12 +280,10 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
 
     edges, xGradient, yGradient
 
-
 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)
 
-
 type Points = HashSet<Point>
 
 let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
@@ -365,15 +358,12 @@ let findExtremum (img: Image<Gray, 'TDepth>) (extremumType: ExtremumType) : IEnu
 
     result.Select(fun l -> Points(l))
 
-
 let findMaxima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
     findExtremum img ExtremumType.Maxima
 
-
 let findMinima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
     findExtremum img ExtremumType.Minima
 
-
 type PriorityQueue () =
     let size = 256
     let q: Points[] = Array.init size (fun i -> Points())
@@ -472,7 +462,6 @@ type PriorityQueue () =
         highest <- -1
         lowest <- size
 
-
 type private AreaState =
     | Removed = 1
     | Unprocessed = 2
@@ -556,9 +545,9 @@ let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperatio
                         nextElements.Add(p) |> ignore
 
                 else
-                    let m' = pixels.[p.Y, p.X]
-                    if m' <> null
-                    then
+                    match pixels.[p.Y, p.X] with
+                    | null -> ()
+                    | m' ->
                         if m'.Elements.Count + m.Elements.Count <= area
                         then
                             m'.State <- AreaState.Removed
@@ -614,7 +603,6 @@ let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperatio
             | _ -> ()
     ()
 
-
 let areaOpen (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Opening
 
@@ -627,8 +615,7 @@ type Island (cmp: IComparer<float32>) =
     member val Level = 0.f with get, set
     member val Surface = 0 with get, set
 
-
-let private areaOperationF (img: Image<Gray, float32>) (area: int) (op: AreaOperation) =
+let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) =
     let w = img.Width
     let h = img.Height
     let earth = img.Data
@@ -655,82 +642,96 @@ let private areaOperationF (img: Image<Gray, float32>) (area: int) (op: AreaOper
                 let ni = i + p.Y
                 let nj = j + p.X
                 let neighbor = Point(nj, ni)
-                if ni >= 0 && ni < h && nj >= 0 && nj < w && ownership.[ni, nj] = null && not (shorePoints.Contains(neighbor))
+                if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor))
                 then
                     shorePoints.Add(neighbor) |> ignore
                     island.Shore.Add earth.[ni, nj, 0] neighbor
 
-    for island in islands do
-        let mutable stop = island.Shore.IsEmpty
-
-        // 'true' if 'p' is owned or adjacent to 'island'.
-        let ownedOrAdjacent (p: Point) : bool =
-            ownership.[p.Y, p.X] = island ||
-            (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
-            (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
-            (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
-            (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
-
-        while not stop && island.Surface < area do
-            let level, next = island.Shore.Max
-            let other = ownership.[next.Y, next.X]
-            if other = island // During merging, some points on the shore may be owned by the island itself -> ignored.
-            then
-                island.Shore.RemoveNext ()
-            else
-                if other <> null
-                then // We touching another island.
-                    if island.Surface + other.Surface >= area
-                    then
-                        stop <- true
-                    else // We can merge 'other' into 'surface'.
-                        island.Surface <- island.Surface + other.Surface
-                        island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level
-                        for l, p in other.Shore do
-                            let mutable currentY = p.Y + 1
-                            while currentY < h && ownership.[currentY, p.X] = other do
-                                ownership.[currentY, p.X] <- island
-                                currentY <- currentY + 1
-                            island.Shore.Add l p
-                        other.Shore.Clear()
-
-                elif comparer.Compare(level, island.Level) > 0
+    for area, obj in areas do
+        for island in islands do
+            let mutable stop = island.Shore.IsEmpty
+
+            // 'true' if 'p' is owned or adjacent to 'island'.
+            let inline ownedOrAdjacent (p: Point) : bool =
+                ownership.[p.Y, p.X] = island ||
+                (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
+                (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
+                (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
+                (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
+
+            while not stop && island.Surface < area do
+                let level, next = island.Shore.Max
+                let other = ownership.[next.Y, next.X]
+                if other = island // During merging, some points on the shore may be owned by the island itself -> ignored.
                 then
-                    stop <- true
-                else
                     island.Shore.RemoveNext ()
-                    for i, j in se do
-                        let ni = i + next.Y
-                        let nj = j + next.X
-                        if ni < 0 || ni >= h || nj < 0 || nj >= w
+                else
+                    if not <| Object.ReferenceEquals(other, null)
+                    then // We touching another island.
+                        if island.Surface + other.Surface >= area
                         then
-                            island.Surface <- Int32.MaxValue
                             stop <- true
-                        else
-                            let neighbor = Point(nj, ni)
-                            if not <| ownedOrAdjacent neighbor
-                            then
-                                island.Shore.Add earth.[ni, nj, 0] neighbor
-                    if not stop
+                        else // We can merge 'other' into 'surface'.
+                            island.Surface <- island.Surface + other.Surface
+                            island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level
+                            for l, p in other.Shore do
+                                let mutable currentY = p.Y + 1
+                                while currentY < h && ownership.[currentY, p.X] = other do
+                                    ownership.[currentY, p.X] <- island
+                                    currentY <- currentY + 1
+                                island.Shore.Add l p
+                            other.Shore.Clear()
+
+                    elif comparer.Compare(level, island.Level) > 0
                     then
-                        ownership.[next.Y, next.X] <- island
-                        island.Level <- level
-                        island.Surface <- island.Surface + 1
-
-    for i in 0 .. h - 1 do
-        for j in 0 .. w - 1 do
-            let island = ownership.[i, j]
-            if island <> null
-            then
-                earth.[i, j, 0] <- island.Level
+                        stop <- true
+                    else
+                        island.Shore.RemoveNext ()
+                        for i, j in se do
+                            let ni = i + next.Y
+                            let nj = j + next.X
+                            if ni < 0 || ni >= h || nj < 0 || nj >= w
+                            then
+                                island.Surface <- Int32.MaxValue
+                                stop <- true
+                            else
+                                let neighbor = Point(nj, ni)
+                                if not <| ownedOrAdjacent neighbor
+                                then
+                                    island.Shore.Add earth.[ni, nj, 0] neighbor
+                        if not stop
+                        then
+                            ownership.[next.Y, next.X] <- island
+                            island.Level <- level
+                            island.Surface <- island.Surface + 1
+
+        let mutable diff = 0.f
+
+        for i in 0 .. h - 1 do
+            for j in 0 .. w - 1 do
+                match ownership.[i, j] with
+                | null -> ()
+                | island ->
+                    let l = island.Level
+                    diff <- diff + l - earth.[i, j, 0]
+                    earth.[i, j, 0] <- l
+
+        match f with
+        | Some f' -> f' obj diff
+        | _ -> ()
     ()
 
-
 let areaOpenF (img: Image<Gray, float32>) (area: int) =
-    areaOperationF img area AreaOperation.Opening
+    areaOperationF img [ area, () ] None AreaOperation.Opening
 
 let areaCloseF (img: Image<Gray, float32>) (area: int) =
-    areaOperationF img area AreaOperation.Closing
+    areaOperationF img [ area, () ] None AreaOperation.Closing
+
+let areaOpenFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
+    areaOperationF img areas (Some f) AreaOperation.Opening
+
+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) =
@@ -787,7 +788,6 @@ let areaOpen2 (img: Image<Gray, byte>) (area: int) =
                             for p in pointsChecked do
                                 imgData.[p.Y, p.X, 0] <- maxNeighborValue
 
-
 // Zhang and Suen algorithm.
 // Modify 'mat' in place.
 let thin (mat: Matrix<byte>) =
@@ -838,7 +838,6 @@ 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.
 let removeArea (mat: Matrix<byte>) (areaSize: int) =
@@ -907,17 +906,13 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
 
     List<Point>(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);
 
-
 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) =
-
     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)
@@ -938,11 +933,9 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo
             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) =
     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) =
     if drawCellContent