{- 	
--- Programming competition entry - Karol Pogonowski, matriculation number: 0838671 ---
	This is the raytracer I wrote for the competition. It supports quite a lot of 
	cool techniques, including:
	-Phong shading
	-Soft shadows
	-Adaptive supersampling antialiasing
	-Multi-depth refraction and reflection
	-Transparency shadows
	-High Dynamic Range rendering with bloom
	-Parallel computing
	
	In general, this is not optimized a lot - it should have a BSP tree construction
	to reduce amount of time spend on raytracing. However, it works good without it too,
	just need some time to render whole screen ;)
	
	To compile it, run > ghc -O2 -threaded --make raytracer.hs <. Do not try to run this 
	in ghci, as it will probably crash very quickly. There are two renderers - default one
	based on lazy lists and a HDR one, which works on Data.Map. Although the first one runs 
	on very low memory usage, the second is implemented quite poorly (should be using mutable
	unboxed arrays instead) and can eat quite a lot of RAM.
	
	To invoke the HDR version, run > ./raytracer hdr +RTS -N<number of CPU cores> <. Be sure 
	you run this on computer with at least 2GB RAM. If the program crashes due to being out of
	memory in the middle of computation, add '-c' flag to RTS options, which will allow
	compaction GC algorithm to be run and optimize the RAM usage, sacrificing a bit of performance.
	
	Non-HDR renderer is invoked by > ./raytracer +RTS -N<number of CPU cores> <. It doesn't require
	that much RAM, in fact it eats no more than 10 MB of RAM due to lazy evaluation.
	
	My competition entry would be the HDR rendered one, however if you want to see how an ordinary
	one looks (or HDR renderer would fail, which should not happen though) feel free to run it and
	compare the results.
	
	This was tested and working on GHC 6.10.1 on Windows XP SP3 and GHC 6.8.1 on student.compute
	servers. The render time is around 20 minutes without compression algorithm on my Core 2 Duo 
	2.3 Ghz.
-}


import System.IO
import System.Environment
import Control.Parallel.Strategies
import qualified Data.List as List
import qualified Data.Map as Map
import Data.Word
import Data.Time
import Text.Printf

{----- Datatype definitions -----}

-- Vector class, double precision
data VecT = VecT { vecX, vecY, vecZ :: Double }
-- Colours are lower precision than vectors
data ColourT = ColT { r, g, b :: Float } deriving (Eq)
-- Ray data
data RayT = RayT { rayOr, rayDir :: VecT }
-- Object primitives
data ObjectT 	= Sphere { objOr :: VecT, rad :: Double, mat :: MatT, isLight :: Bool, name :: String }
				| Plane { objOr :: VecT, distance :: Double, mat :: MatT, isLight :: Bool, name :: String }-- objOr is plane's normal
-- Material type
data MatT = MatT { diffuse, reflection, specular, refraction, transparency :: Double, clr :: ColourT }
-- Hit return
data HitT = MISS | HIT | INOBJ deriving (Enum,Eq)

instance NFData ColourT where
	rnf (ColT r g b) = r `seq` g `seq` b `seq` ()
	
{----- Global render defines -----}

-- This sets field of view from the origin. 
-- Note: keep aspect ratio here
renderBoundX = -16.0
renderBoundY = -10.0

-- Get delta x for interpolation	between pixels
getDX = (renderBoundX * (-2)) / fromIntegral screenWidth
-- Get delta y for interpolation	between pixels			
getDY = (renderBoundY * (-2)) / fromIntegral screenHeight

-- Origin of the camera. Keep this relatively far away.
origin = (VecT 0.0 (0.0) (-50.0))

-- Resolution
screenWidth = 1440 :: Word16
screenHeight = 900 :: Word16
-- Resolution for HDR - 1
hdrHeight = screenHeight `div` 2 - 1
hdrWidth = screenWidth `div` 2 - 1	

-- Max trace distance
maxDist = 1000.0 :: Double
-- Maximum recursive depth for traces
maxDepth = 6 ::Int
-- Small value used for finding a point just above surface of the object
epsilon = 0.0001 :: Double
-- Maximum depth of transparency shading trace
shadowDepth = 2 :: Int
-- Sets number of passes for x-axis for soft shadows. To look smooth this has to be at least 1/100 of resolution
shadowPasses = 16
-- Number of supersample anti-aliasing passes (Always multiple of 2!)
superPasses = 8
-- Intensity of bloom in HDR renderer
bloomInt = 0.2
-- Pixel bloom passes. Note this happens on 1/4 of original screen
bloomPasses = 8

-- Zero vector
zeroVector = (VecT 0.0 0.0 0.0)
-- Background colour
bgColor = black
-- Background object
background = (Plane zeroVector maxDist (MatT 0.0 0.0 0.0 0.0 0.0 bgColor) False "background")

--  Predefined colours for convience
red    	 	= (ColT 1.0 0.0 0.0)
green  		= (ColT 0.0 1.0 0.0)
blue   	 	= (ColT 0.0 0.0 1.0)
white  	 	= (ColT 1.0 1.0 1.0)
black   	= (ColT 0.0 0.0 0.0)
midgrey 	= (ColT 0.5 0.5 0.5)
nearlywhite = (ColT 0.8 0.8 0.8)
nearlyblack = (ColT 0.05 0.05 0.05)

