Diffraction effects detection for HDR image-based measurements

Modern imaging techniques have proved to be very efficient to recover a scene with high dynamic range (HDR) values. However, this high dynamic range can introduce star-burst patterns around highlights arising from the diffraction of the camera aperture. The spatial extent of this effect can be very wide and alters pixels values, which, in a measurement context, are not reliable anymore. To address this problem, we introduce a novel algorithm that, utilizing a closed-form PSF, predicts where the diffraction will affect the pixels of an HDR image, making it possible to discard them from the measurement. Our approach gives better results than common deconvolution techniques and the uncertainty values (convolution kernel and noise) of the algorithm output are recovered.


Introduction
In a wide variety of applications, the camera dynamic range does not permit to capture the whole dynamic range of the scene. High dynamic range (HDR) imaging [1] is therefore necessary in order to fully recover the whole scene dynamic range. HDR photography merges photographs of a scene taken at different levels of exposure, in order to increase the native camera dynamic range. HDR images are very useful because they speed up the acquisition process when using imaging devices.
A common artifact arising from the high dynamic range are star burst patterns that appear around highlights. This effect is due to light diffraction through the lens diaphragm and cannot be avoided. From a metrology perspective, these diffraction patterns pollute lots of pixels that cannot be taken as reliable measurements. Since the camera diffraction pattern has a very high dynamic range, the higher the image dynamic range is, the more prominent is the pollution by diffraction. More generally, even if the effect becomes less visible, high value pixels can always affect the lower value pixels through diffraction because the spatial range of diffraction is not bounded. One has to be very careful when considering a low value pixel as a reliable measurement. Mathematically, this diffraction effect is described by a convolution, to which is added a classical measurement noise. Our proposed algorithm aims to detect and discard the pixels polluted by the light diffraction, while keeping the rest as reliable measurements.
Recovering a noisy measurement blurred by a convolution kernel (impulse response) is an issue of main interest since it focuses on removing the impact of the measuring instrument on the acquired data. The main difficulty is that it is an ill-posed mathematical problem (cf. [2], p.7): a solution is not unique, may not exist, and may not be stable. In fact, if the deconvolved solution is not stable, a slight error in the data may lead to a very large error in the solution. It means that for measurement purposes, where noise is always present, recovering the true unbiased data is mathematically impossible. Yet, a wide variety of deconvolution techniques have been developed, divided into 4 major categories: Fourier based techniques, constrained iterative algorithms, entropy maximization, and maximum likelihood estimation (Bayesian methods).
Fourier techniques, such as inverse filtering [3], Wiener filtering [4], CLEAN [5], or the Eldar algorithm [6], suffer from the lack of a priori information and the uniqueness of solution [7]. Noise amplification is also a great issue, even if some of these methods, such as the Eldar algorithm, aim to minimize it . Constrained iterative algorithms try to recover the original measurement by iterative processes under constraints, such as Jansson Van-Cittert algorithms [3,8] and the Gold algorithm [3] (non-negativity of the solution), or the Combettes and Trussel algorithm [9] (bounded noise). These iterative algorithms converge to the inverse filtering but they benefit from less suffering from the noise amplification problem. In this category, modelfitting techniques try to describe the measurement by a set of constrained parameters or basis functions and to find the best match [10][11][12]. However, these fitting methods generally lead to unstable solutions [7], which is a well-known problem, in spectrometry for instance [12,13].
Compared to all these previous techniques, the entropy maximizations algorithm [14,15] performs better in reducing the RMS error on the solution [16].
The most commonly used algorithms [17] are all Bayesian methods [18], based on the maximum likelihood estimation [17]. They are very flexible in the constraint definition, and the output is robust to noise, even more when a regularization function is used [3]. A lot of variation exists: ICTM [19,20], blind deconvolution [21-23], Pixon [24], and a wide variety of algorithms implementing different input constraints. Yet, the famous Richardson-Lucy deconvolution algorithm [18] is the one that performs best compared to its concurrents in MSE error minimization [25]. Thanks to the included constraints, these methods inject a wide variety of a priori information about the noise and the impulse response, therefore lead to a better solution [7].
None of these algorithms guarantee any uncertainty value for the deconvolution output because it depends on the problem unknowns [6,18]. In his original paper [18], Richardson writes that the value of his process is that "it can give intelligible results in some cases where the Fourier process cannot", highlighting the fact that the deconvolution techniques are not aimed at guaranteeing a measurement value.
All in all, the main issue is that deconvolution algorithms are not able to guarantee any boundaries for the recovered pixel value, in spite of a good shape of the reconstructed image. However, when doing metrology-grade measurements, uncertainties are necessary. In this paper, we propose to tackle the problem differently by predicting and identifying the pixels in the image that are polluted by diffraction and then discard them. Since our technique classifies pixels instead of recovering their original value, no pixel value is modified, and therefore, we can keep track of the measurement uncertainty.

Overview of the Method
The first step is to precompute the optical impulse response (also called the point spread function, PSF) of the camera for a given setup. This computation is based on the diaphragm fitting with a proposed model, which is general enough to cover a wide variety of apertures, but also gives a closed-form solution of the PSF. Therefore, our algorithm predicts the amount of diffraction present in the HDR image-based measurement. The algorithm is based on an incremental prediction of the effect of diffraction, from the highest to the lowest pixel values. At first, we discuss the Fourier optics basics needed for evaluating the PSF and the validity conditions implied. Then, we set the parameters used for a real camera lens system, deriving a closed-form solution for the PSF equation. We describe the diffraction detection algorithm that finds the pixels to be discarded from the image. The method is then confronted to simulations and real photographs, and finally the validity of the algorithm is discussed. Fig. 1: Left: Principle of diffraction through a thin-lens camera, composed of a finite aperture lens and a camera sensor separated by a distance D. Even if the object point at d were in the focal plane at d * , the image of the sensor would not be a point as expected, but a pattern due to light behaving as a wave. Right: Mathematical model of a standard n-bladed camera aperture. The full pattern can be divided into similar geometries, themselves sub-divided into two elementary parts : a triangle OAB (blue), and a section of parabola whose axis of symmetry passes through the middle point M (red).

Point Spread Function for a Thin-Lens Camera Model
A lot of projective cameras can be well-described by a thin-lens model, composed of a finite aperture lens of focal length f shaped by a pupil function P(x, y) and a sensor behind it at a distance D (cf. Fig. 1). Because of the wave nature of light, the image of a point through a finite aperture is not a point, but a Point Spread Function (PSF), which depends on the wavelength λ of the light, the camera settings, and the distance d on the optical axis to the camera aperture. From Fourier optics [26], the PSF function is given by where F [·] represents the Fourier transform operator, and S pup the lens surface. Neglecting the lens thickness and aberrations, this formulation is valid under two approximations. Noting f # the lens f-number, the first is the well-known Fresnel approximation such that any wave entering the camera can be considered as a parabolic wave. The second one, called the in-focus approximation, necessitates the scene to be comprised within a distance from the object plane (conjugate of the sensor plane by the lens, cf. Fig. 1) verifying , where d * = D f /(D − f ) is the distance between the lens and the object plane. Verified in most of the real case scenarios, the image formation is considered incoherent. Therefore, each point in the object space contributes to the image formation by adding its own intensity. The PSF then is the function that applies a blur on the perfect image I * , such that the captured image I is given by the following image formation formula: where ⊗ is the convolution operator.

Standard Diaphragm Modeling
The great majority of lens diaphragms are designed with blades, also known as bladed apertures. An aperture with n blades will be referred as n-bladed in the following. In the case of a circular diaphragm, the resulting PSF is the well-known Airy pattern. Shung-Wu and Mittra [27] have studied diaphragms with polygonal shapes but only for straight edges. However, by construction, each blade is an arc of a circle, thus its shape is of constant curvature, giving a good description for any edge of the diaphragm. Some manufacturers change the curvature in a discrete way along the blade edge, yet, for each camera lens f-number, the edge of the blade on the contour of the aperture is in fact an arc of a circle, but of different curvature. Generally, if two consecutive blades cross each other in a certain point, referred as a vertex in the following, the shape described by the set formed by these points is an irregular polygon (cf. Fig. 1 Right). Although one might assume that an aperture is designed to fit a regular polygon, it is not the case because of mechanical constraints between the blades, mostly when they are tight at high f-numbers. Fourier transform being linear, the diaphragm shape has to be mathematically described with independent elements. The origin O is chosen to be geometrical center of the diaphragm surface. Then, for each two neighboring vertices, for instance A and B in the right part of Figure 1, the sub-shape is separated into a triangle (blue) and a section of parabola (red). This sub-shape is described in its local frame defined by a rotation of the main frame from an angle α. In this local frame, points A and B can be designated by coordinates (x, y) = (r AB , z A ) and (r AB , z B ), respectively. The parabola is the only parabola of curvature C at origin, passing by A and B, whose axis of symmetry is the x-axis passing by the middle point M = (A + B)/2 of ordinate z M . Its equation is given by For simplicity, the width of the parabola will be referred as h = 1 8 C(z A − z B ) 2 , and its height L = z A − z B .

