module Classifier
open System
+open System.Collections.Generic
open System.Drawing
open Emgu.CV
open Emgu.CV.Structure
+open Types
+open Utils
-type CellClass = HealthyRBC | InfectedRBC | Peculiar
-type Cell = {
- cellClass: CellClass
- center: Point
- elements: Matrix<byte> }
+type private EllipseFlaggedKd (e: Ellipse) =
+ inherit Ellipse (e.Cx, e.Cy, e.A, e.B, e.Alpha)
-type KdEllipse (e: Types.Ellipse) =
- inherit Types.Ellipse (e.Cx, e.Cy, e.A, e.B, e.Alpha)
-
- interface KdTree.I2DCoords with
+ member val Removed = false with get, set
+
+ interface KdTree.I2DCoords with
member this.X = this.Cx
member this.Y = this.Cy
-
-
-let findCells (ellipses: Types.Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image<Gray, byte>) : Cell list =
+
+
+let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (fg: Image<Gray, byte>) (config: Config.Config) : Cell list =
if ellipses.IsEmpty
then
[]
else
- let tree = KdTree.Tree.buildTree (ellipses |> List.map KdEllipse)
- []
\ No newline at end of file
+ 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.maxRBCSize) * config.scale
+ KdTree.maxX = e.Cx + (e.A + config.maxRBCSize) * config.scale
+ KdTree.minY = e.Cy - (e.A + config.maxRBCSize) * config.scale
+ KdTree.maxY = e.Cy + (e.A + config.maxRBCSize) * config.scale }
+
+ // The minimum window to contain a given ellipse.
+ let ellipseWindow (e: Ellipse) =
+ let cx, cy = roundInt e.Cx, roundInt e.Cy
+ let a = int (e.A + 0.5)
+ cx - a, cy - a, cx + a, cy + a
+
+ // 1) Remove ellipses touching the edges.
+ let ellipsesInside = ellipses |> List.map (fun e ->
+ EllipseFlaggedKd (e, Removed = e.isOutside (float fg.Width) (float fg.Height)))
+
+ // 2) Associate touching ellipses with each ellipses.
+ let tree = KdTree.Tree.BuildTree ellipsesInside
+ 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 ellipsesWithNeigbors = ellipsesInside |> List.choose (fun e -> if e.Removed then None else Some (e, neighbors e))
+
+ // 3) Remove ellipse with a lower percentage of foreground.
+ let fgData = fg.Data
+ for e, neighbors in ellipsesWithNeigbors do
+ let minX, minY, maxX, maxY = ellipseWindow e
+
+ let mutable totalElement = 0
+ let mutable fgElement = 0
+
+ for y in minY .. maxY do
+ for x in minX .. maxX do
+ let yf, xf = float y, float x
+ if e.Contains xf yf && neighbors |> List.forall (fun (otherE, _, _) -> not <| otherE.Contains xf yf)
+ then
+ totalElement <- totalElement + 1
+ if fgData.[y, x, 0] > 0uy
+ then
+ fgElement <- fgElement + 1
+
+ if totalElement < config.minimumCellArea || (float fgElement) / (float totalElement) < config.percentageOfFgValidCell
+ then
+ e.Removed <- true
+
+ // 3) Define pixels associated to each ellipse and create the cells.
+
+ // 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) =
+ e.Contains p.X p.Y &&
+ seq {
+ let c = PointD(e.Cx, e.Cy)
+ for d1 in lines do
+ let d2 = Utils.lineFromTwoPoints c p
+ 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.
+ } |> Seq.forall id
+
+
+ ellipsesWithNeigbors
+ |> List.choose (fun (e, neighbors) ->
+ if e.Removed
+ then
+ None
+ else
+ let minX, minY, maxX, maxY = ellipseWindow e
+
+ let infectedPixels = List<Point>()
+ let mutable stainPixels = 0
+ let mutable darkStainPixels = 0
+ let mutable nbElement = 0
+
+ let mutable stain_x = 0.0
+ let mutable stain_x2 = 0.0
+ let mutable stain_y = 0.0
+ let mutable stain_y2 = 0.0
+
+ let mutable sumCoords_y = 0
+ let mutable sumCoords_x = 0
+
+ let elements = new Matrix<byte>(maxY - minY + 1, maxX - minX + 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))) &&
+ fg.Data.[y, x, 0] > 0uy
+ then
+ elements.[y-minY, x-minX] <- 1uy
+ nbElement <- nbElement + 1
+ sumCoords_y <- sumCoords_y + y
+ sumCoords_x <- sumCoords_x + x
+
+ if infection.Data.[y, x, 0] > 0uy
+ then
+ infectedPixels.Add(Point(x, y))
+
+ if parasites.stain.Data.[y, x, 0] > 0uy
+ then
+ stainPixels <- stainPixels + 1
+ stain_x <- stain_x + (float x)
+ stain_x2 <- stain_x2 + (float x) ** 2.0
+ stain_y <- stain_y + (float y)
+ stain_y2 <- stain_y2 + (float y) ** 2.0
+
+ if parasites.darkStain.Data.[y, x, 0] > 0uy
+ then
+ darkStainPixels <- darkStainPixels + 1
+
+ let cellClass =
+ if float darkStainPixels > config.MaxDarkStainRatio * (float nbElement) (* ||
+ sqrt (((float sumCoords_x) / (float nbElement) - e.Cx) ** 2.0 + ((float sumCoords_y) / (float nbElement) - e.Cy) ** 2.0) > e.A * config.maxOffcenter *)
+ then
+ Peculiar
+ elif infectedPixels.Count > config.infectionPixelsRequired
+ then
+ let infectionToRemove = ImgTools.connectedComponents parasites.stain infectedPixels
+ for p in infectionToRemove do
+ infection.Data.[p.Y, p.X, 0] <- 0uy
+ InfectedRBC
+ else
+ HealthyRBC
+
+ Some { cellClass = cellClass
+ center = Point(roundInt e.Cx, roundInt e.Cy)
+ elements = elements })
module Config
+type Debug =
+ | DebugOff
+ | DebugOn of string // Output directory.
+
type Config = {
+ debug: Debug
+
scale: float
+ minRBCSize: float
+ maxRBCSize: float
+
doGSigma1: float
doGSigma2: float
doGLowFreqPercentageReduction: float
stainSigma: float
stainLevel: float
stainSpreadRequired: float
-
+
infectionSigma: float
infectionLevel: float
infectionPixelsRequired: int
+
+ percentageOfFgValidCell: float
+
+ MaxDarkStainRatio: float
+
+ minimumCellArea: int
+ maxOffcenter: float
}
\ No newline at end of file
| _ -> -1.0
if nintpts = 0
then Some (area, [||], [||])
- else Some (area, xint.[..nintpts-1], yint.[..nintpts-1])
\ No newline at end of file
+ else
+ let xTransform = Array.zeroCreate nintpts
+ let yTransform = Array.zeroCreate nintpts
+ for i in 0 .. (nintpts - 1) do
+ xTransform.[i] <- cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
+ yTransform.[i] <- sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
+ Some (area, xTransform, yTransform)
\ No newline at end of file
let x = (b + a) / 2.0
x, f x
+// Ellipse.A is always equal or greater than Ellipse.B.
+// Ellipse.Alpha is between 0 and Pi.
let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option =
let accuracy_extremum_search_1 = 4
let accuracy_extremum_search_2 = 3
let alpha1 = atan m1
let alpha2 = atan m2
- let r1 = sqrt (p1x**2.0 + p1y**2.0)
+ let r1 = sqrt (p1x ** 2.0 + p1y ** 2.0)
let theta1 = atan2 p1y p1x
let r2 = sqrt (p2x**2.0 + p2y**2.0)
let valid =
4.0 * sin (alpha1 - theta1) * (-r1 * sin (alpha1 - theta1) + r2 * sin (alpha1 - theta2)) *
sin (alpha2 - theta2) * (-r1 * sin (alpha2 - theta1) + r2 * sin (alpha2 - theta2)) +
- r1 * r2 * sin (alpha1 - alpha2) **2.0 * sin (theta1 - theta2) **2.0 < 0.0
+ r1 * r2 * sin (alpha1 - alpha2) ** 2.0 * sin (theta1 - theta2) ** 2.0 < 0.0
if valid
then
let r theta =
(r1 * r2 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0)
let rabs = r >> abs
let rellipse theta =
sqrt (
- rc**2.0 + (r1**2.0 * r2**2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2))**2.0 * sin (theta1 - theta2)**2.0) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)**2.0 -
+ rc ** 2.0 + (r1 ** 2.0 * r2 ** 2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) ** 2.0 * sin (theta1 - theta2) ** 2.0) /
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0) ** 2.0 -
(2.0 * r1 * r2 * rc * cos (theta - psi) * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) /
- (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0))
+ (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.0))
// We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval.
let r1eTheta, r1e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Maximum // Pi/2 and not pi because the period is Pi.
else
Some (m1, m2)
-
let find (edges: Matrix<byte>)
(xDir: Image<Gray, float>)
(yDir: Image<Gray, float>)
+++ /dev/null
-module ImageAnalysis
-
-open System
-open System.Drawing
-
-open Emgu.CV
-open Emgu.CV.Structure
-
-open Utils
-open ImgTools
-open Config
-open Types
-
-(*type Result = {
- RBCPositions : Point list
- infectedRBCPositions : Point list
- img: Image<Bgr, byte>
-}*)
-
-let doAnalysis (img: Image<Bgr, byte>) (config: Config) : Classifier.Cell list =
-
- let imgFloat = img.Convert<Bgr, float32>()
- use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area)
-
- (*use scaledImg =
- if config.scale = 1.0
- then
- img
- else
- let m = new Mat()
- CvInvoke.Resize(img, m, Size(roundInt (float img.Size.Width * config.scale), roundInt (float img.Size.Height * config.scale)))
- m*)
-
- use green = scaledImg.Item(1)
-
- //use green = new Matrix<byte>(scaledImg.Size)
- //CvInvoke.MixChannels(scaledImg, green, [| 1; 0 |])
-
- //let greenMatrix = new Matrix<byte>(green.Height, green.Width, green.DataPointer)
-
- //let test = greenMatrix.[10, 10]
-
- use filteredGreen = (gaussianFilter green config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter green config.doGSigma2)
-
- 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 = filteredGreen.Convolution(sobelKernel).Convert<Gray, float>()
- use yEdges = filteredGreen.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 minLocation = ref <| Point()
- let max = ref 0.0
- let maxLocation = ref <| Point()
- magnitudes.MinMax(min, max, minLocation, maxLocation)
-
- use magnitudesByte = ((magnitudes / !max) * 255.0).Convert<byte>() // Otsu from OpenCV only support 'byte'.
- 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 "Removing small connected components" (fun () ->
- removeArea edges 12)
-
- saveMat (edges * 255.0) "edges.png"
-
-
- let kmediansResults = KMedians.kmedians filteredGreen 1.0
-
- let parasites = ParasitesMarker.find green filteredGreen kmediansResults config
-
- saveImg parasites.darkStain "parasites_dark_stain.png"
- saveImg parasites.stain "parasites_stain.png"
- saveImg parasites.infection "parasites_infection.png"
-
- let radiusRange = config.scale * 20.0, config.scale * 40.0
- let windowSize = roundInt (1.6 * (snd radiusRange))
- let factorNbPick = 1.5
- let ellipses = logTime "Finding ellipses" (fun () ->
- Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick) |> List.filter (fun e -> not (e.CutAVericalLine 0.0) &&
- not (e.CutAVericalLine (float img.Width)) &&
- not (e.CutAnHorizontalLine 0.0) &&
- not (e.CutAnHorizontalLine (float img.Height)))
-
- drawEllipses img ellipses (Bgr(0.0, 255.0, 255.0))
- //saveImg img "ellipses.png"
-
- Classifier.findCells ellipses parasites kmediansResults.fg
-
- //
-
- (*use imageHSV = scaledImage.Convert<Hsv, uint8>()
- let H, S = match imageHSV.Split() with // Warning: H is from 0 to 179°.
- | [| H; S; _|] -> H, S
- | _ -> failwith "unable to split the HSV channels"
-
- let hueShiftValue = 175
- // Modulo operator doesn't exist on matrix thus we have to apply a function to every pixels.
- let correctedH : Image<Gray, byte> = H.Convert(fun b _ _ ->
- (255 - int(b) * 255 / 179 + hueShiftValue) % 256 |> byte
- )
-
- let correctedS : Image<Gray, byte> = S.Not()
-
- let filteredH = correctedH.SmoothMedian(5)
- let filteredS = correctedS.SmoothMedian(5)*)
-
- //let greenChannel = scaledImage.Item(1)
-
- //let filteredImage = (gaussianFilter greenChannel config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter greenChannel config.doGSigma2)
-
- // let filteredImage = greenChannel.ThresholdAdaptive(Gray(255.), CvEnum.AdaptiveThresholdType.GaussianC, CvEnum.ThresholdType.Binary, 61, Gray(5.0))
- // let thresholdedImage = filteredImage.CopyBlank()
-
- // CvInvoke.Threshold(filteredImage, thresholdedImage, 0., 255., CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.BinaryInv) |> ignore
-
- // filteredImage <|
\ No newline at end of file
open Utils
// Normalize image values between 0uy and 255uy.
-let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
+let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
let min = ref [| 0.0 |]
let minLocation = ref <| [| Point() |]
let max = ref [| 0.0 |]
let size = 2 * int (ceil (4.0 * standardDeviation)) + 1
img.SmoothGaussian(size, size, standardDeviation, standardDeviation)
-// Zhang and Suen algorithm.
+// Zhang and Suen algorithm.
// Modify 'mat' in place.
-let thin (mat: Matrix<byte>) =
- let neighbors = [|
+let thin (mat: Matrix<byte>) =
+ let neighbors = [|
(-1, 0) // p2
(-1, 1) // p3
( 0, 1) // p4
let mutable data2 = Array2D.zeroCreate<byte> h w
// Return the list of neighbor values.
- let neighborsValues (p1i, p1j) =
- Array.map (fun (ni, nj) ->
+ 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 pattern01 (values: byte[]) =
let mutable nb = 0
let mutable lastValue = 255uy
for v in values do
then
nb <- nb + 1
nb
-
+
let mutable pixelChanged = true
let mutable oddIteration = true
while pixelChanged do
else
data2.[i, j] <- 0uy
- oddIteration <- not oddIteration
+ oddIteration <- not oddIteration
let tmp = data1
data1 <- data2
data2 <- tmp
-
+
+
+
+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<byte>) (areaSize: int) =
- let neighbors = [|
+ let neighbors = [|
(-1, 0) // p2
(-1, 1) // p3
( 0, 1) // p4
( 1, -1) // p7
( 0, -1) // p8
(-1, -1) |] // p9
-
+
let mat' = new Matrix<byte>(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))
for (ni, nj) in neighborhood do
data.[ni, nj] <- 0uy
+let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
+ let w = img.Width
+ let h = img.Height
+
+ let pointChecked = HashSet<Point>()
+ let pointToCheck = List<Point>(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<Point>(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<Gray, 'TDeph>(mat.Size)
mat.CopyTo(img)
- saveImg img name
-
+ saveImg img filepath
+
(*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)
-
+let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) =
img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1);
-let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) =
+let drawLineF (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)
+ drawLine img color x0 y0 x1 y1
+
+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
if not first_iteration
then
- drawLine img color x0 y0 x y
+ drawLineF img color x0 y0 x y
else
first_iteration <- false
y0 <- y
let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Types.Ellipse list) (color: 'TColor) =
- List.iter (fun e -> drawEllipse img e color) ellipses
\ No newline at end of file
+ List.iter (fun e -> drawEllipse img e color) ellipses
+
+
+let rngCell = System.Random()
+let drawCell (img: Image<Bgr, byte>) (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 = match c.cellClass with
+ | Types.HealthyRBC -> Bgr(255.0, 0.0, 0.0)
+ | Types.InfectedRBC -> Bgr(0.0, 0.0, 255.0)
+ | Types.Peculiar -> Bgr(0.0, 0.0, 0.0)
+
+ drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y
+ drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3)
+
+let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Types.Cell list) =
+ List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file
let max = ref [| 0.0 |]
let maxLocation = ref <| [| Point() |]
img.MinMax(min, max, minLocation, maxLocation)
-
+
let mutable median_bg = (!max).[0] - ((!max).[0] - (!min).[0]) / 4.0
let mutable median_fg = (!min).[0] + ((!max).[0] - (!min).[0]) / 4.0
let mutable d_bg = new Image<Gray, float32>(img.Size)
let mutable d_fg = new Image<Gray, float32>(img.Size)
let mutable fg = new Image<Gray, byte>(img.Size)
-
+
for i in 1..nbIteration do
CvInvoke.Pow(img - median_bg, 2.0, d_bg)
CvInvoke.Pow(img - median_fg, 2.0, d_fg)
fg <- (d_fg * fgFactor).Cmp(d_bg, CvEnum.CmpType.LessThan)
-
+
median_fg <- MathNet.Numerics.Statistics.Statistics.Median(seq {
for i in 0 .. h - 1 do
for j in 0 .. w - 1 do
if fg.Data.[i, j, 0] > 0uy then yield img.Data.[i, j, 0] |> float })
-
+
median_bg <- MathNet.Numerics.Statistics.Statistics.Median(seq {
for i in 0 .. h - 1 do
for j in 0 .. w - 1 do
if fg.Data.[i, j, 0] = 0uy then yield img.Data.[i, j, 0] |> float })
-
+
CvInvoke.Sqrt(d_fg, d_fg)
{ fg = fg; median_bg = median_bg; median_fg = median_fg; d_fg = d_fg }
| Node of float * Tree<'a> * Tree<'a>
| Leaf of 'a
- static member buildTree (l: 'a list) : Tree<'a> =
+ static member BuildTree (l: 'a list) : Tree<'a> =
let xSorted = List.toArray l
let ySorted = List.toArray l
Array.sortInPlaceWith cmpX xSorted
buildTreeFromSortedArray xSorted ySorted 1
- static member search (tree: Tree<'a>) (searchRegion: Region) : 'a list =
+ member this.Search (searchRegion: Region) : 'a list =
let rec valuesFrom (tree: Tree<'a>) : 'a list =
match tree with
| Leaf v -> [v]
let upRegion = { currentRegion with minY = splitValue }
(valuesInRegion downRegion part1) @ (valuesInRegion upRegion part2)
- searchWithRegion tree { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1
+ searchWithRegion this { minX = Double.MinValue; maxX = Double.MaxValue; minY = Double.MinValue; maxY = Double.MaxValue } 1
///// Tests. TODO: to put in a unit test.
Point(3.5, 1.5)
Point(2.5, 0.5) ]
- let tree = Tree.buildTree pts
+ let tree = Tree.BuildTree pts
Utils.dprintfn "Tree: %A" tree
- let s1 = Tree.search tree { minX = 0.0; maxX = 5.0; minY = 0.0; maxY = 5.0 } // All points.
+ let s1 = tree.Search { minX = 0.0; maxX = 5.0; minY = 0.0; maxY = 5.0 } // All points.
Utils.dprintfn "s1: %A" s1
- let s2 = Tree.search tree { minX = 2.8; maxX = 4.5; minY = 3.0; maxY = 4.5 }
+ let s2 = tree.Search { minX = 2.8; maxX = 4.5; minY = 3.0; maxY = 4.5 }
Utils.dprintfn "s2: %A" s2
- let s3 = Tree.search tree { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 }
+ let s3 = tree.Search { minX = 2.0; maxX = 2.0; minY = 2.0; maxY = 2.0 }
Utils.dprintfn "s3: %A" s3
let test2 () =
Point(1.0, 2.0)
Point(1.0, 3.0) ]
- let tree = Tree.buildTree pts
+ let tree = Tree.BuildTree pts
Utils.dprintfn "Tree: %A" tree
- let s1 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 1.0; maxY = 1.0 }
+ let s1 = tree.Search { minX = 1.0; maxX = 1.0; minY = 1.0; maxY = 1.0 }
Utils.dprintfn "s1: %A" s1
- let s2 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 2.0; maxY = 2.0 }
+ let s2 = tree.Search { minX = 1.0; maxX = 1.0; minY = 2.0; maxY = 2.0 }
Utils.dprintfn "s2: %A" s2
// This case result is wrong: FIXME
- let s3 = Tree.search tree { minX = 1.0; maxX = 1.0; minY = 3.0; maxY = 3.0 }
+ let s3 = tree.Search { minX = 1.0; maxX = 1.0; minY = 3.0; maxY = 3.0 }
Utils.dprintfn "s3: %A" s3
- let s4 = Tree.search tree { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 }
+ let s4 = tree.Search { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 }
Utils.dprintfn "s4: %A" s4
--- /dev/null
+module ImageAnalysis
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open Utils
+open ImgTools
+open Config
+open Types
+
+
+let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) : Cell list =
+
+ let imgFloat = img.Convert<Bgr, float32>()
+ use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area)
+
+ use green = scaledImg.Item(1)
+
+ use filteredGreen = (gaussianFilter green config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter green config.doGSigma2)
+
+ 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 = filteredGreen.Convolution(sobelKernel).Convert<Gray, float>()
+ use yEdges = filteredGreen.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 minLocation = ref <| Point()
+ let max = ref 0.0
+ let maxLocation = ref <| Point()
+ magnitudes.MinMax(min, max, minLocation, maxLocation)
+
+ use magnitudesByte = ((magnitudes / !max) * 255.0).Convert<byte>() // Otsu from OpenCV only support 'byte'.
+ 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 "Removing small connected components from thinning" (fun () -> removeArea edges 12)
+
+ let kmediansResults = KMedians.kmedians filteredGreen 1.0
+
+ let parasites = ParasitesMarker.find green filteredGreen kmediansResults config
+
+ let radiusRange = config.scale * config.minRBCSize, config.scale * config.maxRBCSize
+ let windowSize = roundInt (1.6 * (snd radiusRange))
+ let factorNbPick = 1.5
+ let ellipses = logTime "Finding ellipses" (fun () ->
+ Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick)
+
+ let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites kmediansResults.fg config)
+
+ match config.debug with
+ | DebugOn output ->
+ let buildFileName postfix = System.IO.Path.Combine(output, name + postfix)
+ saveMat (edges * 255.0) (buildFileName " - edges.png")
+ saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
+ saveImg parasites.stain (buildFileName " - parasites - stain.png")
+ saveImg parasites.infection (buildFileName " - parasites - infection.png")
+
+ let imgEllipses = img.Copy()
+ drawEllipses imgEllipses ellipses (Bgr(0.0, 240.0, 240.0))
+ saveImg imgEllipses (buildFileName " - ellipses.png")
+
+ saveImg (kmediansResults.fg * 255.0) (buildFileName " - foreground.png")
+
+ let imgCells = img.Copy()
+ drawCells imgCells false cells
+ saveImg imgCells (buildFileName " - cells.png")
+
+ let imgCells' = img.Copy()
+ drawCells imgCells' true cells
+ saveImg imgCells' (buildFileName " - cells - full.png")
+ | _ -> ()
+
+ cells
let matchingScoreThreshold1 = 0.6
let matchingScoreThreshold2 = 1.0
-type EllipseScore (matchingScore: float, e: Ellipse) =
+type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) =
let mutable matchingScore = matchingScore
- member this.MatchingScore = matchingScore
member this.Ellipse = e
- member val Processed = false with get, set
- member val Removed = false with get, set
+ member this.MatchingScore = matchingScore
member this.AddMatchingScore(score: float) =
matchingScore <- matchingScore + score
-
+
+ member val Processed = false with get, set
+ member val Removed = false with get, set
+
interface KdTree.I2DCoords with
- member this.X = e.Cx
- member this.Y = e.Cy
+ member this.X = this.Ellipse.Cx
+ member this.Y = this.Ellipse.Cy
+
type MatchingEllipses (radiusMin: float) =
- let ellipses = List<EllipseScore>()
+ let ellipses = List<EllipseScoreFlaggedKd>()
member this.Add (e: Ellipse) =
- ellipses.Add(EllipseScore(0.0, e))
+ ellipses.Add(EllipseScoreFlaggedKd(0.0, e))
member this.Ellipses : Ellipse list =
- // 1) Create a kd-tree from the ellipses list.
- let tree = KdTree.Tree.buildTree (List.ofSeq ellipses)
-
- // 2) Compute the matching score of each ellipses.
- let windowSize = radiusMin
- for e in ellipses do
- e.Processed <- true
- let areaE = e.Ellipse.Area
- let window = { KdTree.minX = e.Ellipse.Cx - windowSize / 2.0
- KdTree.maxX = e.Ellipse.Cx + windowSize / 2.0
- KdTree.minY = e.Ellipse.Cy - windowSize / 2.0
- KdTree.maxY = e.Ellipse.Cy + windowSize / 2.0 }
- for other in KdTree.Tree.search tree window do
- if not other.Processed
- then
- let areaOther = other.Ellipse.Area
- match EEOver.EEOverlapArea e.Ellipse other.Ellipse with
- | Some (commonArea, _, _) ->
- let matchingScore = 2.0 * commonArea / (areaE + areaOther)
- if matchingScore >= matchingScoreThreshold1
- then
- other.AddMatchingScore(matchingScore)
- e.AddMatchingScore(matchingScore)
- | _ -> ()
-
- // 3) Sort ellipses by their score.
- ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore))
-
- // 4) Remove ellipses wich have a low score.
- let i = ellipses.BinarySearch(EllipseScore(matchingScoreThreshold2, Ellipse(0.0, 0.0, 0.0, 0.0, 0.0)),
- { new IComparer<EllipseScore> with
- member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs
- let nbToRemove = ellipses.Count - i
- if nbToRemove > 0
+ if ellipses.Count = 0
then
- for j in i .. nbToRemove - 1 do
- ellipses.[j].Removed <- true
- ellipses.RemoveRange(i, nbToRemove)
+ []
+ else
+ // 1) Create a kd-tree from the ellipses list.
+ let tree = KdTree.Tree.BuildTree (List.ofSeq ellipses)
+
+ // 2) Compute the matching score of each ellipses.
+ let windowSize = radiusMin
+ for e in ellipses do
+ e.Processed <- true
+ let areaE = e.Ellipse.Area
+ let window = { KdTree.minX = e.Ellipse.Cx - windowSize / 2.0
+ KdTree.maxX = e.Ellipse.Cx + windowSize / 2.0
+ KdTree.minY = e.Ellipse.Cy - windowSize / 2.0
+ KdTree.maxY = e.Ellipse.Cy + windowSize / 2.0 }
+ for other in tree.Search window do
+ if not other.Processed
+ then
+ let areaOther = other.Ellipse.Area
+ match EEOver.EEOverlapArea e.Ellipse other.Ellipse with
+ | Some (commonArea, _, _) ->
+ let matchingScore = 2.0 * commonArea / (areaE + areaOther)
+ if matchingScore >= matchingScoreThreshold1
+ then
+ other.AddMatchingScore(matchingScore)
+ e.AddMatchingScore(matchingScore)
+ | _ -> ()
- // 5) Remove ellipses whose center is into an ellipse with a better score
- for e in ellipses do
- if not e.Removed
+ // 3) Sort ellipses by their score.
+ ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore))
+
+ // 4) Remove ellipses wich have 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
+ let nbToRemove = ellipses.Count - i
+ if nbToRemove > 0
then
- let window = { KdTree.minX = e.Ellipse.Cx - e.Ellipse.A
- KdTree.maxX = e.Ellipse.Cx + e.Ellipse.A
- KdTree.minY = e.Ellipse.Cy - e.Ellipse.A
- KdTree.maxY = e.Ellipse.Cy + e.Ellipse.A }
- for other in KdTree.Tree.search tree window do
- if not other.Removed && other.MatchingScore < e.MatchingScore
- then
- if e.Ellipse.Contains other.Ellipse.Cx other.Ellipse.Cy
- then
- other.Removed <- true
- ellipses.RemoveAll(fun e -> e.Removed) |> ignore
-
- List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
+ for j in i .. nbToRemove - 1 do
+ ellipses.[j].Removed <- true
+ ellipses.RemoveRange(i, nbToRemove)
+
+ // 5) Remove ellipses whose center is into an ellipse with a better score
+ for e in ellipses do
+ if not e.Removed
+ then
+ let window = { KdTree.minX = e.Ellipse.Cx - e.Ellipse.A
+ KdTree.maxX = e.Ellipse.Cx + e.Ellipse.A
+ KdTree.minY = e.Ellipse.Cy - e.Ellipse.A
+ KdTree.maxY = e.Ellipse.Cy + e.Ellipse.A }
+ for other in tree.Search window do
+ if not other.Removed && other.MatchingScore < e.MatchingScore
+ then
+ if e.Ellipse.Contains other.Ellipse.Cx other.Ellipse.Cy
+ then
+ other.Removed <- true
+ ellipses.RemoveAll(fun e -> e.Removed) |> ignore
+
+ List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
-
-
<PlatformTarget>x64</PlatformTarget>
<DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
<Prefer32Bit>false</Prefer32Bit>
+ <StartArguments>--folder "../../../../Images/08.10.2015/test" --output "../../../Images/output" --debug</StartArguments>
</PropertyGroup>
<PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
<DebugType>pdbonly</DebugType>
<PlatformTarget>x64</PlatformTarget>
<DocumentationFile>bin\Release\Parasitemia.XML</DocumentationFile>
<Prefer32Bit>false</Prefer32Bit>
+ <StartArguments>--folder "../../../../Images/08.10.2015/Without scale" --output "../../../Images/output" --debug</StartArguments>
</PropertyGroup>
<PropertyGroup>
<MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
<Compile Include="NumericUpDown.xaml.fs" />
<Resource Include="MainWindow.xaml" />
<Compile Include="MainWindow.xaml.fs" />
- <Compile Include="Utils.fs" />
<Compile Include="Types.fs" />
+ <Compile Include="Utils.fs" />
<Compile Include="ImgTools.fs" />
<Compile Include="Config.fs" />
<Compile Include="KMedians.fs" />
<Compile Include="EEOver.fs" />
<Compile Include="ParasitesMarker.fs" />
<Compile Include="KdTree.fs" />
- <Compile Include="Classifier.fs" />
<Compile Include="MatchingEllipses.fs" />
+ <Compile Include="Classifier.fs" />
<Compile Include="Ellipse.fs" />
- <Compile Include="ImageAnalysis.fs" />
+ <Compile Include="MainAnalysis.fs" />
<Compile Include="Program.fs" />
<None Include="App.config" />
<Content Include="packages.config" />
stain: Image<Gray, byte>
infection: Image<Gray, byte> }
+// Create three binary markers :
+// * 'Dark stain' corresponds to the colored pixel, it's independent of the size of the areas.
+// * 'Stain' corresponds to the stain around the parasites.
+// * 'Infection' corresponds to the parasite. It shouldn't contain thrombocytes.
let find (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (kmediansResult: KMedians.Result) (config: Config.Config) : Result =
+ let green = ImgTools.gaussianFilter green 1.0
+
// We use the filtered image to find the dark stain.
let { KMedians.fg = fg; KMedians.median_bg = median_bg; KMedians.median_fg = median_fg; KMedians.d_fg = d_fg } = kmediansResult
let darkStain = d_fg.Cmp(median_bg * config.darkStainLevel, CvEnum.CmpType.GreaterThan)
darkStain._And(filteredGreen.Cmp(median_fg, CvEnum.CmpType.LessThan))
darkStain._And(fg)
-
+
let fgFloat = (fg / 255.0).Convert<Gray, float32>()
let greenWithoutBg = green.Copy()
greenWithoutBg.SetValue(Gray(0.0), fg.Not())
- let findSmears (sigma: float) (level: float) : Image<Gray, byte> =
+ let findSmears (sigma: float) (level: float) : Image<Gray, byte> =
let greenWithoutBgSmoothed = ImgTools.gaussianFilter greenWithoutBg sigma
let fgSmoothed = ImgTools.gaussianFilter fgFloat sigma
let smears = (greenWithoutBg.Mul(fgSmoothed)).Cmp(greenWithoutBgSmoothed.Mul(level), CvEnum.CmpType.LessThan)
smears._And(fg)
- smears
+ smears
{ darkStain = darkStain;
stain = findSmears config.stainSigma config.stainLevel
infection = findSmears config.infectionSigma config.infectionLevel }
-
+
module Parasitemia.Main
+open System
open System.IO
open System.Windows
open System.Windows.Media
open Config
-let display (window : Views.MainWindow) (img : IImage) =
+let display (window : Views.MainWindow) (img : IImage) =
let imgControl = window.Root.FindName("img") :?> Controls.Image
imgControl.Source <- BitmapSourceConvert.ToBitmapSource(img)
let log (window : Views.MainWindow) (mess : string) =
let txtLog = window.Root.FindName("txtLog") :?> Controls.TextBlock
txtLog.Text <- txtLog.Text + mess + "\n"
-
-[<System.STAThreadAttribute>]
-do
- Utils.dprintfn "OpenCV test"
-
- let app = new Application()
- let mainWindow = Views.MainWindow()
-
- Utils.log <- (fun m -> log mainWindow m)
-
- let config = {
- scale = 1.0
-
- doGSigma1 = 1.5
- doGSigma2 = 20.0
- doGLowFreqPercentageReduction = 0.9
-
- darkStainLevel = 0.839
-
- stainSigma = 10.0
- stainLevel = 0.9
- stainSpreadRequired = 3.0
-
- infectionSigma = 2.0
- infectionLevel = 0.762
- infectionPixelsRequired = 1 }
-
- ///// ELLIPSES /////
- //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png")
- //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single.png")
- //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single_oblong_4.png")
- //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/strange_rbc_1.png")
- //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single_blurred.png")
- use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/lot.png")
-
-
- ///// PARASITES /////
- //use img = new Image<Bgr, byte>("../../../../src/Parasites/images/1.png")
-
-// KdTree.test3 ()
-// Utils.dprintfn "area: %A" area
-
- let result = Utils.logTime "Total" (fun () ->
- ImageAnalysis.doAnalysis img config )
-
- display mainWindow img
- mainWindow.Root.Show()
- app.Run() |> ignore
+
+
+type Input =
+ | File of string
+ | Dir of string
+
+type RunningMode =
+ | CmdLine of Input * string
+ | Window
+
+type Arguments = RunningMode * bool
+
+let parseArgs (args: string[]) : Arguments =
+
+ let output = Array.tryFindIndex ((=) "--output") args
+
+ let runningMode =
+ match Array.tryFindIndex ((=) "--folder") args, output with
+ | Some i, Some i_output when i < args.Length - 2 && i_output < args.Length - 2 ->
+ CmdLine ((Dir args.[i+1]), args.[i_output + 1])
+ | _ ->
+ match Array.tryFindIndex ((=) "--file") args, output with
+ | Some i, Some i_output when i < args.Length - 2 && i_output < args.Length - 2 ->
+ CmdLine ((File args.[i+1]), args.[i_output + 1])
+ |_ ->
+ Window
+
+ runningMode, Array.exists ((=) "--debug") args
+
+
+[<EntryPoint>]
+let main args =
+ match parseArgs args with
+ | mode, debug ->
+ let scale = 1.
+ let config = {
+ debug = if debug then DebugOn "." else DebugOff
+
+ scale = scale
+
+ // RBC size range in px at scale = 1.0.
+ minRBCSize = 20.
+ maxRBCSize = 40.
+
+ doGSigma1 = 1.5
+ doGSigma2 = 20.
+ doGLowFreqPercentageReduction = 0.75
+
+ darkStainLevel = 0.5
+
+ stainSigma = 10.
+ stainLevel = 0.9
+ stainSpreadRequired = 3.0
+
+ infectionSigma = 2.2
+ infectionLevel = 0.85
+ infectionPixelsRequired = 1
+
+ percentageOfFgValidCell = 0.4
+
+ MaxDarkStainRatio = 0.1
+
+ minimumCellArea = 600. * scale ** 2. |> int
+ maxOffcenter = 0.5 }
+
+ match mode with
+ | CmdLine (input, output) ->
+ let config = { config with debug = DebugOn output }
+ Utils.log <- (fun m -> Console.WriteLine m)
+
+ Directory.CreateDirectory output |> ignore
+
+ let files = match input with
+ | File file -> [ file ]
+ | Dir dir -> Directory.EnumerateFiles dir |> List.ofSeq
+
+ use resultFile = new StreamWriter(new FileStream(Path.Combine(output, "results.txt"), FileMode.Create, FileAccess.Write))
+
+ for file in files do
+ try
+ let fileInfo = FileInfo(file)
+ use img = new Image<Bgr, byte>(file)
+ let cells = Utils.logTime "Whole analyze" (fun () -> ImageAnalysis.doAnalysis img fileInfo.Name config)
+ let total, infected = Utils.countCells cells
+ fprintf resultFile "File: %s %d %d %.2f\n" file total infected (100. * (float infected) / (float total))
+ with
+ | _ as ex -> Utils.log (sprintf "Unable to open the image '%A': %A" file ex)
+ 0
+
+ | Window ->
+ let app = new Application()
+ let mainWindow = Views.MainWindow()
+
+ Utils.log <- (fun m -> log mainWindow m)
+
+ //display mainWindow img
+ mainWindow.Root.Show()
+ app.Run()
module Types
open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) =
member this.Cx = cx
member this.Area = a * b * Math.PI
// Does the ellipse contain the point (x, y)?.
- member this.Contains x y =
+ member this.Contains x y =
((x - cx) * cos alpha + (y - cy) * sin alpha) ** 2.0 / a ** 2.0 + ((x - cx) * sin alpha - (y - cy) * cos alpha) ** 2.0 / b ** 2.0 <= 1.0
// A line is defined as : y = mx + l
member this.CutALine (m: float) (l: float) : bool =
- -2.0 * l ** 2.0 + a ** 2.0 + m ** 2.0 * a ** 2.0 + b ** 2.0 + m ** 2.0 * b ** 2.0 -
- 4.0 * m * l * cx - 2.0 * m ** 2.0 * cx ** 2.0 + 4.0 * l * cy + 4.0 * m * cx * cy -
+ -2.0 * l ** 2.0 + a ** 2.0 + m ** 2.0 * a ** 2.0 + b ** 2.0 + m ** 2.0 * b ** 2.0 -
+ 4.0 * m * l * cx - 2.0 * m ** 2.0 * cx ** 2.0 + 4.0 * l * cy + 4.0 * m * cx * cy -
2.0 * cy ** 2.0 + (-1.0 + m ** 2.0) * (a ** 2.0 - b ** 2.0) * cos (2.0 * alpha) - 2.0 * m * (a ** 2.0 - b ** 2.0) * sin (2.0 * alpha) > 0.0
- member this.CutAVericalLine (y: float) : bool =
+ member this.CutAVericalLine (y: float) : bool =
a ** 2.0 + b ** 2.0 - 2.0 * y ** 2.0 + 4.0 * y * cx - 2.0 * cx ** 2.0 + a ** 2.0 * cos (2.0 * alpha) - b ** 2.0 * cos (2.0 * alpha) > 0.0
- member this.CutAnHorizontalLine (x: float) : bool =
- a ** 2.0 + b ** 2.0 - 2.0 * x ** 2.0 + 4.0 * x * cy - 2.0 * cy ** 2.0 - a ** 2.0 * cos (2.0 * alpha) + b ** 2.0 * cos (2.0 * alpha) > 0.0
\ No newline at end of file
+ member this.CutAnHorizontalLine (x: float) : bool =
+ a ** 2.0 + b ** 2.0 - 2.0 * x ** 2.0 + 4.0 * x * cy - 2.0 * cy ** 2.0 - a ** 2.0 * cos (2.0 * alpha) + b ** 2.0 * cos (2.0 * alpha) > 0.0
+
+ member this.isOutside (width: float) (height: float) =
+ this.Cx <= 0.0 || this.Cx >= width ||
+ this.Cy <= 0.0 || this.Cy >= height ||
+ this.CutAVericalLine 0.0 || this.CutAVericalLine width ||
+ this.CutAnHorizontalLine 0.0 || this.CutAnHorizontalLine height
+
+
+type CellClass = HealthyRBC | InfectedRBC | Peculiar
+
+type Cell = {
+ cellClass: CellClass
+ center: Point
+ elements: Matrix<byte> }
+
+
+type Line (a: float, b: float) =
+ member this.A = a
+ member this.B = b
+
+type PointD (x: float, y: float) =
+ member this.X = x
+ member this.Y = y
+
open System.Diagnostics
+open Types
+
let roundInt = int << round
let inline dprintfn fmt =
let res = f ()
sw.Stop()
log <| sprintf "%A (time: %A ms)" m sw.ElapsedMilliseconds
- res
\ No newline at end of file
+ res
+
+let inline lineFromTwoPoints (p1: PointD) (p2: PointD) : Line =
+ let a = (p1.Y - p2.Y) / (p1.X - p2.X)
+ let b = -(p2.X * p1.Y - p1.X * p2.Y) / (p1.X - p2.X)
+ Line(a, b)
+
+let inline pointFromTwoLines (l1: Line) (l2: Line) : PointD =
+ let x = -(l1.B - l2.B) / (l1.A - l2.A)
+ let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A)
+ PointD(x, y)
+
+let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) =
+ (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0
+
+let distanceTwoPoints (p1: PointD) (p2: PointD) =
+ squaredDistanceTwoPoints p1 p2 |> sqrt
+
+let countCells (cells: Cell list) : int * int =
+ cells |> List.fold (fun (total, infected) { cellClass = cellClass } ->
+ match cellClass with
+ | HealthyRBC -> (total + 1, infected)
+ | InfectedRBC -> (total + 1, infected + 1)
+ | Peculiar -> (total, infected)) (0, 0)
\ No newline at end of file
projects / master-thesis.git / commitdiff