[libcamera-devel,1/3] utils: raspberrypi: ctt: Improved color matrix fitting

18800 diff mbox series

Added code which optimises the color matrices based off
delta E values for the calibration images. Working in LAB
color space.

Signed-off-by Ben Benson <ben.benson@raspberrypi.com>
---
 utils/raspberrypi/ctt/colors.py        |  30 +++
 utils/raspberrypi/ctt/ctt_ccm.py       | 258 +++++++++++++++++++++----
 utils/raspberrypi/ctt/ctt_visualise.py |  38 ++++
 3 files changed, 291 insertions(+), 35 deletions(-)
 create mode 100644 utils/raspberrypi/ctt/colors.py
 create mode 100644 utils/raspberrypi/ctt/ctt_visualise.py

Hi Ben

Thanks very much for this work!

I'm going to make a few more minor or style related comments, but I'm
happy to regard most of them as being for future reference. We could
go on tweaking the Python forever, but with just a couple of
exceptions, I'd rather not!

On Fri, 7 Jul 2023 at 14:41, Ben Benson via libcamera-devel
<libcamera-devel@lists.libcamera.org> wrote:
>
> Added code which optimises the color matrices based off
> delta E values for the calibration images. Working in LAB
> color space.
>
> Signed-off-by Ben Benson <ben.benson@raspberrypi.com>

s/by/by:/
(We probably need to fix this one.)

As I've said, minor observations notwithstanding:

Reviewed-by: David Plowman <david.plowman@raspberrypi.com>

> ---
>  utils/raspberrypi/ctt/colors.py        |  30 +++
>  utils/raspberrypi/ctt/ctt_ccm.py       | 258 +++++++++++++++++++++----
>  utils/raspberrypi/ctt/ctt_visualise.py |  38 ++++
>  3 files changed, 291 insertions(+), 35 deletions(-)
>  create mode 100644 utils/raspberrypi/ctt/colors.py
>  create mode 100644 utils/raspberrypi/ctt/ctt_visualise.py
>
> diff --git a/utils/raspberrypi/ctt/colors.py b/utils/raspberrypi/ctt/colors.py
> new file mode 100644
> index 00000000..bcd87e35
> --- /dev/null
> +++ b/utils/raspberrypi/ctt/colors.py
> @@ -0,0 +1,30 @@
> +# SIMPLE CODE TO CONVERT RGB VALUES TO LAB

No need to shout!  :)

> +def RGB_to_LAB(RGB):  # where RGB is a 1x3 array.   e.g RGB = [100, 255, 230]
> +    num = 0
> +    XYZ = [0, 0, 0]
> +    # converted all the three R, G, B to X, Y, Z
> +    X = RGB[0] * 0.4124 + RGB[1] * 0.3576 + RGB[2] * 0.1805
> +    Y = RGB[0] * 0.2126 + RGB[1] * 0.7152 + RGB[2] * 0.0722
> +    Z = RGB[0] * 0.0193 + RGB[1] * 0.1192 + RGB[2] * 0.9505
> +
> +    XYZ[0] = X / 255 * 100
> +    XYZ[1] = Y / 255 * 100  # XYZ Must be in range 0 -> 100, so scale down from 255
> +    XYZ[2] = Z / 255 * 100
> +    XYZ[0] = XYZ[0] / 95.047  # ref_X =  95.047   Observer= 2°, Illuminant= D65
> +    XYZ[1] = XYZ[1] / 100.0  # ref_Y = 100.000
> +    XYZ[2] = XYZ[2] / 108.883  # ref_Z = 108.883
> +    num = 0
> +    for value in XYZ:
> +        if value > 0.008856:
> +            value = value ** (0.3333333333333333)
> +        else:
> +            value = (7.787 * value) + (16 / 116)
> +        XYZ[num] = value
> +        num = num + 1
> +
> +    # L, A, B, values calculated below
> +    L = (116 * XYZ[1]) - 16
> +    a = 500 * (XYZ[0] - XYZ[1])
> +    b = 200 * (XYZ[1] - XYZ[2])
> +
> +    return [L, a, b]
> diff --git a/utils/raspberrypi/ctt/ctt_ccm.py b/utils/raspberrypi/ctt/ctt_ccm.py
> index 376cc712..bd44b4d8 100644
> --- a/utils/raspberrypi/ctt/ctt_ccm.py
> +++ b/utils/raspberrypi/ctt/ctt_ccm.py
> @@ -6,27 +6,66 @@
>
>  from ctt_image_load import *
>  from ctt_awb import get_alsc_patches
> -
> -
> +import colors
> +from scipy.optimize import minimize
> +from ctt_visualise import visualise_macbeth_chart
> +import numpy as np
>  """
>  takes 8-bit macbeth chart values, degammas and returns 16 bit
>  """
> +
> +'''
> +This program has many options from which to derive the color matrix from.
> +The first is average. This minimises the average delta E across all patches of
> +the macbeth chart. Testing across all cameras yeilded this as the most color
> +accurate and vivid. Other options are avalible however.
> +Maximum minimises the maximum Delta E of the patches. It iterates through till
> +a minimum maximum is found (so that there is
> +not one patch that deviates wildly.)
> +This yeilds generally good results but overall the colors are less accurate

s/yeilds/yields/   (but don't mind)

> +Have a fiddle with maximum and see what you think.
> +The final option allows you to select the patches for which to average across.
> +This means that you can bias certain patches, for instance if you want the
> +reds to be more accurate.
> +'''
> +
> +matrix_selection_types = ["average", "maximum", "patches"]
> +typenum = 0  # select from array above, 0 = average, 1 = maximum, 2 = patches
> +test_patches = [1, 2, 5, 8, 9, 12, 14]
> +
> +'''
> +Enter patches to test for. Can also be entered twice if you
> +would like twice as much bias on one patch.
> +'''
> +
> +
>  def degamma(x):
> -    x = x / ((2**8)-1)
> -    x = np.where(x < 0.04045, x/12.92, ((x+0.055)/1.055)**2.4)
> -    x = x * ((2**16)-1)
> +    x = x / ((2 ** 8) - 1)  # takes 255 and scales it down to one
> +    x = np.where(x < 0.04045, x / 12.92, ((x + 0.055) / 1.055) ** 2.4)
> +    x = x * ((2 ** 16) - 1)  # takes one and scales up to 255

