Raytracer minimaliste en Haskell : Affichage d'une sphère

Pour commencer, nos objets auront des couleurs. Il nous faut donc quelque chose pour les représenter. Allons-y pour un type représentant la matière de nos objets.

module Materiel (
  Red,
  Green,
  Blue,
  Color,
  Material (Material),
  ) where

--Couleur, où 0.0 signifit "pas de couleur" et 1.0 "toute la couleur".
type Red = RealRep
type Green = RealRep
type Blue = RealRep
type Color = (Red, Green, Blue)

data Material = Material Color

Nous allons placer nos objets dans un module « Objects ». Nous allons définir un objet comme constitué de trois éléments : des informations sur sa position dans l'espace, des informations sur la matière dont il est constitué, et enfin une façon de calculer si un rayon intersecte l'objet, ou non.

odule Objects (
  Object (Object),
  Camera (Camera),
  Scene (Scene),
  Intersection,
  ) where

import Tools
import Material

--Représentation de l'observateur
data Camera = Camera Location Distance PlanSize

--Une scène
data Scene = Scene Camera [Object]

--Une intersection
data Intersection = Intersection

--On object raytracable
data Object = Object Location Material (Vector -> Scene -> Intersection)

On en profite pour déplacer le type Object dans le module Objects, et donc ses dépendances, et on pensera bien à ajouter l'import dans le module Raytracing.

Les types PlanSize et PixelLocation sont déplacés dans le module Tools.

Une sphère est définie par l'équation kitxmlcodeinlinelatexdvpx^2 + y^2 + z^2 = r^2finkitxmlcodeinlinelatexdvp où kitxmlcodeinlinelatexdvprfinkitxmlcodeinlinelatexdvp est le rayon de la sphère. Si l'on se place sur la position de la caméra, et que l'on se déplace le long de notre vecteur, construit dans la section précédente, il est possible qu'en avançant assez on « tombe » sur la sphère, c'est-à-dire que pour une bonne valeur de kitxmlcodeinlinelatexdvpkfinkitxmlcodeinlinelatexdvp, on aura les coordonnées du point kitxmlcodeinlinelatexdvp(o_x + k . r_x, o_y + k . r_y, o_z + k r_z) := \text{camera} + k . \text{rayon}finkitxmlcodeinlinelatexdvp solution de l'équation. Si l'on injecte les valeurs de la dernière équation dans celle de la sphère, on trouve un polynôme en kitxmlcodeinlinelatexdvpkfinkitxmlcodeinlinelatexdvp de degré 2.

Cette méthode a quelque chose de général. On définit un objet par une équation implicite, comme l'équation de la sphère, puis l'on injecte le point d'arrivée du rayon (qui nous est pour le moment inconnu). De cette façon, on se retrouve avec une équation à une variable. Si l'on sait la résoudre, les solutions, si elles existent, sont les points d'intersection du rayon et de l'objet. S'il n'y a pas de solution, c'est tout simplement que l'on ne peut pas voir l'objet.

Par bonheur, la sphère a la chance de nous donner une équation polynomiale de degré 2. On sait tous, depuis le lycée, que l'on peut calculer les solutions grâce au coefficient du polynôme. Ajoutons donc le nécessaire pour calculer les racines d'un polynôme de degré 2. Comme il peut y en avoir aucune, une (dite double), ou encore deux, nous stockerons les solutions dans une liste.

On en profite pour définir une intersection comme un point, la distance, et l'objet. On construit l'intersection à partir des distances. N'oublions pas que grâce à la paresse, elles ne seront calculées que si l'on a vraiment besoin de leur valeur.

--Une intersection
type Intersection = (Location, Distance, Object)

--On object raytracable
data Object = Object Location Material (Vector -> Scene -> [Distance])

type Radius = RealRep

