# SPDX-License-Identifier: BSD-2-Clause
#
# Copyright (C) 2019, Raspberry Pi (Trading) Limited
#
# ctt_alsc.py - camera tuning tool for ALSC (auto lens shading correction)
from ctt_image_load import *
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
"""
preform alsc calibration on a set of images
"""
def alsc_all(Cam, do_alsc_colour, plot):
imgs_alsc = Cam.imgs_alsc
"""
create list of colour temperatures and associated calibration tables
"""
list_col = []
list_cr = []
list_cb = []
list_cg = []
for Img in imgs_alsc:
col, cr, cb, cg, size = alsc(Cam, Img, do_alsc_colour, plot)
list_col.append(col)
list_cr.append(cr)
list_cb.append(cb)
list_cg.append(cg)
Cam.log += '\n'
Cam.log += '\nFinished processing images'
w, h, dx, dy = size
Cam.log += '\nChannel dimensions: w = {} h = {}'.format(int(w), int(h))
Cam.log += '\n16x12 grid rectangle size: w = {} h = {}'.format(dx, dy)
"""
convert to numpy array for data manipulation
"""
list_col = np.array(list_col)
list_cr = np.array(list_cr)
list_cb = np.array(list_cb)
list_cg = np.array(list_cg)
cal_cr_list = []
cal_cb_list = []
"""
only do colour calculations if required
"""
if do_alsc_colour:
Cam.log += '\nALSC colour tables'
for ct in sorted(set(list_col)):
Cam.log += '\nColour temperature: {} K'.format(ct)
"""
average tables for the same colour temperature
"""
indices = np.where(list_col == ct)
ct = int(ct)
t_r = np.mean(list_cr[indices], axis=0)
t_b = np.mean(list_cb[indices], axis=0)
"""
force numbers to be stored to 3dp.... :(
"""
t_r = np.where((100*t_r) % 1 <= 0.05, t_r+0.001, t_r)
t_b = np.where((100*t_b) % 1 <= 0.05, t_b+0.001, t_b)
t_r = np.where((100*t_r) % 1 >= 0.95, t_r-0.001, t_r)
t_b = np.where((100*t_b) % 1 >= 0.95, t_b-0.001, t_b)
t_r = np.round(t_r, 3)
t_b = np.round(t_b, 3)
r_corners = (t_r[0], t_r[15], t_r[-1], t_r[-16])
b_corners = (t_b[0], t_b[15], t_b[-1], t_b[-16])
r_cen = t_r[5*16+7]+t_r[5*16+8]+t_r[6*16+7]+t_r[6*16+8]
r_cen = round(r_cen/4, 3)
b_cen = t_b[5*16+7]+t_b[5*16+8]+t_b[6*16+7]+t_b[6*16+8]
b_cen = round(b_cen/4, 3)
Cam.log += '\nRed table corners: {}'.format(r_corners)
Cam.log += '\nRed table centre: {}'.format(r_cen)
Cam.log += '\nBlue table corners: {}'.format(b_corners)
Cam.log += '\nBlue table centre: {}'.format(b_cen)
cr_dict = {
'ct': ct,
'table': list(t_r)
}
cb_dict = {
'ct': ct,
'table': list(t_b)
}
cal_cr_list.append(cr_dict)
cal_cb_list.append(cb_dict)
Cam.log += '\n'
else:
cal_cr_list, cal_cb_list = None, None
"""
average all values for luminance shading and return one table for all temperatures
"""
lum_lut = np.mean(list_cg, axis=0)
lum_lut = np.where((100*lum_lut) % 1 <= 0.05, lum_lut+0.001, lum_lut)
lum_lut = np.where((100*lum_lut) % 1 >= 0.95, lum_lut-0.001, lum_lut)
lum_lut = list(np.round(lum_lut, 3))
"""
calculate average corner for lsc gain calculation further on
"""
corners = (lum_lut[0], lum_lut[15], lum_lut[-1], lum_lut[-16])
Cam.log += '\nLuminance table corners: {}'.format(corners)
l_cen = lum_lut[5*16+7]+lum_lut[5*16+8]+lum_lut[6*16+7]+lum_lut[6*16+8]
l_cen = round(l_cen/4, 3)
Cam.log += '\nLuminance table centre: {}'.format(l_cen)
av_corn = np.sum(corners)/4
return cal_cr_list, cal_cb_list, lum_lut, av_corn
"""
calculate g/r and g/b for 32x32 points arranged in a grid for a single image
"""
def alsc(Cam, Img, do_alsc_colour, plot=False):
Cam.log += '\nProcessing image: ' + Img.name
"""
get channel in correct order
"""
channels = [Img.channels[i] for i in Img.order]
"""
calculate size of single rectangle.
-(-(w-1)//32) is a ceiling division. w-1 is to deal robustly with the case
where w is a multiple of 32.
"""
w, h = Img.w/2, Img.h/2
dx, dy = int(-(-(w-1)//16)), int(-(-(h-1)//12))
"""
average the green channels into one
"""
av_ch_g = np.mean((channels[1: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 b# SPDX-License-Identifier: BSD-2-Clause
#
# Copyright (C) 2019, Raspberry Pi (Trading) Limited
#
# ctt_noise.py - camera tuning tool noise calibration
from ctt_image_load import *
import matplotlib.pyplot as plt
"""
Find noise standard deviation and fit to model:
noise std = a + b*sqrt(pixel mean)
"""
def noise(Cam, Img, plot):
Cam.log += '\nProcessing image: {}'.format(Img.name)
stds = []
means = []
"""
iterate through macbeth square patches
"""
for ch_patches in Img.patches:
for patch in ch_patches:
"""
renormalise patch
"""
patch = np.array(patch)
patch = (patch-Img.blacklevel_16)/Img.againQ8_norm
std = np.std(patch)
mean = np.mean(patch)
stds.append(std)
means.append(mean)
"""
clean data and ensure all means are above 0
"""
stds = np.array(stds)
means = np.array(means)
means = np.clip(np.array(means), 0, None)
sq_means = np.sqrt(means)
"""
least squares fit model
"""
fit = np.polyfit(sq_means, stds, 1)
Cam.log += '\nBlack level = {}'.format(Img.blacklevel_16)
Cam.log += '\nNoise profile: offset = {}'.format(int(fit[1]))
Cam.log += ' slope = {:.3f}'.format(fit[0])
"""
remove any values further than std from the fit
anomalies most likely caused by:
> ucharacteristically noisy white patch
> saturation in the white patch
"""
fit_score = np.abs(stds - fit[0]*sq_means - fit[1])
fit_std = np.std(stds)
fit_score_norm = fit_score - fit_std
anom_ind = np.where(fit_score_norm > 1)
fit_score_norm.sort()
sq_means_clean = np.delete(sq_means, anom_ind)
stds_clean = np.delete(stds, anom_ind)
removed = len(stds) - len(stds_clean)
if removed != 0:
Cam.log += '\nIdentified and removed {} anomalies.'.format(removed)
Cam.log += '\nRecalculating fit'
"""
recalculate fit with outliers removed
"""
fit = np.polyfit(sq_means_clean, stds_clean, 1)
Cam.log += '\nNoise profile: offset = {}'.format(int(fit[1]))
Cam.log += ' slope = {:.3f}'.format(fit[0])
"""
if fit const is < 0 then force through 0 by
dividing by sq_means and fitting poly order 0
"""
corrected = 0
if fit[1] < 0:
corrected = 1
ones = np.ones(len(means))
y_data = stds/sq_means
fit2 = np.polyfit(ones, y_data, 0)
Cam.log += '\nOffset below zero. Fit recalculated with zero offset'
Cam.log += '\nNoise profile: offset = 0'
Cam.log += ' slope = {:.3f}'.format(fit2[0])
# print('new fit')
# print(fit2)
"""
plot fit for debug
"""
if plot:
x = np.arange(sq_means.max()//0.88)
fit_plot = x*fit[0] + fit[1]
plt.scatter(sq_means, stds, label='data', color='blue')
plt.scatter(sq_means[anom_ind], stds[anom_ind], color='orange', label='anomalies')
plt.plot(x, fit_plot, label='fit', color='red', ls=':')
if fit[1] < 0:
fit_plot_2 = x*fit2[0]
plt.plot(x, fit_plot_2, label='fit 0 intercept', color='green', ls='--')
plt.plot(0, 0)
plt.title('Noise Plot\nImg: {}'.format(Img.str))
plt.legend(loc='upper left')
plt.xlabel('Sqrt Pixel Value')
plt.ylabel('Noise Standard Deviation')
plt.grid()
plt.show()
"""
End of plotting code
"""
"""
format output to include forced 0 constant
"""
Cam.log += '\n'
if corrected:
fit = [fit2[0], 0]
return fit
else:
return fit