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 math 16import os.path 17import textwrap 18 19import its.caps 20import its.device 21import its.image 22import its.objects 23 24import matplotlib 25from matplotlib import pylab 26import matplotlib.pyplot as plt 27import numpy as np 28import scipy.signal 29import scipy.stats 30 31BAYER_LIST = ['R', 'GR', 'GB', 'B'] 32NAME = os.path.basename(__file__).split('.')[0] 33RTOL_EXP_GAIN = 0.97 34 35 36def tile(a, tile_size): 37 """Convert a 2D array to 4D w/ dims [tile_size, tile_size, row, col] where row, col are tile indices. 38 """ 39 tile_rows, tile_cols = a.shape[0]/tile_size, a.shape[1]/tile_size 40 a = a.reshape([tile_rows, tile_size, tile_cols, tile_size]) 41 a = a.transpose([1, 3, 0, 2]) 42 return a 43 44 45def main(): 46 """Capture a set of raw images with increasing analog gains and measure the noise. 47 """ 48 49 # How many sensitivities per stop to sample. 50 steps_per_stop = 2 51 # How large of tiles to use to compute mean/variance. 52 tile_size = 64 53 # Exposure bracketing range in stops 54 bracket_stops = 4 55 # How high to allow the mean of the tiles to go. 56 max_signal_level = 0.25 57 # Colors used for plotting the data for each exposure. 58 colors = 'rygcbm' 59 60 # Define a first order high pass filter to eliminate low frequency 61 # signal content when computing variance. 62 f = np.array([-1, 1]).astype('float32') 63 # Make it a higher order filter by convolving the first order 64 # filter with itself a few times. 65 f = np.convolve(f, f) 66 f = np.convolve(f, f) 67 68 # Compute the normalization of the filter to preserve noise 69 # power. Let a be the normalization factor we're looking for, and 70 # Let X and X' be the random variables representing the noise 71 # before and after filtering, respectively. First, compute 72 # Var[a*X']: 73 # 74 # Var[a*X'] = a^2*Var[X*f_0 + X*f_1 + ... + X*f_N-1] 75 # = a^2*(f_0^2*Var[X] + f_1^2*Var[X] + ... + (f_N-1)^2*Var[X]) 76 # = sum(f_i^2)*a^2*Var[X] 77 # 78 # We want Var[a*X'] to be equal to Var[X]: 79 # 80 # sum(f_i^2)*a^2*Var[X] = Var[X] -> a = sqrt(1/sum(f_i^2)) 81 # 82 # We can just bake this normalization factor into the high pass 83 # filter kernel. 84 f /= math.sqrt(np.dot(f, f)) 85 86 bracket_factor = math.pow(2, bracket_stops) 87 88 with its.device.ItsSession() as cam: 89 props = cam.get_camera_properties() 90 props = cam.override_with_hidden_physical_camera_props(props) 91 92 # Get basic properties we need. 93 sens_min, sens_max = props['android.sensor.info.sensitivityRange'] 94 sens_max_analog = props['android.sensor.maxAnalogSensitivity'] 95 sens_max_meas = sens_max_analog 96 white_level = props['android.sensor.info.whiteLevel'] 97 98 print "Sensitivity range: [%f, %f]" % (sens_min, sens_max) 99 print "Max analog sensitivity: %f" % (sens_max_analog) 100 101 # Do AE to get a rough idea of where we are. 102 s_ae, e_ae, _, _, _ = \ 103 cam.do_3a(get_results=True, do_awb=False, do_af=False) 104 # Underexpose to get more data for low signal levels. 105 auto_e = s_ae*e_ae/bracket_factor 106 # Focus at zero to intentionally blur the scene as much as possible. 107 f_dist = 0.0 108 109 # If the auto-exposure result is too bright for the highest 110 # sensitivity or too dark for the lowest sensitivity, report 111 # an error. 112 min_exposure_ns, max_exposure_ns = \ 113 props['android.sensor.info.exposureTimeRange'] 114 if auto_e < min_exposure_ns*sens_max_meas: 115 raise its.error.Error("Scene is too bright to properly expose \ 116 at the highest sensitivity") 117 if auto_e*bracket_factor > max_exposure_ns*sens_min: 118 raise its.error.Error("Scene is too dark to properly expose \ 119 at the lowest sensitivity") 120 121 # Start the sensitivities at the minimum. 122 s = sens_min 123 124 samples = [[], [], [], []] 125 plots = [] 126 measured_models = [[], [], [], []] 127 color_plane_plots = {} 128 while int(round(s)) <= sens_max_meas: 129 s_int = int(round(s)) 130 print 'ISO %d' % s_int 131 fig, [[plt_r, plt_gr], [plt_gb, plt_b]] = plt.subplots(2, 2, figsize=(11, 11)) 132 fig.gca() 133 color_plane_plots[s_int] = [plt_r, plt_gr, plt_gb, plt_b] 134 fig.suptitle('ISO %d' % s_int, x=0.54, y=0.99) 135 for i, plot in enumerate(color_plane_plots[s_int]): 136 plot.set_title('%s' % BAYER_LIST[i]) 137 plot.set_xlabel('Mean signal level') 138 plot.set_ylabel('Variance') 139 140 samples_s = [[], [], [], []] 141 for b in range(0, bracket_stops): 142 # Get the exposure for this sensitivity and exposure time. 143 e = int(math.pow(2, b)*auto_e/float(s)) 144 print 'exp %.3fms' % round(e*1.0E-6, 3) 145 req = its.objects.manual_capture_request(s_int, e, f_dist) 146 cap = cam.do_capture(req, cam.CAP_RAW) 147 planes = its.image.convert_capture_to_planes(cap, props) 148 e_read = cap['metadata']['android.sensor.exposureTime'] 149 s_read = cap['metadata']['android.sensor.sensitivity'] 150 s_err = 's_write: %d, s_read: %d, RTOL: %.2f' % ( 151 s, s_read, RTOL_EXP_GAIN) 152 assert (1.0 >= s_read/float(s_int) >= RTOL_EXP_GAIN), s_err 153 print 'ISO_write: %d, ISO_read: %d' % (s_int, s_read) 154 155 for (pidx, p) in enumerate(planes): 156 plot = color_plane_plots[s_int][pidx] 157 158 p = p.squeeze() 159 160 # Crop the plane to be a multiple of the tile size. 161 p = p[0:p.shape[0] - p.shape[0]%tile_size, 162 0:p.shape[1] - p.shape[1]%tile_size] 163 164 # convert_capture_to_planes normalizes the range 165 # to [0, 1], but without subtracting the black 166 # level. 167 black_level = its.image.get_black_level( 168 pidx, props, cap["metadata"]) 169 p *= white_level 170 p = (p - black_level)/(white_level - black_level) 171 172 # Use our high pass filter to filter this plane. 173 hp = scipy.signal.sepfir2d(p, f, f).astype('float32') 174 175 means_tiled = \ 176 np.mean(tile(p, tile_size), axis=(0, 1)).flatten() 177 vars_tiled = \ 178 np.var(tile(hp, tile_size), axis=(0, 1)).flatten() 179 180 samples_e = [] 181 for (mean, var) in zip(means_tiled, vars_tiled): 182 # Don't include the tile if it has samples that might 183 # be clipped. 184 if mean + 2*math.sqrt(var) < max_signal_level: 185 samples_e.append([mean, var]) 186 187 if samples_e: 188 means_e, vars_e = zip(*samples_e) 189 color_plane_plots[s_int][pidx].plot( 190 means_e, vars_e, colors[b%len(colors)] + '.', 191 alpha=0.5) 192 samples_s[pidx].extend(samples_e) 193 194 for (pidx, p) in enumerate(samples_s): 195 [S, O, R, _, _] = scipy.stats.linregress(samples_s[pidx]) 196 measured_models[pidx].append([s_int, S, O]) 197 print "Sensitivity %d: %e*y + %e (R=%f)" % (s_int, S, O, R) 198 199 # Add the samples for this sensitivity to the global samples list. 200 samples[pidx].extend([(s_int, mean, var) for (mean, var) in samples_s[pidx]]) 201 202 # Add the linear fit to subplot for this sensitivity. 203 # pylab.subplot(2, 2, pidx+1) 204 #pylab.plot([0, max_signal_level], [O, O + S*max_signal_level], 'rgkb'[pidx]+'--', 205 #label="Linear fit") 206 color_plane_plots[s_int][pidx].plot([0, max_signal_level], 207 [O, O + S*max_signal_level], 'rgkb'[pidx]+'--', 208 label="Linear fit") 209 210 xmax = max([max([x for (x, _) in p]) for p in samples_s])*1.25 211 ymax = max([max([y for (_, y) in p]) for p in samples_s])*1.25 212 color_plane_plots[s_int][pidx].set_xlim(xmin=0, xmax=xmax) 213 color_plane_plots[s_int][pidx].set_ylim(ymin=0, ymax=ymax) 214 color_plane_plots[s_int][pidx].legend() 215 pylab.tight_layout() 216 217 fig.savefig('%s_samples_iso%04d.png' % (NAME, s_int)) 218 plots.append([s_int, fig]) 219 220 # Move to the next sensitivity. 221 s *= math.pow(2, 1.0/steps_per_stop) 222 223 # do model plots 224 (fig, (plt_S, plt_O)) = plt.subplots(2, 1, figsize=(11, 8.5)) 225 plt_S.set_title("Noise model") 226 plt_S.set_ylabel("S") 227 plt_O.set_xlabel("ISO") 228 plt_O.set_ylabel("O") 229 230 A = [] 231 B = [] 232 C = [] 233 D = [] 234 for (pidx, p) in enumerate(measured_models): 235 # Grab the sensitivities and line parameters from each sensitivity. 236 S_measured = [e[1] for e in measured_models[pidx]] 237 O_measured = [e[2] for e in measured_models[pidx]] 238 sens = np.asarray([e[0] for e in measured_models[pidx]]) 239 sens_sq = np.square(sens) 240 241 # Use a global linear optimization to fit the noise model. 242 gains = np.asarray([s[0] for s in samples[pidx]]) 243 means = np.asarray([s[1] for s in samples[pidx]]) 244 vars_ = np.asarray([s[2] for s in samples[pidx]]) 245 246 # Define digital gain as the gain above the max analog gain 247 # per the Camera2 spec. Also, define a corresponding C 248 # expression snippet to use in the generated model code. 249 digital_gains = np.maximum(gains/sens_max_analog, 1) 250 digital_gain_cdef = "(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)" % \ 251 (sens_max_analog, sens_max_analog) 252 253 # Find the noise model parameters via least squares fit. 254 ad = gains*means 255 bd = means 256 cd = gains*gains 257 dd = digital_gains*digital_gains 258 a = np.asarray([ad, bd, cd, dd]).T 259 b = vars_ 260 261 # To avoid overfitting to high ISOs (high variances), divide the system 262 # by the gains. 263 a /= (np.tile(gains, (a.shape[1], 1)).T) 264 b /= gains 265 266 [A_p, B_p, C_p, D_p], _, _, _ = np.linalg.lstsq(a, b) 267 A.append(A_p) 268 B.append(B_p) 269 C.append(C_p) 270 D.append(D_p) 271 272 # Plot the noise model components with the values predicted by the 273 # noise model. 274 S_model = A_p*sens + B_p 275 O_model = \ 276 C_p*sens_sq + D_p*np.square(np.maximum(sens/sens_max_analog, 1)) 277 278 plt_S.loglog(sens, S_measured, 'rgkb'[pidx]+'+', basex=10, basey=10, 279 label="Measured") 280 plt_S.loglog(sens, S_model, 'rgkb'[pidx]+'x', basex=10, basey=10, label="Model") 281 282 plt_O.loglog(sens, O_measured, 'rgkb'[pidx]+'+', basex=10, basey=10, 283 label="Measured") 284 plt_O.loglog(sens, O_model, 'rgkb'[pidx]+'x', basex=10, basey=10, label="Model") 285 plt_S.legend() 286 plt_O.legend() 287 288 fig.savefig("%s.png" % (NAME)) 289 290 # add models to subplots and re-save 291 for [s, fig] in plots: # re-step through figs... 292 dg = max(s/sens_max_analog, 1) 293 fig.gca() 294 for (pidx, p) in enumerate(measured_models): 295 S = A[pidx]*s + B[pidx] 296 O = C[pidx]*s*s + D[pidx]*dg*dg 297 color_plane_plots[s][pidx].plot([0, max_signal_level], 298 [O, O + S*max_signal_level], 'rgkb'[pidx]+'-', 299 label="Model", alpha=0.5) 300 color_plane_plots[s][pidx].legend(loc='upper left') 301 fig.savefig('%s_samples_iso%04d.png' % (NAME, s)) 302 303 # Generate the noise model implementation. 304 A_array = ",".join([str(i) for i in A]) 305 B_array = ",".join([str(i) for i in B]) 306 C_array = ",".join([str(i) for i in C]) 307 D_array = ",".join([str(i) for i in D]) 308 noise_model_code = textwrap.dedent("""\ 309 /* Generated test code to dump a table of data for external validation 310 * of the noise model parameters. 311 */ 312 #include <stdio.h> 313 #include <assert.h> 314 double compute_noise_model_entry_S(int plane, int sens); 315 double compute_noise_model_entry_O(int plane, int sens); 316 int main(void) { 317 for (int plane = 0; plane < %d; plane++) { 318 for (int sens = %d; sens <= %d; sens += 100) { 319 double o = compute_noise_model_entry_O(plane, sens); 320 double s = compute_noise_model_entry_S(plane, sens); 321 printf("%%d,%%d,%%lf,%%lf\\n", plane, sens, o, s); 322 } 323 } 324 return 0; 325 } 326 327 /* Generated functions to map a given sensitivity to the O and S noise 328 * model parameters in the DNG noise model. The planes are in 329 * R, Gr, Gb, B order. 330 */ 331 double compute_noise_model_entry_S(int plane, int sens) { 332 static double noise_model_A[] = { %s }; 333 static double noise_model_B[] = { %s }; 334 double A = noise_model_A[plane]; 335 double B = noise_model_B[plane]; 336 double s = A * sens + B; 337 return s < 0.0 ? 0.0 : s; 338 } 339 340 double compute_noise_model_entry_O(int plane, int sens) { 341 static double noise_model_C[] = { %s }; 342 static double noise_model_D[] = { %s }; 343 double digital_gain = %s; 344 double C = noise_model_C[plane]; 345 double D = noise_model_D[plane]; 346 double o = C * sens * sens + D * digital_gain * digital_gain; 347 return o < 0.0 ? 0.0 : o; 348 } 349 """ % (len(A), sens_min, sens_max, A_array, B_array, C_array, D_array, digital_gain_cdef)) 350 print noise_model_code 351 for i, _ in enumerate(BAYER_LIST): 352 read_noise = C[i] * sens_min * sens_min + D[i] 353 e_msg = '%s model min ISO noise < 0! C: %.4e, D: %.4e, rn: %.4e' % ( 354 BAYER_LIST[i], C[i], D[i], read_noise) 355 assert read_noise > 0, e_msg 356 assert C[i] > 0, '%s model slope is negative. slope=%.4e' % ( 357 BAYER_LIST[i], C[i]) 358 text_file = open("noise_model.c", "w") 359 text_file.write("%s" % noise_model_code) 360 text_file.close() 361 362if __name__ == '__main__': 363 main() 364