*3.2.2. Treatment efficiency analysis*

26 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

Furthermore, we illustrated how real modulus of vector of intensity of electric field at phase 0° is interpreted using Matlab (second image in Fig. 4a.). This is the most basic interpretation

Note that this kind of results interpretation is much more flexible than the interpretation allowed by post processing tools in commercial EM simulators. In the following example we shall show how to work with time dependency of phasors. Since the results of EM field simulator are extracted when the steady state is reached time dependency is reduced to angle of phasors depicting the field of vectors. Through the following method we can alter phase of those phasors and show real part and imaginary part through one period. The results can be seen in Fig. 6. Figure 7. shows example of data processing to achieve this.

**Figure 6.** Phase Shifted Data – real part of vector *E* [dB] (left phase = 0°, right phase = 90°)

Note: In many EM field simulators you may encounter various errors. Pay special attention to the data structure of your exported data since it may not be useful in the way we have shown here (e.g. real and imaginary parts are exported as absolute values so the vital

Note: In the part of the script (Fig. 7.) where lowerThan variable is used we are changing the range of values. Since there are parts of model where values of intensity of electric or

**Figure 5.** Colormap Editor Window

of obtained data we can do.

information about phase is lost).

In the following example we will go through one of many useful applications of EM simulators today – evaluation of hyperthermia cancer treatment. Generally, hyperthermia is a method through which tissue is overheated (usually using microwave energy) and cells die (principle of this method is that energy is focused into the cancerous tissue which is less perfused and thus it is more heated – temperature in treated area rises above levels that trigger cell apoptosis). For more information on microwave hyperthermia see for example (Vrba & Oppl., 2008).

Model of simulated experiment can be seen in Figure 8. In this example we use waveguide applicator which is fed by coaxial line ending in protruding inner wire which is located near the shorted end of the waveguide section. There is a horn aperture which helps focus microwave energy to the desired area.

Relatively complex and complicated structure of human body is replaced with so called "phantom" that represents simple muscle tissue. In the model there is a tumor located 1 cm below the surface of phantom. This tumor has the same dielectric parameters as the surrounding muscle tissue (usually the only difference in simulations between muscle tissue and tumorous tissue is in their perfusion, heat transfer rate and heat generation rate – generated heat in tumorous tissue is usually transferred slower than in physiological surrounding tissue).

Additionally there is a water bolus which serves as a coolant body protecting the surface tissue of patient and moving the maximum of temperature to the lower layers. In this example we use source at 434 MHz and the applicator is filled with water (required dimensions of the applicator are effectively lower and impedance matching between waveguide-bolus and body are much better). For more information on waveguide hyperthermia applicators please refer to (Vera et al., 2006).

In this simulation we again extract rough data from EM simulator and we process them further using Matlab. We again have intensity of electric and magnetic field defined in every element of the model. For effective analysis of the treatment we need to evaluate SAR (Specific Absorption Ratio) which describes how much power is absorbed in a weight unit [W/kg]. For more information on SAR see http://www.ets-lindgren.com/pdf/sar\_lo.pdf.

**Figure 7.** Method of phase shifting results

To determine how SAR is distributed we need to use this formula.

$$SAR = \frac{\sigma}{2\rho} \left| E \right|^2 \tag{4}$$

As we see, in this case we can very well utilize *RMS|E|*. From formula (3) we can say that SAR is defined by following equation.

$$SAR = \frac{\sigma}{\rho}RMS \left| E \right|^2 \tag{5}$$

SAR is thus depending on RMS value of intensity of electric field, on conductivity of a material and on its density. We need to obtain those values somehow. In our case whole phantom has homogeneous density and electric conductivity thus we can only mask matrix of *RMS|E|* and multiply each element by coefficient produced by ratio of σ and ρ.

**Figure 8.** Model of microwave hyperthermia at 434 MHz

**Figure 7.** Method of phase shifting results

SAR is defined by following equation.

To determine how SAR is distributed we need to use this formula.

2

(4)

(5)

2 *SAR E*

As we see, in this case we can very well utilize *RMS|E|*. From formula (3) we can say that

<sup>2</sup> *SAR RMS E*

