function [B] = loreta_lap3d(grid)

% LORETA_LAP3D computes the discrete spatial laplacian operator for a
% regular 3D grid.
%
% R.D. Pascual-Marqui (1999) Review of Methods for Solving the EEG Inverse Problem,
% International Journal of Bioelectromagnetism, 1(1):75-86.
% available online from http://www.ijbem.org/


if ~isfield(grid, 'xgrid') || ~isfield(grid, 'ygrid') || ~isfield(grid, 'zgrid')
  error('this only works for regular 3d dipole grids');
end

x = grid.xgrid;
y = grid.ygrid;
z = grid.zgrid;

nx = length(x);
ny = length(y);
nz = length(z);

dx = diff(x);
dy = diff(y);
dz = diff(z);

if any(dx~=dx(1)) || any(dy~=dy(1)) || any(dz~=dz(1))
  error('grid should be regularly spaced');
end

if dy(1)~=dx(1) || dz(1)~=dx(1)
  error('grid should have the same spacing in each dimension');
end

% recompute the position of each dipole to check for consistency with the 3d grid
[X, Y, Z] = ndgrid(x, y, z);
pos = [X(:) Y(:) Z(:)];
if any(pos(:) ~= grid.pos(:))
  error('the dipole positions are not consistent with the expected grid');
end

% determine the grid resolution, i.e. the spacing between grid points
d = dx(1);

% determine the number of grid locations
M = nx*ny*nz;

% this will hold the discrete spatial laplacian for the complete grid
% but only for a single orientation of the estimated current (i.e. Jx, Jy or Jz)
B = spalloc(M*3,M*3,7*M*3);

Jxsel = (1:3:(M*3));
Jysel = (1:3:(M*3)) + 1;
Jzsel = (1:3:(M*3)) + 2;

k = 0;
for iz=2:(nz-1)
  for iy=2:(ny-1)
    for ix=2:(nx-1)
      % determine the linear index of this grid point
      i = ix + (iy-1)*nx  + (iz-1)*nx*ny;
      if any(i==grid.inside)
        % determine the linear indices of its neighbours
        dpx = i + 1;
        dmx = i - 1;
        dpy = i + nx;
        dmy = i - nx;
        dpz = i + nx*ny;
        dmz = i - nx*ny;
        % the current distribution contains three components per grid point [Jx Jy Jz].
        B(Jxsel(i), Jxsel([i dmx dmy dmz dpx dpy dpz])) = [-6 1 1 1 1 1 1];
        B(Jysel(i), Jysel([i dmx dmy dmz dpx dpy dpz])) = [-6 1 1 1 1 1 1];
        B(Jzsel(i), Jzsel([i dmx dmy dmz dpx dpy dpz])) = [-6 1 1 1 1 1 1];
        k = k + 1;
      else
        % do not compute the laplacian for grid points that are not inside the brain
      end
    end
  end
end

fprintf('laplacian matrix is filled for %.0f%%\n', 100*k./M);

% scale the second order derivative with the distance squared
B = B ./ (d^2);