A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分
代码介绍原文：
1.代码结构
1.1 data数据
- net
在这里插入图片描述
2. scripts_preset.m 代码
2.1 代码整体介绍
这段代码的作用是加载PSFs数据，并进行一系列参数设置。
-
首先，通过设置
data_dir变量为数据目录的路径。然后，根据场景名来选择特定于场景的参数。根据不同的场景名，设置d1和d2的值。 -
net 场景包括一个距离相机约30mm的网状物和一个距离网状物后方约100mm的橡皮鸭。
%% 加载 psfs 数据
data_dir = '../data'; % 设置数据目录
scene_name = 'net'; % 设置场景名称
% -- SCENE
% 选择特定于场景的参数
switch scene_name
% 若场景名为“toys”（玩具），设置d1为80，d2为600
case 'toys'
d1 = 80; d2 = 600;
% 若场景名为“two_cards”（两张卡片），设置d1为35，d2为200
case 'two_cards'
d1 = 35; d2 = 200;
% 若场景名为“net”（网），设置d1为20，d2为300
case 'net'
d1 = 20; d2 = 300;
% 若场景名为“carrot”（萝卜），设置d1为15，d2为300
case 'carrot'
d1 = 15; d2 = 300;
end
-
设置
masktype变量来选择掩膜类型，可选的类型有"random"（随机）、“mls”(最小二乘)、“opt”(优化)和"shiftXY"(XY移位)。根据选择的掩膜类型不同，代码会设置不同的掩膜序列。 -
设置
nmask和ndepth的值，分别表示掩膜数量和深度数量。% -- MASKTYPE % 设置掩膜类型，可选项有 "random"(随机)、"mls"(最少二乘)、"opt"(优化) 和 "shiftXY"(XY移位) % masktype = 'random' % masktype = 'mls' masktype = 'opt' % masktype = 'shiftXY' % -- NUM OF MASKS AND DEPTHS % 设置掩膜数量为8，深度数量也为8 nmask = 8; ndepth = 8; -
接下来，设置
recon_method变量来选择重建方法，可选的方法有"adap_l2"（自适应L2）和"t_svd"（基于是否有GPU设备设置为张量奇异值分解）。 -
如果没有检测到
GPU设备，则将重建方法设置为"t_svd"。% -- METHOD % 设置重建方法为 'adap_l2'（自适应L2），或者基于是否有GPU设备设置为 't_svd'（张量奇异值分解） recon_method = 'adap_l2' % recon_method = 't_svd' % 如果没有检测到GPU设备，则将重建方法设置为't_svd' if gpuDeviceCount < 1 recon_method = 't_svd'; end -
加载文件
% 根据设置的数据目录加载PSFs数据文件 load(sprintf('%s/psfs-fullsize-cg-bin8.mat', data_dir)); disp(dists)-
字典内包含三部分
-
dists
-
mask_names
-
psfs
-
-
-
相机参数设置
% 相机参数预设 Ho = 3648 / 2; %图像的高度 Wo = 5472 / 2; %图像的宽度 mask_depth = 10.51; %掩膜距离传感器的距离 bin_size = 8; cam_pitch = 2.4 * 1e-3; %相机的像素距离，以毫米为单位 cam_chnl_pitch = cam_pitch * 2 * bin_size; %是通道间的距离，通过将相机像素间距乘以2和二进制大小来计算。 d = mask_depth; z2a = @(z) 1-d./z; % 深度到a的转换函数 a2z = @(a) d./(1-a); % a到深度的转换函数 % 计算深度范围 depth_range = a2z(linspace(z2a(d1), z2a(d2), ndepth)); -
根据选择的掩膜类型，选择原始测量数据和
PSFs（点扩散函数）进行恢复。%% 选择原始测量数据和PSFs进行恢复 % 根据掩膜类型选择随机二值掩膜 switch masktype % 如果掩膜类型为'random'，选择序列中的第6至第25个掩膜 case 'random' mask_seq = 6:25; ind1 = mask_seq(1:nmask); ind2 = mask_seq(nmask+1)*ones(length(ind1),1); % 如果掩膜类型为'mls'，选择序列中的第47至第66个掩膜 case 'mls' mask_seq = 47:66; ind1 = 47:2:47+nmask*2-1; ind2 = 48:2:47+nmask*2-1; % 如果掩膜类型为'opt'，选择序列中的第27至第46个掩膜 case 'opt' mask_seq =27:46; ind1 = 27:2:27+nmask*2-1; ind2 = 28:2:27+nmask*2-1; % 如果掩膜类型为'shiftXY'，选择序列中的第67至第102个掩膜 case 'shiftXY' mask_seq = 67:102; ind1 = 67:2:67+nmask*2-1; ind2 = 68:2:67+nmask*2-1; end
3 demo_block_diagonal.m 代码
3.1 代码介绍
-
参数初始化
%定义一个函数，该函数接受一个图像数组和图像名称，将其转化为RGB格式并播放 %H×W×C×N ===> H×W×N×C play_rggb_stack = @(im, names)play_stack(rggb2rgb(permute(im, [1 3 4 2])), names); % 定义一个函数，用于显示归一化的图像 imshowsm = @(im)imshow(mat2gray(squeeze(im))); % 定义一个函数，用于找到最大值并规范化 max1 = @(x) x/max(abs(x(:))); run('scripts_preset'); %设置实验ID expid = 1; %% % 初始化图像大小和点扩散函数（PSF）大小 ims_size = [228, 342]; psf_size = [228, 342]; % 获取掩膜数量 nmask = length(ind2); % 初始化存储图像序列和点扩散函数校准序列的数组 ims_stack = zeros([nmask, 4, ims_size]); psf_calib_stack = zeros([nmask,length(dists), 4, psf_size]); mask_stack = zeros(nmask, 769, 1024); % 定义一个函数，用于将数组填充到最大尺寸 padmax = @(x) padarray(x, (psf_size - size(x))/2, 0, 'both'); -
读取数据
-
为什么要存储掩膜数据的差值？
GPT: 存储差值的目的是为了在后续的重建过程中，使用这些差值来恢复被遮挡的信息。差值代表了两个不同掩膜之间的区别，这 些差值可以用于去除遮挡区域或恢复遮挡区域的信息。
通过存储差值，可以在后续的重建过程中利用它们来改善图像重建的质量和准确性。
-
data_dir = ‘…/data’
scene_name = ‘net’
-
mask_seq =27:46;
ind1 = 27:2:27+nmask2-1;
ind2 = 28:2:27+nmask2-1;
% 循环读取不同掩膜的图像和PSF校准数据 for ii = 1:nmask % 读取图像数据 im2 = im2double(imread(sprintf('%s/%s/%s.png',data_dir, scene_name, mask_names{ind2(ii)}))); im1 = im2double(imread(sprintf('%s/%s/%s.png',data_dir, scene_name, mask_names{ind1(ii)}))); % 读取掩膜数据并存储差值 mask1 = im2double(imread(sprintf('%s/masks/%s.png',data_dir, mask_names{ind1(ii)}))); mask2 = im2double(imread(sprintf('%s/masks/%s.png',data_dir, mask_names{ind2(ii)}))); mask_stack(ii, :,:) = mask1 - mask2; % 对图像数据进行去马赛克，并对像素进行合并 im_rggb = debayer(im1 - im2); % (height, width, channels) 重排列为 (channels, height, width)，然后均值池化 im_rggb = bin_pixels( reshape(permute(im_rggb, [3 1 2]), [1, 4, Ho, Wo]), bin_size); ims_stack(ii,:,:,:) = im_rggb; % 循环不同的景深处理点扩散函数（PSF） for use_d = 1:5 psf = padmax(psfs{ind1(ii), use_d}) - padmax(psfs{ind2(ii), use_d}); psf = repmat(reshape(psf, 1, 1, psf_size(1), psf_size(2)), 1, 4, 1, 1); psf_calib_stack(ii,use_d,:,:,:) = psf; end fprintf('\n read in data for mask %02d', ii); end**附：**bin_pixels函数，这个函数是一个用于将图像进行均值池化的函数。
输入参数:
- im: 输入图像，维度为(N, C, H, W)，其中N为图像的数量，C为图像的通道数，H为图像的高度，W为图像的宽度。
- bin_size: 池化操作的窗口大小。
function out = bin_pixels(im, bin_size) %im: (N, C, H, W) %bin_size:int kernel = ones(bin_size); kernel = kernel / sum(kernel(:)); half_bin = round(bin_size/2); ys = 1:bin_size:(size(im,3)-bin_size+1); xs = 1:bin_size:(size(im,4)-bin_size+1); out = zeros(size(im,1), size(im,2), length(ys), length(xs)); for i = 1:size(im,1), for j = 1:size(im,2), tmp = conv2(squeeze(im(i,j,:,:)), kernel, 'valid'); out(i,j,:,:) = tmp(ys, xs); end; end end -
-
计算每个掩膜对应的多深度PSF，创建PSF和图像的FFT堆栈
% 设置深度范围 im_z = depth_range; ndepth = length(im_z); zs = im_z; % 坐标转换 cam_size = size(ims_stack); cam_size = cam_size(3:end); C = size(ims_stack,2); %初始化PSF堆栈 psf_stack = zeros([nmask, ndepth, 4, psf_size]); % 循环计算每个掩膜对应的多深度PSF for ii = 1:nmask for d = 1:ndepth z = zs(d); z_ind = 2; % Compute beta, psf at z, by interpolating from nearest captured. % 通过插值计算z处的PSF beta = zeros([C psf_size]); % scale_fun = @(z, z0)((z0-mask_depth) * z)/((z-mask_depth) * z0); %d in mm scale_fun = @(z, z0)(z0 * (z + mask_depth)) / (z * (z0 + mask_depth)); % 深度以毫米为单位 a1_over_a0 = scale_fun(z, dists(z_ind)); % direct interp for right size直接插值到正确的尺寸 [x, y] = meshgrid(linspace(-.5, .5, psf_size(2)), linspace(-.5, .5, psf_size(1))); for c = 1:C alpha = anti_alias(squeeze(psf_calib_stack(ii,z_ind , c, :, :)), a1_over_a0); beta(c, :, :) = interp2(x, y, alpha, ... x / a1_over_a0, y / a1_over_a0, ... 'linear', 0); end psf_stack(ii, d, :,:,:) = beta; end end %% create fft stacks of PSF and ims%% 创建PSF和图像的FFT堆栈 full_size = ims_size + psf_size - 1; % 定义完整尺寸 psf_fft_stack = zeros([nmask, ndepth, 4, full_size]); % 初始化PSF的FFT堆栈 ims_fft_stack = zeros([nmask, 4, full_size]); % 初始化图像的FFT堆栈 -
傅里叶变换
for ii = 1:nmask for di = 1:ndepth for c = 1:C % 对PSF和图像进行FFT变换 psf_fft_stack(ii, di, c, :,:) = fft2(squeeze(psf_stack(ii, di, c, :,:)), full_size(1),full_size(2)); ims_fft_stack(ii,c,:,:) = fft2( pad_and_center(squeeze(ims_stack(ii, c, :,:)), full_size), full_size(1), full_size(2)); end end end -
重建
计算了正则化项
lambda_mat，然后将其添加到AtA_gpu中。正则项的计算是通过将最大的特征值乘以0.05来实现的。这样，代码实现了在最小二乘问题中添加正则项的功能。lambda_mat = repmat(eye(ndepth),[1,1, size(Ac_piece, 3) ]) .* reshape( max(reshape(squeeze(sqrt(sum(sum(abs(AtA_gpu).^2,1),2))), [], 1) * 0.05,200) ,1,1,size(Ac_piece,3)) ; AtA_gpu = AtA_gpu + lambda_mat;for c = 1:C Ac = squeeze(psf_fft_stack(:,:,c,:,:)); % 提取c通道的PSF FFT bc = ims_fft_stack(:,c,:,:); % 提取c通道的图像FFT switch recon_method % 根据重建方法进行选择 case 'adap_l2' % reconstruction with GPU npiece = 35; % 将数据划分为多个块以适应GPU运算 Ac_piece = reshape(Ac, nmask, ndepth, [], npiece); bc_piece = reshape(bc, nmask, 1 , [], npiece); xhat_fft_gpu = zeros(ndepth, 1, size(Ac_piece,3), npiece); for pp = 1:npiece % 将数据移至GPU进行计算 % move Ac and bc onto gpu Ac_gpu = gpuArray(Ac_piece(:,:,:,pp)); bc_gpu = gpuArray(bc_piece(:,:,:,pp)); % L2 regularization based on frobenius norm of AtA基于AtA的Frobenius范数的L2正则化 AtA_gpu = pagefun( @mtimes, pagefun(@ctranspose, Ac_gpu), Ac_gpu); %lambda_mat = repmat(eye(ndepth),[1,1, size(Ac_piece, 3) ]) .* reshape( max(vec(squeeze(sqrt(sum(sum(abs(AtA_gpu).^2,1),2)))) * 0.05,200) ,1,1,size(Ac_piece,3)) ; lambda_mat = repmat(eye(ndepth),[1,1, size(Ac_piece, 3) ]) .* reshape( max(reshape(squeeze(sqrt(sum(sum(abs(AtA_gpu).^2,1),2))), [], 1) * 0.05,200) ,1,1,size(Ac_piece,3)) ; AtA_gpu = AtA_gpu + lambda_mat; clear lambda_mat %计算重建的FFT xhat_fft_gpu(:,:,:,pp) = gather(pagefun(@mtimes, pagefun(@mldivide, ... AtA_gpu, ... pagefun(@ctranspose, Ac_gpu)), bc_gpu)); end xhat_fft_gpu = reshape(xhat_fft_gpu, ndepth, 1, []); xhat_FFT_stack = permute(gather(reshape(xhat_fft_gpu,ndepth, full_size(1),full_size(2))),[2,3,1]); %通过IFFT换回空间域并取实部 xhat_full(:,:,c,:) = real(ifft2(xhat_FFT_stack)); case 't_svd' % reconstruction without GPU (省略) end fprintf('\n%02d/%02d color channels completed', c, C); end -
其他处理
% 将RGGB转换为RGB格式 xhat_full_rgb = rggb2rgb(permute(xhat_full,[4,1,2,3])); %% rotate/shift % 创建保存文件夹 % 根据实验ID、掩码类型和深度索引创建一个保存文件夹路径 save_folder = fullfile('../results',sprintf('%03d',expid),sprintf('masktype_%s',masktype),sprintf('bd_nd%02d_nm%02d',ndepth,length(ind2))); % 如果文件夹不存在，则创建该文件夹 if ~exist(save_folder,'dir') mkdir(save_folder); end % 图像预处理 % 将xhat_full_rgb转换为灰度图像格式 xhat_full_rgb = mat2gray(xhat_full_rgb); % 初始化一个新的图像数组xhat_rgb xhat_rgb = zeros([ ndepth ,cam_size, 3 ]); % 旋转/移位图像 % 对于每个深度层，将图像旋转180度并根据点扩散函数尺寸进行移位裁剪 for di = 1:ndepth xhat_tmp = squeeze(xhat_full_rgb(di,:,:,:)); xhat_tmp = circshift(rot90(xhat_tmp,2), -(psf_size-1)); xhat_tmp = xhat_tmp(1:cam_size(1), 1:cam_size(2),:); xhat_rgb(di, :,:,:) = xhat_tmp; end % % 将处理后的RGB图像再次转换为灰度图像格式 xhat_rgb = mat2gray(xhat_rgb);
-
去噪
%% BM3D denoising % tic; % xhat_rgb_denoised = zeros(size(xhat_rgb)); % for di = 1:ndepth % for ci = 1:3 % xhat_rgb_denoised(di,:,:,ci) = BM3D(squeeze(xhat_rgb(di,:,:,ci)),0.035); % end % fprintf('\n depth %02d/%02d is denoised', di, ndepth); % end % xhat_rgb = xhat_rgb_denoised; % t_bm3d = toc; %% TV-L1 denoising (optional denoising method) tic; xhat_rgb_denoised = zeros(size(xhat_rgb)); for di = 1:ndepth for ci = 1:3 xhat_rgb_denoised(di,:,:,ci) = TVL1denoise(squeeze(xhat_rgb(di,:,:,ci)),1.3,100); end fprintf('\n depth %02d/%02d is denoised', di, ndepth); end xhat_rgb = xhat_rgb_denoised; t_tvl1 = toc; -
其他操作
%% crop images % 裁剪图像 % 对图像进行裁剪，只保留感兴趣的区域 xhat_rgb_crop = permute((xhat_rgb(:,40:187,36:309,:)),[2,3,4,1]); % 创建存储裁剪后原始结果的文件夹 save_folder1 = fullfile(save_folder,'raw_results'); if ~exist(save_folder1, 'dir') mkdir(save_folder1); end % 将裁剪后的图像保存为图像文件并创建动画gif文件 for di = 1:ndepth xhat_rgb_crop(:,:,:,di) = (xhat_rgb_crop(:,:,:,di)); imagesc(xhat_rgb_crop(:,:,:,di)); axis image off; % title(sprintf('%3.5g mm', depth_range(di))); export_fig(fullfile(save_folder1,sprintf('final_result_scene_%s_depth_%02d',scene_name,di))); frame = export_fig; ims{di} = frame; end filename = fullfile(save_folder1,sprintf('anime_result_scene_%s.gif',scene_name)); % Specify the output file name for idx = 1:ndepth [A,map] = rgb2ind(ims{idx},256); if idx == 1 imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',1*(8/ndepth)); else imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',1*(8/ndepth)); end end
% 调整对比度以便输出 % 创建存储调整对比度后的结果的文件夹 %% adjust contrast for exporting save_folder2 = fullfile(save_folder,'color_adjusted_results'); if ~exist(save_folder2, 'dir') mkdir(save_folder2); end % 设置对比度调整后处理参数 % color adjustment post-processing parameter if strcmp(recon_method,'adap_l2') r = 0.005; else r = 0.001; end % 调整图像对比度并保存为图像文件 xhat_tmp = permute(xhat_rgb_crop, [1,2,4,3]); xhat_tmp_size = size(xhat_tmp); xhat_tmp_vec = reshape(xhat_tmp, [],1,3); xhat_tmp_adjust = imadjust(xhat_tmp_vec, stretchlim(xhat_tmp_vec, [r,1-r]),[]); xhat_adjust = permute(reshape(xhat_tmp_adjust, xhat_tmp_size),[1,2,4,3]); ims = {}; for di = 1:ndepth imagesc(xhat_adjust(:,:,:,di)); axis image off; export_fig(fullfile(save_folder2,sprintf('blockDiagonal_scene_%s_depth_%02d',scene_name,di))); frame = export_fig; ims{di} = frame; end filename = fullfile(save_folder2,sprintf('anime_result_scene_%s_adjusted.gif',scene_name)); % Specify the output file name for idx = 1:ndepth [A,map] = rgb2ind(ims{idx},256); if idx == 1 imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',1*(8/ndepth)); else imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',1*(8/ndepth)); end end % 保存所有处理后的数据 % 将所有相关的变量和处理后的图像数据保存到.mat文件中 %% save all the processed data %xhat_full_rgb: 完整的RGB图像数据 %xhat_rgb: RGB图像数据 %xhat_rgb_crop: 裁剪后的RGB图像数据 %xhat_adjust: 调整后的RGB图像数据 %depth_range: 深度范围 %scene_name: 场景名称 %masktype: 掩模类型 save(fullfile(save_folder,sprintf('results_data_%s.mat', scene_name ) ),'xhat_full_rgb','xhat_rgb','xhat_rgb_crop','xhat_adjust','depth_range','scene_name','masktype');
文章来源：https://www.toymoban.com/news/detail-847076.html
4 generate_all_in_focus_rggb函数 生成全在焦图
%清除变量和关闭所有打开的图像界面
clear; close all;
%将当前工作目录添加到路径中
addpath(genpath(pwd));
run('demo_block_diagonal');
% 定义图像处理参数
fstack_params = [];
fstack_params.filtersize = 5; % 滤波器大小
fstack_params.logstd = 2 ; % 标准差
fstack_params.dilatesize = 13; % 膨胀大小
fstack_params.blendsize = 13; % 混合大小
fstack_params.blendstd = 5; % 混合标准差
fstack_params.logthreshold = 0; % 对数阈值
% 生成全焦深度图像和深度信息
[image_rggb, depth_cone] = generate_all_in_focus_rggb(permute(xhat_rgb_crop,[4,3,1,2]) , fstack_params, depth_range, 0);
%% 显示生成的全焦深度图像和深度信息
subplot(121);
imshow(mat2gray(squeeze((permute(image_rggb,[1,3,4,2])))))
subplot(122);
imshow(depth_cone,[]);colormap parula
%%
[image_rggb, depth_cone] = generate_all_in_focus_rggb(permute(xhat_rgb,[1,4,2,3]) , fstack_params, depth_range, 0);
image_rggb = image_rggb(:,:,40:197,36:309);
depth_cone = depth_cone(40:197,36:309);
%
subplot(121);
imshow(mat2gray(squeeze((permute(image_rggb,[1,3,4,2])))))
subplot(122);
imshow(depth_cone,[]);colormap parula
%%
% % scene_name='ruler'
% % scene_name = 'toys'
% scene_name = 'two_cards'
%
% % -- MASKTYPE
% % masktype = 'random'
% masktype = 'opt'
% % masktype = 'scaleXY'
% % masktype = 'shiftXY'
%%% 定义保存结果的文件夹路径
save_folder = fullfile('../progMask_results','1202_final_wiener');
if ~exist(save_folder,'dir')
mkdir(save_folder);
end
% 保存生成的全焦深度图像
figure(1);
imshow(mat2gray(squeeze((permute(image_rggb,[1,3,4,2])))))
export_fig(fullfile(save_folder, sprintf('image_%s_%s',scene_name,masktype)));
% 保存生成的深度信息
figure(1);
imshow(depth_cone,[]);colormap parula
export_fig(fullfile(save_folder, sprintf('depth_%s_%s',scene_name,masktype)));
文章来源地址https://www.toymoban.com/news/detail-847076.html
到了这里，关于A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
![[论文阅读]Ghost-free High Dynamic Range Imaging with Context-aware Transformer](https://imgs.yssmx.com/Uploads/2024/02/733710-1.jpg)