Closed-form Point Spread Function
In this subsection, we derive a closed-form solution for the Fourier transform of each elementary shape. Indeed, it would be possible to do it numerically with a discrete Fourier transform algorithm of the pupil function. The issue is that the sampling needed, the zero-padding for periodicity breaking, and the high dynamic range of the PSF, would make this algorithm very memory consuming and inaccurate. Furthermore, in a Fourier transform, any resulting value depends on the whole range of input values. For all these reasons, a closed-form solution would give each value of the Fourier transform accurately, without a need for extra memory.
In order to have a closed-form solution of the PSF equation (1), one has to derive the Fourier transform of the basis shapes: triangle and parabola. The shape function of the triangle is then referred as P tri , and the parabola as P par . The two closed-form Fourier transforms are given by Step 1 (resp. 2) aims to set the front ring (resp. diaphragm) of the lens in focus, giving the l 1 (resp. l 2 ) distance between the front ring and the mirror. These measurements then allow us to retrieve the focal plane distance d, the shift ∆ between the aperture and the front of the camera, and the diaphragm parameters with a picture taken at Step 2.
with η A/B = π z A/B ν y + r AB ν x , and γ tri = − r AB 2iπν y ; and where Thanks to the closed-form solutions for the elementary shapes and due to the linearity of the operator, the complete Fourier transform of a n-bladed diaphragm is given by: with k the index of the k-th sub-shape, and R(α k ) the rotation matrix of angle α k . Finally, to compute the PSF (cf. Eq. 1), the diaphragm area S pup is required:

Measurement of Camera Parameters
A complete description of our camera model with the shape of its diaphragm needs a lot of parameters. In this subsection, we present a very simple calibration procedure to retrieve these parameters. The goal of our calibration procedure (cf. Fig. 2) is to determine each camera parameter: the focal plane distance d * (which replaces d in the equations due to condition (3)), the sensor distance D, and the pupil settings: the curvature C and the diaphragm vertices. The exact position of the diaphragm within the lens is not known and is not accessible from a direct measurement. To address this issue, a real object can be taken as reference, in this case the front ring of the camera. A mirror is placed in front of the camera in order to allow the camera to image itself. Moreover, this self-imaging technique has the advantage that a good alignment on the optical axis can be granted.
The procedure then consists in measuring the two distances l 1 and l 2 between the front ring and the mirror in two steps where the front ring of the camera, then the diaphragm is in focus, respectively (cf. Fig. 2, steps 1 and 2). To have a more precise measurement of these two distances, a good approach is to let the diaphragm wide open in order to minimize the depth of field. Then, while the targeted object crosses the focal plane, one can look for the minimum autocorrelation width in order to have a better estimation of the focal plane position. Therefore, simple geometry rules give the following results where ∆ is the distance between the front ring and the diaphragm. With these parameters, the D = d f /(d − f ) sensor distance to the diaphragm can be deduced. In order to find the pupil parameters, we use the image taken when the diaphragm is in focus (cf. Fig. 2, step 2). From this picture, by using the knowledge of the image pixel size (d/D)d x with d x the sensor pixel size, we fit the diaphragm and measure the C curvature and the positions of the vertices.
A good way to check for the consistency of the diaphragm parameters is to compare the indicated f-number with its measured equivalent f # depending on S pup (cf. Eq. 9)

High Dynamic Range Imaging Technique
In imagery, the scene I to acquire has its values ranged within [min(I), max(I)]. The ratio of the two boundaries is called the "dynamic range" of the scene, noted D I and therefore defined by The problem is that the scene dynamic range can have an arbitrary value whereas the camera has a limited one, noted D c . Therefore, we have in general D c < D I ; consequently, taking a picture with a camera corresponds to extracting a band of intensity values from the scene (cf. Fig. 3, left).
This band of values, noted b m 1 , can be chosen by setting the level of exposure of the camera via the exposure time for instance). In any case, pixels of such a picture fall into 3 categories: 1. "underexposed pixels", which are pixels with a value below the minimum camera pixel value, they appear black 2. "overexposed pixels", which are pixels with a value beyond the maximum camera pixel value, they appear bright because of the sensor saturation 3. "well-exposed pixels" are the pixels falling within the camera pixel range of value, they are the only metrologically reliable pixels of the picture.
This means that the underexposed and overexposed pixels are discarded from the measurement since their true values are unknown. To address this issue, the common technique is to use High Dynamic Range (HDR) imaging, the goal of which is to increase the effective dynamic range of the camera. The idea is to take several pictures at different levels of exposure such that the union of each individual band permits to recover the full scene dynamic range (cf. Fig. 3, right). Principle of a picture of a scene with a higher dynamic range than the camera. The level of exposure of the camera sets a band of well-exposed values within the whole scene dynamic range. The under-and over-exposed pixels are discarded from the measurement. Center: HDR imaging principle of a scene. In order to increase the dynamic range of a single picture, multiple pictures can be merged to cover the whole intensity distribution of the scene. Right: Cut of the HDR image values into separate non-overlapping bands of value b k such that the whole image dynamic range is covered.

Diffraction and High Dynamic Range Images
The measured HDR image I hdr is described by a convolution of the perfect original image I * hdr with the PSF (cf. Eq. 4), but with an additive noise term B, so that The most general case is the one with images of very high dynamic range, the effective spatial extent of the PSF (noted PSF radius ) can cover up to the entire image, such that every pixel can be affected by the value of any other. In practice, it means that if the original image size is (n x , n y ), the PSF must be given in an image of size (2n x − 1, 2n y − 1). If the measured image I hdr is taken as the final retained measurement, without any post-computation, the uncertainty on each value is given by the noise B, but also up to a convolution kernel of very big radius (depending on the image dynamic range, up to hundreds of pixels extent) implying a great loss of effective resolution. Therefore, the I hdr image should not be taken as a reliable measurement of I * hdr because of this kernel that makes every pixel of the measurement interdependent with one another.

Diffraction Detection Algorithm
In this section, we present our algorithm, the goal of which is to identify pixels polluted by diffraction.

Algorithm Overview
Our analytical PSF function permits to predict the effects of diffraction. Based on this known PSF, our algorithm simulates a second diffraction on the acquired image (the perfect image is then diffracted once by the physical diaphragm, then through simulation). Our method relies on two ideas: (i) if a pixel is not modified during our simulated diffraction, it was not the case during the physical diffraction, either; and (ii) diffraction pollution on a pixel is always originating from pixels of higher values. Even though these assumptions are not true in general, they become valid if we accept a residual diffraction kernel. The idea of this residual kernel is that, within a certain extent, diffraction makes pixels too interdependent, so that our method cannot separate the effects of diffraction from the true pixel values.
Following these considerations, our diffraction detection algorithm is divided into three parts: 1. The HDR image is cut into non-overlapping bands of values of same dynamic range (cf. Fig. 3).
2. A residual convolution kernel K is removed from the diffraction prediction (cf. Algo. 2).
3. Diffraction is progressively predicted, by iterating from the band of highest values toward the lowest and applying a user thresholding criterion to discard pixels affected by diffraction (cf. Algo. 1).
At first (cf. Fig. 3), the HDR image is cut into non-overlapping bands of value (b 1 , .., b k , ..b N ) of identical dynamic range D b . A binary mask function 1 k describes the domain where I hdr lies in the b k band. Each sub-picture is then referred as I k , noting that I k = I hdr × 1 k . For multiple bands, from k 1 to k 2 , the quantities are subscripted The key idea is that for most lenses, the dynamic range in which the PSF is very similar to a Dirac function is big, between a factor of 10 to 1000. Each sub-picture I k is therefore composed of two separate contributions: its inner value I * k that is considered diffractionfree and a diffraction term coming from the higher bands. The exact definition and implication of this "diffraction-freeness" is discussed in Subsection 6.3.
Our algorithm (cf. Algo. 1) essentially consists of a sequence through the bands, from the highest (b 1 ) to the lowest (b N ). In each iteration, a partial HDR image I 1→k−1 is convolved with the PSF, and the diffracted values present in the 1 k mask are extracted. These values are compared to the original I k picture, and a thresholding criterion ρ is applied to distinguish clean pixels from the ones affected by diffraction. This method is then iteratively applied until the full image dynamic range has been covered.

Conditions on the HDR image.
For the algorithm to predict the effect of diffraction on an image correctly, two conditions are required during the HDR image measurement: (i) an overlap exists for consecutive exposure bands, and (ii) the highest band must not have any over-exposed pixels. For instance, these conditions are met in Figure 3 (right part): each band presents an overlap with its neighboring bands, and the highest band (red) contains the highest pixel value. Our algorithm can only be applied on input HDR images respecting these conditions.

