content: {-# LANGUAGE CPP #-}
  ---------------------------------------------------------------------------
  -- |
  -- Copyright   :  (C) 2015 Edward Kmett
  -- License     :  BSD-style (see the file LICENSE)
  --
  -- Maintainer  :  Edward Kmett <ekmett@gmail.com>
  -- Stability   :  experimental
  -- Portability :  non-portable
  --
  -- Common projection matrices: e.g. perspective/orthographic transformation
  -- matrices.
  --
  -- Analytically derived inverses are also supplied, because they can be
  -- much more accurate in practice than computing them through general
  -- purpose means
  ---------------------------------------------------------------------------
  module Linear.Projection
    ( lookAt
    , perspective, inversePerspective
    , infinitePerspective, inverseInfinitePerspective
    , frustum, inverseFrustum
    , ortho, inverseOrtho
    ) where
  
  import Control.Lens hiding (index)
  import Linear.V3
  import Linear.V4
  import Linear.Matrix
  import Linear.Epsilon
  import Linear.Metric
  
  -- $setup
  -- >>> import Linear.Matrix
  -- >>> import Linear.V2
  -- >>> import Linear.V4
  
  -- | Build a look at view matrix
  lookAt
    :: (Epsilon a, Floating a)
    => V3 a -- ^ Eye
    -> V3 a -- ^ Center
    -> V3 a -- ^ Up
    -> M44 a
  lookAt eye center up =
    V4 (V4 (xa^._x)  (xa^._y)  (xa^._z)  xd)
       (V4 (ya^._x)  (ya^._y)  (ya^._z)  yd)
       (V4 (-za^._x) (-za^._y) (-za^._z) zd)
       (V4 0         0         0          1)
    where za = normalize $ center - eye
          xa = normalize $ cross za up
          ya = cross xa za
          xd = -dot xa eye
          yd = -dot ya eye
          zd = dot za eye
  
  -- | Build a matrix for a symmetric perspective-view frustum
  perspective
    :: Floating a
    => a -- ^ FOV (y direction, in radians)
    -> a -- ^ Aspect ratio
    -> a -- ^ Near plane
    -> a -- ^ Far plane
    -> M44 a
  perspective fovy aspect near far =
    V4 (V4 x 0 0    0)
       (V4 0 y 0    0)
       (V4 0 0 z    w)
       (V4 0 0 (-1) 0)
    where tanHalfFovy = tan $ fovy / 2
          x = 1 / (aspect * tanHalfFovy)
          y = 1 / tanHalfFovy
          fpn = far + near
          fmn = far - near
          oon = 0.5/near
          oof = 0.5/far
          -- z = 1 / (near/fpn - far/fpn) -- would be better by .5 bits
          z = -fpn/fmn
          w = 1/(oof-oon) -- 13 bits error reduced to 0.17
          -- w = -(2 * far * near) / fmn
  
  #ifdef HERBIE
  {-# ANN perspective "NoHerbie" #-}
  #endif
  
  -- | Build an inverse perspective matrix
  inversePerspective
    :: Floating a
    => a -- ^ FOV (y direction, in radians)
    -> a -- ^ Aspect ratio
    -> a -- ^ Near plane
    -> a -- ^ Far plane
    -> M44 a
  inversePerspective fovy aspect near far =
    V4 (V4 a 0 0 0   )
       (V4 0 b 0 0   )
       (V4 0 0 0 (-1))
       (V4 0 0 c d   )
    where tanHalfFovy = tan $ fovy / 2
          a = aspect * tanHalfFovy
          b = tanHalfFovy
          c = oon - oof
          d = oon + oof
          oon = 0.5/near
          oof = 0.5/far
  
  
  -- | Build a perspective matrix per the classic @glFrustum@ arguments.
  frustum
    :: Floating a
    => a -- ^ Left
    -> a -- ^ Right
    -> a -- ^ Bottom
    -> a -- ^ Top
    -> a -- ^ Near
    -> a -- ^ Far
    -> M44 a
  frustum l r b t n f =
    V4 (V4 x 0 a    0)
       (V4 0 y e    0)
       (V4 0 0 c    d)
       (V4 0 0 (-1) 0)
    where
      rml = r-l
      tmb = t-b
      fmn = f-n
      x = 2*n/rml
      y = 2*n/tmb
      a = (r+l)/rml
      e = (t+b)/tmb
      c = negate (f+n)/fmn
      d = (-2*f*n)/fmn
  
  inverseFrustum
    :: Floating a
    => a -- ^ Left
    -> a -- ^ Right
    -> a -- ^ Bottom
    -> a -- ^ Top
    -> a -- ^ Near
    -> a -- ^ Far
    -> M44 a
  inverseFrustum l r b t n f =
    V4 (V4 rx 0 0 ax)
       (V4 0 ry 0 by)
       (V4 0 0 0 (-1))
       (V4 0 0 rd cd)
    where
      hrn  = 0.5/n
      hrnf = 0.5/(n*f)
      rx = (r-l)*hrn
      ry = (t-b)*hrn
      ax = (r+l)*hrn
      by = (t+b)*hrn
      cd = (f+n)*hrnf
      rd = (n-f)*hrnf
  
  -- | Build a matrix for a symmetric perspective-view frustum with a far plane at infinite
  infinitePerspective
    :: Floating a
    => a -- ^ FOV (y direction, in radians)
    -> a -- ^ Aspect Ratio
    -> a -- ^ Near plane
    -> M44 a
  infinitePerspective fovy a n =
    V4 (V4 x 0 0    0)
       (V4 0 y 0    0)
       (V4 0 0 (-1) w)
       (V4 0 0 (-1) 0)
    where
      t = n*tan(fovy/2)
      b = -t
      l = b*a
      r = t*a
      x = (2*n)/(r-l)
      y = (2*n)/(t-b)
      w = -2*n
  
  inverseInfinitePerspective
    :: Floating a
    => a -- ^ FOV (y direction, in radians)
    -> a -- ^ Aspect Ratio
    -> a -- ^ Near plane
    -> M44 a
  inverseInfinitePerspective fovy a n =
    V4 (V4 rx 0 0  0)
       (V4 0 ry 0  0)
       (V4 0 0  0  (-1))
       (V4 0 0  rw (-rw))
    where
      t = n*tan(fovy/2)
      b = -t
      l = b*a
      r = t*a
      hrn = 0.5/n
      rx = (r-l)*hrn
      ry = (t-b)*hrn
      rw = -hrn
  
  -- | Build an orthographic perspective matrix from 6 clipping planes.
  -- This matrix takes the region delimited by these planes and maps it
  -- to normalized device coordinates between [-1,1]
  --
  -- This call is designed to mimic the parameters to the OpenGL @glOrtho@
  -- call, so it has a slightly strange convention: Notably: the near and
  -- far planes are negated.
  --
  -- Consequently:
  --
  -- @
  -- 'ortho' l r b t n f !* 'V4' l b (-n) 1 = 'V4' (-1) (-1) (-1) 1
  -- 'ortho' l r b t n f !* 'V4' r t (-f) 1 = 'V4' 1 1 1 1
  -- @
  --
  -- Examples:
  --
  -- >>> ortho 1 2 3 4 5 6 !* V4 1 3 (-5) 1
  -- V4 (-1.0) (-1.0) (-1.0) 1.0
  --
  -- >>> ortho 1 2 3 4 5 6 !* V4 2 4 (-6) 1
  -- V4 1.0 1.0 1.0 1.0
  ortho
    :: Fractional a
    => a -- ^ Left
    -> a -- ^ Right
    -> a -- ^ Bottom
    -> a -- ^ Top
    -> a -- ^ Near
    -> a -- ^ Far
    -> M44 a
  ortho l r b t n f =
    V4 (V4 (-2*x) 0      0     ((r+l)*x))
       (V4 0      (-2*y) 0     ((t+b)*y))
       (V4 0      0      (2*z) ((f+n)*z))
       (V4 0      0      0     1)
    where x = recip(l-r)
          y = recip(b-t)
          z = recip(n-f)
  
  -- | Build an inverse orthographic perspective matrix from 6 clipping planes
  inverseOrtho
    :: Fractional a
    => a -- ^ Left
    -> a -- ^ Right
    -> a -- ^ Bottom
    -> a -- ^ Top
    -> a -- ^ Near
    -> a -- ^ Far
    -> M44 a
  inverseOrtho l r b t n f =
    V4 (V4 x 0 0 c)
       (V4 0 y 0 d)
       (V4 0 0 z e)
       (V4 0 0 0 1)
    where x = 0.5*(r-l)
          y = 0.5*(t-b)
          z = 0.5*(n-f)
          c = 0.5*(l+r)
          d = 0.5*(b+t)
          e = -0.5*(n+f)
  