Home > bioelectromagnetism > PVC_eqW.m

PVC_eqW

PURPOSE ^

Partial Volume Correction

SYNOPSIS ^

function [pvc]=PVC_eqW(b,R,PSFx,PSFy,PSFz)

DESCRIPTION ^

 Partial Volume Correction 
 
 Usage : [pvc]=PVC(b,R,PSFx,PSFy,PSFz)

 Inputs
 b    = Measured PET Data (3d matrix) b(1:xdim,1:ydim,1:zdim) (b must be in counts/voxel)
 R    = ROI Definitions (4d binary matrix) R(1:xdim,1:ydim,1:zdim,1:no_regions)
 PSFx = Matrix Point Spread Function in x Direction. The ith row is the kernel for the ith voxel.
 PSFy = Matrix Point Spread Function in y Direction. The ith row is the kernel for the ith voxel.
 PSFz = Matrix Point Spread Function in z Direction. The ith row is the kernel for the ith voxel.

 Outputs
 These are contained in the structure pvc
 pvc.xm is the mean parameter estimate for the true activity concentrations in the ROIs
 pvc.xv is contains the associated variance estimates  
 pvc.ROIsize is size of each region.

 Tensor Implementation of the weighted least squares Partial Volume Correction Problem.
 PSF values are for reconstructed voxel size of b. Each number represents the recovery 
 coefficient in the next voxel for a point source in position 1 of the vector. 

 This method applies equal weighting to each voxel
 
 (c) Roger Gunn & John Aston, 26-02-2000

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [pvc]=PVC_eqW(b,R,PSFx,PSFy,PSFz)
0002 
0003 % Partial Volume Correction
0004 %
0005 % Usage : [pvc]=PVC(b,R,PSFx,PSFy,PSFz)
0006 %
0007 % Inputs
0008 % b    = Measured PET Data (3d matrix) b(1:xdim,1:ydim,1:zdim) (b must be in counts/voxel)
0009 % R    = ROI Definitions (4d binary matrix) R(1:xdim,1:ydim,1:zdim,1:no_regions)
0010 % PSFx = Matrix Point Spread Function in x Direction. The ith row is the kernel for the ith voxel.
0011 % PSFy = Matrix Point Spread Function in y Direction. The ith row is the kernel for the ith voxel.
0012 % PSFz = Matrix Point Spread Function in z Direction. The ith row is the kernel for the ith voxel.
0013 %
0014 % Outputs
0015 % These are contained in the structure pvc
0016 % pvc.xm is the mean parameter estimate for the true activity concentrations in the ROIs
0017 % pvc.xv is contains the associated variance estimates
0018 % pvc.ROIsize is size of each region.
0019 %
0020 % Tensor Implementation of the weighted least squares Partial Volume Correction Problem.
0021 % PSF values are for reconstructed voxel size of b. Each number represents the recovery
0022 % coefficient in the next voxel for a point source in position 1 of the vector.
0023 %
0024 % This method applies equal weighting to each voxel
0025 %
0026 % (c) Roger Gunn & John Aston, 26-02-2000
0027 
0028 
0029 pvc.xm=[];pvc.xv=[];
0030 if nargin~=5 disp('Usage : [xm,xv]=PVC(b,R,PSFx,PSFy,PSFz)');return;end
0031 if length(size(R))~=4;disp('Error: ROI matrix should be 4D');return;end
0032 if length(size(b))~=3;disp('Error: b matrix should be 3D');return;end
0033 if size(R,4)<2;disp('Error: There must be at least 1 ROI');return;end
0034 if sum(size(R)==[size(b) 0])~=3;disp('Error: ROI and Image volume have different dimensions');return;end
0035 if sum(size(PSFx)~=size(b,1));disp('Error: Check PSFx');return;end
0036 if sum(size(PSFy)~=size(b,2));disp('Error: Check PSFy');return;end
0037 if sum(size(PSFz)~=size(b,3));disp('Error: Check PSFz');return;end
0038 if (size(b,1)==1|size(b,2)==1|size(b,3)==1);disp('Try a 3D Volume, It''s much nicer !');return;end
0039 
0040 [xdim ydim zdim]     = size(b);
0041 no_regions             = size(R,4);
0042 PSFz                     = PSFz';
0043 
0044 
0045 % Blur ROI's
0046 for r=1:no_regions
0047    fprintf('.');
0048     for x=1:size(b,3);R(:,:,x,r)=PSFx(1:xdim,1:xdim)*squeeze(R(:,:,x,r));end
0049     for y=1:size(b,1);R(y,:,:,r)=PSFy(1:ydim,1:ydim)*squeeze(R(y,:,:,r));end
0050    for z=1:size(b,2);R(:,z,:,r)=squeeze(R(:,z,:,r))*PSFz(1:zdim,1:zdim);end
0051 end
0052 
0053 % xm = inv(R'P'PR)R'Pb
0054 PR         = reshape(R,xdim*ydim*zdim,no_regions);
0055 pvc.xm    = inv(PR'*PR)*PR'*reshape(b,xdim*ydim*zdim,1);
0056 pvc.xv    = inv(PR'*PR)*sum((reshape(b,xdim*ydim*zdim,1)-PR*pvc.xm).^2)/(xdim*ydim*zdim-no_regions);
0057 
0058 pvc.xm=pvc.xm';
0059 return;

Generated on Mon 15-Aug-2005 15:36:19 by m2html © 2003