ITS: update tools/dng_noise_model.py with improvements
Use rawStats instead of raw captures to speed things up.
Use a smaller tile size to reduce low freq image changes
Increase exposure bracketing stops from 4 to 8
Minor gpylint cleanups for alignment and spacing.
bug: 158028170
Change-Id: I6e2007bb8f350071b719af68a051858dc7c80af1
diff --git a/apps/CameraITS/tools/dng_noise_model.py b/apps/CameraITS/tools/dng_noise_model.py
index b490d19..db37deb 100644
--- a/apps/CameraITS/tools/dng_noise_model.py
+++ b/apps/CameraITS/tools/dng_noise_model.py
@@ -1,4 +1,4 @@
-# Copyright 2014 The Android Open Source Project
+# Copyright 2014 The Android Open Source Project. Lint as python2
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
@@ -21,7 +21,6 @@
import its.image
import its.objects
-import matplotlib
from matplotlib import pylab
import matplotlib.pyplot as plt
import numpy as np
@@ -29,61 +28,24 @@
import scipy.stats
BAYER_LIST = ['R', 'GR', 'GB', 'B']
+BRACKET_MAX = 8 # Exposure bracketing range in stops
+COLORS = 'rygcbm' # Colors used for plotting the data for each exposure.
+MAX_SCALE_FUDGE = 1.1
+MAX_SIGNAL_VALUE = 0.25 # Maximum value to allow mean of the tiles to go.
NAME = os.path.basename(__file__).split('.')[0]
RTOL_EXP_GAIN = 0.97
-
-
-def tile(a, tile_size):
- """Convert a 2D array to 4D w/ dims [tile_size, tile_size, row, col] where row, col are tile indices.
- """
- tile_rows, tile_cols = a.shape[0]/tile_size, a.shape[1]/tile_size
- a = a.reshape([tile_rows, tile_size, tile_cols, tile_size])
- a = a.transpose([1, 3, 0, 2])
- return a
+STEPS_PER_STOP = 2 # How many sensitivities per stop to sample.
+# How large of tiles to use to compute mean/variance.
+# Large tiles may have their variance corrupted by low frequency
+# image changes (lens shading, scene illumination).
+TILE_SIZE = 32
def main():
"""Capture a set of raw images with increasing analog gains and measure the noise.
"""
- # How many sensitivities per stop to sample.
- steps_per_stop = 2
- # How large of tiles to use to compute mean/variance.
- tile_size = 64
- # Exposure bracketing range in stops
- bracket_stops = 4
- # How high to allow the mean of the tiles to go.
- max_signal_level = 0.25
- # Colors used for plotting the data for each exposure.
- colors = 'rygcbm'
-
- # Define a first order high pass filter to eliminate low frequency
- # signal content when computing variance.
- f = np.array([-1, 1]).astype('float32')
- # Make it a higher order filter by convolving the first order
- # filter with itself a few times.
- f = np.convolve(f, f)
- f = np.convolve(f, f)
-
- # Compute the normalization of the filter to preserve noise
- # power. Let a be the normalization factor we're looking for, and
- # Let X and X' be the random variables representing the noise
- # before and after filtering, respectively. First, compute
- # Var[a*X']:
- #
- # Var[a*X'] = a^2*Var[X*f_0 + X*f_1 + ... + X*f_N-1]
- # = a^2*(f_0^2*Var[X] + f_1^2*Var[X] + ... + (f_N-1)^2*Var[X])
- # = sum(f_i^2)*a^2*Var[X]
- #
- # We want Var[a*X'] to be equal to Var[X]:
- #
- # sum(f_i^2)*a^2*Var[X] = Var[X] -> a = sqrt(1/sum(f_i^2))
- #
- # We can just bake this normalization factor into the high pass
- # filter kernel.
- f /= math.sqrt(np.dot(f, f))
-
- bracket_factor = math.pow(2, bracket_stops)
+ bracket_factor = math.pow(2, BRACKET_MAX)
with its.device.ItsSession() as cam:
props = cam.get_camera_properties()
@@ -95,28 +57,28 @@
sens_max_meas = sens_max_analog
white_level = props['android.sensor.info.whiteLevel']
- print "Sensitivity range: [%f, %f]" % (sens_min, sens_max)
- print "Max analog sensitivity: %f" % (sens_max_analog)
+ print 'Sensitivity range: [%f, %f]' % (sens_min, sens_max)
+ print 'Max analog sensitivity: %f' % (sens_max_analog)
# Do AE to get a rough idea of where we are.
- s_ae, e_ae, _, _, _ = \
- cam.do_3a(get_results=True, do_awb=False, do_af=False)
+ s_ae, e_ae, _, _, _ = cam.do_3a(
+ get_results=True, do_awb=False, do_af=False)
# Underexpose to get more data for low signal levels.
- auto_e = s_ae*e_ae/bracket_factor
+ auto_e = s_ae * e_ae / bracket_factor
# Focus at zero to intentionally blur the scene as much as possible.
f_dist = 0.0
# If the auto-exposure result is too bright for the highest
# sensitivity or too dark for the lowest sensitivity, report
# an error.
- min_exposure_ns, max_exposure_ns = \
- props['android.sensor.info.exposureTimeRange']
+ min_exposure_ns, max_exposure_ns = props[
+ 'android.sensor.info.exposureTimeRange']
if auto_e < min_exposure_ns*sens_max_meas:
- raise its.error.Error("Scene is too bright to properly expose \
- at the highest sensitivity")
+ raise its.error.Error('Scene is too bright to properly expose '
+ 'at the highest sensitivity')
if auto_e*bracket_factor > max_exposure_ns*sens_min:
- raise its.error.Error("Scene is too dark to properly expose \
- at the lowest sensitivity")
+ raise its.error.Error('Scene is too dark to properly expose '
+ 'at the lowest sensitivity')
# Start the sensitivities at the minimum.
s = sens_min
@@ -138,57 +100,56 @@
plot.set_ylabel('Variance')
samples_s = [[], [], [], []]
- for b in range(0, bracket_stops):
+ for b in range(BRACKET_MAX):
# Get the exposure for this sensitivity and exposure time.
e = int(math.pow(2, b)*auto_e/float(s))
print 'exp %.3fms' % round(e*1.0E-6, 3)
req = its.objects.manual_capture_request(s_int, e, f_dist)
- cap = cam.do_capture(req, cam.CAP_RAW)
- planes = its.image.convert_capture_to_planes(cap, props)
- e_read = cap['metadata']['android.sensor.exposureTime']
+ fmt_raw = {'format': 'rawStats',
+ 'gridWidth': TILE_SIZE,
+ 'gridHeight': TILE_SIZE}
+ cap = cam.do_capture(req, fmt_raw)
+ mean_image, var_image = its.image.unpack_rawstats_capture(cap)
+ idxs = its.image.get_canonical_cfa_order(props)
+ means = [mean_image[:, :, i] for i in idxs]
+ vars_ = [var_image[:, :, i] for i in idxs]
+
s_read = cap['metadata']['android.sensor.sensitivity']
s_err = 's_write: %d, s_read: %d, RTOL: %.2f' % (
s, s_read, RTOL_EXP_GAIN)
assert (1.0 >= s_read/float(s_int) >= RTOL_EXP_GAIN), s_err
print 'ISO_write: %d, ISO_read: %d' % (s_int, s_read)
- for (pidx, p) in enumerate(planes):
+ for pidx in range(len(means)):
plot = color_plane_plots[s_int][pidx]
- p = p.squeeze()
-
- # Crop the plane to be a multiple of the tile size.
- p = p[0:p.shape[0] - p.shape[0]%tile_size,
- 0:p.shape[1] - p.shape[1]%tile_size]
-
# convert_capture_to_planes normalizes the range
# to [0, 1], but without subtracting the black
# level.
black_level = its.image.get_black_level(
- pidx, props, cap["metadata"])
- p *= white_level
- p = (p - black_level)/(white_level - black_level)
+ pidx, props, cap['metadata'])
+ means_p = (means[pidx] - black_level)/(white_level - black_level)
+ vars_p = vars_[pidx]/((white_level - black_level)**2)
- # Use our high pass filter to filter this plane.
- hp = scipy.signal.sepfir2d(p, f, f).astype('float32')
+ # TODO(dsharlet): It should be possible to account for low
+ # frequency variation by looking at neighboring means, but I
+ # have not been able to make this work.
- means_tiled = \
- np.mean(tile(p, tile_size), axis=(0, 1)).flatten()
- vars_tiled = \
- np.var(tile(hp, tile_size), axis=(0, 1)).flatten()
+ means_p = np.asarray(means_p).flatten()
+ vars_p = np.asarray(vars_p).flatten()
samples_e = []
- for (mean, var) in zip(means_tiled, vars_tiled):
+ for (mean, var) in zip(means_p, vars_p):
# Don't include the tile if it has samples that might
# be clipped.
- if mean + 2*math.sqrt(var) < max_signal_level:
+ if mean + 2*math.sqrt(max(var, 0)) < MAX_SIGNAL_VALUE:
samples_e.append([mean, var])
if samples_e:
means_e, vars_e = zip(*samples_e)
color_plane_plots[s_int][pidx].plot(
- means_e, vars_e, colors[b%len(colors)] + '.',
- alpha=0.5)
+ means_e, vars_e, COLORS[b%len(COLORS)] + '.',
+ alpha=0.5, markersize=1)
samples_s[pidx].extend(samples_e)
for (pidx, p) in enumerate(samples_s):
@@ -200,15 +161,12 @@
samples[pidx].extend([(s_int, mean, var) for (mean, var) in samples_s[pidx]])
# Add the linear fit to subplot for this sensitivity.
- # pylab.subplot(2, 2, pidx+1)
- #pylab.plot([0, max_signal_level], [O, O + S*max_signal_level], 'rgkb'[pidx]+'--',
- #label="Linear fit")
- color_plane_plots[s_int][pidx].plot([0, max_signal_level],
- [O, O + S*max_signal_level], 'rgkb'[pidx]+'--',
- label="Linear fit")
+ color_plane_plots[s_int][pidx].plot(
+ [0, MAX_SIGNAL_VALUE], [O, O + S*MAX_SIGNAL_VALUE],
+ 'rgkb'[pidx]+'--', label='Linear fit')
- xmax = max([max([x for (x, _) in p]) for p in samples_s])*1.25
- ymax = max([max([y for (_, y) in p]) for p in samples_s])*1.25
+ xmax = max([max([x for (x, _) in p]) for p in samples_s])*MAX_SCALE_FUDGE
+ ymax = (O + S*xmax)*MAX_SCALE_FUDGE
color_plane_plots[s_int][pidx].set_xlim(xmin=0, xmax=xmax)
color_plane_plots[s_int][pidx].set_ylim(ymin=0, ymax=ymax)
color_plane_plots[s_int][pidx].legend()
@@ -218,7 +176,7 @@
plots.append([s_int, fig])
# Move to the next sensitivity.
- s *= math.pow(2, 1.0/steps_per_stop)
+ s *= math.pow(2, 1.0/STEPS_PER_STOP)
# do model plots
(fig, (plt_S, plt_O)) = plt.subplots(2, 1, figsize=(11, 8.5))
@@ -242,28 +200,23 @@
gains = np.asarray([s[0] for s in samples[pidx]])
means = np.asarray([s[1] for s in samples[pidx]])
vars_ = np.asarray([s[2] for s in samples[pidx]])
+ gains = gains.flatten()
+ means = means.flatten()
+ vars_ = vars_.flatten()
# Define digital gain as the gain above the max analog gain
# per the Camera2 spec. Also, define a corresponding C
# expression snippet to use in the generated model code.
digital_gains = np.maximum(gains/sens_max_analog, 1)
- digital_gain_cdef = "(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)" % \
+ assert np.all(digital_gains == 1)
+ digital_gain_cdef = '(sens / %d.0) < 1.0 ? 1.0 : (sens / %d.0)' % \
(sens_max_analog, sens_max_analog)
- # Find the noise model parameters via least squares fit.
- ad = gains*means
- bd = means
- cd = gains*gains
- dd = digital_gains*digital_gains
- a = np.asarray([ad, bd, cd, dd]).T
- b = vars_
+ # Divide the whole system by gains*means.
+ f = lambda x, a, b, c, d: (c*(x[0]**2) + d + (x[1])*a*x[0] + (x[1])*b)/(x[0])
+ [result, _] = scipy.optimize.curve_fit(f, (gains, means), vars_/(gains))
- # To avoid overfitting to high ISOs (high variances), divide the system
- # by the gains.
- a /= (np.tile(gains, (a.shape[1], 1)).T)
- b /= gains
-
- [A_p, B_p, C_p, D_p], _, _, _ = np.linalg.lstsq(a, b)
+ [A_p, B_p, C_p, D_p] = result[0:4]
A.append(A_p)
B.append(B_p)
C.append(C_p)
@@ -276,16 +229,17 @@
C_p*sens_sq + D_p*np.square(np.maximum(sens/sens_max_analog, 1))
plt_S.loglog(sens, S_measured, 'rgkb'[pidx]+'+', basex=10, basey=10,
- label="Measured")
- plt_S.loglog(sens, S_model, 'rgkb'[pidx]+'x', basex=10, basey=10, label="Model")
-
+ label='Measured')
+ plt_S.loglog(sens, S_model, 'rgkb'[pidx]+'x', basex=10, basey=10,
+ label='Model')
plt_O.loglog(sens, O_measured, 'rgkb'[pidx]+'+', basex=10, basey=10,
- label="Measured")
- plt_O.loglog(sens, O_model, 'rgkb'[pidx]+'x', basex=10, basey=10, label="Model")
+ label='Measured')
+ plt_O.loglog(sens, O_model, 'rgkb'[pidx]+'x', basex=10, basey=10,
+ label='Model')
plt_S.legend()
plt_O.legend()
- fig.savefig("%s.png" % (NAME))
+ fig.savefig('%s.png' % (NAME))
# add models to subplots and re-save
for [s, fig] in plots: # re-step through figs...
@@ -294,9 +248,9 @@
for (pidx, p) in enumerate(measured_models):
S = A[pidx]*s + B[pidx]
O = C[pidx]*s*s + D[pidx]*dg*dg
- color_plane_plots[s][pidx].plot([0, max_signal_level],
- [O, O + S*max_signal_level], 'rgkb'[pidx]+'-',
- label="Model", alpha=0.5)
+ color_plane_plots[s][pidx].plot(
+ [0, MAX_SIGNAL_VALUE], [O, O + S*MAX_SIGNAL_VALUE],
+ 'rgkb'[pidx]+'-', label='Model', alpha=0.5)
color_plane_plots[s][pidx].legend(loc='upper left')
fig.savefig('%s_samples_iso%04d.png' % (NAME, s))
