Fix homogeneous coordinates (and the affine transformations on them)

I'm a little bit confused about how Affine is supposed to work in the linear package.

In particular, affine matrix transformations. The api encourages you to create a 4x4 matrix for your affine transformation:

mkTransformation :: Num a => Quaternion a -> V3 a -> M44 a
mkTransformationMat :: Num a => M33 a -> V3 a -> M44 a

The M44 result presumably needs to be multiplied with a V4 representation in homogeneous coordinates. This is pretty standard fare in computer graphics libraries:

let affineMat = mkTransformationMat rotationScale translation :: M44 Double
     homogeneousVec = V4 x y z 1 :: V4 Double
in affineMat !* homogeneousVec

Another convention that I believe is reasonably standard in computer graphics packages is for Point to represent homogeneous coordinates (since points can be translated and projected while vectors typically represent a scale and direction only). Linear.Affine exposes a Point type which presumably serves this function, but...

1. `Point` versus `point`

In linear Point is just a newtype over the underlying vector type. Confusingly neither point nor normalizePoint has anything to do with Point:

point :: Num a => V3 a -> V4 a
normalizePoint :: Fractional a => V4 a -> V3 a

2. Affine transformations on points

Furthermore, I can't quite make out how one would go about multiplying the before mentioned M44 transformation matrix with my Point:

let affineMat = mkTransformationMat rotationScale translation :: M44 Double
     homogeneousVec = P (V3 x y z) :: Point V3 Double
in affineMat !* homogeneousVec  --    Couldn't match type ‘Point f’ with ‘V3’
                                --      Expected type: V3 Double
                                --      Actual type: Point f Double
                                --      In the second argument of ‘(!*)’, namely ‘P v’

In addition, it is also sometimes useful to use a non-square matrix when there is no perspective projection/shear present in the transform. This is especially true with 2D graphics where most transformations are rotation/scale/translation, conveniently representable as M23:

let scaleAndTranslate = V2
                                         (V3 sx 0 tx)
                                         (V3 0 sy ty)
     homogeneousVec = P (V2 x y) -- Alternatively: V3 x y 1
in affineMat !* homogeneousVec -- Couldn't match type ‘Point f’ with ‘V2’

How should I deal with my shear / translate / perspective transformations?

There's the temptation to add conversion functions:

homogeneous2 :: Point V2 Double -> V3 Double
homogeneous2 (P (V2 x y)) = V3 x y 1

homogeneous3 :: Point V3 Double -> V4 Double
homogeneous3 (P (V3 x y z)) = V4 x y z 1

But this is far from ideal: if you started out with Point then you'd expect to have a Point result after applying a transformation to it.

What to do?

In conclusion, looking at this from the application programmer side, the Point type looks useless since you inevitably need to drop it off of your type to do anything useful with them.

I'm not sure what the best approach would be to fixing this would be or whether I've misunderstood the intended interface for Linear.Affine: E.g. Affine is nice, but I'm used to using matrix multiplication in order to compose long chains of transformations. Furthermore, shear and perspective projection is not available, which is weird considering wikipedia's description:

The above-mentioned augmented matrix is called affine transformation matrix, or projective transformation matrix (as it can also be used to perform projective transformations).

Thanks for any help!

I don't think I've used Point much if at all, so I can't speak from experience. ~~Perhaps the matrix-vector multiplication functions should be phrased in terms of R4 to let the user benefit from Point's intended type safety.~~ (EDIT: I should look at the code before guessing what the problem is.)

The objection then is that it's a bit heavy to write, P (mat !* pt^._xyzw). If the isomorphism between Point f a and f a were automatically plumbed by a type class, we could write a (!*) in terms of that. Or the Affine module could export matrix-vector multiplication adapted to Point. Is either of those appealing?

(For clarity, this was the expression I typed into GHCi to see how things hang together: P (mkTransformation (axisAngle (V3 0 0 1) (pi / 2)) (V3 1 2 0) !* (P (point (V3 1 1 0)))^._xyzw))

The objection then is that it's a bit heavy to write, P (mat !* pt^._xyzw)

