change on histogram for extern value dot affect the definition

[CharlesAuthier]

Context

A problem comes out, during the emva test, with the logarithmic histogram when we had "dead pixel" or "defect pixel", the histogram became useless. The graph lose is resolution and the model doesn't fit the distribution anymore. Here an image generated with nine defects pixels when generated the emva test.

Zoom on the distribution.

The origine of this came form emva1288/emva1288/process/routines.py in the Histogram1288 function. When the image contains defects pixels, that change the maximum and/or minimum values($y_{max}$ and $y_{min}$ of the image $y$). Resulting in this situation the Q(number of bins) and W(binwidth) will change. The model who is determinate with the mean and the standard deviation of the distribution appears to be inaccurate due to the presence of extremes values.

y = np.ravel(img)
ymin = np.min(y)
ymax = np.max(y)

# Because we are working with integers, minimum binwidth is 1
W = 1
q = ymax - ymin
Q = q + 1

The standard already put a $Q_{max}$ at 256 bins so when the image has extremes values like the first picture, the binwidth change and became larger reducing a low resolution for the histogram like the second picture show.

# When too many bins, create a new integer binwidth
if Q > Qmax:
    # We want the number of bins as close as possible to Qmax (256)
    W = int(np.ceil(1. * q / (Qmax - 1)))
    Q = int(np.floor(1. * q / W)) + 1

This pull request has for goal to have the same resolution than if the camera doesn't have defect pixel, regain a model who fit the distribution and without losing the values of the defect pixel.

Solution

For this case, nine defects pixels are simulated for a camera in camera.py. To do that a value of 255 is affected at nine pixels in the DSNU. It's simple to simulate defect pixel in the DSNU because is a kind of nonuniformity and is an addition of the dark signal who is directly adding to the image of electrons.

if dsnu is None:
    self._dsnu = np.zeros(self._shape)
    self._dsnu[0, 0] = 255
    self._dsnu[300, 600] = 255
    self._dsnu[400, 460] = 255
    self._dsnu[460, 300] = 255
    self._dsnu[100, 120] = 255
    self._dsnu[300, 90] = 255
    self._dsnu[20, 50] = 255
    self._dsnu[250, 250] = 255
    self._dsnu[450, 450] = 255

This simulation gives a logarithmic histogram DSNU represented the two first picture.

Starting with the fact of the model resulting from the distribution need to be a normal distribution and 99.9999426697% of the values in a normal distribution are cover by five sigmas. So all values outside the deviation of five sigmas must be a defect pixel (extremes outliers).

To compute the distribution and not the useless outliers, the mean of the image $y$ and the standard deviation are calculated.

# Full distribution sigma, and median
full_median = np.mean(y)
full_sigma = np.std(y)
# with 5 sigma, 99.9999426697% values are in the Normal distribution
# Source: Wikipedia, Normal distribution, section on standart diviation
sigma5 = full_sigma * 5
# Core, main distribution
core = y[(y > full_mean - sigma5) & (y < full_mean + sigma5)]

The region cover by 5 sigmas will be called core. This core has a new minimum and maximum who doesn't include the extremes outliers. A mean and a new sigma are calculated for the model who comes later. In the case where they have no outliers, nothing will change, because the core will cover all the image.

# Normal distribution with the original sigma, and mean
    mu = np.mean(core)
    sigma = np.std(core)

    # Core ymin and ymax
    ymin_core = np.min(core)
    ymax_core = np.max(core)

To change the resolution only a new binwidth is required. The covered region needs to be the same(with the outliers) so the $Q$ will upgrade and probably by over the $Q_{max}$.

# Because we are working with integers, minimum binwidth is 1
W = 1
q = ymax - ymin
# Same for the Core
q_core = ymax_core - ymin_core
Q_core = q_core + 1

# When too many bins, create a new integer binwidth
if Q_core > Qmax:
    # We want the number of bins as close as possible to Qmax (256)
    W = int(np.ceil(1. * q_core / (Qmax - 1)))

Q = int(np.floor(1. * q / W)) + 1

The result of this implementation:

Zoom on the distribution.

Zoom on the defects pixels.

It's hard to see the outliers with the modification because the resolution increases, but they are still there at different values of deviation of the mean with a number of pixels of one.

Conclusion

More complex code, increase the time to process but not very significative. In result, the histogram with a better definition and a model who can be used.

Change to the 1288 EMVA Standard Compliant

These changes require some modifications in the 1288 EMVA Standard Compliant. Based on the release 3.1 of December 30 2016, this are the suggested changes:

Section 4.3.1 Logarithmic Histograms.

