SurfStatResels

Resels of surface or volume data inside a mask.

 Usage: [resels, reselspvert, edg] = SurfStatResels( slm [, mask] )

 slm.resl = e x k matrix of sum over observations of squares of
            differences of normalized residuals along each edge.
 slm.tri  = t x 3 matrix of triangle indices, 1-based, t=#triangles.
 or
 slm.lat  = 3D logical array, 1=in, 0=out.
 mask     = 1 x v, 1=inside, 0=outside, v=#vertices,
          = ones(1,v), i.e. the whole surface, by default.

 resels      = 1 x (D+1) vector of 0,...,D dimensional resels of the mask,
             = EC of the mask if slm.resl is not given.
 reselspvert = 1 x v vector of D-dimensional resels per mask vertex.
 edg         = e x 2 matrix of edge indices, 1-based, e=#edges.

0001 function [resels, reselspvert, edg] = SurfStatResels( slm, mask );
0002 
0003 %Resels of surface or volume data inside a mask.
0004 %
0005 % Usage: [resels, reselspvert, edg] = SurfStatResels( slm [, mask] )
0006 %
0007 % slm.resl = e x k matrix of sum over observations of squares of
0008 %            differences of normalized residuals along each edge.
0009 % slm.tri  = t x 3 matrix of triangle indices, 1-based, t=#triangles.
0010 % or
0011 % slm.lat  = 3D logical array, 1=in, 0=out.
0012 % mask     = 1 x v, 1=inside, 0=outside, v=#vertices,
0013 %          = ones(1,v), i.e. the whole surface, by default.
0014 %
0015 % resels      = 1 x (D+1) vector of 0,...,D dimensional resels of the mask,
0016 %             = EC of the mask if slm.resl is not given.
0017 % reselspvert = 1 x v vector of D-dimensional resels per mask vertex.
0018 % edg         = e x 2 matrix of edge indices, 1-based, e=#edges.
0019 
0020 % The first version of this function (below) was much more elegant.
0021 % There were two reasons why I had to abandon it:
0022 % 1) for volumetric data, the edg, tri and tet matrices become huge, often
0023 % exceeding the memory capacity - they require 62 integers per voxel, as
0024 % opposed to 12 integers per vertex for surface data.
0025 % 2) even if memory is not exceeded, the "unique" and "ismember" functions
0026 % take a huge amount of time; the original version is ~7 times slower.
0027 % To overcome 1), I had to resort to a slice-wise implementation for
0028 % volumetric data. To cut down storage, the edg, tri and tet matrices are
0029 % formed only for two adjacent slices in the volume. The slm.lat data is
0030 % then copied into one slice, then the other, in an alternating fashion.
0031 % Since the tetrahedral filling is also alternating, there is no need to
0032 % recompute the edg, tri and tet matrices. However the ordering of the edg
0033 % matrix is crucial - it must match the SurfStatEdg function so that
0034 % slm.resl is synchronised. The ordering of the tri and tet matrices is not
0035 % important.
0036 % To avoid the "unique" function in 2), the edg and tri matrices were
0037 % created by hand, rather than derived from the tet matrix. To avoid the
0038 % "ismember" function, edge look-up was implemented by a sparse matrix.
0039 % This reduces execution time to ~18%. However it will not work for large
0040 % dimensions, since the largest index of a sparse matrix is 2^31. If this
0041 % is exceeded, then "interp1" is used which is slower but still reduces the
0042 % time to ~45% that of "ismember". Traingle area look-up was avoided by
0043 % simply re-computing trangle areas from the edges.
0044 %
0045 % Here is the original SurfStatResels function:
0046 %
0047 % if isfield(surf,'tri')
0048 %     tri=sort(surf.tri,2);
0049 %     edg=unique([tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])],'rows');
0050 %     tet=[];
0051 % end
0052 % if isfield(surf,'lat')
0053 %     % Fill lattice with 5 tetrahedra per cube
0054 %     [I,J,K]=size(surf.lat);
0055 %     [i,j,k]=ndgrid(1:int32(I),1:int32(J),1:int32(K));
0056 %     i=i(:);
0057 %     j=j(:);
0058 %     k=k(:);
0059 %     c1=find(rem(i+j+k,2)==0 & i<I & j<J & k<K);
0060 %     c2=find(rem(i+j+k,2)==0 & i>1 & j<J & k<K);
0061 %     clear i j k
0062 %     IJ=I*J;
0063 %     tet=[c1     c1+1    c1+1+I    c1+1+IJ;
0064 %          c1     c1+I    c1+1+I    c1+I+IJ;
0065 %          c1     c1+1+I  c1+1+IJ   c1+I+IJ;
0066 %          c1     c1+IJ   c1+1+IJ   c1+I+IJ;
0067 %          c1+1+I c1+1+IJ c1+I+IJ   c1+1+I+IJ;
0068 %          c2-1   c2      c2-1+I    c2-1+IJ;
0069 %          c2     c2-1+I  c2+I      c2+I+IJ;
0070 %          c2     c2-1+I  c2-1+IJ   c2+I+IJ;
0071 %          c2     c2-1+IJ c2+IJ     c2+I+IJ;
0072 %          c2-1+I c2-1+IJ c2-1+I+IJ c2+I+IJ];
0073 %     clear c1 c2
0074 %     % Find triangles and edges
0075 %     tri=unique([tet(:,[1 2 3]); tet(:,[1 2 4]); ...
0076 %                 tet(:,[1 3 4]); tet(:,[2 3 4])],'rows');
0077 %     edg=unique([tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])],'rows');
0078 %     % index by voxels in the lat
0079 %     vid=cumsum(surf.lat(:)).*surf.lat(:);
0080 %     % only inside the lat
0081 %     edg=vid(edg(all(surf.lat(edg),2),:));
0082 %     tri=vid(tri(all(surf.lat(tri),2),:));
0083 %     tet=vid(tet(all(surf.lat(tet),2),:));
0084 % end
0085 %
0086 % m=sum(mask(:));
0087 % lkc(1,1)=m;
0088 % %% LKC of edges
0089 % maskedg=all(mask(edg),2);
0090 % lkc(1,2)=sum(maskedg);
0091 % if isfield(slm,'resl')
0092 %     r1=mean(sqrt(slm.resl(maskedg,:)),2);
0093 %     lkc(2,2)=sum(r1);
0094 % end
0095 %
0096 % %% LKC of triangles
0097 % masktri=all(mask(tri),2);
0098 % lkc(1,3)=sum(masktri);
0099 % if isfield(slm,'resl')
0100 %     [tf,loc]=ismember(tri(masktri,[1 2]),edg,'rows');  l12=slm.resl(loc,:);
0101 %     [tf,loc]=ismember(tri(masktri,[1 3]),edg,'rows');  l13=slm.resl(loc,:);
0102 %     [tf,loc]=ismember(tri(masktri,[2 3]),edg,'rows');  l23=slm.resl(loc,:);
0103 %     a=max(4*l12.*l13-(l12+l13-l23).^2,0);
0104 %     r2=mean(sqrt(a),2)/4;
0105 %     lkc(2,3)=sum(mean(sqrt(l12)+sqrt(l13)+sqrt(l23),2))/2;
0106 %     lkc(3,3)=sum(r2);
0107 % end
0108 %
0109 % if isfield(slm,'lat')
0110 %     %% LKC of tetrahedra
0111 %     masktet=all(mask(tet),2);
0112 %     trimasktri=tri(masktri,:);
0113 %     lkc(1,4)=sum(masktet);
0114 %     if isfield(slm,'resl')
0115 %         [tf,loc]=ismember(tet(masktet,[1 2 3]),trimasktri,'rows');  a4=a(loc,:);
0116 %         [tf,loc]=ismember(tet(masktet,[1 2 4]),trimasktri,'rows');  a3=a(loc,:);
0117 %         [tf,loc]=ismember(tet(masktet,[1 3 4]),trimasktri,'rows');  a2=a(loc,:);
0118 %         [tf,loc]=ismember(tet(masktet,[2 3 4]),trimasktri,'rows');  a1=a(loc,:);
0119 %
0120 %         [tf,loc]=ismember(tet(masktet,[1 2]),edg,'rows');  l12=slm.resl(loc,:);
0121 %         [tf,loc]=ismember(tet(masktet,[1 3]),edg,'rows');  l13=slm.resl(loc,:);
0122 %         [tf,loc]=ismember(tet(masktet,[2 3]),edg,'rows');  l23=slm.resl(loc,:);
0123 %         [tf,loc]=ismember(tet(masktet,[1 4]),edg,'rows');  l14=slm.resl(loc,:);
0124 %         [tf,loc]=ismember(tet(masktet,[2 4]),edg,'rows');  l24=slm.resl(loc,:);
0125 %         [tf,loc]=ismember(tet(masktet,[3 4]),edg,'rows');  l34=slm.resl(loc,:);
0126 %
0127 %         d12=4*l12.*l34-(l13+l24-l23-l14).^2;
0128 %         d13=4*l13.*l24-(l12+l34-l23-l14).^2;
0129 %         d14=4*l14.*l23-(l12+l34-l24-l13).^2;
0130 %
0131 %         h=(a1<=0)|(a2<=0);
0132 %         delta12=sum(mean(sqrt(l34).*pacos((d12-a1-a2)./sqrt(a1.*a2+h)/2.*(1-h)+h),2));
0133 %         h=(a1<=0)|(a3<=0);
0134 %         delta13=sum(mean(sqrt(l24).*pacos((d13-a1-a3)./sqrt(a1.*a3+h)/2.*(1-h)+h),2));
0135 %         h=(a1<=0)|(a4<=0);
0136 %         delta14=sum(mean(sqrt(l23).*pacos((d14-a1-a4)./sqrt(a1.*a4+h)/2.*(1-h)+h),2));
0137 %         h=(a2<=0)|(a3<=0);
0138 %         delta23=sum(mean(sqrt(l14).*pacos((d14-a2-a3)./sqrt(a2.*a3+h)/2.*(1-h)+h),2));
0139 %         h=(a2<=0)|(a4<=0);
0140 %         delta24=sum(mean(sqrt(l13).*pacos((d13-a2-a4)./sqrt(a2.*a4+h)/2.*(1-h)+h),2));
0141 %         h=(a3<=0)|(a4<=0);
0142 %         delta34=sum(mean(sqrt(l12).*pacos((d12-a3-a4)./sqrt(a3.*a4+h)/2.*(1-h)+h),2));
0143 %
0144 %         r3=mean(sqrt(max((4*a1.*a2-(a1+a2-d12).^2)./(l34+(l34<=0)).*(l34>0),0)),2)/48;
0145 %
0146 %         lkc(2,4)=(delta12+delta13+delta14+delta23+delta24+delta34)/(2*pi);
0147 %         lkc(3,4)=sum(mean(sqrt(a1)+sqrt(a2)+sqrt(a3)+sqrt(a4),2))/8;
0148 %         lkc(4,4)=sum(r3);
0149 %     end
0150 % end
0151 
0152 if isfield(slm,'tri')
0153     tri=sort(slm.tri,2);
0154     edg=unique([tri(:,[1 2]); tri(:,[1 3]); tri(:,[2 3])],'rows');
0155 
0156     if nargin<2
0157         v=max(edg(:));
0158         mask=logical(zeros(1,v));
0159         mask(edg)=1;
0160     else
0161         if size(mask,1)>1
0162             mask=mask';
0163         end
0164         v=size(mask,2);
0165     end
0166         
0167     m=sum(mask(:));
0168     lkc(1,1)=m;
0169     %% LKC of edges
0170     maskedg=all(mask(edg),2);
0171     lkc(1,2)=sum(maskedg);
0172     if isfield(slm,'resl')
0173         r1=mean(sqrt(slm.resl(maskedg,:)),2);
0174         lkc(2,2)=sum(r1);
0175     end
0176 
0177     %% LKC of triangles
0178     masktri=all(mask(tri),2);
0179     lkc(1,3)=sum(masktri);
0180     if isfield(slm,'resl')
0181         [tf,loc]=ismember(tri(masktri,[1 2]),edg,'rows');  l12=slm.resl(loc,:);
0182         [tf,loc]=ismember(tri(masktri,[1 3]),edg,'rows');  l13=slm.resl(loc,:);
0183         [tf,loc]=ismember(tri(masktri,[2 3]),edg,'rows');  l23=slm.resl(loc,:);
0184         a=max(4*l12.*l13-(l12+l13-l23).^2,0);
0185         r2=mean(sqrt(a),2)/4;
0186         lkc(2,3)=sum(mean(sqrt(l12)+sqrt(l13)+sqrt(l23),2))/2;
0187         lkc(3,3)=sum(r2);
0188     end
0189     if nargout>=2
0190         reselspvert=zeros(v,1);
0191         for j=1:3
0192             reselspvert=reselspvert+accumarray(tri(masktri,j),r2,[v 1]);
0193         end
0194         D=2;
0195         reselspvert=reselspvert'/(D+1)/sqrt(4*log(2))^D;
0196     end
0197 end
0198 %%
0199 if isfield(slm,'lat')
0200     edg=SurfStatEdg(slm);
0201     % The lattice is filled with 5 alternating tetrahedra per cube
0202     [I,J,K]=size(slm.lat);
0203     IJ=I*J;
0204     [i,j]=ndgrid(1:int32(I),1:int32(J));
0205     i=i(:);
0206     j=j(:);
0207 
0208     c1=int32(find(rem(i+j,2)==0 & i<I & j<J));
0209     c2=int32(find(rem(i+j,2)==0 & i>1 & j<J));
0210     c11=int32(find(rem(i+j,2)==0 & i==I & j<J));
0211     c21=int32(find(rem(i+j,2)==0 & i==I & j>1));
0212     c12=int32(find(rem(i+j,2)==0 & i<I & j==J));
0213     c22=int32(find(rem(i+j,2)==0 & i>1 & j==J));
0214 
0215     d1=find(rem(i+j,2)==1 & i<I & j<J)+IJ;
0216     d2=find(rem(i+j,2)==1 & i>1 & j<J)+IJ;
0217 
0218     tri1=[c1     c1+1    c1+1+I;  % bottom slice
0219         c1     c1+I    c1+1+I;
0220         c2-1   c2      c2-1+I;
0221         c2     c2-1+I  c2+I];
0222     tri2=[c1     c1+1    c1+1+IJ;  % between slices
0223         c1     c1+IJ   c1+1+IJ;
0224         c1     c1+I    c1+I+IJ;
0225         c1     c1+IJ   c1+I+IJ;
0226         c1     c1+1+I  c1+1+IJ;
0227         c1     c1+1+I  c1+I+IJ;
0228         c1     c1+1+IJ c1+I+IJ;
0229         c1+1+I c1+1+IJ c1+I+IJ;
0230         c2-1   c2      c2-1+IJ;
0231         c2     c2-1+IJ c2+IJ;
0232         c2-1   c2-1+I  c2-1+IJ;
0233         c2-1+I c2-1+IJ c2-1+I+IJ;
0234         c2     c2-1+I  c2+I+IJ;
0235         c2     c2-1+IJ c2+I+IJ;
0236         c2     c2-1+I  c2-1+IJ;
0237         c2-1+I c2-1+IJ c2+I+IJ;
0238         c11    c11+I    c11+I+IJ;
0239         c11    c11+IJ   c11+I+IJ;
0240         c21-I  c21      c21-I+IJ;
0241         c21    c21-I+IJ c21+IJ;
0242         c12    c12+1    c12+1+IJ;
0243         c12    c12+IJ   c12+1+IJ;
0244         c22-1  c22      c22-1+IJ;
0245         c22    c22-1+IJ c22+IJ];
0246     tri3=[d1     d1+1    d1+1+I;  % top slice
0247         d1     d1+I    d1+1+I;
0248         d2-1   d2      d2-1+I;
0249         d2     d2-1+I  d2+I];
0250     tet1=[c1     c1+1    c1+1+I    c1+1+IJ; % between slices
0251         c1     c1+I    c1+1+I    c1+I+IJ;
0252         c1     c1+1+I  c1+1+IJ   c1+I+IJ;
0253         c1     c1+IJ   c1+1+IJ   c1+I+IJ;
0254         c1+1+I c1+1+IJ c1+I+IJ   c1+1+I+IJ;
0255         c2-1   c2      c2-1+I    c2-1+IJ;
0256         c2     c2-1+I  c2+I      c2+I+IJ;
0257         c2     c2-1+I  c2-1+IJ   c2+I+IJ;
0258         c2     c2-1+IJ c2+IJ     c2+I+IJ;
0259         c2-1+I c2-1+IJ c2-1+I+IJ c2+I+IJ];
0260 
0261     v=sum(slm.lat(:));
0262     if nargin<2
0263         mask=logical(ones(1,v));
0264     else
0265         if size(mask,1)>1
0266             mask=mask';
0267         end
0268     end
0269     if nargout>=2
0270         reselspvert=zeros(v,1);
0271     end
0272 
0273     vs=cumsum(squeeze(sum(sum(slm.lat))));
0274     vs=int32([0; vs; vs(K)]);
0275     es=0;
0276     lat=logical(zeros(I,J,2));
0277     lat(:,:,1)=slm.lat(:,:,1);
0278     lkc=zeros(4);
0279     fprintf(1,'%s',[num2str(K) ' slices to resel, % remaining: 100 ']);
0280     n10=floor(K/10);
0281     for k=1:K
0282         if rem(k,n10)==0
0283             fprintf(1,'%s',[num2str(round(100*(1-k/K))) ' ']);
0284         end
0285         f=rem(k,2);
0286         if k<K
0287             lat(:,:,f+1)=slm.lat(:,:,k+1);
0288         else
0289             lat(:,:,f+1)=zeros(I,J);
0290         end
0291         vid=int32(cumsum(lat(:)).*lat(:))';
0292         if f
0293             edg1=edg(edg(:,1)>vs(k) & edg(:,1)<=vs(k+1),:)-vs(k);
0294             edg2=edg(edg(:,1)>vs(k) & edg(:,2)<=vs(k+2),:)-vs(k);
0295             tri=[vid(tri1(all(lat(tri1),2),:)); vid(tri2(all(lat(tri2),2),:))];
0296             mask1=mask((vs(k)+1):vs(k+2));
0297         else
0298             edg1=[edg(edg(:,1)>vs(k) & edg(:,2)<=vs(k+1),:)-vs(k)+vs(k+2)-vs(k+1);
0299                 edg(edg(:,1)<=vs(k+1) & edg(:,2)>vs(k+1),2)-vs(k+1) ...
0300                 edg(edg(:,1)<=vs(k+1) & edg(:,2)>vs(k+1),1)-vs(k)+vs(k+2)-vs(k+1)];
0301             edg2=[edg1; edg(edg(:,1)>vs(k+1) & edg(:,2)<=vs(k+2),:)-vs(k+1)];
0302             tri=[vid(tri3(all(lat(tri3),2),:)); vid(tri2(all(lat(tri2),2),:))];
0303             mask1=[mask((vs(k+1)+1):vs(k+2)) mask((vs(k)+1):vs(k+1))];
0304         end
0305         tet=vid(tet1(all(lat(tet1),2),:));
0306 
0307         m1=max(double(edg2(:,1)));
0308         ue=double(edg2(:,1))+m1*(double(edg2(:,2))-1);
0309         e=size(edg2,1);
0310         ae=1:e;
0311         if e<2^31
0312             sparsedg=sparse(ue,1,ae);
0313         end
0314         %%
0315         lkc1=zeros(4);
0316         lkc1(1,1)=sum(mask((vs(k)+1):vs(k+1)));
0317         
0318         %% LKC of edges
0319         maskedg=all(mask1(edg1),2);
0320         lkc1(1,2)=sum(maskedg);
0321         if isfield(slm,'resl')
0322             r1=mean(sqrt(slm.resl(find(maskedg)+es,:)),2);
0323             lkc1(2,2)=sum(r1);
0324         end
0325         
0326         %% LKC of triangles
0327         masktri=all(mask1(tri),2);
0328         lkc1(1,3)=sum(masktri);
0329         if isfield(slm,'resl')
0330             if e<2^31
0331                 l12=slm.resl(sparsedg(tri(masktri,1)+m1*(tri(masktri,2)-1),1)+es,:);
0332                 l13=slm.resl(sparsedg(tri(masktri,1)+m1*(tri(masktri,3)-1),1)+es,:);
0333                 l23=slm.resl(sparsedg(tri(masktri,2)+m1*(tri(masktri,3)-1),1)+es,:);
0334             else
0335                 l12=slm.resl(interp1(ue,ae,tri(masktri,1)+m1*(tri(masktri,2)-1),'nearest')+es,:);
0336                 l13=slm.resl(interp1(ue,ae,tri(masktri,1)+m1*(tri(masktri,3)-1),'nearest')+es,:);
0337                 l23=slm.resl(interp1(ue,ae,tri(masktri,2)+m1*(tri(masktri,3)-1),'nearest')+es,:);
0338             end
0339             a=max(4*l12.*l13-(l12+l13-l23).^2,0);
0340             r2=mean(sqrt(a),2)/4;
0341             lkc1(2,3)=sum(mean(sqrt(l12)+sqrt(l13)+sqrt(l23),2))/2;
0342             lkc1(3,3)=sum(r2);
0343             
0344             if nargout>=2 & K==1
0345                 for j=1:3
0346                     if f
0347                         v1=tri(masktri,j)+vs(k);
0348                     else
0349                         v1=tri(masktri,j)+vs(k+1);
0350                         v1=v1-int32(v1>vs(k+2))*(vs(k+2)-vs(k));
0351                     end
0352                     reselspvert=reselspvert+accumarray(v1,r2,[v 1]);
0353                 end
0354             end
0355         end
0356         
0357         %% LKC of tetrahedra
0358         masktet=all(mask1(tet),2);
0359         lkc1(1,4)=sum(masktet);
0360         if isfield(slm,'resl') & k<K
0361             if e<2^31
0362                 l12=slm.resl(sparsedg(tet(masktet,1)+m1*(tet(masktet,2)-1),1)+es,:);
0363                 l13=slm.resl(sparsedg(tet(masktet,1)+m1*(tet(masktet,3)-1),1)+es,:);
0364                 l23=slm.resl(sparsedg(tet(masktet,2)+m1*(tet(masktet,3)-1),1)+es,:);
0365                 l14=slm.resl(sparsedg(tet(masktet,1)+m1*(tet(masktet,4)-1),1)+es,:);
0366                 l24=slm.resl(sparsedg(tet(masktet,2)+m1*(tet(masktet,4)-1),1)+es,:);
0367                 l34=slm.resl(sparsedg(tet(masktet,3)+m1*(tet(masktet,4)-1),1)+es,:);
0368             else
0369                 l12=slm.resl(interp1(ue,ae,tet(masktet,1)+m1*(tet(masktet,2)-1),'nearest')+es,:);
0370                 l13=slm.resl(interp1(ue,ae,tet(masktet,1)+m1*(tet(masktet,3)-1),'nearest')+es,:);
0371                 l23=slm.resl(interp1(ue,ae,tet(masktet,2)+m1*(tet(masktet,3)-1),'nearest')+es,:);
0372                 l14=slm.resl(interp1(ue,ae,tet(masktet,1)+m1*(tet(masktet,4)-1),'nearest')+es,:);
0373                 l24=slm.resl(interp1(ue,ae,tet(masktet,2)+m1*(tet(masktet,4)-1),'nearest')+es,:);
0374                 l34=slm.resl(interp1(ue,ae,tet(masktet,3)+m1*(tet(masktet,4)-1),'nearest')+es,:);
0375             end
0376             a4=max(4*l12.*l13-(l12+l13-l23).^2,0);
0377             a3=max(4*l12.*l14-(l12+l14-l24).^2,0);
0378             a2=max(4*l13.*l14-(l13+l14-l34).^2,0);
0379             a1=max(4*l23.*l24-(l23+l24-l34).^2,0);
0380 
0381             d12=4*l12.*l34-(l13+l24-l23-l14).^2;
0382             d13=4*l13.*l24-(l12+l34-l23-l14).^2;
0383             d14=4*l14.*l23-(l12+l34-l24-l13).^2;
0384 
0385             h=(a1<=0)|(a2<=0);
0386             delta12=sum(mean(sqrt(l34).*pacos((d12-a1-a2)./sqrt(a1.*a2+h)/2.*(1-h)+h),2));
0387             h=(a1<=0)|(a3<=0);
0388             delta13=sum(mean(sqrt(l24).*pacos((d13-a1-a3)./sqrt(a1.*a3+h)/2.*(1-h)+h),2));
0389             h=(a1<=0)|(a4<=0);
0390             delta14=sum(mean(sqrt(l23).*pacos((d14-a1-a4)./sqrt(a1.*a4+h)/2.*(1-h)+h),2));
0391             h=(a2<=0)|(a3<=0);
0392             delta23=sum(mean(sqrt(l14).*pacos((d14-a2-a3)./sqrt(a2.*a3+h)/2.*(1-h)+h),2));
0393             h=(a2<=0)|(a4<=0);
0394             delta24=sum(mean(sqrt(l13).*pacos((d13-a2-a4)./sqrt(a2.*a4+h)/2.*(1-h)+h),2));
0395             h=(a3<=0)|(a4<=0);
0396             delta34=sum(mean(sqrt(l12).*pacos((d12-a3-a4)./sqrt(a3.*a4+h)/2.*(1-h)+h),2));
0397 
0398             r3=mean(sqrt(max((4*a1.*a2-(a1+a2-d12).^2)./(l34+(l34<=0)).*(l34>0),0)),2)/48;
0399 
0400             lkc1(2,4)=(delta12+delta13+delta14+delta23+delta24+delta34)/(2*pi);
0401             lkc1(3,4)=sum(mean(sqrt(a1)+sqrt(a2)+sqrt(a3)+sqrt(a4),2))/8;
0402             lkc1(4,4)=sum(r3);
0403             
0404             if nargout>=2
0405                 for j=1:4
0406                     if f
0407                         v1=tet(masktet,j)+vs(k);
0408                     else
0409                         v1=tet(masktet,j)+vs(k+1);
0410                         v1=v1-int32(v1>vs(k+2))*(vs(k+2)-vs(k));
0411                     end
0412                     reselspvert=reselspvert+accumarray(v1,r3,[v 1]);
0413                 end
0414             end
0415         end
0416         lkc=lkc+lkc1;
0417         es=es+size(edg1,1);
0418     end
0419     if nargout>=2
0420         D=2+(K>1);
0421         reselspvert=reselspvert'/(D+1)/sqrt(4*log(2))^D;
0422     end
0423     fprintf(1,'%s\n','Done');
0424 end
0425 %% resels
0426 D=size(lkc,1)-1;
0427 D2=size(lkc,2)-1;
0428 tpltz=toeplitz((-1).^(0:D),(-1).^(0:D2));
0429 lkcs=sum(tpltz.*lkc,2)';
0430 lkcs=lkcs(1:max(find(abs(lkcs))));
0431 resels=lkcs./sqrt(4*log(2)).^(0:D);
0432 
0433 return
0434 end
0435 
0436 %%
0437 function y=pacos(x)
0438 y=acos(min(abs(x),1).*sign(x));
0439 return
0440 end