X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=2f021c8ae04a8a60e429a52b4da7f192b24c2546;hp=cc09405b3477005646cf8f560eea32016e54dbcf;hb=84fdf7404133803fdf0dc867a4da68a144975191;hpb=ba64921fb9a0c36cd8cf802cbf1b2c0f79bc25f6 diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index cc09405..2f021c8 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -3,88 +3,475 @@ open System open System.Drawing open System.Collections.Generic +open System.Linq open Emgu.CV open Emgu.CV.Structure open Utils +open Heap + +// 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 - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() + 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) -// Zhang and Suen algorithm. -// Modify 'mat' in place. -let thin (mat: Matrix) = - 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 +type Points = HashSet + +let drawPoints (img: Image) (points: Points) (intensity: byte) = + 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 + 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 <- 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 + +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.zeroCreate h w - - // Return the list of neighbor values. - let neighborsValues (p1i, p1j) = - Array.map (fun (ni, nj) -> - let pi = p1i + ni - let pj = p1j + nj - if pi < 0 || pi >= h || pj < 0 || pj >= w then 0uy else data1.[pi, pj] - ) neighbors - - // Return the number of 01 pattern in 'values' in a circular way. - let pattern01 (values: byte[]) = - let mutable nb = 0 - let mutable lastValue = 255uy - for v in values do - if lastValue = 0uy && v = 1uy - then - nb <- nb + 1 - lastValue <- v - if lastValue = 0uy && values.[0] = 1uy - then - nb <- nb + 1 - nb - + 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 values = neighborsValues (i, j) - let s = Array.reduce (+) values - if s >= 2uy && s <= 6uy && - pattern01 values = 1 && - (not oddIteration || (values.[0] * values.[2] * values.[4] = 0uy && values.[2] * values.[4] * values.[6] = 0uy)) && // Odd iteration. - (oddIteration || (values.[0] * values.[2] * values.[6] = 0uy && values.[0] * values.[4] * values.[6] = 0uy)) // Even iterations. + 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] <- 1uy else data2.[i, j] <- 0uy - oddIteration <- not oddIteration + oddIteration <- not oddIteration let tmp = data1 data1 <- data2 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) = - let neighbors = [| + let neighbors = [| (-1, 0) // p2 (-1, 1) // p3 ( 0, 1) // p4 @@ -93,29 +480,24 @@ let removeArea (mat: Matrix) (areaSize: int) = ( 1, -1) // p7 ( 0, -1) // p8 (-1, -1) |] // p9 - + let mat' = new Matrix(mat.Size) let w = mat'.Width let h = mat'.Height - mat.CopyTo(mat') - + 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 + if data'.[i, j] = 1uy then let neighborhood = List<(int*int)>() let neighborsToCheck = List<(int*int)>() - neighborsToCheck.Add((i, j)) + neighborsToCheck.Add((i, j)) data'.[i, j] <- 0uy - let pop (l: List<'a>) : 'a = - let n = l.[l.Count - 1] - l.RemoveAt(l.Count - 1) - n - while neighborsToCheck.Count > 0 do let (ci, cj) = pop neighborsToCheck neighborhood.Add((ci, cj)) @@ -131,47 +513,104 @@ let removeArea (mat: Matrix) (areaSize: int) = for (ni, nj) in neighborhood do data.[ni, nj] <- 0uy +let connectedComponents (img: Image) (startPoints: List) : List = + let w = img.Width + let h = img.Height + + let pointChecked = Points() + let pointToCheck = List(startPoints); + + let data = img.Data + + while pointToCheck.Count > 0 do + let next = pop pointToCheck + 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.Add(p) + + List(pointChecked) + -let saveImg (img: Image<'TColor, 'TDepth>) (name: string) = - img.Save("output/" + name) - +let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = + img.Save(filepath) -let saveMat (mat: Matrix<'TDepth>) (name: string) = + +let saveMat (mat: Matrix<'TDepth>) (filepath: string) = use img = new Image(mat.Size) mat.CopyTo(img) - saveImg img name - -(*let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = - let e' = Ellipse(PointF(float32 e.cx, float32 e.cy), SizeF(2.0f * float32 e.a, 2.0f * float32 e.b), float32 e.alpha) - img.Draw(e', color)*) - -let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) = - let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1) - - img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1); - -let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = - let cosAlpha = cos e.alpha - let sinAlpha = sin e.alpha - - let mutable x0 = 0.0 - let mutable y0 = 0.0 - let mutable first_iteration = true - - let n = 40 - let thetaIncrement = 2.0 * Math.PI / (float n) - - for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do - let cosTheta = cos theta - let sinTheta = sin theta - let x = e.cx + cosAlpha * e.a * cosTheta - sinAlpha * e.b * sinTheta - let y = e.cy + sinAlpha * e.a * cosTheta + cosAlpha * e.b * sinTheta - - if not first_iteration + 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) + else + let windowPosX = e.Cx - e.A - 5.0 + let gapX = windowPosX - (float (int windowPosX)) + + let windowPosY = e.Cy - e.A - 5.0 + let gapY = windowPosY - (float (int windowPosY)) + + let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int) + + img.ROI <- roi + if roi = img.ROI // We do not display ellipses touching the edges (FIXME) then - drawLine img color x0 y0 x y - else - first_iteration <- false + 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) + 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 + - x0 <- x - y0 <- y \ No newline at end of file +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