function test_suite = test_surficial_neighborhood
% tests for cosmo_surficial_neighborhood
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;
function test_surficial_neighborhood_surface_dijkstra
if cosmo_skip_test_if_no_external('surfing')
return;
end
warning_state=warning();
warning_resetter=onCleanup(@()warning(warning_state));
warning('off');
opt=struct();
opt.progress=false;
ds=cosmo_synthetic_dataset('type','surface');
[vertices,faces]=get_synthetic_surface();
% dijkstra neighborhood fixed number of voxels
args={{vertices,faces},'count',4,'metric','dijkstra',opt};
nh1=cosmo_surficial_neighborhood(ds,args{:});
assertFalse(isfield(nh1.a,'vol'));
assert_equal_cell(nh1.neighbors,{ [ 1 2 4 3 ]
[ 2 1 3 5 ]
[ 3 2 6 5 ]
[ 4 1 5 2 ]
[ 5 2 4 6 ]
[ 6 3 5 2 ] });
assertEqual(nh1.fa.radius,[2 1 sqrt(2) sqrt(2) 1 2]);
assertEqual(nh1.fa.node_indices,1:6);
check_area(vertices,faces,nh1);
args={{vertices,faces},'radius',2.5,'metric','dijkstra',opt};
nh2=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh2.neighbors,{[1 2 3 4 5];
[1 2 3 4 5 6];
[1 2 3 4 5 6];
[1 2 3 4 5 6];
[1 2 3 4 5 6];
[2 3 4 5 6]; });
assertEqual(nh2.fa.radius,[2,2,1+sqrt(2),1+sqrt(2),2,2]);
assertEqual(nh2.fa.node_indices,1:6);
check_partial_neighborhood(ds,nh2,args);
check_area(vertices,faces,nh2);
args{1}{1}([2,6],:)=NaN;
nh3=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh3.neighbors,{[1 4 5];
[];
[3 4 5];
[1 3 4 5];
[1 3 4 5];
[]; });
assertEqual(nh3.fa.radius,[2,NaN,1+sqrt(2),1+sqrt(2),2,NaN]);
check_partial_neighborhood(ds,nh3,args);
check_area(args{1}{:},nh3);
args{2}='count';
args{3}=3;
nh4=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh4.neighbors,{[1 4 5];
[];
[3 5 4];
[4 1 5];
[5 4 3];
[]; });
assertEqual(nh4.fa.radius,[2,NaN,1+sqrt(2),1,sqrt(2),NaN]);
check_area(args{1}{:},nh4);
args{1}{1}=vertices;
args{1}{1}([2,5],:)=NaN; % split in two surfaces
args{3}=2;
nh5=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh5.neighbors,{[1 4];
[];
[3 6];
[4 1];
[];
[6 3]; });
assertEqual(nh5.fa.radius,[1,NaN,1,1,NaN,1]);
check_area(args{1}{:},nh5);
% throw error when too many nodes asked for
aet=@(varargin)assertExceptionThrown(@()...
cosmo_surficial_neighborhood(varargin{:}),'');
args{2}='count';
args{3}=3;
aet(ds,args{:})
function test_surficial_neighborhood_surface_direct
if cosmo_skip_test_if_no_external('surfing')
return;
end
warning_state=warning();
warning_resetter=onCleanup(@()warning(warning_state));
warning('off');
opt=struct();
opt.progress=false;
ds=cosmo_synthetic_dataset('type','surface');
[vertices,faces]=get_synthetic_surface();
% direct neighborhood
args={{vertices,faces},'direct',true,opt};
nh3=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh3.neighbors,{ [ 1 2 4 ]
[ 2 3 1 4 5 ]
[ 3 2 5 6 ]
[ 4 1 2 5 ]
[ 5 3 6 2 4 ]
[ 6 3 5 ] });
assertElementsAlmostEqual(nh3.fa.radius,sqrt([1 2 2 2 2 1]));
check_partial_neighborhood(ds,nh3,args);
args{1}{1}([2 5],:)=NaN;
nh4=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh4.neighbors,{ [ 1 4 ]
[]
[ 3 6 ]
[ 1 4 ]
[]
[6 3] });
check_partial_neighborhood(ds,nh4,args);
function test_surficial_neighborhood_surface_geodesic
if cosmo_skip_test_if_no_external('fast_marching') || ...
cosmo_skip_test_if_no_external('surfing')
return;
end
warning_state=warning();
warning_resetter=onCleanup(@()warning(warning_state));
warning('off');
opt=struct();
opt.progress=false;
ds=cosmo_synthetic_dataset('type','surface');%,'size','normal');
[vertices,faces]=get_synthetic_surface();
args={{vertices,faces},'count',4,opt};
nh=cosmo_surficial_neighborhood(ds,args{:});
assert_equal_cell(nh.neighbors,{[ 1 2 4 5 ]
[ 2 1 3 5 ]
[ 3 2 6 5 ]
[ 4 1 5 2 ]
[ 5 2 4 6 ]
[ 6 3 5 2 ] });
assertEqual(nh.fa.node_indices,1:6);
assertEqual(nh.fa.radius,[sqrt(.5)+1 1 sqrt(2) sqrt(2) 1 sqrt(.5)+1]);
vertices2=[NaN NaN NaN; vertices;NaN NaN NaN];
faces2=[faces+1; 1 1 8];
args={{vertices2,faces2},'count',4,opt};
nh2=cosmo_surficial_neighborhood(ds,args{:});
assertEqual(nh2.neighbors,{ zeros(1,0)
[ 2 3 5 6 ]
[ 3 2 4 6 ]
[ 4 3 6 2 ]
[ 5 2 6 3 ]
[ 6 3 5 4 ] });
function test_surficial_neighborhood_volume_geodesic
if cosmo_skip_test_if_no_external('fast_marching') || ...
cosmo_skip_test_if_no_external('surfing')
return;
end
warning_state=warning();
warning_resetter=onCleanup(@()warning(warning_state));
warning('off');
opt=struct();
opt.progress=false;
ds=cosmo_synthetic_dataset();
vertices=[-2 0 2 -2 0 2;...
-1 -1 -1 1 1 1;...
-1 -1 -1 -1 -1 -1]';
faces= [ 3 2 3 2;...
2 1 5 4;...
5 4 6 5 ]';
pial=vertices;
pial(:,3)=pial(:,3)+1;
white=vertices;
white(:,3)=white(:,3)-1;
nh1=cosmo_surficial_neighborhood(ds,{vertices,[-1 1],faces},...
'count',4,opt);
nh2=cosmo_surficial_neighborhood(ds,{pial,white,faces},...
'count',4,opt);
assert_equal_cell(nh1.neighbors,{[ 1 2 4 5 ]
[ 1 2 3 5 ]
[ 2 3 5 6 ]
[ 4 1 5 2 ]
[ 5 2 4 6 ]
[ 6 3 5 2 ] });
assertEqual(nh1.fa.node_indices,1:6);
assertEqual(nh1,nh2);
function test_surficial_neighborhood_volume_dijkstra
if cosmo_skip_test_if_no_external('surfing')
return;
end
warning_state=warning();
warning_resetter=onCleanup(@()warning(warning_state));
warning('off');
opt=struct();
opt.progress=false;
ds=cosmo_synthetic_dataset();
vertices=[-2 0 2 -2 0 2;
-1 -1 -1 1 1 1
-1 -1 -1 -1 -1 -1]';
faces= [ 3 2 3 2
2 1 5 4
5 4 6 5 ]';
pial=vertices;
pial(:,3)=pial(:,3)+1;
white=vertices;
white(:,3)=white(:,3)-1;
args3={{vertices,[-1 1],faces},'metric','dijkstra','count',4,opt};
args4={{pial,white,faces},'metric','dijkstra','count',4,opt};
nh3=cosmo_surficial_neighborhood(ds,args3{:});
nh4=cosmo_surficial_neighborhood(ds,args4{:});
assert_equal_cell(nh3.neighbors,{ [ 1 2 4 3 ]
[ 2 1 3 5 ]
[ 3 2 6 5 ]
[ 4 1 5 2 ]
[ 5 2 4 6 ]
[ 6 3 5 2 ] });
assertEqual(nh3.fa.node_indices,1:6);
assert_equal_cell(nh4.neighbors,nh3.neighbors);
assertFalse(isfield(nh3.a,'vol'));
% TODO
% check_partial_neighborhood(ds,nh3,args3);
% check_partial_neighborhood(ds,nh3,args3);
function test_surficial_neighborhood_exceptions
if cosmo_skip_test_if_no_external('surfing')
return
end
ds=cosmo_synthetic_dataset('type','surface');%,'size','normal');
[vertices,faces]=get_synthetic_surface();
aet=@(varargin)assertExceptionThrown(@()...
cosmo_surficial_neighborhood(varargin{:},...
'progress',false),'');
aet(ds,{vertices,faces});
% need surfaces
aet(ds,{},'radius',2);
% center_ids not supported for surface dataset
aet(ds,{vertices,faces},'radius',2,'center_ids',1);
% cannot have duplicate feature ids
ds_double=cosmo_stack({ds,ds},2);
ds_double.a.fdim.values{1}=[1:6,1:6];
aet(ds_double,{vertices,faces},'radius',2);
% outside range
ds_double.a.fdim.values{1}=[1:5,1:5];
aet(ds_double,{vertices,faces},'radius',2);
% cannot have fmri and surface dataset combined
ds2=cosmo_synthetic_dataset();
ds2.a.fdim.values{end+1}=ds.a.fdim.values{1};
ds2.a.fdim.labels{end+1}=ds.a.fdim.labels{1};
ds2.fa.node_indices=ds.fa.node_indices;
aet(ds2,{vertices,faces},'radius',2);
% cannot have MEEG dataset
ds_meeg=cosmo_synthetic_dataset('type','meeg');
aet(ds_meeg,{vertices,faces},'radius',2);
% need positive scalar radius
aet(ds,{vertices,faces},'radius',-1);
aet(ds,{vertices,faces},'radius',eye(2));
function check_partial_neighborhood(ds,nh,args)
% see if when have a partial dataset, the neighborbood reflects
% that too
nf=size(ds.samples,2);
rp=randperm(nf);
keep_count=round(nf*.7);
keep_sel=rp(1:keep_count);
keep_all=[keep_sel keep_sel keep_sel];
ds_sel=cosmo_slice(ds,keep_all,2);
fdim=ds_sel.a.fdim.values{1};
rp_fdim=randperm(numel(fdim));
ds_sel.a.fdim.values{1}=fdim(rp_fdim);
nh_sel=cosmo_surficial_neighborhood(ds_sel,args{:});
assertEqual(numel(nh_sel.neighbors), numel(keep_sel));
assertEqual(nh_sel.a.fdim.labels,nh.a.fdim.labels);
assertEqual(nh_sel.a.fdim.values{1},nh.a.fdim.values{1}(rp_fdim));
assertEqual(numel(nh_sel.a.fdim.values),numel(nh.a.fdim.values));
assertEqual(ds_sel.a,nh_sel.a);
opt=cosmo_structjoin(args(2:end));
if isfield(opt,'radius')
metric=opt.metric;
metric_arg=opt.radius;
elseif isfield(opt,'count')
metric=opt.metric;
metric_arg=[10 opt.count];
elseif isfield(opt,'direct')
metric='direct';
if opt.direct
metric_arg=NaN;
else
metric_arg=0;
end
else
assert(false);
end
faces=args{1}{2};
n2f=surfing_invertmapping(faces);
nodes_ds_sel=ds_sel.a.fdim.values{1}(ds_sel.fa.node_indices);
nodes_ds=ds.a.fdim.values{1}(ds.fa.node_indices);
nodes_nh_sel=nh_sel.a.fdim.values{1}(nh_sel.fa.node_indices);
vertices=args{1}{1};
nvertices=size(vertices,1);
nodes_kept=cosmo_match(1:nvertices,nodes_ds_sel);
vertices(~nodes_kept,:)=NaN;
node_mask=all(isfinite(vertices),2);
nodes_removed=setdiff(nodes_ds(:)',nodes_ds_sel(:)');
assertEqual(setxor(nodes_removed,nodes_ds_sel),1:nf);
nb_sel=nh_sel.neighbors;
for k=1:numel(nh_sel.neighbors)
sel_center_node=nodes_nh_sel(k);
idx=find(nodes_ds==sel_center_node);
assert(numel(idx)==1);
center_node=nodes_ds(idx);
assertEqual(sel_center_node, center_node);
switch metric
case 'direct'
if node_mask(sel_center_node)
direct_neighbors=surfing_surface_nbrs(faces',...
vertices');
around_nodes=direct_neighbors(sel_center_node,:);
msk=cosmo_match(around_nodes, ...
find(isfinite(vertices(:,1))));
% add node itself
around_nodes=[sel_center_node,...
around_nodes(msk & around_nodes>0)];
else
around_nodes=[];
end
otherwise
around_nodes=surfing_circleROI(vertices',faces',...
sel_center_node,metric_arg,metric,n2f);
end
sel_around_nodes=nodes_ds_sel(nb_sel{k});
if isempty(sel_around_nodes)
assertTrue(isempty(around_nodes));
else
assertEqual(unique(sel_around_nodes),...
setdiff(around_nodes, nodes_removed))
end
end
function test_surface_subsampling
if cosmo_skip_test_if_no_external('surfing')
return;
end
vertices=[0 -1 -2 -1 1 2 1 3 4 3;
0 -2 0 2 2 0 -2 2 0 -2;
0 0 0 0 0 0 0 0 0 0]';
faces=[1 1 1 1 1 1 5 8 6 6;
2 3 4 5 6 7 8 9 9 10;
3 4 5 6 7 2 6 6 10 7]';
% make custom volume
ds=cosmo_synthetic_dataset('size','small');
cp=cosmo_cartprod({1:7,1:3})';
ds.fa.i=cp(1,:);
ds.fa.j=cp(2,:);
ds.fa.k=ds.fa.i*0+1;
ds.a.fdim.values={1:7;1:3;1};
ds.samples=zeros(numel(ds.sa.targets),numel(ds.fa.i));
ds.a.vol.dim=cellfun(@numel,ds.a.fdim.values)';
ds.a.vol.mat(2,4)=-4;
ds.a.vol.mat(3,4)=-1;
ds.a.vol.mat(1,1)=1;
ds.a.vol.mat(2,2)=2;
ds.a.vol.mat(3,3)=1;
surfs={vertices,faces,[-4 5]};
opt=struct();
opt.progress=false;
opt.radius=3;
opt.metric='euclidean';
nh=cosmo_surficial_neighborhood(ds,surfs,opt);
assertEqual(nh.neighbors,{ [ 2 4 8 10 12 16 18 ]
[ 2 4 8 10 ]
[ 2 8 10 16 ]
[ 8 10 16 18 ]
[ 10 12 16 18 20 ]
[ 4 6 10 12 14 18 20 ]
[ 2 4 6 10 12 ]
[ 12 14 18 20 ]
[ 6 12 14 20 ]
[ 4 6 12 14 ] });
assertEqual(nh.origin.fa,ds.fa);
assertEqual(nh.origin.a,ds.a);
% test subsampling
subsample=2;
surfs={vertices,faces,[-4 5],subsample};
nh2=cosmo_surficial_neighborhood(ds,surfs,opt);
assertEqual(nh2.neighbors,{ [ 2 4 8 10 ]
[ 2 8 10 16 ]
[ 8 10 16 18 ]
[ 10 12 16 18 20 ]
[ 2 4 6 10 12 ]
[ 12 14 18 20 ]
[ 6 12 14 20 ]
[ 4 6 12 14 ] });
assertEqual(nh2.origin.fa,ds.fa);
assertEqual(nh2.origin.a,ds.a);
% subsampling with pial surface
pial=bsxfun(@plus,vertices,[0 0 1]);
white=bsxfun(@plus,vertices,[0 0 -1]);
[vo,fo]=surfing_subsample_surface(vertices,faces,2,.2,0);
surfs={pial,white,faces,vo,fo};
nh3=cosmo_surficial_neighborhood(ds,surfs,opt);
assertEqual(nh2,nh3);
% check center ids options
slice_ids=[5 3 2];
nh4=cosmo_surficial_neighborhood(ds,surfs,opt,'center_ids',slice_ids);
nh4_sl=struct();
nh4_sl.neighbors=nh3.neighbors(slice_ids);
nh4_sl.fa=cosmo_slice(nh3.fa,slice_ids,2,'struct');
nh4_sl.a=nh3.a;
assertEqual(nh4.a.fdim.values{1}(nh4.fa.node_indices),...
nh4_sl.a.fdim.values{1}(nh4_sl.fa.node_indices));
assertEqual(nh4.neighbors,nh4_sl.neighbors);
% try with file names
fn_pial=cosmo_make_temp_filename('pial','.asc');
fn_white=cosmo_make_temp_filename('white','.asc');
fn_tiny=cosmo_make_temp_filename('tiny','.asc');
cleaner1=onCleanup(@()delete(fn_pial));
cleaner2=onCleanup(@()delete(fn_white));
cleaner3=onCleanup(@()delete(fn_tiny));
surfing_write(fn_pial, pial, faces);
surfing_write(fn_white, white, faces);
surfing_write(fn_tiny, vo, fo);
surfs={fn_pial,fn_white,fn_tiny};
nh5=cosmo_surficial_neighborhood(ds,surfs,opt);
assertEqual(nh2,nh5);
% should work with alternative voldef
ds_bad_vol=ds;
ds_bad_vol.a.vol.mat(:)=NaN;
ds_bad_vol.a.vol.dim(:)=NaN;
nh6=cosmo_surficial_neighborhood(ds_bad_vol,surfs,opt,...
'vol_def',ds.a.vol);
nh6.origin.a.vol=ds.a.vol;
assertEqual(nh5,nh6);
% check exceptions
aet=@(varargin)assertExceptionThrown(@()...
cosmo_surficial_neighborhood(varargin{:}),'');
surfs={fn_pial,fn_pial,fn_tiny};
aet(ds,surfs,opt);
white_bad=white;
white_bad=white_bad(2:end,:);
aet(ds,{pial,white_bad,faces},opt);
% missing faces for output surface
aet(ds,{fn_pial,fn_white,vo},opt);
% face mismatch
faces_bad=faces;
faces_bad=faces_bad(end:-1:1,:);
surfing_write(fn_white, white, faces_bad);
aet(ds,{fn_pial,fn_white},opt);
% too many surf arguments
aet(ds,{fn_pial,fn_white,fn_tiny,fn_tiny},opt);
aet(ds,{pial,white,faces,pial,white,white},opt);
aet(ds,{pial,white,faces,fn_pial,white},opt);
% surfs are not a cell
aet(ds,struct,opt);
aet(ds,{pial,white,{}});
function check_area(vertices,faces,nh)
assert(isfield(nh.fa,'area'));
area=surfing_surfacearea(vertices,faces);
node_idxs=nh.fa.node_indices;
n_nodes=numel(node_idxs);
for k=1:n_nodes
nbr_idxs=node_idxs(nh.neighbors{k});
expected_area=sum(area(nbr_idxs));
if isnan(expected_area)
assertEqual(expected_area,nh.fa.area(k))
else
assertElementsAlmostEqual(expected_area,nh.fa.area(k));
end
end
function [vertices,faces]=get_synthetic_surface()
% return the following surface (face indices in [brackets])
%
% 1-----2-----3
% | /| /|
% |[2]/ |[1]/ |
% | / | / |
% | /[4]| /[3]|
% |/ |/ |
% 4-----5-----6
vertices=[0 0 0 1 1 1;
1 2 3 1 2 3;
0 0 0 0 0 0]';
faces= [ 3 2 3 2
2 1 5 4
5 4 6 5 ]';
function assert_equal_cell(x,y)
% small helper
assertEqual(size(x),size(y))
for k=1:numel(x)
xk=x{k};
yk=y{k};
if isempty(xk)
assertTrue(isempty(yk));
else
assertEqual(sort(xk),sort(yk));
end
end
