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 its.device 16import its.objects 17import its.image 18import pprint 19import pylab 20import os.path 21import matplotlib 22import matplotlib.pyplot 23import numpy 24import math 25 26def main(): 27 """Compute the DNG noise model from a color checker chart. 28 29 TODO: Make this more robust; some manual futzing may be needed. 30 """ 31 NAME = os.path.basename(__file__).split(".")[0] 32 33 with its.device.ItsSession() as cam: 34 35 props = cam.get_camera_properties() 36 37 white_level = float(props['android.sensor.info.whiteLevel']) 38 black_levels = props['android.sensor.blackLevelPattern'] 39 idxs = its.image.get_canonical_cfa_order(props) 40 black_levels = [black_levels[i] for i in idxs] 41 42 # Expose for the scene with min sensitivity 43 sens_min, sens_max = props['android.sensor.info.sensitivityRange'] 44 s_ae,e_ae,awb_gains,awb_ccm,_ = cam.do_3a(get_results=True) 45 s_e_prod = s_ae * e_ae 46 47 # Make the image brighter since the script looks at linear Bayer 48 # raw patches rather than gamma-encoded YUV patches (and the AE 49 # probably under-exposes a little for this use-case). 50 s_e_prod *= 2 51 52 # Capture raw frames across the full sensitivity range. 53 NUM_SENS_STEPS = 15 54 sens_step = int((sens_max - sens_min - 1) / float(NUM_SENS_STEPS)) 55 reqs = [] 56 sens = [] 57 for s in range(sens_min, sens_max, sens_step): 58 e = int(s_e_prod / float(s)) 59 req = its.objects.manual_capture_request(s, e) 60 req["android.colorCorrection.transform"] = \ 61 its.objects.float_to_rational(awb_ccm) 62 req["android.colorCorrection.gains"] = awb_gains 63 reqs.append(req) 64 sens.append(s) 65 66 caps = cam.do_capture(reqs, cam.CAP_RAW) 67 68 # A list of the (x,y) coords of the center pixel of a collection of 69 # patches of a color checker chart. Each patch should be uniform, 70 # however the actual color doesn't matter. Note that the coords are 71 # relative to the *converted* RGB image, which is 1/2 x 1/2 of the 72 # full size; convert back to full. 73 img = its.image.convert_capture_to_rgb_image(caps[0], props=props) 74 patches = its.image.get_color_checker_chart_patches(img, NAME+"_debug") 75 patches = [(2*x,2*y) for (x,y) in sum(patches,[])] 76 77 lines = [] 78 for (s,cap) in zip(sens,caps): 79 # For each capture, compute the mean value in each patch, for each 80 # Bayer plane; discard patches where pixels are close to clamped. 81 # Also compute the variance. 82 CLAMP_THRESH = 0.2 83 planes = its.image.convert_capture_to_planes(cap, props) 84 points = [] 85 for i,plane in enumerate(planes): 86 plane = (plane * white_level - black_levels[i]) / ( 87 white_level - black_levels[i]) 88 for j,(x,y) in enumerate(patches): 89 tile = plane[y/2-16:y/2+16:,x/2-16:x/2+16:,::] 90 mean = its.image.compute_image_means(tile)[0] 91 var = its.image.compute_image_variances(tile)[0] 92 if (mean > CLAMP_THRESH and mean < 1.0-CLAMP_THRESH): 93 # Each point is a (mean,variance) tuple for a patch; 94 # for a given ISO, there should be a linear 95 # relationship between these values. 96 points.append((mean,var)) 97 98 # Fit a line to the points, with a line equation: y = mx + b. 99 # This line is the relationship between mean and variance (i.e.) 100 # between signal level and noise, for this particular sensor. 101 # In the DNG noise model, the gradient (m) is "S", and the offset 102 # (b) is "O". 103 points.sort() 104 xs = [x for (x,y) in points] 105 ys = [y for (x,y) in points] 106 m,b = numpy.polyfit(xs, ys, 1) 107 lines.append((s,m,b)) 108 print s, "->", m, b 109 110 # TODO: Clean up these checks (which currently fail in some cases). 111 # Some sanity checks: 112 # * Noise levels should increase with brightness. 113 # * Extrapolating to a black image, the noise should be positive. 114 # Basically, the "b" value should correspnd to the read noise, 115 # which is the noise level if the sensor was operating in zero 116 # light. 117 #assert(m > 0) 118 #assert(b >= 0) 119 120 # Draw a plot. 121 pylab.plot(xs, ys, 'r') 122 pylab.plot([0,xs[-1]],[b,m*xs[-1]+b],'b') 123 matplotlib.pyplot.savefig("%s_plot_mean_vs_variance.png" % (NAME)) 124 125 # Now fit a line across the (m,b) line parameters for each sensitivity. 126 # The gradient (m) params are fit to the "S" line, and the offset (b) 127 # params are fit to the "O" line, both as a function of sensitivity. 128 gains = [d[0] for d in lines] 129 Ss = [d[1] for d in lines] 130 Os = [d[2] for d in lines] 131 mS,bS = numpy.polyfit(gains, Ss, 1) 132 mO,bO = numpy.polyfit(gains, Os, 1) 133 134 # Plot curve "O" as 10x, so it fits in the same scale as curve "S". 135 pylab.plot(gains, [10*o for o in Os], 'r') 136 pylab.plot([gains[0],gains[-1]], 137 [10*mO*gains[0]+10*bO, 10*mO*gains[-1]+10*bO], 'b') 138 pylab.plot(gains, Ss, 'r') 139 pylab.plot([gains[0],gains[-1]], [mS*gains[0]+bS, mS*gains[-1]+bS], 'b') 140 matplotlib.pyplot.savefig("%s_plot_S_O.png" % (NAME)) 141 142 print """ 143 /* Generated test code to dump a table of data for external validation 144 * of the noise model parameters. 145 */ 146 #include <stdio.h> 147 #include <assert.h> 148 double compute_noise_model_entry_S(int sens); 149 double compute_noise_model_entry_O(int sens); 150 int main(void) { 151 int sens; 152 for (sens = %d; sens <= %d; sens += 100) { 153 double o = compute_noise_model_entry_O(sens); 154 double s = compute_noise_model_entry_S(sens); 155 printf("%%d,%%lf,%%lf\\n", sens, o, s); 156 } 157 return 0; 158 } 159 160 /* Generated functions to map a given sensitivity to the O and S noise 161 * model parameters in the DNG noise model. 162 */ 163 double compute_noise_model_entry_S(int sens) { 164 double s = %e * sens + %e; 165 return s < 0.0 ? 0.0 : s; 166 } 167 double compute_noise_model_entry_O(int sens) { 168 double o = %e * sens + %e; 169 return o < 0.0 ? 0.0 : o; 170 } 171 """%(sens_min,sens_max,mS,bS,mO,bO) 172 173if __name__ == '__main__': 174 main() 175 176