X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=cee21c77d5c60fa26c8e25f4c2afa34c15056543;hp=f567e682decb86c97b81aefeb8c7a44eeab1f347;hb=d0c85068bb98a7999ed994f02669befa70edd5f9;hpb=bef2e9f0bf1bba21d4c988fdf654c2dc303ec34a diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index f567e68..cee21c7 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -36,11 +36,13 @@ let suppressMConnections (img: Matrix) = 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 + 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 + 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 = @@ -70,7 +72,8 @@ let findEdges (img: Image) : Matrix * Image * yGradientData.[h - 1, c, 0] <- 0.0 use magnitudes = new Matrix(xGradient.Size) - CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles). + use angles = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles). let thresholdHigh, thresholdLow = let sensibility = 0.1 @@ -80,7 +83,6 @@ let findEdges (img: Image) : Matrix * Image * // Non-maximum suppression. use nms = new Matrix(xGradient.Size) - nms.SetValue(1.0) for i in 0 .. h - 1 do nms.Data.[i, 0] <- 0uy @@ -94,25 +96,45 @@ let findEdges (img: Image) : Matrix * Image * for j in 1 .. w - 2 do let vx = xGradient.Data.[i, j, 0] let vy = yGradient.Data.[i, j, 0] - let angle = - let a = atan2 vy vx - if a < 0.0 then 2. * Math.PI + a else a - - let mNeigbors (sign: int) : float = - if angle < Math.PI / 8. || angle >= 15.0 * Math.PI / 8. then magnitudes.Data.[i, j + sign] - elif angle < 3.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j + sign] - elif angle < 5.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j] - elif angle < 7.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j - sign] - elif angle < 9.0 * Math.PI / 8. then magnitudes.Data.[i, j - sign] - elif angle < 11.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j - sign] - elif angle < 13.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j] - else magnitudes.Data.[i - sign, j + sign] - - let m = magnitudes.Data.[i, j] - if m < mNeigbors 1 || m < mNeigbors -1 || m < thresholdLow then - nms.Data.[i, j] <- 0uy + if vx <> 0. || vy <> 0. + then + let angle = angles.[i, j] - // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection). + let vx', vy' = abs vx, abs vy + let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy' + let ratio1 = 1. - ratio2 + + let mNeigbors (sign: int) : float = + if angle < Math.PI / 4. + then + ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign] + elif angle < Math.PI / 2. + then + ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j] + elif angle < 3.0 * Math.PI / 4. + then + ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign] + elif angle < Math.PI + then + ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign] + elif angle < 5. * Math.PI / 4. + then + ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign] + elif angle < 3. * Math.PI / 2. + then + ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j] + elif angle < 7. * Math.PI / 4. + then + ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign] + else + ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign] + + let m = magnitudes.Data.[i, j] + if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 + then + nms.Data.[i, j] <- 1uy + + // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection). let edges = new Matrix(xGradient.Size) @@ -120,7 +142,8 @@ let findEdges (img: Image) : Matrix * Image * let toVisit = Stack() for i in 0 .. h - 1 do for j in 0 .. w - 1 do - if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh then + if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh + then nms.Data.[i, j] <- 0uy toVisit.Push(Point(j, i)) while toVisit.Count > 0 do @@ -128,14 +151,15 @@ let findEdges (img: Image) : Matrix * Image * edges.Data.[p.Y, p.X] <- 1uy for i' in -1 .. 1 do for j' in -1 .. 1 do - if i' <> 0 || j' <> 0 then + 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 && nms.Data.[ni, nj] = 1uy then + if ni >= 0 && ni < h && nj >= 0 && nj < w && nms.Data.[ni, nj] = 1uy + then nms.Data.[ni, nj] <- 0uy toVisit.Push(Point(nj, ni)) - edges, xGradient, yGradient @@ -146,7 +170,7 @@ let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : type Points = HashSet -let drawPoints (img: Image) (points: Points) (intensity: byte) = +let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = for p in points do img.Data.[p.Y, p.X, 0] <- intensity @@ -154,7 +178,7 @@ type ExtremumType = | Maxima = 1 | Minima = 2 -let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumerable = +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 |] @@ -219,10 +243,11 @@ let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumer result.Select(fun l -> Points(l)) -let findMaxima (img: Image) : IEnumerable = +let findMaxima (img: Image) : IEnumerable = findExtremum img ExtremumType.Maxima -let findMinima (img: Image) : IEnumerable = + +let findMinima (img: Image) : IEnumerable = findExtremum img ExtremumType.Minima @@ -473,6 +498,117 @@ let areaOpen (img: Image) (area: int) = let areaClose (img: Image) (area: int) = areaOperation img area AreaOperation.Closing +[] +type Island (cmp: IComparer) = + member val Shore = Heap.Heap(cmp) with get + member val Level = 0.f with get, set + member val Surface = 0 with get, set + + +let private areaOperationF (img: Image) (area: int) (op: AreaOperation) = + let w = img.Width + let h = img.Height + let earth = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let comparer = if op = AreaOperation.Opening + then { new IComparer with member this.Compare(v1, v2) = v1.CompareTo(v2) } + else { new IComparer with member this.Compare(v1, v2) = v2.CompareTo(v1) } + + let ownership: Island[,] = Array2D.create h w null + + // Initialize islands with their shore. + let islands = List() + let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima + for e in extremum do + let island = + let p = e.First() + Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count) + islands.Add(island) + let shorePoints = Points() + for p in e do + ownership.[p.Y, p.X] <- island + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let neighbor = Point(nj, ni) + if ni >= 0 && ni < h && nj >= 0 && nj < w && ownership.[ni, nj] = null && not (shorePoints.Contains(neighbor)) + then + shorePoints.Add(neighbor) |> ignore + island.Shore.Add earth.[ni, nj, 0] neighbor + + for island in islands do + let mutable stop = island.Shore.IsEmpty + + // 'true' if 'p' is owned or adjacent to 'island'. + let ownedOrAdjacent (p: Point) : bool = + ownership.[p.Y, p.X] = island || + (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) || + (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) || + (p.X > 0 && ownership.[p.Y, p.X - 1] = island) || + (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island) + + while not stop && island.Surface < area do + let level, next = island.Shore.Max + let other = ownership.[next.Y, next.X] + if other = island // During merging, some points on the shore may be owned by the island itself -> ignored. + then + island.Shore.RemoveNext () + else + if other <> null + then // We touching another island. + if island.Surface + other.Surface >= area + then + stop <- true + else // We can merge 'other' into 'surface'. + island.Surface <- island.Surface + other.Surface + island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level + for l, p in other.Shore do + let mutable currentY = p.Y + 1 + while currentY < h && ownership.[currentY, p.X] = other do + ownership.[currentY, p.X] <- island + currentY <- currentY + 1 + island.Shore.Add l p + other.Shore.Clear() + + elif comparer.Compare(level, island.Level) > 0 + then + stop <- true + else + island.Shore.RemoveNext () + for i, j in se do + let ni = i + next.Y + let nj = j + next.X + if ni < 0 || ni >= h || nj < 0 || nj >= w + then + island.Surface <- Int32.MaxValue + stop <- true + else + let neighbor = Point(nj, ni) + if not <| ownedOrAdjacent neighbor + then + island.Shore.Add earth.[ni, nj, 0] neighbor + if not stop + then + ownership.[next.Y, next.X] <- island + island.Level <- level + island.Surface <- island.Surface + 1 + + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + let island = ownership.[i, j] + if island <> null + then + earth.[i, j, 0] <- island.Level + () + + +let areaOpenF (img: Image) (area: int) = + areaOperationF img area AreaOperation.Opening + +let areaCloseF (img: Image) (area: int) = + areaOperationF img area AreaOperation.Closing + // A simpler algorithm than 'areaOpen' but slower. let areaOpen2 (img: Image) (area: int) = let w = img.Width @@ -580,12 +716,6 @@ let thin (mat: Matrix) = data2 <- tmp -// FIXME: replace by a queue or stack. -let pop (l: List<'a>) : 'a = - let n = l.[l.Count - 1] - l.RemoveAt(l.Count - 1) - n - // Remove all 8-connected pixels with an area equal or greater than 'areaSize'. // Modify 'mat' in place. let removeArea (mat: Matrix) (areaSize: int) = @@ -611,37 +741,37 @@ let removeArea (mat: Matrix) (areaSize: int) = for j in 0..w-1 do if data'.[i, j] = 1uy then - let neighborhood = List<(int*int)>() - let neighborsToCheck = List<(int*int)>() - neighborsToCheck.Add((i, j)) + let neighborhood = List() + let neighborsToCheck = Stack() + neighborsToCheck.Push(Point(j, i)) data'.[i, j] <- 0uy while neighborsToCheck.Count > 0 do - let (ci, cj) = pop neighborsToCheck - neighborhood.Add((ci, cj)) + let n = neighborsToCheck.Pop() + neighborhood.Add(n) for (ni, nj) in neighbors do - let pi = ci + ni - let pj = cj + nj + let pi = n.Y + ni + let pj = n.X + nj if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy then - neighborsToCheck.Add((pi, pj)) + neighborsToCheck.Push(Point(pj, pi)) data'.[pi, pj] <- 0uy if neighborhood.Count <= areaSize then - for (ni, nj) in neighborhood do - data.[ni, nj] <- 0uy + for n in neighborhood do + data.[n.Y, n.X] <- 0uy let connectedComponents (img: Image) (startPoints: List) : List = let w = img.Width let h = img.Height let pointChecked = Points() - let pointToCheck = List(startPoints); + let pointToCheck = Stack(startPoints); let data = img.Data while pointToCheck.Count > 0 do - let next = pop pointToCheck + let next = pointToCheck.Pop() pointChecked.Add(next) |> ignore for ny in -1 .. 1 do for nx in -1 .. 1 do @@ -650,7 +780,7 @@ let connectedComponents (img: Image) (startPoints: List) : Li let p = Point(next.X + nx, next.Y + ny) if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p) then - pointToCheck.Add(p) + pointToCheck.Push(p) List(pointChecked)