{----- Object list definitions -----}

-- List of all objects
objects = [light1,light2,backPlane] ++ lambda 0.0 2.0 0.0 ++
		(h (-9.0) (-2.0)) ++ (a (-5.5) (-2.0)) ++ (s (-1.6) (-2.0)) ++ (k 1.0 (-2.0)) ++ (e 4.0 (-2.0)) ++ (l 7.0 (-2.0)) ++ (l 10.0 (-2.0)) ++
		(p (-6.9) 3.0) ++ (o (-2.9) 3.0) ++ (w (1.6) 3.0) ++ (e 4.1 3.0) ++ (r_ 7.3 3.0) 
lights = filter isLight objects

--	Object type	Origin	Radius	Diffuse Reflection Specular Refraction index Transparency Colour	isLight	Name --

-- Lights
light1 = (Sphere (VecT (10.0) (-6.5) (-4.0)) 	0.2 (MatT 0.0 0.0 0.0 0.0 0.0 (ColT 2.2 2.2 2.2)) True "light1")
light2 = (Sphere (VecT (-10.0) (-6.5) (-4.0)) 	0.2 (MatT 0.0 0.0 0.0 0.0 0.0 (ColT 2.2 2.2 2.2)) True "light2")


-- back plane
backPlane = (Plane (VecT 0.0 (0.0) (-1.0)) 	 5.0 (MatT 0.8 0.5 0.0 0.0 0.0 nearlyblack) False "backPlane")

