PentHertz
/
deep-tempest
mirror of https://github.com/PentHertz/deep-tempest.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
deep-tempest/gr-tempest/python/DTutils.py

735 lines
22 KiB
Python
Raw Permalink Blame History

import numpy as np
from matplotlib import pyplot as plt
from scipy import signal
import cv2 as cv
from scipy.spatial import distance_matrix
def autocorr(x):
"""Compute autocorrelation function of 1-D array
Input:
x: 1-D array
Output:
autocorr: autocorrelation function of x
"""
# Use FFT method, which has more computing efectiveness for 1-D numpy arrays
autocorr = signal.correlate(x,x,mode='full', method= 'fft')
# Fix some shifts due FFT
half_idx =int(autocorr.size/2)
max_ind = np.argmax(autocorr[half_idx:])+half_idx
autocorr = autocorr[max_ind:]
# Normalise output
return autocorr/autocorr[0]
def uint8_to_binarray(integer):
"""Convert integer into fixed-length 8-bit binary array. LSB in [0].
Extended and modified code from https://github.com/projf/display_controller/blob/master/model/tmds.py
"""
b_array = [int(i) for i in reversed(bin(integer)[2:])]
b_array += [0]*(8-len(b_array))
return b_array
def uint16_to_binarray(integer):
"""Convert integer into fixed-length 16-bit binary array. LSB in [0].
Extended and modified code from https://github.com/projf/display_controller/blob/master/model/tmds.py
"""
b_array = [int(i) for i in reversed(bin(integer)[2:])]
b_array += [0]*(16-len(b_array))
return b_array
def binarray_to_uint(binarray):
array = binarray[::-1]
num = array[0]
for n in range(1,len(binarray)):
num = (num << 1) + array[n]
return num
def TMDS_pixel (pix,cnt=0):
"""8bit pixel TMDS coding
Inputs:
- pix: 8-bit pixel
- cnt: 0's and 1's balance. Default in 0 (balanced)
Outputs:
- pix_out: TDMS coded 16-bit pixel (only 10 useful)
- cnt: 0's and 1's balance updated with new pixel coding
"""
# Convert 8-bit pixel to binary list D
D = uint8_to_binarray(pix)
# Initialize output q
qm = [D[0]]
# 1's unbalanced condition at current pixel
N1_D = np.sum(D)
if N1_D>4 or (N1_D==4 and not(D[0])):
# XNOR of consecutive bits
for k in range(1,8):
qm.append( not(qm[k-1] ^ D[k]) )
qm.append(0)
else:
# XOR of consecutive bits
for k in range(1,8):
qm.append( qm[k-1] ^ D[k] )
qm.append(1)
# Initialize output qout
qout = qm.copy()
# Unbalanced condition with previous and current pixels
N1_qm = np.sum(qm[:8])
N0_qm = 8 - N1_qm
if cnt==0 or N1_qm==4:
qout.append(not(qm[8]))
qout[8] = qm[8]
qout[:8]=qm[:8] if qm[8] else np.logical_not(qm[:8])
if not(qm[8]):
cnt += N0_qm - N1_qm
else:
cnt += N1_qm - N0_qm
else:
if (cnt>0 and N1_qm>4) or (cnt<0 and N1_qm<4):
qout.append(1)
qout[8] = qm[8]
qout[:8] = np.logical_not(qm[:8])
cnt += 2*qm[8] + N0_qm - N1_qm
else:
qout.append(0)
qout[8] = qm[8]
qout[:8] = qm[:8]
cnt += -2*(not(qm[8])) + N1_qm - N0_qm
# Return the TMDS coded pixel as uint and 0's y 1's balance
return binarray_to_uint(qout), cnt
def TMDS_encoding_original (I, blanking = False):
"""TMDS image coding
Inputs:
- I: 2-D image array
- blanking: Boolean that specifies if horizontal and vertical blanking is applied
Output:
- I_c: TDMS coded 16-bit image (only 10 useful)
"""
# Create "ghost dimension" if I is gray-scale image (not RGB)
if len(I.shape)!= 3:
I = np.repeat(I[:, :, np.newaxis], 3, axis=2).astype('uint8')
chs = 3
# Get image resolution
v_in, h_in = I.shape[:2]
if blanking:
# Get blanking resolution for input image
v = (v_in==1080)*1125 + (v_in==720)*750 + (v_in==600)*628 + (v_in==480)*525
h = (h_in==1920)*2200 + (h_in==1280)*1650 + (h_in==800)*1056 + (h_in==640)*800
vdiff = v - v_in
hdiff = h - h_in
# Create image with blanking and change type to uint16
# Assuming the blanking corresponds to 10bit number [0, 0, 1, 0, 1, 0, 1, 0, 1, 1] (LSB first)
I_c = 852*np.ones((v,h,chs)).astype('uint16')
else:
v_diff = 0
h_diff = 0
I_c = np.zeros((v_in,h_in,chs)).astype('uint16')
# Iterate over channels and pixels
for c in range(chs):
for i in range(v_in):
cnt=[0,0,0]
for j in range(h_in):
# Get pixel and code it TMDS between blanking
pix = I[i,j,c]
I_c[i + v_diff//2 , j + h_diff//2, c], cnt[c] = TMDS_pixel (pix,cnt[c])
return I_c
def TMDS_pixel_cntdiff (pix,cnt=0):
"""8bit pixel TMDS coding
Inputs:
- pix: 8-bit pixel
- cnt: 0's and 1's balance. Default in 0 (balanced)
Outputs:
- pix_out: TDMS coded 16-bit pixel (only 10 useful)
- cntdiff: balance difference given by the actual coded pixel
"""
# Convert 8-bit pixel to binary list D
D = uint8_to_binarray(pix)
# Initialize output q
qm = [D[0]]
# 1's unbalanced condition at current pixelo
N1_D = np.sum(D)
if N1_D>4 or (N1_D==4 and not(D[0])):
# XNOR of consecutive bits
for k in range(1,8):
qm.append( not(qm[k-1] ^ D[k]) )
qm.append(0)
else:
# XOR of consecutive bits
for k in range(1,8):
qm.append( qm[k-1] ^ D[k] )
qm.append(1)
# Initialize output qout
qout = qm.copy()
# Unbalanced condition with previous and current pixels
N1_qm = np.sum(qm[:8])
N0_qm = 8 - N1_qm
if cnt==0 or N1_qm==4:
qout.append(not(qm[8]))
qout[8]=qm[8]
qout[:8]=qm[:8] if qm[8] else [not(val) for val in qm[:8]]
if not(qm[8]):
cnt_diff = N0_qm - N1_qm
else:
cnt_diff = N1_qm - N0_qm
else:
if (cnt>0 and N1_qm>4) or (cnt<0 and N1_qm<4):
qout.append(1)
qout[8]=qm[8]
qout[:8] = [not(val) for val in qm[:8]]
cnt_diff = 2*qm[8] +N0_qm -N1_qm
else:
qout.append(0)
qout[8]=qm[8]
qout[:8] = qm[:8]
cnt_diff = -2*(not(qm[8])) + N1_qm - N0_qm
# Return the TMDS coded pixel as uint and 0's y 1's balance difference
uint_out = binarray_to_uint(qout)
return uint_out, cnt_diff
### Create TMDS LookUp Tables for fast encoding (3 times faster than the other implementation)
byte_range = np.arange(256)
# Initialize pixel coding and cnt-difference arrays
TMDS_pix_table = np.zeros((256,3),dtype='uint16')
TMDS_cntdiff_table = np.zeros((256,3),dtype='int8')
for byte in byte_range:
p0,p_null, p1 = TMDS_pixel_cntdiff(byte,-1),TMDS_pixel_cntdiff(byte,0),TMDS_pixel_cntdiff(byte,1) # 0's and 1's unbalance respect.
TMDS_pix_table[byte,0] = p0[0]
TMDS_pix_table[byte,1] = p_null[0]
TMDS_pix_table[byte,2] = p1[0]
TMDS_cntdiff_table[byte,0] = p0[1]
TMDS_cntdiff_table[byte,1] = p_null[1]
TMDS_cntdiff_table[byte,2] = p1[1]
def pixel_fastencoding(pix,cnt_prev=0):
"""8bit pixel TMDS fast coding
Inputs:
- pix: 8-bit pixel
- cnt: 0's and 1's balance. Default in 0 (balanced)
Outputs:
- pix_out: TDMS coded 16-bit pixel (only 10 useful)
- cnt: 0's and 1's balance updated with new pixel coding
"""
balance_idx = int(np.sign(cnt_prev))+1
pix_out = TMDS_pix_table[pix,balance_idx]
cnt = cnt_prev + TMDS_cntdiff_table[pix,balance_idx]
return pix_out, cnt
def TMDS_blanking (h_total, v_total, h_active, v_active, h_front_porch, v_front_porch, h_back_porch, v_back_porch):
# Initialize blanking image
img_blank = np.zeros((v_total,h_total))
# Get the total blanking on vertical an horizontal axis
h_blank = h_total - h_active
v_blank = v_total - v_active
# (C1,C0)=(0,0) region
img_blank[:v_front_porch,:h_front_porch] = 0b1101010100
img_blank[:v_front_porch,h_blank-h_back_porch:] = 0b1101010100
img_blank[v_blank-v_back_porch:v_blank,:h_front_porch] = 0b1101010100
img_blank[v_blank-v_back_porch:v_blank,h_blank-h_back_porch:] = 0b1101010100
img_blank[v_blank:,:h_blank] = 0b1101010100
# (C1,C0)=(0,1) region
img_blank[:v_front_porch,h_front_porch:h_blank-h_back_porch] = 0b0010101011
img_blank[v_blank-v_back_porch:,h_front_porch:h_blank-h_back_porch] = 0b0010101011
# (C1,C0)=(1,0) region
img_blank[v_front_porch:v_blank-v_back_porch,:h_front_porch] = 0b0101010100
img_blank[v_front_porch:v_blank-v_back_porch,h_blank-h_back_porch:] = 0b0101010100
# (C1,C0)=(1,1) region
img_blank[v_front_porch:v_blank-v_back_porch,v_front_porch:h_blank-h_back_porch] = 0b1010101011
return(img_blank)
def TMDS_encoding (I, blanking = False):
"""TMDS image coding
Inputs:
- I: 2D/3D image array (v_size, h_size, channels)
- blanking: Boolean that specifies if horizontal and vertical blanking is applied or not
Output:
- I_c: 3D TDMS coded 16-bit (only 10 useful) image array
"""
# Create "ghost dimension" if I is gray-scale image (not RGB)
if len(I.shape)!= 3:
# Gray-scale image
I = np.repeat(I[:, :, np.newaxis], 3, axis=2).astype('uint8')
chs = 1
else:
# RGB image
chs = 3
# Get image resolution
v_in, h_in = I.shape[:2]
if blanking:
# Get blanking resolution for input image
v = (v_in==1080)*1125 + (v_in==900)*1000 + (v_in==720)*750 + (v_in==600)*628 + (v_in==480)*525
h = (h_in==1920)*2200 + (h_in==1600)*1800 + (h_in==1280)*1650 + (h_in==800)*1056 + (h_in==640)*800
v_diff = v - v_in
h_diff = h - h_in
v_front_porch = (v_in==1080)*4 + (v_in==900)*1 + (v_in==720)*5 + (v_in==600)*1 + (v_in==480)*2
v_back_porch = (v_in==1080)*36 + (v_in==900)*96 + (v_in==720)*20 + (v_in==600)*23 + (v_in==480)*25
h_front_porch = (h_in==1920)*88 + (h_in==1600)*24 + (h_in==1280)*110 + (h_in==800)*40 + (h_in==640)*8
h_back_porch = (h_in==1920)*148 + (h_in==1600)*96 + (h_in==1280)*220 + (h_in==800)*88 + (h_in==640)*40
# Create image with blanking and change type to uint16
# Assuming the blanking corresponds to 10bit number
# [0, 0, 1, 0, 1, 0, 1, 0, 1, 1] (LSB first) for channels R and G
I_c = 852*np.ones((v,h,chs)).astype('uint16')
I_c[:,:,2] = TMDS_blanking(h_total=h, v_total=v, h_active=h_in, v_active=v_in,
h_front_porch=h_front_porch, v_front_porch=v_front_porch, h_back_porch=h_back_porch, v_back_porch=v_back_porch)
else:
v_diff = 0
h_diff = 0
I_c = np.zeros((v_in,h_in,chs)).astype('uint16')
# Iterate over channels and pixels
for c in range(chs):
for i in range(v_in):
cnt = [0,0,0]
for j in range(h_in):
# Get pixel and code it TMDS between blanking
pix = I[i,j,c]
I_c[i + v_diff, j + h_diff, c], cnt[c] = pixel_fastencoding (pix,cnt[c])
return I_c
def DecTMDS_pixel (pix):
"""10-bit pixel TMDS decoding
Inputs:
- pix: 16-bit pixel (only 10 first bits useful)
Output:
- pix_out: 8-bit TMDS decoded pixel
"""
D = uint16_to_binarray(pix)[:10]
if D[9]:
D[:8] = np.logical_not(D[:8])
Q = D.copy()[:8]
if D[8]:
for k in range(1,8):
Q[k] = D[k] ^ D[k-1]
else:
for k in range(1,8):
Q[k] = not(D[k] ^ D[k-1])
# Return pixel as uint
return binarray_to_uint(Q)
def TMDS_decoding (Ic):
"""Image TMDS decoding
Inputs:
- Ic: TMDS coded image
Output:
- Idec: 8-bit decoded image
"""
# Create "ghost dimension" if gray-scale image (not RGB)
if len(Ic.shape)!= 3:
Ic = Ic.reshape(Ic.shape[0],Ic.shape[1],1)
Idec = Ic.copy()
# Get image dimensions
Nx, Ny, Nz = Ic.shape
# Iterate over channels and pixels
for c in np.arange(Nz):
for i in np.arange(Nx):
for j in np.arange(Ny):
# Get pixel and use TMDS decoding
pix = Ic[i,j,c]
Idec[i,j,c] = DecTMDS_pixel (pix)
return Idec
def TMDS_serial(I):
'''
Serialize an image as an 1D binary array given a 10bit pixel value.
Inputs:
- I: TMDS image to serialize. Pixel values must be between 0 and 1023
Output:
- Iserials: 1D binary array per image channel which represents
the voltage value to be transmitted
'''
assert np.min(I)>=0 and np.max(I)<= 1023, "Pixel values must be between 0 and 1023"
# Initialize lists per channel
Iserials = []
n_rows,n_columns, n_channels = I.shape
# Iterate over pixel
for c in range(n_channels):
channel_list = []
for i in range(n_rows):
for j in range(n_columns):
# Get pixel value and cast it as binary string
binstring = bin(I[i,j,c])[2:]
# Fill string with 0's to get length 10
binstring = '0'*(10-len(binstring))+binstring
# Re-order string for LSB first
binstring = binstring[::-1]
binarray = list(binstring)
# Extend the bit stream
channel_list.extend(binarray)
Iserials.append(channel_list)
# Digital to analog value mapping: [0,1]-->[-A,A] (A=1)
Iserials = np.sum(2*np.array(Iserials,dtype='int32') - 1, axis=0)
del(channel_list)
return Iserials
def remove_outliers(I, radius=3, threshold=20):
"""
Replaces a pixel by the median of the pixels in the surrounding if it deviates from the median by more than a certain value (the threshold).
"""
# Copy input
I_output = I.copy()
# Apply median filter
I_median = cv.medianBlur(I,radius)
# Replace with median value where difference with median exceedes threshold
where_replace = np.abs(I-I_median) > threshold
I_output[where_replace] = I_median[where_replace]
return I_output
def adjust_dynamic_range(I):
I_output = I.astype('float32')
I_output = (I_output - I_output.min()) / (I_output.max() - I_output.min())
I_output = 255 * I_output
return I_output.astype('uint8')
def find_intersection(line1, line2):
"""Finds the intersection of two lines given in Hesse normal form.
Returns closest integer pixel locations.
See https://stackoverflow.com/a/383527/5087436
"""
rho1, theta1 = line1
rho2, theta2 = line2
A = np.array([
[np.cos(theta1), np.sin(theta1)],
[np.cos(theta2), np.sin(theta2)]
])
b = np.array([[rho1], [rho2]])
x0, y0 = np.linalg.solve(A, b)
x0, y0 = int(np.round(x0)), int(np.round(y0))
# Return as row-column coordinates
return [y0, x0]
def apply_blanking_shift(I, h_active=1600, v_active=900,
h_blanking=200,v_blanking=100,
pad_len=300, debug=False):
"""
Correct capture shift to center image using VESA blanking information.
The image must not have geometric distortion (generaly caused by sampling error)
"""
if debug:
# Show original image
plt.figure(figsize=(12,10))
plt.title(f'Original image ')
plt.imshow(I)
plt.axis('off')
plt.show()
# Color to RGB
I_gray = cv.cvtColor(I,cv.COLOR_BGR2GRAY)
# I_gray = cv.medianBlur(I_gray,3)
I_gray = cv.GaussianBlur(I_gray,ksize=(5,5),sigmaX=3,sigmaY=3)
# I_gray = cv.blur(I_gray,(3,3))
# Wrap-padding image to get whole blanking pattern
pad_size_gray = ((pad_len,0),(pad_len,0))
pad_size = ((pad_len,0),(pad_len,0),(0,0))
I_gray = np.pad(I_gray,pad_size_gray,'wrap')
I_pad = np.pad(I,pad_size,'wrap')
# Edge-detector with Canny filter with empirical parameters
I_edge = cv.Canny(I_gray,40,50,apertureSize = 3) # sin blur
if debug:
# Show padded image and edges
plt.figure()
plt.title(f'Original image wrap-padded {pad_len} pixels up and left sided')
plt.imshow(I_pad)
plt.axis('off')
plt.show()
plt.figure()
plt.title('Canny edged image')
plt.imshow(I_edge,cmap='gray')
plt.axis('off')
plt.show()
# Hough Transform resolution and minimum votes threshold
theta_step = np.pi/2
rho_step = 1
# votes_thrs = 255 # sin blur
votes_thrs = 80
rho_max = np.sqrt(I_gray.shape[0]**2 + I_gray.shape[1]**2)
# Find Hough Transform
lines = cv.HoughLines(I_edge, rho_step, theta_step, votes_thrs, None, 0, 0)
if debug:
# Show detected lines
I_lines = I_gray.copy()
for line in lines:
rho = line[0][0]
theta = line[0][1]
a = np.cos(theta)
b = np.sin(theta)
x0 = a * rho
y0 = b * rho
pt1 = (int(x0 + rho_max*(-b)), int(y0 + rho_max*(a)))
pt2 = (int(x0 - rho_max*(-b)), int(y0 - rho_max*(a)))
cv.line(I_lines, pt1, pt2, (255,0,0), 1, cv.LINE_AA)
plt.figure()
plt.title('All detected lines')
plt.imshow(I_lines)
plt.axis('off')
plt.show()
# Angle and rho arrays
try:
lines_angles = lines[:,0,1]
except:
print("No lines detected. Returning original image")
return I
# Find unique lines angles detected
unique_angles = np.unique(lines_angles)
blankings = [h_blanking, v_blanking]
# Initiate blanking lines variable
blanking_lines = []
for angle, blanking in zip(unique_angles, blankings):
# Keep lines with certain angle
angle_lines = lines[lines[:,0,1]==angle]
angle_lines = angle_lines[:,0,:]
# Keep the pair with rho distance that equals blanking
# First, compute the distance matrix over all rhos:
rho_angle_lines = np.array(angle_lines[:,0]).reshape(-1,1)
rho_distances = np.abs(np.abs(distance_matrix(rho_angle_lines, rho_angle_lines)) - blanking)
# Find minimum index
pair_lines_idx = np.unravel_index(rho_distances.argmin(), rho_distances.shape)
# Add blanking limit line
for line in pair_lines_idx:
blanking_lines.append(angle_lines[line])
# List to array
blanking_lines = np.array(blanking_lines)
# Show blanking limit lines
if debug:
I_lines = I_gray.copy()
for line in blanking_lines:
rho = line[0]
theta = line[1]
a = np.cos(theta)
b = np.sin(theta)
x0 = a * rho
y0 = b * rho
pt1 = (int(x0 + rho_max*(-b)), int(y0 + rho_max*(a)))
pt2 = (int(x0 - rho_max*(-b)), int(y0 - rho_max*(a)))
cv.line(I_lines, pt1, pt2, (255,0,0), 3, cv.LINE_AA)
plt.figure(figsize=(12,10))
plt.title('Blanking limit lines')
plt.imshow(I_lines)
plt.axis('off')
plt.show()
# Initialize top-left corner blanking lines to find intersection
# These lines statisfies to be the ones with bigger rho value
blanking_start = []
unique_angles = np.unique(blanking_lines[:,1])
for angulo in unique_angles:
angle_lines = blanking_lines[blanking_lines[:,1]==angulo]
max_rho_line = angle_lines[np.argmax(angle_lines[:,0])]
blanking_start.append(max_rho_line)
if len(blanking_start)!=2:
print("Not enough lines detected. Returning original image")
return I
if debug:
I_lines = I_gray.copy()
for line in blanking_start:
rho = line[0]
theta = line[1]
a = np.cos(theta)
b = np.sin(theta)
x0 = a * rho
y0 = b * rho
pt1 = (int(x0 + rho_max*(-b)), int(y0 + rho_max*(a)))
pt2 = (int(x0 - rho_max*(-b)), int(y0 - rho_max*(a)))
cv.line(I_lines, pt1, pt2, (255,0,0), 3, cv.LINE_AA)
plt.figure(figsize=(12,10))
plt.title('Blanking start lines')
plt.imshow(I_lines)
plt.axis('off')
plt.show()
# Find blanking start lines intersection coordinates
x_shift, y_shift = find_intersection(blanking_start[0],blanking_start[1])
# Adjust to active image only (remove all blanking)
I_shift = I_pad[pad_len:,pad_len:]
I_shift = np.roll(I_shift, -x_shift+pad_len, axis=0)
I_shift = np.roll(I_shift, -y_shift+pad_len, axis=1)
I_shift = I_shift[:v_active,:h_active]
if debug:
plt.figure(figsize=(12,10))
plt.title('Blanking removed image')
plt.imshow(I_shift)
plt.axis('off')
plt.show()
return I_shift
def preprocess_raw_capture(I, h_active, v_active,
h_blanking, v_blanking, debug=False):
"""
Center raw captured image, filter noise and adjust the contrast
"""
# Center image. Returns input image if no sufficient lines where detected
# In the latter case, use image as is without centering
is_centered = True
I_shift_fix = apply_blanking_shift(I,
h_active=h_active, v_active=v_active,
h_blanking=h_blanking, v_blanking=v_blanking,
debug=debug
)
# Remove outliers with median thresholding heuristic
# Default: radius=3, threshold=20
I_no_outliers = remove_outliers(I_shift_fix)
# Stretch dynamic range to [0,255]
I_out = adjust_dynamic_range(I_no_outliers)
if debug:
is_centered = (I == I_shift_fix)
plt.figure(figsize=(12,10))
ax0 = plt.subplot(3,1,1)
ax0.imshow(I_shift_fix, interpolation='none')
ax0.set_title('Centered image'*is_centered + 'Image'*(~is_centered))
ax0.axis('off')
ax1 = plt.subplot(3,1,2, sharex=ax0, sharey=ax0)
ax1.imshow(I_no_outliers, interpolation='none')
ax1.set_title('Outliers removed')
ax1.axis('off')
ax1 = plt.subplot(3,1,3, sharex=ax0, sharey=ax0)
ax1.imshow(I_out, interpolation='none')
ax1.set_title('Contrast adjusted')
ax1.axis('off')
plt.show()
return I_out
Reference in New Issue View Git Blame Copy Permalink