ComputeResels2

 Compute FWHM of the 3D Autocorrelation Function
 and use to Estimate the Number of Resolvable Elements
 in the Image.

0001 function [Sx, Sy, Sz, U, Resels]=ComputeResels2(X, Mask, XDIM, YDIM, ZDIM)
0002 % Compute FWHM of the 3D Autocorrelation Function
0003 % and use to Estimate the Number of Resolvable Elements
0004 % in the Image.
0005 
0006 DIMS=[XDIM YDIM ZDIM];
0007 
0008 %fid=fopen(ImageName,'r');
0009 %  X=fread(fid,'float32');
0010 %fclose(fid);
0011 
0012 %fid=fopen('Mask.img','r');
0013 %  Mask=fread(fid);
0014 %fclose(fid);
0015 M=sum(Mask);
0016 
0017 Z=reshape(X,XDIM,YDIM,ZDIM);
0018 
0019 % Determine 3D Autocorrelation Function of the Image
0020 z=fftn(Z);
0021 w=z.*conj(z);
0022 W=real(ifftn(w));
0023 W=W./W(1,1,1);  % So that the max is unity
0024 W=fftshift(W);  % So that the DC offset is a the center of the volume
0025 
0026 S=zeros(3,3);
0027 SumW=0;
0028 
0029 % Measure Upper Triangle Cross-Products of the 3D Autocorrelation Function
0030 for x=1:DIMS(1),
0031     for y=1:DIMS(2),
0032         for z=1:DIMS(3),
0033             S(1,1)=S(1,1)+abs(W(x,y,z)).*(x-(DIMS(1)+1)./2).^2;
0034 end
0035 end
0036 end
0037 
0038 for x=1:DIMS(1),
0039     for y=1:DIMS(2),
0040         for z=1:DIMS(3),
0041             S(2,2)=S(2,2)+abs(W(x,y,z)).*(y-(DIMS(2)+1)./2).^2;
0042 end
0043 end
0044 end
0045 
0046 for x=1:DIMS(1),
0047     for y=1:DIMS(2),
0048         for z=1:DIMS(3),
0049             S(3,3)=S(3,3)+abs(W(x,y,z)).*(z-(DIMS(2)+1)./2).^2;
0050 end
0051 end
0052 end
0053 
0054 for x=1:DIMS(1),
0055     for y=1:DIMS(2),
0056         for z=1:DIMS(3),
0057             S(1,2)=S(1,2)+abs(W(x,y,z)).*(x-(DIMS(1)+1)./2).*(y-(DIMS(2)+1)./2);
0058 end
0059 end
0060 end
0061 
0062 for x=1:DIMS(1),
0063     for y=1:DIMS(2),
0064         for z=1:DIMS(3),
0065             S(1,3)=S(1,3)+abs(W(x,y,z)).*(x-(DIMS(1)+1)./2).*(z-(DIMS(3)+1)./2);
0066 end
0067 end
0068 end
0069 
0070 for x=1:DIMS(1),
0071     for y=1:DIMS(2),
0072         for z=1:DIMS(3),
0073             S(2,3)=S(2,3)+abs(W(x,y,z)).*(y-(DIMS(2)+1)./2).*(z-(DIMS(3)+1)./2);
0074 end
0075 end
0076 end
0077 
0078 % Compute Sum of Absolute Autocorrelation Values
0079 SumW=sum(sum(sum(abs(W))));
0080 
0081 % Reflect and Scale the Autocorrelation Cross-Product Matrix by 2*log(2)./SumW
0082 %S=2*log(2)*(S+S'.*(1-eye(3,3)))./SumW;
0083 S=(S+S'.*(1-eye(3,3)))./SumW;
0084 
0085 % Obtain SVD Decomposition of S
0086 [U V P]=svd(S);
0087 
0088 % Estimated FWHM in Each Axis
0089 Sx=sqrt(V(1,1));
0090 Sy=sqrt(V(2,2));
0091 Sz=sqrt(V(3,3));
0092 
0093 % Compute the Area of an Ellipsoid With These Major Axes
0094 A=(S(1,1)*S(2,2)*S(3,3))^(-1/6);
0095 B=(4/3)*pi*sqrt(det(S./2));
0096 
0097 % Contour Method: Find ISO contour at W==0.5 in X and Y
0098 %C=contourc(W(:,:,(ZDIM+1)./2),[0.5 0.5]);
0099 %F(1,:)=C(1,:)-(XDIM+1)./2;
0100 %F(2,:)=C(2,:)-(YDIM+1)./2;
0101 %b=F(:,2:end)*F(:,2:end)';
0102 %Sx=sqrt(b(1,1));
0103 %Sy=sqrt(b(2,2));
0104 
0105 % Find ISO contour at W==0.5 in Z
0106 %clear C F;
0107 %a=permute(W,[3 1 2]);
0108 %C=contourc(a(:,:,ceil((YDIM+1)./2)),[0.5 0.5]);
0109 %F(1,:)=C(1,:)-(ZDIM+1)./2;
0110 %b=F(:,2:end)*F(:,2:end)';
0111 %Sz=sqrt(b(1,1));
0112 
0113 % Compute Resolveable Elements
0114 %Resels=M.*A;
0115 Resels=M.*B^(-1/3);
0116 
0117 % Clear Temp Variables
0118 clear X Z z w Mask M V a F;
0119 
0120 return;
0121 
0122 
0123 
0124 
0125 
0126 
0127 
0128 
0129 
0130 
0131 
0132 
0133 
0134 
0135 
0136 
0137 
0138 
0139 
0140 
0141 
0142 
0143 
0144 
0145 
0146 
0147 
0148 
0149 
0150 
0151 
0152 
0153 
0154 
0155 
0156 
0157 
0158 
0159