From 81d1b86719a1ebaf649c1de4c1364603155a53e1 Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Sun, 17 Jan 2016 18:00:46 +0100 Subject: [PATCH] * Add an exact method to compute an ellipse from three points and two tangents. * Simplification of the matching ellipses process. --- Parasitemia/Parasitemia/Config.fs | 1 - Parasitemia/Parasitemia/Const.fs | 1 - Parasitemia/Parasitemia/EEOver.fs | 10 +-- Parasitemia/Parasitemia/Ellipse.fs | 63 +++++++++++++++++- .../Parasitemia/GUI/AnalysisWindow.xaml | 2 +- Parasitemia/Parasitemia/GUI/GUI.fs | 24 +++++-- Parasitemia/Parasitemia/GUI/MainWindow.xaml | 4 +- Parasitemia/Parasitemia/Granulometry.fs | 1 - Parasitemia/Parasitemia/ImgTools.fs | 20 ------ Parasitemia/Parasitemia/KMeans.fs | 1 - Parasitemia/Parasitemia/KdTree.fs | 1 - Parasitemia/Parasitemia/MainAnalysis.fs | 22 +++---- Parasitemia/Parasitemia/MatchingEllipses.fs | 64 ++++++------------- Parasitemia/Parasitemia/Parasitemia.fsproj | 6 ++ Parasitemia/Parasitemia/ParasitesMarker.fs | 2 - Parasitemia/Parasitemia/Program.fs | 17 +++-- Parasitemia/Parasitemia/Types.fs | 1 - Parasitemia/Parasitemia/Utils.fs | 2 +- 18 files changed, 133 insertions(+), 109 deletions(-) diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index cd8f659..e7deb61 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -85,7 +85,6 @@ type RBCRadius (radius: float32, parameters: Parameters) = override this.ToString() = sprintf "%d px (%.1f μm)" (Utils.roundInt <| 2.f * radius) (2. * this.μm) - type Config (param: Parameters) = let RBCadiusInPixels (rbcDiameter: float<μm>) (resolution: float) : float32 = let rbcRadiusInch: float = (μmToInch rbcDiameter) / 2. diff --git a/Parasitemia/Parasitemia/Const.fs b/Parasitemia/Parasitemia/Const.fs index 0d2e52e..304d3a9 100644 --- a/Parasitemia/Parasitemia/Const.fs +++ b/Parasitemia/Parasitemia/Const.fs @@ -1,4 +1,3 @@ module Const -// TODO: try with a literal and check performance. let PI = float32 System.Math.PI \ No newline at end of file diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index 5a3c9e4..2bd13b2 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -69,7 +69,6 @@ let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (b then TANGENT_POINT else INTERSECTION_POINT - let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) = if abs x.[0] > a1 then x.[0] <- if x.[0] < 0.0 then -a1 else a1 @@ -181,7 +180,6 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: area1 + area2 - let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) (phi_1: float) (a2: float) (b2: float) (h2_tr: float) (k2_tr: float) (phi_2: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float = let mutable tanpts = 0 let mutable tanindex = 0 @@ -218,7 +216,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( let theta = Array.zeroCreate 4 - for i in 0..3 do + for i in 0 .. 3 do if abs xint.[i] > a1 then xint.[i] <- if xint.[i] < 0.0 then -a1 else a1 @@ -229,7 +227,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] #endif - for j in 1..3 do + for j in 1 .. 3 do let tmp0 = theta.[j] let tmp1 = xint.[j] let tmp2 = yint.[j] @@ -334,7 +332,6 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( area1 + area2 + area3 + area4 + area5 - let private quadroots (p: float[]) (r: float[,]) = let mutable b = -p.[1] / (2.0 * p.[0]) let c = p.[2] / p.[0] @@ -429,7 +426,6 @@ let private cubicroots (p: float[]) (r: float[,]) = for k in 1..3 do r.[2, k] <- 0.0 - let private biquadroots (p: float[]) (r: float[,]) = if p.[0] <> 1.0 then @@ -446,7 +442,7 @@ let private biquadroots (p: float[]) (r: float[,]) = a <- a - d let mutable quadExecuted = false - let quad () = + let inline quad () = if not quadExecuted then p.[2] <- c / b diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index be35940..ac3c041 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -4,6 +4,8 @@ open System open System.Collections.Generic open System.Drawing +open MathNet.Numerics.LinearAlgebra + open Emgu.CV open Emgu.CV.Structure @@ -147,6 +149,61 @@ let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: else None +let ellipse2 (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option = + let p0 = pointFromTwoLines (Types.Line(float32 m1, float32 (p1y - m1 * p1x))) (Types.Line(float32 m2, float32(p2y - m2 * p2x))) + let p0x, p0y = float p0.X, float p0.Y + + let s = matrix [[ 1.; 0.; 0. ] + [ 0.; 0.; -0.5 ] + [ 0.; -0.5; 0. ]] + + let v0 = matrix [[ 1.; p0x; p0y ]] + let v1 = matrix [[ 1.; p1x; p1y ]] + let v2 = matrix [[ 1.; p2x; p2y ]] + let v3 = matrix [[ 1.; p3x; p3y ]] + + let p = (v3.Stack(v1).Stack(v2).Determinant() * v0).Stack(v0.Stack(v3).Stack(v2).Determinant() * v1).Stack(v0.Stack(v1).Stack(v3).Determinant() * v2).Transpose() + let conicMat = p * s.Inverse() * p.Transpose() + let a = conicMat.[0, 0] + let b = conicMat.[0, 1] + let c = conicMat.[1, 1] + let d = conicMat.[0, 2] + let e = conicMat.[1, 2] + let f = conicMat.[2, 2] + + // Center. + let cx = b / a + let cy = d / a + + let at = c * f - e ** 2. + (e * d - b * f) * cx + (b * e - c * d) * cy + if at = 0. + then + None + else + let q = (-1. / at) * (matrix [[ a * f - d ** 2.0; b * d - a * e ]; [ b * d - a * e; a * c - b ** 2.0 ]]) + let eigen = q.Evd() + let eigenValues = eigen.EigenValues + let lambda = eigenValues.[1].Real + let mu = eigenValues.[0].Real + + if lambda <= 0. || mu <= 0. + then + None + else + let r1, r2 = 1. / (sqrt lambda), 1. / (sqrt mu) + + let eigenVectors = eigen.EigenVectors + let v_a = eigenVectors.[0, 0] + let v_b = eigenVectors.[1, 0] // [0, 1] + + // Angle against the longest axis. + let phi = (if r2 > r1 then atan (v_b / v_a) else atan (v_a / v_b)) + + let phi' = if phi < 0. then phi + Math.PI else phi + let majorAxis, minorAxis = if r1 > r2 then r1, r2 else r2, r1 + + Some (Types.Ellipse(float32 cx, float32 cy, float32 majorAxis, float32 minorAxis, float32 phi')) + let private vectorRotation (p1x: float32) (p1y: float32) (v1x: float32) (v1y: float32) (px: float32) (py: float32) : float32 = let mutable rotation = 1.f @@ -172,7 +229,6 @@ let private vectorRotation (p1x: float32) (p1y: float32) (v1x: float32) (v1y: fl rotation <- -1.f rotation - 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 @@ -238,7 +294,7 @@ let find (edges: Matrix) let rng = Random(42) - let ellipses = MatchingEllipses(r1) + let ellipses = MatchingEllipses(config.RBCRadius.Pixel) for window_i in -windowSize + increment .. increment .. h - increment do for window_j in -windowSize + increment .. increment .. w - increment do @@ -280,7 +336,8 @@ let find (edges: Matrix) then match areVectorsValid (float32 p1xf) (float32 p1yf) (float32 p2xf) (float32 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 (float m1) p2xf p2yf (float m2) p3xf p3yf with + //let pouet = ellipse2 p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf + match ellipse2 p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with | 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 -> ellipses.Add e diff --git a/Parasitemia/Parasitemia/GUI/AnalysisWindow.xaml b/Parasitemia/Parasitemia/GUI/AnalysisWindow.xaml index 7684ad9..599d581 100644 --- a/Parasitemia/Parasitemia/GUI/AnalysisWindow.xaml +++ b/Parasitemia/Parasitemia/GUI/AnalysisWindow.xaml @@ -3,7 +3,7 @@ xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" mc:Ignorable="d" - x:Name="AnalysisWindow" Height="453" Width="515" MinHeight="100" MinWidth="100" Title="Analysis" Icon="pack://application:,,,/Resources/logo_256.png"> + x:Name="AnalysisWindow" Height="453" Width="515" MinHeight="100" MinWidth="100" Title="Analysis" Icon="pack://application:,,,/Resources/icon.ico"> diff --git a/Parasitemia/Parasitemia/GUI/GUI.fs b/Parasitemia/Parasitemia/GUI/GUI.fs index 65c3cdc..8c57e20 100644 --- a/Parasitemia/Parasitemia/GUI/GUI.fs +++ b/Parasitemia/Parasitemia/GUI/GUI.fs @@ -16,7 +16,7 @@ open Emgu.CV.WPF open Config open Types -let run (defaultConfig: Config) = +let run (defaultConfig: Config) (fileToOpen: string option) = let app = new Application() let mainWindow = Views.MainWindow() let ctrl (name: string): 'a = @@ -34,7 +34,9 @@ let run (defaultConfig: Config) = let menuLoadFile: MenuItem = ctrl "menuOpen" let menuNewFile: MenuItem = ctrl "menuNew" let menuAddSourceImage: MenuItem = ctrl "menuAddSourceImage" + let menuAnalysis: MenuItem = ctrl "menuAnalysis" let menuStartAnalysis: MenuItem = ctrl "menuStartAnalysis" + let menuView: MenuItem = ctrl "menuView" let menuHightlightRBC: MenuItem = ctrl "menuHightlightRBC" let txtPatient: TextBox = ctrl "txtPatient" @@ -334,6 +336,12 @@ let run (defaultConfig: Config) = updatePreviews () updateGlobalParasitemia () + let loadFile (filepath: string) = + askSaveCurrent () + state.FilePath <- filepath + state.Load() + updateGUI () + txtPatient.LostFocus.AddHandler(fun obj args -> state.PatientID <- txtPatient.Text) menuExit.Click.AddHandler(fun obj args -> @@ -347,11 +355,7 @@ let run (defaultConfig: Config) = let dialog = OpenFileDialog(Filter = PiaZ.filter) let res = dialog.ShowDialog() if res.HasValue && res.Value - then - askSaveCurrent () - state.FilePath <- dialog.FileName - state.Load() - updateGUI ()) + then loadFile dialog.FileName) menuNewFile.Click.AddHandler(fun obj args -> askSaveCurrent () @@ -370,6 +374,9 @@ let run (defaultConfig: Config) = then updateCurrentImage ()) + + menuAnalysis.SubmenuOpened.AddHandler(fun obj args -> menuStartAnalysis.IsEnabled <- state.SourceImages.Count() > 0) + menuStartAnalysis.Click.AddHandler(fun obj args -> if Analysis.showWindow mainWindow.Root state then @@ -467,4 +474,9 @@ let run (defaultConfig: Config) = scrollViewCurrentImage.ScrollChanged.AddHandler(fun obj args -> updateViewportPreview ()) mainWindow.Root.Show() + + match fileToOpen with + | Some filepath -> loadFile filepath + | None -> () + app.Run() \ No newline at end of file diff --git a/Parasitemia/Parasitemia/GUI/MainWindow.xaml b/Parasitemia/Parasitemia/GUI/MainWindow.xaml index 314b9dd..a2151d7 100644 --- a/Parasitemia/Parasitemia/GUI/MainWindow.xaml +++ b/Parasitemia/Parasitemia/GUI/MainWindow.xaml @@ -16,10 +16,10 @@ - + - + diff --git a/Parasitemia/Parasitemia/Granulometry.fs b/Parasitemia/Parasitemia/Granulometry.fs index 7886043..a1c2fd4 100644 --- a/Parasitemia/Parasitemia/Granulometry.fs +++ b/Parasitemia/Parasitemia/Granulometry.fs @@ -52,7 +52,6 @@ let findRadiusByClosing (img: Image) (range: int * int) (scale: f float (max + r1') / scale |> roundInt - let findRadiusByAreaClosing (img: Image) (range: int * int) : int = let r1, r2 = range diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 23c9aee..096fd94 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -31,7 +31,6 @@ let saveMat (mat: Matrix<'TDepth>) (filepath: string) = mat.CopyTo(img) saveImg img filepath - type Histogram = { data: int[]; total: int; sum: int; min: float32; max: float32 } let histogramImg (img: Image) (nbSamples: int) : Histogram = @@ -148,7 +147,6 @@ let otsu (hist: Histogram) : float32 * float32 * float32 = toFloat level, toFloat mean1, toFloat mean2 - let suppressMConnections (img: Matrix) = let w = img.Width let h = img.Height @@ -163,7 +161,6 @@ let suppressMConnections (img: Matrix) = then img.[i, j] <- 0uy - let findEdges (img: Image) : Matrix * Image * Image = let w = img.Width let h = img.Height @@ -283,12 +280,10 @@ let findEdges (img: Image) : Matrix * Image edges, xGradient, yGradient - let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> = let size = 2 * int (ceil (4.0 * standardDeviation)) + 1 img.SmoothGaussian(size, size, standardDeviation, standardDeviation) - type Points = HashSet let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = @@ -363,15 +358,12 @@ let findExtremum (img: Image) (extremumType: ExtremumType) : IEnu result.Select(fun l -> Points(l)) - let findMaxima (img: Image) : IEnumerable = findExtremum img ExtremumType.Maxima - let findMinima (img: Image) : IEnumerable = findExtremum img ExtremumType.Minima - type PriorityQueue () = let size = 256 let q: Points[] = Array.init size (fun i -> Points()) @@ -470,7 +462,6 @@ type PriorityQueue () = highest <- -1 lowest <- size - type private AreaState = | Removed = 1 | Unprocessed = 2 @@ -612,7 +603,6 @@ let private areaOperation (img: Image) (area: int) (op: AreaOperatio | _ -> () () - let areaOpen (img: Image) (area: int) = areaOperation img area AreaOperation.Opening @@ -625,7 +615,6 @@ type Island (cmp: IComparer) = member val Level = 0.f with get, set member val Surface = 0 with get, set - let private areaOperationF (img: Image) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) = let w = img.Width let h = img.Height @@ -732,7 +721,6 @@ let private areaOperationF (img: Image) (areas: (int * 'a) list) | _ -> () () - let areaOpenF (img: Image) (area: int) = areaOperationF img [ area, () ] None AreaOperation.Opening @@ -800,7 +788,6 @@ let areaOpen2 (img: Image) (area: int) = for p in pointsChecked do imgData.[p.Y, p.X, 0] <- maxNeighborValue - // Zhang and Suen algorithm. // Modify 'mat' in place. let thin (mat: Matrix) = @@ -851,7 +838,6 @@ let thin (mat: Matrix) = data1 <- data2 data2 <- tmp - // Remove all 8-connected pixels with an area equal or greater than 'areaSize'. // Modify 'mat' in place. let removeArea (mat: Matrix) (areaSize: int) = @@ -920,17 +906,13 @@ let connectedComponents (img: Image) (startPoints: List) : Li List(pointChecked) - let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) = img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness); - let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) = img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias); - let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) (alpha: float) = - if alpha >= 1.0 then img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) @@ -951,11 +933,9 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img) img.ROI <- Rectangle.Empty - let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) (alpha: float) = List.iter (fun e -> drawEllipse img e color alpha) ellipses - let rngCell = System.Random() let drawCell (img: Image) (drawCellContent: bool) (c: Types.Cell) = if drawCellContent diff --git a/Parasitemia/Parasitemia/KMeans.fs b/Parasitemia/Parasitemia/KMeans.fs index afb6daf..98352fe 100644 --- a/Parasitemia/Parasitemia/KMeans.fs +++ b/Parasitemia/Parasitemia/KMeans.fs @@ -6,7 +6,6 @@ open System.Drawing open Emgu.CV open Emgu.CV.Structure - type Result = { fg: Image mean_bg: float32 diff --git a/Parasitemia/Parasitemia/KdTree.fs b/Parasitemia/Parasitemia/KdTree.fs index 7da0255..04a1152 100644 --- a/Parasitemia/Parasitemia/KdTree.fs +++ b/Parasitemia/Parasitemia/KdTree.fs @@ -96,7 +96,6 @@ type Tree<'a when 'a :> I2DCoords> = 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: float32, y: float32) = diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index 8f59514..53429b2 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -36,18 +36,17 @@ let doAnalysis (img: Image) (name: string) (config: Config) (reportPr let initialAreaOpening = int <| config.RBCRadiusByResolution.Area * config.Parameters.ratioAreaPaleCenter * 1.2f // We do an area opening a little larger to avoid to do a second one in the case the radius found is near the initial one. logTimeWithName "Area opening number one" (fun () -> ImgTools.areaOpenF filteredGreen initialAreaOpening) - report 8 + report 10 let range = let delta = config.Parameters.granulometryRange * config.RBCRadiusByResolution.Pixel int <| config.RBCRadiusByResolution.Pixel - delta, int <| config.RBCRadiusByResolution.Pixel + delta //let r1 = log "Granulometry (morpho)" (fun() -> Granulometry.findRadiusByClosing (filteredGreen.Convert()) range 1.0 |> float32) config.SetRBCRadius <| logTimeWithName "Granulometry (area)" (fun() -> Granulometry.findRadiusByAreaClosing filteredGreen range |> float32) - // log (sprintf "r1: %A, r2: %A" r1 r2) logWithName (sprintf "Found erytrocyte diameter: %A" config.RBCRadius) - report 24 + report 20 let secondAreaOpening = int <| config.RBCRadius.Area * config.Parameters.ratioAreaPaleCenter if secondAreaOpening > initialAreaOpening @@ -62,13 +61,15 @@ let doAnalysis (img: Image) (name: string) (config: Config) (reportPr removeArea edges (config.RBCRadius.Pixel ** 2.f / 50.f |> int) edges, xGradient, yGradient) - let allEllipses, ellipses = logTimeWithName "Finding ellipses" (fun () -> - let matchingEllipses = Ellipse.find edges xGradient yGradient config - matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses) + let matchingEllipses = logTimeWithName "Finding ellipses" (fun () -> Ellipse.find edges xGradient yGradient config) + + report 60 + + let prunedEllipses = logTimeWithName "Ellipses pruning" (fun () -> matchingEllipses.PrunedEllipses) report 80 - let cells = logTimeWithName "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config) + let cells = logTimeWithName "Classifier" (fun () -> Classifier.findCells prunedEllipses parasites filteredGreenWhitoutStain config) report 100 @@ -89,11 +90,11 @@ let doAnalysis (img: Image) (name: string) (config: Config) (reportPr saveImg parasites.infection (buildFileName " - parasites - infection.png") let imgAllEllipses = img.Copy() - drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.05 + drawEllipses imgAllEllipses matchingEllipses.Ellipses (Bgr(255.0, 255.0, 255.0)) 0.04 saveImg imgAllEllipses (buildFileName " - ellipses - all.png") let imgEllipses = filteredGreenWhitoutStain.Convert() - drawEllipses imgEllipses ellipses (Bgr(0.0, 240.0, 240.0)) 1.0 + drawEllipses imgEllipses prunedEllipses (Bgr(0.0, 240.0, 240.0)) 1.0 saveImg imgEllipses (buildFileName " - ellipses.png") let imgCells = img.Copy() @@ -138,10 +139,9 @@ let doMultipleAnalysis (imgs: (string * Config * Image) list) (report progressPerAnalysis.AddOrUpdate(id, progress, (fun _ _ -> progress)) |> ignore report (progressPerAnalysis.Values.Sum() / nbImgs) - let nbConcurrentTaskLimit = 4 // To reduce the memory taken. let n = Environment.ProcessorCount imgs |> PSeq.map (fun (id, config, img) -> id, doAnalysis img id config (Some (fun p -> reportProgressImg id p))) - |> PSeq.withDegreeOfParallelism (if n > nbConcurrentTaskLimit then nbConcurrentTaskLimit else n) + |> PSeq.withDegreeOfParallelism n |> PSeq.toList diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 6e57199..3599bb5 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -8,14 +8,6 @@ open System.Collections.Generic open Types open Utils - -// Do not take in account matching score below this when two ellipses are matched. -[] -let matchingScoreThreshold1 = 0.6f - -[] -let scaleOverlapTest = 0.8f - type private EllipseScoreFlaggedKd (matchingScore: float32, e: Ellipse) = let mutable matchingScore = matchingScore @@ -23,7 +15,7 @@ type private EllipseScoreFlaggedKd (matchingScore: float32, e: Ellipse) = member this.MatchingScore = matchingScore - member this.AddMatchingScore(score: float32) = + member this.AddMatchingScore (score: float32) = matchingScore <- matchingScore + score member val Processed = false with get, set @@ -33,12 +25,11 @@ type private EllipseScoreFlaggedKd (matchingScore: float32, e: Ellipse) = member this.X = this.Ellipse.Cx member this.Y = this.Ellipse.Cy - -type MatchingEllipses (radiusMin: float32) = +type MatchingEllipses (radius: float32) = let ellipses = List() // All ellipses with a score below this are removed. - let matchingScoreThreshold2 = 20.f * radiusMin + let matchingScoreThreshold = 1.f member this.Add (e: Ellipse) = ellipses.Add(EllipseScoreFlaggedKd(0.f, e)) @@ -46,7 +37,6 @@ type MatchingEllipses (radiusMin: float32) = member this.Ellipses : Ellipse list = List.ofSeq ellipses |> List.map (fun e -> e.Ellipse) - // Process all ellipses and return ellipses ordered from the best score to the worst. member this.PrunedEllipses : Ellipse list = if ellipses.Count = 0 then @@ -56,7 +46,7 @@ type MatchingEllipses (radiusMin: float32) = let tree = KdTree.Tree.BuildTree (List.ofSeq ellipses) // 2) Compute the matching score of each ellipses. - let windowSize = radiusMin + let windowSize = radius / 2.f for e in ellipses do e.Processed <- true let areaE = e.Ellipse.Area @@ -69,45 +59,33 @@ type MatchingEllipses (radiusMin: float32) = then let areaOther = other.Ellipse.Area match EEOver.EEOverlapArea e.Ellipse other.Ellipse with - | Some (commonArea, _, _) -> - let matchingScore = 2.f * commonArea / (areaE + areaOther) - if matchingScore >= matchingScoreThreshold1 + | Some (overlapArea, _, _) -> + let matchingScore = (2.f * overlapArea / (areaE + areaOther)) ** 20.f + if matchingScore <= 1.f // For approximation error. then - other.AddMatchingScore(matchingScore * e.Ellipse.Perimeter) - e.AddMatchingScore(matchingScore * other.Ellipse.Perimeter) + other.AddMatchingScore(matchingScore) + e.AddMatchingScore(matchingScore) | _ -> () - // 3) Sort ellipses by their score. - 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.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 - if nbToRemove > 0 - then - for j in i .. nbToRemove - 1 do - ellipses.[j].Removed <- true - ellipses.RemoveRange(i, nbToRemove) - - // 5) Remove ellipses whose center is into an ellipse with a better score + // 3) Remove ellipses whose center is near the center of another ellipse with a better score. for e in ellipses do - if not e.Removed + if e.MatchingScore < matchingScoreThreshold then + e.Removed <- true + else let window = { KdTree.minX = e.Ellipse.Cx - e.Ellipse.A KdTree.maxX = e.Ellipse.Cx + e.Ellipse.A KdTree.minY = e.Ellipse.Cy - e.Ellipse.A KdTree.maxY = e.Ellipse.Cy + e.Ellipse.A } for other in tree.Search window do - if not other.Removed && other.MatchingScore < e.MatchingScore + if not other.Removed && e.MatchingScore > other.MatchingScore && + distanceTwoPoints (PointD(e.Ellipse.Cx, e.Ellipse.Cy)) (PointD(other.Ellipse.Cx, other.Ellipse.Cy)) < 0.3f * e.Ellipse.B then - if e.Ellipse.Scale(scaleOverlapTest).Contains other.Ellipse.Cx other.Ellipse.Cy - then - other.Removed <- true - ellipses.RemoveAll(fun e -> e.Removed) |> ignore - - List.ofSeq ellipses |> List.map (fun e -> e.Ellipse) - + other.Removed <- true + ellipses + |> List.ofSeq + |> List.filter (fun e -> not e.Removed) + |> List.sortWith (fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore)) + |> List.map (fun e -> e.Ellipse) diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj index bb9ea52..610d802 100644 --- a/Parasitemia/Parasitemia/Parasitemia.fsproj +++ b/Parasitemia/Parasitemia/Parasitemia.fsproj @@ -112,6 +112,12 @@ + + + + + + diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs index a6fce83..671755c 100644 --- a/Parasitemia/Parasitemia/ParasitesMarker.fs +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -12,7 +12,6 @@ type Result = { infection: Image stain: Image } - // Create three binary markers : // * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas. // * 'Stain' corresponds to the stain around the parasites. @@ -44,7 +43,6 @@ let findMa (green: Image) (filteredGreen: Image) ( tmp, tmp - // Create three binary markers : // * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas. // * 'Stain' corresponds to the stain around the parasites. diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 2dcd2ff..06059f6 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -14,8 +14,8 @@ type Input = | Dir of string type RunningMode = - | CmdLine of Input * string - | Window + | CmdLine of Input * string // A file or a directory to process and the output directory. + | Window of string option // An optional path to a file to open can be given in window mode. type Arguments = RunningMode * bool @@ -32,14 +32,16 @@ let parseArgs (args: string[]) : Arguments = | Some i, Some i_output when i < args.Length - 2 && i_output < args.Length - 2 -> CmdLine ((File args.[i+1]), args.[i_output + 1]) |_ -> - Window + Window (if args.Length > 0 && not (args.[0].StartsWith("--")) then Some args.[0] else None) runningMode, Array.exists ((=) "--debug") args - [] [] let main args = + + let e = Ellipse.ellipse2 -11.4 -7.8 -0.169811 -23.75 0.8 -3.885714 -19. 1.5 + match parseArgs args with | mode, debug -> let config = Config(defaultParameters) @@ -70,15 +72,16 @@ let main args = let results = ImageAnalysis.doMultipleAnalysis images None for id, cells in results do + let config = images |> List.pick (fun (id', config', _) -> if id' = id then Some config' else None) let total, infected = Utils.countCells cells - fprintf resultFile "File: %s %d %d %.2f\n" id total infected (100. * (float infected) / (float total))) + fprintf resultFile "File: %s %d %d %.2f (diameter: %A)\n" id total infected (100. * (float infected) / (float total)) config.RBCRadius) //Utils.log (sprintf "== File: %A" file) //with //| :? IOException as ex -> Utils.log (sprintf "Unable to open the image '%A': %A" file ex) 0 - | Window -> + | Window fileToOpen -> (*let display (window : Views.MainWindow) (img : IImage) = let imgControl = window.Root.FindName("img") :?> Controls.Image imgControl.Source <- BitmapSourceConvert.ToBitmapSource(img) @@ -88,4 +91,4 @@ let main args = txtLog.Text <- txtLog.Text + mess + "\n"*) if debug then config.Debug <- DebugOn "." - GUI.Main.run config + GUI.Main.run config fileToOpen diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs index 5275b3a..5db6bb6 100644 --- a/Parasitemia/Parasitemia/Types.fs +++ b/Parasitemia/Parasitemia/Types.fs @@ -42,7 +42,6 @@ type Ellipse (cx: float32, cy: float32, a: float32, b: float32, alpha: float32) override this.ToString () = sprintf "(cx: %A, cy: %A, a: %A, b: %A, alpha: %A)" this.Cx this.Cy this.A this.B this.Alpha - type CellClass = HealthyRBC | InfectedRBC | Peculiar type Cell = { diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs index c2440b9..638f9f5 100644 --- a/Parasitemia/Parasitemia/Utils.fs +++ b/Parasitemia/Parasitemia/Utils.fs @@ -37,7 +37,7 @@ let inline linePassThroughSegment (l: Line) (p1: PointD) (p2: PointD) : bool = let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) = (p1.X - p2.X) ** 2.f + (p1.Y - p2.Y) ** 2.f -let distanceTwoPoints (p1: PointD) (p2: PointD) = +let inline distanceTwoPoints (p1: PointD) (p2: PointD) = squaredDistanceTwoPoints p1 p2 |> sqrt let countCells (cells: Cell list) : int * int = -- 2.43.0