Skip to content

Commit

Permalink
initial
Browse files Browse the repository at this point in the history
  • Loading branch information
Beliy Nikita committed Jul 4, 2024
0 parents commit bca86d1
Show file tree
Hide file tree
Showing 52 changed files with 5,078 additions and 0 deletions.
191 changes: 191 additions & 0 deletions +pet_IDIF/IFextractionPET.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,191 @@
function IDIF = IFextractionPET(pet_files, ...
SubsetCaro, MaskCaro, ...
tFrames, ...
varargin)

% IFextractionPET is a function implementing the Pairwise correlation
% method to compute the image derived input function based on the PET
% following the method devlopped by M.Schain [Schain, J. Cereb Blood
% Flow Metab, 2013].
%
% Input:
% pet_files: Dynamic PET volumes, either in list (3D images) or one 4D
% SubsetCaro: binary mask of regions including carodtid arteries
% MaskCaro: binary mask of seeds representing hot part of carotid arteries
% tFrames: cell array of mid-times for each frame of pet_files
% varargin:
% PSF: PSF of the medical imaging system for the partial volume correction,
% (default: 6)
% Tcorr: Time in seconds during wich the PET data are corrected
% for spillover correction (default: 120)
% thresh: Threshold for correlation, used in Shain's method
% (default: 0.97)
%
% Output:
% Extracted carotid: specified_name.nii
% Estimated IF: specified_name.txt
%
% Variable:
% D = data from the PET
% Dsub = subset of the data from the PET
% DsubN = centered and normed subset of the data from the PET
% Mask1 = mask from the lower part of the brain
% Mask2 = mask of the carotid (from coregistred MRI)
% Caro = mask of the carotid from Schain's algorithm
% Caro3D = 3D binary mask of the carotid from Schain's algorithm
% Caro3Dsmooth = smoothed image of Caro3D
% tFrame = time of the frame
% IF = input function
%
% =========================================================================
%
% M.A. Bahri, 2021.
% Cyclotron Research Centre, University of Liege, Belgium
% =========================================================================

args = inputParser;
args.addParameter('PSF', 6);
args.addParameter('Tcorr', 120);
args.addParameter('treshold', 0.97);
args.parse(varargin{:});


PSF = args.Results.PSF;
Tcorr = args.Results.Tcorr;
thresh = args.Results.treshold;

% Loading volumes
V = crc_read_spm_vol(pet_files);
in_pth = fileparts(V(1).fname);
bids_name = bids.File(V(1).fname);

if size(V, 1) < 5
error('IF extraction requires at least %d images', 5);
end

if numel(tFrames) ~= numel(V)
error('Number of mid-times must be same as number of volumes');
end

% Loading masks
% Old name Mask1
if ischar(SubsetCaro)
V_mask_subset = spm_vol(SubsetCaro);
mask_subset = spm_read_vols(V_mask_subset);
else
mask_subset = SubsetCaro;
end

% Old name Mask2
if ischar(MaskCaro)
V_mask_seed = spm_vol(MaskCaro);
mask_seed = spm_read_vols(V_mask_seed);
else
mask_seed = MaskCaro;
end

% generating output path and basename
out_pth = fullfile(in_pth, '..', 'IDIF');
if ~exist(out_pth, 'dir')
mkdir(out_pth);
end

% Variable attribution and Data initialization
for ind = 1:numel(V)
dat = spm_read_vols(V(ind));
D(ind, :) = dat(:);
end

% Pixel spacing of the medical imaging system for the partial volume
% correction
pixelspacing = [abs(V(1).mat(1, 1)) V(1).mat(2,2) V(1).mat(3,3)];

% =========================================================================
%% Begining of the Schain's algorithm
fprintf('-> Applying partial volume correction...\n');
% Extraction of the image subset out of the PET data
% id_subset == idx1
% id_seed == idx2
id_subset = find(mask_subset);
id_seed = find(mask_seed(id_subset));
Dsub = D(:, id_subset);

% Calculation of the correlation matrix
[t, nbv]= size(Dsub);
DsubC = Dsub - ones(t, 1) * mean(Dsub); % recentering the value
DsubN = zeros(t, nbv); %initialisation of the normalized matrix
for i = 1:nbv
DsubN(:, i) = DsubC(:, i) ./ norm(DsubC(:, i));
end

