avw_shrinkwrap.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV, intensityAtR, R] = locate_val(FV,vol,origin,fit),

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

Nvert = size(FV.vertices,1);
progress = round(.1 * Nvert);

% Initialise difference intensity & distance arrays
intensityAtR = zeros(Nvert,1);
R = intensityAtR;

% Find distance and direction cosines for line from 
% origin to all vertices
XV = FV.vertices(:,1);
YV = FV.vertices(:,2);
ZV = FV.vertices(:,3);
RV = sqrt( (XV-xo).^2 + (YV-yo).^2 + (ZV-zo).^2 );
LV = (XV-xo)./RV; % cos alpha
MV = (YV-yo)./RV; % cos beta
NV = (ZV-zo)./RV; % cos gamma

% Check for binary volume data, if empty, binary
binvol = find(vol > 1);

% Locate each vertex at a given fit value
tic
for v = 1:Nvert,
    
    if v > progress,
        fprintf('...interp3 processed %4d of %4d vertices',progress,Nvert);
        t = toc; fprintf(' (%5.2f sec)\n',t);
        progress = progress + progress;
    end
    
    % Find direction cosines for line from origin to vertex
    x = XV(v);
    y = YV(v);
    z = ZV(v);
    d = RV(v);
    l = LV(v); % cos alpha
    m = MV(v); % cos beta
    n = NV(v); % cos gamma
    
    % find discrete points between origin
    % and vertex + 20% of vertex distance
    POINTS = 250;
    
    Rarray = linspace(0,(d + .2 * d),POINTS);
    
    L = repmat(l,1,POINTS);
    M = repmat(m,1,POINTS);
    N = repmat(n,1,POINTS);
    
    XI = (L .* Rarray) + xo;
    YI = (M .* Rarray) + yo;
    ZI = (N .* Rarray) + zo;
    
    % interpolate volume values at these points
    % ( not sure why have to swap XI,YI here )
    VI = interp3(vol,YI,XI,ZI,'*linear');
    
    % do we have a binary volume (no values > 1)
    if isempty(binvol),
        maxindex = max(find(VI>0));
        if maxindex,
            R(v) = Rarray(maxindex);
        end
    else
        % find the finite values of VI
        index = max(find(VI(isfinite(VI))));
        if index,
            
            % Find nearest volume value to the required fit value
            nearest = max(find(VI >= fit.val));
            
            %[ nearest, value ] = NearestArrayPoint( VI, fit.val );
            
            % Check this nearest index against a differential
            % negative peak value
            %diffVI = diff(VI);
            %if max(VI) > 1,
            %    diffindex = find(diffVI < -20);
            %else
            % probably a binary volume
            %    diffindex = find(diffVI < 0);
            %end
            
            % now set d
            if nearest,
                R(v) = Rarray(nearest);
            end
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Constrain relocation by fit.vattr, 
    % some % of distance from neighbours
    
    vi = find(FV.edge(v,:));  % the neighbours' indices
    X = FV.vertices(vi,1);    % the neighbours' vertices
    Y = FV.vertices(vi,2);
    Z = FV.vertices(vi,3);
    
    % Find neighbour distances
    RN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
    % Find mean distance of neighbours
    neighbour_distances_mean = mean(RN);
    
    minattr = fit.vattr;
    maxattr = 1 + fit.vattr;
    
    if R(v) < (minattr * neighbour_distances_mean),
        R(v) = minattr * neighbour_distances_mean;
    end
    if R(v) > (maxattr * neighbour_distances_mean),
        R(v) = maxattr * neighbour_distances_mean;
    end
    if R(v) == 0, R(v) = neighbour_distances_mean; end
    
    % relocate vertex to new distance
    x = (l * R(v)) + xo;
    y = (m * R(v)) + yo;
    z = (n * R(v)) + zo;
    
    FV.vertices(v,:) = [ x y z ];
    
    % Find intensity value at this distance
    intensityAtR(v) = interp1(Rarray,VI,R(v),'linear');
    
end

if isempty(binvol),
    % check outliers and smooth twice for binary volumes
    FV = vertex_outliers(FV, R, origin);
    FV = mesh_smooth(FV,origin,fit.vattr);
end
FV = mesh_smooth(FV,origin,fit.vattr);

return









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV, intensityAtR, R] = shrink_wrap(FV,vol,origin,fit),

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

Nvert = size(FV.vertices,1);
progress = round(.1 * Nvert);

% Initialise difference intensity & distance arrays
R = zeros(Nvert,1);
intensityAtR = R;

% Check for binary volume data, if empty, binary
binvol = find(vol > 1);

Imin = 0;
Imax = max(max(max(vol)));

