X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=f64b1c3318d687d6dc314fe98f4dbdd125446d6f;hp=e9cbbe9ee0f59f7a4f54c557ff345ca2f20be6df;hb=044b0ae69df3ac565432545b2fa934589016f9bd;hpb=b070295cf67b2025164a34b6594e84f0d771cdc9 diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index e9cbbe9..f64b1c3 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -8,17 +8,163 @@ open System.Linq open Emgu.CV open Emgu.CV.Structure -open Utils open Heap +open Const +open Utils // Normalize image values between 0uy and 255uy. -let normalizeAndConvert (img: Image) : Image = +let normalizeAndConvert (img: Image) : Image = let min = ref [| 0.0 |] let minLocation = ref <| [| Point() |] let max = ref [| 0.0 |] let maxLocation = ref <| [| Point() |] img.MinMax(min, max, minLocation, maxLocation) - ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() + ((img.Convert() - (!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) + let yGradient = img.Convolution(sobelKernel.Transpose()) + + let xGradientData = xGradient.Data + let yGradientData = yGradient.Data + for r in 0 .. h - 1 do + xGradientData.[r, 0, 0] <- 0.f + xGradientData.[r, w - 1, 0] <- 0.f + yGradientData.[r, 0, 0] <- 0.f + yGradientData.[r, w - 1, 0] <- 0.f + + for c in 0 .. w - 1 do + xGradientData.[0, c, 0] <- 0.f + xGradientData.[h - 1, c, 0] <- 0.f + yGradientData.[0, c, 0] <- 0.f + yGradientData.[h - 1, c, 0] <- 0.f + + use magnitudes = new Matrix(xGradient.Size) + use angles = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles). + + let thresholdHigh, thresholdLow = + let sensibilityHigh = 0.1f + let sensibilityLow = 0.1f + use magnitudesByte = magnitudes.Convert() + let threshold = float32 <| CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) + threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold) + + // Non-maximum suppression. + use nms = new Matrix(xGradient.Size) + + let nmsData = nms.Data + let anglesData = angles.Data + let magnitudesData = magnitudes.Data + let xGradientData = xGradient.Data + let yGradientData = yGradient.Data + + let PI = float32 Math.PI + + for i in 0 .. h - 1 do + nmsData.[i, 0] <- 0uy + nmsData.[i, w - 1] <- 0uy + + for j in 0 .. w - 1 do + nmsData.[0, j] <- 0uy + nmsData.[h - 1, j] <- 0uy + + for i in 1 .. h - 2 do + for j in 1 .. w - 2 do + let vx = xGradientData.[i, j, 0] + let vy = yGradientData.[i, j, 0] + if vx <> 0.f || vy <> 0.f + then + let angle = anglesData.[i, j] + + let vx', vy' = abs vx, abs vy + let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy' + let ratio1 = 1.f - ratio2 + + let mNeigbors (sign: int) : float32 = + if angle < PI / 4.f + then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign] + elif angle < PI / 2.f + then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j] + elif angle < 3.f * PI / 4.f + then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign] + elif angle < PI + then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign] + elif angle < 5.f * PI / 4.f + then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign] + elif angle < 3.f * PI / 2.f + then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j] + elif angle < 7.f * PI / 4.f + then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign] + else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign] + + let m = magnitudesData.[i, j] + if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 + then + nmsData.[i, j] <- 1uy + + // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection). + + let edges = new Matrix(xGradient.Size) + let edgesData = edges.Data + + // Hysteresis thresholding. + let toVisit = Stack() + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh + then + nmsData.[i, j] <- 0uy + toVisit.Push(Point(j, i)) + while toVisit.Count > 0 do + let p = toVisit.Pop() + edgesData.[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 && nmsData.[ni, nj] = 1uy + then + nmsData.[ni, nj] <- 0uy + toVisit.Push(Point(nj, ni)) + + edges, xGradient, yGradient let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> = @@ -26,13 +172,23 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : img.SmoothGaussian(size, size, standardDeviation, standardDeviation) -let findMaxima (img: Image) : IEnumerable> = - use suppress = new Image(img.Size) +type Points = HashSet + +let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = + for p in points do + img.Data.[p.Y, p.X, 0] <- intensity + +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>() @@ -45,9 +201,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() @@ -57,27 +213,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 @@ -88,198 +241,434 @@ 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 - result.Select(fun l -> HashSet(l)) + +let findMinima (img: Image) : IEnumerable = + findExtremum img ExtremumType.Minima type PriorityQueue () = - let q = List>() // TODO: Check performance with an HasSet + let size = 256 + 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" else - let l = q.[0] + let l = q.[highest] let next = l.First() l.Remove(next) |> ignore let value = byte highest + if l.Count = 0 then - q.RemoveAt(0) highest <- highest - 1 - while q.Count > 0 && q.[0] = null do - q.RemoveAt(0) + while highest > lowest && q.[highest].Count = 0 do highest <- highest - 1 + if highest = lowest + then + highest <- -1 + lowest <- size + + 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.UnionWith (other: PriorityQueue) = - while not other.IsEmpty do - let p, v = other.Next - this.Add p v*) + member this.Min = + lowest + 1 |> byte member this.Add (value: byte) (p: Point) = let vi = int value - if this.IsEmpty + if vi > highest then - highest <- int value - q.Insert(0, null) - elif vi > highest - then - for i in highest .. vi - 1 do - q.Insert(0, null) highest <- vi - elif highest - vi >= q.Count + if vi <= lowest then - for i in 0 .. highest - vi - q.Count do - q.Add(null) + lowest <- vi - 1 + + q.[vi].Add(p) |> ignore - let pos = highest - vi - if q.[pos] = null + member this.Remove (value: byte) (p: Point) = + let vi = int value + if q.[vi].Remove(p) && q.[vi].Count = 0 then - q.[pos] <- HashSet([p]) - else - q.[pos].Add(p) |> ignore + if vi = highest + then + highest <- highest - 1 + while highest > lowest && q.[highest].Count = 0 do + highest <- highest - 1 + elif vi - 1 = lowest + then + lowest <- lowest + 1 + while lowest < highest && q.[lowest + 1].Count = 0 do + lowest <- lowest + 1 + if highest = lowest // The queue is now empty. + then + highest <- -1 + lowest <- size member this.IsEmpty = - q.Count = 0 + highest = -1 member this.Clear () = - while highest >= 0 do - q.[highest] <- null + while highest > lowest do + q.[highest].Clear() highest <- highest - 1 + highest <- -1 + lowest <- size +type private AreaState = + | Removed = 1 + | Unprocessed = 2 + | Validated = 3 -type MaximaState = Uncertain | Validated | TooBig -type Maxima = { - elements : HashSet - mutable intensity: byte option - mutable state: MaximaState } +type private AreaOperation = + | Opening = 1 + | Closing = 2 +[] +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 maxima = findMaxima img |> Seq.map (fun m -> { elements = m; intensity = None; state = Uncertain }) |> List.ofSeq - let toValidated = Stack(maxima) + let areas = List((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area) - while toValidated.Count > 0 do - let m = toValidated.Pop() - if m.elements.Count <= area + let pixels: Area[,] = Array2D.create h w null + for m in areas do + for e in m.Elements do + pixels.[e.Y, e.X] <- m + + let queue = PriorityQueue() + + let addEdgeToQueue (elements: Points) = + for p in elements do + 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 + let m = areas.[i] + if m.Elements.Count <= area && m.State <> AreaState.Removed then - let queue = - let q = PriorityQueue() - let firstElements = HashSet() - for p in m.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 (m.elements.Contains(p')) && not (firstElements.Contains(p')) - then - firstElements.Add(p') |> ignore - q.Add (img.Data.[ni, nj, 0]) p' - q + queue.Clear() + addEdgeToQueue m.Elements - let mutable intensity = queue.Max - let nextElements = HashSet() + 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. then - if m.elements.Count + nextElements.Count + 1 > area + if m.Elements.Count + nextElements.Count + 1 > area then - m.state <- Validated - m.intensity <- Some intensity + m.State <- AreaState.Validated + m.Intensity <- Some intensity stop <- true else nextElements.Add(p) |> ignore - elif intensity' < intensity + + elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity then - m.elements.UnionWith(nextElements) - if m.elements.Count = area + m.Elements.UnionWith(nextElements) + for e in nextElements do + pixels.[e.Y, e.X] <- m + + if m.Elements.Count = area then - m.state <- Validated - m.intensity <- Some (intensity') + m.State <- AreaState.Validated + m.Intensity <- Some (intensity') stop <- true else intensity <- intensity' nextElements.Clear() nextElements.Add(p) |> ignore - else // i' > i - seq { - for m' in maxima do - if m' <> m && m'.elements.Contains(p) then - if m'.elements.Count + m.elements.Count <= area - then - m'.state <- Uncertain - m'.elements.UnionWith(m.elements) - if not <| toValidated.Contains m' // FIXME: Maybe use state instead of scanning the whole list. - then - toValidated.Push(m') - stop <- true - yield false - } |> Seq.forall id |> ignore - if not stop + else + let m' = pixels.[p.Y, p.X] + if m' <> null + then + if m'.Elements.Count + m.Elements.Count <= area + then + m'.State <- AreaState.Removed + for e in m'.Elements do + pixels.[e.Y, e.X] <- m + queue.Remove imgData.[e.Y, e.X, 0] e + addEdgeToQueue m'.Elements + m.Elements.UnionWith(m'.Elements) + let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min + if intensityMax <> intensity + then + intensity <- intensityMax + nextElements.Clear() + merged <- true + + if not merged then - m.state <- Validated - m.intensity <- Some (intensity) + m.State <- AreaState.Validated + m.Intensity <- Some (intensity) stop <- true - if not stop + 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 <- Validated - m.intensity <- Some (intensity) - stop <- true - elif not (m.elements.Contains(p')) && not (nextElements.Contains(p')) - then - queue.Add (img.Data.[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 - if m.elements.Count + nextElements.Count <= area + if m.Elements.Count + nextElements.Count <= area then - m.state <- Validated - m.intensity <- Some intensity' - m.elements.UnionWith(nextElements) + m.State <- AreaState.Validated + m.Intensity <- Some intensity' + m.Elements.UnionWith(nextElements) stop <- true - for m in maxima do - if m.state = Validated + for m in areas do + if m.State = AreaState.Validated then - match m.intensity with + match m.Intensity with | Some i -> - for p in m.elements do - img.Data.[p.Y, p.X, 0] <- i + for p in m.Elements do + imgData.[p.Y, p.X, 0] <- i | _ -> () () +let areaOpen (img: Image) (area: int) = + areaOperation img area AreaOperation.Opening + +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 + 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) = @@ -331,12 +720,6 @@ let thin (mat: Matrix) = data2 <- tmp -// FIXME: replace by a queue or stack. -let pop (l: List<'a>) : 'a = - let n = l.[l.Count - 1] - l.RemoveAt(l.Count - 1) - n - // Remove all 8-connected pixels with an area equal or greater than 'areaSize'. // Modify 'mat' in place. let removeArea (mat: Matrix) (areaSize: int) = @@ -350,7 +733,7 @@ let removeArea (mat: Matrix) (areaSize: int) = ( 0, -1) // p8 (-1, -1) |] // p9 - let mat' = new Matrix(mat.Size) + use mat' = new Matrix(mat.Size) let w = mat'.Width let h = mat'.Height mat.CopyTo(mat') @@ -362,37 +745,37 @@ let removeArea (mat: Matrix) (areaSize: int) = for j in 0..w-1 do if data'.[i, j] = 1uy then - let neighborhood = List<(int*int)>() - let neighborsToCheck = List<(int*int)>() - neighborsToCheck.Add((i, j)) + let neighborhood = List() + let neighborsToCheck = Stack() + neighborsToCheck.Push(Point(j, i)) data'.[i, j] <- 0uy while neighborsToCheck.Count > 0 do - let (ci, cj) = pop neighborsToCheck - neighborhood.Add((ci, cj)) + let n = neighborsToCheck.Pop() + neighborhood.Add(n) for (ni, nj) in neighbors do - let pi = ci + ni - let pj = cj + nj + let pi = n.Y + ni + let pj = n.X + nj if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy then - neighborsToCheck.Add((pi, pj)) + neighborsToCheck.Push(Point(pj, pi)) data'.[pi, pj] <- 0uy if neighborhood.Count <= areaSize then - for (ni, nj) in neighborhood do - data.[ni, nj] <- 0uy + for n in neighborhood do + data.[n.Y, n.X] <- 0uy let connectedComponents (img: Image) (startPoints: List) : List = let w = img.Width let h = img.Height - let pointChecked = HashSet() - let pointToCheck = List(startPoints); + let pointChecked = Points() + let pointToCheck = Stack(startPoints); let data = img.Data while pointToCheck.Count > 0 do - let next = pop pointToCheck + let next = pointToCheck.Pop() pointChecked.Add(next) |> ignore for ny in -1 .. 1 do for nx in -1 .. 1 do @@ -401,45 +784,38 @@ let connectedComponents (img: Image) (startPoints: List) : Li let p = Point(next.X + nx, next.Y + ny) if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p) then - pointToCheck.Add(p) + pointToCheck.Push(p) 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 then - img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias) + 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) else - let windowPosX = e.Cx - e.A - 5.0 - let gapX = windowPosX - (float (int windowPosX)) + let windowPosX = e.Cx - e.A - 5.f + let gapX = windowPosX - (float32 (int windowPosX)) - let windowPosY = e.Cy - e.A - 5.0 - let gapY = windowPosY - (float (int windowPosY)) + let windowPosY = e.Cy - e.A - 5.f + let gapY = windowPosY - (float32 (int windowPosY)) - let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int) + let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int) img.ROI <- roi if roi = img.ROI // We do not display ellipses touching the edges (FIXME) then use i = new Image<'TColor, 'TDepth>(img.ROI.Size) - i.Draw(Ellipse(PointF(float32 <| (e.A + 5. + gapX) , float32 <| (e.A + 5. + gapY)), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias) + i.Draw(Ellipse(PointF(float32 <| (e.A + 5.f + gapX) , float32 <| (e.A + 5.f + gapY)), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img) img.ROI <- Rectangle.Empty