atlasblobs/make_atlas_blobs.m at master · kjamison/atlasblobs · GitHub

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
function atlasblobs = make_atlas_blobs(atlasfile,varargin)
% atlasblobs = make_atlas_blobs(atlasfile,'param',value,...)
%
% Required inputs:
% atlasfile = MRI volume file (eg nii.gz) with labeled voxels
%
% Optional inputs:
% backgroundfile = MRI volume file to extract a midsagittal slice from to
%   display in background (must be in same coordinate space as atlasfile)
% roinames = list of ROI names for later reference
% roihemis = list of 'lh' or 'rh' that will be used to decide which views
%   display each ROI on (default = determine automatically from ROI
%   centers)
% volumesmoothing = Box size for voxel smoothing on
%   each ROI before extracting surface? (helpful especially for small ROIs,
%   but downside = ROIs get larger/puffier and might overlap)
% surfsmoothing = iterations of surface vertex smoothing after extracting
%   from volume. Less necessary if volume smoothing performed (0=none)
%
% Output:
% atlasblobs = struct with all the info needed for displaying surface
% rendering later

args = inputParser;
args.addParameter('backgroundfile',[]);
args.addParameter('roinames',[]);
args.addParameter('roihemis',[]);
args.addParameter('volumesmoothing',false);
args.addParameter('surfacesmoothing',0);

args.parse(varargin{:});
args = args.Results;


atlasblobs=struct();

atlasblobs.backgroundslice=[];
atlasblobs.backgroundposition=[];
if(~isempty(args.backgroundfile) && exist(args.backgroundfile))

    [Vbg,~,bgscales,~,~,bghdr] = read_avw_and_header(args.backgroundfile);
    Vbg_x=squeeze(Vbg(round(end/2),:,:));

    atlasblobs.backgroundslice=Vbg_x;

    bgijk=[size(Vbg,1)/2 1 1;
        size(Vbg,1)/2 size(Vbg,2) 1;
        size(Vbg,1)/2 size(Vbg,2) size(Vbg,3);
        size(Vbg,1)/2 1 size(Vbg,3)];
    bgxyz=affine_transform(bghdr.mat,bgijk);

    atlasblobs.backgroundposition=bgxyz;
end

if(ischar(atlasfile) && exist(atlasfile,'file'))
    [V,~,scales,~,~,hdr] =read_avw_and_header(atlasfile);
elseif(isstruct(atlasfile) && isfield(atlasfile,'hdr'))
    V=atlasfile.V;
    hdr=atlasfile.hdr;
end
uval=unique(V(V>0));
FV_all={};
fprintf('Generating %d isosurfaces:\n',numel(uval));
for i = 1:numel(uval)
    fprintf('%d ',i);
    if(mod(i,10)==0)
        fprintf('\n');
    end
    Vtmp=V==uval(i);
    %isoval=0.5;
    if(args.volumesmoothing > 0)
        %Vtmp=smooth3(Vtmp);
        Vtmp=smooth3(Vtmp,'box',max(3,args.volumesmoothing));
        %isoval=.5*max(Vtmp(:))
    end
    FV = isosurface(Vtmp);
    %FV = isosurface(Vtmp,isoval);


    FV.vertices=FV.vertices(:,[2 1 3]); %isosurface seems to swap X and Y for some reason

    FV.vertices=affine_transform(hdr.mat,FV.vertices);

    if(args.surfacesmoothing>0)
        FV.conn = vertex_neighbours(FV);
        FV.vertices = mesh_smooth_vertices(FV.vertices, FV.conn,[],args.surfacesmoothing);
    end
    FV_all{i}=FV;
end


fprintf('Done!\n');

atlasblobs.FV=FV_all;

atlasblobs.roilabels=uval;
if(~isempty(args.roinames))
    atlasblobs.roinames=args.roinames;
else
    atlasblobs.roinames={};
end

roimean=cellfun(@(x)mean(x.vertices,1),FV_all,'uniformoutput',false);
roimean=cat(1,roimean{:});
atlasblobs.roicenters=roimean;

if(~isempty(args.roihemis))
    atlasblobs.hemisphere=args.roihemis;
else

    roimin=cellfun(@(x)min(x.vertices,[],1),FV_all,'uniformoutput',false);
    roimin=cat(1,roimin{:});
    roimax=cellfun(@(x)max(x.vertices,[],1),FV_all,'uniformoutput',false);
    roimax=cat(1,roimax{:});

    %atlasblobs.hemisphere=repmat({'both'},numel(FV_all),1);
    %atlasblobs.hemisphere(roimax(:,1)<0)={'lh'};
    %atlasblobs.hemisphere(roimin(:,1)>0)={'rh'};
    atlasblobs.hemisphere(roimean(:,1)<0)={'lh'};
    atlasblobs.hemisphere(roimean(:,1)>=0)={'rh'};

end

seqidx=(1:max(uval))';
if(~isequal(uval,seqidx))
    if(numel(uval)~=numel(atlasblobs.roinames) && numel(atlasblobs.roinames)==numel(seqidx))
        atlasblobs.roinames=atlasblobs.roinames(ismember(seqidx,uval));
    end
    if(numel(uval)~=numel(atlasblobs.hemisphere) && numel(atlasblobs.hemisphere)==numel(seqidx))
        atlasblobs.hemisphere=atlasblobs.hemisphere(ismember(seqidx,uval));
    end
end