Stochastic photographic grain synthesis from crystallographic structure simulation

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 soft-light or hard-light blend modes ),
• synthesizing grain from an a-priori 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 Monte-Carlo simulation¹, producing very good visual results but with runtimes unsuited for real-time 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 physically-realistic while staying computationnaly-affordable and producing good-looking 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

Inter-layers 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 low-speed emulsions have a typical size of 0.5 to 0.7 µm ². Kodak tabular grains have a the property of being large but shallow (non-cubic), with a depth as small as 0.12 μm ³. The typical depth of photographic emulsions is 15 μm, which gives us 21 to 30 elementary layers for non-tabular 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 ⁴.

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 scene-linear 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 single-channel 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}^{n-1}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 odd-sized kernels), and $\epsilon$ is the discretization error, equal to the inverse of the size. For $n = 5$, $\phi = 0$ and an 11-taps kernel, we get:

 1[[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
 2 [0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
 3 [0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
 4 [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
 5 [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
 6 [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
 7 [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
 8 [0., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.],
 9 [0., 0., 1., 1., 1., 1., 1., 1., 1., 0., 0.],
10 [0., 0., 0., 1., 1., 1., 1., 1., 0., 0., 0.],
11 [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 log-normal distribution of an user-defined standard deviation for the size. Indeed, distributions of particles sizes are commonly found to follow the log-normal distribution and this will make our results comparable with Newson & al. ¹.

The number of vertices should follow a gaussian distribution centered in $\mu = 6$ since most AgX crystals seem to be either cubes or octahedra.

Intra-layer 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 FFT-based convolution) and efficient (contiguous) memory access.

Since the whole crystal should have the same optical density over its surface, and since we build non-normalized 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 2-3 %) 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 ². Remains now to assign a density to seeds (which will give their density to crystals grown with un-normalized kernels).

On negative film, the density is the level of “blackness”, expressed as the logarithm of the optical transmittance $T$⁴:

$$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 sub-region 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 sub-region $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 user-defined 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 non-normalized 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 :

1. 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 log-normal for the size),
2. 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.
3. 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 log-normal 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 non-tabular (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 full-resolution 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 time-consuming), 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 post-process 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 sub-pixel 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 well-known scheme in computer graphics when compositing non-square 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 resolution-dependent, and I don’t think it is a real problem. The medium-format 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⁵ 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 locally-averaged 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. ¹ is an over-accurate 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 locally-averaged 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 Monte-Carlo 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 (un-parallelized with a full array I/O at each matrix operation) runs in less than 8 s on single-core CPU for 30 layers of grains at an average of 3 px, over a 1361×2048 px image. This means near-real-time 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

  1import numpy as np
  2from PIL import Image
  3from scipy.signal import convolve2d
  4from scipy.ndimage import gaussian_filter
  5from scipy.special import erfinv, erf
  6from scipy.optimize import curve_fit
  7import matplotlib.pyplot as plt
  8
  9GAMMA = 2.4
 10
 11def show_image(array):
 12    # Re-gamma and convert to 8 bits
 13    result = Image.fromarray(np.uint8(array**(1. / GAMMA) * 255))
 14    result.show()
 15    return result
 16
 17def open_image(path):
 18    # Open, un-gamma and convert to monochrome using luminance
 19    image = Image.open(path)
 20    im = (np.array(image) / 255.)**GAMMA
 21    return 0.20 * im[:,:,0] + 0.75 * im[:,:,1] + 0.05 * im[:,:,2]
 22
 23def create_crystal(width, n, rotation):
 24    # w is the width of the kernel, aka "grain". Needs to be odd.
 25    # n is number of vertices: 3 for triangle, 4 for square, 5 for pentagon, etc.
 26    # see https://math.stackexchange.com/a/4160104/498090
 27    # rotation is an angle in radians
 28    # blur controls how fuzzy the crystal edges are
 29
 30    # Spatial coordinates rounding error
 31    eps = 1. / width;
 32    radius = max(int((width - 1) / 2.), 1) # exact if width is odd
 33    kernel = np.empty((width, width))
 34
 35    for i in range(width):
 36        for j in range(width):
 37            # get normalized kernel coordinates from center in [-1 ; 1]
 38            x = i / radius - 1.
 39            y = j / radius - 1.
 40
 41            # get current radial distance from kernel center
 42            r = np.hypot(x, y)
 43
 44            # get the radial distance at current angle of the shape envelope
 45            M = np.cos(np.pi / n) / np.cos((2. * np.arcsin(np.cos(n * (np.arctan2(y, x) + rotation))) + np.pi) / (2. * n))
 46
 47            # write 1 if we are inside the envelope of the shape, else 0
 48            kernel[i, j] = (M >= r - eps)
 49
 50    return kernel
 51
 52
 53def pick_crystal(size, std):
 54    # Pick one crystal size, shape and orientation based on random dice rolls
 55    rng = np.random.default_rng()
 56
 57    dice_shape = np.clip(rng.normal(6, 1.5), 3, 10)
 58    dice_rotation = np.random.rand() * 2. * np.pi
 59
 60    log_normal_var = np.random.lognormal(mean=np.log(size), sigma=std)
 61    random_size = int(min(max(log_normal_var, 1), 3 * size))
 62
 63    # ensure kernel size is odd
 64    if random_size % 2 == 0.:
 65        random_size += 1
 66
 67    return create_crystal(random_size, dice_shape, dice_rotation)
 68
 69def distribution_to_variable(population, sigma):
 70    # Find t so the normal centered unit gaussian distribution G(x, mu=0, sigma=1)
 71    # contains the population ratio below x = t.
 72    # population is in[0; 1]
 73    return erfinv(2. * population - 1.) * np.sqrt(2.) * sigma
 74
 75def variable_to_distribution(x, sigma=1.):
 76    # Find the ratio of the population of the normal centered unit
 77    # gaussian distribution G(x, mu=0, sigma=1)
 78    # contained below x
 79    return (1. + erf(x / (sigma * np.sqrt(2)))) / 2.
 80
 81def filling_to_rand_variable(x):
 82    # Map a filling ratio with a gaussian distribution variable
 83    # fitting in [0; 1] : 
 84    # return 1.3790633 * x*1.34123704 / np.abs(x - 1)**0.27288849
 85    # fitting in [0; 0.8] : 
 86    return 1.12139063 * x*1.05577624 / np.abs(x - 1)**0.34443235
 87
 88def grainify(image, filling=0.9, size=9, n=20, std=1):
 89    """
 90    Params:
 91        - image (array): B&W image or single color channel
 92        - filling (float): ratio of surface filling with AgX crystals.
 93        - size (int): average pixel size of the grain (odd values only)
 94        - n (int): number of layers of silver halide.
 95        - std (float): the standard deviation of the log-normal distribution of crystal sizes.
 96
 97    Photo emulsions are 10 to 15 μm thick, [1][2]
 98        - non-tabular grains have a typical diameter between 0.5 and 0.7 μm (up to 3μm) [1][2]
 99        - tabular grains have a thickness of 0.12 micrometer but are used for color emulsions (×3) [4]
100    That gives us 14 to 40 crystal layers in practice.
101
102    The surface filling ratio was 15% of Ilford emulsions in the 1960's. [1]
103    Visual results of the simulation indicate Ilford Delta 100 could be around 25 %.
104
105    References : 
106         [1] https://cds.cern.ch/record/870005/files/p129.pdf
107         [2] https://www.epfl.ch/labs/gdp/wp-content/uploads/2019/09/PC2_Lesson_10.pdf
108         [3] https://pubs.acs.org/doi/pdf/10.1021/bk-1982-0200.ch001
109         [4] https://www.jstage.jst.go.jp/article/photogrst1964/49/6/49_6_499/_pdf
110
111    """
112    result = np.zeros_like(image)
113    crystals = np.zeros_like(image)
114
115    final_filling = 0.
116    sigma = filling_to_rand_variable(filling)
117    layer_density = image / n
118
119    for i in range(n):
120        # Init a random crystal shape and size for the current layer
121        crystal_area = 0
122        while crystal_area == 0.:
123            # Some shapes applied on small kernel sizes lead to null kernels.
124            # In that case, try again.
125            crystal = pick_crystal(size, std)
126            crystal_area = crystal.sum()
127
128        # Ensure a constant surface filling over layers, no matter the crystal size
129        n_a = filling * image.size / crystal_area
130        value = filling if crystal_area == 1 else sigma / crystal_area
131        bound = distribution_to_variable(value, np.sqrt(n_a))
132
133        # Lightness value for seeds
134        layer = np.clip(layer_density, 0, (image - result))
135
136        # Init random seeds over the surface of the layer
137        seeds = np.where(np.random.normal(n_a, np.sqrt(n_a), image.shape) < n_a + bound, layer, 0.)
138
139        # Grow the crystals by convolving seeds with non-normalized crystal kernel.
140        grains = convolve2d(seeds, crystal, mode='same', boundary='symm')
141
142        # Crystal kernel is not normalized (to enlarge seeds), meaning when 2 crystals
143        # overlap, energy is added as if AgX captured more photons than available.
144        # Normalize it now for energy conservation.
145        grains = np.where(grains > layer, layer, grains)
146
147        result += grains
148
149    # Adjust global exposure to match the original
150    coef = np.mean(image) / np.mean(result)
151    print("exposure coef:", coef)
152    grainy = np.clip(result * coef, 0., 1.)
153
154    # Printing model of opacity : fully white on paper is fully opaque on negative
155    # so grains should not appear there. Use alpha blending to mix grain with the original smooth pic.
156    mask = 1. - image # assuming white = 1, otherwise normalize first
157    return np.clip(mask * grainy + (1. - mask) * image, 0., 1.), final_filling / n

Execution over a picture : 

1image = open_image("/path/to/image.jpg")
2grainy, fill_ratio = grainify(image, filling=0.25, size=7, n=30, std=0.25)
3to_save = show_image(grainy)
4to_save.save("result.jpg") # optional