一种图像的星迹去除方法（Part 2）

欢迎大家来到IT世界,在知识的湖畔探索吧!

import numpy as np import matplotlib.pyplot as plt from scipy.signal import convolve2d from scipy.ndimage import binary_closing, binary_dilation, binary_erosion, binary_opening from numpy.fft import fft2, fftshift, ifft2, ifftshift, rfft2, irfft2 import matplotlib.patches as patches import skimage.restoration as restore from skimage.metrics import peak_signal_noise_ratio, structural_similarity import cv2 from scipy.ndimage import map_coordinates from IPython import display import time import os import tifffile as tiff

欢迎大家来到IT世界,在知识的湖畔探索吧!

欢迎大家来到IT世界,在知识的湖畔探索吧!def 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()

欢迎大家来到IT世界,在知识的湖畔探索吧!

Polar Transforms

欢迎大家来到IT世界,在知识的湖畔探索吧!def polar_transform(img, center, output=None): r = (np.array(img.shape[:2])2).sum()0.5/2 if 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.shape r_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 output

def 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.shape y_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.5 ts = np.arccos(x_samples / (rs + 1e-8)) ts[y_samples<0] = np.pi * 2 - ts[y_samples<0] ts *= (img.shape[1]-1) / (np.pi * 2) map_coordinates(img, (rs, ts), order=1, output=output) return output

Utils

欢迎大家来到IT世界,在知识的湖畔探索吧!def generate_star_trails(img, kernel): return convolve2d(img, kernel, mode='same')

def contrast_stetch(img): max = img.max() min = img.min() return (img - min) / (max - min)

欢迎大家来到IT世界,在知识的湖畔探索吧!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)

欢迎大家来到IT世界,在知识的湖畔探索吧!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.close()

Center

def find_center(img, threshold, ax=None): # Binarize Image img = contrast_stetch(img) bin_img = (img > threshold).astype('int') H, W = img.shape[:2] factor = 1 if (H > 1000 or W > 1000): factor = max(H, W) // 1000 # Dilating and subsampling bin_img = binary_dilation(bin_img, np.ones((factor, factor))) bin_img = bin_img[::factor, ::factor].astype(np.uint8) * 255 # To remove noise blur = 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 circular if ax is None: fig, ax = plt.subplots(1) ax.set_aspect('equal') # Show the image ax.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.show() circ = circles[0] circ[1::-1] = circ[:2] print("Predicted Center: ", circ[:2] * factor) return circ[:2] * factor 

Deconv

欢迎大家来到IT世界,在知识的湖畔探索吧!def 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/2 total_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_theta return estimated_kernel

Wiener Decovolution

def wiener_deconvolution(image, psf, nsr=0.01): """ image (np.array) : MxN psf (np.array) : MxN nsr (noise to signal ratio - 1/SNR from above derivation) : float """ # YOUR CODE HERE deblurred_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_image

Vanilla GD

欢迎大家来到IT世界,在知识的湖畔探索吧!def F(x): return rfft2(x,axes=(0,1)) def Fh(x): return irfft2(x,axes=(0,1))

def crop(x): M,N = x.shape rl = M//4 ru = 3*np.ceil(M/4) cl = N//4 cu = 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.shape pad_psf = np.zeros((2*H, 2*W)) pad_psf[M//2:M//2+psf.shape[0], N//2:N//2+psf.shape[1]] = psf return pad_psf

欢迎大家来到IT世界,在知识的湖畔探索吧!def 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)) / 100 psf = np.zeros_like(img) kh, kw = kernel.shape H, W = psf.shape psf[H//2-kh//2:H//2-kh//2+kh, W//2-kw//2:W//2-kw//2+kw] = kernel [M,N] = img.shape x = 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('<A(x),y> =',Ax_y) print('<x,Ah(y)>=',x_Ahy) 

欢迎大家来到IT世界,在知识的湖畔探索吧!(300, 400) <A(x),y> = 29904.2 <x,Ah(y)>= 29904.

def vanilla_gd(meas, pad_psf, n_iters=500): M, N = meas.shape x_k = np.zeros((2*M, 2*N)) #Initialize with zeros mu = 1 #Step size fig = plt.figure() #Initialize figure im_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 progress if k%10 == 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_k im_obj.set_data(np.maximum(crop(x_k)/np.max(x_k),0)) plt.title('Gradient Descent, Iteration {}'.format(k)) plt.show() return x_gd

FISTA

欢迎大家来到IT世界,在知识的湖畔探索吧!def gd_fista(meas, pad_psf, n_iters=500, mu=1): M, N = meas.shape # Initialize variables t_k = 1 #Momentum parameter x_k = np.zeros((2*M, 2*N)) # Image y_k = np.zeros((2*M, 2*N)) z_k = np.zeros((2*M, 2*N)) # Initialize plotting fig = plt.figure() im_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 parameters t_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2 beta = (t_k - 1) / t_k1 y_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_k1 z_k = z_k1 #print("Beta:", beta, " t_k: ", t_k) if k%10 == 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.show() return x_fista

FISTA + Regularization

def 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 variables t_k = 1 #Momentum parameter x_k = np.zeros((2*M, 2*N)) # Image y_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 regularizer delta_x = np.zeros_like(x_k) # Initialize plotting fig = plt.figure() im_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 parameters t_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2 beta = (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_k1 z_k = z_k1 if k%10 == 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.show() return x_fista_reg

欢迎大家来到IT世界,在知识的湖畔探索吧!def 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 grad M, N = meas.shape # Initialize variables t_k = 1 #Momentum parameter x_k = np.zeros((2*M, 2*N)) # Image y_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 regularizer delta_x = np.zeros_like(x_k) # Initialize plotting fig = plt.figure() im_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 parameters t_k1 = (1 + np.sqrt(1+4*t_k*t_k)) / 2 beta = (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_k1 z_k = z_k1 if k%10 == 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.show() return x_fista_tv_reg

Deconvoltion Manager

def 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-5 rl_out = restore.richardson_lucy(polar_img, kernel, iterations=iterations) # 8. Display results f = 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 too recon_images = {} recon_images['wiener_recon'] = wiener_recon recon_images['rl_recon'] = rl_recon recon_images['gd_recon'] = gd_recon recon_images['fista_recon'] = fista_recon recon_images['gd_fista_tikhonov_recon'] = gd_fista_tikhonov_recon recon_images['gd_fista_tv_recon'] = gd_fista_tv_recon return recon_images

Tests

Synthetic

欢迎大家来到IT世界,在知识的湖畔探索吧!res_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.0000000000000002 Predicted Center: [420. 312.] True Center: (400, 320)

欢迎大家来到IT世界,在知识的湖畔探索吧!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))

