First, the new definitions:

data Region = UnitCircle
| Polygon [Coordinate]
| AffineTransform Matrix3x3 Region
| Empty
deriving Show
type Vector3 = (Float, Float, Float)
type Matrix3x3 = (Vector3, Vector3, Vector3)
type Matrix2x2 = ((Float, Float), (Float, Float)) |

data Region = UnitCircle
| Polygon [Coordinate]
| AffineTransform Matrix3x3 Region
| Empty
deriving Show
type Vector3 = (Float, Float, Float)
type Matrix3x3 = (Vector3, Vector3, Vector3)
type Matrix2x2 = ((Float, Float), (Float, Float))

And some old ones that will suffice unchanged:

type Coordinate = (Float, Float)
type Ray = (Coordinate, Coordinate)
isLeftOf :: Coordinate -> Ray -> Bool
(px,py) `isLeftOf` ((ax,ay), (bx,by))
= let (s,t) = (px-ax, py-ay)
(u,v) = (px-bx, py-by)
in s * v >= t * u
isRightOf :: Coordinate -> Ray -> Bool
(px,py) `isRightOf` ((ax,ay), (bx,by))
= let (s,t) = (px-ax, py-ay)
(u,v) = (px-bx, py-by)
in s * v <= t * u
containsR :: Region -> Coordinate -> Bool
(Polygon pts) `containsR` p
= let leftOfList = map isLeftOfp (zip pts (tail pts ++ [head pts]))
isLeftOfp p' = isLeftOf p p'
rightOfList = map isRightOfp (zip pts (tail pts ++ [head pts]))
isRightOfp p' = isRightOf p p'
in and leftOfList || and rightOfList
Empty `containsR` p = False |

type Coordinate = (Float, Float)
type Ray = (Coordinate, Coordinate)
isLeftOf :: Coordinate -> Ray -> Bool
(px,py) `isLeftOf` ((ax,ay), (bx,by))
= let (s,t) = (px-ax, py-ay)
(u,v) = (px-bx, py-by)
in s * v >= t * u
isRightOf :: Coordinate -> Ray -> Bool
(px,py) `isRightOf` ((ax,ay), (bx,by))
= let (s,t) = (px-ax, py-ay)
(u,v) = (px-bx, py-by)
in s * v <= t * u
containsR :: Region -> Coordinate -> Bool
(Polygon pts) `containsR` p
= let leftOfList = map isLeftOfp (zip pts (tail pts ++ [head pts]))
isLeftOfp p' = isLeftOf p p'
rightOfList = map isRightOfp (zip pts (tail pts ++ [head pts]))
isRightOfp p' = isRightOf p p'
in and leftOfList || and rightOfList
Empty `containsR` p = False

Now we have the new forms of `containsR` to write. The `UnitCircle` is easy:

UnitCircle `containsR` (x,y) = x^2 + y^2 <= 1 |

UnitCircle `containsR` (x,y) = x^2 + y^2 <= 1

And in general, the transform of a region contains a point if the region contains the inverse transform of that point.

(AffineTransform m r) `containsR` (x,y)
= if determinant3 m == 0
then singularContains m (x,y)
else let m' = inverse m
(x', y', _) = matrixMul m' (x,y,1)
in r `containsR` (x', y') |

(AffineTransform m r) `containsR` (x,y)
= if determinant3 m == 0
then singularContains m (x,y)
else let m' = inverse m
(x', y', _) = matrixMul m' (x,y,1)
in r `containsR` (x', y')

Now some standard code for multiplying and inverting matrices:

matrixMul :: Matrix3x3 -> Vector3 -> Vector3
matrixMul (r1, r2, r3) v
= (dotProduct r1 v,
dotProduct r2 v,
dotProduct r3 v)
dotProduct :: Vector3 -> Vector3 -> Float
dotProduct (a,b,c) (x,y,z) = a*x + b*y + c*z
inverse :: Matrix3x3 -> Matrix3x3
inverse ((a,b,c), (d,e,f), (g,h,i))
= let det = determinant3 ((a,b,c), (d,e,f), (g,h,i))
a' = determinant2 ((e,f), (h,i)) / det
b' = determinant2 ((c,b), (i,h)) / det
c' = determinant2 ((b,c), (e,f)) / det
d' = determinant2 ((f,d), (i,g)) / det
e' = determinant2 ((a,c), (g,i)) / det
f' = determinant2 ((c,a), (f,d)) / det
g' = determinant2 ((d,e), (g,h)) / det
h' = determinant2 ((b,a), (h,g)) / det
i' = determinant2 ((a,b), (d,e)) / det
in ((a',b',c'), (d',e',f'), (g',h',i'))
determinant3 :: Matrix3x3 -> Float
determinant3 ((a,b,c), (d,e,f), (g,h,i))
= let aei = a * e * i
afh = a * f * h
bdi = b * d * i
cdh = c * d * h
bfg = b * f * g
ceg = c * e * g
in aei - afh - bdi + cdh + bfg - ceg
determinant2 :: Matrix2x2 -> Float
determinant2 ((a,b), (c,d))
= a * d - b * c |

matrixMul :: Matrix3x3 -> Vector3 -> Vector3
matrixMul (r1, r2, r3) v
= (dotProduct r1 v,
dotProduct r2 v,
dotProduct r3 v)
dotProduct :: Vector3 -> Vector3 -> Float
dotProduct (a,b,c) (x,y,z) = a*x + b*y + c*z
inverse :: Matrix3x3 -> Matrix3x3
inverse ((a,b,c), (d,e,f), (g,h,i))
= let det = determinant3 ((a,b,c), (d,e,f), (g,h,i))
a' = determinant2 ((e,f), (h,i)) / det
b' = determinant2 ((c,b), (i,h)) / det
c' = determinant2 ((b,c), (e,f)) / det
d' = determinant2 ((f,d), (i,g)) / det
e' = determinant2 ((a,c), (g,i)) / det
f' = determinant2 ((c,a), (f,d)) / det
g' = determinant2 ((d,e), (g,h)) / det
h' = determinant2 ((b,a), (h,g)) / det
i' = determinant2 ((a,b), (d,e)) / det
in ((a',b',c'), (d',e',f'), (g',h',i'))
determinant3 :: Matrix3x3 -> Float
determinant3 ((a,b,c), (d,e,f), (g,h,i))
= let aei = a * e * i
afh = a * f * h
bdi = b * d * i
cdh = c * d * h
bfg = b * f * g
ceg = c * e * g
in aei - afh - bdi + cdh + bfg - ceg
determinant2 :: Matrix2x2 -> Float
determinant2 ((a,b), (c,d))
= a * d - b * c

The only thing left is the question of how to deal with a singular (non-invertible) matrix. If an affine matrix is non-invertible, that means that AE – BD = 0. Either all four coefficients are 0, or AE = BD. If all 4 are 0, then every point in the region will be collapsed into the single point (C,F), and we need to check that this is the given point. If AE = BD otherwise, we have for any point `(x,y)` in the region:

x' = Ax + By + C
y' = Dx + Ey + F
Dx' = ADx + BDy + CD
Ay' = ADx + AEy + AF
= ADx + BDy + AF [AE = BD]
∴ Dx' - Ay' = CD - AF
Dx' - Ay' + AF - CD = 0

so any point in the region will be collapsed into the line Dx’ – Ay’ + AF – CD = 0, and we need to check that the given point satisfies this equation.

singularContains :: Matrix3x3 -> Coordinate -> Bool
singularContains ((a,b,c), (d,e,f), _) (x,y)
= if (a == 0 && b == 0 && d == 0 && e == 0)
then (x == c && y == f)
else (d*x - a*y + a*f - c*d == 0) |

singularContains :: Matrix3x3 -> Coordinate -> Bool
singularContains ((a,b,c), (d,e,f), _) (x,y)
= if (a == 0 && b == 0 && d == 0 && e == 0)
then (x == c && y == f)
else (d*x - a*y + a*f - c*d == 0)