* Add an exact method to compute an ellipse from three points and two tangents.
authorGreg Burri <greg.burri@gmail.com>
Sun, 17 Jan 2016 17:00:46 +0000 (18:00 +0100)
committerGreg Burri <greg.burri@gmail.com>
Sun, 17 Jan 2016 17:00:46 +0000 (18:00 +0100)
* Simplification of the matching ellipses process.

18 files changed:
Parasitemia/Parasitemia/Config.fs
Parasitemia/Parasitemia/Const.fs
Parasitemia/Parasitemia/EEOver.fs
Parasitemia/Parasitemia/Ellipse.fs
Parasitemia/Parasitemia/GUI/AnalysisWindow.xaml
Parasitemia/Parasitemia/GUI/GUI.fs
Parasitemia/Parasitemia/GUI/MainWindow.xaml
Parasitemia/Parasitemia/Granulometry.fs
Parasitemia/Parasitemia/ImgTools.fs
Parasitemia/Parasitemia/KMeans.fs
Parasitemia/Parasitemia/KdTree.fs
Parasitemia/Parasitemia/MainAnalysis.fs
Parasitemia/Parasitemia/MatchingEllipses.fs
Parasitemia/Parasitemia/Parasitemia.fsproj
Parasitemia/Parasitemia/ParasitesMarker.fs
Parasitemia/Parasitemia/Program.fs
Parasitemia/Parasitemia/Types.fs
Parasitemia/Parasitemia/Utils.fs

index cd8f659..e7deb61 100644 (file)
@@ -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<ppi>) : float32 =
         let rbcRadiusInch: float<inch> = (μmToInch rbcDiameter) / 2.
index 0d2e52e..304d3a9 100644 (file)
@@ -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
index 5a3c9e4..2bd13b2 100644 (file)
@@ -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
index be35940..ac3c041 100644 (file)
@@ -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<byte>)
 
     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<byte>)
                         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
index 7684ad9..599d581 100644 (file)
@@ -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">
    <Grid>
       <Grid.RowDefinitions>
          <RowDefinition Height="50*"/>
index 65c3cdc..8c57e20 100644 (file)
@@ -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
index 314b9dd..a2151d7 100644 (file)
          <MenuItem Header="_Images">
             <MenuItem x:Name="menuAddSourceImage" Header="_Add a source image" />
          </MenuItem>
-         <MenuItem Header="_Analysis">
+         <MenuItem x:Name="menuAnalysis" Header="_Analysis">
             <MenuItem x:Name="menuStartAnalysis" Header="_Show analysis window" />
          </MenuItem>
-         <MenuItem Header="_View">
+         <MenuItem x:Name="menuView" Header="_View">
             <MenuItem x:Name="menuHightlightRBC" Header="_Highlight healthy erytrocytes" IsCheckable="True" />
          </MenuItem>
       </Menu>
index 7886043..a1c2fd4 100644 (file)
@@ -52,7 +52,6 @@ let findRadiusByClosing (img: Image<Gray, 'TDepth>) (range: int * int) (scale: f
 
     float (max + r1') / scale |> roundInt
 
-
 let findRadiusByAreaClosing (img: Image<Gray, float32>) (range: int * int) : int =
     let r1, r2 = range
 
index 23c9aee..096fd94 100644 (file)
@@ -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<Gray, float32>) (nbSamples: int) : Histogram =
@@ -148,7 +147,6 @@ let otsu (hist: Histogram) : float32 * float32 * float32 =
 
     toFloat level, toFloat mean1, toFloat mean2
 
-
 let suppressMConnections (img: Matrix<byte>) =
     let w = img.Width
     let h = img.Height
@@ -163,7 +161,6 @@ let suppressMConnections (img: Matrix<byte>) =
             then
                 img.[i, j] <- 0uy
 
-
 let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32> * Image<Gray, float32> =
     let w = img.Width
     let h = img.Height
@@ -283,12 +280,10 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
 
     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<Point>
 
 let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
@@ -363,15 +358,12 @@ let findExtremum (img: Image<Gray, 'TDepth>) (extremumType: ExtremumType) : IEnu
 
     result.Select(fun l -> Points(l))
 
-
 let findMaxima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
     findExtremum img ExtremumType.Maxima
 
-
 let findMinima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
     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<Gray, byte>) (area: int) (op: AreaOperatio
             | _ -> ()
     ()
 
-
 let areaOpen (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Opening
 
@@ -625,7 +615,6 @@ type Island (cmp: IComparer<float32>) =
     member val Level = 0.f with get, set
     member val Surface = 0 with get, set
 
-
 let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) =
     let w = img.Width
     let h = img.Height
@@ -732,7 +721,6 @@ let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list)
         | _ -> ()
     ()
 
-
 let areaOpenF (img: Image<Gray, float32>) (area: int) =
     areaOperationF img [ area, () ] None AreaOperation.Opening
 
@@ -800,7 +788,6 @@ let areaOpen2 (img: Image<Gray, byte>) (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<byte>) =
@@ -851,7 +838,6 @@ let thin (mat: Matrix<byte>) =
         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<byte>) (areaSize: int) =
@@ -920,17 +906,13 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
 
     List<Point>(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<Bgr, byte>) (drawCellContent: bool) (c: Types.Cell) =
     if drawCellContent
index afb6daf..98352fe 100644 (file)
@@ -6,7 +6,6 @@ open System.Drawing
 open Emgu.CV
 open Emgu.CV.Structure
 
-
 type Result = {
     fg: Image<Gray, byte>
     mean_bg: float32
index 7da0255..04a1152 100644 (file)
@@ -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) =
index 8f59514..53429b2 100644 (file)
@@ -36,18 +36,17 @@ let doAnalysis (img: Image<Bgr, byte>) (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<Gray, byte>()) 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<Bgr, byte>) (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<Bgr, byte>) (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<Bgr, byte>()
-        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<Bgr, byte>) 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
index 6e57199..3599bb5 100644 (file)
@@ -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.
-[<Literal>]
-let matchingScoreThreshold1 = 0.6f
-
-[<Literal>]
-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<EllipseScoreFlaggedKd>()
 
     // 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<EllipseScoreFlaggedKd> 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)
 
index bb9ea52..610d802 100644 (file)
     <None Include="App.config" />
     <Content Include="packages.config" />
     <Resource Include="Resources\icon.ico" />
+    <None Include="Scripts\load-references-debug.fsx" />
+    <None Include="Scripts\load-project-debug.fsx" />
+    <None Include="Scripts\load-references-release.fsx" />
+    <None Include="Scripts\load-project-release.fsx" />
+    <None Include="Scripts\load-references-debuggui.fsx" />
+    <None Include="Scripts\load-project-debuggui.fsx" />
   </ItemGroup>
   <ItemGroup>
     <Reference Include="Castle.Core">
index a6fce83..671755c 100644 (file)
@@ -12,7 +12,6 @@ type Result = {
     infection: Image<Gray, byte>
     stain: Image<Gray, byte> }
 
-
 // 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<Gray, float32>) (filteredGreen: Image<Gray, float32>) (
     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.
index 2dcd2ff..06059f6 100644 (file)
@@ -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
 
-
 [<EntryPoint>]
 [<STAThread()>]
 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
index 5275b3a..5db6bb6 100644 (file)
@@ -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 = {
index c2440b9..638f9f5 100644 (file)
@@ -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 =