{----- Letter constructors -----}
-- Build lambda sign from small spheres
lambda x y z = [lam_l1,lam_l2,lam_l3,lam_l4,lam_l5,lam_l6,lam_l7,lam_l8,lam_l9,lam_l10,lam_l11,lam_l12,lam_l13,lam_l14,lam_l15,lam_l16,lam_l17,lam_l18,lam_l19,lam_r1,lam_r2,lam_r3,lam_r4,lam_r5,lam_r6,lam_r7]
	where
	lam_l1 = (Sphere (VecT (x-8.2) (y-6.7) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l1")
	lam_l2 = (Sphere (VecT (x-8.6) (y-5.5) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l2")
	lam_l3 = (Sphere (VecT (x-7) (y-7) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l3")
	lam_l4 = (Sphere (VecT (x-6.2) (y-6) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l4")
	lam_l5 = (Sphere (VecT (x-5.4) (y-5) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l5")
	lam_l6 = (Sphere (VecT (x-4.6) (y-4) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l6")
	lam_l7 = (Sphere (VecT (x-3.8) (y-3) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l7")
	lam_l8 = (Sphere (VecT (x-3) (y-2) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l8")
	lam_l9 = (Sphere (VecT (x-2.2) (y-1) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l9")
	lam_l10 = (Sphere (VecT (x-1.4) (y) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l10")
	lam_l11 = (Sphere (VecT (x-0.6) (y+1) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l11")
	lam_l12 = (Sphere (VecT (x+0.2) (y+2) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l12")
	lam_l13 = (Sphere (VecT (x+1) (y+3) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l13")
	lam_l14 = (Sphere (VecT (x+1.8) (y+4) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l14")
	lam_l15 = (Sphere (VecT (x+2.6) (y+5) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l15")
	lam_l16 = (Sphere (VecT (x+3.4) (y+6) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l16")
	lam_l17 = (Sphere (VecT (x+4.2) (y+7) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l17")
	lam_l18 = (Sphere (VecT (x+5.4) (y+6.7) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l18")
	lam_l19 = (Sphere (VecT (x+5.8) (y+5.5) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_l19")

	lam_r1 = (Sphere (VecT (x-2.0) (y+1) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r1")
	lam_r2 = (Sphere (VecT (x-2.8) (y+2) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r2")
	lam_r3 = (Sphere (VecT (x-3.6) (y+3) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r3")
	lam_r4 = (Sphere (VecT (x-4.4) (y+4) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r4")
	lam_r5 = (Sphere (VecT (x-5.2) (y+5) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r5")
	lam_r6 = (Sphere (VecT (x-6.0) (y+6) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r6")
	lam_r7 = (Sphere (VecT (x-6.8) (y+7) z) 0.6 (MatT 0.6 0.5 0.8 0.0 0.0 (ColT 0.01 0.01 0.01)) False "lam_r7")

-- H
h x y = [h_l1,h_l2,h_l3,h_l4,h_l5,h_r1,h_r2,h_r3,h_r4,h_r5,h_b1,h_b2]
	where
	h_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_l1")
	h_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_l2")
	h_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_l3")
	h_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_l4")
	h_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_l5")

	h_b1 = (Sphere (VecT (x-0.6) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_b1")
	h_b2 = (Sphere (VecT (x+0.2) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_b2")

	h_r1 = (Sphere (VecT (x+1.0) (y+0) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_r1")
	h_r2 = (Sphere (VecT (x+1.0) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_r2")
	h_r3 = (Sphere (VecT (x+1.0) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_r3")
	h_r4 = (Sphere (VecT (x+1.0) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_r4")
	h_r5 = (Sphere (VecT (x+1.0) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.58 0.0 0.82)) False "h_r5")

-- A
a x y = [a_r1,a_r2,a_r3,a_r4,a_l1,a_l2,a_l3,a_l4,a_l5,a_m1,a_m2]
	where
	a_l1 = (Sphere (VecT (x) (y+0.1) 	 (-2.0)) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_l1")
	a_l2 = (Sphere (VecT (x-0.35) (y+0.9) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_l2")
	a_l3 = (Sphere (VecT (x-0.7) (y+1.65) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_l3")
	a_l4 = (Sphere (VecT (x-1.05) (y+2.5) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_l4")
	a_l5 = (Sphere (VecT (x-1.4) (y+3.25) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_l5")
	
	a_m1 = (Sphere (VecT (x-0.35) (y+2.35) (-2.0)) 0.35 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_m1")
	a_m2 = (Sphere (VecT (x+0.35) (y+2.35) (-2.0)) 0.35 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_m2")
	
	a_r1 = (Sphere (VecT (x+0.35) (y+0.9) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_r1")
	a_r2 = (Sphere (VecT (x+0.7) (y+1.65) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_r2")
	a_r3 = (Sphere (VecT (x+1.05) (y+2.5) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_r3")
	a_r4 = (Sphere (VecT (x+1.4) (y+3.25) (-2.0)) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.94 0.64 0.1)) False "a_r4")

-- S
s x y = [s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8]
	where
	s_1 = (Sphere (VecT (x) (y+0.05) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_1")
	s_2 = (Sphere (VecT (x-0.8) (y+0.1) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_2")
	s_3 = (Sphere (VecT (x-1.4) (y+0.65) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_3")
	s_4 = (Sphere (VecT (x-0.9) (y+1.35) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_4")
	s_5 = (Sphere (VecT (x-0.4) (y+2.05) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_5")
	s_6 = (Sphere (VecT (x) 	(y+2.8) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_6")
	s_7 = (Sphere (VecT (x-0.65) (y+3.15) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_7")
	s_8 = (Sphere (VecT (x-1.4) (y+3.35) 	 (-2.0)) 0.37 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 0.0 0.0 0.93)) False "s_8")
	
-- K
k x y = [k_l1,k_l2,k_l3,k_l4,k_l5,k_r1,k_r2,k_r3,k_r4,k_r5]
	where
	k_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_l1")
	k_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_l2")
	k_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_l3")
	k_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_l4")
	k_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_l5")
	
	k_r1 = (Sphere (VecT (x+0.6) (y+0.1) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_r1")
	k_r2 = (Sphere (VecT (x+0.0) (y+0.8) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_r2")
	k_r3 = (Sphere (VecT (x-0.6) (y+1.6) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_r3")	
	k_r4 = (Sphere (VecT (x+0.0) (y+2.35) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_r4")
	k_r5 = (Sphere (VecT (x+0.6) (y+3.1) ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.75 0.24 1.0)) False "k_r5")
	
-- E
e x y = [e_l1,e_l2,e_l3,e_l4,e_l5,e_t1,e_t2,e_m1,e_m2,e_b1,e_b2]
	where
	e_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_l1")
	e_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_l2")
	e_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_l3")
	e_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_l4")
	e_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_l5")

	e_t1 = (Sphere (VecT (x-0.6) (y) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_t1")
	e_t2 = (Sphere (VecT (x+0.2) (y) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_t2")
	
	e_m1 = (Sphere (VecT (x-0.6) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_m1")
	e_m2 = (Sphere (VecT (x+0.2) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_m2")
	
	e_b1 = (Sphere (VecT (x-0.6) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_b1")
	e_b2 = (Sphere (VecT (x+0.2) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.9 0.0 0.0)) False "e_b2")
	
-- L
l x y = [l_l1,l_l2,l_l3,l_l4,l_l5,l_b1,l_b2]
	where
	l_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_l1")
	l_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_l2")
	l_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_l3")
	l_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_l4")
	l_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_l5")
	
	l_b1 = (Sphere (VecT (x-0.6) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_b1")
	l_b2 = (Sphere (VecT (x+0.2) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.5 1.0 0.62)) False "l_b2")

-- P
p x y = [p_o1,p_o2,p_o3,p_l1,p_l2,p_l3,p_l4,p_l5]
	where
	p_o1 = (Sphere (VecT (x-0.6) (y) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_o1")
	p_o2 = (Sphere (VecT (x+0.1) (y+0.8) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_o2")
	p_o3 = (Sphere (VecT (x-0.6) (y+1.6) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_o3")
	
	p_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_l1")
	p_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_l2")
	p_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_l3")
	p_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_l4")
	p_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 1.0 0.2 0.0)) False "p_l5")

-- O
o x y = [(Sphere (VecT (x-1.6) (y+1.6) 	 ((-2.0))) 1.8 (MatT 0.4 0.3 1.0 1.3 0.9 (ColT 0.1 0.1 0.1)) False "r_l1")]

-- W
w x y = [w_r1,w_r2,w_r3,w_r4,w_m,w_l1,w_l2,w_l3,w_l4]
	where
	w_r1 = (Sphere (VecT (x) (y) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_r1")
	w_r2 = (Sphere (VecT (x-0.4) (y+1.0) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_r3")
	w_r3 = (Sphere (VecT (x-0.8) (y+2) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_r3")
	w_r4 = (Sphere (VecT (x-1.2) (y+3) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_r4")

	w_m = (Sphere (VecT (x-2.0) (y+2.3) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_r4")
	
	w_l1 = (Sphere (VecT (x-2.8) (y+3) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_l1")
	w_l2 = (Sphere (VecT (x-3.2) (y+2.0) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_l3")
	w_l3 = (Sphere (VecT (x-3.6) (y+1) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_l3")
	w_l4 = (Sphere (VecT (x-4.0) (y) 	 ((-2.0))) 0.5 (MatT 0.4 0.3 1.0 1.05 0.7 (ColT 1.0 0.61 0.71)) False "w_l4")
	
-- R
r_ x y = [r_l1,r_l2,r_l3,r_l4,r_l5,r_r1,r_r2,r_o1,r_o2,r_o3]
	where
	r_l1 = (Sphere (VecT (x-1.4) (y) 	 ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_l1")
	r_l2 = (Sphere (VecT (x-1.4) (y+0.8) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_l2")
	r_l3 = (Sphere (VecT (x-1.4) (y+1.6) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_l3")
	r_l4 = (Sphere (VecT (x-1.4) (y+2.4) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_l4")
	r_l5 = (Sphere (VecT (x-1.4) (y+3.2) ((-2.0))) 0.375 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_l5")
	
	r_o1 = (Sphere (VecT (x-0.6) (y) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_o1")
	r_o2 = (Sphere (VecT (x+0.1) (y+0.8) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_o2")
	r_o3 = (Sphere (VecT (x-0.6) (y+1.6) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_o3")
			
	r_r1 = (Sphere (VecT (x-0.55) (y+2.45) 	 ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_r1")
	r_r2 = (Sphere (VecT (x+0.1) (y+3.15) ((-2.0))) 0.4 (MatT 0.4 0.3 1.0 0.0 0.0 (ColT 0.0 0.6 0.0)) False "r_r2")

{----- Math functions -----}
-- We have a right handed coordinate system.  If x increases to your right, and Y increases downwards then
-- you are looking in the direction of increasing Z.

-- Colour addition
colAdd :: ColourT -> ColourT -> ColourT
colAdd (ColT a b c) (ColT x y z) = ColT (a + x) (b + y) (c + z)

-- Colour subtraction
colSub :: ColourT -> ColourT -> ColourT
colSub (ColT a b c) (ColT x y z) = ColT (a - x) (b - y) (c - z)

-- Product of two colours
colProduct :: ColourT -> ColourT -> ColourT
colProduct (ColT a b c) (ColT x y z) = ColT (a*x) (b*y) (c*z)

-- Colour dot product
colDot :: ColourT -> ColourT -> Float
colDot (ColT a b c) (ColT x y z) = (a*x) + (b*y) + (c*z)

-- Colour scalar multiplication
colScalar :: ColourT -> Double -> ColourT
colScalar (ColT a b c) d 
	= ColT (a*d') (b*d') (c*d')
	where
	d' = realToFrac d

-- Returns normal given incident vector
normalVec :: ObjectT -> VecT -> VecT
normalVec (Sphere objOr rad mat isLight name) inc = (inc `vecSub` objOr) `scalarMult` (1.0/rad)
normalVec (Plane normal dist mat isLight name) inc = normal
				
-- Returns vector length
vecLength :: VecT -> Double
vecLength (VecT x y z) = sqrt((x*x) + (y*y) + (z*z))

-- Vector addition
vecAdd :: VecT -> VecT -> VecT
vecAdd (VecT a b c) (VecT x y z) = VecT (a + x) (b + y) (c + z)

-- Vector subtraction
vecSub :: VecT -> VecT -> VecT
vecSub (VecT a b c) (VecT x y z) = VecT (a - x) (b - y) (c - z)


{-
-- Vector cross product
crossProduct :: VecT -> VecT -> VecT
crossProduct (VecT a b c) (VecT x y z) = VecT (b*z + c*y) (-(a*z + c*x)) (a*y - b*a)

-- Negates vector
negVec :: VecT -> VecT
negVec (VecT x y z) = VecT (-x) (-y) (-z)
-}

-- Vector dot product
dotProduct :: VecT -> VecT -> Double
dotProduct (VecT a b c) (VecT x y z) = (a*x) + (b*y) + (c*z)

-- Vector scalar multiplication
scalarMult :: VecT -> Double -> VecT 
scalarMult (VecT a b c) d = VecT (a*d) (b*d) (c*d)

-- Returns unit vector
normalize :: VecT -> VecT
normalize v 	
	| len > 0 = v `scalarMult` (1 / len)
	| otherwise = zeroVector
	where
	len = vecLength v
	
{----- Core raytracing functions -----}
									
-- Naive raytracing. Intersects the ray with every object, then checks which was is closest and return it with its colour.
-- Input: Ray, depth, refraction index
-- Output: (Hit object, colour, distance)
raytrace :: RayT -> Int -> Double -> (ObjectT, ColourT, Double)
raytrace ray depth rIdx	
	| isLight hitObject == True = (hitObject, matClr, dist)
	| otherwise = (hitObject, colour, dist)
	where
	-- Get closest object
	(hitObject, dist, hitType) = rayObjects (background, maxDist, MISS) ray objects
	pi = rayOr ray `vecAdd` (rayDir ray `scalarMult` dist)	
	-- Calculate colour of the pixel
	colour = diffSpecColour `colAdd` reflColour `colAdd` refrColour
	diffSpecColour = calcColour pi hitObject ray matClr lights
	reflColour = calcReflections pi nPi rDir (reflection hitMat) matClr rIdx depth
	refrColour = calcRefraction pi nPi rDir rIdx (refraction hitMat) matClr (hitToDouble hitType) trans depth
	hitMat = mat hitObject
	matClr = clr hitMat
	rDir = rayDir ray
	nPi = normalVec hitObject pi
	trans = transparency hitMat

-- Checks whether specified ray intersect with the object
intersect :: RayT -> ObjectT -> (HitT, Double)
intersect ray (Sphere objOr rad mat isLight name) 	
	| det > 0 && i2 > 0 = if i1 < 0 then (INOBJ, i2) else (HIT, i1)
	| otherwise = (MISS, 0.0)
	where
	dist = rayOr ray `vecSub` objOr
	dot = -(dist `dotProduct` rayDir ray)
	det = (dot*dot) - (dist `dotProduct` dist) + (rad * rad)
	detroot = sqrt det
	i1 = dot - detroot
	i2 = dot + detroot
								
intersect ray (Plane normal distance mat isLight name) 	
	| det /= 0 && dist > 0 = (HIT, dist)
	| otherwise = (MISS, 0.0)	
	where
	det = normal `dotProduct` rayDir ray
	dist = -((normal `dotProduct` rayOr ray) + distance)/det

	
-- Returns numeric value of provided HitT
-- Input: hit data
-- Output: Double value
hitToDouble :: HitT -> Double
hitToDouble (INOBJ) = (-1.0)
hitToDouble (MISS) 	= 0.0
hitToDouble (HIT) 	= 1.0

-- Calculates colour by taking into account all lights illuminating the point
-- Uses diffuse shading and specular lighting
-- Input: point at intersection, object to be illuminated, ray, first colour (for recursion), list of all objects
-- Output: colour
calcColour :: VecT -> ObjectT -> RayT -> ColourT -> [ObjectT] ->  ColourT
calcColour pi hit ray col [] 			= col
calcColour pi hit ray col (lgt:lgts) 	
	= calcColour pi hit ray hitColour lgts
	where
	lightDir = objOr lgt `vecSub` pi
	normalPi = normalVec hit pi
	diffColour = if diff > 0 && shade > 0 then (calcDiffuse lightDir normalPi diff col lgtColour) `colScalar` shade else black
	specColour = if spec > 0 && shade > 0 then (calcSpecular lightDir lgtColour (rayDir ray) normalPi spec) `colScalar` shade else black
	shade = calcShadowSphere pi lgt
	spec = specular (mat hit)
	diff = diffuse (mat hit)
	hitColour = col `colAdd` diffColour `colAdd` specColour
	lgtColour = clr (mat lgt)
							
-- Calculates diffuse shading
-- Input: light direction from PoI, normal at point of intersection, diffuse index, material colour, light colour	
-- Output: colour after diffusion				
calcDiffuse :: VecT -> VecT -> Double -> ColourT -> ColourT -> ColourT
calcDiffuse lDir nPi diffIndex col lClr 
	| dot > 0 = dCol
	| otherwise = black
	where
	l = normalize lDir
	dot = nPi `dotProduct` l
	diff = dot * diffIndex
	dCol = (col `colScalar` diff) `colProduct` lClr

					
-- Calculates specular component
-- Input: light direction from PoI, light colour, ray direction,  normal at point of intersection, specular index				
-- Output: colour after specular
calcSpecular :: VecT -> ColourT -> VecT -> VecT -> Double -> ColourT
calcSpecular lDir lClr rdir n specIndex 	
	| dot > 0 = sCol
	| otherwise = black
	where
	l = normalize (lDir)
	q = l `vecSub` (n `scalarMult` (2 * (l `dotProduct` n)))
	dot = rdir `dotProduct` q
	spec = specIndex * (dot ^ 20)
	sCol = lClr `colScalar` spec
	
-- Calculates whether specified object is in shadow from sphere light source	
-- This outputs soft shadows, although they are only softened along x-axis.
-- Input: point of intersection, light
-- Output: magnitude of shade
calcShadowSphere :: VecT -> ObjectT -> Double
calcShadowSphere pi lgt 	
	= (shadowPass 0) / shadowPasses
	where		
	radius = rad lgt
	step = (VecT 1.0 0.0 0.0) `scalarMult` ((radius * 2) / (shadowPasses-1))	
	leftSide = objOr lgt `vecSub` ((VecT 1.0 0.0 0.0) `scalarMult` radius)							
	shadowPass pass	
		| pass < shadowPasses = passShade + shadowPass (pass+1)
		| otherwise = 0
		where
		passOr = leftSide `vecAdd` (step `scalarMult` pass)
		rayDir = passOr `vecSub` pi
		nDir = rayDir `scalarMult` (1/dirLen)
		dirLen = vecLength rayDir
		shadowRay = RayT (pi `vecAdd` (nDir `scalarMult` epsilon)) nDir
		dist = dirLen * (1.0 - (4*epsilon))
		passShade = rayShadow shadowRay dist objectsNoLights 0
		objectsNoLights = filter (\f -> not (isLight f)) objects
		-- Check if shadow ray hits any of the objects that are in our way to the light. Return 1 for fully visible objects and 0 for fully shaded ones
		-- Take transparency of occluding object into account. Note this takes into account only first hit object for now.
		rayShadow ray dist [] dep = 1.0
		rayShadow ray dist (obj:objs) dep 	
			| dep < shadowDepth && rslt == HIT && dst < dist = transparency (mat obj) * rayShadow ray dist objs (dep+1)
			| otherwise = rayShadow ray dist objs dep
			where
			(rslt,dst) = intersect ray obj		
			
-- Calculates refracted colour recursively.
-- Input:  point of intersection, normal at PoI, ray direction, previous refraction index, refraction index, material colour, hit enum, transparency, depth
-- Output: colour after refraction
calcRefraction :: VecT -> VecT -> VecT -> Double -> Double -> ColourT -> Double -> Double -> Int -> ColourT
calcRefraction pi nPi rDir oldIndex refrIndex clr hitT t depth	
	| t > 0 && refrIndex > 0 && depth <= maxDepth && cosT2 > 0 && refrDst < maxDist = (refrColour `colProduct` transparency `colScalar` t) 
	| otherwise = black
	where
	n = oldIndex / refrIndex
	rPi = nPi `scalarMult` hitT
	cosI = -(rPi `dotProduct` rDir)
	cosT2 = 1.0 - (n * n * (1.0 - (cosI * cosI)))
	refrDir = (rDir `scalarMult` n) `vecAdd` (rPi `scalarMult` (n * cosI - sqrt cosT2))
	refrRay = RayT (pi `vecAdd` (refrDir `scalarMult` epsilon)) refrDir
	(refrObj,refrColour,refrDst) = raytrace refrRay (depth+1) refrIndex
	-- Apply Beer's law
	absorbance = clr `colScalar` (0.15 * (-refrDst))
	transparency = (ColT (exp (r absorbance)) (exp (g absorbance)) (exp (b absorbance)))													
																			
-- Calculates reflected colour recursively
-- Input: normal at point of intersection, point of intersection, ray direction, reflection index, material colour, refraction index, depth
-- Output: colour after reflection
calcReflections :: VecT -> VecT -> VecT -> Double -> ColourT -> Double -> Int -> ColourT
calcReflections pi nPi rDir refl col rIdx depth	
	| refl > 0 && depth <= maxDepth && reflDist < maxDist = rCol -- Don't reflect the background
	| otherwise = black
	where		
	(reflObj,reflColour,reflDist) = raytrace reflRay (depth+1) rIdx
	reflVec = rDir `vecSub` (nPi `scalarMult` (2.0 * (rDir `dotProduct` nPi)))
	reflRay = RayT (pi `vecAdd` (reflVec `scalarMult` epsilon)) reflVec				
	rCol = col `colProduct` reflColour `colScalar` refl
		
-- Find nearest intersection by recursing through every object
-- Input: (old object, old distance, old hittype), ray, list of all objects
-- Output: object hit, distance, hit type
rayObjects :: (ObjectT, Double, HitT) -> RayT -> [ObjectT] -> (ObjectT, Double, HitT)
rayObjects (oldObj, oldDist, oldHit) ray [] = (oldObj, oldDist, oldHit)
-- Check if ray hits specified object. Then check whether the specified object is closer than old one.
-- If yes, return it. If not, return to old one. Either way, recurse.
rayObjects (oldObj, oldDist, oldHit) ray (obj:objs) 
	| rslt /= MISS && dst < oldDist = rayObjects (obj, dst, rslt) ray objs
	| otherwise = rayObjects (oldObj, oldDist, oldHit) ray objs
	where
	(rslt,dst) = intersect ray obj

{----- Main rendering functions -----}
-- Runs raytraces for every line, starting with input integer	
-- This runs every row computation in parallel, which is then dynamically
-- allocated to idle cores by the compiler. Note this is anti-aliased 
-- by multisampling, however only in row, as parallelism doesnt allow
-- passing last hit object through rows
-- Output: map with (y,x) coordinates paired with pixel colour
render :: [ColourT]
render =  concat (parMap rnf (renderRow 0 origin background) [0..(screenHeight-1)])

-- Code below is for adaptive multisampling anti-aliasing. 
-- Returns a list of colours of all pixels in the line X		
-- Input: X line number, view origin, Y line number,  last hit primitive
-- Output: lsit of succeeding colours			
renderRow :: Word16 -> VecT -> ObjectT -> Word16 -> [ColourT]
renderRow lineX or lstObj lineY 
	| lineX < screenWidth = colour : (renderRow (lineX+1) or obj lineY)
	| otherwise = []
	where
	sX = renderBoundX + (fromIntegral lineX * getDX)
	sY = renderBoundY + (fromIntegral lineY * getDY)
	rayDir = normalize ((VecT sX sY 0.0) `vecSub` or)
	ray = RayT or rayDir
	(obj,clr,dst) = raytrace ray 0 1.0
	-- Do supersampling only when new primitive is hit
	colour = if (name lstObj) /= (name obj) then renderSupersample sX sY or else clr
												
-- Fires multiple rays for supersampling
-- Returns averaged colours based on number of rays
-- Input: sX sY origin
-- Output: averaged colour, hit object
renderSupersample :: Double -> Double -> VecT -> ColourT
renderSupersample sx sy or 	
	= (foldl1 colAdd (superSampleX (-(superPasses / 2)))) `colScalar` (1/((superPasses+1)^2))
	where
	superSampleX passX 
		| passX < (superPasses / 2) = superSampleY (-(superPasses / 2)) ++ superSampleX (passX+1)
		| otherwise = superSampleY (-(superPasses / 2))
		where
		superSampleY passY 
			| passY < (superPasses / 2) = col : (superSampleY (passY+1))
			| otherwise = [col]
			where
			superDir = normalize((VecT (sx + dx) (sy + dy) 0.0) `vecSub` or)
			dx = getDX * passX / superPasses
			dy = getDY * passY / superPasses
			superRay = RayT or superDir
			(obj,col,dst) = raytrace superRay 0 1.0
				
{---- HDR renderer -----}
-- This is alternative version of the renderer which uses Data.Map to store pixels with their coordinates.
-- It then applies a HDR+bloom filter to the screen (although it's not fully HDR, the end effect is nice ;))
-- First, it downsamples the image by factor of 1/4.Then, it suppresses all LDR values in the map. Next, 
-- a horizontal Gaussian-like filter is applied to bleed the pixels by defined offset. Eventually, vertical 
-- Gaussian-like filter is also applied and the two images are merged. The reason I didn't use it by default 
-- is its enormous memory usage. It seems the GC can't handle this code well, and I don't want to spend another 
-- week trying to rewrite it using arrays. This is probably a great example of how you should not write progams
-- with Data.Map ;)

-- Bleed all the pixels in the map and blend them together
-- Input: map of pixels
-- Output: map of blended bled pixels
bleedPixels :: Int -> (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)
bleedPixels offset pixels 
	= foldl1 (Map.unionWith colAdd) (parMap rnf (\f -> bleedColumn offset (Map.filterWithKey (\(_,x) _ -> x == f ) allRows)) [0..hdrWidth])
	where
	allRows = foldl1 (Map.unionWith colAdd) (parMap rnf (\f -> bleedRow offset (Map.filterWithKey (\(y,_) _ -> y == f) pixels)) [0..hdrHeight])
	
-- Bleed a column of pixels
-- Input: bleed intensity, map of one column
-- Output: bled map of one column
bleedColumn :: Int -> (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)										
bleedColumn offset pixelsInColumn 
	= if bleedList /= [] then bleedFromList bleedList Map.empty else Map.empty
	where
	offWord = fromIntegral offset
	bleedList = Map.toList pixelsInColumn 
	bleedFromList (((y,x),col):xs) mm 	
		| xs /= [] =  bleedFromList xs (aux (-offset) mm)
		| otherwise = aux (-offset) mm
		where
		aux off m 
			| off == 0 = aux (off+1) (Map.insertWith colAdd (y,x) (col `colScalar` bloomInt ) m)
			| off < offWord && (y') >= 0 && (y') <= hdrHeight = aux (off+1) (Map.insertWith colAdd ((y'),x) (col `colScalar` (bloomInt / (fromIntegral (abs off)))) m)
			| off == offWord && (y') >= 0 && (y') <= hdrHeight = Map.insertWith colAdd ((y'),x) (col `colScalar` (bloomInt / (fromIntegral (abs off)))) m
			| otherwise = Map.empty
			where
			y' = y + fromIntegral off
			
-- Bleed a row of pixels. NOTE: remember the coordinates are reversed, ie. (y,x) - easier cast to list
-- Input: bleed intensity, map of one row
-- Output: bled map of one row
bleedRow :: Int -> (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)
bleedRow offset pixelsInLine 		
	= if bleedList /= [] then bleedFromList bleedList Map.empty else Map.empty
	where
	offWord = fromIntegral offset
	bleedList = Map.toList pixelsInLine
	bleedFromList (((y,x),col):xs) mm 	
		| xs /= [] =  bleedFromList xs (aux (-offset) mm)
		| otherwise = aux (-offset) mm
		where
		aux off m 
			| off == 0 = aux (off+1) (Map.insertWith colAdd (y,x) (col `colScalar` bloomInt ) m)
			| off < offWord && (x') >= 0 && (x') <= hdrWidth = aux (off+1) (Map.insertWith colAdd (y,(x')) (col `colScalar` (bloomInt / (fromIntegral (abs off)))) m)
			| off == offWord && (x') >= 0 && (x') <= hdrWidth = Map.insertWith colAdd (y,(x')) (col `colScalar` (bloomInt / (fromIntegral (abs off)))) m
			| otherwise = Map.empty
			where
			x' = x + fromIntegral off
																																	
				
-- Remove low dynamic range colours from the list (sets them to 0)													
suppressLDR :: (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)
suppressLDR col = Map.filter (\f -> r f > 1 || g f > 1 || b f > 1) col
	
-- Downsample the screen by factor of 4. This speeds up HDR rendering and adds up some blur
downsampleScreen :: (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)
downsampleScreen clr = Map.mapKeysWith (\f y -> f `colAdd` y `colScalar` 0.25)	(\(y,x) -> (y `div` 2, x `div` 2)) (clr)

-- Transforms the scene into high dynamic range image
-- Then applies bloom and finally blends HDR and original images					
renderHDR :: (Map.Map (Word16,Word16) ColourT) -> (Map.Map (Word16,Word16) ColourT)
renderHDR clr 	= resample
	where
	resample = Map.mapWithKey (\(y,x) col -> col `colAdd` (Map.findWithDefault black (y `div` 2, x `div` 2) bled)) clr
	bled = bleedPixels bloomPasses (downsampleScreen (suppressLDR clr))

-- Code below is for supersampling anti-aliasing. 
-- Runs raytraces for every line, starting with input integer	
-- This runs every row computation in parallel, which is then dynamically
-- allocated to idle cores by the compiler. Note this is anti-aliased 
-- by multisampling, however only in row, as parallelism doesnt allow
-- passing last hit object through rows
-- Output: map with (y,x) coordinates paired with pixel colour
hdr_render :: Map.Map (Word16,Word16) ColourT
hdr_render = foldl1 (Map.unionWith colAdd) (parMap rnf (\f -> hdr_renderRow 0 f origin) [0..(screenHeight-1)])

-- Returns a list of colours of all pixels in the line X		
-- Input: X line number, view origin, Y line number,  last hit primitive
-- Output: list of map singletons with ((y,x) colour)			
hdr_renderRow :: Word16 -> Word16 -> VecT -> Map.Map (Word16,Word16) ColourT
hdr_renderRow x y o 
	= renderRow' x y o background Map.empty
	where
	renderRow' lineX lineY or lstObj m 
		| lineX < (screenWidth-1) = renderRow' (lineX+1) lineY or obj (Map.insertWith colAdd (lineY,lineX) colour m) 
		| otherwise = Map.insertWith colAdd (lineY,lineX) colour m
		where
		sX = renderBoundX + (fromIntegral lineX * getDX)
		sY = renderBoundY + (fromIntegral lineY * getDY)
		rayDir = normalize ((VecT sX sY 0.0) `vecSub` or)
		ray = (RayT or rayDir)
		(obj,clr,dst) = raytrace ray 0 1.0
		-- Do supersampling only when new primitive is hit
		colour = if (name lstObj) /= (name obj) then renderSupersample sX sY or else clr
												
-- Converts list of list of vectors into strings, to allow it to be saved to file														
convertToPPM :: [ColourT] -> String
convertToPPM [] = []
convertToPPM (y:ys) 	
	= convertToPPM' y ++ convertToPPM ys
	where
	convertToPPM' col	 = (printf "%u %u %u\n" ( r' :: Word16) (g' :: Word16) (b' :: Word16))
		where
		r' = clamp 0 255 (round(255 * r col))
		g' = clamp 0 255 (round(255 * g col))
		b' = clamp 0 255 (round(255 * b col))
												
-- Truncates values beyond the range
clamp :: Word16 -> Word16 -> Word16 -> Word16
clamp min max val 	
	| val > max = max
	| val < min = min
	| otherwise = val

-- Opens file and returns its handle
getAndOpenFile :: String -> IOMode -> IO Handle
getAndOpenFile name mode = catch (openFile name mode)
             (\_ -> do putStrLn ("Cannot open "++ name ++ "\n")
                       getAndOpenFile name mode)

-- Run tracer and save its result to file					   
main = do 
		args <- getArgs
		putStrLn "Starting raytracing..."
		startTime <- getCurrentTime
		writeFile "raytrace.ppm" ("P3\n" ++ show screenWidth ++ "\n" ++ show screenHeight ++ "\n255\n" ++ convertToPPM (if filter (List.isInfixOf "hdr") args /= [] then Map.elems (renderHDR hdr_render) else render))
		endTime <- getCurrentTime
		putStrLn ("Render time: " ++ show (diffUTCTime endTime startTime))
		getChar