File: brain1020.m

package info (click to toggle)
octave-brain2mesh 0.7.9-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, sid
  • size: 11,736 kB
  • sloc: makefile: 54
file content (373 lines) | stat: -rw-r--r-- 15,268 bytes parent folder | download
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
function [landmarks, curves, initpoints]=brain1020(node, face, initpoints, perc1, perc2, varargin)
%
% landmarks=brain1020(node, face)
%   or
% landmarks=brain1020(node, face, [], perc1, perc2)
% landmarks=brain1020(node, face, initpoints)
% landmarks=brain1020(node, face, initpoints, perc1, perc2)
% [landmarks, curves, initpoints]=brain1020(node, face, initpoints, perc1, perc2, options)
%
% compute 10-20-like scalp landmarks with user-specified density on a head mesh
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% == Input ==
%    node: full head mesh node list
%    face: full head mesh element list- a 3-column array defines face list
%          for the exterior (scalp) surface; a 4-column array defines the
%          tetrahedral mesh of the full head.
%    initpoints:(optional) one can provide the 3-D coordinates of the below
%          5 landmarks: nz, iz, lpa, rpa, cz0 (cz0 is the initial guess of cz)
%          initpoints can be a struct with the above landmark names
%          as subfield, or a 5x3 array definining these points in the above
%          mentioned order (one can use the output landmarks as initpoints)
%    perc1:(optional) the percentage of geodesic distance towards the rim of 
%          the landmarks; this is the first number of the 10-20 or 10-10 or
%          10-5 systems, in this case, it is 10 (for 10%). default is 10.
%    perc2:(optional) the percentage of geodesic distance twoards the center 
%          of the landmarks; this is the 2nd number of the 10-20 or 10-10 or
%          10-5 systems, which are 20, 10, 5, respectively, default is 20
%    options: one can add additional 'name',value pairs to the function
%          call to provide additional control. Supported optional names
%          include
%           'display' : [1] or 0, if set to 1, plot landmarks and curves
%           'cztol' : [1e-6], the tolerance for searching cz that bisects
%                   saggital and coronal reference curves
%           'maxcziter' : [10] the maximum number of iterations to update
%                   cz to bisect both cm and sm curves
%           'baseplane' : [1] or 0, if set to 1, create the reference
%                   curves along the primary control points (nz,iz,lpa,rpa)
%           'minangle' : [0] if set to a positive number, this specifies
%                   the minimum angle (radian) between adjacent segments in
%                   the reference curves to avoid sharp turns (such as the
%                   dips near ear canals), this parameter will be
%                   passed to polylinesimplify to simplify the curve first.
%                   Please be noted that the landmarks generated with
%                   simplified curves may not land exactly on the surface.
%
% == Output ==
%    landmarks: a structure storing all computed landmarks. The subfields
%          include two sections: 
%          1) 'nz','iz','lpa','rpa','cz': individual 3D positions defining
%             the 5 principle reference points: nasion (nz), inion (in),
%             left-pre-auricular-point (lpa), right-pre-auricular-point 
%             (rpa) and vertex (cz) - cz is updated from initpoints to bisect
%             the saggital and coronal ref. curves.
%          2) landmarks along specific cross-sections, each cross section
%             may contain more than 1 position. The cross-sections are
%             named in the below format:
%             2.1: a fieldname starting from 'c', 's' and 'a' indicates
%                  the cut is along coronal, saggital and axial directions,
%                  respectively;
%             2.2: a name starting from 'pa' indicates the cut is along the
%                  axial plane acrossing principle reference points
%             2.3: the following letter 'm', 'a','p' suggests the 'medial',
%                  'anterior' and 'posterior', respectively
%             2.4: the last letter 'l' or 'r' suggests the 'left' and
%                  'right' side, respectively
%             2.5: non-medial coronal cuts are divided into two groups, the
%                  anterior group (ca{l,r}) and the posterior group
%                  (cp{lr}), with a number indicates the node spacing 
%                  stepping away from the medial plane.
%
%             for example, landmarks.cm refers to the landmarks along the
%                  medial-coronal plane, anterior-to-posteior order
%             similarly, landmarks.cpl_3 refers to the landmarks along the
%                  coronal (c) cut plane located in the posterior-left side
%                  of the head, with 3 saggital landmark spacing from the
%                  medial-coronal reference curve.
%    curves: a structure storing all computed cross-section curves. The
%             subfields are named similarly to landmarks, except that
%             landmarks stores the 10-? points, and curves stores the
%             detailed cross-sectional curves
%    initpoints: a 5x3 array storing the principle reference points in the
%             orders of 'nz','iz','lpa','rpa','cz'
%
% 
% == Example ==
% See brain2mesh/examples/SPM_example_brain.m for an example
% https://github.com/fangq/brain2mesh/blob/master/examples/SPM_example_brain.m
%
% == Dependency ==
% This function requires a pre-installed Iso2Mesh Toolbox
%  Download URL: http://github.com/fangq/iso2mesh
%  Website: http://iso2mesh.sf.net
%
% == Reference ==
% If you use this function in your publication, the authors of this toolbox
% apprecitate if you can cite the below paper
%
%  Anh Phong Tran, Shijie Yan and Qianqian Fang, "Improving model-based
%  fNIRS analysis using mesh-based anatomical and light-transport models,"
%  Neurophotonics, 7(1), 015008, 2020 
%  URL: http://dx.doi.org/10.1117/1.NPh.7.1.015008
%
%
% -- this function is part of brain2mesh toolbox (http://mcx.space/brain2mesh)
%    License: GPL v3 or later, see LICENSE.txt for details
%

if(nargin<2)
    error('one must provide a head-mesh to call this function');
end

if(isempty(node) || isempty(face) || size(face,2)<=2 || size(node,2)<3)
    error('input node must have 3 columns, face must have at least 3 columns');
end

if(nargin<5)
    perc2=20;
    if(nargin<4)
        perc1=10;
        if(nargin<3)
            initpoints=[];
        end
    end
end

% parse user options
opt=varargin2struct(varargin{:});

showplot=jsonopt('display',1,opt);
baseplane=jsonopt('baseplane',1,opt);
tol=jsonopt('cztol',1e-6,opt);
dosimplify=jsonopt('minangle',0,opt);
maxcziter=jsonopt('maxcziter',10,opt);

% convert initpoints input to a 5x3 array
if(isstruct(initpoints))
    initpoints=struct('nz', initpoints.nz(:).','iz', initpoints.iz(:).',...
                     'lpa',initpoints.lpa(:).','rpa',initpoints.rpa(:).',...
		     'cz', initpoints.cz(:).');
    landmarks=initpoints;
    if(exist('struct2array','file'))
        initpoints=struct2array(initpoints);
        initpoints=reshape(initpoints(:),3,length(initpoints(:))/3)';
    else
        initpoints=[initpoints.nz(:).'; initpoints.iz(:).';initpoints.lpa(:).';initpoints.rpa(:).';initpoints.cz(:).'];
    end
end

% convert tetrahedral mesh into a surface mesh
if(size(face,2)>=4)
    face=volface(face(:,1:4));
end

% remove nodes not located in the surface
[node,face]=removeisolatednode(node,face);

% if initpoints is not sufficient, ask user to interactively select nz, iz, lpa, rpa and cz first
if(isempty(initpoints) || size(initpoints,1)<5)
    hf=figure;
    plotmesh(node,face);
    set(hf,'userdata',initpoints);
    if(~isempty(initpoints))
        hold on;
        plotmesh(initpoints,'gs', 'LineWidth',4);
    end
    idx=size(initpoints,1)+1;
    landmarkname={'Nasion','Inion','Left-pre-auricular-point','Right-pre-auricular-point','Vertex/Cz','Done'};
    title(sprintf('Rotate the mesh, select data cursor, click on P%d: %s',idx, landmarkname{idx}));
    rotate3d('on');
    set(datacursormode(hf),'UpdateFcn',@myupdatefcn);
end

% wait until all 5 points are defined
if(exist('hf','var'))
    try
        while(size(get(hf,'userdata'),1)<5)
            pause(0.1);
        end
    catch
        error('user aborted');
    end
    datacursormode(hf,'off');
    initpoints=get(hf,'userdata');
    close(hf);
end

if(showplot)
    disp(initpoints);
end

% save input initpoints to landmarks output, cz is not finalized
if(size(initpoints,1)>=5)
    landmarks=struct('nz', initpoints(1,:),'iz', initpoints(2,:),...
                     'lpa',initpoints(3,:),'rpa',initpoints(4,:),...
		     'cz', initpoints(5,:));
end

% at this point, initpoints contains {nz, iz, lpa, rpa, cz0}
% plot the head mesh
if(showplot)
    figure;
    hp=plotmesh(node,face,'facealpha',0.6,'facecolor',[1 0.8 0.7]);
    if(~isoctavemesh)
        set(hp,'linestyle','none');
    end
    camlight;
    lighting gouraud
    hold on;
end

lastcz=[1 1 1]*inf;
cziter=0;

%% Find cz that bisects cm and sm curves within a tolerance, using UI 10-10 approach
while(norm(initpoints(5,:)-lastcz)>tol && cziter<maxcziter)
    %% Step 1: nz, iz and cz0 to determine saggital reference curve
    
    nsagg=slicesurf(node, face, initpoints([1,2,5],:));
    %% Step 1.1: get cz1 as the mid-point between iz and nz
    [slen, nsagg]=polylinelen(nsagg, initpoints(1,:), initpoints(2,:), initpoints(5,:));
    if(dosimplify)
        [nsagg, slen]=polylinesimplify(nsagg,dosimplify);
    end
    [idx, weight, cz]=polylineinterp(slen, sum(slen)*0.5, nsagg);
    initpoints(5,:)=cz(1,:);

    %% Step 1.2: lpa, rpa and cz1 to determine coronal reference curve, update cz1
    curves.cm=slicesurf(node, face, initpoints([3,4,5],:));
    [len, curves.cm]=polylinelen(curves.cm, initpoints(3,:), initpoints(4,:), initpoints(5,:));
    if(dosimplify)
        [curves.cm, len]=polylinesimplify(curves.cm,dosimplify);
    end
    [idx, weight, coro]=polylineinterp(len, sum(len)*0.5, curves.cm);
    lastcz=initpoints(5,:);
    initpoints(5,:)=coro(1,:);
    cziter=cziter+1;
    if(showplot)
        fprintf('cz iteration %d error %e\n',cziter, norm(initpoints(5,:)-lastcz));
    end
end

% set the finalized cz to output
landmarks.cz=initpoints(5,:);

if(showplot)
    disp(initpoints);
end

%% Step 2: subdivide saggital (sm) and coronal (cm) ref curves

[idx, weight, coro]=polylineinterp(len, sum(len)*(perc1:perc2:(100-perc1))*0.01, curves.cm);
landmarks.cm=coro;                 % t7, c3, cz, c4, t8

curves.sm=slicesurf(node, face, initpoints([1,2,5],:));
[slen, curves.sm]=polylinelen(curves.sm, initpoints(1,:), initpoints(2,:), initpoints(5,:));
if(dosimplify)
    [curves.sm, slen]=polylinesimplify(curves.sm,dosimplify);
end
[idx, weight, sagg]=polylineinterp(slen, sum(slen)*(perc1:perc2:(100-perc1))*0.01, curves.sm);
landmarks.sm=sagg;           % fpz, fz, cz, pz, oz

%% Step 3: fpz, t7 and oz to determine left 10% axial reference curve

[landmarks.aal, curves.aal, landmarks.apl, curves.apl]=slicesurf3(node,face,landmarks.sm(1,:), landmarks.cm(1,:), landmarks.sm(end,:),perc2*2);

%% Step 4: fpz, t8 and oz to determine right 10% axial reference curve

[landmarks.aar, curves.aar, landmarks.apr, curves.apr]=slicesurf3(node,face,landmarks.sm(1,:), landmarks.cm(end,:),landmarks.sm(end,:), perc2*2);


%% show plots of the landmarks
if(showplot)
    plotmesh(curves.sm,'r-','LineWidth',1);
    plotmesh(curves.cm,'g-','LineWidth',1);
    plotmesh(curves.aal,'k-','LineWidth',1);
    plotmesh(curves.aar,'k-','LineWidth',1);
    plotmesh(curves.apl,'b-','LineWidth',1);
    plotmesh(curves.apr,'b-','LineWidth',1);

    plotmesh(landmarks.sm,'ro','LineWidth',2);
    plotmesh(landmarks.cm,'go','LineWidth',2);
    plotmesh(landmarks.aal,'ko','LineWidth',2);
    plotmesh(landmarks.aar,'mo','LineWidth',2);
    plotmesh(landmarks.apl,'ko','LineWidth',2);
    plotmesh(landmarks.apr,'mo','LineWidth',2);
end

%% Step 5: computing all anterior coronal cuts, moving away from the medial cut (cm) toward frontal

idxcz=closestnode(landmarks.sm,landmarks.cz);

skipcount=max(floor(10/perc2),1);

for i=1:size(landmarks.aal,1)-skipcount
    step=(perc2*25)*0.1*(1+((perc2<20 + perc2<10) && i==size(landmarks.aal,1)-skipcount));
    [landmarks.(sprintf('cal_%d',i)), leftpart, landmarks.(sprintf('car_%d',i)), rightpart]=slicesurf3(node,face,landmarks.aal(i,:), landmarks.sm(idxcz-i,:), landmarks.aar(i,:),step);
    if(showplot)
        plotmesh(leftpart,'k-','LineWidth',1);
        plotmesh(rightpart,'k-','LineWidth',1);

        plotmesh(landmarks.(sprintf('cal_%d',i)),'yo','LineWidth',2);
        plotmesh(landmarks.(sprintf('car_%d',i)),'co','LineWidth',2);
    end
end

%% Step 6: computing all posterior coronal cuts, moving away from the medial cut (cm) toward occipital

for i=1:size(landmarks.apl,1)-skipcount
    step=(perc2*25)*0.1*(1+((perc2<20 + perc2<10) && i==size(landmarks.apl,1)-skipcount));
    [landmarks.(sprintf('cpl_%d',i)), leftpart, landmarks.(sprintf('cpr_%d',i)), rightpart]=slicesurf3(node,face,landmarks.apl(i,:), landmarks.sm(idxcz+i,:), landmarks.apr(i,:),step);
    if(showplot)
        plotmesh(leftpart,'k-','LineWidth',1);
        plotmesh(rightpart,'k-','LineWidth',1);

        plotmesh(landmarks.(sprintf('cpl_%d',i)),'yo','LineWidth',2);
        plotmesh(landmarks.(sprintf('cpr_%d',i)),'co','LineWidth',2);
    end
end

%% Step 7: create the axial cuts across priciple ref. points: left: nz, lpa, iz, right: nz, rpa, iz

if(baseplane && perc2<=10)
    [landmarks.paal, curves.paal, landmarks.papl, curves.papl]=slicesurf3(node,face,landmarks.nz, landmarks.lpa, landmarks.iz, perc2*2);
    [landmarks.paar, curves.paar, landmarks.papr, curves.papr]=slicesurf3(node,face,landmarks.nz, landmarks.rpa, landmarks.iz, perc2*2);
    if(showplot)
        plotmesh(curves.paal,'k-','LineWidth',1);
        plotmesh(curves.paar,'k-','LineWidth',1);
        plotmesh(curves.papl,'k-','LineWidth',1);
        plotmesh(curves.papr,'k-','LineWidth',1);

        plotmesh(landmarks.paal,'yo','LineWidth',2);
        plotmesh(landmarks.papl,'co','LineWidth',2);
        plotmesh(landmarks.paar,'yo','LineWidth',2);
        plotmesh(landmarks.papr,'co','LineWidth',2);
    end
end

%%-------------------------------------------------------------------------------------
% helper functions
%--------------------------------------------------------------------------------------

% the respond function when a data-cursor tip to popup

function txt=myupdatefcn(empt,event_obj)
pt=get(gcf,'userdata');
pos = get(event_obj,'Position');
%idx=  get(event_obj,'DataIndex');
txt = {['x: ',num2str(pos(1))],...
       ['y: ',num2str(pos(2))],['z: ',num2str(pos(3))]};
if(~isempty(pt) && ismember(pos,pt,'rows'))
    return;
end
targetup=get(get(event_obj,'Target'),'parent');
idx=size(pt,1)+2;
landmarkname={'Nasion','Inion','Left-pre-auricular-point','Right-pre-auricular-point','Vertex/Cz','Done'};
title(sprintf('Rotate the mesh, select data cursor, click on P%d: %s',idx, landmarkname{idx}));
set(targetup,'userdata',struct('pos',pos));
pt=[pt;pos];
if(size(pt,1)<6)
     set(gcf,'userdata',pt);
end
hpt=findobj(gcf,'type','line');
if(~isempty(hpt))
    set(hpt,'xdata',pt(:,1),'ydata',pt(:,2),'zdata',pt(:,3));
else
    hold on;
    plotmesh([pt;pos],'gs', 'LineWidth',4);
end
disp(['Adding landmark ' landmarkname{idx-1} ':' txt]);
datacursormode(gcf,'off');
rotate3d('on');