Skip to content
Snippets Groups Projects
Commit 5674b984 authored by vpustova's avatar vpustova
Browse files

Silico Code

parent 949357b5
Branches main
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 776 additions and 0 deletions
AugmentedInSilico/Fig3b_silicoaug_RSULM.m 0 → 100644
+ 167
− 0
View file @ 5674b984
%% 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');
AugmentedInSilico/Fig3c_silicoaug_SRPCA.m 0 → 100644
+ 173
− 0
View file @ 5674b984
%% 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');
AugmentedInSilico/MatOutGT.mat 0 → 100644
+ 0
− 0
View file @ 5674b984
File added
AugmentedInSilico/RFbUpscaled.mat 0 → 100644
+ 0
− 0
View file @ 5674b984
File added
AugmentedInSilico/SRPCA.m 0 → 100644
+ 119
− 0
View file @ 5674b984
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
AugmentedInSilico/SilicoAugMetrics.m 0 → 100644
+ 140
− 0
View file @ 5674b984
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
AugmentedInSilico/ULM_localization2DSR.m 0 → 100644
+ 177
− 0
View file @ 5674b984
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment