使用matlab实现衍射受限成像系统2

题目1(扩展)

题干相关

  根据你学过的衍射受限系统成像的计算方法，以自己的图像作为输入$I_{input}$，模拟计算出经过衍射受限系统后的输出图像$I_{image}$
  衍射受限系统表达如下：
$$I_{image}=I_{input}\otimes PSF$$
  $\otimes$ 表示什么运算很好解答–>$\otimes$ 可以理解为卷积操作，即有一个卷积核并使用该卷积核处理输入图像后得输出图像
  PSF(point spread function),点扩散函数。描述了成像系统对点源或点对象的响应。–>将这个理解为和卷积核相似得矩阵即可

不同掩模板的判别标准

本博客主要针对不同点扩散函数的掩模板对结果进行扩充：
  光波经过小孔的衍射分布均满足基尔霍夫衍射公式，根据成像平面与孔的距离，采取不同的近似处理方法，主要分为三个区域：几何投影区，菲涅尔衍射区和夫琅禾费衍射区–>本质区别：成像距离的远近

• 几何投影区：光强分布就是衍射孔的几何形状
• 菲涅尔衍射区：衍射图样锐变得边缘消失，衍射图像与孔得相似性逐渐消失，而且衍射图像的中心产生明暗变化
• 夫琅禾费衍射区：衍射图像只有大小变化，形状基本不变，光强分布具有相似性

  在此基础上，本文实现了夫琅禾费圆孔衍射，夫琅禾费衍矩孔射，菲涅尔圆孔衍射和菲涅尔方孔衍射点扩散函数下的衍射受限系统成像计算方法，下面将给出具体的代码实现和结果图：

夫琅禾费圆孔衍射代码及结果图

%%
% 2024.08.06
% StarsZong
% main.m

clear, clc;
close all;

img = imread("lenna.jpg");
img_gray = rgb2gray(img);

% 夫琅禾费圆孔衍射
lambda = 500e-9;
r = 5e-4;    
z = 2;         
k = 255;
I = F(lambda, r, z, k);

img_out = imfilter(double(img_gray), I);
img_out1 = img_out ./ max(max(img_out));
img_out1 = uint8(255 * img_out1);

figure;
imshow(I);
colormap(hot(255));
figure;
imshow([img_gray, img_out1]);

% 2024.08.06
% StarsZong

function I = fun1(lambda, r, z, size)
% 夫琅禾费圆孔衍射掩模版
% lambda：   光的波长
% r：        圆孔的半径
% z：        光源与圆孔的距离
% size：     掩模板的大小（必须为奇数）

N = 1e-2;
[x, y] = meshgrid(linspace(-N, N, size));
I_0 = ((pi * r.^2) / (lambda * z)).^2;
s = 2 * pi * r * sqrt(x.^2+y.^2) ./ (lambda * z);
I = 4 * I_0 * (besselj(1,s) ./ (s+eps)).^2;
A = max(max(I));
for i = 1: length(x)
    for j =  1: length(y)
        if I(i, j) == 0
            I(i, j) = A;
        end
    end
end
I = 255 * I ./ max(max(I));
end

结果图：

夫琅禾费矩孔衍射代码及结果图

%%
% 2024.08.06
% StarsZong
% main.m

clear, clc;
close all;

img = imread("lenna.jpg");
img_gray = rgb2gray(img);

% 夫琅禾费矩孔衍射
lambda=500e-9;
a = 2e-3;
b = 1e-3;
z = 2;
size = 256;
I = fun2(lambda, a, b, z, size);

img_out = imfilter(double(img_gray), I);
img_out1 = img_out ./ max(max(img_out));
img_out1 = uint8(255 * img_out1);

figure;
imshow(I);
colormap(hot(255));
figure;
imshow([img_gray, img_out1]);

% 2024.08.06
% StarsZong

