[Master Index] [Index for Toolbox]

# 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

This function calls

    vec                                C:\BrainStorm_2001\Toolbox\vec.m

This function is called by

    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

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON