Skip to content

Commit 0098053

Browse files
delucaaldelucaal
committed
Last stable update before major changes
Fixes to the automatic detection of gradient directions Added a function to automatically determine the correct spatial orientation Several bug fixes
1 parent 692357a commit 0098053

File tree

22 files changed

+2101
-10859
lines changed

22 files changed

+2101
-10859
lines changed

‎ExploreDTIInterface/Auxiliaries/CSD.m‎

Lines changed: 1 addition & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -1,140 +1 @@
1-
%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%%%%% Alberto De Luca - alberto@isi.uu.nl $%%%%%% Distributed under the terms of LGPLv3 %%%
2-
3-
4-
5-
classdef CSD
6-
% Originally written from Ben Jeurissen (ben.jeurissen@uantwerpen.be)
7-
% under the supervision of Alexander Leemans (a.leemans@umcutrecht.nl)
8-
%
9-
% Constrained Spherical Deconvolution (CSD) class
10-
%
11-
% Usage:
12-
% r_sh = SD.response(dwi, grad, bval, minFA, lmax);
13-
% csd = CSD(r_sh, target_lmax);
14-
% fod_sh = csd.deconv(dwi_sh);
15-
%
16-
% for super-resolution: make sure target_lmax is higher than SH order of r_sh
17-
%
18-
% default parameters can be changed after construction by:
19-
% csd.setLambda(lambda);
20-
% csd.setTau(tau)
21-
% csd.setIter(iter)
22-
% csd.setInitLmax(init_lmax)
23-
% csd.setNegativeDirs(dirs)
24-
%
25-
% see also SD
26-
%
27-
% Copyright Ben Jeurissen (ben.jeurissen@ua.ac.be)
28-
%
29-
properties (GetAccess = public, SetAccess = private)
30-
lmax;
31-
end
32-
33-
properties (Access = private)
34-
n;
35-
convmat;
36-
hr_sh;
37-
lambda_;
38-
r_rh;
39-
my_fod;
40-
41-
lambda = 1;
42-
tau = 0.1;
43-
dirs = textread('dir300.txt');
44-
niter = 50;
45-
init_lmax = 4;
46-
end
47-
48-
methods (Access = public)
49-
function this = CSD(r_sh, target_lmax)
50-
if ~exist('target_lmax','var') || isempty(target_lmax)
51-
target_lmax = SH.n2lmax(size(r_sh,1));
52-
end
53-
this.lmax = target_lmax;
54-
r_sh_n = size(r_sh,1);
55-
r_sh_lmax = SH.n2lmax(r_sh_n);
56-
target_n = SH.lmax2n(target_lmax);
57-
58-
lmax_ = min(r_sh_lmax,target_lmax);
59-
this.n = min(r_sh_n,target_n);
60-
61-
d_sh = SH.eval(lmax_,[0 0])';
62-
k = find(d_sh);
63-
this.r_rh = r_sh(k)./d_sh(k);
64-
m = [];
65-
for l = 0:2:lmax_
66-
m = [m; ones(2*l+1,1)*this.r_rh(l/2+1)];
67-
end
68-
this.convmat = diag(m);
69-
70-
this.hr_sh = SH.eval(this.lmax,c2s(this.dirs));
71-
this.convmat = [this.convmat zeros(size(this.convmat,1),size(this.hr_sh,2)-size(this.convmat,2)) ];
72-
this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1);
73-
end
74-
75-
function this = setLambda(this, lambda)
76-
this.lambda = lambda;
77-
this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1);
78-
end
79-
80-
function this = setTau(this, tau)
81-
this.tau = tau;
82-
end
83-
84-
function this = setIter(this, iter)
85-
this.niter = iter;
86-
end
87-
88-
function this = setInitLmax(this, init_lmax)
89-
this.init_lmax = init_lmax;
90-
end
91-
92-
function this = setNegativeDirs(this, dirs)
93-
this.dirs = dirs;
94-
this.hr_sh = SH.eval(this.lmax,c2s(this.dirs));
95-
this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1);
96-
end
97-
98-
function fod_sh = process(this, dwi_sh)
99-
fod_sh = this.deconv(dwi_sh);
100-
end
101-
102-
function dwi_sh = conv(this, fod_sh)
103-
if ndims(fod_sh) ~= 2; [fod_sh, mask] = vec(fod_sh); end;
104-
dwi_sh = this.convmat*fod_sh;
105-
if exist('mask','var'); dwi_sh = unvec(dwi_sh,mask); end;
106-
end
107-
108-
function fod_sh = deconv(this, dwi_sh)
109-
if ndims(dwi_sh) ~= 2; [dwi_sh, mask] = vec(dwi_sh); end;
110-
fod_sh = this.convmat\dwi_sh(1:this.n,:); fod_sh(SH.lmax2n(this.init_lmax)+1:end,:) = 0;
111-
threshold = this.tau*mean(this.hr_sh*fod_sh,1);
112-
fod_sh = CSD.actual_deconv(fod_sh, dwi_sh(1:this.n,:), this.convmat, this.hr_sh, threshold, this.lambda_, this.niter);
113-
if exist('mask','var'); fod_sh = unvec(fod_sh,mask); end;
114-
end
115-
end
116-
117-
methods (Access = private, Static = true)
118-
function f_sh = actual_deconv(f_sh,dwi_sh,fconv,hr_sh,threshold,lambda,niter)
119-
parfor i = 1:size(f_sh,2) % for each voxel
120-
f_sh_i = f_sh(:,i); dwi_sh_i = dwi_sh(:,i); threshold_i = threshold(i);
121-
for it = 1:niter % maximum of 50 CSD iterations
122-
f_hr = hr_sh*f_sh_i; % calculate high-res FOD
123-
neg = find (f_hr < threshold_i); % find negative (and small) directions in high-res FOD
124-
if size(neg,1) + size(fconv,1) < size(hr_sh,2) % if there aren't enough negative directions to allow super resolution
125-
disp('Warning: too few negative directions identified - failed to converge'); % stop
126-
break;
127-
end
128-
m = [fconv; lambda*hr_sh(neg,:)]; % assume the negative directions
129-
s = [dwi_sh_i; zeros(size(neg,1),1)]; % are zero
130-
f_sh_i = m\s; % re-estimate f_sh as if more directions are available
131-
ssd = sum((f_sh_i-f_sh(:,i)).^2);
132-
f_sh(:,i) = f_sh_i;
133-
if (ssd == 0.0) % if there is no improvement of f_sh over the previous iteration
134-
break; % we have converged
135-
end
136-
end
137-
end
138-
end
139-
end
140-
end
1+
%%%$ Included in MRIToolkit (https://github.com/delucaal/MRIToolkit) %%%%%% Alberto De Luca - alberto@isi.uu.nl $%%%%%% Distributed under the terms of LGPLv3 %%%classdef CSD% Originally written from Ben Jeurissen (ben.jeurissen@uantwerpen.be)% under the supervision of Alexander Leemans (a.leemans@umcutrecht.nl)% % Constrained Spherical Deconvolution (CSD) class % % Usage: % r_sh = SD.response(dwi, grad, bval, minFA, lmax); % csd = CSD(r_sh, target_lmax); % fod_sh = csd.deconv(dwi_sh); % % for super-resolution: make sure target_lmax is higher than SH order of r_sh % % default parameters can be changed after construction by: % csd.setLambda(lambda); % csd.setTau(tau) % csd.setIter(iter) % csd.setInitLmax(init_lmax) % csd.setNegativeDirs(dirs) % % see also SD % % Copyright Ben Jeurissen (ben.jeurissen@ua.ac.be) % properties (GetAccess = public, SetAccess = private) lmax; end properties (Access = private) n; convmat; hr_sh; lambda_; r_rh; my_fod; lambda = 1; tau = 0.1; dirs = textread('dir300.txt'); niter = 50; init_lmax = 4; end methods (Access = public) function this = CSD(r_sh, target_lmax) if ~exist('target_lmax','var') || isempty(target_lmax) target_lmax = SH.n2lmax(size(r_sh,1)); end this.lmax = target_lmax; r_sh_n = size(r_sh,1); r_sh_lmax = SH.n2lmax(r_sh_n); target_n = SH.lmax2n(target_lmax); lmax_ = min(r_sh_lmax,target_lmax); this.n = min(r_sh_n,target_n); d_sh = SH.eval(lmax_,[0 0])'; k = find(d_sh); this.r_rh = r_sh(k)./d_sh(k); m = []; for l = 0:2:lmax_ m = [m; ones(2*l+1,1)*this.r_rh(l/2+1)]; end this.convmat = diag(m); this.hr_sh = SH.eval(this.lmax,c2s(this.dirs)); this.convmat = [this.convmat zeros(size(this.convmat,1),size(this.hr_sh,2)-size(this.convmat,2)) ]; this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1); end function this = setLambda(this, lambda) this.lambda = lambda; this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1); end function this = setTau(this, tau) this.tau = tau; end function this = setIter(this, iter) this.niter = iter; end function this = setInitLmax(this, init_lmax) this.init_lmax = init_lmax; end function this = setNegativeDirs(this, dirs) this.dirs = dirs; this.hr_sh = SH.eval(this.lmax,c2s(this.dirs)); this.lambda_ = this.lambda*size(this.convmat,1)*this.r_rh(1)/size(this.hr_sh,1); end function fod_sh = process(this, dwi_sh) fod_sh = this.deconv(dwi_sh); end function dwi_sh = conv(this, fod_sh) if ndims(fod_sh) ~= 2; [fod_sh, mask] = vec(fod_sh); end; dwi_sh = this.convmat*fod_sh; if exist('mask','var'); dwi_sh = unvec(dwi_sh,mask); end; end function fod_sh = deconv(this, dwi_sh) if ndims(dwi_sh) ~= 2; [dwi_sh, mask] = vec(dwi_sh); end; fod_sh = this.convmat\dwi_sh(1:this.n,:); fod_sh(SH.lmax2n(this.init_lmax)+1:end,:) = 0; threshold = this.tau*mean(this.hr_sh*fod_sh,1); fod_sh = CSD.actual_deconv(fod_sh, dwi_sh(1:this.n,:), this.convmat, this.hr_sh, threshold, this.lambda_, this.niter); if exist('mask','var'); fod_sh = unvec(fod_sh,mask); end; end end methods (Access = private, Static = true) function f_sh = actual_deconv(f_sh,dwi_sh,fconv,hr_sh,threshold,lambda,niter) parfor i = 1:size(f_sh,2) % for each voxel f_sh_i = f_sh(:,i); dwi_sh_i = dwi_sh(:,i); threshold_i = threshold(i); for it = 1:niter % maximum of 50 CSD iterations f_hr = hr_sh*f_sh_i; % calculate high-res FOD neg = find (f_hr < threshold_i); % find negative (and small) directions in high-res FOD if size(neg,1) + size(fconv,1) < size(hr_sh,2) % if there aren't enough negative directions to allow super resolution disp('Warning: too few negative directions identified - failed to converge'); % stop break; end m = [fconv; lambda*hr_sh(neg,:)]; % assume the negative directions s = [dwi_sh_i; zeros(size(neg,1),1)]; % are zero f_sh_i = m\s; % re-estimate f_sh as if more directions are available ssd = sum((f_sh_i-f_sh(:,i)).^2); f_sh(:,i) = f_sh_i; if (ssd == 0.0) % if there is no improvement of f_sh over the previous iteration break; % we have converged end end end end endend

0 commit comments

Comments
 (0)