function I = fun2(lambda, a, b, z, size)
% 夫琅禾费矩孔衍射
% lambda：   光的波长
% a：        矩孔的长
% b：        矩孔的宽
% z：        光源与矩孔的距离
% size：     掩模板的大小（必须为偶数）

N = 1e-2;
x = linspace(-N, N, size);
for i = 1: size
    sinphi1 = x(i) / sqrt(x(i)^2 + z^2);
    sinphi2 = x ./ sqrt(x.^2 + z^2);
    u1 = pi * a * sinphi1 / lambda;
    u2 = pi * b * sinphi2 ./ lambda;
    I(:,i) = (sin(u1) ./ u1).^2 .* (sin(u2) ./ u2).^2;
end
I = 255 * (I ./ max(max(I)));
end

结果图：

菲涅尔圆孔衍射代码及结果图

%%
% 2024.08.06
% StarsZong
% main.m

clear, clc;
close all;

img = imread("lenna.jpg");
img_gray = rgb2gray(img);

% 菲涅尔圆孔衍射
lambda=500e-9;
r = 1e-3;
z = 2; 
size = 256;        
I = fun3(lambda, r, z, size);

img_out = imfilter(double(img_gray), I);
img_out1 = img_out ./ max(max(img_out));
img_out1 = uint8(255 * img_out1);

figure;
imshow(I);
colormap(hot(255));
figure;
imshow([img_gray, img_out1]);

% 2024.08.06
% StarsZong

function I = fun3(lambda, r, z, size)
% 菲涅尔圆孔衍射
% lambda：   光的波长
% r：        圆孔的半径
% z：        光源与圆孔的距离
% size：     掩模板的大小（奇偶均可）

N= 1e-2;
[x, y] = meshgrid(linspace(-N, N, size));
z1 = zeros(size, size);
z2 = zeros(size, size);
for i = 1: size
    for j = 1: size
        if x(i, j)^2 + y(i, j)^2 <= r^2
            z1(i, j) = 1;
            z2(i, j) = x(i, j)^2 + y(i, j)^2;
        end
    end
end
I = fftshift(abs((fft2(z1.*exp(1i*2*pi/lambda/2/z*z2))).^2));
I = 4 * (I / max(max(I)));
end

结果图

菲涅尔方孔衍射代码及结果图

%%
% 2024.08.06
% StarsZong
% main.m

clear, clc;
close all;

img = imread("lenna.jpg");
img_gray = rgb2gray(img);

% 菲涅尔方孔衍射
lambda=500e-9;
a = 2e-3;
z = 20; 
size = 255;        
I = fun4(lambda, a, z, size);

img_out = imfilter(double(img_gray), I);
img_out1 = img_out ./ max(max(img_out));
img_out1 = uint8(255 * img_out1);

figure;
imshow(I);
colormap(hot(255));
figure;
imshow([img_gray, img_out1]);

% 2024.08.06
% StarsZong

function I = fun4(lambda, a, z, size)
% 菲涅尔方孔衍射
% lambda：   光的波长
% a：        方孔的边长
% z:         光源与方孔的距离
% size：     掩模板的大小（奇偶均可）

N = 1e-2;
[x, y] = meshgrid(linspace(-N, N, size));
z1 = zeros(size, size);
z2 = zeros(size, size);
for i = 1: size
    for j = 1: size
        if abs(x(i, j)) <= a && abs(y(i, j)) <= a
            z1(i, j) = 1;
            z2(i, j) = x(i, j)^2 + y(i, j)^2;
        end
    end
end
I = fftshift(abs((fft2(z1.*exp(1i*2*pi/lambda/2/z*z2))).^2));
I = 255 * (I / max(max(I)));
end

结果图

代码相关

imread --> 读取图像
rgb2gray --> 将rgb彩色图转为灰度图
imfilter --> 将原图和卷积核输入，输出卷积操作后的图像
imshow --> 展示图象
meshfird --> 将一维数据关联为二维数据
besselj --> 贝塞尔函数

还有问题可在matlab内通过："help 函数名" 进行函数查看