From 84fdf7404133803fdf0dc867a4da68a144975191 Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Mon, 21 Dec 2015 23:34:17 +0100 Subject: [PATCH] Change the parasite detection method. --- Parasitemia/Parasitemia/Classifier.fs | 40 ++-- Parasitemia/Parasitemia/Config.fs | 39 ++-- Parasitemia/Parasitemia/Ellipse.fs | 10 +- Parasitemia/Parasitemia/Heap.fs | 41 ++++ Parasitemia/Parasitemia/ImgTools.fs | 217 ++++++++++++++------ Parasitemia/Parasitemia/MainAnalysis.fs | 50 +++-- Parasitemia/Parasitemia/MatchingEllipses.fs | 6 +- Parasitemia/Parasitemia/Parasitemia.fsproj | 3 +- Parasitemia/Parasitemia/ParasitesMarker.fs | 8 +- Parasitemia/Parasitemia/ParasitesMarker2.fs | 50 +++++ Parasitemia/Parasitemia/Program.fs | 47 +++-- 11 files changed, 359 insertions(+), 152 deletions(-) create mode 100644 Parasitemia/Parasitemia/Heap.fs create mode 100644 Parasitemia/Parasitemia/ParasitesMarker2.fs diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index 3ed42f0..4c6cc56 100644 --- a/Parasitemia/Parasitemia/Classifier.fs +++ b/Parasitemia/Parasitemia/Classifier.fs @@ -21,7 +21,7 @@ type private EllipseFlaggedKd (e: Ellipse) = member this.Y = this.Cy -let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image) (config: Config.Config) : Cell list = +let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker2.Result) (img: Image) (config: Config.Config) : Cell list = if ellipses.IsEmpty then [] @@ -29,10 +29,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: let infection = parasites.infection.Copy() // To avoid to modify the parameter. // This is the minimum window size to check if other ellipses touch 'e'. - let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.maxRBCSize) * config.scale - KdTree.maxX = e.Cx + (e.A + config.maxRBCSize) * config.scale - KdTree.minY = e.Cy - (e.A + config.maxRBCSize) * config.scale - KdTree.maxY = e.Cy + (e.A + config.maxRBCSize) * config.scale } + let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMax) * config.Parameters.scale + KdTree.maxX = e.Cx + (e.A + config.RBCMax) * config.Parameters.scale + KdTree.minY = e.Cy - (e.A + config.RBCMax) * config.Parameters.scale + KdTree.maxY = e.Cy + (e.A + config.RBCMax) * config.Parameters.scale } // The minimum window to contain a given ellipse. let ellipseWindow (e: Ellipse) = @@ -40,8 +40,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: let a = int (e.A + 0.5) cx - a, cy - a, cx + a, cy + a - let w = float fg.Width - let h = float fg.Height + let w = img.Width + let w_f = float w + let h = img.Height + let h_f = float h let ellipses = ellipses |> List.map EllipseFlaggedKd @@ -60,7 +62,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: let ellipsesWithNeigbors = ellipses |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e)) // 2) Remove ellipses with a lower percentage of foreground. (taken from the lower score to the highest). - for e, neighbors in List.rev ellipsesWithNeigbors do + (*for e, neighbors in List.rev ellipsesWithNeigbors do let minX, minY, maxX, maxY = ellipseWindow e let mutable totalElement = 0 @@ -75,13 +77,13 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: then fgElement <- fgElement + 1 - if (float fgElement) / (float totalElement) < config.percentageOfFgValidCell + if (float fgElement) / (float totalElement) < config.Parameters.percentageOfFgValidCell then - e.Removed <- true + e.Removed <- true*) // 3) Remove ellipses touching the edges. for e in ellipses do - if e.isOutside w h then e.Removed <- true + if e.isOutside w_f h_f then e.Removed <- true // 4) Remove ellipses with little area. for e, neighbors in ellipsesWithNeigbors do @@ -90,16 +92,15 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: let minX, minY, maxX, maxY = ellipseWindow e let mutable area = 0 - for y in (if minY < 0 then 0 else minY) .. (if maxY >= fg.Height then fg.Height - 1 else maxY) do - for x in (if minX < 0 then 0 else minX) .. (if maxX >= fg.Width then fg.Width - 1 else maxX) do + 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 yf, xf = float y, float x - if fg.Data.[y, x, 0] > 0uy && - e.Contains xf yf && + if e.Contains xf yf && neighbors |> List.forall (fun (otherE, _, _) -> otherE.Removed || not <| otherE.Contains xf yf) then area <- area + 1 - if area < config.minimumCellArea + if area < int config.RBCMinArea then e.Removed <- true @@ -142,8 +143,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: for y in minY .. maxY do for x in minX .. maxX do let p = PointD(float x, float y) - if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) && - fg.Data.[y, x, 0] > 0uy + if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) then elements.[y-minY, x-minX] <- 1uy nbElement <- nbElement + 1 @@ -167,11 +167,11 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: darkStainPixels <- darkStainPixels + 1 let cellClass = - if float darkStainPixels > config.MaxDarkStainRatio * (float nbElement) (* || + if float darkStainPixels > config.Parameters.maxDarkStainRatio * (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 *) then Peculiar - elif infectedPixels.Count > config.infectionPixelsRequired + elif infectedPixels.Count > config.Parameters.infectionPixelsRequired then let infectionToRemove = ImgTools.connectedComponents parasites.stain infectedPixels for p in infectionToRemove do diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index 48a6334..55e142c 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -1,16 +1,16 @@ module Config +open System + type Debug = | DebugOff | DebugOn of string // Output directory. -type Config = { - debug: Debug - +type Parameters = { scale: float - minRBCSize: float - maxRBCSize: float + minRbcRadius: float + maxRbcRadius: float preFilterSigma: float @@ -21,18 +21,27 @@ type Config = { // Parasites detection. darkStainLevel: float - stainSigma: float - stainLevel: float - stainSpreadRequired: float + stainArea: float // Factor of a RBC area. 0.5 means the half of RBC area. + stainLevel: float // [0, 1] - infectionSigma: float - infectionLevel: float + infectionArea: float // Factor of a RBC area. 0.5 means the half of RBC area. + infectionLevel: float // [0, 1] infectionPixelsRequired: int - percentageOfFgValidCell: float + maxDarkStainRatio: float + + minimumCellArea: float // Factor of RBC area. +} + +type Config (param: Parameters) = + member this.Parameters = param + member val Debug = DebugOff with get, set + member val RBCSize = 30. with get, set + + member this.RBCMin = this.RBCSize + param.minRbcRadius * this.RBCSize + member this.RBCMax = this.RBCSize + param.maxRbcRadius * this.RBCSize - MaxDarkStainRatio: float + member this.RBCMinArea = param.minimumCellArea * Math.PI * this.RBCSize ** 2.0 - minimumCellArea: int - maxOffcenter: float -} \ No newline at end of file + member this.ParasiteArea = param.infectionArea * Math.PI * this.RBCSize ** 2.0 + member this.StainArea = param.stainArea * Math.PI * this.RBCSize ** 2.0 diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index 139b7b5..33333bf 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -46,8 +46,8 @@ let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) // 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 accuracy_extremum_search_1 = 7 // 3 - let accuracy_extremum_search_2 = 7 // 4 + let accuracy_extremum_search_1 = 8 // 3 + let accuracy_extremum_search_2 = 8 // 4 // p3 as the referencial. let p1x = p1x - p3x @@ -208,9 +208,9 @@ let find (edges: Matrix) (yDir: Image) (config: Config) : MatchingEllipses = - let r1, r2 = config.scale * config.minRBCSize, config.scale * config.maxRBCSize - let windowSize = roundInt (config.factorWindowSize * r2) - let factorNbPick = config.factorNbPick + let r1, r2 = config.Parameters.scale * config.RBCMin, config.Parameters.scale * config.RBCMax // FIXME: scale factor should be applied in Config!? + let windowSize = roundInt (config.Parameters.factorWindowSize * r2) + let factorNbPick = config.Parameters.factorNbPick let increment = windowSize / 4 diff --git a/Parasitemia/Parasitemia/Heap.fs b/Parasitemia/Parasitemia/Heap.fs new file mode 100644 index 0000000..c2ab69f --- /dev/null +++ b/Parasitemia/Parasitemia/Heap.fs @@ -0,0 +1,41 @@ +module Heap + +open System.Collections.Generic + +let inline private parent (i: int) : int = (i - 1) / 2 +let inline private left (i: int) : int = 2 * (i + 1) - 1 +let inline private right (i: int) : int = 2 * (i + 1) + +type private Node<'k, 'v> = { key: 'k; value : 'v } with + override this.ToString () = sprintf "%A -> %A" this.key this.value + +type Heap<'k, 'v when 'k : comparison> () = + let a = List>() + + let rec heapUp (i: int) = + let l, r = left i, right i + let mutable max = if l < a.Count && a.[l].key > a.[i].key then l else i + if r < a.Count && a.[r].key > a.[max].key + then + max <- r + if max <> i + then + let tmp = a.[i] + a.[i] <- a.[max] + a.[max] <- tmp + heapUp max + + member this.Next () : 'k * 'v = + let { key = key; value = value} = a.[0] + a.RemoveAt(0) + key, value + + member this.Add (key: 'k) (value: 'v) = + a.Insert(0, { key = key; value = value}) + heapUp 0 + + member this.IsEmpty = a.Count = 0 + member this.Max : 'k = a.[0].key + member this.Clear () = a.Clear() + + diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 24c6998..2f021c8 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -32,13 +32,18 @@ let drawPoints (img: Image) (points: Points) (intensity: byte) = for p in points do img.Data.[p.Y, p.X, 0] <- intensity -let findMaxima (img: Image) : IEnumerable = - use suppress = new Image(img.Size) + +type ExtremumType = + | Maxima = 1 + | Minima = 2 + +let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumerable = let w = img.Width let h = img.Height + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] let imgData = img.Data - let suppressData = suppress.Data + let suppress: bool[,] = Array2D.zeroCreate h w let result = List>() @@ -51,9 +56,9 @@ let findMaxima (img: Image) : IEnumerable = while betterLevelToCheck.Count > 0 do let p = betterLevelToCheck.Pop() - if suppressData.[p.Y, p.X, 0] = 0uy + if not suppress.[p.Y, p.X] then - suppressData.[p.Y, p.X, 0] <- 1uy + suppress.[p.Y, p.X] <- true sameLevelToCheck.Push(p) let current = List() @@ -63,27 +68,24 @@ let findMaxima (img: Image) : IEnumerable = let p' = sameLevelToCheck.Pop() let currentLevel = imgData.[p'.Y, p'.X, 0] current.Add(p') |> ignore - for i in -1 .. 1 do - for j in -1 .. 1 do - if i <> 0 || j <> 0 + for i, j in se do + let ni = i + p'.Y + let nj = j + p'.X + if ni >= 0 && ni < h && nj >= 0 && nj < w + then + let level = imgData.[ni, nj, 0] + let notSuppressed = not suppress.[ni, nj] + + if level = currentLevel && notSuppressed then - let ni = i + p'.Y - let nj = j + p'.X - if ni >= 0 && ni < h && nj >= 0 && nj < w + suppress.[ni, nj] <- true + sameLevelToCheck.Push(Point(nj, ni)) + elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel + then + betterExists <- true + if notSuppressed then - let level = imgData.[ni, nj, 0] - let notSuppressed = suppressData.[ni, nj, 0] = 0uy - - if level = currentLevel && notSuppressed - then - suppressData.[ni, nj, 0] <- 1uy - sameLevelToCheck.Push(Point(nj, ni)) - elif level > currentLevel - then - betterExists <- true - if notSuppressed - then - betterLevelToCheck.Push(Point(nj, ni)) + betterLevelToCheck.Push(Point(nj, ni)) if not betterExists then @@ -94,18 +96,26 @@ let findMaxima (img: Image) : IEnumerable = for j in 0 .. w - 1 do let maxima = flood (Point(j, i)) if maxima.Count > 0 - then result.AddRange(maxima) + then + result.AddRange(maxima) 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()) // TODO: Check performance with an HasSet + let q: Points[] = Array.init size (fun i -> Points()) let mutable highest = -1 // Value of the first elements of 'q'. let mutable lowest = size - member this.Next () : byte * Point = + member this.NextMax () : byte * Point = if this.IsEmpty then invalidOp "Queue is empty" @@ -127,9 +137,34 @@ type PriorityQueue () = value, next + member this.NextMin () : byte * Point = + if this.IsEmpty + then + invalidOp "Queue is empty" + else + let l = q.[lowest + 1] + let next = l.First() + l.Remove(next) |> ignore + let value = byte (lowest + 1) + + if l.Count = 0 + then + lowest <- lowest + 1 + while lowest < highest && q.[lowest + 1].Count = 0 do + lowest <- lowest + 1 + if highest = lowest + then + highest <- -1 + lowest <- size + + value, next + member this.Max = highest |> byte + member this.Min = + lowest + 1 |> byte + member this.Add (value: byte) (p: Point) = let vi = int value @@ -173,23 +208,28 @@ type PriorityQueue () = lowest <- size -type AreaState = +type private AreaState = | Removed = 1 | Unprocessed = 2 | Validated = 3 +type private AreaOperation = + | Opening = 1 + | Closing = 2 + [] -type Area (elements: Points) = +type private Area (elements: Points) = member this.Elements = elements member val Intensity = None with get, set member val State = AreaState.Unprocessed with get, set -let areaOpen (img: Image) (area: int) = +let private areaOperation (img: Image) (area: int) (op: AreaOperation) = let w = img.Width let h = img.Height let imgData = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] - let areas = List(findMaxima img |> Seq.map Area) + let areas = List((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area) let pixels: Area[,] = Array2D.create h w null for m in areas do @@ -200,16 +240,13 @@ let areaOpen (img: Image) (area: int) = let addEdgeToQueue (elements: Points) = for p in elements do - for i in -1 .. 1 do - for j in -1 .. 1 do - if i <> 0 || j <> 0 - then - let ni = i + p.Y - let nj = j + p.X - let p' = Point(nj, ni) - if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p')) - then - queue.Add (imgData.[ni, nj, 0]) p' + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let p' = Point(nj, ni) + if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p')) + then + queue.Add (imgData.[ni, nj, 0]) p' // Reverse order is quicker. for i in areas.Count - 1 .. -1 .. 0 do @@ -219,12 +256,12 @@ let areaOpen (img: Image) (area: int) = queue.Clear() addEdgeToQueue m.Elements - let mutable intensity = queue.Max + let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min let nextElements = Points() let mutable stop = false while not stop do - let intensity', p = queue.Next () + let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin () let mutable merged = false if intensity' = intensity // The intensity doesn't change. @@ -237,7 +274,7 @@ let areaOpen (img: Image) (area: int) = else nextElements.Add(p) |> ignore - elif intensity' < intensity + elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity then m.Elements.UnionWith(nextElements) for e in nextElements do @@ -253,7 +290,7 @@ let areaOpen (img: Image) (area: int) = nextElements.Clear() nextElements.Add(p) |> ignore - else // i' > i + else let m' = pixels.[p.Y, p.X] if m' <> null then @@ -265,7 +302,7 @@ let areaOpen (img: Image) (area: int) = queue.Remove imgData.[e.Y, e.X, 0] e addEdgeToQueue m'.Elements m.Elements.UnionWith(m'.Elements) - let intensityMax = queue.Max + let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min if intensityMax <> intensity then intensity <- intensityMax @@ -280,21 +317,18 @@ let areaOpen (img: Image) (area: int) = if not stop && not merged then - for i in -1 .. 1 do - for j in -1 .. 1 do - if i <> 0 || j <> 0 - then - let ni = i + p.Y - let nj = j + p.X - let p' = Point(nj, ni) - if ni < 0 || ni >= h || nj < 0 || nj >= w - then - m.State <- AreaState.Validated - m.Intensity <- Some (intensity) - stop <- true - elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) - then - queue.Add (imgData.[ni, nj, 0]) p' + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let p' = Point(nj, ni) + if ni < 0 || ni >= h || nj < 0 || nj >= w + then + m.State <- AreaState.Validated + m.Intensity <- Some (intensity) + stop <- true + elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) + then + queue.Add (imgData.[ni, nj, 0]) p' if queue.IsEmpty then @@ -316,6 +350,67 @@ let areaOpen (img: Image) (area: int) = () +let areaOpen (img: Image) (area: int) = + areaOperation img area AreaOperation.Opening + +let areaClose (img: Image) (area: int) = + areaOperation img area AreaOperation.Closing + +let areaOpen2 (img: Image) (area: int) = + let w = img.Width + let h = img.Height + let imgData = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let histogram = Array.zeroCreate 256 + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + let v = imgData.[i, j, 0] |> int + histogram.[v] <- histogram.[v] + 1 + + let flooded : bool[,] = Array2D.zeroCreate h w + + let pointsChecked = HashSet() + let pointsToCheck = Stack() + + for level in 255 .. -1 .. 0 do + let mutable n = histogram.[level] + if n > 0 + then + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + if not flooded.[i, j] && imgData.[i, j, 0] = byte level + then + let mutable maxNeighborValue = 0uy + pointsChecked.Clear() + pointsToCheck.Clear() + pointsToCheck.Push(Point(j, i)) + + while pointsToCheck.Count > 0 do + let next = pointsToCheck.Pop() + pointsChecked.Add(next) |> ignore + flooded.[next.Y, next.X] <- true + + for nx, ny in se do + let p = Point(next.X + nx, next.Y + ny) + if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h + then + let v = imgData.[p.Y, p.X, 0] + if v = byte level + then + if not (pointsChecked.Contains(p)) + then + pointsToCheck.Push(p) + elif v > maxNeighborValue + then + maxNeighborValue <- v + + if int maxNeighborValue < level && pointsChecked.Count <= area + then + for p in pointsChecked do + imgData.[p.Y, p.X, 0] <- maxNeighborValue + + // Zhang and Suen algorithm. // Modify 'mat' in place. let thin (mat: Matrix) = diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index 39383fc..d959624 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -14,31 +14,40 @@ open Types let doAnalysis (img: Image) (name: string) (config: Config) : Cell list = - use scaledImg = if config.scale = 1.0 then img else img.Resize(config.scale, CvEnum.Inter.Area) + use scaledImg = if config.Parameters.scale = 1.0 then img else img.Resize(config.Parameters.scale, CvEnum.Inter.Area) - use blue = scaledImg.Item(0) use green = scaledImg.Item(1) - use red = scaledImg.Item(2) let greenFloat = green.Convert() - let filteredGreen = gaussianFilter green config.preFilterSigma - logTime "areaOpen" (fun () -> ImgTools.areaOpen filteredGreen 2000) + let filteredGreen = gaussianFilter green config.Parameters.preFilterSigma + + (*let maximaImg = filteredGreen.Copy() + let maxima = logTime "maxima" (fun () -> ImgTools.findMaxima maximaImg) + for m in maxima do + ImgTools.drawPoints maximaImg m 255uy - let RBCSize = Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float - let RBCSizeDelta = RBCSize / 3.0 + let greenOpen1 = filteredGreen.Copy() + logTime "areaOpen1" (fun () -> ImgTools.areaOpen greenOpen1 2000)*) + + logTime "areaOpen" (fun () -> ImgTools.areaOpen filteredGreen 2000) - let config = { config with minRBCSize = RBCSize - RBCSizeDelta ; maxRBCSize = RBCSize + RBCSizeDelta } + config.RBCSize <- Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float let filteredGreenFloat = filteredGreen.Convert() // Is it neccessary? + let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreenFloat 1.0) + + let parasites, filteredGreenWhitoutParasites, filteredGreenWhitoutStain = ParasitesMarker2.find filteredGreen filteredGreenFloat kmediansResults config + let filteredGreenWhitoutParasitesFloat = filteredGreenWhitoutParasites.Convert() + use sobelKernel = new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] [ 2.0f; 0.0f; -2.0f ] [ 1.0f; 0.0f; -1.0f ]], Point(0, 0)) - use xEdges = filteredGreenFloat.Convolution(sobelKernel).Convert() - use yEdges = filteredGreenFloat.Convolution(sobelKernel.Transpose()).Convert() + use xEdges = filteredGreenWhitoutParasitesFloat.Convolution(sobelKernel).Convert() + use yEdges = filteredGreenWhitoutParasitesFloat.Convolution(sobelKernel.Transpose()).Convert() let xEdgesData = xEdges.Data let yEdgesData = yEdges.Data @@ -70,18 +79,14 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li logTime "Finding edges" (fun() -> thin edges) logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12) - let kmediansResults = logTime "Finding foreground (k-medians)" (fun () -> KMedians.kmedians filteredGreenFloat 1.0) - - let parasites = ParasitesMarker.find greenFloat filteredGreenFloat kmediansResults config - let allEllipses, ellipses = logTime "Finding ellipses" (fun () -> let matchingEllipses = Ellipse.find edges xEdges yEdges config matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses ) - let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites kmediansResults.fg config) + let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutParasites config) // Output pictures if debug flag is set. - match config.debug with + match config.Debug with | DebugOn output -> let dirPath = System.IO.Path.Combine(output, name) System.IO.Directory.CreateDirectory dirPath |> ignore @@ -112,14 +117,23 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li drawCells imgCells' true cells saveImg imgCells' (buildFileName " - cells - full.png") - let filteredGreenMaxima = gaussianFilter green config.preFilterSigma + let filteredGreenMaxima = gaussianFilter green config.Parameters.preFilterSigma for m in ImgTools.findMaxima filteredGreenMaxima do ImgTools.drawPoints filteredGreenMaxima m 255uy saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png") saveImg filteredGreen (buildFileName " - filtered.png") - saveImg blue (buildFileName " - blue.png") + saveImg filteredGreenWhitoutParasites (buildFileName " - filtered closed.png") + + (*saveImg parasitesMarker (buildFileName " - parasites (area closing).png") + saveImg stainMarker (buildFileName " - stain (area closing).png")*) + saveImg green (buildFileName " - green.png") + + use blue = scaledImg.Item(0) + saveImg blue (buildFileName " - blue.png") + + use red = scaledImg.Item(2) saveImg red (buildFileName " - red.png") | _ -> () diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 99e3cb5..2196122 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -13,19 +13,17 @@ open Utils let matchingScoreThreshold1 = 0.6 // All ellipsee with a score below this is removed. -let matchingScoreThreshold2 = 1. / 100. +let matchingScoreThreshold2 = 2. type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) = let mutable matchingScore = matchingScore - let perimeter = e.Perimeter member this.Ellipse = e member this.MatchingScore = matchingScore - // The score is proportional to the perimeter because large ellipse will receive more votes. member this.AddMatchingScore(score: float) = - matchingScore <- matchingScore + score / perimeter + matchingScore <- matchingScore + score member val Processed = false with get, set member val Removed = false with get, set diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj index 3e6609a..fe81290 100644 --- a/Parasitemia/Parasitemia/Parasitemia.fsproj +++ b/Parasitemia/Parasitemia/Parasitemia.fsproj @@ -28,7 +28,7 @@ x64 bin\Debug\Parasitemia.XML false - --folder "../../../Images/areaopen" --output "../../../Images/output" --debug + --folder "../../../Images/debug" --output "../../../Images/output" --debug pdbonly @@ -72,6 +72,7 @@ + diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs index 9adf185..d1a51fa 100644 --- a/Parasitemia/Parasitemia/ParasitesMarker.fs +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -14,11 +14,11 @@ type Result = { // * '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 (green: Image) (filteredGreen: Image) (kmediansResult: KMedians.Result) (config: Config.Config) : Result = +(*let find (green: Image) (filteredGreen: Image) (kmediansResult: KMedians.Result) (config: Config.Config) : Result = // We use the filtered image to find the dark stain. let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResult - let darkStain = d_fg.Cmp(median_bg * config.darkStainLevel, CvEnum.CmpType.GreaterThan) + let darkStain = d_fg.Cmp(median_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan) darkStain._And(filteredGreen.Cmp(median_fg, CvEnum.CmpType.LessThan)) darkStain._And(fg) @@ -34,8 +34,8 @@ let find (green: Image) (filteredGreen: Image) (km smears { darkStain = darkStain; - stain = findSmears config.stainSigma config.stainLevel - infection = findSmears config.infectionSigma config.infectionLevel } + stain = findSmears config.Parameters.stainSigma config.Parameters.stainLevel + infection = findSmears config.Parameters.infectionSigma config.Parameters.infectionLevel }*) diff --git a/Parasitemia/Parasitemia/ParasitesMarker2.fs b/Parasitemia/Parasitemia/ParasitesMarker2.fs new file mode 100644 index 0000000..03af2b6 --- /dev/null +++ b/Parasitemia/Parasitemia/ParasitesMarker2.fs @@ -0,0 +1,50 @@ +module ParasitesMarker2 + +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +type Result = { + darkStain: Image + stain: Image + infection: 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. +// * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes. +let find (filteredGreen: Image) (filteredGreenFloat: Image) (kmediansResult: KMedians.Result) (config: Config.Config) : Result * Image * Image = + + // We use the filtered image to find the dark stain. + let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResult + let darkStain = d_fg.Cmp(median_bg * config.Parameters.darkStainLevel, CvEnum.CmpType.GreaterThan) + darkStain._And(filteredGreenFloat.Cmp(median_fg, CvEnum.CmpType.LessThan)) + + let marker (area: int) (threshold: float) : Image * Image = + let closed = filteredGreen.Copy() + ImgTools.areaClose closed area + let diff = closed - filteredGreen + + let min = ref [| 0. |] + let minLocation = ref <| [| Point() |] + let max = ref [| 0. |] + let maxLocation = ref <| [| Point() |] + diff.MinMax(min, max, minLocation, maxLocation) + + diff._ThresholdBinary((!max).[0] * threshold |> Gray, Gray(255.)) + diff, closed + + let parasitesMarker, filteredGreenWithoutParasites = marker (int config.ParasiteArea) config.Parameters.infectionLevel + let stainMarker, filteredGreenWithoutStain = marker (int config.StainArea) config.Parameters.stainLevel + + { darkStain = darkStain + stain = parasitesMarker + infection = stainMarker }, + filteredGreenWithoutParasites, + filteredGreenWithoutStain + + + + + diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 86de435..103596a 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -57,41 +57,36 @@ let parseArgs (args: string[]) : Arguments = let main args = match parseArgs args with | mode, debug -> - let scale = 1. - let config = { - debug = if debug then DebugOn "." else DebugOff + let config = + Config( + { scale = 1. - scale = scale + minRbcRadius = -0.3 + maxRbcRadius = 0.3 - // RBC size range in px at scale = 1.0. - minRBCSize = 20. - maxRBCSize = 42. + preFilterSigma = 1.5 - preFilterSigma = 1.5 + factorNbPick = 1.0 + factorWindowSize = 2.0 - factorNbPick = 1.0 // 1.0 - factorWindowSize = 2.0 + darkStainLevel = 0.4 // Lower -> more sensitive. - darkStainLevel = 0.4 // Lower -> more sensitive. + stainArea = 0.02 + stainLevel = 0.2 - stainSigma = 10. - stainLevel = 0.9 - stainSpreadRequired = 3.0 + infectionArea = 0.06 + infectionLevel = 0.2 + infectionPixelsRequired = 1 - infectionSigma = 2.2 - infectionLevel = 0.87 - infectionPixelsRequired = 1 + maxDarkStainRatio = 0.1 - percentageOfFgValidCell = 0.4 - - MaxDarkStainRatio = 0.1 - - minimumCellArea = 1000. * scale ** 2. |> int - maxOffcenter = 0.5 } + minimumCellArea = 0.3 }) match mode with | CmdLine (input, output) -> - let config = if debug then { config with debug = DebugOn output } else config + if debug + then + config.Debug <- DebugOn output Directory.CreateDirectory output |> ignore @@ -121,6 +116,10 @@ let main args = let app = new Application() let mainWindow = Views.MainWindow() + if debug + then + config.Debug <- DebugOn "." + Utils.log <- (fun m -> log mainWindow m) //display mainWindow img -- 2.43.0