avw_shrinkwrap.m

mri_toolbox是一个工具用来MRI. 来自于SourceForge, 我上传这个软件,希望能结识对医疗软件感兴趣的兄弟.

function [FV, Edges] = avw_shrinkwrap(avw,FV,smooth,vthresh,interpVal,...
    fitval,fittol,fititer,fitchange,fitvattr)

% avw_shrinkwrap - Tesselate the surface of a 3D Analyze 7.5 avw struct
% 
% [FV, Edges] = avw_shrinkwrap(avw,FV,smooth,vthresh,interpVal,...
%               fitval,fittol,fititer,fitchange,fitvattr)
% 
% avw       - an Analyze 7.5 data struct (see avw_read)
% FV        - input tesselation; if empty, sphere tesselation 
%             is created.  FV has fields FV.vertices, FV.faces
% smooth    - Gaussian smoothing of avw (5x5x5 kernel), 0|1 (default 0)
% vthresh   - binarise avw at fitval, 0|1 (default 0)
% interpVal - use radial interpolation (faster) or incremental
%             radial shrink (slower), 0|1 (default 1, faster)
% fitval    - image intensity to shrink wrap (default 20)
% fittol    - image intensity tolerance (default 5)
% fititer   - max number of iterations to fit (default 200)
% fitchange - least sig. change in intensity values 
%             between fit iterations (default 2)
% fitvattr  - vertex attraction (constraint), 0:1, smaller
%             values are less constraint; close to 0 for 
%             no constraint is useful when dealing with 
%             binary volumes, otherwise 0.4 (40%) seems OK
% 
% FV        - a struct with 2562 vertices and 5120 faces
% Edges     - a [2562,2562] matrix of edge connectivity for FV
% 
% This function has been developed to find the scalp surface 
% for MRI of the human head.
% 
% It starts with a sphere tesselation (large radius) and moves
% each vertex point toward the center of the volume until it
% lies at or near the fitval.  The vertices are constrained to
% move only along the radial projection from the origin and they
% are also required to stay within a small distance of their
% neighbours.  The function is not optimised for speed, but 
% it should produce reasonable results.
% 
% Example of creating a scalp tesselation for SPM T1 MRI template:
% 
%   avw = avw_read('T1');
%   FV = avw_shrinkwrap(avw,[],0,0,[],intensity,5.0,50,0.5,0.4);
%   patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[.6 .6 .6]);
% 
% The output vertex coordinates are in meters with an origin at (0,0,0), 
% which lies at the center of the MRI volume.  This function uses the
% avw.hdr.dime.pixdim values to convert from voxel coordinates into
% meters.
% 
% See also: ISOSURFACE, SPHERE_TRI, MESH_REFINE_TRI4,
%           MESH_BEM_SHELLS_FUNC, MESH_BEM_SHELLS_SCRIPT
% 

% $Revision: 1.2 $ $Date: 2004/02/07 01:41:51 $

% Licence:  GNU GPL, no implied or express warranties
% History:  08/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - adapted from mesh_shrinkwrap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parse arguments

if ~exist('avw','var'), error('...no input volume\n');
elseif isempty(avw),    error('...empty input volume\n');
end

% Find approximate center of volume
center = avw_center(avw);
origin = center.abs.voxels;

version = '[$Revision: 1.2 $]';
fprintf('AVW_SHRINKWRAP [v%s]\n',version(12:16));  tic;

if ~exist('smooth','var'),   smooth = 0;
elseif isempty(smooth),      smooth = 0;
end

if ~exist('vthresh','var'),  vthresh = 0;
elseif isempty(vthresh),     vthresh = 0;
end

if ~exist('interpVal','var'), interpVal = 0;
elseif isempty(interpVal),    interpVal = 0;
end

if ~exist('fitval','var'),   fit.val = 20;
elseif isempty(fitval),      fit.val = 20;
else                         fit.val = fitval;
end

if ~exist('fittol','var'),   fit.tol = 5;
elseif isempty(fittol),      fit.tol = 5;
else                         fit.tol = fittol;
end

if fit.val <= fit.tol,
    error('...must use fit tolerance < fit value\n');
end

if ~exist('fititer','var'),  fit.iter = 200;
elseif isempty(fititer),     fit.iter = 200;
else                         fit.iter = fititer;
end

if ~exist('fitchange','var'),fit.change = 2;
elseif isempty(fitchange),   fit.change = 2;
else                         fit.change = fitchange;
end

