雑記帳
僕用勉強ノート 「レイトレーシング」の巻

Haskell でレイトレーシングのチュートリアルを追いかける その2 - 法線ベクトル

引き続きこのサイトのチュートリアルに則って、レイトレーシングによる画像の生成に挑戦。
進捗状況としては、ひとまず「section 6」まで完了。
コードの実行結果
リンク先の結果と一緒なのでそちらを是非見てほしい
ソースコード
{-# LANGUAGE TypeOperators #-}

import Data.Char
import Data.Functor
import Control.Monad
import Control.Lens
import Linear.Vector
import Linear.Metric
import Linear.V3
import Linear.Quaternion

-- https://raytracing.github.io/books/RayTracingInOneWeekend.html section 6 with Haskell--


main = do
  let
    -- Image
    aspect_ratio = 16.0 / 9
    image_width = 256
    image_height = round $ fromInteger image_width / aspect_ratio
    -- World
    world = (([]
      `add` RT_Sphere{center = V3 0 0 (-1), radius = 0.5})
      `add` RT_Sphere{center = V3 0 (-100.5) (-1), radius = 100})
    -- Camera
    viewport_height = 2.0
    viewport_width = aspect_ratio * viewport_height
    focal_length = 1.0
    origin = zero
    horizontal = viewport_width *^ unit _x
    vertical = viewport_height *^ unit _y
    lower_left_corner = origin - horizontal ^/2 - vertical ^/2 - focal_length *^ unit _z


  -- Render

  img_data <- return $ "P3\n" ++ show image_width ++ " " ++ show image_height ++ "\n255\n"
  putStr $ img_data

  foldr (>>) (return ()) $ (fmap $ ($) $ \(j, i) ->
    let
      u = fromInteger i / (fromInteger image_width - 1.0)
      v = fromInteger j / (fromInteger image_height - 1.0)
      r = Ray {orig = origin, dir = lower_left_corner + u *^ horizontal + v *^ vertical - origin}
      pixcel_color = ray_color r world
    in
      write_color $ pixcel_color) $
      (,) <$> [image_height - 1, image_height - 2 .. 0] <*> [0, 1 .. image_width - 1]



---------------------
-- A Ray Data Type --
---------------------

data Ray = Ray {
  orig :: V3 Double,  -- Origin of this ray (As a position in 3D Euclidean space)
  dir :: V3 Double    -- Direction of this ray (As a direction vector in 3D Euclidean space)
} deriving (Show)


at' :: Ray -> Double -> V3 Double
at' r t = (orig r) + t *^ (dir r)


----------------------
-- Hittable objects --
----------------------

data HitRecord = HitRecord {
  p :: V3 Double,
  normal :: V3 Double,
  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)
  }

data RT_Sphere = RT_Sphere {
  center :: V3 Double,
  radius :: Double
} deriving (Show)


type HittableObjects = (RT_Sphere + RT_Sphere) + RT_Sphere -- Second and third RT_Sphere are just dummies

class Hittable a where
  toSum :: a -> HittableObjects
  hit :: a -> Ray -> Double -> Double -> Maybe HitRecord

instance (Hittable a, Hittable b) => Hittable (a + b) where
  toSum = coPair(toSum, toSum)
  hit = coPair(hit, hit)

instance Hittable RT_Sphere where
  toSum = Inj1 -: Inj1
  hit obj r t_min t_max =
    let
      oc = orig r - center obj
      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
                    } r ((at' r x - center obj) ^/ radius obj)
                else
                  f $ xs

              [] ->
                Nothing

        in
          f $ [(-half_b - root) / a, (-half_b + root) / a]
      else
        Nothing

add :: Hittable a => [HittableObjects] -> a -> [HittableObjects]
add list obj = (toSum obj) : list


hitSomething :: [HittableObjects] -> 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 t_max
        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)

---------------
-- Utilities --
---------------

write_color :: RealFrac a => V3 a -> IO ()
write_color v =
  let
    f = show.floor.(255.999*)
  in
    do
      tmp <- return $ f(v ^._x) ++ " " ++ f(v ^._y) ++ " " ++ f(v ^._z) ++ "\n"
      putStr $ tmp


ray_color :: Ray -> [HittableObjects] -> V3 Double
ray_color r objects =
  let
    record = hitSomething objects r 0 10000
  in
    case record of
      Just a ->
        0.5 *^ ((normal a) + (V3 1 1 1))

      Nothing ->
        let
          unit_direction = normalize $ (dir r)
          s = 0.5 * (unit_direction ^._y + 1.0)
        in
          lerp s (V3 0.5 0.7 1.0) (V3 1.0 1.0 1.0)


-- Diagrammatic-order composition
(-:) = flip (.)

-- Sum objects and injections
data (+) a b = Inj1 a | Inj2 b

instance (Show a, Show b) => Show (a + b) where
  show = coPair(show, show)

-- Dual concept of pairs
coPair :: (a1 -> b, a2 -> b) -> (a1 + a2 -> b)
coPair (f, g) x = case x of
  Inj1 x -> f x
  Inj2 x -> g x