% Extraction of the correlation matrix only for the voxel of interrest
M = DsubN(:, id_seed)' * DsubN;
% All autocorrelation coefficients are set to zero
for i = 1:length(id_seed)
M(i, id_seed(i)) = 0;
end

% Extraction of voxels part of the carotid
[~, B] = find(M > thresh);
C = union(B, B);
fprintf('Voxels limit: %f\n', 1.1 * length(id_seed));
while size(C, 1) < 1.1 * length(id_seed)
thresh = thresh - 0.025;
[~, B] = find(M > thresh);
C = union(B, B);
fprintf('treshold = %f; size(C) = %d\n', thresh, size(C, 1));
if thresh < 0
warning(['Reached 0 treshold, ',...
'number of carotide voxels may be too small']);
break;
end
end

Caro = zeros(size(D, 2), 1);
Caro(id_subset(C)) = 1;

% =========================================================================
%% Partial Volume correction

% Transfromation of the extracted carotid from vector to 3D matrix
Caro3D = zeros(V(1).dim);
Caro3D(:)=Caro(:);

% Uncomment to see the extracted carotid with an other viewer than SPM
% figure;
% imshow3D(Caro3D)

% Writing tbe .nii file to verify if the extracted carotid is correct
bids_name.entities.desc = 'IDIF';
bids_name.suffix = 'mask';
write_image(V(1), Caro3D, fullfile(out_pth, bids_name.filename()));

% Computation of the GTM with one compartiment for the Tcorr first
% minutes
% Smoothing of the carotid mask
sig = fwhm2sigma(PSF); %sigma
sigma = [sig; sig; sig];
% Correction of the splillover effect and computation of IF
% idx3 = find(Caro3D);
Caro3Dsmooth = gauss3filter(Caro3D, sigma, pixelspacing);
idx3 = find(Caro3D);
%Caro3Dsmooth = Caro3Dsmooth(idx3);
tcorr = find(tFrames(:) > Tcorr, 1);

IF = zeros(numel(tFrames), 1);
for i = 1:numel(tFrames)
if i < tcorr
tempIF = D(i, :) ./ Caro3Dsmooth(:)';
IF(i) = mean(tempIF(idx3));
else
IF(i) = mean(D(i, idx3));
end
end

bids_name.entities.desc = '';
bids_name.suffix = 'IDIF';
bids_name.extension = '.tsv';
IDIF.fname = fullfile(out_pth, bids_name.filename());
IDIF.mid_time = tFrames(:);
IDIF.input_function = IF;
end
56 changes: 56 additions & 0 deletions +pet_IDIF/caro_correlate.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
function corr_mask = caro_correlate(data, core_mask, exp_mask, params)
def_params.threshold = 0.97;
def_params.expansion = 1.1;

fprintf('-> Selecting carotides by correlation\n');
params = crc_update_config(params, def_params);

% Reorienting data
dim = size(data);
n_vols = dim(4);
dim = dim(1:3);
D = permute(data, [4, 1, 2, 3]);

id_expanded = find(exp_mask);
id_core = find(core_mask(id_expanded));
Dexp = D(:, id_expanded);

% Calculation of the correlation matrix
[t, nbv]= size(Dexp);
DexpC = Dexp - ones(t, 1) * mean(Dexp); % recentering the value
DexpN = zeros(t, nbv);
for i = 1:nbv
DexpN(:, i) = DexpC(:, i) ./ norm(DexpC(:, i));
end

% Extraction of the correlation matrix only for the voxel of interrest
M = DexpN(:, id_core)' * DexpN;
% All autocorrelation coefficients are set to zero
for i = 1:length(id_core)
M(i, id_core(i)) = 0;
end

% Extraction of voxels part of the carotid
thresh = 1;
[~, B] = find(M > thresh);
C = union(B, B);
fprintf('---> Expansion limit: %.0f\n', params.expansion * length(id_core));
while true
[~, B] = find(M > thresh);
C = union(B, B);
fprintf('---> treshold = %f; size(C) = %d\n', thresh, size(C, 1));

thresh = thresh - 0.025;
if thresh < params.threshold
fprintf('--> Reached correlation threshold\n');
break;
end
if size(C, 1) >= params.expansion * length(id_core)
fprintf('--> Reached expansion threshold\n');
break;
end
end

