// We reverse the list to get the lower score ellipses first.
let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev
- // 2) Remove ellipses with a high standard deviation (high contrast).
+
+ // 2) Remove ellipses touching the edges.
+ for e in ellipses do
+ if e.isOutside w_f h_f then e.Removed <- true
+
+ // 3) Remove ellipses with a high standard deviation (high contrast).
// CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed.
yield float img.Data.[y, x, 0] })
for e in ellipses do
- let minX, minY, maxX, maxY = ellipseWindow e
-
- let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {
- for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
- for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
- if e.Contains (float x) (float y)
- then
- yield float img.Data.[y, x, 0] })
+ if not e.Removed
+ then
+ let shrinkedE = e.Scale 0.9
+ let minX, minY, maxX, maxY = ellipseWindow shrinkedE
- if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
- e.Removed <- true
+ let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {
+ for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
+ for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
+ if shrinkedE.Contains (float x) (float y)
+ then
+ yield float img.Data.[y, x, 0] })
- // 3) Remove ellipses touching the edges.
- for e in ellipses do
- if e.isOutside w_f h_f then e.Removed <- true
+ if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
+ e.Removed <- true
// 4) Remove ellipses with little area.
let minArea = config.RBCMinArea
p.[k] <- p.[k] / p.[0]
p.[0] <- 1.0
let e = 0.25 * p.[1]
- let b = ref (2.0 * e)
- let c = ref (!b ** 2.0)
- let mutable d = 0.75 * !c
- b := p.[3] + !b *(!c - p.[2])
+ let mutable b = 2.0 * e
+ let mutable c = b ** 2.0
+ let mutable d = 0.75 * c
+ b <- p.[3] + b *(c - p.[2])
let mutable a = p.[2] - d
- c := p.[4] + e * (e * a - p.[3])
+ c <- p.[4] + e * (e * a - p.[3])
a <- a - d
- let quadExecuted = ref false
+ let mutable quadExecuted = false
let quad () =
- if not !quadExecuted
+ if not quadExecuted
then
- p.[2] <- !c / !b
+ p.[2] <- c / b
quadroots p r
for k in 1..2 do
for j in 1..2 do
r.[j, k+2] <- r.[j, k]
p.[1] <- -p.[1]
- p.[2] <- !b
+ p.[2] <- b
quadroots p r
for k in 1..4 do
r.[1,k] <- r.[1,k] - e
- quadExecuted := true
+ quadExecuted <- true
p.[1] <- 0.5 * a
- p.[2] <- (p.[1] * p.[1] - !c) * 0.25
- p.[3] <- !b * !b / -64.0
+ p.[2] <- (p.[1] * p.[1] - c) * 0.25
+ p.[3] <- b * b / -64.0
if p.[3] < 0.0
then
cubicroots p r
then
d <- r.[1, k] * 4.0
a <- a + d
- if a >= 0.0 && !b >= 0.0
+ if a >= 0.0 && b >= 0.0
then
p.[1] <- sqrt d
- elif a <= 0.0 && !b <= 0.0
+ elif a <= 0.0 && b <= 0.0
then
p.[1] <- sqrt d
else
p.[1] <- -(sqrt d)
- b := 0.5 * (a + !b / p.[1])
+ b <- 0.5 * (a + b / p.[1])
quad ()
k <- 4
k <- k + 1
- if not !quadExecuted && p.[2] < 0.0
+ if not quadExecuted && p.[2] < 0.0
then
- b := sqrt !c
- d <- !b + !b - a
+ b <- sqrt c
+ d <- b + b - a
p.[1] <- 0.0
if d > 0.0
then
p.[1] <- sqrt d
- elif not !quadExecuted
+ elif not quadExecuted
then
if p.[1] > 0.0
then
- b := (sqrt p.[2]) * 2.0 + p.[1]
+ b <- (sqrt p.[2]) * 2.0 + p.[1]
else
- b := -(sqrt p.[2]) * 2.0 + p.[1]
+ b <- -(sqrt p.[2]) * 2.0 + p.[1]
- if !b <> 0.0
+ if b <> 0.0
then
p.[1] <- 0.0
else
for k in 1..4 do
r.[1, k] <- -e
r.[2, k] <- 0.0
- quadExecuted := true
+ quadExecuted <- true
quad ()
printf "nroots = %d\n" nroots
#endif
- let ychk = [|
- for i in 1 .. nroots do
- if abs r.[2, i] < EPS
- then
- yield r.[1, i] * b1
+ let ychk = Array.init nroots (fun _ -> Double.MaxValue)
+ let mutable nychk = 0
+ for i in 1 .. nroots do
+ if abs r.[2, i] < EPS
+ then
+ ychk.[nychk] <- r.[1, i] * b1
+ nychk <- nychk + 1
#if DEBUG_LOG
- printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+ printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
#endif
- |]
Array.sortInPlace ychk
#if DEBUG_LOG
printf "\t j=%d, ychk=%f\n" j ychk.[j]
#endif
- let nychk = Array.length ychk
let mutable nintpts = 0
let xint = Array.zeroCreate 4
let inline private left (i: int) : int = 2 * (i + 1) - 1
let inline private right (i: int) : int = 2 * (i + 1)
-// TODO: measure performance with a struct.
-type private Node<'k, 'v> = { mutable key: 'k; value : 'v } with
+[<Struct>]
+type private Node<'k, 'v> =
+ val mutable key: 'k
+ val value: 'v
+ new (k, v) = { key = k; value = v }
override this.ToString () = sprintf "%A -> %A" this.key this.value
type Heap<'k, 'v> (kComparer : IComparer<'k>) =
a.[max] <- tmp
heapUp max
+ // Check the integrity of the heap, use for debug only.
let rec checkIntegrity (i: int) : bool =
let l, r = left i, right i
-
let leftIntegrity =
if l < a.Count
then
(this :> IEnumerable<'k * 'v>).GetEnumerator() :> System.Collections.IEnumerator
member this.Next () : 'k * 'v =
- let { key = key; value = value } = a.[0]
+ let node = a.[0]
a.[0] <- a.[a.Count - 1]
a.RemoveAt(a.Count - 1)
heapUp 0
- key, value
+ node.key, node.value
member this.RemoveNext () =
a.[0] <- a.[a.Count - 1]
heapUp 0
member this.Add (key: 'k) (value: 'v) =
- a.Add({ key = key; value = value })
+ a.Add(Node(key, value))
let mutable i = a.Count - 1
while i > 0 && kComparer.Compare(a.[parent i].key, a.[i].key) < 0 do
then
img.[i, j] <- 0uy
+
let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
let w = img.Width
let h = img.Height
CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles).
let thresholdHigh, thresholdLow =
- let sensibility = 0.1
+ let sensibilityHigh = 0.1
+ let sensibilityLow = 0.1
use magnitudesByte = magnitudes.Convert<byte>()
let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
- threshold + (sensibility * threshold), threshold - (sensibility * threshold)
+ threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
// Non-maximum suppression.
use nms = new Matrix<byte>(xGradient.Size)
let mNeigbors (sign: int) : float =
if angle < Math.PI / 4.
- then
- ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign]
+ then ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign]
elif angle < Math.PI / 2.
- then
- ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j]
+ then ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j]
elif angle < 3.0 * Math.PI / 4.
- then
- ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign]
+ then ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign]
elif angle < Math.PI
- then
- ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign]
+ then ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign]
elif angle < 5. * Math.PI / 4.
- then
- ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign]
+ then ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign]
elif angle < 3. * Math.PI / 2.
- then
- ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j]
+ then ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j]
elif angle < 7. * Math.PI / 4.
- then
- ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign]
- else
- ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign]
+ then ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign]
+ else ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign]
let m = magnitudes.Data.[i, j]
if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
logTime "areaOpen 1" (fun () -> ImgTools.areaOpenF filteredGreen config.Parameters.initialAreaOpen)
- config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.5 |> float)
+ config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.3 |> float)
let secondAreaOpen = int <| config.RBCArea / 3.
if secondAreaOpen > config.Parameters.initialAreaOpen
let ellipses = List<EllipseScoreFlaggedKd>()
// All ellipses with a score below this are removed.
- let matchingScoreThreshold2 = 25. * radiusMin // 600.
+ let matchingScoreThreshold2 = 20. * radiusMin
member this.Add (e: Ellipse) =
ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
factorWindowSize = 2.0
darkStainLevel = 0.25 // Lower -> more sensitive. 0.3
- maxDarkStainRatio = 0.1
+ maxDarkStainRatio = 0.1 // 10 %
infectionArea = 0.012 // 1.2 %
infectionLevel = 1.12 // Lower -> more sensitive.
stainLevel = 1.1 // Lower -> more sensitive.
maxStainRatio = 0.12 // 12 %
- standardDeviationMaxRatio = 0.58 // 0.55
+ standardDeviationMaxRatio = 0.5 // 0.55
minimumCellArea = 0.5 })
match mode with