X-Git-Url: http://git.euphorik.ch/?a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=3cfdc89fe9f3264168c97c1ec0e72969c9d2b7da;hb=6b550c3faf4dea77738fa5c27cd9af277f45549c;hp=cee21c77d5c60fa26c8e25f4c2afa34c15056543;hpb=d0c85068bb98a7999ed994f02669befa70edd5f9;p=master-thesis.git diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index cee21c7..3cfdc89 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -8,17 +8,18 @@ open System.Linq open Emgu.CV open Emgu.CV.Structure -open Utils open Heap +open Const +open Utils // Normalize image values between 0uy and 255uy. -let normalizeAndConvert (img: Image) : Image = +let normalizeAndConvert (img: Image) : 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) - ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() + ((img.Convert() - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).Convert() let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) = @@ -31,6 +32,123 @@ let saveMat (mat: Matrix<'TDepth>) (filepath: string) = 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 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 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 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) |> Array.sum |> float + + 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 + + let suppressMConnections (img: Matrix) = let w = img.Width let h = img.Height @@ -45,7 +163,8 @@ let suppressMConnections (img: Matrix) = then img.[i, j] <- 0uy -let findEdges (img: Image) : Matrix * Image * Image = + +let findEdges (img: Image) : Matrix * Image * Image = let w = img.Width let h = img.Height @@ -54,110 +173,114 @@ let findEdges (img: Image) : Matrix * Image * [ 2.0f; 0.0f; -2.0f ] [ 1.0f; 0.0f; -1.0f ]], Point(1, 1)) - let xGradient = img.Convolution(sobelKernel).Convert() - let yGradient = img.Convolution(sobelKernel.Transpose()).Convert() + let xGradient = img.Convolution(sobelKernel) + let yGradient = img.Convolution(sobelKernel.Transpose()) let xGradientData = xGradient.Data let yGradientData = yGradient.Data for r in 0 .. h - 1 do - xGradientData.[r, 0, 0] <- 0.0 - xGradientData.[r, w - 1, 0] <- 0.0 - yGradientData.[r, 0, 0] <- 0.0 - yGradientData.[r, w - 1, 0] <- 0.0 + xGradientData.[r, 0, 0] <- 0.f + xGradientData.[r, w - 1, 0] <- 0.f + yGradientData.[r, 0, 0] <- 0.f + yGradientData.[r, w - 1, 0] <- 0.f for c in 0 .. w - 1 do - xGradientData.[0, c, 0] <- 0.0 - xGradientData.[h - 1, c, 0] <- 0.0 - yGradientData.[0, c, 0] <- 0.0 - yGradientData.[h - 1, c, 0] <- 0.0 + xGradientData.[0, c, 0] <- 0.f + xGradientData.[h - 1, c, 0] <- 0.f + yGradientData.[0, c, 0] <- 0.f + yGradientData.[h - 1, c, 0] <- 0.f - use magnitudes = new Matrix(xGradient.Size) - use angles = new Matrix(xGradient.Size) + use magnitudes = new Matrix(xGradient.Size) + use angles = new Matrix(xGradient.Size) CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes (without angles). let thresholdHigh, thresholdLow = - let sensibility = 0.1 + let sensibilityHigh = 0.1f + let sensibilityLow = 0.0f use magnitudesByte = magnitudes.Convert() - let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) - threshold + (sensibility * threshold), threshold - (sensibility * threshold) + let threshold = float32 <| CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) + 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 + + let PI = float32 Math.PI + for i in 0 .. h - 1 do - nms.Data.[i, 0] <- 0uy - nms.Data.[i, w - 1] <- 0uy + nmsData.[i, 0] <- 0uy + nmsData.[i, w - 1] <- 0uy for j in 0 .. w - 1 do - nms.Data.[0, j] <- 0uy - nms.Data.[h - 1, j] <- 0uy + 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 = xGradient.Data.[i, j, 0] - let vy = yGradient.Data.[i, j, 0] - if vx <> 0. || vy <> 0. + let vx = xGradientData.[i, j, 0] + let vy = yGradientData.[i, j, 0] + if vx <> 0.f || vy <> 0.f then - let angle = angles.[i, j] + 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. - ratio2 - - let mNeigbors (sign: int) : float = - if angle < Math.PI / 4. - then - ratio1 * magnitudes.Data.[i, j + sign] + ratio2 * magnitudes.Data.[i + sign, j + sign] - elif angle < Math.PI / 2. - then - ratio2 * magnitudes.Data.[i + sign, j + sign] + ratio1 * magnitudes.Data.[i + sign, j] - elif angle < 3.0 * Math.PI / 4. - then - ratio1 * magnitudes.Data.[i + sign, j] + ratio2 * magnitudes.Data.[i + sign, j - sign] - elif angle < Math.PI - then - ratio2 * magnitudes.Data.[i + sign, j - sign] + ratio1 * magnitudes.Data.[i, j - sign] - elif angle < 5. * Math.PI / 4. - then - ratio1 * magnitudes.Data.[i, j - sign] + ratio2 * magnitudes.Data.[i - sign, j - sign] - elif angle < 3. * Math.PI / 2. - then - ratio2 * magnitudes.Data.[i - sign, j - sign] + ratio1 * magnitudes.Data.[i - sign, j] - elif angle < 7. * Math.PI / 4. - then - ratio1 * magnitudes.Data.[i - sign, j] + ratio2 * magnitudes.Data.[i - sign, j + sign] - else - ratio2 * magnitudes.Data.[i - sign, j + sign] + ratio1 * magnitudes.Data.[i, j + sign] - - let m = magnitudes.Data.[i, j] + 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 - nms.Data.[i, j] <- 1uy + 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 - // Histeresis thresholding. + // Hysteresis thresholding. let toVisit = Stack() for i in 0 .. h - 1 do for j in 0 .. w - 1 do - if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh + if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh then - nms.Data.[i, j] <- 0uy + nmsData.[i, j] <- 0uy toVisit.Push(Point(j, i)) while toVisit.Count > 0 do let p = toVisit.Pop() - edges.Data.[p.Y, p.X] <- 1uy + 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 && nms.Data.[ni, nj] = 1uy + if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy then - nms.Data.[ni, nj] <- 0uy + nmsData.[ni, nj] <- 0uy toVisit.Push(Point(nj, ni)) edges, xGradient, yGradient @@ -729,7 +852,7 @@ let removeArea (mat: Matrix) (areaSize: int) = ( 0, -1) // p8 (-1, -1) |] // p9 - let mat' = new Matrix(mat.Size) + use mat' = new Matrix(mat.Size) let w = mat'.Width let h = mat'.Height mat.CopyTo(mat') @@ -797,21 +920,21 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColo if alpha >= 1.0 then - img.Draw(Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias) + img.Draw(Ellipse(PointF(float32 e.Cx, float32 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.0 - let gapX = windowPosX - (float (int windowPosX)) + let windowPosX = e.Cx - e.A - 5.f + let gapX = windowPosX - (float32 (int windowPosX)) - let windowPosY = e.Cy - e.A - 5.0 - let gapY = windowPosY - (float (int windowPosY)) + let windowPosY = e.Cy - e.A - 5.f + let gapY = windowPosY - (float32 (int windowPosY)) - let roi = Rectangle(int windowPosX, int windowPosY, 2. * (e.A + 5.0) |> int, 2.* (e.A + 5.0) |> int) + 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(Ellipse(PointF(float32 <| (e.A + 5. + gapX) , float32 <| (e.A + 5. + gapY)), SizeF(2. * e.B |> float32, 2. * e.A |> float32), float32 <| e.Alpha / Math.PI * 180.), color, 1, CvEnum.LineType.AntiAlias) + i.Draw(Ellipse(PointF(float32 <| (e.A + 5.f + gapX) , float32 <| (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