corr_mask = zeros(dim, 'uint8');
corr_mask(id_expanded(C)) = 1;
end
88 changes: 88 additions & 0 deletions +pet_IDIF/caro_search.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
function [seed_mask, caro_mask, exp_mask] = caro_search(seed_candidates, caro_candidates,...
p_seed, p_expansion)
% Performs carotides search based on seeds candidates,
% and expanded carotides candidates
def_seed.min_value = 1000;
def_seed.candidate = 0.9;

def_expansion.hot_carotides = 0.8;
def_expansion.cold_carotides = 0.5;
def_expansion.min_size = 20;
def_expansion.max_size = 30;

fprintf('--> Analyzing seeds...\n');
p_seed = crc_update_config(p_seed, def_seed);
def_expansion = crc_update_config(p_expansion, def_expansion);

caro_mask = zeros(size(seed_candidates), 'uint8');
exp_mask = zeros(size(seed_candidates), 'uint8');
seed_mask = zeros(size(seed_candidates), 'uint8');

[val_cand, index_cand] = max(seed_candidates(:));
val_first = val_cand;

while (seed_candidates(index_cand) > 0)
if caro_candidates(index_cand) < p_seed.min_value
seed_candidates(index_cand) = 0;
[val_cand, index_cand] = max(seed_candidates(:));
continue;
end

[cand_x, cand_y, cand_z] = ind2sub(size(seed_candidates), index_cand);

% fprintf('---> Candidate %f at (%d, %d, %d): ',...
% val_cand, cand_x, cand_y, cand_z);
[caro, caro_exp] = expand_carotide(caro_candidates, index_cand,...
p_expansion.hot_carotides,...
p_expansion.cold_carotides,...
p_expansion.max_size);

% Masking tested seed
seed_candidates(index_cand) = 0;
% Masking expanded hot carotide
seed_candidates(caro_exp > 0) = 0;
% Masking expanded carotide
caro_candidates(caro_exp > 0) = 0;
if nnz(caro_exp) < p_expansion.min_size
% fprintf('---> Rejected due to size\n');
else
caro_mask = caro_mask | caro;
exp_mask = exp_mask | caro_exp;
seed_mask(index_cand) = 1;
end
[val_cand, index_cand] = max(seed_candidates(:));

end

n_caro = nnz(seed_mask);
if n_caro == 0
error('No valid carotides candidates found')
end
fprintf('Found %d carotides candidates of %d (core)/ %d (expanded)\n', ...
n_caro, nnz(caro_mask), nnz(exp_mask));

end

function [caro, caro_exp] = expand_carotide(img, seed,...
cutoff_core, cutoff_exp,...
max_size)
caro = zeros(size(img), 'uint8');
cutoff_mask = (img > (img(seed) * cutoff_core));
caro(seed) = 1;
SE = strel('sphere', 1);
caro_size = nnz(caro);
d_size = caro_size;
while d_size > 0 && caro_size < max_size
caro = imdilate(caro, SE) & cutoff_mask;
d_size = nnz(caro) - caro_size;
caro_size = caro_size + d_size;
end

cutoff_mask = (img > (img(seed) * cutoff_exp));
caro_exp = caro;
for i = 1:5
caro_exp = imdilate(caro_exp, SE) & cutoff_mask;
end

% fprintf('Expanded to %d (core) / %d (expanded) voxels\n', caro_size, nnz(caro_exp));
end
20 changes: 20 additions & 0 deletions +pet_IDIF/cilindrical_mask.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
function mask = get_geom_mask(mean_im, params)
% Generates cilindrical mask around COM of image

fprintf('--> Creating geometrical mask\n')
% default values
def_params.radius = 15; % Radius, in voxels of cilindrical mask
def_params.z_rejection = 10; % Number of top and bottom layer to ignore

params = crc_update_config(params, def_params);

[dim_x, dim_y, dim_z] = size(mean_im);
mask = zeros(dim_x, dim_y);
mask(com(mean_im, 1), com(mean_im, 2)) = 1;
mask = imdilate(mask, strel("disk", params.radius, 0));
mask = repmat(mask, 1, 1, dim_z);
% mask(:, :, 1:params.z_rejection) = 0;
mask(:, :, 1:com(mean_im, 3)) = 0;
mask(:, :, end-params.z_rejection: end) = 0;
mask = logical(mask);
end
Loading

0 comments on commit bca86d1

Please sign in to comment.