/
ft_volumerealign.m
2131 lines (1880 loc) · 84.1 KB
/
ft_volumerealign.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
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
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [realign, snap] = ft_volumerealign(cfg, mri, target)
% FT_VOLUMEREALIGN spatially aligns an anatomical MRI with head coordinates based on
% external fiducials or anatomical landmarks. This function typically does not change
% the anatomical MRI volume itself, but only adjusts the homogeneous transformation
% matrix that describes the mapping from voxels to the coordinate system. It also
% appends a coordsys-field to the output data, or it updates it. This field specifies
% how the x/y/z-axes of the coordinate system should be interpreted. Occasionally,
% the orientation and handedness of the output volume may be different from the orientation
% and handedness of the input volume. This is determined by the cfg.flip
% argument. See the code for more details.
%
% For spatial normalisation and deformation (i.e. warping) an MRI to a template brain
% you should use the FT_VOLUMENORMALISE function.
%
% Different methods for aligning the anatomical MRI to a coordinate system are
% implemented, which are described in detail below:
%
% INTERACTIVE - Use a graphical user interface to click on the location of anatomical
% landmarks or fiducials. The anatomical data can be displayed as three orthogonal
% MRI slices or as a rendering of the head surface. The coordinate system is updated
% according to the definition of the coordinates of these fiducials.
%
% FIDUCIAL - The coordinate system is updated according to the definition of the
% coordinates of anatomical landmarks or fiducials that are specified in the
% configuration.
%
% HEADSHAPE - Match the head surface from the MRI with a measured head surface using
% an iterative closest point procedure. The MRI will be updated to match the measured
% head surface. You can optionally do an initial manual coregistration of the two head
% surfaces.
%
% SPM - Align the individual MRI to the coordinate system of a target or template MRI
% by matching the two volumes.
%
% FSL - Align the individual MRI to the coordinate system of a target or template MRI
% by matching the two volumes.
%
% Use as
% [mri] = ft_volumerealign(cfg, mri)
% or
% [mri] = ft_volumerealign(cfg, mri, target)
% where the first input is the configuration structure, the second input is an
% anatomical or functional MRI volume and the third (optional) input is the the
% target anatomical MRI for SPM or FSL.
%
% The configuration can contain the following options
% cfg.method = string representing the method for aligning
% 'interactive' use the GUI to specify the fiducials
% 'fiducial' use pre-specified fiducials
% 'headshape' match the MRI surface to a headshape
% 'spm' match to template anatomical MRI
% 'fsl' match to template anatomical MRI
% cfg.coordsys = string specifying the origin and the axes of the coordinate
% system. Supported coordinate systems are 'ctf', '4d', 'bti',
% 'eeglab', 'neuromag', 'itab', 'yokogawa', 'asa', 'acpc',
% and 'paxinos'. See http://tinyurl.com/ojkuhqz
% cfg.clim = [min max], scaling of the anatomy color (default is automatic)
% cfg.parameter = 'anatomy' the parameter which is used for the visualization
% cfg.viewresult = string, 'yes' or 'no', whether or not to visualize aligned volume(s)
% after realignment (default = 'no')
% cfg.flip = string, 'yes' or 'no', to realign the volume approximately to the
% input coordinate axes, this may reorient the output volume relative
% to the input (default = 'yes', when cfg.method = 'interactive', and 'no' otherwise)
%
% When cfg.method = 'interactive', a user interface allows for the specification of
% the fiducials or landmarks using the mouse, cursor keys and keyboard. The fiducials
% can be specified by pressing the corresponding key on the keyboard (n/l/r or
% a/p/z). When pressing q the interactive mode will stop and the transformation
% matrix is computed. This method supports the following options:
% cfg.viewmode = 'ortho' or 'surface', visualize the anatomical MRI as three
% slices or visualize the extracted head surface (default = 'ortho')
% cfg.snapshot = 'no' ('yes'), making a snapshot of the image once a
% fiducial or landmark location is selected. The optional second
% output argument to the function will contain the handles to these
% figures.
% cfg.snapshotfile = 'ft_volumerealign_snapshot' or string, the root of
% the filename for the snapshots, including the path. If no path
% is given the files are saved to the pwd. The consecutive
% figures will be numbered and saved as png-file.
%
% When cfg.method = 'fiducial' and cfg.coordsys is based on external anatomical
% landmarks, as is common for EEG and MEG, the following is required to specify the
% voxel indices of the fiducials:
% cfg.fiducial.nas = [i j k], position of nasion
% cfg.fiducial.lpa = [i j k], position of LPA
% cfg.fiducial.rpa = [i j k], position of RPA
% cfg.fiducial.zpoint = [i j k], a point on the positive z-axis. This is
% an optional 'fiducial', and can be used to determine
% whether the input voxel coordinate axes are left-handed
% (i.e. flipped in one of the dimensions). If this additional
% point is specified, and the voxel coordinate axes are left
% handed, the volume is flipped to yield right handed voxel
% axes.
%
% When cfg.method = 'fiducial' and cfg.coordsys = 'acpc', as is common for fMRI,
% the following is required to specify the voxel indices of the fiducials:
% cfg.fiducial.ac = [i j k], position of anterior commissure
% cfg.fiducial.pc = [i j k], position of posterior commissure
% cfg.fiducial.xzpoint = [i j k], point on the midsagittal-plane with a
% positive Z-coordinate, i.e. an interhemispheric
% point above ac and pc
% The coordinate system will be according to the RAS_Tal convention, i.e.
% the origin corresponds with the anterior commissure the Y-axis is along
% the line from the posterior commissure to the anterior commissure the
% Z-axis is towards the vertex, in between the hemispheres the X-axis is
% orthogonal to the YZ-plane, positive to the right.
%
% When cfg.method = 'fiducial' and cfg.coordsys = 'paxinos' for a mouse brain,
% the following is required to specify the voxel indices of the fiducials:
% cfg.fiducial.bregma = [i j k], position of bregma
% cfg.fiducial.lambda = [i j k], position of lambda
% cfg.fiducial.yzpoint = [i j k], point on the midsagittal-plane
%
% With the 'interactive' and 'fiducial' methods it is possible to define an
% additional point (with the key 'z'), which should be a point on the positive side
% of the xy-plane, i.e. with a positive z-coordinate in world coordinates. This point
% will subsequently be used to check whether the input coordinate system is left or
% right-handed. For the 'interactive' method you can also specify an additional
% control point (with the key 'r'), that should be a point with a positive coordinate
% on the left-right axis, i.e.', a point on the right of the head.
%
% When cfg.method = 'headshape', the function extracts the scalp surface from the
% anatomical MRI, and aligns this surface with the user-supplied headshape.
% Additional options pertaining to this method should be defined in the subcfg
% cfg.headshape. The following option is required:
% cfg.headshape.headshape = string pointing to a headshape structure or a
% file containing headshape, see FT_READ_HEADSHAPE
%
% Additional options pertaining to the headshape method should be specified in
% the sub-structure cfg.headshape and can include:
% cfg.headshape.scalpsmooth = scalar, smoothing parameter for the scalp
% extraction (default = 2)
% cfg.headshape.scalpthreshold = scalar, threshold parameter for the scalp
% extraction (default = 0.1)
% cfg.headshape.interactive = 'yes' or 'no', use interactive realignment to
% align headshape with scalp surface (default = 'yes')
% cfg.headshape.icp = 'yes' or 'no', use automatic realignment
% based on the icp-algorithm. If both 'interactive'
% and 'icp' are executed, the icp step follows the
% interactive realignment step (default = 'yes')
%
% When cfg.method = 'spm', a third input argument is required. The input volume is
% coregistered to this target volume, using SPM. You can specify the version of
% the SPM toolbox to use with
% cfg.spmversion = string, 'spm2', 'spm8', 'spm12' (default = 'spm12')
%
% Additional options pertaining to SPM2 and SPM8 should be defined in the
% sub-structure cfg.spm and can include:
% cfg.spm.regtype = 'subj', 'rigid'
% cfg.spm.smosrc = scalar value
% cfg.spm.smoref = scalar value
%
% Additional options pertaining to SPM12 should be defined in the
% sub-structure cfg.spm and can include:
% cfg.spm.sep = optimisation sampling steps (mm), default: [4 2]
% cfg.spm.params = starting estimates (6 elements), default: [0 0 0 0 0 0]
% cfg.spm.cost_fun = cost function string:
% 'mi' - Mutual Information (default)
% 'nmi' - Normalised Mutual Information
% 'ecc' - Entropy Correlation Coefficient
% 'ncc' - Normalised Cross Correlation
% cfg.spm.tol = tolerences for accuracy of each param, default: [0.02 0.02 0.02 0.001 0.001 0.001]
% cfg.spm.fwhm = smoothing to apply to 256x256 joint histogram, default: [7 7]
%
% When cfg.method is 'fsl', a third input argument is required. The input volume is
% coregistered to this target volume, using FSL-flirt. Additional options pertaining
% to the FSL method should be defined in the sub-structure cfg.fsl and can include:
% cfg.fsl.path = string, specifying the path to fsl
% cfg.fsl.costfun = string, specifying the cost-function used for
% coregistration
% cfg.fsl.interpmethod = string, specifying the interpolation method, can be
% 'trilinear', 'nearestneighbour', or 'sinc'
% cfg.fsl.dof = scalar, specifying the number of parameters for the
% affine transformation. 6 (rigid body), 7 (global
% rescale), 9 (traditional) or 12.
% cfg.fsl.reslice = string, specifying whether the output image will be
% resliced conform the target image (default = 'yes')
%
% To facilitate data-handling and distributed computing you can use
% cfg.inputfile = ...
% cfg.outputfile = ...
% If you specify one of these (or both) the input data will be read from a
% *.mat file on disk and/or the output data will be written to a *.mat
% file. These mat files should contain only a single variable,
% corresponding with the input/output structure.
%
% See also FT_READ_MRI, FT_VOLUMERESLICE, FT_INTERACTIVEREALIGN, FT_ELECTRODEREALIGN,
% FT_DETERMINE_COORDSYS, SPM_AFFREG, SPM_NORMALISE, SPM_COREG
% Undocumented options:
%
% cfg.weights = vector of weights that is used to weight the individual headshape
% points in the icp algorithm. Used optionally in cfg.method = 'headshape'. If not
% specified, weights are put on points with z-coordinate<0 (assuming those to be eye
% rims and nose ridges, i.e. important points.
% Copyright (C) 2006-2023, Robert Oostenveld, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation, either version 3 of the License, or (at your
% option) any later version.
%
% FieldTrip is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$
% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin = nargin;
ft_nargout = nargout;
% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar mri
ft_preamble provenance mri
% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
return
end
% the data can be passed as input argument or can be read from disk
hastarget = exist('target', 'var');
% check if the input data is valid for this function
mri = ft_checkdata(mri, 'datatype', 'volume', 'feedback', 'yes');
% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'renamedval', {'method', 'realignfiducial', 'fiducial'});
cfg = ft_checkconfig(cfg, 'renamed', {'landmark', 'fiducial'}); % cfg.landmark -> cfg.fiducial
% mni/spm/tal are to be interpreted as acpc with native scaling, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3304
cfg = ft_checkconfig(cfg, 'renamedval', {'coordsys', 'mni', 'acpc'});
cfg = ft_checkconfig(cfg, 'renamedval', {'coordsys', 'spm', 'acpc'});
cfg = ft_checkconfig(cfg, 'renamedval', {'coordsys', 'tal', 'acpc'});
% see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2837
cfg = ft_checkconfig(cfg, 'renamed', {'viewdim', 'axisratio'});
% set the defaults
cfg.coordsys = ft_getopt(cfg, 'coordsys', []);
cfg.method = ft_getopt(cfg, 'method', []); % the default is set below
cfg.flip = ft_getopt(cfg, 'flip', []); % the default is set below
cfg.fiducial = ft_getopt(cfg, 'fiducial', []);
cfg.parameter = ft_getopt(cfg, 'parameter', 'anatomy');
cfg.clim = ft_getopt(cfg, 'clim', []);
cfg.viewmode = ft_getopt(cfg, 'viewmode', 'ortho'); % for method=interactive
cfg.snapshot = ft_getopt(cfg, 'snapshot', false);
cfg.snapshotfile = ft_getopt(cfg, 'snapshotfile', fullfile(pwd, 'ft_volumerealign_snapshot'));
cfg.spmversion = ft_getopt(cfg, 'spmversion', 'spm12');
cfg.voxelratio = ft_getopt(cfg, 'voxelratio', 'data'); % display size of the voxel, 'data' or 'square'
cfg.axisratio = ft_getopt(cfg, 'axisratio', 'data'); % size of the axes of the three orthoplots, 'square', 'voxel', or 'data'
cfg.viewresult = ft_getopt(cfg, 'viewresult', 'no');
viewresult = istrue(cfg.viewresult);
if isempty(cfg.method)
if isempty(cfg.fiducial)
% fiducials have not yet been specified
cfg.method = 'interactive';
else
% fiducials have already been specified
cfg.method = 'fiducial';
end
end
if isempty(cfg.flip)
if strcmp(cfg.method, 'interactive')
cfg.flip = 'yes';
else
cfg.flip = 'no';
end
end
if isempty(cfg.coordsys)
if isstruct(cfg.fiducial) && all(ismember(fieldnames(cfg.fiducial), {'lpa', 'rpa', 'nas', 'zpoint'}))
cfg.coordsys = 'ctf';
elseif isstruct(cfg.fiducial) && all(ismember(fieldnames(cfg.fiducial), {'ac', 'pc', 'xzpoint', 'right'}))
cfg.coordsys = 'acpc';
elseif isstruct(cfg.fiducial) && all(ismember(fieldnames(cfg.fiducial), {'bregma', 'lambda', 'yzpoint'}))
cfg.coordsys = 'paxinos';
elseif strcmp(cfg.method, 'interactive')
cfg.coordsys = 'ctf';
end
ft_warning('defaulting to "%s" coordinate system', cfg.coordsys);
end
% these two have to be simultaneously true for a snapshot to be taken
dosnapshot = istrue(cfg.snapshot);
if dosnapshot
% create an empty array of handles
snap = [];
end
% select the parameter that should be displayed
cfg.parameter = parameterselection(cfg.parameter, mri);
if iscell(cfg.parameter) && ~isempty(cfg.parameter)
cfg.parameter = cfg.parameter{1};
elseif iscell(cfg.parameter) && isempty(cfg.parameter)
% cfg.parameter has been cleared by parameterselection due to a
% dimensionality mismatch. Most probable cause here is the fact that a 4D
% volume (e.g. DTI data) is in the input. This needs to be patched in a
% more structural way at some point, but for the time being we'll use a
% workaround here.
% assume anatomy to be the parameter of interest
siz = size(mri.anatomy);
if all(siz(1:3)==mri.dim) && numel(siz)==4
% it's OK
cfg.parameter= 'anatomy';
else
ft_error('there''s an unexpected dimension mismatch');
end
end
if strcmp(cfg.flip, 'yes')
if strcmp(cfg.method, 'interactive')
% align the anatomical volume approximately to coordinate system, this puts it upright
origmethod = cfg.method;
tmpcfg = [];
tmpcfg.method = 'flip';
tmpcfg.trackcallinfo = 'no';
tmpcfg.showcallinfo = 'no';
mri = ft_volumereslice(tmpcfg, mri);
[cfg, mri] = rollback_provenance(cfg, mri);
cfg.method = origmethod;
else
ft_error('flipping is only supported for method="interactive", please use FT_VOLUMERESLICE');
end
end
% start with an empty transform and coordsys
transform = [];
coordsys = [];
if any(strcmp(cfg.method, {'fiducial', 'interactive'}))
switch cfg.coordsys
case {'ctf', '4d', 'bti', 'eeglab', 'neuromag', 'itab', 'yokogawa', 'asa'}
fidlabel = {'nas', 'lpa', 'rpa', 'zpoint'};
fidletter = {'n', 'l', 'r', 'z'};
fidexplanation1 = ' press n for nas, l for lpa, r for rpa';
fidexplanation2 = ' press z for an extra control point that should have a positive z-value';
case 'acpc'
fidlabel = {'ac', 'pc', 'xzpoint', 'right'};
fidletter = {'a', 'p', 'z', 'r'};
fidexplanation1 = ' press a for ac, p for pc, z for xzpoint';
fidexplanation2 = ' press r for an extra control point that should be on the right side';
case 'paxinos'
fidlabel = {'bregma', 'lambda', 'yzpoint'};
fidletter = {'b', 'l', 'z'};
fidexplanation1 = ' press b for bregma, l for lambda, z for yzpoint';
fidexplanation2 = '';
otherwise
ft_error('unknown coordinate system "%s"', cfg.coordsys);
end
for i=1:length(fidlabel)
if ~isfield(cfg.fiducial, fidlabel{i}) || isempty(cfg.fiducial.(fidlabel{i}))
cfg.fiducial.(fidlabel{i}) = [nan nan nan];
end
end
end % interactive or fiducial
switch cfg.method
case 'fiducial'
% the actual coordinate transformation will be done further down
case 'landmark'
% the actual coordinate transformation will be done further down
case 'interactive'
% this requires the user to click the anatomical landmarks or fiducials
switch cfg.viewmode
case 'ortho'
% start building the figure
h = figure;
set(h, 'visible', 'on');
% axes settings
if strcmp(cfg.axisratio, 'voxel')
% determine the number of voxels to be plotted along each axis
axlen1 = mri.dim(1);
axlen2 = mri.dim(2);
axlen3 = mri.dim(3);
elseif strcmp(cfg.axisratio, 'data')
% determine the length of the edges along each axis
[cp_voxel, cp_head] = cornerpoints(mri.dim, mri.transform);
axlen1 = norm(cp_head(2,:)-cp_head(1,:));
axlen2 = norm(cp_head(4,:)-cp_head(1,:));
axlen3 = norm(cp_head(5,:)-cp_head(1,:));
elseif strcmp(cfg.axisratio, 'square')
% the length of the axes should be equal
axlen1 = 1;
axlen2 = 1;
axlen3 = 1;
end
% this is the size reserved for subplot h1, h2 and h3
h1size(1) = 0.82*axlen1/(axlen1 + axlen2);
h1size(2) = 0.82*axlen3/(axlen2 + axlen3);
h2size(1) = 0.82*axlen2/(axlen1 + axlen2);
h2size(2) = 0.82*axlen3/(axlen2 + axlen3);
h3size(1) = 0.82*axlen1/(axlen1 + axlen2);
h3size(2) = 0.82*axlen2/(axlen2 + axlen3);
if strcmp(cfg.voxelratio, 'square')
voxlen1 = 1;
voxlen2 = 1;
voxlen3 = 1;
elseif strcmp(cfg.voxelratio, 'data')
% the size of the voxel is scaled with the data
[cp_voxel, cp_head] = cornerpoints(mri.dim, mri.transform);
voxlen1 = norm(cp_head(2,:)-cp_head(1,:))/norm(cp_voxel(2,:)-cp_voxel(1,:));
voxlen2 = norm(cp_head(4,:)-cp_head(1,:))/norm(cp_voxel(4,:)-cp_voxel(1,:));
voxlen3 = norm(cp_head(5,:)-cp_head(1,:))/norm(cp_voxel(5,:)-cp_voxel(1,:));
end
%% the figure is interactive, add callbacks
set(h, 'windowbuttondownfcn', @cb_buttonpress);
set(h, 'windowbuttonupfcn', @cb_buttonrelease);
set(h, 'windowkeypressfcn', @cb_keyboard);
set(h, 'CloseRequestFcn', @cb_quit);
% axis handles will hold the anatomical functional if present, along with labels etc.
h1 = axes('position', [0.06 0.06+0.06+h3size(2) h1size(1) h1size(2)]);
h2 = axes('position', [0.06+0.06+h1size(1) 0.06+0.06+h3size(2) h2size(1) h2size(2)]);
h3 = axes('position', [0.06 0.06 h3size(1) h3size(2)]);
set(h1, 'Tag', 'ik', 'Visible', 'off', 'XAxisLocation', 'top');
set(h2, 'Tag', 'jk', 'Visible', 'off', 'YAxisLocation', 'right'); % after rotating in ft_plot_ortho this becomes top
set(h3, 'Tag', 'ij', 'Visible', 'off');
set(h1, 'DataAspectRatio', 1./[voxlen1 voxlen2 voxlen3]);
set(h2, 'DataAspectRatio', 1./[voxlen1 voxlen2 voxlen3]);
set(h3, 'DataAspectRatio', 1./[voxlen1 voxlen2 voxlen3]);
xc = round(mri.dim(1)/2); % start with center view
yc = round(mri.dim(2)/2);
zc = round(mri.dim(3)/2);
% enhance the contrast of the volumetric data, see also FT_DEFACEVOLUME
dat = double(mri.(cfg.parameter));
dmin = min(dat(:));
dmax = max(dat(:));
dat = (dat-dmin)./(dmax-dmin);
if isempty(cfg.clim)
cfg.clim = [0 1];
else
% apply the same scaling as to the data
cfg.clim = (cfg.clim-dmin)./(dmax-dmin);
end
if isfield(cfg, 'pnt')
pnt = cfg.pnt;
else
pnt = zeros(0,3);
end
markerpos = zeros(0,3);
markerlabel = {};
markercolor = {};
% determine the apprioriate [left bottom width height] position of the intensity range sliders
posbase = [];
posbase(1) = h1size(1) + h2size(1)/2 + 0.06*2; % horizontal center of the second plot
posbase(2) = h3size(2)/2 + 0.06; % vertical center of the third plot
posbase(3) = 0.01; % width of the sliders is not so important, if it falls below a certain value, it's a vertical slider, otherwise a horizontal one
posbase(4) = h3size(2)/3 + 0.06; % one-third of the height of the third plot
%
posh45text = [posbase(1)-posbase(3)*5 posbase(2)-.1 posbase(3)*10 posbase(4)+0.07];
posh4text = [posbase(1)-.04-posbase(3)*2 posbase(2)-.1 posbase(3)*5 posbase(4)+0.035];
posh5text = [posbase(1)+.04-posbase(3)*2 posbase(2)-.1 posbase(3)*5 posbase(4)+0.035];
posh4slid = [posbase(1)-.04 posbase(2)-.1 posbase(3) posbase(4)];
posh5slid = [posbase(1)+.04 posbase(2)-.1 posbase(3) posbase(4)];
% intensity range sliders
uicontrol('Style', 'text',...
'String', 'Intensity', ...
'Units', 'normalized', ...
'Position', posh45text, ...
'HandleVisibility', 'on');
h4text = uicontrol('Style', 'text',...
'String', 'Min', ...
'Units', 'normalized', ...
'Position', posh4text, ... % text is centered, so height adjust vertical position
'HandleVisibility', 'on');
h4 = uicontrol('Style', 'slider', ...
'Parent', h, ...
'Min', 0, 'Max', 1, ...
'Value', cfg.clim(1), ...
'Units', 'normalized', ...
'Position', posh4slid, ...
'Callback', @cb_minslider);
h5text = uicontrol('Style', 'text',...
'String', 'Max',...
'Units', 'normalized', ...
'Position',posh5text, ... % text is centered, so height adjust vertical position
'HandleVisibility', 'on');
h5 = uicontrol('Style', 'slider', ...
'Parent', h, ...
'Min', 0, 'Max', 1, ...
'Value', cfg.clim(2), ...
'Units', 'normalized', ...
'Position', posh5slid, ...
'Callback', @cb_maxslider);
% create structure to be passed to gui
opt = [];
opt.viewmode = 'ortho';
opt.viewresult = false; % flag to use for certain keyboard/redraw calls
opt.init = true;
opt.quit = false;
opt.twovol = false; % flag to use for certain options of viewresult
opt.dim = mri.dim;
opt.ijk = [xc yc zc];
opt.h1size = h1size;
opt.h2size = h2size;
opt.h3size = h3size;
opt.h4 = h4;
opt.h5 = h5;
opt.handlesaxes = [h1 h2 h3];
opt.handlesfigure = h;
opt.ana = dat;
opt.update = [1 1 1];
opt.tag = 'ik';
opt.mri = mri;
opt.showcrosshair = true;
opt.showmarkers = false;
opt.markers = {markerpos markerlabel markercolor};
opt.clim = cfg.clim;
opt.fiducial = cfg.fiducial;
opt.fidlabel = fidlabel;
opt.fidletter = fidletter;
opt.fidexplanation1 = fidexplanation1;
opt.fidexplanation2 = fidexplanation2;
opt.pnt = pnt;
if isfield(mri, 'unit') && ~strcmp(mri.unit, 'unknown')
opt.unit = mri.unit; % this is shown in the feedback on screen
else
opt.unit = ''; % this is not shown
end
setappdata(h, 'opt', opt);
cb_redraw(h);
cb_help(h);
case 'surface'
% make a mesh from the scalp surface
cfg.headshape = ft_getopt(cfg, 'headshape');
cfg.headshape.scalpsmooth = ft_getopt(cfg.headshape, 'scalpsmooth', 2, 1); % empty is OK
cfg.headshape.scalpthreshold = ft_getopt(cfg.headshape, 'scalpthreshold', 0.1);
if ~isfield(mri, 'scalp') || ~islogical(mri.scalp)
% extract the scalp surface from the anatomical image
tmpcfg = [];
tmpcfg.output = 'scalp';
tmpcfg.scalpsmooth = cfg.headshape.scalpsmooth;
tmpcfg.scalpthreshold = cfg.headshape.scalpthreshold;
if isfield(cfg, 'template')
tmpcfg.template = cfg.template;
end
seg = ft_volumesegment(tmpcfg, mri);
else
% use the scalp segmentation that is provided
seg = mri;
end
tmpcfg = [];
tmpcfg.tissue = 'scalp';
tmpcfg.method = 'isosurface';
tmpcfg.spmversion = cfg.spmversion;
tmpcfg.numvertices = inf;
scalp = ft_prepare_mesh(tmpcfg, seg);
scalp = ft_convert_units(scalp, 'mm');
% start building the figure
h = figure;
set(h, 'color', [1 1 1]);
set(h, 'visible', 'on');
% add callbacks
set(h, 'windowkeypressfcn', @cb_keyboard);
set(h, 'CloseRequestFcn', @cb_quit);
% create figure handles
h1 = axes;
% create structure to be passed to gui
opt = [];
opt.viewmode = 'surface';
opt.viewresult = false; % flag to use for certain keyboard/redraw calls
opt.init = true;
opt.quit = false;
opt.handlesfigure = h;
opt.handlesaxes = h1;
opt.handlesfigure = h;
opt.handlesmarker = [];
opt.camlighthandle = [];
opt.scalp = scalp;
opt.showmarkers = false;
opt.mri = mri;
opt.fiducial = cfg.fiducial;
opt.fidlabel = fidlabel;
opt.fidletter = fidletter;
opt.fidexplanation1 = fidexplanation1;
opt.fidexplanation2 = fidexplanation2;
if isfield(scalp, 'unit') && ~strcmp(scalp.unit, 'unknown')
opt.unit = scalp.unit; % this is shown in the feedback on screen
else
opt.unit = ''; % this is not shown
end
setappdata(h, 'opt', opt);
cb_redraw(h);
otherwise
ft_error('unsupported viewmode "%s"', cfg.viewmode);
end % switch viewmode
while(opt.quit==0)
uiwait(h);
opt = getappdata(h, 'opt');
end
delete(h);
% store the interactively determined fiducials in the configuration
% the actual coordinate transformation will be done further down
cfg.fiducial = opt.fiducial;
case 'headshape'
if ischar(cfg.headshape)
% old-style specification, convert cfg into new representation
cfg.headshape = struct('headshape', cfg.headshape);
if isfield(cfg, 'scalpsmooth')
cfg.headshape.scalpsmooth = cfg.scalpsmooth;
cfg = rmfield(cfg, 'scalpsmooth');
end
if isfield(cfg, 'scalpthreshold')
cfg.headshape.scalpthreshold = cfg.scalpthreshold;
cfg = rmfield(cfg, 'scalpthreshold');
end
elseif isstruct(cfg.headshape) && isfield(cfg.headshape, 'pos')
% old-style specification, convert into new representation
cfg.headshape = struct('headshape', cfg.headshape);
if isfield(cfg, 'scalpsmooth')
cfg.headshape.scalpsmooth = cfg.scalpsmooth;
cfg = rmfield(cfg, 'scalpsmooth');
end
if isfield(cfg, 'scalpthreshold')
cfg.headshape.scalpthreshold = cfg.scalpthreshold;
cfg = rmfield(cfg, 'scalpthreshold');
end
elseif isstruct(cfg.headshape)
% new-style specification, do nothing
else
ft_error('incorrect specification of cfg.headshape');
end
if ischar(cfg.headshape.headshape)
shape = ft_read_headshape(cfg.headshape.headshape);
else
shape = cfg.headshape.headshape;
end
shape = ft_convert_units(shape, mri.unit); % make the units of the headshape consistent with the MRI
cfg.headshape.interactive = ft_getopt(cfg.headshape, 'interactive', true);
cfg.headshape.icp = ft_getopt(cfg.headshape, 'icp', true);
cfg.headshape.scalpsmooth = ft_getopt(cfg.headshape, 'scalpsmooth', 2, 1); % empty is OK
cfg.headshape.scalpthreshold = ft_getopt(cfg.headshape, 'scalpthreshold', 0.1);
dointeractive = istrue(cfg.headshape.interactive);
doicp = istrue(cfg.headshape.icp);
if ~isfield(mri, 'scalp') || ~islogical(mri.scalp)
% extract the scalp surface from the anatomical image
tmpcfg = [];
tmpcfg.output = 'scalp';
tmpcfg.spmversion = cfg.spmversion;
tmpcfg.scalpsmooth = cfg.headshape.scalpsmooth;
tmpcfg.scalpthreshold = cfg.headshape.scalpthreshold;
if isfield(cfg, 'template')
tmpcfg.template = cfg.template;
end
seg = ft_volumesegment(tmpcfg, mri);
else
% use the scalp segmentation that is provided
seg = mri;
end
tmpcfg = [];
tmpcfg.tissue = 'scalp';
tmpcfg.method = 'projectmesh'; %'isosurface';
tmpcfg.spmversion = cfg.spmversion;
tmpcfg.numvertices = 20000;
scalp = ft_prepare_mesh(tmpcfg, seg);
if dointeractive
fprintf('doing interactive realignment with headshape\n');
tmpcfg = [];
tmpcfg.template.headshape = shape; % this is the Polhemus recorded headshape
tmpcfg.template.headshapestyle = 'vertex';
tmpcfg.individual.headshape = scalp; % this is the headshape extracted from the anatomical MRI
tmpcfg.individual.headshapestyle = 'surface';
tmpcfg = ft_interactiverealign(tmpcfg);
M = tmpcfg.m;
cfg.transform_interactive = M;
% update the relevant geometrical info
scalp = ft_transform_geometry(M, scalp);
end % dointeractive
% always perform an icp-step, because this will give an estimate of the
% initial distance of the corresponding points. depending on the value
% for doicp, deal with the output differently
if doicp
numiter = 50;
else
numiter = 1;
end
if ~isfield(cfg, 'weights')
w = ones(size(shape.pos,1),1);
else
w = cfg.weights(:);
if numel(w)~=size(shape.pos,1)
ft_error('number of weights should be equal to the number of points in the headshape');
end
end
% the icp function wants this as a function handle.
weights = @(x)assignweights(x,w);
ft_hastoolbox('fileexchange',1);
% construct the coregistration matrix
nrm = surface_normals(scalp.pos, scalp.tri, 'vertex');
[R, t, err, dummy, info] = icp(scalp.pos', shape.pos', numiter, 'Minimize', 'plane', 'Normals', nrm', 'Weight', weights, 'Extrapolation', true, 'WorstRejection', 0.05);
if doicp
fprintf('doing iterative closest points realignment with headshape\n');
% create the additional transformation matrix and compute the
% distance between the corresponding points, both prior and after icp
% this one transforms from scalp 'headspace' to shape 'headspace'
M2 = inv([R t;0 0 0 1]);
% warp the extracted scalp points to the new positions
scalp.pos = ft_warp_apply(M2, scalp.pos);
target = scalp;
target.pos = target.pos;
target.inside = (1:size(target.pos,1))';
functional = rmfield(shape, 'pos');
functional.distance = info.distanceout(:);
functional.pos = info.qout';
tmpcfg = [];
tmpcfg.parameter = 'distance';
tmpcfg.interpmethod = 'sphere_avg';
tmpcfg.sphereradius = 10;
tmpcfg.feedback = 'none';
smoothdist = ft_sourceinterpolate(tmpcfg, functional, target);
scalp.distance = smoothdist.distance(:);
functional.distance = info.distancein(:);
smoothdist = ft_sourceinterpolate(tmpcfg, functional, target);
scalp.distancein = smoothdist.distance(:);
cfg.icpinfo = info;
cfg.transform_icp = M2;
else
% compute the distance between the corresponding points, prior to icp:
% this corresponds to the final result after interactive only
M2 = eye(4); % this is needed later on
target = scalp;
target.pos = target.pos;
target.inside = (1:size(target.pos,1))';
functional = rmfield(shape, 'pos');
functional.pow = info.distancein(:);
functional.pos = info.qout';
tmpcfg = [];
tmpcfg.parameter = 'pow';
tmpcfg.interpmethod = 'sphere_avg';
tmpcfg.sphereradius = 10;
smoothdist = ft_sourceinterpolate(tmpcfg, functional, target);
scalp.distance = smoothdist.pow(:);
end % doicp
% create headshape structure for mri-based surface point cloud
if isfield(mri, 'coordsys')
scalp.coordsys = mri.coordsys;
% coordsys is the same as input mri
coordsys = mri.coordsys;
else
coordsys = 'unknown';
end
% update the cfg
cfg.headshape.headshape = shape;
cfg.headshape.headshapemri = scalp;
if doicp && dointeractive
transform = M2*M;
elseif doicp
transform = M2;
elseif dointeractive
transform = M;
end
case 'fsl'
if ~isfield(cfg, 'fsl'), cfg.fsl = []; end
cfg.fsl.path = ft_getopt(cfg.fsl, 'path', '');
cfg.fsl.costfun = ft_getopt(cfg.fsl, 'costfun', 'corratio');
cfg.fsl.interpmethod = ft_getopt(cfg.fsl, 'interpmethod', 'trilinear');
cfg.fsl.dof = ft_getopt(cfg.fsl, 'dof', 6);
cfg.fsl.reslice = ft_getopt(cfg.fsl, 'reslice', 'yes');
cfg.fsl.searchrange = ft_getopt(cfg.fsl, 'searchrange', [-180 180]);
% write the input and target to a temporary file
% and create some additional temporary file names to contain the output
filename_mri = tempname;
filename_target = tempname;
filename_output = tempname;
filename_mat = tempname;
tmpcfg = [];
tmpcfg.parameter = 'anatomy';
tmpcfg.filename = filename_mri;
tmpcfg.filetype = 'nifti';
fprintf('writing the input volume to a temporary file: %s\n', [filename_mri, '.nii']);
ft_volumewrite(tmpcfg, mri);
tmpcfg.filename = filename_target;
fprintf('writing the target volume to a temporary file: %s\n', [filename_target, '.nii']);
ft_volumewrite(tmpcfg, target);
% create the command to call flirt
fprintf('using flirt for the coregistration\n');
r1 = num2str(cfg.fsl.searchrange(1));
r2 = num2str(cfg.fsl.searchrange(2));
str = sprintf('%s/flirt -in %s -ref %s -out %s -omat %s -bins 256 -cost %s -searchrx %s %s -searchry %s %s -searchrz %s %s -dof %s -interp %s',...
cfg.fsl.path, filename_mri, filename_target, filename_output, filename_mat, cfg.fsl.costfun, r1, r2, r1, r2, r1, r2, num2str(cfg.fsl.dof), cfg.fsl.interpmethod);
if isempty(cfg.fsl.path), str = str(2:end); end % remove the first filesep, assume path to flirt to be known
% system call
system(str);
% process the output
if ~istrue(cfg.fsl.reslice)
% get the transformation that corresponds to the coregistration and
% reconstruct the mapping from the target's world coordinate system
% to the input's voxel coordinate system
fid = fopen(filename_mat);
flirtmat = textscan(fid, '%f');
fclose(fid);
% this contains the coregistration information in FSL convention
flirtmat = reshape(flirtmat{1},4,4)';
% The following chunck of code is from Ged Ridgway's
% flirtmat2worldmat code
% src = inv(flirtmat) * trg
% srcvox = src.mat \ inv(flirtmat) * trg.mat * trgvox
% BUT, flirt doesn't use src.mat, only absolute values of the
% scaling elements from it,
% AND, if images are not radiological, the x-axis is flipped, see:
% https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind0810&L=FSL&P=185638
% https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind0903&L=FSL&P=R93775
scl_target = diag([sqrt(sum(target.transform(1:3,1:3).^2)) 1]);
if det(target.transform(1:3,1:3)) > 0
% neurological, x-axis is flipped, such that [3 2 1 0] and [0 1 2 3]
% have the same *scaled* coordinates:
xflip = diag([-1 1 1 1]);
xflip(1, 4) = target.dim(1)-1; % reflect about centre
scl_target = scl_target * xflip;
end
scl_mri = diag([sqrt(sum(mri.transform(1:3,1:3).^2)) 1]);
if det(mri.transform(1:3,1:3)) > 0
% neurological, x-axis is flipped, such that [3 2 1 0] and [0 1 2 3]
% have the same *scaled* coordinates:
xflip = diag([-1 1 1 1]);
xflip(1, 4) = mri.dim(1)-1; % reflect about centre
scl_mri = scl_mri * xflip;
end
% AND, Flirt's voxels are zero-based, while SPM's are one-based...
addone = eye(4);
addone(:, 4) = 1;
fslvoxmat = inv(scl_mri) * inv(flirtmat) * scl_target;
spmvoxmat = addone * (fslvoxmat / addone);
target2mri = mri.transform * (spmvoxmat / target.transform);
transform = inv(target2mri);
if isfield(target, 'coordsys')
coordsys = target.coordsys;
else
coordsys = 'unknown';
end
else
% get the updated anatomy
mrinew = ft_read_mri([filename_output, '.nii.gz']);
mri.anatomy = mrinew.anatomy;
mri.transform = mrinew.transform;
mri.dim = mrinew.dim;
transform = eye(4);
if isfield(target, 'coordsys')
coordsys = target.coordsys;
else
coordsys = 'unknown';
end
end
delete([filename_mri, '.nii']);
delete([filename_target, '.nii']);
delete([filename_output, '.nii.gz']);
delete(filename_mat);
case 'spm'
% check that the preferred SPM version is on the path
ft_hastoolbox(cfg.spmversion, 1);
if strcmpi(cfg.spmversion, 'spm2') || strcmpi(cfg.spmversion, 'spm8')
if ~isfield(cfg, 'spm'), cfg.spm = []; end
cfg.spm.regtype = ft_getopt(cfg.spm, 'regtype', 'subj');
cfg.spm.smosrc = ft_getopt(cfg.spm, 'smosrc', 2);
cfg.spm.smoref = ft_getopt(cfg.spm, 'smoref', 2);
if ~isfield(mri, 'coordsys')
mri = ft_determine_coordsys(mri);
else
fprintf('Input volume has coordinate system ''%s''\n', mri.coordsys);
end
if ~isfield(target, 'coordsys')
target = ft_determine_coordsys(target);
else
fprintf('Target volume has coordinate system ''%s''\n', target.coordsys);
end
if strcmp(mri.coordsys, target.coordsys)
% this should hopefully work
else
% only works when it is possible to approximately align the input to the target coordsys
if strcmp(target.coordsys, 'acpc')
mri = ft_convert_coordsys(mri, 'acpc');
else
ft_error('The coordinate systems of the input and target volumes are different, coregistration is not possible');
end
end
% flip and permute the 3D volume itself, so that the voxel and
% headcoordinates approximately correspond
[tmp, pvec_mri, flip_mri, T] = align_ijk2xyz(mri);
[target] = align_ijk2xyz(target);
tname1 = [tempname, '.img'];
tname2 = [tempname, '.img'];
V1 = ft_write_mri(tname1, mri.anatomy, 'transform', mri.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
V2 = ft_write_mri(tname2, target.anatomy, 'transform', target.transform, 'spmversion', spm('ver'), 'dataformat', 'nifti_spm');
flags = cfg.spm;
flags.nits = 0; %set number of non-linear iterations to zero
params = spm_normalise(V2,V1, [], [], [],flags);
%mri.transform = (target.transform/params.Affine)/T;
transform = (target.transform/params.Affine)/T/mri.transform;
% transform = eye(4);
elseif strcmpi(cfg.spmversion, 'spm12')
if ~isfield(cfg, 'spm'), cfg.spm = []; end
tname1 = [tempname, '.nii'];
tname2 = [tempname, '.nii'];