module Solid where
import Vec
import Clr
import Data.List hiding (group)
--COMMON DATATYPES AND FUNCTIONS--
data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show
data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple
--union of two bounding boxes
bbjoin :: Bbox -> Bbox -> Bbox
bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =
(Bbox (vmin p1a p1b) (vmax p2a p2b))
--overlap of two bounding boxes
bboverlap :: Bbox -> Bbox -> Bbox
bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =
(Bbox (vmax p1a p1b) (vmin p2a p2b))
--split a bounding box into two
bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)
bbsplit (Bbox p1 p2) axis offset =
if (offset < (va p1 axis)) || (offset > (va p2 axis))
then error "degenerate bounding box split"
else ((Bbox p1 (vset p2 axis offset)),
(Bbox (vset p1 axis offset) p2))
-- generate a bounding box from a list of points
bbpts :: [Vec] -> Bbox
bbpts [] = empty_bbox
bbpts ((Vec x y z):[]) =
Bbox (Vec (x-delta) (y-delta) (z-delta))
(Vec (x+delta) (y+delta) (z+delta))
bbpts ((Vec x y z):pts) =
let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts
minx = fmin (x-delta) p1x
miny = fmin (y-delta) p1y
minz = fmin (z-delta) p1z
maxx = fmax (x+delta) p2x
maxy = fmax (y+delta) p2y
maxz = fmax (z+delta) p2z in
Bbox (Vec minx miny minz) (Vec maxx maxy maxz)
-- surface area, volume
bbsa :: Bbox -> Flt
bbsa (Bbox p1 p2) =
let Vec dx dy dz = vsub p2 p1
in dx*dy + dx*dz + dy*dz
bbvol :: Bbox -> Flt
bbvol (Bbox p1 p2) =
let (Vec dx dy dz) = vsub p2 p1
in dx*dy*dz
empty_bbox =
Bbox (Vec infinity infinity infinity)
(Vec (-infinity) (-infinity) (-infinity))
everything_bbox =
Bbox (Vec (-infinity) (-infinity) (-infinity))
(Vec infinity infinity infinity)
bbclip :: Ray -> Bbox -> Interval
bbclip (Ray (Vec ox oy oz) (Vec dx dy dz))
(Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =
let dxrcp = 1/dx
dyrcp = 1/dy
dzrcp = 1/dz
Interval inx outx = if dx > 0
then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)
else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)
Interval iny outy = if dy > 0
then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)
else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)
Interval inz outz = if dz > 0
then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)
else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)
in
Interval (fmax3 inx iny inz) (fmin3 outx outy outz)
data Rayint = RayHit {
depth :: !Flt,
pos :: !Vec,
norm :: !Vec,
texture :: !Texture
} | RayMiss deriving Show
nearest :: Rayint -> Rayint -> Rayint
nearest a RayMiss = a
nearest RayMiss b = b
nearest (RayHit da pa na ta) (RayHit db pb nb tb) =
if da < db
then RayHit da pa na ta
else RayHit db pb nb tb
furthest :: Rayint -> Rayint -> Rayint
furthest a RayMiss = RayMiss
furthest RayMiss b = RayMiss
furthest (RayHit da pa na ta) (RayHit db pb nb tb) =
if da > db
then RayHit da pa na ta
else RayHit db pb nb tb
hit :: Rayint -> Bool
hit (RayHit _ _ _ _) = True
hit RayMiss = False
dist :: Rayint -> Flt
dist RayMiss = infinity
dist (RayHit d _ _ _) = d
--LIGHTS--
data Light = Light {litpos :: !Vec,
litcol :: !Color} deriving Show
--MATERIALS--
data Material = Material {clr :: Color,
reflect, refract, ior,
kd, shine :: !Flt} deriving Show
type Texture = Rayint -> Material
-- this is sort of a no-op; we don't have a
-- good way to show an arbitrary function
showTexture :: Texture -> String
showTexture t = show $ t RayMiss
instance Show Texture where
show = showTexture
m_white = (Material c_white 0 0 0 1 2)
t_white ri = m_white
t_uniform :: Material -> Texture
t_uniform m = \x -> m
interp :: Flt -> Flt -> Flt -> Flt
interp scale a b =
scale*a + (1-scale)*b
--not really correct, but we'll go with it for now
m_interp :: Material -> Material -> Flt -> Material
m_interp m1 m2 scale =
let (Material m1c m1refl m1refr m1ior m1kd m1shine) = m1
(Material m2c m2refl m2refr m2ior m2kd m2shine) = m2
intp = interp scale
c = cadd (cscale m1c scale) (cscale m2c (1-scale))
refl = intp m1refl m2refl
refr = intp m1refr m2refr
ior = intp m1ior m2ior
kd = intp m1kd m2kd
shine = intp m1shine m2shine
in (Material c refl refr ior kd shine)
--SOLID TYPES--
data Solid = Sphere {center :: !Vec,
radius, invradius :: !Flt}
| Triangle {v1, v2, v3 :: Vec}
| TriangleNorm {v1, v2, v3, n1, n2, n3 :: Vec}
| Disc !Vec !Vec !Flt -- position, normal, r*r
| Cylinder !Flt !Flt !Flt -- radius height1 height2
| Cone !Flt !Flt !Flt !Flt -- r clip1 clip2 height
| Plane Vec Flt -- normal, offset from origin
| Box !Bbox
| Group ![Solid]
| Intersection ![Solid]
| Bevel !Solid !Flt
| Bound Solid Solid
| Difference !Solid !Solid
| Bih {bihbb :: !Bbox, bihroot :: !BihNode}
| Instance !Solid !Xfm
| Tex !Solid Texture
| Portal Solid Solid -- if we hit a, intersect with b
| Photon Vec Vec Color Flt -- pos, incident ray, color, radius
| Nothing deriving Show -- conflicts with
-- Nothing :: Maybe a from prelude
data BihNode = BihLeaf !Solid
| BihBranch {lmax :: !Flt, rmin :: !Flt, ax :: !Int,
l :: BihNode, r :: BihNode} deriving Show
--CONSTRUCTORS--
sphere :: Vec -> Flt -> Solid
sphere c r =
Sphere c r (1.0/r)
triangle :: Vec -> Vec -> Vec -> Solid
triangle v1 v2 v3 =
Triangle v1 v2 v3
--simple tesselation
triangles :: [Vec] -> [Solid]
triangles (v1:vs) =
zipWith (\v2 v3 -> triangle v1 v2 v3) vs (tail vs)
trianglenorm v1 v2 v3 n1 n2 n3 =
-- Triangle v1 v2 v3
TriangleNorm v1 v2 v3 n1 n2 n3
trianglesnorms :: [(Vec,Vec)] -> [Solid]
trianglesnorms (vn1:vns) =
zipWith (\vn2 vn3 -> trianglenorm (fst vn1) (fst vn2) (fst vn3)
(snd vn1) (snd vn2) (snd vn3))
vns (tail vns)
box :: Vec -> Vec -> Solid
box p1 p2 =
Box (Bbox p1 p2)
disc :: Vec -> Vec -> Flt -> Solid
disc pos norm r =
Disc pos norm (r*r)
cylinder_z :: Flt -> Flt -> Flt -> Solid
cylinder_z r h1 h2 = Cylinder r h1 h2
cone_z :: Flt -> Flt -> Flt -> Flt -> Solid
cone_z r h1 h2 height = Cone r h1 h2 height
-- construct a general cylinder from p1 to p2 with radius r
cylinder :: Vec -> Vec -> Flt -> Solid
cylinder p1 p2 r =
let axis = vsub p2 p1
len = vlen axis
ax1 = vscale axis (1/len)
(ax2,ax3) = orth ax1
in Instance (cylinder_z r 0 len)
(compose [ (xyz_to_uvw ax2 ax3 ax1),
(translate p1) ])
-- similar for cone
cone :: Vec -> Flt -> Vec -> Flt -> Solid
cone p1 r1 p2 r2 =
if r1 < r2
then cone p2 r2 p1 r1
else if r1-r2 < delta
then cylinder p1 p2 r2
else
let axis = vsub p2 p1
len = vlen axis
ax1 = vscale axis (1/len)
(ax2,ax3) = orth ax1
height = (r1*len)/(r1-r2) -- distance to end point
in
Instance (cone_z r1 0 len height)
(compose [ (xyz_to_uvw ax2 ax3 ax1),
(translate p1) ])
plane :: Vec -> Vec -> Solid
plane orig norm_ = Plane norm d
where norm = vnorm norm_
d = vdot orig norm
-- flatten tree of groups into a single group
flatten_group :: [Solid] -> [Solid]
flatten_group ((Group slds):xs) =
(flatten_group slds) ++ xs
flatten_group ((Solid.Nothing):xs) = xs
flatten_group x = x
group :: [Solid] -> Solid
group [] = Solid.Nothing
group (sld:[]) = sld
group slds =
Group (flatten_group slds)
transform :: Solid -> [Xfm] -> Solid
transform (Instance s xfm2) xfm1 =
transform s [compose ([xfm2] ++ xfm1)]
transform s xfm =
Instance s (compose xfm)
-- push all the transforms out to the leaves
-- and throw away pre-existing bounding volumes
-- so we can run the bih constructor on the
-- resulting group
flatten_transform :: Solid -> [Solid]
flatten_transform (Group slds) =
flatten_group $ concat (map flatten_transform slds)
flatten_transform (Instance s xfm) =
case s of
Group slds -> flatten_transform $ group (map (\x -> transform x [xfm]) slds)
Bound sa sb -> flatten_transform (transform sb [xfm])
Instance sa xfm2 -> flatten_transform (transform s [xfm])
_ -> [transform s [xfm]]
flatten_transform (Bound sa sb) = flatten_transform sb
-- bih construction
build_leaf objs =
BihLeaf (group (map snd objs))
max_bih_sa = 0.3 :: Flt
build_rec :: [(Bbox,Solid)] -> Bbox -> Bbox -> Int -> BihNode
build_rec objs nodebox splitbox depth =
-- if (null objs) || (null $ tail objs) ||
-- (null $ tail $ tail objs)
if length objs < 2
then build_leaf objs
else
let (Bbox nodeboxp1 nodeboxp2) = nodebox
(Bbox splitboxp1 splitboxp2) = splitbox
axis = vmaxaxis (vsub splitboxp2 splitboxp1)
bbmin = va splitboxp1 axis
bbmax = va splitboxp2 axis
candidate = (bbmin + bbmax) * 0.5
in
if candidate > (va nodeboxp2 axis) then
build_rec objs nodebox
(Bbox splitboxp1 (vset splitboxp2 axis candidate))
depth
else
if candidate < (va nodeboxp1 axis) then
build_rec objs nodebox (
Bbox (vset splitboxp1 axis candidate) splitboxp2)
depth
else
-- not sure if this is a big win
let nbsa = bbsa nodebox
(big,small) = partition (\ (bb,_) ->
(bbsa bb) > (nbsa * max_bih_sa)) objs
in
if (not $ null big) && ((length big) < ((length small)*2))
then (BihBranch (va nodeboxp2 0) (va nodeboxp1 0) 0
(build_rec big nodebox splitbox (depth+1))
(build_rec small nodebox splitbox (depth+1)) )
else
let (l,r) = partition (\((Bbox bbp1 bbp2),_)->
(((va bbp1 axis)+(va bbp2 axis))*0.5)
< candidate ) objs
lmax = foldl fmax (-infinity) (map (\((Bbox _ p2),_) -> va p2 axis) l)
rmin = foldl fmin infinity (map (\((Bbox p1 _),_) -> va p1 axis) r)
(lsplit,rsplit) = bbsplit splitbox axis candidate
lnb = (Bbox nodeboxp1 (vset nodeboxp2 axis lmax))
rnb = (Bbox (vset nodeboxp1 axis rmin) nodeboxp2)
in
-- stop if there's no progress being made
if ((null l) && (rmin <= bbmin)) ||
((null r) && (lmax >= bbmax))
then build_leaf objs
else
(BihBranch (lmax+delta) (rmin-delta) axis
(build_rec l lnb lsplit (depth+1))
(build_rec r rnb rsplit (depth+1)) )
bih :: [Solid] -> Solid
bih [] = Solid.Nothing
-- bih (sld:[]) = sld -- sometimes we'd like to be able to use a
-- single object bih just for its aabb
bih slds =
let objs = map (\x -> ((bound x),x)) (flatten_group slds)
bb = foldl bbjoin empty_bbox (map (\(b,_)->b) objs)
root = build_rec objs bb bb 0
(Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bb
in
if p1x == (-infinity) || p1y == (-infinity) || p1z == (-infinity) ||
p2x == infinity || p2y == infinity || p2z == infinity
then
error $ "bih: infinite bounding box " ++ (show objs)
else
(Bih bb root)
--INTERSECTION TESTS--
rayint :: Solid -> Ray -> Flt -> Texture -> Rayint
--Basic Primitives--
--Triangle--
-- adaptation of Moller and Trumbore from pbrt page 127
rayint (Triangle p1 p2 p3) (Ray o dir) dist tex =
let e1 = vsub p2 p1
e2 = vsub p3 p1
s1 = vcross dir e2
divisor = vdot s1 e1
in
if (divisor == 0)
then RayMiss
else
let invdivisor = 1.0 / divisor
d = vsub o p1
b1 = (vdot d s1) * invdivisor
in
if (b1 < 0) || (b1 > 1)
then RayMiss
else
let s2 = vcross d e1
b2 = (vdot dir s2) * invdivisor
in
if (b2 < 0) || (b1 + b2 > 1)
then RayMiss
else
let t = (vdot e2 s2) * invdivisor
in
if (t < 0) || (t > dist)
then RayMiss
else
RayHit t (vscaleadd o dir t) (vnorm (vcross e1 e2)) tex
rayint (TriangleNorm p1 p2 p3 n1 n2 n3) (Ray o dir) dist tex =
let e1 = vsub p2 p1
e2 = vsub p3 p1
s1 = vcross dir e2
divisor = vdot s1 e1
in
if (divisor == 0)
then RayMiss
else
let invdivisor = 1.0 / divisor
d = vsub o p1
b1 = (vdot d s1) * invdivisor
in
if (b1 < 0) || (b1 > 1)
then RayMiss
else
let s2 = vcross d e1
b2 = (vdot dir s2) * invdivisor
in
if (b2 < 0) || (b1 + b2 > 1)
then RayMiss
else
let t = (vdot e2 s2) * invdivisor
in
if (t < 0) || (t > dist)
then RayMiss
else
let n1scaled = (vscale n1 (1-(b1+b2)))
n2scaled = (vscale n2 b1)
n3scaled = (vscale n3 b2)
norm = vnorm $ vadd3 n1scaled n2scaled n3scaled
in RayHit t (vscaleadd o dir t) norm tex
--Sphere--
-- adapted from graphics gems volume 1
rayint (Sphere center r invr) (Ray e dir) dist t =
let eo = vsub center e
v = vdot eo dir
in
if (dist >= (v - r)) && (v > 0.0)
then
let vsqr = v*v
csqr = vdot eo eo
rsqr = r*r
disc = rsqr - (csqr - vsqr) in
if disc < 0.0 then
RayMiss
else
let d = sqrt disc
hitdist = if (v-d) > 0 then (v-d) else (v+d)
in if (hitdist < 0) || (hitdist > dist)
then RayMiss
else
let p = vscaleadd e dir hitdist
-- n = vscale (vsub p center) invr in
-- n = vsub (vscale p invr) (vscale center invr) in
n = vnorm (vsub p center)
in RayHit hitdist p n t
else
RayMiss
-- nice and simple
rayint (Disc point norm radius_sqr) r d t =
let (Ray orig dir) = r
dist = plane_int_dist r point norm
in if dist < 0 || dist > d
then RayMiss
else let pos = vscaleadd orig dir dist
offset = vsub pos point
in
if (vdot offset offset) > radius_sqr
then RayMiss
else RayHit dist pos norm t
-- cylinder aligned to z axis
-- no end caps
-- adapted from pbrt
rayint (Cylinder r h1 h2) (Ray orig dir) d t =
let Vec ox oy oz = orig
Vec dx dy dz = dir
a = dx*dx + dy*dy
b = 2*(dx*ox + dy*oy)
c = ox*ox + oy*oy - r*r
disc = b*b - 4*a*c
in if disc < 0
then RayMiss
else
let discsqrt = sqrt disc
q = if b < 0
then (b-discsqrt)*(-0.5)
else (b+discsqrt)*(-0.5)
t0' = q/a
t1' = c/q
t0 = fmin t0' t1'
t1 = fmax t0' t1'
in if t1 < 0 || t0 > d
then RayMiss
else let dist = if t0 < 0
then t1
else t0
in if dist < 0 || dist > d
then RayMiss
else let pos = vscaleadd orig dir dist
Vec posx posy posz = pos
in if posz > h1 && posz < h2
then RayHit dist pos (Vec (posx/r) (posy/r) 0) t
else if dz > 0 -- ray pointing up from bottom
then if oz < h1
then rayint (Disc (Vec 0 0 h1) nvz (r*r)) (Ray orig dir) d t
--then rayint_aadisc h1 r (Ray orig dir) d t
else RayMiss
else if oz > h2
then rayint (Disc (Vec 0 0 h2) vz (r*r)) (Ray orig dir) d t
--rayint_aadisc h2 r (Ray orig dir) d t -- todo: fix normal
else RayMiss
-- cone centered on z axiz, height of hp, clipped at h1 and h2
rayint (Cone r clip1 clip2 height) (Ray orig dir) d t =
let Vec ox oy oz = orig
Vec dx dy dz = dir
k' = (r/height)
k = k'*k'
a = dx*dx + dy*dy - k*dz*dz
b = 2*(dx*ox + dy*oy - k*dz*(oz-height))
c = ox*ox + oy*oy - k*(oz-height)*(oz-height)
disc = b*b - 4*a*c
in if disc < 0
then RayMiss
else
let discsqrt = sqrt disc
q = if b < 0
then (b-discsqrt)*(-0.5)
else (b+discsqrt)*(-0.5)
t0' = q/a
t1' = c/q
t0 = fmin t0' t1'
t1 = fmax t0' t1'
in if t1 < 0 || t0 > d
then RayMiss
else let dist = if t0 < 0
then t1
else t0
in if dist < 0 || dist > d
then RayMiss
else
let pos = vscaleadd orig dir dist
Vec posx posy posz = pos
in if posz > clip1 && posz < clip2
then let invhyp = 1 / (sqrt (height*height + r*r))
up = r * invhyp
out = height * invhyp
r_ = sqrt (posx*posx + posy*posy)
correction = (out)/(r_)
in RayHit dist pos (Vec (posx*correction) (posy*correction) up) t
else
if dz > 0 -- ray pointing up from bottom
then if oz < clip1
then rayint (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d t
else RayMiss
else if oz > clip2
then let r2 = r*(1-((clip2-clip1)/(height)))
in rayint (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d t
--rayint_aadisc clip2 r2 (Ray orig dir) d t
else RayMiss
-- then rayint_aadisc clip1 r (Ray orig dir) d t
-- else RayMiss -- rayint_aadisc clip2
-- (r*((clip2-clip1)/height))
-- (Ray orig dir) d t -- todo: fix normal
--Plane--
rayint (Plane norm offset) (Ray orig dir) d t =
let hit = -(((vdot norm orig)-offset) / (vdot norm dir))
in if hit < 0 || hit > d
then RayMiss
else RayHit hit (vscaleadd orig dir hit) norm t
--Box--
-- this could be optimized a bit more
rayint (Box (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z))) r d t =
let (Ray orig dir) = r
(Vec ox oy oz) = orig
(Vec dx dy dz) = dir
dxrcp = 1/dx
dyrcp = 1/dy
dzrcp = 1/dz
Interval inx outx = if dx > 0
then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)
else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)
Interval iny outy = if dy > 0
then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)
else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)
Interval inz outz = if dz > 0
then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)
else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)
lastin = (fmax3 inx iny inz)
firstout = (fmin3 outx outy outz)
in if lastin > firstout || firstout < 0 || lastin > d
then RayMiss
else
let n = if inx == lastin
then if dx > 0 then nvx else vx
else if iny == lastin
then if dy > 0 then nvy else vy
else if dz > 0 then nvz else vz
norm = if lastin > 0 then n else vinvert n
hitdepth = fmax 0 lastin
in
RayHit hitdepth (vscaleadd orig dir hitdepth) norm t
--Composite objects--
--Instance--
-- transforming the distance is a little awkward
-- the normal shouldn't have to be re-normalized, should it?
rayint (Instance sld xfm) (Ray orig dir) d t =
let newdir = invxfm_vec xfm dir
neworig = invxfm_point xfm orig
lenscale = vlen newdir
invlenscale = 1/lenscale
in
case (rayint_check sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale) t) of
RayMiss -> RayMiss
RayHit depth pos n tex -> RayHit (depth*invlenscale) (xfm_point xfm pos) (vnorm (invxfm_norm xfm n)) tex
--Group--
-- rayint (Group lst) r d t = foldl nearest RayMiss (map (\x -> rayint x r d t) lst)
rayint (Group xs) r d t =
let rig [] = RayMiss
rig (x:xs) = nearest (rayint_check x r d t) (rig xs)
in rig xs
--Difference--
-- csg of object a - object b --
rayint (Difference sa sb) r d t =
let dif = Difference sa sb
Ray orig dir = r
ria = rayint sa r d t
in
case ria of
RayMiss -> RayMiss
RayHit ad ap an at ->
if inside sb orig
then
case rayint sb r d t of
RayMiss -> RayMiss
RayHit bd bp bn bt ->
if bd < ad
then if inside sa bp
then RayHit bd bp (vinvert bn) bt
else rayint_advance dif r d t bd
else rayint_advance dif r d t bd
else
if inside sb ap
then rayint_advance dif r d t ad
else RayHit ad ap an at
--Intersection--
-- fixme: there's some numerical instability near edges
rayint (Intersection slds) r d t =
let (Ray orig dir) = r
in
if null slds || d < 0
then RayMiss
else
let s = head slds
in case tail slds of
[] -> rayint_check s r d t
ss -> if inside s orig
then case rayint s r d t of
RayMiss -> rayint (Intersection ss) r d t
RayHit sd sp sn st ->
case rayint_check (Intersection ss) r sd t of
RayMiss -> rayint_advance (Intersection slds)
r d t sd
hit -> hit
else case rayint s r d t of
RayMiss -> RayMiss
RayHit sd sp sn st ->
if inside (Intersection ss) sp
then RayHit sd sp sn st
else rayint_advance (Intersection slds)
r d t sd
rayint (Bound sa sb) r d t =
let (Ray orig _) = r
in if inside sa orig || shadow sa r d
then rayint sb r d t
else RayMiss
--Bih--
rayint (Bih bb root) r d t =
let Ray orig dir = r
dir_rcp = vrcp dir
Interval near far = bbclip r bb
traverse (BihLeaf s) near far = rayint_check s r (fmin d far) t
traverse (BihBranch lsplit rsplit axis l r) near far =
let dirr = va dir_rcp axis
o = va orig axis
dl = (lsplit - o) * dirr
dr = (rsplit - o) * dirr
in
if near > far
then RayMiss
else
-- this is ugly and verbose,
-- but it does what it needs to
if dirr > 0
then
(nearest
(if near < dl
then traverse l near (fmin dl far)
else RayMiss)
(if dr < far
then traverse r (fmax dr near) far
else RayMiss))
else
(nearest
(if near < dr
then traverse r near (fmin dr far)
else RayMiss)
(if dl < far
then traverse l (fmax dl near) far
else RayMiss)) --}
in
traverse root near far
-- anything that hits it gets teleported
-- one-way door to another world
-- this is dangerous and doesn't really work right
rayint (Portal port world) r d t =
case rayint port r d t of
RayMiss -> RayMiss
RayHit depth _ _ _ ->
rayint_advance world r d t depth
--Tex--
-- this is a little odd; rather than associate
-- a texture with each primitive, we use a
-- container object; everything inside has that
-- texture, unless it's overridden by a nested
-- texture
rayint (Tex s tex) r d t = rayint_check s r d tex
--Nothing--
rayint (Solid.Nothing) _ _ _ = RayMiss
-- default case: miss
-- rayint _ _ _ _ = RayMiss
-- various specialized ray intersections, used as helper functions
-- used by cylinder / cone code
-- broken, do not use
rayint_aadisc :: Flt -> Flt-> Ray -> Flt -> Texture -> Rayint
rayint_aadisc height radius r d t =
let Ray orig dir = r
-- Vec _ _ oz = orig
-- Vec _ _ dz = dir
-- dist = (height-oz)/dz
dist = plane_int_dist r (Vec 0 0 height) vx
in if dist <= 0 || dist >= d || isNaN dist
then RayMiss
else let pos = vscaleadd orig dir dist
(Vec px py _) = pos
in
if (px*px + py+py) > radius*radius
then RayMiss
else RayHit dist pos vz t
-- move ray forward, intersect, fix result
-- useful in csg
rayint_advance :: Solid -> Ray -> Flt -> Texture -> Flt -> Rayint
rayint_advance s r d t adv =
let a = adv+delta
in
case (rayint s (ray_move r a) (d-a) t) of
RayMiss -> RayMiss
RayHit depth pos norm tex -> RayHit (depth+a) pos norm tex
-- check results of a ray-intersection test
rayint_check s r d t =
let Ray orig dir = r in
case rayint s r d t of
RayMiss -> RayMiss
RayHit depth pos norm tex ->
if depth < 0
then error $ "rayint depth < 0 " ++ (show depth) ++ " " ++ (show s)
else
if depth > d
then error $ "rayint depth (" ++ (show depth) ++ ") > d (" ++ (show d) ++ ") " ++ (show s)
else
if not $ veq pos (vscaleadd orig dir depth)
then error $ "rayint position doesn't match depth" ++
(show pos) ++ (show $ vscaleadd orig dir depth) ++ (show depth) ++ (show s)
else
if (vlen norm) < 1-delta
then error "normal too short"
else
if (vlen norm) > 1+delta
then error $ "normal too long " ++ (show norm) ++ " " ++ (show s)
else RayHit depth pos norm tex
--SHADOW--
shadow :: Solid -> Ray -> Flt -> Bool
--Sphere--
shadow (Sphere center r invr) (Ray e dir) dist =
let eo = vsub center e
v = vdot eo dir
in
if (dist >= (v - r)) && (v > 0.0)
then
let vsqr = v*v
csqr = vdot eo eo
rsqr = r*r
disc = rsqr - (csqr - vsqr) in
if disc < 0.0 then
False
else
let d = sqrt disc
hitdist = if (v-d) > 0 then (v-d) else (v+d)
in if (hitdist < 0) || (hitdist > dist)
then False
else True
else
False
shadow (Triangle p1 p2 p3) (Ray o dir) dist =
let e1 = vsub p2 p1
e2 = vsub p3 p1
s1 = vcross dir e2
divisor = vdot s1 e1
in
if (divisor == 0)
then False
else
let invdivisor = 1.0 / divisor
d = vsub o p1
b1 = (vdot d s1) * invdivisor
in
if (b1 < 0) || (b1 > 1)
then False
else
let s2 = vcross d e1
b2 = (vdot dir s2) * invdivisor
in
if (b2 < 0) || (b1 + b2 > 1)
then False
else
let t = (vdot e2 s2) * invdivisor
in
if (t < 0) || (t > dist)
then False
else True
shadow (TriangleNorm p1 p2 p3 n1 n2 n3) r d =
shadow (Triangle p1 p2 p3) r d
shadow (Disc point norm radius_sqr) r d =
let (Ray orig dir) = r
dist = plane_int_dist r point norm
in if dist < 0 || dist > d
then False
else let pos = vscaleadd orig dir dist
offset = vsub pos point
in
if (vdot offset offset) > radius_sqr
then False
else True
shadow (Box box) r d =
let Interval near far = bbclip r box
in
if (near > far) || far <= 0 || far > d
then False
else True
-- cone centered on z axiz, height of hp, clipped at h1 and h2
shadow (Cone r clip1 clip2 height) (Ray orig dir) d =
let Vec ox oy oz = orig
Vec dx dy dz = dir
k' = (r/height)
k = k'*k'
a = dx*dx + dy*dy - k*dz*dz
b = 2*(dx*ox + dy*oy - k*dz*(oz-height))
c = ox*ox + oy*oy - k*(oz-height)*(oz-height)
disc = b*b - 4*a*c
in if disc < 0
then False
else
let discsqrt = sqrt disc
q = if b < 0
then (b-discsqrt)*(-0.5)
else (b+discsqrt)*(-0.5)
t0' = q/a
t1' = c/q
t0 = fmin t0' t1'
t1 = fmax t0' t1'
in if t1 < 0 || t0 > d
then False
else let dist = if t0 < 0
then t1
else t0
in if dist < 0 || dist > d
then False
else
let pos = vscaleadd orig dir dist
Vec posx posy posz = pos
in if posz > clip1 && posz < clip2
then True
else
if dz > 0 -- ray pointing up from bottom
then if oz < clip1
then shadow (Disc (Vec 0 0 clip1) nvz (r*r)) (Ray orig dir) d
else False
else if oz > clip2
then let r2 = r*(1-((clip2-clip1)/(height)))
in shadow (Disc (Vec 0 0 clip2) vz (r2*r2)) (Ray orig dir) d
else False
shadow (Instance sld xfm) (Ray orig dir) d =
let newdir = invxfm_vec xfm dir
neworig = invxfm_point xfm orig
lenscale = vlen newdir
invlenscale = 1/lenscale
in
shadow sld (Ray neworig (vscale newdir invlenscale)) (d*lenscale)
shadow (Tex s t) r d = shadow s r d
shadow (Group xs) r d =
let sg [] = False
sg (x:xs) = (shadow x r d) || (sg xs)
in sg xs
shadow (Bound sa sb) r d =
let (Ray orig _ ) = r
in if inside sa orig || shadow sa r d
then shadow sb r d
else False
shadow (Bih bb root) r d =
let (Ray orig dir) = r
dir_rcp = vrcp dir
Interval near far = bbclip r bb
traverse (BihLeaf s) near far = shadow s r (fmin d far)
traverse (BihBranch lsplit rsplit axis l r) near far =
let dirr = va dir_rcp axis
o = va orig axis
dl = (lsplit - o) * dirr
dr = (rsplit - o) * dirr
in
if near > far
then False
else
if dirr > 0
then
((if near < dl
then traverse l near (fmin dl far)
else False)
||
(if dr < far
then traverse r (fmax dr near) far
else False))
else
((if near < dr
then traverse r near (fmin dr far)
else False)
||
(if dl < far
then traverse l (fmax dl near) far
else False))
in traverse root near far
-- default shadow test, in case an optimized test
-- isn't implemented yet, we use the regular
-- trace function; we have to be careful with
-- container objects, though; using the slow
-- test on the container means doing the slow
-- test against everything it contains
shadow s r d =
case (rayint s r d t_white) of
RayHit _ _ _ _ -> True
RayMiss -> False
--INSIDE--
inside :: Solid -> Vec -> Bool
inside (Sphere center r invr) pt =
let offset = vsub center pt
in (vdot offset offset) < r*r
inside (Group slds) pt =
foldl' (||) False (map (\x -> inside x pt) slds)
inside (Difference sa sb) pt =
(inside sa pt) && (not $ inside sb pt)
-- note: inside is True for an empty intersection.
-- this is actually the preferred semantics in
-- some cases, strange as it may seem.
inside (Intersection slds) pt =
foldl' (&&) True (map (\x -> inside x pt) slds)
inside (Instance s xfm) pt =
inside s (xfm_point xfm pt)
inside (Plane norm offset) pt =
let onplane = (vscale norm offset)
newvec = vsub onplane pt
in vdot newvec norm > 0
inside (Box (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2))) (Vec x y z) =
x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2
inside (Tex s t) pt = inside s pt
inside (Bih (Bbox (Vec x1 y1 z1) (Vec x2 y2 z2)) root) pt =
let (Vec x y z) = pt
traverse (BihLeaf s) = inside s pt
traverse (BihBranch lsplit rsplit axis l r) =
let o = va pt axis
in (if o < lsplit
then (traverse l)
else False)
||
(if o > rsplit
then (traverse r)
else False)
in
(x > x1) && (x < x2) && (y > y1) && (y < y2) && (z > z1) && (z < z2) && (traverse root)
inside (Cylinder r h1 h2) (Vec x y z) =
z > h1 && z < h2 && x*x + y*y < r*r
inside (Cone rbase h1 h2 height) (Vec x y z) =
let r = rbase*(1-(((z-h1)/height)))
in z > h1 && z < h2 && x*x + y*y < r*r
inside (Bound sa sb) pt = inside sa pt && inside sb pt
inside _ _ = False
-- return distance to surface, positive if inside, negative if outside
-- this isn't used for anything in particular yet
power :: Solid -> Vec -> Flt
power (Sphere center r invr) pt =
let offset = vsub center pt
in r - (vlen offset)
power (Group slds) pt =
foldl' (max) (-infinity) (map (\x -> power x pt) slds)
-- not accurate
power (Instance s xfm) pt =
power s (xfm_point xfm pt)
power (Plane norm offset) pt =
let onplane = (vscale norm offset)
newvec = vsub onplane pt
in vdot newvec norm
--BOUND--
bound :: Solid -> Bbox
bound (Sphere center r invr) =
let offset = (vec r r r) in
(Bbox (vsub center offset) (vadd center offset))
bound (Triangle (Vec v1x v1y v1z)
(Vec v2x v2y v2z)
(Vec v3x v3y v3z)) =
Bbox
(Vec ((fmin (fmin v1x v2x) v3x) - delta)
((fmin (fmin v1y v2y) v3y) - delta)
((fmin (fmin v1z v2z) v3z) - delta) )
(Vec ((fmax (fmax v1x v2x) v3x) + delta)
((fmax (fmax v1y v2y) v3y) + delta)
((fmax (fmax v1z v2z) v3z) + delta) )
bound (TriangleNorm v1 v2 v3 n1 n2 n3) =
bound (Triangle v1 v2 v3)
--dangerous to use inside a bih
bound (Plane norm offset) = everything_bbox
bound (Box box) = box
-- this could be a tighter fit
bound (Disc pos norm rsqr) =
bound (sphere pos (sqrt rsqr))
bound (Cylinder r h1 h2) =
Bbox (Vec (-r) (-r) h1) (Vec r r h2)
bound (Cone r h1 h2 height) =
Bbox (Vec (-r) (-r) h1) (Vec r r h2)
bound (Group slds) =
foldl' bbjoin empty_bbox (map bound slds)
bound (Difference sa sb) = bound sa
bound (Intersection slds) =
if null slds
then empty_bbox
else foldl' bboverlap everything_bbox (map bound slds)
-- the reason the following doesn't work:
-- bound (Bound sa sb) = bboverlap (bound sa) (bound sb)
-- is that the ray has to hit solid a, and it might not
-- if the bound is smaller than sa.
-- update: this is probably no longer true, as I added an "inside"
-- test to the bounds check as well
bound (Bound sa sb) = bound sa
-- not optimal, but it does the job
bound (Instance sld xfm) =
let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bound sld
pxfm = xfm_point xfm
in
bbpts [(pxfm (Vec x y z)) | x <- [p1x,p2x],
y <- [p1y,p2y],
z <- [p1z,p2z]]
bound (Bih bb root) = bb
bound (Portal port world) = bound port
bound (Tex s t) = bound s
bound (Photon p indir clr r) = bound (sphere p r)
bound Solid.Nothing = empty_bbox
-- no default case: we want an exception
-- if there is no match
--CAMERA--
data Camera = Camera {campos, fwd, up, right :: !Vec}
deriving Show
default_cam = (Camera (vec 0.0 0.0 (-3.0))
(vec 0.0 0.0 1.0)
(vec 0.0 1.0 0.0)
(vec 1.0 0.0 0.0) )
camera :: Vec -> Vec -> Vec -> Flt -> Camera
camera pos at up angle =
let fwd = vnorm $ vsub at pos
right = vnorm $ vcross up fwd
up_ = vnorm $ vcross fwd right
cam_scale = tan ((pi/180)*(angle/2))
in
Camera pos fwd
(vscale up_ cam_scale)
(vscale right cam_scale)
--SCENE--
data Scene = Scene {sld :: !Solid,
lights :: ![Light],
cam :: !Camera,
dtex :: !Texture,
bground :: !Color} deriving Show
default_scene = (Scene (sphere (vec 0.0 0.0 0.0) 1.0)
[] default_cam t_white c_white)