A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分-Toy模板网

这篇具有很好参考价值的文章主要介绍了A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分

代码介绍原文：

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

1.代码结构

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

1.1 data数据

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

在这里插入图片描述 A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

2. scripts_preset.m 代码

2.1 代码整体介绍

这段代码的作用是加载PSFs数据，并进行一系列参数设置。

首先，通过设置data_dir变量为数据目录的路径。然后，根据场景名来选择特定于场景的参数。根据不同的场景名，设置d1和d2的值。
net 场景包括一个距离相机约30mm的网状物和一个距离网状物后方约100mm的橡皮鸭。

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

%% 加载 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

相机参数设置

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

% 相机参数预设
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;

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

% 循环读取不同掩膜的图像和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

重建

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

计算了正则化项 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;

其他操作

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

%% 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

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

% 调整对比度以便输出
% 创建存储调整对比度后的结果的文件夹
%% 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');

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab

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)));

A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分,matlab 文章来源地址https://www.toymoban.com/news/detail-847076.html

到了这里，关于A Simple Framework for 3D Lensless Imaging with Programmablle Masks 论文代码部分的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！