Core Algorithm
Firstly, the HDR image is cut into non-overlapping bands of values b k (cf. Fig. 3). Without loss of generality, I hdr can be normalized such that its maximum becomes 1 (cf. Algo. 1, line 2). Therefore, the band cut (cf. Algo. 1, lines 3 & 6-9) is defined by: As stated previously, the I hdr HDR image is already subject to diffraction effects, since it has been measured. Yet we propose to numerically diffract it a second time, by computing I hdr ⊗ PSF (cf. Algo. 1, line 10). From this computation, the method is based for k ← 2, N do 6: I k ← I hdr * 1 k 9: I 1→k−1 ← I hdr * 1 1→k−1 10: end for 13: return Discarded, K 14: end procedure on the following principle: if a pixel is unchanged from I hdr to I hdr ⊗ PSF, it also remains unaltered from I * hdr to I hdr . This principle is applied band by band, from the highest b 1 to the lowest b N .
In order to justify the principle and find the conditions of validity, we can derive the effect of diffraction of the HDR image over one band b k . Indeed, this operation is formally described as follows Since the values of the lower bands are smaller than the current k-th band and the PSF is rapidly decreasing, the third term is considered negligible. In fact, this assumption is not always true, yet it is only valid under the bottom-up influence condition (cf. Subsection 6.3.2). Therefore, the previous equation can be simplified to: Furthermore, as described before, if the dynamic range of a band D b is small enough, the PSF acts as a Dirac function. Here again, this assumption is not generally true, but it is valid under the within-band influence condition (discussed in Subsection 6.3.1). Then, for the diffraction of the I k value, one can neglect the non-Dirac term of the PSF such as and normalize it so that the approximation remains energy conservative. Therefore, equation (15) can be simplified to: Finally, from equation (17), computing I hdr ⊗ PSF essentially consists of the sum of two terms: I k , which is the inner value without convolution, and the other term, which is the diffraction from the upper bands impacting the I k picture. Therefore, the pixels where this diffraction term is negligible are considered "diffraction-free". A threshold criterion is then chosen as input of the algorithm, noted ρ, which defines that any pixel affected by diffraction verifies Consequently, for each band, this threshold determines the pixels to discard from the measurement. For RGB color images, the algorithm can be applied separately on each color channel. The PSFs have to be recomputed for each color because of the different wavelengths λ. Generally, the scene illumination spectrum is not known, thus a good estimate of the λ value is the maximum spectral sensitivity per channel. In this band of values, our algorithm principle states that if a pixel is unchanged from I hdr to I hdr ⊗ PSF, then it is the same case from I * hdr to I hdr . For this principle to be applicable, two conditions are required: (i) diffraction effects are negligible within a single b k band (within-band influence), and (ii) the b k band is not affected by diffraction coming from the lower bands (bottom-up influence). However, these assumptions are not true in general, and this enables to quantify to what extent the algorithm is capable of detecting diffraction.
To this end, we have to define a kernel K (cf. Algo. 2) in which interdependencies between pixels are too strong, such that our algorithm cannot separate the inner value from the diffraction contribution. This kernel comes from the two conditions described previously, and is thus defined as their combination with K wb and K bu two binary functions (cf. Eqs. 22 and 25).
Since we cannot sort out pixels that are so strongly interdependent, when predicting the influence of diffraction on the I k picture, it is necessary to remove K from the prediction. Concretely, the sorting condition (18) becomes more flexible, Fig. 4: Left: Within-band influence effect. In this worst-case scenario, pixels within a band can be linked through diffraction, while we assume this is not the case. Thus, the effect of diffraction can be removed up to a residual convolution kernel K wb . Right: Bottom-up influence effect. In this worst-case scenario, pixels from the lower bands should never be able to diffract more than ρ % of the pixel values in the current band. Therefore, a band is strongly interdependent with lower bands by a residual kernel K bu .
After executing the algorithm, the remaining pixels can be characterized with their uncertainty from noise B but also up to a residual convolution kernel given by the function K . Therefore, the remaining (i.e., non-discarded) pixels I output are metrologically characterized by Even if the function K can be arbitrarily shaped, it can be characterized by its maximum outer radius, thereafter noted K radius .
Our algorithm aims to be conservative regarding equation (21): any pixel that does not fit in this equation is rejected. However, a lot of pixels that fit it can be discarded by our algorithm, implying that many more pixels than intended are lost.

Within-band Influence
Let us consider the two extreme cases of an I k picture for the effect of convolution: I k is a constant function, and I k is a Dirac function. In the constant case, the PSF being normalized, it is easy to conclude that a convolution by the PSF does not affect I k . So, if I k is a constant, the diffraction effect is always negligible. In the case of a Dirac function (cf. Fig. 4, left), the I k ⊗ PSF becomes the PSF function itself. So when the band of value b k needs to be considered "diffraction-free", a small convolution kernel still remains. This remaining kernel then have to be removed from the diffraction prediction since our method cannot separate strongly related pixels. Therefore, the following binary mask defines the inseparable diffraction kernel, caused by this within-band influence condition.

