-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsurf_recon.m
176 lines (153 loc) · 4.93 KB
/
surf_recon.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
% WIP. I am still not happy with these surface reconstruction methods. It may have to do with imperfect Laplace solutions.
clear; close all;
Laplace_radial = niftiread('Laplace_radial.nii.gz');
Laplace_tangential = niftiread('Laplace_tangential.nii.gz');
LBL = niftiread('labelmap-postProc.nii.gz');
hdr = niftiinfo('labelmap-postProc.nii.gz');
i = find(LBL>0);
sz = size(LBL);
[x,y,z] = ind2sub(sz,i);
Laplace_radial = Laplace_radial(i);
Laplace_tangential = Laplace_tangential(i);
LBL = LBL(i);
%% use cosine basis functions to describe concentric rings
%
% interval = linspace(0,1,100);
%
% figure('units','normalized','outerposition',[0 0 1 1])
% for band = 1:length(interval)-1
% d = find(Laplace_radial>interval(band) & Laplace_radial<interval(band+1));
% scatter3(x(d),y(d),z(d),1,Laplace_tangential(d),'.');
% xlim([748, 1464]);
% ylim([185, 1960]);
% zlim([284, 1605]);
% view(-90,0);
% drawnow;
%
% %extract this 'band' and sort by tangential coords
% [LP,i] = sort(Laplace_tangential(d));
% v = [x(d(i)) y(d(i)) z(d(i))];
% % make and fit a cosine basis design matrix with k cosine functions
% for k = 1:10
% DX(:,k) = cos(linspace(-pi*k,pi*k,length(v)));
% end
% beta = pinv(DX'*DX)*DX'*v;
% % reconstruct pointcloud
% vrec = DX*beta;
% plot3(vrec(:,1),vrec(:,2),vrec(:,3));
% end
%% reduce points (averaging bins) to a meshgrid
longitudebins = 100;
latittudebins = 100;
longlin = linspace(0,1,longitudebins);
latlin = linspace(0,1,latittudebins);
grid = nan(longitudebins,latittudebins,3);
for long = 1:longitudebins-1
for lat = 1:latittudebins-1
vertices = find(Laplace_radial>longlin(long) & Laplace_radial<longlin(long+1) & ...
Laplace_tangential>latlin(lat) & Laplace_tangential<latlin(lat+1));
try
grid(long,lat,1) = mean(x(vertices));
grid(long,lat,2) = mean(y(vertices));
grid(long,lat,3) = mean(z(vertices));
end
end
end
grid(:,end,:) = grid(:,1,:);
for dir = 1:3
grid(:,:,dir) = fillmissing(grid(:,:,dir),'linear');
end
% get face connectivity
tri = [1:(longitudebins*latittudebins)]';
F = [tri,tri+1,tri+(longitudebins) ; tri,tri-1,tri-(longitudebins)];
F = reshape(F',[3,longitudebins,latittudebins,2]);
F(:,longitudebins,:,1) = nan;
F(:,1,:,2) = nan;
F(:,:,latittudebins,1) = nan;
F(:,:,1,2) = nan;
F(isnan(F)) = [];
F=reshape(F,[3,(longitudebins-1)*(latittudebins-1)*2])';
% save as gifti
gii = gifti();
gii.faces = F;
gii.vertices = reshape(grid,[100*100,3]);
gii.vertices(:,4) = 0;
gii.vertices = gii.vertices * hdr.Transform.T;
gii.vertices(:,4) = [];
gii.vertices = gii.vertices+hdr.Transform.T(4,1:3);
save(gii,'meshgrid-100x100.surf.gii','Base64Binary');
figure;
vertices = patch('faces',F,'vertices',reshape(grid,[100*100,3]));
vertices.LineStyle = 'none';
%p.FaceColor = 'b';
material dull;
light;
hold on;
c = jet(latittudebins);
for lat = 1:latittudebins-1
plot3(grid(lat,:,1),grid(lat,:,2),grid(lat,:,3),'color',c(lat,:));
end
axis equal tight;
view(90,0);
saveas(gcf,'longitudinalLines-lat.png')
view(-90,0);
saveas(gcf,'longitudinalLines-med.png')
% try smoothing with cosineRep
% Vrecon = CosineRep_2Dsurf(grid,5,0.001);
% p = patch('faces',F,'vertices',Vrecon);
% p.LineStyle = 'none';
% p.FaceColor = 'b';
% material dull;
% light;
% hold on;
% axis equal tight;
%% reduce points to a meshgrid using interpolant
% see https://github.com/jordandekraker/Hippocampal_AutoTop/blob/master/tools/coords_SurfMap.m
lt = Laplace_tangential*2*pi - pi;
lr = 1-Laplace_radial;
[xx,yy] = pol2cart(lt,lr);
fd=@(p) sqrt(sum(p.^2,2)) -1;
[vertices,tri] = distmesh2d(fd,@huniform,0.01,[-1,-1;1,1],[]); % 36k vertices
scattInterp = scatteredInterpolant(xx,yy,x,'nearest','nearest');
xxx = scattInterp(vertices(:,1),vertices(:,2));
scattInterp = scatteredInterpolant(xx,yy,y,'nearest','nearest');
yyy = scattInterp(vertices(:,1),vertices(:,2));
scattInterp = scatteredInterpolant(xx,yy,z,'nearest','nearest');
zzz = scattInterp(vertices(:,1),vertices(:,2));
clear scattInterp;
Vxyz = [xxx yyy zzz];
figure;
p = patch('faces',tri,'vertices',Vxyz);
p.LineStyle = 'none';
p.FaceColor = 'b';
material dull;
light;
axis equal tight;
gii = gifti();
gii.faces = tri;
gii.vertices = Vxyz;
gii.vertices(:,4) = 0;
gii.vertices = gii.vertices * hdr.Transform.T;
gii.vertices(:,4) = [];
gii.vertices = gii.vertices+hdr.Transform.T(4,1:3);
save(gii,'distmesh-36k.surf.gii','Base64Binary');
gii.vertices = [vertices(:,1),vertices(:,2), zeros(length(vertices),1)];
save(gii,'distmesh-36k_unfolded.surf.gii','Base64Binary');
% smooth
FV = struct();
FV.faces = tri;
FV.vertices = Vxyz;
FV2 = smoothpatch(FV);
figure;
p = patch(FV2);
p.LineStyle = 'none';
p.FaceColor = 'b';
material dull;
light;
axis equal tight;
gii.vertices = FV2.vertices;
gii.vertices(:,4) = 0;
gii.vertices = gii.vertices * hdr.Transform.T;
gii.vertices(:,4) = [];
gii.vertices = gii.vertices+hdr.Transform.T(4,1:3);
save(gii,'distmesh-36k_smoothed.surf.gii','Base64Binary');