From 044b0ae69df3ac565432545b2fa934589016f9bd Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Tue, 5 Jan 2016 11:42:12 +0100 Subject: [PATCH] Use float32 to reduce memory footprint. --- Parasitemia/Parasitemia/Classifier.fs | 41 +++--- Parasitemia/Parasitemia/Config.fs | 20 +-- Parasitemia/Parasitemia/EEOver.fs | 18 +-- Parasitemia/Parasitemia/Ellipse.fs | 134 ++++++++++---------- Parasitemia/Parasitemia/Granulometry.fs | 2 +- Parasitemia/Parasitemia/ImgTools.fs | 134 +++++++++++--------- Parasitemia/Parasitemia/KMeans.fs | 46 ++++--- Parasitemia/Parasitemia/KMedians.fs | 2 +- Parasitemia/Parasitemia/KdTree.fs | 98 +++++++------- Parasitemia/Parasitemia/MainAnalysis.fs | 19 ++- Parasitemia/Parasitemia/MatchingEllipses.fs | 28 ++-- Parasitemia/Parasitemia/Parasitemia.fsproj | 19 ++- Parasitemia/Parasitemia/ParasitesMarker.fs | 21 +-- Parasitemia/Parasitemia/Program.fs | 41 +++--- Parasitemia/Parasitemia/Types.fs | 36 +++--- Parasitemia/Parasitemia/Utils.fs | 4 +- Parasitemia/Parasitemia/packages.config | 9 +- 17 files changed, 364 insertions(+), 308 deletions(-) diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index 0d21375..f0c9d01 100644 --- a/Parasitemia/Parasitemia/Classifier.fs +++ b/Parasitemia/Parasitemia/Classifier.fs @@ -37,13 +37,13 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: // 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. @@ -79,13 +79,17 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: 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 [] @@ -101,24 +105,24 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: // 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 @@ -133,7 +137,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: 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 @@ -159,7 +163,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: let elements = new Matrix(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 @@ -179,8 +183,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: 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 diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index 7f841c4..aa8f0e5 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -2,6 +2,8 @@ open System +open Const + type Debug = | DebugOff | DebugOn of string // Output directory. @@ -11,40 +13,42 @@ type Parameters = { 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 diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index 528632a..5a3c9e4 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -514,9 +514,9 @@ let private biquadroots (p: float[]) (r: float[,]) = 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 @@ -718,11 +718,11 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float * float[] * f | 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 diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index 73771be..ef167b2 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -2,6 +2,7 @@ open System open System.Collections.Generic +open System.Drawing open Emgu.CV open Emgu.CV.Structure @@ -9,12 +10,12 @@ 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) @@ -40,12 +41,12 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) 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 @@ -60,27 +61,27 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: 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 @@ -92,7 +93,7 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: 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) @@ -101,11 +102,11 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: 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) @@ -116,27 +117,27 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: 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 @@ -147,32 +148,32 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: 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 @@ -191,12 +192,12 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) 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 @@ -204,27 +205,32 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) let find (edges: Matrix) - (xGradient: Image) - (yGradient: Image) + (xGradient: Image) + (yGradient: Image) (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() let edgesData = edges.Data let xDirData = xGradient.Data @@ -243,40 +249,40 @@ let find (edges: Matrix) 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 | _ -> () | _ -> () diff --git a/Parasitemia/Parasitemia/Granulometry.fs b/Parasitemia/Parasitemia/Granulometry.fs index b18f53c..0ef862b 100644 --- a/Parasitemia/Parasitemia/Granulometry.fs +++ b/Parasitemia/Parasitemia/Granulometry.fs @@ -11,7 +11,7 @@ open Utils // 'range': a minimum and maximum radius. // 'scale': <= 1.0, to speed up the process. let findRadius (img: Image) (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) diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index b9a31f8..f64b1c3 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -8,17 +8,18 @@ open System.Linq 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) : Image = +let normalizeAndConvert (img: Image) : Image = 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() + ((img.Convert() - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = @@ -46,7 +47,7 @@ let suppressMConnections (img: Matrix) = img.[i, j] <- 0uy -let findEdges (img: Image) : Matrix * Image * Image = +let findEdges (img: Image) : Matrix * Image * Image = let w = img.Width let h = img.Height @@ -55,103 +56,112 @@ let findEdges (img: Image) : Matrix * Image * [ 2.0f; 0.0f; -2.0f ] [ 1.0f; 0.0f; -1.0f ]], Point(1, 1)) - let xGradient = img.Convolution(sobelKernel).Convert() - let yGradient = img.Convolution(sobelKernel.Transpose()).Convert() + 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(xGradient.Size) - use angles = new Matrix(xGradient.Size) + use magnitudes = new Matrix(xGradient.Size) + use angles = new Matrix(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() - 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(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(xGradient.Size) + let edgesData = edges.Data - // Histeresis thresholding. + // Hysteresis thresholding. let toVisit = Stack() 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 @@ -723,7 +733,7 @@ let removeArea (mat: Matrix) (areaSize: int) = ( 0, -1) // p8 (-1, -1) |] // p9 - let mat' = new Matrix(mat.Size) + use mat' = new Matrix(mat.Size) let w = mat'.Width let h = mat'.Height mat.CopyTo(mat') @@ -791,21 +801,21 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo 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 diff --git a/Parasitemia/Parasitemia/KMeans.fs b/Parasitemia/Parasitemia/KMeans.fs index 7c18120..15651ae 100644 --- a/Parasitemia/Parasitemia/KMeans.fs +++ b/Parasitemia/Parasitemia/KMeans.fs @@ -8,11 +8,11 @@ open Emgu.CV.Structure type Result = { fg: Image - mean_bg: float - mean_fg: float + mean_bg: float32 + mean_fg: float32 d_fg: Image } // Euclidean distances of the foreground to mean_fg. -let kmeans (img: Image) (fgFactor: float) : Result = +let kmeans (img: Image) : Result = let nbIteration = 3 let w = img.Width let h = img.Height @@ -23,35 +23,47 @@ let kmeans (img: Image) (fgFactor: float) : Result = 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(img.Size) - let mutable d_fg = new Image(img.Size) - let mutable fg = new Image(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 = null + let mutable d_fg : Image = null + let fg = new Image(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 diff --git a/Parasitemia/Parasitemia/KMedians.fs b/Parasitemia/Parasitemia/KMedians.fs index 8a3fd2b..f7f2e54 100644 --- a/Parasitemia/Parasitemia/KMedians.fs +++ b/Parasitemia/Parasitemia/KMedians.fs @@ -12,7 +12,7 @@ type Result = { median_fg: float d_fg: Image } // Euclidean distances of the foreground to median_fg. -let kmedians (img: Image) (fgFactor: float) : Result = +let kmedians (img: Image) : Result = let nbIteration = 3 let w = img.Width let h = img.Height diff --git a/Parasitemia/Parasitemia/KdTree.fs b/Parasitemia/Parasitemia/KdTree.fs index 3e714ac..7da0255 100644 --- a/Parasitemia/Parasitemia/KdTree.fs +++ b/Parasitemia/Parasitemia/KdTree.fs @@ -3,40 +3,40 @@ 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> = @@ -47,11 +47,11 @@ type Tree<'a when 'a :> I2DCoords> = 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. @@ -74,32 +74,32 @@ type Tree<'a when 'a :> I2DCoords> = 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 @@ -108,47 +108,47 @@ type Point (x: float, y: float) = 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 - + diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index fdae997..ec2aabb 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -13,31 +13,30 @@ open Types let doAnalysis (img: Image) (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() - 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()) (10, 100) 0.3 |> float) + config.RBCRadius <- logTime "Granulometry" (fun() -> Granulometry.findRadius (filteredGreen.Convert()) (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) @@ -49,8 +48,6 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li 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") @@ -80,7 +77,7 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li 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") diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 2022c26..6e57199 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -10,18 +10,20 @@ open Utils // Do not take in account matching score below this when two ellipses are matched. -let matchingScoreThreshold1 = 0.6 +[] +let matchingScoreThreshold1 = 0.6f -let scaleOverlapTest = 0.8 +[] +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 @@ -32,14 +34,14 @@ type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) = member this.Y = this.Ellipse.Cy -type MatchingEllipses (radiusMin: float) = +type MatchingEllipses (radiusMin: float32) = let ellipses = List() // 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) @@ -58,17 +60,17 @@ type MatchingEllipses (radiusMin: float) = 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) @@ -79,7 +81,7 @@ type MatchingEllipses (radiusMin: float) = 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 with member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs let nbToRemove = ellipses.Count - i diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj index b4e3645..14e1c11 100644 --- a/Parasitemia/Parasitemia/Parasitemia.fsproj +++ b/Parasitemia/Parasitemia/Parasitemia.fsproj @@ -65,14 +65,15 @@ + + - @@ -90,10 +91,16 @@ ..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll - + + ..\packages\FSharp.Collections.ParallelSeq.1.0.2\lib\net40\FSharp.Collections.ParallelSeq.dll + True + + + ..\packages\FSharp.Core.4.0.0.1\lib\net40\FSharp.Core.dll + True - - ..\packages\FSharp.ViewModule.Core.0.9.9.1\lib\net45\FSharp.ViewModule.Core.Wpf.dll + + ..\packages\FSharp.ViewModule.Core.0.9.9.2\lib\portable-net45+netcore45+wpa81+wp8+MonoAndroid1+MonoTouch1\FSharp.ViewModule.dll True @@ -109,11 +116,11 @@ True - ..\packages\MathNet.Numerics.3.9.0\lib\net40\MathNet.Numerics.dll + ..\packages\MathNet.Numerics.3.10.0\lib\net40\MathNet.Numerics.dll True - ..\packages\MathNet.Numerics.FSharp.3.9.0\lib\net40\MathNet.Numerics.FSharp.dll + ..\packages\MathNet.Numerics.FSharp.3.10.0\lib\net40\MathNet.Numerics.FSharp.dll True diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs index 71c15a9..f1a94e1 100644 --- a/Parasitemia/Parasitemia/ParasitesMarker.fs +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -20,9 +20,9 @@ type Result = { let findMa (green: Image) (filteredGreen: Image) (config: Config.Config) : Result * Image * Image = // 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) @@ -50,9 +50,9 @@ let findMa (green: Image) (filteredGreen: Image) ( // * '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) (config: Config.Config) : Result * Image * Image = +let find (filteredGreen: Image) (config: Config.Config) : Result * Image = - let filteredGreenWithoutInfection = filteredGreen.Copy() + use filteredGreenWithoutInfection = filteredGreen.Copy() ImgTools.areaCloseF filteredGreenWithoutInfection (int config.InfectionArea) let filteredGreenWithoutStain = filteredGreenWithoutInfection.Copy() @@ -61,18 +61,20 @@ let find (filteredGreen: Image) (config: Config.Config) : Result // 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) (closed: Image) (level: float) : Image = - 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() @@ -82,7 +84,6 @@ let find (filteredGreen: Image) (config: Config.Config) : Result { darkStain = darkStain infection = infectionMarker stain = stainMarker }, - filteredGreenWithoutInfection, filteredGreenWithoutStain diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 6531da8..09f4312 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -9,6 +9,9 @@ open System.Windows.Shapes open System.Windows.Controls open System.Drawing open System.Diagnostics +open System.Threading + +open FSharp.Collections.ParallelSeq open Emgu.CV open Emgu.CV.Structure @@ -63,26 +66,25 @@ let main args = 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) -> @@ -102,16 +104,25 @@ let main args = 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(file) - Utils.log (sprintf "== File: %A" file) + //try + let images = seq { for file in files -> Path.GetFileNameWithoutExtension(FileInfo(file).Name), new Image(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 -> diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs index e27f7d0..6240519 100644 --- a/Parasitemia/Parasitemia/Types.fs +++ b/Parasitemia/Parasitemia/Types.fs @@ -6,36 +6,38 @@ open System.Drawing 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 @@ -49,13 +51,13 @@ type Cell = { elements: Matrix } [] -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) [] -type PointD (x: float, y: float) = +type PointD (x: float32, y: float32) = member this.X = x member this.Y = y diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs index 58e914f..b30e9ac 100644 --- a/Parasitemia/Parasitemia/Utils.fs +++ b/Parasitemia/Parasitemia/Utils.fs @@ -4,7 +4,7 @@ open System.Diagnostics open Types -let roundInt = int << round +let inline roundInt v = v |> round |> int let inline dprintfn fmt = Printf.ksprintf System.Diagnostics.Debug.WriteLine fmt @@ -35,7 +35,7 @@ let inline linePassThroughSegment (l: Line) (p1: PointD) (p2: PointD) : bool = 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 diff --git a/Parasitemia/Parasitemia/packages.config b/Parasitemia/Parasitemia/packages.config index 20de8bc..d6eabef 100644 --- a/Parasitemia/Parasitemia/packages.config +++ b/Parasitemia/Parasitemia/packages.config @@ -2,10 +2,11 @@ - - + + + - - + + \ No newline at end of file -- 2.43.0