mesh_laplacian_interp - Computes the zero Laplacian interpolation matrix Useage: [int, keepindex, repindex] = mesh_laplacian_interp(lap, index) This function calculates an interpolation matrix that provides the coefficients for the calculation of potential values at all unknown vertices of a mesh, given known potential values at a subset of the mesh vertices (at 'index'). The interpolation solution is constrained by a minimal norm of the Laplacian of the mesh. See the reference below for details. 'lap' is the laplacian matrix for the full mesh (see mesh_laplacian) 'int' is the matrix which interpolates from the points in 'index' to the full mesh. 'index' is a row vector of indices into a subset of the vertices used to calculate 'lap'. This subset is where the electric potential is known and usually corresponds to the given electrode vertices, eg: index = dsearchn(scalpvert,elecvert)'; If 'index' contains repeated indices, this is a serious problem but this function will try to continue with only the unique indices. The 'keepindex' array can be used to select these. The 'repindex' array is the repeated indices. Interpolations can be done using matrix 'int', eg: [int, keepindex, repindex] = mesh_laplacian_interp(lap,index); if isempty(repindex), Vint = int * Vknown; else Vint = int * Vknown(keepindex); end This attempts to implement interpolation method B (p. 336) of Oostendorp T, Oosterom A & Huiskamp G (1989), Interpolation on a triangulated 3D surface. Journal of Computational Physics, 80: 331-343.
0001 function [int, keepindex, repindex] = mesh_laplacian_interp(lap, index); 0002 0003 % mesh_laplacian_interp - Computes the zero Laplacian interpolation matrix 0004 % 0005 % Useage: [int, keepindex, repindex] = mesh_laplacian_interp(lap, index) 0006 % 0007 % This function calculates an interpolation matrix that provides 0008 % the coefficients for the calculation of potential values at all 0009 % unknown vertices of a mesh, given known potential values at 0010 % a subset of the mesh vertices (at 'index'). The interpolation 0011 % solution is constrained by a minimal norm of the Laplacian 0012 % of the mesh. See the reference below for details. 0013 % 0014 % 'lap' is the laplacian matrix for the full mesh (see mesh_laplacian) 0015 % 'int' is the matrix which interpolates from the points in 'index' 0016 % to the full mesh. 'index' is a row vector of indices into a 0017 % subset of the vertices used to calculate 'lap'. This subset 0018 % is where the electric potential is known and usually corresponds 0019 % to the given electrode vertices, eg: 0020 % 0021 % index = dsearchn(scalpvert,elecvert)'; 0022 % 0023 % If 'index' contains repeated indices, this is a serious problem 0024 % but this function will try to continue with only the unique indices. 0025 % The 'keepindex' array can be used to select these. 0026 % The 'repindex' array is the repeated indices. 0027 % 0028 % Interpolations can be done using matrix 'int', eg: 0029 % 0030 % [int, keepindex, repindex] = mesh_laplacian_interp(lap,index); 0031 % if isempty(repindex), 0032 % Vint = int * Vknown; 0033 % else 0034 % Vint = int * Vknown(keepindex); 0035 % end 0036 % 0037 % This attempts to implement interpolation method B (p. 336) of 0038 % Oostendorp T, Oosterom A & Huiskamp G (1989), 0039 % Interpolation on a triangulated 3D surface. 0040 % Journal of Computational Physics, 80: 331-343. 0041 % 0042 0043 % $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $ 0044 0045 % Licence: GNU GPL, no implied or express warranties 0046 % History: (c) 04/2002 Robert Oostenveld 0047 % - agreed to release 'lapint' under GNU GPL 0048 % 04/2002, Darren.Weber_at_radiology.ucsf.edu 0049 % - introduced check for index replications and 0050 % adjusted calculations/output accordingly 0051 % - converted lap to sparse matrix and solution 0052 % of interpolation matrix with \ operator 0053 % - accepts sparse lap input and returns sparse int 0054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0055 0056 if size(lap,1)~=size(lap,2), error('MESH_LAPLACIAN_INTERP: lap matrix is not square'); end 0057 0058 if issparse(lap), lap = full(lap); end 0059 0060 tic 0061 0062 % -- Check for any duplicate values in the known index 0063 0064 if size(index,1)~=1, index = index'; end 0065 0066 % Remove any replicate indices from 'index' 0067 [KnownIndex, i, i] = union(index,index); 0068 if ~isequal(length(KnownIndex),length(index)), 0069 fprintf('\a\n\nMESH_LAPLACIAN_INTERP: Warning:\n\nTrimming duplicate values from index, use keepindex!\n\n'); 0070 end 0071 keepindex = sort(i); 0072 repindex = setdiff(1:length(index),sort(i)); 0073 0074 KnownIndex = index(keepindex); % unsort KnownIndex 0075 clear index 0076 0077 k = length(KnownIndex); 0078 n = length(lap); 0079 0080 fprintf('MESH_LAPLACIAN_INTERP: Calc Interpolation matrix for %5d to %5d vertices...',k,n); 0081 0082 % find 'unknown' indices of lap matrix 0083 UnknownIndex = setdiff(1:n, KnownIndex); 0084 0085 % reshuffle rows & columns of lap matrix 0086 lapi = [KnownIndex, UnknownIndex]; 0087 lap = lap(lapi, :); % rows 0088 lap = lap(:, lapi); % columns 0089 0090 % Segregate known/unknown portions of lap 0091 k = length(KnownIndex); 0092 n = length(lap); 0093 0094 L11 = lap(1:k ,1:k ); 0095 L12 = lap(1:k ,(k+1):n); 0096 L21 = lap((k+1):n,1:k ); 0097 L22 = lap((k+1):n,(k+1):n); 0098 0099 clear lap lapi; % tidy up some memory 0100 0101 % Convert to sparse for quicker computation 0102 A = sparse([L12; L22]); 0103 B = sparse([L11; L21]); 0104 0105 clear L11 L12 L21 L22; % tidy up some memory 0106 0107 %int = -pinv(A) * B; % cannot use pinv with sparse matrix 0108 int = -A \ B; 0109 0110 % Convert result back to full matrix 0111 int = full(int); 0112 0113 % append the interpolating piece to the identity matrix 0114 % these take care of the known potentials 0115 int = [eye(k); int]; 0116 0117 % reshuffle the columns of the interpolating matrix 0118 [tmp, order] = sort(KnownIndex); 0119 int = int(:,order); 0120 0121 % reshuffle the rows of the interpolating matrix 0122 [tmp, order] = sort([KnownIndex, UnknownIndex]); 0123 int = int(order, :); 0124 0125 int = sparse(int); 0126 0127 t = toc; 0128 fprintf('done (%6.2f sec).\n',t); 0129 0130 return