mesh_shrinkwrap.m

Matlab下的EEG处理程序库

    YV = FV.vertices(:,2);
    ZV = FV.vertices(:,3);
    DV = sqrt( (XV-xo).^2 + (YV-yo).^2 + (ZV-zo).^2 );
    LV = (XV-xo)./DV; % cos alpha
    MV = (YV-yo)./DV; % cos beta
    NV = (ZV-zo)./DV; % 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 = DV(v);
        l = LV(v); % cos alpha
        m = MV(v); % cos beta
        n = NV(v); % cos gamma
        
        % find discrete points between origin
        % and vertex + 50% of vertex distance
        points = 250;
        
        Darray = linspace(0,(d + .2 * d),points);
        
        L = repmat(l,1,points);
        M = repmat(m,1,points);
        N = repmat(n,1,points);
        
        XI = (L .* Darray) + xo;
        YI = (M .* Darray) + yo;
        ZI = (N .* Darray) + zo;
        
        % interpolate volume values at these points
        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,
                D(v) = Darray(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,
                    D(v) = Darray(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
        DN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
        % Find mean distance of neighbours
        DNmean = mean(DN);
        
        minattr = fit.vattr;
        maxattr = 1 + fit.vattr;
        
        if D(v) < (minattr * DNmean),
            D(v) = minattr * DNmean;
        end
        if D(v) > (maxattr * DNmean),
            D(v) = maxattr * DNmean;
        end
        if D(v) == 0, D(v) = DNmean; end
        
        % relocate vertex to new distance
        x = (l * D(v)) + xo;
        y = (m * D(v)) + yo;
        z = (n * D(v)) + zo;
        
        FV.vertices(v,:) = [ x y z ];
        
        % Find intensity value at this distance
        Dval(v) = interp1(Darray,VI,D(v),'linear');
        
    end
    
    if isempty(binvol),
        % check outliers and smooth twice for binary volumes
        FV = vertex_outliers(FV, D, origin);
        FV = mesh_smooth(FV,origin,fit.vattr);
    end
    FV = mesh_smooth(FV,origin,fit.vattr);
    
return







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Currently not calling this function (Oct 02)
function [FV, Dval, D] = shrink_wrap(FV,vol,origin,fit),
    
    xo = origin(1); yo = origin(2); zo = origin(3);
    
    Nvert = size(FV.vertices,1);
    
    Dval = zeros(Nvert,1);
    D    = Dval;
    
    for v = 1:Nvert,
        
        x = FV.vertices(v,1);
        y = FV.vertices(v,2);
        z = FV.vertices(v,3);
        
        % Find distance of vertex from origin
        d = sqrt( (x-xo)^2 + (y-yo)^2 + (z-zo)^2 );
        
        % check whether vertex is already at fit.val in vol
        
        volval = vol_val(vol,x,y,z);
        
        if isnan(volval) | (volval < fit.val - fit.tol) | (volval > fit.val + fit.tol),
            
            % Find direction cosines for line from centre to vertex
            l = (x-xo)/d; % cos alpha
            m = (y-yo)/d; % cos beta
            n = (z-zo)/d; % cos gamma
            
            % now modify d by fit.dist
            if isnan(volval) | (volval < fit.val - fit.tol),
                d = d - (d * fit.vattr);
            else
                d = d + (d * fit.vattr);
            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
            DN = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2 );
            % Find mean distance of neighbours
            DNmean = mean(DN);
            
            minattr = fit.vattr;
            maxattr = 1 + fit.vattr;
            
            if d < (minattr * DNmean),
                d = minattr * DNmean;
            elseif d > (maxattr * DNmean),
                d = maxattr * DNmean;
            end
            
            % locate vertex at this new distance
            x = (l * d) + xo;
            y = (m * d) + yo;
            z = (n * d) + zo;
            
            FV.vertices(v,:) = [ x y z ];
        end
        
        Dval(v) = vol_val(vol,x,y,z);
        D(v)    = d;
    end
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Currently not calling this function (Oct 02)
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
        Xv = size(vol,1);
        Yv = size(vol,2);
        Zv = size(vol,3);
        
        if x <= Xv & y <= Yv & z <= Zv,
            % OK return volume value at xyz
            val = vol(x,y,z);
        end
    end
    
return



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV] = vertex_outliers(FV, D, 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(D);
	DistStDev = std(D);
	DistMax = DistMean + 2 * DistStDev;
	DistMin = DistMean - 2 * DistStDev;
	
	for v = 1:size(FV.vertices,1),
        
        if D(v) >= DistMax,
            D(v) = DistMean;
            relocate = 1;
        elseif D(v) <= DistMin,
            D(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 centre 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 * D(v)) + xo;
            y = (m * D(v)) + yo;
            z = (n * D(v)) + zo;
            
            FV.vertices(v,:) = [ x y z ];
        end
	end
return