--Cherche les intersections avec un sphère de rayon "radius" centré en (0, 0, 0)
sphereIntersect :: Radius -> Ray -> Scene -> [Distance]
sphereIntersect radius ray@(rx, ry, rz) scene = solvePoly2 (a, b, c)
  where
    Scene camera _ = scene
    Camera cameraLocation _ _ = camera
    (cx, cy, cz) = cameraLocation
    a = rx^2 + ry^2 + rz^2
    b = 2 * (rx * cx + ry * cy + rz * cz)
    c = cx^2 + cy^2 + cz^2 - radius^2

buildSphere :: Radius -> Location -> Material -> Object
buildSphere radius loc mat = Object loc mat (sphereIntersect radius)

module Tools (
  Location,
  Distance,
  RealRep,
  Vector,
  Ray,
  PlanSize,
  PlanWidth,
  PlanHeight,
  PixelLocation,
  normalise,
  int2RealRep,
  (+!),
  (-!),
  (*!),
  ) where

-- ...

(+!) :: Vector -> Vector -> Vector
(+!) (a, b, c) (a', b', c') = (a + a', b + b', c + c')
(-!) :: Vector -> Vector -> Vector
(-!) (a, b, c) (a', b', c') = (a - a', b - b', c - c')
(*!) :: RealRep -> Vector -> Vector
(*!) l (a, b, c) = (l * a, l* b, l * c)

--Les coefficients (a, b, c) d'un polynôme de la forme ax^2 + bx + c
type Polynom2 = (RealRep, RealRep, RealRep)
solvePoly2 :: Polynom2 -> [RealRep]
solvePoly2 (a, b, c)
  | discriminant < 0 = []
  | discriminant > 0 = [(-b + sqrtDiscriminant) / den, (-b - sqrtDiscriminant) / den]
  | otherwise = [-b / den]
  where
    discriminant = b^2 - 4 * a * c
    sqrtDiscriminant = sqrt discriminant
    den = 2 * a

On construit maintenant une fonction pour transformer la liste de distances en liste d'intersections.

distanceToIntersection :: Object -> Ray -> Distance -> Scene -> Intersection
distanceToIntersection obj ray objectDistance scene = (cameraLocation +! (objectDistance *! ray), objectDistance, obj)
  where
    Scene camera _ = scene
    Camera cameraLocation _ _ = camera

Construisons une vraie scène vide

emptyScene :: Camera -> Scene
emptyScene camera = Scene camera []

On modifie donc l'appel dans le main.

    pixels = [(x, y) | y <- [0..height - 1], x <- [0..width - 1]] >>= (getPixel . emptyScene $ buildDefaultCamera 300 width height)