SAR is thus depending on RMS value of intensity of electric field, on conductivity of a material and on its density. We need to obtain those values somehow. In our case whole phantom has homogeneous density and electric conductivity thus we can only mask matrix

of *RMS|E|* and multiply each element by coefficient produced by ratio of σ and ρ.

Note: There are other options how to do that. For example some EM simulator allow users to export matrix of dielectric parameters of a model (this will produce similar matrix to our masking matrix).

Now we should look into how to produce masking matrix efficiently. From our model we know that every element of the voxelized model is phantom if its value of y axis is higher than 120 mm (we can find out this in CAD part of EM simulator). We can prepare our three dimensional matrix of *RMS|E|* according to previous sections. We also have actual axis of this matrix which we obtained from exported file as well. We simply need to identify which element is at the 120 mm value and then we shall mask the elements with lower value of y axis.

An elegant way to find the index at y axis of element which has the value nearest to 120 mm (note that the value is not usually exactly 120 so if we looked for the exact match we would probably fail) is to use function min. Example of this method is as follows.

**Figure 9.** Example of finding the closest element position

In this example, we can illustrate one other useful feature of Matlab processing of EM simulator results. We can instantly normalize obtained results. Usually it is possible to extract actual power in the source which was simulated (e.g. 0.002496 W) and we can utilize this to obtain normalization coefficient. Power normalization coefficient is given by the following formula.

**Figure 10.** Covering Field of *RMS|E|* with Mask

$$\text{Coef } f\_{power} = \frac{\text{Power to which we want to normalize}}{\text{Actual Simulated Power}} \tag{6}$$

Since this is the coefficient of power ratio we need to use coefficient of intensity of electric field ratio which is the square root of coefficient of power (see the equation below).

$$\text{coef } f\_{power} \approx \frac{E^2\_{\text{norm}}}{E^2\_{\text{actual}}} \Longrightarrow \text{coef } f\_E = \sqrt{\text{coef } f\_{power}} \tag{7}$$

Simply by multiplying the matrix of data by this coefficient we obtain normalized data to the desired value of delivered power.

Now we can show several interpretations of results we obtained (Z section in the middle of the model). Thanks to Matlab's flexibility and versatility we can highlight the results in the way that is far more efficient than by using post processing of the EM field simulator. In the following figures we show intensity of electric field and SAR in the treated area.

As you can see in figure 11. each colour scheme highlights different thing (as for the first three). In the fourth graph we used custom colormap to show some kind of a critical zone. In this example we expected that higher values of SAR than 1000 W/kg are considered to be dangerous in some regions (note that these critical limits differ significantly in every application – many factors contribute, e.g. vital or sensitive organs near the treated area etc.) and thus we highlighted the zone with higher values using bright red colour. This may be very important but commercial simulation software generally omits such option.

It may be also very useful to show obtained data only in the region of the tumour or to show values in form of several slices, semi-transparent layers or iso-surface view. Now we shall show how to prepare masking matrix for some simple geometrical objects representing tumours. For more information on the mentioned advanced techniques of results representation please see section 4. Advanced Viewing Techniques.

**Figure 11.** SAR [W/kg] in Different Colour Schemes (from left upper corner to right lower corner: hot, jet, lines, custom)

#### *3.2.3. Preparation of masking matrix*

30 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

**Figure 10.** Covering Field of *RMS|E|* with Mask

the desired value of delivered power.

such option.

*power*

*Powertowhichwe wanttonormalize coef f ActualSimulated Power* (6)

(7)

Since this is the coefficient of power ratio we need to use coefficient of intensity of electric

*power E power*

Simply by multiplying the matrix of data by this coefficient we obtain normalized data to

Now we can show several interpretations of results we obtained (Z section in the middle of the model). Thanks to Matlab's flexibility and versatility we can highlight the results in the way that is far more efficient than by using post processing of the EM field simulator. In the

As you can see in figure 11. each colour scheme highlights different thing (as for the first three). In the fourth graph we used custom colormap to show some kind of a critical zone. In this example we expected that higher values of SAR than 1000 W/kg are considered to be dangerous in some regions (note that these critical limits differ significantly in every application – many factors contribute, e.g. vital or sensitive organs near the treated area etc.) and thus we highlighted the zone with higher values using bright red colour. This may be very important but commercial simulation software generally omits

It may be also very useful to show obtained data only in the region of the tumour or to show values in form of several slices, semi-transparent layers or iso-surface view. Now we shall show how to prepare masking matrix for some simple geometrical objects representing tumours. For more information on the mentioned advanced techniques of results

representation please see section 4. Advanced Viewing Techniques.

field ratio which is the square root of coefficient of power (see the equation below).

*acctual <sup>E</sup> coef f coef f coef f <sup>E</sup>*

following figures we show intensity of electric field and SAR in the treated area.

2 2 *norm*

> Very simple masking (i.e. when the region that needs to be masked is in the form of a cuboid) was shown in previous section. Now we should look into masking regions round in shape (i.e. spheres, ellipsoids and other common shapes that parts of model can be represented by).

> Spheres are case of ellipsoid and we shall treat them as such. So our main aim now is to mask region generally in ellipsoidal form. Equation defining an ellipsoid (with its origin at *xc, yc, zc*) in the three dimensional coordinate system is shown below.

$$\frac{(\mathbf{x} - \mathbf{x}\_c)^2}{a^2} + \frac{(y - y\_c)^2}{b^2} + \frac{(z - z\_c)^2}{c^2} \le 1\tag{8}$$

The graphical representation of a general ellipsoid can be seen in figure 12. There we can see that *a*, *b* and *c* are semi-principal axes of ellipsoid and define its dimensions.

Albeit there is a built-in function designed to generate ellipsoid (see the Product Help of Matlab) we shall show you an easy way which allows you to generate desired masking matrix (e.g. with values 0 or 1).

**Figure 12.** Representation of Ellipsoid (from left: scheme, tri-axial, oblate, prolate)

As you can see in Fig. 13. we simply utilized the formula representing ellipsoid and all elements whose coordinates meet the given restrictions are filled with ones. The other elements remain zero. Constants elementsX, elementsY or elementsZ define the length of axis of computational domain (i.e. actualAxisX and so on). Variables xc, yc and zc define centre of our ellipsoid.

**Figure 13.** Defining Ellipsoid in the Masking Matrix

Now we can simply use generated masking matrix and multiply every element of result matrix by according element of masking matrix. Then we can look how our results are interpreted in detailed view in the tumour.

From masked matrix of SAR we can also easily extract total power radiated to the tumorous region. Total power lost in the healthy tissue is than given by the difference between total power (in the case of perfectly matched source) and power lost in the tumour.

**Figure 14.** SAR [W/kg] in Tumour (from left: Y-slice, X-slice, Z-slice)

**Figure 12.** Representation of Ellipsoid (from left: scheme, tri-axial, oblate, prolate)

matrix (e.g. with values 0 or 1).

and zc define centre of our ellipsoid.

**Figure 13.** Defining Ellipsoid in the Masking Matrix

interpreted in detailed view in the tumour.

Albeit there is a built-in function designed to generate ellipsoid (see the Product Help of Matlab) we shall show you an easy way which allows you to generate desired masking

As you can see in Fig. 13. we simply utilized the formula representing ellipsoid and all elements whose coordinates meet the given restrictions are filled with ones. The other elements remain zero. Constants elementsX, elementsY or elementsZ define the length of axis of computational domain (i.e. actualAxisX and so on). Variables xc, yc

Now we can simply use generated masking matrix and multiply every element of result matrix by according element of masking matrix. Then we can look how our results are

From masked matrix of SAR we can also easily extract total power radiated to the tumorous region. Total power lost in the healthy tissue is than given by the difference between total

power (in the case of perfectly matched source) and power lost in the tumour.

Note: In some cases it might be beneficial to set the colorbar range the same in all pictures (as shown in Fig. 14.). The simplest way to do this is to set value of [1,1] element of displayed matrix (i.e. a slice) to the maximum value of the whole three dimensional matrix (if needed, the minimum value should be set to element [max,max] of displayed matrix).

Additional information that we can obtain from masked result matrix is total absorbed power in the tumour. Since we know the dimensions of the tumour we can easily calculate its volume and use following formula.

