SurfStatView1

Viewer for surface and/or volume data or P-value or Q-value structures.
 
 Usage: SurfStatView( struct, surf, ... );
 
 struct        = 1 x v vector of data, v=#vertices,
                 or one of the following structures:
 P-value structure:
 struct.P      = 1 x v vector of corrected P-values for vertices.
 struct.C      = 1 x v vector of corrected P-values for clusters.
 struct.mask   = 1 x v, 1=inside, 0=outside, v=#vertices.
 struct.thresh = P-value threshold for plot, 0.05 by default.
 Q-value structure:
 struct.Q      = 1 x v matrix of Q-values. 
 struct.mask   = 1 x v, 1=inside, 0=outside, v=#vertices.
 struct.thresh = Q-value threshold for plot, 0.05 by default.

 surf is either a surface or volume structure:
 Surface structure:
 surf.coord   = 3 x v matrix of surface coordinates.
 surf.tri     = t x 3 matrix of triangle indices, 1-based, t=#triangles.
 Volume structure:
 surf.lat     = 3D logical array, 1=in, 0=out.
 surf.vox     = 1 x 3 vector of voxel sizes in mm, [1 1 1] by default.
 surf.origin  = position in mm of the first voxel, [0 0 0] by default.
 
 The following can be added as extra pairs of parameters in any order: 

 For slices of volume structures:
 'x', vector of slice positions in x direction, default mean(x coord).
 'y', vector of slice positions in y direction, default mean(y coord).
 'z', vector of slice positions in z direction, default mean(z coord).

 For isosurfaces of volume structures: 
 'maskthresh' mask threshold, default 0.5.
 'datathresh' data threshold, default [].
 'maskalpha', mask    alpha, default 0.2.
 'dataalpha', data    alpha, default 1.
 'clusalpha', cluster alpha, default 0.2.
 'maskcolor', mask    colour, any matlab ColorSpec, default 'green'.
 'datacolor', data    colour, any matlab ColorSpec, default 'red'.
 'cluscolor', cluster colour, any matlab ColorSpec, default 'blue'.

 'title', any string, data name by default.
 'background', background colour, any matlab ColorSpec, such as 
   'white' (default), 'black'=='k', 'r'==[1 0 0], [1 0.4 0.6] (pink) etc.
   Letter and line colours are inverted if background is dark (mean<0.5). 

 E.g. SurfStatView( struct, surf, 'x', 10, 'datathresh', 3, ...
          'title', 'One x slice at 10mm, thresholded at 3' );

 To change the colour limits to e.g. [2 4], use set( gca, 'CLim', [2 4] );
 To change the colour map, use e.g. colormap( 'jet' ). 
 Surfaces can be edited in the figure window by clicking e.g. "Rotate 3D".

0001 function SurfStatView1(struct, surf, varargin );
0002 
0003 %Viewer for surface and/or volume data or P-value or Q-value structures.
0004 %
0005 % Usage: SurfStatView( struct, surf, ... );
0006 %
0007 % struct        = 1 x v vector of data, v=#vertices,
0008 %                 or one of the following structures:
0009 % P-value structure:
0010 % struct.P      = 1 x v vector of corrected P-values for vertices.
0011 % struct.C      = 1 x v vector of corrected P-values for clusters.
0012 % struct.mask   = 1 x v, 1=inside, 0=outside, v=#vertices.
0013 % struct.thresh = P-value threshold for plot, 0.05 by default.
0014 % Q-value structure:
0015 % struct.Q      = 1 x v matrix of Q-values.
0016 % struct.mask   = 1 x v, 1=inside, 0=outside, v=#vertices.
0017 % struct.thresh = Q-value threshold for plot, 0.05 by default.
0018 %
0019 % surf is either a surface or volume structure:
0020 % Surface structure:
0021 % surf.coord   = 3 x v matrix of surface coordinates.
0022 % surf.tri     = t x 3 matrix of triangle indices, 1-based, t=#triangles.
0023 % Volume structure:
0024 % surf.lat     = 3D logical array, 1=in, 0=out.
0025 % surf.vox     = 1 x 3 vector of voxel sizes in mm, [1 1 1] by default.
0026 % surf.origin  = position in mm of the first voxel, [0 0 0] by default.
0027 %
0028 % The following can be added as extra pairs of parameters in any order:
0029 %
0030 % For slices of volume structures:
0031 % 'x', vector of slice positions in x direction, default mean(x coord).
0032 % 'y', vector of slice positions in y direction, default mean(y coord).
0033 % 'z', vector of slice positions in z direction, default mean(z coord).
0034 %
0035 % For isosurfaces of volume structures:
0036 % 'maskthresh' mask threshold, default 0.5.
0037 % 'datathresh' data threshold, default [].
0038 % 'maskalpha', mask    alpha, default 0.2.
0039 % 'dataalpha', data    alpha, default 1.
0040 % 'clusalpha', cluster alpha, default 0.2.
0041 % 'maskcolor', mask    colour, any matlab ColorSpec, default 'green'.
0042 % 'datacolor', data    colour, any matlab ColorSpec, default 'red'.
0043 % 'cluscolor', cluster colour, any matlab ColorSpec, default 'blue'.
0044 %
0045 % 'title', any string, data name by default.
0046 % 'background', background colour, any matlab ColorSpec, such as
0047 %   'white' (default), 'black'=='k', 'r'==[1 0 0], [1 0.4 0.6] (pink) etc.
0048 %   Letter and line colours are inverted if background is dark (mean<0.5).
0049 %
0050 % E.g. SurfStatView( struct, surf, 'x', 10, 'datathresh', 3, ...
0051 %          'title', 'One x slice at 10mm, thresholded at 3' );
0052 %
0053 % To change the colour limits to e.g. [2 4], use set( gca, 'CLim', [2 4] );
0054 % To change the colour map, use e.g. colormap( 'jet' ).
0055 % Surfaces can be edited in the figure window by clicking e.g. "Rotate 3D".
0056 
0057 mask=[];
0058 maskthresh=0.5;
0059 datathresh=[];
0060 maskalpha=0.2;
0061 dataalpha=1;
0062 clusalpha=0.2;
0063 maskcolor='green';
0064 datacolor='red';
0065 cluscolor='blue';
0066 xs=NaN;
0067 ys=NaN;
0068 zs=NaN;
0069 title1=inputname(1);
0070 background='white';
0071 
0072 optargin=size(varargin,2);
0073 for narg=1:2:optargin
0074     switch varargin{narg}
0075         case 'mask'
0076             mask=varargin{narg+1};
0077         case 'maskthresh'
0078             maskthresh=varargin{narg+1};
0079         case 'datathresh'
0080             datathresh=varargin{narg+1};
0081         case 'maskalpha'
0082             maskalpha=varargin{narg+1};
0083         case 'dataalpha'
0084             dataalpha=varargin{narg+1};
0085         case 'clusalpha'
0086             clusalpha=varargin{narg+1};
0087         case 'maskcolor'
0088             maskcolor=varargin{narg+1};
0089         case 'datacolor'
0090             datacolor=varargin{narg+1};
0091         case 'cluscolor'
0092             cluscolor=varargin{narg+1};
0093         case 'x'
0094             xs=varargin{narg+1};
0095         case 'y'
0096             ys=varargin{narg+1};
0097         case 'z'
0098             zs=varargin{narg+1};
0099         case 'title'
0100             title1=varargin{narg+1};
0101         case 'background'
0102             background=varargin{narg+1};
0103         otherwise
0104             warning(['Unrecognized optional argument ' varargin{narg} ' ignored :-(']);
0105     end
0106 end
0107 
0108 hold on;
0109 
0110 if isfield(surf,'coord')
0111     if isempty(struct)
0112         struct=zeros(1,size(surf.coord,2));
0113     end
0114     
0115     if ~isstruct(struct)
0116         clim=[min(struct),max(struct)];
0117         if clim(1)==clim(2)
0118             clim=clim(1)+[-1 0];
0119         end
0120         a=trisurf(surf.tri,surf.coord(1,:),surf.coord(2,:),surf.coord(3,:),...
0121             double(struct),'EdgeColor','none','FaceColor','interp');  
0122         colorbar; colormap(spectral(256));
0123         set(gca,'CLim',clim);
0124         set(gca,'Tag','SurfStatView 1 0');
0125         dcm_obj=datacursormode(gcf);
0126         set(dcm_obj,'UpdateFcn',@SurfStatDataCursor,'DisplayStyle','window');
0127     else
0128         if isfield(struct,'P')
0129             if ~isfield(struct,'thresh')
0130                 struct.thresh=0.05;
0131             end
0132 
0133             signifpeak=struct.P<struct.thresh;
0134             if isfield(struct,'C')
0135                 signifclus=struct.C<struct.thresh;
0136                 t1=signifclus.*(1-signifpeak).*(127-struct.C/struct.thresh*126);
0137             else
0138                 signifclus=0;
0139                 t1=0;
0140             end
0141             t2=signifpeak.*(255-struct.P/struct.thresh*126);
0142             t3=(1-signifpeak).*(1-signifclus)*128;
0143             tt=(t1+t2+t3).*struct.mask*struct.thresh;
0144 
0145             trisurf(surf.tri,surf.coord(1,:),surf.coord(2,:),surf.coord(3,:),...
0146                 double(tt),'EdgeColor','none');
0147             set(gca,'Tag','SurfStatView 1 0');
0148             
0149             colorbar; colormap(spectral(256));
0150             cm=[zeros(1,3);
0151                 zeros(127,1)   (0:126)'/127   ones(127,1); ones(1,3)*0.8;
0152                 ones(127,1)    (0:126)'/127   zeros(127,1)];
0153             SurfStatColormap(cm);
0154             cb = SurfStatColLim( [0 255]*struct.thresh );
0155 
0156             set(cb,'YLim',[0 255]*struct.thresh);
0157             h=get(cb,'Children');
0158             set(h,'YData',[0 255]*struct.thresh);
0159             set(cb,'YTick',[1+(0:5)/5*126 129+(0:5)/5*126]*struct.thresh);
0160             pstr=strvcat(num2str(struct.thresh*(5:-1:1)'/5),'       0');
0161             set(cb,'YTickLabel',strvcat(pstr,pstr));
0162             xl=get(cb,'YLabel');
0163             set(xl,'String','P Cluster               P Vertex');
0164             
0165             dcm_obj=datacursormode(gcf);
0166             set(dcm_obj,'UpdateFcn',@SurfStatDataCursorQ,'DisplayStyle','window');
0167         end
0168         if isfield(struct,'Q')
0169             if ~isfield(struct,'thresh')
0170                 struct.thresh=0.05;
0171             end
0172 
0173             t1=(struct.Q<struct.thresh).*(255-struct.Q/struct.thresh*253);
0174             t2=(struct.Q>=struct.thresh);
0175             tt=(t1+t2).*struct.mask*struct.thresh;
0176 
0177             trisurf(surf.tri,surf.coord(1,:),surf.coord(2,:),surf.coord(3,:),...
0178                 double(tt),'EdgeColor','none');
0179             set(gca,'Tag','SurfStatView 1 0');
0180 
0181             colorbar; colormap(spectral(256));
0182             cm=[zeros(1,3); ones(1,3)*0.8; ones(254,1) (0:253)'/254 zeros(254,1)];
0183             SurfStatColormap(cm);
0184             cb = SurfStatColLim( [0 255]*struct.thresh );
0185 
0186             set(cb,'YLim',[0 255]*struct.thresh);
0187             h=get(cb,'Children');
0188             set(h,'YData',[0 255]*struct.thresh);
0189             set(cb,'YTick',(2+(0:5)/5*253)*struct.thresh);
0190             set(cb,'YTickLabel',num2str(struct.thresh*(5:-1:0)'/5));
0191             xl=get(cb,'YLabel');
0192             set(xl,'String','Q');
0193             
0194             dcm_obj=datacursormode(gcf);
0195             set(dcm_obj,'UpdateFcn',@SurfStatDataCursorP,'DisplayStyle','window');
0196         end
0197     end
0198 end
0199 
0200 if isfield(surf,'lat')
0201     if ~isfield(surf,'vox')
0202         surf.vox=ones(1,3);
0203     end
0204     if ~isfield(surf,'origin');
0205         surf.origin=zeros(1,3);
0206     end    
0207     dim=size(surf.lat);
0208     x1=surf.origin(1)+((1:dim(1))-1)*surf.vox(1);
0209     y1=surf.origin(2)+((1:dim(2))-1)*surf.vox(2);
0210     z1=surf.origin(3)+((1:dim(3))-1)*surf.vox(3);
0211     [x,y,z]=meshgrid(x1,y1,z1);
0212 
0213     % Do a little bit of smoothing with FWHM=1.15 voxels;
0214     % since max(filt)>0.5, then mask>0.5 should still be the same as mask:
0215     [i,j,k]=ndgrid(-1:1,-1:1,-1:1);
0216     FWHM=1.15;
0217     filt=2.^(-4*(i.^2+j.^2+k.^2)/FWHM^2);
0218     filt=filt/sum(filt(:));
0219 
0220     if ~isstruct(struct)
0221         T=zeros(dim)+min(struct(:));
0222         T(surf.lat)=struct;
0223         T=permute(T,[2 1 3]);
0224         if (isempty(datathresh) & isempty(mask)) | ~isnan(xs) | ~isnan(ys) | ~isnan(zs) 
0225             if isnan(xs)
0226                 xs=mean(x1);
0227             end
0228             if isnan(ys)
0229                 ys=mean(y1);
0230             end
0231             if isnan(zs)
0232                 zs=mean(z1);
0233             end
0234             s=slice(x,y,z,T,xs,ys,zs);
0235             set(s,'EdgeColor','none','FaceColor','interp');
0236             colorbar; colormap(spectral(256));
0237         end
0238         if ~isempty(datathresh)
0239             p1=patch(isosurface(x,y,z,T,datathresh));
0240             isonormals(x,y,z,T,p1);
0241             set(p1,'FaceColor',datacolor,'EdgeColor','none','FaceAlpha',dataalpha);
0242         end
0243     end
0244 
0245     if isfield(struct,'mask') | ~isempty(mask)
0246         M=zeros(dim);
0247         if isfield(struct,'mask') 
0248             M(surf.lat)=struct.mask;
0249         else
0250             M(surf.lat)=mask;
0251         end            
0252         M=convn(M,filt,'same');
0253         M=permute(M,[2 1 3]);
0254         p1=patch(isosurface(x,y,z,M,maskthresh));
0255         isonormals(x,y,z,M,p1);
0256         set(p1,'FaceColor',maskcolor,'EdgeColor','none','FaceAlpha',maskalpha);
0257     end
0258 
0259     if isfield(struct,'P')
0260         if ~isfield(struct,'thresh')
0261             struct.thresh=0.05;
0262         end
0263         P=ones(dim);
0264         P(surf.lat)=struct.P.*struct.mask+(1-struct.mask);
0265         P=permute(P,[2 1 3]);
0266         p2=patch(isosurface(x,y,z,P,struct.thresh));
0267         isonormals(x,y,z,P,p2)
0268         set(p2,'FaceColor',datacolor,'EdgeColor','none','FaceAlpha',dataalpha);
0269 
0270         if isfield(struct,'C')
0271             C=zeros(dim);
0272             C(surf.lat)=(struct.C<=struct.thresh).*struct.mask;
0273             C=convn(C,filt,'same');
0274             C=permute(C,[2 1 3]);
0275             p3=patch(isosurface(x,y,z,C,0.5));
0276             isonormals(x,y,z,C,p3)
0277             set(p3,'FaceColor',cluscolor,'EdgeColor','none','FaceAlpha',clusalpha);
0278         end
0279     end
0280 
0281     if isfield(struct,'Q')
0282         if ~isfield(struct,'thresh')
0283             struct.thresh=0.05;
0284         end
0285         Q=ones(dim);
0286         Q(surf.lat)=struct.Q.*struct.mask+(1-struct.mask);
0287         Q=permute(Q,[2 1 3]);
0288         p2=patch(isosurface(x,y,z,Q,struct.thresh));
0289         isonormals(x,y,z,Q,p2)
0290         set(p2,'FaceColor',datacolor,'EdgeColor','none','FaceAlpha',dataalpha);
0291     end
0292 end
0293 
0294 xlabel('x'); ylabel('y'); zlabel('z');
0295 view(35,25);
0296 daspect([1 1 1]); axis tight; axis vis3d; grid on;
0297 lighting phong; material shiny; camlight; 
0298 
0299 hold off;
0300 
0301 title(title1);
0302 
0303 whitebg(gcf,background);
0304 set(gcf,'Color',background,'InvertHardcopy','off');
0305 
0306 set(gcf,'PaperPosition',[0.25 2.5 6 4.5]);
0307 
0308 return
0309 end