* Add area granulometry (not used for the moment)
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 3cfdc89..23c9aee 100644 (file)
@@ -198,9 +198,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.
@@ -556,9 +554,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
@@ -628,7 +626,7 @@ type Island (cmp: IComparer<float32>) =
     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 +653,97 @@ 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) =