{-# LANGUAGE GeneralizedNewtypeDeriving, ExistentialQuantification, Rank2Types #-}
module RSAGL.Math.Curve
    (Curve,
     zipCurve,
     iterateCurve,
     transposeCurve,
     curve,
     Surface,
     surface,
     wrapSurface,
     unwrapSurface,
     transposeSurface,
     zipSurface,
     iterateSurface,
     halfIterateSurface,
     flipTransposeSurface,
     translateCurve,
     scaleCurve,
     clampCurve,
     loopCurve,
     controlCurve,
     transformCurve2,
     uv_identity,
     translateSurface,
     scaleSurface,
     transformSurface,
     transformSurface2,
     surfaceDerivative,
     curveDerivative,
     orientationLoops,
     newellCurve,
     surfaceNormals3D,
     SamplingAlgorithm,
     IntervalSample,
     intervalRange,
     intervalSize,
     intervalSample,
     intervalSingleIntegral,
     linearSamples,
     adaptiveMagnitudeSamples,
     integrateCurve)
    where

import RSAGL.Math.ListUtils
import RSAGL.Math.Vector
import RSAGL.Math.Angle
import RSAGL.Math.Affine
import Data.List
import Data.Maybe
import Control.DeepSeq
import Control.Parallel.Strategies
import Control.Applicative
import RSAGL.Math.AbstractVector
import Debug.Trace
import RSAGL.Math.Interpolation
import RSAGL.Math.FMod
import RSAGL.Math.Types
import RSAGL.Math.BoundingBox

-- | A parametric function that is aware of it's own sampling interval.  The first parameter is the sampling interval, while the second is the curve input parameter.
type CurveF a = (RSdouble,RSdouble) -> a
-- | A surface is a curve of curves.
type SurfaceF a = CurveF (CurveF a)
-- | A 'Curve' is a parametric function from a one-dimensional space into a space of an arbitrary datatype.  The key feature of a 'Curve' is that it is aware of it's own
-- sampling interval.  Using this information and appropriate arithmetic and scalar multiplication functions provided by RSAGL.AbstractVector, a 'Curve' can be differentiated or integrated.
newtype Curve a = Curve { fromCurve :: CurveF a }

instance Functor Curve where
   fmap g (Curve f) = Curve $ g . f

instance Applicative Curve where
    pure a = Curve $ const a
    f <*> a = zipCurve ($) f a

instance (AffineTransformable a) => AffineTransformable (Curve a) where
    scale v = fmap (scale v)
    translate v = fmap (translate v)
    rotate vector angle = fmap (rotate vector angle)
    transform m = fmap (transform m)
    rotateX angle = fmap (rotateX angle)
    rotateY angle = fmap (rotateY angle)
    rotateZ angle = fmap (rotateZ angle)

instance NFData (Curve a) where
    rnf (Curve f) = seq f ()

-- | Sample a specific point on a curve given the sampling interval (first parameter) and the point on the curve (second parameter).
sampleCurve :: Curve a -> RSdouble -> RSdouble -> a
sampleCurve (Curve f) = curry f

-- | Sample a curve at regular intervals in the range 0..1 inclusive.
iterateCurve :: Integer -> Curve x -> [x]
iterateCurve n c = map f $ zeroToOne n
    where f = sampleCurve c (0.25/fromInteger n)

-- | Combine two curves using an arbitrary function.
zipCurve :: (x -> y -> z) -> Curve x -> Curve y -> Curve z
zipCurve f (Curve x) (Curve y) = Curve $ \hu -> f (x hu) (y hu)

-- | Arbitrarily transform a 'Curve'.
mapCurve :: (CurveF a -> CurveF a) -> Curve a -> Curve a
mapCurve f = Curve . f . fromCurve

-- | Arbitrarily transform a surface 'Curve'.
mapCurve2 :: (SurfaceF a -> SurfaceF a) -> Curve (Curve a) -> Curve (Curve a)
mapCurve2 f = Curve . (Curve .) . f . (fromCurve .) . fromCurve

-- | Translate a curve along the axis of the input parameter.
translateCurve :: RSdouble -> Curve a -> Curve a
translateCurve x = mapCurve $ \f (h,u) -> f (h,u-x)

-- | Scale a curve along the axis of the input parameter.  Factors greater than one have a "zoom in" effect, while
-- factors less than one have a "zoom out" effect.
scaleCurve :: RSdouble -> Curve a -> Curve a
scaleCurve s = mapCurve (\f (h,u) -> f (h/s,u/s))

-- | Clamp lower and upper bounds of a curve along the axis of the input parameter.
clampCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a
clampCurve (a,b) | b < a = clampCurve (b,a)
clampCurve (a,b) = mapCurve $ \f (h,u) -> f (h,min b $ max a u)

-- | Loop a curve onto itself at the specified bounds.
loopCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a
loopCurve (a,b) | b < a = loopCurve (b,a)
loopCurve (a,b) = mapCurve $ \f (h,u) -> f (h,(u-a) `fmod` (b-a) + a)

-- | Transform a curve by manipulating control points.
controlCurve :: (RSdouble,RSdouble) -> (RSdouble,RSdouble) -> Curve a -> Curve a
controlCurve (u,v) (u',v') = translateCurve u' . scaleCurve ((u'-v') / (u-v)) . translateCurve (negate u)

-- | Transpose the inner and outer components of a curve.
transposeCurve :: Curve (Curve a) -> Curve (Curve a)
transposeCurve = mapCurve2 flip

-- | Lift two curve transformations onto each axis of a second order curve.
transformCurve2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Curve (Curve a) -> Curve (Curve a)
transformCurve2 fu fv = transposeCurve . fu . transposeCurve . fv

-- | Define a simple curve.
curve :: (RSdouble -> a) -> Curve a
curve = Curve . uncurry . const

-- | A 'Surface' is a based on a 'Curve' with an output of another 'Curve'.
newtype Surface a = Surface (Curve (Curve a)) deriving (NFData,AffineTransformable)

-- | Define a simple surface.
surface :: (RSdouble -> RSdouble -> a) -> Surface a
surface f = Surface $ curve (\u -> curve $ flip f u)

wrapSurface :: Curve (Curve a) -> Surface a
wrapSurface = Surface

unwrapSurface :: Surface a -> Curve (Curve a)
unwrapSurface (Surface s) = s

-- | Transpose the axes of a 'Surface'.
transposeSurface :: Surface a -> Surface a
transposeSurface (Surface s) = Surface $ transposeCurve s

-- | Sample a surface at regularly spaced lattice points in the range 0..1 inclusive.
iterateSurface :: (Integer,Integer) -> Surface a -> [[a]]
iterateSurface (u,v) (Surface s) = map (iterateCurve u) $ iterateCurve v s

-- | Sample the outer 'Curve' of a 'Surface' at regularly spaced intervals.
halfIterateSurface :: Integer -> Surface a -> [Curve a]
halfIterateSurface u = iterateCurve u . unwrapSurface

instance Functor Surface where
    fmap f (Surface x) = Surface $ fmap (fmap f) x

instance Applicative Surface where
    pure a = surface (const $ const a)
    f <*> a = zipSurface ($) f a

-- | Combine two surfaces using an arbitrary function.
zipSurface :: (x -> y -> z) -> Surface x -> Surface y -> Surface z
zipSurface f (Surface x) (Surface y) = Surface $ zipCurve (zipCurve f) x y

-- | Lift a transformation on a second order 'Curve' onto a Surface.
transformSurface :: (Curve (Curve a) -> Curve (Curve a)) -> Surface a -> Surface a
transformSurface f = Surface . f . unwrapSurface

-- | Lift two curve transformations onto each axis of a Surface.
transformSurface2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Surface a -> Surface a
transformSurface2 fu fv = transformSurface (transformCurve2 fu fv)

-- | Translate a surface over each of its input parameter axes, as translateCurve.
translateSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a
translateSurface (u,v) = transformSurface2 (translateCurve u) (translateCurve v)

-- | Scale a surface along each of its input parameter axes, as scaleCurve.
scaleSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a
scaleSurface (u,v) = transformSurface2 (scaleCurve u) (scaleCurve v)

-- | Transpose a surface while flipping the inner curve, so that that orientable surfaces retain their original orientation.
flipTransposeSurface :: Surface a -> Surface a
flipTransposeSurface = transformSurface (mapCurve2 $ \f (hu,u) (hv,v) -> f (hu,u) (hv,1-v)) . transposeSurface

-- | An identity 'Surface'.
uv_identity :: Surface (RSdouble,RSdouble)
uv_identity = surface (curry id)

-- | Take the derivative of a 'Curve'.
curveDerivative :: (AbstractSubtract p v,AbstractScale v) => Curve p -> Curve v
curveDerivative (Curve f) = Curve $ \(h,u) -> scalarMultiply (recip $ 2 * h) $ f (h/2,u+h) `sub` f (h/2,u-h)

-- | Take the piecewise derivative of a 'Surface' along the inner and outer curves.
surfaceDerivative :: (AbstractSubtract p v,AbstractScale v) => Surface p -> Surface (v,v)
surfaceDerivative s = zipSurface (,) (curvewiseDerivative s) (transposeSurface $ curvewiseDerivative $ transposeSurface s)
    where curvewiseDerivative (Surface t) = Surface $ fmap curveDerivative t

-- | Determine the orientation of a 'Surface' by passing very small circles centered on each sampled point as the parametric input.
--
-- A gotchya with this operation is that as much as 3/4ths of the orientation loop may lie outside of the 0..1 range that is normally
-- sampled.  Depending on how the surface is constructed, this may produce unexpected results.  The solution is to clamp the
-- the problematic parametric inputs at 0 and 1 using 'clampSurface'.
--
-- As a rule, do clamp longitudinal axes that come to a singularity at each end.
-- Do not clamp latitudinal axes that are connected at each end.
--
orientationLoops :: Surface p -> Surface (Curve p)
orientationLoops (Surface s) = Surface $ Curve $ \(uh,u) -> Curve $ \(vh,v) ->
                                     curve $ \t -> f (uh/2,u + uh*(sine $ fromRotations t)) 
				                     (vh/2,v + vh*(cosine $ fromRotations t))
   where f = fromCurve . fromCurve s

-- | Try to determine the normal vector to a curve.
newellCurve :: Curve Point3D -> Maybe Vector3D
newellCurve c = newell $ iterateCurve 16 c

-- | Try to determine the normal vectors of a surface using orientation loops.  This is usually slower but more successful than 'surfaceNormals3DByPartialDerivatives'.
-- This generate a warning message if it cannot determine the normal vector at a sampled point.
-- See also 'orientationLoops'.
surfaceNormals3DByOrientationLoops :: Surface Point3D -> Surface SurfaceVertex3D
surfaceNormals3DByOrientationLoops s = SurfaceVertex3D <$> s <*> ((\c -> errmsg c (newellCurve c)) <$> orientationLoops s)
    where errmsg c = fromMaybe (trace ("surfaceNormals3DByOrientationLoops: zero normal gave up: " ++ show (iterateCurve 16 c)) (Vector3D 0 0 0))

-- | Try to determine the normal vectors of a surface using partial derivatives.
surfaceNormals3DByPartialDerivatives :: Surface Point3D -> Surface (Maybe Vector3D)
surfaceNormals3DByPartialDerivatives s = safeCrossProduct <$> surfaceDerivative s
    where x = snd $ boundingCenterRadius $ boundingBox $ concat $ iterateSurface (8,8) s
          safeCrossProduct (u_,v_) =
              do u <- aLargeVector (x/100) u_
	         v <- aLargeVector (x/100) v_
		 return $ vectorNormalize $ crossProduct u v

-- | Try to determine the normal vectors of a surface using multiple techniques.
surfaceNormals3D :: Surface Point3D -> Surface SurfaceVertex3D
surfaceNormals3D s = (\p by_pd by_newell -> case by_pd of
                           Just v -> SurfaceVertex3D p v
			   Nothing -> by_newell) <$>
                     s <*> (surfaceNormals3DByPartialDerivatives s) <*> (surfaceNormals3DByOrientationLoops s)

-- | An interval of a curve, including the curve, lower and upper bounds of the interval, and an instantaneous sample value for that interval.
data IntervalSample a = IntervalSample (Curve a) RSdouble RSdouble a

intervalSample :: Curve a -> RSdouble -> RSdouble -> IntervalSample a
intervalSample c l h = IntervalSample c l h $ sampleCurve c ((abs $ l - h) / 2) ((l+h) / 2) 

-- | Lower and upper bounds of an 'IntervalSample'.
intervalRange :: IntervalSample a -> (RSdouble,RSdouble)
intervalRange (IntervalSample _ l h _) = (l,h)

-- | Size of the range of an 'IntervalSample'.
intervalSize :: IntervalSample a -> RSdouble
intervalSize (IntervalSample _ l h _) = abs $ h - l

-- | Instantaneous sample value of an 'Interval'.
intervalValue :: IntervalSample a -> a
intervalValue (IntervalSample _ _ _ a) = a

-- | Integral of the sample value over the range of the 'IntervalSample'.
intervalSingleIntegral :: (AbstractScale a) => IntervalSample a -> a
intervalSingleIntegral x = scalarMultiply (intervalSize x) $ intervalValue x

-- | Split an interval into three equal parts.
splitInterval :: IntervalSample a -> [IntervalSample a]
splitInterval (IntervalSample c l h a) = [intervalSample c l l',IntervalSample c l' h' a,intervalSample c h' h]
    where l' = lerp (1/3) (l,h)
          h' = lerp (2/3) (l,h)

type SamplingAlgorithm a = Curve a -> [IntervalSample a]

-- | Definite integral of a curve.
integrateCurve :: (AbstractAdd p v,AbstractScale v,AbstractZero p) => SamplingAlgorithm v -> Curve v -> p -> p
integrateCurve samplingAlgorithm c initial_value = foldl' add initial_value $ map intervalSingleIntegral $ samplingAlgorithm c

-- | Sampling algorithm that takes a fixed count of samples.
linearSamples :: Integer -> SamplingAlgorithm a
linearSamples n c = map (\(l,h) -> intervalSample c l h) $ doubles $ zeroToOne (n+1)

-- | Sampling algorithm that takes increasing numbers of samples over intervals where the magnitude of the sample is large.
adaptiveMagnitudeSamples :: (AbstractMagnitude a) => Integer -> SamplingAlgorithm a
adaptiveMagnitudeSamples n c = resampleLoop (\xs -> if genericLength xs > n then Nothing else Just $ newSamples xs) $ linearSamples (max 1 $ n `div` 10) c
    where newSamples xs = let a = abstractAverage $ map intervalMagnitude xs
                              in flip concatMap xs $ \x -> if intervalMagnitude x >= a then splitInterval x else [x]
          intervalMagnitude :: (AbstractMagnitude a) => IntervalSample a -> RSdouble
          intervalMagnitude (IntervalSample _ l h a) = magnitude a * (abs $ h-l)

-- | Loop to keep generating samples until finished.
resampleLoop :: (b -> Maybe b) -> b -> b
resampleLoop nextPass initial_value = f $ initial_value
    where f x = maybe x f $ nextPass x