Nous allons décomposer l'environnement en trois ensembles:
Pour réaliser une image contenant le faisceau lumineux et les ombres portées des objets nous allons procéder en suivant les étapes suivantes:
On va supposer que le faisceau lumineux a une étendue d'au plus 180° et que les objets placés dans l'environnement n'ont pas des formes "trop tordues" (notion que l'on précisera plus tard :) ).
Le code utilise les librairies numpy et matplotlib ainsi que le fichier geom_utils.py défini plus haut et contenant des fonctions utiles pour la géométrie.
In [1]:
import numpy as np
from matplotlib.patches import PathPatch
from matplotlib.path import Path
from matplotlib import pyplot as plt
from matplotlib import animation
from geom_utils import *
Ray et Beam:On commence par définir une classe Ray représentant un rayon lumineux issu d'une source et s'arrêtant en un point limite. Ce rayon est caractérisé par un l'angle qu'il fait avec:
u représentant l'origine des angles (typiquement le vecteur coupant le faisceau lumineux en deux faisceaux de même taille),x si aucun vecteur origine u n'est précisé.Dans la suite, on comparera deux rayons lumineux en considérant qu'un rayon r est supérieur à un rayon s si l'angle de r est plus grand que celui de s.
In [2]:
class Ray(object):
"""A class representing a ray of light which is nothing but a segment
characterized by a start, an end and an angle computed with regards to
an origin vector `u`."""
def __init__(self, source, limit, u=None):
self._start = source
self._end = limit
ray = limit - source
if u is None:
self._angle = vector_to_angle(ray)
else:
origin = vector_to_angle(u)
rot = rotation(-origin, two_dim=True).dot(ray)
self._angle = vector_to_angle(rot)
@property
def angle(self):
return self._angle
@property
def end(self):
return self._end
def __cmp__(self, ray):
return -1 if self.angle < ray.angle else \
0 if self.angle == ray.angle else 1
Ensuite on définit une classe Beam qui est représente les limites d'un faisceau lumineux. Elle est caractérisée par une source de lumière et des vecteurs limites aux extrémités du faisceau:
In [3]:
class Beam(object):
"""A class representing a beam of light.
Attributes:
source : a numpy 1D array with two elements representing
the origin of the beam.
start : a numpy 1D array with two elements representing
the limit vector where the beam starts.
end : a numpy 1D array with two elements representing
the limit vector where the beam ends."""
def __init__(self, source, start, end):
self._source = source
self._start = start
self._end = end
@property
def start(self):
return self._start
@property
def end(self):
return self._end
@property
def source(self):
return self._source
In [4]:
# A rect with bottom-left vertex (0,0) and top-right (10, 10)
WORLD = rectangle_vertices(0, 0, 10, 10)
world_patch = get_polygon_patch(WORLD, facecolor='black', lw=3)
Pour colorier le faisceau, il faut déterminer le polygone formé par le faisceau et les bords du monde virtuel (le faisceau ne peut pas déborder au delà des limites du monde). On utilise pour cela la fonction world_beam_intersection et world_intersection. Cette fonction utilise la fonction lines_intersection de geom_utils.py pour trouver l'intersection entre deux droites données chacune par un point et un vecteur directeur.
Le principe est le suivant:
vertex_1 == vertex_2), et dans ce cas le polygone recherché est le triangle formé par les deux rayons limites du faisceau et un bout du côté du monde virtuel.
In [5]:
def world_beam_intersection(beam):
"""Returns a list of polygon vertices representing the limits of the beam of
light when it is stopped by the `WORLD` boundaries."""
source, start, end = beam.source, beam.start, beam.end
intersection_1, vertex_1 = world_intersection(source, start)
intersection_2, vertex_2 = world_intersection(source, end)
l = len(WORLD)
if vertex_1 == vertex_2:
v = WORLD[(vertex_1+1)%l] - source
if v.dot(start) < v.dot(end):
return [source, intersection_1, intersection_2]
else:
world = [WORLD[(vertex_1+1+i)%l] for i in xrange(l)]
return [source, intersection_1] + world + [intersection_2]
else:
return [source, intersection_1] + \
[WORLD[(vertex_1+1+i)%l] for i in xrange((vertex_2-vertex_1)%l)] + \
[intersection_2]
def world_intersection(point, vector):
"""Returns the point which is the intersection of the line starting
from `point` with direction `vector` and the `WORLD`, as well as the
index of the `WORLD` vertex of the intersection."""
world = WORLD + [WORLD[0]]
intersections = [(lines_intersection(point, vector,
world[i], world[i+1] - world[i]),
i)
for i in xrange(len(WORLD))]
# filter only the intersections in the direction of vector
intersections = [(p, i) for p, i in intersections
if p is not None
if (p - point).dot(vector) >= 0]
distances = [(p - point).dot(p - point) for p, i in intersections]
# a trick to get the index of the minimum:
min_distance_index = min(xrange(len(distances)),
key=distances.__getitem__)
return intersections[min_distance_index]
SOURCE = np.r_[5,5]
START = np.r_[1., 1.5]
END = np.r_[-2, -0.5]
START = START/(START.dot(START))**0.5
END = END/(END.dot(END))**0.5
BEAM = Beam(SOURCE, START, END)
light_patch = get_polygon_patch(world_beam_intersection(BEAM),
facecolor='yellow', alpha=0.7)
shade_of_obstacles est la fonction la plus grosse de notre programme. Elle cherche à déterminer un polygone représentant l'ombre portée d'un obstacle éclairé par un faisceau.
On procède de la façon suivante:
In [6]:
def shade_of_obstacle(obstacle, beam):
"""Returns a list of polygon vertices representing the shade of an `obstacle`
lightened by a `beam`."""
source, start, end = beam.source, beam.start, beam.end
# we place the origin for the angles at the middle of the beam:
origin = start + end
theta_origin = vector_to_angle(origin)
theta_min = vector_to_angle(rotation(-theta_origin, two_dim=True).dot(start))
theta_max = vector_to_angle(rotation(-theta_origin, two_dim=True).dot(end))
# We consider the rays which are intersecting a vertex of the obstacle:
candidates = [Ray(source, vertex, origin) for vertex in obstacle]
# Here we suppose that the range of the beam is maximum pi so we restrict
# our search to -pi/2, +pi/2 and take the rays with min and max angle
if not any([c for c in candidates if -np.pi/2 <= c.angle <= np.pi/2]):
return None
ray_min, ray_max = min(candidates), max(candidates)
# If the obstacle is outside the scope of vision of the source: return None
if (ray_max.angle < theta_min) or (ray_min.angle > theta_max) or \
abs(ray_max.angle - ray_min.angle) > np.pi:
return None
# If the obstacle is partially in the scope of vision of the source:
# adjust by finding a point P inside the obstacle intersecting the beam
if ray_min.angle < theta_min:
P = lines_intersection(source, start, ray_min.end,
ray_max.end - ray_min.end)
if P is not None:
ray_min = Ray(source, P, origin)
if ray_max.angle > theta_max:
P = lines_intersection(source, end, ray_min.end,
ray_max.end - ray_min.end)
if P is not None:
ray_max = Ray(source, P, origin)
return world_shade_intersection(ray_min.end,
angle_to_vector(ray_min.angle + theta_origin),
ray_max.end,
angle_to_vector(ray_max.angle + theta_origin))
def world_shade_intersection(P1, u1, P2, u2):
"""Returns a list of polygon vertices representing the intersection
of the semi-lines `(P1, u1)` and `(P2, u2)` with the `WORLD` boudaries."""
intersection_1, vertex_1 = world_intersection(P1, u1)
intersection_2, vertex_2 = world_intersection(P2, u2)
l = len(WORLD)
if vertex_1 == vertex_2: return [P1, P2, intersection_2, intersection_1]
else:
return [P1, intersection_1] + \
[WORLD[(vertex_1+1+i)%l] for i in xrange((vertex_2-vertex_1)%l)] + \
[intersection_2, P2]
In [7]:
OBSTACLES = [
[np.r_[4.5, 7], np.r_[3, 8], np.r_[1.5, 7], np.r_[3, 6]],
[np.r_[8.5, 7], np.r_[7, 8], np.r_[5.5, 7], np.r_[7, 6]],
[np.r_[4.5, 3], np.r_[5.5, 3], np.r_[5.5, 4], np.r_[4.5, 4]],
[np.r_[3, 1], np.r_[4, 1], np.r_[5, 1.5], np.r_[6, 1],
np.r_[7, 1], np.r_[8, 2], np.r_[2, 2]],
]
shades = [shade_of_obstacle(o, BEAM) for o in OBSTACLES]
shades = [shade for shade in shades if shade is not None]
shades_patches = [get_polygon_patch(shade, facecolor='#171408',
ec='#171408', lw=2)
for shade in shades]
In [8]:
obstacles_patches = [get_polygon_patch(o, facecolor='#5e5bb2')
for o in OBSTACLES]
In [9]:
%pylab inline
fig = plt.figure()
ax = plt.axes(xlim=(0, 10), ylim=(0, 10))
ax.set_axis_off()
plt.axis('equal')
plt.xlim((0, 10))
plt.ylim((0, 10))
plt.grid(False)
# On affiche le résultat de l'étape 1:
ax.add_patch(world_patch)
# Celui de l'étape 2:
ax.add_patch(light_patch)
# Celui de l'étape 3:
for sp in shades_patches:
ax.add_patch(sp)
# Et enfin celui de l'étape 4:
for op in obstacles_patches:
ax.add_patch(op)
Pour faire quelque chose de plus joli, on peut réaliser une animation qui fait tourner le faisceau lumineux... On utilise le module défini ici qui permet d'embarquer directement une animation dans un notebook ipython.
In [10]:
from mpl_animation_html import display_animation
NB_FRAMES = 400
def animate(i):
ax.patches = []
u = angle_to_vector(vector_to_angle(START) + i*2./NB_FRAMES*np.pi)
v = angle_to_vector(vector_to_angle(END) + i*2./NB_FRAMES*np.pi)
world_patch = get_polygon_patch(WORLD, facecolor='#171408')
obstacles_patches = [get_polygon_patch(o, facecolor='#5e5bb2') for o in OBSTACLES]
ax.add_patch(world_patch)
beam = Beam(SOURCE, u, v)
light = get_polygon_patch(world_beam_intersection(beam), facecolor='#ffefa2', lw=0)
ax.add_patch(light)
for o in OBSTACLES:
shade = shade_of_obstacle(o, beam)
if shade:
p = get_polygon_patch(shade, facecolor='#171408', ec='#171408', lw=2)
ax.add_patch(p)
for o in obstacles_patches:
ax.add_patch(o)
anim = animation.FuncAnimation(fig, animate, frames=NB_FRAMES, interval=20, blit=False)
display_animation(anim)
Out[10]: