Lancer de rayons (sujet Centrale 2021)

Sujet : https://www.concours-centrale-supelec.fr/CentraleSupelec/2021/Multi/sujets/I007.pdf

import numpy as np
import plotly
import matplotlib.pyplot as plt

plotly.offline.init_notebook_mode(connected=True)  # pour export html

def point(x, y, z):
    return np.array([x, y, z])

I Géométrie

# Q1
def vec(A, B):
    return np.array(B) - np.array(A)

# Q2 
def ps(v1, v2):
    return np.sum(v1 * v2)

# Q3
def norme(v):
    return np.sqrt(np.sum(v**2))

# Q4
def unitaire(v):
    return v / norme(v)

# Q5
def pt(r, t):
    S, u = r
    return S + t * u

def dir(A, B):
    return unitaire(vec(A, B))

def ra(A, B):
    return A, dir(A, B)

# Q6
def sp(A, B):
    return (A, norme(B - A))

# Q8
def intersection(r, s):
    S, u = r
    C, r_s = s
    CS = vec(S, C)
    b = 2*ps(u, CS)
    delta = b**2 - 4*(norme(CS)**2 - r_s**2)
    if delta < 0:
        return None
    t1, t2 = -b - delta**0.5, -b + delta**0.5
    t = min(abs(t1), abs(t2))/2
    return S + t*u, t

II Optique

Q 9. P est le premier point d'intersection d'un rayon partant de S et allant en P. (ou : Produit scalaire positif avec la normale au plan tangent)

# Q10
def au_dessus(s, P, src):
    return ps(vec(P, src), vec(s[0], P)) > 0

# Q11
def visible(obj, j, P, src):
    for i, o in enumerate(obj):
        inter = intersection(ra(src, P), o)
        if i != j and inter != None and inter[1] < norme(vec(src, P)):
            return False
    return au_dessus(obj[j], P, src)

# Q12
def couleur_diffusée(r, Cs, N, kd):
    cos = ps(-r[1], N)
    return cos*(Cs*kd)

# Q13
def rayon_réfléchi(s, P, src):
    C, r = s
    u = dir(src, P)
    N = dir(C, P)
    return P, 2*N + u

Digression : affichage d'un rayon avec Plotly

import plotly.graph_objects as go

def line(P, Q):
    return go.Scatter3d(x=[P[0], Q[0]], y=[P[1], Q[1]], z=[P[2], Q[2]])

theta, phi = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j]
s = (point(0, 0, 0), 1)
src = point(3, 3, 2)
P = point(1, 0, 0)
_, u = rayon_réfléchi(s, P, src)

fig = go.Figure(data=[ 
    line(src, P),
    go.Surface(x=np.cos(theta)*np.cos(phi), y=np.sin(theta)*np.cos(phi), z=np.sin(phi)),
    line(P, P + 3*u),
])
fig.show()

img_bytes = fig.to_image(format="png")
from IPython.display import Image
Image(img_bytes)  # image statique

Lancer de rayons

# Q18
def grille(i, j):
    return np.array([(Delta/N)*(j - N/2), (Delta/N)*(N/2 - i), 0])

# Q19
def rayon_écran(omega, i, j):
    return omega, dir(omega, grille(i, j))

# Q20
def interception(r, eviter=-1):
    x, minim = None, np.inf
    for i in range(len(Objet)):
        if i != eviter:
            t = intersection(r, Objet[i])
            if t != None and t[1] < minim:
                x, minim = (t[0], i), t[1]
    return x

# Q21
def couleur_diffusion(P, j):
    coul = np.array([.0, .0, .0])
    for i in range(len(Source)):
        if visible(Objet, j, P, Source[i]):
            coul += couleur_diffusée(ra(Source[i], Objet[j][0]), ColSrc[i], dir(Objet[j][0], P), KdObj[j])
    return coul

# Q22
def lancer(omega, fond):
    img = np.full((N, N, 3), fond, dtype=np.float64)
    for i in range(N):
        for j in range(N):
            inter = interception(rayon_écran(omega, i, j))
            if inter is not None:
                P, k = inter
                img[i, j] = couleur_diffusion(P, k)
    return img

V Améliorations

# Q25
def réflexions(r, rmax):
    res = []
    I = interception(r)
    while len(res) < rmax + 1 and I != None:
        res.append(I)
        r = rayon_réfléchi(Objet[I[1]], I[0], r[0])
        I = interception(r, I[1])
    return res

# Q26
def couleur_perçue(r, rmax, fond):
    C = np.array([.0, .0, .0])  # noir
    refl = réflexions(r, rmax)[::-1]
    if not refl:
        return fond
    for P, i in refl:
        C = couleur_diffusion(P, i) + KrObj[i]*C
    return C

# Q27
def lancer_complet(omega, fond, rmax):
    img = np.full((N, N, 3), fond, dtype=np.float64)
    for i in range(N):
        for j in range(N):
            img[i, j] = couleur_perçue(rayon_écran(omega, i, j), rmax, fond)
    return np.clip(img, .0, 1.)

Objet = [
    (point(-2.5, 0, -3), 1),
    (point(2.5, 0, -3), 1.5),
    (point(0, 2, -2), 1.5),
    (point(.5, -1.5, -1), 1.2),
]
KdObj = [
    point(.0, 1., .0),
    point(.2, .7, 1.),
    point(1., .0, .0),
    point(1.0, 1.0, 1.0),
]    

Source = [point(-2, 0, 3), point(3, 0, 3)]
ColSrc = [point(1, 1, 1), point(0, 0, 1)]

N = 1000
Delta = 5
KrObj = [.5, .5, .5, .5]
im = lancer_complet(point(0, 0, 4), point(.7, .7, .7), 2)
plt.figure(figsize = (20,10))
plt.imshow(im)
plt.imsave("ray.png", im)