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