Bottom-up Influence
In order to check for the validity of neglecting the contribution of I k+1→N ⊗ PSF over the I k term, let us consider the I k picture in the worst-case scenario (cf. Fig. 4, right): I k is composed of a single pixel of value v − + η, and I k+1 an image of constant value v − − η (other than the single I k pixel), with η > 0 an infinitesimal quantity.
In the limit η → 0, this situation describes a constant image of value v − , where one pixel is considered in the b k band, and all the others in the b k−1 band. Our method is supposed to discard a pixel if it predicts a relative amount of the diffraction contribution superior to ρ. In this situation, referring to the rejection condition (18), diffraction would be neglected if Indeed, in this situation, since the pixels are of equal intensity, this condition may not always be satisfied, and this effect is noted as bottom-up influence. The solution is, as for the within-band influence condition, to consider that our algorithm is not capable of separating the diffraction effect in this worst-case situation.
Hereafter, a residual kernel has to be accepted and removed from the diffraction prediction. This kernel is defined such that if we remove this residual kernel from the PSF function, condition (23) has to be respected. Among multiple possible solutions, we chose to keep the one that minimizes the area of this residual kernel. This method is to find the threshold s * to the PSF that best fits condition (23) : With this best threshold, the following binary mask defines the second inseparable diffraction kernel, caused by this bottom-up influence condition. Fig. 6: Comparison of the PSF resulting from the fitted diaphragm against a real HDR photograph of a quasi-point light source. Some slight differences can be observed in the repartition of light within the widened star branches of the PSF, which is explained by the random variations along the diaphragm edges that we do not take into account.

Results
In this section we present our results for the PSF computation as well as the diffraction detection algorithm. Firstly, comparisons between analytical PSF and the true recorded PSF with a fitted diaphragm are provided, showing the high accuracy of our model. Secondly, the algorithm is applied on real case scenarios, providing a good understanding on the different limitations. Finally, the algorithm is applied on simulated HDR images. The contribution of diffraction is natively known for each pixel, that makes it possible to assess the efficiency of our algorithm to separate highly polluted pixels and compare it to state-of-the-art deconvolution techniques.

Real Aperture Fitting and Point Spread Function
The aperture model composed of an irregular polygon with curved edges is assessed to be general enough to cover a wide range of camera lenses. We tested it on our available camera lenses: one scientific-class lens of focal 50mm from Linos and two consumer Canon lenses of 50 and 100mm focal length. The goal is to compare how the diaphragm model fits a real aperture and to demonstrate that the resulting theoretical PSF also fits well a true PSF image. The variety of diaphragms in Figure 5 highlights the need to have an elaborated enough mathematical modeling. Our model allows a very good fit of a wide range of common diaphragms and its Fourier transform is analytical, as well as the resulting PSF. As shown in Figure 5, the irregular polygon and the curved edges features have their importance. For the Canon 100mm lens at f/11, it is sufficient to fit an irregular polygonal shape, with no need for a curvature term. In contrary, the Linos 50mm at f/4 could not have been described with a regular polygon, as the curvature of the edges really needs to be taken into account.
Even if the aperture is well fit by our diaphragm model, the theory driving the PSF also needs to fit well a real photograph. Concerning the resulting PSF, the simulation is compared to a real PSF measured in HDR as shown in Figure 6. In our case, the diaphragm fit is man-made, so the PSF is subject to some uncertainty. A certain roughness is present on the edges of the diaphragm, which is not taken into account in our model. Due to Fourier transform properties, roughness has the effect to change the light repartition by widening the star branches of the PSF. This effect is clearly visible in Figure 6, in the bottom left star stripe of the Canon 100mm PSF, where the prediction under-estimates the widening of the stripe. In order to include this effect, the distribution of normals of each edge is needed, which is far off the camera resolution available Fig. 7: Results of the algorithm applied on real HDR images (tonemapped with Drago et al. [1]) for various camera configurations, with input parameters D b = 10 and ρ = 5%. The wavelengths used for each color channel are [λ R , λ G , λ B ] = [600 nm, 540 nm, 470 nm]. The segmentation images show the discarded pixels (red), the valid pixels (green), and the under-exposed ones (black). If the HDR images exhibits obvious star shaped patterns, the algorithm detects it, and they are finally removed. Such result is qualitative in nature, because there is no reference HDR image without diffraction. False predictions are present in the first two cases (l), where the diffraction prediction seems rotated from the real one. This problem emerges from the misfit of the lens diaphragm, as discussed in subsection 7.7.1.
during the diaphragm fitting.
However, the main problem comes from the fact that for strongly closed down apertures the diaphragm shape is not repeatable (depending on the quality of the diaphragm). In fact, we observed that the polygon of the diaphragm seems to be rotated by a few degrees. This effect is directly emerging from the fact that closing down an aperture essentially consists in reducing the size of the polygon while rotating it. If the diaphragm is not in the exact same configuration for each user setting of the fnumber, we mainly observe that the diaphragm shape obtained is a rotated version of the measured one. For instance, with our Canon 100mm at f/27, it results in a 3°tilt deviation from our prediction.
As a consequence, since a good description of the PSF function implies a good repeatability of the diaphragm closure, our method is more suitable for scientific-grade cameras, as well as for fixed and toothed manual apertures.

