* Some adjustments.
let infection = parasites.infection.Copy() // To avoid to modify the parameter.
// This is the minimum window size to check if other ellipses touch 'e'.
- let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMaxRadius) * config.Parameters.scale
- KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius) * config.Parameters.scale
- KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius) * config.Parameters.scale
- KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) * config.Parameters.scale }
+ let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMaxRadius)
+ KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius)
+ KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius)
+ KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) }
// The minimum window to contain a given ellipse.
let ellipseWindow (e: Ellipse) =
// Return 'true' if the point 'p' is owned by e.
// The lines represents all intersections with other ellipses.
- let pixelOwnedByE (p: PointD) (e: Ellipse) (lines: Line list) =
+ let pixelOwnedByE (p: PointD) (e: Ellipse) (others: (Ellipse * Line) list) =
e.Contains p.X p.Y &&
seq {
let c = PointD(e.Cx, e.Cy)
- for d1 in lines do
+ for e', d1 in others do
let d2 = Utils.lineFromTwoPoints c p
+ let c' = PointD(e'.Cx, e'.Cy)
+ let v = pointFromTwoLines d1 (lineFromTwoPoints c c')
+ let case1 = sign (v.X - c.X) <> sign (v.X - c'.X) || Utils.squaredDistanceTwoPoints v c > Utils.squaredDistanceTwoPoints v c'
if d2.Valid
then
let p' = Utils.pointFromTwoLines d1 d2
- yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p // 'false' -> the point is owned by another ellipse.
+ // Yield 'false' when the point is owned by another ellipse.
+ if case1
+ then
+ yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p
+ else
+ yield sign (c.X - p.X) = sign (c.X - p'.X) && Utils.squaredDistanceTwoPoints c p' < Utils.squaredDistanceTwoPoints c p
else
- yield true
+ yield case1
} |> Seq.forall id
let ellipses = ellipses |> List.map EllipseFlaggedKd
- // 1) Associate touching ellipses with each ellipses.
+ // 1) Associate touching ellipses with each ellipses and remove ellipse with more than two intersections.
let tree = KdTree.Tree.BuildTree ellipses
- let neighbors (e: Ellipse) : (EllipseFlaggedKd * PointD * PointD) list =
- tree.Search (searchRegion e)
- // We only keep the ellipses touching 'e'.
- |> List.choose (fun otherE ->
- match EEOver.EEOverlapArea e otherE with
- | Some (area, px, py) when area > 0.0 && px.Length >= 2 && py.Length >= 2 ->
- Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1]))
- | _ ->
- None )
+ let neighbors (e: EllipseFlaggedKd) : (EllipseFlaggedKd * PointD * PointD) list =
+ if not e.Removed
+ then
+ tree.Search (searchRegion e)
+ // We only keep the ellipses touching 'e'.
+ |> List.choose (fun otherE ->
+ match EEOver.EEOverlapArea e otherE with
+ | Some (_, px, _) when px.Length > 2 ->
+ otherE.Removed <- true
+ None
+ | Some (area, px, py) when area > 0.0 && px.Length = 2 ->
+ Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1]))
+ | _ ->
+ None )
+ else
+ []
// We reverse the list to get the lower score ellipses first.
let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev
// 2) Remove ellipses with a high standard deviation (high contrast).
- let globalStdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq {
+
+ // CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed.
+
+ let globalStdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation(seq {
for y in 0 .. h - 1 do
for x in 0 .. w - 1 do
yield img.Data.[y, x, 0] |> float })
for e in ellipses do
let minX, minY, maxX, maxY = ellipseWindow e
- let stdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq {
+ let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq {
for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
if e.Contains (float x) (float y)
then
- yield img.Data.[y, x, 0] |> float })
+ yield float img.Data.[y, x, 0] })
- if stdDiviation > globalStdDiviation * config.Parameters.standardDeviationMaxRatio
- then
+ if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then
e.Removed <- true
// 3) Remove ellipses touching the edges.
for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do
for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do
let p = PointD(float x, float y)
- if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2)))
+ if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
then
area <- area + 1
for y in minY .. maxY do
for x in minX .. maxX do
let p = PointD(float x, float y)
- if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2)))
+ if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2)))
then
elements.[y-minY, x-minX] <- 1uy
nbElement <- nbElement + 1
infectionLevel: float // > 1
standardDeviationMaxRatio: float // The standard deviation of the pixel values of a cell can't be greater than standardDeviationMaxRatio * global standard deviation
- minimumCellArea: float // Factor of RBC area.
+ minimumCellArea: float // Factor of the nominal RBC area.
}
type Config (param: Parameters) =
let find (edges: Matrix<byte>)
- (xDir: Image<Gray, float>)
- (yDir: Image<Gray, float>)
+ (xGradient: Image<Gray, float>)
+ (yGradient: Image<Gray, float>)
(config: Config) : MatchingEllipses =
- let r1, r2 = config.Parameters.scale * config.RBCMinRadius, config.Parameters.scale * config.RBCMaxRadius // FIXME: scale factor should be applied in Config!?
+ let r1, r2 = config.RBCMinRadius, config.RBCMaxRadius
let windowSize = roundInt (config.Parameters.factorWindowSize * r2)
let factorNbPick = config.Parameters.factorNbPick
let currentElements = List<(int * int)>()
let edgesData = edges.Data
- let xDirData = xDir.Data
- let yDirData = yDir.Data
+ let xDirData = xGradient.Data
+ let yDirData = yGradient.Data
let rng = Random(42)
((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert<Gray, byte>()
+let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
+ img.Save(filepath)
+
+
+let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
+ use img = new Image<Gray, 'TDeph>(mat.Size)
+ mat.CopyTo(img)
+ saveImg img filepath
+
+
+let suppressMConnections (img: Matrix<byte>) =
+ 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<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
+ 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).Convert<Gray, float>()
+ let yGradient = img.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+
+ let xGradientData = xGradient.Data
+ let yGradientData = yGradient.Data
+ for r in 0 .. h - 1 do
+ xGradientData.[r, 0, 0] <- 0.0
+ xGradientData.[r, w - 1, 0] <- 0.0
+ yGradientData.[r, 0, 0] <- 0.0
+ yGradientData.[r, w - 1, 0] <- 0.0
+
+ for c in 0 .. w - 1 do
+ xGradientData.[0, c, 0] <- 0.0
+ xGradientData.[h - 1, c, 0] <- 0.0
+ yGradientData.[0, c, 0] <- 0.0
+ yGradientData.[h - 1, c, 0] <- 0.0
+
+ use magnitudes = new Matrix<float>(xGradient.Size)
+ CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles).
+
+ let thresholdHigh, thresholdLow =
+ let sensibility = 0.1
+ use magnitudesByte = magnitudes.Convert<byte>()
+ let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+ threshold + (sensibility * threshold), threshold - (sensibility * threshold)
+
+ // Non-maximum suppression.
+ use nms = new Matrix<byte>(xGradient.Size)
+ nms.SetValue(1.0)
+
+ for i in 0 .. h - 1 do
+ nms.Data.[i, 0] <- 0uy
+ nms.Data.[i, w - 1] <- 0uy
+
+ for j in 0 .. w - 1 do
+ nms.Data.[0, j] <- 0uy
+ nms.Data.[h - 1, j] <- 0uy
+
+ for i in 1 .. h - 2 do
+ 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
+
+ // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection).
+
+ let edges = new Matrix<byte>(xGradient.Size)
+
+ // Histeresis thresholding.
+ let toVisit = Stack<Point>()
+ 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
+ nms.Data.[i, j] <- 0uy
+ toVisit.Push(Point(j, i))
+ while toVisit.Count > 0 do
+ let p = toVisit.Pop()
+ 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
+ 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
+ nms.Data.[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)
for p in points do
img.Data.[p.Y, p.X, 0] <- intensity
-
type ExtremumType =
| Maxima = 1
| Minima = 2
let areaClose (img: Image<Gray, byte>) (area: int) =
areaOperation img area AreaOperation.Closing
+// A simpler algorithm than 'areaOpen' but slower.
let areaOpen2 (img: Image<Gray, byte>) (area: int) =
let w = img.Width
let h = img.Height
List<Point>(pointChecked)
-let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
- img.Save(filepath)
-
-
-let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
- use img = new Image<Gray, 'TDeph>(mat.Size)
- mat.CopyTo(img)
- 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
let filteredGreenWhitoutInfectionFloat = filteredGreenWhitoutInfection.Convert<Gray, float32>()
let filteredGreenWhitoutStainFloat = filteredGreenWhitoutStain.Convert<Gray, float32>()
- use sobelKernel =
- new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
- [ 2.0f; 0.0f; -2.0f ]
- [ 1.0f; 0.0f; -1.0f ]], Point(0, 0))
-
- use xEdges = filteredGreenWhitoutStainFloat.Convolution(sobelKernel).Convert<Gray, float>()
- use yEdges = filteredGreenWhitoutStainFloat.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
-
- let xEdgesData = xEdges.Data
- let yEdgesData = yEdges.Data
- for r in 0..xEdges.Rows-1 do
- xEdgesData.[r, 0, 0] <- 0.0
- xEdgesData.[r, xEdges.Cols-1, 0] <- 0.0
- yEdgesData.[r, 0, 0] <- 0.0
- yEdgesData.[r, xEdges.Cols-1, 0] <- 0.0
-
- for c in 0..xEdges.Cols-1 do
- xEdgesData.[0, c, 0] <- 0.0
- xEdgesData.[xEdges.Rows-1, c, 0] <- 0.0
- yEdgesData.[0, c, 0] <- 0.0
- yEdgesData.[xEdges.Rows-1, c, 0] <- 0.0
-
- use magnitudes = new Matrix<float>(xEdges.Size)
- CvInvoke.CartToPolar(xEdges, yEdges, magnitudes, new Mat()) // Compute the magnitudes (without angles).
-
- let min = ref 0.0
+ (*let min = ref 0.0
let minLocation = ref <| Point()
let max = ref 0.0
let maxLocation = ref <| Point()
use edges = new Matrix<byte>(xEdges.Size)
CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore
- logTime "Finding edges" (fun() -> thin edges)
+ logTime "Finding edges" (fun() -> thin edges)*)
+
+ let edges, xGradient, yGradient = ImgTools.findEdges filteredGreenWhitoutStainFloat
+
logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12)
let allEllipses, ellipses = logTime "Finding ellipses" (fun () ->
- let matchingEllipses = Ellipse.find edges xEdges yEdges config
+ let matchingEllipses = Ellipse.find edges xGradient yGradient config
matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses )
let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config)
saveImg parasites.infection (buildFileName " - parasites - infection.png")
let imgAllEllipses = img.Copy()
- drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.1
+ drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.05
saveImg imgAllEllipses (buildFileName " - ellipses - all.png")
let imgEllipses = img.Copy()
// Do not take in account matching score below this when two ellipses are matched.
let matchingScoreThreshold1 = 0.6
-// All ellipses with a score below this are removed.
-let matchingScoreThreshold2 = 600.
+let scaleOverlapTest = 0.8
type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
let mutable matchingScore = matchingScore
type MatchingEllipses (radiusMin: float) =
let ellipses = List<EllipseScoreFlaggedKd>()
+ // All ellipses with a score below this are removed.
+ let matchingScoreThreshold2 = 25. * radiusMin // 600.
+
member this.Add (e: Ellipse) =
ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
// 3) Sort ellipses by their score.
ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore))
- // 4) Remove ellipses wich have a low score.
+ // 4) Remove ellipses with a low score.
let i = ellipses.BinarySearch(EllipseScoreFlaggedKd(matchingScoreThreshold2, Ellipse(0.0, 0.0, 0.0, 0.0, 0.0)),
{ new IComparer<EllipseScoreFlaggedKd> with
member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs
for other in tree.Search window do
if not other.Removed && other.MatchingScore < e.MatchingScore
then
- if e.Ellipse.Scale(0.8).Contains other.Ellipse.Cx other.Ellipse.Cy
+ if e.Ellipse.Scale(scaleOverlapTest).Contains other.Ellipse.Cx other.Ellipse.Cy
then
other.Removed <- true
ellipses.RemoveAll(fun e -> e.Removed) |> ignore
factorNbPick = 1.0
factorWindowSize = 2.0
- darkStainLevel = 0.3 // Lower -> more sensitive.
+ darkStainLevel = 0.25 // Lower -> more sensitive. 0.3
maxDarkStainRatio = 0.1
infectionArea = 0.012 // 1.2 %
stainLevel = 1.1 // Lower -> more sensitive.
maxStainRatio = 0.12 // 12 %
- standardDeviationMaxRatio = 0.55
+ standardDeviationMaxRatio = 0.58 // 0.55
minimumCellArea = 0.5 })
match mode with
let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A)
PointD(x, y)
+let inline linePassThroughSegment (l: Line) (p1: PointD) (p2: PointD) : bool =
+ let p = pointFromTwoLines l (lineFromTwoPoints p1 p2)
+ sign (p.X - p1.X) <> sign (p.X - p2.X)
+
let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) =
(p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0
projects / master-thesis.git / commitdiff