summaryrefslogtreecommitdiff
path: root/src/py/cam/helpers.py
blob: 2d906667eba0333fe44fc95e5f89248d56e7a271 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright (C) 2022, Tomi Valkeinen <tomi.valkeinen@ideasonboard.com>
#
# Debayering code from PiCamera documentation

from numpy.lib.stride_tricks import as_strided
import libcamera as libcam
import libcamera.utils
import numpy as np


def demosaic(data, r0, g0, g1, b0):
    # Separate the components from the Bayer data to RGB planes

    rgb = np.zeros(data.shape + (3,), dtype=data.dtype)
    rgb[r0[1]::2, r0[0]::2, 0] = data[r0[1]::2, r0[0]::2]  # Red
    rgb[g0[1]::2, g0[0]::2, 1] = data[g0[1]::2, g0[0]::2]  # Green
    rgb[g1[1]::2, g1[0]::2, 1] = data[g1[1]::2, g1[0]::2]  # Green
    rgb[b0[1]::2, b0[0]::2, 2] = data[b0[1]::2, b0[0]::2]  # Blue

    # Below we present a fairly naive de-mosaic method that simply
    # calculates the weighted average of a pixel based on the pixels
    # surrounding it. The weighting is provided by a byte representation of
    # the Bayer filter which we construct first:

    bayer = np.zeros(rgb.shape, dtype=np.uint8)
    bayer[r0[1]::2, r0[0]::2, 0] = 1  # Red
    bayer[g0[1]::2, g0[0]::2, 1] = 1  # Green
    bayer[g1[1]::2, g1[0]::2, 1] = 1  # Green
    bayer[b0[1]::2, b0[0]::2, 2] = 1  # Blue

    # Allocate an array to hold our output with the same shape as the input
    # data. After this we define the size of window that will be used to
    # calculate each weighted average (3x3). Then we pad out the rgb and
    # bayer arrays, adding blank pixels at their edges to compensate for the
    # size of the window when calculating averages for edge pixels.

    output = np.empty(rgb.shape, dtype=rgb.dtype)
    window = (3, 3)
    borders = (window[0] - 1, window[1] - 1)
    border = (borders[0] // 2, borders[1] // 2)

    rgb = np.pad(rgb, [
        (border[0], border[0]),
        (border[1], border[1]),
        (0, 0),
    ], 'constant')
    bayer = np.pad(bayer, [
        (border[0], border[0]),
        (border[1], border[1]),
        (0, 0),
    ], 'constant')

    # For each plane in the RGB data, we use a nifty numpy trick
    # (as_strided) to construct a view over the plane of 3x3 matrices. We do
    # the same for the bayer array, then use Einstein summation on each
    # (np.sum is simpler, but copies the data so it's slower), and divide
    # the results to get our weighted average:

    for plane in range(3):
        p = rgb[..., plane]
        b = bayer[..., plane]
        pview = as_strided(p, shape=(
            p.shape[0] - borders[0],
            p.shape[1] - borders[1]) + window, strides=p.strides * 2)
        bview = as_strided(b, shape=(
            b.shape[0] - borders[0],
            b.shape[1] - borders[1]) + window, strides=b.strides * 2)
        psum = np.einsum('ijkl->ij', pview)
        bsum = np.einsum('ijkl->ij', bview)
        output[..., plane] = psum // bsum

    return output


def to_rgb(fmt, size, data):
    w = size.width
    h = size.height

    if fmt == libcam.formats.YUYV:
        # YUV422
        yuyv = data.reshape((h, w // 2 * 4))

        # YUV444
        yuv = np.empty((h, w, 3), dtype=np.uint8)
        yuv[:, :, 0] = yuyv[:, 0::2]                    # Y
        yuv[:, :, 1] = yuyv[:, 1::4].repeat(2, axis=1)  # U
        yuv[:, :, 2] = yuyv[:, 3::4].repeat(2, axis=1)  # V

        m = np.array([
            [1.0, 1.0, 1.0],
            [-0.000007154783816076815, -0.3441331386566162, 1.7720025777816772],
            [1.4019975662231445, -0.7141380310058594, 0.00001542569043522235]
        ])

        rgb = np.dot(yuv, m)
        rgb[:, :, 0] -= 179.45477266423404
        rgb[:, :, 1] += 135.45870971679688
        rgb[:, :, 2] -= 226.8183044444304
        rgb = rgb.astype(np.uint8)

    elif fmt == libcam.formats.RGB888:
        rgb = data.reshape((h, w, 3))
        rgb[:, :, [0, 1, 2]] = rgb[:, :, [2, 1, 0]]

    elif fmt == libcam.formats.BGR888:
        rgb = data.reshape((h, w, 3))

    elif fmt in [libcam.formats.ARGB8888, libcam.formats.XRGB8888]:
        rgb = data.reshape((h, w, 4))
        rgb = np.flip(rgb, axis=2)
        # drop alpha component
        rgb = np.delete(rgb, np.s_[0::4], axis=2)

    elif str(fmt).startswith('S'):
        fmt = str(fmt)
        bayer_pattern = fmt[1:5]
        bitspp = int(fmt[5:])

        if bitspp == 8:
            data = data.reshape((h, w))
            data = data.astype(np.uint16)
        elif bitspp in [10, 12]:
            data = data.view(np.uint16)
            data = data.reshape((h, w))
        else:
            raise Exception('Bad bitspp:' + str(bitspp))

        idx = bayer_pattern.find('R')
        assert(idx != -1)
        r0 = (idx % 2, idx // 2)

        idx = bayer_pattern.find('G')
        assert(idx != -1)
        g0 = (idx % 2, idx // 2)

        idx = bayer_pattern.find('G', idx + 1)
        assert(idx != -1)
        g1 = (idx % 2, idx // 2)

        idx = bayer_pattern.find('B')
        assert(idx != -1)
        b0 = (idx % 2, idx // 2)

        rgb = demosaic(data, r0, g0, g1, b0)
        rgb = (rgb >> (bitspp - 8)).astype(np.uint8)

    else:
        rgb = None

    return rgb


# A naive format conversion to 24-bit RGB
def mfb_to_rgb(mfb: libcamera.utils.MappedFrameBuffer, cfg: libcam.StreamConfiguration):
    data = np.array(mfb.planes[0], dtype=np.uint8)
    rgb = to_rgb(cfg.pixel_format, cfg.size, data)
    return rgb
hl opt">:2]), axis=0) if do_alsc_colour: """ obtain 16x12 grid of intensities for each channel and subtract black level """ g = get_16x12_grid(av_ch_g, dx, dy) - Img.blacklevel_16 r = get_16x12_grid(channels[0], dx, dy) - Img.blacklevel_16 b = get_16x12_grid(channels[3], dx, dy) - Img.blacklevel_16 """ calculate ratios as 32 bit in order to be supported by medianBlur function """ cr = np.reshape(g/r, (12, 16)).astype('float32') cb = np.reshape(g/b, (12, 16)).astype('float32') cg = np.reshape(1/g, (12, 16)).astype('float32') """ median blur to remove peaks and save as float 64 """ cr = cv2.medianBlur(cr, 3).astype('float64') cb = cv2.medianBlur(cb, 3).astype('float64') cg = cv2.medianBlur(cg, 3).astype('float64') cg = cg/np.min(cg) """ debugging code showing 2D surface plot of vignetting. Quite useful for for sanity check """ if plot: hf = plt.figure(figsize=(8, 8)) ha = hf.add_subplot(311, projection='3d') """ note Y is plotted as -Y so plot has same axes as image """ X, Y = np.meshgrid(range(16), range(12)) ha.plot_surface(X, -Y, cr, cmap=cm.coolwarm, linewidth=0) ha.set_title('ALSC Plot\nImg: {}\n\ncr'.format(Img.str)) hb = hf.add_subplot(312, projection='3d') hb.plot_surface(X, -Y, cb, cmap=cm.coolwarm, linewidth=0) hb.set_title('cb') hc = hf.add_subplot(313, projection='3d') hc.plot_surface(X, -Y, cg, cmap=cm.coolwarm, linewidth=0) hc.set_title('g') # print(Img.str) plt.show() return Img.col, cr.flatten(), cb.flatten(), cg.flatten(), (w, h, dx, dy) else: """ only perform calculations for luminance shading """ g = get_16x12_grid(av_ch_g, dx, dy) - Img.blacklevel_16 cg = np.reshape(1/g, (12, 16)).astype('float32') cg = cv2.medianBlur(cg, 3).astype('float64') cg = cg/np.min(cg) if plot: hf = plt.figure(figssize=(8, 8)) ha = hf.add_subplot(1, 1, 1, projection='3d') X, Y = np.meashgrid(range(16), range(12)) ha.plot_surface(X, -Y, cg, cmap=cm.coolwarm, linewidth=0) ha.set_title('ALSC Plot (Luminance only!)\nImg: {}\n\ncg').format(Img.str) plt.show() return Img.col, None, None, cg.flatten(), (w, h, dx, dy) """ Compresses channel down to a 16x12 grid """ def get_16x12_grid(chan, dx, dy): grid = [] """ since left and bottom border will not necessarily have rectangles of dimension dx x dy, the 32nd iteration has to be handled separately. """ for i in range(11): for j in range(15): grid.append(np.mean(chan[dy*i:dy*(1+i), dx*j:dx*(1+j)])) grid.append(np.mean(chan[dy*i:dy*(1+i), 15*dx:])) for j in range(15): grid.append(np.mean(chan[11*dy:, dx*j:dx*(1+j)])) grid.append(np.mean(chan[11*dy:, 15*dx:])) """ return as np.array, ready for further manipulation """ return np.array(grid) """ obtains sigmas for red and blue, effectively a measure of the 'error' """ def get_sigma(Cam, cal_cr_list, cal_cb_list): Cam.log += '\nCalculating sigmas' """ provided colour alsc tables were generated for two different colour temperatures sigma is calculated by comparing two calibration temperatures adjacent in colour space """ """ create list of colour temperatures """ cts = [cal['ct'] for cal in cal_cr_list] # print(cts) """ calculate sigmas for each adjacent cts and return worst one """ sigma_rs = [] sigma_bs = [] for i in range(len(cts)-1): sigma_rs.append(calc_sigma(cal_cr_list[i]['table'], cal_cr_list[i+1]['table'])) sigma_bs.append(calc_sigma(cal_cb_list[i]['table'], cal_cb_list[i+1]['table'])) Cam.log += '\nColour temperature interval {} - {} K'.format(cts[i], cts[i+1]) Cam.log += '\nSigma red: {}'.format(sigma_rs[-1]) Cam.log += '\nSigma blue: {}'.format(sigma_bs[-1]) """ return maximum sigmas, not necessarily from the same colour temperature interval """ sigma_r = max(sigma_rs) if sigma_rs else 0.005 sigma_b = max(sigma_bs) if sigma_bs else 0.005 Cam.log += '\nMaximum sigmas: Red = {} Blue = {}'.format(sigma_r, sigma_b) # print(sigma_rs, sigma_bs) # print(sigma_r, sigma_b) return sigma_r, sigma_b """ calculate sigma from two adjacent gain tables """ def calc_sigma(g1, g2): """ reshape into 16x12 matrix """ g1 = np.reshape(g1, (12, 16)) g2 = np.reshape(g2, (12, 16)) """ apply gains to gain table """ gg = g1/g2 if np.mean(gg) < 1: gg = 1/gg """ for each internal patch, compute average difference between it and its 4 neighbours, then append to list """ diffs = [] for i in range(10): for j in range(14): """ note indexing is incremented by 1 since all patches on borders are not counted """ diff = np.abs(gg[i+1][j+1]-gg[i][j+1]) diff += np.abs(gg[i+1][j+1]-gg[i+2][j+1])