diff --git a/InVivoPalaSup/Fig7db_InVivoPala_SRPCA.m b/InVivoPalaSup/Fig7db_InVivoPala_SRPCA.m
new file mode 100644
index 0000000000000000000000000000000000000000..f872638d2cae1c7dcd3a8af45be6e57058acb119
--- /dev/null
+++ b/InVivoPalaSup/Fig7db_InVivoPala_SRPCA.m
@@ -0,0 +1,311 @@
+%% InVivoPala_SRPCA.m 
+% 
+%
+% Based on the code created by Arthur Chavignon
+% Adapted for the Pala supplementary in vivo data: https://zenodo.org/records/7883227
+%
+% 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
+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 = 15.625;
+UF.FrameRateUF = 1000;
+
+% Here you put the size of your data
+SizeOfBloc = [SR*180 SR*256 800];
+
+
+% 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 = 250;          % 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 = 1540/(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 = 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', 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',[25 SizeOfBloc(3)],...  % SVD filtering, to be adapted to your clutter/SNR levels
+    'max_linking_distance',3,...        % 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]*3,...                  % 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= @SRPCAPala;
+Track_tot = {};
+blksz = [360 512 400];
+
+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 (or other kind of images without compression)
+    % tmp = load([IQfiles(hhh).folder filesep IQfiles(hhh).name],'IQ');
+
+    tmp = load(strcat(PALA_data_folder,'/IQRBS/IQRBS_',int2str(hhh),'.mat'),'IQ');
+
+    % Temporal filtering
+    tmp.IQ = filter(but_b,but_a,tmp.IQ,[],3); %(optional)
+    tmp.IQ(~isfinite(tmp.IQ))=0;
+
+    IQ = imresize(tmp.IQ,SR,"bicubic");tmp = [];
+    IQ = IQ/max(abs(IQ(:)));
+
+    IQ =  circshift(IQ, [-pixOffset -pixOffset]);
+    IQ_filt = blockproc3(IQ, blksz, myfun, [40 40 20]); IQ = [];
+    IQ_filt =  circshift(IQ_filt, [pixOffset pixOffset]);
+
+    % Detection process (return a list of super resolved coordinates in pixel)
+    [MatTracking] = ULM_localization2DSR_Pala(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_MatOutPala10.mat'),'MatOut','MatOut_zdir','MatOut_vel','ULM','lambda','llx','llz')
+
+%% Crop empty columns and lines
+
+MatOut=MatOut(41:end,81+60:end-80-60);
+MatOut_vel=MatOut_vel(41:end,81+60:end-80-60);
+MatOut_zdir=MatOut_zdir(41:end,81+60:end-80-60);
+
+llz = llz(1:1761);
+llx = llx(1:2281);
+
+%% Figures
+ULM.SRscale=0.50731/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
+
+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;
+
+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_MatOutPala10.mat'),'MatOut')
+% MatOut=MatOut(41:end,81+60:end-80-60);
+Saturation  = nnz(MatOut()>0)/numel(MatOut());  %Calculated on Cropped Matrix
+
+
+% CR
+% load(strcat(pwd,'/Result/SRPCA_MatOutPala10.mat'),'MatOut')
+% MatOut=MatOut(41:end,81+60:end-80-60);
+CR = ContrastRatio(MatOut);
+
+% Image for FRC estimation
+% load(strcat(pwd,'/Result/SRPCA_MatOutPala10.mat'),'MatOut')
+% MatOut=MatOut(41:end,81+60:end-80-60);
+
+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/SRPCA_MatOutPala10.tif'),'caxis',caxis,'Overwrite',1)
+
+% Gridding index
+load(strcat(pwd,'/Result/SRPCA_MatOutPala10_Gridding.mat'),'MatOut')
+MatOut=MatOut(41:end,81+60:end-80-60);
+Gridding_index = PALA_ComputeGriddingIndex({MatOut},1);  %Calculated on Cropped Matrix
+Gridding = 1-Gridding_index/30; % [100 at 0; 0 at -30]
\ No newline at end of file
diff --git a/InVivoPalaSup/H10.mat b/InVivoPalaSup/H10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..ec64f162eb7dfa568b569d5134a276c6b84acc46
Binary files /dev/null and b/InVivoPalaSup/H10.mat differ
diff --git a/InVivoPalaSup/Result/RSULM_MatOutPala10.mat b/InVivoPalaSup/Result/RSULM_MatOutPala10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..99fd1e77712648144fa83bb3397814b96d4346ba
Binary files /dev/null and b/InVivoPalaSup/Result/RSULM_MatOutPala10.mat differ
diff --git a/InVivoPalaSup/Result/RSULM_MatOutPala10.tif b/InVivoPalaSup/Result/RSULM_MatOutPala10.tif
new file mode 100644
index 0000000000000000000000000000000000000000..0bb2f3e3a49122e0cec55dee6ee99272b085f65b
Binary files /dev/null and b/InVivoPalaSup/Result/RSULM_MatOutPala10.tif differ
diff --git a/InVivoPalaSup/Result/RSULM_MatOutPala10_Gridding.mat b/InVivoPalaSup/Result/RSULM_MatOutPala10_Gridding.mat
new file mode 100644
index 0000000000000000000000000000000000000000..62b5d70f6342db6f32369ac5ca849143aa565278
Binary files /dev/null and b/InVivoPalaSup/Result/RSULM_MatOutPala10_Gridding.mat differ
diff --git a/InVivoPalaSup/Result/SRPCANonInterpTracks.mat b/InVivoPalaSup/Result/SRPCANonInterpTracks.mat
new file mode 100644
index 0000000000000000000000000000000000000000..581b07e7d3a332952083575a01c4a1150840861c
Binary files /dev/null and b/InVivoPalaSup/Result/SRPCANonInterpTracks.mat differ
diff --git a/InVivoPalaSup/Result/SRPCA_MatOutPala10.mat b/InVivoPalaSup/Result/SRPCA_MatOutPala10.mat
new file mode 100644
index 0000000000000000000000000000000000000000..117f8aaafa47fea89e33bda9c3dca3661ca5f934
Binary files /dev/null and b/InVivoPalaSup/Result/SRPCA_MatOutPala10.mat differ
diff --git a/InVivoPalaSup/Result/SRPCA_MatOutPala10.tif b/InVivoPalaSup/Result/SRPCA_MatOutPala10.tif
new file mode 100644
index 0000000000000000000000000000000000000000..b32585d6ced0cb9bbabbfefefae1e37fc9953bb1
Binary files /dev/null and b/InVivoPalaSup/Result/SRPCA_MatOutPala10.tif differ
diff --git a/InVivoPalaSup/Result/SRPCA_MatOutPala10_Gridding.mat b/InVivoPalaSup/Result/SRPCA_MatOutPala10_Gridding.mat
new file mode 100644
index 0000000000000000000000000000000000000000..e9acbda1e4ec4d736f8f0d959a6516c24667f104
Binary files /dev/null and b/InVivoPalaSup/Result/SRPCA_MatOutPala10_Gridding.mat differ
diff --git a/InVivoPalaSup/Result/SRPCA_MatOutPala4.mat b/InVivoPalaSup/Result/SRPCA_MatOutPala4.mat
new file mode 100644
index 0000000000000000000000000000000000000000..01257e56b080c132d6bce1f30c155f9fc645475b
Binary files /dev/null and b/InVivoPalaSup/Result/SRPCA_MatOutPala4.mat differ
diff --git a/InVivoPalaSup/SRPCAPala.m b/InVivoPalaSup/SRPCAPala.m
new file mode 100644
index 0000000000000000000000000000000000000000..02af8c7a91d92fdc355f20ec8fbecb544a97e9b4
--- /dev/null
+++ b/InVivoPalaSup/SRPCAPala.m
@@ -0,0 +1,103 @@
+function x = SRPCAPala(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('H10.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.0116; %lambda = 2.5*SR/sqrt(m*n);
+
+    mu = 10; 
+    
+    rho = 1e-6;
+
+    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/InVivoPalaSup/ULM_localization2DSR_Pala.m b/InVivoPalaSup/ULM_localization2DSR_Pala.m
new file mode 100644
index 0000000000000000000000000000000000000000..0db7f7afef32b3601ed121c5f90d7b0ffa7a4e68
--- /dev/null
+++ b/InVivoPalaSup/ULM_localization2DSR_Pala.m
@@ -0,0 +1,182 @@
+function MatTracking = ULM_localization2DSR_Pala(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), [512 513 1024 1025 1536 1537 2048 2049]+pixOffset) ||  ismember(index_mask_z(iscat), [360 361 720 721 1080 1081 1440 1441]+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/InVivoPalaSup/psfSR4.mat b/InVivoPalaSup/psfSR4.mat
new file mode 100644
index 0000000000000000000000000000000000000000..443743a7bf33a442430c4cfcb6e828f62acd7508
Binary files /dev/null and b/InVivoPalaSup/psfSR4.mat differ