一个有趣的光学图案绘制程序

最新推荐文章于 2022-07-08 15:05:48 发布

Matrix_11

最新推荐文章于 2022-07-08 15:05:48 发布

阅读量650

点赞数 1

分类专栏：计算摄影与图像处理 Photoshop 图像处理算法

本文链接：https://blog.csdn.net/matrix_space/article/details/105499667

版权

Photoshop 图像处理算法同时被 2 个专栏收录

92 篇文章 180 订阅

订阅专栏

计算摄影与图像处理

71 篇文章 74 订阅

订阅专栏

一个有趣的光学图案绘制程序

$\lambda = 600 nm$
$\alpha = 0.5$
$\lambda$
$\pi / \lambda$
$w_0 = 5 \lambda$
$z_0 = k \cdot w_0 \cdot w_0 / 2$
$\alpha \cdot z_0$
$\sqrt {2 \pi p!/(p+s)!}$
$w_0 \sqrt{1+(z/z_0)^2}$
$z_1 = 50$

$\frac{1}{w(z)^2}(1 + j\frac{z}{z_0}), \quad b^{*} = \frac{1}{w(z)^2}(1 - j\frac{z}{z_0})$

$C_1 = \frac{1}{2^s} \left( \frac{A \pi k}{z} \right)^2 \frac{(b^{*} w(z))^{2p} }{w_0^{2(s+p)}}$

$r_n = n \frac{a}{N}, \quad \varphi_m = (m-1)\frac{2 \pi}{M}$

$C_2 = exp(-j\frac{k}{z_1}(x\cdot r_n \cdot cos(\varphi_m) + y \cdot r_n \cdot sin(\varphi_m)))$

$C_3 = (\frac{k r_n}{z_1})^s \cdot exp(-\frac{1}{4b}(\frac{kr_n}{z_1})^2)$

$C_4 = L_p^s (\frac{1}{4|b|^2}(\frac{kr_n}{z_1})^2) \cdot exp(-js\varphi_m)$

$C_s = \sum_{m=1}^{M} \sum_{n=1}^{N} C_2 \cdot C_3 \cdot C_4$

$z_1) = C_1 \cdot |C_s|^2$

Python 程序代码如下：

import numpy as np
import matplotlib.pyplot as plt
import math
import cmath

PI = math.pi

def Lps(p, s, x):
    
    if p == 0:
        y = 1
    elif p == 1:
        y = -x + s + 1
    elif p == 2:
        y = 0.5 * x * x - (s+2) * x + (s+1)*(s+2)/2 
    else:
        y = 0
    
    return y


def fac(x):
    return math.factorial(x)

# p 取值 0,1,2
# s 取值 0,1,2
# alpha 取值 0-1
    
p = 0
s = 1
alpha = 0.5

lam = 550 * 1e-9
Num = 200
a = 5 * lam
k = 2 * PI / lam
w0 = 5 * lam
z0 = k * w0 * w0 / 2
z = alpha * z0

A = (2 * PI * fac(p) / fac(p + abs(s))) ** 0.5
w_z = w0 * math.sqrt(1+(z/z0)**2)

b= complex(1.0/(w_z **2), (z/z0)/(w_z **2))
con_b = complex(1.0/(w_z **2), -(z/z0)/(w_z **2))
z1 = 50

I = np.zeros((Num * 2, Num * 2))
C1_1 = 1 / (2 **s)
C1_2 = (A * PI * k / z1) ** 2
C1_3 = (con_b * w_z) ** (2 * p) / (w0 ** (2*(p+s)))
C1 = C1_1 * C1_2 * C1_3 

M = 8
N = 3

Cen_x = Num + 0.5
Cen_y = Num + 0.5

for yy in range(Num*2):
    for xx in range(Num*2):
        xi = xx - Cen_x
        yi = Cen_y - yy
        Sum_xy = 0 
        for n in range(1, N+1):
            r_n = n * a / N
            for m in range(M):
                phi_m = m * 2 * PI / M
                c2_1 = yi * r_n * math.sin(phi_m) + xi * r_n * math.cos(phi_m)
                c2_j = complex(0, -k/z1 * c2_1)
                # c2_j = complex(y0*r_n*math.sin(phi_m), -k/z1*(x0*r_n* math.cos(phi_m)))
                c2 = cmath.exp(c2_j)
                c3 = ((k*r_n/z1) **p) * (cmath.exp(-(k*r_n/z)**2/(4*b)))
                c4_1 = (k*r_n/z1) ** 2/(4 * abs(b) **2)
                c4 = Lps(p, s, c4_1) * cmath.exp(complex(0, -s * phi_m))
                Sum_xy = Sum_xy + c2 * c3 * c4
        
        I[yy, xx] = abs(C1 * (abs(Sum_xy) ** 2))


I_image=I/np.max(I)

cv2.imwrite('Img_out_3.jpg', I_image * 255)

plt.figure()
plt.imshow(I_image, cmap='gray')
plt.show()