Finds peaks (local maxima) and clusters for surface data.
Usage: [ peak, clus, clusid ] = SurfStatPeakClus( slm, mask, thresh ...
[, reselspvert [, edg ] ] );
slm.t = l x v matrix of data, v=#vertices; the first row
slm.t(1,:) is used for the clusters, and the other
rows are used to calculate cluster resels if k>1. See
SurfStatF for the precise definition of the extra rows.
slm.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
or
slm.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume).
mask = 1 x v vector, 1=inside, 0=outside, v=#vertices.
thresh = clusters are vertices where slm.t(1,mask)>=thresh.
reselspvert = 1 x v vector of resels per vertex, default ones(1,v).
edg = e x 2 matrix of edge indices, 1-based, e=#edges.
The following are optional:
slm.df = degrees of freedom - only the length (1 or 2) is used
to determine if slm.t is Hotelling's T or T^2 when k>1.
slm.k = k=#variates, 1 by default.
peak.t = np x 1 vector of peaks (local maxima).
peak.vertid = np x 1 vector of vertex id's (1-based).
peak.clusid = np x 1 vector of cluster id's that contain the peak.
clus.clusid = nc x 1 vector of cluster id numbers
clus.nverts = nc x 1 vector of number of vertices in the cluster.
clus.resels = nc x 1 vector of resels in the cluster.
clusid = 1 x v vector of cluster id's for each vertex.0001 function [ peak, clus, clusid ] = SurfStatPeakClus( slm, mask, thresh, ... 0002 reselspvert, edg ); 0003 0004 %Finds peaks (local maxima) and clusters for surface data. 0005 % 0006 % Usage: [ peak, clus, clusid ] = SurfStatPeakClus( slm, mask, thresh ... 0007 % [, reselspvert [, edg ] ] ); 0008 % 0009 % slm.t = l x v matrix of data, v=#vertices; the first row 0010 % slm.t(1,:) is used for the clusters, and the other 0011 % rows are used to calculate cluster resels if k>1. See 0012 % SurfStatF for the precise definition of the extra rows. 0013 % slm.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles. 0014 % or 0015 % slm.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume). 0016 % mask = 1 x v vector, 1=inside, 0=outside, v=#vertices. 0017 % thresh = clusters are vertices where slm.t(1,mask)>=thresh. 0018 % reselspvert = 1 x v vector of resels per vertex, default ones(1,v). 0019 % edg = e x 2 matrix of edge indices, 1-based, e=#edges. 0020 % The following are optional: 0021 % slm.df = degrees of freedom - only the length (1 or 2) is used 0022 % to determine if slm.t is Hotelling's T or T^2 when k>1. 0023 % slm.k = k=#variates, 1 by default. 0024 % 0025 % peak.t = np x 1 vector of peaks (local maxima). 0026 % peak.vertid = np x 1 vector of vertex id's (1-based). 0027 % peak.clusid = np x 1 vector of cluster id's that contain the peak. 0028 % clus.clusid = nc x 1 vector of cluster id numbers 0029 % clus.nverts = nc x 1 vector of number of vertices in the cluster. 0030 % clus.resels = nc x 1 vector of resels in the cluster. 0031 % clusid = 1 x v vector of cluster id's for each vertex. 0032 0033 if nargin<5 0034 edg=SurfStatEdg(slm); 0035 end 0036 0037 [l,v]=size(slm.t); 0038 0039 slm.t(1,~mask)=min(slm.t(1,:)); 0040 t1=slm.t(1,edg(:,1)); 0041 t2=slm.t(1,edg(:,2)); 0042 islm=ones(1,v); 0043 islm(edg(t1<t2,1))=0; 0044 islm(edg(t2<t1,2))=0; 0045 lmvox=find(islm); 0046 lmt=slm.t(1,lmvox); 0047 0048 excurset=(slm.t(1,:)>=thresh); 0049 n=sum(excurset); 0050 0051 if n<1 0052 peak=[]; 0053 clus=[]; 0054 clusid=[]; 0055 return 0056 end 0057 0058 voxid=cumsum(excurset); 0059 edg=voxid(edg(all(excurset(edg),2),:)); 0060 0061 % Find cluster id's in nf (from Numerical Recipes in C, page 346): 0062 nf=1:n; 0063 for el=1:size(edg,1) 0064 j=edg(el,1); 0065 k=edg(el,2); 0066 while nf(j)~=j j=nf(j); end 0067 while nf(k)~=k k=nf(k); end 0068 if j~=k nf(j)=k; end 0069 end 0070 for j=1:n 0071 while nf(j)~=nf(nf(j)) nf(j)=nf(nf(j)); end 0072 end 0073 0074 % find the unique cluster id's corresponding to the local maxima: 0075 vox=find(excurset); 0076 ivox=find(ismember(vox,lmvox)); 0077 clmid=nf(ivox); 0078 [uclmid,iclmid,jclmid]=unique(clmid); 0079 0080 % find their volumes: 0081 ucid=unique(nf); 0082 nclus=length(ucid); 0083 ucvol=histc(nf,ucid); 0084 0085 % find their resels: 0086 if nargin<4 0087 reselsvox=ones(size(vox)); 0088 else 0089 reselsvox=reselspvert(vox); 0090 end 0091 nf1=interp1([0 ucid],0:nclus,nf,'nearest'); 0092 0093 % if k>1, find volume of cluster in added sphere: 0094 if ~isfield(slm,'k') | slm.k==1 0095 ucrsl=accumarray(nf1',reselsvox)'; 0096 end 0097 if isfield(slm,'k') & slm.k==2 0098 if l==1 0099 ndf=length(slm.df); 0100 r=2*acos((thresh./slm.t(1,vox)).^(1/ndf)); 0101 else 0102 r=2*acos(sqrt((thresh-slm.t(2,vox)).*(thresh>=slm.t(2,vox))./ ... 0103 (slm.t(1,vox)-slm.t(2,vox)))); 0104 end 0105 ucrsl=accumarray(nf1',r'.*reselsvox')'; 0106 end 0107 if isfield(slm,'k') & slm.k==3 0108 if l==1 0109 ndf=length(slm.df); 0110 r=2*pi*(1-(thresh./slm.t(1,vox)).^(1/ndf)); 0111 else 0112 nt=20; 0113 theta=((1:nt)'-1/2)/nt*pi/2; 0114 s=(cos(theta).^2)*slm.t(2,vox); 0115 if l==3 0116 s=s+(sin(theta).^2)*slm.t(3,vox); 0117 end 0118 r=2*pi*mean(1-sqrt((thresh-s).*(thresh>=s)./ ... 0119 (ones(nt,1)*slm.t(1,vox)-s))); 0120 end 0121 ucrsl=accumarray(nf1',r'.*reselsvox')'; 0122 end 0123 0124 % and their ranks (in ascending order): 0125 [sortucrsl,iucrsl]=sort(ucrsl); 0126 rankrsl=zeros(1,nclus); 0127 rankrsl(iucrsl)=nclus:-1:1; 0128 0129 % add these to lm as extra columns: 0130 lmid=lmvox(ismember(lmvox,vox)); 0131 lm=flipud(sortrows([slm.t(1,lmid)' lmid' rankrsl(jclmid)'],1)); 0132 cl=sortrows([rankrsl' ucvol' ucrsl'],1); 0133 0134 clusid=zeros(1,v); 0135 clusid(vox)=interp1([0 ucid],[0 rankrsl],nf,'nearest'); 0136 0137 peak.t=lm(:,1); 0138 peak.vertid=lm(:,2); 0139 peak.clusid=lm(:,3); 0140 clus.clusid=cl(:,1); 0141 clus.nverts=cl(:,2); 0142 clus.resels=cl(:,3); 0143 0144 return 0145 end 0146 0147