I guess the comment here is not correct, it should be 65535 and not
255. (This might be another change to make as it's a bit confusing for
the reader otherwise.)

>      return x
>
>
> +def gamma(x):
> +    # return (x *  * (1 / 2.4)  *  1.055 - 0.055)
> +    e = []
> +    for i in range(len(x)):
> +        e.append(((x[i] / 255) ** (1 / 2.4) * 1.055 - 0.055) * 255)
> +    return e
> +
> +
>  """
>  FInds colour correction matrices for list of images
>  """
> +
> +
>  def ccm(Cam, cal_cr_list, cal_cb_list):
> +    global matrix_selection_types, typenum
>      imgs = Cam.imgs
>      """
>      standard macbeth chart colour values
>      """
> -    m_rgb = np.array([  # these are in sRGB
> +    m_rgb = np.array([  # these are in RGB
>          [116, 81, 67],    # dark skin
>          [199, 147, 129],  # light skin
>          [91, 122, 156],   # blue sky
> @@ -34,7 +73,7 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          [130, 128, 176],  # blue flower
>          [92, 190, 172],   # bluish green
>          [224, 124, 47],   # orange
> -        [68, 91, 170],     # purplish blue
> +        [68, 91, 170],    # purplish blue
>          [198, 82, 97],    # moderate red
>          [94, 58, 106],    # purple
>          [159, 189, 63],   # yellow green
> @@ -52,16 +91,22 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          [82, 84, 86],     # neutral 3.5
>          [49, 49, 51]      # black 2
>      ])
> -
>      """
>      convert reference colours from srgb to rgb
>      """
> -    m_srgb = degamma(m_rgb)
> +    m_srgb = degamma(m_rgb)  # now in 16 bit color.
> +
> +    m_lab = []
> +    for col in m_srgb:
> +        m_lab.append(colors.RGB_to_LAB(col / 256))
> +    # This produces matrix of LAB values for ideal color chart)

I guess we could have looked the Lab values up explicitly, but this is OK too.

> +
>      """
>      reorder reference values to match how patches are ordered
>      """
>      m_srgb = np.array([m_srgb[i::6] for i in range(6)]).reshape((24, 3))
> -
> +    m_lab = np.array([m_lab[i::6] for i in range(6)]).reshape((24, 3))
> +    m_rgb = np.array([m_rgb[i::6] for i in range(6)]).reshape((24, 3))
>      """
>      reformat alsc correction tables or set colour_cals to None if alsc is
>      deactivated
> @@ -76,8 +121,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>              """
>              normalise tables so min value is 1
>              """
> -            cr_tab = cr_tab/np.min(cr_tab)
> -            cb_tab = cb_tab/np.min(cb_tab)
> +            cr_tab = cr_tab / np.min(cr_tab)
> +            cb_tab = cb_tab / np.min(cb_tab)
>              colour_cals[cr['ct']] = [cr_tab, cb_tab]
>
>      """
> @@ -94,6 +139,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          the function will simply return the macbeth patches
>          """
>          r, b, g = get_alsc_patches(Img, colour_cals, grey=False)
> +        # 256 values for each patch of sRGB values
> +
>          """
>          do awb
>          Note: awb is done by measuring the macbeth chart in the image, rather
> @@ -101,34 +148,123 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          and the ccm matrices will be more accurate.
>          """
>          r_greys, b_greys, g_greys = r[3::4], b[3::4], g[3::4]
> -        r_g = np.mean(r_greys/g_greys)
> -        b_g = np.mean(b_greys/g_greys)
> +        r_g = np.mean(r_greys / g_greys)
> +        b_g = np.mean(b_greys / g_greys)
>          r = r / r_g
>          b = b / b_g
> -
>          """
>          normalise brightness wrt reference macbeth colours and then average
>          each channel for each patch
>          """
> -        gain = np.mean(m_srgb)/np.mean((r, g, b))
> +        gain = np.mean(m_srgb) / np.mean((r, g, b))
>          Cam.log += '\nGain with respect to standard colours: {:.3f}'.format(gain)
> -        r = np.mean(gain*r, axis=1)
> -        b = np.mean(gain*b, axis=1)
> -        g = np.mean(gain*g, axis=1)
> -
> +        r = np.mean(gain * r, axis=1)
> +        b = np.mean(gain * b, axis=1)
> +        g = np.mean(gain * g, axis=1)
>          """
>          calculate ccm matrix
>          """
> +        # ==== All of below should in sRGB ===##
> +        sumde = 0
>          ccm = do_ccm(r, g, b, m_srgb)
> +        # This is the initial guess that our optimisation code works with.
> +
> +        r1 = ccm[0]
> +        r2 = ccm[1]
> +        g1 = ccm[3]
> +        g2 = ccm[4]
> +        b1 = ccm[6]
> +        b2 = ccm[7]
> +        '''
> +        COLOR MATRIX LOOKS AS BELOW
> +        R1 R2 R3   Rval     Outr
> +        G1 G2 G3  *  Gval  =  G
> +        B1 B2 B3   Bval     B
> +        Will be optimising 6 elements and working out the third element using 1-r1-r2 = r3
> +        '''
> +
> +        x0 = [r1, r2, g1, g2, b1, b2]
> +        '''
> +        We use our old CCM as the initial guess for the program to find the
> +        optimised matrix
> +        '''
> +        result = minimize(guess, x0, args=(r, g, b, m_lab), tol=0.0000000001)
> +        '''
> +        This produces a color matrix which has the lowest delta E possible,
> +        based off the input data. Note it is impossible for this to reach
> +        zero since inperfect data.

s/inperfect/imperfect/       (though I'm not overly fussed)
or even "the input data will be imperfect".

> +        '''
> +
> +        Cam.log += ("\n \n Optimised Matrix Below: \n \n")
> +        [r1, r2, g1, g2, b1, b2] = result.x
> +        # The new, optimised color correction matrix values
> +        optimised_ccm = [r1, r2, (1 - r1 - r2), g1, g2, (1 - g1 - g2), b1, b2, (1 - b1 - b2)]
> +        # This is the optimised Color Matrix (preserving greys by summing rows up to 1)
> +        Cam.log += str(optimised_ccm)
> +        Cam.log += "\n Old Color Correction Matrix Below \n"
> +        Cam.log += str(ccm)
> +
> +        formatted_ccm = np.array(ccm).reshape((3, 3))
> +
> +        '''
> +        below is a whole load of code that then applies the latest color
> +        matrix, and returns LAB values for color. This can then be used
> +        to calculate the final delta E
> +        '''
> +        optimised_ccm_rgb = []  # Original Color Corrected Matrix RGB / LAB
> +        optimised_ccm_lab = []
> +        for w in range(24):
> +            RGB = np.array([r[w], g[w], b[w]])
> +            ccm_applied_rgb = np.dot(formatted_ccm, (RGB / 256))
> +            optimised_ccm_rgb.append(gamma(ccm_applied_rgb))
> +            optimised_ccm_lab.append(colors.RGB_to_LAB(ccm_applied_rgb))
> +
> +        formatted_optimised_ccm = np.array(ccm).reshape((3, 3))
> +        after_gamma_rgb = []
> +        after_gamma_lab = []
> +        for w in range(24):
> +            RGB = np.array([r[w], g[w], b[w]])
> +            optimised_ccm_applied_rgb = np.dot(formatted_optimised_ccm, RGB / 256)
> +            after_gamma_rgb.append(gamma(optimised_ccm_applied_rgb))
> +            after_gamma_lab.append(colors.RGB_to_LAB(optimised_ccm_applied_rgb))
> +        '''
> +        Gamma After RGB / LAB
> +        We now want to spit out some data that shows
> +        how the optimisation has improved the color matricies

s/matricies/matrices/   but again, I don't really mind

> +        '''
> +        Cam.log += "Here are the Improvements"
> +
> +        # CALCULATE WORST CASE delta e
> +        old_worst_delta_e = 0
> +        bavg = transform_and_evaluate(formatted_ccm, r, g, b, m_lab)
> +        new_worst_delta_e = 0
> +        aavg = transform_and_evaluate(formatted_optimised_ccm, r, g, b, m_lab)

One might consider renaming "bavg" as "before_avg" and "aavg" as
"after_avg" to make this a bit more readable? But not urgent.

> +        for i in range(24):
> +            old_delta_e = deltae(optimised_ccm_lab[i], m_lab[i])  # Current Old Delta E
> +            new_delta_e = deltae(after_gamma_lab[i], m_lab[i])  # Current New Delta E
> +            if old_delta_e > old_worst_delta_e:
> +                old_worst_delta_e = old_delta_e
> +            if new_delta_e > new_worst_delta_e:
> +                new_worst_delta_e = new_delta_e

There's probably some 1 or 2-line Python idiom for the above... though
I don't particularly feel there's a need to go and find it!

> +
> +        Cam.log += "Before color correction matrix was optimised, we got an average delta E of " + str(bavg) + " and a maximum delta E of " + str(old_worst_delta_e)
> +        Cam.log += "After color correction matrix was optimised, we got an average delta E of " + str(aavg) + " and a maximum delta E of " + str(new_worst_delta_e)
> +
> +        visualise_macbeth_chart(m_rgb, optimised_ccm_rgb, after_gamma_rgb, str(Img.col) + str(matrix_selection_types[typenum]))
> +        '''
> +        The program will also save some visualisations of improvements.
> +        Very pretty to look at. Top rectangle is ideal, Left square is
> +        before optimisation, right square is after.
> +        '''
>
>          """
>          if a ccm has already been calculated for that temperature then don't
>          overwrite but save both. They will then be averaged later on
> -        """
> +        """  # NOW GOING TO USE OPTIMISED COLOR MATRIX, optimised_ccm
>          if Img.col in ccm_tab.keys():
> -            ccm_tab[Img.col].append(ccm)
> +            ccm_tab[Img.col].append(optimised_ccm)
>          else:
> -            ccm_tab[Img.col] = [ccm]
> +            ccm_tab[Img.col] = [optimised_ccm]
>          Cam.log += '\n'
>
>      Cam.log += '\nFinished processing images'
> @@ -137,8 +273,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>      """
>      for k, v in ccm_tab.items():
>          tab = np.mean(v, axis=0)
> -        tab = np.where((10000*tab) % 1 <= 0.05, tab+0.00001, tab)
> -        tab = np.where((10000*tab) % 1 >= 0.95, tab-0.00001, tab)
> +        tab = np.where((10000 * tab) % 1 <= 0.05, tab + 0.00001, tab)
> +        tab = np.where((10000 * tab) % 1 >= 0.95, tab - 0.00001, tab)
>          ccm_tab[k] = list(np.round(tab, 5))
>          Cam.log += '\nMatrix calculated for colour temperature of {} K'.format(k)
>
> @@ -156,6 +292,50 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>      return ccms
>
>
> +def guess(x0, r, g, b, m_lab):       # provides a method of numerical feedback for the optimisation code
> +    [r1, r2, g1, g2, b1, b2] = x0
> +    ccm = np.array([r1, r2, (1 - r1 - r2),
> +                    g1, g2, (1 - g1 - g2),
> +                    b1, b2, (1 - b1 - b2)]).reshape((3, 3))  # format the matrix correctly
> +    return transform_and_evaluate(ccm, r, g, b, m_lab)
> +
> +
> +def transform_and_evaluate(ccm, r, g, b, m_lab):  # Transforms colors to LAB and applies the correction matrix
> +    # create list of matrix changed colors
> +    realrgb = []
> +    for i in range(len(r)):
> +        RGB = np.array([r[i], g[i], b[i]])
> +        rgb_post_ccm = np.dot(ccm, RGB)  # This is RGB values after the color correction matrix has been applied
> +        realrgb.append(colors.RGB_to_LAB(rgb_post_ccm))
> +    # now compare that with m_lab and return numeric result, averaged for each patch
> +    return (sumde(realrgb, m_lab) / 24)  # returns an average result of delta E
> +
> +
> +def sumde(listA, listB):
> +    global typenum, test_patches
> +    sumde = 0
> +    maxde = 0
> +    patchde = []
> +    for i in range(len(listA)):
> +        if maxde < (deltae(listA[i], listB[i])):
> +            maxde = deltae(listA[i], listB[i])
> +        patchde.append(deltae(listA[i], listB[i]))
> +        sumde += deltae(listA[i], listB[i])
> +    '''
> +    The different options specified at the start allow for
> +    the maximum to be returned, average or specific patches
> +    '''
> +    if typenum == 0:
> +        return sumde
> +    if typenum == 1:
> +        return maxde
> +    if typenum == 2:
> +        output = 0
> +        for y in range(len(test_patches)):
> +            output += patchde[test_patches[y]]   # grabs the specific patches (no need for averaging here)
> +        return output
> +
> +
>  """
>  calculates the ccm for an individual image.
>  ccms are calculate in rgb space, and are fit by hand. Although it is a 3x3

s/calculate/calculated/

> @@ -164,12 +344,14 @@ calculation.
>  Should you want to fit them in another space (e.g. LAB) we wish you the best of
>  luck and send us the code when you are done! :-)
>  """
> +
> +
>  def do_ccm(r, g, b, m_srgb):
>      rb = r-b
>      gb = g-b
> -    rb_2s = (rb*rb)
> -    rb_gbs = (rb*gb)
> -    gb_2s = (gb*gb)
> +    rb_2s = (rb * rb)
> +    rb_gbs = (rb * gb)
> +    gb_2s = (gb * gb)
>
>      r_rbs = rb * (m_srgb[..., 0] - b)
>      r_gbs = gb * (m_srgb[..., 0] - b)
> @@ -191,7 +373,7 @@ def do_ccm(r, g, b, m_srgb):
>      b_rb = np.sum(b_rbs)
>      b_gb = np.sum(b_gbs)
>
> -    det = rb_2*gb_2 - rb_gb*rb_gb
> +    det = rb_2 * gb_2 - rb_gb * rb_gb
>
>      """
>      Raise error if matrix is singular...
> @@ -201,19 +383,19 @@ def do_ccm(r, g, b, m_srgb):
>      if det < 0.001:
>          raise ArithmeticError
>
> -    r_a = (gb_2*r_rb - rb_gb*r_gb)/det
> -    r_b = (rb_2*r_gb - rb_gb*r_rb)/det
> +    r_a = (gb_2 * r_rb - rb_gb * r_gb) / det
> +    r_b = (rb_2 * r_gb - rb_gb * r_rb) / det
>      """
>      Last row can be calculated by knowing the sum must be 1
>      """
>      r_c = 1 - r_a - r_b
>
> -    g_a = (gb_2*g_rb - rb_gb*g_gb)/det
> -    g_b = (rb_2*g_gb - rb_gb*g_rb)/det
> +    g_a = (gb_2 * g_rb - rb_gb * g_gb) / det
> +    g_b = (rb_2 * g_gb - rb_gb * g_rb) / det
>      g_c = 1 - g_a - g_b
>
> -    b_a = (gb_2*b_rb - rb_gb*b_gb)/det
> -    b_b = (rb_2*b_gb - rb_gb*b_rb)/det
> +    b_a = (gb_2 * b_rb - rb_gb * b_gb) / det
> +    b_b = (rb_2 * b_gb - rb_gb * b_rb) / det
>      b_c = 1 - b_a - b_b
>
>      """
> @@ -222,3 +404,9 @@ def do_ccm(r, g, b, m_srgb):
>      ccm = [r_a, r_b, r_c, g_a, g_b, g_c, b_a, b_b, b_c]
>
>      return ccm
> +
> +
> +def deltae(colorA, colorB):
> +    return ((colorA[0] - colorB[0]) ** 2 + (colorA[1] - colorB[1]) ** 2 + (colorA[2] - colorB[2]) ** 2) ** 0.5
> +    # return ((colorA[1]-colorB[1]) *  * 2 + (colorA[2]-colorB[2]) *  * 2) *  * 0.5
> +    # UNCOMMENT IF YOU WANT TO NEGLECT LUMINANCE FROM CALCULATION OF DELTA E
> diff --git a/utils/raspberrypi/ctt/ctt_visualise.py b/utils/raspberrypi/ctt/ctt_visualise.py
> new file mode 100644
> index 00000000..fe5381c6
> --- /dev/null
> +++ b/utils/raspberrypi/ctt/ctt_visualise.py
> @@ -0,0 +1,38 @@
> +# Some code that will save virtual macbeth charts that show the difference between optimised matricies and non optimised matricies

As above ("matrices")

> +import numpy as np
> +from PIL import Image
> +'''make patch 200 by 200 pixels:
> +200x100 will be correct color
> +100x100 will be new color, 100x100 old color
> +if 50px gap between pixels, image will be 7 gaps + 6 patches wide, and 5 gaps + 4 patches high.
> +  y --->
> +x
> +|
> +|
> +'''
> +
> +
> +def visualise_macbeth_chart(macbeth_rgb, original_rgb, new_rgb, output_filename):

Maybe worth a docstring to say what it's doing?

> +    image = np.zeros((1050, 1550, 3), dtype=np.uint8)
> +    colorindex = -1
> +    for y in range(6):
> +        for x in range(4):  # Creates 6 x 4 grid of macbeth chart
> +            colorindex += 1
> +            xlocation = 50 + 250 * x  # Means there is 50px of black gap between each square, more like the real macbeth chart.
> +            ylocation = 50 + 250 * y
> +            for g in range(200):
> +                for i in range(100):
> +                    image[xlocation + i, ylocation + g] = macbeth_rgb[colorindex]
> +            xlocation = 150 + 250 * x
> +            ylocation = 50 + 250 * y
> +            for i in range(100):
> +                for g in range(100):
> +                    image[xlocation + i, ylocation + g] = original_rgb[colorindex]  # Smaller squares below to compare the old colors with the new ones
> +            xlocation = 150 + 250 * x
> +            ylocation = 150 + 250 * y
> +            for i in range(100):
> +                for g in range(100):
> +                    image[xlocation + i, ylocation + g] = new_rgb[colorindex]
> +
> +    img = Image.fromarray(image, 'RGB')
> +    img.save(str(output_filename) + 'Generated Macbeth Chart.png')
> --
> 2.34.1
>

Thanks very much!

David

Hi Ben,

Thank you for this work!

On Fri, 7 Jul 2023 at 14:41, Ben Benson via libcamera-devel
<libcamera-devel@lists.libcamera.org> wrote:
>
> Added code which optimises the color matrices based off
> delta E values for the calibration images. Working in LAB
> color space.
>
> Signed-off-by Ben Benson <ben.benson@raspberrypi.com>

With David's minor style suggestions applied:

Reviewed-by: Naushir Patuck <naush@raspberrypi.com>

> ---
>  utils/raspberrypi/ctt/colors.py        |  30 +++
>  utils/raspberrypi/ctt/ctt_ccm.py       | 258 +++++++++++++++++++++----
>  utils/raspberrypi/ctt/ctt_visualise.py |  38 ++++
>  3 files changed, 291 insertions(+), 35 deletions(-)
>  create mode 100644 utils/raspberrypi/ctt/colors.py
>  create mode 100644 utils/raspberrypi/ctt/ctt_visualise.py
>
> diff --git a/utils/raspberrypi/ctt/colors.py b/utils/raspberrypi/ctt/colors.py
> new file mode 100644
> index 00000000..bcd87e35
> --- /dev/null
> +++ b/utils/raspberrypi/ctt/colors.py
> @@ -0,0 +1,30 @@
> +# SIMPLE CODE TO CONVERT RGB VALUES TO LAB
> +def RGB_to_LAB(RGB):  # where RGB is a 1x3 array.   e.g RGB = [100, 255, 230]
> +    num = 0
> +    XYZ = [0, 0, 0]
> +    # converted all the three R, G, B to X, Y, Z
> +    X = RGB[0] * 0.4124 + RGB[1] * 0.3576 + RGB[2] * 0.1805
> +    Y = RGB[0] * 0.2126 + RGB[1] * 0.7152 + RGB[2] * 0.0722
> +    Z = RGB[0] * 0.0193 + RGB[1] * 0.1192 + RGB[2] * 0.9505
> +
> +    XYZ[0] = X / 255 * 100
> +    XYZ[1] = Y / 255 * 100  # XYZ Must be in range 0 -> 100, so scale down from 255
> +    XYZ[2] = Z / 255 * 100
> +    XYZ[0] = XYZ[0] / 95.047  # ref_X =  95.047   Observer= 2°, Illuminant= D65
> +    XYZ[1] = XYZ[1] / 100.0  # ref_Y = 100.000
> +    XYZ[2] = XYZ[2] / 108.883  # ref_Z = 108.883
> +    num = 0
> +    for value in XYZ:
> +        if value > 0.008856:
> +            value = value ** (0.3333333333333333)
> +        else:
> +            value = (7.787 * value) + (16 / 116)
> +        XYZ[num] = value
> +        num = num + 1
> +
> +    # L, A, B, values calculated below
> +    L = (116 * XYZ[1]) - 16
> +    a = 500 * (XYZ[0] - XYZ[1])
> +    b = 200 * (XYZ[1] - XYZ[2])
> +
> +    return [L, a, b]
> diff --git a/utils/raspberrypi/ctt/ctt_ccm.py b/utils/raspberrypi/ctt/ctt_ccm.py
> index 376cc712..bd44b4d8 100644
> --- a/utils/raspberrypi/ctt/ctt_ccm.py
> +++ b/utils/raspberrypi/ctt/ctt_ccm.py
> @@ -6,27 +6,66 @@
>
>  from ctt_image_load import *
>  from ctt_awb import get_alsc_patches
> -
> -
> +import colors
> +from scipy.optimize import minimize
> +from ctt_visualise import visualise_macbeth_chart
> +import numpy as np
>  """
>  takes 8-bit macbeth chart values, degammas and returns 16 bit
>  """
> +
> +'''
> +This program has many options from which to derive the color matrix from.
> +The first is average. This minimises the average delta E across all patches of
> +the macbeth chart. Testing across all cameras yeilded this as the most color
> +accurate and vivid. Other options are avalible however.
> +Maximum minimises the maximum Delta E of the patches. It iterates through till
> +a minimum maximum is found (so that there is
> +not one patch that deviates wildly.)
> +This yeilds generally good results but overall the colors are less accurate
> +Have a fiddle with maximum and see what you think.
> +The final option allows you to select the patches for which to average across.
> +This means that you can bias certain patches, for instance if you want the
> +reds to be more accurate.
> +'''
> +
> +matrix_selection_types = ["average", "maximum", "patches"]
> +typenum = 0  # select from array above, 0 = average, 1 = maximum, 2 = patches
> +test_patches = [1, 2, 5, 8, 9, 12, 14]
> +
> +'''
> +Enter patches to test for. Can also be entered twice if you
> +would like twice as much bias on one patch.
> +'''
> +
> +
>  def degamma(x):
> -    x = x / ((2**8)-1)
> -    x = np.where(x < 0.04045, x/12.92, ((x+0.055)/1.055)**2.4)
> -    x = x * ((2**16)-1)
> +    x = x / ((2 ** 8) - 1)  # takes 255 and scales it down to one
> +    x = np.where(x < 0.04045, x / 12.92, ((x + 0.055) / 1.055) ** 2.4)
> +    x = x * ((2 ** 16) - 1)  # takes one and scales up to 255
>      return x
>
>
> +def gamma(x):
> +    # return (x *  * (1 / 2.4)  *  1.055 - 0.055)
> +    e = []
> +    for i in range(len(x)):
> +        e.append(((x[i] / 255) ** (1 / 2.4) * 1.055 - 0.055) * 255)
> +    return e
> +
> +
>  """
>  FInds colour correction matrices for list of images
>  """
> +
> +
>  def ccm(Cam, cal_cr_list, cal_cb_list):
> +    global matrix_selection_types, typenum
>      imgs = Cam.imgs
>      """
>      standard macbeth chart colour values
>      """
> -    m_rgb = np.array([  # these are in sRGB
> +    m_rgb = np.array([  # these are in RGB
>          [116, 81, 67],    # dark skin
>          [199, 147, 129],  # light skin
>          [91, 122, 156],   # blue sky
> @@ -34,7 +73,7 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          [130, 128, 176],  # blue flower
>          [92, 190, 172],   # bluish green
>          [224, 124, 47],   # orange
> -        [68, 91, 170],     # purplish blue
> +        [68, 91, 170],    # purplish blue
>          [198, 82, 97],    # moderate red
>          [94, 58, 106],    # purple
>          [159, 189, 63],   # yellow green
> @@ -52,16 +91,22 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          [82, 84, 86],     # neutral 3.5
>          [49, 49, 51]      # black 2
>      ])
> -
>      """
>      convert reference colours from srgb to rgb
>      """
> -    m_srgb = degamma(m_rgb)
> +    m_srgb = degamma(m_rgb)  # now in 16 bit color.
> +
> +    m_lab = []
> +    for col in m_srgb:
> +        m_lab.append(colors.RGB_to_LAB(col / 256))
> +    # This produces matrix of LAB values for ideal color chart)
> +
>      """
>      reorder reference values to match how patches are ordered
>      """
>      m_srgb = np.array([m_srgb[i::6] for i in range(6)]).reshape((24, 3))
> -
> +    m_lab = np.array([m_lab[i::6] for i in range(6)]).reshape((24, 3))
> +    m_rgb = np.array([m_rgb[i::6] for i in range(6)]).reshape((24, 3))
>      """
>      reformat alsc correction tables or set colour_cals to None if alsc is
>      deactivated
> @@ -76,8 +121,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>              """
>              normalise tables so min value is 1
>              """
> -            cr_tab = cr_tab/np.min(cr_tab)
> -            cb_tab = cb_tab/np.min(cb_tab)
> +            cr_tab = cr_tab / np.min(cr_tab)
> +            cb_tab = cb_tab / np.min(cb_tab)
>              colour_cals[cr['ct']] = [cr_tab, cb_tab]
>
>      """
> @@ -94,6 +139,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          the function will simply return the macbeth patches
>          """
>          r, b, g = get_alsc_patches(Img, colour_cals, grey=False)
> +        # 256 values for each patch of sRGB values
> +
>          """
>          do awb
>          Note: awb is done by measuring the macbeth chart in the image, rather
> @@ -101,34 +148,123 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>          and the ccm matrices will be more accurate.
>          """
>          r_greys, b_greys, g_greys = r[3::4], b[3::4], g[3::4]
> -        r_g = np.mean(r_greys/g_greys)
> -        b_g = np.mean(b_greys/g_greys)
> +        r_g = np.mean(r_greys / g_greys)
> +        b_g = np.mean(b_greys / g_greys)
>          r = r / r_g
>          b = b / b_g
> -
>          """
>          normalise brightness wrt reference macbeth colours and then average
>          each channel for each patch
>          """
> -        gain = np.mean(m_srgb)/np.mean((r, g, b))
> +        gain = np.mean(m_srgb) / np.mean((r, g, b))
>          Cam.log += '\nGain with respect to standard colours: {:.3f}'.format(gain)
> -        r = np.mean(gain*r, axis=1)
> -        b = np.mean(gain*b, axis=1)
> -        g = np.mean(gain*g, axis=1)
> -
> +        r = np.mean(gain * r, axis=1)
> +        b = np.mean(gain * b, axis=1)
> +        g = np.mean(gain * g, axis=1)
>          """
>          calculate ccm matrix
>          """
> +        # ==== All of below should in sRGB ===##
> +        sumde = 0
>          ccm = do_ccm(r, g, b, m_srgb)
> +        # This is the initial guess that our optimisation code works with.
> +
> +        r1 = ccm[0]
> +        r2 = ccm[1]
> +        g1 = ccm[3]
> +        g2 = ccm[4]
> +        b1 = ccm[6]
> +        b2 = ccm[7]
> +        '''
> +        COLOR MATRIX LOOKS AS BELOW
> +        R1 R2 R3   Rval     Outr
> +        G1 G2 G3  *  Gval  =  G
> +        B1 B2 B3   Bval     B
> +        Will be optimising 6 elements and working out the third element using 1-r1-r2 = r3
> +        '''
> +
> +        x0 = [r1, r2, g1, g2, b1, b2]
> +        '''
> +        We use our old CCM as the initial guess for the program to find the
> +        optimised matrix
> +        '''
> +        result = minimize(guess, x0, args=(r, g, b, m_lab), tol=0.0000000001)
> +        '''
> +        This produces a color matrix which has the lowest delta E possible,
> +        based off the input data. Note it is impossible for this to reach
> +        zero since inperfect data.
> +        '''
> +
> +        Cam.log += ("\n \n Optimised Matrix Below: \n \n")
> +        [r1, r2, g1, g2, b1, b2] = result.x
> +        # The new, optimised color correction matrix values
> +        optimised_ccm = [r1, r2, (1 - r1 - r2), g1, g2, (1 - g1 - g2), b1, b2, (1 - b1 - b2)]
> +        # This is the optimised Color Matrix (preserving greys by summing rows up to 1)
> +        Cam.log += str(optimised_ccm)
> +        Cam.log += "\n Old Color Correction Matrix Below \n"
> +        Cam.log += str(ccm)
> +
> +        formatted_ccm = np.array(ccm).reshape((3, 3))
> +
> +        '''
> +        below is a whole load of code that then applies the latest color
> +        matrix, and returns LAB values for color. This can then be used
> +        to calculate the final delta E
> +        '''
> +        optimised_ccm_rgb = []  # Original Color Corrected Matrix RGB / LAB
> +        optimised_ccm_lab = []
> +        for w in range(24):
> +            RGB = np.array([r[w], g[w], b[w]])
> +            ccm_applied_rgb = np.dot(formatted_ccm, (RGB / 256))
> +            optimised_ccm_rgb.append(gamma(ccm_applied_rgb))
> +            optimised_ccm_lab.append(colors.RGB_to_LAB(ccm_applied_rgb))
> +
> +        formatted_optimised_ccm = np.array(ccm).reshape((3, 3))
> +        after_gamma_rgb = []
> +        after_gamma_lab = []
> +        for w in range(24):
> +            RGB = np.array([r[w], g[w], b[w]])
> +            optimised_ccm_applied_rgb = np.dot(formatted_optimised_ccm, RGB / 256)
> +            after_gamma_rgb.append(gamma(optimised_ccm_applied_rgb))
> +            after_gamma_lab.append(colors.RGB_to_LAB(optimised_ccm_applied_rgb))
> +        '''
> +        Gamma After RGB / LAB
> +        We now want to spit out some data that shows
> +        how the optimisation has improved the color matricies
> +        '''
> +        Cam.log += "Here are the Improvements"
> +
> +        # CALCULATE WORST CASE delta e
> +        old_worst_delta_e = 0
> +        bavg = transform_and_evaluate(formatted_ccm, r, g, b, m_lab)
> +        new_worst_delta_e = 0
> +        aavg = transform_and_evaluate(formatted_optimised_ccm, r, g, b, m_lab)
> +        for i in range(24):
> +            old_delta_e = deltae(optimised_ccm_lab[i], m_lab[i])  # Current Old Delta E
> +            new_delta_e = deltae(after_gamma_lab[i], m_lab[i])  # Current New Delta E
> +            if old_delta_e > old_worst_delta_e:
> +                old_worst_delta_e = old_delta_e
> +            if new_delta_e > new_worst_delta_e:
> +                new_worst_delta_e = new_delta_e
> +
> +        Cam.log += "Before color correction matrix was optimised, we got an average delta E of " + str(bavg) + " and a maximum delta E of " + str(old_worst_delta_e)
> +        Cam.log += "After color correction matrix was optimised, we got an average delta E of " + str(aavg) + " and a maximum delta E of " + str(new_worst_delta_e)
> +
> +        visualise_macbeth_chart(m_rgb, optimised_ccm_rgb, after_gamma_rgb, str(Img.col) + str(matrix_selection_types[typenum]))
> +        '''
> +        The program will also save some visualisations of improvements.
> +        Very pretty to look at. Top rectangle is ideal, Left square is
> +        before optimisation, right square is after.
> +        '''
>
>          """
>          if a ccm has already been calculated for that temperature then don't
>          overwrite but save both. They will then be averaged later on
> -        """
> +        """  # NOW GOING TO USE OPTIMISED COLOR MATRIX, optimised_ccm
>          if Img.col in ccm_tab.keys():
> -            ccm_tab[Img.col].append(ccm)
> +            ccm_tab[Img.col].append(optimised_ccm)
>          else:
> -            ccm_tab[Img.col] = [ccm]
> +            ccm_tab[Img.col] = [optimised_ccm]
>          Cam.log += '\n'
>
>      Cam.log += '\nFinished processing images'
> @@ -137,8 +273,8 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>      """
>      for k, v in ccm_tab.items():
>          tab = np.mean(v, axis=0)
> -        tab = np.where((10000*tab) % 1 <= 0.05, tab+0.00001, tab)
> -        tab = np.where((10000*tab) % 1 >= 0.95, tab-0.00001, tab)
> +        tab = np.where((10000 * tab) % 1 <= 0.05, tab + 0.00001, tab)
> +        tab = np.where((10000 * tab) % 1 >= 0.95, tab - 0.00001, tab)
>          ccm_tab[k] = list(np.round(tab, 5))
>          Cam.log += '\nMatrix calculated for colour temperature of {} K'.format(k)
>
> @@ -156,6 +292,50 @@ def ccm(Cam, cal_cr_list, cal_cb_list):
>      return ccms
>
>
> +def guess(x0, r, g, b, m_lab):       # provides a method of numerical feedback for the optimisation code
> +    [r1, r2, g1, g2, b1, b2] = x0
> +    ccm = np.array([r1, r2, (1 - r1 - r2),
> +                    g1, g2, (1 - g1 - g2),
> +                    b1, b2, (1 - b1 - b2)]).reshape((3, 3))  # format the matrix correctly
> +    return transform_and_evaluate(ccm, r, g, b, m_lab)
> +
> +
> +def transform_and_evaluate(ccm, r, g, b, m_lab):  # Transforms colors to LAB and applies the correction matrix
> +    # create list of matrix changed colors
> +    realrgb = []
> +    for i in range(len(r)):
> +        RGB = np.array([r[i], g[i], b[i]])
> +        rgb_post_ccm = np.dot(ccm, RGB)  # This is RGB values after the color correction matrix has been applied
> +        realrgb.append(colors.RGB_to_LAB(rgb_post_ccm))
> +    # now compare that with m_lab and return numeric result, averaged for each patch
> +    return (sumde(realrgb, m_lab) / 24)  # returns an average result of delta E
> +
> +
> +def sumde(listA, listB):
> +    global typenum, test_patches
> +    sumde = 0
> +    maxde = 0
> +    patchde = []
> +    for i in range(len(listA)):
> +        if maxde < (deltae(listA[i], listB[i])):
> +            maxde = deltae(listA[i], listB[i])
> +        patchde.append(deltae(listA[i], listB[i]))
> +        sumde += deltae(listA[i], listB[i])
> +    '''
> +    The different options specified at the start allow for
> +    the maximum to be returned, average or specific patches
> +    '''
> +    if typenum == 0:
> +        return sumde
> +    if typenum == 1:
> +        return maxde
> +    if typenum == 2:
> +        output = 0
> +        for y in range(len(test_patches)):
> +            output += patchde[test_patches[y]]   # grabs the specific patches (no need for averaging here)
> +        return output
> +
> +
>  """
>  calculates the ccm for an individual image.
>  ccms are calculate in rgb space, and are fit by hand. Although it is a 3x3
> @@ -164,12 +344,14 @@ calculation.
>  Should you want to fit them in another space (e.g. LAB) we wish you the best of
>  luck and send us the code when you are done! :-)
>  """
> +
> +
>  def do_ccm(r, g, b, m_srgb):
>      rb = r-b
>      gb = g-b
> -    rb_2s = (rb*rb)
> -    rb_gbs = (rb*gb)
> -    gb_2s = (gb*gb)
> +    rb_2s = (rb * rb)
> +    rb_gbs = (rb * gb)
> +    gb_2s = (gb * gb)
>
>      r_rbs = rb * (m_srgb[..., 0] - b)
>      r_gbs = gb * (m_srgb[..., 0] - b)
> @@ -191,7 +373,7 @@ def do_ccm(r, g, b, m_srgb):
>      b_rb = np.sum(b_rbs)
>      b_gb = np.sum(b_gbs)
>
> -    det = rb_2*gb_2 - rb_gb*rb_gb
> +    det = rb_2 * gb_2 - rb_gb * rb_gb
>
>      """
>      Raise error if matrix is singular...
> @@ -201,19 +383,19 @@ def do_ccm(r, g, b, m_srgb):
>      if det < 0.001:
>          raise ArithmeticError
>
> -    r_a = (gb_2*r_rb - rb_gb*r_gb)/det
> -    r_b = (rb_2*r_gb - rb_gb*r_rb)/det
> +    r_a = (gb_2 * r_rb - rb_gb * r_gb) / det
> +    r_b = (rb_2 * r_gb - rb_gb * r_rb) / det
>      """
>      Last row can be calculated by knowing the sum must be 1
>      """
>      r_c = 1 - r_a - r_b
>
> -    g_a = (gb_2*g_rb - rb_gb*g_gb)/det
> -    g_b = (rb_2*g_gb - rb_gb*g_rb)/det
> +    g_a = (gb_2 * g_rb - rb_gb * g_gb) / det
> +    g_b = (rb_2 * g_gb - rb_gb * g_rb) / det
>      g_c = 1 - g_a - g_b
>
> -    b_a = (gb_2*b_rb - rb_gb*b_gb)/det
> -    b_b = (rb_2*b_gb - rb_gb*b_rb)/det
> +    b_a = (gb_2 * b_rb - rb_gb * b_gb) / det
> +    b_b = (rb_2 * b_gb - rb_gb * b_rb) / det
>      b_c = 1 - b_a - b_b
>
>      """
> @@ -222,3 +404,9 @@ def do_ccm(r, g, b, m_srgb):
>      ccm = [r_a, r_b, r_c, g_a, g_b, g_c, b_a, b_b, b_c]
>
>      return ccm
> +
> +
> +def deltae(colorA, colorB):
> +    return ((colorA[0] - colorB[0]) ** 2 + (colorA[1] - colorB[1]) ** 2 + (colorA[2] - colorB[2]) ** 2) ** 0.5
> +    # return ((colorA[1]-colorB[1]) *  * 2 + (colorA[2]-colorB[2]) *  * 2) *  * 0.5
> +    # UNCOMMENT IF YOU WANT TO NEGLECT LUMINANCE FROM CALCULATION OF DELTA E
> diff --git a/utils/raspberrypi/ctt/ctt_visualise.py b/utils/raspberrypi/ctt/ctt_visualise.py
> new file mode 100644
> index 00000000..fe5381c6
> --- /dev/null
> +++ b/utils/raspberrypi/ctt/ctt_visualise.py
> @@ -0,0 +1,38 @@
> +# Some code that will save virtual macbeth charts that show the difference between optimised matricies and non optimised matricies
> +import numpy as np
> +from PIL import Image
> +'''make patch 200 by 200 pixels:
> +200x100 will be correct color
> +100x100 will be new color, 100x100 old color
> +if 50px gap between pixels, image will be 7 gaps + 6 patches wide, and 5 gaps + 4 patches high.
> +  y --->
> +x
> +|
> +|
> +'''
> +
> +
> +def visualise_macbeth_chart(macbeth_rgb, original_rgb, new_rgb, output_filename):
> +    image = np.zeros((1050, 1550, 3), dtype=np.uint8)
> +    colorindex = -1
> +    for y in range(6):
> +        for x in range(4):  # Creates 6 x 4 grid of macbeth chart
> +            colorindex += 1
> +            xlocation = 50 + 250 * x  # Means there is 50px of black gap between each square, more like the real macbeth chart.
> +            ylocation = 50 + 250 * y
> +            for g in range(200):
> +                for i in range(100):
> +                    image[xlocation + i, ylocation + g] = macbeth_rgb[colorindex]
> +            xlocation = 150 + 250 * x
> +            ylocation = 50 + 250 * y
> +            for i in range(100):
> +                for g in range(100):
> +                    image[xlocation + i, ylocation + g] = original_rgb[colorindex]  # Smaller squares below to compare the old colors with the new ones
> +            xlocation = 150 + 250 * x
> +            ylocation = 150 + 250 * y
> +            for i in range(100):
> +                for g in range(100):
> +                    image[xlocation + i, ylocation + g] = new_rgb[colorindex]
> +
> +    img = Image.fromarray(image, 'RGB')
> +    img.save(str(output_filename) + 'Generated Macbeth Chart.png')
> --
> 2.34.1
>