From d0c85068bb98a7999ed994f02669befa70edd5f9 Mon Sep 17 00:00:00 2001 From: Greg Burri Date: Sun, 3 Jan 2016 21:19:44 +0100 Subject: [PATCH] Use float32 images instead of byte to improve the edge detection precision. --- Parasitemia/Parasitemia/Classifier.fs | 4 +- Parasitemia/Parasitemia/Granulometry.fs | 2 +- Parasitemia/Parasitemia/Heap.fs | 71 ++++++-- Parasitemia/Parasitemia/ImgTools.fs | 180 ++++++++++++++++++--- Parasitemia/Parasitemia/MainAnalysis.fs | 39 +---- Parasitemia/Parasitemia/ParasitesMarker.fs | 16 +- Parasitemia/Parasitemia/Program.fs | 2 +- 7 files changed, 234 insertions(+), 80 deletions(-) diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index 7f11e13..f17c0b9 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) (img: Image) (config: Config.Config) : Cell list = +let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: Image) (config: Config.Config) : Cell list = if ellipses.IsEmpty then [] @@ -100,7 +100,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: let globalStdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation(seq { for y in 0 .. h - 1 do for x in 0 .. w - 1 do - yield img.Data.[y, x, 0] |> float }) + yield float img.Data.[y, x, 0] }) for e in ellipses do let minX, minY, maxX, maxY = ellipseWindow e diff --git a/Parasitemia/Parasitemia/Granulometry.fs b/Parasitemia/Parasitemia/Granulometry.fs index db299dc..b18f53c 100644 --- a/Parasitemia/Parasitemia/Granulometry.fs +++ b/Parasitemia/Parasitemia/Granulometry.fs @@ -10,7 +10,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 = +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) let r1, r2 = range diff --git a/Parasitemia/Parasitemia/Heap.fs b/Parasitemia/Parasitemia/Heap.fs index 39c6a4b..82298e2 100644 --- a/Parasitemia/Parasitemia/Heap.fs +++ b/Parasitemia/Parasitemia/Heap.fs @@ -6,19 +6,25 @@ 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 +// TODO: measure performance with a struct. +type private Node<'k, 'v> = { mutable key: 'k; value : 'v } with override this.ToString () = sprintf "%A -> %A" this.key this.value -type Heap<'k, 'v when 'k : comparison> () = +type Heap<'k, 'v> (kComparer : IComparer<'k>) = 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 + + // Is the left child greater than the parent? + let mutable max = if l < a.Count && kComparer.Compare(a.[l].key, a.[i].key) > 0 then l else i + + // Is the right child greater than the parent and the left child? + if r < a.Count && kComparer.Compare(a.[r].key, a.[max].key) > 0 then max <- r + // If a child is greater than the parent. if max <> i then let tmp = a.[i] @@ -26,17 +32,64 @@ type Heap<'k, 'v when 'k : comparison> () = a.[max] <- tmp heapUp max + let rec checkIntegrity (i: int) : bool = + let l, r = left i, right i + + let leftIntegrity = + if l < a.Count + then + if kComparer.Compare(a.[l].key, a.[i].key) > 0 + then false + else checkIntegrity l + else + true + let rightIntegrity = + if r < a.Count + then + if kComparer.Compare(a.[r].key, a.[i].key) > 0 + then false + else checkIntegrity r + else + true + leftIntegrity && rightIntegrity + + interface IEnumerable<'k * 'v> with + member this.GetEnumerator () : IEnumerator<'k * 'v> = + (seq { for e in a -> e.key, e.value }).GetEnumerator() + + interface System.Collections.IEnumerable with + member this.GetEnumerator () : System.Collections.IEnumerator = + (this :> IEnumerable<'k * 'v>).GetEnumerator() :> System.Collections.IEnumerator + member this.Next () : 'k * 'v = - let { key = key; value = value} = a.[0] - a.RemoveAt(0) + let { key = key; value = value } = a.[0] + a.[0] <- a.[a.Count - 1] + a.RemoveAt(a.Count - 1) + heapUp 0 key, value - member this.Add (key: 'k) (value: 'v) = - a.Insert(0, { key = key; value = value}) + member this.RemoveNext () = + a.[0] <- a.[a.Count - 1] + a.RemoveAt(a.Count - 1) heapUp 0 + member this.Add (key: 'k) (value: 'v) = + a.Add({ key = key; value = value }) + + let mutable i = a.Count - 1 + while i > 0 && kComparer.Compare(a.[parent i].key, a.[i].key) < 0 do + let tmp = a.[parent i] + a.[parent i] <- a.[i] + a.[i] <- tmp + i <- parent i + member this.IsEmpty = a.Count = 0 - member this.Max : 'k = a.[0].key + member this.Count = a.Count + + member this.Max : 'k * 'v = + let max = a.[0] + max.key, max.value + member this.Clear () = a.Clear() diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 5e32ccb..cee21c7 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -72,7 +72,8 @@ let findEdges (img: Image) : Matrix * Image * yGradientData.[h - 1, c, 0] <- 0.0 use magnitudes = new Matrix(xGradient.Size) - CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles). + use angles = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles). let thresholdHigh, thresholdLow = let sensibility = 0.1 @@ -82,7 +83,6 @@ let findEdges (img: Image) : Matrix * Image * // Non-maximum suppression. use nms = new Matrix(xGradient.Size) - nms.SetValue(1.0) for i in 0 .. h - 1 do nms.Data.[i, 0] <- 0uy @@ -96,26 +96,45 @@ let findEdges (img: Image) : Matrix * Image * for j in 1 .. w - 2 do let vx = xGradient.Data.[i, j, 0] let vy = yGradient.Data.[i, j, 0] - let angle = - let a = atan2 vy vx - if a < 0.0 then 2. * Math.PI + a else a - - let mNeigbors (sign: int) : float = - if angle < Math.PI / 8. || angle >= 15.0 * Math.PI / 8. then magnitudes.Data.[i, j + sign] - elif angle < 3.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j + sign] - elif angle < 5.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j] - elif angle < 7.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j - sign] - elif angle < 9.0 * Math.PI / 8. then magnitudes.Data.[i, j - sign] - elif angle < 11.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j - sign] - elif angle < 13.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j] - else magnitudes.Data.[i - sign, j + sign] - - let m = magnitudes.Data.[i, j] - if m < mNeigbors 1 || m < mNeigbors -1 || m < thresholdLow + if vx <> 0. || vy <> 0. then - nms.Data.[i, j] <- 0uy + let angle = angles.[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] + if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 + then + nms.Data.[i, j] <- 1uy - // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection). + // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection). let edges = new Matrix(xGradient.Size) @@ -141,7 +160,6 @@ let findEdges (img: Image) : Matrix * Image * nms.Data.[ni, nj] <- 0uy toVisit.Push(Point(nj, ni)) - edges, xGradient, yGradient @@ -152,7 +170,7 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : type Points = HashSet -let drawPoints (img: Image) (points: Points) (intensity: byte) = +let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = for p in points do img.Data.[p.Y, p.X, 0] <- intensity @@ -160,7 +178,7 @@ type ExtremumType = | Maxima = 1 | Minima = 2 -let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumerable = +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 |] @@ -225,10 +243,11 @@ let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumer result.Select(fun l -> Points(l)) -let findMaxima (img: Image) : IEnumerable = +let findMaxima (img: Image) : IEnumerable = findExtremum img ExtremumType.Maxima -let findMinima (img: Image) : IEnumerable = + +let findMinima (img: Image) : IEnumerable = findExtremum img ExtremumType.Minima @@ -479,6 +498,117 @@ let areaOpen (img: Image) (area: int) = let areaClose (img: Image) (area: int) = areaOperation img area AreaOperation.Closing +[] +type Island (cmp: IComparer) = + member val Shore = Heap.Heap(cmp) with get + member val Level = 0.f with get, set + member val Surface = 0 with get, set + + +let private areaOperationF (img: Image) (area: int) (op: AreaOperation) = + let w = img.Width + let h = img.Height + let earth = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let comparer = if op = AreaOperation.Opening + then { new IComparer with member this.Compare(v1, v2) = v1.CompareTo(v2) } + else { new IComparer with member this.Compare(v1, v2) = v2.CompareTo(v1) } + + let ownership: Island[,] = Array2D.create h w null + + // Initialize islands with their shore. + let islands = List() + let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima + for e in extremum do + let island = + let p = e.First() + Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count) + islands.Add(island) + let shorePoints = Points() + for p in e do + ownership.[p.Y, p.X] <- island + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let neighbor = Point(nj, ni) + if ni >= 0 && ni < h && nj >= 0 && nj < w && ownership.[ni, nj] = null && not (shorePoints.Contains(neighbor)) + then + shorePoints.Add(neighbor) |> ignore + island.Shore.Add earth.[ni, nj, 0] neighbor + + for island in islands do + let mutable stop = island.Shore.IsEmpty + + // 'true' if 'p' is owned or adjacent to 'island'. + let ownedOrAdjacent (p: Point) : bool = + ownership.[p.Y, p.X] = island || + (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) || + (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) || + (p.X > 0 && ownership.[p.Y, p.X - 1] = island) || + (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island) + + while not stop && island.Surface < area do + let level, next = island.Shore.Max + let other = ownership.[next.Y, next.X] + if other = island // During merging, some points on the shore may be owned by the island itself -> ignored. + then + island.Shore.RemoveNext () + else + if other <> null + then // We touching another island. + if island.Surface + other.Surface >= area + then + stop <- true + else // We can merge 'other' into 'surface'. + island.Surface <- island.Surface + other.Surface + island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level + for l, p in other.Shore do + let mutable currentY = p.Y + 1 + while currentY < h && ownership.[currentY, p.X] = other do + ownership.[currentY, p.X] <- island + currentY <- currentY + 1 + island.Shore.Add l p + other.Shore.Clear() + + elif comparer.Compare(level, island.Level) > 0 + then + stop <- true + else + island.Shore.RemoveNext () + for i, j in se do + let ni = i + next.Y + let nj = j + next.X + if ni < 0 || ni >= h || nj < 0 || nj >= w + then + island.Surface <- Int32.MaxValue + stop <- true + else + let neighbor = Point(nj, ni) + if not <| ownedOrAdjacent neighbor + then + island.Shore.Add earth.[ni, nj, 0] neighbor + if not stop + then + ownership.[next.Y, next.X] <- island + island.Level <- level + island.Surface <- island.Surface + 1 + + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + let island = ownership.[i, j] + if island <> null + then + earth.[i, j, 0] <- island.Level + () + + +let areaOpenF (img: Image) (area: int) = + areaOperationF img area AreaOperation.Opening + +let areaCloseF (img: Image) (area: int) = + areaOperationF img area AreaOperation.Closing + // A simpler algorithm than 'areaOpen' but slower. let areaOpen2 (img: Image) (area: int) = let w = img.Width diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index ff82b8b..a240878 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -20,7 +20,7 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li let greenFloat = green.Convert() - let filteredGreen = gaussianFilter green config.Parameters.preFilterSigma + let filteredGreen = gaussianFilter greenFloat config.Parameters.preFilterSigma (*let maximaImg = filteredGreen.Copy() let maxima = logTime "maxima" (fun () -> ImgTools.findMaxima maximaImg) @@ -31,37 +31,19 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li logTime "areaOpen1" (fun () -> ImgTools.areaOpen greenOpen1 2000)*) let initialAreaOpen = 2000 - logTime "areaOpen 1" (fun () -> ImgTools.areaOpen filteredGreen 2000) + logTime "areaOpen 1" (fun () -> ImgTools.areaOpenF filteredGreen initialAreaOpen) config.RBCRadius <- Granulometry.findRadius filteredGreen (10, 100) 0.5 |> float let secondAreaOpen = int <| config.RBCArea / 3. if secondAreaOpen > initialAreaOpen then - logTime "areaOpen 2" (fun () -> ImgTools.areaOpen filteredGreen secondAreaOpen) + logTime "areaOpen 2" (fun () -> ImgTools.areaOpenF filteredGreen secondAreaOpen) - let filteredGreenFloat = filteredGreen.Convert() // Is it neccessary? - - let parasites, filteredGreenWhitoutInfection, filteredGreenWhitoutStain = ParasitesMarker.find filteredGreen filteredGreenFloat config + let parasites, filteredGreenWhitoutInfection, filteredGreenWhitoutStain = ParasitesMarker.find filteredGreen config //let parasites, filteredGreenWhitoutInfection, filteredGreenWhitoutStain = ParasitesMarker.findMa greenFloat filteredGreenFloat config - let filteredGreenWhitoutInfectionFloat = filteredGreenWhitoutInfection.Convert() - let filteredGreenWhitoutStainFloat = filteredGreenWhitoutStain.Convert() - - let canny = green.Canny(60., 40.) - - (*let min = ref 0.0 - let minLocation = ref <| Point() - let max = ref 0.0 - let maxLocation = ref <| Point() - magnitudes.MinMax(min, max, minLocation, maxLocation) - - use magnitudesByte = ((magnitudes / !max) * 255.0).Convert() // Otsu from OpenCV only support 'byte'. - use edges = new Matrix(xEdges.Size) - CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore - logTime "Finding edges" (fun() -> thin edges)*) - - let edges, xGradient, yGradient = ImgTools.findEdges filteredGreenWhitoutStainFloat + let edges, xGradient, yGradient = ImgTools.findEdges filteredGreenWhitoutStain logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12) @@ -79,7 +61,7 @@ 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") + //saveImg canny (buildFileName " - canny.png") saveMat (edges * 255.0) (buildFileName " - edges.png") @@ -95,8 +77,6 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li drawEllipses imgEllipses ellipses (Bgr(0.0, 240.0, 240.0)) 1.0 saveImg imgEllipses (buildFileName " - ellipses.png") - // saveImg (kmediansResults.fg * 255.0) (buildFileName " - foreground.png") - let imgCells = img.Copy() drawCells imgCells false cells saveImg imgCells (buildFileName " - cells.png") @@ -105,18 +85,15 @@ 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.Parameters.preFilterSigma + let filteredGreenMaxima = gaussianFilter greenFloat config.Parameters.preFilterSigma for m in ImgTools.findMaxima filteredGreenMaxima do - ImgTools.drawPoints filteredGreenMaxima m 255uy + ImgTools.drawPoints filteredGreenMaxima m 255.f saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png") saveImg filteredGreen (buildFileName " - filtered.png") saveImg filteredGreenWhitoutStain (buildFileName " - filtered closed stain.png") saveImg filteredGreenWhitoutInfection (buildFileName " - filtered closed infection.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) diff --git a/Parasitemia/Parasitemia/ParasitesMarker.fs b/Parasitemia/Parasitemia/ParasitesMarker.fs index b50f7e9..9e46c23 100644 --- a/Parasitemia/Parasitemia/ParasitesMarker.fs +++ b/Parasitemia/Parasitemia/ParasitesMarker.fs @@ -50,13 +50,13 @@ 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) (filteredGreenFloat: Image) (config: Config.Config) : Result * Image * Image = +let find (filteredGreen: Image) (config: Config.Config) : Result * Image * Image = let filteredGreenWithoutInfection = filteredGreen.Copy() - ImgTools.areaClose filteredGreenWithoutInfection (int config.InfectionArea) + ImgTools.areaCloseF filteredGreenWithoutInfection (int config.InfectionArea) let filteredGreenWithoutStain = filteredGreenWithoutInfection.Copy() - ImgTools.areaClose filteredGreenWithoutStain (int config.StainArea) + ImgTools.areaCloseF filteredGreenWithoutStain (int config.StainArea) // We use the filtered image to find the dark stain. @@ -71,10 +71,10 @@ let find (filteredGreen: Image) (filteredGreenFloat: Image) (closed: Image) (level: float) : Image = + let marker (img: Image) (closed: Image) (level: float) : Image = let diff = closed - (img * level) diff._ThresholdBinary(Gray(0.0), Gray(255.)) - diff + diff.Convert() let infectionMarker = marker filteredGreen filteredGreenWithoutInfection config.Parameters.infectionLevel let stainMarker = marker filteredGreenWithoutInfection filteredGreenWithoutStain config.Parameters.stainLevel @@ -86,9 +86,3 @@ let find (filteredGreen: Image) (filteredGreenFloat: Image