CameraITS: Added DNG noise model tool.

Also added method to its.image to extract center coords of patches
in a color checker chart.

Change-Id: I2a4a91861234b71347d552d63be9baa29b08c429

2 files changed, 345 insertions, 0 deletions

diff --git a/apps/CameraITS/pymodules/its/image.py b/apps/CameraITS/pymodules/its/image.py
index b96896a..14bd248 100644
--- a/apps/CameraITS/pymodules/its/image.py
+++ b/apps/CameraITS/pymodules/its/image.py
@@ -23,6 +23,7 @@ import numpy
 import math
 import unittest
 import cStringIO
+import scipy.stats
 
 DEFAULT_YUV_TO_RGB_CCM = numpy.matrix([
                                 [1.000,  0.000,  1.402],
@@ -193,6 +194,7 @@ def convert_raw_to_rgb_image(r_plane, gr_plane, gb_plane, b_plane,
         RGB float-3 image array, with pixel values in [0.0, 1.0]
     """
     # Values required for the RAW to RGB conversion.
+    assert(props is not None)
     white_level = float(props['android.sensor.info.whiteLevel'])
     black_levels = props['android.sensor.blackLevelPattern']
     gains = cap_res['android.colorCorrection.gains']
@@ -479,6 +481,173 @@ def downscale_image(img, f):
     img = numpy.vstack(chs).T.reshape(h/f,w/f,chans)
     return img
 
+def __measure_color_checker_patch(img, xc,yc, patch_size):
+    r = patch_size/2
+    tile = img[yc-r:yc+r+1:, xc-r:xc+r+1:, ::]
+    means = tile.mean(1).mean(0)
+    return means
+
+def get_color_checker_chart_patches(img, debug_fname_prefix=None):
+    """Return the center coords of each patch in a color checker chart.
+
+    Assumptions:
+    * Chart is vertical or horizontal w.r.t. camera, but not diagonal.
+    * Chart is (roughly) planar-parallel to the camera.
+    * Chart is centered in frame (roughly).
+    * Around/behind chart is white/grey background.
+    * The only black pixels in the image are from the chart.
+    * Chart is 100% visible and contained within image.
+    * No other objects within image.
+    * Image is well-exposed.
+    * Standard color checker chart with standard-sized black borders.
+
+    The values returned are in the coordinate system of the chart; that is,
+    the "origin" patch is the brown patch that is in the chart's top-left
+    corner when it is in the normal upright/horizontal orientation. (The chart
+    may be any of the four main orientations in the image.)
+
+    The chart is 6x4 patches in the normal upright orientation. The return
+    values of this function are the center coordinate of the top-left patch,
+    and the displacement vectors to the next patches to the right and below
+    the top-left patch. From these pieces of data, the center coordinates of
+    any of the patches can be computed.
+
+    Args:
+        img: Input image, as a numpy array with pixels in [0,1].
+        debug_fname_prefix: If not None, the (string) name of a file prefix to
+            use to save a number of debug images for visulaizing the output of
+            this function; can be used to see if the patches are being found
+            successfully.
+
+    Returns:
+        Three tuples:
+        (1) The integer (x,y) coords of the center of the chart's patch (0,0).
+        (2) The float (dx,dy) displacement to the chart's (0,1) patch, which
+            is the pink patch to the right of the top-left brown patch.
+        (3) The float (dx,dt) displacement to the chart's (1,0) patch, which
+            is the orange patch below the top-left brown patch.
+    """
+
+    # Shrink the original image.
+    DOWNSCALE_FACTOR = 4
+    img_small = downscale_image(img, DOWNSCALE_FACTOR)
+
+    # Make a threshold image, which is 1.0 where the image is black,
+    # and 0.0 elsewhere.
+    BLACK_PIXEL_THRESH = 0.2
+    mask_img = scipy.stats.threshold(
+                img_small.max(2), BLACK_PIXEL_THRESH, 1.1, 0.0)
+    mask_img = 1.0 - scipy.stats.threshold(mask_img, -0.1, 0.1, 1.0)
+
+    if debug_fname_prefix is not None:
+        h,w = mask_img.shape
+        write_image(img, debug_fname_prefix+"_0.jpg")
+        write_image(mask_img.repeat(3).reshape(h,w,3),
+                debug_fname_prefix+"_1.jpg")
+
+    # Mask image flattened to a single row or column (by averaging).
+    # Also apply a threshold to these arrays.
+    FLAT_PIXEL_THRESH = 0.05
+    flat_row = mask_img.mean(0)
+    flat_col = mask_img.mean(1)
+    flat_row = [0 if v < FLAT_PIXEL_THRESH else 1 for v in flat_row]
+    flat_col = [0 if v < FLAT_PIXEL_THRESH else 1 for v in flat_col]
+
+    # Start and end of the non-zero region of the flattened row/column.
+    flat_row_nonzero = [i for i in range(len(flat_row)) if flat_row[i]>0]
+    flat_col_nonzero = [i for i in range(len(flat_col)) if flat_col[i]>0]
+    flat_row_min, flat_row_max = min(flat_row_nonzero), max(flat_row_nonzero)
+    flat_col_min, flat_col_max = min(flat_col_nonzero), max(flat_col_nonzero)
+
+    # Orientation of chart, and number of grid cells horz. and vertically.
+    orient = "h" if flat_row_max-flat_row_min>flat_col_max-flat_col_min else "v"
+    xgrids = 6 if orient=="h" else 4
+    ygrids = 6 if orient=="v" else 4
+
+    # Get better bounds on the patches region, lopping off some of the excess
+    # black border.
+    HRZ_BORDER_PAD_FRAC = 0.0138
+    VERT_BORDER_PAD_FRAC = 0.0395
+    xpad = HRZ_BORDER_PAD_FRAC if orient=="h" else VERT_BORDER_PAD_FRAC
+    ypad = HRZ_BORDER_PAD_FRAC if orient=="v" else VERT_BORDER_PAD_FRAC
+    xchart = flat_row_min + (flat_row_max - flat_row_min) * xpad
+    ychart = flat_col_min + (flat_col_max - flat_col_min) * ypad
+    wchart = (flat_row_max - flat_row_min) * (1 - 2*xpad)
+    hchart = (flat_col_max - flat_col_min) * (1 - 2*ypad)
+
+    # Get the colors of the 4 corner patches, in clockwise order, by measuring
+    # the average value of a small patch at each of the 4 patch centers.
+    colors = []
+    centers = []
+    for (x,y) in [(0,0), (xgrids-1,0), (xgrids-1,ygrids-1), (0,ygrids-1)]:
+        xc = xchart + (x + 0.5)*wchart/xgrids
+        yc = ychart + (y + 0.5)*hchart/ygrids
+        xc = int(xc * DOWNSCALE_FACTOR + 0.5)
+        yc = int(yc * DOWNSCALE_FACTOR + 0.5)
+        centers.append((xc,yc))
+        chan_means = __measure_color_checker_patch(img, xc,yc, 32)
+        colors.append(sum(chan_means) / len(chan_means))
+
+    # The brightest corner is the white patch, the darkest is the black patch.
+    # The black patch should be counter-clockwise from the white patch.
+    white_patch_index = None
+    for i in range(4):
+        if colors[i] == max(colors) and \
+                colors[(i-1+4)%4] == min(colors):
+            white_patch_index = i%4
+    assert(white_patch_index is not None)
+
+    # Return the coords of the origin (top-left when the chart is in the normal
+    # upright orientation) patch's center, and the vector displacement to the
+    # center of the second patch on the first row of the chart (when in the
+    # normal upright orienation).
+    origin_index = (white_patch_index+1)%4
+    prev_index = (origin_index-1+4)%4
+    next_index = (origin_index+1)%4
+    origin_center = centers[origin_index]
+    prev_center = centers[prev_index]
+    next_center = centers[next_index]
+    vec_across = tuple([(next_center[i]-origin_center[i])/5.0 for i in [0,1]])
+    vec_down = tuple([(prev_center[i]-origin_center[i])/3.0 for i in [0,1]])
+
+    # Sanity check: test that the R,G,B,black,white patches are correct.
+    patch_info = [("r",2,2), ("g",2,1), ("b",2,0), ("w",3,0), ("k",3,5)]
+    for i in range(len(patch_info)):
+        color,yi,xi = patch_info[i]
+        xc = int(origin_center[0] + vec_across[0]*xi + vec_down[0]*yi)
+        yc = int(origin_center[1] + vec_across[1]*xi + vec_down[1]*yi)
+        means = __measure_color_checker_patch(img, xc,yc, 64)
+        if color == "r":
+            high_chans = [0]
+            low_chans = [1,2]
+        elif color == "g":
+            high_chans = [1]
+            low_chans = [0,2]
+        elif color == "b":
+            high_chans = [2]
+            low_chans = [0,1]
+        elif color == "w":
+            high_chans = [0,1,2]
+            low_chans = []
+        elif color == "k":
+            high_chans = []
+            low_chans = [0,1,2]
+        else:
+            assert(False)
+
+        # If the debug info is requested, then don't assert that the patches
+        # are matched, to allow the caller to see the output.
+        if debug_fname_prefix is not None:
+            img[int(yc),int(xc)] = 1.0
+        else:
+            assert(min([means[i] for i in high_chans]+[1]) > \
+                   max([means[i] for i in low_chans]+[0]))
+
+    if debug_fname_prefix is not None:
+        write_image(img, debug_fname_prefix+"_2.jpg")
+
+    return origin_center, vec_across, vec_down
+
 class __UnitTest(unittest.TestCase):
     """Run a suite of unit tests on this module.
     """
diff --git a/apps/CameraITS/tools/compute_dng_noise_model.py b/apps/CameraITS/tools/compute_dng_noise_model.py
new file mode 100644
index 0000000..e528332
--- /dev/null
+++ b/apps/CameraITS/tools/compute_dng_noise_model.py
@@ -0,0 +1,176 @@
+# Copyright 2014 The Android Open Source Project
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import its.device
+import its.objects
+import its.image
+import pprint
+import pylab
+import os.path
+import matplotlib
+import matplotlib.pyplot
+import numpy
+import math
+
+def main():
+    """Compute the DNG noise model from a color checker chart.
+
+    TODO: Make this more robust; some manual futzing may be needed to get
+    this to work on a new device.
+    """
+    NAME = os.path.basename(__file__).split(".")[0]
+
+    with its.device.ItsSession() as cam:
+
+        props = cam.get_camera_properties()
+
+        white_level = float(props['android.sensor.info.whiteLevel'])
+        black_levels = props['android.sensor.blackLevelPattern']
+        idxs = its.image.get_canonical_cfa_order(props)
+        black_levels = [black_levels[i] for i in idxs]
+
+        # Expose for the scene with min sensitivity
+        sens_min, sens_max = props['android.sensor.info.sensitivityRange']
+        s_ae,e_ae,awb_gains,awb_ccm,_  = cam.do_3a()
+        s_e_prod = s_ae * e_ae
+
+        # Make the image brighter since the script looks at linear Bayer
+        # raw patches rather than gamma-encoded YUV patches (and the AE
+        # probably under-exposes a little for this use-case).
+        s_e_prod *= 2
+
+        # Capture raw frames across the full sensitivity range.
+        SENSITIVITY_STEP = 400
+        reqs = []
+        sens = []
+        for s in range(sens_min, sens_max, SENSITIVITY_STEP):
+            e = int(s_e_prod / float(s))
+            req = its.objects.manual_capture_request(s, e)
+            req["android.colorCorrection.transform"] = \
+                    its.objects.float_to_rational(awb_ccm)
+            req["android.colorCorrection.gains"] = awb_gains
+            reqs.append(req)
+            sens.append(s)
+
+        caps = cam.do_capture(reqs, cam.CAP_RAW)
+
+        # A list of the (x,y) coords of the top-left pixel of a collection of
+        # 64x64 pixel patches of a color checker chart. Each patch should be
+        # uniform, however the actual color doesn't matter.
+        img = its.image.convert_capture_to_rgb_image(caps[0], props=props)
+        (x0,y0),(dxh,dyh),(dxv,dyv) = \
+                its.image.get_color_checker_chart_patches(img, NAME+"_debug")
+        patches = []
+        for xi in range(6):
+            for yi in range(4):
+                xc = int(x0 + dxh*xi + dxv*yi)
+                yc = int(y0 + dyh*xi + dyv*yi)
+                patches.append((xc-32,yc-32))
+
+        lines = []
+        for (s,cap) in zip(sens,caps):
+            # For each capture, compute the mean value in each patch, for each
+            # Bayer plane; discard patches where pixels are close to clamped.
+            # Also compute the variance.
+            CLAMP_THRESH = 0.2
+            planes = its.image.convert_capture_to_planes(cap, props)
+            points = []
+            for i,plane in enumerate(planes):
+                plane = (plane * white_level - black_levels[i]) / (
+                        white_level - black_levels[i])
+                for (x,y) in patches:
+                    tile = plane[y/2:y/2+32,x/2:x/2+32,:]
+                    mean = its.image.compute_image_means(tile)[0]
+                    var = its.image.compute_image_variances(tile)[0]
+                    if (mean > CLAMP_THRESH and mean < 1.0-CLAMP_THRESH):
+                        # Each point is a (mean,variance) tuple for a patch;
+                        # for a given ISO, there should be a linear
+                        # relationship between these values.
+                        points.append((mean,var))
+
+            # Fit a line to the points, with a line equation: y = mx + b.
+            # This line is the relationship between mean and variance (i.e.)
+            # between signal level and noise, for this particular sensor.
+            # In the DNG noise model, the gradient (m) is "S", and the offset
+            # (b) is "O".
+            points.sort()
+            xs = [x for (x,y) in points]
+            ys = [y for (x,y) in points]
+            m,b = numpy.polyfit(xs, ys, 1)
+            lines.append((s,m,b))
+            print s, "->", m, b
+
+            # Some sanity checks:
+            # * Noise levels should increase with brightness.
+            # * Extrapolating to a black image, the noise should be positive.
+            # Basically, the "b" value should correspnd to the read noise,
+            # which is the noise level if the sensor was operating in zero
+            # light.
+            #assert(m > 0)
+            #assert(b >= 0)
+
+            # Draw a plot.
+            pylab.plot(xs, ys, 'r')
+            pylab.plot([0,xs[-1]],[b,m*xs[-1]+b],'b')
+            matplotlib.pyplot.savefig("%s_plot_mean_vs_variance.png" % (NAME))
+
+        # Now fit a line across the (m,b) line parameters for each sensitivity.
+        # The gradient (m) params are fit to the "S" line, and the offset (b)
+        # params are fit to the "O" line, both as a function of sensitivity.
+        gains = [d[0] for d in lines]
+        Ss = [d[1] for d in lines]
+        Os = [d[2] for d in lines]
+        mS,bS = numpy.polyfit(gains, Ss, 1)
+        mO,bO = numpy.polyfit(gains, Os, 1)
+
+        # Plot curve "O" as 10x, so it fits in the same scale as curve "S".
+        pylab.plot(gains, [10*o for o in Os], 'r')
+        pylab.plot([gains[0],gains[-1]],
+                [10*mO*gains[0]+10*bO, 10*mO*gains[-1]+10*bO], 'b')
+        pylab.plot(gains, Ss, 'r')
+        pylab.plot([gains[0],gains[-1]], [mS*gains[0]+bS, mS*gains[-1]+bS], 'b')
+        matplotlib.pyplot.savefig("%s_plot_S_O.png" % (NAME))
+
+        print """
+        /* Generated test code to dump a table of data for external validation
+         * of the noise model parameters.
+         */
+        #include <stdio.h>
+        #include <assert.h>
+        void compute_noise_model_entries(int sens, double *o, double *s);
+        int main(void) {
+            int sens;
+            double o, s;
+            for (sens = %d; sens <= %d; sens += 100) {
+                compute_noise_model_entries(sens, &o, &s);
+                printf("%%d,%%lf,%%lf\\n", sens, o, s);
+            }
+            return 0;
+        }
+
+        /* Generated function to map a given sensitivity to the O and S noise
+         * model parameters in the DNG noise model.
+         */
+        void compute_noise_model_entries(int sens, double *o, double *s) {
+            assert(sens >= %d && sens <= %d && o && s);
+            *s = %e * sens + %e;
+            *o = %e * sens + %e;
+            *s = *s < 0.0 ? 0.0 : *s;
+            *o = *o < 0.0 ? 0.0 : *o;
+        }
+        """%(sens_min,sens_max,sens_min,sens_max,mS,bS,mO,bO)
+
+if __name__ == '__main__':
+    main()
+