$$P\_{\text{absorbed}} = Volume \sum\_{\text{all voxels}} RMS \left| E \right|^2 \sigma \tag{9}$$

The volume of standard ellipsoid is given as follows.

$$V = \frac{3}{4}\pi abc$$

For example in our case tumour has volume only 8.37766e-6 m3. Information on total absorbed power may be very efficient way of preliminary evaluation for the treatment planning. Through this method we can determine power lost in any part of a simulated model. There may be some intricate volumes which are not as easily described as ellipsoid. Then we need to export their volume or masking matrix from the EM simulator (if possible).

We can determine the volume of such complicated shapes by summing elementary volumes of voxels representing those shapes. This may be unnecessary in some more advanced EM field simulators since they allow us to export model data in many suitable forms for Matlab processing. But the method presented in the following text is universal and can serve for better understanding of the model, its grid and working principals of EM simulations.

First of all we need to define which voxels are occupied by the model we want to evaluate. For this example we shall use our previously generated ellipsoid. As you can see in Fig. 14. generally grid of a simulation does not have to be symmetrical (i.e. voxels are not cubes but they are in the form of general cuboids). This means that each element may be representing different volume. We have our matrix of zeros and ones which we generated in Fig. 13. Now we need to apply this matrix to the actual coordinate system.

We can use built-in function meshgrid to produce three dimensional matrices which contain coordinates for specified element (for more information on function meshgrid please refer to the Product Help of Matlab). Now we can easily determine exact coordinates of elements occupied by our model (see Fig. 15.).

**Figure 15.** Determination of Coordinates of Model Elements

At this moment we need to get better understanding of how models are defined. Since the values (in this case ones) are located in the nodes of the grid not in the centers we need to determine how the model actually looks like. See the following figure for better understanding (2D-plane is used to illustrate).

As you can see there will be an error caused by this node to cell transformation. When the simulation grid is defined appropriately (i.e. it is fine enough) the error will be only marginal. If such error is unacceptable more advanced techniques of node to cell transformation should be used. But for purposes of this example this method is more than sufficient. For example, volume of previously defined ellipsoid determined voxel by voxel is 8.5302e-6 m3. This means that the relative error is 1.8 % (and this error includes error caused by voxelization itself – this means that node to cell transformation really brings only marginal error in simulations with fine grid).

In the Fig. 16. you can see that node to cell transformation can be done in a few (precisely 8) ways. We can easily determine which way is the most precise (i.e. [(A+1) – A] or [(A-1) – A] and its combinations with B and C). Since our model is voxelized symmetrically there is almost no difference between volumes computed accordingly to various node to cell transformations (ranging between 8.5302e-6 and 8.5304e-6).

**Figure 16.** Actual Model Dimensions

34 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

we need to apply this matrix to the actual coordinate system.

of elements occupied by our model (see Fig. 15.).

**Figure 15.** Determination of Coordinates of Model Elements

understanding (2D-plane is used to illustrate).

marginal error in simulations with fine grid).

different volume. We have our matrix of zeros and ones which we generated in Fig. 13. Now

We can use built-in function meshgrid to produce three dimensional matrices which contain coordinates for specified element (for more information on function meshgrid please refer to the Product Help of Matlab). Now we can easily determine exact coordinates

At this moment we need to get better understanding of how models are defined. Since the values (in this case ones) are located in the nodes of the grid not in the centers we need to determine how the model actually looks like. See the following figure for better

As you can see there will be an error caused by this node to cell transformation. When the simulation grid is defined appropriately (i.e. it is fine enough) the error will be only marginal. If such error is unacceptable more advanced techniques of node to cell transformation should be used. But for purposes of this example this method is more than sufficient. For example, volume of previously defined ellipsoid determined voxel by voxel is 8.5302e-6 m3. This means that the relative error is 1.8 % (and this error includes error caused by voxelization itself – this means that node to cell transformation really brings only

In the Fig. 16. you can see that node to cell transformation can be done in a few (precisely 8) ways. We can easily determine which way is the most precise (i.e. [(A+1) – A] or [(A-1) – A] and its combinations with B and C). Since our model is voxelized symmetrically there is
