diff options
Diffstat (limited to 'src/Geometry.hs')
| -rw-r--r-- | src/Geometry.hs | 230 |
1 files changed, 230 insertions, 0 deletions
diff --git a/src/Geometry.hs b/src/Geometry.hs new file mode 100644 index 0000000..dd8d8e1 --- /dev/null +++ b/src/Geometry.hs @@ -0,0 +1,230 @@ +module Geometry + ( Line (..), + lineThrough2Points, + lineThroughPointWithAngle, + LineSide (..), + sideOfLine, + proj, + perp', + orthoproj, + distPointLine, + distPointSegment, + rotate, + -- Rectangle types + GeometryRectangle (..), + rectangleToGeometry, + geometryToRectangle, + -- Intersection + IntersectionResult (..), + getIntersectingRectangles, + ) +where + +import Data.Int (Int32) +import Data.List (sortOn) +import Data.Maybe (catMaybes, listToMaybe) +import Graphics.X11 (Rectangle (..)) +import Prelude hiding (subtract) + +-- | A line in parametric form: (x + t*dx, y + t*dy) for t ∈ ℝ +data Line a = Line + { x :: !a, + y :: !a, + dx :: !a, + dy :: !a + } + deriving (Eq, Show, Read) + +-- | Construct a line passing through two points +lineThrough2Points :: (Num a, Ord a) => (a, a) -> (a, a) -> Line a +lineThrough2Points (x1, y1) (x2, y2) = + Line + { x = x1, + y = y1, + dx = x2 - x1, + dy = y2 - y1 + } + +-- | Construct a line through a point with a given angle (in radians) +-- Angle is measured counterclockwise from the positive x-axis +lineThroughPointWithAngle :: (Floating a, Ord a) => (a, a) -> a -> Line a +lineThroughPointWithAngle (x, y) angle = + Line + { x = x, + y = y, + dx = cos angle, + dy = sin angle + } + +-- | Which side of the line is a point on? +data LineSide + = LeftSide + | RightSide + | OnLine + deriving (Eq, Show) + +-- | Determine which side of the line a point lies on +-- Uses the cross product to determine orientation +sideOfLine :: (Fractional a, Ord a) => Line a -> (a, a) -> LineSide +sideOfLine line (px, py) = + case compare cross 0 of + GT -> LeftSide + LT -> RightSide + EQ -> OnLine + where + -- Vector from line origin to point + relX = px - x line + relY = py - y line + -- Cross product in 2D: (dx, dy) × (px, py) = dx*py - dy*px + cross = dx line * relY - dy line * relX + +-- | Project point onto line +proj :: (Fractional a, Ord a) => Line a -> (a, a) -> (a, a) +proj line (px, py) = + (x line + t * dx line, y line + t * dy line) + where + -- Vector from line origin to point + relX = px - x line + relY = py - y line + -- t = (rel · direction) / (direction · direction) + dotProd = relX * dx line + relY * dy line + dirSq = dx line * dx line + dy line * dy line + t = dotProd / dirSq + +-- | Perpendicular vector from point to line +perp' :: (Fractional a, Ord a) => Line a -> (a, a) -> (a, a) +perp' line point = + (px - projX, py - projY) + where + (projX, projY) = proj line point + (px, py) = point + +-- | Orthogonal projection (same as proj, but returns the projected point) +orthoproj :: (Fractional a, Ord a) => Line a -> (a, a) -> (a, a) +orthoproj = proj + +-- | Distance from point to line +distPointLine :: (Floating a, Ord a) => Line a -> (a, a) -> a +distPointLine line point = + sqrt (dx * dx + dy * dy) + where + (dx, dy) = perp' line point + +-- | Distance from point to line segment (not infinite line) +distPointSegment :: (Floating a, Ord a) => (a, a) -> (a, a) -> (a, a) -> a +distPointSegment (ax, ay) (bx, by) (px, py) + | t <= 0 = sqrt (dx1 * dx1 + dy1 * dy1) + | t >= 1 = sqrt (dx2 * dx2 + dy2 * dy2) + | otherwise = sqrt (dx * dx + dy * dy) + where + abX = bx - ax + abY = by - ay + apX = px - ax + apY = py - ay + t = (apX * abX + apY * abY) / (abX * abX + abY * abY) + projX = ax + t * abX + projY = ay + t * abY + dx = px - projX + dy = py - projY + dx1 = px - ax + dy1 = py - ay + dx2 = px - bx + dy2 = py - by + +-- | Rotate a point around the origin by an angle (in radians) +rotate :: (Floating a, Ord a) => a -> (a, a) -> (a, a) +rotate angle (x, y) = + (x * cos angle - y * sin angle, x * sin angle + y * cos angle) + +-- | A simple Rectangle type for geometry operations +-- Uses (x, y) as the top-left corner, with width and height +data GeometryRectangle a = GeometryRectangle + { geoX :: !a, + geoY :: !a, + geoWidth :: !a, + geoHeight :: !a + } + deriving (Eq, Show, Read) + +-- | Convert X11's Rectangle to Geometry's Rectangle +rectangleToGeometry :: Num a => Graphics.X11.Rectangle -> GeometryRectangle a +rectangleToGeometry (Graphics.X11.Rectangle x y w h) = + GeometryRectangle + { geoX = fromIntegral x, + geoY = fromIntegral y, + geoWidth = fromIntegral w, + geoHeight = fromIntegral h + } + +-- | Convert Geometry's Rectangle to X11's Rectangle +geometryToRectangle :: Integral a => GeometryRectangle a -> Graphics.X11.Rectangle +geometryToRectangle (GeometryRectangle x y w h) = + Graphics.X11.Rectangle (fromIntegral x) (fromIntegral y) (fromIntegral w) (fromIntegral h) + +-- | Alias for GeometryRectangle for backwards compatibility +type Rectangle a = GeometryRectangle a + +-- | Result of intersecting a line with a rectangle +data IntersectionResult a b = IntersectionResult + { tag :: a, + rectangle :: GeometryRectangle b, + intersectionPoint :: (b, b) + } + deriving (Eq, Show) + +-- | Find all rectangles that intersect the given line, sorted by distance from origin +-- Returns tags of intersected rectangles in order of their intersection point +-- closest to (0,0) +getIntersectingRectangles :: (Floating a, Ord a) => Line a -> [(GeometryRectangle a, b)] -> [IntersectionResult b a] +getIntersectingRectangles line rects = catMaybes $ map (checkIntersection line) rects + +-- | Check if a line intersects a rectangle, and if so, return the intersection point +-- closest to the origin +checkIntersection :: (Floating a, Ord a) => Line a -> (GeometryRectangle a, b) -> Maybe (IntersectionResult b a) +checkIntersection line (rect, tag) = IntersectionResult tag rect <$> findClosestIntersection line rect + +-- | Find the intersection point closest to the origin +findClosestIntersection :: (Floating a, Ord a) => Line a -> GeometryRectangle a -> Maybe (a, a) +findClosestIntersection line rect = + listToMaybe $ sortOn distanceFromOrigin $ filter (pointInRectangle rect) candidates + where + -- Find intersections with each edge of the rectangle + topEdge = lineThrough2Points (geoX rect, geoY rect) (geoX rect + geoWidth rect, geoY rect) + bottomEdge = lineThrough2Points (geoX rect, geoY rect + geoHeight rect) (geoX rect + geoWidth rect, geoY rect + geoHeight rect) + leftEdge = lineThrough2Points (geoX rect, geoY rect) (geoX rect, geoY rect + geoHeight rect) + rightEdge = lineThrough2Points (geoX rect + geoWidth rect, geoY rect) (geoX rect + geoWidth rect, geoY rect + geoHeight rect) + + candidates = + catMaybes + [ intersectLines line topEdge, + intersectLines line bottomEdge, + intersectLines line leftEdge, + intersectLines line rightEdge + ] + +-- | Check if a point is inside a rectangle (including edges) +pointInRectangle :: (Num a, Ord a) => GeometryRectangle a -> (a, a) -> Bool +pointInRectangle rect (px, py) = + px >= geoX rect + && px <= geoX rect + geoWidth rect + && py >= geoY rect + && py <= geoY rect + geoHeight rect + +-- | Find the intersection point of two lines, if it exists +intersectLines :: (Fractional a, Eq a) => Line a -> Line a -> Maybe (a, a) +intersectLines (Line x1 y1 dx1 dy1) (Line x2 y2 dx2 dy2) = + -- Solve: x1 + t*dx1 = x2 + s*dx2 and y1 + t*dy1 = y2 + s*dy2 + -- Using Cramer's rule + let det = dx1 * (- dy2) - dy1 * (- dx2) + detT = (x2 - x1) * (- dy2) - (y2 - y1) * (- dx2) + in if det == 0 + then Nothing -- Parallel lines + else + let t = detT / det + intersectionX = x1 + t * dx1 + intersectionY = y1 + t * dy1 + in Just (intersectionX, intersectionY) + +-- | Euclidean distance from origin +distanceFromOrigin :: (Num a, Floating a) => (a, a) -> a +distanceFromOrigin (x, y) = sqrt (x * x + y * y) |