欢迎大家来到IT世界,在知识的湖畔探索吧!Kernel Length: 89 1.0 Predicted Center: [398. 320.] 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)

欢迎大家来到IT世界,在知识的湖畔探索吧!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.0 Predicted Center: [418. 324.] True Center: (400, 320)

欢迎大家来到IT世界,在知识的湖畔探索吧!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_30deg/" polar_img, psf, pad_psf, center, img = preprocessing_synthetic_images("real/stars3.jpg", th=0.2, deg=30, center=(400, 320))

欢迎大家来到IT世界,在知识的湖畔探索吧!Kernel Length: 268 0.99999 Predicted Center: [416. 328.] True Center: (400, 320)

Real Images

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) 

欢迎大家来到IT世界,在知识的湖畔探索吧!Predicted Center: [436. 358.] Kernel Shape: (1, 76) 1.0

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)

欢迎大家来到IT世界,在知识的湖畔探索吧!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) 

知乎学术咨询：

欢迎大家来到IT世界,在知识的湖畔探索吧!https://www.zhihu.com/consult/people/?isMe=1

担任《Mechanical System and Signal Processing》《中国电机工程学报》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测。

分割线分割线分割线

基于小波分析的时间序列降噪（Python,ipynb文件）

完整代码：

https://mbd.pub/o/bread/mbd-Z5yZkppv

时间序列的最大重叠离散小波分解（Python）

完整代码：

欢迎大家来到IT世界,在知识的湖畔探索吧!https://mbd.pub/o/bread/mbd-Z5yZmZ1s

基于连续小波变换的信号滤波方法（Python，ipynb文件）

完整代码：

https://mbd.pub/o/bread/mbd-Z5yZlJpx

信号的时域、频域和时频域特征提取（Python）

完整代码：

欢迎大家来到IT世界,在知识的湖畔探索吧!https://mbd.pub/o/bread/mbd-Z5yalJ1s

不同小波族的优缺点

完整代码：

https://mbd.pub/o/bread/mbd-Z5yak59x

同步压缩变换的一些应用（时频分析，盲源分离等，MATLAB）

完整代码和数据可通过知乎学术咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）

欢迎大家来到IT世界,在知识的湖畔探索吧!https://www.zhihu.com/consult/people/

引力波信号的同步压缩变换，高阶同步压缩变换，脊线提取与信号重建（MATLAB R2021B）

完整代码和数据可通过知乎学术咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）

https://www.zhihu.com/consult/people/

基于同步压缩变换（重分配算子机制）的非平稳信号瞬时频率估计（MATLAB R2018A）

完整代码和数据可通过知乎学术咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）

欢迎大家来到IT世界,在知识的湖畔探索吧!https://www.zhihu.com/consult/people/

基于优化Morlet小波的一维信号瞬态特征提取方法（MATLAB ）

程序运行环境为MATLAB R2018A，利用优化Morlet小波对一维信号进行瞬态特征提取。程序测试了模拟信号，地震信号，发动机销子活塞故障振动信号，发动机气门正常振动信号，发动机排气门故障振动信号，结果如下。

完整代码可通过知乎付费咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）

https://www.zhihu.com/consult/people/

MATLAB环境下基于高分辨时频分析与阈值法的一维信号降噪

完整代码可通过知乎付费咨询获得（注意：一次知乎学术付费咨询只能获得一套代码＋学术指导）

欢迎大家来到IT世界,在知识的湖畔探索吧!https://www.zhihu.com/consult/people/

基于经验模态分解、离散小波变换和AR模型的肌电图信号降噪（Python，ipynb文件）

完整代码：

https://mbd.pub/o/bread/mbd-ZpyTm55x

采用8种方法对一维信号进行降噪（Python）

采用8种方法对一维信号进行降噪（Python）

• NS: noisy signal
• S: original siganl
• mean filter: ws = window size
• median 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 = threshold
• std filter:
• NN: neural network

完整代码

欢迎大家来到IT世界,在知识的湖畔探索吧!https://mbd.pub/o/bread/mbd-ZpyTmJZy