<?xml version="1.0" encoding="utf-8"?>
<configuration>
<startup>
- <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.6" />
+ <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.6.1" />
</startup>
<runtime>
<assemblyBinding xmlns="urn:schemas-microsoft-com:asm.v1">
module Config
type Config = {
+ scale: float
+
doGSigma1: float
doGSigma2: float
doGLowFreqPercentageReduction: float
- scale: float
+ // Parasites detection.
+ darkStainLevel: float
+
+ stainSigma: float
+ stainLevel: float
+ stainSpreadRequired: float
+
+ infectionSigma: float
+ infectionLevel: float
+ infectionPixelsRequired: int
}
\ No newline at end of file
let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff
let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff
+#if DEBUG_LOG
+ printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1
+ printf "theta=%f\n" theta
+ printf "eps_Radian=%f\n" eps_radian
+ printf "(x1, y1)=(%f, %f)\n" x1 y1
+ printf "(x2, y2)=(%f, %f)\n" x2 y2
+ printf "test1=%f\n" test1
+ printf "test2=%f\n" test2
+#endif
+
if test1 * test2 > 0.0
then TANGENT_POINT
else INTERSECTION_POINT
if area1 < 0.0
then
+#if DEBUG_LOG
+ printf "TWO area1=%f\n" area1
+#endif
area1 <- area1 + a1 * b1
let cosphi = cos (phi_1 - phi_2)
let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
if area2 < 0.0
then
+#if DEBUG_LOG
+ printf "TWO area2=%f\n" area2
+#endif
area2 <- area2 + a2 * b2
area1 + area2
for i in 0..2 do
if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
then
- tanpts <- tanpts + 1;
- tanindex <- i;
-
+ tanpts <- tanpts + 1
+ tanindex <- i
+#if DEBUG_LOG
+ printf "tanindex=%d\n" tanindex
+#endif
+
if tanpts <> 1
then
-1.0
xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
+#if DEBUG_LOG
+ for k in 0..3 do
+ printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
+#endif
+
for j in 1..3 do
let tmp0 = theta.[j]
let tmp1 = xint.[j]
let tmp2 = yint.[j]
let mutable k = j - 1
+ let mutable k2 = 0
while k >= 0 do
if theta.[k] <= tmp0
then
+ k2 <- k + 1
k <- -1
else
theta.[k+1] <- theta.[k]
xint.[k+1] <- xint.[k]
yint.[k+1] <- yint.[k]
- k <- k - 1
+ k <- k - 1
+ k2 <- k + 1
+
+ theta.[k2] <- tmp0
+ xint.[k2] <- tmp1
+ yint.[k2] <- tmp2
- theta.[k+1] <- tmp0
- xint.[k+1] <- tmp1
- yint.[k+1] <- tmp2
+
+#if DEBUG_LOG
+ printf "AFTER sorting\n"
+ for k in 0..3 do
+ printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
+#endif
let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
if area5 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
+#endif
area5 <- area5 + area_2
if area4 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
+#endif
area4 <- area4 + area_2
if area3 < 0.0
then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
+#endif
area3 <- area3 + area_1
if area2 < 0.0
- then
+ then
+#if DEBUG_LOG
+ printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
+#endif
area2 <- area2 + area_1
+#if DEBUG_LOG
+ printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
+#endif
+
area1 + area2 + area3 + area4 + area5
quad ()
-open Types
-
-let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float =
- let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1
- let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2
+let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float =
+ let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha
+ let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha
if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
then
let sinphi = sin phi_1
(h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
+#if DEBUG_LOG
+ printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
+#endif
+
let cosphi = cos phi_2r
let cosphi2 = cosphi ** 2.0
let sinphi = sin phi_2r
a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
|]
+#if DEBUG_LOG
+ for i in 0..4 do
+ printf "cy[%d]=%f\n" i cy.[i]
+#endif
+
let py = Array.zeroCreate<float> 5
let r = Array2D.zeroCreate<float> 3 5
for i in 0 .. 3 do
py.[4-i] <- cy.[i] / cy.[4]
py.[0] <- 1.0
+#if DEBUG_LOG
+ for i in 0..4 do
+ printf "py[%d]=%f\n" i py.[i]
+#endif
biquadroots py r
4
else
0
+#if DEBUG_LOG
+ printf "nroots = %d\n" nroots
+#endif
+
let ychk = [|
for i in 1 .. nroots do
if abs r.[2, i] < EPS
then
yield r.[1, i] * b1
+#if DEBUG_LOG
+ printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
+#endif
|]
Array.sortInPlace ychk
+#if DEBUG_LOG
+ printf "nychk=%d\n" ychk.Length
+ for j in 0 .. ychk.Length - 1 do
+ printf "\t j=%d, ychk=%f\n" j ychk.[j]
+#endif
+
let nychk = Array.length ychk
let mutable nintpts = 0
let mutable i = 0
while returnValue = 0.0 && i < nychk do
+#if DEBUG_LOG
+ printf "------------->i=%d (nychk=%d)\n" i nychk
+#endif
+
if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
then
+#if DEBUG_LOG
+ printf "check intersecting points. nintps is %d" nintpts
+#endif
+
let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
let x2 = -x1
+#if DEBUG_LOG
+ printf "\tx1=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x1)= %f\n" x1 ychk.[i] a1 b1 (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff)
+ printf "\tx2=%f, y1=%f, A=%f. B=%f ---> ellipse2tr(x2) %f\n" x2 ychk.[i] a1 b1 (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff)
+#endif
+
if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
then
nintpts <- nintpts + 1
+#if DEBUG_LOG
+ printf "first if x1. acc nintps=%d\n" nintpts
+#endif
if nintpts > 4
then
returnValue <- -1.0
else
xint.[nintpts-1] <- x1
yint.[nintpts-1] <- ychk.[i]
+#if DEBUG_LOG
+ printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
+#endif
if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
then
nintpts <- nintpts + 1
+#if DEBUG_LOG
+ printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
+#endif
if nintpts > 4
then
returnValue <- -1.0
else
xint.[nintpts-1] <- x2
yint.[nintpts-1] <- ychk.[i]
+
+#if DEBUG_LOG
+ printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
+#endif
+
+#if DEBUG_LOG
+ else
+ printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
+#endif
i <- i + 1
match nintpts with
| 0 | 1 -> nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
| 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
- | TANGENT_POINT -> nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
- | INTERSECTION_POINT -> twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
+ | TANGENT_POINT ->
+#if DEBUG_LOG
+ printf "one point is tangent\n"
+#endif
+ nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
+
+ | INTERSECTION_POINT ->
+#if DEBUG_LOG
+ printf "check twointpts\n"
+#endif
+ twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
| _ -> -1.0
-
if fc < fd
then
- b <- d;
- d <- c;
- c <- b - gr * (b - a);
+ b <- d
+ d <- c
+ c <- b - gr * (b - a)
else
- a <- c;
- c <- d;
- d <- a + gr * (b - a);
+ a <- c
+ c <- d
+ d <- a + gr * (b - a)
let x = (b + a) / 2.0
x, f x
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 accuracy_extremum_search_1 = 4
+ let accuracy_extremum_search_2 = 3
// p3 as the referencial.
let p1x = p1x - p3x
let cx = cx + p3x
let cy = cy + p3y
- Some { cx = cx; cy = cy; a = r1e; b = r2e; alpha = alpha }
+ Some (Types.Ellipse(cx, cy, r1e, r2e, alpha))
else
None
(windowSize: int)
(factorNbPick: float) : Types.Ellipse list =
- let increment = windowSize / 4;
+ let increment = windowSize / 4
let r1, r2 = radiusRange
let radiusTolerance = (r2 - r1) * 0.2
- let minimumDistance = (r2 / 1.5) ** 2.0;
- let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0;
+ let minimumDistance = (r2 / 1.5) ** 2.0
+ let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0
let h = edges.Height
let w = edges.Width
let xDirData = xDir.Data
let yDirData = yDir.Data
- let rng = Random()
+ let rng = Random(42)
- let ellipses = MatchingEllipses ()
+ let ellipses = MatchingEllipses(r1)
for window_i in -windowSize + increment .. increment .. h - increment do
for window_j in -windowSize + increment .. increment .. w - increment do
match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0] with
| Some (m1, m2) ->
match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with
- | Some e when e.cx > 0.0 && e.cx < (float w) - 1.0 && e.cy > 0.0 && e.cy < (float h) - 1.0 &&
- e.a >= r1 - radiusTolerance && e.a <= r2 + radiusTolerance && e.b >= r1 - radiusTolerance && e.b <= r2 + radiusTolerance ->
+ | Some e when e.Cx > 0.0 && e.Cx < (float w) - 1.0 && e.Cy > 0.0 && e.Cy < (float h) - 1.0 &&
+ e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
ellipses.Add e
| _ -> ()
| _ -> ()
//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 ]
use magnitudesByte = ((magnitudes / !max) * 255.0).Convert<byte>() // Otsu from OpenCV only support 'byte'.
use edges = new Matrix<byte>(xEdges.Size)
- let threshold = CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
- thin edges
- removeArea edges 12
+ let threshold = CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+
+// let filteredGreenMat = new Matrix<float32>(filteredGreen.Size)
+// filteredGreen.CopyTo(filteredGreenMat)
+ let parasites = ParasitesMarker.find green filteredGreen config
+
+ saveImg parasites.darkStain "parasites_dark_stain.png"
+ saveImg parasites.stain "parasites_stain.png"
+ saveImg parasites.infection "parasites_infection.png"
+
+ logTime "Finding edges" (fun() ->
+ thin edges)
+
+ logTime "Removing small connected components" (fun () ->
+ removeArea edges 12)
saveMat (edges * 255.0) "edges.png"
let radiusRange = config.scale * 20.0, config.scale * 40.0
let windowSize = roundInt (1.6 * (snd radiusRange))
- let factorNbPick = 1.0;
- let ellipses = Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick
-
- drawEllipse img (List.head ellipses) (Bgr(0.0, 255.0, 255.0))
- saveImg img "ellipses.png"
+ let factorNbPick = 1.5
+ let ellipses = logTime "Finding ellipses" (fun () ->
+ Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick)
+
+ drawEllipses img ellipses (Bgr(0.0, 255.0, 255.0))
+ //saveImg img "ellipses.png"
{ RBCPositions = []; infectedRBCPositions = []; img = img }
open Utils
+// Normalize image values between 0uy and 255uy.
+let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
+ 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<Gray, byte>()
+
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)
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 cosAlpha = cos e.Alpha
+ let sinAlpha = sin e.Alpha
let mutable x0 = 0.0
let mutable y0 = 0.0
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
+ 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
then
first_iteration <- false
x0 <- x
- y0 <- y
\ No newline at end of file
+ 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
module KMedians
+open System.Collections.Generic
+open System.Drawing
+
open Emgu.CV
open Emgu.CV.Structure
-open System.Drawing
+type Result = {
+ fg: Image<Gray, byte>
+ median_bg: float
+ median_fg: float
+ d_fg: Image<Gray, float32> } // Distances to median_fg.
-let kmedians (mat: Matrix<float32>) (fgFactor: float) : Matrix<byte> =
+let kmedians (img: Image<Gray, float32>) (fgFactor: float) : Result =
let nbIteration = 3
-
- let min = ref 0.0
- let minLocation = ref <| Point()
- let max = ref 0.0
- let maxLocation = ref <| Point()
- mat.MinMax(min, max, minLocation, maxLocation)
+ let w = img.Width
+ let h = img.Height
+
+ let min = ref [| 0.0 |]
+ let minLocation = ref <| [| Point() |]
+ let max = ref [| 0.0 |]
+ let maxLocation = ref <| [| Point() |]
+ img.MinMax(min, max, minLocation, maxLocation)
- let mutable bgValue = !max - (!max - !min) / 4.0
- let mutable fgValue = !min + (!max - !min) / 4.0
- let mutable d_bg = new Matrix<float32>(mat.Size)
- let mutable d_fg = new Matrix<float32>(mat.Size)
-
- for i in 1..3 do
- CvInvoke.Pow(mat - bgValue, 2.0, d_bg)
- CvInvoke.Pow(mat - fgValue, 2.0, d_fg)
- let fg = (d_fg * fgFactor).Cmp(d_bg, CvEnum.CmpType.LessThan)
-
- printfn "test"
-
-
- (*d_bg = (imgFloat - color_bg) .^ 2.0;
- d_fg = (imgFloat - color_fg) .^ 2.0;
- fg = d_fg * Background_weight < d_bg;
- imgFilteredFg = imgFloat;
- imgFilteredFg(~fg) = nan;
- color_fg = median(reshape(imgFilteredFg, [], 1), 'omitnan');
- imgFilteredBg = imgFloat;
- imgFilteredBg(fg) = nan;
- color_bg = median(reshape(imgFilteredBg, [], 1), 'omitnan');
- *)
-
+ 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)
- new Matrix<byte>(mat.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 }
--- /dev/null
+module KdTree
+
+open System
+
+type I2DCoords =
+ abstract X : float
+ abstract Y : float
+
+// Compare 'e1' and 'e2' by X.
+let cmpX (e1: I2DCoords) (e2: I2DCoords) : int =
+ match e1.X.CompareTo(e2.X) with
+ | 0 -> match e1.Y.CompareTo(e2.Y) with
+ | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
+ | v -> v
+ | v -> v
+
+// Compare 'e1' and 'e2' by Y.
+let cmpY (e1: I2DCoords) (e2: I2DCoords) : int =
+ match e1.Y.CompareTo(e2.Y) with
+ | 0 -> match e1.X.CompareTo(e2.X) with
+ | 0 -> e1.GetHashCode().CompareTo(e2.GetHashCode())
+ | v -> v
+ | v -> v
+
+type Region = { minX: float; maxX: float; minY: float; maxY: float } with
+ member this.Contains px py : bool =
+ px >= this.minX && px <= this.maxX &&
+ py >= this.minY && py <= this.maxY
+
+ member this.IsSub otherRegion : bool =
+ this.minX >= otherRegion.minX && this.maxX <= otherRegion.maxX &&
+ this.minY >= otherRegion.minY && this.maxY <= otherRegion.maxY
+
+ member this.Intersects otherRegion : bool =
+ this.minX < otherRegion.maxX && this.maxX >= otherRegion.minX &&
+ this.minY < otherRegion.maxY && this.maxY >= otherRegion.minY
+
+type Tree<'a when 'a :> I2DCoords> =
+ | Node of float * Tree<'a> * Tree<'a>
+ | Leaf of 'a
+
+ static member buildTree (l: 'a list) : Tree<'a> =
+ let xSorted = List.toArray l
+ let ySorted = List.toArray l
+ Array.sortInPlaceWith cmpX xSorted
+ Array.sortInPlaceWith cmpY ySorted
+
+ let rec buildTreeFromSortedArray (pXSorted: 'a[]) (pYSorted: 'a[]) (depth: int) : Tree<'a> =
+ if pXSorted.Length = 1
+ then
+ Leaf pXSorted.[0]
+ else
+ if depth % 2 = 1 // 'depth' is odd -> vertical splitting else horizontal splitting.
+ then
+ let leftX, rightX = Array.splitAt ((pXSorted.Length + 1) / 2) pXSorted
+ let splitElement = Array.last leftX
+ let leftY, rightY = Array.partition (fun (e: 'a) -> cmpX e splitElement <= 0) pYSorted // FIXME: Maybe this operation can be optimized.
+ Node (splitElement.X, buildTreeFromSortedArray leftX leftY (depth + 1), buildTreeFromSortedArray rightX rightY (depth + 1))
+ else
+ let downY, upY = Array.splitAt ((pYSorted.Length + 1) / 2) pYSorted
+ let splitElement = Array.last downY
+ let downX, upX = Array.partition (fun (e: 'a) -> cmpY e splitElement <= 0) pXSorted // FIXME: Maybe this operation can be optimized.
+ Node (splitElement.Y, buildTreeFromSortedArray downX downY (depth + 1), buildTreeFromSortedArray upX upY (depth + 1))
+
+ buildTreeFromSortedArray xSorted ySorted 1
+
+ static member search (tree: Tree<'a>) (searchRegion: Region) : 'a list =
+ let rec valuesFrom (tree: Tree<'a>) : 'a list =
+ match tree with
+ | Leaf v -> [v]
+ | Node (_, part1, part2) -> (valuesFrom part1) @ (valuesFrom part2)
+
+ let rec searchWithRegion (tree: Tree<'a>) (currentRegion: Region) (depth: int) : 'a list =
+ match tree with
+ | Leaf v -> if searchRegion.Contains v.X v.Y then [v] else []
+ | Node (splitValue, part1, part2) ->
+ let valuesInRegion (region: Region) (treeRegion: Tree<'a>) =
+ if region.IsSub searchRegion
+ then
+ valuesFrom treeRegion
+ elif region.Intersects searchRegion
+ then
+ searchWithRegion treeRegion region (depth + 1)
+ else
+ []
+
+ if depth % 2 = 1 // Vertical splitting.
+ then
+ let leftRegion = { currentRegion with maxX = splitValue }
+ let rightRegion = { currentRegion with minX = splitValue }
+ (valuesInRegion leftRegion part1) @ (valuesInRegion rightRegion part2)
+ else // Horizontal splitting.
+ let downRegion = { currentRegion with maxY = splitValue }
+ 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
+
+
+///// Tests. TODO: to put in a unit test.
+
+type Point (x: float, y: float) =
+ interface I2DCoords with
+ member this.X = x
+ member this.Y = y
+
+ override this.ToString () =
+ sprintf "(%.1f, %.1f)" x y
+
+// TODO: test with identical X or Y coords
+let test () =
+ let pts = [
+ Point(1.0, 1.0)
+ Point(2.0, 2.0)
+ Point(1.5, 3.6)
+ Point(3.0, 3.2)
+ Point(4.0, 4.0)
+ Point(3.5, 1.5)
+ Point(2.5, 0.5) ]
+
+ 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.
+ Utils.dprintfn "s1: %A" s1
+
+ let s2 = Tree.search tree { 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 }
+ Utils.dprintfn "s3: %A" s3
+
+let test2 () =
+ let pts = [
+ Point(1.0, 1.0)
+ Point(1.0, 2.0)
+ Point(1.0, 3.0) ]
+
+ 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 }
+ Utils.dprintfn "s1: %A" s1
+
+ let s2 = Tree.search tree { 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 }
+ Utils.dprintfn "s3: %A" s3
+
+ let s4 = Tree.search tree { minX = 0.0; maxX = 2.0; minY = 0.0; maxY = 4.0 }
+ Utils.dprintfn "s4: %A" s4
+
open System
open System.Linq
+open System.Collections
open System.Collections.Generic
open Types
open Utils
-type EllipseScore = { mutable matchingScore: float; e: Ellipse }
+let matchingScoreThreshold1 = 0.6
+let matchingScoreThreshold2 = 1.0
-let matchingScoreThreshold1 = 0.5
-let matchingScoreThreshold2 = 2.0
+type EllipseScore (matchingScore: float, e: Ellipse) =
+ let mutable matchingScore = matchingScore
-let ellipseArea e = e.a * e.b * Math.PI
+ member this.MatchingScore = matchingScore
+ member this.Ellipse = e
+ member val Processed = false with get, set
+ member val Removed = false with get, set
-type MatchingEllipses () =
+ member this.AddMatchingScore(score: float) =
+ matchingScore <- matchingScore + score
+
+ interface KdTree.I2DCoords with
+ member this.X = e.Cx
+ member this.Y = e.Cy
+
+type MatchingEllipses (radiusMin: float) =
let ellipses = List<EllipseScore>()
member this.Add (e: Ellipse) =
-
- // dprintfn "new ellipse: %A, nb: %A" e ellipses.Count
-
- let areaE = ellipseArea e
-
- let mutable matchingScoreE = 0.0
+ ellipses.Add(EllipseScore(0.0, e))
- for other in ellipses do
- let areaOther = ellipseArea other.e
- let commonArea = EEOver.EEOverlapArea e other.e
- let matchingScore = (2.0 * commonArea / (areaE + areaOther)) ** 2.0
- if matchingScore >= matchingScoreThreshold1
+ 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
+ let commonArea = EEOver.EEOverlapArea e.Ellipse other.Ellipse
+ 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
+ then
+ 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
- other.matchingScore <- other.matchingScore + matchingScore
- matchingScoreE <- matchingScoreE + matchingScore
- printfn "Score"
-
- ellipses.Add({ matchingScore = matchingScoreE; e = e })
+ 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
- member this.Ellipses : Ellipse list =
dprintfn "Number of ellipse: %A" ellipses.Count
- (*let sortedEllipses =
- List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses)
- |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore))*)
-
- List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses)
- |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore))
- |> List.map (fun { e = e } -> e)
+ List.ofSeq ellipses |> List.map (fun e -> e.Ellipse)
- // ellipses.Where(fun e -> e.matchingScore >= matchingScoreThreshold2)
- // ([], fun acc { matchingScore = score; e = e } -> if score >= matchingScoreThreshold2 then e :: acc else acc)
<OutputType>WinExe</OutputType>
<RootNamespace>Parasitemia</RootNamespace>
<AssemblyName>Parasitemia</AssemblyName>
- <TargetFrameworkVersion>v4.6</TargetFrameworkVersion>
+ <TargetFrameworkVersion>v4.6.1</TargetFrameworkVersion>
<AutoGenerateBindingRedirects>true</AutoGenerateBindingRedirects>
<TargetFSharpCoreVersion>4.4.0.0</TargetFSharpCoreVersion>
<Name>Parasitemia</Name>
<Compile Include="ImgTools.fs" />
<Compile Include="Config.fs" />
<Compile Include="EEOver.fs" />
+ <Compile Include="KdTree.fs" />
<Compile Include="MatchingEllipses.fs" />
<Compile Include="Ellipse.fs" />
<Compile Include="KMedians.fs" />
+ <Compile Include="ParasitesMarker.fs" />
<Compile Include="ImageAnalysis.fs" />
<Compile Include="Program.fs" />
<None Include="App.config" />
<Reference Include="Emgu.Util">
<HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll</HintPath>
</Reference>
- <Reference Include="FSharp.Data.TypeProviders" />
+ <Reference Include="FSharp.Core, Version=$(TargetFSharpCoreVersion), Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a">
+ </Reference>
<Reference Include="FSharp.ViewModule.Core.Wpf">
<HintPath>..\packages\FSharp.ViewModule.Core.0.9.9.1\lib\net45\FSharp.ViewModule.Core.Wpf.dll</HintPath>
<Private>True</Private>
<HintPath>..\packages\log4net.1.2.10\lib\2.0\log4net.dll</HintPath>
<Private>True</Private>
</Reference>
- <Reference Include="mscorlib" />
- <Reference Include="FSharp.Core, Version=$(TargetFSharpCoreVersion), Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a">
+ <Reference Include="MathNet.Numerics">
+ <HintPath>..\packages\MathNet.Numerics.3.9.0\lib\net40\MathNet.Numerics.dll</HintPath>
<Private>True</Private>
</Reference>
+ <Reference Include="MathNet.Numerics.FSharp">
+ <HintPath>..\packages\MathNet.Numerics.FSharp.3.9.0\lib\net40\MathNet.Numerics.FSharp.dll</HintPath>
+ <Private>True</Private>
+ </Reference>
+ <Reference Include="mscorlib" />
<Reference Include="PresentationCore" />
<Reference Include="PresentationFramework" />
<Reference Include="System" />
<Reference Include="System.Core" />
<Reference Include="System.Data" />
- <Reference Include="System.Data.Linq" />
+ <Reference Include="System.Data.DataSetExtensions" />
<Reference Include="System.Drawing" />
<Reference Include="System.Numerics" />
<Reference Include="System.Windows.Interactivity">
--- /dev/null
+module ParasitesMarker
+
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+type Result = {
+ darkStain: Image<Gray, byte>
+ stain: Image<Gray, byte>
+ infection: Image<Gray, byte> }
+
+let find (green: Image<Gray, float32>) (filteredGreen: Image<Gray, float32>) (config: Config.Config) : Result =
+
+ // 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 } = KMedians.kmedians filteredGreen 1.0
+ 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 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
+
+ { darkStain = darkStain;
+ stain = findSmears config.stainSigma config.stainLevel
+ infection = findSmears config.infectionSigma config.infectionLevel }
+
+
+
+
Utils.log <- (fun m -> log mainWindow m)
let config = {
+ scale = 1.0
+
doGSigma1 = 1.5
doGSigma2 = 20.0
doGLowFreqPercentageReduction = 0.9
- scale = 1.0
- }
+
+ darkStainLevel = 0.839
+
+ stainSigma = 10.0
+ stainLevel = 0.9
+ stainSpreadRequired = 3.0
+
+ infectionSigma = 2.0
+ infectionLevel = 0.762
+ infectionPixelsRequired = 1 }
- //use img = new Image<Rgb, byte>("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png")
+ ///// 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_2.png")
- use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/two_rbc_1.png")
-
- let result = ImageAnalysis.doAnalysis img config
+ //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 = ImageAnalysis.doAnalysis img config
+
display mainWindow result.img
mainWindow.Root.Show()
-
app.Run() |> ignore
open System
-type Ellipse = { cx: float; cy: float; a: float; b: float; alpha: float }
-
+type Ellipse (cx: float, cy: float, a: float, b: float, alpha: float) =
+ member this.Cx = cx
+ member this.Cy = cy
+ member this.A = a
+ member this.B = b
+ member this.Alpha = alpha
+ member this.Area = a * b * Math.PI
-//type PointImg = { x: int; y: int }
\ No newline at end of file
+ // Does the ellipse contain the point (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
\ No newline at end of file
<packages>
<package id="Castle.Core" version="1.1.0" targetFramework="net46" />
<package id="Expression.Blend.Sdk" version="1.0.2" targetFramework="net46" />
+ <package id="FSharp.Core" version="3.1.2.5" targetFramework="net46" />
<package id="FSharp.ViewModule.Core" version="0.9.9.1" targetFramework="net46" />
<package id="FsXaml.Wpf" version="0.9.9" targetFramework="net46" />
<package id="log4net" version="1.2.10" targetFramework="net46" />
+ <package id="MathNet.Numerics" version="3.9.0" targetFramework="net461" />
+ <package id="MathNet.Numerics.FSharp" version="3.9.0" targetFramework="net461" />
</packages>
\ No newline at end of file
projects / master-thesis.git / commitdiff