if ~exist('fitvattr','var'), fit.vattr = 0.4;
elseif isempty(fitvattr),    fit.vattr = 0.4;
else                         fit.vattr = fitvattr;
end
if fit.vattr > 1,
    fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 1\n');
    fit.vattr = 1;
end
if fit.vattr < 0,
    fprintf('...fit vertattr (v) must be 0 <= v <= 1, setting v = 0.\n');
    fit.vattr = 0;
end


% get size of volume, in voxels
xdim = size(avw.img,1);
ydim = size(avw.img,2);
zdim = size(avw.img,3);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check whether to create a sphere tesselation
% or use an input tesselation as the start point

sphere = 0;
if ~exist('FV','var'),
    sphere = 1;
elseif ~isfield(FV,'vertices'),
    sphere = 1;
elseif ~isfield(FV,'faces'),
    sphere = 1;
elseif isempty(FV.vertices),
    sphere = 1;
elseif isempty(FV.faces),
    sphere = 1;
end
if sphere,
    % Create a sphere tesselation to encompass the volume
    diameter = max([xdim ydim zdim]);
    radius = diameter / 1.5;
    FV = sphere_tri('ico',4,radius); % 2562 vertices
    % Shift the center of the sphere to the center of the volume
    FV.vertices = FV.vertices + repmat(origin,size(FV.vertices,1),1);
else
    fprintf('...using input FV tesselation...\n');
end

% the 'edge' matrix is the connectivity of all vertices,
% used to find neighbours during movement of vertices,
% including smoothing the tesselation
FV.edge = mesh_edges(FV);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% theshold the volume based on background noise estimate
slice_x0 = avw.img(1,:,:);
slice_x1 = avw.img(end,:,:);
slice_y0 = avw.img(:,1,:);
slice_y1 = avw.img(:,end,:);
slice_z0 = avw.img(:,:,1);
slice_z1 = avw.img(:,:,end);

slice_x0_mean = mean(mean(slice_x0));
slice_x1_mean = mean(mean(slice_x1));
slice_y0_mean = mean(mean(slice_y0));
slice_y1_mean = mean(mean(slice_y1));
slice_z0_mean = mean(mean(slice_z0));
slice_z1_mean = mean(mean(slice_z1));

background_mean = mean([slice_x0_mean,slice_x1_mean,slice_y0_mean,slice_y1_mean,slice_z0_mean,slice_z1_mean]);

% threshold the volume by this background value
index = find(avw.img < background_mean);
if index,
  avw.img(index) = 0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate a binary version of the volume, zeroing all 
% values below a threshold, setting all others to threshold
if vthresh,
    fprintf('...binarising volume at fit.val threshold...'); tic;
    Vindex = find(avw.img < fit.val);
    avw.img(Vindex) = 0;
    Vindex = find(avw.img >= fit.val);
    avw.img(Vindex) = fit.val;
    t = toc; fprintf('done (%5.2f sec)\n',t);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% smooth the volume
if smooth,
    fprintf('...gaussian smoothing (5-10 minutes)...'); tic;
    avw.img = smooth3(avw.img,'gaussian',5,.8);
    t = toc; fprintf('done (%5.2f sec)\n',t);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate scalp intensity
fit = estimate_scalp(FV,avw.img,origin,fit);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now begin recursion
fprintf('...fitting...\n');	tic;

i = 1;
Fminima = 0;
intensityAtRMean = [0 0];

