Ce document traite le sujet de centrale d’informatique de 2021.

import math, numpy as np

from typing import Tuple

vecteur = point = np.ndarray
rayon = Tuple[point, vecteur]
sphère = Tuple[point, float]

noir = np.array([0., 0., 0.])
blanc = np.array([1., 1., 1.])
couleur = image = np.ndarray

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

I Géométrie

Q 1.

def vec(A: point, B:point) -> vecteur:
    return B - A

Q 2.

def ps(v1:vecteur, v2:vecteur) -> float:
    return np.sum(v1 * v2)

Q 3.

def norme(v:vecteur) -> float:
    return np.sqrt(ps(v, v))

Q 4.

def unitaire(v:vecteur) -> vecteur:
    return (1. / norme(v)) * v

Q 5.

def pt(r:rayon, t:float) -> point:
    """Calcule les coordonées du point M du rayon r=(S, u) tel que SM = t*u"""
    assert t >= 0
    (S, u) = r
    return S + t * u

def dir(A:point, B:point) -> vecteur:
    """Construit le vecteur unitaire associé à la direction du vecteur AB"""
    return unitaire(vec(A, B))

def ra(A:point, B:point) -> rayon:
    """Construit le rayon lumineux issu du point A et de direction AB"""
    return A, dir(A, B)

Q 6.

def sp(A:point, B:point) -> sphère:
    return A, norme(vec(A, B))

Q 7.

Soit $M \in (A, \vec{u})$ alors $\exists t \in \mathbb{R}^+$ tel que $\vec{AM} = t \vec{u}$

Alors

$\begin{align} \vec{AC} + \vec{CM} = t \vec{u} & \iff (\vec{AC} + \vec{CM})^2 = t^2 \vec{u}^2 = t^2\\ & \iff \lVert\vec{CA}\rVert^2 + 2 \vec{AC} \cdot \vec{CM} + \lVert\vec{CM}\rVert^2 = t^2\\ & \iff \lVert\vec{CA}\rVert^2 + 2 \vec{AC} \cdot (\vec{CA} + \vec{AM}) + \lVert\vec{CM}\rVert^2 = t^2\\ & \iff \lVert\vec{CA}\rVert^2 - 2 \lVert\vec{CA}\rVert^2 + 2t (\vec{u}\cdot\vec{AC}) + \lVert\vec{CM}\rVert^2 = t^2\\ & \iff -\lVert\vec{CA}\rVert^2 - 2t (\vec{u}\cdot\vec{CA}) + \lVert\vec{CM}\rVert^2 = t^2\\ & \iff t^2 + 2t (\vec{u}\cdot\vec{CA}) + \lVert\vec{CA}\rVert^2 - \lVert\vec{CM}\rVert^2 = 0 \end{align}$

i.e. $\forall M \in (A, \vec{u}) \,\, \exists t \in \mathbb{R}^+$ tel que $t^2 + 2t (\vec{u}\cdot\vec{CA}) + \lVert\vec{CA}\rVert^2 - \lVert\vec{CM}\rVert^2 = 0$ (1)

1) Supposons que $(A, \vec{u})$ et $(C, r)$ sont séquentes

L’intersection entre un cercle et une droite ne peut résulter qu’en un ou deux point que l’on note $M$ ou $M_1$ et $M_2$.

$M \in (C, r)$ implique que $\lVert\vec{CM}\rVert^2 = r^2$

donc d’après (1) on obtient : $t^2 + 2t (\vec{u}\cdot\vec{CA}) + \lVert\vec{CA}\rVert^2 - r^2 = 0$ et cette équation n’admet de solutions qu’en un ou deux points i.e. n’admet qu’une ou deux solutions réelles.

2) Supposons que $t^2 + 2t (\vec{u}\cdot\vec{CA}) + \lVert\vec{CA}\rVert^2 - r^2 = 0$ possède deux solutions réelles

