-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathprepareForAbaqus.m
153 lines (143 loc) · 7.12 KB
/
prepareForAbaqus.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
% This function creates rnge_elec and rnge_gel files for post scanIP
% processing. Apply this script in batch using batch_prepare_forAbaqus.m
% where each subject's data is stored in separate folders. This script must
% be applied after running script batch_roast_forAbaqus.m
function prepareForAbaqus(Image,elecgellabelsMat,optMat)
[dirname,~,~]=fileparts(Image);
%%
vol=load_untouch_nii(Image);
data=vol.img;
data_backup=data;
for v=1:6
data(data~=v)=0;
data(data==v)=1;
vol.img=data;
%save_untouch_nii(vol,sprintf('class_%d_%s',v,n));
save_untouch_nii(vol,fullfile(dirname,sprintf('class_%d',v)));
data=data_backup;
end
%%
load(elecgellabelsMat)
load(optMat)
%hdrInfo= electrodePlacement(Image,Image,[],elecConfig,opt,[]);
%%
template = load_untouch_nii(fullfile(dirname,'class_5.nii'));
scalp = template.img;
nasion=landmarks_original(1,:);inion=landmarks_original(2,:); right=landmarks_original(3,:);
left=landmarks_original(4,:);front_neck=landmarks_original(5,:);back_neck=landmarks_original(6,:);
e1 = right-left; e1 = e1/norm(e1);
e2 = nasion-inion; e2 = e2/norm(e2);
e3 = nasion-front_neck; e3 = e3/norm(e3); % detect the orientation of the head based on the anatomical landmarks
[~,perm1] = max(abs(e1)); [~,perm2] = max(abs(e2)); [~,perm3] = max(abs(e3));
isFlip = [sign(e1(perm1)) sign(e2(perm2)) sign(e3(perm3))]; % detect if the head is flipped or not in each direction compared to RAS system
perm = [perm1,perm2,perm3]; % permutation order into RAS
[~,iperm] = sort(perm); % inverse permutation order back to original orientation of the head
scalp = permute(scalp,perm); % permute the head into RAS
[Nx, Ny, Nz]=size(scalp); % size of head in RAS orientation
%% for elec rnge mat
elec_C_final = cell(1,size(electrode_coord,1));
rnge = cell(1,size(electrode_coord,1));
isOut = zeros(length(rnge),1);
for i = 1:size(electrode_coord,1)
elec_C_final{i} = zeros(size(elec_C{i},1),3);
rng = zeros(2,3);
for j=1:size(elec_C{i},1)
if elec_C{i}(j,1)>0 && elec_C{i}(j,1)<=Nx && elec_C{i}(j,2)>0 && elec_C{i}(j,2)<=Ny && elec_C{i}(j,3)>0 && elec_C{i}(j,3)<=Nz
elec_C_final{i}(j,1) = elec_C{i}(j,iperm(1));
elec_C_final{i}(j,2) = elec_C{i}(j,iperm(2));
elec_C_final{i}(j,3) = elec_C{i}(j,iperm(3));
end % permute back electrode coordinates to match the original orientation of the head
end
if all(sum(elec_C_final{i},2))
rng(1,:) = max(elec_C_final{i});
rng(2,:) = min(elec_C_final{i});
rnge{i} = rng;
else if any(sum(elec_C_final{i},2)) % in case of some electrode voxels go out of image boundary then elec_C_final has 0
elec_C_temp = elec_C_final{i}(sum(elec_C_final{i},2)>0,:);
rng(1,:) = max(elec_C_temp);
rng(2,:) = min(elec_C_temp);
rnge{i} = rng;
else rnge{i} = []; % this is the case that the electrode completely goes out of image boundary, thus this electrode actually does not exist
warning(['Electrode #' num2str(i) ' goes out of image boundary!']);
isOut(i) = 1;
end
end
end
if any(isOut)
warning('Some of the electrodes go out of image boundary, the program can continue, but it is highly recommended that you expand the image by adding empty slices on the boundaries.');
end
for i=1:length(rnge)
if ~isempty(rnge{i})
rnge{i} = rnge{i}.*repmat(template.hdr.dime.pixdim(2:4),2,1);
end
end % use NIFTI header info to convert range info into world coordinates for subsequent electrode labeling ANDY 2014-08-12
% It's in fact a hack, i.e., only applies the scaling to the range
% information, to match the pseudo-world space of the mesh coordinates
% (generated by ScanIP)
% Get the range of coordinates for each electrode for subsequent electrode labelling
% (so that we can use a script to automatically specify anode and cathode
% when solving current flow in Abaqus, and also calculate EXACT energized area for
% each electrode)
save(fullfile(dirname,'rnge_elec.mat'),'rnge');
%% for gel rnge mat
gel_C_final = cell(1,size(electrode_coord,1));
for i = 1:size(electrode_coord,1)
gel_C_final{i} = zeros(size(gel_C{i},1),3);
for j=1:size(gel_C{i},1)
if gel_C{i}(j,1)>0 && gel_C{i}(j,1)<=Nx && gel_C{i}(j,2)>0 && gel_C{i}(j,2)<=Ny && gel_C{i}(j,3)>0 && gel_C{i}(j,3)<=Nz
gel_C_final{i}(j,1) = gel_C{i}(j,iperm(1));
gel_C_final{i}(j,2) = gel_C{i}(j,iperm(2));
gel_C_final{i}(j,3) = gel_C{i}(j,iperm(3));
end % permute back gel coordinates to match the original orientation of the head
end
end % Get ready for the generation of the range of coordinates for gel in the following codes
% ANDY 2014-03-04
disp('constructing electrode and gel volume to be exported...')
for i = 1:size(electrode_coord,1)
for j=1:size(gel_C{i},1)
if gel_C{i}(j,1)>0 && gel_C{i}(j,1)<=Nx && gel_C{i}(j,2)>0 && gel_C{i}(j,2)<=Ny && gel_C{i}(j,3)>0 && gel_C{i}(j,3)<=Nz
volume_gel(gel_C{i}(j,1), gel_C{i}(j,2), gel_C{i}(j,3)) = 1;
end
end
for j=1:size(elec_C{i},1)
if elec_C{i}(j,1)>0 && elec_C{i}(j,1)<=Nx && elec_C{i}(j,2)>0 && elec_C{i}(j,2)<=Ny && elec_C{i}(j,3)>0 && elec_C{i}(j,3)<=Nz
volume_elec(elec_C{i}(j,1), elec_C{i}(j,2), elec_C{i}(j,3)) = 1;
end
end
end
rnge = cell(1,size(electrode_coord,1));
for i = 1:size(electrode_coord,1)
isGel = zeros(size(gel_C_final{i},1),1);
for j=1:length(isGel)
if all(gel_C_final{i}(j,:)) && volume_gel(gel_C_final{i}(j,1),gel_C_final{i}(j,2),gel_C_final{i}(j,3))==1
isGel(j)=1;
end
end
gel_C_final{i} = gel_C_final{i}(isGel==1,:); % remove those gel coordinates that go into scalp/electrode/bone
rng = zeros(2,3);
if all(sum(gel_C_final{i},2))
rng(1,:) = max(gel_C_final{i});
rng(2,:) = min(gel_C_final{i});
rnge{i} = rng;
else if any(sum(gel_C_final{i},2)) % in case of some gel voxels go out of image boundary then gel_C_final has 0
gel_C_temp = gel_C_final{i}(sum(gel_C_final{i},2)>0,:);
rng(1,:) = max(gel_C_temp);
rng(2,:) = min(gel_C_temp);
rnge{i} = rng;
else rnge{i} = []; % this is the case that the gel completely goes out of image boundary, thus this gel actually does not exist
warning(['Gel of electrode #' num2str(i) ' goes out of image boundary!']);
end
end
end
for i=1:length(rnge)
if ~isempty(rnge{i})
rnge{i} = rnge{i}.*repmat(template.hdr.dime.pixdim(2:4),2,1);
end
end % use NIFTI header info to convert range info into world coordinates for subsequent gel labeling ANDY 2014-08-12
% It's in fact a hack, i.e., only applies the scaling to the range
% information, to match the pseudo-world space of the mesh coordinates
% (generated by ScanIP)
% Get the range of coordinates for each gel for subsequent gel labelling
% (so that we can use a script to automatically calculate EXACT energized area for each electrode)
save(fullfile(dirname,'rnge_gel.mat'),'rnge');
end