X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemiaCore%2FImgTools%2FMorpho.fs;fp=Parasitemia%2FParasitemiaCore%2FImgTools%2FMorpho.fs;h=b3a2e759a7e92315c3c23133f6e4e54d464528f6;hp=0000000000000000000000000000000000000000;hb=3f8b0d281b3058faf23dbd0363de440bd04c6574;hpb=e3842630f4d36c5ea8c8a0c3d4762684e1c510f4 diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs b/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs new file mode 100644 index 0000000..b3a2e75 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs @@ -0,0 +1,670 @@ +module ParasitemiaCore.Morpho + +open System +open System.Drawing +open System.Collections.Generic +open System.Linq + +open Emgu.CV +open Emgu.CV.Structure + +open Types + +/// +/// Remove M-adjacent pixels. It may be used after thinning. +/// +let suppressMAdjacency (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 + +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 + | _ -> () + () + +/// +/// Area opening on byte image. +/// +let areaOpen (img: Image) (area: int) = + areaOperation img area AreaOperation.Opening + +/// +/// Area closing on byte image. +/// +let areaClose (img: Image) (area: int) = + areaOperation img area AreaOperation.Closing + +// A simpler algorithm than 'areaOpen' on byte image 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 + +[] +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 + member this.IsInfinite = this.Surface = Int32.MaxValue + +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.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0 + then + stop <- true + else // We can merge 'other' into 'surface'. + island.Surface <- island.Surface + other.Surface + island.Level <- other.Level + // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.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 + | _ -> () + () + +/// +/// Area opening on float image. +/// +let areaOpenF (img: Image) (area: int) = + areaOperationF img [ area, () ] None AreaOperation.Opening + +/// +/// Area closing on float image. +/// +let areaCloseF (img: Image) (area: int) = + areaOperationF img [ area, () ] None AreaOperation.Closing + +/// +/// Area closing on float image with different areas. Given areas must be sorted increasingly. +/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference +/// Between the previous and the current closing. +/// +let areaOpenFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Opening + +/// +/// Same as 'areaOpenFWithFun' for closing operation. +/// +let areaCloseFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Closing + +/// +/// Zhang and Suen thinning 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) : Points = + 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) + + pointChecked