Raytracer minimaliste en Haskell : Gestion des lumières

Cochoy Jérémy

Deuxième article d'une série visant la réalisation d'un raytracer minimaliste en haskell. Quelques connaissances de base sur les monades et la manipulation des listes seront nécessaires, ainsi que quelques notions de mathématiques comme le produit scalaire, l'angle formé par deux vecteurs, etc.

Votre avis et vos suggestions nous intéressent !
Alors après votre lecture, n'hésitez pas  : Commentez

- Titre : Raytracer minimaliste en Haskell : Gestion des lumières

- Auteur : Cochoy Jérémy

- Parution : 3 juin 2013

- Mise à jour : 7 juin 2013

- Public visé :

- Licence : Les sources présentées sur cette page sont libres de droits et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une œuvre intellectuelle protégée par les droits d'auteur. Copyright ® 2013 Zenol. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts. Droits de diffusion permanents accordés à Developpez LLC.

Sommaire

- I. Introduction

- II. Lumières : Révision du calcul de luminosité

- II-1. La lumière ambiante

- II-2. Calcul de la lumière

- II-3. Exposition

- III. Lumières : les sources lumineuses

- III-1. Les normales

- III-2. Lumière diffuse

- III-3. Lumière spéculaire

- IV. Les ombres

- V. Anti-aliasing : supersampling stochastique

- VI. Remerciements

I - Introduction

Cet article fait suite à l'affichage d'une sphère. Nous allons repartir de ces sources pour ajouter des sources de lumières, des ombres, et de l'anti-aliasing.

Vous pouvez trouver les sources de cet article en cliquant sur ce lien.

II - Lumières : Révision du calcul de luminosité

II-1 - La lumière ambiante

Nous allons maintenant introduire les sources de lumières. Il y a diverses sources possibles : projecteurs, boules de feu, lumière éloignée (soleil)... On peut aussi souhaiter que la lumière décroisse lorsque l'on s'éloigne de la source.

Nous allons, de façon tout à fait arbitraire, introduire les sources lumineuses suivantes :

- Une lumière ambiante, qui sera appliquée à tous les objets. On fait le choix que la partie réfléchie reçue sera multipliée par le coefficient où les coefficients a, b et c pourront être librement choisis, et d est la distance de l'observateur à l'objet. Cela signifie que plus l'observateur sera loin, moins la luminosité sera intense. Une lumière ambiante n'a pas de position précise dans l'espace, donc pas de coordonnées.

- Les lumières directionnelles. Tous les rayons sont parallèles, et « proviennent de l'infini ». C'est une bonne approximation des rayons du soleil. On dira qu'une telle lumière ne s'atténue pas. Une lumière directionnelle n'a pas non plus de position dans l'espace.

