ImaGIN_spm_mesh_to_grid.m
4.06 KB
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
function R = ImaGIN_spm_mesh_to_grid(M, V, T)
% Non-linear interpolation of surface-based data onto a regular grid
% FORMAT R = spm_mesh_to_grid(M, V, T)
% M - a patch structure with fields 'faces' and 'vertices'
% V - an spm_vol structure with fields 'dim' and 'mat'
% T - array of data to be interpolated
%
% R - interpolated data on grid defined by V
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging
% Karl Friston, Guillaume Flandin
% $Id: spm_mesh_to_grid.m 4079 2010-10-07 11:41:54Z guillaume $
%-Precompute interpolation kernel and voxel grid
%==========================================================================
%-Compute a densely sampled triangular mask
%--------------------------------------------------------------------------
[tx, ty] = meshgrid(0.05:0.1:0.95, 0.05:0.1:0.95);
tx = tx(:);
ty = ty(:);
ti = find(sum([tx ty],2) <= 0.9);
tx = tx(ti);
ty = ty(ti);
tz = 1 - tx - ty;
%-Map the dense template triangle onto each face of the surface mesh
%--------------------------------------------------------------------------
P1 = M.vertices(M.faces(:,1),:);
P2 = M.vertices(M.faces(:,2),:);
P3 = M.vertices(M.faces(:,3),:);
If = speye(size(M.faces,1));
vertDense = kron(If,tx)*P2 + kron(If,ty)*P3 + kron(If,tz)*P1;
kernel = [tz tx ty];
%-Get voxel indices of all vertices from dense mesh within the image
%--------------------------------------------------------------------------
voxDense = V.mat\[vertDense';ones(1,size(vertDense,1))];
voxDense = round(voxDense(1:3,:)');
% Remove = unique([find(voxDense(:,1)<=0); find(voxDense(:,2)<=0);find(voxDense(:,3)<=0)]); %OD=in case vo
% Keep=setdiff(1:size(voxDense,1),Remove);
% voxDense=voxDense(Keep,:);%OD=in case vo
voxDense(find(voxDense<=0))=1; %OD, quick fix, don't find the solution
voxDense(find(voxDense(:,1)>V.dim(1)),1)=V.dim(1); %OD, quick fix, don't find the solution
voxDense(find(voxDense(:,2)>V.dim(2)),2)=V.dim(2); %OD, quick fix, don't find the solution
voxDense(find(voxDense(:,3)>V.dim(3)),3)=V.dim(3); %OD, quick fix, don't find the solution
voxInd = sub2ind(V.dim,voxDense(:,1),voxDense(:,2),voxDense(:,3));
voxInd = reshape(voxInd,size(kernel,1),[]);
%-Interpolation
%==========================================================================
integralConservation = true;
integralConservation = false; %OD
%-Normalise vertex data by the number of faces they are in
%--------------------------------------------------------------------------
if integralConservation
A = spm_mesh_adjacency(M);
T = spdiags(1./sum(A,2),0,size(A,1),size(A,1)) * T;
end
%-Interpolate each dense face data in the voxel grid
%--------------------------------------------------------------------------
R = zeros([V.dim size(T,2)]);
if ~integralConservation, countR = zeros(size(R)); end
for i = 1:size(M.faces,1)
%- Values at the vertices of the face
faceVal = T(M.faces(i,:),:);
% isData = any(faceVal);
if sum(isnan(T)>0) %OD
isData = sum(isnan(faceVal))<3;
else
isData = any(faceVal);
end
if any(isData)
%- Voxels corresponding to those vertices
[iV,I,J] = unique(voxInd(:,i));
for j=find(isData)
%- Values at the vertices of the dense face
faceValDense = kernel*faceVal(:,j);
if integralConservation
faceValDense = faceValDense * sum(faceVal(:,j))/sum(faceValDense);
end
%- Integrating values within voxels
k = (j-1) * prod(V.dim);
if integralConservation
R(iV+k) = R(iV+k) + accumarray(J,faceValDense);
else
R(iV+k) = R(iV+k) + accumarray(J,faceValDense,[],@mean);
countR(iV+k) = countR(iV+k) + 1;
end
end
end
end
if ~integralConservation
% countR(~countR) = 1; %OD
R = R ./ countR;
end
if size(T,2) == 1, R = R(:,:,:,1); end