-
Notifications
You must be signed in to change notification settings - Fork 0
/
contrack_score.m
112 lines (97 loc) · 3.96 KB
/
contrack_score.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
function [ scores ] = contrack_score(fg, dt6, fib2voxXform, dwiSeg, dwiROI )
% Runs the Contrack score algorithm to rate a set of fiber tracts
% [ scores ] = contrack_score( fg, dt6, dwiSeg, dwiROI )
%
% Inputs:
% fg : A fiber group, which essentially summarizes a pdb file.
% The pdb file is imported as fg.fibers (struct with 3xn-paths)
% NOTE : This is typically loaded from a pdb file
% >> fg = dtiLoadFiberGroup('fname.pdb')
%
% dt6 : The diffusion weighted data
%
% dwiSeg : A co-registered segmented volume. Could be higher res than the
% diffusion data. Contains the ROI.
%
% dwiROI : An ROI to avoid (passing this will give the fiber a zero
% score). A matrix of size dwiSeg, where 0 indicates region to
% avoid and 1 indicates acceptable region.
%
% Outuputs:
% score : A vector of n * 1, with scores for each fiber, ordered using
% the same numerical index as the fg.fibers input.
% HISTORY:
% 2012.12.05 SM: wrote it.
% Count the fibers. Each will be scored independently
n_fibers = length(fg.fibers);
% Size the scores
scores = zeros(n_fibers,1);
% Compute the Bingham distribution constants
% NOTE TODO : Get these from somewhere.
C =1; % Fix this
CL = 1; % Fix this
sigmaM = pi*14/180; % User param. From paper (pg. 7 col. 2, para 1)
eta = .175; % User param. From paper (pg. 7 col. 2, para 2)
% Compute the Watson distribution constants
sigmaC = 1; % Angular dispersion
CW = 1; % Normalizing constant
lambda = exp(-2); % User length scoring param. From paper (pg. 7 col. 2, para 3)
angleCutoff = 2.26;% radians = 129.488462 degrees
for f_ctr=1:n_fibers,
% Algorithm to compute the score :
% Q(s) = p(D|s) p(s)
% s = Estimated pathway
% D = Raw diffusion data
% Get the raw fiber data (xyz tangent vectors along a trajectory)
fgf = fg.fibers{1}; % size = 3xfiber-len
fgftan = diff(fg.fibers{1}')'; % size = 3xfiber-len - 1
% Get the diffusion tensors along the fiber path
tensors = ctrExtractDWITensorsAlongPath(fgf,dt6,fib2voxXform);
% Stage 1: Compute p(D|s) = Π_{i=1:n} [ p(Di | ti) ]
% Di are diffusion tensors at each point along the pathway
% ti are tangent vectors at each point along the pathway
pds = 1;
%For each tangent vector along the length of the fiber
for j=1:size(fgftan,2),
% Compute the score for the tangent
pdt = ctrBinghamScore(fgftan(:,j), tensors{j+1}.D, C, CL, sigmaM, eta );
% Multiply the estimates for this point
pds = pds * pdt;
% No need to loop if the score goes to zero somewhere in the middle
if(pds == 0) scores(f_ctr) = 0; break; end;
end
if(pds == 0) scores(f_ctr) = 0; continue; end;
% Stage 2: Compute p(s) = pend(s1)*pend(sn) Π_{i=1:n} [ p(Di | ti) ]
% Get the roi values along the fiber path (used in computing ps).
roi = ctrExtractROIAlongPath(fgf, dwiROI, fib2voxXform);
% Make sure the fiber starts and ends in the ROI.
pends1 = roi(1);
pendsn = roi(end);
% Initialize ps
ps = pends1 * pendsn;
% No need to loop if the score goes to zero somewhere in the middle
if(ps == 0) scores(f_ctr) = 0; continue; end;
%For each edge of the fiber
for j=2:size(fgf,2)-1,
% Local angle between two edges of the tract
segPre = (fgf(:,j) - fgf(:,j-1));
segPre = segPre./norm(segPre);
segPost = (fgf(:,j+1) - fgf(:,j));
segPost = segPost./norm(segPost);
% Inv cos of the dot product
thetaSeg = acos(segPost'*segPre);
% Apply a cut-off angle
if( thetaSeg > angleCutoff ) ps = 0; break; end;
% Compute the watson score
pcurve = ctrWatsonScore(CW, sigmaC, thetaSeg);
% Compute the manual knob and wm mask
plen = lambda * roi(j); % If the ROI is zero, this will exit.
% Multiply the estimates for this point
ps = ps * pcurve * plen;
% No need to loop if the score goes to zero somewhere in the middle
if(ps == 0) break; end;
end
% Stage 3: Compute the score Q(s)
scores(f_ctr) = pds * ps;
end
end