A point source of light¶

Each lamp's irradiance at every cell it reaches, falling off with the shell and shaded by what stands between, and the shade, its shadow map.

Requires

each cell's Luminous power Phi (W; zero where nothing shines), its Transmittance t (1, 0..1: the share of light a cell lets through), the Reach r of a lamp (cells), the lattice's spacing (m); the lattice a plane or a space.

Provides

the Light at each cell (W/m on a plane: the power over the shell's length; W/m^2 in space), and its Shade (1, 0..1: the light that arrived over the light that would have in the open; 1 where no lamp reaches, nothing there being shadowed).

Law

from each lamp, E = Phi / shell(d) times the product of the transmittances of the cells strictly between the lamp and the cell, on its ray from one to the other (rays.fan within), for every cell within the reach; the lamps' light adds. On a plane the shell is the circle's length, 2 pi d; in space the sphere's area, 4 pi d^2. A cell is lit on its face whatever its own transmittance: a wall is bright and what is behind it dark.

Source

Pharr, Jakob & Humphreys, Physically Based Rendering, 4th ed., 12.1 (point lights) and 4.1 (radiometry: the inverse-square law); Born & Wolf, Principles of Optics, 7th ed., 1.1 (the energy of a spherical wave).

Proof

curve, invariant, limit.

No Trait form: waits on the middle: a lamp's rays are rays.fan's samples, and the light a cell receives a product over `ordered(it.ray, it.step)` with `@ before`, which draft 5 does not write yet; until then a reference.

light.point_source(luminous, transmittance, reach, spacing)

Each lamp's light at every cell it reaches: the power over the shell, shaded by what stands between, and the shade, the light that arrived over the light in the open (a point source lights and shadows in one pass).

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import library
from lusine import dmath as D
from lusine.kernels.light.point_source import point_source

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import library from lusine import dmath as D from lusine.kernels.light.point_source import point_source

The reference¶

The kernel's reference, in numpy, over rays.fan's samples, one ray to each cell within the reach, folded with @ before; its Trait form is not written yet.

In [2]:

library.reference(point_source)

library.reference(point_source)

Out[2]:

the reference, point_source.py

def point_source(luminous, transmittance, reach, spacing=1.0):
    """A POINT SOURCE OF LIGHT: each lamp's irradiance at every cell it
    reaches, falling off with the shell and shaded by what stands between,
    and the shade, its shadow map.

    Requires: each cell's `Luminous` power Phi (W; zero where nothing shines),
    its `Transmittance` t (1, 0..1: the share of light a cell lets through),
    the `Reach` r of a lamp (cells), the lattice's spacing (m); the lattice a
    plane or a space.
    Provides: the `Light` at each cell (W/m on a plane: the power over the
    shell's length; W/m^2 in space), and its `Shade` (1, 0..1: the light that
    arrived over the light that would have in the open; 1 where no lamp
    reaches, nothing there being shadowed).
    Law: from each lamp, E = Phi / shell(d) times the product of the
    transmittances of the cells strictly between the lamp and the cell, on
    its ray from one to the other (rays.fan `within`), for every cell within
    the reach; the
    lamps' light adds. On a plane the shell is the circle's length, 2 pi d;
    in space the sphere's area, 4 pi d^2. A cell is lit on its face whatever
    its own transmittance: a wall is bright and what is behind it dark.
    Source: Pharr, Jakob & Humphreys, Physically Based Rendering, 4th ed.,
    12.1 (point lights) and 4.1 (radiometry: the inverse-square law); Born &
    Wolf, Principles of Optics, 7th ed., 1.1 (the energy of a spherical wave).
    Proof: curve, invariant, limit."""
    Phi = np.asarray(luminous, dtype=np.float64)
    t = np.asarray(transmittance, dtype=np.float64)
    dims = Phi.ndim
    light = np.zeros_like(Phi)
    open_ = np.zeros_like(Phi)
    offsets = within(reach, dims)
    for at in sources(Phi):
        for off, cell, between in _received(Phi.shape, at, t, reach, offsets):
            d = float(np.sqrt(sum(o * o for o in off))) * float(spacing)
            e = Phi[at] / _shell(d, dims)
            light[cell] += e * between
            open_[cell] += e
    shade = np.ones_like(Phi)
    lit = open_ > 0
    shade[lit] = light[lit] / open_[lit]
    return light, shade

Curve: the shell's falloff¶

In the open the light falls as the shell grows: as 1/d on a plane, the power over the circle's length 2 pi d, and as 1/d^2 in space, over the sphere's area 4 pi d^2, the inverse-square law.

In [3]:

n, c = 61, 30
P = np.zeros((n, n)); P[c, c] = 1.0
E, S = point_source(P, np.ones((n, n)), 25)
d = np.arange(1, 26)
plane = np.array([E[c, c + k] for k in d])
assert np.allclose(plane * 2 * D.PI * d, 1.0, rtol=0, atol=1e-12) and np.all(S == 1.0)
Q = np.zeros((25, 25, 25)); Q[12, 12, 12] = 1.0
E3, _ = point_source(Q, np.ones((25, 25, 25)), 11)
d3 = np.arange(1, 12)
space = np.array([E3[12, 12, 12 + k] for k in d3])
assert np.allclose(space * 4 * D.PI * d3 ** 2, 1.0, rtol=0, atol=1e-12)
fig, ax = plt.subplots(figsize=(7, 3.2))
ax.loglog(d, plane, 'o', label='kernel, a plane')
ax.loglog(d, 1 / (2 * D.PI * d), label='Phi / (2 pi d)')
ax.loglog(d3, space, 's', label='kernel, a space')
ax.loglog(d3, 1 / (4 * D.PI * d3 ** 2), label='Phi / (4 pi d^2)')
ax.set_xlabel('distance (cells)')
ax.set_ylabel('light (W/m, W/m2)')
ax.legend()
plt.show()

n, c = 61, 30 P = np.zeros((n, n)); P[c, c] = 1.0 E, S = point_source(P, np.ones((n, n)), 25) d = np.arange(1, 26) plane = np.array([E[c, c + k] for k in d]) assert np.allclose(plane * 2 * D.PI * d, 1.0, rtol=0, atol=1e-12) and np.all(S == 1.0) Q = np.zeros((25, 25, 25)); Q[12, 12, 12] = 1.0 E3, _ = point_source(Q, np.ones((25, 25, 25)), 11) d3 = np.arange(1, 12) space = np.array([E3[12, 12, 12 + k] for k in d3]) assert np.allclose(space * 4 * D.PI * d3 ** 2, 1.0, rtol=0, atol=1e-12) fig, ax = plt.subplots(figsize=(7, 3.2)) ax.loglog(d, plane, 'o', label='kernel, a plane') ax.loglog(d, 1 / (2 * D.PI * d), label='Phi / (2 pi d)') ax.loglog(d3, space, 's', label='kernel, a space') ax.loglog(d3, 1 / (4 * D.PI * d3 ** 2), label='Phi / (4 pi d^2)') ax.set_xlabel('distance (cells)') ax.set_ylabel('light (W/m, W/m2)') ax.legend() plt.show()

Invariant: lamps add, and the shade is the light's share of the open's¶

The law is linear in the power: two lamps' light is the sum of each one's, walls and all. The shade is a share, between 0 and 1, and the light is the shade times the light in the open.

In [4]:

rng = np.random.default_rng(1)
t = rng.uniform(0.3, 1.0, (31, 31))
a = np.zeros((31, 31)); a[10, 10] = 1.0
b = np.zeros((31, 31)); b[20, 22] = 3.0
both, S = point_source(a + b, t, 12)
assert np.allclose(both, point_source(a, t, 12)[0] + point_source(b, t, 12)[0], rtol=0, atol=1e-14)
open_ = point_source(a + b, np.ones((31, 31)), 12)[0]
assert np.all((S >= 0) & (S <= 1)) and np.allclose(both, S * open_, rtol=0, atol=1e-14)
print('superposition holds; the light is the shade times the open light')

rng = np.random.default_rng(1) t = rng.uniform(0.3, 1.0, (31, 31)) a = np.zeros((31, 31)); a[10, 10] = 1.0 b = np.zeros((31, 31)); b[20, 22] = 3.0 both, S = point_source(a + b, t, 12) assert np.allclose(both, point_source(a, t, 12)[0] + point_source(b, t, 12)[0], rtol=0, atol=1e-14) open_ = point_source(a + b, np.ones((31, 31)), 12)[0] assert np.all((S >= 0) & (S <= 1)) and np.allclose(both, S * open_, rtol=0, atol=1e-14) print('superposition holds; the light is the shade times the open light')

superposition holds; the light is the shade times the open light

Limit: a wall that lets nothing through¶

An opaque wall is lit on its face and shades everything behind it to nothing; a wall letting half through halves the light straight behind it (an oblique ray crosses a thin wall in more cells, its longer path through it, and is shaded more). The shade is the shadow map a view draws.

In [5]:

n, c = 41, 20
P = np.zeros((n, n)); P[c, c] = 1.0
open_, _ = point_source(P, np.ones((n, n)), 15)
for share in (0.0, 0.5):
    w = np.ones((n, n)); w[:, c + 4] = share
    E, S = point_source(P, w, 15)
    assert np.allclose(E[c, c + 5:], share * open_[c, c + 5:], rtol=0, atol=1e-15)
    assert E[c, c + 4] == open_[c, c + 4]
    if share == 0.0:
        assert float(E[:, c + 5:].sum()) == 0.0
w = np.ones((n, n)); w[5:36, c + 4] = 0.0
E, S = point_source(P, w, 15)
fig, ax = plt.subplots(1, 2, figsize=(8, 3.6))
ax[0].imshow(np.log10(E + 1e-6)); ax[0].set_title('light')
ax[1].imshow(S, cmap='gray', vmin=0, vmax=1); ax[1].set_title('shade')
plt.show()

n, c = 41, 20 P = np.zeros((n, n)); P[c, c] = 1.0 open_, _ = point_source(P, np.ones((n, n)), 15) for share in (0.0, 0.5): w = np.ones((n, n)); w[:, c + 4] = share E, S = point_source(P, w, 15) assert np.allclose(E[c, c + 5:], share * open_[c, c + 5:], rtol=0, atol=1e-15) assert E[c, c + 4] == open_[c, c + 4] if share == 0.0: assert float(E[:, c + 5:].sum()) == 0.0 w = np.ones((n, n)); w[5:36, c + 4] = 0.0 E, S = point_source(P, w, 15) fig, ax = plt.subplots(1, 2, figsize=(8, 3.6)) ax[0].imshow(np.log10(E + 1e-6)); ax[0].set_title('light') ax[1].imshow(S, cmap='gray', vmin=0, vmax=1); ax[1].set_title('shade') plt.show()