"""
submodule focussed completely on images, all functions with prefix "img"
take an image as main input and return a processed image
"""
import copy
from collections.abc import Iterable
import cv2
import matplotlib.pyplot as plt
import numpy as np
from numba import njit, prange
from scipy.special import erf
from skimage import exposure
#%% padding manipulation
[docs]
def img_padding_attenuation(padded_img,pad_width,mode="linear",
parameters=None,plot=False):
"""extending the functionality of np.pad (following same syntax: pad_width)
or cv2.copyMakeBorder allowing custom attenuation (value reduction towards 0)
of added padding (borders)
modes: 'linear', 'exponential', 'inverse_exponential' and 'sigmoid_normalCDF'
"""
if parameters is None:
parameters=dict()
if not isinstance(pad_width, Iterable):
pw=pad_width,pad_width,pad_width,pad_width
elif not isinstance(pad_width[0], Iterable):
pw=pad_width[0],pad_width[1],pad_width[0],pad_width[1]
else:
pw=pad_width[0][0],pad_width[0][1],pad_width[1][0],pad_width[1][1]
imgshape=padded_img.shape
seqs=[]
if mode=="linear":
if "end_value" not in parameters:
parameters["end_value"]=0
for i in range(4):
seqs.append(np.linspace(1,parameters["end_value"],pw[i]))
if mode=="exponential" or mode=="inverse_exponential":
if "end_value" not in parameters:
parameters["end_value"]=0
if "lambda" not in parameters:
parameters["lambda"]=3
for i in range(4):
seq=np.exp(-np.linspace(0,parameters["lambda"],pw[i]))
factor=(1-parameters["end_value"])/(1-seq[-1])
seq*=factor
seq-=seq[-1]-parameters["end_value"]
seqs.append(seq)
if mode=="inverse_exponential":
for i in range(4):
seqs[i]=1+parameters["end_value"]-seqs[i][::-1]
if mode=="sigmoid_normalCDF":
if "end_value" not in parameters:
parameters["end_value"]=0
if "sigma" not in parameters:
parameters["sigma"]=0.5
for i in range(4):
seq=np.linspace(1,-1,pw[i])
seq=0.5+0.5*erf(seq/(np.sqrt(2)*parameters["sigma"]))
factor=(1-parameters["end_value"])/(1-2*seq[-1])
seq*=factor
seq-=seq[-1]-parameters["end_value"]
seqs.append(seq)
if plot:
width=np.min([len(seq) for seq in seqs])
plt.plot(np.arange(len(seqs[0]))[::-1],seqs[0],c="r")
plt.plot(np.arange(len(seqs[2]))[::-1],seqs[2],c="b")
plt.plot(np.arange(width)+len(seqs[0]),np.ones(width),'--',c='r',alpha=0.5)
plt.plot(np.arange(width)+len(seqs[2]),np.ones(width),'-.',c="b",alpha=0.5)
plt.plot(np.arange(len(seqs[1]))+len(seqs[0])+width,seqs[1],c="r",label="vertical")
plt.plot(np.arange(len(seqs[3]))+len(seqs[2])+width,seqs[3],c="b",label="horizontal")
plt.legend()
n_cols=imgshape[1]
b_up = np.tile(seqs[0][::-1, None], (1, n_cols))
b_down= np.tile(seqs[1][:, None], (1, n_cols))
n_rows=imgshape[0]
b_left = np.tile(seqs[2][None,::-1], ( n_rows,1))
b_right = np.tile(seqs[3][None,:], ( n_rows,1))
atten_img=padded_img*1.
atten_img[-pw[1]:,:]*=b_down
atten_img[:,-pw[3]:]*=b_right
atten_img[:pw[0],:]*=b_up
atten_img[:,:pw[2]]*=b_left
return atten_img
#%% autoclipping
[docs]
def img_autoclip(img,ratio=0.001):
"""
auto clipping removing small and big outliers
Args:
img (TYPE): DESCRIPTION.
ratio (TYPE, optional): DESCRIPTION. Defaults to 0.001.
Returns:
TYPE: DESCRIPTION.
"""
imgdata=np.sort(np.ravel(img))
threshold=int(len(imgdata)*ratio)
return np.clip(img,imgdata[threshold],imgdata[-threshold])
#%% morphLaplace
[docs]
def img_morphLaplace(image, kernel):
"""
morphological Laplace Filter
Args:
image (MxN array_like):
np.uint8.
kernel (TYPE):
DESCRIPTION.
Returns:
filtered_image (MxN array_like):
np.uint8.
"""
return cv2.erode(image, kernel) + cv2.dilate(image, kernel) - 2 * image - 128
#%% gammaCorrection
[docs]
def img_gammaCorrection(img, gamma):
"""
transformation adjusting the gamma level
(gamma=1 means no change)
Args:
img (MxN array_like):
np.uint8.
gamma (float):
gamma value between 0 and ~10.
Returns:
gamma_transformed_image (MxN array_like):
np.uint8.
"""
invGamma = 1 / gamma
table = [((i / 255) ** invGamma) * 255 for i in range(256)]
table = np.array(table, np.uint8)
return cv2.LUT(img, table)
#%% format image dtypes
[docs]
def img_to_uint8(img,newmax=255,imgmin=None,imgmax=None):
"""
transform contrast range to unsigned integer 8bit
Args:
img (MxN array_like):
input image.
newmax (int, optional):
optionally reduce contrast range,
by setting a lower, than datatype given, threshold to the maximum value.
Defaults to 255.
imgmin (float, optional):
optionally set the lower limit,
of the original image contrast,
imgmin=0 preserves the original relative contrast range offset
imgmin=None sets the contrast range offset to 0
Defaults to None.
imgmax (float, optional):
optionally set the upper limit
of the original image contrast,
imgmax=None uses the biggest occuring value of the image
Defaults to None.
Returns:
datatype_conform_image (MxN array_like):
np.uint8.
"""
if imgmax is None:
imgmax = np.max(img)
if imgmin is None:
imgmin = np.min(img)
newimg = img - imgmin
newimgmax=imgmax - imgmin
return ((newimg / newimgmax) * (newmax+0.5)).astype(np.uint8)
[docs]
def img_to_uint8_fast(img):
newimg=np.empty(img.shape,dtype=np.uint8)
img_to_uint8_numba(img, newimg)
return newimg
[docs]
@njit(parallel=True)
def img_to_uint8_numba(img,newimg):
N0,N1=img.shape
imgmin=np.min(img)
imgmax=np.max(img)-imgmin
for i in prange(N0):
for j in range(N1):
newimg[i,j]=(img[i,j]-imgmin)/imgmax *255.5
[docs]
def img_normalize(img):
#if imgmin is None:
imgmin = np.min(img)
newimg = img - imgmin
return newimg / np.max(newimg)
[docs]
def img_normalize_fast(img):
newimg=np.empty(img.shape,dtype=np.float32)
img_normalize_numba(img, newimg)
return newimg
[docs]
@njit(parallel=True)
def img_normalize_numba(img,newimg):
N0,N1=img.shape
imgmin=np.min(img)
imgmax=np.max(img)-imgmin
for i in prange(N0):
for j in range(N1):
newimg[i,j]=(img[i,j]-imgmin)/imgmax
[docs]
def img_to_uint16(img,newmax=65535,imgmin=None):
"""
transform contrast range to unsigned integer 16bit
Args:
img (MxN array_like):
input image.
newmax (int, optional):
optionally reduce contrast range,
by setting a lower, than datatype given, threshold to the maximum value.
Defaults to 65535.
imgmin (int, optional):
optionally set the lower limit,
of the original image contrast,
imgmin=0 preserves the original relative contrast range offset
imgmin=None sets the contrast range offset to 0
Defaults to None.
Returns:
datatype_conform_image (MxN array_like):
np.uint16.
"""
if imgmin is None:
imgmin = np.min(img)
newimg = img - imgmin
return (newimg / np.max(newimg) * (newmax+0.5)).astype(np.uint16)
[docs]
def img_to_int8(img,newmax=127,imgmin=None):
"""
transform contrast range to signed integer 8bit
Args:
img (MxN array_like):
input image.
newmax (int, optional):
optionally reduce contrast range,
by setting a lower, than datatype given, threshold to the maximum value.
Defaults to 127.
imgmin (int, optional):
optionally set the lower limit,
of the original image contrast,
imgmin=0 preserves the original relative contrast range offset
imgmin=None sets the contrast range offset to 0
Defaults to None.
Returns:
datatype_conform_image (MxN array_like):
np.int8.
"""
newmax+=128
if imgmin is None:
imgmin = np.min(img)
newimg = img - imgmin
return (newimg / np.max(newimg) * (newmax+0.5) -128).astype(np.int8)
[docs]
def img_to_int16(img,newmax=32767,imgmin=None):
"""
transform contrast range to signed integer 16bit
Args:
img (MxN array_like):
input image.
newmax (int, optional):
optionally reduce contrast range,
by setting a lower, than datatype given, threshold to the maximum value.
Defaults to 32767.
imgmin (int, optional):
optionally set the lower limit,
of the original image contrast,
imgmin=0 preserves the original relative contrast range offset
imgmin=None sets the contrast range offset to 0
Defaults to None.
Returns:
datatype_conform_image (MxN array_like):
np.int16.
"""
newmax+=32768
if imgmin is None:
imgmin = np.min(img)
newimg = img-imgmin
return (newimg / np.max(newimg) * (newmax+0.5) -32768).astype(np.int16)
[docs]
def img_to_half_int8(img):
"""
transform contrast to half the range of signed integer 8bit,
so values are between -64 and 63 for possible range -128 to 127
Args:
img (MxN array_like):
input image.
Returns:
datatype_conform_image (MxN array_like):
np.int8.
"""
img -= np.min(img)
return (img / np.max(img) * 63.5 -64).astype(np.int8)
[docs]
def img_to_half_int16(img):
"""
transform contrast to half the range of signed integer 16bit,
so values are between -16384 and 16383 for possible range -32768 to 32767
Args:
img (MxN array_like):
input image.
Returns:
datatype_conform_image (MxN array_like):
np.int16.
"""
img -= np.min(img)
return (img / np.max(img) * 16383.5-16384).astype(np.int16)
[docs]
def img_gray_to_rgba(img):
"""
transform grayscale (uint8) to rgba (uint8) image
Args:
img (TYPE): DESCRIPTION.
Returns:
rgba_img (TYPE): DESCRIPTION.
"""
rgba_img=np.dstack((img,img,img,np.zeros(img.shape,dtype=np.uint8)+255))
return rgba_img
[docs]
def img_single_to_double_channel(img):
#gray_alpha
if len(img.shape)==3:
return img
else:
return np.dstack((img,np.ones(img.shape,dtype=img.dtype)))
[docs]
def img_add_weighted_gray_alpha(img1,img2):
newimg=np.empty(img1.shape,dtype=np.uint8)
img_add_weighted_gray_alpha_numba(img1,img2,newimg)
#totalalpha= np.zeros(img1.shape[:2],dtype=np.uint16)
#totalalpha+=img1[:,:,1]
#totalalpha+=img2[:,:,1]
#totalalpha[totalalpha==0]=1
#alpha1=img1[:,:,1]/totalalpha
#alpha2=img2[:,:,1]/totalalpha
#newimg[:,:,0] = alpha1*img1[:,:,0]+alpha2*img2[:,:,0]
#newimg[:,:,1] = (alpha1+alpha2) * 255.5
return newimg
[docs]
@njit(parallel=True)
def img_add_weighted_gray_alpha_numba(img1,img2,newimg):
N0=img1.shape[0]
N1=img1.shape[1]
for i in prange(N0):
for j in range(N1):
totalalpha=img1[i,j,1]+img2[i,j,1]
if totalalpha ==0:
totalalpha=1
alpha1=img1[i,j,1]/totalalpha
alpha2=img2[i,j,1]/totalalpha
newimg[i,j,0]=alpha1*img1[i,j,0]+alpha2*img2[i,j,0]
newimg[i,j,1]=max(img1[i,j,1],img2[i,j,1])#(alpha1+alpha2)*255.5#(img1[i,j,1]+img2[i,j,1])/2#
#totalalpha= numba.uint16[:,:]#np.zeros(img1.shape[:2],dtype=np.uint16)
#totalalpha+=img1[:,:,1]
#totalalpha+=img2[:,:,1]
#totalalpha[totalalpha==0]=1
#alpha1=img1[:,:,1]/totalalpha
#alpha2=img2[:,:,1]/totalalpha
#newimg[:,:,0] = alpha1*img1[:,:,0]+alpha2*img2[:,:,0]
#newimg[:,:,1] = (alpha1+alpha2) * 255.5
#return newimg
[docs]
def img_add_weighted_rgba(img1,img2):
newimg1=np.copy(img1)
newimg2=np.copy(img2)
totalalpha= np.zeros(img1.shape[:2],dtype=np.uint16)
totalalpha+=img1[:,:,3]
totalalpha+=img2[:,:,3]
totalalpha[totalalpha==0]=1
alpha1=img1[:,:,3]/totalalpha
alpha2=img2[:,:,3]/totalalpha
for i in range(3):
newimg1[:,:,i] = alpha1*newimg1[:,:,i]
newimg2[:,:,i] = alpha2*newimg2[:,:,i]
newimg1[:,:,3] = alpha1 * 255.5
newimg2[:,:,3] = alpha2 * 255.5
return cv2.add(newimg1,newimg2)
#%% noise_line_suppression
[docs]
def img_noise_line_suppression(image, ksize):
"""
morphological opening with a horizontal line as structuring element
(first erode, than dilate ; only horizotal lines with length ksize remain)
Args:
image (array_like):
input.
ksize (int):
uneven integer.
Returns:
processed_image (array_like):
output.
"""
erod_img = cv2.erode(image, np.ones([1, ksize]))
return cv2.dilate(erod_img, np.ones([1, ksize]))
#%% rebin
[docs]
def img_rebin_by_mean(image, new_shape):
"""
reduce the resolution of an image MxN to mxn by taking an average,
whereby M and N must be multiples of m and n
Args:
image (MxN array_like):
image.
new_shape (tuple):
containing two integers with the new shape.
Returns:
rebinned_image (mxn array_like):
smaller image.
"""
if image.shape[0]%new_shape[0] !=0 or image.shape[1]%new_shape[1]:
raise ValueError("image shape is not a multiple of new_shape")
shape = (
new_shape[0],
image.shape[0] // new_shape[0],
new_shape[1],
image.shape[1] // new_shape[1],
)
return image.reshape(shape).mean(-1).mean(1)
#%% make_square
[docs]
def img_make_square(image, startindex=None):
"""
crops the largest square image from the original, by default from the center.
The position of the cropped square can be specified via startindex,
moving the frame from the upper left corner at startindex=0
to the lower right corner at startindex=abs(M-N)
Args:
image (MxN array_like):
DESCRIPTION.
startindex (int, optional):
must be within 0 <= startindex <= abs(M-N).
Defaults to None.
Returns:
square_image (array_like):
either MxM or NxN array.
"""
ishape = image.shape
index_small, index_big = np.argsort(ishape)
roi = np.zeros([2, 2], dtype=int)
roi[index_small] = [0, ishape[index_small]]
delta = np.abs(ishape[1] - ishape[0])
if startindex is None:
startindex = np.floor(delta / 2)
else:
if startindex > delta or startindex < 0:
print("Error: Invalid startindex")
print("0 <= startindex <= " + str(delta))
roi[index_big] = startindex, startindex + ishape[index_small]
square_image = image[roi[0, 0] : roi[0, 1], roi[1, 0] : roi[1, 1]]
return square_image
#%% image rotation
# this function is a modified version of the original from
# https://github.com/PyImageSearch/imutils/blob/master/imutils/convenience.py#L41
[docs]
def img_rotate_bound(image, angle, flag="cubic", bm=1):
"""
rotates an image by the given angle clockwise;
The rotated image is given in a rectangular bounding box
without cutting off parts of the original image.
Args:
image (MxN array_like):
np.uint8.
angle (float):
angle given in degrees.
flag (string, optional):
possibilities:"cubic","linear";
sets the method of interplation.
Defaults to "cubic".
bm (int, optional):
sets the border mode,
extrapolating from the borders of the image.
0: continues the image by padding zeros
1: continues the image by repeating the border-pixel values.
Defaults to 1.
Returns:
rotated_image (KxL array_like):
np.uint8.
log (list):
looking like [M,N,inverse_rotation_matrix],
contains the original shape M,N and the matrix needed
to invert the rotation for the function img_rotate_back.
"""
(h, w) = image.shape[:2]
(cX, cY) = (w / 2, h / 2)
# grab the rotation matrix (applying the negative of the
# angle to rotate clockwise), then grab the sine and cosine
# (i.e., the rotation components of the matrix)
M = cv2.getRotationMatrix2D((cX, cY), -angle, 1.0)
cos = np.abs(M[0, 0])
sin = np.abs(M[0, 1])
# compute the new bounding dimensions of the image
nW = int(np.round((h * sin) + (w * cos)))
nH = int(np.round((h * cos) + (w * sin)))
if bm == 0:
bm = cv2.BORDER_CONSTANT
elif bm == 1:
bm = cv2.BORDER_REPLICATE
# adjust the rotation matrix to take into account translation
M[0, 2] += (nW / 2) - cX
M[1, 2] += (nH / 2) - cY
invRotateMatrix = cv2.invertAffineTransform(M)
log = [(h, w), invRotateMatrix]
# perform the actual rotation and return the image
if flag == "cubic":
return (
cv2.warpAffine(image, M, (nW, nH), flags=cv2.INTER_CUBIC, borderMode=bm),
log,
)
else:
return (
cv2.warpAffine(image, M, (nW, nH), flags=cv2.INTER_LINEAR, borderMode=bm),
log,
)
[docs]
def img_rotate_back(image, log, flag="cubic", bm=1):
"""
invert the rotation done by img_rotate_bound returning the image
to its original shape MxN, cutting away padded values for the
bounding box generated by img_rotate_bound
Args:
image (KxL array_like):
np.uint8.
log (list):
[M,N,inverse_rotation_matrix],
contains the original shape M,N and the matrix needed
to invert the rotation. log is given by the function img_rotate_bound.
flag (string, optional):
possibilities:"cubic","linear";
sets the method of interplation.
Defaults to "cubic".
bm (int, optional):
possibilities: 0,1;
sets the border mode, extrapolating from the borders of the image.
0: continues the image by padding zeros
1: continues the image by repeating the border-pixel values.
(bm=1 allows more exact back transformation, avoiding the
decrease of border-pixel values due to interpolation with zeros.)
Defaults to 1.
Returns:
inverse_rotated_image (MxN array_like):
np.uint8.
"""
(h, w), invM = log
if bm == 0:
bm = cv2.BORDER_CONSTANT
elif bm == 1:
bm = cv2.BORDER_REPLICATE
if flag == "cubic":
return cv2.warpAffine(image, invM, (w, h), flags=cv2.INTER_CUBIC, borderMode=bm)
else:
return cv2.warpAffine(
image, invM, (w, h), flags=cv2.INTER_LINEAR, borderMode=bm
)
#%% tiling
[docs]
def img_periodic_tiling(img, tiles=3):
"""
takes an image as a tile and creates a tiling of
tiles x tiles by duplicating it.
Args:
img (MxN array_like):
2d-dataset / image.
tiles (int, optional):
number of tiles in vertical and horizontal direction.
the number of tiles must be uneven.
Defaults to 3.
Returns:
tiled (array_like):
with shape tiles*M x tiles*N.
orig (tuples):
containing the bounding coordinates of the center image
(lower_row_limit,upper_row_limit),(lower_column_limit,upper_column_limit).
"""
s = np.array(img.shape)
tiled = np.zeros(s * tiles, dtype=img.dtype)
for i in range(tiles):
for j in range(tiles):
tiled[s[0] * i : s[0] * (i + 1), s[1] * j : s[1] * (j + 1)] = img
oij = tiles // 2
orig = (s[0] * oij, s[0] * (oij + 1)), (s[1] * oij, s[1] * (oij + 1))
return tiled, orig
#%% special image transformations
#%% asymmetric non maximum supppression
@njit
def _aysmmetric_non_maximum_suppression(
newimg, img, cimg, rimg, mask, thresh_ratio, ksize, asympix, damping
):
ioffs = ksize // 2
joffs = ksize // 2 + asympix // 2
for i in range(ioffs, img.shape[0] - ioffs):
for j in range(joffs, img.shape[1] - joffs):
if not mask[i, j]:
pass
elif (
not mask[i - ioffs, j - joffs]
or not mask[i + ioffs, j + joffs]
or not mask[i - ioffs, j + joffs]
or not mask[i + ioffs, j - joffs]
):
pass
else:
v = max(cimg[i, j - joffs : j + joffs + 1])
h = max(rimg[i - ioffs : i + ioffs + 1, j]) * ksize / (ksize + asympix)
if h > v * thresh_ratio:
newimg[i, j] = img[i, j]
else:
newimg[i, j] = (
img[i, j] / damping
) # np.min(img[i-ioffs:i+ioffs+1,j-joffs:j+joffs+1])
return newimg
[docs]
def img_anms(img, mask, thresh_ratio=1.5, ksize=5, asympix=0, damping=5):
"""
asymmetric non maximum supppression
Args:
img (array_like):
DESCRIPTION.
mask (TYPE):
DESCRIPTION.
thresh_ratio (float, optional):
DESCRIPTION.
Defaults to 1.5.
ksize (int, optional):
uneven integer.
Defaults to 5.
asympix (TYPE, optional):
DESCRIPTION.
Defaults to 0.
damping (TYPE, optional):
DESCRIPTION.
Defaults to 5.
Returns:
processed_image (array_like):
DESCRIPTION.
"""
newimg = copy.deepcopy(img)
cimg = cv2.sepFilter2D(
img, cv2.CV_64F, np.ones(1), np.ones(ksize), borderType=cv2.BORDER_ISOLATED
)
rimg = cv2.sepFilter2D(
img,
cv2.CV_64F,
np.ones(ksize + asympix),
np.ones(1),
borderType=cv2.BORDER_ISOLATED,
)
return _aysmmetric_non_maximum_suppression(
newimg, img, cimg, rimg, mask, thresh_ratio, ksize, asympix, damping
)