Il ne nous reste plus qu'à récupérer la liste d'objets, calculer les intersections pour chacun, et garder la plus proche. On doit translater toute la scène, car les objets sont supposés se trouver en (0, 0, 0). Il est plus simple de raisonner ainsi, que ce soit pour calculer les intersections (les équations implicites sont plus simples pour les objets centrés à l'origine) ou encore l'orientation (plus tard on voudra pouvoir faire tourner des objets). On déplace donc la caméra.

--Calcule les intersections du rayon avec un objet.
computeIntersection :: Ray -> Object -> Scene -> [Intersection]
computeIntersection ray object scene = intersect ray >>= mapM (distanceToIntersection object ray) $ scene'
  where
    Object objectLocation _ intersect = object
    Scene camera objects = scene
    Camera cameraLocation distance planSize = camera
    camera' = Camera (cameraLocation -! objectLocation) distance planSize
    scene' = Scene camera' objects

--Applique le calcul d'intersections sur chacun des objets
computeIntersections :: Ray -> Scene -> [Intersection]
computeIntersections ray scene = concat (contextualIntersections scene)
  where Scene _ objects = scene
  contextualIntersections = sequence $ map (computeIntersection ray) objects

La fonction computeIntersection prend le rayon que nous avons lancé, et l'intersecte avec l'objet qu'elle reçoit. Elle réalise le calcul de l'intersection grâce à la fonction gentiment fournie par l'objet. Cette fonction fournit une liste de distances. On produit donc des intersections (contenant : la distance, la position de l'intersection, et l'objet, car nous aurons besoin de tout ça dans la suite) en mappant la fonction de conversion sur la liste.

Comme le typage est un peu compliqué, voici le détail. L'expression intersect ray est de type Scene -> [Distance], que l'on pensera comme du m [a]. L'expression distanceToIntersection object ray est elle de type Distance -> Scene -> Intersection que l'on pensera a -> m b. On voudrait appliquer map, mais alors on aurait quelque chose de type [a] -> [m b]. En appliquant mapM à cette expression, on la transforme en du [a] -> m [b]. On n'a plus qu'à utiliser l'opérateur >>= pour appliquer cette expression à l'intérieur de intersect ray. On convertit donc bien une liste de distances dépendant de la scène en une liste d'objets, dépendant de la scène.

Dans la seconde fonction, on doit convertir une liste d'objets en une liste d'intersections. Mais comme il y a quelques monades dans l'histoire, on utilise sequence pour convertir du [m [b]] en m [[b]].

Maintenant, on élimine toutes les intersections situées derrière la caméra, et on récupère l'intersection la plus proche de nous, si c'est possible.

--Ne conserve que les intersections devant la caméra, et récupère la plus proche
getClosestIntersection :: [Intersection] -> Scene -> Maybe Intersection
getClosestIntersection intersections = return intersection
  where
    allowedIntersections = filter isAhead intersections
    isAhead intersection@(_, distance, _) = distance > 0
    compareZ a b = compare ad bd
      where
        (_, ad, _) = a
        (_, bd, _) = b
    closestIntersection = minimumBy compareZ allowedIntersections
    intersection = if null allowedIntersections
                   then Nothing
                   else Just closestIntersection

On filtre et compare les intersections d'après leur distance à la caméra. On déconstruit le type Intersection pour récupérer les distances ad, bd et distance.

Il nous faut récupérer la couleur de l'objet depuis l'intersection, afin d'afficher un pixel de la bonne couleur. Comme nous n'avons pas encore de gestion de la lumière, on dira que plus l'objet est loin, moins il est lumineux. Ce n'est pas tout à fait vrai, mais on pourra obtenir une première image de cette façon.

import Material

-- ...

--Construit la couleur à partir d'une intersection
computeColor :: Intersection -> Color
computeColor (_, distance, Object _ material _) = (r / distance, g / distance, b / distance)
  where
    Material (r, g, b) = material

Et pour finir, convertissons cette couleur en pixel, sous la forme de trois Word8. On ajoute donc une fonction utilitaire dans le module Tools que l'on pensera bien à exporter (il faudra importer Data.Word). Il faudra aussi exporter la fonction.

realRep2Word8 :: RealRep -> Word8
realRep2Word8 = truncate

On peut alors ajouter la fonction de conversion dans le module Material. Il faudra importer les modules Tools et Data.Word.

colorToPixel :: Color -> [Word8]
colorToPixel (r, g, b) = map realRep2Word8 [r * 255, g * 255, b * 255, 255]

Nous pouvons maintenant écrire getPixels à partir de tout ce que nous avons produit.

getPixel :: Scene -> PixelLocation -> [Word8]
getPixel scene pixelLocation =  case color of
                                     Nothing -> [0, 0, 0, 255]
                                     Just pixels -> pixels
  where
    intersection = return pixelLocation
                   >>= computeRay
                   >>= computeIntersections
                   >>= getClosestIntersection
    color = fmap (colorToPixel . computeColor) (intersection scene)

On obtient à nouveau une image, mais elle est vide. Il faut rajouter des objets. Testons une simple sphère :

tableData :: Int -> Int -> UArray (Int,Int,Int) Word8
tableData width height = listArray ((0,0,0), (height-1, width-1, 3)) $ pixels
  where
    pixels = [(x, y) | y <- [0..height - 1], x <- [0..width - 1]]
                 >>= getPixel scene
    scene = Scene camera objects
    camera = buildDefaultCamera 300 width height
    objects = [buildSphere 1.0 (0, 0, -2) (Material (0.0, 0.0, 1.0))]

Voici donc le premier rendu :