diff --git a/AugmentedInSilico/Fig3b_silicoaug_RSULM.m b/AugmentedInSilico/Fig3b_silicoaug_RSULM.m new file mode 100644 index 0000000000000000000000000000000000000000..7bebc6c0c5cc6a9821153d8371ecad17a0bdea87 --- /dev/null +++ b/AugmentedInSilico/Fig3b_silicoaug_RSULM.m @@ -0,0 +1,167 @@ +%% silicoaug_RSULM.m +% +% +% Based on the code created by Arthur Chavignon +% Adapted for the Augmented in Silico Data. +% +% DATE 2020.12.17 - VERSION 1.1 +% AUTHORS: Arthur Chavignon, Baptiste Heiles, Vincent Hingot. CNRS, Sorbonne Universite, INSERM. +% Laboratoire d'Imagerie Biomedicale, Team PPM. 15 rue de l'Ecole de Medecine, 75006, Paris +% Code Available under Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (see https://creativecommons.org/licenses/by-nc-sa/4.0/) +% ACADEMIC REFERENCES TO BE CITED +% Details of the code in the article by Heiles, Chavignon, Hingot, Lopez, Teston and Couture. +% Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy, Nature Biomedical Engineering, 2021. +% General description of super-resolution in: Couture et al., Ultrasound localization microscopy and super-resolution: A state of the art, IEEE UFFC 2018 + +clear +%cd ./.. +%cd('PALA_silico'); + +PALA_folder = fileparts(pwd); % location of the PALA folder +% DEFINE THE ADDONS DIRECTORY ON YOUR COMPUTER +addpath(genpath(fullfile(PALA_folder, 'PALA_addons'))); + +%% Data parameters +% A few parameters must be provided by the user depending of your input images (size of pixel, wavelength) +% These parameters will be copied and used later during the creation of the ULM structure. +% in this example, UF.TwFreq = 15MHz, UF.FrameRateUF = 1000Hz; +UF.TwFreq = 14.8; +UF.FrameRateUF = 1000; +SR = 10; +% Here you put the size of your data +SizeOfBloc = [SR*84 SR*143 1000]; + +% Here you put the size of pixel in the prefered unit. It can be um, mm, m, wavelength, or an arbitrary scale. +ScaleOfPixel = [1 1]; % [pixel_size_z, pixel_size_x] + +% The imaging frame rate is required for velocity calculation, and temporal filtering. +framerate = 500; % imaging framerate in [Hz] + +% Number of blocs to process +Nbuffers = 5; % number of bloc to process (used in the parfor) + +% If pixel sizes are in wavelength, lambda must be provided for a velocity maps in mm/s, +lambda = 0.09856; + +%% ULM parameters +% this script can be run using different scaling, it can be wavelength, pixel size, mm, um. +% In this example, input pixel are isotropic and equal to lambda (pixelPitch_x = pixelPitch_y = lambda) +% All size defined later are expressed in lambda + +res = SR; + +ULM = struct('numberOfParticles', 120,... % Number of particles per frame. (30-100) + 'size',[SizeOfBloc(1) SizeOfBloc(2) SizeOfBloc(3)],... % size of input data [nb_pixel_z nb_pixel_x nb_frame_per_bloc] + 'scale',[ScaleOfPixel 1/framerate],...% Scale [z x dt], size of pixel in the scaling unit. (here, pixsize = 1*lambda) + 'res',res,... % Resolution factor. Typically 10 for final image rendering at lambda/10. + 'SVD_cutoff',[10 SizeOfBloc(3)],... % SVD filtering, to be adapted to your clutter/SNR levels + 'max_linking_distance',2,... % Maximum linking distance between two frames to reject pairing, in pixels units (UF.scale(1)). (2-4 pixel). + 'min_length', 15,... % Minimum allowed length of the tracks in time. (5-20 frames) + 'fwhm',[1 1]*3,... % Size [pixel] of the mask for localization. (3x3 for pixel at lambda, 5x5 at lambda/2). [fmwhz fmwhx] + 'max_gap_closing', 0,... % Allowed gap in microbubbles' pairing. (if you want to skip frames 0) + 'interp_factor',1/res,... % Interpfactor (decimation of tracks) + 'LocMethod','Radial'... % Select localization algorithm (WA,Interp,Radial,CurveFitting,NoLocalization) + ); +% ULM.ButterCuttofFreq = [50 249]; % Cutoff frequency (Hz) for additional filter. Typically [20 300] at 1kHz. +ULM.parameters.NLocalMax = 3; % Safeguard on the number of maxLocal in the fwhm*fwhm grid (3 for fwhm=3, 7 for fwhm=5) +% [but_b,but_a] = butter(2,ULM.ButterCuttofFreq/(framerate/2),'bandpass'); +ULM.lambda = lambda; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Load data and localize microbubbles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Track_tot = {}; +fprintf('--- ULM PROCESSING --- \n\n');t1=tic; +SNR = '_20dB' +for hhh = 1:Nbuffers + fprintf('Processing bloc %d/%d\n',hhh,Nbuffers); + + % Load IQ data + tmp = load(strcat(pwd,'/SilicoData',SNR,'/SilicoAug',int2str(hhh),'.mat'),'SilicoAug'); + + % Filtering of IQ to remove clutter (optional) + IQ_filt = SVDfilter(tmp.SilicoAug,ULM.SVD_cutoff);tmp = []; + + % Temporal filtering + % Removes MBs in the in Silico Data ! + % IQ_filt = filter(but_b,but_a,IQ_filt,[],3); %(optional) + % IQ_filt(~isfinite(IQ_filt))=0; + + % Detection and localization process (return a list of coordinates in pixel) + [MatTracking] = ULM_localization2D(abs(IQ_filt),ULM); IQ_filt=[]; + % Convert pixel into isogrid (pixel are not necessary isometric); + MatTracking(:,2:3) = (MatTracking(:,2:3) - [1 1]).*ULM.scale(1:2); + + save(strcat(pwd,'/SilicoData',SNR,'/Result/RSMatTrack',int2str(hhh),'.mat'),'MatTracking') + + % Tracking algorithm (list of tracks) + Track_tot_i = ULM_tracking2D(MatTracking,ULM); + + % Saving part: + %--- if for, you can save tracks at each loop to avoid RAM out of memory + % save([workingdir filesep filename '_tracks' num2str(hhh,'%.3d') '.mat'],'Track_tot_i','ULM') % + %--- if parfor you can cat all tracks in a huge cells variable + Track_tot{hhh} = Track_tot_i; + Track_tot_i={};MatTracking = []; +end +Track_tot = cat(1,Track_tot{:}); +Tend = toc(t1);disp('Done') +fprintf('ULM done in %d hours %.1f minutes.\n', floor(Tend/60/60), rem(Tend/60,60)); + + +%% Create individual variable to save using v6 version. +% By cutting Tracks into different variables small than 2GB, the save v6 is faster than save v7.3 +% CutTracks = round(linspace(1,numel(Track_tot),4)); +% Track_tot_1 = Track_tot(CutTracks(1):CutTracks(2)-1); +% Track_tot_2 = Track_tot(CutTracks(2):CutTracks(3)-1); +% Track_tot_3 = Track_tot(CutTracks(3):end); +% save(['example_tracks.mat'],'Track_tot_1','Track_tot_2','Track_tot_3','Tend','ULM','-v6') +% clear Track_tot_1 Track_tot_2 Track_tot_3 + +% load([workingdir filesep filename 'example_tracks.mat']) +% Track_tot = cat(1,Track_tot_1,Track_tot_2,Track_tot_3);clear Track_tot_1 Track_tot_2 Track_tot_3 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Create MatOut %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% create the MatOut density with interpolated tracks for visual analysis, and with non interpolated tracks for aliasing index calculation. +% Define the size of SRpixel for displaying (default 10) +ULM.SRscale = ULM.scale(1)/ULM.res; +ULM.SRsize = round(ULM.size(1:2).*ULM.scale(1:2)/ULM.SRscale); % + +% Convert tracks into SRpixel +Track_matout = cellfun(@(x) (x(:,[1 2 3 4])+[1 1 0 0]*1)./[ULM.SRscale ULM.SRscale 1 1],Track_tot,'UniformOutput',0); +llz = [0:ULM.SRsize(1)]*ULM.SRscale;llx = [0:ULM.SRsize(2)]*ULM.SRscale; + +%% Accumulate tracks on the final MatOut grid. +fprintf('--- CREATING MATOUTS --- \n\n') +MatOut = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1); %pos in superpix [z x] +% MatOut_zdir = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1,'mode','2D_vel_z'); %pos in superpix [z x] +% MatOut_vel = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1,'mode','2D_velnorm'); %pos in superpix [z x] +% MatOut_vel = MatOut_vel*ULM.lambda; % Convert into [mm/s] + +save(strcat(pwd,'/SilicoData',SNR,'/Result/RSULMSilico10.mat'),'MatOut','ULM','lambda','llx','llz') %,'MatOut_zdir','MatOut_vel', + +MatOutGT = load('MatOutGT.mat', 'MatOut').MatOut; + +crossCorrMatrix = normxcorr2(MatOut,MatOutGT); +[maxCorr, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); + +xOffset = xPeak - size(MatOutGT, 2); +yOffset = yPeak - size(MatOutGT, 1); + +MatOut = circshift(MatOut,[yOffset xOffset]); + +figure('WindowState', 'maximized'); +imagesc(log(MatOutGT+0.5)); +axis image +colormap hot +brighten(0.2); +set(gca,'XColor', 'none','YColor','none') +set(gca, 'color', 'none'); + +figure('WindowState', 'maximized'); +imagesc(log(MatOut+0.5)); +axis image +colormap hot +brighten(0.2); +set(gca,'XColor', 'none','YColor','none') +set(gca, 'color', 'none'); diff --git a/AugmentedInSilico/Fig3c_silicoaug_SRPCA.m b/AugmentedInSilico/Fig3c_silicoaug_SRPCA.m new file mode 100644 index 0000000000000000000000000000000000000000..0e743022ffd1445707cf38da5875077f6e9df351 --- /dev/null +++ b/AugmentedInSilico/Fig3c_silicoaug_SRPCA.m @@ -0,0 +1,173 @@ +%% silicoaug_SRPCA.m +% +% +% Based on the code created by Arthur Chavignon +% Adapted for the Augmented in Silico Data. +% +% DATE 2020.12.17 - VERSION 1.1 +% AUTHORS: Arthur Chavignon, Baptiste Heiles, Vincent Hingot. CNRS, Sorbonne Universite, INSERM. +% Laboratoire d'Imagerie Biomedicale, Team PPM. 15 rue de l'Ecole de Medecine, 75006, Paris +% Code Available under Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (see https://creativecommons.org/licenses/by-nc-sa/4.0/) +% ACADEMIC REFERENCES TO BE CITED +% Details of the code in the article by Heiles, Chavignon, Hingot, Lopez, Teston and Couture. +% Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy, Nature Biomedical Engineering, 2021. +% General description of super-resolution in: Couture et al., Ultrasound localization microscopy and super-resolution: A state of the art, IEEE UFFC 2018 + +clear +%cd ./.. +%cd('PALA_silico'); + +PALA_folder = fileparts(pwd); % location of the PALA folder +% DEFINE THE ADDONS DIRECTORY ON YOUR COMPUTER +addpath(genpath(fullfile(PALA_folder, 'PALA_addons'))); + +%% Data parameters +% A few parameters must be provided by the user depending of your input images (size of pixel, wavelength) +% These parameters will be copied and used later during the creation of the ULM structure. +% in this example, UF.TwFreq = 15MHz, UF.FrameRateUF = 1000Hz; +UF.TwFreq = 14.8; +UF.FrameRateUF = 1000; +SR = 10; +% Here you put the size of your data +SizeOfBloc = [SR*84 SR*143 1000]; + +% Here you put the size of pixel in the prefered unit. It can be um, mm, m, wavelength, or an arbitrary scale. +ScaleOfPixel = [1 1]; % [pixel_size_z, pixel_size_x] + +% The imaging frame rate is required for velocity calculation, and temporal filtering. +framerate = 500; % imaging framerate in [Hz] + +% Number of blocs to process +Nbuffers = 5; % number of bloc to process (used in the parfor) + +% If pixel sizes are in wavelength, lambda must be provided for a velocity maps in mm/s, +lambda = 0.09856; + +%% ULM parameters +% this script can be run using different scaling, it can be wavelength, pixel size, mm, um. +% In this example, input pixel are isotropic and equal to lambda (pixelPitch_x = pixelPitch_y = lambda) +% All size defined later are expressed in lambda + +res = SR; + +ULM = struct('numberOfParticles', 120,... % Number of particles per frame. (30-100) + 'size',[SizeOfBloc(1) SizeOfBloc(2) SizeOfBloc(3)],... % size of input data [nb_pixel_z nb_pixel_x nb_frame_per_bloc] + 'scale',[ScaleOfPixel 1/framerate],...% Scale [z x dt], size of pixel in the scaling unit. (here, pixsize = 1*lambda) + 'res',res,... % Resolution factor. Typically 10 for final image rendering at lambda/10. + 'SVD_cutoff',[10 SizeOfBloc(3)],... % SVD filtering, to be adapted to your clutter/SNR levels + 'max_linking_distance',2,... % Maximum linking distance between two frames to reject pairing, in pixels units (UF.scale(1)). (2-4 pixel). + 'min_length', 15,... % Minimum allowed length of the tracks in time. (5-20 frames) + 'fwhm',[1 1]*3,... % Size [pixel] of the mask for localization. (3x3 for pixel at lambda, 5x5 at lambda/2). [fmwhz fmwhx] + 'max_gap_closing', 0,... % Allowed gap in microbubbles' pairing. (if you want to skip frames 0) + 'interp_factor',1/res,... % Interpfactor (decimation of tracks) + 'LocMethod','Radial'... % Select localization algorithm (WA,Interp,Radial,CurveFitting,NoLocalization) + ); +% ULM.ButterCuttofFreq = [50 249]; % Cutoff frequency (Hz) for additional filter. Typically [20 300] at 1kHz. +ULM.parameters.NLocalMax = 3; % Safeguard on the number of maxLocal in the fwhm*fwhm grid (3 for fwhm=3, 7 for fwhm=5) +% [but_b,but_a] = butter(2,ULM.ButterCuttofFreq/(framerate/2),'bandpass'); +ULM.lambda = lambda; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Load data and localize microbubbles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +load('RFbUpscaled.mat', 'RFb') +SNR = '_20dB' + +t1=tic;Track_tot = {}; +fprintf('--- ULM PROCESSING --- \n\n');t1=tic; + + +for hhh = 1:Nbuffers + fprintf('Processing bloc %d/%d\n',hhh,Nbuffers); + + % Load IQ data + tmp = load(strcat(pwd,'/SilicoData',SNR,'/SilicoAug',int2str(hhh),'.mat'),'SilicoAug'); + + tmp.SilicoAug = imresize(tmp.SilicoAug,SR,"bicubic"); + + % Temporal filtering + % Removes MBs in the in Silico Data ! + % tmp.SilicoAug = filter(but_b,but_a,tmp.SilicoAug,[],3); %(optional) + % tmp.SilicoAug(~isfinite(tmp.SilicoAug))=0; + + IQ = tmp.SilicoAug/max(abs(tmp.SilicoAug(:)));tmp = []; + IQ_filt = SRPCA(IQ, RFb, 1, 20); IQ = []; + + % Detection and localization process (return a list of coordinates in pixel) + [MatTracking] = ULM_localization2DSR(abs(IQ_filt),ULM); IQ_filt=[]; + MatTracking(:,2:3) = MatTracking(:,2:3)/SR; + % Convert pixel into isogrid (pixel are not necessary isometric); + MatTracking(:,2:3) = (MatTracking(:,2:3) - [1 1]).*ULM.scale(1:2); + %save(strcat(pwd,'/SilicoData',SNR,'/Result/SRPCAMatTrack',int2str(hhh),'.mat'),'MatTracking') + + % Tracking algorithm (list of tracks) + Track_tot_i = ULM_tracking2D(MatTracking,ULM); + + % Saving part: + %--- if for, you can save tracks at each loop to avoid RAM out of memory + % save([workingdir filesep filename '_tracks' num2str(hhh,'%.3d') '.mat'],'Track_tot_i','ULM') % + %--- if parfor you can cat all tracks in a huge cells variable + Track_tot{hhh} = Track_tot_i; + Track_tot_i={};MatTracking = []; +end +Track_tot = cat(1,Track_tot{:}); +Tend = toc(t1);disp('Done') +fprintf('ULM done in %d hours %.1f minutes.\n', floor(Tend/60/60), rem(Tend/60,60)); + + +%% Create individual variable to save using v6 version. +% By cutting Tracks into different variables small than 2GB, the save v6 is faster than save v7.3 +% CutTracks = round(linspace(1,numel(Track_tot),4)); +% Track_tot_1 = Track_tot(CutTracks(1):CutTracks(2)-1); +% Track_tot_2 = Track_tot(CutTracks(2):CutTracks(3)-1); +% Track_tot_3 = Track_tot(CutTracks(3):end); +% save(['example_tracks.mat'],'Track_tot_1','Track_tot_2','Track_tot_3','Tend','ULM','-v6') +% clear Track_tot_1 Track_tot_2 Track_tot_3 + +% load([workingdir filesep filename 'example_tracks.mat']) +% Track_tot = cat(1,Track_tot_1,Track_tot_2,Track_tot_3);clear Track_tot_1 Track_tot_2 Track_tot_3 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Create MatOut %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% create the MatOut density with interpolated tracks for visual analysis, and with non interpolated tracks for aliasing index calculation. +% Define the size of SRpixel for displaying (default 10) +ULM.SRscale = ULM.scale(1)/ULM.res; +ULM.SRsize = round(ULM.size(1:2).*ULM.scale(1:2)); %/ULM.SRscale + +% Convert tracks into SRpixel +Track_matout = cellfun(@(x) (x(:,[1 2 3 4])+[1 1 0 0]*1)./[ULM.SRscale ULM.SRscale 1 1],Track_tot,'UniformOutput',0); +llz = [0:ULM.SRsize(1)]*ULM.SRscale;llx = [0:ULM.SRsize(2)]*ULM.SRscale; + +%% Accumulate tracks on the final MatOut grid. +fprintf('--- CREATING MATOUTS --- \n\n') +MatOut = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1); %pos in superpix [z x] +% MatOut_zdir = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1,'mode','2D_vel_z'); %pos in superpix [z x] +% MatOut_vel = ULM_Track2MatOut(Track_matout,ULM.SRsize+[1 1]*1,'mode','2D_velnorm'); %pos in superpix [z x] +% MatOut_vel = MatOut_vel*ULM.lambda; % Convert into [mm/s] + +save(strcat(pwd,'/SilicoData',SNR,'/Result/SRPCASilico10.mat'),'MatOut','ULM','lambda','llx','llz') %,'MatOut_zdir','MatOut_vel', + +MatOutGT = load('MatOutGT.mat', 'MatOut').MatOut; + +crossCorrMatrix = normxcorr2(MatOut,MatOutGT); +[maxCorr, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); + +xOffset = xPeak - size(MatOutGT, 2); +yOffset = yPeak - size(MatOutGT, 1); + +MatOut = circshift(MatOut,[yOffset xOffset]); + +figure('WindowState', 'maximized'); +imagesc(log(MatOutGT+0.5)); +axis image +colormap hot +brighten(0.2); +set(gca,'XColor', 'none','YColor','none') +set(gca, 'color', 'none'); + +figure('WindowState', 'maximized'); +imagesc(log(MatOut+0.5)); +axis image +colormap hot +brighten(0.2); +set(gca,'XColor', 'none','YColor','none') +set(gca, 'color', 'none'); diff --git a/AugmentedInSilico/MatOutGT.mat b/AugmentedInSilico/MatOutGT.mat new file mode 100644 index 0000000000000000000000000000000000000000..bb4b711b2f97ee1ab54c16219e14708631780f89 Binary files /dev/null and b/AugmentedInSilico/MatOutGT.mat differ diff --git a/AugmentedInSilico/RFbUpscaled.mat b/AugmentedInSilico/RFbUpscaled.mat new file mode 100644 index 0000000000000000000000000000000000000000..d1448b494464b4f89128fbd4033fded8fb3b6a59 Binary files /dev/null and b/AugmentedInSilico/RFbUpscaled.mat differ diff --git a/AugmentedInSilico/SRPCA.m b/AugmentedInSilico/SRPCA.m new file mode 100644 index 0000000000000000000000000000000000000000..8f26557a62404520869a99ebb3ca752be2d4ed3d --- /dev/null +++ b/AugmentedInSilico/SRPCA.m @@ -0,0 +1,119 @@ +function x = SRPCA8(S,H,lambda,mu,rho,max_iter) + +% Computational Super-Resolved RPCA. +% Input +% - S is a 2d+t data matrix (of size m x n x p) to be decomposed +% - H: Point Spread Function (PSF) +% - lambda: Sparse regularization parameter, default = 1/sqrt(max(n*m, p)) +% - mu: Low-Rank regularization parameter, default = 10 +% - rho: the augmented lagrangian parameter (convergence), default = 1e-3 +% - max_iter: maximum number of iterations, default = 20 +% +% Ouput +% - x is the high resolution blood +% - T (not returned) is the tissue + + [m, n, p] = size(S); + unobserved = isnan(S); + S(unobserved) = 0; + + if nargin < 2 + H = 1; + else + H = H / sum(abs(H(:))); + [Mh, Nh] = size(H); + center = round([Mh, Nh] / 2); + H = fft2(circshift(padarray(H, [m - Mh, n - Nh], 'post'), 1 - center)); + end + + + if nargin < 3 + lambda = 1/sqrt(m*n); + else + lambda = lambda/sqrt(m*n); + end + + if nargin < 4 + mu = 10; + end + + if nargin < 5 + rho = 1e-3; + end + + if nargin < 6 + max_iter = 20; + end + + + rho_max = 1; + eps = 1.5; + + %% initial solution + T = single(zeros(m, n, p)); + y = single(zeros(m, n, p)); + x = single(zeros(m, n, p)); + Z = single(zeros(m, n, p)); + N = single(zeros(m, n, p)); + W = single(zeros(m, n, p)); + Hx = single(zeros(m, n, p)); + + Dt = conj(H); + DD = abs(H).^2; + + for iter = (1:max_iter-1) + y = (S-Hx + rho*(T - W))./(1+rho); + + x = ifft2(fft2(ifft2(Dt.*fft2(S-y)) + rho*(Z - N))./(DD + rho)); + + Hx = ifft2(H.*fft2(x)); + + T = fastDo(y + W, mu); + Z = So(lambda, x + N); + + W = W + y - T; %gamma2 + N = N + x - Z; %gamma1 + + rho = min(rho*eps, rho_max); + + end + + y = (S-Hx + rho*(T - W))./(1+rho); + + x = ifft2(fft2(ifft2(Dt.*fft2(S-y)) + rho*(Z - N))./(DD + rho)); + + +end + +function r = So(tau, S) + % shrinkage operator + r = sign(S) .* max(abs(S) - tau, 0); +end + +% Matlab SVD implementation is slow +% function r = Do(tau, S) +% % shrinkage operator for singular values +% [U, D, V] = svd(S, 'econ'); +% r = U*So(tau, D)*V'; +% end + +function SX = fastDo(IQ,tau) + + % Fast SVD implementation based on the code of + % Vipin Vijayan, Fast SVD and PCA, https://www.mathworks.com/matlabcentral/fileexchange/47132-fast-svd-and-pca + + initsize = size(IQ); + X = reshape(IQ,prod(initsize(1:end-1)),initsize(end)); + + C = X'*X; + [V,D] = eig(double(C)); + clear C; + U = X*V; % convert evecs from X'*X to X*X'. the evals are the same. + s = sqrt(abs(diag(D))); + U = U./s'; + s = max(s - tau, 0); + S = diag(s); + r=U*S*V'; + SX = reshape(r,initsize); +end + diff --git a/AugmentedInSilico/SilicoAugMetrics.m b/AugmentedInSilico/SilicoAugMetrics.m new file mode 100644 index 0000000000000000000000000000000000000000..1b5973dfa2e547eb29f34f63a96b76a084665422 --- /dev/null +++ b/AugmentedInSilico/SilicoAugMetrics.m @@ -0,0 +1,140 @@ +clear +numblocks = 5; +NoiseLevels = {'_40dB','_20dB','_10dB','_0dB'}; +SR = 2; +MatOutGT = load('MatOutGT.mat', 'MatOut').MatOut; + + +%% 'SRPCA' + +for ijk = 1:length(NoiseLevels) +SNR = NoiseLevels{ijk}; +formatSpec = ['SRPCA' SNR ' prec %4.2f, sens %4.2f, jacc %4.2f, ssim %4.2f, MSE %4.2f \n']; + +TP = zeros(numblocks,1); +FP = zeros(numblocks,1); +FN = zeros(numblocks,1); + +for hhh = 1:numblocks + + load(strcat('BubbleGT/BubbleGT',int2str(SR),'_Silico',int2str(hhh),'.mat')) + sumBubbleGT=sum(BubbleGT,3); + load(strcat('SilicoData',SNR,'/Result/SRPCAMatTrack',int2str(hhh),'.mat')) + +Est = zeros(size(BubbleGT)); + +for kk = 1:length(MatTracking) + Est(max(1, round( MatTracking(kk,2)*SR)),max(1, round( MatTracking(kk,3)*SR)),MatTracking(kk,4))=Est(max(1, round( MatTracking(kk,2)*SR)),max(1, round( MatTracking(kk,3)*SR)),MatTracking(kk,4))+1; +end + +crossCorrMatrix = normxcorr2(sum(Est,3),sumBubbleGT); +[~, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); +yOffset = yPeak - size(sumBubbleGT, 1); +xOffset = xPeak - size(sumBubbleGT, 2); + +Est = circshift(Est,[yOffset xOffset 0]); + +temp = max(BubbleGT - Est, 0) ; +NumBullesGT = sum(BubbleGT(:)); +NumUndetected = sum(temp(:)); + +FN(hhh) = NumUndetected; + +TP(hhh)= NumBullesGT - NumUndetected; + +temp = min(BubbleGT - Est, 0) ; + +FP(hhh)= - sum(temp(:)); + + +end + + + J_p = mean(TP./(FP+TP)*100); % precision + J_r = mean(TP./(FN+TP)*100); % sensitivity + J_ac = mean(TP./(FP+FN+TP)*100); % Jaccard + +load(['SilicoData',SNR,'/Result/SRPCASilico10.mat'],'MatOut') + +crossCorrMatrix = normxcorr2(MatOut,MatOutGT); +[maxCorr, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); + +xOffset = xPeak - size(MatOutGT, 2); +yOffset = yPeak - size(MatOutGT, 1); + +MatOut = circshift(MatOut,[yOffset xOffset]); + +ssimval = ssim(MatOut,MatOutGT)*100; +cNRMSE = US_ADM_calc_PSNR(MatOutGT,MatOut).NRMSE*100; +fprintf(formatSpec,J_p,J_r,J_ac,ssimval,cNRMSE) +end + +%% 'RS ULM' +for ijk = 1:length(NoiseLevels) +SNR = NoiseLevels{ijk}; + +formatSpec = ['RSULM' SNR ' prec %4.2f, sens %4.2f, jacc %4.2f, ssim %4.2f, MSE %4.2f \n']; + + +TP = zeros(numblocks,1); +FP = zeros(numblocks,1); +FN = zeros(numblocks,1); + +for hhh = 1:numblocks + + load(strcat('BubbleGT/BubbleGT',int2str(SR),'_Silico',int2str(hhh),'.mat')) + sumBubbleGT=sum(BubbleGT,3); + load(strcat('SilicoData',SNR,'/Result/RSMatTrack',int2str(hhh),'.mat')) + +Est = zeros(size(BubbleGT)); + +for kk = 1:length(MatTracking) + Est(max(1, round( MatTracking(kk,2)*SR)),max(1, round( MatTracking(kk,3)*SR)),MatTracking(kk,4))=Est(max(1, round( MatTracking(kk,2)*SR)),max(1, round( MatTracking(kk,3)*SR)),MatTracking(kk,4))+1; +end + +crossCorrMatrix = normxcorr2(sum(Est,3),sumBubbleGT); +[~, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); +yOffset = yPeak - size(sumBubbleGT, 1); +xOffset = xPeak - size(sumBubbleGT, 2); + +Est = circshift(Est,[yOffset xOffset 0]); + +temp = max(BubbleGT - Est, 0) ; +NumBullesGT = sum(BubbleGT(:)); +NumUndetected = sum(temp(:)); + +FN(hhh) = NumUndetected; + +TP(hhh)= NumBullesGT - NumUndetected; + +temp = min(BubbleGT - Est, 0) ; + +FP(hhh)= - sum(temp(:)); + + +end + + + J_p = mean(TP./(FP+TP)*100); % precision + J_r = mean(TP./(FN+TP)*100); % sensitivity + J_ac = mean(TP./(FP+FN+TP)*100); % Jaccard + +load(['SilicoData',SNR,'/Result/RSULMSilico10.mat'],'MatOut') + +crossCorrMatrix = normxcorr2(MatOut,MatOutGT); +[maxCorr, maxIndex] = max(abs(crossCorrMatrix(:))); +[yPeak, xPeak] = ind2sub(size(crossCorrMatrix), maxIndex); + +xOffset = xPeak - size(MatOutGT, 2); +yOffset = yPeak - size(MatOutGT, 1); + +MatOut = circshift(MatOut,[yOffset xOffset]); + +ssimval = ssim(MatOut,MatOutGT)*100; +cNRMSE = US_ADM_calc_PSNR(MatOutGT,MatOut).NRMSE*100; + +fprintf(formatSpec,J_p,J_r,J_ac,ssimval,cNRMSE) +end \ No newline at end of file diff --git a/AugmentedInSilico/ULM_localization2DSR.m b/AugmentedInSilico/ULM_localization2DSR.m new file mode 100644 index 0000000000000000000000000000000000000000..6059bcd7cd9edcc38b2c45c855ae640f23e177d0 --- /dev/null +++ b/AugmentedInSilico/ULM_localization2DSR.m @@ -0,0 +1,177 @@ +function MatTracking = ULM_localization2DSR(MatIn,ULM) +% Adapted from ULM_localization2D for regional max detection only without using any localization method +%% function MatTracking = ULM_localization2D(MatIn,ULM) +% This function performs the detection, selection and sub-pixel localization of bubbles on +% a list of input images (MatIn). +% +% - The detection step is performed with the imregionalmax function. It returns the list of local maxima. +% - The selection step consists of sorting intensities, in each frames, and keeping the highest maxima. +% - The localization steps consists of applying a sub-wavelength localization kernel +% (weighted average, interpolation, radial symmetry...) to a cropped image centered on a +% local maxima (a 5x5 or 3x3 large image). Localization kernel are discussed in the cited article. +% +% This function can be easily adapt by anyone who wants to try a new localization kernel, +% keeping all the framework unchanged. +% +% INPUTS: +% - MatIn is the sequence containing all the images +% - ULM structure must contains all parameters for localization method +% - NLocalMax : number of local maxima +% - LocMethod : localisation mehtod {'wa','interp','radial','curvefitting','nolocalization'}; +% - InterpMethod : {'bicubic','lanczos3','spline'} (methods avaliable for imresize) +% - numberOfParticles is an estimation of the number of particles per image +% - fwhm is the fwhm of a bubble (usually 3 if pixel at \lambda, 5 if pixel at \lambda/2) +% OUTPUT: +% - MatTracking is the table that stores the particles values and position in [pixel] +% MatInReduced is the input matrix within a zeros frame of FWHM/2 +% +% This function was created by Baptiste Heiles 07/04/17, last modifications Arthur Chavignon, 18/11/19 +% +% DATE 2020.07.22 - VERSION 1.1 +% AUTHORS: Baptiste Heiles, Arthur Chavignon, Vincent Hingot. CNRS, Sorbonne Universite, INSERM. +% Laboratoire d'Imagerie Biomedicale, Team PPM. 15 rue de l'Ecole de Medecine, 75006, Paris +% Code Available under Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International (see https://creativecommons.org/licenses/by-nc-sa/4.0/) +% ACADEMIC REFERENCES TO BE CITED +% Details of the code in the article by Heiles, Chavignon, Hingot, Lopez, Teston and Couture. +% Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy, Nature Biomedical Engineering, 2021. +% General description of super-resolution in: Couture et al., Ultrasound localization microscopy and super-resolution: A state of the art, IEEE UFFC 2018 + +%% Get input data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +fwhmz = ULM.fwhm(2); +fwhmx = ULM.fwhm(1); + +% Building a vector from -FWHM to FWHM, this vector will be used for the mask's shifting +vectfwhmz = -1*round(fwhmz/2):round(fwhmz/2); +vectfwhmx = -1*round(fwhmx/2):round(fwhmx/2); + +[height,width,numberOfFrames]=size(MatIn);% Get sizes of the Matrix, height denotes number of rows/depth of imaging, width denotes number of lines/width of imaging, +% numberOfFrames denotes the number of elements in the third dimension/number of Frames +MatIn = abs(MatIn);% Make sure you work with the intensity matrix +info = whos('MatIn');typename = info.class; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Initialize structures +% if fields are missing, they will be set with default values. +if ~isfield(ULM,'LocMethod'),ULM.LocMethod = 'radial';end + +if ~isfield(ULM,'parameters')% Create an empty structure for parameters hosting + ULM.parameters = struct(); +end + +if strcmp(ULM.LocMethod,'interp') + if ~isfield(ULM.parameters,'InterpMethod') + ULM.parameters.InterpMethod = 'spline'; + end + if sum(strcmp(ULM.parameters.InterpMethod,{'bilinear','bicubic'})) + warning('Faster but pixelated, Weighted Average will be faster and smoother.') + end +end + +if ~isfield(ULM.parameters,'NLocalMax') + if fwhmz==3,ULM.parameters.NLocalMax = 2; + else,ULM.parameters.NLocalMax = 3; + end +end + +%% 1 PREPARE INTENSITY MATRIX +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% MATRIX CROPPING %% +% Creates smaller matrix MatInReduced to avoid boundaries, where microbubbles cannot be localized. We avoided padding because padding would +% result in erroneous localization in the boundaries. +MatInReduced = zeros(height,width,numberOfFrames,typename); +MatInReduced(1+round(fwhmz/2)+1:height-round(fwhmz/2)-1,1+round(fwhmx/2)+1:width-round(fwhmx/2)-1,:) = MatIn(1+round(fwhmz/2)+1:height-round(fwhmz/2)-1, 1+round(fwhmx/2)+1:width-round(fwhmx/2)-1,:); +[height,width,numberOfFrames] = size(MatInReduced); + + +%% 2 DETECTION AND SELECTION OF MICROBUBBLES +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% DETECTION OF LOCAL MAXIMA %% +% Concatenates MatInReduced into a 2D matrix with time and space in its rows and only space in columns to apply imregionalmax function +% The imregionalmax connectivity is set to default (8), to consider the 8 adjacent pixels in the horizontal(2), vertical(2), and diagonal directions(4). +% Generates an IntensityMatrix (3D) with only local maximal pixels with associated values. +Mat2D = permute(MatInReduced, [1,3,2]); %so that all the frames are in columns +Mat2D = reshape(Mat2D,height*numberOfFrames,width);% Concatenate Matrix +mask2D = imregionalmax(Mat2D); clear Mat2D % Perform imregionalmax +mask = reshape(mask2D,height,numberOfFrames,width);clear mask2D % reshape concatenated mask +mask = permute(mask,[1,3,2]); % so that we restore (z,x,t) table + +IntensityMatrix = MatInReduced.*mask; %Values of intensities at regional maxima + +% SELECTION OF MICROBUBBLES %% +% Only the first numberOfParticles highest local max will be kept for localization. +% Other local max will be considered as noise. +% Sort intensites in each frames, and store pixel coordinates +% At the end of this section, spatial and temporal coordinates microbubbles are +% stored into: index_mask_z, index_mask_x, index_numberOfFrames +[tempMatrix,~] = sort(reshape(IntensityMatrix,[],size(IntensityMatrix,3)),1,'descend'); + +% Remove the last kept intensity values to each frame. This means that you cannot fix an intensity threshold, +% we rely on number of particles. This is key for transparency/parallelization. +IntensityFinal = IntensityMatrix - ones(size(IntensityMatrix)) .* reshape(tempMatrix( ULM.numberOfParticles+1,:),[1 1 numberOfFrames]); +clear tempMatrix +% Construction of the final mask with only the kept microbubbles low resolved and their associated intensity +MaskFinal = (mask.*IntensityFinal)>0; +MaskFinal(isnan(MaskFinal))=0; +MaskFinal = (MaskFinal>0).*IntensityMatrix; + +% Preparing intensities and coordinates for further calculation of average, intensities etc... +index_mask = find(MaskFinal); +[index_mask_z,index_mask_x,index_numberOfFrames]=ind2sub([height, width, numberOfFrames], index_mask); +clear mask IntensityFinal MaskFinal IntensityMatrix +clear index_mask + + +%% 3 SUBWALENGTH LOCALIZATION OF MICROBUBBLES +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% LOCALIZATION OF MICROBUBBLES %% +% The position of each microbubble will be localized with a subwavelength precision. +% For example a microbubble at pixel coordinates [34 67] will be localized at [34.4 66.8] + +% Initialize variables averageXc, averageZc which are the super-resolved position of microbubbles +averageXc = nan(1,size(index_mask_z,1),typename); +averageZc = nan(1,size(index_mask_z,1),typename); + +for iscat=1:size(index_mask_z,1) + % For each microbubble, create a 2D intensity matrix of the Region of interest defined by fwhm + IntensityRoi = MatIn(index_mask_z(iscat)+vectfwhmz,index_mask_x(iscat)+vectfwhmx,index_numberOfFrames(iscat)); + + % NLocal max + % If there are too many localmax in the region of interest, the microbubble shape will be affected and the localization distorted. + % In that case, we set averageZc, averageXc to NaN value. + if nnz(imregionalmax(IntensityRoi))>ULM.parameters.NLocalMax + continue + end + + % Apply the localization method selected + % functions are detailed at the end of the code (excepted LocRadialSym which requires an additional function) + + % Store the final super-resolved position of the microbubble as its pixel position and an axial/lateral sub-pixel shift. + averageZc(iscat) = index_mask_z(iscat); + averageXc(iscat) = index_mask_x(iscat); + + % Additional safeguards + % sigma evaluates the size of the microbubble. If it appears to be too large, the microbubble can be removed (optional) + + + % If the final axial/lateral shift is higher that the fwhmz, + % localization has diverged and the microbubble is ignored. +end +keepIndex = ~isnan(averageXc); + +ind = sub2ind([height,width,numberOfFrames],index_mask_z(keepIndex),index_mask_x(keepIndex),index_numberOfFrames(keepIndex)); +clear index_mask_z index_mask_x IntensityRoi + +%% BUILD MATTRACKING %% +% Creating the table which stores the high resolved microbubbles coordinates and the density value +MatTracking = zeros(nnz(keepIndex),4,typename); + +MatTracking(:,1) = MatInReduced(ind); % Initial intensity of the microbubble +MatTracking(:,2) = averageZc(keepIndex); % Super-resolved axial coordinate +MatTracking(:,3) = averageXc(keepIndex); % Super-resolved lateral coordinate +MatTracking(:,4) = index_numberOfFrames(keepIndex); % Frame number of the microbubble +clear averageXc averageZc index_numberOfFrames MatInReduced +end +