Diffraction Prediction in Real Case Scenarios
Using the same camera lenses as described previously, HDR images have been taken in laboratory but also in real uncontrolled conditions (night pictures).
The algorithm seems to discard a lot more pixels than one would expect, highlighting the fact that the method does not pretend to discard only pixels affected by diffrac- Fig. 8: Histograms of the error of magnitude against a virtual reference of the remaining valid pixels for various methods and three different SNRs. The PSF function used is given by our Linos 50mm lens at f/11. The E max factor measures the maximum error remaining after applying our method (red curve). The resulting histogram is much more concentrated towards smaller errors than compared to all deconvolution algorithms (blue curves). Of course, the quality of the original image (green curve) is not reached because of the residual kernel contribution, but our output error matches very well with the achieved output (brown curve) prediction.
tion, but also diffraction-free pixels. Since the algorithm can be too conservative, the percentage of discarded measurements can significantly decrease the efficiency of an HDR image-based measurement. The K kernel is also much smaller than the PSF kernel, falling into a range of few pixels, which guarantees that the long-range blurring effect of the PSF has been removed.
In laboratory conditions, where we used our Linos lens, the scene is perfectly stable and controlled, and the camera response is also very stationary. In this situation, shown in Figure 7, our diffraction removal algorithm completely removes the widened star shaped pattern making it very useful for measurements. In an uncontrolled scenario (e.g., with outdoor imaging) the illumination conditions are not stable wrt. time, and HDR values can be shifted up or down because of the intensity variation of lamps. Moreover, as shown for our Canon lens (cf. Subsection 7.1), the diaphragm fitting can be incorrect because of the lack of repeatability of the lens diaphragm setting. Then, the PSF prediction is biased, so are the discarded pixels. This is visible on the two left cases of Figure 7, where the removed pixels seem tilted with respect to the star shaped pattern.

