1# Copyright 2014 The Android Open Source Project
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7#      http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
14
15import its.device
16import its.caps
17import its.objects
18import its.image
19import os.path
20from matplotlib import pylab
21import matplotlib
22import matplotlib.pyplot as plt
23import math
24import textwrap
25import time
26import numpy as np
27import scipy.stats
28import scipy.signal
29
30
31# Convert a 2D array a to a 4D array with dimensions [tile_size,
32# tile_size, row, col] where row, col are tile indices.
33def tile(a, tile_size):
34    tile_rows, tile_cols = a.shape[0]/tile_size, a.shape[1]/tile_size
35    a = a.reshape([tile_rows, tile_size, tile_cols, tile_size])
36    a = a.transpose([1, 3, 0, 2])
37    return a
38
39
40def main():
41    """Capture a set of raw images with increasing gains and measure the noise.
42    """
43    NAME = os.path.basename(__file__).split(".")[0]
44
45    # How many sensitivities per stop to sample.
46    steps_per_stop = 2
47    # How large of tiles to use to compute mean/variance.
48    tile_size = 64
49    # Exposure bracketing range in stops
50    bracket_stops = 4
51    # How high to allow the mean of the tiles to go.
52    max_signal_level = 0.5
53    # Colors used for plotting the data for each exposure.
54    colors = 'rygcbm'
55
56    # Define a first order high pass filter to eliminate low frequency
57    # signal content when computing variance.
58    f = np.array([-1, 1]).astype('float32')
59    # Make it a higher order filter by convolving the first order
60    # filter with itself a few times.
61    f = np.convolve(f, f)
62    f = np.convolve(f, f)
63
64    # Compute the normalization of the filter to preserve noise
65    # power. Let a be the normalization factor we're looking for, and
66    # Let X and X' be the random variables representing the noise
67    # before and after filtering, respectively. First, compute
68    # Var[a*X']:
69    #
70    #   Var[a*X'] = a^2*Var[X*f_0 + X*f_1 + ... + X*f_N-1]
71    #             = a^2*(f_0^2*Var[X] + f_1^2*Var[X] + ... + (f_N-1)^2*Var[X])
72    #             = sum(f_i^2)*a^2*Var[X]
73    #
74    # We want Var[a*X'] to be equal to Var[X]:
75    #
76    #    sum(f_i^2)*a^2*Var[X] = Var[X] -> a = sqrt(1/sum(f_i^2))
77    #
78    # We can just bake this normalization factor into the high pass
79    # filter kernel.
80    f /= math.sqrt(np.dot(f, f))
81
82    bracket_factor = math.pow(2, bracket_stops)
83
84    with its.device.ItsSession() as cam:
85        props = cam.get_camera_properties()
86
87        # Get basic properties we need.
88        sens_min, sens_max = props['android.sensor.info.sensitivityRange']
89        sens_max_analog = props['android.sensor.maxAnalogSensitivity']
90        white_level = props['android.sensor.info.whiteLevel']
91
92        print "Sensitivity range: [%f, %f]" % (sens_min, sens_max)
93        print "Max analog sensitivity: %f" % (sens_max_analog)
94
95        # Do AE to get a rough idea of where we are.
96        s_ae, e_ae, _, _, _  = \
97            cam.do_3a(get_results=True, do_awb=False, do_af=False)
98        # Underexpose to get more data for low signal levels.
99        auto_e = s_ae*e_ae/bracket_factor
100
101        # If the auto-exposure result is too bright for the highest
102        # sensitivity or too dark for the lowest sensitivity, report
103        # an error.
104        min_exposure_ns, max_exposure_ns = \
105            props['android.sensor.info.exposureTimeRange']
106        if auto_e < min_exposure_ns*sens_max:
107            raise its.error.Error("Scene is too bright to properly expose \
108                                  at the highest sensitivity")
109        if auto_e*bracket_factor > max_exposure_ns*sens_min:
110            raise its.error.Error("Scene is too dark to properly expose \
111                                  at the lowest sensitivity")
112
113        # Start the sensitivities at the minimum.
114        s = sens_min
115
116        samples = []
117        plots = []
118        measured_models = []
119        while s <= sens_max + 1:
120            print "ISO %d" % round(s)
121            fig = plt.figure()
122            plt_s = fig.gca()
123            plt_s.set_title("ISO %d" % round(s))
124            plt_s.set_xlabel("Mean signal level")
125            plt_s.set_ylabel("Variance")
126
127            samples_s = []
128            for b in range(0, bracket_stops + 1):
129                # Get the exposure for this sensitivity and exposure time.
130                e = int(math.pow(2, b)*auto_e/float(s))
131                req = its.objects.manual_capture_request(round(s), e)
132                cap = cam.do_capture(req, cam.CAP_RAW)
133                planes = its.image.convert_capture_to_planes(cap, props)
134
135                samples_e = []
136                for (pidx, p) in enumerate(planes):
137                    p = p.squeeze()
138
139                    # Crop the plane to be a multiple of the tile size.
140                    p = p[0:p.shape[0] - p.shape[0]%tile_size,
141                          0:p.shape[1] - p.shape[1]%tile_size]
142
143                    # convert_capture_to_planes normalizes the range
144                    # to [0, 1], but without subtracting the black
145                    # level.
146                    black_level = its.image.get_black_level(
147                        pidx, props, cap["metadata"])
148                    p *= white_level
149                    p = (p - black_level)/(white_level - black_level)
150
151                    # Use our high pass filter to filter this plane.
152                    hp = scipy.signal.sepfir2d(p, f, f).astype('float32')
153
154                    means_tiled = \
155                        np.mean(tile(p, tile_size), axis=(0, 1)).flatten()
156                    vars_tiled = \
157                        np.var(tile(hp, tile_size), axis=(0, 1)).flatten()
158
159                    for (mean, var) in zip(means_tiled, vars_tiled):
160                        # Don't include the tile if it has samples that might
161                        # be clipped.
162                        if mean + 2*math.sqrt(var) < max_signal_level:
163                            samples_e.append([mean, var])
164
165                    means_e, vars_e = zip(*samples_e)
166                    plt_s.plot(means_e, vars_e, colors[b%len(colors)] + ',')
167
168                    samples_s.extend(samples_e)
169
170            [S, O, R, p, stderr] = scipy.stats.linregress(samples_s)
171            measured_models.append([round(s), S, O])
172            print "Sensitivity %d: %e*y + %e (R=%f)" % (round(s), S, O, R)
173
174            # Add the samples for this sensitivity to the global samples list.
175            samples.extend([(round(s), mean, var) for (mean, var) in samples_s])
176
177            # Add the linear fit to the plot for this sensitivity.
178            plt_s.plot([0, max_signal_level], [O, O + S*max_signal_level], 'r-',
179                       label="Linear fit")
180            xmax = max([x for (x, _) in samples_s])*1.25
181            plt_s.set_xlim(xmin=0, xmax=xmax)
182            plt_s.set_ylim(ymin=0, ymax=(O + S*xmax)*1.25)
183            fig.savefig("%s_samples_iso%04d.png" % (NAME, round(s)))
184            plots.append([round(s), fig])
185
186            # Move to the next sensitivity.
187            s *= math.pow(2, 1.0/steps_per_stop)
188
189        # Grab the sensitivities and line parameters from each sensitivity.
190        S_measured = [e[1] for e in measured_models]
191        O_measured = [e[2] for e in measured_models]
192        sens = np.asarray([e[0] for e in measured_models])
193        sens_sq = np.square(sens)
194
195        # Use a global linear optimization to fit the noise model.
196        gains = np.asarray([s[0] for s in samples])
197        means = np.asarray([s[1] for s in samples])
198        vars_ = np.asarray([s[2] for s in samples])
199
200        # Define digital gain as the gain above the max analog gain
201        # per the Camera2 spec. Also, define a corresponding C
202        # expression snippet to use in the generated model code.
203        digital_gains = np.maximum(gains/sens_max_analog, 1)
204        digital_gain_cdef = "(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)" % \
205            (sens_max_analog, sens_max_analog)
206
207        # Find the noise model parameters via least squares fit.
208        ad = gains*means
209        bd = means
210        cd = gains*gains
211        dd = digital_gains*digital_gains
212        a = np.asarray([ad, bd, cd, dd]).T
213        b = vars_
214
215        # To avoid overfitting to high ISOs (high variances), divide the system
216        # by the gains.
217        a /= (np.tile(gains, (a.shape[1], 1)).T)
218        b /= gains
219
220        [A, B, C, D], _, _, _ = np.linalg.lstsq(a, b)
221
222        # Plot the noise model components with the values predicted by the
223        # noise model.
224        S_model = A*sens + B
225        O_model = \
226            C*sens_sq + D*np.square(np.maximum(sens/sens_max_analog, 1))
227
228        (fig, (plt_S, plt_O)) = plt.subplots(2, 1)
229        plt_S.set_title("Noise model")
230        plt_S.set_ylabel("S")
231        plt_S.loglog(sens, S_measured, 'r+', basex=10, basey=10,
232                     label="Measured")
233        plt_S.loglog(sens, S_model, 'bx', basex=10, basey=10, label="Model")
234        plt_S.legend(loc=2)
235
236        plt_O.set_xlabel("ISO")
237        plt_O.set_ylabel("O")
238        plt_O.loglog(sens, O_measured, 'r+', basex=10, basey=10,
239                     label="Measured")
240        plt_O.loglog(sens, O_model, 'bx', basex=10, basey=10, label="Model")
241        fig.savefig("%s.png" % (NAME))
242
243        for [s, fig] in plots:
244            plt_s = fig.gca()
245
246            dg = max(s/sens_max_analog, 1)
247            S = A*s + B
248            O = C*s*s + D*dg*dg
249            plt_s.plot([0, max_signal_level], [O, O + S*max_signal_level], 'b-',
250                       label="Model")
251            plt_s.legend(loc=2)
252
253            plt.figure(fig.number)
254
255            # Re-save the plot with the global model.
256            fig.savefig("%s_samples_iso%04d.png" % (NAME, round(s)))
257
258        # Generate the noise model implementation.
259        noise_model_code = textwrap.dedent("""\
260            /* Generated test code to dump a table of data for external validation
261             * of the noise model parameters.
262             */
263            #include <stdio.h>
264            #include <assert.h>
265            double compute_noise_model_entry_S(int sens);
266            double compute_noise_model_entry_O(int sens);
267            int main(void) {
268                int sens;
269                for (sens = %d; sens <= %d; sens += 100) {
270                    double o = compute_noise_model_entry_O(sens);
271                    double s = compute_noise_model_entry_S(sens);
272                    printf("%%d,%%lf,%%lf\\n", sens, o, s);
273                }
274                return 0;
275            }
276
277            /* Generated functions to map a given sensitivity to the O and S noise
278             * model parameters in the DNG noise model.
279             */
280            double compute_noise_model_entry_S(int sens) {
281                double s = %e * sens + %e;
282                return s < 0.0 ? 0.0 : s;
283            }
284
285            double compute_noise_model_entry_O(int sens) {
286                double digital_gain = %s;
287                double o = %e * sens * sens + %e * digital_gain * digital_gain;
288                return o < 0.0 ? 0.0 : o;
289            }
290            """ % (sens_min, sens_max, A, B, digital_gain_cdef, C, D))
291        print noise_model_code
292        text_file = open("noise_model.c", "w")
293        text_file.write("%s" % noise_model_code)
294        text_file.close()
295
296if __name__ == '__main__':
297    main()
298
299