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