具有傅立叶域正则化的实时单像素视频成像附matlab代码.zip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% MATLAB script to reproduce figure 4 of "Real-time THz imaging with a single-pixel detector" %%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Script written by Rayko Stantchev %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % % % % % % % % % % NOTE % % % % % % % % % % % % % % % % %%%% %%%% %%%% The fdri.m was originally written by others! See DOI:10.1364/OE.26.020009 %%%% %%%% "Real-time single-pixel video imaging with Fourier domain regularization" %%%% %%%% Optics Express Vol. 26, Issue 16, pp. 20009-20022 (2018). %%%% %%%% The script is downloaded from %%%% %%%% https://www.igf.fuw.edu.pl/en/articles/popularyzacja-l/fdri-x2018-08-06/ %%%% %%%% But I modified it so it works with matrices of sizes other than 2^n with n integer! %%%% %%%% %%%% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %%%% %Get where script is located (does not work in with Run Section!) scriptlocation = fileparts(which(mfilename() )); cd( scriptlocation ) %Set current folder to where script is located %% Load data Paley_raw = disp_imgload( 'undersamp32_paley_dispimg' ); %% Process the raw signals %%%%% This is done by the function named bias_Texp. It does section Software Optimization of the manuscript raw_sigs = zeros(9, 2048); for i=1:9 processed_data{i} = bias_Texp(Paley_raw, i); raw_sigs(i, :) = processed_data{i}.Amps; end hline = Paley_raw.hline; n= Paley_raw.n; %% Plot the fitted signals i = 1; y = processed_data{i}.y; %Fitted signals x = processed_data{i}.x; %Raw signals t = processed_data{i}.t; %Fitted signals pp = 1:2500; %number of points to plot below figure plot(t(pp), x(pp), '.', t(pp), y(pp)); xlabel('Time (\mus)'); ylabel('Voltage (V)'); legend('Raw signals', 'Fitted signal') %% Fully sampled images %Ground Truth image ground_T = img_recon( mean( raw_sigs(:, 1:2:end) - raw_sigs(:, 2:2:end),1), hline, n); ground_T = flipud(ground_T'); %This is get images to be displayed correctly, (ie. up is up) %Single meassurement i = 1; sigsF = mean( raw_sigs(i, 1:2:end) - raw_sigs(i, 2:2:end) ,1); single_T = img_recon( sigsF, hline, n); single_T = flipud(single_T'); %1 - 0 meassurement sigs = mean(raw_sigs(i, 1:2:end), 1); sig01 = 2*(sigs(1:end) - mean(sigs(2:end)) ); sig01(1) = sigs(1); single_01 = img_recon( sig01, hline, n); single_01 = flipud(single_01'); figure subplot(1, 3, 1) imagesc(ground_T); axis image; colorbar; title('Ground Truth') subplot(1, 3, 2) imagesc(single_T); axis image; colorbar; title('Single image, [1 -1] masks') subplot(1, 3, 3) imagesc(single_01); axis image; colorbar; title('Single image, [1 0] masks') %% Undersampling, this is section Undersampling cc = linspace( 0.01, 1, 135 ); %Compression ratios pics10 = cell(1, length(cc)); pics_pn= cell(1, length(cc)); clear snr nh = length(hline); sigsF(end:nh+1) =0; sig01(end:nh+1) =0; % dispstat('', 'init') for j=1:length(cc) % dispstat( num2str(j)) %This function is useful for displaying evaluation number. Get it from the MATLAB % file exchange! %Sampling vector smp = round( linspace(1, nh+1, cc(j)*n^2) ); M = zeros( length(smp) , nh+1); for i=1:length(smp) if i==1 hline_t = ones(nh+1,1); else hline_t = [1 circshift(hline ,[1, smp(i)-2])]; end M(i, :) =hline_t(:); end y01 = sig01(smp)'; %1 and 0 Signals y_pn = sigsF(smp)'; %1 and -1 Signals [P, DD]=fdri(M, n, n, 0.95, 5e-9, 1); % get the reconsturction matirx x01=P*y01; x01=flipud( reshape(x01(1:n^2), 32, 32)'*n^2 ); x_pn=P*y_pn; x_pn=flipud( reshape(x_pn(1:n^2), 32, 32)'*n^2 ); pics10{j} = x01; pics_pn{j} = x_pn; snr.sim_10(j) = ssim(pics10{j}, ground_T); snr.sim_pn(j) = ssim(pics_pn{j}, ground_T); end rec_x = pinv(M)*y_pn; rec_x = rec_x(1:n^2); rec_x=flipud( reshape(rec_x, 32, 32)'*n^2 ); %% This plots the results like in the manuscript figure set(gcf, 'Units', 'centimeters', 'Position', [5 9 18.7 6.4] ); xx = 10.1; %Font size subplot(2, 4 , [1 5]) imagesc( ground_T ); axis image off; title('Ground Truth', 'FontSize', xx); set(gca,'FontSize', xx) g = colorbar('SouthOutside'); ylabel(g, 'Field strength (arb. unit)', 'FontSize', xx); set(g, 'Position', [-0.074 +0.01 0 -0.01]+get(g, 'Position') ) set(gca, 'Position', [-0.075 -0.08 0 0]+get(gca, 'Position')) text(-0.16, 1.12, 'a', 'FontSize', xx+1, 'Color', [0 0 0], 'Units', 'normalized',... 'VerticalAlignment', 'cap', 'FontWeight', 'bold') sigs_10 = snr.sim_10; %Signals to plot sigs_p = snr.sim_pn; ax = subplot(2, 4 , [2:4, 6:8]); plt_1 = plot(cc, sigs_10, '-', 2*cc, sigs_p, '-.'); grid on xlim([0 cc(end)]); ylim([-0.475 1.25]) ylabel('Structure similary index', 'FontSize', xx); xlabel('N0 of measurements / N0 of image pixels', 'FontSize', xx) set(gca,'FontSize', xx) set(gca, 'Position', [-0.0173 0.073 0.07 -0.072]+get(gca, 'Position')) text(0.01, 0.98, 'b', 'FontSize', xx+1, 'Color', [0 0 0], 'Units', 'normalized',... 'VerticalAlignment', 'cap', 'FontWeight', 'bold') % Adding the images ss = 0.179; %Image size plts_01 = round(linspace(4, length(cc)-1, 9)); %Images to plot for [1 0] masks plts_pn = round(0.5*linspace(4, length(cc)-1, 9)); %Images to plot for [1 -1] masks xxa = -[0.08 0.085 0.087 0.092 0.095 0.095 0.095 0.1 0.1]; %X shift values for images lxx = [0.01 0.02 0.00 -0.01 -0.02 -0.025 -0.025 -0.02 -0.02 ]; %X shift values for lines for i=1:length(plts_01) [xa1, ya] = ds2nfu(ax, cc(plts_01(i)) , sigs_10(plts_01(i)) ); %See ds2nfu function axes('Position', [xa1+xxa(i) ya*0.7+0.3 ss ss]) img01 = pics10{ plts_01(i) }; %image to plot imagesc(img01); axis off square; set(gca, 'xticklabel', [], 'yticklabel', [], 'FontSize', xx); text(16, 28, num2str( sigs_10(plts_01(i)), 3 ), 'Color', [1 1 1], 'FontSize', xx, 'HorizontalAlignment', 'center' ) line(ax, cc(plts_01(i)).*[1 1]+[0 lxx(i)], [sigs_10(plts_01(i)) 1.25*sigs_10(plts_01(i))+0.375], 'LineStyle', '-', 'Color', 'Blue') [xa2, ya] = ds2nfu(ax, 2*cc(plts_pn(i)) , sigs_p(plts_pn(i)) ); axes('Position', [xa1+xxa(i) 0.95*ya-0.24 ss ss]) img_pn = pics_pn{ plts_pn(i) }; %image to plot imagesc(img_pn); axis off square; set(gca, 'xticklabel', [], 'yticklabel', [], 'FontSize', xx); text(16, 28, num2str( sigs_p(plts_pn(i)) , 3 ), 'Color', [1 1 1], 'FontSize', xx, 'HorizontalAlignment', 'center' ) line(ax, 2*cc(plts_pn(i)).*[1 1]+[0 lxx(i)], [sigs_p(plts_pn(i)) sigs_p(plts_pn(i))-0.375 ], 'LineStyle', '-', 'Color', [1 0.64453125 0]) end legend(plt_1, {'[1 0] masks', '[1 -1] masks'}, 'Location', 'SouthEast', 'FontSize', xx) set(gcf, 'PaperSize', [17.82 6.2442])