From: Greg Burri Date: Sat, 2 Jan 2016 16:20:22 +0000 (+0100) Subject: * Treat some special cases when ellipses intersecting. X-Git-Tag: 1.0.11~69 X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=commitdiff_plain;h=bef2e9f0bf1bba21d4c988fdf654c2dc303ec34a * Treat some special cases when ellipses intersecting. * Some adjustments. --- diff --git a/Parasitemia/Parasitemia/Classifier.fs b/Parasitemia/Parasitemia/Classifier.fs index 99b6edc..7f11e13 100644 --- a/Parasitemia/Parasitemia/Classifier.fs +++ b/Parasitemia/Parasitemia/Classifier.fs @@ -29,10 +29,10 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: let infection = parasites.infection.Copy() // To avoid to modify the parameter. // This is the minimum window size to check if other ellipses touch 'e'. - let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMaxRadius) * config.Parameters.scale - KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius) * config.Parameters.scale - KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius) * config.Parameters.scale - KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) * config.Parameters.scale } + let searchRegion (e: Ellipse) = { KdTree.minX = e.Cx - (e.A + config.RBCMaxRadius) + KdTree.maxX = e.Cx + (e.A + config.RBCMaxRadius) + KdTree.minY = e.Cy - (e.A + config.RBCMaxRadius) + KdTree.maxY = e.Cy + (e.A + config.RBCMaxRadius) } // The minimum window to contain a given ellipse. let ellipseWindow (e: Ellipse) = @@ -47,39 +47,57 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: // Return 'true' if the point 'p' is owned by e. // The lines represents all intersections with other ellipses. - let pixelOwnedByE (p: PointD) (e: Ellipse) (lines: Line list) = + let pixelOwnedByE (p: PointD) (e: Ellipse) (others: (Ellipse * Line) list) = e.Contains p.X p.Y && seq { let c = PointD(e.Cx, e.Cy) - for d1 in lines do + for e', d1 in others do let d2 = Utils.lineFromTwoPoints c p + let c' = PointD(e'.Cx, e'.Cy) + let v = pointFromTwoLines d1 (lineFromTwoPoints c c') + let case1 = sign (v.X - c.X) <> sign (v.X - c'.X) || Utils.squaredDistanceTwoPoints v c > Utils.squaredDistanceTwoPoints v c' if d2.Valid then let p' = Utils.pointFromTwoLines d1 d2 - yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p // 'false' -> the point is owned by another ellipse. + // Yield 'false' when the point is owned by another ellipse. + if case1 + then + yield sign (c.X - p.X) <> sign (c.X - p'.X) || Utils.squaredDistanceTwoPoints c p' > Utils.squaredDistanceTwoPoints c p + else + yield sign (c.X - p.X) = sign (c.X - p'.X) && Utils.squaredDistanceTwoPoints c p' < Utils.squaredDistanceTwoPoints c p else - yield true + yield case1 } |> Seq.forall id let ellipses = ellipses |> List.map EllipseFlaggedKd - // 1) Associate touching ellipses with each ellipses. + // 1) Associate touching ellipses with each ellipses and remove ellipse with more than two intersections. let tree = KdTree.Tree.BuildTree ellipses - let neighbors (e: Ellipse) : (EllipseFlaggedKd * PointD * PointD) list = - tree.Search (searchRegion e) - // We only keep the ellipses touching 'e'. - |> List.choose (fun otherE -> - match EEOver.EEOverlapArea e otherE with - | Some (area, px, py) when area > 0.0 && px.Length >= 2 && py.Length >= 2 -> - Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1])) - | _ -> - None ) + let neighbors (e: EllipseFlaggedKd) : (EllipseFlaggedKd * PointD * PointD) list = + if not e.Removed + then + tree.Search (searchRegion e) + // We only keep the ellipses touching 'e'. + |> List.choose (fun otherE -> + match EEOver.EEOverlapArea e otherE with + | Some (_, px, _) when px.Length > 2 -> + otherE.Removed <- true + None + | Some (area, px, py) when area > 0.0 && px.Length = 2 -> + Some (otherE, PointD(px.[0], py.[0]), PointD(px.[1], py.[1])) + | _ -> + None ) + else + [] // We reverse the list to get the lower score ellipses first. let ellipsesWithNeigbors = ellipses |> List.map (fun e -> e, neighbors e) |> List.rev // 2) Remove ellipses with a high standard deviation (high contrast). - let globalStdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq { + + // CvInvoke.Normalize(img, img, 0.0, 255.0, CvEnum.NormType.MinMax) // Not needed. + + let globalStdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation(seq { for y in 0 .. h - 1 do for x in 0 .. w - 1 do yield img.Data.[y, x, 0] |> float }) @@ -87,15 +105,14 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: for e in ellipses do let minX, minY, maxX, maxY = ellipseWindow e - let stdDiviation = MathNet.Numerics.Statistics.StreamingStatistics.StandardDeviation(seq { + let stdDeviation = MathNet.Numerics.Statistics.Statistics.StandardDeviation (seq { for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do if e.Contains (float x) (float y) then - yield img.Data.[y, x, 0] |> float }) + yield float img.Data.[y, x, 0] }) - if stdDiviation > globalStdDiviation * config.Parameters.standardDeviationMaxRatio - then + if stdDeviation > globalStdDeviation * config.Parameters.standardDeviationMaxRatio then e.Removed <- true // 3) Remove ellipses touching the edges. @@ -113,7 +130,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: for y in (if minY < 0 then 0 else minY) .. (if maxY >= h then h - 1 else maxY) do for x in (if minX < 0 then 0 else minX) .. (if maxX >= w then w - 1 else maxX) do let p = PointD(float x, float y) - if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) + if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2))) then area <- area + 1 @@ -139,7 +156,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img: for y in minY .. maxY do for x in minX .. maxX do let p = PointD(float x, float y) - if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (Utils.lineFromTwoPoints p1 p2))) + if pixelOwnedByE p e (neighbors |> List.choose (fun (otherE, p1, p2) -> if otherE.Removed then None else Some (otherE :> Ellipse, Utils.lineFromTwoPoints p1 p2))) then elements.[y-minY, x-minX] <- 1uy nbElement <- nbElement + 1 diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs index 45b1ca3..5740e8b 100644 --- a/Parasitemia/Parasitemia/Config.fs +++ b/Parasitemia/Parasitemia/Config.fs @@ -30,7 +30,7 @@ type Parameters = { infectionLevel: float // > 1 standardDeviationMaxRatio: float // The standard deviation of the pixel values of a cell can't be greater than standardDeviationMaxRatio * global standard deviation - minimumCellArea: float // Factor of RBC area. + minimumCellArea: float // Factor of the nominal RBC area. } type Config (param: Parameters) = diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs index 48c7f54..73771be 100644 --- a/Parasitemia/Parasitemia/Ellipse.fs +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -204,11 +204,11 @@ let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) let find (edges: Matrix) - (xDir: Image) - (yDir: Image) + (xGradient: Image) + (yGradient: Image) (config: Config) : MatchingEllipses = - let r1, r2 = config.Parameters.scale * config.RBCMinRadius, config.Parameters.scale * config.RBCMaxRadius // FIXME: scale factor should be applied in Config!? + let r1, r2 = config.RBCMinRadius, config.RBCMaxRadius let windowSize = roundInt (config.Parameters.factorWindowSize * r2) let factorNbPick = config.Parameters.factorNbPick @@ -227,8 +227,8 @@ let find (edges: Matrix) let currentElements = List<(int * int)>() let edgesData = edges.Data - let xDirData = xDir.Data - let yDirData = yDir.Data + let xDirData = xGradient.Data + let yDirData = yGradient.Data let rng = Random(42) diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs index 2f021c8..f567e68 100644 --- a/Parasitemia/Parasitemia/ImgTools.fs +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -21,6 +21,124 @@ let normalizeAndConvert (img: Image) : Image = ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).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 + + +let suppressMConnections (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 + +let findEdges (img: Image) : Matrix * Image * Image = + let w = img.Width + let h = img.Height + + use sobelKernel = + new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] + [ 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 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 + + 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 + + use magnitudes = new Matrix(xGradient.Size) + CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles). + + let thresholdHigh, thresholdLow = + let sensibility = 0.1 + 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) + + // Non-maximum suppression. + use nms = new Matrix(xGradient.Size) + nms.SetValue(1.0) + + for i in 0 .. h - 1 do + nms.Data.[i, 0] <- 0uy + nms.Data.[i, w - 1] <- 0uy + + for j in 0 .. w - 1 do + nms.Data.[0, j] <- 0uy + nms.Data.[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] + let angle = + let a = atan2 vy vx + if a < 0.0 then 2. * Math.PI + a else a + + let mNeigbors (sign: int) : float = + if angle < Math.PI / 8. || angle >= 15.0 * Math.PI / 8. then magnitudes.Data.[i, j + sign] + elif angle < 3.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j + sign] + elif angle < 5.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j] + elif angle < 7.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j - sign] + elif angle < 9.0 * Math.PI / 8. then magnitudes.Data.[i, j - sign] + elif angle < 11.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j - sign] + elif angle < 13.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j] + else magnitudes.Data.[i - sign, j + sign] + + let m = magnitudes.Data.[i, j] + if m < mNeigbors 1 || m < mNeigbors -1 || m < thresholdLow then + nms.Data.[i, j] <- 0uy + + // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection). + + let edges = new Matrix(xGradient.Size) + + // Histeresis 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 then + nms.Data.[i, j] <- 0uy + toVisit.Push(Point(j, i)) + while toVisit.Count > 0 do + let p = toVisit.Pop() + edges.Data.[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 then + nms.Data.[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) @@ -32,7 +150,6 @@ let drawPoints (img: Image) (points: Points) (intensity: byte) = for p in points do img.Data.[p.Y, p.X, 0] <- intensity - type ExtremumType = | Maxima = 1 | Minima = 2 @@ -356,6 +473,7 @@ let areaOpen (img: Image) (area: int) = let areaClose (img: Image) (area: int) = areaOperation img area AreaOperation.Closing +// A simpler algorithm than 'areaOpen' but slower. let areaOpen2 (img: Image) (area: int) = let w = img.Width let h = img.Height @@ -537,21 +655,14 @@ let connectedComponents (img: Image) (startPoints: List) : Li List(pointChecked) -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 - 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: Types.Ellipse) (color: 'TColor) (alpha: float) = if alpha >= 1.0 diff --git a/Parasitemia/Parasitemia/MainAnalysis.fs b/Parasitemia/Parasitemia/MainAnalysis.fs index c987e63..6111f06 100644 --- a/Parasitemia/Parasitemia/MainAnalysis.fs +++ b/Parasitemia/Parasitemia/MainAnalysis.fs @@ -47,32 +47,7 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li let filteredGreenWhitoutInfectionFloat = filteredGreenWhitoutInfection.Convert() let filteredGreenWhitoutStainFloat = filteredGreenWhitoutStain.Convert() - use sobelKernel = - new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] - [ 2.0f; 0.0f; -2.0f ] - [ 1.0f; 0.0f; -1.0f ]], Point(0, 0)) - - use xEdges = filteredGreenWhitoutStainFloat.Convolution(sobelKernel).Convert() - use yEdges = filteredGreenWhitoutStainFloat.Convolution(sobelKernel.Transpose()).Convert() - - let xEdgesData = xEdges.Data - let yEdgesData = yEdges.Data - for r in 0..xEdges.Rows-1 do - xEdgesData.[r, 0, 0] <- 0.0 - xEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 - yEdgesData.[r, 0, 0] <- 0.0 - yEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 - - for c in 0..xEdges.Cols-1 do - xEdgesData.[0, c, 0] <- 0.0 - xEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 - yEdgesData.[0, c, 0] <- 0.0 - yEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 - - use magnitudes = new Matrix(xEdges.Size) - CvInvoke.CartToPolar(xEdges, yEdges, magnitudes, new Mat()) // Compute the magnitudes (without angles). - - let min = ref 0.0 + (*let min = ref 0.0 let minLocation = ref <| Point() let max = ref 0.0 let maxLocation = ref <| Point() @@ -82,11 +57,14 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li use edges = new Matrix(xEdges.Size) CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) |> ignore - logTime "Finding edges" (fun() -> thin edges) + logTime "Finding edges" (fun() -> thin edges)*) + + let edges, xGradient, yGradient = ImgTools.findEdges filteredGreenWhitoutStainFloat + logTime "Removing small connected components from thinning" (fun () -> removeArea edges 12) let allEllipses, ellipses = logTime "Finding ellipses" (fun () -> - let matchingEllipses = Ellipse.find edges xEdges yEdges config + let matchingEllipses = Ellipse.find edges xGradient yGradient config matchingEllipses.Ellipses, matchingEllipses.PrunedEllipses ) let cells = logTime "Classifier" (fun () -> Classifier.findCells ellipses parasites filteredGreenWhitoutStain config) @@ -106,7 +84,7 @@ let doAnalysis (img: Image) (name: string) (config: Config) : Cell li saveImg parasites.infection (buildFileName " - parasites - infection.png") let imgAllEllipses = img.Copy() - drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.1 + drawEllipses imgAllEllipses allEllipses (Bgr(0.0, 240.0, 240.0)) 0.05 saveImg imgAllEllipses (buildFileName " - ellipses - all.png") let imgEllipses = img.Copy() diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs index 6e6218d..c39541e 100644 --- a/Parasitemia/Parasitemia/MatchingEllipses.fs +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -12,8 +12,7 @@ open Utils // Do not take in account matching score below this when two ellipses are matched. let matchingScoreThreshold1 = 0.6 -// All ellipses with a score below this are removed. -let matchingScoreThreshold2 = 600. +let scaleOverlapTest = 0.8 type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) = let mutable matchingScore = matchingScore @@ -36,6 +35,9 @@ type private EllipseScoreFlaggedKd (matchingScore: float, e: Ellipse) = type MatchingEllipses (radiusMin: float) = let ellipses = List() + // All ellipses with a score below this are removed. + let matchingScoreThreshold2 = 25. * radiusMin // 600. + member this.Add (e: Ellipse) = ellipses.Add(EllipseScoreFlaggedKd(0.0, e)) @@ -76,7 +78,7 @@ type MatchingEllipses (radiusMin: float) = // 3) Sort ellipses by their score. ellipses.Sort(fun e1 e2 -> e2.MatchingScore.CompareTo(e1.MatchingScore)) - // 4) Remove ellipses wich have a low score. + // 4) Remove ellipses with a low score. let i = ellipses.BinarySearch(EllipseScoreFlaggedKd(matchingScoreThreshold2, Ellipse(0.0, 0.0, 0.0, 0.0, 0.0)), { new IComparer with member this.Compare(e1, e2) = e2.MatchingScore.CompareTo(e1.MatchingScore) }) |> abs @@ -98,7 +100,7 @@ type MatchingEllipses (radiusMin: float) = for other in tree.Search window do if not other.Removed && other.MatchingScore < e.MatchingScore then - if e.Ellipse.Scale(0.8).Contains other.Ellipse.Cx other.Ellipse.Cy + if e.Ellipse.Scale(scaleOverlapTest).Contains other.Ellipse.Cx other.Ellipse.Cy then other.Removed <- true ellipses.RemoveAll(fun e -> e.Removed) |> ignore diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs index 528f2de..dafff56 100644 --- a/Parasitemia/Parasitemia/Program.fs +++ b/Parasitemia/Parasitemia/Program.fs @@ -69,7 +69,7 @@ let main args = factorNbPick = 1.0 factorWindowSize = 2.0 - darkStainLevel = 0.3 // Lower -> more sensitive. + darkStainLevel = 0.25 // Lower -> more sensitive. 0.3 maxDarkStainRatio = 0.1 infectionArea = 0.012 // 1.2 % @@ -79,7 +79,7 @@ let main args = stainLevel = 1.1 // Lower -> more sensitive. maxStainRatio = 0.12 // 12 % - standardDeviationMaxRatio = 0.55 + standardDeviationMaxRatio = 0.58 // 0.55 minimumCellArea = 0.5 }) match mode with diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs index d6fa2f0..58e914f 100644 --- a/Parasitemia/Parasitemia/Utils.fs +++ b/Parasitemia/Parasitemia/Utils.fs @@ -30,6 +30,10 @@ let inline pointFromTwoLines (l1: Line) (l2: Line) : PointD = let y = -(l2.A * l1.B - l1.A * l2.B) / (l1.A - l2.A) PointD(x, y) +let inline linePassThroughSegment (l: Line) (p1: PointD) (p2: PointD) : bool = + let p = pointFromTwoLines l (lineFromTwoPoints p1 p2) + sign (p.X - p1.X) <> sign (p.X - p2.X) + let inline squaredDistanceTwoPoints (p1: PointD) (p2: PointD) = (p1.X - p2.X) ** 2.0 + (p1.Y - p2.Y) ** 2.0