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
