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
20import pylab
21import matplotlib
22import matplotlib.pyplot as plt
23import math
24import Image
25import time
26import numpy as np
27import scipy.stats
28import scipy.signal
29
30# Convert a 2D array a to a 4D array with dimensions [tile_size,
31# tile_size, row, col] where row, col are tile indices.
32def tile(a, tile_size):
33    tile_rows, tile_cols = a.shape[0]/tile_size, a.shape[1]/tile_size
34    a = a.reshape([tile_rows, tile_size, tile_cols, tile_size])
35    a = a.transpose([1, 3, 0, 2])
36    return a
37
38def main():
39    """Capture a set of raw images with increasing gains and measure the noise.
40    """
41    NAME = os.path.basename(__file__).split(".")[0]
42
43    # How many sensitivities per stop to sample.
44    steps_per_stop = 2
45    # How large of tiles to use to compute mean/variance.
46    tile_size = 64
47    # Exposure bracketing range in stops
48    bracket_stops = 4
49    # How high to allow the mean of the tiles to go.
50    max_signal_level = 0.5
51    # Colors used for plotting the data for each exposure.
52    colors = 'rygcbm'
53
54    # Define a first order high pass filter to eliminate low frequency
55    # signal content when computing variance.
56    f = np.array([-1, 1]).astype('float32')
57    # Make it a higher order filter by convolving the first order
58    # filter with itself a few times.
59    f = np.convolve(f, f)
60    f = np.convolve(f, f)
61
62    # Compute the normalization of the filter to preserve noise
63    # power. Let a be the normalization factor we're looking for, and
64    # Let X and X' be the random variables representing the noise
65    # before and after filtering, respectively. First, compute
66    # Var[a*X']:
67    #
68    #   Var[a*X'] = a^2*Var[X*f_0 + X*f_1 + ... + X*f_N-1]
69    #             = a^2*(f_0^2*Var[X] + f_1^2*Var[X] + ... + (f_N-1)^2*Var[X])
70    #             = sum(f_i^2)*a^2*Var[X]
71    #
72    # We want Var[a*X'] to be equal to Var[X]:
73    #
74    #    sum(f_i^2)*a^2*Var[X] = Var[X] -> a = sqrt(1/sum(f_i^2))
75    #
76    # We can just bake this normalization factor into the high pass
77    # filter kernel.
78    f = f/math.sqrt(np.dot(f, f))
79
80    bracket_factor = math.pow(2, bracket_stops)
81
82    with its.device.ItsSession() as cam:
83        props = cam.get_camera_properties()
84
85        # Get basic properties we need.
86        sens_min, sens_max = props['android.sensor.info.sensitivityRange']
87        sens_max_analog = props['android.sensor.maxAnalogSensitivity']
88        white_level = props['android.sensor.info.whiteLevel']
89        black_levels = props['android.sensor.blackLevelPattern']
90        idxs = its.image.get_canonical_cfa_order(props)
91        black_levels = [black_levels[i] for i in idxs]
92
93        print "Sensitivity range: [%f, %f]" % (sens_min, sens_max)
94        print "Max analog sensitivity: %f" % (sens_max_analog)
95
96        # Do AE to get a rough idea of where we are.
97        s_ae,e_ae,_,_,_  = \
98            cam.do_3a(get_results=True, do_awb=False, do_af=False)
99        # Underexpose to get more data for low signal levels.
100        auto_e = s_ae*e_ae/bracket_factor
101
102        # If the auto-exposure result is too bright for the highest
103        # sensitivity or too dark for the lowest sensitivity, report
104        # an error.
105        min_exposure_ns, max_exposure_ns = \
106            props['android.sensor.info.exposureTimeRange']
107        if auto_e < min_exposure_ns*sens_max:
108            raise its.error.Error("Scene is too bright to properly expose \
109                                  at the highest sensitivity")
110        if auto_e*bracket_factor > max_exposure_ns*sens_min:
111            raise its.error.Error("Scene is too dark to properly expose \
112                                  at the lowest sensitivity")
113
114        # Start the sensitivities at the minimum.
115        s = sens_min
116
117        samples = []
118        plots = []
119        measured_models = []
120        while s <= sens_max + 1:
121            print "ISO %d" % round(s)
122            fig = plt.figure()
123            plt_s = fig.gca()
124            plt_s.set_title("ISO %d" % round(s))
125            plt_s.set_xlabel("Mean signal level")
126            plt_s.set_ylabel("Variance")
127
128            samples_s = []
129            for b in range(0, bracket_stops + 1):
130                # Get the exposure for this sensitivity and exposure time.
131                e = int(math.pow(2, b)*auto_e/float(s))
132                req = its.objects.manual_capture_request(round(s), e)
133                cap = cam.do_capture(req, cam.CAP_RAW)
134                planes = its.image.convert_capture_to_planes(cap, props)
135
136                samples_e = []
137                for (pidx, p) in enumerate(planes):
138                    p = p.squeeze()
139
140                    # Crop the plane to be a multiple of the tile size.
141                    p = p[0:p.shape[0] - p.shape[0]%tile_size,
142                          0:p.shape[1] - p.shape[1]%tile_size]
143
144                    # convert_capture_to_planes normalizes the range
145                    # to [0, 1], but without subtracting the black
146                    # level.
147                    black_level = black_levels[pidx]
148                    p = 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 = 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 = a/(np.tile(gains, (a.shape[1], 1)).T)
218        b = 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        print """
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
292if __name__ == '__main__':
293    main()
294
295