X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;fp=Parasitemia%2FParasitemia%2FImgTools.fs;h=0000000000000000000000000000000000000000;hp=096fd94ec7e9736cd73a60200ce95be017ca6c35;hb=4bfa3cbdc6145e6944f02e24829ab2ef3a851ac1;hpb=48ecdfc43001c444eff6ad442986049384674af2 diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs deleted file mode 100644 index 096fd94..0000000 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ /dev/null @@ -1,973 +0,0 @@ -module ImgTools - -open System -open System.Drawing -open System.Collections.Generic -open System.Linq - -open Emgu.CV -open Emgu.CV.Structure - -open Heap -open Const -open Utils - -// Normalize image values between 0uy and 255uy. -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.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 - -type Histogram = { data: int[]; total: int; sum: int; min: float32; max: float32 } - -let histogramImg (img: Image) (nbSamples: int) : Histogram = - let imgData = img.Data - - let min, max = - let min = ref [| 0.0 |] - let minLocation = ref <| [| Point() |] - let max = ref [| 0.0 |] - let maxLocation = ref <| [| Point() |] - img.MinMax(min, max, minLocation, maxLocation) - float32 (!min).[0], float32 (!max).[0] - - let bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for i in 0 .. img.Height - 1 do - for j in 0 .. img.Width - 1 do - let p = bin imgData.[i, j, 0] - data.[p] <- data.[p] + 1 - - { data = data; total = img.Height * img.Width; sum = Array.sum data; min = min; max = max } - -let histogramMat (mat: Matrix) (nbSamples: int) : Histogram = - let matData = mat.Data - - let min, max = - let min = ref 0.0 - let minLocation = ref <| Point() - let max = ref 0.0 - let maxLocation = ref <| Point() - mat.MinMax(min, max, minLocation, maxLocation) - float32 !min, float32 !max - - let bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for i in 0 .. mat.Height - 1 do - for j in 0 .. mat.Width - 1 do - let p = bin matData.[i, j] - data.[p] <- data.[p] + 1 - - { data = data; total = mat.Height * mat.Width; sum = Array.sum data; min = min; max = max } - -let histogram (values: float32 seq) (nbSamples: int) : Histogram = - let mutable min = Single.MaxValue - let mutable max = Single.MinValue - let mutable n = 0 - - for v in values do - n <- n + 1 - if v < min then min <- v - if v > max then max <- v - - let bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for v in values do - let p = bin v - data.[p] <- data.[p] + 1 - - { data = data; total = n; sum = Array.sum data; min = min; max = max } - -let otsu (hist: Histogram) : float32 * float32 * float32 = - let mutable sumB = 0 - let mutable wB = 0 - let mutable maximum = 0.0 - let mutable level = 0 - let sum = hist.data |> Array.mapi (fun i v -> i * v) |> Array.sum |> float - - for i in 0 .. hist.data.Length - 1 do - wB <- wB + hist.data.[i] - if wB <> 0 - then - let wF = hist.total - wB - if wF <> 0 - then - sumB <- sumB + i * hist.data.[i] - let mB = (float sumB) / (float wB) - let mF = (sum - float sumB) / (float wF) - let between = (float wB) * (float wF) * (mB - mF) ** 2.; - if between >= maximum - then - level <- i - maximum <- between - - let mean1 = - let mutable sum = 0 - let mutable nb = 0 - for i in 0 .. level - 1 do - sum <- sum + i * hist.data.[i] - nb <- nb + hist.data.[i] - (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) - - let mean2 = - let mutable sum = 0 - let mutable nb = 0 - for i in level + 1 .. hist.data.Length - 1 do - sum <- sum + i * hist.data.[i] - nb <- nb + hist.data.[i] - (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) - - let toFloat l = - float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min - - toFloat level, toFloat mean1, toFloat mean2 - -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.0f - use magnitudesByte = magnitudes.Convert() - let threshold, _, _ = otsu (histogramMat magnitudes 300) - 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> = - let size = 2 * int (ceil (4.0 * standardDeviation)) + 1 - img.SmoothGaussian(size, size, standardDeviation, standardDeviation) - -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 suppress: bool[,] = Array2D.zeroCreate h w - - let result = List>() - - let flood (start: Point) : List> = - let sameLevelToCheck = Stack() - let betterLevelToCheck = Stack() - betterLevelToCheck.Push(start) - - let result' = List>() - - while betterLevelToCheck.Count > 0 do - let p = betterLevelToCheck.Pop() - if not suppress.[p.Y, p.X] - then - suppress.[p.Y, p.X] <- true - sameLevelToCheck.Push(p) - let current = List() - - let mutable betterExists = false - - while sameLevelToCheck.Count > 0 do - let p' = sameLevelToCheck.Pop() - let currentLevel = imgData.[p'.Y, p'.X, 0] - current.Add(p') |> ignore - 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 - betterExists <- true - if notSuppressed - then - betterLevelToCheck.Push(Point(nj, ni)) - - if not betterExists - then - result'.Add(current) - result' - - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - let maxima = flood (Point(j, i)) - if maxima.Count > 0 - 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()) - let mutable highest = -1 // Value of the first elements of 'q'. - let mutable lowest = size - - member this.NextMax () : byte * Point = - if this.IsEmpty - then - invalidOp "Queue is empty" - else - let l = q.[highest] - let next = l.First() - l.Remove(next) |> ignore - let value = byte highest - - if l.Count = 0 - then - highest <- highest - 1 - 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.Min = - lowest + 1 |> byte - - member this.Add (value: byte) (p: Point) = - let vi = int value - - if vi > highest - then - highest <- vi - if vi <= lowest - then - lowest <- vi - 1 - - q.[vi].Add(p) |> ignore - - member this.Remove (value: byte) (p: Point) = - let vi = int value - if q.[vi].Remove(p) && q.[vi].Count = 0 - then - 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 = - highest = -1 - - member this.Clear () = - while highest > lowest do - q.[highest].Clear() - highest <- highest - 1 - highest <- -1 - lowest <- size - -type private AreaState = - | Removed = 1 - | Unprocessed = 2 - | Validated = 3 - -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 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((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 - 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 - queue.Clear() - addEdgeToQueue m.Elements - - 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 = 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 - then - m.State <- AreaState.Validated - m.Intensity <- Some intensity - stop <- true - else - nextElements.Add(p) |> ignore - - elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity - then - m.Elements.UnionWith(nextElements) - for e in nextElements do - pixels.[e.Y, e.X] <- m - - if m.Elements.Count = area - then - m.State <- AreaState.Validated - m.Intensity <- Some (intensity') - stop <- true - else - intensity <- intensity' - nextElements.Clear() - nextElements.Add(p) |> ignore - - else - match pixels.[p.Y, p.X] with - | null -> () - | m' -> - 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 <- AreaState.Validated - m.Intensity <- Some (intensity) - stop <- true - - if not stop && not merged - then - 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 - then - m.State <- AreaState.Validated - m.Intensity <- Some intensity' - m.Elements.UnionWith(nextElements) - stop <- true - - for m in areas do - if m.State = AreaState.Validated - then - match m.Intensity with - | Some 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) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (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 && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor)) - then - shorePoints.Add(neighbor) |> ignore - island.Shore.Add earth.[ni, nj, 0] neighbor - - for area, obj in areas do - for island in islands do - let mutable stop = island.Shore.IsEmpty - - // 'true' if 'p' is owned or adjacent to 'island'. - let inline 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 not <| Object.ReferenceEquals(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 - - let mutable diff = 0.f - - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - match ownership.[i, j] with - | null -> () - | island -> - let l = island.Level - diff <- diff + l - earth.[i, j, 0] - earth.[i, j, 0] <- l - - match f with - | Some f' -> f' obj diff - | _ -> () - () - -let areaOpenF (img: Image) (area: int) = - areaOperationF img [ area, () ] None AreaOperation.Opening - -let areaCloseF (img: Image) (area: int) = - areaOperationF img [ area, () ] None AreaOperation.Closing - -let areaOpenFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = - areaOperationF img areas (Some f) AreaOperation.Opening - -let areaCloseFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = - areaOperationF img areas (Some f) 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) = - let w = mat.Width - let h = mat.Height - let mutable data1 = mat.Data - let mutable data2 = Array2D.copy data1 - - let mutable pixelChanged = true - let mutable oddIteration = true - - while pixelChanged do - pixelChanged <- false - for i in 0..h-1 do - for j in 0..w-1 do - if data1.[i, j] = 1uy - then - let p2 = if i = 0 then 0uy else data1.[i-1, j] - let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1] - let p4 = if j = w-1 then 0uy else data1.[i, j+1] - let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1] - let p6 = if i = h-1 then 0uy else data1.[i+1, j] - let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1] - let p8 = if j = 0 then 0uy else data1.[i, j-1] - let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1] - - let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 - if sumNeighbors >= 2uy && sumNeighbors <= 6uy && - (if p2 = 0uy && p3 = 1uy then 1 else 0) + - (if p3 = 0uy && p4 = 1uy then 1 else 0) + - (if p4 = 0uy && p5 = 1uy then 1 else 0) + - (if p5 = 0uy && p6 = 1uy then 1 else 0) + - (if p6 = 0uy && p7 = 1uy then 1 else 0) + - (if p7 = 0uy && p8 = 1uy then 1 else 0) + - (if p8 = 0uy && p9 = 1uy then 1 else 0) + - (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 && - if oddIteration - then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy - else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy - then - data2.[i, j] <- 0uy - pixelChanged <- true - else - data2.[i, j] <- 0uy - - oddIteration <- not oddIteration - let tmp = data1 - 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) (areaSize: int) = - let neighbors = [| - (-1, 0) // p2 - (-1, 1) // p3 - ( 0, 1) // p4 - ( 1, 1) // p5 - ( 1, 0) // p6 - ( 1, -1) // p7 - ( 0, -1) // p8 - (-1, -1) |] // p9 - - use mat' = new Matrix(mat.Size) - let w = mat'.Width - let h = mat'.Height - mat.CopyTo(mat') - - let data = mat.Data - let data' = mat'.Data - - for i in 0..h-1 do - for j in 0..w-1 do - if data'.[i, j] = 1uy - then - let neighborhood = List() - let neighborsToCheck = Stack() - neighborsToCheck.Push(Point(j, i)) - data'.[i, j] <- 0uy - - while neighborsToCheck.Count > 0 do - let n = neighborsToCheck.Pop() - neighborhood.Add(n) - for (ni, nj) in neighbors do - 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.Push(Point(pj, pi)) - data'.[pi, pj] <- 0uy - if neighborhood.Count <= areaSize - then - 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 = Points() - let pointToCheck = Stack(startPoints); - - let data = img.Data - - while pointToCheck.Count > 0 do - let next = pointToCheck.Pop() - pointChecked.Add(next) |> ignore - for ny in -1 .. 1 do - for nx in -1 .. 1 do - if ny <> 0 && nx <> 0 - then - 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.Push(p) - - List(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) - else - let windowPosX = e.Cx - e.A - 5.f - let gapX = windowPosX - (float32 (int windowPosX)) - - let windowPosY = e.Cy - e.A - 5.f - let gapY = windowPosY - (float32 (int windowPosY)) - - 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.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 - -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) (drawCellContent: bool) (c: Types.Cell) = - if drawCellContent - then - let colorB = rngCell.Next(20, 70) - let colorG = rngCell.Next(20, 70) - let colorR = rngCell.Next(20, 70) - - for y in 0 .. c.elements.Height - 1 do - for x in 0 .. c.elements.Width - 1 do - if c.elements.[y, x] > 0uy - then - let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2 - let b = img.Data.[y + dy, x + dx, 0] |> int - let g = img.Data.[y + dy, x + dx, 1] |> int - let r = img.Data.[y + dy, x + dx, 2] |> int - img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB) - img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG) - img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR) - - let crossColor, crossColor2 = - match c.cellClass with - | Types.HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.) - | Types.InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.) - | Types.Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.) - - drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2 - drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2 - - drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1 - drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1 - - -let drawCells (img: Image) (drawCellContent: bool) (cells: Types.Cell list) = - List.iter (fun c -> drawCell img drawCellContent c) cells \ No newline at end of file