-
Notifications
You must be signed in to change notification settings - Fork 24
/
prepareForGetDP.m
129 lines (105 loc) · 5.57 KB
/
prepareForGetDP.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
function prepareForGetDP(subj,node,elem,elecNeeded,uniTag)
% prepareForGetDP(subj,node,elem,elecNeeded,uniTag)
%
% Prepare to solve in getDP
%
% (c) Yu (Andy) Huang, Parra Lab at CCNY
% October 2017
[dirname,subjName] = fileparts(subj);
if isempty(dirname), dirname = pwd; end
% node = node + 0.5; already done right after mesh
% indNode_elecElm = elem(find(elem(:,5) == 8),1:4);
% X = zeros(size(indNode_elecElm,1),3);
% for e = 1:size(indNode_elecElm,1), X(e,:) = mean ( node(indNode_elecElm(e,:),1:3) ); end
% % figure; plot3(X(:,1),X(:,2),X(:,3),'r.');
%
% Xt = round(X);
% label_elec = volume_elecLabel(sub2ind(size(volume_elecLabel),Xt(:,1),Xt(:,2),Xt(:,3)));
% indBad = find(label_elec==0);
% indGood = find(label_elec>0);
% [~,indOnGood] = map2Points(X(indBad,:),X(indGood,:),'closest'); % using nearest neighbor to fix bad labeling
% label_elec(indBad) = label_elec(indGood(indOnGood));
%
% indNode_gelElm = elem(find(elem(:,5) == 7),1:4);
% X = zeros(size(indNode_gelElm,1),3);
% for e = 1:size(indNode_gelElm,1), X(e,:) = mean ( node(indNode_gelElm(e,:),1:3) ); end
% % figure; plot3(X(:,1),X(:,2),X(:,3),'r.');
%
% Xt = round(X);
% label_gel = volume_gelLabel(sub2ind(size(volume_gelLabel),Xt(:,1),Xt(:,2),Xt(:,3)));
% indBad = find(label_gel==0);
% indGood = find(label_gel>0);
% [~,indOnGood] = map2Points(X(indBad,:),X(indGood,:),'closest'); % using nearest neighbor to fix bad labeling
% label_gel(indBad) = label_gel(indGood(indOnGood));
%
% save([dirname filesep subjName '_' uniTag '_elecMeshLabels.mat'],'label_elec','label_gel');
numOfTissue = 6; % hard coded across ROAST.
numOfElec = length(elecNeeded);
element_elecNeeded = cell(numOfElec,1);
area_elecNeeded = zeros(numOfElec,1);
% resolution = mean(hdrInfo.pixdim);
% % mean() here to handle anisotropic rmesolution; ugly. Maybe just
% % resample MRI to isotropic in the very beginning?
% Not needed now when mesh coordinates are in pseudo-world space % ANDY 2019-03-13
warning('off','MATLAB:TriRep:PtsNotInTriWarnId');
for i=1:numOfElec
% if isempty(indNode_elecElm(label_elec==i,:))
% error(['Electrode ' elecNeeded{i} ' was not meshed properly. Reasons may be: 1) electrode size is too small so the mesher cannot capture it; 2) mesh resolution is not high enough. Consider using bigger electrodes or increasing the mesh resolution by specifying the mesh options.']);
% end
%
% [faces_elec,verts_elec] = freeBoundary(TriRep(indNode_elecElm(label_elec==i,:),node(:,1:3)));
% [faces_gel,verts_gel] = freeBoundary(TriRep(indNode_gelElm(label_gel==i,:),node(:,1:3)));
indNode_gelElm = elem(find(elem(:,5) == numOfTissue+i),1:4);
indNode_elecElm = elem(find(elem(:,5) == numOfTissue+numOfElec+i),1:4);
if isempty(indNode_gelElm)
error(['Gel under electrode ' elecNeeded{i} ' was not meshed properly. Reasons may be: 1) electrode size is too small so the mesher cannot capture it; 2) mesh resolution is not high enough. Consider using bigger electrodes or increasing the mesh resolution by specifying the mesh options.']);
end
if isempty(indNode_elecElm)
error(['Electrode ' elecNeeded{i} ' was not meshed properly. Reasons may be: 1) electrode size is too small so the mesher cannot capture it; 2) mesh resolution is not high enough. Consider using bigger electrodes or increasing the mesh resolution by specifying the mesh options.']);
end
[faces_gel,verts_gel] = freeBoundary(TriRep(indNode_gelElm,node(:,1:3)));
[faces_elec,verts_elec] = freeBoundary(TriRep(indNode_elecElm,node(:,1:3)));
[~,iE,iG] = intersect(verts_elec,verts_gel,'rows');
tempTag = ismember(faces_elec,iE);
% faces_overlap = faces_elec(sum(tempTag,2)==3,:);
faces_elecOuter = faces_elec(~(sum(tempTag,2)==3),:);
[~,Loc] = ismember(verts_elec,node(:,1:3),'rows');
element_elecNeeded{i} = Loc(faces_elecOuter);
% calculate the surface area
a = (verts_elec(faces_elecOuter(:, 2),:) - verts_elec(faces_elecOuter(:, 1),:)); %*resolution;
b = (verts_elec(faces_elecOuter(:, 3),:) - verts_elec(faces_elecOuter(:, 1),:)); %*resolution;
c = cross(a, b, 2);
area_elecNeeded(i) = sum(0.5*sqrt(sum(c.^2, 2)));
end
if ~exist([dirname filesep subjName '_' uniTag '_usedElecArea.mat'],'file')
save([dirname filesep subjName '_' uniTag '_usedElecArea.mat'],'area_elecNeeded');
end
if ~exist([dirname filesep subjName '_' uniTag '_ready.msh'],'file')
disp('setting up boundary conditions...');
fid_in = fopen([dirname filesep subjName '_' uniTag '.msh']);
fid_out = fopen([dirname filesep subjName '_' uniTag '_ready.msh'],'w');
numOfPart = length(unique(elem(:,5)));
while ~feof(fid_in)
s = fgetl(fid_in);
if strcmp(s,'$Elements')
fprintf(fid_out,'%s\n',s);
s = fgetl(fid_in);
numOfElem = str2num(s);
fprintf(fid_out,'%s\n',num2str(numOfElem+size(cell2mat(element_elecNeeded),1)));
elseif strcmp(s,'$EndElements')
ii = 0;
for j=1:numOfElec
for i=1:size(element_elecNeeded{j},1)
fprintf(fid_out,'%s \n',[num2str(numOfElem+i+ii) ' 2 2 ' num2str(numOfPart+j) ' ' num2str(numOfPart+j) ' ' num2str(element_elecNeeded{j}(i,1)) ' ' num2str(element_elecNeeded{j}(i,2)) ' ' num2str(element_elecNeeded{j}(i,3))]);
end
ii = ii + i;
end
fprintf(fid_out,'%s\n',s);
else
fprintf(fid_out,'%s\n',s);
end
end
fclose(fid_in);
fclose(fid_out);
end