雑記帳
Haskell でレイトレーシングのチュートリアルを追いかける その14 - カメラ
進捗
カメラの実装
今回は、
- カメラの位置、視野、向きの設定
- 焦点距離とピンボケの表現
といったカメラに関する内容を一挙に進めた。
中空ガラス描画テクニックをトーラスにも適用してみた
2次元トーラスの方程式は2次元球面の方程式の場合と同様、半径の符号は交点を求める過程には全く影響を与えないが、コード内の
return $ set_face_normal HitRecord {
p = c1 + x,
normal = zero,
t = k,
front_face = False,
mat = mat_Torus obj
} r ((x - (r1 *^ (normalize $ x - (n `dot` x) *^ n))) ^/ r2)
という部分を見ての通り、衝突面の法線方向を求める際には影響してくるといった感じになっている。
この場合、小半径 r2 の符号を負にすることで、チュートリアルで説明されているような中空ガラス (hollow glass) の描画テクを応用できる状況を作れるので、今回それもやってみた。(想像以上に綺麗に描画された！)
r2 の符号を負にすることで、チュートリアルで説明されているような中空ガラス (hollow glass) の描画テクを応用できる状況を作れるので、今回それもやってみた。(想像以上に綺麗に描画された！)前回からの変更点
y-軸の向きの変更
カメラの設定を行うにあたって、y-軸の方向をこれまで想定していた方向に対して真逆の方向に変更した。
理由としては、これまで用いていた基底に関するベクトルの成分表示だと、クロス積を取る操作について適切に振舞わなかったため。
一応クロス積について補足しておくと、クロス積の意味合いとしては
- 2つの空間ベクトル \(\boldsymbol{v},\boldsymbol{w}\) が与えられたとき、「始点を共有させた形で得られるそれら2つのベクトルを2辺に取るような空間内の平行四辺形」を考えることができるが、この時その平行四辺形の面積を大きさに持ち、\(\boldsymbol{v}\) から \(\boldsymbol{w}\) へと (平行四辺形内を通りながら) 向かう方向に右ねじを回した際に、その右ねじが進む向きを方向として持つようなベクトル \(\boldsymbol{v}\times \boldsymbol{w}\)
といった感じになる。
ここで問題になるのが、通常「クロス積の計算式」として与えられる式を用いて計算される空間ベクトルというのは、必ずしも上で述べた意味を持つベクトルにならないという点であり、実際、使用している成分表示を与えている座標系が右手座標系 (right-handed coordinate system) と呼ばれるものでないと、想定している方向を向いてくれない。
別段、左手座標系用に別途、前述した意味を持つベクトルを算出する関数を作ってしまうという手もあるのだけど、それだとややこしくなってしまうかもしれないということで、この判断に至った。
疑似乱数の発生させ方の変更
前回・前々回は各ピクセルの色を求める際に使用する疑似乱数の幾つかを先にリストとして作り、そのリストから値を引っぱってくるという方法を用いていたが、今回のピンボケ実装のために get_ray 関数の中での疑似乱数の使用が求められたため、その辺のやり方を変更した。
get_ray 関数の中での疑似乱数の使用が求められたため、その辺のやり方を変更した。余談
なんで乱数ごときにそんなに大騒ぎする必要があるのかを疑問に思う人は、一度 Haskell で乱数を使ってみよう。参照透過性を保障する Haskell で疑似乱数を用いた (IO とは無関係な) 関数を書きたいとなったとき、疑似乱数発生に使用する状態データを明示的に持ち歩く必要があるので、なかなか骨が折れる。
コードの実行結果
ソースコード
{-# LANGUAGE TypeOperators #-}
module Main where
import Data.Complex
import Control.Monad.Fix
import Control.Lens
import System.Random
import Linear.Vector
import Linear.Metric
import Linear.V3
-- https://raytracing.github.io/books/RayTracingInOneWeekend.html
main :: IO ()
main = do
let
-- Image
aspect_ratio = 16.0 / 9
image_width = 512
image_height = round $ fromIntegral image_width / aspect_ratio
samples_per_pixel = 100
max_depth = 100
-- World
material_ground = make_shared $ MAT_Lambertian {albedo_Lamb = V3 0.8 0.8 0.0}
material_red = make_shared $ MAT_Lambertian {albedo_Lamb = V3 0.7 0.3 0.3}
material_blue = make_shared $ MAT_Lambertian {albedo_Lamb = V3 0.3 0.5 0.7}
material_silver = make_shared $ MAT_Metal {albedo_Metal = V3 0.7 0.7 0.7, fuzz = 0.0}
material_gold = make_shared $ MAT_Metal {albedo_Metal = V3 0.8 0.6 0.2, fuzz = 0.0}
material_glass = make_shared $ MAT_Dielectric {ref_idx = 1.7}
material_water = make_shared $ MAT_Dielectric {ref_idx = 1.33}
world = []
`add` RT_Torus{
centerOfTorus = V3 (-1) 1 0,
majorRadius = 0.3,
minorRadius = 0.2,
orientationOfTorus = normalize $ V3 0.5 0.3 1.0,
mat_Torus = material_blue
}
`add` RT_Torus{
centerOfTorus = V3 0 1 0,
majorRadius = 0.35,
minorRadius = 0.15,
orientationOfTorus = normalize $ V3 0.4 (-1.5) 1.0,
mat_Torus = material_water
}
`add` RT_Torus{
centerOfTorus = V3 0 1 0,
majorRadius = 0.35,
minorRadius = -0.1,
orientationOfTorus = normalize $ V3 0.4 (-1.5) 1.0,
mat_Torus = material_water
}
`add` RT_Torus{
centerOfTorus = V3 1 1 0,
majorRadius = 0.3,
minorRadius = 0.2,
orientationOfTorus = normalize $ V3 (-0.7) (-0.7) 0.9,
mat_Torus = material_silver
}
`add` RT_Sphere{center = V3 0 1 (-100000.5), radius = 100000, mat_Sphere = material_ground}
-- Camera
lookfrom = V3 2.7 (-1.7) 3
lookat = V3 0 1 0
dist_to_focus = norm $ lookat - lookfrom
aperture = 0.1
camera = cam(lookfrom, lookat, V3 0 0 1, 20, aspect_ratio, aperture, dist_to_focus)
img_data = "P3\n" ++ show image_width ++ " " ++ show image_height ++ "\n255\n"
putStr $ img_data
foldr (>>) (return ()) $ do
let
indices = [image_height - 1, image_height - 2 .. 0] `prod` [0 .. image_width - 1]
seeds = (randomRs (0, 536870912) (mkStdGen 21) :: [Int])
((j,i), seed) <- zip indices seeds
return $ do
let
pixcel_color = getColor(0, 0, mkStdGen seed)
getColor = fix $ \rec (s, pixcel_color', gen_current) -> do
let
(randNum1, gen1) = random gen_current
(randNum2, gen2) = random gen1
u = (fromIntegral i + randNum1) / (fromIntegral image_width - 1.0)
v = (fromIntegral j + randNum2) / (fromIntegral image_height - 1.0)
(r,gen3) = get_ray camera (u, v) gen2
(d,gen4) = ray_color r world max_depth gen3
if s < samples_per_pixel then
rec(succ s, pixcel_color' + d, gen4)
else
pixcel_color'
write_color pixcel_color samples_per_pixel
data Ray = Ray {
orig :: V3 Double,
dir :: V3 Double
} deriving (Show)
at' :: Ray -> Double -> V3 Double
at' r t = (orig r) + t *^ (dir r)
data Camera = Camera {
viewport_height :: Double,
viewport_width :: Double,
focal_length :: Double,
origin :: V3 Double,
horizontal :: V3 Double,
vertical :: V3 Double,
lower_left_corner :: V3 Double,
axis_u :: V3 Double,
axis_v :: V3 Double,
axis_w :: V3 Double,
lens_radius :: Double
} deriving (Show)
cam (lookfrom,lookat,vup,vfov,aspect_ratio,aperture,focus_dist) =
Camera {
viewport_height = viewport_height',
viewport_width = viewport_width',
origin = origin',
horizontal = horizontal',
vertical = vertical',
lower_left_corner = lower_left_corner',
axis_u = u,
axis_v = v,
axis_w = w,
lens_radius = lens_radius'
}
where
theta = deg2rad(vfov)
h = tan(theta/2)
viewport_height' = 2.0 * h
viewport_width' = aspect_ratio * viewport_height'
w = normalize $ lookat - lookfrom
u = normalize $ w `cross` vup
v = u `cross` w
origin' = lookfrom
horizontal' = (focus_dist * viewport_width') *^ u
vertical' = (focus_dist * viewport_height') *^ v
lower_left_corner' =
origin' - horizontal' ^/2 - vertical' ^/2
+ focus_dist *^ w
lens_radius' = aperture / 2
get_ray :: Camera -> (Double, Double) -> StdGen -> (Ray, StdGen)
get_ray this (s, t) gen0 = (r,gen1)
where
(in_unit_disk, gen1) = random_in_unit_sphere gen0
rd = lens_radius this *^ in_unit_disk
offset = (rd ^._x) *^ axis_u this + (rd ^._y) *^ axis_v this
r =
Ray {
orig = origin this + offset,
dir = lower_left_corner this + s *^ horizontal this + t *^ vertical this - (origin this + offset)
}
type HittableData = (RT_Sphere + RT_Torus) + RT_Sphere
class Hittable a where
toSum :: a -> HittableData
hit :: a -> Ray -> Double -> Double -> Maybe HitRecord
instance (Hittable a, Hittable b) => Hittable (Either a b) where
toSum = coPair(toSum, toSum)
hit = coPair(hit, hit)
add :: Hittable a => [HittableData] -> a -> [HittableData]
add list obj = (toSum obj) : list
data HitRecord = HitRecord {
p :: V3 Double,
normal :: V3 Double,
mat :: MaterialData,
t :: Double,
front_face :: Bool
} deriving (Show)
set_face_normal :: HitRecord -> Ray -> V3 Double -> HitRecord
set_face_normal this r outward_normal = HitRecord {
p = p this,
normal = if dir r `dot` outward_normal < 0 then outward_normal else -outward_normal,
t = t this,
front_face = (dir r `dot` outward_normal < 0),
mat = mat this
}
hitSomething :: [HittableData] -> Ray -> Double -> Double -> Maybe HitRecord
hitSomething list r t_min t_max =
let
f (list', r', closest_so_far, currentRecord) =
case list' of
x:xs ->
let temp = hit x r' t_min closest_so_far in
case temp of
Just a ->
f $ (xs, r', t a, temp)
Nothing ->
f $ (xs, r', closest_so_far, currentRecord)
[] ->
currentRecord
in
f $ (list, r, t_max, Nothing)
----------------------
-- Hittable Objects --
----------------------
-- Sphere
data RT_Sphere = RT_Sphere {
center :: V3 Double,
radius :: Double,
mat_Sphere :: MaterialData
} deriving (Show)
instance Hittable RT_Sphere where
toSum = inj1 -: inj1
hit obj r t_min t_max =
let
p0 = orig r
c1 = center obj
r1 = radius obj
oc = p0 - c1
a = quadrance (dir r)
half_b = oc `dot` dir r
c = quadrance oc - (radius obj) ^ 2
discriminant = half_b ^ 2 - a*c in
if discriminant > 0 then
let
root = sqrt discriminant
f k =
case k of
x:xs ->
if t_min < x && x < t_max then
return $ set_face_normal HitRecord {
p = at' r x,
normal = zero,
t = x,
front_face = False,
mat = mat_Sphere obj
} r ((at' r x - c1) ^/ r1)
else
f $ xs
[] ->
Nothing
in
f $ [(-half_b - root) / a, (-half_b + root) / a]
else
Nothing
-- Torus
data RT_Torus = RT_Torus {
centerOfTorus :: V3 Double,
majorRadius :: Double,
minorRadius :: Double,
orientationOfTorus :: V3 Double,
mat_Torus :: MaterialData
} deriving (Show)
instance Hittable RT_Torus where
toSum = inj2 -: inj1
hit obj r t_min t_max =
let
p0 = orig r
a = dir r
a_norm = norm a
c1 = centerOfTorus obj
r1 = majorRadius obj
r2 = minorRadius obj
n = orientationOfTorus obj
s = getIntersection_forTorus (p0,a,c1,r1,r2,n)
oc = p0 - c1
a_sq = quadrance (dir r)
half_b = oc `dot` dir r
c = quadrance oc - (r1 + r2 + 0.01) ^ 2
discriminant = half_b ^ 2 - a_sq*c
in
if discriminant > 0 then
if null s then
Nothing
else
let
k = minimum s
x = at' r k - c1
in
if t_min < k && k < t_max then
return $ set_face_normal HitRecord {
p = c1 + x,
normal = zero,
t = k,
front_face = False,
mat = mat_Torus obj
} r ((x - (r1 *^ (normalize $ x - (n `dot` x) *^ n))) ^/ r2)
else
Nothing
else
Nothing
write_color :: V3 Double -> Int -> IO ()
write_color (V3 r g b) spp =
let
v' = V3 (sqrt $ r / fromIntegral spp) (sqrt $ g / fromIntegral spp) (sqrt $ b / fromIntegral spp)
f = show.floor.(256*).(clamp 0 0.999)
in
putStr $ f(v' ^._x) ++ " " ++ f(v' ^._y) ++ " " ++ f(v' ^._z) ++ "\n"
ray_color :: Ray -> [HittableData] -> Int -> StdGen -> (V3 Double, StdGen)
ray_color r objects depth gen =
if depth <= 0 then
(zero, gen)
else
let
record = hitSomething objects r 0.001 infinity
in
case record of
Just record' ->
let
(ret, gen1) = scatter (mat record') r record' gen
in
case ret of
Just (scattered, attenuation) ->
let
(col_t , gen_t) = ray_color scattered objects (pred depth) gen1
in
(attenuation * col_t , gen_t)
Nothing ->
(zero, gen1)
Nothing ->
let
unit_direction = normalize $ (dir r)
s = 0.5 * (unit_direction ^._z + 1.0)
in
(lerp s (V3 0.5 0.7 1.0) (V3 1.0 1.0 1.0), gen)
--------------------
-- Random numbers --
--------------------
random_in_unit_sphere :: StdGen -> (V3 Double, StdGen)
random_in_unit_sphere gen0 =
let
(rand1,gen1) = randomR (-1, 1) gen0 :: (Double, StdGen)
(rand2,gen2) = randomR (-1, 1) gen1 :: (Double, StdGen)
(rand3,gen3) = randomR (-1, 1) gen2 :: (Double, StdGen)
v = V3 rand1 rand2 rand3
in
if quadrance v >= 1 then
random_in_unit_sphere gen3
else
(v, gen3)
random_unit_vector :: StdGen -> (V3 Double, StdGen)
random_unit_vector gen0 =
let
(a, gen1) = randomR (0, 2*pi) gen0 :: (Double, StdGen)
(z, gen2) = randomR (-1, 1) gen1 :: (Double, StdGen)
r = sqrt $ 1 - z^2
in
(V3 (r*cos(a)) (r*sin(a)) z, gen2)
random_in_hemisphere :: V3 Double -> StdGen -> (V3 Double, StdGen)
random_in_hemisphere normal gen0 =
let
(in_unit_sphere, gen1) = random_in_unit_sphere gen0
in
if in_unit_sphere `dot` normal > 0 then
(in_unit_sphere, gen1)
else
(-in_unit_sphere, gen1)
random_in_unit_disk :: StdGen -> (V3 Double, StdGen)
random_in_unit_disk gen0 =
let
(rand1,gen1) = randomR (-1, 1) gen0 :: (Double, StdGen)
(rand2,gen2) = randomR (-1, 1) gen1 :: (Double, StdGen)
v = V3 rand1 rand2 0
in
if quadrance v >= 1 then
random_in_unit_sphere gen2
else
(v, gen2)
---------------
-- Utilities --
---------------
infinity :: RealFloat a => a
infinity = encodeFloat (floatRadix 0 - 1) (snd $ floatRange 0)
deg2rad :: Floating a => a -> a
deg2rad degrees = degrees * pi / 180
clamp :: (Ord a, Num a) => a -> a -> a -> a
clamp x y val = (max x).(min y) $ val
reflect :: V3 Double -> V3 Double -> V3 Double
reflect v n = v - (2 * (n `dot` v)) *^ n
refract :: V3 Double -> V3 Double -> Double -> V3 Double
refract uv n eta_over_eta' =
r_out_perp + r_out_parallel
where
cos_theta = min 1 ((-uv) `dot` n)
r_out_perp = eta_over_eta' *^ (uv + cos_theta *^ n)
r_out_parallel = (sqrt $ abs (1 - quadrance r_out_perp)) *^ (-n)
reflectance :: Double -> Double -> Double
reflectance cosine ref_idx =
let
r0 = ((1 - ref_idx) / (1 + ref_idx)) ^ 2
in
r0 + (1-r0) * (1 - cosine) ^ 5
near_zero (V3 r1 r2 r3) =
(abs(r1) < s) && (abs(r2) < s) && (abs(r3) < s)
where
s = 1.0E-7
--------------------
-- Material Class --
--------------------
type MaterialData = (MAT_Lambertian + MAT_Metal) + MAT_Dielectric
class Material a where
make_shared :: a -> MaterialData
scatter :: a -> Ray -> HitRecord -> StdGen -> (Maybe (Ray, V3 Double), StdGen)
instance (Material a, Material b) => Material (Either a b) where
make_shared = coPair(make_shared, make_shared)
scatter = coPair(scatter, scatter)
-- Lambertian
data MAT_Lambertian = MAT_Lambertian {
albedo_Lamb :: V3 Double
} deriving (Show)
instance Material MAT_Lambertian where
make_shared = inj1 -: inj1
scatter this r_in record gen =
let
(rand1, gen1) = random_unit_vector gen
scattered_direction =
if near_zero(normal record + rand1) then
normal record
else
normal record + rand1
scattered = Ray{orig = p record, dir = scattered_direction}
attenuation = albedo_Lamb this
in
(Just (scattered, attenuation), gen1)
-- Metal
data MAT_Metal = MAT_Metal {
albedo_Metal :: V3 Double,
fuzz :: Double
} deriving (Show)
instance Material MAT_Metal where
make_shared = inj2 -: inj1
scatter this r_in record gen =
let
(rand1, gen1) = random_in_unit_sphere gen
reflected = reflect (normalize $ dir r_in) (normal record)
scattered = Ray{orig = p record, dir = reflected + (fuzz this) *^ rand1}
attenuation = albedo_Metal this
in
if (dir scattered `dot` normal record) > 0 then
(Just (scattered, attenuation), gen1)
else
(Nothing, gen1)
-- Dielectric
data MAT_Dielectric = MAT_Dielectric {
ref_idx :: Double
} deriving (Show)
instance Material MAT_Dielectric where
make_shared = inj2
scatter this r_in record gen =
let
(rand1, gen1) = random gen
attenuation = V3 1 1 1
eta_over_eta' = if front_face record then 1 / ref_idx this else ref_idx this
unit_direction = normalize $ dir r_in
cos_theta = min 1 (-unit_direction `dot` normal record)
sin_theta = sqrt $ 1 - cos_theta ^ 2
cannot_refract = eta_over_eta' * sin_theta > 1
in
if cannot_refract || (reflectance cos_theta eta_over_eta' > rand1) then
let
reflected = reflect unit_direction (normal record)
scattered = Ray{orig = p record, dir = reflected}
in
(Just (scattered, attenuation), gen1)
else
let
refracted = refract unit_direction (normal record) eta_over_eta'
scattered = Ray{orig = p record, dir = refracted}
in
(Just (scattered, attenuation), gen)
getIntersection_forTorus :: (V3 Double, V3 Double, V3 Double, Double, Double, V3 Double) -> [Double]
getIntersection_forTorus = solveQuarticEq . genCoefficients
genCoefficients (x0,a,c,r1,r2,n) = (b4,b3,b2,b1,b0)
where
d0 = x0 - c
k = (r1^2) - (r2^2)
a_sq = quadrance a
d0_sq = quadrance d0
b4 = a_sq^2
b3 = 4*(d0 `dot` a)*a_sq
b2 = 2*d0_sq*a_sq+4*((d0 `dot` a)^2) + 2*k*a_sq - 4*(r1^2)*a_sq + 4*(r1^2)*(n `dot` a)^2
b1 = 4*d0_sq*(d0 `dot` a)+4*k*(d0 `dot` a) - 8*(r1^2)*(d0 `dot` a) + 8*(r1^2)*(n `dot` d0)*(n `dot` a)
b0 = d0_sq*d0_sq+2*k*d0_sq+k^2 - 4*(r1^2)*d0_sq + 4*(r1^2)*(n `dot` d0)^2
solveQuarticEq (a4,a3,a2,a1,a0) =
let
sol = do
(x_Re :+ x_Im) <- [x1,x2,x3,x4]
if (abs(x_Im) < 1.0E-8) && (1.0E-8 <= x_Re) then
return x_Re
else
[]
in
sol
where
l1 = (toCmp $ k3/4)/sqrt(k4)
l2 = (toCmp $ (cbrt(2)*k5)/(3*a4))/k8 + k8/(toCmp $ 3*cbrt(2)*a4)
l3 = (toCmp $ (a3^2)/(2*a4^2) - (4*a2)/(3*a4)) - l2
k1 = l1 + l3
k2 = -l1 + l3
k3 = -((a3/a4)^3) + (4*a2*a3)/(a4^2) - (8*a1)/a4
k4 = (toCmp $ ((a3/(2*a4))^2) - (2*a2)/(3*a4)) + l2
k5 = a2^2 - 3*a1*a3 + 12*a0*a4
k6 = 2*a2^3 - 9*a1*a2*a3 + 27*a0*a3^2 + 27*a1^2*a4 - 72*a0*a2*a4
k7 = -4*k5^3 + k6^2
k8 = cbrt((toCmp $ k6) + sqrt(toCmp $ k7))
l4 = toCmp $ -a3/(4*a4)
l5 = sqrt(k2)/2
l6 = sqrt(k1)/2
l7 = sqrt(k4)/2
x1 = l4 - l5 - l7
x2 = l4 + l5 - l7
x3 = l4 - l6 + l7
x4 = l4 + l6 + l7
cbrt x = x ** (1/3)
toCmp x = x :+ 0
prod x y = x >>= (\u -> zip (repeat u) y)
(-:) = flip (.)
type (+) a b = Either a b
inj1 :: a -> a + b
inj1 = Left
inj2 :: b -> a + b
inj2 = Right
coPair :: (a1 -> b, a2 -> b) -> (a1 + a2 -> b)
coPair = uncurry eitherタグ一覧: