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 logging 16import math 17import os.path 18import textwrap 19from matplotlib import pylab 20import matplotlib.pyplot as plt 21from mobly import test_runner 22import numpy as np 23import scipy.signal 24import scipy.stats 25 26import its_base_test 27import capture_request_utils 28import image_processing_utils 29import its_session_utils 30 31 32_BAYER_LIST = ('R', 'GR', 'GB', 'B') 33_BRACKET_MAX = 8 # Exposure bracketing range in stops 34_BRACKET_FACTOR = math.pow(2, _BRACKET_MAX) 35_PLOT_COLORS = 'rygcbm' # Colors used for plotting the data for each exposure. 36_MAX_SCALE_FUDGE = 1.1 37_MAX_SIGNAL_VALUE = 0.25 # Maximum value to allow mean of the tiles to go. 38_NAME = os.path.basename(__file__).split('.')[0] 39_RTOL_EXP_GAIN = 0.97 40_STEPS_PER_STOP = 3 # How many sensitivities per stop to sample. 41_TILE_SIZE = 32 # Tile size to compute mean/variance. Large tiles may have 42 # their variance corrupted by low freq image changes. 43_TILE_CROP_N = 0 # Number of tiles to crop from edge of image. Usually 0. 44 45 46def check_auto_exposure_targets(auto_e, sens_min, sens_max, props): 47 """Checks if AE too bright for highest gain & too dark for lowest gain.""" 48 49 min_exposure_ns, max_exposure_ns = props[ 50 'android.sensor.info.exposureTimeRange'] 51 if auto_e < min_exposure_ns*sens_max: 52 raise AssertionError('Scene is too bright to properly expose at highest ' 53 f'sensitivity: {sens_max}') 54 if auto_e*_BRACKET_FACTOR > max_exposure_ns*sens_min: 55 raise AssertionError('Scene is too dark to properly expose at lowest ' 56 f'sensitivity: {sens_min}') 57 58 59def create_noise_model_code(noise_model_a, noise_model_b, 60 noise_model_c, noise_model_d, 61 sens_min, sens_max, digital_gain_cdef, log_path): 62 """Creates the c file for the noise model.""" 63 64 noise_model_a_array = ','.join([str(i) for i in noise_model_a]) 65 noise_model_b_array = ','.join([str(i) for i in noise_model_b]) 66 noise_model_c_array = ','.join([str(i) for i in noise_model_c]) 67 noise_model_d_array = ','.join([str(i) for i in noise_model_d]) 68 code = textwrap.dedent(f"""\ 69 /* Generated test code to dump a table of data for external validation 70 * of the noise model parameters. 71 */ 72 #include <stdio.h> 73 #include <assert.h> 74 double compute_noise_model_entry_S(int plane, int sens); 75 double compute_noise_model_entry_O(int plane, int sens); 76 int main(void) {{ 77 for (int plane = 0; plane < {len(noise_model_a)}; plane++) {{ 78 for (int sens = {sens_min}; sens <= {sens_max}; sens += 100) {{ 79 double o = compute_noise_model_entry_O(plane, sens); 80 double s = compute_noise_model_entry_S(plane, sens); 81 printf("%d,%d,%lf,%lf\\n", plane, sens, o, s); 82 }} 83 }} 84 return 0; 85 }} 86 87 /* Generated functions to map a given sensitivity to the O and S noise 88 * model parameters in the DNG noise model. The planes are in 89 * R, Gr, Gb, B order. 90 */ 91 double compute_noise_model_entry_S(int plane, int sens) {{ 92 static double noise_model_A[] = {{ {noise_model_a_array:s} }}; 93 static double noise_model_B[] = {{ {noise_model_b_array:s} }}; 94 double A = noise_model_A[plane]; 95 double B = noise_model_B[plane]; 96 double s = A * sens + B; 97 return s < 0.0 ? 0.0 : s; 98 }} 99 100 double compute_noise_model_entry_O(int plane, int sens) {{ 101 static double noise_model_C[] = {{ {noise_model_c_array:s} }}; 102 static double noise_model_D[] = {{ {noise_model_d_array:s} }}; 103 double digital_gain = {digital_gain_cdef:s}; 104 double C = noise_model_C[plane]; 105 double D = noise_model_D[plane]; 106 double o = C * sens * sens + D * digital_gain * digital_gain; 107 return o < 0.0 ? 0.0 : o; 108 }} 109 """) 110 text_file = open(os.path.join(log_path, 'noise_model.c'), 'w') 111 text_file.write('%s' % code) 112 text_file.close() 113 114 # Creates the noise profile C++ file 115 code = textwrap.dedent(f"""\ 116 /* noise_profile.cc 117 Note: gradient_slope --> gradient of API slope parameter 118 offset_slope --> offset of API slope parameter 119 gradient_intercept--> gradient of API intercept parameter 120 offset_intercept --> offset of API intercept parameter 121 Note: SENSOR_NOISE_PROFILE in Android Developers doc uses 122 N(x) = sqrt(Sx + O), where 'S' is 'slope' & 'O' is 'intercept' 123 */ 124 .noise_profile = 125 {{.noise_coefficients_r = {{.gradient_slope = {noise_model_a[0]}, 126 .offset_slope = {noise_model_b[0]}, 127 .gradient_intercept = {noise_model_c[0]}, 128 .offset_intercept = {noise_model_d[0]}}}, 129 .noise_coefficients_gr = {{.gradient_slope = {noise_model_a[1]}, 130 .offset_slope = {noise_model_b[1]}, 131 .gradient_intercept = {noise_model_c[1]}, 132 .offset_intercept = {noise_model_d[1]}}}, 133 .noise_coefficients_gb = {{.gradient_slope = {noise_model_a[2]}, 134 .offset_slope = {noise_model_b[2]}, 135 .gradient_intercept = {noise_model_c[2]}, 136 .offset_intercept = {noise_model_d[2]}}}, 137 .noise_coefficients_b = {{.gradient_slope = {noise_model_a[3]}, 138 .offset_slope = {noise_model_b[3]}, 139 .gradient_intercept = {noise_model_c[3]}, 140 .offset_intercept = {noise_model_d[3]}}}}}, 141 """) 142 text_file = open(os.path.join(log_path, 'noise_profile.cc'), 'w') 143 text_file.write('%s' % code) 144 text_file.close() 145 146 147class DngNoiseModel(its_base_test.ItsBaseTest): 148 """Create DNG noise model. 149 150 Captures RAW images with increasing analog gains to create the model. 151 def requires 'test' in name to actually run. 152 """ 153 154 def test_dng_noise_model_generation(self): 155 with its_session_utils.ItsSession( 156 device_id=self.dut.serial, 157 camera_id=self.camera_id, 158 hidden_physical_id=self.hidden_physical_id) as cam: 159 props = cam.get_camera_properties() 160 props = cam.override_with_hidden_physical_camera_props(props) 161 log_path = self.log_path 162 name_with_log_path = os.path.join(log_path, _NAME) 163 if self.hidden_physical_id: 164 camera_name = f'{self.camera_id}.{self.hidden_physical_id}' 165 else: 166 camera_name = self.camera_id 167 logging.info('Starting %s for camera %s', _NAME, camera_name) 168 169 # Get basic properties we need. 170 sens_min, sens_max = props['android.sensor.info.sensitivityRange'] 171 sens_max_analog = props['android.sensor.maxAnalogSensitivity'] 172 sens_max_meas = sens_max_analog 173 white_level = props['android.sensor.info.whiteLevel'] 174 175 logging.info('Sensitivity range: [%d, %d]', sens_min, sens_max) 176 logging.info('Max analog sensitivity: %d', sens_max_analog) 177 178 # Do AE to get a rough idea of where we are. 179 iso_ae, exp_ae, _, _, _ = cam.do_3a( 180 get_results=True, do_awb=False, do_af=False) 181 182 # Underexpose to get more data for low signal levels. 183 auto_e = iso_ae * exp_ae / _BRACKET_FACTOR 184 check_auto_exposure_targets(auto_e, sens_min, sens_max_meas, props) 185 186 # Focus at zero to intentionally blur the scene as much as possible. 187 f_dist = 0.0 188 189 # Start the sensitivities at the minimum. 190 iso = sens_min 191 samples = [[], [], [], []] 192 plots = [] 193 measured_models = [[], [], [], []] 194 color_plane_plots = {} 195 isos = [] 196 while int(round(iso)) <= sens_max_meas: 197 iso_int = int(round(iso)) 198 isos.append(iso_int) 199 logging.info('ISO %d', iso_int) 200 fig, [[plt_r, plt_gr], [plt_gb, plt_b]] = plt.subplots( 201 2, 2, figsize=(11, 11)) 202 fig.gca() 203 color_plane_plots[iso_int] = [plt_r, plt_gr, plt_gb, plt_b] 204 fig.suptitle('ISO %d' % iso_int, x=0.54, y=0.99) 205 for i, plot in enumerate(color_plane_plots[iso_int]): 206 plot.set_title('%s' % _BAYER_LIST[i]) 207 plot.set_xlabel('Mean signal level') 208 plot.set_ylabel('Variance') 209 210 samples_s = [[], [], [], []] 211 for b in range(_BRACKET_MAX): 212 # Get the exposure for this sensitivity and exposure time. 213 exposure = int(math.pow(2, b)*auto_e/iso) 214 logging.info('exp %.3fms', round(exposure*1.0E-6, 3)) 215 req = capture_request_utils.manual_capture_request(iso_int, exposure, 216 f_dist) 217 fmt_raw = {'format': 'rawStats', 218 'gridWidth': _TILE_SIZE, 219 'gridHeight': _TILE_SIZE} 220 cap = cam.do_capture(req, fmt_raw) 221 if self.debug_mode: 222 img = image_processing_utils.convert_capture_to_rgb_image( 223 cap, props=props) 224 image_processing_utils.write_image( 225 img, f'{name_with_log_path}_{iso_int}_{exposure}ns.jpg', True) 226 227 mean_img, var_img = image_processing_utils.unpack_rawstats_capture( 228 cap) 229 idxs = image_processing_utils.get_canonical_cfa_order(props) 230 raw_stats_size = mean_img.shape 231 means = [mean_img[_TILE_CROP_N:raw_stats_size[0]-_TILE_CROP_N, 232 _TILE_CROP_N:raw_stats_size[1]-_TILE_CROP_N, i] 233 for i in idxs] 234 vars_ = [var_img[_TILE_CROP_N:raw_stats_size[0]-_TILE_CROP_N, 235 _TILE_CROP_N:raw_stats_size[1]-_TILE_CROP_N, i] 236 for i in idxs] 237 if self.debug_mode: 238 logging.info('rawStats image size: %s', str(raw_stats_size)) 239 logging.info('R planes means image size: %s', str(means[0].shape)) 240 logging.info('means min: %.3f, median: %.3f, max: %.3f', 241 np.min(means), np.median(means), np.max(means)) 242 logging.info('vars_ min: %.4f, median: %.4f, max: %.4f', 243 np.min(vars_), np.median(vars_), np.max(vars_)) 244 245 s_read = cap['metadata']['android.sensor.sensitivity'] 246 if not 1.0 >= s_read/float(iso_int) >= _RTOL_EXP_GAIN: 247 raise AssertionError( 248 f's_write: {iso}, s_read: {s_read}, RTOL: {_RTOL_EXP_GAIN}') 249 logging.info('ISO_write: %d, ISO_read: %d', iso_int, s_read) 250 251 for pidx in range(len(means)): 252 plot = color_plane_plots[iso_int][pidx] 253 254 # convert_capture_to_planes normalizes the range to [0, 1], but 255 # without subtracting the black level. 256 black_level = image_processing_utils.get_black_level( 257 pidx, props, cap['metadata']) 258 means_p = (means[pidx] - black_level)/(white_level - black_level) 259 vars_p = vars_[pidx]/((white_level - black_level)**2) 260 261 # TODO(dsharlet): It should be possible to account for low 262 # frequency variation by looking at neighboring means, but I 263 # have not been able to make this work. 264 265 means_p = np.asarray(means_p).flatten() 266 vars_p = np.asarray(vars_p).flatten() 267 268 samples_e = [] 269 for (mean, var) in zip(means_p, vars_p): 270 # Don't include the tile if it has samples that might be clipped. 271 if mean + 2*math.sqrt(max(var, 0)) < _MAX_SIGNAL_VALUE: 272 samples_e.append([mean, var]) 273 274 if samples_e: 275 means_e, vars_e = zip(*samples_e) 276 color_plane_plots[iso_int][pidx].plot( 277 means_e, vars_e, _PLOT_COLORS[b%len(_PLOT_COLORS)] + '.', 278 alpha=0.5, markersize=1) 279 samples_s[pidx].extend(samples_e) 280 281 for (pidx, p) in enumerate(samples_s): 282 [slope, intercept, rvalue, _, _] = scipy.stats.linregress( 283 samples_s[pidx]) 284 measured_models[pidx].append([iso_int, slope, intercept]) 285 logging.info('%s sensitivity %d: %e*y + %e (R=%f)', 286 'RGKB'[pidx], iso_int, slope, intercept, rvalue) 287 288 # Add the samples for this sensitivity to the global samples list. 289 samples[pidx].extend( 290 [(iso_int, mean, var) for (mean, var) in samples_s[pidx]]) 291 292 # Add the linear fit to subplot for this sensitivity. 293 color_plane_plots[iso_int][pidx].plot( 294 [0, _MAX_SIGNAL_VALUE], 295 [intercept, intercept + slope * _MAX_SIGNAL_VALUE], 296 'rgkb'[pidx] + '--', 297 label='Linear fit') 298 299 xmax = max([max([x for (x, _) in p]) for p in samples_s 300 ]) * _MAX_SCALE_FUDGE 301 ymax = (intercept + slope * xmax) * _MAX_SCALE_FUDGE 302 color_plane_plots[iso_int][pidx].set_xlim(xmin=0, xmax=xmax) 303 color_plane_plots[iso_int][pidx].set_ylim(ymin=0, ymax=ymax) 304 color_plane_plots[iso_int][pidx].legend() 305 pylab.tight_layout() 306 307 fig.savefig(f'{name_with_log_path}_samples_iso{iso_int:04d}.png') 308 plots.append([iso_int, fig]) 309 310 # Move to the next sensitivity. 311 iso *= math.pow(2, 1.0/_STEPS_PER_STOP) 312 313 # Do model plots 314 (fig, (plt_s, plt_o)) = plt.subplots(2, 1, figsize=(11, 8.5)) 315 plt_s.set_title('Noise model: N(x) = sqrt(Sx + O)') 316 plt_s.set_ylabel('S') 317 plt_o.set_xlabel('ISO') 318 plt_o.set_ylabel('O') 319 320 noise_model = [] 321 for (pidx, p) in enumerate(measured_models): 322 # Grab the sensitivities and line parameters from each sensitivity. 323 s_measured = [e[1] for e in measured_models[pidx]] 324 o_measured = [e[2] for e in measured_models[pidx]] 325 sens = np.asarray([e[0] for e in measured_models[pidx]]) 326 sens_sq = np.square(sens) 327 328 # Use a global linear optimization to fit the noise model. 329 gains = np.asarray([s[0] for s in samples[pidx]]) 330 means = np.asarray([s[1] for s in samples[pidx]]) 331 vars_ = np.asarray([s[2] for s in samples[pidx]]) 332 gains = gains.flatten() 333 means = means.flatten() 334 vars_ = vars_.flatten() 335 336 # Define digital gain as the gain above the max analog gain 337 # per the Camera2 spec. Also, define a corresponding C 338 # expression snippet to use in the generated model code. 339 digital_gains = np.maximum(gains/sens_max_analog, 1) 340 if not np.all(digital_gains == 1): 341 raise AssertionError(f'Digital gain! gains: {gains}, ' 342 f'Max analog gain: {sens_max_analog}.') 343 digital_gain_cdef = '(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)' % ( 344 sens_max_analog, sens_max_analog) 345 346 # Divide the whole system by gains*means. 347 f = lambda x, a, b, c, d: (c*(x[0]**2)+d+(x[1])*a*x[0]+(x[1])*b)/(x[0]) 348 result, _ = scipy.optimize.curve_fit(f, (gains, means), vars_/(gains)) 349 # result[0:4] = s_gradient, s_offset, o_gradient, o_offset 350 # Note 'S' and 'O' are the API terms for the 2 model params. 351 # The noise_profile.cc uses 'slope' for 'S' and 'intercept' for 'O'. 352 # 'gradient' and 'offset' are used to describe the linear fit 353 # parameters for 'S' and 'O'. 354 noise_model.append(result[0:4]) 355 356 # Plot noise model components with the values predicted by the model. 357 s_model = result[0]*sens + result[1] 358 o_model = result[2]*sens_sq + result[3]*np.square(np.maximum( 359 sens/sens_max_analog, 1)) 360 361 plt_s.loglog(sens, s_measured, 'rgkb'[pidx]+'+', base=10, 362 label='Measured') 363 plt_s.loglog(sens, s_model, 'rgkb'[pidx]+'o', base=10, 364 label='Model', alpha=0.3) 365 plt_o.loglog(sens, o_measured, 'rgkb'[pidx]+'+', base=10, 366 label='Measured') 367 plt_o.loglog(sens, o_model, 'rgkb'[pidx]+'o', base=10, 368 label='Model', alpha=0.3) 369 plt_s.legend() 370 plt_s.set_xticks(isos) 371 plt_s.set_xticklabels(isos) 372 373 plt_o.set_xticks([]) 374 plt_o.set_xticks(isos) 375 plt_o.set_xticklabels(isos) 376 plt_o.legend() 377 fig.savefig(f'{name_with_log_path}.png') 378 379 # Generate individual noise model components 380 noise_model_a, noise_model_b, noise_model_c, noise_model_d = zip( 381 *noise_model) 382 383 # Add models to subplots and re-save 384 for [iso, fig] in plots: # re-step through figs... 385 dig_gain = max(iso/sens_max_analog, 1) 386 fig.gca() 387 for (pidx, p) in enumerate(measured_models): 388 s = noise_model_a[pidx]*iso + noise_model_b[pidx] 389 o = noise_model_c[pidx]*iso**2 + noise_model_d[pidx]*dig_gain**2 390 color_plane_plots[iso][pidx].plot( 391 [0, _MAX_SIGNAL_VALUE], [o, o+s*_MAX_SIGNAL_VALUE], 392 'rgkb'[pidx]+'-', label='Model', alpha=0.5) 393 color_plane_plots[iso][pidx].legend(loc='upper left') 394 fig.savefig(f'{name_with_log_path}_samples_iso{iso:04d}.png') 395 396 # Validity checks on model: read noise > 0, positive intercept gradient. 397 for i, _ in enumerate(_BAYER_LIST): 398 read_noise = noise_model_c[i] * sens_min * sens_min + noise_model_d[i] 399 if read_noise <= 0: 400 raise AssertionError(f'{_BAYER_LIST[i]} model min ISO noise < 0! ' 401 f'API intercept gradient: {noise_model_c[i]:.4e}, ' 402 f'API intercept offset: {noise_model_d[i]:.4e}, ' 403 f'read_noise: {read_noise:.4e}') 404 if noise_model_c[i] <= 0: 405 raise AssertionError(f'{_BAYER_LIST[i]} model API intercept gradient ' 406 f'is negative: {noise_model_c[i]:.4e}') 407 408 # Generate the noise model file. 409 create_noise_model_code( 410 noise_model_a, noise_model_b, noise_model_c, noise_model_d, 411 sens_min, sens_max, digital_gain_cdef, log_path) 412 413if __name__ == '__main__': 414 test_runner.main() 415