- Les spots, ou projecteurs. Ce sont des sources de lumières directionnelles en forme de cône. Plus on est dans l'axe, plus la zone est éclairée. L'angle du cône de lumière sera nécessairement aigu (afin de rendre l'implémentation plus simple). On disposera de trois coefficients a, b et c pour multiplier le résultat par où e est la distance du projecteur à l'objet.

- Les sources omnidirectionnelles. Ce sont des sortes de boules de lumière qui éclairent dans toutes les directions. À nouveau, ces lumières disposeront d'une atténuation fonction de leur distance.

Vous remarquerez que pour toutes les lumières disposant de trois coefficients (a, b, c) pour contrôler l'atténuation, le choix (0, 0, 1) nous donne une lumière non atténuée.

Tout comme pour nos objets, une lumière sera constituée de ses caractéristiques, et d'une fonction de calcul de luminosité. Mais puisque la seule caractéristique commune à nos lumières est leur couleur, et que l'on n'en a besoin que pour calculer la luminosité, on n'a pas besoin de la stocker.

Lights.hs

module Lights (
Light (Light),
Attenuation,
computeAttenuation,
buildAmbiantLight,
buildDirectionalLight,
buildOmnidirectionalLight,
buildSpotLight,
) where

import Tools
import Objects

data Light = Light (Intersection -> Scene -> Color)

-- Les facteurs a, b et c de l'aténuation, qui est un polynome en d
type AConstant = RealRep
type ALinear = RealRep
type ASquare = RealRep

-- Le triplet (a, b, c) utilisé pour calculer 1/(ad^2 + b*d + c)
type Attenuation = (ASquare, ALinear, AConstant)

computeAttenuation :: Attenuation -> Distance -> RealRep
computeAttenuation (a, b, c) d = 1 / (a * d ^ 2 + b * d + c)

buildAmbiantLight = undefined
buildDirectionalLight = undefined
buildOmnidirectionalLight = undefined
buildSpotLight = undefined

Toutes les fonctions exportées serviront à construire les différents types de lumières listés quelques paragraphes plus tôt. Commençons par le plus simple : la lumière ambiante.

Il nous faut pouvoir multiplier deux à deux les coordonnées d'un vecteur. On rajoute donc une fonction utilitaire.

Tools.hs

vectorProduct :: Vector -> Vector -> Vector
vectorProduct (a, b, c) (x, y, z) = (a * x, b * y, c * z)

On calcule maintenant la lumière ambiante réfléchie au niveau de l'intersection.

Lights.hs

buildAmbiantLight :: Attenuation -> Color -> Light
buildAmbiantLight attenuation lightColor = Light lightFunction
where
lightFunction intersection scene = (computeAttenuation attenuation distance) *! light
where
light = lightColor `vectorProduct` objectColor
Material objectColor = material
(_, distance, Object _ material _) = intersection

Afin de pouvoir tester petit à petit nos lumières, et les voir naître sous nos yeux, on va maintenant mettre en place l'application des lumières aux objets de la scène avant de poursuivre.

II-2 - Calcul de la lumière

Nous nous retrouvons confrontés à un souci de dépendances cycliques. On voudrait le type objet dans Object.hs, le type Light dans Light.hs, et le type Scene dans Raytracer.hs. Seulement Scene dépend de Light et Object. Comme ghc ne supporte pas les dépendances cycliques, nous avons deux solutions. L'une utilise un système de « bootstraping », l'autre consiste à placer tous les types dans un Type.hs. Pour des raisons de simplicité et maintenabilité nous choisissons la seconde.

Types.hs

module Types (
Red,
Green,
Blue,
Color,
Camera (Camera),
Scene (Scene),
Intersection,
Object (Object),
Light (Light),
Material (Material),
) where

import Tools

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

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

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

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

--Une lumière
data Light = Light (Intersection -> Scene -> Color)

--Couleur, 0.0 signifie "pas de couleur" et 1.0 "toute la couleur".
type Red = RealRep
type Green = RealRep
type Blue = RealRep

type Color = (Red, Green, Blue)

-- Une matière
data Material = Material Color

Nous ajoutons maintenant une liste de lumières :

Types.hs

data Scene = Scene Camera [Object] [Light]

On met à jour les diverses apparitions du constructeur Scene, en ajoutant un troisième argument (_ s'il s'agit de décomposer une scène, ou [] s‘il s'agit de construire une scène). Dans computeItersection, on doit récupérer les lumières pour construire scene.

Raytracer.hs

computeIntersection ray object scene =  -- ...
-- ...
Scene camera objects lights = scene
-- ...
scene' = Scene camera' objects lights

Il nous faut maintenant revoir computeColor pour utiliser cette liste de lumières. On commence par ajouter une fonction pour sommer les couleurs provenant d'une liste de couleurs. N'oublions pas de l'exporter. Comme on se doute que l'on va devoir sommer beaucoup de valeurs, on ne veut pas construire une liste de thunk qui seront longs à évaluer, et l'on préfère donc utiliser foldl' plutôt que foldl.

Tools.hs

import Data.List (foldl')
-- ...
vectorSum :: [Vector] -> Vector
vectorSum = foldl' (+!) (0, 0, 0)

On revoit maintenant computeColor pour qu'elle calcule chacune des lumières sur l'intersection, puis somme les couleurs produites. Cela implique de changer son type.

Raytracer.hs

--Construit la couleur à partir d'une intersection
computeColor :: Intersection -> Scene -> Color
computeColor intersection scene = color
where
colors = map (\ (Light f) -> f intersection scene) lights
color = vectorSum colors
Scene _ _ lights = scene

Il nous faut donc modifier les appels. Par chance, il n'y en a qu'un.

Raytracer.hs

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 intersectionToPixel (intersection scene)
intersectionToPixel intersection = colorToPixel . computeColor intersection $ scene
Scene _ _ lights = scene

Si l'on lance maintenant notre raytracer, on n'obtient qu'un écran noir. C'est que l'on doit ajouter des lumières à notre scène pour y voir quelque chose. On importe Lights dans Main.hs, et on modifie donc dans tableData la ligne suivante :

Main.hs

scene = Scene camera objects [buildAmbiantLight (0, 1, 0) (1, 1, 1)]

On retrouve donc exactement la lumière que l'on avait dans l'article précédent (car , , , donc on multiplie par ).

II-3 - Exposition

On peut s'amuser à ajouter plus de lumières. Ajoutons un effet torche avec une lumière ambiante jaune (1, 1, 0) ainsi qu'une lumière ambiante blanche (1, 1, 1) "plate" sur une sphère blanche.

Main.hs

lights = [buildAmbiantLight (0, 0, 1) (1, 1, 1), buildAmbiantLight (0, 1, 0) (1, 1, 0)]
objects = [buildSphere 1.0 (0, 0, -2) (Material (1.0, 1.0, 1.0))]
scene = Scene camera objects lights

Nom d'un phoque ! D'où sort ce bleu ? C'est bien sur un overflow dû au fait que l'on multiplie 255 par une valeur supérieure à 1 dans colorToPixel, elle-même appelée par getPixel. En fait, on voulait garantir que computeColor ne produise jamais des valeurs supérieures à 1.

Il y a deux solutions simples pour gérer cette situation. D'abord, on pourrait diviser par la somme de toutes les lumières de la scène. Ça fonctionnerait, mais on obtiendrait une image très sombre, surtout dans les zones où seulement une partie des lumières participent. Une autre solution, simple et efficace, est d'appliquer une fonction avec une courbe écrasée et tendant vers 1 en plus l'infini. Il faut aussi qu'elle soit linéaire vers les valeurs proches de 0, et qu'elle vaille 0 en 0 (en d'autres termes, on veut qu'elle soit équivalente à en 0). Une fonction qui répond à nos attentes est . Le choix de lambda dépend des scènes, donc on laissera cette information dans la caméra, sous le nom "d'exposition".

On rajoute donc le champ dans le type Camera

Types.hs

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

--L'exposition
type Exposition = RealRep

On s'occupe de modifier tous les appels au constructeur de Camera. Dans buildDefaultCamera, on mettra une exposition de 1. Ensuite, on peut revoir la fonction computeColor.

Types.hs

--Construit la couleur à partir d'une intersection
computeColor :: Intersection -> Scene -> Color
computeColor intersection scene = vectorApply (\x -> 1 - exp(-exposition*x)) color
where
colors = map (\ (Light f) -> f intersection scene) lights
color = vectorSum colors
Scene camera _ lights = scene
Camera _ _ _ exposition = camera

Le résultat obtenu est bien plus convainquant :

III - Lumières : les sources lumineuses

III-1 - Les normales

Pour un bon éclairage, on a besoin de trois lumières : la lumière ambiante, la lumière diffuse, et la lumière spéculaire (cf http://en.wikipedia.org/wiki/Phong_shading). Pour la lumière ambiante, nous avons fait le choix de créer des sources de lumières particulières. Il nous reste à traiter les cas diffus et spéculaire. Ces deux cas dépendent de la direction de la source de lumière, et cette direction varie pour chacune de nos sources lumineuses. Pour les calculer, il faut connaître la normale à la surface. La normale est un vecteur unitaire qui indique la direction extérieure à la surface. Bien sûr, l'extérieur et l'intérieur d'une surface sont bien subjectifs. Par exemple, pour un plan, la normale est le long d'une droite orthogonale à celui-ci. Il faut que la normale à la sphère pointe "vers" la caméra, et donc selon le cas, c'est le vecteur normal que l'on veut avoir. On fera donc un petit ajustement au moment des calculs, pour choisir le bon vecteur. Modifions donc nos objets pour qu'ils fournissent une normale, à partir d'une intersection.

On commence par modifier nos intersections. Il faut bien différencier les coordonnées dans le repère local (centré sur) de l'objet (celles que l'on calculait jusque-là) et les coordonnées dans notre scène (qui pourront nous être utiles pour calculer la lumière).

Types.hs

--Représentation d'une position dans un repère local
type LocalLoc = Location

Types.hs

--Une intersection
type Normal = vector
type Intersection = Intersection Location LocalLoc Distance Object Normal

L'ajout de la normale nous servira plus tard. Pour le moment, on placera le vecteur nul dans les constructeurs. On utilise les messages du compilateur pour retrouver les utilisations du type Intersection et les modifier, et en exportant bien le constructeur de type.

Object.hs

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

localIntersectionToGlobal :: Intersection -> Intersection
localIntersectionToGlobal intersection = Intersection (objectLocation +! localLoc) localLoc distance object normal
where
Intersection location localLoc distance object normal = intersection
Object objectLocation _ _ _ = object

En effet, à l'appel de distanceToIntersection, le repère local est le repère utilisé. On ajoute la possibilité, avec l'ancienne caméra, de recalculer la position de l'intersection.

Raytracer.hs

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

On modifie le type Object pour ajouter le calcul de normal (et l'on modifie les appels du constructeur) :

Type.hs

--Un objet raytraçable
type IntersectionFunction = Ray -> Scene -> [Distance]
type NormalFunction = Object -> LocalLoc -> Vector
data Object = Object Location Material IntersectionFunction NormalFunction

Calculons les normales de la sphère :

Objects.hs

--Détermine les normales de la sphère
sphereNormal :: Object -> LocalLoc -> Vector
sphereNormal object localLoc = normalise localLoc

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

En règle générale, quand un objet (comme la sphère) est défini de façon implicite, les normales s'obtiennent facilement en calculant le gradient. C'est-à-dire que l'on dérive l'équation par rapport à X, puis par rapport à Y, et enfin par rapport à Z. On obtient donc trois formules. Chacune nous donne une coordonnée du vecteur normal (respectivement la coordonnée en X, en Y et en Z). Il se trouve que pour la sphère, on trouve le vecteur . Ce sont donc deux fois le vecteur coordonnées, et l'on peut alors directement normaliser l'intersection (en coordonnées locales).

Pour finir, on utilise la fonction de calcul de normales dans le calcul d'une intersection :

Objects.hs

--Choisit la bonne direction de la normale, dans le repère local de l'objet
selectNormalSign :: Camera -> Normal -> Normal
selectNormalSign camera normal = if aligned then normal else (-1) *! normal
where
aligned = (cameraLocation .! normal) >= 0
Camera cameraLocation _ _ _ = camera

--Compute intersection's location, from local rep
distanceToIntersection :: Object -> Ray -> Distance -> Scene -> Intersection
distanceToIntersection object ray objectDistance scene = Intersection intersectionLocation intersectionLocation objectDistance object (selectNormalSign camera normal)
where
normal = (normalFunction object intersectionLocation)
Scene camera _ _ = scene
Camera cameraLocation _ _ _ = camera
intersectionLocation = cameraLocation +! (objectDistance *! ray)
Object _ _ _ normalFunction = object

Comme annoncé, on veut que la normale pointe vers la caméra, et on obtient cela grâce à selectNormalSign. De cette façon, pas de bug de lumière parce que la normale est dans le mauvais sens, et qu'elle vous fait dos alors qu'elle devrait vous faire face. La normale sera toujours dans le bon sens.

III-2 - Lumière diffuse

Commençons par la fonction de calcul de la lumière diffuse. Celle-ci repose sur l'approximation selon laquelle la lumière est réfléchie équitablement dans toutes les directions, selon l'angle entre la normale à la surface et la direction de la lumière.

Tools.hs

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

Lights.hs

type UnaryLightVector = Vector
type UnaryNormalVector = Vector
-- Le type "type DiffuseFactor = Vector" sera placé dans Types.hs
-- Calcule la lumière diffuse à partir de deux rayons provenant de l'intersection et pointant respectivement
-- vers la source de lumière et vers "l'extérieur" de l'objet.
computeDiffuse :: DiffuseFactor -> UnaryLightVector -> UnaryNormalVector -> RealRep
computeDiffuse diffuseFactor unaryLightVector unaryNormalVector = diffuseValue
where
diffuseValue = diffuseFactor * diffuseFactor
diffuseFactor = max unaryLightVector .! unaryNormalVector 0

Le calcul de lumière tient donc en une ligne : on pondère par le cosinus de l'angle entre le vecteur Objet -> Lumière, nommé unaryLightVector, et la normale de l'objet, nommé unaryNormalVector.

On ajoute un caractère diffus à nos objets. Remarquez que l'on ne garde le facteur diffuseValue (correspondant au cosinus de l'angle entre la source de lumière et la normale) que s‘il est positif (une lumière négative n'a pas de sens). Dans le main, on mettra la diffusion et la lumière spéculaire (les deux champs) de notre sphère à 1, histoire de pouvoir compiler.

Types.hs

-- Une matière
type DiffuseFactor = RealRep
type SpecularFactor = RealRep
data Material = Material Color DiffuseFactor SpecularFactor SpecularExponent

Ajoutons maintenant la source de lumière directionnelle.

Lights.hs

type LightDirection = Vector
buildDirectionalLight :: LightDirection -> Color -> Light
buildDirectionalLight lightDirection lightColor = Light lightFunction
where
unaryLightVector = normalise $ (-1) *! lightDirection
lightFunction intersection scene = (diffuseLightFactor + specularLightFactor) *! color
where
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = 0
Intersection _ _ _ object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

On a préparé l'emplacement pour le calcul de la lumière spéculaire, que l'on fixe pour le moment à 0.

Je vous propose une nouvelle scène, certes pas très esthétique, mais parfaite pour vérifier que nos calculs et nos lumières fonctionnent bien.

Raytracer.hs

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
camera = buildDefaultCamera 600 width height
lights = [buildAmbiantLight (0, 0, 2) (1, 1, 1), buildDirectionalLight (1, -1, -3) (1, 1, 1)]
objects = [buildSphere 1.0 (0, 1.2, -7) (Material (1.0, 1.0, 1.0) 1 1 1),
buildSphere 1.0 (0, -1.2, -7) (Material (1.0, 0.0, 1.0) 1 1 1),
buildSphere 1.0 (1.2, 0, -7) (Material (0.0, 1.0, 1.0) 1 1 1),
buildSphere 1.0 (-1.5, 0, -7) (Material (0.0, 0.0, 1.0) 1 1 1),
buildSphere 0.5 (0, 2, -6.6) (Material (1.0, 1.0, 0.0) 1 1 1)]
scene = Scene camera objects lights

On passe maintenant à la lumière omnidirectionnelle.

Main.hs

buildOmnidirectionalLight :: Attenuation -> Location -> Color -> Light
buildOmnidirectionalLight attenuation lightLocation lightColor = Light lightFunction
where
lightFunction intersection scene = factor *! color
where
factor = attenuationFactor * (diffuseLightFactor + specularLightFactor)
unaryLightVector = normalise (lightLocation -! (intersectionLocation))
attenuationFactor = computeAttenuation attenuation distance
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = 0
Intersection intersectionLocation _ distance object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

Voici les nouvelles lumières, pour tester :

Main.hs

lights = [buildAmbiantLight (0, 0, 2) (1, 1, 1), buildOmnidirectionalLight (0, 0, 0.8) (0, 0, -5) (1, 1, 1)]

Ajoutons maintenant notre spot lumineux :

Types.hs

--Un angle
type RadianAngle = RealRep

Lights.hs

buildSpotLight :: LightDirection -> RadianAngle -> Attenuation -> Location -> Color -> Light
buildSpotLight lightDirection coneAngle attenuation lightLocation lightColor = Light lightFunction
where
unaryLightDirection = normalise $ (-1) *! lightDirection
lightFunction intersection scene = factor *! color
where
unaryLightVector = normalise (lightLocation -! intersectionLocation)
spotAngle = acos(unaryLightVector .! unaryLightDirection)
halfConeAngle = coneAngle / 2
spotFactor = if spotAngle < halfConeAngle
then cos(pi * spotAngle / coneAngle) -- = cos(2 * pi * spotAngle / halfConeAngle)
else 0
factor = spotFactor * attenuationFactor * (diffuseLightFactor + specularLightFactor)
attenuationFactor = computeAttenuation attenuation distance
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = 0
Intersection intersectionLocation _ distance object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

Pour tester, vous pouvez utiliser les lumières suivantes :

Main.hs

lights = [buildAmbiantLight (0, 0, 15) (1, 1, 1), buildSpotLight (1, -1, -13) (pi/6) (0, 0, 1) (-1, 1, 0) (1, 1, 1)]

Pour calculer un "beau cône", on commence par récupérer l'angle entre la direction du spot et la position de l'objet. Si l'angle est inférieur à l'angle maximal (comme on regarde le long de la direction du spot, on doit diviser l'angle d'ouverture par deux) alors on calcule un "dégradé". Pour ce faire, on utilise simplement le cosinus de l'angle "si jamais le cône avait une ouverture de 180 degrés". D'où le pi/2 * spotAngle / halfConeAngle. On fait donc le cosinus d'une interpolation linéaire.

III-3 - Lumière spéculaire

Il ne nous reste plus grand-chose à faire pour gérer aussi la lumière spéculaire.

La lumière spéculaire est un petit halo qui se dessine sur les objets et leur donne un aspect « poli ». Pour produire cet effet, on va calculer la réflexion du rayon de lumière incident (normalisé) sur l'objet. Si sont respectivement le rayon réfléchi que l'on cherche, moins le rayon de lumière incident, et la normale, alors on obtient par la formule . Ensuite, on prend le cosinus de l'angle entre le rayon réfléchi et le vecteur Objet -> Caméra normalisé. On obtient alors un terme représentant la lumière spéculaire. Il nous faut ensuite l'adapter légèrement selon les caractéristiques de la matière. On utilise donc les deux coefficients specularFactor et specularExponent représentant respectivement l'intensité lumineuse spéculaire et la « rugosité » de la matière. Pour plus d'informations vous avez la page wikipedia.

Lights.hs

type UnaryCameraVector = Vector
computeSpecular :: SpecularFactor -> SpecularExponent -> UnaryCameraVector -> UnaryLightVector -> UnaryNormalVector -> RealRep
computeSpecular specularFactor specularExponent unaryCameraVector unaryLightVector unaryNormalVector = specularValue
where
specularValue = specularFactor * (phongTerm ** specularExponent)
reflect = 2 * (unaryLightVector .! unaryNormalVector)
phongVector = (reflect *! unaryNormalVector) -! unaryLightVector
phongTerm = max (phongVector .! unaryCameraVector) 0

On implémente alors les lumières manquantes :

Lights.hs

type LightDirection = Vector
buildDirectionalLight :: LightDirection -> Color -> Light
buildDirectionalLight lightDirection lightColor = Light lightFunction
where
unaryLightVector = normalise $ (-1) *! lightDirection
lightFunction intersection scene = (diffuseLightFactor + specularLightFactor) *! color
where
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = computeSpecular specularFactor specularExponent unaryCameraVector unaryLightVector normal
unaryCameraVector = normalise $ (cameraLocation -! intersectionLocation)
Scene camera _ _ = scene
Camera cameraLocation _ _ _ = camera
Intersection intersectionLocation _ _ object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

Lights.hs

buildOmnidirectionalLight :: Attenuation -> Location -> Color -> Light
buildOmnidirectionalLight attenuation lightLocation lightColor = Light lightFunction
where
lightFunction intersection scene = factor *! color
where
factor = attenuationFactor * (diffuseLightFactor + specularLightFactor)
unaryLightVector = normalise (lightLocation -! (intersectionLocation))
attenuationFactor = computeAttenuation attenuation distance
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = computeSpecular specularFactor specularExponent unaryCameraVector unaryLightVector normal
unaryCameraVector = normalise $ (cameraLocation -! intersectionLocation)
Scene camera _ _ = scene
Camera cameraLocation _ _ _ = camera
Intersection intersectionLocation _ distance object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

Lights.hs

buildSpotLight :: LightDirection -> RadianAngle -> Attenuation -> Location -> Color -> Light
buildSpotLight lightDirection coneAngle attenuation lightLocation lightColor = Light lightFunction
where
unaryLightDirection = normalise $ (-1) *! lightDirection
lightFunction intersection scene = factor *! color
where
unaryLightVector = normalise (lightLocation -! intersectionLocation)
spotAngle = acos(unaryLightVector .! unaryLightDirection)
halfConeAngle = coneAngle / 2
spotFactor = if spotAngle < halfConeAngle
then cos(pi * spotAngle / coneAngle) -- = cos(2 * pi * spotAngle / halfConeAngle)
else 0
factor = spotFactor * attenuationFactor * (diffuseLightFactor + specularLightFactor)
attenuationFactor = computeAttenuation attenuation distance
color = lightColor `vectorProduct` objectColor
diffuseLightFactor = computeDiffuse objectDiffuseFactor unaryLightVector normal
specularLightFactor = computeSpecular specularFactor specularExponent unaryCameraVector unaryLightVector normal
unaryCameraVector = normalise $ (cameraLocation -! intersectionLocation)
Scene camera _ _ = scene
Camera cameraLocation _ _ _ = camera
Intersection intersectionLocation _ distance object normal = intersection
Object _ material _ normalFunction = object
Material objectColor objectDiffuseFactor specularFactor specularExponent = material

Pour illustrer le résultat, je vous propose de revoir les réglages de la scène :

Main.hs

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
camera = buildDefaultCamera 600 width height
lights = [buildAmbiantLight (0, 0, 15) (1, 1, 1), buildSpotLight (1, -1, -13) (pi/6) (0, 0, 1) (-1, 1, 0) (1, 1, 1)]
objects = [buildSphere 1.0 (0, 1.2, -7) (Material (1.0, 1.0, 1.0) 1 2 6),
buildSphere 1.0 (0, -1.2, -7) (Material (1.0, 0.0, 1.0) 1 2 50),
buildSphere 1.0 (1.2, 0, -7) (Material (0.0, 1.0, 1.0) 1 2 12),
buildSphere 1.0 (-1.5, 0, -7) (Material (0.0, 0.0, 1.0) 1 1 12),
buildSphere 0.5 (0, 2, -6.6) (Material (1.0, 1.0, 0.0) 1 1 12)]
scene = Scene camera objects lights

Et voici le résultat :

IV - Les ombres

Pour commencer, nous déplaçons les fonctions de calcul d'intersection pour qu'elles puissent être importées dans le module Lights.

Main.hs

module Intersections (
computeIntersection,
computeIntersections,
) where

import Types
import Tools
import Objects

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

--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

Commençons avec la lumière correctionnelle. Tout objet sur la trajectoire obtenue en remontant le rayon de lumière crée une ombre. Il faut donc détecter si un objet se trouve dans cette direction. Si oui, on ne doit produire aucune lumière. On va donc lancer un rayon depuis l'intersection.

On commence par exporter une constante, signifiant que toute distance inférieure à celle-ci est considérée comme égale à 0. Cela va nous permettre d'éviter de ré-intersecter la position depuis laquelle on lance notre rayon.

Tools.hs

epsilon :: RealRep
epsilon = 1e-12

Maintenant, on veut pouvoir relancer un rayon depuis une intersection. On ajoute donc une fonction spécifique dans le module Intersections. On constate ici que cela consiste simplement à modifier la position de la caméra, et de filtrer les intersections obtenues.

Intersections.hs

--Recherche les intersections depuis la position location vers la direction ray
lookFrom :: Location -> Ray -> Scene -> [Intersection]
lookFrom location ray scene = filter (\(Intersection _ _ distance _ _) -> distance > epsilon ) $ computeIntersections ray scene'
where
Scene camera objects light = scene
Camera cameraLocation cameraDistance cameraPlanSize cameraExposition = camera
camera' = Camera location cameraDistance cameraPlanSize cameraExposition
scene' = Scene camera' objects light

On modifie maintenant le calcul de lumière en rajoutant une instruction conditionnelle. Si la liste obtenue par lookFrom en regardant dans la direction de la source lumineuse est vide, on peut calculer la lumière, sinon on renvoie (0, 0, 0).

Light.hs

buildDirectionalLight lightDirection lightColor = Light lightFunction
where
unaryLightVector = normalise $ (-1) *! lightDirection
lightFunction intersection scene = if hidden then (0, 0, 0) else computedLight
where
hidden = not . null $ lookFrom intersectionLocation unaryLightVector scene
computedLight = (diffuseLightFactor + specularLightFactor) *! color
-- ...

Au tour maintenant de l'ombre du spot. Pour celle-ci, on ne doit pas compter les objets derrière le spot. On filtre donc pour exclure les intersections trop éloignées. Puisqu'on fera exactement la même chose pour la lumière omnidirectionnelle, on ajoute une nouvelle fonction.

Lights.hs

isHidden :: Distance -> Location -> UnaryLightVector -> Scene -> Bool
isHidden lightDistance location unaryLightVector scene = not . null $ filteredIntersections
where
intersections = lookFrom location unaryLightVector scene
filteredIntersections = filter (\ (Intersection _ _ distance _ _) -> distance < lightDistance) intersections

buildSpotLight :: LightDirection -> RadianAngle -> Attenuation -> Location -> Color -> Light
buildSpotLight lightDirection coneAngle attenuation lightLocation lightColor = Light lightFunction
where
lightDistance = norm lightDirection
unaryLightDirection = normalise $ (-1) *! lightDirection
lightFunction intersection scene = if hidden then (0, 0, 0) else computedLight
where
hidden = isHidden lightDistance intersectionLocation unaryLightVector scene
computedLight = factor *! color

Pour ce qui est de la source omnidirectionnelle, c'est identique au précédent. Il nous faut juste obtenir le vecteur pointant vers la source et la distance de la source à l'intersection.

Lights.hs

buildOmnidirectionalLight :: Attenuation -> Location -> Color -> Light
buildOmnidirectionalLight attenuation lightLocation lightColor = Light lightFunction
where
lightFunction intersection scene = if hidden then (0, 0, 0) else computedLight
where
lightDistance = norm lightDirection
hidden = isHidden lightDistance intersectionLocation unaryLightVector scene
computedLight = factor *! color
lightDirection = lightLocation -! intersectionLocation
unaryLightVector = normalise lightDirection
-- ...

Pour tester les ombres, j'ai rajouté quelques objets à notre scène précédente :

Main.hs

objects = [buildSphere 1.0 (0, 1.2, -7) (Material (1.0, 1.0, 1.0) 1 2 6),
buildSphere 1.0 (0, -1.2, -7) (Material (1.0, 0.0, 1.0) 1 2 50),
buildSphere 1.0 (1.2, 0, -7) (Material (0.0, 1.0, 1.0) 1 2 12),
buildSphere 1.0 (-1.5, 0, -7) (Material (0.0, 0.0, 1.0) 1 1 12),
buildSphere 0.5 (0, 2, -6.6) (Material (1.0, 1.0, 0.0) 1 1 12),
buildSphere 1.0 (-2, 2, -4) (Material (1.0, 0.0, 0.0) 1 1 12),
buildSphere 0.1 (-0.5, 0.7, -4.5) (Material (1.0, 0.0, 0.0) 1 1 12)]

V - Anti-aliasing : supersampling stochastique

L'aliasing est un phénomène inhérent aux calculs de nature discrète. Il regroupe les deux phénomènes communs suivants :

- Le phénomène de crénelage, qui se caractérise par l'apparition d'un motif en escalier le long du bord des objets. Il est fortement visible sur toutes les images que nous avons générées. Vous trouverez d'autres exemples sur la page wikipedia.

- Un second phénomène très courant, que nous n'avons pas encore rencontré (faute de disposer d'assez d'objets), mais qui se produit trop souvent, par exemple dès qu'il est question d'appliquer une texture sur un plan, est l'apparition de motifs d'interférences. Cela est dû au fait que l'on superpose notre grille de pixels (que l'on peut considérer comme un réseau) sur un réseau formé par des objets / un motif qui se répète. Des exemples sont disponibles sur cette page ainsi que sur wikipedia.

Il existe plusieurs façons de réduire l'effet d'aliasing. L'une des plus simples consiste à lancer 4 rayons sur les sommets du pixel (vue comme un carré) plutôt qu'un unique rayon central, et de faire la moyenne. Cela revient à générer une image de Largeur + 1 par Hauteur +1 pixels. On est donc dans la famille du sur-échantillonnage, ou encore supersampling.

La solution que je vous propose est du même genre, mais légèrement plus raffinée (et malheureusement, coûteuse en termes de performances). Nous allons subdiviser chaque pixel en sous-pixels, où n est le taux de sur-échantillonnage. Pour chacun de ces sous-pixels, on lancera un rayon dont la position est aléatoire, mais reste à l'intérieur de ce sous-pixel. De cette façon, on construit une "grille irrégulière", ce qui permet de réduire les interférences. Ensuite, pour obtenir la couleur du pixel, on moyenne tous les sous-pixels qui le composent.

Plutôt que de moyenner tous les sous-pixels de façon aveugle, on pourrait les pondérer selon une courbe gaussienne, ou encore courbe en cloche.

On commence par une fonction utilitaire, qui va nous permettre de subdiviser une liste infinie de nombres aléatoires en sous-listes de taille fixe.

Tools.hs

--Coupe une liste en une liste de morceaux de longueur length
spliter :: Int -> [a] -> [[a]]
spliter n [] = []
spliter n list = first : (spliter second)
where
(first, second) = splitAt n list

On ajoute aussi la possibilité de produire une liste de valeurs aléatoires qui restent identique à chaque exécution (on ne souhaite pas que deux exécutions de notre programme avec la même scène nous donne deux valeurs différentes).

Tools.hs

import System.Random

-- ...

--Produit des valeurs entre -offset et +offset
randomOffsets :: RealRep -> [RealRep]
randomOffsets offset = randomRs (0, offset) g
where
g = mkStdGen 42

On va vouloir manipuler des pixels aux coordonnées non entières. Puisqu'au final, nos pixels ne sont que des coordonnées sur une grille, pourquoi ne pas les stocker dans des RealReps ?

Tools.hs

--Une coordonnée dans l'image
type PixelLocation = (RealRep, RealRep)

Les conversions de computeRay deviennent alors inutiles et doivent être modifiées.

Raytracer.hs

--Construit un rayon à partir d'un pixel
computeRay :: PixelLocation -> Scene -> Ray
computeRay (x, y) (Scene camera objects _) = normalise (vx, vy, vz)
where Camera _ distance (planWidth, planHeight) _ = camera
halfPlanWidth = int2RealRep (planWidth `div` 2)
halfPlanHeight = int2RealRep (planHeight `div` 2)
vx = x - halfPlanWidth
vy = y - halfPlanHeight
vz = -distance

Nous allons modifier la fonction getPixel pour qu'elle prenne un paramètre entier supplémentaire : le taux de sur-échantillonnage. La valeur 1 correspondra à lancer un seul rayon, alors que n signifiera n*n rayons lancés. Dans cette première version, c'est un sur-échantillonnage où le rayon sera toujours lancé dans l'angle inférieur gauche (x minimal, y minimal) du sous pixel.

Raytracer.hs

type SamplingRate = Int
getPixel :: SamplingRate -> Scene -> PixelLocation -> [Word8]
getPixel samplingRate scene pixelLocation@(px, py) = colorToPixel color
where
--Conversion en nombre à virgule
samplingRate' = int2RealRep samplingRate
--Produit l'intersection
pixelIntersection pixelLocation = return pixelLocation
>>= computeRay
>>= computeIntersections
>>= getClosestIntersection
--Construit la couleur du pixel
Scene _ _ lights = scene
intersectionToPixel intersection = computeColor intersection $ scene
maybePixelColor pixelLocation = fmap intersectionToPixel (pixelIntersection pixelLocation scene)
pixelColor maybeColor = case maybeColor of
Nothing -> (0, 0, 0)
Just color -> color
--Calcule l'écart entre deux sous-pixels
offset = 1 / samplingRate'
--Construit une liste de couleurs à partir des sous-pixels
colorList = map (pixelColor.maybePixelColor) (pixels samplingRate)
--Construit la liste des sous-pixels
pixels 1 = [pixelLocation]
pixels samplingRate = [(px + x, py + y) | x <- subCells, y <- subCells]
where
subCells = map (\x -> (int2RealRep x) * offset - 1/2) [0..samplingRate-1]
--Fait la moyenne de toutes les sources lumineuses
color = (offset*offset) *! (vectorSum colorList)

Vous pouvez maintenant tester cette anti-aliasing avec le main (main si ça veut le coup).

Pour rajouter un comportement aléatoire, nous allons légèrement modifier la fonction getPixel pour qu'elle puisse prendre en paramètre une liste d'offsets à appliquer à chaque sous-pixel (elle se voit donc renommée).

Raytracer.hs

type SamplingRate = Int
getPixel :: SamplingRate -> Scene -> PixelLocation -> [Word8]
getPixel sr sc px = stochasticGetPixel sr sc px []

type RandList = [RealRep]
stochasticGetPixel :: SamplingRate -> Scene -> PixelLocation -> RandList -> [Word8]
stochasticGetPixel samplingRate scene pixelLocation@(px, py) randList = colorToPixel color
where
-- ...
--Construit la liste des sous pixels
pixels 1 = [pixelLocation]
pixels samplingRate = [(px + x, py + y) | x <- subCells, y <- subCells]
where
subCells = map (\x -> (int2RealRep x) * offset - 1/2) [0..samplingRate-1]
stocasticCells = if null randList then subCells else newSubCells
where
newSubCells = zipWith (+) subCells randList
-- ...

On ajoute quelques fonctions utilitaires pour appeler plus simplement stochasticGetPixel :

Raytracer.hs

getPixelMap :: SamplingRate -> Scene -> [PixelLocation] -> [Word8]
getPixelMap samplingRate scene locations = pixels
where
stochasticPixels = map
(stochasticGetPixel samplingRate scene)
locations
pixels = concat $ zipWith ($) stochasticPixels (generateRandLists samplingRate)

generateRandLists :: SamplingRate -> [RandList]
generateRandLists samplingRate = spliter (samplingRate*samplingRate) $ randomOffsets (1 / samplingRate')
where
samplingRate' = int2RealRep samplingRate

On remplace maintenant l'appel dans le main :

Main.hs

tableData width height = listArray ((0,0,0), (height-1, width-1, 3)) $ pixels
where
pixels = getPixelMap 2 scene [(int2RealRep x, int2RealRep y)
| y <- [0..height - 1], x <- [0..width - 1]]

Et c'est terminé! On se quitte sur une scène illustrant tout ce que nous avons ajouté au cours de cet article.

VI - Remerciements

Un grand merci à Cédric Duprez pour sa relecture orthographique.