SurfStatSmooth

Smooths surface data by repeatedly averaging over edges.

 Usage: Y = SurfStatSmooth( Y, sv, FWHM );
 
 Y        = n x v or n x v x k matrix of surface data, v=#vertices;
            n=#observations; k=#variates, or memory map of same.
 surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
 or
 surf.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume). 
 FWHM     = approximate FWHM of Gaussian smoothing filter, in mesh units.

 Note that if the data is memory mapped, then the data is overwriten by
 the smoothed data.

0001 function Y = SurfStatSmooth( Y, surf, FWHM );
0002 
0003 %Smooths surface data by repeatedly averaging over edges.
0004 %
0005 % Usage: Y = SurfStatSmooth( Y, sv, FWHM );
0006 %
0007 % Y        = n x v or n x v x k matrix of surface data, v=#vertices;
0008 %            n=#observations; k=#variates, or memory map of same.
0009 % surf.tri = t x 3 matrix of triangle indices, 1-based, t=#triangles.
0010 % or
0011 % surf.lat = nx x ny x nz matrix, 1=in, 0=out, [nx,ny,nz]=size(volume).
0012 % FWHM     = approximate FWHM of Gaussian smoothing filter, in mesh units.
0013 %
0014 % Note that if the data is memory mapped, then the data is overwriten by
0015 % the smoothed data.
0016 
0017 niter=ceil(FWHM^2/(2*log(2)));
0018 
0019 if isnumeric(Y)
0020     [n,v,k]=size(Y);
0021     isnum=true;
0022 else
0023     Ym=Y;
0024     s=Ym.Format{2};
0025     if length(s)==2
0026         s=s([2 1]);
0027         k=1;
0028     else
0029         s=s([3 1 2]);
0030         k=s(3);
0031     end
0032     n=s(1);
0033     v=s(2);
0034     isnum=false;
0035 end
0036 
0037 edg=SurfStatEdg(surf);
0038 
0039 Y1=accumarray(edg(:,1),2,[v 1])'+accumarray(edg(:,2),2,[v 1])';
0040 
0041 if n>1
0042     fprintf(1,'%s',[num2str(n) ' x ' num2str(k) ' surfaces to smooth, % remaining: 100 ']);
0043 end
0044 n10=floor(n/10);
0045 for i=1:n
0046     if rem(i,n10)==0
0047         fprintf(1,'%s',[num2str(100-i/n10*10) ' ']);
0048     end
0049     for j=1:k
0050         if isnum
0051             Ys=squeeze(Y(i,:,j));
0052             for iter=1:niter
0053                 Yedg=Ys(edg(:,1))+Ys(edg(:,2));
0054                 Ys=(accumarray(edg(:,1),Yedg',[v 1]) + ...
0055                     accumarray(edg(:,2),Yedg',[v 1]))'./Y1;
0056             end
0057             Y(i,:,j)=Ys;
0058         else
0059             if length(s)==2
0060                 Y=Ym.Data(1).Data(:,i);
0061             else
0062                 Y=Ym.Data(1).Data(:,j,i);
0063             end            
0064             for iter=1:niter
0065                 Yedg=Y(edg(:,1))+Y(edg(:,2));
0066                 Y=(accumarray(edg(:,1),Yedg',[v 1]) + ...
0067                     accumarray(edg(:,2),Yedg',[v 1]))'./Y1;
0068             end
0069             if length(s)==2
0070                 Ym.Data(1).Data(:,i)=Y;
0071             else
0072                 Ym.Data(1).Data(:,j,i)=Y;
0073             end            
0074         end
0075     end
0076 end
0077 if n>1
0078     fprintf(1,'%s\n','Done');
0079 end
0080 if ~isnum
0081     Y=Ym;
0082 end
0083 
0084 return
0085 end