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.Vector
import RSAGL.Math.Angle
import RSAGL.Auxiliary.Auxiliary
import RSAGL.Math.Affine
import Data.List
import Data.Maybe
import Control.Parallel.Strategies
import Control.Applicative
import RSAGL.Math.AbstractVector
import Debug.Trace
import RSAGL.Modeling.BoundingBox
import RSAGL.Math.Interpolation
import RSAGL.Math.FMod
import RSAGL.Types
type CurveF a = (RSdouble,RSdouble) -> a
type SurfaceF a = CurveF (CurveF a)
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 ()
sampleCurve :: Curve a -> RSdouble -> RSdouble -> a
sampleCurve (Curve f) = curry f
iterateCurve :: Integer -> Curve x -> [x]
iterateCurve n c = map f $ zeroToOne n
where f = sampleCurve c (0.25/fromInteger n)
zipCurve :: (x -> y -> z) -> Curve x -> Curve y -> Curve z
zipCurve f (Curve x) (Curve y) = Curve $ \hu -> f (x hu) (y hu)
mapCurve :: (CurveF a -> CurveF a) -> Curve a -> Curve a
mapCurve f = Curve . f . fromCurve
mapCurve2 :: (SurfaceF a -> SurfaceF a) -> Curve (Curve a) -> Curve (Curve a)
mapCurve2 f = Curve . (Curve .) . f . (fromCurve .) . fromCurve
translateCurve :: RSdouble -> Curve a -> Curve a
translateCurve x = mapCurve $ \f (h,u) -> f (h,ux)
scaleCurve :: RSdouble -> Curve a -> Curve a
scaleCurve s = mapCurve (\f (h,u) -> f (h/s,u/s))
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)
loopCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a
loopCurve (a,b) | b < a = loopCurve (b,a)
loopCurve (a,b) = mapCurve $ \f (h,u) -> f (h,(ua) `fmod` (ba) + a)
controlCurve :: (RSdouble,RSdouble) -> (RSdouble,RSdouble) -> Curve a -> Curve a
controlCurve (u,v) (u',v') = translateCurve u' . scaleCurve ((u'v') / (uv)) . translateCurve (negate u)
transposeCurve :: Curve (Curve a) -> Curve (Curve a)
transposeCurve = mapCurve2 flip
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
curve :: (RSdouble -> a) -> Curve a
curve = Curve . uncurry . const
newtype Surface a = Surface (Curve (Curve a)) deriving (NFData,AffineTransformable)
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
transposeSurface :: Surface a -> Surface a
transposeSurface (Surface s) = Surface $ transposeCurve s
iterateSurface :: (Integer,Integer) -> Surface a -> [[a]]
iterateSurface (u,v) (Surface s) = map (iterateCurve u) $ iterateCurve v s
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
zipSurface :: (x -> y -> z) -> Surface x -> Surface y -> Surface z
zipSurface f (Surface x) (Surface y) = Surface $ zipCurve (zipCurve f) x y
transformSurface :: (Curve (Curve a) -> Curve (Curve a)) -> Surface a -> Surface a
transformSurface f = Surface . f . unwrapSurface
transformSurface2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Surface a -> Surface a
transformSurface2 fu fv = transformSurface (transformCurve2 fu fv)
translateSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a
translateSurface (u,v) = transformSurface2 (translateCurve u) (translateCurve v)
scaleSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a
scaleSurface (u,v) = transformSurface2 (scaleCurve u) (scaleCurve v)
flipTransposeSurface :: Surface a -> Surface a
flipTransposeSurface = transformSurface (mapCurve2 $ \f (hu,u) (hv,v) -> f (hu,u) (hv,1v)) . transposeSurface
uv_identity :: Surface (RSdouble,RSdouble)
uv_identity = surface (curry id)
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,uh)
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
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
newellCurve :: Curve Point3D -> Maybe Vector3D
newellCurve c = newell $ iterateCurve 16 c
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))
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
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)
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)
intervalRange :: IntervalSample a -> (RSdouble,RSdouble)
intervalRange (IntervalSample _ l h _) = (l,h)
intervalSize :: IntervalSample a -> RSdouble
intervalSize (IntervalSample _ l h _) = abs $ h l
intervalValue :: IntervalSample a -> a
intervalValue (IntervalSample _ _ _ a) = a
intervalSingleIntegral :: (AbstractScale a) => IntervalSample a -> a
intervalSingleIntegral x = scalarMultiply (intervalSize x) $ intervalValue x
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]
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
linearSamples :: Integer -> SamplingAlgorithm a
linearSamples n c = map (\(l,h) -> intervalSample c l h) $ doubles $ zeroToOne (n+1)
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 $ hl)
resampleLoop :: (b -> Maybe b) -> b -> b
resampleLoop nextPass initial_value = f $ initial_value
where f x = maybe x f $ nextPass x