module ParasitemiaCore.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