Matmap Documentation
[Documentation home]

The Boundary Element Method

Currently the files in the BEM directory, solve the forward problem for the case where the source is defined as a surface with potentials, i.e. the epicardial surface. There are two functions that perform the boundary element calculation:

bemCheckModel(model)   - This function checks the integrity of a model and tries to correct for mistakes in the model definition
bemMatrixPP(model) - The forward computation, i.e. the computation of the forward matrix.

These two functions form the main function in computing the forward solutions. The model is setup as follows:
model itself is a structured matrix, currently containing only one field, called 'surface'. Surface is a cell array with each cell containing the definition of a surface. In this model definition the surfaces are assumed to be indexed from the outer one towards the most inner boundary, i.e., the first cell of the surfaces array contains the outermost surface and the last one contains the inner most surface, e.g. the epicardial one.

In order for the boundary element method to be used all surfaces need to be closed.


Example of defining a model

Assume there are two geometry files, cage.geom and torso.geom:
The geometries can be loaded as following
>> indextorso = ioReadGEOM('torso.geom');
>> indexcage = ioReadGEOM('cage.geom');

Now a model can be constructed
>> model.surface{1} = GEOM{indextorso};
>> model.surface{2} = GEOM{indexcage};

To complete a model, the conductivities need to be defined. The conductivities are defined for each surface, the conductivity outside of this surface and the conductivity on the inside. In our example the conductivity of the torso is set to 0.2 S/m:

>> model.surface{1}.sigma = [0 0.2];
>> model.surface{2}.sigma = [0.2 0];
We defined the conductivity inside the inner surface as zero, actually any value can be chosen here as this value will fall out of the equation. It is used in the  sub-functions of numerical computation, but in the end it should fall out of the equation

Now we are ready to do some boundary element computation


Checking the correctness of a model

Before computing the transfer function it is good to test the model for its correctness. The function bemCheckModel(model) will check the integrity of the model; if it encounters some minor mistakes like a surface is defined Clock Wise instead of Counter Clock Wise it will correct these mistakes and return a corrected version of the model. Hence we use the following function:

>> model = bemCheckModel(model)

If this function returns the model, it is assumed to be OK.
The following checks were performed by the function:
 The following checks are not performed:
Calculating the forward matrix

The following function call will now compute the forward matrix:
>> T = bemMatrixPP(model);

Here T is the matrix that takes the nodes at the inner surface as a source and computes the potentials at all the other nodes in the matrix assuming there is no current flowing out of the model over the outer most boundary. If the number of nodes on the inner surface is N and the combined number of nodes on the other surfaces is M, T will be a M times N matrix.

Assume Ui is the potential distribution on the inner boundary then the other ones can be computed as follows:
>> Uo = T*Ui;

In case of three or more surfaces the Uo vector needs to be distributed over multiple surfaces starting with second most inner one and going to the outermost one.
The potentials are ordered according to the node ordering in the pts array of that surface.


Validation of the BEM method

Since the code may be altered of the course of some further development, the function bemValidate() checks whether the BEM code still produces accurate solutions by comparing the solutions of a spherical boundary element model versus an analytical solution:
 
>> bemValidate;
This should give a RDM value in the order of 7.4E-4 and a magnification error in the order of 1.0007. In case the error is smaller, the boundary element method apparently works even better than when this function was written. In case the errors are far bigger, something broke in the code and needs an urgent fix.

Of course this is only one case. For the case of a more thorough investigation in the accuracy of the BEM method, two functions are provided to test the BEM forward code with spherical models:

bemGenerateSphere() - Generates a multi layer spherical model
anaSolveSphere()          - Analytical solution (expansion in spherical harmonics) of the multi layer concentric spherical model

for example:

>> model = bemGenerateSphere([1 0.5],[0 0.2 0.2],0.15)
This function will generate a model with two surfaces, both are spherical with radii 1 and 0.5 respectively. The compartment has a conductivity of 0.2 S/m and the mesh used  in the spherical tesselation has a mean edge length of 0.15*radius of the mesh.

In the analytical model we assume that the inner volume of the model has a conductivity of 0.2 as well and as a source we use a dipole. Then we calculate the potential on both the inner and the outer surface. Now to test the boundary element method we use the potential at the inner surface as source to compute the one on the outer surface:

>>Ui = anaSolveSphere('U',model.surface{2}.pts,[0 0 0],[0 0 1],[1],[0.2 0],40);
The first parameter specifies that we want to compute potentials ('U'), similarly the function will compute the magnetic field ('H') and the current density field ('J') as well. The second parameter are the positions at which the potentials need to be computed; the third one is the location of the dipole, in this case the origin of the sphere, the fourth one is the direction of the current dipole, in this case point up along the z-axis; the fifth one is the radii in the model (note that here the opposite order is used in comparison to bemGeneratesphere (due to historical reasons)) and the sixth parameter is the corresponding conductivity vector. The last one is the number spherical harmonics used to evaluate the analytical solution. Normally 40 should be enough, going to really high high numbers e.g. 100 may cause instabilities in the algorithm computing the Legendre Polynomials, and hence may not be more accurate than the one using 40 harmonics.

 
Similarly the potentials at the outer bounday can be computed
>>Uo = anaSolveSphere('U',model.surface{1}.pts,[0 0 0],[0 0 1],[1],[0.2 0],40);

The potential at the outer boundary can be computed as well using the boundary element method:
>> Uobem = T*Ui;

Comparison of Uo and Uobem will approximate the accuracy of the BEM method.


Displaying data directly in matlab

The function bemPlotSurface() can be used to plot data at a boundary.

Plotting a surface:
>> figure; bemPlotSurface(model.surfaces{1});
This will plot the first surface

Plotting a surface with a potential distribution:
>> figure; bemPlotSurface(model.surfaces{1},Uo);
This will also add the solution to the graph.

For more on the functionality on bemPlotSurface, see the help in the reference section


Error norms

To assist in the evaluation of the accuracy of forward solutions the following error norms have been implemented:

errRMS   - RMS error between two solutions
errRDM   - Relative Difference Measure
errRDMS - Relative Difference Measure *, the magnification error has been removed from this measure and it is only sensitive to patterns
errMAG   - Magnification error (1 is a perfect match between solutions)



[Documentation home]