diff --git a/InVivoTours/Fig6ac_InVivoTours_RSULM.m b/InVivoTours/Fig6ac_InVivoTours_RSULM.m
new file mode 100644
index 0000000000000000000000000000000000000000..53e0e70337d7a4f3b0886161582af75be7492378
--- /dev/null
+++ b/InVivoTours/Fig6ac_InVivoTours_RSULM.m
@@ -0,0 +1,297 @@
+%% InVivoTours_RSULM.m
+%
+%
+% Based on the code created by Arthur Chavignon
+% Adapted for the WKY in vivo 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 ./..
+
+PALA_folder = fileparts(pwd); % location of the PALA folder
+PALA_data_folder = fullfile(PALA_folder, 'PALA_data_folder');
+% DEFINE THE ADDONS DIRECTORY ON YOUR COMPUTER
+addpath(genpath(fullfile(PALA_folder, 'PALA_addons')));
+
+%% Adapt parameters to your data
+% 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;
+
+% Here you put the size of your data
+SizeOfBloc = [160 128 500];
+% 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]
+% In that example, the size of pixels is lambda x lambda. The localization process will be
+% performs in wavelength. For velocity rendering, velocities in [wavelength/s] will be converted in [mm/s].
+
+% The imaging frame rate is required for velocity calculation, and temporal filtering.
+framerate = 1000; % imaging framerate in [Hz]
+
+% Number of blocs to process
+Nbuffers = 80; % 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 = 1480/(UF.TwFreq*1e3);
+
+%% 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 = 10; % final ratio of localization rendering, it's approximately resolution factor of localization in scale(1) units.
+% for a pixel size of 100um, we can assume that ULM algorithm provides precision 10 (res) times
+% smaller than pixel size. Final rendering will be reconstructed on a 10x10um grid.
+
+ULM = struct('numberOfParticles', 90,... % 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',[5 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;
+
+for hhh = 1:Nbuffers
+
+ fprintf('Processing bloc %d/%d\n',hhh,Nbuffers);
+ % Load IQ data
+ tmp = load(strcat(PALA_data_folder,'/IQTours/IQTours_',int2str(hhh),'.mat'),'IQ');
+ tmp.IQ = tmp.IQ(5:164,:,:);
+
+ % Filtering of IQ to remove clutter (optional)
+ IQ_filt = SVDfilter(tmp.IQ,ULM.SVD_cutoff);tmp = [];
+
+ % Temporal filtering
+ 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);
+
+ % Tracking algorithm (list of tracks)
+ Track_tot_i = ULM_tracking2D(MatTracking,ULM);
+
+
+ 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,'/Result/RSULM_MatOutTours10.mat'),'MatOut','MatOut_zdir','MatOut_vel','ULM','lambda','llx','llz')
+
+%% Figures
+ULM.SRscale=0.625/ULM.res; %for velocity maps in mm/s and scale-bars
+maxSRPCA = 64;
+
+IntPower = 1/2.3;
+
+cmin = min(MatOut(:).^IntPower);
+cmax = max(MatOut(:).^IntPower);
+numTicks = 7; % Adjust as needed
+
+h1=figure('WindowState', 'maximized');
+im=imagesc((MatOut).^IntPower);
+axis image
+colormap hot
+clbar = colorbar;
+set(clbar, 'YTick', linspace(cmin, cmax, numTicks));
+
+clbar.Label.String = 'Normalized localization density';
+%clbar.TickLabels = clbar.Ticks.^2.3; %Not normalized
+clbar.TickLabels = round(clbar.Ticks.^2.3 ./ maxSRPCA, 2);
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:2],60+[0:BarWidth])=max(clim);
+text(60+BarWidth/2, size(MatOut,1)-50+3 +15, '1 mm', ...
+ 'Color', 'w', 'FontSize', 12, 'FontWeight', 'bold', ...
+ 'HorizontalAlignment', 'center');
+
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2)+clbar.Ticks(end)/2, x(3)];
+
+%%
+IntPower = 1/2.3;
+
+OriginalData = MatOut.^IntPower.*sign(imgaussfilt(MatOut_zdir,.8));
+cmin = min(OriginalData(:));
+cmax = max(OriginalData(:));
+clear OriginalData;
+numTicks = 9; % Adjust as needed
+
+figure('WindowState', 'maximized');
+velColormap = cat(1,flip(flip(hot(128),1),2),hot(128));
+velColormap = velColormap(5:end-5,:);
+im=imagesc(llx,llz,(MatOut).^IntPower.*sign(imgaussfilt(MatOut_zdir,.8)))
+im.CData = im.CData - sign(im.CData)/2;axis image
+colormap(gca,velColormap)
+
+
+clbar = colorbar;
+set(clbar, 'YTick', linspace(cmin*0.90, cmax*0.90, numTicks));
+
+clbar.Label.String = 'Normalized localization density';
+clbar.Ticks = clbar.Ticks;%customTicks;
+TempTickLabels = round(sign(clbar.Ticks).*abs(clbar.Ticks/0.90).^2.3);
+clbar.TickLabels = round(TempTickLabels ./ maxSRPCA,2);
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+caxis([-1 1]*max(caxis)*.75)
+
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:3],60+[0:BarWidth])=max(caxis);
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2), x(3)];
+
+%%
+
+
+%vmax_disp = ceil(quantile(MatOut_vel(abs(MatOut_vel)>0),.98)/10)*10;
+vmax_disp = 80;
+
+h1=figure('WindowState', 'maximized'); clf,set(gcf,'Position',[652 393 941 585]);
+clbsize = [180,50];
+Mvel_rgb = MatOut_vel/vmax_disp; % normalization
+% Mvel_rgb(1:clbsize(1),1:clbsize(2)) = repmat(linspace(1,0,clbsize(1))',1,clbsize(2)); % add velocity colorbar
+Mvel_rgb = Mvel_rgb.^(1/2);Mvel_rgb(Mvel_rgb>1)=1;
+Mvel_rgb = imgaussfilt(Mvel_rgb,.5);
+Mvel_rgb = ind2rgb(round(Mvel_rgb*256),jet(256)); % convert ind into RGB
+IntPower = 1/2.3;
+MatShadow = MatOut;MatShadow = MatShadow./max(MatShadow(:)*.3);MatShadow(MatShadow>1)=1;
+%MatShadow(1:clbsize(1),1:clbsize(2))=repmat(linspace(0,1,clbsize(2)),clbsize(1),1);
+Mvel_rgb = Mvel_rgb.*(MatShadow.^IntPower);
+Mvel_rgb = brighten(Mvel_rgb,.4);
+im=imshow(Mvel_rgb,'XData',llx,'YData',llz);axis on
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:3],60+[0:BarWidth],1:3)=1;
+
+
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+colormap(jet);
+clbar = colorbar;
+clim([0 vmax_disp]);
+clbar.Label.String = 'Velocity magnitude (mm/s)';
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2)+clbar.Ticks(end)/2, x(3)];
+
+%% Metrics
+% Saturation Value
+% load(strcat(pwd,'/Result/RSULM_MatOutTours10.mat'),'MatOut')
+Saturation = nnz(MatOut()>0)/numel(MatOut());
+
+
+% CR
+% load(strcat(pwd,'/Result/RSULM_MatOutTours10.mat'),'MatOut')
+CR = ContrastRatio(MatOut);
+
+% Image for FRC estimation
+% load(strcat(pwd,'/Result/RSULM_MatOutTours10.mat'),'MatOut')
+
+figure(1);
+IntPower = 1/2.3;
+im=imagesc(llx,llz,MatOut.^IntPower);axis image
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+im.CData = max(im.CData-.5,0);
+
+colormap(gca,gray(128))
+caxis(caxis*.7) % add saturation in image
+ca = gca;
+
+WriteTif(im.CData,ca.Colormap,strcat(pwd,'/Result/RSULM_MatOutTours10.tif'),'caxis',caxis,'Overwrite',1)
+
+% Gridding index
+load(strcat(pwd,'/Result/RSULM_MatOutTours10_Gridding.mat'),'MatOut')
+Gridding_index = PALA_ComputeGriddingIndex({MatOut},1);
+Gridding = 1-Gridding_index/30; % [100 at 0; 0 at -30]
diff --git a/InVivoTours/Fig6db_InVivoTours_SRPCA.m b/InVivoTours/Fig6db_InVivoTours_SRPCA.m
new file mode 100644
index 0000000000000000000000000000000000000000..0574bd0e595e6ff45264b5cdbdca53ce2074ec82
--- /dev/null
+++ b/InVivoTours/Fig6db_InVivoTours_SRPCA.m
@@ -0,0 +1,301 @@
+%% InVivoTours_SRPCA.m
+%
+%
+% Based on the code created by Arthur Chavignon
+% Adapted for the WKY in vivo 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 ./..
+
+PALA_folder = fileparts(pwd); % location of the PALA folder
+PALA_data_folder = fullfile(PALA_folder, 'PALA_data_folder');
+% DEFINE THE ADDONS DIRECTORY ON YOUR COMPUTER
+addpath(genpath(fullfile(PALA_folder, 'PALA_addons')));
+
+
+%% Adapt parameters to your data
+% 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;
+
+SR = 10; %Super Resolution factor
+UF.TwFreq = 14.8;
+UF.FrameRateUF = 1000;
+
+% Here you put the size of your data
+SizeOfBloc = [SR*160 SR*128 500];
+
+
+% 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]
+% In that example, the size of pixels is lambda x lambda. The localization process will be
+% performs in wavelength. For velocity rendering, velocities in [wavelength/s] will be converted in [mm/s].
+
+% The imaging frame rate is required for velocity calculation, and temporal filtering.
+framerate = 1000; % imaging framerate in [Hz]
+
+% Number of blocs to process
+Nbuffers = 80; % 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 = 1480/(UF.TwFreq*1e3); % lambda = 1.002031144211239e-1
+
+%% 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; % final ratio of localization rendering, it's approximately resolution factor of localization in scale(1) units.
+% for a pixel size of 100um, we can assume that ULM algorithm provides precision 10 (res) times
+% smaller than pixel size. Final rendering will be reconstructed on a 10x10um grid.
+
+ULM = struct('numberOfParticles', 70,... % 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',[25 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', 10,... % Minimum allowed length of the tracks in time. (5-20 frames)
+ 'fwhm',[1 1]*7,... % Size [pixel] of the mask for localization. (3x3 for pixel at lambda, 5x5 at lambda/2). [fmwhz fmwhx]
+ 'max_gap_closing', 1,... % 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; % Safeguard on the number of maxLocal in the fwhm*fwhm grid (3 for fwhm=3, 7 for fwhm=5)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Load data and localize microbubbles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+myfun= @SRPCATours;
+blksz = [336 272 500];
+Track_tot = {};
+
+fprintf('--- ULM PROCESSING --- \n\n');t1=tic;
+for hhh = 1:Nbuffers
+ pixOffset = 2*floor((hhh-1) / 5) ;
+ fprintf('Processing bloc %d/%d\n',hhh,Nbuffers);
+
+ % Load IQ data
+ tmp = load(strcat(PALA_data_folder,'/IQTours/IQTours_',int2str(hhh),'.mat'),'IQ');
+
+ tmp.IQ = tmp.IQ(5:164,:,:);
+ tmp.IQ = imresize(tmp.IQ,SR,"bicubic");
+
+ % Temporal filter
+ tmp.IQ = filter(but_b,but_a,tmp.IQ,[],3); %(optional)
+ tmp.IQ(~isfinite(tmp.IQ))=0;
+
+ IQ = tmp.IQ/max(abs(tmp.IQ(:)));tmp = [];
+ IQ = padarray(IQ, [40 40]);
+
+ IQ = circshift(IQ, [-pixOffset -pixOffset]);
+ IQ_filt = blockproc3(IQ, blksz, myfun, [40 40 0]); IQ = [];
+ IQ_filt = circshift(IQ_filt, [pixOffset pixOffset]);
+
+ IQ_filt=IQ_filt(41:end-40,41:end-40,:);
+
+ % Detection process (return a list of super resolved coordinates in pixel)
+ [MatTracking] = ULM_localization2DSR_Tours(abs(IQ_filt),ULM,pixOffset); IQ_filt=[];
+
+ % Normalize for tracking
+ 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([pwd,'/Result/MatTracking',int2str(hhh),'.mat'],'MatTracking')
+
+ % Tracking algorithm (list of tracks)
+ Track_tot_i = ULM_tracking2D(MatTracking,ULM);
+
+ 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));
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ULM.SRscale = ULM.scale(1)/ULM.res;
+ULM.SRsize = round(ULM.size(1:2).*ULM.scale(1:2));
+%
+% % 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,'/Result/SRPCA_MatOutTours10.mat'),'MatOut','MatOut_zdir','MatOut_vel','ULM','lambda','llx','llz')
+
+%% Figures
+ULM.SRscale=0.625/ULM.res; %for velocity maps in mm/s and scale-bars
+
+IntPower = 1/2.3;
+
+cmin = min(MatOut(:).^IntPower);
+cmax = max(MatOut(:).^IntPower);
+numTicks = 7; % Adjust as needed
+
+h1=figure('WindowState', 'maximized');
+im=imagesc((MatOut).^IntPower);
+axis image
+colormap hot
+clbar = colorbar;
+set(clbar, 'YTick', linspace(cmin, cmax, numTicks));
+
+clbar.Label.String = 'Normalized localization density';
+%clbar.TickLabels = clbar.Ticks.^2.3; %Not normalized
+clbar.TickLabels = round(clbar.Ticks.^2.3 ./ max(clbar.Ticks(:).^2.3),2);
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:2],60+[0:BarWidth])=max(clim);
+text(60+BarWidth/2, size(MatOut,1)-50+3 +15, '1 mm', ...
+ 'Color', 'w', 'FontSize', 12, 'FontWeight', 'bold', ...
+ 'HorizontalAlignment', 'center');
+
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2)+clbar.Ticks(end)/2, x(3)];
+
+%%
+IntPower = 1/2.3;
+
+OriginalData = MatOut.^IntPower.*sign(imgaussfilt(MatOut_zdir,.8));
+cmin = min(OriginalData(:));
+cmax = max(OriginalData(:));
+clear OriginalData;
+numTicks = 9; % Adjust as needed
+
+h1=figure('WindowState', 'maximized');
+velColormap = cat(1,flip(flip(hot(128),1),2),hot(128));
+velColormap = velColormap(5:end-5,:);
+im=imagesc(llx,llz,(MatOut).^IntPower.*sign(imgaussfilt(MatOut_zdir,.8)))
+im.CData = im.CData - sign(im.CData)/2;axis image
+colormap(gca,velColormap)
+
+
+clbar = colorbar;
+set(clbar, 'YTick', linspace(cmin*0.90, cmax*0.90, numTicks));
+
+clbar.Label.String = 'Normalized localization density';
+clbar.Ticks = clbar.Ticks;%customTicks;
+TempTickLabels = round(sign(clbar.Ticks).*abs(clbar.Ticks/0.90).^2.3);
+clbar.TickLabels = round(TempTickLabels ./ max(TempTickLabels(:)),2);
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+caxis([-1 1]*max(caxis)*.75)
+
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:3],60+[0:BarWidth])=max(caxis);
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2), x(3)];
+
+%%
+
+
+
+%vmax_disp = ceil(quantile(MatOut_vel(abs(MatOut_vel)>0),.98)/10)*10;
+vmax_disp = 80;
+
+h1=figure('WindowState', 'maximized'); clf,set(gcf,'Position',[652 393 941 585]);
+clbsize = [180,50];
+Mvel_rgb = MatOut_vel/vmax_disp; % normalization
+% Mvel_rgb(1:clbsize(1),1:clbsize(2)) = repmat(linspace(1,0,clbsize(1))',1,clbsize(2)); % add velocity colorbar
+Mvel_rgb = Mvel_rgb.^(1/2);Mvel_rgb(Mvel_rgb>1)=1;
+Mvel_rgb = imgaussfilt(Mvel_rgb,.5);
+Mvel_rgb = ind2rgb(round(Mvel_rgb*256),jet(256)); % convert ind into RGB
+IntPower = 1/2.3;
+MatShadow = MatOut;MatShadow = MatShadow./max(MatShadow(:)*.3);MatShadow(MatShadow>1)=1;
+%MatShadow(1:clbsize(1),1:clbsize(2))=repmat(linspace(0,1,clbsize(2)),clbsize(1),1);
+Mvel_rgb = Mvel_rgb.*(MatShadow.^IntPower);
+Mvel_rgb = brighten(Mvel_rgb,.4);
+im=imshow(Mvel_rgb,'XData',llx,'YData',llz);axis on
+BarWidth = round(1./(ULM.SRscale*lambda)); % 1 mm
+im.CData(size(MatOut,1)-50+[0:3],60+[0:BarWidth],1:3)=1;
+
+
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+
+colormap(jet);
+clbar = colorbar;
+clim([0 vmax_disp]);
+clbar.Label.String = 'Velocity magnitude (mm/s)';
+
+% Smaller Colorbar
+% If there is a display issue, try executing the code below separately
+% after you created the figure h1 as it depends on your screen size
+% x1=get(gca,'position');
+% x=get(clbar,'Position');
+% x(1)=x(1)-0.003;
+% x(3)=0.005;
+% set(clbar,'Position',x)
+% set(gca,'position',x1)
+% clbar.Label.Position = [x(1)+2.5, x(2)+clbar.Ticks(end)/2, x(3)];
+
+%% Metrics
+% Saturation Value
+% load(strcat(pwd,'/Result/SRPCA_MatOutTours10.mat'),'MatOut')
+Saturation = nnz(MatOut()>0)/numel(MatOut());
+
+
+% CR
+% load(strcat(pwd,'/Result/SRPCA_MatOutTours10.mat'),'MatOut')
+CR = ContrastRatio(MatOut);
+
+% Image for FRC estimation
+% load(strcat(pwd,'/Result/SRPCA_MatOutTours10.mat'),'MatOut')
+
+figure(1);
+IntPower = 1/2.3;
+im=imagesc(llx,llz,MatOut.^IntPower);axis image
+% im=imagesc(llx,llz,PowDop2);axis image
+set(gca,'XColor', 'none','YColor','none')
+set(gca, 'color', 'none');
+im.CData = max(im.CData-.5,0);
+
+colormap(gca,gray(128))
+caxis(caxis*.7) % add saturation in image
+ca = gca;
+
+WriteTif(im.CData,ca.Colormap,strcat(pwd,'/Result/SRPCA_MatOutTours10.tif'),'caxis',caxis,'Overwrite',1)
+
+% Gridding index
+load(strcat(pwd,'/Result/SRPCA_MatOutTours10_Gridding.mat'),'MatOut')
+Gridding_index = PALA_ComputeGriddingIndex({MatOut},1);
+Gridding = 1-Gridding_index/30; % [100 at 0; 0 at -30]
diff --git a/InVivoTours/Result/MatTracks/MatTracking1.mat b/InVivoTours/Result/MatTracks/MatTracking1.mat
new file mode 100644
index 0000000000000000000000000000000000000000..90025d20ca3bad8882efcc7f785f5103e9835b93
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking1.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking10.mat b/InVivoTours/Result/MatTracks/MatTracking10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..139e436ace88fb6b810125d0c8f372bf267af410
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking10.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking11.mat b/InVivoTours/Result/MatTracks/MatTracking11.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5f5af2dde3277e192ed71c4b11418f8a64d674ff
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking11.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking12.mat b/InVivoTours/Result/MatTracks/MatTracking12.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c0e79b2c71ea9e8ed117b7bcdd83b1862b74d6b0
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking12.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking13.mat b/InVivoTours/Result/MatTracks/MatTracking13.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b2b878507f893469b84ff004e8f5b9a0420d5c9c
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking13.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking14.mat b/InVivoTours/Result/MatTracks/MatTracking14.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3504f078fe2337f56a6b5ac58495d8d5500d6754
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking14.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking15.mat b/InVivoTours/Result/MatTracks/MatTracking15.mat
new file mode 100644
index 0000000000000000000000000000000000000000..90230ba1482776133c7180ee373ff8f057ba69d9
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking15.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking16.mat b/InVivoTours/Result/MatTracks/MatTracking16.mat
new file mode 100644
index 0000000000000000000000000000000000000000..11ad938b6fb822d44e5258a580a0fc1da4dde9b5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking16.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking17.mat b/InVivoTours/Result/MatTracks/MatTracking17.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2a58a1e2f67a27b12e281f150db53fcecc0b096b
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking17.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking18.mat b/InVivoTours/Result/MatTracks/MatTracking18.mat
new file mode 100644
index 0000000000000000000000000000000000000000..62cb6e27fcb270eafa19338c649a338e3b1b8d06
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking18.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking19.mat b/InVivoTours/Result/MatTracks/MatTracking19.mat
new file mode 100644
index 0000000000000000000000000000000000000000..43f3a8a9a9a87f99f97b26c6d9c7421aa6dc8a42
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking19.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking2.mat b/InVivoTours/Result/MatTracks/MatTracking2.mat
new file mode 100644
index 0000000000000000000000000000000000000000..989d7caae1e57e7d2608b2956f1efc2c3ef9e795
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking2.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking20.mat b/InVivoTours/Result/MatTracks/MatTracking20.mat
new file mode 100644
index 0000000000000000000000000000000000000000..036c0c445776257ea02700a5261f08db0b7becb5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking20.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking21.mat b/InVivoTours/Result/MatTracks/MatTracking21.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2b5c28c00cf04aab1572108b0ba5da00213c0d66
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking21.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking22.mat b/InVivoTours/Result/MatTracks/MatTracking22.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0bcd54ba940d36e40076061da82eb3f0152ccb18
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking22.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking23.mat b/InVivoTours/Result/MatTracks/MatTracking23.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1ced8ec14d3c95f8ada4f1d2b1151db60162bd7b
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking23.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking24.mat b/InVivoTours/Result/MatTracks/MatTracking24.mat
new file mode 100644
index 0000000000000000000000000000000000000000..30ad268384aac4f49a0e516e3d6694663d79e809
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking24.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking25.mat b/InVivoTours/Result/MatTracks/MatTracking25.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8633c36de0c8bbe752d433a2e6c3f2a05005c231
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking25.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking26.mat b/InVivoTours/Result/MatTracks/MatTracking26.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3cf7cf657ed9406382f56e577783c19a603d3e4c
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking26.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking27.mat b/InVivoTours/Result/MatTracks/MatTracking27.mat
new file mode 100644
index 0000000000000000000000000000000000000000..552ae5c527990ea5d695f79bf1eb825a9f5eb097
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking27.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking28.mat b/InVivoTours/Result/MatTracks/MatTracking28.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4a8f103c5517902cab5126e7311ed0955ca6e1be
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking28.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking29.mat b/InVivoTours/Result/MatTracks/MatTracking29.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6f4b19121019041ae915be7fe77c798cad16b4d1
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking29.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking3.mat b/InVivoTours/Result/MatTracks/MatTracking3.mat
new file mode 100644
index 0000000000000000000000000000000000000000..110060438dd331e60f2fce0e3a1388a5132a03c6
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking3.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking30.mat b/InVivoTours/Result/MatTracks/MatTracking30.mat
new file mode 100644
index 0000000000000000000000000000000000000000..af6e2f41222c33f58c7de27d249798994bcc94e3
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking30.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking31.mat b/InVivoTours/Result/MatTracks/MatTracking31.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c40d8b862b9c1159b128333e42989e6b1fd03e90
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking31.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking32.mat b/InVivoTours/Result/MatTracks/MatTracking32.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a8f620ed925db82be8eed043343f353dbbacc108
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking32.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking33.mat b/InVivoTours/Result/MatTracks/MatTracking33.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9f718764f780cd3cfdd63695221b33d297c28d96
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking33.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking34.mat b/InVivoTours/Result/MatTracks/MatTracking34.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9605edaa93e89ff6bb7518b8cbb582b3b79a204d
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking34.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking35.mat b/InVivoTours/Result/MatTracks/MatTracking35.mat
new file mode 100644
index 0000000000000000000000000000000000000000..ba5b02bc17b2f1c70b7eb6a73c2745fb86dc9b40
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking35.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking36.mat b/InVivoTours/Result/MatTracks/MatTracking36.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6bd86983be818779b13bdef78b85761df04a2a59
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking36.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking37.mat b/InVivoTours/Result/MatTracks/MatTracking37.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b6176953b06d1be9fe824ce248fd9a5361cb5185
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking37.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking38.mat b/InVivoTours/Result/MatTracks/MatTracking38.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9859010eebd255e2b84ee545461a37c32738016a
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking38.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking39.mat b/InVivoTours/Result/MatTracks/MatTracking39.mat
new file mode 100644
index 0000000000000000000000000000000000000000..cf72906a8d294252cfb2fcafbe89b0ffbc3c665c
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking39.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking4.mat b/InVivoTours/Result/MatTracks/MatTracking4.mat
new file mode 100644
index 0000000000000000000000000000000000000000..72edecc619c6ef4f0f63322c00e98ac42f116900
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking4.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking40.mat b/InVivoTours/Result/MatTracks/MatTracking40.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b911900347169772eb39f6edbea13fc956b1a239
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking40.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking41.mat b/InVivoTours/Result/MatTracks/MatTracking41.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2657359bde72f8dc490ebab0eedaa821c9444936
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking41.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking42.mat b/InVivoTours/Result/MatTracks/MatTracking42.mat
new file mode 100644
index 0000000000000000000000000000000000000000..cca4985e41ce802e052fdb2904c3e037ee7c62c5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking42.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking43.mat b/InVivoTours/Result/MatTracks/MatTracking43.mat
new file mode 100644
index 0000000000000000000000000000000000000000..65672988fb5ad84a5217bedb2bd2e1f65cc96587
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking43.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking44.mat b/InVivoTours/Result/MatTracks/MatTracking44.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d50a7446022ce6f70c2f2a3f0f7e15330309cef0
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking44.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking45.mat b/InVivoTours/Result/MatTracks/MatTracking45.mat
new file mode 100644
index 0000000000000000000000000000000000000000..abe840ea2a36c836bd93e892952950d2200071ae
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking45.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking46.mat b/InVivoTours/Result/MatTracks/MatTracking46.mat
new file mode 100644
index 0000000000000000000000000000000000000000..5fb93ce2b0ff98d000ce057775cc04cf9f3ae1d6
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking46.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking47.mat b/InVivoTours/Result/MatTracks/MatTracking47.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b0ca1402937c5d62bc1ae185be0cfa3ff2f77751
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking47.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking48.mat b/InVivoTours/Result/MatTracks/MatTracking48.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3cb20ecbc25496ca6d5509dc68bf63ef0d80c0ab
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking48.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking49.mat b/InVivoTours/Result/MatTracks/MatTracking49.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8804abc5bf41f80e7f7d05f72728b746778d032e
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking49.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking5.mat b/InVivoTours/Result/MatTracks/MatTracking5.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9db382afadc389766bba41a79c7aef0c51c1a31c
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking5.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking50.mat b/InVivoTours/Result/MatTracks/MatTracking50.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d04fe4c6a4df5d15dcf1c0eedd4378e948f0f6e8
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking50.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking51.mat b/InVivoTours/Result/MatTracks/MatTracking51.mat
new file mode 100644
index 0000000000000000000000000000000000000000..61ca5627155ef6c45f6976cd88ae71423318c8c5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking51.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking52.mat b/InVivoTours/Result/MatTracks/MatTracking52.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a859c03bc228277fe254f9e52e8396e782d7aa3f
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking52.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking53.mat b/InVivoTours/Result/MatTracks/MatTracking53.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8618caddc7f4cda624170469b31149bc5bb143c0
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking53.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking54.mat b/InVivoTours/Result/MatTracks/MatTracking54.mat
new file mode 100644
index 0000000000000000000000000000000000000000..47116ae1f5da4419eba773ea2ae363bde696c906
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking54.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking55.mat b/InVivoTours/Result/MatTracks/MatTracking55.mat
new file mode 100644
index 0000000000000000000000000000000000000000..83b128eb3ae523638d160067774ac4b3e6e83ec5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking55.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking56.mat b/InVivoTours/Result/MatTracks/MatTracking56.mat
new file mode 100644
index 0000000000000000000000000000000000000000..4cde3e94bbbd0c23d3594ff6e18f965696883530
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking56.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking57.mat b/InVivoTours/Result/MatTracks/MatTracking57.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7958c2a304dbfeece0575ac80275dd92ff7eec54
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking57.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking58.mat b/InVivoTours/Result/MatTracks/MatTracking58.mat
new file mode 100644
index 0000000000000000000000000000000000000000..90dcef7a9e6bd905becab04d20823891c774de8d
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking58.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking59.mat b/InVivoTours/Result/MatTracks/MatTracking59.mat
new file mode 100644
index 0000000000000000000000000000000000000000..372026120a2358d6a2b22996dd270eb35066c79b
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking59.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking6.mat b/InVivoTours/Result/MatTracks/MatTracking6.mat
new file mode 100644
index 0000000000000000000000000000000000000000..717d79bf17680518e9d08160c4ae8b2472229982
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking6.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking60.mat b/InVivoTours/Result/MatTracks/MatTracking60.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e2b986e9eb1755611a4f40d78457b767eaba9dc9
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking60.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking61.mat b/InVivoTours/Result/MatTracks/MatTracking61.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9a18ec4892041ff43692b2977ae189a1c42758aa
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking61.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking62.mat b/InVivoTours/Result/MatTracks/MatTracking62.mat
new file mode 100644
index 0000000000000000000000000000000000000000..6e077940720fc714c5ac9ff41bfa3447ff342933
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking62.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking63.mat b/InVivoTours/Result/MatTracks/MatTracking63.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d708fb88d8a7da6865101ca6094ef5db766f03f5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking63.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking64.mat b/InVivoTours/Result/MatTracks/MatTracking64.mat
new file mode 100644
index 0000000000000000000000000000000000000000..86c81fcae03a418cac287c08786d31897639ade2
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking64.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking65.mat b/InVivoTours/Result/MatTracks/MatTracking65.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e42583dcf20545b14a6fae1b9b41267d1c5e241f
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking65.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking66.mat b/InVivoTours/Result/MatTracks/MatTracking66.mat
new file mode 100644
index 0000000000000000000000000000000000000000..78cf768685a6f04f23bf06dc7b35a3ce8fb37a80
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking66.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking67.mat b/InVivoTours/Result/MatTracks/MatTracking67.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a91035d2c669954d648c971ef1e53af159e9fcba
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking67.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking68.mat b/InVivoTours/Result/MatTracks/MatTracking68.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3b8014100a530c236923b5aff08966045fcd283d
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking68.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking69.mat b/InVivoTours/Result/MatTracks/MatTracking69.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f24af09378e7b9f97a3e58e6c97d799b90abe5ac
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking69.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking7.mat b/InVivoTours/Result/MatTracks/MatTracking7.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8f0afa9461fbd04b403ee53b90fe7280103c9f66
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking7.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking70.mat b/InVivoTours/Result/MatTracks/MatTracking70.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7eedcb881c6934f6ab8d80d4727b7c8c111d6a77
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking70.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking71.mat b/InVivoTours/Result/MatTracks/MatTracking71.mat
new file mode 100644
index 0000000000000000000000000000000000000000..2eb8ae43eec12b7a67ced2e125b62b0d4167d9e9
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking71.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking72.mat b/InVivoTours/Result/MatTracks/MatTracking72.mat
new file mode 100644
index 0000000000000000000000000000000000000000..850c4462f2319b65161f08ce2b587aebb23c40fa
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking72.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking73.mat b/InVivoTours/Result/MatTracks/MatTracking73.mat
new file mode 100644
index 0000000000000000000000000000000000000000..9cc04ae645c49067ed5b0af238ab38cb3d24e404
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking73.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking74.mat b/InVivoTours/Result/MatTracks/MatTracking74.mat
new file mode 100644
index 0000000000000000000000000000000000000000..ef6c0d8b7c59aa3c8e837a21cfc87f5f13cadf15
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking74.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking75.mat b/InVivoTours/Result/MatTracks/MatTracking75.mat
new file mode 100644
index 0000000000000000000000000000000000000000..f6c09bb4e3b210d0c2519e2d7d071dd778e63910
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking75.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking76.mat b/InVivoTours/Result/MatTracks/MatTracking76.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1fa72500c62680f8e9830e164002216bd3e1855f
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking76.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking77.mat b/InVivoTours/Result/MatTracks/MatTracking77.mat
new file mode 100644
index 0000000000000000000000000000000000000000..842bcd2d32929136d9877189599f1892bd46cde7
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking77.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking78.mat b/InVivoTours/Result/MatTracks/MatTracking78.mat
new file mode 100644
index 0000000000000000000000000000000000000000..114349a62996fc38f51d0defeef8cff384c9e857
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking78.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking79.mat b/InVivoTours/Result/MatTracks/MatTracking79.mat
new file mode 100644
index 0000000000000000000000000000000000000000..3811531a4484701d8d2f7a261334166033041517
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking79.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking8.mat b/InVivoTours/Result/MatTracks/MatTracking8.mat
new file mode 100644
index 0000000000000000000000000000000000000000..1b44a4c789876972ee7edfc66b0ea06fe29bc73a
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking8.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking80.mat b/InVivoTours/Result/MatTracks/MatTracking80.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e0891a0fdb9f82a22b768689d40a73bd4d7c53c5
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking80.mat differ
diff --git a/InVivoTours/Result/MatTracks/MatTracking9.mat b/InVivoTours/Result/MatTracks/MatTracking9.mat
new file mode 100644
index 0000000000000000000000000000000000000000..fb03b9db41e7352bb85ca48ae400a8347cfd0af0
Binary files /dev/null and b/InVivoTours/Result/MatTracks/MatTracking9.mat differ
diff --git a/InVivoTours/Result/RSULM_MatOutTours10.mat b/InVivoTours/Result/RSULM_MatOutTours10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..aa0a15b3bc5a22c8ef9c3cbad690171a73d3933a
Binary files /dev/null and b/InVivoTours/Result/RSULM_MatOutTours10.mat differ
diff --git a/InVivoTours/Result/RSULM_MatOutTours10.tif b/InVivoTours/Result/RSULM_MatOutTours10.tif
new file mode 100644
index 0000000000000000000000000000000000000000..96219df64a1b5ea489cdd9fda782bdca8c377142
Binary files /dev/null and b/InVivoTours/Result/RSULM_MatOutTours10.tif differ
diff --git a/InVivoTours/Result/RSULM_MatOutTours10_Gridding.mat b/InVivoTours/Result/RSULM_MatOutTours10_Gridding.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7db855b2c6df3e84a9e70a159528434ed5af411c
Binary files /dev/null and b/InVivoTours/Result/RSULM_MatOutTours10_Gridding.mat differ
diff --git a/InVivoTours/Result/SRPCANonInterpTracks.mat b/InVivoTours/Result/SRPCANonInterpTracks.mat
new file mode 100644
index 0000000000000000000000000000000000000000..01f688f540da29b6a612335478672daf18272f85
Binary files /dev/null and b/InVivoTours/Result/SRPCANonInterpTracks.mat differ
diff --git a/InVivoTours/Result/SRPCA_MatOutTours10.mat b/InVivoTours/Result/SRPCA_MatOutTours10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..95a8da7cb733ac53aa7b06cb51b0ba5d748ab9f9
Binary files /dev/null and b/InVivoTours/Result/SRPCA_MatOutTours10.mat differ
diff --git a/InVivoTours/Result/SRPCA_MatOutTours10.tif b/InVivoTours/Result/SRPCA_MatOutTours10.tif
new file mode 100644
index 0000000000000000000000000000000000000000..6f37aae519d892ccb5541b4e05e2fd76027edeae
Binary files /dev/null and b/InVivoTours/Result/SRPCA_MatOutTours10.tif differ
diff --git a/InVivoTours/Result/SRPCA_MatOutTours10_Gridding.mat b/InVivoTours/Result/SRPCA_MatOutTours10_Gridding.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c1d064a1b0a1989322d330fca6830408ac1e264b
Binary files /dev/null and b/InVivoTours/Result/SRPCA_MatOutTours10_Gridding.mat differ
diff --git a/InVivoTours/Result/SRPCA_MatOutTours4.mat b/InVivoTours/Result/SRPCA_MatOutTours4.mat
new file mode 100644
index 0000000000000000000000000000000000000000..964a2e067ed7440399e08563052e36110ffdad8e
Binary files /dev/null and b/InVivoTours/Result/SRPCA_MatOutTours4.mat differ
diff --git a/InVivoTours/SRPCATours.m b/InVivoTours/SRPCATours.m
new file mode 100644
index 0000000000000000000000000000000000000000..238c16c309a82a21eb9c000cd401a7f4ba9a9fad
--- /dev/null
+++ b/InVivoTours/SRPCATours.m
@@ -0,0 +1,104 @@
+function x = SRPCATours(S)
+
+% 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;
+
+ load('psf10Tours.mat', 'H')
+ [Mh, Nh] = size(H);
+ center = round([Mh, Nh] / 2);
+ H = fft2(circshift(padarray(H, [m - Mh, n - Nh], 'post'), 1 - center));
+
+
+ lambda = 0.1333; %lambda = 20*SR/sqrt(m*n);
+
+ mu = 20;
+
+ rho = 1e-3;
+
+ eps = 1.5;
+
+ rho_max = 1;
+
+ max_iter = 20;
+
+ %% 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/InVivoTours/ULM_localization2DSR_Tours.m b/InVivoTours/ULM_localization2DSR_Tours.m
new file mode 100644
index 0000000000000000000000000000000000000000..7d10271346633f39ae4cc6d16fac80aca7efab42
--- /dev/null
+++ b/InVivoTours/ULM_localization2DSR_Tours.m
@@ -0,0 +1,182 @@
+function MatTracking = ULM_localization2DSR_Tours(MatIn,ULM,pixOffset)
+%% 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 ismember(index_mask_x(iscat), [232 233 504 505 776 777 1048 1049]+pixOffset) || ismember(index_mask_z(iscat), [296 297 632 633 968 969 1304 1305]+pixOffset)
+ continue
+ end
+
+ 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
+
diff --git a/InVivoTours/psf10Tours.mat b/InVivoTours/psf10Tours.mat
new file mode 100644
index 0000000000000000000000000000000000000000..c08d7cb5c5beaea4ec951ca9230b3a2d0dbf911b
Binary files /dev/null and b/InVivoTours/psf10Tours.mat differ
diff --git a/InVivoTours/psfSR4_Tours.mat b/InVivoTours/psfSR4_Tours.mat
new file mode 100644
index 0000000000000000000000000000000000000000..62259618ba33dabfc343177113301249ef1a5594
Binary files /dev/null and b/InVivoTours/psfSR4_Tours.mat differ