Donc $\Delta = 4 (\vec{u}\cdot\vec{CA})^2 - 4 (\lVert\vec{CA}\rVert^2 - r^2) \ge 0$

et $t_1 = - \vec{u}\cdot\vec{CA} + \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)}$

et $t_2 = - \vec{u}\cdot\vec{CA} - \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)}$

sont les expressions des deux solutions de l’équation.

Soit $M_1 \in (A, \vec{u})$ tel que $\vec{AM_1} = t_1 \vec{u}$ et $M_2 \in (A, \vec{u})$ tel que $\vec{AM_2} = t_2 \vec{u}$

On a $\vec{AM_1} = \vec{AC} + \vec{CM_1} = t_1 \vec{u} = - (\vec{u}\cdot\vec{CA})\,\, \vec{u} + \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)}\,\, \vec{u}$

i.e. $\vec{CM_1} = - (\vec{u}\cdot\vec{CA})\,\, \vec{u} + \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)}\,\, \vec{u} + \vec{CA}$

donc $\vec{CM_1}\cdot\vec{u} = - \vec{u}\cdot\vec{CA} + \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)} + \vec{CA}\cdot\vec{u} = \sqrt{(\vec{u}\cdot\vec{CA})^2 - (\lVert\vec{CA}\rVert^2 - r^2)}$

Soit $\vec{v}$ le vecteur unitaire orthogonal à $\vec{u}$ obtenu en lui applicant une rotation de $\frac{\pi}{2}$

alors $\vec{u}\cdot\vec{v} = 0$ et $\vec{u},\vec{v} = \frac{\pi}{2} \iff \vec{u},\vec{CA} + \vec{CA},\vec{v} = \frac{\pi}{2} \iff \vec{CA},\vec{v} = \frac{\pi}{2} - \vec{u},\vec{CA}$

donc $\vec{CM_1}\cdot\vec{v} = \vec{CA}\cdot\vec{v} = \lVert\vec{CA}\rVert \cos\left(\vec{CA},\vec{v}\right) = \lVert\vec{CA}\rVert \cos\left(\frac{\pi}{2} - \vec{u},\vec{CA}\right) = \lVert\vec{CA}\rVert \sin\left(\vec{u},\vec{CA}\right)$

ainsi

$\begin{align} \lVert\vec{CM_1}\rVert^2 &= (\vec{CM_1}\cdot\vec{u})^2 + (\vec{CM_1}\cdot\vec{v})^2\\ &= (\vec{u}\cdot\vec{CA})^2 - \lVert\vec{CA}\rVert^2 + r^2 + \lVert\vec{CA}\rVert^2 \sin^2\left(\vec{u},\vec{CA}\right)\\ &= (\vec{u}\cdot\vec{CA})^2 - \lVert\vec{CA}\rVert^2 + r^2 + \lVert\vec{CA}\rVert^2 \left(1 - \cos^2\left(\vec{u},\vec{CA}\right)\right)\\ &= (\vec{u}\cdot\vec{CA})^2 - \lVert\vec{CA}\rVert^2 + r^2 + \lVert\vec{CA}\rVert^2 - (\vec{u}\cdot\vec{CA})^2\\ &= r^2 \end{align}$

ainsi $M_1 \in (C, r)$, de même pour $M_2$ donc $(A, \vec{u})$ et $(C, r)$ sont séquentes en un ou deux points respectivement $M$ ou $M_1$ et $M_2$

Q 8.

def intersection(r:rayon, s:sphère) -> (point, float):
    A, u = r
    C, R = s
    
    d = ps(u, vec(C, A))**2 - (norme(vec(C, A))**2 - R**2)
    
    if (d < 0) or (norme(vec(C, A)) < R) or (ps(u, vec(C, A)) > 0):
        # je le savais !!! il manque à l'énoncé la réciproque que c'est des solutions réelles positives pour que les solutions correspondent bien à une intersection du bon côté du rayon c.a.d côté direction de u
        # donc il n'y a même pas besoin de cette longue preuve, la plus courte ou CM = .... = r marche !!
        return np.array([-np.inf, -np.inf, -np.inf]), -np.inf
    
    t = - ps(u, vec(C, A)) - np.sqrt(d)
    
    return pt(r, t), t

Vérification

A = np.array([0., 0.5, 0.])
C = np.array([3., 0., 0.])

u = unitaire(vec(A, C + np.array([0., -0.5, 0.])))
r = 1.

M, t = intersection((A, u), (C, r))

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_aspect(1)
ax.set_xlim(-1, 5)
ax.set_ylim(-2, 2)

ax.add_artist(plt.Circle(C[:-1], r, fc='white', ec='green'))
ax.add_artist(plt.Line2D(*np.transpose([A, pt((A, u), 100)])[:-1]))
if (t >= 0):
    ax.add_artist(plt.Circle(M[:-1], 0.05, fc='red', ec=None))

png

II Optique

II.A - Visibilité

Q 9.

Les directions des vecteurs $\vec{SP}$ et $\vec{CP}$ doivent être opposés i.e. $\vec{SP}\cdot\vec{CP} < 0$

Q 10.

def au_dessus(s:sphère, P:point, src:point) -> bool:
    C = s[0]
    S = src
    return ps(dir(S, P), dir(C, P)) < 0

Q 11.

def visible(obj:[sphère], j:int, P:point, src:point) -> bool:
    if not au_dessus(obj[j], P, src):
        return False
    
    for i, s in enumerate(obj):
        if i != j:
            if intersection(ra(P, src), s)[1] >= 0:
                return False
    
    return True

II.B - Diffusion

Q 12.

def couleur_diffusée(r:rayon, Cs:couleur, N:vecteur, kd:couleur) -> couleur:
    A, u = r
    return kd * Cs * ps(N, -u)

II.C - Réflexion

Q 13.

def pv(u:vecteur, v:vecteur) -> vecteur:
    """Produit vectoriel de u et v"""
    return np.array([
        u[1]*v[2] - u[2]*v[1],
        u[2]*v[0] - u[0]*v[2],
        u[0]*v[1] - u[1]*v[0]
    ])

def rotation(u:vecteur, N:vecteur, phi:float) -> vecteur:
    """Formule de rotation d'Euler - Rodrigues"""
    return np.cos(phi) * u + (1 - np.cos(phi)) * ps(u, N) * N + np.sin(phi) * pv(N, u)

def rayon_réfléchi(s:sphère, P:point, src:point) -> rayon:
    u = unitaire(vec(P, src))
    C, r = s
    N = unitaire(vec(C, P))
    return P, rotation(u, N, np.pi)

Vérification

A = np.array([0., 0.5, 0.])
C = np.array([3., 0., 0.])

u = unitaire(vec(A, C + np.array([0., -0.5, 0.])))
r = 1.

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_aspect(1)
ax.set_xlim(-1, 5)
ax.set_ylim(-2, 2)

# Cercle
ax.add_artist(plt.Circle(C[:-1], r, fc='white', ec='green'))

# Rayon
ax.plot(*np.transpose([A, pt((A, u), 100)])[:-1])
ax.arrow(*A[:2], *u[:2], head_width=.05, head_length=.07)

M, t = intersection((A, u), (C, r))

r = rayon_réfléchi((C, r), M, A)
A, u = r
ax.plot(*np.transpose([A, pt((A, u), 100)])[:-1])
ax.arrow(*A[:2], *u[:2], head_width=.05, head_length=.07)

N = unitaire(vec(C, M))
ax.arrow(*M[:2], *N[:2], head_width=.05, head_length=.07)

<matplotlib.patches.FancyArrow at 0x7ff85c5cfb10>

png

IV Lancer de rayons

IV.A - Écran

Q 18.

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

Q 19.

def rayon_écran(omega:point, i:int, j:int) -> rayon:
    return (omega, unitaire(vec(omega, grille(i, j))))

