From: Greg Burri Date: Tue, 26 Jan 2016 13:55:21 +0000 (+0100) Subject: Split the module 'ImgTools' in many modules. X-Git-Tag: 1.0.11~30 X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=commitdiff_plain;h=3f8b0d281b3058faf23dbd0363de440bd04c6574 Split the module 'ImgTools' in many modules. --- diff --git a/Parasitemia/ParasitemiaCore/Classifier.fs b/Parasitemia/ParasitemiaCore/Classifier.fs index 615840c..f4b0844 100644 --- a/Parasitemia/ParasitemiaCore/Classifier.fs +++ b/Parasitemia/ParasitemiaCore/Classifier.fs @@ -196,7 +196,6 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: then parasiteArea <- parasiteArea + 1 - let cellClass = if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement) then @@ -204,7 +203,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: elif nucleusPixels.Count > 0 && parasiteArea >= minimumParasiteArea then - let infectionToRemove = ImgTools.connectedComponents parasites.parasite nucleusPixels + let infectionToRemove = Morpho.connectedComponents parasites.parasite nucleusPixels for p in infectionToRemove do nucleusData.[p.Y, p.X, 0] <- 0uy InfectedRBC diff --git a/Parasitemia/ParasitemiaCore/EEOver.fs b/Parasitemia/ParasitemiaCore/EEOver.fs index eb8ac0b..6862dc5 100644 --- a/Parasitemia/ParasitemiaCore/EEOver.fs +++ b/Parasitemia/ParasitemiaCore/EEOver.fs @@ -1,4 +1,5 @@ -module ParasitemiaCore.EEOver +// Translation from https://github.com/chraibi/EEOver. +module ParasitemiaCore.EEOver open System @@ -508,7 +509,9 @@ let private biquadroots (p: float[]) (r: float[,]) = quad () -// Return a tuple (area, x intersections, y intersections) +/// +/// Return a tuple (area, x intersections, y intersections). +/// let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option = let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha diff --git a/Parasitemia/ParasitemiaCore/Ellipse.fs b/Parasitemia/ParasitemiaCore/Ellipse.fs index 1059e3b..4e48230 100644 --- a/Parasitemia/ParasitemiaCore/Ellipse.fs +++ b/Parasitemia/ParasitemiaCore/Ellipse.fs @@ -14,142 +14,12 @@ open Config open MatchingEllipses open Const -type private SearchExtremum = Minimum | Maximum - -let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) (xmax: float) (searchExtremum: SearchExtremum) : (float * float) = - let gr = 1. / 1.6180339887498948482 - let mutable a = xmin - let mutable b = xmax - let mutable c = b - gr * (b - a) - let mutable d = a + gr * (b - a) - - for i in 1 .. nbIter do - let mutable fc = f c - let mutable fd = f d - - if searchExtremum = Maximum - then - let tmp = fc - fc <- fd - fd <- tmp - - if fc < fd - then - b <- d - d <- c - c <- b - gr * (b - a) - else - a <- c - c <- d - d <- a + gr * (b - a) - - let x = (b + a) / 2. - x, f x - -// Ellipse.A is always equal or greater than Ellipse.B. -// Ellipse.Alpha is between 0 and Pi. +/// +/// Try to build an ellipse from three points and two tangents passing by the first and the second point. +/// '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 = 10 // 3 - let accuracy_extremum_search_2 = 10 // 4 - - // p3 as the referencial. - let p1x = p1x - p3x - let p1y = p1y - p3y - - let p2x = p2x - p3x - let p2y = p2y - p3y - - // Convert to polar coordinates. - let alpha1 = atan m1 - let alpha2 = atan m2 - - let r1 = sqrt (p1x ** 2. + p1y ** 2.) - let theta1 = atan2 p1y p1x - - let r2 = sqrt (p2x ** 2. + p2y ** 2.) - let theta2 = atan2 p2y p2x - - let valid = - 4. * 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. * sin (theta1 - theta2) ** 2. < 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. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.) - - let rabs = r >> abs - - // We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval. - let thetaTan, _ = goldenSectionSearch rabs accuracy_extremum_search_1 0. Math.PI Maximum - let rTan = r thetaTan - - let PTanx = rTan * cos thetaTan - let PTany = rTan * sin thetaTan - - let d1a = tan alpha1 - let d1b = -d1a * p1x + p1y - - let d2a = tan alpha2 - let d2b = -d2a * p2x + p2y - - let d3a = -1. / tan thetaTan - let d3b = -d3a * PTanx + PTany - - let Ux = -(d1b - d2b) / (d1a - d2a) - let Uy = -(d2a * d1b - d1a * d2b) / (d1a - d2a) - - let Vx = -(d1b - d3b) / (d1a - d3a) - let Vy = -(d3a * d1b - d1a * d3b) / (d1a - d3a) - - let Wx = p1x + (p2x - p1x) / 2. - let Wy = p1y + (p2y - p1y) / 2. - - let Zx = p1x + (PTanx - p1x) / 2. - let Zy = p1y + (PTany - p1y) / 2. - - let va = -(-Vy + Zy) / (Vx - Zx) - let vb = -(Zx * Vy - Vx * Zy) / (Vx - Zx) - - let ua = -(-Uy + Wy) / (Ux - Wx) - let ub = -(Wx * Uy - Ux * Wy) / (Ux - Wx) - - let cx = -(vb - ub) / (va - ua) - let cy = -(ua * vb - va * ub) / (va - ua) - - let rc = sqrt (cx ** 2. + cy ** 2.) - let psi = atan2 cy cx - - let rellipse theta = - sqrt ( - rc ** 2. + (r1 ** 2. * r2 ** 2. * (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. * sin (theta1 - theta2) ** 2.) / - (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.) ** 2. - - (2. * 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. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.)) - - // 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. (Math.PI / 2.) Maximum // Pi/2 and not pi because the period is Pi. - let r2eTheta, r2e = goldenSectionSearch rellipse accuracy_extremum_search_2 0. (Math.PI / 2.) Minimum - - let rr1e = r r1eTheta - let r1ex = rr1e * cos r1eTheta - let r1ey = rr1e * sin r1eTheta - let mutable alpha = atan ((r1ey - cy) / (r1ex - cx)) - if alpha < 0. - then - alpha <- alpha + Math.PI - - // Ride off the p3 referential. - let cx = cx + p3x - let cy = cy + p3y - - Some (Types.Ellipse(float32 cx, float32 cy, float32 r1e, float32 r2e, float32 alpha)) - else - None - -let ellipse2 (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option = let p0 = pointFromTwoLines (Types.Line(float32 m1, float32 (p1y - m1 * p1x))) (Types.Line(float32 m2, float32(p2y - m2 * p2x))) let p0x, p0y = float p0.X, float p0.Y @@ -259,7 +129,6 @@ let private areVectorsValid (p1x: float32) (p1y: float32) (p2x: float32) (p2y: f else Some (m1, m2) - let find (edges: Matrix) (xGradient: Matrix) (yGradient: Matrix) @@ -336,7 +205,7 @@ let find (edges: Matrix) then match areVectorsValid (float32 p1xf) (float32 p1yf) (float32 p2xf) (float32 p2yf) -xDirData.[p1.Y, p1.X] -yDirData.[p1.Y, p1.X] -xDirData.[p2.Y, p2.X] -yDirData.[p2.Y, p2.X] with | Some (m1, m2) -> - match ellipse2 p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with + match ellipse p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with | Some e when e.Cx > 0.f && e.Cx < w_f - 1.f && e.Cy > 0.f && e.Cy < h_f - 1.f && e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance -> ellipses.Add e diff --git a/Parasitemia/ParasitemiaCore/Granulometry.fs b/Parasitemia/ParasitemiaCore/Granulometry.fs index 0afe5cc..b325534 100644 --- a/Parasitemia/ParasitemiaCore/Granulometry.fs +++ b/Parasitemia/ParasitemiaCore/Granulometry.fs @@ -71,7 +71,7 @@ let findRadiusByAreaClosing (img: Image) (radiusRange: int * int) let mutable maxDiff = 0.f let mutable max_r = r1 - ImgTools.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (fun r diff -> + Morpho.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (fun r diff -> if r <> r1 && diff > maxDiff then maxDiff <- diff diff --git a/Parasitemia/ParasitemiaCore/ImgTools.fs b/Parasitemia/ParasitemiaCore/ImgTools.fs deleted file mode 100644 index cbb65b4..0000000 --- a/Parasitemia/ParasitemiaCore/ImgTools.fs +++ /dev/null @@ -1,1018 +0,0 @@ -module ParasitemiaCore.ImgTools - -open System -open System.Drawing -open System.Collections.Generic -open System.Linq - -open Emgu.CV -open Emgu.CV.Structure - -open Heap -open Const -open Types -open Utils - -let normalize (img: Image) (upperLimit: float) : Image = - let min = ref [| 0.0 |] - let minLocation = ref <| [| Point() |] - let max = ref [| 0.0 |] - let maxLocation = ref <| [| Point() |] - img.MinMax(min, max, minLocation, maxLocation) - let normalized = (img - (!min).[0]) / ((!max).[0] - (!min).[0]) - if upperLimit = 1.0 - then normalized - else upperLimit * normalized - -let mergeChannels (img: Image) (rgbWeights: float * float * float) : Image = - match rgbWeights with - | 1., 0., 0. -> img.[2] - | 0., 1., 0. -> img.[1] - | 0., 0., 1. -> img.[0] - | redFactor, greenFactor, blueFactor -> - let result = new Image(img.Size) - CvInvoke.AddWeighted(result, 1., img.[2], redFactor, 0., result) - CvInvoke.AddWeighted(result, 1., img.[1], greenFactor, 0., result) - CvInvoke.AddWeighted(result, 1., img.[0], blueFactor, 0., result) - result - -let mergeChannelsWithProjection (img: Image) (v1r: float32, v1g: float32, v1b: float32) (v2r: float32, v2g: float32, v2b: float32) (upperLimit: float) : Image = - let vr, vg, vb = v2r - v1r, v2g - v1g, v2b - v1b - let vMagnitude = sqrt (vr ** 2.f + vg ** 2.f + vb ** 2.f) - let project (r: float32) (g: float32) (b: float32) = ((r - v1r) * vr + (g - v1g) * vg + (b - v1b) * vb) / vMagnitude - let result = new Image(img.Size) - // TODO: Essayer en bindant Data pour gagner du temps - for i in 0 .. img.Height - 1 do - for j in 0 .. img.Width - 1 do - result.Data.[i, j, 0] <- project img.Data.[i, j, 2] img.Data.[i, j, 1] img.Data.[i, j, 0] - normalize result upperLimit - -// Normalize image values between 0uy and 255uy. -let normalizeAndConvert (img: Image) : Image = - (normalize (img.Convert()) 255.).Convert() - -let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = - img.Save(filepath) - -let saveMat (mat: Matrix<'TDepth>) (filepath: string) = - use img = new Image(mat.Size) - mat.CopyTo(img) - saveImg img filepath - -type Histogram = { data: int[]; total: int; sum: int; min: float32; max: float32 } - -let histogramImg (img: Image) (nbSamples: int) : Histogram = - let imgData = img.Data - - let min, max = - let min = ref [| 0.0 |] - let minLocation = ref <| [| Point() |] - let max = ref [| 0.0 |] - let maxLocation = ref <| [| Point() |] - img.MinMax(min, max, minLocation, maxLocation) - float32 (!min).[0], float32 (!max).[0] - - let inline bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for i in 0 .. img.Height - 1 do - for j in 0 .. img.Width - 1 do - let p = bin imgData.[i, j, 0] - data.[p] <- data.[p] + 1 - - { data = data; total = img.Height * img.Width; sum = Array.sum data; min = min; max = max } - -let histogramMat (mat: Matrix) (nbSamples: int) : Histogram = - let matData = mat.Data - - let min, max = - let min = ref 0.0 - let minLocation = ref <| Point() - let max = ref 0.0 - let maxLocation = ref <| Point() - mat.MinMax(min, max, minLocation, maxLocation) - float32 !min, float32 !max - - let inline bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for i in 0 .. mat.Height - 1 do - for j in 0 .. mat.Width - 1 do - let p = bin matData.[i, j] - data.[p] <- data.[p] + 1 - - { data = data; total = mat.Height * mat.Width; sum = Array.sum data; min = min; max = max } - -let histogram (values: float32 seq) (nbSamples: int) : Histogram = - let mutable min = Single.MaxValue - let mutable max = Single.MinValue - let mutable n = 0 - - for v in values do - n <- n + 1 - if v < min then min <- v - if v > max then max <- v - - let inline bin (x: float32) : int = - let p = int ((x - min) / (max - min) * float32 nbSamples) - if p >= nbSamples then nbSamples - 1 else p - - let data = Array.zeroCreate nbSamples - - for v in values do - let p = bin v - data.[p] <- data.[p] + 1 - - { data = data; total = n; sum = Array.sum data; min = min; max = max } - -let otsu (hist: Histogram) : float32 * float32 * float32 = - let mutable sumB = 0 - let mutable wB = 0 - let mutable maximum = 0.0 - let mutable level = 0 - let sum = hist.data |> Array.mapi (fun i v -> i * v |> float) |> Array.sum - - for i in 0 .. hist.data.Length - 1 do - wB <- wB + hist.data.[i] - if wB <> 0 - then - let wF = hist.total - wB - if wF <> 0 - then - sumB <- sumB + i * hist.data.[i] - let mB = (float sumB) / (float wB) - let mF = (sum - float sumB) / (float wF) - let between = (float wB) * (float wF) * (mB - mF) ** 2.; - if between >= maximum - then - level <- i - maximum <- between - - let mean1 = - let mutable sum = 0 - let mutable nb = 0 - for i in 0 .. level - 1 do - sum <- sum + i * hist.data.[i] - nb <- nb + hist.data.[i] - (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) - - let mean2 = - let mutable sum = 0 - let mutable nb = 0 - for i in level + 1 .. hist.data.Length - 1 do - sum <- sum + i * hist.data.[i] - nb <- nb + hist.data.[i] - (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) - - let toFloat l = - float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min - - toFloat level, toFloat mean1, toFloat mean2 - -/// -/// Remove M-adjacent pixels. It may be used after thinning. -/// -let suppressMAdjacency (img: Matrix) = - 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 - -/// -/// Find edges of an image by using the Canny approach. -/// The thresholds are automatically defined with otsu on gradient magnitudes. -/// -/// -let findEdges (img: Image) : Matrix * Matrix * Matrix = - let w = img.Width - let h = img.Height - - use sobelKernel = - new Matrix(array2D [[ 1.0f; 0.0f; -1.0f ] - [ 2.0f; 0.0f; -2.0f ] - [ 1.0f; 0.0f; -1.0f ]]) - - let xGradient = new Matrix(img.Size) - let yGradient = new Matrix(img.Size) - CvInvoke.Filter2D(img, xGradient, sobelKernel, Point(1, 1)) - CvInvoke.Filter2D(img, yGradient, sobelKernel.Transpose(), Point(1, 1)) - - use magnitudes = new Matrix(xGradient.Size) - use angles = new Matrix(xGradient.Size) - CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. - - let thresholdHigh, thresholdLow = - let sensibilityHigh = 0.1f - let sensibilityLow = 0.0f - let threshold, _, _ = otsu (histogramMat magnitudes 300) - threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold) - - // Non-maximum suppression. - use nms = new Matrix(xGradient.Size) - - let nmsData = nms.Data - let anglesData = angles.Data - let magnitudesData = magnitudes.Data - let xGradientData = xGradient.Data - let yGradientData = yGradient.Data - - for i in 0 .. h - 1 do - nmsData.[i, 0] <- 0uy - nmsData.[i, w - 1] <- 0uy - - for j in 0 .. w - 1 do - nmsData.[0, j] <- 0uy - nmsData.[h - 1, j] <- 0uy - - for i in 1 .. h - 2 do - for j in 1 .. w - 2 do - let vx = xGradientData.[i, j] - let vy = yGradientData.[i, j] - if vx <> 0.f || vy <> 0.f - then - let angle = anglesData.[i, j] - - let vx', vy' = abs vx, abs vy - let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy' - let ratio1 = 1.f - ratio2 - - let mNeigbors (sign: int) : float32 = - if angle < PI / 4.f - then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign] - elif angle < PI / 2.f - then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j] - elif angle < 3.f * PI / 4.f - then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign] - elif angle < PI - then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign] - elif angle < 5.f * PI / 4.f - then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign] - elif angle < 3.f * PI / 2.f - then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j] - elif angle < 7.f * PI / 4.f - then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign] - else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign] - - let m = magnitudesData.[i, j] - if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 - then - nmsData.[i, j] <- 1uy - - // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection). - - let edges = new Matrix(xGradient.Size) - let edgesData = edges.Data - - // Hysteresis thresholding. - let toVisit = Stack() - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh - then - nmsData.[i, j] <- 0uy - toVisit.Push(Point(j, i)) - while toVisit.Count > 0 do - let p = toVisit.Pop() - edgesData.[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 && nmsData.[ni, nj] = 1uy - then - nmsData.[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) - -let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = - for p in points do - img.Data.[p.Y, p.X, 0] <- intensity - -type ExtremumType = - | Maxima = 1 - | Minima = 2 - -let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumerable = - let w = img.Width - let h = img.Height - let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] - - let imgData = img.Data - let suppress: bool[,] = Array2D.zeroCreate h w - - let result = List>() - - let flood (start: Point) : List> = - let sameLevelToCheck = Stack() - let betterLevelToCheck = Stack() - betterLevelToCheck.Push(start) - - let result' = List>() - - while betterLevelToCheck.Count > 0 do - let p = betterLevelToCheck.Pop() - if not suppress.[p.Y, p.X] - then - suppress.[p.Y, p.X] <- true - sameLevelToCheck.Push(p) - let current = List() - - let mutable betterExists = false - - while sameLevelToCheck.Count > 0 do - let p' = sameLevelToCheck.Pop() - let currentLevel = imgData.[p'.Y, p'.X, 0] - current.Add(p') |> ignore - for i, j in se do - let ni = i + p'.Y - let nj = j + p'.X - if ni >= 0 && ni < h && nj >= 0 && nj < w - then - let level = imgData.[ni, nj, 0] - let notSuppressed = not suppress.[ni, nj] - - if level = currentLevel && notSuppressed - then - suppress.[ni, nj] <- true - sameLevelToCheck.Push(Point(nj, ni)) - elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel - then - betterExists <- true - if notSuppressed - then - betterLevelToCheck.Push(Point(nj, ni)) - - if not betterExists - then - result'.Add(current) - result' - - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - let maxima = flood (Point(j, i)) - if maxima.Count > 0 - then - result.AddRange(maxima) - - result.Select(fun l -> Points(l)) - -let findMaxima (img: Image) : IEnumerable = - findExtremum img ExtremumType.Maxima - -let findMinima (img: Image) : IEnumerable = - findExtremum img ExtremumType.Minima - -type PriorityQueue () = - let size = 256 - let q: Points[] = Array.init size (fun i -> Points()) - let mutable highest = -1 // Value of the first elements of 'q'. - let mutable lowest = size - - member this.NextMax () : byte * Point = - if this.IsEmpty - then - invalidOp "Queue is empty" - else - let l = q.[highest] - let next = l.First() - l.Remove(next) |> ignore - let value = byte highest - - if l.Count = 0 - then - highest <- highest - 1 - while highest > lowest && q.[highest].Count = 0 do - highest <- highest - 1 - if highest = lowest - then - highest <- -1 - lowest <- size - - value, next - - member this.NextMin () : byte * Point = - if this.IsEmpty - then - invalidOp "Queue is empty" - else - let l = q.[lowest + 1] - let next = l.First() - l.Remove(next) |> ignore - let value = byte (lowest + 1) - - if l.Count = 0 - then - lowest <- lowest + 1 - while lowest < highest && q.[lowest + 1].Count = 0 do - lowest <- lowest + 1 - if highest = lowest - then - highest <- -1 - lowest <- size - - value, next - - member this.Max = - highest |> byte - - member this.Min = - lowest + 1 |> byte - - member this.Add (value: byte) (p: Point) = - let vi = int value - - if vi > highest - then - highest <- vi - if vi <= lowest - then - lowest <- vi - 1 - - q.[vi].Add(p) |> ignore - - member this.Remove (value: byte) (p: Point) = - let vi = int value - if q.[vi].Remove(p) && q.[vi].Count = 0 - then - if vi = highest - then - highest <- highest - 1 - while highest > lowest && q.[highest].Count = 0 do - highest <- highest - 1 - elif vi - 1 = lowest - then - lowest <- lowest + 1 - while lowest < highest && q.[lowest + 1].Count = 0 do - lowest <- lowest + 1 - - if highest = lowest // The queue is now empty. - then - highest <- -1 - lowest <- size - - member this.IsEmpty = - highest = -1 - - member this.Clear () = - while highest > lowest do - q.[highest].Clear() - highest <- highest - 1 - highest <- -1 - lowest <- size - -type private AreaState = - | Removed = 1 - | Unprocessed = 2 - | Validated = 3 - -type private AreaOperation = - | Opening = 1 - | Closing = 2 - -[] -type private Area (elements: Points) = - member this.Elements = elements - member val Intensity = None with get, set - member val State = AreaState.Unprocessed with get, set - -let private areaOperation (img: Image) (area: int) (op: AreaOperation) = - let w = img.Width - let h = img.Height - let imgData = img.Data - let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] - - let areas = List((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area) - - let pixels: Area[,] = Array2D.create h w null - for m in areas do - for e in m.Elements do - pixels.[e.Y, e.X] <- m - - let queue = PriorityQueue() - - let addEdgeToQueue (elements: Points) = - for p in elements do - for i, j in se do - let ni = i + p.Y - let nj = j + p.X - let p' = Point(nj, ni) - if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p')) - then - queue.Add (imgData.[ni, nj, 0]) p' - - // Reverse order is quicker. - for i in areas.Count - 1 .. -1 .. 0 do - let m = areas.[i] - if m.Elements.Count <= area && m.State <> AreaState.Removed - then - queue.Clear() - addEdgeToQueue m.Elements - - let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min - let nextElements = Points() - - let mutable stop = false - while not stop do - let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin () - let mutable merged = false - - if intensity' = intensity // The intensity doesn't change. - then - if m.Elements.Count + nextElements.Count + 1 > area - then - m.State <- AreaState.Validated - m.Intensity <- Some intensity - stop <- true - else - nextElements.Add(p) |> ignore - - elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity - then - m.Elements.UnionWith(nextElements) - for e in nextElements do - pixels.[e.Y, e.X] <- m - - if m.Elements.Count = area - then - m.State <- AreaState.Validated - m.Intensity <- Some (intensity') - stop <- true - else - intensity <- intensity' - nextElements.Clear() - nextElements.Add(p) |> ignore - - else - match pixels.[p.Y, p.X] with - | null -> () - | m' -> - if m'.Elements.Count + m.Elements.Count <= area - then - m'.State <- AreaState.Removed - for e in m'.Elements do - pixels.[e.Y, e.X] <- m - queue.Remove imgData.[e.Y, e.X, 0] e - addEdgeToQueue m'.Elements - m.Elements.UnionWith(m'.Elements) - let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min - if intensityMax <> intensity - then - intensity <- intensityMax - nextElements.Clear() - merged <- true - - if not merged - then - m.State <- AreaState.Validated - m.Intensity <- Some (intensity) - stop <- true - - if not stop && not merged - then - for i, j in se do - let ni = i + p.Y - let nj = j + p.X - let p' = Point(nj, ni) - if ni < 0 || ni >= h || nj < 0 || nj >= w - then - m.State <- AreaState.Validated - m.Intensity <- Some (intensity) - stop <- true - elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) - then - queue.Add (imgData.[ni, nj, 0]) p' - - if queue.IsEmpty - then - if m.Elements.Count + nextElements.Count <= area - then - m.State <- AreaState.Validated - m.Intensity <- Some intensity' - m.Elements.UnionWith(nextElements) - stop <- true - - for m in areas do - if m.State = AreaState.Validated - then - match m.Intensity with - | Some i -> - for p in m.Elements do - imgData.[p.Y, p.X, 0] <- i - | _ -> () - () - -/// -/// Area opening on byte image. -/// -let areaOpen (img: Image) (area: int) = - areaOperation img area AreaOperation.Opening - -/// -/// Area closing on byte image. -/// -let areaClose (img: Image) (area: int) = - areaOperation img area AreaOperation.Closing - -// A simpler algorithm than 'areaOpen' on byte image but slower. -let areaOpen2 (img: Image) (area: int) = - let w = img.Width - let h = img.Height - let imgData = img.Data - let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] - - let histogram = Array.zeroCreate 256 - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - let v = imgData.[i, j, 0] |> int - histogram.[v] <- histogram.[v] + 1 - - let flooded : bool[,] = Array2D.zeroCreate h w - - let pointsChecked = HashSet() - let pointsToCheck = Stack() - - for level in 255 .. -1 .. 0 do - let mutable n = histogram.[level] - if n > 0 - then - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - if not flooded.[i, j] && imgData.[i, j, 0] = byte level - then - let mutable maxNeighborValue = 0uy - pointsChecked.Clear() - pointsToCheck.Clear() - pointsToCheck.Push(Point(j, i)) - - while pointsToCheck.Count > 0 do - let next = pointsToCheck.Pop() - pointsChecked.Add(next) |> ignore - flooded.[next.Y, next.X] <- true - - for nx, ny in se do - let p = Point(next.X + nx, next.Y + ny) - if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h - then - let v = imgData.[p.Y, p.X, 0] - if v = byte level - then - if not (pointsChecked.Contains(p)) - then - pointsToCheck.Push(p) - elif v > maxNeighborValue - then - maxNeighborValue <- v - - if int maxNeighborValue < level && pointsChecked.Count <= area - then - for p in pointsChecked do - imgData.[p.Y, p.X, 0] <- maxNeighborValue - -[] -type Island (cmp: IComparer) = - member val Shore = Heap.Heap(cmp) with get - member val Level = 0.f with get, set - member val Surface = 0 with get, set - member this.IsInfinite = this.Surface = Int32.MaxValue - -let private areaOperationF (img: Image) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) = - let w = img.Width - let h = img.Height - let earth = img.Data - let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] - - let comparer = if op = AreaOperation.Opening - then { new IComparer with member this.Compare(v1, v2) = v1.CompareTo(v2) } - else { new IComparer with member this.Compare(v1, v2) = v2.CompareTo(v1) } - - let ownership: Island[,] = Array2D.create h w null - - // Initialize islands with their shore. - let islands = List() - let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima - for e in extremum do - let island = - let p = e.First() - Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count) - islands.Add(island) - let shorePoints = Points() - for p in e do - ownership.[p.Y, p.X] <- island - for i, j in se do - let ni = i + p.Y - let nj = j + p.X - let neighbor = Point(nj, ni) - if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor)) - then - shorePoints.Add(neighbor) |> ignore - island.Shore.Add earth.[ni, nj, 0] neighbor - - for area, obj in areas do - for island in islands do - let mutable stop = island.Shore.IsEmpty - - // 'true' if 'p' is owned or adjacent to 'island'. - let inline ownedOrAdjacent (p: Point) : bool = - ownership.[p.Y, p.X] = island || - (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) || - (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) || - (p.X > 0 && ownership.[p.Y, p.X - 1] = island) || - (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island) - - while not stop && island.Surface < area do - let level, next = island.Shore.Max - let other = ownership.[next.Y, next.X] - if other = island // During merging, some points on the shore may be owned by the island itself -> ignored. - then - island.Shore.RemoveNext () - else - if not <| Object.ReferenceEquals(other, null) - then // We touching another island. - if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0 - then - stop <- true - else // We can merge 'other' into 'surface'. - island.Surface <- island.Surface + other.Surface - island.Level <- other.Level - // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.Level - for l, p in other.Shore do - let mutable currentY = p.Y + 1 - while currentY < h && ownership.[currentY, p.X] = other do - ownership.[currentY, p.X] <- island - currentY <- currentY + 1 - island.Shore.Add l p - other.Shore.Clear() - - elif comparer.Compare(level, island.Level) > 0 - then - stop <- true - else - island.Shore.RemoveNext () - for i, j in se do - let ni = i + next.Y - let nj = j + next.X - if ni < 0 || ni >= h || nj < 0 || nj >= w - then - island.Surface <- Int32.MaxValue - stop <- true - else - let neighbor = Point(nj, ni) - if not <| ownedOrAdjacent neighbor - then - island.Shore.Add earth.[ni, nj, 0] neighbor - if not stop - then - ownership.[next.Y, next.X] <- island - island.Level <- level - island.Surface <- island.Surface + 1 - - let mutable diff = 0.f - - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - match ownership.[i, j] with - | null -> () - | island -> - let l = island.Level - diff <- diff + l - earth.[i, j, 0] - earth.[i, j, 0] <- l - - match f with - | Some f' -> f' obj diff - | _ -> () - () - -/// -/// Area opening on float image. -/// -let areaOpenF (img: Image) (area: int) = - areaOperationF img [ area, () ] None AreaOperation.Opening - -/// -/// Area closing on float image. -/// -let areaCloseF (img: Image) (area: int) = - areaOperationF img [ area, () ] None AreaOperation.Closing - -/// -/// Area closing on float image with different areas. Given areas must be sorted increasingly. -/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference -/// Between the previous and the current closing. -/// -let areaOpenFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = - areaOperationF img areas (Some f) AreaOperation.Opening - -/// -/// Same as 'areaOpenFWithFun' for closing operation. -/// -let areaCloseFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = - areaOperationF img areas (Some f) AreaOperation.Closing - -/// -/// Zhang and Suen thinning algorithm. -/// Modify 'mat' in place. -/// -let thin (mat: Matrix) = - let w = mat.Width - let h = mat.Height - let mutable data1 = mat.Data - let mutable data2 = Array2D.copy data1 - - let mutable pixelChanged = true - let mutable oddIteration = true - - while pixelChanged do - pixelChanged <- false - for i in 0..h-1 do - for j in 0..w-1 do - if data1.[i, j] = 1uy - then - let p2 = if i = 0 then 0uy else data1.[i-1, j] - let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1] - let p4 = if j = w-1 then 0uy else data1.[i, j+1] - let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1] - let p6 = if i = h-1 then 0uy else data1.[i+1, j] - let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1] - let p8 = if j = 0 then 0uy else data1.[i, j-1] - let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1] - - let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 - if sumNeighbors >= 2uy && sumNeighbors <= 6uy && - (if p2 = 0uy && p3 = 1uy then 1 else 0) + - (if p3 = 0uy && p4 = 1uy then 1 else 0) + - (if p4 = 0uy && p5 = 1uy then 1 else 0) + - (if p5 = 0uy && p6 = 1uy then 1 else 0) + - (if p6 = 0uy && p7 = 1uy then 1 else 0) + - (if p7 = 0uy && p8 = 1uy then 1 else 0) + - (if p8 = 0uy && p9 = 1uy then 1 else 0) + - (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 && - if oddIteration - then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy - else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy - then - data2.[i, j] <- 0uy - pixelChanged <- true - else - data2.[i, j] <- 0uy - - oddIteration <- not oddIteration - let tmp = data1 - data1 <- data2 - data2 <- tmp - -/// -/// Remove all 8-connected pixels with an area equal or greater than 'areaSize'. -/// Modify 'mat' in place. -/// -let removeArea (mat: Matrix) (areaSize: int) = - let neighbors = [| - (-1, 0) // p2 - (-1, 1) // p3 - ( 0, 1) // p4 - ( 1, 1) // p5 - ( 1, 0) // p6 - ( 1, -1) // p7 - ( 0, -1) // p8 - (-1, -1) |] // p9 - - use mat' = new Matrix(mat.Size) - let w = mat'.Width - let h = mat'.Height - 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 - then - let neighborhood = List() - let neighborsToCheck = Stack() - neighborsToCheck.Push(Point(j, i)) - data'.[i, j] <- 0uy - - while neighborsToCheck.Count > 0 do - let n = neighborsToCheck.Pop() - neighborhood.Add(n) - for (ni, nj) in neighbors do - let pi = n.Y + ni - let pj = n.X + nj - if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy - then - neighborsToCheck.Push(Point(pj, pi)) - data'.[pi, pj] <- 0uy - if neighborhood.Count <= areaSize - then - for n in neighborhood do - data.[n.Y, n.X] <- 0uy - -let connectedComponents (img: Image) (startPoints: List) : Points = - let w = img.Width - let h = img.Height - - let pointChecked = Points() - let pointToCheck = Stack(startPoints); - - let data = img.Data - - while pointToCheck.Count > 0 do - let next = pointToCheck.Pop() - 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.Push(p) - - pointChecked - -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: Ellipse) (color: 'TColor) (alpha: float) = - if alpha >= 1.0 - then - img.Draw(Emgu.CV.Structure.Ellipse(PointF(e.Cx, e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) - else - let windowPosX = e.Cx - e.A - 5.f - let gapX = windowPosX - (float32 (int windowPosX)) - - let windowPosY = e.Cy - e.A - 5.f - let gapY = windowPosY - (float32 (int windowPosY)) - - let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int) - - img.ROI <- roi - if roi = img.ROI // We do not display ellipses touching the edges (FIXME) - then - use i = new Image<'TColor, 'TDepth>(img.ROI.Size) - i.Draw(Emgu.CV.Structure.Ellipse(PointF(e.A + 5.f + gapX, e.A + 5.f + gapY), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) - CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img) - img.ROI <- Rectangle.Empty - -let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Ellipse list) (color: 'TColor) (alpha: float) = - List.iter (fun e -> drawEllipse img e color alpha) ellipses - -let rngCell = System.Random() -let drawCell (img: Image) (drawCellContent: bool) (c: Cell) = - if drawCellContent - then - let colorB = rngCell.Next(20, 70) - let colorG = rngCell.Next(20, 70) - let colorR = rngCell.Next(20, 70) - - for y in 0 .. c.elements.Height - 1 do - for x in 0 .. c.elements.Width - 1 do - if c.elements.[y, x] > 0uy - then - let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2 - let b = img.Data.[y + dy, x + dx, 0] |> int - let g = img.Data.[y + dy, x + dx, 1] |> int - let r = img.Data.[y + dy, x + dx, 2] |> int - img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB) - img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG) - img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR) - - let crossColor, crossColor2 = - match c.cellClass with - | HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.) - | InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.) - | Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.) - - drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2 - drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2 - - drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1 - drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1 - - -let drawCells (img: Image) (drawCellContent: bool) (cells: Cell list) = - List.iter (fun c -> drawCell img drawCellContent c) cells \ No newline at end of file diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs b/Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs new file mode 100644 index 0000000..363acda --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs @@ -0,0 +1,79 @@ +module ParasitemiaCore.Drawing + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +open Const +open Types + +let drawPoints (img: Image) (points: Points) (intensity: 'TDepth) = + for p in points do + img.Data.[p.Y, p.X, 0] <- intensity + +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: Ellipse) (color: 'TColor) (alpha: float) = + if alpha >= 1.0 + then + img.Draw(Emgu.CV.Structure.Ellipse(PointF(e.Cx, e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) + else + let windowPosX = e.Cx - e.A - 5.f + let gapX = windowPosX - (float32 (int windowPosX)) + + let windowPosY = e.Cy - e.A - 5.f + let gapY = windowPosY - (float32 (int windowPosY)) + + let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int) + + img.ROI <- roi + if roi = img.ROI // We do not display ellipses touching the edges (FIXME) + then + use i = new Image<'TColor, 'TDepth>(img.ROI.Size) + i.Draw(Emgu.CV.Structure.Ellipse(PointF(e.A + 5.f + gapX, e.A + 5.f + gapY), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias) + CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img) + img.ROI <- Rectangle.Empty + +let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Ellipse list) (color: 'TColor) (alpha: float) = + List.iter (fun e -> drawEllipse img e color alpha) ellipses + +let rngCell = System.Random() +let drawCell (img: Image) (drawCellContent: bool) (c: Cell) = + if drawCellContent + then + let colorB = rngCell.Next(20, 70) + let colorG = rngCell.Next(20, 70) + let colorR = rngCell.Next(20, 70) + + for y in 0 .. c.elements.Height - 1 do + for x in 0 .. c.elements.Width - 1 do + if c.elements.[y, x] > 0uy + then + let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2 + let b = img.Data.[y + dy, x + dx, 0] |> int + let g = img.Data.[y + dy, x + dx, 1] |> int + let r = img.Data.[y + dy, x + dx, 2] |> int + img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB) + img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG) + img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR) + + let crossColor, crossColor2 = + match c.cellClass with + | HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.) + | InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.) + | Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.) + + drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2 + drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2 + + drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1 + drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1 + +let drawCells (img: Image) (drawCellContent: bool) (cells: Cell list) = + List.iter (fun c -> drawCell img drawCellContent c) cells \ No newline at end of file diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Edges.fs b/Parasitemia/ParasitemiaCore/ImgTools/Edges.fs new file mode 100644 index 0000000..b174ee9 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Edges.fs @@ -0,0 +1,121 @@ +module ParasitemiaCore.Edges + +open System +open System.Drawing +open System.Collections.Generic + +open Emgu.CV +open Emgu.CV.Structure + +open Const +open Histogram +open Otsu + +/// +/// Find edges of an image by using the Canny approach. +/// The thresholds are automatically defined with otsu on gradient magnitudes. +/// +/// +let find (img: Image) : Matrix * Matrix * Matrix = + let w = img.Width + let h = img.Height + + use sobelKernel = + new Matrix(array2D [[ 1.0f; 0.0f; -1.0f ] + [ 2.0f; 0.0f; -2.0f ] + [ 1.0f; 0.0f; -1.0f ]]) + + let xGradient = new Matrix(img.Size) + let yGradient = new Matrix(img.Size) + CvInvoke.Filter2D(img, xGradient, sobelKernel, Point(1, 1)) + CvInvoke.Filter2D(img, yGradient, sobelKernel.Transpose(), Point(1, 1)) + + use magnitudes = new Matrix(xGradient.Size) + use angles = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. + + let thresholdHigh, thresholdLow = + let sensibilityHigh = 0.1f + let sensibilityLow = 0.0f + let threshold, _, _ = otsu (histogramMat magnitudes 300) + threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold) + + // Non-maximum suppression. + use nms = new Matrix(xGradient.Size) + + let nmsData = nms.Data + let anglesData = angles.Data + let magnitudesData = magnitudes.Data + let xGradientData = xGradient.Data + let yGradientData = yGradient.Data + + for i in 0 .. h - 1 do + nmsData.[i, 0] <- 0uy + nmsData.[i, w - 1] <- 0uy + + for j in 0 .. w - 1 do + nmsData.[0, j] <- 0uy + nmsData.[h - 1, j] <- 0uy + + for i in 1 .. h - 2 do + for j in 1 .. w - 2 do + let vx = xGradientData.[i, j] + let vy = yGradientData.[i, j] + if vx <> 0.f || vy <> 0.f + then + let angle = anglesData.[i, j] + + let vx', vy' = abs vx, abs vy + let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy' + let ratio1 = 1.f - ratio2 + + let mNeigbors (sign: int) : float32 = + if angle < PI / 4.f + then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign] + elif angle < PI / 2.f + then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j] + elif angle < 3.f * PI / 4.f + then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign] + elif angle < PI + then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign] + elif angle < 5.f * PI / 4.f + then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign] + elif angle < 3.f * PI / 2.f + then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j] + elif angle < 7.f * PI / 4.f + then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign] + else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign] + + let m = magnitudesData.[i, j] + if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 + then + nmsData.[i, j] <- 1uy + + // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection). + + let edges = new Matrix(xGradient.Size) + let edgesData = edges.Data + + // Hysteresis thresholding. + let toVisit = Stack() + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh + then + nmsData.[i, j] <- 0uy + toVisit.Push(Point(j, i)) + while toVisit.Count > 0 do + let p = toVisit.Pop() + edgesData.[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 && nmsData.[ni, nj] = 1uy + then + nmsData.[ni, nj] <- 0uy + toVisit.Push(Point(nj, ni)) + + edges, xGradient, yGradient \ No newline at end of file diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs b/Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs new file mode 100644 index 0000000..92b73f6 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs @@ -0,0 +1,84 @@ +module ParasitemiaCore.Histogram + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +type Histogram = { + data: int[] + total: int // Number of elements. + sum: int // Sum of all intensity. + min: float32 + max: float32 } + +let histogramImg (img: Image) (nbSamples: int) : Histogram = + let imgData = img.Data + + let min, max = + let min = ref [| 0.0 |] + let minLocation = ref <| [| Point() |] + let max = ref [| 0.0 |] + let maxLocation = ref <| [| Point() |] + img.MinMax(min, max, minLocation, maxLocation) + float32 (!min).[0], float32 (!max).[0] + + let inline bin (x: float32) : int = + let p = int ((x - min) / (max - min) * float32 nbSamples) + if p >= nbSamples then nbSamples - 1 else p + + let data = Array.zeroCreate nbSamples + + for i in 0 .. img.Height - 1 do + for j in 0 .. img.Width - 1 do + let p = bin imgData.[i, j, 0] + data.[p] <- data.[p] + 1 + + { data = data; total = img.Height * img.Width; sum = Array.sum data; min = min; max = max } + +let histogramMat (mat: Matrix) (nbSamples: int) : Histogram = + let matData = mat.Data + + let min, max = + let min = ref 0.0 + let minLocation = ref <| Point() + let max = ref 0.0 + let maxLocation = ref <| Point() + mat.MinMax(min, max, minLocation, maxLocation) + float32 !min, float32 !max + + let inline bin (x: float32) : int = + let p = int ((x - min) / (max - min) * float32 nbSamples) + if p >= nbSamples then nbSamples - 1 else p + + let data = Array.zeroCreate nbSamples + + for i in 0 .. mat.Height - 1 do + for j in 0 .. mat.Width - 1 do + let p = bin matData.[i, j] + data.[p] <- data.[p] + 1 + + { data = data; total = mat.Height * mat.Width; sum = Array.sum data; min = min; max = max } + +let histogram (values: float32 seq) (nbSamples: int) : Histogram = + let mutable min = Single.MaxValue + let mutable max = Single.MinValue + let mutable n = 0 + + for v in values do + n <- n + 1 + if v < min then min <- v + if v > max then max <- v + + let inline bin (x: float32) : int = + let p = int ((x - min) / (max - min) * float32 nbSamples) + if p >= nbSamples then nbSamples - 1 else p + + let data = Array.zeroCreate nbSamples + + for v in values do + let p = bin v + data.[p] <- data.[p] + 1 + + { data = data; total = n; sum = Array.sum data; min = min; max = max } \ No newline at end of file diff --git a/Parasitemia/ParasitemiaCore/ImgTools/IO.fs b/Parasitemia/ParasitemiaCore/ImgTools/IO.fs new file mode 100644 index 0000000..219eab6 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/IO.fs @@ -0,0 +1,15 @@ +module ParasitemiaCore.IO + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = + img.Save(filepath) + +let saveMat (mat: Matrix<'TDepth>) (filepath: string) = + use img = new Image(mat.Size) + mat.CopyTo(img) + saveImg img filepath diff --git a/Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs b/Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs new file mode 100644 index 0000000..0510bc8 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs @@ -0,0 +1,49 @@ +module ParasitemiaCore.ImgTools + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +let normalize (img: Image) (upperLimit: float) : Image = + let min = ref [| 0.0 |] + let minLocation = ref <| [| Point() |] + let max = ref [| 0.0 |] + let maxLocation = ref <| [| Point() |] + img.MinMax(min, max, minLocation, maxLocation) + let normalized = (img - (!min).[0]) / ((!max).[0] - (!min).[0]) + if upperLimit = 1.0 + then normalized + else upperLimit * normalized + +let mergeChannels (img: Image) (rgbWeights: float * float * float) : Image = + match rgbWeights with + | 1., 0., 0. -> img.[2] + | 0., 1., 0. -> img.[1] + | 0., 0., 1. -> img.[0] + | redFactor, greenFactor, blueFactor -> + let result = new Image(img.Size) + CvInvoke.AddWeighted(result, 1., img.[2], redFactor, 0., result) + CvInvoke.AddWeighted(result, 1., img.[1], greenFactor, 0., result) + CvInvoke.AddWeighted(result, 1., img.[0], blueFactor, 0., result) + result + +let mergeChannelsWithProjection (img: Image) (v1r: float32, v1g: float32, v1b: float32) (v2r: float32, v2g: float32, v2b: float32) (upperLimit: float) : Image = + let vr, vg, vb = v2r - v1r, v2g - v1g, v2b - v1b + let vMagnitude = sqrt (vr ** 2.f + vg ** 2.f + vb ** 2.f) + let project (r: float32) (g: float32) (b: float32) = ((r - v1r) * vr + (g - v1g) * vg + (b - v1b) * vb) / vMagnitude + let result = new Image(img.Size) + // TODO: Essayer en bindant Data pour gagner du temps + for i in 0 .. img.Height - 1 do + for j in 0 .. img.Width - 1 do + result.Data.[i, j, 0] <- project img.Data.[i, j, 2] img.Data.[i, j, 1] img.Data.[i, j, 0] + normalize result upperLimit + +// Normalize image values between 0uy and 255uy. +let normalizeAndConvert (img: Image) : Image = + (normalize (img.Convert()) 255.).Convert() + +let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> = + let size = 2 * int (ceil (4.0 * standardDeviation)) + 1 + img.SmoothGaussian(size, size, standardDeviation, standardDeviation) diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs b/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs new file mode 100644 index 0000000..b3a2e75 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs @@ -0,0 +1,670 @@ +module ParasitemiaCore.Morpho + +open System +open System.Drawing +open System.Collections.Generic +open System.Linq + +open Emgu.CV +open Emgu.CV.Structure + +open Types + +/// +/// Remove M-adjacent pixels. It may be used after thinning. +/// +let suppressMAdjacency (img: Matrix) = + 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 + +type ExtremumType = + | Maxima = 1 + | Minima = 2 + +let findExtremum (img: Image) (extremumType: ExtremumType) : IEnumerable = + let w = img.Width + let h = img.Height + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let imgData = img.Data + let suppress: bool[,] = Array2D.zeroCreate h w + + let result = List>() + + let flood (start: Point) : List> = + let sameLevelToCheck = Stack() + let betterLevelToCheck = Stack() + betterLevelToCheck.Push(start) + + let result' = List>() + + while betterLevelToCheck.Count > 0 do + let p = betterLevelToCheck.Pop() + if not suppress.[p.Y, p.X] + then + suppress.[p.Y, p.X] <- true + sameLevelToCheck.Push(p) + let current = List() + + let mutable betterExists = false + + while sameLevelToCheck.Count > 0 do + let p' = sameLevelToCheck.Pop() + let currentLevel = imgData.[p'.Y, p'.X, 0] + current.Add(p') |> ignore + for i, j in se do + let ni = i + p'.Y + let nj = j + p'.X + if ni >= 0 && ni < h && nj >= 0 && nj < w + then + let level = imgData.[ni, nj, 0] + let notSuppressed = not suppress.[ni, nj] + + if level = currentLevel && notSuppressed + then + suppress.[ni, nj] <- true + sameLevelToCheck.Push(Point(nj, ni)) + elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel + then + betterExists <- true + if notSuppressed + then + betterLevelToCheck.Push(Point(nj, ni)) + + if not betterExists + then + result'.Add(current) + result' + + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + let maxima = flood (Point(j, i)) + if maxima.Count > 0 + then + result.AddRange(maxima) + + result.Select(fun l -> Points(l)) + +let findMaxima (img: Image) : IEnumerable = + findExtremum img ExtremumType.Maxima + +let findMinima (img: Image) : IEnumerable = + findExtremum img ExtremumType.Minima + +type PriorityQueue () = + let size = 256 + let q: Points[] = Array.init size (fun i -> Points()) + let mutable highest = -1 // Value of the first elements of 'q'. + let mutable lowest = size + + member this.NextMax () : byte * Point = + if this.IsEmpty + then + invalidOp "Queue is empty" + else + let l = q.[highest] + let next = l.First() + l.Remove(next) |> ignore + let value = byte highest + + if l.Count = 0 + then + highest <- highest - 1 + while highest > lowest && q.[highest].Count = 0 do + highest <- highest - 1 + if highest = lowest + then + highest <- -1 + lowest <- size + + value, next + + member this.NextMin () : byte * Point = + if this.IsEmpty + then + invalidOp "Queue is empty" + else + let l = q.[lowest + 1] + let next = l.First() + l.Remove(next) |> ignore + let value = byte (lowest + 1) + + if l.Count = 0 + then + lowest <- lowest + 1 + while lowest < highest && q.[lowest + 1].Count = 0 do + lowest <- lowest + 1 + if highest = lowest + then + highest <- -1 + lowest <- size + + value, next + + member this.Max = + highest |> byte + + member this.Min = + lowest + 1 |> byte + + member this.Add (value: byte) (p: Point) = + let vi = int value + + if vi > highest + then + highest <- vi + if vi <= lowest + then + lowest <- vi - 1 + + q.[vi].Add(p) |> ignore + + member this.Remove (value: byte) (p: Point) = + let vi = int value + if q.[vi].Remove(p) && q.[vi].Count = 0 + then + if vi = highest + then + highest <- highest - 1 + while highest > lowest && q.[highest].Count = 0 do + highest <- highest - 1 + elif vi - 1 = lowest + then + lowest <- lowest + 1 + while lowest < highest && q.[lowest + 1].Count = 0 do + lowest <- lowest + 1 + + if highest = lowest // The queue is now empty. + then + highest <- -1 + lowest <- size + + member this.IsEmpty = + highest = -1 + + member this.Clear () = + while highest > lowest do + q.[highest].Clear() + highest <- highest - 1 + highest <- -1 + lowest <- size + +type private AreaState = + | Removed = 1 + | Unprocessed = 2 + | Validated = 3 + +type private AreaOperation = + | Opening = 1 + | Closing = 2 + +[] +type private Area (elements: Points) = + member this.Elements = elements + member val Intensity = None with get, set + member val State = AreaState.Unprocessed with get, set + +let private areaOperation (img: Image) (area: int) (op: AreaOperation) = + let w = img.Width + let h = img.Height + let imgData = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let areas = List((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area) + + let pixels: Area[,] = Array2D.create h w null + for m in areas do + for e in m.Elements do + pixels.[e.Y, e.X] <- m + + let queue = PriorityQueue() + + let addEdgeToQueue (elements: Points) = + for p in elements do + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let p' = Point(nj, ni) + if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p')) + then + queue.Add (imgData.[ni, nj, 0]) p' + + // Reverse order is quicker. + for i in areas.Count - 1 .. -1 .. 0 do + let m = areas.[i] + if m.Elements.Count <= area && m.State <> AreaState.Removed + then + queue.Clear() + addEdgeToQueue m.Elements + + let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min + let nextElements = Points() + + let mutable stop = false + while not stop do + let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin () + let mutable merged = false + + if intensity' = intensity // The intensity doesn't change. + then + if m.Elements.Count + nextElements.Count + 1 > area + then + m.State <- AreaState.Validated + m.Intensity <- Some intensity + stop <- true + else + nextElements.Add(p) |> ignore + + elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity + then + m.Elements.UnionWith(nextElements) + for e in nextElements do + pixels.[e.Y, e.X] <- m + + if m.Elements.Count = area + then + m.State <- AreaState.Validated + m.Intensity <- Some (intensity') + stop <- true + else + intensity <- intensity' + nextElements.Clear() + nextElements.Add(p) |> ignore + + else + match pixels.[p.Y, p.X] with + | null -> () + | m' -> + if m'.Elements.Count + m.Elements.Count <= area + then + m'.State <- AreaState.Removed + for e in m'.Elements do + pixels.[e.Y, e.X] <- m + queue.Remove imgData.[e.Y, e.X, 0] e + addEdgeToQueue m'.Elements + m.Elements.UnionWith(m'.Elements) + let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min + if intensityMax <> intensity + then + intensity <- intensityMax + nextElements.Clear() + merged <- true + + if not merged + then + m.State <- AreaState.Validated + m.Intensity <- Some (intensity) + stop <- true + + if not stop && not merged + then + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let p' = Point(nj, ni) + if ni < 0 || ni >= h || nj < 0 || nj >= w + then + m.State <- AreaState.Validated + m.Intensity <- Some (intensity) + stop <- true + elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p')) + then + queue.Add (imgData.[ni, nj, 0]) p' + + if queue.IsEmpty + then + if m.Elements.Count + nextElements.Count <= area + then + m.State <- AreaState.Validated + m.Intensity <- Some intensity' + m.Elements.UnionWith(nextElements) + stop <- true + + for m in areas do + if m.State = AreaState.Validated + then + match m.Intensity with + | Some i -> + for p in m.Elements do + imgData.[p.Y, p.X, 0] <- i + | _ -> () + () + +/// +/// Area opening on byte image. +/// +let areaOpen (img: Image) (area: int) = + areaOperation img area AreaOperation.Opening + +/// +/// Area closing on byte image. +/// +let areaClose (img: Image) (area: int) = + areaOperation img area AreaOperation.Closing + +// A simpler algorithm than 'areaOpen' on byte image but slower. +let areaOpen2 (img: Image) (area: int) = + let w = img.Width + let h = img.Height + let imgData = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let histogram = Array.zeroCreate 256 + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + let v = imgData.[i, j, 0] |> int + histogram.[v] <- histogram.[v] + 1 + + let flooded : bool[,] = Array2D.zeroCreate h w + + let pointsChecked = HashSet() + let pointsToCheck = Stack() + + for level in 255 .. -1 .. 0 do + let mutable n = histogram.[level] + if n > 0 + then + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + if not flooded.[i, j] && imgData.[i, j, 0] = byte level + then + let mutable maxNeighborValue = 0uy + pointsChecked.Clear() + pointsToCheck.Clear() + pointsToCheck.Push(Point(j, i)) + + while pointsToCheck.Count > 0 do + let next = pointsToCheck.Pop() + pointsChecked.Add(next) |> ignore + flooded.[next.Y, next.X] <- true + + for nx, ny in se do + let p = Point(next.X + nx, next.Y + ny) + if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h + then + let v = imgData.[p.Y, p.X, 0] + if v = byte level + then + if not (pointsChecked.Contains(p)) + then + pointsToCheck.Push(p) + elif v > maxNeighborValue + then + maxNeighborValue <- v + + if int maxNeighborValue < level && pointsChecked.Count <= area + then + for p in pointsChecked do + imgData.[p.Y, p.X, 0] <- maxNeighborValue + +[] +type Island (cmp: IComparer) = + member val Shore = Heap.Heap(cmp) with get + member val Level = 0.f with get, set + member val Surface = 0 with get, set + member this.IsInfinite = this.Surface = Int32.MaxValue + +let private areaOperationF (img: Image) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) = + let w = img.Width + let h = img.Height + let earth = img.Data + let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |] + + let comparer = if op = AreaOperation.Opening + then { new IComparer with member this.Compare(v1, v2) = v1.CompareTo(v2) } + else { new IComparer with member this.Compare(v1, v2) = v2.CompareTo(v1) } + + let ownership: Island[,] = Array2D.create h w null + + // Initialize islands with their shore. + let islands = List() + let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima + for e in extremum do + let island = + let p = e.First() + Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count) + islands.Add(island) + let shorePoints = Points() + for p in e do + ownership.[p.Y, p.X] <- island + for i, j in se do + let ni = i + p.Y + let nj = j + p.X + let neighbor = Point(nj, ni) + if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor)) + then + shorePoints.Add(neighbor) |> ignore + island.Shore.Add earth.[ni, nj, 0] neighbor + + for area, obj in areas do + for island in islands do + let mutable stop = island.Shore.IsEmpty + + // 'true' if 'p' is owned or adjacent to 'island'. + let inline ownedOrAdjacent (p: Point) : bool = + ownership.[p.Y, p.X] = island || + (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) || + (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) || + (p.X > 0 && ownership.[p.Y, p.X - 1] = island) || + (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island) + + while not stop && island.Surface < area do + let level, next = island.Shore.Max + let other = ownership.[next.Y, next.X] + if other = island // During merging, some points on the shore may be owned by the island itself -> ignored. + then + island.Shore.RemoveNext () + else + if not <| Object.ReferenceEquals(other, null) + then // We touching another island. + if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0 + then + stop <- true + else // We can merge 'other' into 'surface'. + island.Surface <- island.Surface + other.Surface + island.Level <- other.Level + // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.Level + for l, p in other.Shore do + let mutable currentY = p.Y + 1 + while currentY < h && ownership.[currentY, p.X] = other do + ownership.[currentY, p.X] <- island + currentY <- currentY + 1 + island.Shore.Add l p + other.Shore.Clear() + + elif comparer.Compare(level, island.Level) > 0 + then + stop <- true + else + island.Shore.RemoveNext () + for i, j in se do + let ni = i + next.Y + let nj = j + next.X + if ni < 0 || ni >= h || nj < 0 || nj >= w + then + island.Surface <- Int32.MaxValue + stop <- true + else + let neighbor = Point(nj, ni) + if not <| ownedOrAdjacent neighbor + then + island.Shore.Add earth.[ni, nj, 0] neighbor + if not stop + then + ownership.[next.Y, next.X] <- island + island.Level <- level + island.Surface <- island.Surface + 1 + + let mutable diff = 0.f + + for i in 0 .. h - 1 do + for j in 0 .. w - 1 do + match ownership.[i, j] with + | null -> () + | island -> + let l = island.Level + diff <- diff + l - earth.[i, j, 0] + earth.[i, j, 0] <- l + + match f with + | Some f' -> f' obj diff + | _ -> () + () + +/// +/// Area opening on float image. +/// +let areaOpenF (img: Image) (area: int) = + areaOperationF img [ area, () ] None AreaOperation.Opening + +/// +/// Area closing on float image. +/// +let areaCloseF (img: Image) (area: int) = + areaOperationF img [ area, () ] None AreaOperation.Closing + +/// +/// Area closing on float image with different areas. Given areas must be sorted increasingly. +/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference +/// Between the previous and the current closing. +/// +let areaOpenFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Opening + +/// +/// Same as 'areaOpenFWithFun' for closing operation. +/// +let areaCloseFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Closing + +/// +/// Zhang and Suen thinning algorithm. +/// Modify 'mat' in place. +/// +let thin (mat: Matrix) = + let w = mat.Width + let h = mat.Height + let mutable data1 = mat.Data + let mutable data2 = Array2D.copy data1 + + let mutable pixelChanged = true + let mutable oddIteration = true + + while pixelChanged do + pixelChanged <- false + for i in 0..h-1 do + for j in 0..w-1 do + if data1.[i, j] = 1uy + then + let p2 = if i = 0 then 0uy else data1.[i-1, j] + let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1] + let p4 = if j = w-1 then 0uy else data1.[i, j+1] + let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1] + let p6 = if i = h-1 then 0uy else data1.[i+1, j] + let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1] + let p8 = if j = 0 then 0uy else data1.[i, j-1] + let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1] + + let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + if sumNeighbors >= 2uy && sumNeighbors <= 6uy && + (if p2 = 0uy && p3 = 1uy then 1 else 0) + + (if p3 = 0uy && p4 = 1uy then 1 else 0) + + (if p4 = 0uy && p5 = 1uy then 1 else 0) + + (if p5 = 0uy && p6 = 1uy then 1 else 0) + + (if p6 = 0uy && p7 = 1uy then 1 else 0) + + (if p7 = 0uy && p8 = 1uy then 1 else 0) + + (if p8 = 0uy && p9 = 1uy then 1 else 0) + + (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 && + if oddIteration + then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy + else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy + then + data2.[i, j] <- 0uy + pixelChanged <- true + else + data2.[i, j] <- 0uy + + oddIteration <- not oddIteration + let tmp = data1 + data1 <- data2 + data2 <- tmp + +/// +/// Remove all 8-connected pixels with an area equal or greater than 'areaSize'. +/// Modify 'mat' in place. +/// +let removeArea (mat: Matrix) (areaSize: int) = + let neighbors = [| + (-1, 0) // p2 + (-1, 1) // p3 + ( 0, 1) // p4 + ( 1, 1) // p5 + ( 1, 0) // p6 + ( 1, -1) // p7 + ( 0, -1) // p8 + (-1, -1) |] // p9 + + use mat' = new Matrix(mat.Size) + let w = mat'.Width + let h = mat'.Height + 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 + then + let neighborhood = List() + let neighborsToCheck = Stack() + neighborsToCheck.Push(Point(j, i)) + data'.[i, j] <- 0uy + + while neighborsToCheck.Count > 0 do + let n = neighborsToCheck.Pop() + neighborhood.Add(n) + for (ni, nj) in neighbors do + let pi = n.Y + ni + let pj = n.X + nj + if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy + then + neighborsToCheck.Push(Point(pj, pi)) + data'.[pi, pj] <- 0uy + if neighborhood.Count <= areaSize + then + for n in neighborhood do + data.[n.Y, n.X] <- 0uy + +let connectedComponents (img: Image) (startPoints: List) : Points = + let w = img.Width + let h = img.Height + + let pointChecked = Points() + let pointToCheck = Stack(startPoints); + + let data = img.Data + + while pointToCheck.Count > 0 do + let next = pointToCheck.Pop() + 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.Push(p) + + pointChecked diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs b/Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs new file mode 100644 index 0000000..6b8ee53 --- /dev/null +++ b/Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs @@ -0,0 +1,47 @@ +module ParasitemiaCore.Otsu + +open Histogram + +let otsu (hist: Histogram) : float32 * float32 * float32 = + let mutable sumB = 0 + let mutable wB = 0 + let mutable maximum = 0.0 + let mutable level = 0 + let sum = hist.data |> Array.mapi (fun i v -> i * v |> float) |> Array.sum + + for i in 0 .. hist.data.Length - 1 do + wB <- wB + hist.data.[i] + if wB <> 0 + then + let wF = hist.total - wB + if wF <> 0 + then + sumB <- sumB + i * hist.data.[i] + let mB = (float sumB) / (float wB) + let mF = (sum - float sumB) / (float wF) + let between = (float wB) * (float wF) * (mB - mF) ** 2.; + if between >= maximum + then + level <- i + maximum <- between + + let mean1 = + let mutable sum = 0 + let mutable nb = 0 + for i in 0 .. level - 1 do + sum <- sum + i * hist.data.[i] + nb <- nb + hist.data.[i] + (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) + + let mean2 = + let mutable sum = 0 + let mutable nb = 0 + for i in level + 1 .. hist.data.Length - 1 do + sum <- sum + i * hist.data.[i] + nb <- nb + hist.data.[i] + (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2) + + let toFloat l = + float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min + + toFloat level, toFloat mean1, toFloat mean2 \ No newline at end of file diff --git a/Parasitemia/ParasitemiaCore/MainAnalysis.fs b/Parasitemia/ParasitemiaCore/MainAnalysis.fs index 9c85c7b..1d2ced3 100644 --- a/Parasitemia/ParasitemiaCore/MainAnalysis.fs +++ b/Parasitemia/ParasitemiaCore/MainAnalysis.fs @@ -12,6 +12,7 @@ open Emgu.CV.Structure open Logger open Utils +open Morpho open ImgTools open Config open Types @@ -88,7 +89,7 @@ let doAnalysis (img: Image) (name: string) (config: Config) (reportPr logTimeWithName "Parasites segmentation" (fun () -> reportWithVal 40 (ParasitesMarker.find img_parasites_filtered config)) let! edges, xGradient, yGradient = logTimeWithName "Finding edges" (fun () -> - let edges, xGradient, yGradient = findEdges img_RBC_filtered + let edges, xGradient, yGradient = Edges.find img_RBC_filtered removeArea edges (config.RBCRadius.Pixel ** 2.f / 50.f |> int) reportWithVal 50 (edges, xGradient, yGradient)) @@ -109,43 +110,43 @@ let doAnalysis (img: Image) (name: string) (config: Config) (reportPr let buildFileName postfix = System.IO.Path.Combine(dirPath, name + postfix) - saveMat (edges * 255.0) (buildFileName " - edges.png") + IO.saveMat (edges * 255.0) (buildFileName " - edges.png") - saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png") - saveImg parasites.parasite (buildFileName " - parasites - stain.png") - saveImg parasites.nucleus (buildFileName " - parasites - infection.png") + IO.saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png") + IO.saveImg parasites.parasite (buildFileName " - parasites - stain.png") + IO.saveImg parasites.nucleus (buildFileName " - parasites - infection.png") let imgAllEllipses = img.Copy() - drawEllipses imgAllEllipses matchingEllipses.Ellipses (Bgr(255.0, 255.0, 255.0)) 0.04 - saveImg imgAllEllipses (buildFileName " - ellipses - all.png") + Drawing.drawEllipses imgAllEllipses matchingEllipses.Ellipses (Bgr(255.0, 255.0, 255.0)) 0.04 + IO.saveImg imgAllEllipses (buildFileName " - ellipses - all.png") let imgEllipses = img_RBC_filtered.Convert() - drawEllipses imgEllipses prunedEllipses (Bgr(0.0, 240.0, 240.0)) 1.0 - saveImg imgEllipses (buildFileName " - ellipses.png") + Drawing.drawEllipses imgEllipses prunedEllipses (Bgr(0.0, 240.0, 240.0)) 1.0 + IO.saveImg imgEllipses (buildFileName " - ellipses.png") let imgCells = img.Copy() - drawCells imgCells false cells - saveImg imgCells (buildFileName " - cells.png") + Drawing.drawCells imgCells false cells + IO.saveImg imgCells (buildFileName " - cells.png") let imgCells' = img.Copy() - drawCells imgCells' true cells - saveImg imgCells' (buildFileName " - cells - full.png") + Drawing.drawCells imgCells' true cells + IO.saveImg imgCells' (buildFileName " - cells - full.png") let filteredGreenMaxima = gaussianFilter img_RBC config.LPFStandardDeviationRBC for m in findMaxima filteredGreenMaxima do - drawPoints filteredGreenMaxima m 255.f - saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png") + Drawing.drawPoints filteredGreenMaxima m 255.f + IO.saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png") - saveImg img_RBC_filtered (buildFileName " - filtered.png") - saveImg imgWhitoutParasite (buildFileName " - filtered closed stain.png") - saveImg imgWithoutNucleus (buildFileName " - filtered closed infection.png") + IO.saveImg img_RBC_filtered (buildFileName " - filtered.png") + IO.saveImg imgWhitoutParasite (buildFileName " - filtered closed stain.png") + IO.saveImg imgWithoutNucleus (buildFileName " - filtered closed infection.png") - saveImg img_RBC (buildFileName " - source - RBC.png") - saveImg img_parasites (buildFileName " - source - parasites.png") + IO.saveImg img_RBC (buildFileName " - source - RBC.png") + IO.saveImg img_parasites (buildFileName " - source - parasites.png") - saveImg (normalize img_float.[2] 255.) (buildFileName " - source - red.png") - saveImg (normalize img_float.[1] 255.) (buildFileName " - source - green.png") - saveImg (normalize img_float.[0] 255.) (buildFileName " - source - blue.png") + IO.saveImg (normalize img_float.[2] 255.) (buildFileName " - source - red.png") + IO.saveImg (normalize img_float.[1] 255.) (buildFileName " - source - green.png") + IO.saveImg (normalize img_float.[0] 255.) (buildFileName " - source - blue.png") | _ -> () return cells } diff --git a/Parasitemia/ParasitemiaCore/ParasitemiaCore.fsproj b/Parasitemia/ParasitemiaCore/ParasitemiaCore.fsproj index 58886c0..6a9eb75 100644 --- a/Parasitemia/ParasitemiaCore/ParasitemiaCore.fsproj +++ b/Parasitemia/ParasitemiaCore/ParasitemiaCore.fsproj @@ -51,7 +51,13 @@ - + + + + + + + diff --git a/Parasitemia/ParasitemiaCore/ParasitesMarker.fs b/Parasitemia/ParasitemiaCore/ParasitesMarker.fs index 6c25d04..0d7c7ae 100644 --- a/Parasitemia/ParasitemiaCore/ParasitesMarker.fs +++ b/Parasitemia/ParasitemiaCore/ParasitesMarker.fs @@ -7,6 +7,9 @@ open Emgu.CV open Emgu.CV.Structure open Utils +open Histogram +open Otsu +open Morpho open ImgTools type Result = { diff --git a/Parasitemia/ParasitemiaCore/UnitsOfMeasure.fs b/Parasitemia/ParasitemiaCore/UnitsOfMeasure.fs index 0d84d16..36c9495 100644 --- a/Parasitemia/ParasitemiaCore/UnitsOfMeasure.fs +++ b/Parasitemia/ParasitemiaCore/UnitsOfMeasure.fs @@ -1,9 +1,9 @@ module ParasitemiaCore.UnitsOfMeasure -[] type px +[] type px // Pixel. [] type μm [] type inch -[] type ppi = px / inch +[] type ppi = px / inch // Pixel per inch. let μmInchRatio = 25.4e3<μm/inch> diff --git a/Parasitemia/ParasitemiaUI/XAML/AboutWindow.xaml b/Parasitemia/ParasitemiaUI/XAML/AboutWindow.xaml index 3cb0de4..e35d234 100644 --- a/Parasitemia/ParasitemiaUI/XAML/AboutWindow.xaml +++ b/Parasitemia/ParasitemiaUI/XAML/AboutWindow.xaml @@ -3,7 +3,7 @@ xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" mc:Ignorable="d" - x:Name="AboutWindow" Height="200.969" Width="282.313" MinHeight="100" MinWidth="100" Title="About" Icon="pack://application:,,,/Resources/icon.ico"> + x:Name="AboutWindow" Height="220" Width="280" MinHeight="220" MinWidth="280" Title="About" Icon="pack://application:,,,/Resources/icon.ico" ResizeMode="NoResize">