Error Analysis
A good way to quantify the quality of the separation between polluted and non-polluted pixels by diffraction is to test the algorithm on a great variety of generated HDR images. Given one image, its "real" measurement is simulated by convolving it with the precomputed PSF and by adding an additive Gaussian white noise. Our algorithm is then applied to this resulting image.
In order to remain as general as possible, our HDR test images are tuned by their bandwidth limit (Gaussian speckle pattern), their histogram of magnitude, and their HDR dynamic range (D hdr ). It is possible to generate a wide variety of such images. Since the different features and conclusions do not seem to be altered whatever the input image, by default, the chosen generated image is a HDR image with a flat histogram, D hdr = 10 10 and a speckle size of 20 pixels.
Since our method focuses on guaranteeing no diffraction pollution on the remaining pixels, the data of interest is the histogram of relative errors between the "true" image and the "measured" one. One relevant metric to be used is the "maximum error of Fig. 9: Variation of the residual kernel range K radius (in pixels) depending on the input parameters. According to its definition, the residual kernel size arises from two contributions which are easily separable (green curve): when one effect is dominant, the other effect does not interfere with it.
This metric allows sorting the different methods, comparing our method to the ones from the state of the art. In Figure 8 are plotted relative histograms of the E error for various SNR values. The PSF used to simulate a measurement is that of the 50mm Linos lens at f/11 and the noise is a Gaussian white noise, which power is given by a signal-to-noise ratio (SNR).
Since convolution problems depend on the image frequency content, the algorithm has been tested on different SNR values and generated images: high values, well-ranged values, low values, large sized speckle and small sized speckle. The conclusion does not depend on the image content: the maximum error E max resulting from our algorithm (with D b = 10 and ρ = 5%) is always better than any other tested deconvolution method (cf. Fig. 8, blue curves), and the result histogram (red curve) fits very well what we expect to recover (cf. Eq. 21), which corresponds to a measurement quality up to a K kernel convolution (brown curve). Figure 8 also makes it evident that not considering diffraction may lead to a very inaccurate measurement: the quality of the ground truth (green curve) is far off the real initial measurement (black curve).
The residual kernel K represents the incapacity of the algorithm to separate diffraction interdependence between certain pixels. Its relationship to the input parameters is given by its definition (cf. Subsection 6.3), therefore we observe the two following phenomena (cf. Fig. 9): 1. when the within-band influence effect is predominant, the size of K only increases with the D b input parameter 2. when the bottom-up influence effect is predominant, the size of K only increases with the thresholding criterion ρ input parameter.
Thus, when one parameter shapes the residual kernel, the other's parameter variation tunes the amount of discarded pixels, along with the quality of the remaining ones (E max ). However, the relationship between E max and the input parameters is not analytically known. Therefore, from a measurement perspective, there are two ways of characterizing the algorithm output: Fig. 10: The output dependencies of the algorithm on the input parameters. The generated input HDR image is a well-distributed HDR image with a dynamic range of 10 10 and a speckle size of 20 pixels. In the region of minimal maximum error (dashed green square), the extent of the residual kernel and the percentage of discarded pixels go in opposite trends.
1. According to Figure 8, the output fits well the analytical prediction of equation (21). Thus, it is possible to characterize the output by: a spatial resolution given by the convolution with K, and the values are given with an uncertainty that directly equals to B.
2. Secondly, it is also possible to omit the residual kernel and to characterize directly the remaining pixels values : the maximum remaining error, E max , is taken as the global relative error on all the output pixels. Yet the E max is not known since the reference image without diffraction is also unknown. In this scenario, the only proposed solution consists in creating an image with similar content to determine a good estimate of the E max value.
Despite the good characteristics of the algorithm output, an inconvenient issue is that for one measurement the loss in terms of number of pixels can be huge (up to 95% of the whole image). The ρ parameter may be too strict and discard too many pixels compared to the user tolerance. However, even if we cannot mathematically link the E max value with the ρ criterion, it is possible to understand the existing trade-off between the input parameters (ρ and D b ) and the result of our algorithm (E max , the percentage of discarded pixels and K).
The three graphs in Figure 10 map the evolution of the outputs of the algorithm with the two inputs. This figure, computed in the generic case previously described, exhibits general features that are present in every test case. The most important feature is that an area of minimal E max error exists in the input space (cf. Fig. 10, green dashed box). The existence of such area stems logically from the definition of the residual kernel: • if ρ is too high, too many pixels are accepted while their diffraction amount is high; • if ρ is too low, the bottom-up influence implies that the K kernel becomes wider and wider, thus the pixels within the K kernel are accepted even though they can be highly affected by diffraction; • as D b increases, the within-band influence also forces the algorithm to accept more and more pixels.
Within this minimal error area, as required by the user, some flexibility exists to choose the best option between minimizing the number of rejected pixels, and minimizing the residual kernel radius.
Performance. We tested the performance of our algorithm by implementing it using Matlab with the GPU toolbox on a computer with a CPU Intel ® Core ™ i7-4790K CPU @ 4.00GHz and an Nvidia ® GPU GeForce GTX 980 Ti. For a 1000×1000 grayscale image, setting the 6-bladed Linos lens camera parameters to: d = f = 50mm, f # = 11, pixel size of 6.45×6.45 µm 2 , and λ = 555nm, the timings are: 202s to precompute the PSF, 0.8s to compute the residual kernel K and 2s for the main iterative algorithm (7 bands). The PSF precomputation time is long because of the over-sampling needed to respect the Shannon criterion as well as the calls to erfi functions that are very time consuming. This could be further optimized by using more intensively the GPU.

Conclusion and Future Work
Since it is not possible to recover values of a HDR measurement polluted by diffraction within uncertainty boundaries, our algorithm focuses on separating between pixels that are affected and unaffected by diffraction, respectively. Our algorithm exploits that diffraction mainly implies high value pixels affecting low value pixels. Predicting diffraction effects needs a good fitting of the diaphragm which is provided by our model, however a bad repeatability of aperture closing may lead to inaccuracies. After applying our algorithm, the remaining "clean" pixels are not modified, their uncertainties are then those given by a direct calibration. The resulting convolution kernel is also greatly reduced, so is the effective spatial image resolution. The result of the algorithm ensures a good quality of the measurement, yet the link between the algorithm parameters and the resulting image characteristics is not known, despite clues on their dependence.
As future work, we intend to focus on the precise analysis of the impact of the input image on the result. The histogram, the frequency content and the spatial coherence of the HDR image should give more insight on how to predict well the resulting error from any measurement; at the moment we still have to infer it from a generated contentequivalent image. The PSF model can also be improved, by improving the diaphragm edge description. In particular, a roughness term may be added for the edges, a method that could be inspired from the prediction of radio wave propagation above rough landscapes [28].