IV.B - Couleur d’un pixel

Q 20.

def interception(r:rayon) -> (point, int):
    tmin = np.inf
    imin = -1
    pmin = np.array([np.inf, np.inf, np.inf])
    
    for i, s in enumerate(Objet):
        M, t = intersection(r, s)
        if 0 <= t < tmin:
            tmin = t
            imin = i
            pmin = M
    
    return pmin, imin

Q 21.

def couleur_diffusion(P:point, j:int) -> couleur:
    
    c = noir
    
    for i, src in enumerate(Source):
        if visible(Objet, j, P, src):
            r = (src, unitaire(vec(src, P)))
            
            C = Objet[j][0]
            N = unitaire(vec(C, P))
            
            Cs = ColSrc[i]
            kd = KdObj[j]
            
            c = c + couleur_diffusée(r, Cs, N, kd)
    
    return c

IV.C - Constitution de l’image

Q 22.

def lancer(omega:point, fond:couleur) -> image:
    image = np.zeros((N, N, 3))
    
    for i in range(N):
        for j in range(N):
            
            P, k = interception(rayon_écran(omega, i, j))
            
            if k >= 0:
                image[i][j] = couleur_diffusion(P, k)
            else:
                image[i][j] = fond
            
    return image

IV.D - Compléxités

Q 23. Q 24.

Compléxité de interception : $\mathcal{O}(n_o)$

Compléxité de visible : $\mathcal{O}(n_o)$

Compléxité de couleur_diffusion : $\mathcal{O}(n_on_s)$

Compléxité de lancer : $\mathcal{O}(N^2n_on_s)$

Vérification

%%time

Objet = [
    (np.array([2., 0., -7]), 3),
    (np.array([-.7, 1.5, -2]), 1)
]
KdObj = [
    np.array([1.0, 0.0, 0.0]),
    np.array([0.0, 0.0, 1.0])
]
Source = [
    np.array([-5, 5, 1]),
]
ColSrc = [
    blanc,
]

Delta = 10
N = 500

omega = np.array([0., 0., 20])
gris_bleuté = np.array([0.8, 0.8, 0.8])

img = lancer(omega, gris_bleuté)

fig, ax = plt.subplots(figsize=(6, 6), dpi=120)
ax.imshow(img)
_ = plt.axis('off')

CPU times: user 24.7 s, sys: 41 ms, total: 24.8 s
Wall time: 24.8 s

png

%%time

Objet = [
    (np.array([0., 0., -4]), 3),
]
KdObj = [
    np.array([1.0, 0.0, 0.0]),
]
Source = [
    np.array([-5, 5, -1]),
    np.array([5, 5, -1]),
]
ColSrc = [
    blanc,
    blanc,
]

Delta = 10
N = 500

omega = np.array([0., 0., 20])
gris_bleuté = np.array([0.8, 0.8, 0.8])

img = lancer(omega, gris_bleuté)

fig, ax = plt.subplots(figsize=(6, 6), dpi=120)
ax.imshow(img)
_ = plt.axis('off')

CPU times: user 16.3 s, sys: 40.4 ms, total: 16.3 s
Wall time: 16.3 s

png

V Améliorations

V.A - Prise en compte de la réflexion

Q 25.

def réflexions(r:rayon, rmax:int) -> [(point, int)]:
    pts = []
    
    for k in range(rmax):
        P, i = interception(r)
        if i >= 0:
            r = rayon_réfléchi(Objet[i], P, r[0])
            pts.append((P, i))
    
    return pts

Vérification

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

O = np.array([0, 0, 0])
x = np.array([1, 0, 0])
y = np.array([0, 1, 0])
z = np.array([0, 0, 1])

L = np.array([-5, 5, 1])

# base 3d

ax.scatter(*O, color='y')
ax.quiver(*O, *x, color='r')
ax.quiver(*O, *y, color='b')
ax.quiver(*O, *z, color='g')


