656 lines
20 KiB
Python
656 lines
20 KiB
Python
# -*- coding: utf-8 -*-
|
|
import numpy as np
|
|
import scipy
|
|
from scipy import fftpack
|
|
import torch
|
|
|
|
from math import cos, sin
|
|
from numpy import zeros, ones, prod, array, pi, log, min, mod, arange, sum, mgrid, exp, pad, round
|
|
from numpy.random import randn, rand
|
|
from scipy.signal import convolve2d
|
|
import cv2
|
|
import random
|
|
# import utils_image as util
|
|
|
|
'''
|
|
modified by Kai Zhang (github: https://github.com/cszn)
|
|
03/03/2019
|
|
'''
|
|
|
|
|
|
def get_uperleft_denominator(img, kernel):
|
|
'''
|
|
img: HxWxC
|
|
kernel: hxw
|
|
denominator: HxWx1
|
|
upperleft: HxWxC
|
|
'''
|
|
V = psf2otf(kernel, img.shape[:2])
|
|
denominator = np.expand_dims(np.abs(V)**2, axis=2)
|
|
upperleft = np.expand_dims(np.conj(V), axis=2) * np.fft.fft2(img, axes=[0, 1])
|
|
return upperleft, denominator
|
|
|
|
|
|
def get_uperleft_denominator_pytorch(img, kernel):
|
|
'''
|
|
img: NxCxHxW
|
|
kernel: Nx1xhxw
|
|
denominator: Nx1xHxW
|
|
upperleft: NxCxHxWx2
|
|
'''
|
|
V = p2o(kernel, img.shape[-2:]) # Nx1xHxWx2
|
|
denominator = V[..., 0]**2+V[..., 1]**2 # Nx1xHxW
|
|
upperleft = cmul(cconj(V), rfft(img)) # Nx1xHxWx2 * NxCxHxWx2
|
|
return upperleft, denominator
|
|
|
|
|
|
def c2c(x):
|
|
return torch.from_numpy(np.stack([np.float32(x.real), np.float32(x.imag)], axis=-1))
|
|
|
|
|
|
def r2c(x):
|
|
return torch.stack([x, torch.zeros_like(x)], -1)
|
|
|
|
|
|
def cdiv(x, y):
|
|
a, b = x[..., 0], x[..., 1]
|
|
c, d = y[..., 0], y[..., 1]
|
|
cd2 = c**2 + d**2
|
|
return torch.stack([(a*c+b*d)/cd2, (b*c-a*d)/cd2], -1)
|
|
|
|
|
|
def cabs(x):
|
|
return torch.pow(x[..., 0]**2+x[..., 1]**2, 0.5)
|
|
|
|
|
|
def cmul(t1, t2):
|
|
'''
|
|
complex multiplication
|
|
t1: NxCxHxWx2
|
|
output: NxCxHxWx2
|
|
'''
|
|
real1, imag1 = t1[..., 0], t1[..., 1]
|
|
real2, imag2 = t2[..., 0], t2[..., 1]
|
|
return torch.stack([real1 * real2 - imag1 * imag2, real1 * imag2 + imag1 * real2], dim=-1)
|
|
|
|
|
|
def cconj(t, inplace=False):
|
|
'''
|
|
# complex's conjugation
|
|
t: NxCxHxWx2
|
|
output: NxCxHxWx2
|
|
'''
|
|
c = t.clone() if not inplace else t
|
|
c[..., 1] *= -1
|
|
return c
|
|
|
|
|
|
def rfft(t):
|
|
return torch.rfft(t, 2, onesided=False)
|
|
|
|
|
|
def irfft(t):
|
|
return torch.irfft(t, 2, onesided=False)
|
|
|
|
|
|
def fft(t):
|
|
return torch.fft(t, 2)
|
|
|
|
|
|
def ifft(t):
|
|
return torch.ifft(t, 2)
|
|
|
|
|
|
def p2o(psf, shape):
|
|
'''
|
|
# psf: NxCxhxw
|
|
# shape: [H,W]
|
|
# otf: NxCxHxWx2
|
|
'''
|
|
otf = torch.zeros(psf.shape[:-2] + shape).type_as(psf)
|
|
otf[...,:psf.shape[2],:psf.shape[3]].copy_(psf)
|
|
for axis, axis_size in enumerate(psf.shape[2:]):
|
|
otf = torch.roll(otf, -int(axis_size / 2), dims=axis+2)
|
|
otf = torch.rfft(otf, 2, onesided=False)
|
|
n_ops = torch.sum(torch.tensor(psf.shape).type_as(psf) * torch.log2(torch.tensor(psf.shape).type_as(psf)))
|
|
otf[...,1][torch.abs(otf[...,1])<n_ops*2.22e-16] = torch.tensor(0).type_as(psf)
|
|
return otf
|
|
|
|
|
|
|
|
# otf2psf: not sure where I got this one from. Maybe translated from Octave source code or whatever. It's just math.
|
|
def otf2psf(otf, outsize=None):
|
|
insize = np.array(otf.shape)
|
|
psf = np.fft.ifftn(otf, axes=(0, 1))
|
|
for axis, axis_size in enumerate(insize):
|
|
psf = np.roll(psf, np.floor(axis_size / 2).astype(int), axis=axis)
|
|
if type(outsize) != type(None):
|
|
insize = np.array(otf.shape)
|
|
outsize = np.array(outsize)
|
|
n = max(np.size(outsize), np.size(insize))
|
|
# outsize = postpad(outsize(:), n, 1);
|
|
# insize = postpad(insize(:) , n, 1);
|
|
colvec_out = outsize.flatten().reshape((np.size(outsize), 1))
|
|
colvec_in = insize.flatten().reshape((np.size(insize), 1))
|
|
outsize = np.pad(colvec_out, ((0, max(0, n - np.size(colvec_out))), (0, 0)), mode="constant")
|
|
insize = np.pad(colvec_in, ((0, max(0, n - np.size(colvec_in))), (0, 0)), mode="constant")
|
|
|
|
pad = (insize - outsize) / 2
|
|
if np.any(pad < 0):
|
|
print("otf2psf error: OUTSIZE must be smaller than or equal than OTF size")
|
|
prepad = np.floor(pad)
|
|
postpad = np.ceil(pad)
|
|
dims_start = prepad.astype(int)
|
|
dims_end = (insize - postpad).astype(int)
|
|
for i in range(len(dims_start.shape)):
|
|
psf = np.take(psf, range(dims_start[i][0], dims_end[i][0]), axis=i)
|
|
n_ops = np.sum(otf.size * np.log2(otf.shape))
|
|
psf = np.real_if_close(psf, tol=n_ops)
|
|
return psf
|
|
|
|
|
|
# psf2otf copied/modified from https://github.com/aboucaud/pypher/blob/master/pypher/pypher.py
|
|
def psf2otf(psf, shape=None):
|
|
"""
|
|
Convert point-spread function to optical transfer function.
|
|
Compute the Fast Fourier Transform (FFT) of the point-spread
|
|
function (PSF) array and creates the optical transfer function (OTF)
|
|
array that is not influenced by the PSF off-centering.
|
|
By default, the OTF array is the same size as the PSF array.
|
|
To ensure that the OTF is not altered due to PSF off-centering, PSF2OTF
|
|
post-pads the PSF array (down or to the right) with zeros to match
|
|
dimensions specified in OUTSIZE, then circularly shifts the values of
|
|
the PSF array up (or to the left) until the central pixel reaches (1,1)
|
|
position.
|
|
Parameters
|
|
----------
|
|
psf : `numpy.ndarray`
|
|
PSF array
|
|
shape : int
|
|
Output shape of the OTF array
|
|
Returns
|
|
-------
|
|
otf : `numpy.ndarray`
|
|
OTF array
|
|
Notes
|
|
-----
|
|
Adapted from MATLAB psf2otf function
|
|
"""
|
|
if type(shape) == type(None):
|
|
shape = psf.shape
|
|
shape = np.array(shape)
|
|
if np.all(psf == 0):
|
|
# return np.zeros_like(psf)
|
|
return np.zeros(shape)
|
|
if len(psf.shape) == 1:
|
|
psf = psf.reshape((1, psf.shape[0]))
|
|
inshape = psf.shape
|
|
psf = zero_pad(psf, shape, position='corner')
|
|
for axis, axis_size in enumerate(inshape):
|
|
psf = np.roll(psf, -int(axis_size / 2), axis=axis)
|
|
# Compute the OTF
|
|
otf = np.fft.fft2(psf, axes=(0, 1))
|
|
# Estimate the rough number of operations involved in the FFT
|
|
# and discard the PSF imaginary part if within roundoff error
|
|
# roundoff error = machine epsilon = sys.float_info.epsilon
|
|
# or np.finfo().eps
|
|
n_ops = np.sum(psf.size * np.log2(psf.shape))
|
|
otf = np.real_if_close(otf, tol=n_ops)
|
|
return otf
|
|
|
|
|
|
def zero_pad(image, shape, position='corner'):
|
|
"""
|
|
Extends image to a certain size with zeros
|
|
Parameters
|
|
----------
|
|
image: real 2d `numpy.ndarray`
|
|
Input image
|
|
shape: tuple of int
|
|
Desired output shape of the image
|
|
position : str, optional
|
|
The position of the input image in the output one:
|
|
* 'corner'
|
|
top-left corner (default)
|
|
* 'center'
|
|
centered
|
|
Returns
|
|
-------
|
|
padded_img: real `numpy.ndarray`
|
|
The zero-padded image
|
|
"""
|
|
shape = np.asarray(shape, dtype=int)
|
|
imshape = np.asarray(image.shape, dtype=int)
|
|
if np.alltrue(imshape == shape):
|
|
return image
|
|
if np.any(shape <= 0):
|
|
raise ValueError("ZERO_PAD: null or negative shape given")
|
|
dshape = shape - imshape
|
|
if np.any(dshape < 0):
|
|
raise ValueError("ZERO_PAD: target size smaller than source one")
|
|
pad_img = np.zeros(shape, dtype=image.dtype)
|
|
idx, idy = np.indices(imshape)
|
|
if position == 'center':
|
|
if np.any(dshape % 2 != 0):
|
|
raise ValueError("ZERO_PAD: source and target shapes "
|
|
"have different parity.")
|
|
offx, offy = dshape // 2
|
|
else:
|
|
offx, offy = (0, 0)
|
|
pad_img[idx + offx, idy + offy] = image
|
|
return pad_img
|
|
|
|
|
|
'''
|
|
Reducing boundary artifacts
|
|
'''
|
|
|
|
|
|
def opt_fft_size(n):
|
|
'''
|
|
Kai Zhang (github: https://github.com/cszn)
|
|
03/03/2019
|
|
# opt_fft_size.m
|
|
# compute an optimal data length for Fourier transforms
|
|
# written by Sunghyun Cho (sodomau@postech.ac.kr)
|
|
# persistent opt_fft_size_LUT;
|
|
'''
|
|
|
|
LUT_size = 2048
|
|
# print("generate opt_fft_size_LUT")
|
|
opt_fft_size_LUT = np.zeros(LUT_size)
|
|
|
|
e2 = 1
|
|
while e2 <= LUT_size:
|
|
e3 = e2
|
|
while e3 <= LUT_size:
|
|
e5 = e3
|
|
while e5 <= LUT_size:
|
|
e7 = e5
|
|
while e7 <= LUT_size:
|
|
if e7 <= LUT_size:
|
|
opt_fft_size_LUT[e7-1] = e7
|
|
if e7*11 <= LUT_size:
|
|
opt_fft_size_LUT[e7*11-1] = e7*11
|
|
if e7*13 <= LUT_size:
|
|
opt_fft_size_LUT[e7*13-1] = e7*13
|
|
e7 = e7 * 7
|
|
e5 = e5 * 5
|
|
e3 = e3 * 3
|
|
e2 = e2 * 2
|
|
|
|
nn = 0
|
|
for i in range(LUT_size, 0, -1):
|
|
if opt_fft_size_LUT[i-1] != 0:
|
|
nn = i-1
|
|
else:
|
|
opt_fft_size_LUT[i-1] = nn+1
|
|
|
|
m = np.zeros(len(n))
|
|
for c in range(len(n)):
|
|
nn = n[c]
|
|
if nn <= LUT_size:
|
|
m[c] = opt_fft_size_LUT[nn-1]
|
|
else:
|
|
m[c] = -1
|
|
return m
|
|
|
|
|
|
def wrap_boundary_liu(img, img_size):
|
|
|
|
"""
|
|
Reducing boundary artifacts in image deconvolution
|
|
Renting Liu, Jiaya Jia
|
|
ICIP 2008
|
|
"""
|
|
if img.ndim == 2:
|
|
ret = wrap_boundary(img, img_size)
|
|
elif img.ndim == 3:
|
|
ret = [wrap_boundary(img[:, :, i], img_size) for i in range(3)]
|
|
ret = np.stack(ret, 2)
|
|
return ret
|
|
|
|
|
|
def wrap_boundary(img, img_size):
|
|
|
|
"""
|
|
python code from:
|
|
https://github.com/ys-koshelev/nla_deblur/blob/90fe0ab98c26c791dcbdf231fe6f938fca80e2a0/boundaries.py
|
|
Reducing boundary artifacts in image deconvolution
|
|
Renting Liu, Jiaya Jia
|
|
ICIP 2008
|
|
"""
|
|
(H, W) = np.shape(img)
|
|
H_w = int(img_size[0]) - H
|
|
W_w = int(img_size[1]) - W
|
|
|
|
# ret = np.zeros((img_size[0], img_size[1]));
|
|
alpha = 1
|
|
HG = img[:, :]
|
|
|
|
r_A = np.zeros((alpha*2+H_w, W))
|
|
r_A[:alpha, :] = HG[-alpha:, :]
|
|
r_A[-alpha:, :] = HG[:alpha, :]
|
|
a = np.arange(H_w)/(H_w-1)
|
|
# r_A(alpha+1:end-alpha, 1) = (1-a)*r_A(alpha,1) + a*r_A(end-alpha+1,1)
|
|
r_A[alpha:-alpha, 0] = (1-a)*r_A[alpha-1, 0] + a*r_A[-alpha, 0]
|
|
# r_A(alpha+1:end-alpha, end) = (1-a)*r_A(alpha,end) + a*r_A(end-alpha+1,end)
|
|
r_A[alpha:-alpha, -1] = (1-a)*r_A[alpha-1, -1] + a*r_A[-alpha, -1]
|
|
|
|
r_B = np.zeros((H, alpha*2+W_w))
|
|
r_B[:, :alpha] = HG[:, -alpha:]
|
|
r_B[:, -alpha:] = HG[:, :alpha]
|
|
a = np.arange(W_w)/(W_w-1)
|
|
r_B[0, alpha:-alpha] = (1-a)*r_B[0, alpha-1] + a*r_B[0, -alpha]
|
|
r_B[-1, alpha:-alpha] = (1-a)*r_B[-1, alpha-1] + a*r_B[-1, -alpha]
|
|
|
|
if alpha == 1:
|
|
A2 = solve_min_laplacian(r_A[alpha-1:, :])
|
|
B2 = solve_min_laplacian(r_B[:, alpha-1:])
|
|
r_A[alpha-1:, :] = A2
|
|
r_B[:, alpha-1:] = B2
|
|
else:
|
|
A2 = solve_min_laplacian(r_A[alpha-1:-alpha+1, :])
|
|
r_A[alpha-1:-alpha+1, :] = A2
|
|
B2 = solve_min_laplacian(r_B[:, alpha-1:-alpha+1])
|
|
r_B[:, alpha-1:-alpha+1] = B2
|
|
A = r_A
|
|
B = r_B
|
|
|
|
r_C = np.zeros((alpha*2+H_w, alpha*2+W_w))
|
|
r_C[:alpha, :] = B[-alpha:, :]
|
|
r_C[-alpha:, :] = B[:alpha, :]
|
|
r_C[:, :alpha] = A[:, -alpha:]
|
|
r_C[:, -alpha:] = A[:, :alpha]
|
|
|
|
if alpha == 1:
|
|
C2 = C2 = solve_min_laplacian(r_C[alpha-1:, alpha-1:])
|
|
r_C[alpha-1:, alpha-1:] = C2
|
|
else:
|
|
C2 = solve_min_laplacian(r_C[alpha-1:-alpha+1, alpha-1:-alpha+1])
|
|
r_C[alpha-1:-alpha+1, alpha-1:-alpha+1] = C2
|
|
C = r_C
|
|
# return C
|
|
A = A[alpha-1:-alpha-1, :]
|
|
B = B[:, alpha:-alpha]
|
|
C = C[alpha:-alpha, alpha:-alpha]
|
|
ret = np.vstack((np.hstack((img, B)), np.hstack((A, C))))
|
|
return ret
|
|
|
|
|
|
def solve_min_laplacian(boundary_image):
|
|
(H, W) = np.shape(boundary_image)
|
|
|
|
# Laplacian
|
|
f = np.zeros((H, W))
|
|
# boundary image contains image intensities at boundaries
|
|
boundary_image[1:-1, 1:-1] = 0
|
|
j = np.arange(2, H)-1
|
|
k = np.arange(2, W)-1
|
|
f_bp = np.zeros((H, W))
|
|
f_bp[np.ix_(j, k)] = -4*boundary_image[np.ix_(j, k)] + boundary_image[np.ix_(j, k+1)] + boundary_image[np.ix_(j, k-1)] + boundary_image[np.ix_(j-1, k)] + boundary_image[np.ix_(j+1, k)]
|
|
|
|
del(j, k)
|
|
f1 = f - f_bp # subtract boundary points contribution
|
|
del(f_bp, f)
|
|
|
|
# DST Sine Transform algo starts here
|
|
f2 = f1[1:-1,1:-1]
|
|
del(f1)
|
|
|
|
# compute sine tranform
|
|
if f2.shape[1] == 1:
|
|
tt = fftpack.dst(f2, type=1, axis=0)/2
|
|
else:
|
|
tt = fftpack.dst(f2, type=1)/2
|
|
|
|
if tt.shape[0] == 1:
|
|
f2sin = np.transpose(fftpack.dst(np.transpose(tt), type=1, axis=0)/2)
|
|
else:
|
|
f2sin = np.transpose(fftpack.dst(np.transpose(tt), type=1)/2)
|
|
del(f2)
|
|
|
|
# compute Eigen Values
|
|
[x, y] = np.meshgrid(np.arange(1, W-1), np.arange(1, H-1))
|
|
denom = (2*np.cos(np.pi*x/(W-1))-2) + (2*np.cos(np.pi*y/(H-1)) - 2)
|
|
|
|
# divide
|
|
f3 = f2sin/denom
|
|
del(f2sin, x, y)
|
|
|
|
# compute Inverse Sine Transform
|
|
if f3.shape[0] == 1:
|
|
tt = fftpack.idst(f3*2, type=1, axis=1)/(2*(f3.shape[1]+1))
|
|
else:
|
|
tt = fftpack.idst(f3*2, type=1, axis=0)/(2*(f3.shape[0]+1))
|
|
del(f3)
|
|
if tt.shape[1] == 1:
|
|
img_tt = np.transpose(fftpack.idst(np.transpose(tt)*2, type=1)/(2*(tt.shape[0]+1)))
|
|
else:
|
|
img_tt = np.transpose(fftpack.idst(np.transpose(tt)*2, type=1, axis=0)/(2*(tt.shape[1]+1)))
|
|
del(tt)
|
|
|
|
# put solution in inner points; outer points obtained from boundary image
|
|
img_direct = boundary_image
|
|
img_direct[1:-1, 1:-1] = 0
|
|
img_direct[1:-1, 1:-1] = img_tt
|
|
return img_direct
|
|
|
|
|
|
"""
|
|
Created on Thu Jan 18 15:36:32 2018
|
|
@author: italo
|
|
https://github.com/ronaldosena/imagens-medicas-2/blob/40171a6c259edec7827a6693a93955de2bd39e76/Aulas/aula_2_-_uniform_filter/matlab_fspecial.py
|
|
"""
|
|
|
|
"""
|
|
Syntax
|
|
h = fspecial(type)
|
|
h = fspecial('average',hsize)
|
|
h = fspecial('disk',radius)
|
|
h = fspecial('gaussian',hsize,sigma)
|
|
h = fspecial('laplacian',alpha)
|
|
h = fspecial('log',hsize,sigma)
|
|
h = fspecial('motion',len,theta)
|
|
h = fspecial('prewitt')
|
|
h = fspecial('sobel')
|
|
"""
|
|
|
|
|
|
def fspecial_average(hsize=3):
|
|
"""Smoothing filter"""
|
|
return np.ones((hsize, hsize))/hsize**2
|
|
|
|
|
|
def fspecial_disk(radius):
|
|
"""Disk filter"""
|
|
raise(NotImplemented)
|
|
rad = 0.6
|
|
crad = np.ceil(rad-0.5)
|
|
[x, y] = np.meshgrid(np.arange(-crad, crad+1), np.arange(-crad, crad+1))
|
|
maxxy = np.zeros(x.shape)
|
|
maxxy[abs(x) >= abs(y)] = abs(x)[abs(x) >= abs(y)]
|
|
maxxy[abs(y) >= abs(x)] = abs(y)[abs(y) >= abs(x)]
|
|
minxy = np.zeros(x.shape)
|
|
minxy[abs(x) <= abs(y)] = abs(x)[abs(x) <= abs(y)]
|
|
minxy[abs(y) <= abs(x)] = abs(y)[abs(y) <= abs(x)]
|
|
m1 = (rad**2 < (maxxy+0.5)**2 + (minxy-0.5)**2)*(minxy-0.5) +\
|
|
(rad**2 >= (maxxy+0.5)**2 + (minxy-0.5)**2)*\
|
|
np.sqrt((rad**2 + 0j) - (maxxy + 0.5)**2)
|
|
m2 = (rad**2 > (maxxy-0.5)**2 + (minxy+0.5)**2)*(minxy+0.5) +\
|
|
(rad**2 <= (maxxy-0.5)**2 + (minxy+0.5)**2)*\
|
|
np.sqrt((rad**2 + 0j) - (maxxy - 0.5)**2)
|
|
h = None
|
|
return h
|
|
|
|
|
|
def fspecial_gaussian(hsize, sigma):
|
|
hsize = [hsize, hsize]
|
|
siz = [(hsize[0]-1.0)/2.0, (hsize[1]-1.0)/2.0]
|
|
std = sigma
|
|
[x, y] = np.meshgrid(np.arange(-siz[1], siz[1]+1), np.arange(-siz[0], siz[0]+1))
|
|
arg = -(x*x + y*y)/(2*std*std)
|
|
h = np.exp(arg)
|
|
h[h < scipy.finfo(float).eps * h.max()] = 0
|
|
sumh = h.sum()
|
|
if sumh != 0:
|
|
h = h/sumh
|
|
return h
|
|
|
|
|
|
def fspecial_laplacian(alpha):
|
|
alpha = max([0, min([alpha,1])])
|
|
h1 = alpha/(alpha+1)
|
|
h2 = (1-alpha)/(alpha+1)
|
|
h = [[h1, h2, h1], [h2, -4/(alpha+1), h2], [h1, h2, h1]]
|
|
h = np.array(h)
|
|
return h
|
|
|
|
|
|
def fspecial_log(hsize, sigma):
|
|
raise(NotImplemented)
|
|
|
|
|
|
def fspecial_motion(motion_len, theta):
|
|
raise(NotImplemented)
|
|
|
|
|
|
def fspecial_prewitt():
|
|
return np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]])
|
|
|
|
|
|
def fspecial_sobel():
|
|
return np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]])
|
|
|
|
|
|
def fspecial(filter_type, *args, **kwargs):
|
|
'''
|
|
python code from:
|
|
https://github.com/ronaldosena/imagens-medicas-2/blob/40171a6c259edec7827a6693a93955de2bd39e76/Aulas/aula_2_-_uniform_filter/matlab_fspecial.py
|
|
'''
|
|
if filter_type == 'average':
|
|
return fspecial_average(*args, **kwargs)
|
|
if filter_type == 'disk':
|
|
return fspecial_disk(*args, **kwargs)
|
|
if filter_type == 'gaussian':
|
|
return fspecial_gaussian(*args, **kwargs)
|
|
if filter_type == 'laplacian':
|
|
return fspecial_laplacian(*args, **kwargs)
|
|
if filter_type == 'log':
|
|
return fspecial_log(*args, **kwargs)
|
|
if filter_type == 'motion':
|
|
return fspecial_motion(*args, **kwargs)
|
|
if filter_type == 'prewitt':
|
|
return fspecial_prewitt(*args, **kwargs)
|
|
if filter_type == 'sobel':
|
|
return fspecial_sobel(*args, **kwargs)
|
|
|
|
|
|
def fspecial_gauss(size, sigma):
|
|
x, y = mgrid[-size // 2 + 1 : size // 2 + 1, -size // 2 + 1 : size // 2 + 1]
|
|
g = exp(-((x ** 2 + y ** 2) / (2.0 * sigma ** 2)))
|
|
return g / g.sum()
|
|
|
|
|
|
def blurkernel_synthesis(h=37, w=None):
|
|
# https://github.com/tkkcc/prior/blob/879a0b6c117c810776d8cc6b63720bf29f7d0cc4/util/gen_kernel.py
|
|
w = h if w is None else w
|
|
kdims = [h, w]
|
|
x = randomTrajectory(250)
|
|
k = None
|
|
while k is None:
|
|
k = kernelFromTrajectory(x)
|
|
|
|
# center pad to kdims
|
|
pad_width = ((kdims[0] - k.shape[0]) // 2, (kdims[1] - k.shape[1]) // 2)
|
|
pad_width = [(pad_width[0],), (pad_width[1],)]
|
|
|
|
if pad_width[0][0]<0 or pad_width[1][0]<0:
|
|
k = k[0:h, 0:h]
|
|
else:
|
|
k = pad(k, pad_width, "constant")
|
|
x1,x2 = k.shape
|
|
if np.random.randint(0, 4) == 1:
|
|
k = cv2.resize(k, (random.randint(x1, 5*x1), random.randint(x2, 5*x2)), interpolation=cv2.INTER_LINEAR)
|
|
y1, y2 = k.shape
|
|
k = k[(y1-x1)//2: (y1-x1)//2+x1, (y2-x2)//2: (y2-x2)//2+x2]
|
|
|
|
if sum(k)<0.1:
|
|
k = fspecial_gaussian(h, 0.1+6*np.random.rand(1))
|
|
k = k / sum(k)
|
|
# import matplotlib.pyplot as plt
|
|
# plt.imshow(k, interpolation="nearest", cmap="gray")
|
|
# plt.show()
|
|
return k
|
|
|
|
|
|
def kernelFromTrajectory(x):
|
|
h = 5 - log(rand()) / 0.15
|
|
h = round(min([h, 27])).astype(int)
|
|
h = h + 1 - h % 2
|
|
w = h
|
|
k = zeros((h, w))
|
|
|
|
xmin = min(x[0])
|
|
xmax = max(x[0])
|
|
ymin = min(x[1])
|
|
ymax = max(x[1])
|
|
xthr = arange(xmin, xmax, (xmax - xmin) / w)
|
|
ythr = arange(ymin, ymax, (ymax - ymin) / h)
|
|
|
|
for i in range(1, xthr.size):
|
|
for j in range(1, ythr.size):
|
|
idx = (
|
|
(x[0, :] >= xthr[i - 1])
|
|
& (x[0, :] < xthr[i])
|
|
& (x[1, :] >= ythr[j - 1])
|
|
& (x[1, :] < ythr[j])
|
|
)
|
|
k[i - 1, j - 1] = sum(idx)
|
|
if sum(k) == 0:
|
|
return
|
|
k = k / sum(k)
|
|
k = convolve2d(k, fspecial_gauss(3, 1), "same")
|
|
k = k / sum(k)
|
|
return k
|
|
|
|
|
|
def randomTrajectory(T):
|
|
x = zeros((3, T))
|
|
v = randn(3, T)
|
|
r = zeros((3, T))
|
|
trv = 1 / 1
|
|
trr = 2 * pi / T
|
|
for t in range(1, T):
|
|
F_rot = randn(3) / (t + 1) + r[:, t - 1]
|
|
F_trans = randn(3) / (t + 1)
|
|
r[:, t] = r[:, t - 1] + trr * F_rot
|
|
v[:, t] = v[:, t - 1] + trv * F_trans
|
|
st = v[:, t]
|
|
st = rot3D(st, r[:, t])
|
|
x[:, t] = x[:, t - 1] + st
|
|
return x
|
|
|
|
|
|
def rot3D(x, r):
|
|
Rx = array([[1, 0, 0], [0, cos(r[0]), -sin(r[0])], [0, sin(r[0]), cos(r[0])]])
|
|
Ry = array([[cos(r[1]), 0, sin(r[1])], [0, 1, 0], [-sin(r[1]), 0, cos(r[1])]])
|
|
Rz = array([[cos(r[2]), -sin(r[2]), 0], [sin(r[2]), cos(r[2]), 0], [0, 0, 1]])
|
|
R = Rz @ Ry @ Rx
|
|
x = R @ x
|
|
return x
|
|
|
|
|
|
if __name__ == '__main__':
|
|
a = opt_fft_size([111])
|
|
print(a)
|
|
|
|
print(fspecial('gaussian', 5, 1))
|
|
|
|
print(p2o(torch.zeros(1,1,4,4).float(),(14,14)).shape)
|
|
|
|
k = blurkernel_synthesis(11)
|
|
import matplotlib.pyplot as plt
|
|
plt.imshow(k, interpolation="nearest", cmap="gray")
|
|
plt.show()
|