Skip to content
Snippets Groups Projects
Commit 5bff4fd4 authored by vpustova's avatar vpustova
Browse files

Update 94 files

- /InVivoTours/Result/RSULM_MatOutTours10_Gridding.mat
- /InVivoTours/Result/RSULM_MatOutTours10.tif
- /InVivoTours/Result/SRPCA_MatOutTours10.tif
- /InVivoTours/Result/SRPCA_MatOutTours4.mat
- /InVivoTours/Result/SRPCA_MatOutTours10_Gridding.mat
- /InVivoTours/Result/RSULM_MatOutTours10.mat
- /InVivoTours/Result/SRPCANonInterpTracks.mat
- /InVivoTours/Result/SRPCA_MatOutTours10.mat
- /InVivoTours/Result/MatTracks/MatTracking1.mat
- /InVivoTours/Result/MatTracks/MatTracking10.mat
- /InVivoTours/Result/MatTracks/MatTracking11.mat
- /InVivoTours/Result/MatTracks/MatTracking12.mat
- /InVivoTours/Result/MatTracks/MatTracking13.mat
- /InVivoTours/Result/MatTracks/MatTracking14.mat
- /InVivoTours/Result/MatTracks/MatTracking15.mat
- /InVivoTours/Result/MatTracks/MatTracking16.mat
- /InVivoTours/Result/MatTracks/MatTracking17.mat
- /InVivoTours/Result/MatTracks/MatTracking18.mat
- /InVivoTours/Result/MatTracks/MatTracking19.mat
- /InVivoTours/Result/MatTracks/MatTracking2.mat
- /InVivoTours/Result/MatTracks/MatTracking20.mat
- /InVivoTours/Result/MatTracks/MatTracking21.mat
- /InVivoTours/Result/MatTracks/MatTracking22.mat
- /InVivoTours/Result/MatTracks/MatTracking23.mat
- /InVivoTours/Result/MatTracks/MatTracking24.mat
- /InVivoTours/Result/MatTracks/MatTracking25.mat
- /InVivoTours/Result/MatTracks/MatTracking26.mat
- /InVivoTours/Result/MatTracks/MatTracking27.mat
- /InVivoTours/Result/MatTracks/MatTracking28.mat
- /InVivoTours/Result/MatTracks/MatTracking3.mat
- /InVivoTours/Result/MatTracks/MatTracking29.mat
- /InVivoTours/Result/MatTracks/MatTracking30.mat
- /InVivoTours/Result/MatTracks/MatTracking31.mat
- /InVivoTours/Result/MatTracks/MatTracking32.mat
- /InVivoTours/Result/MatTracks/MatTracking33.mat
- /InVivoTours/Result/MatTracks/MatTracking34.mat
- /InVivoTours/Result/MatTracks/MatTracking35.mat
- /InVivoTours/Result/MatTracks/MatTracking36.mat
- /InVivoTours/Result/MatTracks/MatTracking37.mat
- /InVivoTours/Result/MatTracks/MatTracking38.mat
- /InVivoTours/Result/MatTracks/MatTracking39.mat
- /InVivoTours/Result/MatTracks/MatTracking4.mat
- /InVivoTours/Result/MatTracks/MatTracking40.mat
- /InVivoTours/Result/MatTracks/MatTracking42.mat
- /InVivoTours/Result/MatTracks/MatTracking45.mat
- /InVivoTours/Result/MatTracks/MatTracking41.mat
- /InVivoTours/Result/MatTracks/MatTracking44.mat
- /InVivoTours/Result/MatTracks/MatTracking43.mat
- /InVivoTours/Result/MatTracks/MatTracking46.mat
- /InVivoTours/Result/MatTracks/MatTracking47.mat
- /InVivoTours/Result/MatTracks/MatTracking48.mat
- /InVivoTours/Result/MatTracks/MatTracking49.mat
- /InVivoTours/Result/MatTracks/MatTracking5.mat
- /InVivoTours/Result/MatTracks/MatTracking50.mat
- /InVivoTours/Result/MatTracks/MatTracking51.mat
- /InVivoTours/Result/MatTracks/MatTracking52.mat
- /InVivoTours/Result/MatTracks/MatTracking53.mat
- /InVivoTours/Result/MatTracks/MatTracking54.mat
- /InVivoTours/Result/MatTracks/MatTracking55.mat
- /InVivoTours/Result/MatTracks/MatTracking56.mat
- /InVivoTours/Result/MatTracks/MatTracking57.mat
- /InVivoTours/Result/MatTracks/MatTracking58.mat
- /InVivoTours/Result/MatTracks/MatTracking59.mat
- /InVivoTours/Result/MatTracks/MatTracking6.mat
- /InVivoTours/Result/MatTracks/MatTracking60.mat
- /InVivoTours/Result/MatTracks/MatTracking62.mat
- /InVivoTours/Result/MatTracks/MatTracking61.mat
- /InVivoTours/Result/MatTracks/MatTracking63.mat
- /InVivoTours/Result/MatTracks/MatTracking64.mat
- /InVivoTours/Result/MatTracks/MatTracking65.mat
- /InVivoTours/Result/MatTracks/MatTracking66.mat
- /InVivoTours/Result/MatTracks/MatTracking67.mat
- /InVivoTours/Result/MatTracks/MatTracking68.mat
- /InVivoTours/Result/MatTracks/MatTracking69.mat
- /InVivoTours/Result/MatTracks/MatTracking7.mat
- /InVivoTours/Result/MatTracks/MatTracking70.mat
- /InVivoTours/Result/MatTracks/MatTracking71.mat
- /InVivoTours/Result/MatTracks/MatTracking72.mat
- /InVivoTours/Result/MatTracks/MatTracking73.mat
- /InVivoTours/Result/MatTracks/MatTracking74.mat
- /InVivoTours/Result/MatTracks/MatTracking75.mat
- /InVivoTours/Result/MatTracks/MatTracking76.mat
- /InVivoTours/Result/MatTracks/MatTracking77.mat
- /InVivoTours/Result/MatTracks/MatTracking78.mat
- /InVivoTours/Result/MatTracks/MatTracking79.mat
- /InVivoTours/Result/MatTracks/MatTracking8.mat
- /InVivoTours/Result/MatTracks/MatTracking80.mat
- /InVivoTours/Result/MatTracks/MatTracking9.mat
- /InVivoTours/Fig6ac_InVivoTours_RSULM.m
- /InVivoTours/Fig6db_InVivoTours_SRPCA.m
- /InVivoTours/psf10Tours.mat
- /InVivoTours/psfSR4_Tours.mat
- /InVivoTours/SRPCATours.m
- /InVivoTours/ULM_localization2DSR_Tours.m
parent f98fb140
No related branches found
No related tags found
No related merge requests found
Showing
with 286 additions and 0 deletions
File added
File added
File added
File added
File added
File added
File added
File added
File added
File added
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
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
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment