tri_interp

(Toolbox/tri_interp.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[H] = tri_interp(dims,verts,Faces)  ;

Help Text

TRI_INTERP - Compute the interpolation matrix from a cortical tessellation to the MRI volume
 function   [H] = tri_interp(dims,verts,Faces)  ;
 TRI_INTERP : Compute the interpolation matrix from a cortical tessellation to the MRI volume

 dims - 3-element vector : size of the orginal (MR) volume
 verts - Mx3 : Coordinates (integers) of the cortical vertices in the MR cube
 Faces - Nx1 : Connectivity matrix of the vertices to form the triangles

 H : interpolation matrix from values on triangular tessellation to voxels

 How to use this ? ----------------------------------------------------------------------------
 H is a one-time computed large but sparse matrix (nvoxels X nvertices) and is used as follows.
 Values of currents at each MRI voxel are obtained by
 MRIcurrent = H x SurfaceCurrent
 SurfaceCurrent is a nvertices-tall vector of estimated currents on the tessellated surface
 MRIcurrent is a nvoxels-tall vector of 3D-interpolated current values at each MR voxels.
 The index of each MRIcurrent entry corresponds to the index in the 3D cube of the MR.
 MRs in BrainStorm are stored as n1 X n2 X n3 data arrays.
 If my_MRI is one of these, values of current at voxel i is MRIcurrent(i). i is the index
 in the original volume taken in lexicographical (i.e. Matlab's natural) order. Don't bother too much about this;
 the only concern one should have is to have the vertex coordinates expressed in the proper MR indices. Use PCS2MRI in that respect.

Cross-Reference Information

vec                                C:\BrainStorm_2001\Toolbox\vec.m

mritool_cb                         C:\BrainStorm_2001\Toolbox\mritool_cb.m
surf2vol_interp                    C:\BrainStorm_2001\Toolbox\surf2vol_interp.m

Listing of function C:\BrainStorm_2001\Toolbox\tri_interp.m

function   [H] = tri_interp(dims,verts,Faces)  ;
%TRI_INTERP - Compute the interpolation matrix from a cortical tessellation to the MRI volume
% function   [H] = tri_interp(dims,verts,Faces)  ;
% TRI_INTERP : Compute the interpolation matrix from a cortical tessellation to the MRI volume
%
% dims - 3-element vector : size of the orginal (MR) volume
% verts - Mx3 : Coordinates (integers) of the cortical vertices in the MR cube
% Faces - Nx1 : Connectivity matrix of the vertices to form the triangles
%
% H : interpolation matrix from values on triangular tessellation to voxels
%
% How to use this ? ----------------------------------------------------------------------------
% H is a one-time computed large but sparse matrix (nvoxels X nvertices) and is used as follows.
% Values of currents at each MRI voxel are obtained by
% MRIcurrent = H x SurfaceCurrent
% SurfaceCurrent is a nvertices-tall vector of estimated currents on the tessellated surface
% MRIcurrent is a nvoxels-tall vector of 3D-interpolated current values at each MR voxels.
% The index of each MRIcurrent entry corresponds to the index in the 3D cube of the MR.
% MRs in BrainStorm are stored as n1 X n2 X n3 data arrays.
% If my_MRI is one of these, values of current at voxel i is MRIcurrent(i). i is the index
% in the original volume taken in lexicographical (i.e. Matlab's natural) order. Don't bother too much about this;
% the only concern one should have is to have the vertex coordinates expressed in the proper MR indices. Use PCS2MRI in that respect.

%<autobegin> ---------------------- 26-May-2004 11:34:41 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Subfunctions in this file, in order of occurrence in file:
%   c = cross(a,b);
%
% At Check-in: $Author: Mosher $  $Revision: 16 $  $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:34:41 -----------------------


 %<autobegin> ---------------------- 18-Nov-2003 22:02:32 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Subfunctions in this file, in order of occurrence in file:
%   c = cross(a,b);
%
% At Check-in: $Author: Mosher $  $Revision: 16 $  $Date: 5/26/04 10:02a $
%
% Overall BrainStorm authors are:
% ** Dr. John C. Mosher,  Los Alamos National Laboratory
% ** Dr. Sylvain Baillet, CNRS Cognitive Neuroscience & Brain Imaging Laboratory
%
% Copyright (c) 2003 BrainStorm MMIII by the University of Southern California
% Principal Investigator:
% ** Professor Richard M. Leahy, USC Signal & Image Processing Institute
%
% Source code may not be distributed in original or altered form.
% For license and copyright notices, type 'bst_splashscreen' in Matlab,
%   or see BrainStorm website at http://neuroimage.usc.edu,
%   or email Professor Richard M. Leahy at leahy@sipi.usc.edu
%   for further information and permissions.
%<autoend> ------------------------ 18-Nov-2003 22:02:32 -----------------------

% /---Script Author--------------------------------------\
% |                                                      |
% | *** Sylvain Baillet Ph.D.                            |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   | 
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% |                                                      |
% \------------------------------------------------------/
%  
% Date of creation: October 2002

Faces = Faces(1:1:end,:); % CHEAT - testing
if 0
    
    dims = [256 256 256];% Size of the Cube
    %% Triangle
    X = [3 3 8;3 8 5]';
    Y = [3 8 5; 3 3 8]';
    Z = [3 3 8; 1 2 4]' ;
    CData = ([1 .5 0; 0 1 .2]');
    figure
    colormap hot
    hsp = patch(X,Y,Z,'r','Cdata',CData,'Facecolor','interp','edgecolor','none');
    hold on
    
    %% Get Patch properties
    Faces = get(hsp,'Faces');
    verts = get(hsp,'Vertices');
    V = get(hsp,'FaceVertexCdata');
    
end

%Flag = 1-spones(dims(1),dims(2),dims(3));
dMAX = 2/sqrt(3); % Distance Max Authorized between any voxel and the plane

if 0
    minn = min(verts); 
    maxx = max(verts); % Where is the brain in the MR volume
    u = [minn(1):maxx(1)];
    v = [minn(2):maxx(2)];
    w = [minn(3):maxx(3)];
end

u = [1:dims(1)];
v = [1:dims(2)];
w = [1:dims(3)];

nFaces = size(Faces,1);


hwait = waitbar(0,'Computing the interpolation matrix...');
steps = linspace(1,nFaces,100); % Just for the waitbar
istep = 1;

X = cell(size(Faces,1),1);
MAT = X;

%H = sparse([],[],[],prod(dims),size(verts,1),0);
%iFlag  = spalloc(prod(dims),1,20*size(verts,1));
%profile on
for tri = 1:nFaces% For every triangle
    faces = (Faces(tri,:));
    bary = mean(verts(faces,:),1); % Barycenter
    %hold on
    %scatter3(bary(1),bary(2),bary(3),'r')
    R = max([abs(verts(faces,1)-bary(1)),abs(verts(faces,2)-bary(2)),...
            abs(verts(faces,3)-bary(3))]); % Vertex block around the current triangle    
    
    R(R<1) = 1;
    
    % Find the voxels within the triangle
    I = u(abs(bary(1)-u)<=R(1));
    J = v(abs(bary(2)-v)<=R(2));
    K = w(abs(bary(3)-w)<=R(3));
    
    if 0
        % Check if OK for barycenter
        I(end+1) = bary(1);
        J(end+1) = bary(2);
        K(end+1) = bary(3);
    end
    
    [II,JJ,KK] = meshgrid(I,J,K);
    II = squeeze(II);
    JJ = squeeze(JJ);
    KK = squeeze(KK);
    
    if 0
        newind = sub2ind(dims,II(:),JJ(:),KK(:));
        newind = setdiff(newind,find(iFlag)); % Don't compute what has already been done
        [II JJ KK] = ind2sub(dims,newind);
    end
    
    % Triangle vertex coordinates, with regards to 1st vertex
    r13 = verts(faces(3),:)' - verts(faces(1),:)'; % negative of r31
    r12 = verts(faces(2),:)' -  verts(faces(1),:)';            % from 1 to 2
    
    N = cross(r12+1e-10,r13+eps)'; % eps : it already occured that some triangles are NULL (vertex locations are colinear): must be an effect of the dowsampling on integer coordinate values
    dArea = norm(N); % Area of the triangle x 2
    N = N/dArea; % Normal to the triangle
    
    % Find the voxels "in" the triangle plane
    bary_p = [II(:) - bary(1), JJ(:) - bary(2),KK(:) - bary(3)];
    Iplane = bary_p * N';
    Iplane = find(abs(Iplane)<dMAX);
    II = II(Iplane);
    JJ = JJ(Iplane);
    KK = KK(Iplane);
    
    % Flag the voxels if they belong to the triangle
    r1p = [II(:) - verts(faces(1),1),JJ(:)- verts(faces(1),2),KK(:)-verts(faces(1),3)]';
    
    % Barycentric Coordinates
    s = (N*cross(r1p,repmat(r13,1,size(r1p,2))))/dArea;
    is = find(s>=0 & s<=1);
    
    t = (N*cross(repmat(r12,1,size(r1p(:,is),2)),r1p(:,is)))/dArea;
    it = find(t>=0 & t<=1);
    
    %itmp = intersect(is,it);
    itmp = is(it);
    
    %r = 1-s(itmp)-t(itmp);
    r = 1-s(itmp)-t(it);
    ir = find(r>=0 & r<=1);
    
    ind = itmp(ir);
    
    
    MAT{tri} =  [r(ir);s(ind);t(it(ir))];
    
    tmp = repmat(sub2ind(dims,II((ind)),JJ((ind)),KK((ind)))',3,1);
    X{tri} = tmp(:)'; 
    
    Y{tri} = repmat(faces,1,length(ind)); 
    
    %Flag(ind2) = 1;
    
    %Ctmp = H*GridAmp(Faces(tri,:)); 
    clear ind ind2 itmp ir is it r s t Iplane 
    
    %CInterp(sub2ind(dims,II(Iplane(ind)),JJ(Iplane(ind)),KK(Iplane(ind)))) = Ctmp;
    %IND(Iplane(ind)) = Iplane(ind); % Redundant ?
    
    %hold on
    % scatter3(II(Iplane(ind)),JJ(Iplane(ind)),KK(Iplane(ind)),50,Cinterp,'filled');
    % disp('')
    
    
    if ~rem(tri,500)
        waitbar(tri/nFaces,hwait)
        istep = istep+1;
    end
    
end
%profile report
close(hwait)
drawnow

clear faces

MAT = [MAT{:}];
MAT = MAT(:);

X = [X{:}];
X = X';

Y = [Y{:}];
Y = Y';

H = sparse(X, Y, MAT, prod(dims),size(verts,1)); 
clear X Y MAT verts 

%save interp_mat H

%X =unique(X);

tmp= sum(H,2);
ind = find(tmp>1+10*eps);
tmp = full(1./tmp(ind));

tempo = sparse(length(ind),size(H,2));

hwait = waitbar(0,'Fixing the interpolation matrix...');

blocks = unique(round(linspace(10,length(ind),20)));

for k = 1:length(blocks)
    waitbar(blocks(k)/blocks(end),hwait)
    if k > 1
        vec = blocks(k-1)+1:blocks(k);
    else
        vec = 1:blocks(1);
    end
    
    H(ind(vec),:) = spdiags(tmp(vec),0,length(vec),length(vec))*H(ind(vec),:);
    
end
clear tmp
delete(hwait)
disp(' ')

function c = cross(a,b);
c = [a(2,:).*b(3,:)-a(3,:).*b(2,:)
    a(3,:).*b(1,:)-a(1,:).*b(3,:)
    a(1,:).*b(2,:)-a(2,:).*b(1,:)];
return