while i <= fit.iter,
    
    if interpVal,
        % use radial interpolation method, moving directly
        % to the intensity value nearest correct intensity
        [FV, intensityAtR, R] = locate_val(FV,avw.img,origin,fit);
    else
        % use incremental method, moving along radial line
        % gradually until finding correct intensity
        [FV, intensityAtR, R] = shrink_wrap(FV,avw.img,origin,fit);
    end
    
    intensityAtRMean = [ intensityAtRMean(2), mean(intensityAtR) ] ;
    
    fprintf('...distance:  mean = %8.4f voxels, std = %8.4f voxels\n',mean(R),std(R));
    fprintf('...intensity: mean = %8.4f,        std = %8.4f\n',...
        mean(intensityAtR),std(intensityAtR));
    fprintf('...real iteration: %3d\n',i);
    
    % Is the mean distance reasonable?
    if mean(R) < 0.5,
        error('...mean distance < 0.5 voxel!\n');
    end
    
    % MDifVal is the mean of the absolute difference
    % between the vertex intensity and the fit intensity
    MDifVal = abs(intensityAtRMean(2) - fit.val);
    
    % Is the mean difference within the tolerance range?
    if MDifVal < fit.tol,
        fprintf('...mean intensity difference < tolerance (%5.2f +/- %5.2f)\n',...
            fit.val,fit.tol);
        break;
    else
        fprintf('...mean intensity difference > tolerance (%5.2f +/- %5.2f)\n',...
            fit.val,fit.tol);
    end
    
    % How much has the intensity values changed?
    if (i > 1) & intensityAtRMean(2) > 0,
        if intensityAtRMean(2) - intensityAtRMean(1) < fit.change,
            fprintf('...no significant intensity change (< %5.2f) in this iteration\n',...
                fit.change);
            Fminima = Fminima + 1;
            if Fminima >= 5,
                fprintf('...no significant intensity change in last 5 iterations\n');
                break;
            end
        else
            Fminima = 0;
        end
    end
    
    % Ensure that iterations begin when MDifVal is truly initialised
    if isnan(MDifVal),
        i = 1;
    else,
        i = i + 1;
    end
    
end

FV = mesh_smooth(FV,origin,fit.vattr);

% Remove large edges matrix from FV
Edges = FV.edge;
FV = struct('vertices',FV.vertices,'faces',FV.faces);

% Now center the output vertices at 0,0,0 by subtracting
% the volume centroid
FV.vertices(:,1) = FV.vertices(:,1) - origin(1);
FV.vertices(:,2) = FV.vertices(:,2) - origin(2);
FV.vertices(:,3) = FV.vertices(:,3) - origin(3);

% Scale the vertices by the pixdim values, after
% converting them from mm to meters
xpixdim = double(avw.hdr.dime.pixdim(2)) / 1000;
ypixdim = double(avw.hdr.dime.pixdim(3)) / 1000;
zpixdim = double(avw.hdr.dime.pixdim(4)) / 1000;

FV.vertices(:,1) = FV.vertices(:,1) .* xpixdim;
FV.vertices(:,2) = FV.vertices(:,2) .* ypixdim;
FV.vertices(:,3) = FV.vertices(:,3) .* zpixdim;


t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return










%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunctions...






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fit] = estimate_scalp(FV,vol,origin,fit),

xo = origin(1); yo = origin(2); zo = origin(3);

Nvert = size(FV.vertices,1);

% Estimate the scalp intensity, using 10% of vertices
N = round(0.05 * Nvert);
scalp_intensity = zeros(N,1);

fprintf('...estimating scalp intensity from %d vertices\n', N);
fprintf('...starting scalp intensity = %d\n', fit.val);

indices = round(rand(1,N) * Nvert);

i = 0;
for v = indices,
  
  x = FV.vertices(v,1);
  y = FV.vertices(v,2);
  z = FV.vertices(v,3);
  r = sqrt( (x-xo).^2 + (y-yo).^2 + (z-zo).^2 );
  l = (x-xo)/r; % cos alpha
  m = (y-yo)/r; % cos beta
  n = (z-zo)/r; % cos gamma
  
  % find discrete points from zero to the vertex
  POINTS = 250;
  radial_distances = linspace(0,r,POINTS);
  
  L = repmat(l,1,POINTS);
  M = repmat(m,1,POINTS);
  N = repmat(n,1,POINTS);
  
  XI = (L .* radial_distances) + xo;
  YI = (M .* radial_distances) + yo;
  ZI = (N .* radial_distances) + zo;
  
  % interpolate volume values at these points
  % ( not sure why have to swap XI,YI here )
  VI = interp3(vol,YI,XI,ZI,'*nearest');
  
  % find location in VI where intensity gradient is steep
  half_max = 0.5 * max(VI);
  index_maxima = find(VI > half_max);
  % use the largest index value to locate the maxima intensity that lie
  % furthest toward the outside of the head
  index_max = index_maxima(end);
  
  i = i + 1;
  scalp_intensity(i,1) = VI(index_max);
  
  %figure; plot(radial_distances,VI)
  
end

fit.val = mean(scalp_intensity);
fit.tol = std(scalp_intensity);

fprintf('...estimated scalp intensity = %g\n', fit.val);
fprintf('...estimated tolerance intensity = %g\n', fit.tol);

return