// The minimum window to contain a given ellipse.
let ellipseWindow (e: Ellipse) =
let cx, cy = roundInt e.Cx, roundInt e.Cy
- let a = int (e.A + 0.5)
+ let a = int (e.A + 0.5f)
cx - a, cy - a, cx + a, cy + a
let w = img.Width
- let w_f = float w
+ let w_f = float32 w
let h = img.Height
- let h_f = float h
+ let h_f = float32 h
// Return 'true' if the point 'p' is owned by e.
// The lines represents all intersections with other ellipses.
tree.Search (searchRegion e)
// We only keep the ellipses touching 'e'.
|> List.choose (fun otherE ->
- match EEOver.EEOverlapArea e otherE with
- | Some (_, px, _) when px.Length > 2 ->
- otherE.Removed <- true
- None
- | Some (area, px, py) when area > 0.0 && px.Length = 2 ->
- Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1]))
- | _ ->
+ if e <> otherE
+ then
+ match EEOver.EEOverlapArea e otherE with
+ | Some (_, px, _) when px.Length > 2 ->
+ otherE.Removed <- true
+ None
+ | Some (area, px, py) when area > 0.f && px.Length = 2 ->
+ Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1]))
+ | _ ->
+ None
+ else
None )
else
[]
// 3) Remove ellipses with a high standard deviation (high contrast).
// CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed.
-
+ let imgData = img.Data
let globalStdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation(seq {
for y in 0 .. h - 1 do
for x in 0 .. w - 1 do
- yield float img.Data.[y, x, 0] })
+ yield float imgData.[y, x, 0] })
for e in ellipses do
if not e.Removed
then
- let shrinkedE = e.Scale 0.9
+ let shrinkedE = e.Scale 0.9f
let minX, minY, maxX, maxY = ellipseWindow shrinkedE
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)
+ if shrinkedE.Contains (float32 x) (float32 y)
then
- yield float img.Data.[y, x, 0] })
+ yield float imgData.[y, x, 0] })
if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
e.Removed <- true
let mutable area = 0
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
- let p = PointD(float x, float y)
+ let p = PointD(float32 x, float32 y)
if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
then
area <- area + 1
let elements = new Matrix<byte>(maxY - minY + 1, maxX - minX + 1)
for y in minY .. maxY do
for x in minX .. maxX do
- let p = PointD(float x, float y)
+ let p = PointD(float32 x, float32 y)
if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
then
elements.[y-minY, x-minX] <- 1uy
let cellClass =
if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement) ||
- float stainPixels > config.Parameters.maxStainRatio * (float nbElement) (* ||
- sqrt (((float sumCoords_x) / (float nbElement) - e.Cx) ** 2.0 + ((float sumCoords_y) / (float nbElement) - e.Cy) ** 2.0) > e.A * config.maxOffcenter *)
+ float stainPixels > config.Parameters.maxStainRatio * (float nbElement)
then
Peculiar
elif infectedPixels.Count >= 1
open System
+open Const
+
type Debug =
| DebugOff
| DebugOn of string // Output directory.
initialAreaOpen: int
- minRbcRadius: float
- maxRbcRadius: float
+ minRbcRadius: float32
+ maxRbcRadius: float32
preFilterSigma: float
// Ellipse.
factorNbPick: float
- factorWindowSize: float // factor of 'maxRBCSize'.
// Parasites detection.
darkStainLevel: float
maxDarkStainRatio: float
- stainArea: float // Factor of a RBC area. 0.5 means the half of RBC area.
+ stainArea: float32 // Factor of a RBC area. 0.5 means the half of RBC area.
stainLevel: float // > 1
maxStainRatio: float // [0, 1]
- infectionArea: float // Factor of a RBC area. 0.5 means the half of RBC area.
+ infectionArea: float32 // Factor of a RBC area. 0.5 means the half of RBC area.
infectionLevel: float // > 1
standardDeviationMaxRatio: float // The standard deviation of the pixel values of a cell can't be greater than standardDeviationMaxRatio * global standard deviation
- minimumCellArea: float // Factor of the nominal RBC area.
+ minimumCellArea: float32 // Factor of the nominal RBC area.
}
type Config (param: Parameters) =
member this.Parameters = param
member val Debug = DebugOff with get, set
- member val RBCRadius = 30. with get, set
+ member val RBCRadius = 30.f with get, set
member this.RBCMinRadius = this.RBCRadius + param.minRbcRadius * this.RBCRadius
member this.RBCMaxRadius = this.RBCRadius + param.maxRbcRadius * this.RBCRadius
- member this.RBCArea = Math.PI * this.RBCRadius ** 2.0
+ member this.RBCArea = PI * this.RBCRadius ** 2.f
member this.RBCMinArea = param.minimumCellArea * this.RBCArea
member this.InfectionArea = param.infectionArea * this.RBCArea
member this.StainArea = param.stainArea * this.RBCArea
+
+ member this.Copy () =
+ this.MemberwiseClone() :?> Config
quad ()
// Return a tuple (area, x intersections, y intersections)
-let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * float[]) option =
- let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
- let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
+let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
+ let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
+ let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
then
| 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| _ -> -1.0
if nintpts = 0
- then Some (area, [||], [||])
+ then Some (float32 area, [||], [||])
else
- let xTransform = Array.zeroCreate nintpts
- let yTransform = Array.zeroCreate nintpts
+ let xTransform : float32[] = Array.zeroCreate nintpts
+ let yTransform : float32[] = Array.zeroCreate nintpts
for i in 0 .. (nintpts - 1) do
- xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
- yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
- Some (area, xTransform, yTransform)
\ No newline at end of file
+ xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
+ yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
+ Some (float32 area, xTransform, yTransform)
\ No newline at end of file
open System
open System.Collections.Generic
+open System.Drawing
open Emgu.CV
open Emgu.CV.Structure
open Utils
open Config
open MatchingEllipses
-
+open Const
type private SearchExtremum = Minimum | Maximum
-let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) (xmax: float) (searchExtremum: SearchExtremum) : (float * float) =
- let gr = 1.0 / 1.6180339887498948482
+let private goldenSectionSearch (f: float32 -> float32) (nbIter: int) (xmin: float32) (xmax: float32) (searchExtremum: SearchExtremum) : (float32 * float32) =
+ let gr = 1.f / 1.6180339887498948482f
let mutable a = xmin
let mutable b = xmax
let mutable c = b - gr * (b - a)
c <- d
d <- a + gr * (b - a)
- let x = (b + a) / 2.0
+ let x = (b + a) / 2.f
x, f x
// Ellipse.A is always equal or greater than Ellipse.B.
// Ellipse.Alpha is between 0 and Pi.
-let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option =
+let ellipse (p1x: float32) (p1y: float32) (m1: float32) (p2x: float32) (p2y: float32) (m2: float32) (p3x: float32) (p3y: float32) : Types.Ellipse option =
let accuracy_extremum_search_1 = 8 // 3
let accuracy_extremum_search_2 = 8 // 4
let alpha1 = atan m1
let alpha2 = atan m2
- let r1 = sqrt (p1x ** 2.0 + p1y ** 2.0)
+ let r1 = sqrt (p1x ** 2.f + p1y ** 2.f)
let theta1 = atan2 p1y p1x
- let r2 = sqrt (p2x**2.0 + p2y**2.0)
+ let r2 = sqrt (p2x ** 2.f + p2y ** 2.f)
let theta2 = atan2 p2y p2x
let valid =
- 4.0 * sin (alpha1 - theta1) * (-r1 * sin (alpha1 - theta1) + r2 * sin (alpha1 - theta2)) *
+ 4.f * sin (alpha1 - theta1) * (-r1 * sin (alpha1 - theta1) + r2 * sin (alpha1 - theta2)) *
sin (alpha2 - theta2) * (-r1 * sin (alpha2 - theta1) + r2 * sin (alpha2 - theta2)) +
- r1 * r2 * sin (alpha1 - alpha2) ** 2.0 * sin (theta1 - theta2) ** 2.0 < 0.0
+ r1 * r2 * sin (alpha1 - alpha2) ** 2.f * sin (theta1 - theta2) ** 2.f < 0.f
if valid
then
let r theta =
(r1 * r2 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0)
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.f - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.f)
let rabs = r >> abs
// We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval.
- let thetaTan, _ = goldenSectionSearch rabs accuracy_extremum_search_1 0.0 Math.PI Maximum
+ let thetaTan, _ = goldenSectionSearch rabs accuracy_extremum_search_1 0.f PI Maximum
let rTan = r thetaTan
let PTanx = rTan * cos thetaTan
let d2a = tan alpha2
let d2b = -d2a * p2x + p2y
- let d3a = -1.0 / tan thetaTan
+ let d3a = -1.f / tan thetaTan
let d3b = -d3a * PTanx + PTany
let Ux = -(d1b - d2b) / (d1a - d2a)
let Vx = -(d1b - d3b) / (d1a - d3a)
let Vy = -(d3a * d1b - d1a * d3b) / (d1a - d3a)
- let Wx = p1x + (p2x - p1x) / 2.0
- let Wy = p1y + (p2y - p1y) / 2.0
+ let Wx = p1x + (p2x - p1x) / 2.f
+ let Wy = p1y + (p2y - p1y) / 2.f
- let Zx = p1x + (PTanx - p1x) / 2.0
- let Zy = p1y + (PTany - p1y) / 2.0
+ let Zx = p1x + (PTanx - p1x) / 2.f
+ let Zy = p1y + (PTany - p1y) / 2.f
let va = -(-Vy + Zy) / (Vx - Zx)
let vb = -(Zx * Vy - Vx * Zy) / (Vx - Zx)
let cx = -(vb - ub) / (va - ua)
let cy = -(ua * vb - va * ub) / (va - ua)
- let rc = sqrt (cx**2.0 + cy**2.0)
+ let rc = sqrt (cx ** 2.f + cy ** 2.f)
let psi = atan2 cy cx
let rellipse theta =
sqrt (
- rc ** 2.0 + (r1 ** 2.0 * r2 ** 2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) ** 2.0 * sin (theta1 - theta2) ** 2.0) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0) ** 2.0 -
- (2.0 * r1 * r2 * rc * cos (theta - psi) * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0))
+ rc ** 2.f + (r1 ** 2.f * r2 ** 2.f * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) ** 2.f * sin (theta1 - theta2) ** 2.f) /
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.f - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.f) ** 2.f -
+ (2.f * r1 * r2 * rc * cos (theta - psi) * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) /
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.f - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.f))
// We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval.
- let r1eTheta, r1e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Maximum // Pi/2 and not pi because the period is Pi.
- let r2eTheta, r2e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Minimum
+ let r1eTheta, r1e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.f (PI / 2.f) Maximum // Pi/2 and not pi because the period is Pi.
+ let r2eTheta, r2e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.f (PI / 2.f) Minimum
let rr1e = r r1eTheta
let r1ex = rr1e * cos r1eTheta
let r1ey = rr1e * sin r1eTheta
let mutable alpha = atan ((r1ey - cy) / (r1ex - cx))
- if alpha < 0.0
+ if alpha < 0.f
then
- alpha <- alpha + Math.PI
+ alpha <- alpha + PI
// Ride off the p3 referential.
let cx = cx + p3x
None
-let private vectorRotation (p1x: float) (p1y: float) (v1x: float) (v1y: float) (px: float) (py: float) : float =
- let mutable rotation = 1.0
+let private vectorRotation (p1x: float32) (p1y: float32) (v1x: float32) (v1y: float32) (px: float32) (py: float32) : float32 =
+ let mutable rotation = 1.f
if p1y > py
then
- if v1x > 0.0
+ if v1x > 0.f
then
- rotation <- -1.0
+ rotation <- -1.f
elif p1y < py
then
- if v1x < 0.0
+ if v1x < 0.f
then
- rotation <- -1.0
+ rotation <- -1.f
elif p1x > px
then
- if v1y < 0.0
+ if v1y < 0.f
then
- rotation <- -1.0
+ rotation <- -1.f
elif p1x < px
then
- if v1y > 0.0
+ if v1y > 0.f
then
- rotation <- -1.0
+ rotation <- -1.f
rotation
-let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) (v1x: float) (v1y: float) (v2x: float) (v2y: float) : (float * float) option =
+let private areVectorsValid (p1x: float32) (p1y: float32) (p2x: float32) (p2y: float32) (v1x: float32) (v1y: float32) (v2x: float32) (v2y: float32) : (float32 * float32) option =
let m1 = -v1x / v1y
let m2 = -v2x / v2y
let alpha1 = atan2 (p1y - py) (p1x - px)
let alpha2 = atan2 (p2y - py) (p2x - px)
- let alpha1' = if alpha1 < 0.0 then 2.0 * Math.PI + alpha1 else alpha1
- let alpha2' = if alpha2 < 0.0 then 2.0 * Math.PI + alpha2 else alpha2
+ let alpha1' = if alpha1 < 0.f then 2.f * PI + alpha1 else alpha1
+ let alpha2' = if alpha2 < 0.f then 2.f * PI + alpha2 else alpha2
let diff = rot1 * alpha1' + rot2 * alpha2'
- if diff > Math.PI || (diff < 0.0 && diff > -Math.PI)
+ if diff > PI || (diff < 0.f && diff > -PI)
then
None
else
let find (edges: Matrix<byte>)
- (xGradient: Image<Gray, float>)
- (yGradient: Image<Gray, float>)
+ (xGradient: Image<Gray, float32>)
+ (yGradient: Image<Gray, float32>)
(config: Config) : MatchingEllipses =
let r1, r2 = config.RBCMinRadius, config.RBCMaxRadius
- let windowSize = roundInt (config.Parameters.factorWindowSize * r2)
+ let incrementWindowDivisor = 4.f
+
+ // We choose a window size for which the biggest ellipse can always be fitted in.
+ let windowSize = roundInt (2.f * r2 / (incrementWindowDivisor - 1.f) * incrementWindowDivisor)
let factorNbPick = config.Parameters.factorNbPick
- let increment = windowSize / 4
+ let increment = windowSize / (int incrementWindowDivisor)
- let radiusTolerance = (r2 - r1) * 0.2
+ let radiusTolerance = (r2 - r1) * 0.2f
- let minimumDistance = (r2 / 1.5) ** 2.0
- let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0
+ let squaredMinimumDistance = (r2 / 1.5f) ** 2.f
+ let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.f + (y1 - y2) ** 2.f
let h = edges.Height
let w = edges.Width
+ let h_f = float32 h
+ let w_f = float32 w
let mutable last_i, last_j = Int32.MaxValue, Int32.MaxValue
- let currentElements = List<(int * int)>()
+ let currentElements = List<Point>()
let edgesData = edges.Data
let xDirData = xGradient.Data
let window_j_end = if window_j + windowSize - 1 >= w then w - 1 else window_j + windowSize - 1
// Remove old elements.
- let indexFirstElement = currentElements.FindIndex(fun (_, pj) -> pj >= window_j)
+ let indexFirstElement = currentElements.FindIndex(fun p -> p.X >= window_j_begin)
if indexFirstElement > 0
then currentElements.RemoveRange(0, indexFirstElement)
// Add the new elements.
- for j in window_j + windowSize - increment .. window_j + windowSize - 1 do
+ let newElemsBegin_j = window_j + windowSize - increment
+ let newElemsEnd_j = window_j + windowSize - 1
+ for j in (if newElemsBegin_j < 0 then 0 else newElemsBegin_j) .. (if newElemsEnd_j >= w then w - 1 else newElemsEnd_j) do
for i in window_i_begin .. window_i_end do
- if j >= 0 && j < w && edgesData.[i, j] = 1uy
- then currentElements.Add((i, j))
+ if edgesData.[i, j] = 1uy
+ then currentElements.Add(Point(j, i))
if currentElements.Count >= 10
then
let mutable nbOfPicks = (float currentElements.Count) * factorNbPick |> int
while nbOfPicks > 0 do
- let (p1y, p1x) as p1 = currentElements.[rng.Next(currentElements.Count)]
- let (p2y, p2x) as p2 = currentElements.[rng.Next(currentElements.Count)]
- let (p3y, p3x) as p3 = currentElements.[rng.Next(currentElements.Count)]
+ let p1 = currentElements.[rng.Next(currentElements.Count)]
+ let p2 = currentElements.[rng.Next(currentElements.Count)]
+ let p3 = currentElements.[rng.Next(currentElements.Count)]
if p1 <> p2 && p1 <> p3 && p2 <> p3
then
nbOfPicks <- nbOfPicks - 1
- let p1yf, p1xf = float p1y, float p1x
- let p2yf, p2xf = float p2y, float p2x
- let p3yf, p3xf = float p3y, float p3x
- if squaredDistance p1xf p1yf p2xf p2yf >= minimumDistance &&
- squaredDistance p1xf p1yf p3xf p3yf >= minimumDistance &&
- squaredDistance p2xf p2yf p3xf p3yf >= minimumDistance
+ let p1yf, p1xf = float32 p1.Y, float32 p1.X
+ let p2yf, p2xf = float32 p2.Y, float32 p2.X
+ let p3yf, p3xf = float32 p3.Y, float32 p3.X
+ if squaredDistance p1xf p1yf p2xf p2yf >= squaredMinimumDistance &&
+ squaredDistance p1xf p1yf p3xf p3yf >= squaredMinimumDistance &&
+ squaredDistance p2xf p2yf p3xf p3yf >= squaredMinimumDistance
then
- match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0] with
+ match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1.Y, p1.X, 0] -yDirData.[p1.Y, p1.X, 0] -xDirData.[p2.Y, p2.X, 0] -yDirData.[p2.Y, p2.X, 0] with
| Some (m1, m2) ->
match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with
- | Some e when e.Cx > 0.0 && e.Cx < (float w) - 1.0 && e.Cy > 0.0 && e.Cy < (float h) - 1.0 &&
+ | Some e when e.Cx > 0.f && e.Cx < w_f - 1.f && e.Cy > 0.f && e.Cy < h_f - 1.f &&
e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
-
- let prout = areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0]
ellipses.Add e
| _ -> ()
| _ -> ()
// 'range': a minimum and maximum radius.
// 'scale': <= 1.0, to speed up the process.
let findRadius (img: Image<Gray, 'TDepth>) (range: int * int) (scale: float) : int =
- use scaledImg = if scale = 1.0 then img else img.Resize(scale, CvEnum.Inter.Area)
+ use scaledImg = if scale = 1. then img else img.Resize(scale, CvEnum.Inter.Area)
let r1, r2 = range
let r1', r2' = roundInt (float r1 * scale), roundInt (float r2 * scale)
open Emgu.CV
open Emgu.CV.Structure
-open Utils
open Heap
+open Const
+open Utils
// Normalize image values between 0uy and 255uy.
-let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
+let normalizeAndConvert (img: Image<Gray, 'TDepth>) : Image<Gray, byte> =
let min = ref [| 0.0 |]
let minLocation = ref <| [| Point() |]
let max = ref [| 0.0 |]
let maxLocation = ref <| [| Point() |]
img.MinMax(min, max, minLocation, maxLocation)
- ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
+ ((img.Convert<Gray, float32>() - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
img.[i, j] <- 0uy
-let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
+let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32> * Image<Gray, float32> =
let w = img.Width
let h = img.Height
[ 2.0f; 0.0f; -2.0f ]
[ 1.0f; 0.0f; -1.0f ]], Point(1, 1))
- let xGradient = img.Convolution(sobelKernel).Convert<Gray, float>()
- let yGradient = img.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+ let xGradient = img.Convolution(sobelKernel)
+ let yGradient = img.Convolution(sobelKernel.Transpose())
let xGradientData = xGradient.Data
let yGradientData = yGradient.Data
for r in 0 .. h - 1 do
- xGradientData.[r, 0, 0] <- 0.0
- xGradientData.[r, w - 1, 0] <- 0.0
- yGradientData.[r, 0, 0] <- 0.0
- yGradientData.[r, w - 1, 0] <- 0.0
+ xGradientData.[r, 0, 0] <- 0.f
+ xGradientData.[r, w - 1, 0] <- 0.f
+ yGradientData.[r, 0, 0] <- 0.f
+ yGradientData.[r, w - 1, 0] <- 0.f
for c in 0 .. w - 1 do
- xGradientData.[0, c, 0] <- 0.0
- xGradientData.[h - 1, c, 0] <- 0.0
- yGradientData.[0, c, 0] <- 0.0
- yGradientData.[h - 1, c, 0] <- 0.0
+ xGradientData.[0, c, 0] <- 0.f
+ xGradientData.[h - 1, c, 0] <- 0.f
+ yGradientData.[0, c, 0] <- 0.f
+ yGradientData.[h - 1, c, 0] <- 0.f
- use magnitudes = new Matrix<float>(xGradient.Size)
- use angles = new Matrix<float>(xGradient.Size)
+ 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).
let thresholdHigh, thresholdLow =
- let sensibilityHigh = 0.1
- let sensibilityLow = 0.1
+ let sensibilityHigh = 0.1f
+ let sensibilityLow = 0.1f
use magnitudesByte = magnitudes.Convert<byte>()
- let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+ let threshold = float32 <| CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
// Non-maximum suppression.
use nms = new Matrix<byte>(xGradient.Size)
+ let nmsData = nms.Data
+ let anglesData = angles.Data
+ let magnitudesData = magnitudes.Data
+ let xGradientData = xGradient.Data
+ let yGradientData = yGradient.Data
+
+ let PI = float32 Math.PI
+
for i in 0 .. h - 1 do
- nms.Data.[i, 0] <- 0uy
- nms.Data.[i, w - 1] <- 0uy
+ nmsData.[i, 0] <- 0uy
+ nmsData.[i, w - 1] <- 0uy
for j in 0 .. w - 1 do
- nms.Data.[0, j] <- 0uy
- nms.Data.[h - 1, j] <- 0uy
+ nmsData.[0, j] <- 0uy
+ nmsData.[h - 1, j] <- 0uy
for i in 1 .. h - 2 do
for j in 1 .. w - 2 do
- let vx = xGradient.Data.[i, j, 0]
- let vy = yGradient.Data.[i, j, 0]
- if vx <> 0. || vy <> 0.
+ let vx = xGradientData.[i, j, 0]
+ let vy = yGradientData.[i, j, 0]
+ if vx <> 0.f || vy <> 0.f
then
- let angle = angles.[i, j]
+ let angle = anglesData.[i, j]
let vx', vy' = abs vx, abs vy
let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
- let ratio1 = 1. - ratio2
-
- let mNeigbors (sign: int) : float =
- if angle < Math.PI / 4.
- 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]
- elif angle < 3.0 * Math.PI / 4.
- 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]
- elif angle < 5. * Math.PI / 4.
- 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]
- 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]
-
- let m = magnitudes.Data.[i, j]
+ let ratio1 = 1.f - ratio2
+
+ let mNeigbors (sign: int) : float32 =
+ if angle < PI / 4.f
+ then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
+ elif angle < PI / 2.f
+ then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
+ elif angle < 3.f * PI / 4.f
+ then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
+ elif angle < PI
+ then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
+ elif angle < 5.f * PI / 4.f
+ then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
+ elif angle < 3.f * PI / 2.f
+ then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
+ elif angle < 7.f * PI / 4.f
+ then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
+ else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
+
+ let m = magnitudesData.[i, j]
if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
then
- nms.Data.[i, j] <- 1uy
+ nmsData.[i, j] <- 1uy
// suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
let edges = new Matrix<byte>(xGradient.Size)
+ let edgesData = edges.Data
- // Histeresis thresholding.
+ // Hysteresis thresholding.
let toVisit = Stack<Point>()
for i in 0 .. h - 1 do
for j in 0 .. w - 1 do
- if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh
+ if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh
then
- nms.Data.[i, j] <- 0uy
+ nmsData.[i, j] <- 0uy
toVisit.Push(Point(j, i))
while toVisit.Count > 0 do
let p = toVisit.Pop()
- edges.Data.[p.Y, p.X] <- 1uy
+ edgesData.[p.Y, p.X] <- 1uy
for i' in -1 .. 1 do
for j' in -1 .. 1 do
if i' <> 0 || j' <> 0
then
let ni = p.Y + i'
let nj = p.X + j'
- if ni >= 0 && ni < h && nj >= 0 && nj < w && nms.Data.[ni, nj] = 1uy
+ if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
then
- nms.Data.[ni, nj] <- 0uy
+ nmsData.[ni, nj] <- 0uy
toVisit.Push(Point(nj, ni))
edges, xGradient, yGradient
( 0, -1) // p8
(-1, -1) |] // p9
- let mat' = new Matrix<byte>(mat.Size)
+ use mat' = new Matrix<byte>(mat.Size)
let w = mat'.Width
let h = mat'.Height
mat.CopyTo(mat')
if alpha >= 1.0
then
- img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
+ 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)
else
- let windowPosX = e.Cx - e.A - 5.0
- let gapX = windowPosX - (float (int windowPosX))
+ let windowPosX = e.Cx - e.A - 5.f
+ let gapX = windowPosX - (float32 (int windowPosX))
- let windowPosY = e.Cy - e.A - 5.0
- let gapY = windowPosY - (float (int windowPosY))
+ let windowPosY = e.Cy - e.A - 5.f
+ let gapY = windowPosY - (float32 (int windowPosY))
- let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int)
+ let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int)
img.ROI <- roi
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. + gapX) , float32 <| (e.A + 5. + gapY)), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias)
+ 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)
CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
img.ROI <- Rectangle.Empty
type Result = {
fg: Image<Gray, byte>
- mean_bg: float
- mean_fg: float
+ mean_bg: float32
+ mean_fg: float32
d_fg: Image<Gray, float32> } // Euclidean distances of the foreground to mean_fg.
-let kmeans (img: Image<Gray, float32>) (fgFactor: float) : Result =
+let kmeans (img: Image<Gray, float32>) : Result =
let nbIteration = 3
let w = img.Width
let h = img.Height
let maxLocation = ref <| [| Point() |]
img.MinMax(min, max, minLocation, maxLocation)
- let mutable mean_bg = (!max).[0] - ((!max).[0] - (!min).[0]) / 4.0
- let mutable mean_fg = (!min).[0] + ((!max).[0] - (!min).[0]) / 4.0
- use mutable d_bg = new Image<Gray, float32>(img.Size)
- let mutable d_fg = new Image<Gray, float32>(img.Size)
- let mutable fg = new Image<Gray, byte>(img.Size)
+ let minf = float32 (!min).[0]
+ let maxf = float32 (!max).[0]
+
+ let mutable mean_bg = maxf - (maxf - minf) / 4.f
+ let mutable mean_fg = minf + (maxf - minf) / 4.f
+ use mutable d_bg : Image<Gray, float32> = null
+ let mutable d_fg : Image<Gray, float32> = null
+ let fg = new Image<Gray, byte>(img.Size)
+
+ let imgData = img.Data
+ let fgData = fg.Data
for i in 1 .. nbIteration do
- d_bg <- img.AbsDiff(Gray(mean_bg))
- d_fg <- img.AbsDiff(Gray(mean_fg))
+ if d_bg <> null
+ then
+ d_bg.Dispose()
+ d_fg.Dispose()
+
+ // EmGu doesn't import the in-place version of 'AbsDiff' so we have to create two images for each iteration.
+ d_bg <- img.AbsDiff(Gray(float mean_bg))
+ d_fg <- img.AbsDiff(Gray(float mean_fg))
CvInvoke.Compare(d_fg, d_bg, fg, CvEnum.CmpType.LessThan)
- let mutable bg_total = 0.0
+ let mutable bg_total = 0.f
let mutable bg_nb = 0
- let mutable fg_total = 0.0
+ let mutable fg_total = 0.f
let mutable fg_nb = 0
for i in 0 .. h - 1 do
for j in 0 .. w - 1 do
- if fg.Data.[i, j, 0] > 0uy
+ if fgData.[i, j, 0] > 0uy
then
- fg_total <- fg_total + float img.Data.[i, j, 0]
+ fg_total <- fg_total + imgData.[i, j, 0]
fg_nb <- fg_nb + 1
else
- bg_total <- bg_total + float img.Data.[i, j, 0]
+ bg_total <- bg_total + imgData.[i, j, 0]
bg_nb <- bg_nb + 1
- mean_bg <- bg_total / float bg_nb
- mean_fg <- fg_total / float fg_nb
+ mean_bg <- bg_total / float32 bg_nb
+ mean_fg <- fg_total / float32 fg_nb
{ fg = fg; mean_bg = mean_bg; mean_fg = mean_fg; d_fg = d_fg }
\ No newline at end of file
median_fg: float
d_fg: Image<Gray, float32> } // Euclidean distances of the foreground to median_fg.
-let kmedians (img: Image<Gray, float32>) (fgFactor: float) : Result =
+let kmedians (img: Image<Gray, float32>) : Result =
let nbIteration = 3
let w = img.Width
let h = img.Height
open System
type I2DCoords =
- abstract X : float
- abstract Y : float
+ abstract X : float32
+ abstract Y : float32
// Compare 'e1' and 'e2' by X.
let cmpX (e1: I2DCoords) (e2: I2DCoords) : int =
match e1.X.CompareTo(e2.X) with
- | 0 -> match e1.Y.CompareTo(e2.Y) with
+ | 0 -> match e1.Y.CompareTo(e2.Y) with
| 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
| v -> v
| v -> v
-
+
// Compare 'e1' and 'e2' by Y.
let cmpY (e1: I2DCoords) (e2: I2DCoords) : int =
match e1.Y.CompareTo(e2.Y) with
- | 0 -> match e1.X.CompareTo(e2.X) with
+ | 0 -> match e1.X.CompareTo(e2.X) with
| 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
| v -> v
| v -> v
-type Region = { minX: float; maxX: float; minY: float; maxY: float } with
+type Region = { minX: float32; maxX: float32; minY: float32; maxY: float32 } with
member this.Contains px py : bool =
- px >= this.minX && px <= this.maxX &&
+ px >= this.minX && px <= this.maxX &&
py >= this.minY && py <= this.maxY
member this.IsSub otherRegion : bool =
- this.minX >= otherRegion.minX && this.maxX <= otherRegion.maxX &&
+ this.minX >= otherRegion.minX && this.maxX <= otherRegion.maxX &&
this.minY >= otherRegion.minY && this.maxY <= otherRegion.maxY
- member this.Intersects otherRegion : bool =
+ member this.Intersects otherRegion : bool =
this.minX < otherRegion.maxX && this.maxX >= otherRegion.minX &&
this.minY < otherRegion.maxY && this.maxY >= otherRegion.minY
-type Tree<'a when 'a :> I2DCoords> =
- | Node of float * Tree<'a> * Tree<'a>
+type Tree<'a when 'a :> I2DCoords> =
+ | Node of float32 * Tree<'a> * Tree<'a>
| Leaf of 'a
static member BuildTree (l: 'a list) : Tree<'a> =
let rec buildTreeFromSortedArray (pXSorted: 'a[]) (pYSorted: 'a[]) (depth: int) : Tree<'a> =
if pXSorted.Length = 1
- then
+ then
Leaf pXSorted.[0]
- else
+ else
if depth % 2 = 1 // 'depth' is odd -> vertical splitting else horizontal splitting.
- then
+ then
let leftX, rightX = Array.splitAt ((pXSorted.Length + 1) / 2) pXSorted
let splitElement = Array.last leftX
let leftY, rightY = Array.partition (fun (e: 'a) -> cmpX e splitElement <= 0) pYSorted // FIXME: Maybe this operation can be optimized.
match tree with
| Leaf v -> if searchRegion.Contains v.X v.Y then [v] else []
| Node (splitValue, part1, part2) ->
- let valuesInRegion (region: Region) (treeRegion: Tree<'a>) =
+ let valuesInRegion (region: Region) (treeRegion: Tree<'a>) =
if region.IsSub searchRegion
then
valuesFrom treeRegion
elif region.Intersects searchRegion
then
searchWithRegion treeRegion region (depth + 1)
- else
+ else
[]
if depth % 2 = 1 // Vertical splitting.
- then
- let leftRegion = { currentRegion with maxX = splitValue }
+ then
+ let leftRegion = { currentRegion with maxX = splitValue }
let rightRegion = { currentRegion with minX = splitValue }
(valuesInRegion leftRegion part1) @ (valuesInRegion rightRegion part2)
else // Horizontal splitting.
- let downRegion = { currentRegion with maxY = splitValue }
+ let downRegion = { currentRegion with maxY = splitValue }
let upRegion = { currentRegion with minY = splitValue }
(valuesInRegion downRegion part1) @ (valuesInRegion upRegion part2)
-
- searchWithRegion this { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1
-
-
+
+ searchWithRegion this { minX = Single.MinValue; maxX = Single.MaxValue; minY = Single.MinValue; maxY = Single.MaxValue } 1
+
+
///// Tests. TODO: to put in a unit test.
-type Point (x: float, y: float) =
+type Point (x: float32, y: float32) =
interface I2DCoords with
member this.X = x
member this.Y = y
sprintf "(%.1f, %.1f)" x y
// TODO: test with identical X or Y coords
-let test () =
+let test () =
let pts = [
- Point(1.0, 1.0)
- Point(2.0, 2.0)
- Point(1.5, 3.6)
- Point(3.0, 3.2)
- Point(4.0, 4.0)
- Point(3.5, 1.5)
- Point(2.5, 0.5) ]
-
+ Point(1.0f, 1.0f)
+ Point(2.0f, 2.0f)
+ Point(1.5f, 3.6f)
+ Point(3.0f, 3.2f)
+ Point(4.0f, 4.0f)
+ Point(3.5f, 1.5f)
+ Point(2.5f, 0.5f) ]
+
let tree = Tree.BuildTree pts
Utils.dprintfn "Tree: %A" tree
- let s1 = tree.Search { minX = 0.0; maxX = 5.0; minY = 0.0; maxY = 5.0 } // All points.
+ let s1 = tree.Search { minX = 0.0f; maxX = 5.0f; minY = 0.0f; maxY = 5.0f } // All points.
Utils.dprintfn "s1: %A" s1
- let s2 = tree.Search { minX = 2.8; maxX = 4.5; minY = 3.0; maxY = 4.5 }
+ let s2 = tree.Search { minX = 2.8f; maxX = 4.5f; minY = 3.0f; maxY = 4.5f }
Utils.dprintfn "s2: %A" s2
-
- let s3 = tree.Search { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 }
+
+ let s3 = tree.Search { minX = 2.0f; maxX = 2.0f; minY = 2.0f; maxY = 2.0f }
Utils.dprintfn "s3: %A" s3
-let test2 () =
+let test2 () =
let pts = [
- Point(1.0, 1.0)
- Point(1.0, 2.0)
- Point(1.0, 3.0) ]
-
+ Point(1.0f, 1.0f)
+ Point(1.0f, 2.0f)
+ Point(1.0f, 3.0f) ]
+
let tree = Tree.BuildTree pts
Utils.dprintfn "Tree: %A" tree
-
- let s1 = tree.Search { minX = 1.0; maxX = 1.0; minY = 1.0; maxY = 1.0 }
+
+ let s1 = tree.Search { minX = 1.0f; maxX = 1.0f; minY = 1.0f; maxY = 1.0f }
Utils.dprintfn "s1: %A" s1
- let s2 = tree.Search { minX = 1.0; maxX = 1.0; minY = 2.0; maxY = 2.0 }
+ let s2 = tree.Search { minX = 1.0f; maxX = 1.0f; minY = 2.0f; maxY = 2.0f }
Utils.dprintfn "s2: %A" s2
-
+
// This case result is wrong: FIXME
- let s3 = tree.Search { minX = 1.0; maxX = 1.0; minY = 3.0; maxY = 3.0 }
+ let s3 = tree.Search { minX = 1.0f; maxX = 1.0f; minY = 3.0f; maxY = 3.0f }
Utils.dprintfn "s3: %A" s3
- let s4 = tree.Search { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 }
+ let s4 = tree.Search { minX = 0.0f; maxX = 2.0f; minY = 0.0f; maxY = 4.0f }
Utils.dprintfn "s4: %A" s4
-
+
let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell list =
-
- use scaledImg = if config.Parameters.scale = 1.0 then img else img.Resize(config.Parameters.scale, CvEnum.Inter.Area)
+ let scaledImg = if config.Parameters.scale = 1.0 then img else img.Resize(config.Parameters.scale, CvEnum.Inter.Area)
use green = scaledImg.Item(1)
let greenFloat = green.Convert<Gray, float32>()
- let filteredGreen = gaussianFilter greenFloat config.Parameters.preFilterSigma
+ let filteredGreen = gaussianFilter greenFloat (float config.Parameters.preFilterSigma)
logTime "areaOpen 1" (fun () -> ImgTools.areaOpenF filteredGreen config.Parameters.initialAreaOpen)
- config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.3 |> float)
+ config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert<Gray, byte>()) (10, 100) 0.3 |> float32)
- let secondAreaOpen = int <| config.RBCArea / 3.
+ let secondAreaOpen = int <| config.RBCArea / 3.f
if secondAreaOpen > config.Parameters.initialAreaOpen
then
logTime "areaOpen 2" (fun () -> ImgTools.areaOpenF filteredGreen secondAreaOpen)
- let parasites, filteredGreenWhitoutInfection, filteredGreenWhitoutStain = ParasitesMarker.find filteredGreen config
+ let parasites, filteredGreenWhitoutStain = ParasitesMarker.find filteredGreen config
//let parasites, filteredGreenWhitoutInfection, filteredGreenWhitoutStain = ParasitesMarker.findMa greenFloat filteredGreenFloat config
let edges, xGradient, yGradient = logTime "Finding edges" (fun () -> ImgTools.findEdges filteredGreenWhitoutStain)
- logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12)
+ logTime "Removing small connected components from thinning" (fun () -> removeArea edges (config.RBCRadius ** 2.f / 50.f |> int))
let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
let matchingEllipses = Ellipse.find edges xGradient yGradient config
- matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses )
+ matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses)
let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config)
let buildFileName postfix = System.IO.Path.Combine(dirPath, name + postfix)
- //saveImg canny (buildFileName " - canny.png")
-
saveMat (edges * 255.0) (buildFileName " - edges.png")
saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
saveImg filteredGreen (buildFileName " - filtered.png")
saveImg filteredGreenWhitoutStain (buildFileName " - filtered closed stain.png")
- saveImg filteredGreenWhitoutInfection (buildFileName " - filtered closed infection.png")
+ //saveImg filteredGreenWhitoutInfection (buildFileName " - filtered closed infection.png")
saveImg green (buildFileName " - green.png")
// Do not take in account matching score below this when two ellipses are matched.
-let matchingScoreThreshold1 = 0.6
+[<Literal>]
+let matchingScoreThreshold1 = 0.6f
-let scaleOverlapTest = 0.8
+[<Literal>]
+let scaleOverlapTest = 0.8f
-type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
+type private EllipseScoreFlaggedKd (matchingScore: float32, e: Ellipse) =
let mutable matchingScore = matchingScore
member this.Ellipse = e
member this.MatchingScore = matchingScore
- member this.AddMatchingScore(score: float) =
+ member this.AddMatchingScore(score: float32) =
matchingScore <- matchingScore + score
member val Processed = false with get, set
member this.Y = this.Ellipse.Cy
-type MatchingEllipses (radiusMin: float) =
+type MatchingEllipses (radiusMin: float32) =
let ellipses = List<EllipseScoreFlaggedKd>()
// All ellipses with a score below this are removed.
- let matchingScoreThreshold2 = 20. * radiusMin
+ let matchingScoreThreshold2 = 20.f * radiusMin
member this.Add (e: Ellipse) =
- ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
+ ellipses.Add(EllipseScoreFlaggedKd(0.f, e))
member this.Ellipses : Ellipse list =
List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
for e in ellipses do
e.Processed <- true
let areaE = e.Ellipse.Area
- let window = { KdTree.minX = e.Ellipse.Cx - windowSize / 2.0
- KdTree.maxX = e.Ellipse.Cx + windowSize / 2.0
- KdTree.minY = e.Ellipse.Cy - windowSize / 2.0
- KdTree.maxY = e.Ellipse.Cy + windowSize / 2.0 }
+ let window = { KdTree.minX = e.Ellipse.Cx - windowSize / 2.f
+ KdTree.maxX = e.Ellipse.Cx + windowSize / 2.f
+ KdTree.minY = e.Ellipse.Cy - windowSize / 2.f
+ KdTree.maxY = e.Ellipse.Cy + windowSize / 2.f }
for other in tree.Search window do
if not other.Processed
then
let areaOther = other.Ellipse.Area
match EEOver.EEOverlapArea e.Ellipse other.Ellipse with
| Some (commonArea, _, _) ->
- let matchingScore = 2.0 * commonArea / (areaE + areaOther)
+ let matchingScore = 2.f * commonArea / (areaE + areaOther)
if matchingScore >= matchingScoreThreshold1
then
other.AddMatchingScore(matchingScore * e.Ellipse.Perimeter)
ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore))
// 4) Remove ellipses with a low score.
- let i = ellipses.BinarySearch(EllipseScoreFlaggedKd(matchingScoreThreshold2, Ellipse(0.0, 0.0, 0.0, 0.0, 0.0)),
+ let i = ellipses.BinarySearch(EllipseScoreFlaggedKd(matchingScoreThreshold2, Ellipse(0.f, 0.f, 0.f, 0.f, 0.f)),
{ new IComparer<EllipseScoreFlaggedKd> with
member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs
let nbToRemove = ellipses.Count - i
<Resource Include="MainWindow.xaml" />
<Compile Include="MainWindow.xaml.fs" />
<Compile Include="Heap.fs" />
+ <Compile Include="Const.fs" />
<Compile Include="Types.fs" />
+ <Compile Include="EEOver.fs" />
<Compile Include="Utils.fs" />
<Compile Include="ImgTools.fs" />
<Compile Include="Granulometry.fs" />
<Compile Include="Config.fs" />
<Compile Include="KMedians.fs" />
<Compile Include="KMeans.fs" />
- <Compile Include="EEOver.fs" />
<Compile Include="ParasitesMarker.fs" />
<Compile Include="KdTree.fs" />
<Compile Include="MatchingEllipses.fs" />
<Reference Include="Emgu.Util">
<HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll</HintPath>
</Reference>
- <Reference Include="FSharp.Core, Version=$(TargetFSharpCoreVersion), Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a">
+ <Reference Include="FSharp.Collections.ParallelSeq">
+ <HintPath>..\packages\FSharp.Collections.ParallelSeq.1.0.2\lib\net40\FSharp.Collections.ParallelSeq.dll</HintPath>
+ <Private>True</Private>
+ </Reference>
+ <Reference Include="FSharp.Core">
+ <HintPath>..\packages\FSharp.Core.4.0.0.1\lib\net40\FSharp.Core.dll</HintPath>
+ <Private>True</Private>
</Reference>
- <Reference Include="FSharp.ViewModule.Core.Wpf">
- <HintPath>..\packages\FSharp.ViewModule.Core.0.9.9.1\lib\net45\FSharp.ViewModule.Core.Wpf.dll</HintPath>
+ <Reference Include="FSharp.ViewModule">
+ <HintPath>..\packages\FSharp.ViewModule.Core.0.9.9.2\lib\portable-net45+netcore45+wpa81+wp8+MonoAndroid1+MonoTouch1\FSharp.ViewModule.dll</HintPath>
<Private>True</Private>
</Reference>
<Reference Include="FsXaml.Wpf">
<Private>True</Private>
</Reference>
<Reference Include="MathNet.Numerics">
- <HintPath>..\packages\MathNet.Numerics.3.9.0\lib\net40\MathNet.Numerics.dll</HintPath>
+ <HintPath>..\packages\MathNet.Numerics.3.10.0\lib\net40\MathNet.Numerics.dll</HintPath>
<Private>True</Private>
</Reference>
<Reference Include="MathNet.Numerics.FSharp">
- <HintPath>..\packages\MathNet.Numerics.FSharp.3.9.0\lib\net40\MathNet.Numerics.FSharp.dll</HintPath>
+ <HintPath>..\packages\MathNet.Numerics.FSharp.3.10.0\lib\net40\MathNet.Numerics.FSharp.dll</HintPath>
<Private>True</Private>
</Reference>
<Reference Include="mscorlib" />
let findMa (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (config: Config.Config) : Result * Image<Gray, byte> * Image<Gray, byte> =
// We use the filtered image to find the dark stain.
- let kmediansResults = logTime "Finding fg/bg (k-medians)" (fun () -> KMedians.kmedians filteredGreen 1.0)
+ let kmediansResults = logTime "Finding fg/bg (k-medians)" (fun () -> KMedians.kmedians filteredGreen)
let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResults
- let darkStain = d_fg.Cmp(median_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
+ let darkStain = d_fg.Cmp(median_bg * float config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
darkStain._And(filteredGreen.Cmp(median_fg, CvEnum.CmpType.LessThan))
darkStain._And(fg)
// * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas.
// * 'Stain' corresponds to the stain around the parasites.
// * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes.
-let find (filteredGreen: Image<Gray, float32>) (config: Config.Config) : Result * Image<Gray, float32> * Image<Gray, float32> =
+let find (filteredGreen: Image<Gray, float32>) (config: Config.Config) : Result * Image<Gray, float32> =
- let filteredGreenWithoutInfection = filteredGreen.Copy()
+ use filteredGreenWithoutInfection = filteredGreen.Copy()
ImgTools.areaCloseF filteredGreenWithoutInfection (int config.InfectionArea)
let filteredGreenWithoutStain = filteredGreenWithoutInfection.Copy()
// We use the filtered image to find the dark stain.
// With K-Means.
- let kmeansResults = logTime "Finding fg/bg (k-means)" (fun () -> KMeans.kmeans (filteredGreenWithoutInfection) 1.0)
+ let kmeansResults = logTime "Finding fg/bg (k-means)" (fun () -> KMeans.kmeans (filteredGreenWithoutInfection))
let { KMeans.mean_bg = value_bg; KMeans.mean_fg = value_fg; KMeans.d_fg = d_fg } = kmeansResults
// With K-Medians.
- (* let kmediansResults = logTime "Finding fg/bg (k-medians)" (fun () -> KMedians.kmedians (filteredGreenWithoutInfection) 1.0) // FIXME: avoid converting this again in MainAnalysis
+ (* let kmediansResults = logTime "Finding fg/bg (k-medians)" (fun () -> KMedians.kmedians (filteredGreenWithoutInfection)) // FIXME: avoid converting this again in MainAnalysis
let { KMedians.median_bg = value_bg; KMedians.median_fg = value_fg; KMedians.d_fg = d_fg } = kmediansResults *)
- let darkStain = d_fg.Cmp(value_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
- darkStain._And(filteredGreenWithoutInfection.Cmp(value_fg, CvEnum.CmpType.LessThan))
+ let darkStain = d_fg.Cmp((float value_bg) * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan)
+ darkStain._And(filteredGreenWithoutInfection.Cmp(float value_fg, CvEnum.CmpType.LessThan))
let marker (img: Image<Gray, float32>) (closed: Image<Gray, float32>) (level: float) : Image<Gray, byte> =
- let diff = closed - (img * level)
+ let diff = img.Copy() // closed - (img * level)
+ diff._Mul(level)
+ CvInvoke.Subtract(closed, diff, diff)
diff._ThresholdBinary(Gray(0.0), Gray(255.))
diff.Convert<Gray, byte>()
{ darkStain = darkStain
infection = infectionMarker
stain = stainMarker },
- filteredGreenWithoutInfection,
filteredGreenWithoutStain
open System.Windows.Controls
open System.Drawing
open System.Diagnostics
+open System.Threading
+
+open FSharp.Collections.ParallelSeq
open Emgu.CV
open Emgu.CV.Structure
initialAreaOpen = 2000
- minRbcRadius = -0.32
- maxRbcRadius = 0.32
+ minRbcRadius = -0.32f
+ maxRbcRadius = 0.32f
preFilterSigma = 1.7 // 1.5
factorNbPick = 1.0
- factorWindowSize = 2.0
- darkStainLevel = 0.25 // Lower -> more sensitive. 0.3
+ darkStainLevel = 0.22 // Lower -> more sensitive. 0.3. Careful about illumination on the borders.
maxDarkStainRatio = 0.1 // 10 %
- infectionArea = 0.012 // 1.2 %
+ infectionArea = 0.012f // 1.2 %
infectionLevel = 1.12 // Lower -> more sensitive.
- stainArea = 0.08 // 8 %
+ stainArea = 0.08f // 8 %
stainLevel = 1.1 // Lower -> more sensitive.
maxStainRatio = 0.12 // 12 %
standardDeviationMaxRatio = 0.5 // 0.55
- minimumCellArea = 0.5 })
+ minimumCellArea = 0.5f })
match mode with
| CmdLine (input, output) ->
use resultFile = new StreamWriter(new FileStream(Path.Combine(output, "results.txt"), FileMode.Append, FileAccess.Write))
- for file in files do
- try
- use img = new Image<Bgr, byte>(file)
- Utils.log (sprintf "== File: %A" file)
+ //try
+ let images = seq { for file in files -> Path.GetFileNameWithoutExtension(FileInfo(file).Name), new Image<Bgr, byte>(file) }
+
+ let nbConcurrentTaskLimit = 4
+ let n = Environment.ProcessorCount
- let cells = Utils.logTime "Whole analyze" (fun () -> ImageAnalysis.doAnalysis img (Path.GetFileNameWithoutExtension(FileInfo(file).Name)) config)
+ Utils.logTime "Whole analyze" (fun () ->
+ let results =
+ images
+ |> PSeq.map (fun (id, img) -> id, ImageAnalysis.doAnalysis img id (config.Copy()))
+ |> PSeq.withDegreeOfParallelism (if n > nbConcurrentTaskLimit then nbConcurrentTaskLimit else n)
+
+ for id, cells in results do
let total, infected = Utils.countCells cells
- fprintf resultFile "File: %s %d %d %.2f\n" file total infected (100. * (float infected) / (float total))
- with
- | :? IOException as ex -> Utils.log (sprintf "Unable to open the image '%A': %A" file ex)
+ fprintf resultFile "File: %s %d %d %.2f\n" id total infected (100. * (float infected) / (float total)))
+
+ //Utils.log (sprintf "== File: %A" file)
+ //with
+ //| :? IOException as ex -> Utils.log (sprintf "Unable to open the image '%A': %A" file ex)
0
| Window ->
open Emgu.CV
open Emgu.CV.Structure
-type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) =
+open Const
+
+type Ellipse (cx: float32, cy: float32, a: float32, b: float32, alpha: float32) =
member this.Cx = cx
member this.Cy = cy
member this.A = a
member this.B = b
member this.Alpha = alpha
- member this.Area = a * b * Math.PI
+ member this.Area = a * b * PI
// Does the ellipse contain the point (x, y)?.
member this.Contains x y =
- ((x - cx) * cos alpha + (y - cy) * sin alpha) ** 2.0 / a ** 2.0 + ((x - cx) * sin alpha - (y - cy) * cos alpha) ** 2.0 / b ** 2.0 <= 1.0
+ ((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
- member this.CutAVericalLine (y: float) : bool =
- a ** 2.0 + b ** 2.0 - 2.0 * y ** 2.0 + 4.0 * y * cx - 2.0 * cx ** 2.0 + a ** 2.0 * cos (2.0 * alpha) - b ** 2.0 * cos (2.0 * alpha) > 0.0
+ member this.CutAVericalLine (y: float32) : bool =
+ a ** 2.f + b ** 2.f - 2.f * y ** 2.f + 4.f * y * cx - 2.f * cx ** 2.f + a ** 2.f * cos (2.f * alpha) - b ** 2.f * cos (2.f * alpha) > 0.f
- member this.CutAnHorizontalLine (x: float) : bool =
- a ** 2.0 + b ** 2.0 - 2.0 * x ** 2.0 + 4.0 * x * cy - 2.0 * cy ** 2.0 - a ** 2.0 * cos (2.0 * alpha) + b ** 2.0 * cos (2.0 * alpha) > 0.0
+ member this.CutAnHorizontalLine (x: float32) : bool =
+ a ** 2.f + b ** 2.f - 2.f * x ** 2.f + 4.f * x * cy - 2.f * cy ** 2.f - a ** 2.f * cos (2.f * alpha) + b ** 2.f * cos (2.f * alpha) > 0.f
- member this.isOutside (width: float) (height: float) =
- this.Cx < 0.0 || this.Cx >= width ||
- this.Cy < 0.0 || this.Cy >= height ||
- this.CutAVericalLine 0.0 || this.CutAVericalLine width ||
- this.CutAnHorizontalLine 0.0 || this.CutAnHorizontalLine height
+ member this.isOutside (width: float32) (height: float32) =
+ this.Cx < 0.f || this.Cx >= width ||
+ this.Cy < 0.f || this.Cy >= height ||
+ this.CutAVericalLine 0.f || this.CutAVericalLine width ||
+ this.CutAnHorizontalLine 0.f || this.CutAnHorizontalLine height
- member this.Scale (factor: float) =
+ member this.Scale (factor: float32) =
Ellipse(this.Cx, this.Cy, this.A * factor, this.B * factor, alpha)
// Approximation of Ramanujan.
member this.Perimeter =
- Math.PI * (3.0 * (this.A + this.B) - sqrt ((3.0 * this.A + this.B) * (this.A + 3.0 * this.B)))
+ PI * (3.f * (this.A + this.B) - sqrt ((3.f * this.A + this.B) * (this.A + 3.f * this.B)))
override this.ToString () =
sprintf "(cx: %A, cy: %A, a: %A, b: %A, alpha: %A)" this.Cx this.Cy this.A this.B this.Alpha
elements: Matrix<byte> }
[<Struct>]
-type Line (a: float, b: float) =
+type Line (a: float32, b: float32) =
member this.A = a
member this.B = b
- member this.Valid = not (Double.IsInfinity this.A)
+ member this.Valid = not (Single.IsInfinity this.A)
[<Struct>]
-type PointD (x: float, y: float) =
+type PointD (x: float32, y: float32) =
member this.X = x
member this.Y = y
open Types
-let roundInt = int << round
+let inline roundInt v = v |> round |> int
let inline dprintfn fmt =
Printf.ksprintf System.Diagnostics.Debug.WriteLine fmt
sign (p.X - p1.X) <> sign (p.X - p2.X)
let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) =
- (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0
+ (p1.X - p2.X) ** 2.f + (p1.Y - p2.Y) ** 2.f
let distanceTwoPoints (p1: PointD) (p2: PointD) =
squaredDistanceTwoPoints p1 p2 |> sqrt
<packages>
<package id="Castle.Core" version="1.1.0" targetFramework="net46" />
<package id="Expression.Blend.Sdk" version="1.0.2" targetFramework="net46" />
- <package id="FSharp.Core" version="3.1.2.5" targetFramework="net46" />
- <package id="FSharp.ViewModule.Core" version="0.9.9.1" targetFramework="net46" />
+ <package id="FSharp.Collections.ParallelSeq" version="1.0.2" targetFramework="net461" />
+ <package id="FSharp.Core" version="4.0.0.1" targetFramework="net461" />
+ <package id="FSharp.ViewModule.Core" version="0.9.9.2" targetFramework="net461" />
<package id="FsXaml.Wpf" version="0.9.9" targetFramework="net46" />
<package id="log4net" version="1.2.10" targetFramework="net46" />
- <package id="MathNet.Numerics" version="3.9.0" targetFramework="net461" />
- <package id="MathNet.Numerics.FSharp" version="3.9.0" targetFramework="net461" />
+ <package id="MathNet.Numerics" version="3.10.0" targetFramework="net461" />
+ <package id="MathNet.Numerics.FSharp" version="3.10.0" targetFramework="net461" />
</packages>
\ No newline at end of file
projects / master-thesis.git / commitdiff