摘要：import numpy as npimport matplotlib.pyplot as pltfrom scipy.signal import convolve2dfrom scipy.ndimage import binary_closing, bina
import numpy as npimport matplotlib.pyplot as pltfrom scipy.signal import convolve2dfrom scipy.ndimage import binary_closing, binary_dilation, binary_erosion, binary_openingfrom numpy.fft import fft2, fftshift, ifft2, ifftshift, rfft2, irfft2import matplotlib.patches as patchesimport skimage.restoration as restorefrom skimage.metrics import peak_signal_noise_ratio, structural_similarityimport cv2 from scipy.ndimage import map_coordinatesfrom IPython import displayimport timeimport osimport tifffile as tiffdef drawnow(fig):display.display(fig)display.clear_output(wait=True)time.sleep(.01)img = plt.imread("synthetic/random_1px.png")if len(img.shape) > 2:print("Converting image to grayscale")img = img.mean(axis=2)plt.imshow(img, cmap="gray")plt.title("Input Image")plt.show
def polar_transform(img, center, output=None):r = (np.array(img.shape[:2])**2).sum**0.5/2if output is None:shp = int(round(r)), int(round(r*2*np.pi))output = np.zeros(shp, dtype=img.dtype)elif isinstance(output, tuple):output = np.zeros(output, dtype=img.dtype)out_h, out_w = output.shaper_samples = np.linspace(0, r, out_h)theta_samples = np.linspace(0, np.pi*2, out_w)xs = r_samples[:,None] * np.cos(theta_samples) + center[1]ys = r_samples[:,None] * np.sin(theta_samples) + center[0]map_coordinates(img, (ys, xs), order=1, output=output)return outputdef inverse_polar_transform(img, center, output=None):r = img.shape[0]if output is None:output = np.zeros((r*2, r*2), dtype=img.dtype)elif isinstance(output, tuple):output = np.zeros(output, dtype=img.dtype)out_h, out_w = output.shapey_samples, x_samples = np.mgrid[:out_h, :out_w]y_samples -= center[0]x_samples -= center[1]rs = (y_samples * y_samples + x_samples * x_samples) **0.5ts = np.arccos(x_samples / (rs + 1e-8))ts[y_samplesdef generate_star_trails(img, kernel):return convolve2d(img, kernel, mode='same')def contrast_stetch(img):max = img.maxmin = img.minreturn (img - min) / (max - min)def plot_polar_img(img, mode=0, title=''):if mode == 0:plt.imshow(img, cmap='gray', aspect='auto')else:plt.imshow(img, aspect='auto')plt.title(title)def save_images(recon_images: dict, path: str):if not os.path.exists(path):os.makedirs(path)for key, val in recon_images.items:val = val.astype(np.float32)tiff.imsave(os.path.join(path, key) + '.tif', val)def compute_stats(true_img: np.ndarray, recon_images: dict, path: str):file = open(os.path.join(path, "stats.txt"), 'w')file.write("Method\tPSNR\tSIMM\n")true_img = true_img.astype(np.float32)for key, val in recon_images.items:val = val.astype(np.float32)psnr = peak_signal_noise_ratio(true_img, val)simm = structural_similarity(true_img, val, win_size=11, data_range=1.0)file.write(f"{key}\t{psnr:.05}\t{simm:.05}\n")file.closedef find_center(img, threshold, ax=None):# Binarize Imageimg = contrast_stetch(img)bin_img = (img > threshold).astype('int')H, W = img.shape[:2]factor = 1if (H > 1000 or W > 1000):factor = max(H, W) // 1000# Dilating and subsamplingbin_img = binary_dilation(bin_img, np.ones((factor, factor)))bin_img = bin_img[::factor, ::factor].astype(np.uint8) * 255# To remove noiseblur = cv2.GaussianBlur(bin_img, ksize=(31, 31), sigmaX=1)circles = cv2.HoughCircles(blur, cv2.HOUGH_GRADIENT, 1, minDist=200, param1=50, param2=30, minRadius=10, maxRadius=0)circles = np.around(circles).reshape((-1, 3))# Create a figure. Equal aspect so circles look circularif ax is None:fig, ax = plt.subplots(1)ax.set_aspect('equal')# Show the imageax.imshow(img, cmap="gray")for circ in circles[:1]:ax.add_patch(patches.Circle(circ[:2] * factor, circ[2], edgecolor='r', linewidth=1, fill=False))ax.add_patch(patches.Rectangle(circ[:2] * factor, 20*factor, 20*factor))ax.set_title("Hough Transform")# # Show the image# # plt.showcirc = circles[0]circ[1::-1] = circ[:2]print("Predicted Center: ", circ[:2] * factor)return circ[:2] * factordef estimate_deblurring_kernel(img, time):'''time: (in minutes) the duration of exposure'''H, W = img.shape[:2]r = (np.array(img.shape[:2])**2).sum**0.5/2total_theta = int(round(2 * np.pi * r))estimated_theta = total_theta * (time / (24 * 60))estimated_theta = np.around(estimated_theta).astype('int')estimated_kernel = np.zeros((1, estimated_theta)) # Is it correct though?estimated_kernel[0, :estimated_theta] = 1 / estimated_thetareturn estimated_kerneldef wiener_deconvolution(image, psf, nsr=0.01):"""image (np.array) : MxNpsf (np.array) : MxN nsr (noise to signal ratio - 1/SNR from above derivation) : float"""# YOUR CODE HEREdeblurred_image = np.zeros_like(image)Y = fftshift(fft2(image))H = fftshift(fft2(psf))X = Y * np.conj(H) / (nsr + np.conj(H) * H)deblurred_image = ifftshift(ifft2(X))deblurred_image = np.abs(deblurred_image)return deblurred_imagedef F(x):return rfft2(x,axes=(0,1))def Fh(x):return irfft2(x,axes=(0,1))def crop(x):M,N = x.shaperl = M//4ru = 3*np.ceil(M/4)cl = N//4cu = 4*np.ceil(N/4)return x[M//4:3*M//4, N//4:3*N//4]def pad(psf, x):M, N = x.shape H, W = psf.shapepad_psf = np.zeros((2*H, 2*W))pad_psf[M//2:M//2+psf.shape[0], N//2:N//2+psf.shape[1]] = psfreturn pad_psfdef A(x, pad_psf):H_bar = F(ifftshift(pad_psf))H = np.real(Fh(H_bar * F(x)))return crop(H) def Ah(y, pad_psf):H = F(ifftshift(pad_psf))H_star = np.conj(H)return np.real(Fh(H_star*F(pad(y, y)))) Testing the implemented Functions.img = np.random.randn(300, 400)print(img.shape)kernel = np.ones((1, 100)) / 100psf = np.zeros_like(img)kh, kw = kernel.shapeH, W = psf.shapepsf[H//2-kh//2:H//2-kh//2+kh, W//2-kw//2:W//2-kw//2+kw] = kernel[M,N] = img.shapex = np.random.rand(2*M, 2*N)y = np.random.rand(M,N)pad_psf = pad(psf, x)Ax_y = np.sum(A(x, pad_psf)*y)x_Ahy = np.sum(x*Ah(y, pad_psf))## They should have the same value.print(' =',Ax_y)print('=',x_Ahy)(300, 400) = 29904.650702617542= 29904.65070261754def vanilla_gd(meas, pad_psf, n_iters=500):M, N = meas.shapex_k = np.zeros((2*M, 2*N)) #Initialize with zerosmu = 1 #Step sizefig = plt.figure #Initialize figureim_obj = plt.imshow(x_k, cmap='gray', aspect='auto') title_obj = plt.title('Gradient Descent, Iteration {}'.format(0))for k in range(n_iters):x_k = x_k - mu * Ah(A(x_k, pad_psf) - meas, pad_psf)## Display progressif k == 0:im_obj.set_data(np.maximum(crop(x_k)/np.max(x_k),0))plt.title('Gradient Descent, Iteration {}'.format(k))drawnow(fig)x_gd = x_kim_obj.set_data(np.maximum(crop(x_k)/np.max(x_k),0))plt.title('Gradient Descent, Iteration {}'.format(k))plt.showreturn x_gddef gd_fista(meas, pad_psf, n_iters=500, mu=1):M, N = meas.shape# Initialize variablest_k = 1 #Momentum parameterx_k = np.zeros((2*M, 2*N)) # Imagey_k = np.zeros((2*M, 2*N))z_k = np.zeros((2*M, 2*N))# Initialize plottingfig = plt.figureim_obj = plt.imshow(x_k, cmap='gray', aspect='auto') title_obj = plt.title('FISTA, Iteration {}'.format(0))for k in range(n_iters):# Update momentum parameterst_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2beta = (t_k - 1) / t_k1y_k1 = x_k - mu * Ah(A(x_k, pad_psf) - meas, pad_psf)z_k1 = np.maximum(0, y_k1)x_k = z_k1 + beta * (z_k1 - z_k)t_k = t_k1z_k = z_k1#print("Beta:", beta, " t_k: ", t_k)if k == 0:im_obj.set_data((np.maximum(crop(z_k)/np.max(z_k),0)))plt.title('FISTA, Iteration {}'.format(k))drawnow(fig)x_fista = z_k im_obj.set_data(np.maximum(crop(z_k)/np.max(z_k),0))plt.title('FISTA, Iteration {}'.format(k))plt.showreturn x_fistadef gd_fista_tikhonov_reg(meas, pad_psf, n_iters=500, mu=1, lamda=0.05):M, N = meas.shape## Function to compute gradient for the Tikhonov Regulariation Term.kernel_tv = np.asarray([[0, -1, 0], [-1, 4, -1], [0, -1 , 0]])# raise NotImplementedError# Initialize variablest_k = 1 #Momentum parameterx_k = np.zeros((2*M, 2*N)) # Imagey_k = np.zeros((2*M, 2*N))z_k = np.zeros((2*M, 2*N))mu = 1 #Step size# lamda = 0.00005 # Step Size for regularizerdelta_x = np.zeros_like(x_k)# Initialize plottingfig = plt.figureim_obj = plt.imshow(x_k, cmap='gray', aspect='auto') title_obj = plt.title('Fista + Tikhonov Regulariation, Iteration {}'.format(0))for k in range(n_iters):# Update momentum parameterst_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2beta = (t_k - 1) / t_k1# for i in range(3):# delta_x[:,:,i] = convolve2d(x_k[:,:,i], kernel, mode='same')delta_x = convolve2d(x_k, kernel_tv, mode='same')y_k1 = x_k - mu * (Ah(A(x_k, pad_psf) - meas, pad_psf) + 2 * lamda * delta_x)z_k1 = np.maximum(0, y_k1)x_k = z_k1 + beta * (z_k1 - z_k)t_k = t_k1z_k = z_k1if k == 0:im_obj.set_data((np.maximum(crop(z_k)/np.max(z_k),0)))plt.title('Fista + Tikhonov Regulariation, Iteration {}'.format(k))drawnow(fig)x_fista_reg = z_k im_obj.set_data(np.maximum(crop(z_k)/np.max(z_k),0))plt.title('Fista + Tikhonov Regulariation, Iteration {}'.format(k))plt.showreturn x_fista_regdef gd_fista_total_variation_reg(meas, pad_psf, n_iters=500, mu=1, lamda=0.05):def gradient_tv(x):grad = np.zeros_like(x)grad[:, :-1] += x[:, :-1] - x[:, 1:]grad[:, 1:] += x[:, 1:] - x[:, :-1]grad[:-1, :] += x[:-1, :] - x[1:, :]grad[1:, :] += x[1:, :] - x[:-1, :]return gradM, N = meas.shape# Initialize variablest_k = 1 #Momentum parameterx_k = np.zeros((2*M, 2*N)) # Imagey_k = np.zeros((2*M, 2*N))z_k = np.zeros((2*M, 2*N))mu = 1 #Step size# lamda = 0.00005 # Step Size for regularizerdelta_x = np.zeros_like(x_k)# Initialize plottingfig = plt.figureim_obj = plt.imshow(x_k, cmap='gray', aspect='auto') title_obj = plt.title('Fista + TV Regulariation, Iteration {}'.format(0))for k in range(n_iters):# Update momentum parameterst_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2beta = (t_k - 1) / t_k1# for i in range(3):# delta_x[:,:,i] = convolve2d(x_k[:,:,i], kernel, mode='same')delta_x = gradient_tv(x_k)y_k1 = x_k - mu * (Ah(A(x_k, pad_psf) - meas, pad_psf) + 2 * lamda * delta_x)z_k1 = np.maximum(0, y_k1)x_k = z_k1 + beta * (z_k1 - z_k)t_k = t_k1z_k = z_k1if k == 0:im_obj.set_data((np.maximum(crop(z_k)/np.max(z_k),0)))plt.title('Fista + TV Regulariation, Iteration {}'.format(k))drawnow(fig)x_fista_tv_reg = z_k im_obj.set_data(np.maximum(crop(z_k)/np.max(z_k),0))plt.title('Fista + TV Regulariation, Iteration {}'.format(k))plt.showreturn x_fista_tv_regdef perform_deconvolution(polar_img: np.ndarray, psf:np.ndarray, pad_psf: np.ndarray, center: tuple, img: np.ndarray, iterations: int = 500, nsr: float = 1e-3):wiener_out = wiener_deconvolution(polar_img, psf, nsr=nsr)gd_out = vanilla_gd(polar_img, pad_psf, n_iters=iterations)fista_out = gd_fista(polar_img, pad_psf, n_iters=iterations)gd_fista_tikhonov_out = gd_fista_tikhonov_reg(polar_img, pad_psf, n_iters=iterations)gd_fista_tv = gd_fista_total_variation_reg(polar_img, pad_psf, n_iters=iterations, lamda=0.05)polar_img += polar_img.max * 1e-5rl_out = restore.richardson_lucy(polar_img, kernel, iterations=iterations)# 8. Display resultsf = plt.figure(figsize=(15, 15))plt.subplot(3, 3, 1)plt.imshow(img, cmap='gray')plt.axis('off')plt.title("Original Image")plt.subplot(3, 3, 2)wiener_out1 = contrast_stetch(wiener_out)wiener_recon = inverse_polar_transform(wiener_out1, center=center, output=img.shape)plt.title("Wiener Deconvolution")plt.axis('off')plt.imshow(wiener_recon, cmap='gray')plt.subplot(3, 3, 3)rl_out1 = contrast_stetch(rl_out)rl_recon = inverse_polar_transform(rl_out1, center=center, output=img.shape)plt.title("Richardson Lucy Deconvolution")plt.axis('off')plt.imshow(rl_recon, cmap='gray')plt.subplot(3, 3, 4)gd_out1 = contrast_stetch(crop(gd_out))gd_recon = inverse_polar_transform(gd_out1, center=center, output=img.shape)plt.title("GD Deconvolution")plt.axis('off')plt.imshow(gd_recon, cmap='gray')plt.subplot(3, 3, 5)fista_out1 = contrast_stetch(crop(fista_out))fista_recon = inverse_polar_transform(fista_out1, center=center, output=img.shape)plt.title("FISTA Deconvolution")plt.axis('off')plt.imshow(fista_recon, cmap='gray')plt.subplot(3, 3, 6)gd_fista_tikhonov_out1 = contrast_stetch(crop(gd_fista_tikhonov_out))gd_fista_tikhonov_recon = inverse_polar_transform(gd_fista_tikhonov_out1, center=center, output=img.shape)plt.title("FISTA+Reg Deconvolution")plt.axis('off')plt.imshow(gd_fista_tikhonov_recon, cmap='gray')plt.subplot(3, 3, 7)gd_fista_tv1 = contrast_stetch(crop(gd_fista_tv))gd_fista_tv_recon = inverse_polar_transform(gd_fista_tv1, center=center, output=img.shape)plt.title("FISTA+TV Deconvolution")plt.axis('off')plt.imshow(gd_fista_tv_recon, cmap='gray')f.suptitle(f'Deconvolution Results', fontsize=10)plt.show# Save them toorecon_images = {}recon_images['wiener_recon'] = wiener_reconrecon_images['rl_recon'] = rl_reconrecon_images['gd_recon'] = gd_reconrecon_images['fista_recon'] = fista_reconrecon_images['gd_fista_tikhonov_recon'] = gd_fista_tikhonov_reconrecon_images['gd_fista_tv_recon'] = gd_fista_tv_reconreturn recon_imagesres_path = "results/synth_05deg/"polar_img, psf, pad_psf, center, img = preprocessing_synthetic_images("real/stars3.jpg", th=0.2, deg=5, center=(400, 320))Kernel Length: 45 1.0000000000000002Predicted Center: [420. 312.]True Center: (400, 320)
recon_images = perform_deconvolution(polar_img, psf, pad_psf, center, img, iterations=500)save_images(recon_images, res_path)compute_stats(img, recon_images, res_path)
res_path = "results/synth_10deg/"polar_img, psf, pad_psf, center, img = preprocessing_synthetic_images("real/stars3.jpg", th=0.2, deg=10, center=(400, 320))Kernel Length: 89 1.0Predicted Center: [398. 320.]True Center: (400, 320)
res_path = "results/synth_15deg/"polar_img, psf, pad_psf, center, img = preprocessing_synthetic_images("real/stars3.jpg", th=0.2, deg=15, center=(400, 320))Kernel Length: 134 1.0Predicted Center: [418. 324.]True Center: (400, 320)
res_path = "results/synth_30deg/"polar_img, psf, pad_psf, center, img = preprocessing_synthetic_images("real/stars3.jpg", th=0.2, deg=30, center=(400, 320))Kernel Length: 268 0.9999999999999999Predicted Center: [416. 328.]True Center: (400, 320)
res_path = "results/real_30min/"polar_img, psf, pad_psf, center, img = preprocessing_natural_image("real/star_trails_30min_small.jpg", th=0.2, center=(480, 320), time=30)Predicted Center: [436. 358.]Kernel Shape: (1, 76) 1.0
res_path = "results/real_10min/"polar_img, psf, pad_psf, center, img = preprocessing_natural_image("real/star_trails_10min_small.jpg", th=0.2, center=(200, 350), time=10)
res_path = "results/real_05min/"polar_img, psf, pad_psf, center, img = preprocessing_natural_image("real/star_trails_5min_small.jpg", th=0.2, center=(450, 350), time=5)
知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。
分割线分割线分割线
基于小波分析的时间序列降噪（Python,ipynb文件）
完整代码：
时间序列的最大重叠离散小波分解（Python）
完整代码：
基于连续小波变换的信号滤波方法（Python，ipynb文件）
完整代码：
信号的时域、频域和时频域特征提取（Python）
完整代码：
不同小波族的优缺点
完整代码：
同步压缩变换的一些应用（时频分析，盲源分离等，MATLAB）
完整代码和数据可通过知乎学术咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）
引力波信号的同步压缩变换，高阶同步压缩变换，脊线提取与信号重建（MATLAB R2021B）
基于同步压缩变换（重分配算子机制）的非平稳信号瞬时频率估计（MATLAB R2018A）
基于优化Morlet小波的一维信号瞬态特征提取方法（MATLAB ）
程序运行环境为MATLAB R2018A，利用优化Morlet小波对一维信号进行瞬态特征提取。程序测试了模拟信号，地震信号，发动机销子活塞故障振动信号，发动机气门正常振动信号，发动机排气门故障振动信号，结果如下。
完整代码可通过知乎付费咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）
MATLAB环境下基于高分辨时频分析与阈值法的一维信号降噪
基于经验模态分解、离散小波变换和AR模型的肌电图信号降噪（Python，ipynb文件）
NS: noisy signalS: original siganlmean filter: ws = window sizemedian filter:average filter: ns = number of noisy signal(different)bandpass filter: l = low cut-off frequency, h = high ...threshold filter: r = ratio(max abs(fft) / min ...)wavelet filter: a = thresholdstd filter:NN: neural network完整代码
来源：张跃银看英语
免责声明：本站系转载，并不代表本网赞同其观点和对其真实性负责。如涉及作品内容、版权和其它问题，请在30日内与本站联系，我们将在第一时间删除内容!
