X-Git-Url: http://git.euphorik.ch/?a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=f567e682decb86c97b81aefeb8c7a44eeab1f347;hb=bef2e9f0bf1bba21d4c988fdf654c2dc303ec34a;hp=24c6998ed97cd505204d4b253849d3f35dc56456;hpb=5ac2dedf8ead5275ac216e0b41829ab39c843800;p=master-thesis.git diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 24c6998..f567e68 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -21,6 +21,124 @@ let normalizeAndConvert (img: Image) : Image = ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() +let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = + img.Save(filepath) + + +let saveMat (mat: Matrix<'TDepth>) (filepath: string) = + use img = new Image(mat.Size) + mat.CopyTo(img) + saveImg img filepath + + +let suppressMConnections (img: Matrix) = + let w = img.Width + let h = img.Height + for i in 1 .. h - 2 do + for j in 1 .. w - 2 do + if img.[i, j] > 0uy && img.Data.[i + 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i - 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i - 1, j - 1] = 0uy) then + img.[i, j] <- 0uy + for i in 1 .. h - 2 do + for j in 1 .. w - 2 do + if img.[i, j] > 0uy && img.Data.[i - 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i + 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i + 1, j - 1] = 0uy) then + img.[i, j] <- 0uy + +let findEdges (img: Image) : Matrix * Image * Image = + let w = img.Width + let h = img.Height + + use sobelKernel = + new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] + [ 2.0f; 0.0f; -2.0f ] + [ 1.0f; 0.0f; -1.0f ]], Point(1, 1)) + + let xGradient = img.Convolution(sobelKernel).Convert() + let yGradient = img.Convolution(sobelKernel.Transpose()).Convert() + + let xGradientData = xGradient.Data + let yGradientData = yGradient.Data + for r in 0 .. h - 1 do + xGradientData.[r, 0, 0] <- 0.0 + xGradientData.[r, w - 1, 0] <- 0.0 + yGradientData.[r, 0, 0] <- 0.0 + yGradientData.[r, w - 1, 0] <- 0.0 + + for c in 0 .. w - 1 do + xGradientData.[0, c, 0] <- 0.0 + xGradientData.[h - 1, c, 0] <- 0.0 + yGradientData.[0, c, 0] <- 0.0 + yGradientData.[h - 1, c, 0] <- 0.0 + + use magnitudes = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles). + + let thresholdHigh, thresholdLow = + let sensibility = 0.1 + use magnitudesByte = magnitudes.Convert() + let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) + threshold + (sensibility * threshold), threshold - (sensibility * threshold) + + // 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 + nms.Data.[i, w - 1] <- 0uy + + for j in 0 .. w - 1 do + nms.Data.[0, j] <- 0uy + nms.Data.[h - 1, j] <- 0uy + + for i in 1 .. h - 2 do + for j in 1 .. w - 2 do + let vx = xGradient.Data.[i, j, 0] + let vy = yGradient.Data.[i, j, 0] + 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 then + nms.Data.[i, j] <- 0uy + + // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection). + + let edges = new Matrix(xGradient.Size) + + // Histeresis thresholding. + let toVisit = Stack() + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh then + nms.Data.[i, j] <- 0uy + toVisit.Push(Point(j, i)) + while toVisit.Count > 0 do + let p = toVisit.Pop() + edges.Data.[p.Y, p.X] <- 1uy + for i' in -1 .. 1 do + for j' in -1 .. 1 do + if i' <> 0 || j' <> 0 then + let ni = p.Y + i' + let nj = p.X + j' + if ni >= 0 && ni < h && nj >= 0 && nj < w && nms.Data.[ni, nj] = 1uy then + nms.Data.[ni, nj] <- 0uy + toVisit.Push(Point(nj, ni)) + + + 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) @@ -32,13 +150,17 @@ 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 +173,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 +185,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 + suppress.[ni, nj] <- true + sameLevelToCheck.Push(Point(nj, ni)) + elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel then - let ni = i + p'.Y - let nj = j + p'.X - if ni >= 0 && ni < h && nj >= 0 && nj < w + 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 +213,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 +254,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 +325,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 +357,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 +373,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 +391,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 +407,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 +419,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 +434,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 +467,68 @@ 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 + +// A simpler algorithm than 'areaOpen' but slower. +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) = @@ -442,21 +655,14 @@ let connectedComponents (img: Image) (startPoints: List) : Li List(pointChecked) -let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = - img.Save(filepath) - - -let saveMat (mat: Matrix<'TDepth>) (filepath: string) = - use img = new Image(mat.Size) - mat.CopyTo(img) - saveImg img filepath - 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