Unfortunately I think the problem may run a little a bit deeper than that. For example, assuming that we have an w explicit element for our homogeneous coordinates, then how should .+^ look?

let p = _ :: Point V4 Double
     v = _ :: V3 Double
in p .+^ v :: Point V4 Double   -- This does not conform to the current interface...

alternatively

let p = _ :: Point V4 Double
     v = _ :: V4 Double  -- Having w ≠ 0 doesn't make sense when adding a vector to a point...
in p .+^ v :: Point V4 Double

I think the proper conversion between Point and vector may be point and normalizePoint - and notice the implementation for normalizePoint is this:

normalizePoint (V4 a b c w) = (1/w) *^ V3 a b c

At least, where perspective projection is concerned, this is how it works:

Or, in matrix form using homogeneous coordinates, the system

${\begin{bmatrix}\mathbf {f} _{x}\\\mathbf {f} _{y}\\\mathbf {f} _{z}\\\mathbf {f} _{w}\\\end{bmatrix}}={\begin{bmatrix}1&0&-{\frac {\mathbf {e} _{x}}{\mathbf {e} _{z}}}&0\\0&1&-{\frac {\mathbf {e} _{y}}{\mathbf {e} _{z}}}&0\\0&0&1&0\\0&0&1/\mathbf {e} _{z}&0\\\end{bmatrix}}{\begin{bmatrix}\mathbf {d} _{x}\\\mathbf {d} _{y}\\\mathbf {d} _{z}\\1\\\end{bmatrix}}} {\begin{bmatrix}{\mathbf {f}}_{x}\\{\mathbf {f}}_{y}\\{\mathbf {f}}_{z}\\{\mathbf {f}}_{w}\\\end{bmatrix}}={\begin{bmatrix}1&0&-{\frac {{\mathbf {e}}_{x}}{{\mathbf {e}}_{z}}}&0\\0&1&-{\frac {{\mathbf {e}}_{y}}{{\mathbf {e}}_{z}}}&0\\0&0&1&0\\0&0&1/{\mathbf {e}}_{z}&0\\\end{bmatrix}}{\begin{bmatrix}{\mathbf {d}}_{x}\\{\mathbf {d}}_{y}\\{\mathbf {d}}_{z}\\1\\\end{bmatrix}$

in conjunction with an argument using similar triangles, leads to division by the homogeneous coordinate, giving

${\begin{array}{lcl}\mathbf {b} _{x}&=&\mathbf {f} _{x}/\mathbf {f} _{w}\\\mathbf {b} _{y}&=&\mathbf {f} _{y}/\mathbf {f} _{w}\\\end{array}}.} {\begin{array}{lcl}{\mathbf {b}}_{x}&=&{\mathbf {f}}_{x}/{\mathbf {f}}_{w}\\{\mathbf {b}}_{y}&=&{\mathbf {f}}_{y}/{\mathbf {f}}_{w}\\\end{array}$

There seems to be three options:

1. Implicit homogeneous coordinate

We could leave the w element (I think this may be called the homogeneous coordinate) implicit (and hidden) within P V3 Double and implement some version of !* that does the division automatically

(!*) :: M44 Double -> Point V3 Double -> Point V3 Double

~~But then how can we get both of these matrix multiplications to work? (both of which are useful):~~

(!*.) :: M34 Double -> Point V3 Double -> Point V3 Double

-- and --

(!*.) :: M23 Double -> Point V2 Double -> Point V2 Double

EDIT: On reflection, I suppose that !*. is no different from how !* already works at present (except with the before-mentioned implicit division)... :fireworks:

2. Affine ignores the implicit homogeneous coordinate

Alternatively (with some type-level Nat magic) we could alter Affine so that Diff takes us down one dimension:

(.+^) :: Num a => p a -> Diff p a -> p a
(.-^) :: Num a => p a -> Diff p a -> p a
(.-.) :: Num a => p a -> p a -> Diff p

becomes

(.+^) :: Point V4 Double -> V3 Double -> Point V4 Double
(.-^) :: Point V4 Double -> V3 Double -> Point V4 Double
(.-.) :: Point V4 Double -> Point V4 Double -> V3 Double

but these are a little bit strange because the homogeneous coordinate is sort of just hanging around - I don't remember what happens to the 4th element which is usually 1 in this circumstance.

3. The homogeneous coordinate is explicit, but not part of the underlying type

Or, one could modify Point so that w is not implicit, but not part of the underlying type

data Point f a = P (f a) a

Perhaps Point is a functor such that (.+^) = fmap (^+^).

Please excuse the somewhat disorganized stream of consciousness in my previous comment.

To be clear, my feeling is that 1. Implicit homogeneous coordinate would be the closest thing to what I would expect from Point.

All it would really require is two new operators added to the current interface. Something like this:

(!*.) :: (KnownNat n, r' a ~ V (n + 1) a, r ~ V n a) => m (r' a) -> Point r a -> m a
(.*!) :: (KnownNat n, r' a ~ V (n + 1) a, r ~ V n a) => Point r a -> r' (m a) -> r' a

This has some advantages:

It is still possible to rotate/scale a point using a lower rank matrix.
```
(m :: M33) !* (p :: Point V3 Double)
```
The alternative would be a little bit more cumbersome (and adds unnecessary overhead):
```
m33_to_m44 (m :: M33) !* (p :: Point V3 Double)
```
I think that morally Point might have been intended to abstract away the details of homogeneous coordinates and implies w = 1. This is clear if you look at how Diff is implemented:

type Diff (Point f) = f

In the current interface Point is represented by the same number of components as the underlying euclidean space.

(!*.) as I've presented it has an interesting result: It forces you to handle w explicitly in order to get a Point type out. E.g. Filling in some of the types give you:

(!*.) :: V4 (V4 a) -> Point V3 a -> V4 a
(!*.) :: V3 (V4 a) -> Point V3 a -> V3 a

(.*!) :: Point V3 a -> V4 (V4 a) -> V4 a
(.*!) :: Point V3 a -> V4 (V3 a) -> V3 a

Modifying point and normalizePoint would allow you to turn V4 and V3 back into points:

normalizePoint3 :: V4 a -> Point V3 a
point3 :: V3 a -> Point V3 a

-- You may also want to skip the extra division involved in normalization somehow:
unsafeNormalizePoint3 :: V4 a -> Point V3 a

I threw together a module in our project for type-safe affine transformations. It'd theoretically be possible to create a PR from it, but unfortunately (at least, as far as I was able to make out) it might require a fair amount of surgery on linear to wire in elegantly. Possibly something along the lines of https://github.com/ekmett/linear/compare/DimUnpacked.

In any case, I'll leave this here for posterity:

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}

module Linear.Extra where

import qualified Data.Vector as Vector
import Linear
import Linear.V hiding (Dim)
import Linear.Affine
import qualified Data.Foldable as Foldable
import GHC.TypeLits

type family Vn (n :: Nat) :: * -> * where
  Vn 0 = V0
  Vn 1 = V1
  Vn 2 = V2
  Vn 3 = V3
  Vn 4 = V4
  Vn n = V n

-- Vector dimensions
class VDim f where
  type Dim f :: Nat
  snocV :: f a -> a -> Vn (Dim f + 1) a
  unsnocV :: Vn (Dim f + 1) a -> (f a, a)
  consV :: a -> f a -> Vn (Dim f + 1) a
  unconsV :: Vn (Dim f + 1) a -> (a, f a)
instance VDim V0 where
  type Dim V0 = 0
  snocV V0 x = V1 x
  unsnocV (V1 a) = (V0, a)
  consV x V0 = V1 x
  unconsV (V1 a) = (a, V0)
instance VDim V1 where
  type Dim V1 = 1
  snocV (V1 a) x = V2 a x
  unsnocV (V2 a b) = (V1 a, b)
  consV x (V1 a) = V2 x a
  unconsV (V2 a b) = (a, V1 b)
instance VDim V2 where
  type Dim V2 = 2
  snocV (V2 a b) x = V3 a b x
  unsnocV (V3 a b c) = (V2 a b, c)
  consV x (V2 a b) = V3 x a b
  unconsV (V3 a b c) = (a, V2 b c)
instance VDim V3 where
  type Dim V3 = 3
  snocV (V3 a b c) x = V4 a b c x
  unsnocV (V4 a b c d) = (V3 a b c, d)
  consV x (V3 a b c) = V4 x a b c
  unconsV (V4 a b c d) = (a, V3 b c d)
instance VDim V4 where
  type Dim V4 = 4
  snocV (V4 a b c d) x = V (Vector.fromList [a,b,c,d,x])
  unsnocV (V v) = case Vector.toList v of
                    [a,b,c,d,e] -> (V4 a b c d, e)
                    _ -> error "Incorrect vector dimension in unsnocV"
  consV x (V4 a b c d) = V (Vector.fromList [x,a,b,c,d])
  unconsV (V v) = case Vector.toList v of
                    [a,b,c,d,e] -> (a, V4 b c d e)
                    _ -> error "Incorrect vector dimension in unconsV"

-- type checker needs a helping hand on this one
-- instance VDim (V (n :: Nat)) where
--   type Dim (V n) = n
--   snocV (V v) x = V (Vector.snoc v x)
--   unsnocV (V v) = (V (Vector.init v), Vector.last v)
--   consV x (V v) = V (Vector.cons x v)
--   unconsV (V v) = (Vector.head v, V (Vector.tail v))

-- | Affine transformations on a point
infixl 7 !*.
(!*.) :: (Functor m, VDim v, Num a, Additive (Vn (Dim v + 1)), Foldable (Vn (Dim v + 1))) => m (Vn (Dim v + 1) a) -> Point v a -> m a
m !*. P v = fmap (\r -> Foldable.sum $ liftI2 (*) r (snocV v 1)) m

infixl 7 .*!
(.*!) :: (Num a, Additive f, VDim v, Additive (Vn (Dim v + 1)), Foldable (Vn (Dim v + 1))) => Point v a -> Vn (Dim v + 1) (f a) -> f a
P v .*! g = sumV $ liftI2 (*^) (snocV v 1) g

-- | This is a generic, type-safe version of Linear.V4.normalizePoint
-- https://en.wikipedia.org/wiki/Scaling_(geometry)#Using_homogeneous_coordinates
homogenizePoint :: (Functor f, VDim f, Fractional a) => Vn (Dim f + 1) a -> Point f a
homogenizePoint v =
  case unsnocV v of
    (p,scale) -> P (p ^/ scale)

-- | This is an unsafe version of homogenizePoint which assumes that the point is properly homogenized with the last component = 1
homogeneousPoint :: (Epsilon a, Show a, VDim f) => Vn (Dim f + 1) a -> Point f a
homogeneousPoint v =
  case unsnocV v of
    (p,scale) -> if nearZero (scale - 1)
                 then P p
                 else error ("homogeneousPoint expects the last component to be 1, got " ++ show scale)

-- See also m33_to_m44 :: Num a => M33 a -> M44 a
m22_to_m33 :: Num a => M22 a -> M33 a
m22_to_m33 (V2 v0 v1) = V3 (snocV v0 0) (snocV v1 0) (V3 0 0 1)

-- See also m43_to_m44 :: Num a => M43 a -> M44 a
m23_to_m33 :: Num a => M23 a -> M33 a
m23_to_m33 m = snocV m (V3 0 0 1)

The reasoning behind Linear.Affine is that lots of users of this library actually work with real affine transformations rather than moving to 4x4 homogeneous coordinates. When I'm working with graphics I by-and-large ignore Linear.Affine, and use the usual x,y,z,0 as vector-like homogeneous coordinates (points at infinity) and x,y,z,1 as (normalized) point-like homogeneous coordinates, and just let w float so I can look at things like linear combinations of points as representing lines between them, etc.

However, Linear.Affine lets you work with Point V3 with 3 coefficients, and V3 for differences between points. When you don't want projective transformations included in the mix, and want to use types rather than whether w = 0 to distinguish points from points-at-infinity, some folks prefer to work with affine transformations directly. This is basically like the computer graphics school of using 3x4 matrices for affine transformations knowing the last column (or row depending on handedness, etc.) is 0,0,0,1

The major client of that module is diagrams. If you don't want to subscribe to that view, I'd mostly just ignore that module.

ekmett / linear