module Objects where

import Vectors
import Basics
import Data.List(sort)
import Debug.Trace

newtype ProtoObject = ProtoObject { fromProto :: Int -> Object }
data Object = Object !Int !Shape !Texture deriving Show

object sh tex = ProtoObject (\ id -> Object id sh tex)

--objectTexture (Object _ _ material) = material

data Shape 
    = Sphere !Vector !Scal -- center, radius
    | Plane !Vector !Scal  -- normal, distance to origin
    | Quadric !Vector !Vector !Scal  -- center, director, const

{- Equation of (Quadric c d k) is (norm (x-c) == abs ( d `dotProd` (x-c)) + k)
let e = norm d, then
e = 0 => sphere
0 < e < 1 => ellipsoid
e = 1 => cylinder
e > 1 => hyperboloid (or cone, when k = 0)

-}
    | Meta ![MetaPoint] !Scal -- points, threshold
  deriving Show

data MetaPoint = MetaPoint !Vector !Scal -- center, strength
  deriving Show




eps :: Scal
eps = (1.0e-3)

distance src ray (Object id sh m) = distance' src ray sh
getNormal (Object id sh m) = getNormal' sh

distance' source ray (Sphere center rad) = solution
  -- we solve sq t - 2*b*t + * c = 0 (w.r.t. t)
    where b = ray `dotProd` (center - source)
	  c = sqnorm (center - source) - sq rad
	  d = b * b - c
	  solution
	      | d < 0.0 = infinite
	      | sol1 > 0 = sol1
	      | sol2 > 0 = sol2
	      | otherwise = infinite
	  sol1 = b - sqrt d
	  sol2 = b + sqrt d

distance' source ray (Plane norm dist) = solution
 where solution = if abs v1 > eps && ((v > 0.0) == (v1 > 0.0)) then v / v1 else infinite
       v1 = norm `dotProd` ray
       v = dist - norm `dotProd` source


distance' source ray (Quadric center dir sqk0) = solution
    where 
	  a = 1 - sq (ray .* dir)
	  b = ray .* ((c0 .* dir) `scale` dir - c0)
	  c = sqnorm c0 - sqk0 - sq (c0 .* dir)
	  d = b * b - a * c
	  sol1 = (b - sqrt d)/a
	  sol2 = (b + sqrt d)/a
	  sol = if a > 0 then sol1 else sol2
	  -- the good solution depends on the concavity/convexy of the curve
	  solution = if d >= 0 && sol > 0 then sol else infinite
  	  c0 = source - center

distance' source ray (Meta points threshold) = --trace ("Maxs = " ++ show maximums) $ 
					       --trace ("Dens = " ++ show (map density maximums)) $ 
                                               if ubound > 0 then solution else infinite
    where ubound = sum [(strength/(c-b*b)) | (b, c, strength) <- pointInfo] - threshold
	  pointInfo = [(ray .* (center - source),
			sqnorm (center - source) + 0.001,
			strength) |
		       MetaPoint center strength <- points]
	  
	  maximums = [b | (b, _, _) <- pointInfo, b > 0]
	  
	  density x = sum [(strength/(x*x-2*b*x+c)) | (b, c, strength) <- pointInfo] - threshold
	  forwardFindSolution (x0:x1:xs)
	      | density x1 > 0 = dichoFindSolution x0 x1
	      | otherwise = forwardFindSolution (x1:xs)
	  forwardFindSolution (x0:[]) = infinite
	  dichoFindSolution x0 x1 
	      | x1 - x0 < eps = mid
	      | density mid >  0 = dichoFindSolution x0  mid
	      | density mid <= 0 = dichoFindSolution mid x1
	      where mid = (x1 + x0) / 2
	  solution = forwardFindSolution (0:sort maximums)
	      


sq x = x * x

getNormal' (Sphere center rad) hit = (1.0/rad) `scale` (hit - center)
getNormal' (Plane norm dist) hit = norm

getNormal' (Quadric center dir k) hit = normalized $ 
      v - (c `scale` dir)
    where v = hit - center
	  c = dir `dotProd` v 

getNormal' (Meta points threshold) hit = 
    normalized $ sum [ (strength / sq (sqnorm (hit-center))) `scale` 
			(hit-center)
			| MetaPoint center strength <- points ]