Table of Contents
The silver halide film grain has an unique look that introduces some texture in smooth image areas, where digital imagery may look synthetic and too clean. Emulating this effect in digital photography cannot be done by simply adding gaussian or poissonian noise to the image, since silver halide grain has varying shapes and sizes. There are two techniques to do it :
 overlaying an image of natural grain sampled on scanned negatives, on top of the digital image (typically using softlight or hardlight blend modes),
 synthesizing grain from an apriori model and the image content.
The first method is unsatisfactory because it implicitely uses the additive noise model, inherited from signal processing, which silver halide grain is not : far from being noise, grains are the actual sensors. The second method has been investigated with MonteCarlo simulation^{[1]}, producing very good visual results but with runtimes unsuited for realtime applications and requiring large amounts of memory. It also relies solely on a probabilistic approach that doesn’t take into account crystallographic structures, so the input parameters are probabilistic distribution parameters that won’t speak much to photographers.
I propose here a method to simulate a random crystallographic structure of photographic emulsion as a stack of crystal layers parametrized by the size of the grains and the concentration of silver halide (AgX) in the emulsion. The purpose of this simulation is to be fairly physicallyrealistic while staying computationnalyaffordable and producing goodlooking grain, but also to provide user parameters that can be intuitively understood by the layman.
#Grain formation↑
A photographic emulsion is silver halide (AgX) crystals uniformingly dispersed into a translucent gelatin medium. Those crystals may have varying shapes and sizes depending on the experimental conditions of production. When enough photons hit them, some $Ag^+$
ions from the surface of the crystal are turned into metallic silver, forming surface clusters. When those clusters are large enough, during the developping process, they will catalyse the reduction of the whole crystal by the chemicals. At this stage, the surface of the crystals will act as an amplification parameter.
#Model↑
#Interlayers model↑
Let us posit that the photographic emulsion can be represented as a stack of elementary crystal layers having each the depth of a crystal. Cubic crystals of lowspeed emulsions have a typical size of 0.5 to 0.7 µm ^{[2]}. Kodak tabular grains have a the property of being large but shallow (noncubic), with a depth as small as 0.12 μm ^{[3]}. The typical depth of photographic emulsions is 15 μm, which gives us 21 to 30 elementary layers for nontabular grains, and up to 125 for tabular grains.
Each layer will receive a fraction of the original light emission, which is bounded by the amount of light already captured at higher layers and by the saturation level of the layer itself (that is, the number of photons that would completely reduce the silver halide into metallic silver). This saturation level is a function of the grain surface, meaning larger grains have an higher saturation level and an increased light sensibility, and it appears that there is a minimum in the size of silver cluster below which they can’t catalize a reduction ^{[4]}.
Photons that were not captured on higher layers will move on to the lower ones, with some amount of scattering, resulting in some image blurring. At the end, it is not expected that all input photons will be captured and therefore the conversion efficiency of the photochemical reaction is expected lower than 100 %. In particular, some vertical sections of the emulsion can be completely free from silver halide crystals, producing translucent “noise” in otherwise dense areas, producing “black grain” (actually, grain holes) after inversion on the final print.
If we take the digital scenelinear photograph $I$
as a 2D projection of the input light field in the direction normal to the sensor plane (using a single RGB channel or a singlechannel monochrome conversion), then the amount of light captured on the whole surface, at each layer $n$
of a stack of $N$
layers at pixel $(x, y)$
is modeled by :
$$
L_n(x, y) = I  \sum_{i = 0}^{n1}L_i(x, y), 0 \leq L_n(x, y) \leq S(x, y)
$$
where $S(x, y)$
is the saturation level of the grain covering $(x, y)$
on that layer, or $0$
if that coordinate does not have a grain. The final image can therefore be reconstructed by :
$$
I_{out}(x, y) = \bar{c} \sum_{n=0}^{N}L_n(x, y)
$$
where $\bar{c}$
is a gain compensating for light losses through grain holes, such that the average exposure of $I_{out}$
matches the average exposure of $I$
, so $\bar{c} = \frac{\overline{I}}{\overline{I_{out}}}$
despite the efficiency being lower than 100 %. This is because digital images are assumed to have been already corrected for exposure and brightness before the grain synthesis, so the model should respect the brightness intent.
#Crystal model↑
Silver halide crystals typically have a cubic or octahedric shape, depending on how they were prepared, with some variability in the number of vertices (triangles and hexagons can also be found).
The polar equation for the outer border of any regular polyhedron is :
$$
r(\theta, \phi, n) = \dfrac{\cos\left(\frac{\pi}{n}\right)}{\cos\left(\frac{2 \arcsin(\cos(n \theta + \phi)) + \pi}{2 n}\right)}\label{crystal}
$$
where $r$
is the radial distance, $\theta$
is the angular parameter, $\phi$
is the orientation angle of the polyhedron and $n \geq 3$
is the number of vertices. From this equation, we can create 2D kernels for polyhedric crystals with :
$$
K(x, y) = \begin{cases}1 \text{ if } r(\theta, \phi, n) \geq \sqrt{x^2 + y^2}  \epsilon\\ 0 \text{ if } r(\theta, \phi, n)< \sqrt{x^2 + y^2}  \epsilon\end{cases}, \theta = \arctan(y / x)
$$
where $(x, y)$
are the normalized kernel coordinates expressed from the center, in $[1; 1]$
(which suggests the use of square and oddsized kernels), and $\epsilon$
is the discretization error, equal to the inverse of the size. For $n = 5$
, $\phi = 0$
and an 11taps kernel, we get :
[[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
[0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
[0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
[0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
[0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
[0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
[0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
[0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
[0., 0., 0., 1., 1., 1., 1., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.]]
The parameters $n$
, $\phi$
and the size of the kernel can be each mapped to a generator of random variable following a normal distribution (for $n$
and $\phi$
a lognormal distribution of an userdefined standard deviation for the size. Indeed, distributions of particles sizes are commonly found to follow the lognormal distribution and this will make our results comparable with Newson & al. ^{[1]}.
The number of vertices should follow a gaussian distribution centered in $\mu = 6$
since most AgX crystals seem to be either cubes or octahedra.
#Intralayer model↑
The stochastic aspect of crystallization should happen within each crystal layer. However, since we will simulate and stack several layers on top of each other, for a sufficient number of layers, the grain appearance will result from the vertical averaging of the crystal distribution.
We posit that the vertical arrangement of crystal shapes does not matter as long as there is enough variety of shapes along the vertical axis, at each position on the film plane. Therefore, we will enforce only one crystal shape and size per layer, which will allow us to simulate the crystal growth numerically using a simple 2D convolution product of a seeds map by one single kernel, allowing fast evaluation (possible using FFTbased convolution) and efficient (contiguous) memory access.
Since the whole crystal should have the same optical density over its surface, and since we build nonnormalized kernels that will act as enlargers, the whole problem is reduced to uniformingly initializing the seeds maps with a controlled seed density and initializing each seed with the proper density.
The number of crystals $n_a$
per layer should follow a poissonian distribution with parameter $\lambda = n_a$
. For a large enough $\lambda$
parameter, the poissonian distribution can be approximated with a gaussian distribution of parameters $\mu = \lambda$
and $\sigma = \sqrt{\lambda}$
. Using a random variable generator following a gaussian distribution of parameters $\mu, \sigma$
, evaluating a random number for each pixel of the digital image, we can control the surface density of seeds by thresholding the magnitude of the random variable while still having an uniform random seed map. Indeed, the cumulative distribution function of the gaussian distribution is :
$$
p = \frac{1}{2}\left(1 + \text{erf}\left(\frac{x  \mu}{\sigma\sqrt{2}}\right)\right)
$$
This means that we find a population $p$
(expressed as a ratio of the total sample) below the variable threshold $x$
. To find the threshold $x$
below which we find a target population, we can derive the inverse equation :
$$
x = \sigma \sqrt{2} \, \text{erf}^{1}(2 p  1)
$$
This equation is exact ± the sampling errors and allows to define a surface filling ratio for each layer as the ratio $seeds / pixels$
. However, since we are going to turn seeds into crystals by enlarging them, we need to account for the size of the crystal but also for the probability of overlapping crystals in computing that surface filling ratio. A theoritical probabilistic study would be interesting here, but is beyond my skills, so I performed a parametric study of the filling ratio and used a numerical fitting of the relationship between $p$
and the crystal surface. The empiric fitting found is :
$$
p' = \frac{1.42225247}{A} \frac{p^{1.31315748}}{p  1^{0.26087773}}
$$
with $A$
, the surface of the grain, that can be computed from the kernel as a direct sum of all elements. This fitting is mostly accurate for $p \in[50; 95]\%$
, it tends to underestimate (by an offset of 23 %) the filling ratio below 50 %.
In practice, the most relevant part of the domain is $[0; 80] \%$
, redoing a fitting for that part only provides a better estimation below 50 % and yields :
$$
p'' = \frac{1.10724714}{A} \frac{p^{1.04877389}}{p  1^{0.37207405}}
$$
Crystal seeds will therefore be planted at each pixel where the gaussian random variable generator yields a value $v$
such that :
$$
v < \sqrt{p\Omega} \sqrt{2} \, \text{erf}^{1}(2 p''  1) + p\Omega
$$
where $p$
is the user parameter linked to the target surface filling ratio, $\Omega$
is the number of pixels in the image and $p''$
is the empiric fitting between seed ratio and grain surface ratio above. Ilford emulsions from the 1960’s have been documented to use cubic crystals in concentrations of 15 % (in volume), that directly gives a surface filling ratio of 15 % considering crystal layers which depth has the dimension of a single crystal ^{[2]}. Remains now to assign a density to seeds (which will give their density to crystals grown with unnormalized kernels).
On negative film, the density is the level of “blackness”, expressed as the logarithm of the optical transmittance $T$
^{[4]} :
$$
D = \log(T) = \log\left(1  \frac{n_A \bar{a}}{A}\right) = S(E)\label{density}
$$
where $A$
is the surface of measurement, $n_A$
is the number of crystal observed on the surface $A$
and $\bar{a}$
is the average crystal surface. This density is a sigmoidal function of the exposure (expressed as a metric of the number of photons received) $S(E)$
. When printing the negative on paper, we add a second stage of sigmoidal transmittance to density conversion.
Digital images use an additive light scheme (RGB in, out and through), instead of a subtractive light scheme relying on dyes transmittance. Also, the exposure → density mapping function $S(E)$
is handled in digital images through usual tone curves. All in all, the current model should only account for how the original luminous signal is splitted and discretized into the elementary “light holders” that are grains.
Let us define the complement of the transmittance, $T' = 1  T = n_A \frac{\bar{a}}{A}$
. For any subregion of surface $A$
of the original image $I$
, we should have $T' = I_A$
, where $I_A = \frac{1}{A}\int_A I \mathrm{d}A$
is the local average of $I$
over the subregion $A$
. The $n_A \bar{a}$
member represents the average surface filling ratio. We can rewrite it $n_A \bar{a} = \int_N\int_A \dot{k} \, \mathrm{d}A \, \mathrm{d}N$
where $\dot{k}$
is a boolean random variable modelling whether the elementary surface $\mathrm{d}A$
is filled by a grain or not, on an elementary layer $\mathrm{d}N$
. Overall, this yields :
$$
\int_A I \mathrm{d}A = \int_N\int_A \dot{k} \, \mathrm{d}A \, \mathrm{d}N
$$
This provides the generalized additive light translation of the optical density equation $\eqref{density}$
, which can be numerically solved by discretization technics over 3D layered meshes of arbitrary shapes. For a large enough domain $A$
over a single elementary layer $\mathrm{d}N$
, $\int_A \dot{k} \, \mathrm{d}A$
is known as the userdefined surface filling ratio $p$
. Locally, the probability of $\dot{k} = 1$
is equal to that ratio. Over the vertical section $\mathrm{d}A$
, the average number of crystals is then $\overline{n_A} = \int_N p \, \mathrm{d}N = pN$
and the maximum number of crystals is $N$
.
From this, we deduce that the image formation process over a vertical section of film emulsion of elementary surface $\mathrm{d}A$
decomposed into $N$
layers can be seen as a decomposition process of the original image intensity $I$
where an average of $pN$
grains capture each an equal fraction of $I$
until they have completely depleted the light emission $I$
of its photons.
We will therefore assign an intensity $\frac{I}{N}$
to each seed, which will be propagated to the whole crystal through the convolution by a nonnormalized kernel. This will lead to :
$$
\sum_N \dot{k} \frac{I}{N} < I
$$
in other words, the result will be underexposed compared to the initial light emission, which is consistent with how film works (the initial exposure producing metallic silver cluster over silver halide crystals surface, that will later catalyse the full reduction of the crystals by the developer). This can be fixed by an exposure compensation $\bar{c} = \overline{I} / \overline{\sum_N \dot{k} \frac{I}{N}}$
at the final stage.
#Printing model↑
The areas from the negative film where density is maximal will print white and will fully block light, acting as an opaque mask. On these regions, the grain should not be apparent. However, our model does not honour this rule and grain will appear in those regions. To fix this, we will use an alpha (opacity) compositing model where the grainy version of the image is overlayed on top of the original smooth version, using the normalized original as the alpha mask as $\alpha = 1  I$
(assuming $I$
white point is normalized to 1) :
$$
I_{final} = \alpha I_{grainy} + (1  \alpha) I
$$
#Summary↑
For each layer $L_n$
of a stack of $N$
layers :
 create a random crystal kernel
$K$
using equation$\eqref{crystal}$
where the parameters are set by random variable generators (following normal distribution for number of vertices and angle, and lognormal for the size),  seed the layer pixels of coordinates
$(x, y)$
with the value$I(x, y) / N$
where a random variable generator following a normal distribution$\mu = p\Omega$
and$\sigma = \sqrt{p\Omega}$
yields a value$v$
such that$v < \sqrt{p\Omega} \sqrt{2} \, \text{erf}^{1}(2 p''  1) + p\Omega$
. Wherever the condition is not met on$v$
, set the pixels to 0.  convolve the seed map
$L_n$
with the crystal kernel$K$
. Since the kernel is not normalized, it does not preserve energy and pixels where 2 crystals overlap will be overbrightened. The condition$(L_n * K)(x, y) \leq I(x, y) / N$
will need to be manually enforced after the convolution product.
Then, aggregate all the layers as $I_{temp} = \sum_n L_n$
and correct the global exposure such that $I_{out} = I_{temp} \bar{I}/\overline{I_{temp}}$
.
Finally, apply the printing model using alpha compositing $I_{final} = \alpha I_{grainy} + (1  \alpha) I$
.
#Results↑
#Varying surface filling ratio and number of layers↑
The average grain size is set at 5 px with a lognormal standard deviation of 0.25. The surface filling ratio is the ratio of the total crystal surface by the total film surface over a single layer. Note that real film surface fillings should be between $[15;30] \%$
and nontabular (monochrome) emulsions should have 20 to 30 layers. Tabular grains in color emulsions could have a maximum of 125 layers but splitted between the 3 color channels.
Adding more layers dilute and average grains in a way that make them effectively disappear.
#Varying grain size and distribution variability↑
The number of layers is set at 30 and the surface filling ratio at 25 %. Note that the grain size in real film is linked to its ISO/ASA sensibility and the variability of the grain size is linked to the quality of the manufacturing process and the experimental conditions of the chemicals manipulation.
Adding variability in the grain size produces a more flaky and cloudy grain texture.
#Natural images↑
These are from my personal collection. Grain is added after the “S” tonemapping curve, I have not investigated yet how it would look when applied before the curve. The fullresolution files are given, click on the images to open them.
For the following color example, grain was applied in Adobe RGB color space. For better accuracy, a color profile emulating film color response should be used.
#Film scan samples↑
These 2 samples are again from my personal collection, and given here for reference of how film grain could look.
A word of advice : “film grain” does not exist. First of all, the structure of grain depends on how the chemical reactions happened : how the film was conserved before exposure, how much time it waited before exposure and developement, whether it was pulled or pushed at exposure, whether the developer was fresh or already used, stirred or still during development. But then, the appearance of grain varies depending on the viewing conditions : until the 1970’s, most professional film labs would use collimated light under enlargers, allowing sharper rendition of both the image content and the film grain. Because this is not forgiving to scratches and dust (so more difficult to work with and more timeconsuming), diffuse light started to take over during the 1970’s and was mainstream by the 1990’s. The same thing happens when scanning. Then, scanning software typically postprocess images, adding more or less sharpening and the scanners may add their own noise on top of the grain.
That is to say that modern emulsions are a lot less grainy than what people may expect (or fantasize), especially when seen through diffuse light, unless expired film or alternative developers (vitamin C, caffenol) were used.
#Discussion↑
The vertical pixel averaging over $N$
layers acts as a subpixel occlusion model, allowing to represent grains partially occluding the pixel area without having to oversample the image : the complexity of the algorithm is $\mathcal{o}N$
. This is a wellknown scheme in computer graphics when compositing nonsquare objects on top of each other : the average intensity over the full pixel is the weighted average of all layers, where the weights are either the opacities or the ratio of overlapping, which falls back to the same idea of average opacity.
The limitation of the model is to be resolutiondependent, and I don’t think it is a real problem. The mediumformat film scan samples above show that, even zoomed in at 1:1 over the 58 Mpx file, grains appear at least several pixels wide, even though silver halide crystals have a radius of 0.5 to 0.7 µm, which would lead to an underlying halftone image of equivalent resolution 3.5 gigapixels for 135 film or 19.6 gigapixels^{[5]} for 120 film in square format. This means that, even on the best possible film scan, each pixel is already a local averaging of the underlying halftone image, whith a local density predicted by equation $\eqref{density}$
. This also mean that the perceived grain is not directly the silver crystals, but rather the spatial variability of the crystal clustering appearing within areas of constant locallyaveraged density.
It is interesting to note that my model is actually wrong on that regard : we simulate crystals as polyhedra which diameter is several pixels in the output image, while actual silver halide crystals have the physical dimension of a 1/100th of pixel by today’s sensor resolution. This scale discrepancy and the fact that we allow “crystals” to overlap here means we actually model random aggregates. Nevertheless, results show something consistent with natural grain.
Based on this, I argue that the boolean elementary crystals prior used by Newson and al. ^{[1]} is an overaccurate description of the problem at a wrong scale, which triggers overly complex probability computations that will get averaged anyway at the pixel scale later in the process, but at the cost of much more expensive computations.
#Conclusion↑
I have shown here a model inspired by crystallography allowing to generate grain from an input RGB image, allowing to virtually build a film emulsion from its locallyaveraged physical and chemical properties, and to simulate how it catches the light from the input image by decomposing the input signal into a random stack of virtual crystals. The user parameters of the model are linked to the physical structure of the film emulsion and can be easily (graphically) explained to a photographer. They can also be traced back to real film emulsions properties. Finally, since we separate surface and depth parameters, tabular grains emulsions can be simulated by adding more crystal planes or larger crystals in the simulation.
This method is faster than a fully stochastic approach based on MonteCarlo schemes because, not only can it be parallelized, but the memory accesses are contiguous and the algorithm uses typical 2D convolutions. The Python code below (unparallelized with a full array I/O at each matrix operation) runs in less than 8 s on singlecore CPU for 30 layers of grains at an average of 3 px, over a 1361×2048 px image. This means nearrealtime performance at desktop resolution once translated into C/C++. It doesn’t allow to isolate grains when magnifying images above reason, but this cool trick has no application in a photography workflow.
#Source code↑
Python 3.11
import numpy as np
from PIL import Image
from scipy.signal import convolve2d
from scipy.ndimage import gaussian_filter
from scipy.special import erfinv, erf
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
GAMMA = 2.4
def show_image(array):
# Regamma and convert to 8 bits
result = Image.fromarray(np.uint8(array**(1. / GAMMA) * 255))
result.show()
return result
def open_image(path):
# Open, ungamma and convert to monochrome using luminance
image = Image.open(path)
im = (np.array(image) / 255.)**GAMMA
return 0.20 * im[:,:,0] + 0.75 * im[:,:,1] + 0.05 * im[:,:,2]
def create_crystal(width, n, rotation):
# w is the width of the kernel, aka "grain". Needs to be odd.
# n is number of vertices: 3 for triangle, 4 for square, 5 for pentagon, etc.
# see https://math.stackexchange.com/a/4160104/498090
# rotation is an angle in radians
# blur controls how fuzzy the crystal edges are
# Spatial coordinates rounding error
eps = 1. / width;
radius = max(int((width  1) / 2.), 1) # exact if width is odd
kernel = np.empty((width, width))
for i in range(width):
for j in range(width):
# get normalized kernel coordinates from center in [1 ; 1]
x = i / radius  1.
y = j / radius  1.
# get current radial distance from kernel center
r = np.hypot(x, y)
# get the radial distance at current angle of the shape envelope
M = np.cos(np.pi / n) / np.cos((2. * np.arcsin(np.cos(n * (np.arctan2(y, x) + rotation))) + np.pi) / (2. * n))
# write 1 if we are inside the envelope of the shape, else 0
kernel[i, j] = (M >= r  eps)
return kernel
def pick_crystal(size, std):
# Pick one crystal size, shape and orientation based on random dice rolls
rng = np.random.default_rng()
dice_shape = np.clip(rng.normal(6, 1.5), 3, 10)
dice_rotation = np.random.rand() * 2. * np.pi
log_normal_var = np.random.lognormal(mean=np.log(size), sigma=std)
random_size = int(min(max(log_normal_var, 1), 3 * size))
# ensure kernel size is odd
if random_size % 2 == 0.:
random_size += 1
return create_crystal(random_size, dice_shape, dice_rotation)
def distribution_to_variable(population, sigma):
# Find t so the normal centered unit gaussian distribution G(x, mu=0, sigma=1)
# contains the population ratio below x = t.
# population is in[0; 1]
return erfinv(2. * population  1.) * np.sqrt(2.) * sigma
def variable_to_distribution(x, sigma=1.):
# Find the ratio of the population of the normal centered unit
# gaussian distribution G(x, mu=0, sigma=1)
# contained below x
return (1. + erf(x / (sigma * np.sqrt(2)))) / 2.
def filling_to_rand_variable(x):
# Map a filling ratio with a gaussian distribution variable
# fitting in [0; 1] :
# return 1.3790633 * x*1.34123704 / np.abs(x  1)**0.27288849
# fitting in [0; 0.8] :
return 1.12139063 * x*1.05577624 / np.abs(x  1)**0.34443235
def grainify(image, filling=0.9, size=9, n=20, std=1):
"""
Params:
 image (array): B&W image or single color channel
 filling (float): ratio of surface filling with AgX crystals.
 size (int): average pixel size of the grain (odd values only)
 n (int): number of layers of silver halide.
 std (float): the standard deviation of the lognormal distribution of crystal sizes.
Photo emulsions are 10 to 15 μm thick, [1][2]
 nontabular grains have a typical diameter between 0.5 and 0.7 μm (up to 3μm) [1][2]
 tabular grains have a thickness of 0.12 micrometer but are used for color emulsions (×3) [4]
That gives us 14 to 40 crystal layers in practice.
The surface filling ratio was 15% of Ilford emulsions in the 1960's. [1]
Visual results of the simulation indicate Ilford Delta 100 could be around 25 %.
References :
[1] https://cds.cern.ch/record/870005/files/p129.pdf
[2] https://www.epfl.ch/labs/gdp/wpcontent/uploads/2019/09/PC2_Lesson_10.pdf
[3] https://pubs.acs.org/doi/pdf/10.1021/bk19820200.ch001
[4] https://www.jstage.jst.go.jp/article/photogrst1964/49/6/49_6_499/_pdf
"""
result = np.zeros_like(image)
crystals = np.zeros_like(image)
final_filling = 0.
sigma = filling_to_rand_variable(filling)
layer_density = image / n
for i in range(n):
# Init a random crystal shape and size for the current layer
crystal_area = 0
while crystal_area == 0.:
# Some shapes applied on small kernel sizes lead to null kernels.
# In that case, try again.
crystal = pick_crystal(size, std)
crystal_area = crystal.sum()
# Ensure a constant surface filling over layers, no matter the crystal size
n_a = filling * image.size / crystal_area
value = filling if crystal_area == 1 else sigma / crystal_area
bound = distribution_to_variable(value, np.sqrt(n_a))
# Lightness value for seeds
layer = np.clip(layer_density, 0, (image  result))
# Init random seeds over the surface of the layer
seeds = np.where(np.random.normal(n_a, np.sqrt(n_a), image.shape) < n_a + bound, layer, 0.)
# Grow the crystals by convolving seeds with nonnormalized crystal kernel.
grains = convolve2d(seeds, crystal, mode='same', boundary='symm')
# Crystal kernel is not normalized (to enlarge seeds), meaning when 2 crystals
# overlap, energy is added as if AgX captured more photons than available.
# Normalize it now for energy conservation.
grains = np.where(grains > layer, layer, grains)
result += grains
# Adjust global exposure to match the original
coef = np.mean(image) / np.mean(result)
print("exposure coef:", coef)
grainy = np.clip(result * coef, 0., 1.)
# Printing model of opacity : fully white on paper is fully opaque on negative
# so grains should not appear there. Use alpha blending to mix grain with the original smooth pic.
mask = 1.  image # assuming white = 1, otherwise normalize first
return np.clip(mask * grainy + (1.  mask) * image, 0., 1.), final_filling / n
Execution over a picture :
image = open_image("/path/to/image.jpg")
grainy, fill_ratio = grainify(image, filling=0.25, size=7, n=30, std=0.25)
to_save = show_image(grainy)
to_save.save("result.jpg") # optional

NEWSON, Alasdair, FARAJ, Noura, GALERNE, Bruno, et al. Realistic film grain rendering. Image Processing On Line, 2017, vol. 7, p. 165183. https://www.ipol.im/pub/art/2017/192/article.pdf ↵ ↵ ↵

DAHLJENSEN, Erik. V. 1 The properties of photographic emulsions. 1963. https://cds.cern.ch/record/870005/files/p129.pdf ↵ ↵

KOFRON, J. T. et BOOMS, R. E. Kodak Tgrain emulsions in color films. Journal of The Society of Photographic Science and Technology of Japan, 1986, vol. 49, no 6, p. 499504. https://www.jstage.jst.go.jp/article/photogrst1964/49/6/49_6_499/_pdf ↵

RAMSDEN, J. J. The photolysis of small silver halide particles. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 1984, vol. 392, no 1803, p. 427444. ↵ ↵

“pixel” is here to be understood as “elementary image cell”, regardless of technological means to render images. ↵