X-Git-Url: http://git.euphorik.ch/?a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FImgTools.fs;h=23c9aee38e18c6f41bd684ecaeb1ea44a8dabd70;hb=3ddaf64dc5ba6a7066a279ad75b9a1ee72194639;hp=f64b1c3318d687d6dc314fe98f4dbdd125446d6f;hpb=044b0ae69df3ac565432545b2fa934589016f9bd;p=master-thesis.git diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index f64b1c3..23c9aee 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -32,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 @@ -79,9 +196,9 @@ let findEdges (img: Image) : Matrix * Image let thresholdHigh, thresholdLow = let sensibilityHigh = 0.1f - let sensibilityLow = 0.1f + let sensibilityLow = 0.0f use magnitudesByte = magnitudes.Convert() - 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. @@ -437,9 +554,9 @@ let private areaOperation (img: Image) (area: int) (op: AreaOperatio nextElements.Add(p) |> ignore else - let m' = pixels.[p.Y, p.X] - if m' <> null - then + match pixels.[p.Y, p.X] with + | null -> () + | m' -> if m'.Elements.Count + m.Elements.Count <= area then m'.State <- AreaState.Removed @@ -509,7 +626,7 @@ type Island (cmp: IComparer) = member val Surface = 0 with get, set -let private areaOperationF (img: Image) (area: int) (op: AreaOperation) = +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 @@ -536,82 +653,97 @@ let private areaOperationF (img: Image) (area: int) (op: AreaOper let ni = i + p.Y let nj = j + p.X let neighbor = Point(nj, ni) - if ni >= 0 && ni < h && nj >= 0 && nj < w && ownership.[ni, nj] = null && not (shorePoints.Contains(neighbor)) + 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 island in islands do - let mutable stop = island.Shore.IsEmpty - - // 'true' if 'p' is owned or adjacent to 'island'. - let 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 other <> null - then // We touching another island. - if island.Surface + other.Surface >= area - then - stop <- true - else // We can merge 'other' into 'surface'. - island.Surface <- island.Surface + other.Surface - island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.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 + 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 - 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 + else + if not <| Object.ReferenceEquals(other, null) + then // We touching another island. + if island.Surface + other.Surface >= area 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 + else // We can merge 'other' into 'surface'. + island.Surface <- island.Surface + other.Surface + island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.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 - ownership.[next.Y, next.X] <- island - island.Level <- level - island.Surface <- island.Surface + 1 - - for i in 0 .. h - 1 do - for j in 0 .. w - 1 do - let island = ownership.[i, j] - if island <> null - then - earth.[i, j, 0] <- island.Level + 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 + | _ -> () () let areaOpenF (img: Image) (area: int) = - areaOperationF img area AreaOperation.Opening + areaOperationF img [ area, () ] None AreaOperation.Opening let areaCloseF (img: Image) (area: int) = - areaOperationF img area AreaOperation.Closing + areaOperationF img [ area, () ] None AreaOperation.Closing + +let areaOpenFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Opening + +let areaCloseFWithFun (img: Image) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) = + areaOperationF img areas (Some f) AreaOperation.Closing // A simpler algorithm than 'areaOpen' but slower. let areaOpen2 (img: Image) (area: int) =