1# Copyright 2014 The Android Open Source Project. Lint as python2 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 math 16import os.path 17import textwrap 18 19import its.caps 20import its.device 21import its.image 22import its.objects 23 24from matplotlib import pylab 25import matplotlib.pyplot as plt 26import numpy as np 27import scipy.signal 28import scipy.stats 29 30BAYER_LIST = ['R', 'GR', 'GB', 'B'] 31BRACKET_MAX = 8 # Exposure bracketing range in stops 32COLORS = 'rygcbm' # Colors used for plotting the data for each exposure. 33MAX_SCALE_FUDGE = 1.1 34MAX_SIGNAL_VALUE = 0.25 # Maximum value to allow mean of the tiles to go. 35NAME = os.path.basename(__file__).split('.')[0] 36RTOL_EXP_GAIN = 0.97 37STEPS_PER_STOP = 2 # How many sensitivities per stop to sample. 38# How large of tiles to use to compute mean/variance. 39# Large tiles may have their variance corrupted by low frequency 40# image changes (lens shading, scene illumination). 41TILE_SIZE = 32 42 43 44def main(): 45 """Capture a set of raw images with increasing analog gains and measure the noise. 46 """ 47 48 bracket_factor = math.pow(2, BRACKET_MAX) 49 50 with its.device.ItsSession() as cam: 51 props = cam.get_camera_properties() 52 props = cam.override_with_hidden_physical_camera_props(props) 53 54 # Get basic properties we need. 55 sens_min, sens_max = props['android.sensor.info.sensitivityRange'] 56 sens_max_analog = props['android.sensor.maxAnalogSensitivity'] 57 sens_max_meas = sens_max_analog 58 white_level = props['android.sensor.info.whiteLevel'] 59 60 print 'Sensitivity range: [%f, %f]' % (sens_min, sens_max) 61 print 'Max analog sensitivity: %f' % (sens_max_analog) 62 63 # Do AE to get a rough idea of where we are. 64 s_ae, e_ae, _, _, _ = cam.do_3a( 65 get_results=True, do_awb=False, do_af=False) 66 # Underexpose to get more data for low signal levels. 67 auto_e = s_ae * e_ae / bracket_factor 68 # Focus at zero to intentionally blur the scene as much as possible. 69 f_dist = 0.0 70 71 # If the auto-exposure result is too bright for the highest 72 # sensitivity or too dark for the lowest sensitivity, report 73 # an error. 74 min_exposure_ns, max_exposure_ns = props[ 75 'android.sensor.info.exposureTimeRange'] 76 if auto_e < min_exposure_ns*sens_max_meas: 77 raise its.error.Error('Scene is too bright to properly expose ' 78 'at the highest sensitivity') 79 if auto_e*bracket_factor > max_exposure_ns*sens_min: 80 raise its.error.Error('Scene is too dark to properly expose ' 81 'at the lowest sensitivity') 82 83 # Start the sensitivities at the minimum. 84 s = sens_min 85 86 samples = [[], [], [], []] 87 plots = [] 88 measured_models = [[], [], [], []] 89 color_plane_plots = {} 90 while int(round(s)) <= sens_max_meas: 91 s_int = int(round(s)) 92 print 'ISO %d' % s_int 93 fig, [[plt_r, plt_gr], [plt_gb, plt_b]] = plt.subplots(2, 2, figsize=(11, 11)) 94 fig.gca() 95 color_plane_plots[s_int] = [plt_r, plt_gr, plt_gb, plt_b] 96 fig.suptitle('ISO %d' % s_int, x=0.54, y=0.99) 97 for i, plot in enumerate(color_plane_plots[s_int]): 98 plot.set_title('%s' % BAYER_LIST[i]) 99 plot.set_xlabel('Mean signal level') 100 plot.set_ylabel('Variance') 101 102 samples_s = [[], [], [], []] 103 for b in range(BRACKET_MAX): 104 # Get the exposure for this sensitivity and exposure time. 105 e = int(math.pow(2, b)*auto_e/float(s)) 106 print 'exp %.3fms' % round(e*1.0E-6, 3) 107 req = its.objects.manual_capture_request(s_int, e, f_dist) 108 fmt_raw = {'format': 'rawStats', 109 'gridWidth': TILE_SIZE, 110 'gridHeight': TILE_SIZE} 111 cap = cam.do_capture(req, fmt_raw) 112 mean_image, var_image = its.image.unpack_rawstats_capture(cap) 113 idxs = its.image.get_canonical_cfa_order(props) 114 means = [mean_image[:, :, i] for i in idxs] 115 vars_ = [var_image[:, :, i] for i in idxs] 116 117 s_read = cap['metadata']['android.sensor.sensitivity'] 118 s_err = 's_write: %d, s_read: %d, RTOL: %.2f' % ( 119 s, s_read, RTOL_EXP_GAIN) 120 assert (1.0 >= s_read/float(s_int) >= RTOL_EXP_GAIN), s_err 121 print 'ISO_write: %d, ISO_read: %d' % (s_int, s_read) 122 123 for pidx in range(len(means)): 124 plot = color_plane_plots[s_int][pidx] 125 126 # convert_capture_to_planes normalizes the range 127 # to [0, 1], but without subtracting the black 128 # level. 129 black_level = its.image.get_black_level( 130 pidx, props, cap['metadata']) 131 means_p = (means[pidx] - black_level)/(white_level - black_level) 132 vars_p = vars_[pidx]/((white_level - black_level)**2) 133 134 # TODO(dsharlet): It should be possible to account for low 135 # frequency variation by looking at neighboring means, but I 136 # have not been able to make this work. 137 138 means_p = np.asarray(means_p).flatten() 139 vars_p = np.asarray(vars_p).flatten() 140 141 samples_e = [] 142 for (mean, var) in zip(means_p, vars_p): 143 # Don't include the tile if it has samples that might 144 # be clipped. 145 if mean + 2*math.sqrt(max(var, 0)) < MAX_SIGNAL_VALUE: 146 samples_e.append([mean, var]) 147 148 if samples_e: 149 means_e, vars_e = zip(*samples_e) 150 color_plane_plots[s_int][pidx].plot( 151 means_e, vars_e, COLORS[b%len(COLORS)] + '.', 152 alpha=0.5, markersize=1) 153 samples_s[pidx].extend(samples_e) 154 155 for (pidx, p) in enumerate(samples_s): 156 [S, O, R, _, _] = scipy.stats.linregress(samples_s[pidx]) 157 measured_models[pidx].append([s_int, S, O]) 158 print "Sensitivity %d: %e*y + %e (R=%f)" % (s_int, S, O, R) 159 160 # Add the samples for this sensitivity to the global samples list. 161 samples[pidx].extend([(s_int, mean, var) for (mean, var) in samples_s[pidx]]) 162 163 # Add the linear fit to subplot for this sensitivity. 164 color_plane_plots[s_int][pidx].plot( 165 [0, MAX_SIGNAL_VALUE], [O, O + S*MAX_SIGNAL_VALUE], 166 'rgkb'[pidx]+'--', label='Linear fit') 167 168 xmax = max([max([x for (x, _) in p]) for p in samples_s])*MAX_SCALE_FUDGE 169 ymax = (O + S*xmax)*MAX_SCALE_FUDGE 170 color_plane_plots[s_int][pidx].set_xlim(xmin=0, xmax=xmax) 171 color_plane_plots[s_int][pidx].set_ylim(ymin=0, ymax=ymax) 172 color_plane_plots[s_int][pidx].legend() 173 pylab.tight_layout() 174 175 fig.savefig('%s_samples_iso%04d.png' % (NAME, s_int)) 176 plots.append([s_int, fig]) 177 178 # Move to the next sensitivity. 179 s *= math.pow(2, 1.0/STEPS_PER_STOP) 180 181 # do model plots 182 (fig, (plt_S, plt_O)) = plt.subplots(2, 1, figsize=(11, 8.5)) 183 plt_S.set_title("Noise model") 184 plt_S.set_ylabel("S") 185 plt_O.set_xlabel("ISO") 186 plt_O.set_ylabel("O") 187 188 A = [] 189 B = [] 190 C = [] 191 D = [] 192 for (pidx, p) in enumerate(measured_models): 193 # Grab the sensitivities and line parameters from each sensitivity. 194 S_measured = [e[1] for e in measured_models[pidx]] 195 O_measured = [e[2] for e in measured_models[pidx]] 196 sens = np.asarray([e[0] for e in measured_models[pidx]]) 197 sens_sq = np.square(sens) 198 199 # Use a global linear optimization to fit the noise model. 200 gains = np.asarray([s[0] for s in samples[pidx]]) 201 means = np.asarray([s[1] for s in samples[pidx]]) 202 vars_ = np.asarray([s[2] for s in samples[pidx]]) 203 gains = gains.flatten() 204 means = means.flatten() 205 vars_ = vars_.flatten() 206 207 # Define digital gain as the gain above the max analog gain 208 # per the Camera2 spec. Also, define a corresponding C 209 # expression snippet to use in the generated model code. 210 digital_gains = np.maximum(gains/sens_max_analog, 1) 211 assert np.all(digital_gains == 1) 212 digital_gain_cdef = '(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)' % \ 213 (sens_max_analog, sens_max_analog) 214 215 # Divide the whole system by gains*means. 216 f = lambda x, a, b, c, d: (c*(x[0]**2) + d + (x[1])*a*x[0] + (x[1])*b)/(x[0]) 217 [result, _] = scipy.optimize.curve_fit(f, (gains, means), vars_/(gains)) 218 219 [A_p, B_p, C_p, D_p] = result[0:4] 220 A.append(A_p) 221 B.append(B_p) 222 C.append(C_p) 223 D.append(D_p) 224 225 # Plot the noise model components with the values predicted by the 226 # noise model. 227 S_model = A_p*sens + B_p 228 O_model = \ 229 C_p*sens_sq + D_p*np.square(np.maximum(sens/sens_max_analog, 1)) 230 231 plt_S.loglog(sens, S_measured, 'rgkb'[pidx]+'+', basex=10, basey=10, 232 label='Measured') 233 plt_S.loglog(sens, S_model, 'rgkb'[pidx]+'x', basex=10, basey=10, 234 label='Model') 235 plt_O.loglog(sens, O_measured, 'rgkb'[pidx]+'+', basex=10, basey=10, 236 label='Measured') 237 plt_O.loglog(sens, O_model, 'rgkb'[pidx]+'x', basex=10, basey=10, 238 label='Model') 239 plt_S.legend() 240 plt_O.legend() 241 242 fig.savefig('%s.png' % (NAME)) 243 244 # add models to subplots and re-save 245 for [s, fig] in plots: # re-step through figs... 246 dg = max(s/sens_max_analog, 1) 247 fig.gca() 248 for (pidx, p) in enumerate(measured_models): 249 S = A[pidx]*s + B[pidx] 250 O = C[pidx]*s*s + D[pidx]*dg*dg 251 color_plane_plots[s][pidx].plot( 252 [0, MAX_SIGNAL_VALUE], [O, O + S*MAX_SIGNAL_VALUE], 253 'rgkb'[pidx]+'-', label='Model', alpha=0.5) 254 color_plane_plots[s][pidx].legend(loc='upper left') 255 fig.savefig('%s_samples_iso%04d.png' % (NAME, s)) 256 257 # Generate the noise model implementation. 258 A_array = ",".join([str(i) for i in A]) 259 B_array = ",".join([str(i) for i in B]) 260 C_array = ",".join([str(i) for i in C]) 261 D_array = ",".join([str(i) for i in D]) 262 noise_model_code = textwrap.dedent("""\ 263 /* Generated test code to dump a table of data for external validation 264 * of the noise model parameters. 265 */ 266 #include <stdio.h> 267 #include <assert.h> 268 double compute_noise_model_entry_S(int plane, int sens); 269 double compute_noise_model_entry_O(int plane, int sens); 270 int main(void) { 271 for (int plane = 0; plane < %d; plane++) { 272 for (int sens = %d; sens <= %d; sens += 100) { 273 double o = compute_noise_model_entry_O(plane, sens); 274 double s = compute_noise_model_entry_S(plane, sens); 275 printf("%%d,%%d,%%lf,%%lf\\n", plane, sens, o, s); 276 } 277 } 278 return 0; 279 } 280 281 /* Generated functions to map a given sensitivity to the O and S noise 282 * model parameters in the DNG noise model. The planes are in 283 * R, Gr, Gb, B order. 284 */ 285 double compute_noise_model_entry_S(int plane, int sens) { 286 static double noise_model_A[] = { %s }; 287 static double noise_model_B[] = { %s }; 288 double A = noise_model_A[plane]; 289 double B = noise_model_B[plane]; 290 double s = A * sens + B; 291 return s < 0.0 ? 0.0 : s; 292 } 293 294 double compute_noise_model_entry_O(int plane, int sens) { 295 static double noise_model_C[] = { %s }; 296 static double noise_model_D[] = { %s }; 297 double digital_gain = %s; 298 double C = noise_model_C[plane]; 299 double D = noise_model_D[plane]; 300 double o = C * sens * sens + D * digital_gain * digital_gain; 301 return o < 0.0 ? 0.0 : o; 302 } 303 """ % (len(A), sens_min, sens_max, A_array, B_array, C_array, D_array, digital_gain_cdef)) 304 print noise_model_code 305 for i, _ in enumerate(BAYER_LIST): 306 read_noise = C[i] * sens_min * sens_min + D[i] 307 e_msg = '%s model min ISO noise < 0! C: %.4e, D: %.4e, rn: %.4e' % ( 308 BAYER_LIST[i], C[i], D[i], read_noise) 309 assert read_noise > 0, e_msg 310 assert C[i] > 0, '%s model slope is negative. slope=%.4e' % ( 311 BAYER_LIST[i], C[i]) 312 text_file = open("noise_model.c", "w") 313 text_file.write("%s" % noise_model_code) 314 text_file.close() 315 316if __name__ == '__main__': 317 main() 318