It is useful to plot the histograms with logarithmic y-axis for two reasons (Fig. 2a). Firstly, it is easy to compare the measured histograms with a normal distribution, witch shows up as a negatively shaped parabola in a logarithmic plot. Thus it is easy to see deviations from normal distributions. Secondly, also rare outliers, i.e., a few pixels out of millions of pixels can be seen easily.
In the case of a death pixel, it is generally when a pixel is saturated for no reason. They can be considered like extreme outliers because they are generally very far from the mean. They cannot really be considered in the deviation, they didn't give any useful statistical information. Setting a limit for the standard deviation can be useful, all outliers outside of the limit (example five standard deviation) are death pixel or useless pixel. Like that, the only outliers, who are in the limit, provide statistical information to the distribution and the model of normal distribution.
All histograms have to be computed from pixel values that come from averaging over many images. In this way the histograms only reflect the statistics of the spatial noise and the temporal noise is averaged out. The statistics from a single image is different. It contains the total noise, i.e. the spatial and the temporal noise. It is, however, useful to see in how far the outliers of the averaged image histogram will vanish in the temporal noise (Fig 2b).
It is hard to generally predict in how far a deviation from the model will impact the final applications. Some of them will have human spectators, while others use a variety of algorithms to make use the images. While a human spectator is usually able to work well with pictures in which some pixel show odd behaviors, some algorithms may suffer from it. Some applications will require defect-free images, some will tolerate some outliers, while other still have problems with a large number of pixels slightly deviating. All this information can be read out of the logarithmic histograms.

Section 8.4 Defect Pixel Characterization

Before evaluation of the PRNU image $<y_{50}> - <y_{dark}>$, this image must be highpass-filtered as detailed in Section 4.4 and Appendix C.5. No highpass-filtering is performed with the DSNU image $<y_{dark}>$.
The computation of the logarithmic histogram involves the following computing steps for every averaged and highpass-filtered nonuniformity image $y$. The following procedure is suggested to obtain optimally smooth histograms. It is based on the simple fact that nonuniformity images are averaged over $L$ integer-valued images. Therefore they there values are integer multiples of $1/L$ and the optimum interval width is therefore also an integer multiple of $L$.

1. Compute minimum and maximum values $y_{min}$ and $y_{max}$ of the image $y$.
2. Compute the standard variation $s_{nw}$ of the image $y$ and multiplies by the number of time to have a good cover of the data. An appropriate one will be five times to reach 99.9999426697%[make a reference to the wikipedia page].
3. Part of the $y$ image can be extremes outliers like death pixel. To not interfere with the distribution.
   $core = (\mu_{y} - 5\dot s_{nw}) < y < (\mu_{y} + 5\dot s_{nw})$
   The $\mu_{y}$ is not to far from the real mean of the distribution, because of the large number of good pixels compare to the few defects pixels and in the case where they have no outliers the image $y$ will become the core.
4. Compute a new minimum and maximum values $y_{min.core}$ and $y_{max.core}$ of the image $core$ without the extremes outliers
5. Part the interval between $y_{min}$ and $y_{max}$ into $Q = L(y_{max} - y_{min})/I +1$ bins of equal width with the optimal bin width of $I/L$. Choose an appropriate ($I = 1,2,3,...$) so that the number of bins is lower than or equal to 256. This condition is met by
   $I = floor[\frac{L(y_{max} - y_{min})}{256}] + 1.$
   It is optimal except if the image has extremes outliers, in that case, the outliers can be too far the histogram will lose in resolution because of an $I$ to hight. Changing $y_{min}$ and $y_{max}$ to $y_{min.core}$ and $y_{max.core}$ in the equation (number of the equation just above) will set the optimal bin width for the distribution. Therefore the $Q$ needs to stay the same to cover even the extremes outliers so the number of bins can exceed 256 in that case.
6. Compute a histogram with all values of the image using these bins. The bin $q$ to be incremented for a value y is
   $q = floor[\frac{L(y - y_{min})}{I}]$
   and for the core distribution for a value $core$ is
   $q_{core} = floor[\frac{L(core - y_{min.core})}{I}].$
   In this way the $Q$ bins of the histogram (indices from 0 to $Q$ - 1) cover values from $y_{min}$ to $y_{max} + (I - 1)/L$.
   The value of the center of the bins as a deviation from the mean value are given as:
   $core[q_{core}] = y_{min.core} + \frac{I - 1}{2L} + q_{core}\frac{I}{L}$.
7. Draw the histogram in a semilogarithmic plot. Use a $x$-axis with the values of the bins relative to the mean value. The $y$-axis must start below 1 so that single outlying pixels can be observed.
8. Add the normal probability density distribution corresponding to the non-white variance $s_{nw}^{2}$ as a dashed line to the graph obtained for an $M \cross N$ image with an interval width $I/L$:
   $p_{normal}[q_{core}] = \frac{I}{L}\dot\frac{NM}{\sqrt{2\pi}s_{nw}}\dot exp(-\frac{core[q_{core}]^{2}}{2 s_{nw}^{2}})$.

The accumulated logarithmic histogram [...].