# source
ax.scatter(*L, color='y')


# rayon
def dessine_rayon(r:rayon, ax):
    A, u = r
    ax.quiver(*A, *u)
    ax.scatter(*A, color='g')
    ax.plot3D(*np.transpose([A, pt((A, u), 100)]), lw=0.5)
    
C = np.array([ 2.,  0., -7.])
r = 3.0

C1 = np.array([-1.3, 1.5, -2])
r1 = 1.0

u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
x = C[0] + r * np.cos(u)*np.sin(v)
y = C[1] + r * np.sin(u)*np.sin(v)
z = C[2] + r * np.cos(v)
ax.plot_wireframe(x, y, z, color="r", alpha=.2)

x = C1[0] + r1 * np.cos(u)*np.sin(v)
y = C1[1] + r1 * np.sin(u)*np.sin(v)
z = C1[2] + r1 * np.cos(v)
ax.plot_wireframe(x, y, z, color="b", alpha=.2)

ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([-10, 1])

def réflexions_test(r:rayon, rmax:int) -> [(point, int)]:
    pts = []
    
    for k in range(rmax):
        P, i = interception(r)
        if i >= 0:
            r = rayon_réfléchi(Objet[i], P, r[0])
            pts.append((r, i))
    
    return pts

screenray = rayon_écran(omega, 240, 290)
refray = réflexions_test(screenray, 10)
dessine_rayon(screenray, ax)
for ray in refray:
    dessine_rayon(ray[0], ax)

ax.view_init(20, 45)

png

Q 26.

def couleur_perçue(r:rayon, rmax:int, fond:couleur) -> couleur:
    
    refx = réflexions(r, rmax)
    
    Ck = noir
    
    for P, k in reversed(refx):
        if k >= 0:
            Cdk = couleur_diffusion(P, k)
            Ck = Cdk + KrObj[k] * Ck
    
    return Ck

Q 27.

def lancer_complet(omega:point, fond:couleur, rmax:int) -> image:
    image = np.zeros((N, N , 3))
    
    for i in range(N):
        for j in range(N):
            r = rayon_écran(omega, i, j)
            P, k = interception(r)
            if k < 0:
                image[i][j] = fond
            else:
                image[i][j] = couleur_perçue(r, rmax, fond)
    
    return image

Vérification

%%time

Objet = [
    (np.array([2., 0., -6.5]), 3),
    (np.array([-2, 1.5, -1]), 2)
]

KdObj = [
    np.array([0.8, 0.8, 0.8]),
    np.array([0.0, 0.0, 1.0])
]

KrObj = [
    1.0,
    0.1
]

Source = [
    np.array([4, 4, 1]),
]

ColSrc = [
    blanc,
    blanc,
    blanc
]

Delta = 10
N = 500

omega = np.array([0., 0., 20])
gris_bleuté = np.array([0.8, 0.8, 0.8])

img = lancer_complet(omega, gris_bleuté, 3)

fig, ax = plt.subplots(figsize=(6, 6), dpi=120)
ax.imshow(img)
_ = plt.axis('off')

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).


CPU times: user 42.5 s, sys: 76.1 ms, total: 42.6 s
Wall time: 42.5 s

png

Sujet d'informatique MP Centrale 2021

I Géométrie

Q 1.

Q 2.

Q 3.

Q 4.

Q 5.

Q 6.

Q 7.

Q 8.

II Optique

II.A - Visibilité

Q 9.

Q 10.

Q 11.

II.B - Diffusion

Q 12.

II.C - Réflexion

Q 13.

IV Lancer de rayons

IV.A - Écran

Q 18.

Q 19.

IV.B - Couleur d’un pixel

Q 20.

Q 21.

IV.C - Constitution de l’image

Q 22.

IV.D - Compléxités

Q 23. Q 24.

V Améliorations

V.A - Prise en compte de la réflexion

Q 25.

Q 26.

Q 27.