% Manipulate each vertex
tic
for v = 1:Nvert,
  
  if v > progress,
    fprintf('...shrink-wrap processed %4d of %4d vertices',progress,Nvert);
    t = toc; fprintf(' (%5.2f sec)\n',t);
    progress = progress + progress;
  end
  
  % Find direction cosines for line from origin to vertex
  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
  
  % interpolate volume values at this point ( x,y are swapped here
  % because the Analyze volume is a left handed coordinate system )
  intensity_old = interp3(vol,y,x,z,'*nearest');
  
  % move vertex closer to the origin, in a radial direction
  r_change = r * 0.05;
  r_new = r - r_change;
  
  % Calculate new vertex coordinates
  x = (l * r_new) + xo; % cos alpha
  y = (m * r_new) + yo; % cos beta
  z = (n * r_new) + zo; % cos gamma
  
  % interpolate volume values at this point ( x,y are swapped here
  % because the Analyze volume is a left handed coordinate system )
  intensity_new = interp3(vol,y,x,z,'*nearest');
  
  intensity_dif = intensity_new - intensity_old;
  
  if intensity_dif == 0,
    % relocate vertex to new distance
    FV.vertices(v,1) = x;
    FV.vertices(v,2) = y;
    FV.vertices(v,3) = z;
    intensityAtR(v,1) = intensity_new;
    R(v) = r_new;
  elseif (fit.val - intensity_new) < (fit.val - intensity_old),
    % relocate vertex to new distance
    FV.vertices(v,1) = x;
    FV.vertices(v,2) = y;
    FV.vertices(v,3) = z;
    intensityAtR(v,1) = intensity_new;
    R(v) = r_new;
  else
    intensityAtR(v,1) = intensity_old;
    R(v) = r;
  end
  
  FV = constrain_vertex(FV,v,origin);
  
end

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

return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV] = constrain_vertex(FV,index,origin),

% This function adapts Smith, S. (2002), Fast robust automated brain
% extraction.  Human Brain Mapping, 17: 143-155.  It corresponds to update
% component 2: surface smoothness control.  It assumes that vertices are
% moved along a radial line from an origin, given by direction
% cosines, rather than calculating the surface normal vector.

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

v = FV.vertices(index,:);
x = FV.vertices(index,1);
y = FV.vertices(index,2);
z = FV.vertices(index,3);

% Find radial distance of vertex from origin
r = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );

% Calculate unit vector
v_unit_vector = ( v - origin ) / r;

% Find direction cosines for line from center to vertex
l = (x-xo)/r; % cos alpha
m = (y-yo)/r; % cos beta
n = (z-zo)/r; % cos gamma

% Find neighbouring vertex coordinates
vi = find(FV.edge(index,:));  % the indices
neighbour_vertices = FV.vertices(vi,:);
X = neighbour_vertices(:,1);
Y = neighbour_vertices(:,2);
Z = neighbour_vertices(:,3);

% Find neighbour radial distances
r_neighbours = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
r_neighbours_mean = mean(r_neighbours);

% Find difference in radial distance between the vertex of interest and its
% neighbours; this value approximates the magnitude of sn in 
% Smith (2002, eq. 1 to 4)
r_diff = r - r_neighbours_mean;

% Find the vector sn, in the direction of the vertex of interest, given the
% difference in radial distance between vertex and mean of neighbours
sn = r_diff * v_unit_vector;

% Find distances between vertex and neighbours, using edge lengths.
% The mean value is l in Smith (2002, eq. 4)
edge_distance = FV.edge(index,vi);
edge_distance_mean = mean(edge_distance);

% Calculate radius of local curvature, solve Smith (2002, eq. 4)
if r_diff,
  radius_of_curvature = (edge_distance_mean ^ 2) / (2 * r_diff);
else
  radius_of_curvature = 10000;
end

% Define limits for radius of curvature
radius_min =  3.33; % mm
radius_max = 10.00; % mm

% Sigmoid function parameters,
% "where E and F control the scale and offset of the sigmoid"
E = mean([(1 / radius_min),  (1 / radius_max)]);
F = 6 * ( (1 / radius_min) - (1 / radius_max) );

Fsigmoid = (1 + tanh( F * (1 / radius_of_curvature - E))) / 2;

% multiply sigmoid function by sn
move_vector = Fsigmoid * sn;

FV.vertices(index,:) = v + move_vector;

return






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [val] = vol_val(vol,x,y,z),

% This function just ensures that xyz are
% actually within the volume before trying
% to get a volume value

val = nan; % assume zero value

x = round(x);
y = round(y);
z = round(z);

if x > 0 & y > 0 & z > 0,
    
    % get bounds of vol
    xdim = size(vol,1);
    ydim = size(vol,2);
    zdim = size(vol,3);
    
    if x <= xdim & y <= ydim & z <= zdim,
        % OK return volume value at xyz
        val = vol(x,y,z);
    end
end

return






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV] = vertex_outliers(FV, R, origin),

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

% Screen FV for outlying vertices, using
% mean +/- 2 * stdev of distance from origin
DistMean = mean(R);
DistStDev = std(R);
DistMax = DistMean + 2 * DistStDev;
DistMin = DistMean - 2 * DistStDev;

for v = 1:size(FV.vertices,1),
    
    if R(v) >= DistMax,
        R(v) = DistMean;
        relocate = 1;
    elseif R(v) <= DistMin,
        R(v) = DistMean;
        relocate = 1;
    else
        relocate = 0;
    end
    
    if relocate,
        x = FV.vertices(v,1);
        y = FV.vertices(v,2);
        z = FV.vertices(v,3);
        
        % Find direction cosines for line from center to vertex
        d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
        l = (x-xo)/d; % cos alpha
        m = (y-yo)/d; % cos beta
        n = (z-zo)/d; % cos gamma
        
        % relocate vertex to this new distance
        x = (l * R(v)) + xo;
        y = (m * R(v)) + yo;
        z = (n * R(v)) + zo;
        
        FV.vertices(v,:) = [ x y z ];
    end
end
return