Stroke Monitoring

Abstract

A computer-implemented process for continuous monitoring of a brain stroke during a critical rehabilitation period, the process including the steps of: (i) accessing initial image data representing an initial image of a subjects brain containing a stroke region; (ii) accessing scattering parameter data representing microwaves scattered by the subjects brain and originating from a plurality of antennas disposed around the subjects brain; and (iii) processing the scattering parameter data and the initial image data using a gradient-free optimisation method to generate estimates of spatial dimensions of the stroke region within the subjects brain, wherein the initial image of the subjects brain is used as a priori information to improve the accuracy of the generated estimates, and the spatial dimensions of the stroke region are generated as global parameters of the gradient-free optimisation method.

Claims

1. A computer-implemented process for continuous monitoring of a brain stroke during a critical rehabilitation period, the process including the steps of: (i) accessing initial image data representing an initial image of a subject's brain containing a stroke region; (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed around the subject's brain; and (iii) processing the scattering parameter data and the initial image data using a gradient-free optimisation method to generate estimates of spatial dimensions of the stroke region within the subject's brain, wherein the initial image of the subject's brain is used as a priori information to improve the accuracy of the generated estimates, and the spatial dimensions of the stroke region are generated as global parameters of the gradient-free optimisation method.

2. The process of claim 1, wherein the spatial dimensions of the stroke region are initially determined by optimising the spatial dimensions of a first predetermined permittivity value of the stroke region and a second predetermined permittivity value for non-stroke regions of the subject's brain.

3. The process of claim 1, wherein the shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimensions of the stroke region are determined by determining the spatial dimensions of the overlapping ellipses.

4. The process of claim 3, wherein the overlapping ellipses have minor axes with fixed spatial dimensions, and the spatial dimensions of the overlapping ellipses are determined as two parameters corresponding to major axes of the overlapping ellipses.

5. The process of claim 1, wherein the spatial dimensions of the stroke region are determined by determining four geometrical parameters.

6. The process of claim 1, including repeating steps (ii) and (iii) at successive times to monitor growth or shrinkage of the stroke region over time.

7. The process of claim 1, wherein the gradient-free optimisation method is a Nelder-Mead gradient-free optimisation method.

8. The process of claim 1, wherein the spatial dimensions and relative permittivity of the stroke region are generated as global parameters of the gradient-free optimisation method.

9. The process of claim 1, wherein the initial image of the subject's brain is generated by magnetic resonance imaging or x-ray imaging or electromagnetic tomography imaging.

10. An apparatus for continuous monitoring of a brain stroke during a critical rehabilitation period, the apparatus including: a memory; at least one processor; and at least one computer-readable storage medium having stored thereon instructions that, when executed by the at least one processor, cause the at least one processor to execute the steps of: (i) accessing initial image data representing an initial image of a subject's brain containing a stroke region; (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed around the subject's brain; and (iii) processing the scattering parameter data and the initial image data to estimate spatial dimensions of the stroke region within the subject's brain, wherein the initial image of the subject's brain is used as a priori information to improve the accuracy of the determination, and the spatial dimensions of the stroke region are determined as global parameters of a gradient-free optimisation method.

11. The apparatus of claim 10, wherein the spatial dimensions of the stroke region are initially determined by optimising the spatial dimensions of a first predetermined permittivity value of the stroke region and a second predetermined permittivity value for non-stroke regions of the subject's brain.

12. The apparatus of claim 10, wherein the shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimensions of the stroke region are determined by determining the spatial dimensions of the overlapping ellipses.

13. The apparatus of claim 12, wherein the spatial dimensions of each of the overlapping ellipses are determined as two geometrical parameters.

14. The apparatus of claim 10, wherein the spatial dimensions of the stroke region are determined by determining four geometrical parameters.

15. The apparatus of claim 10, including repeating steps (ii) and (iii) at successive times to monitor growth or shrinkage of the stroke region over time.

16. The apparatus of claim 10, wherein the gradient-free optimisation method is a Nelder-Mead gradient-free optimisation method.

17. The apparatus of claim 10, wherein the spatial dimensions and relative permittivity of the stroke region are generated as global parameters of the gradient-free optimisation method.

18. The apparatus of claim 10, wherein the initial image of the subject's brain is generated by magnetic resonance imaging or x-ray imaging or electromagnetic tomography imaging.

19. At least one computer-readable storage medium having stored thereon instructions that, when executed by at least one processor of a brain monitoring apparatus, cause the at least one processor to execute the steps of claim 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Some embodiments of the present invention are hereinafter described, by way of example only, with reference to the accompanying drawings, wherein:

[0040] FIG. 1A is derived from a cross-sectional plan view MRI image of a subject's head containing a stroke region, and shows the stroke model information utilized to initiate EMT, with the corresponding tissue dielectric properties (ε.sub.r, σ) at 1.5 GHz taken from Gabriel;

[0041] FIGS. 1B and 1C show two corresponding possible critical responses of the stroke of FIG. 1A during CRP as per Schellinger;

[0042] FIG. 1D is an image of the subject's head generated by a prior art EMT process based on robust gradient-based optimization and without being initiated by the stroke-MRI image, and is clearly not able to accurately identify the stroke region;

[0043] FIGS. 1E and 1F are corresponding EMT images generated by the same process as for FIG. 1D, but initiated by the stroke-MRI image of FIG. 1A, FIG. 1E accurately identifying the stroke region, and FIG. 1F showing its subsequent response to medication (thrombolysis) at every pixel. Nevertheless, this improved EMT retrieval process takes more than an hour, making the approach impractical for stroke monitoring during CRP;

[0044] FIG. 2 shows the different optimization operations of Nelder-Mead (“NM”) optimisation represented as geometrical operations in a 3D space of parameters for the three stroke parameters a.sub.0, b.sub.0, and ε.sub.r (the geometric parameters a.sub.0 and b.sub.0 being shown in FIG. 1A);

[0045] FIG. 3 is a flowchart of an NM optimization process in accordance with the described embodiments of the present invention; after constructing the vertices, the process firstly checks whether the “best” vertex satisfies the truncation condition or not. If not, the first error reduction operation i.e. “reflection”, is performed. Accordingly, if the reflection vertex has a lower error with respect to the best vertex, then the expansion operation is performed with the hope of finding an even better vertex. If not, depending on the error level of the reflection vertex with respect to the worst and second worst vertices, forward/backward contractions and shrinking are performed to finally find a new vertex that is, at minimum, better (has a lower error) than the “worst” vertex in each iteration;

[0046] FIGS. 4A to 4C are respective graphs representing the retrieval process and convergence behavior of the NM optimization process described herein;

[0047] FIGS. 5A to 5D are respective graphs showing the effects of some influential factors on the accuracy and computational time of NM optimization in stroke monitoring during the CRP with 8 antennas, and using the two stroke shape parameters a.sub.0, b.sub.0; FIG. 5A shows the effect of the total number of antennas when SNR=15 dB, and FIG. 5B shows the effect of SNR, FIG. 5C shows the effect of stroke size when SNR=15 dB, and FIG. 5D shows the effect of difference between each individual tissue and the database in Gabriel when SNR=15 dB;

[0048] FIG. 6: Retrieval process of NM optimization for a 3D MRI-derived stroke mode, wherein:

[0049] FIG. 6A is a representation of the three-dimensional EMT geometry with a circular array of antennas disposed around the patient's head;

[0050] FIG. 6B shows a close up view of a stroke region within the patient's head as imaged by the three-dimensional EMT process described herein;

[0051] FIG. 6C shows the imaged stroke region with four geometrical parameters that define the three-dimensional stroke region;

[0052] FIG. 6D shows the imaged three-dimensional stroke region during clot breakdown;

[0053] FIG. 6E is a graph of the four geometrical stroke region parameters and the relative permittivity determined by the process as a function of iteration number;

[0054] FIG. 6F is a graph showing convergence of the process in dB as a function of iteration number;

[0055] FIG. 7 is a schematic diagram of a stroke monitoring apparatus in accordance with an embodiment of the present invention; and

[0056] FIG. 8 is a flow diagram of a stroke monitoring process in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION

[0057] The potential of EMT as a monitoring modality motivated the inventors to develop a new EMT process that is suitable for stroke monitoring during the CRP. Indeed, while some of the intrinsic limitations of EMT described above, namely the diffraction effect and the presence of evanescent waves at UHF and S bands, are unavoidable, the inventors determined that the prohibitively long computational time of prior art EMT processes is due to the numerical formulation of the gradient-based optimizations utilized in every prior art EMT system. As these optimizations find the optimum values of “variables” at every pixel of the resulting images, a high computational time is typically required to find these optimum values. To put it another way, if the spatial distribution and dielectric properties of the stroke region within the brain are unknown, then they are considered as variables whose optimum values must be retrieved at every “pixel” of the image (such as that shown in FIG. 1A). Since having a clinically acceptable spatial resolution requires a relatively high number of pixels (even when considering the limited maximum achievable spatial resolution of EMT), prior art EMT computational processes are prohibitively inefficient for stroke monitoring during the CRP. Moreover, the resulting images calculated by these inefficient processes are insufficiently accurate in any case, as is apparent from FIG. 1D.

[0058] With this in mind, the inventors identified that, as the patient is immediately transferred to the imaging unit upon arrival, the initial spatial distribution of the stroke region will be available as a priori information from the early diagnosis stage using a high-resolution imaging modality such as MRI (or X-ray or even EMT), as shown in FIG. 1A, before transferring the patient to the ICU (it takes less than an hour to generate the initial image of the stroke region using the stroke-MRI imaging modality). Although MRI (or X-ray or EMT) imaging cannot practically be used for continuous monitoring for the reasons described above (safety, cost and computational time), the stroke-MRI (or X-ray or EMT) image generated upon patient arrival can be aligned or ‘registered’ as described in G. Boverman, C. E. L. Davis, S. D. Geimer and P. M. Meaney, Image registration for microwave tomography of the breast using priors from nonsimultaneous previous magnetic resonance images, IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, vol. 2, no. 1, pp. 2-9, March 2018, and in R. L. Leijsen, W. M. Brink, C. A. T. van den Berg, A. G. Webb and R. F. Remis, Three-dimensional contrast source inversion-electrical properties tomography, IEEE Trans. Med. Imag., vol. 37, no. 9, September 2018 (“Leijsen”), where the dielectric properties of the different tissues are given in Leijsen.

[0059] The registered image (of FIG. 1A in this example) can be adopted in an otherwise conventional EMT apparatus installed in the ICU (indeed, the MRI is performed only once), and then used as a priori information to generate images of the stroke responses (shown schematically in FIGS. 1B and 1C), as shown in FIGS. 1E and 1F, demonstrating a significant improvement in retrieval accuracy. The reason for the improvement is that the registered Stroke-MRI image is used to confine the space of possible retrievable images generated by the gradient-based optimizations, leading the EMT to converge to a more accurate image. However, due to the pixel-based optimization nature of conventional (and thus gradient-based) EMT, the computational time of the method is still prohibitively high, around one hour for this example. This is just a slight reduction in time compared to the example of FIG. 1D, where no MRI image supplements the EMT (76 minutes). Therefore, gradient-based EMT might only be useful for non-emergency scenarios such as tumor screening, if initiated by MRI, X-Ray or CT-Scan images (for 3D problems, the reported computational time is around 11 hours on a standard computer, as described in Leijsen).

[0060] Described herein are an apparatus and process for continuous monitoring of a brain stroke during a critical rehabilitation period (CRP), also referred to herein for convenience as a stroke monitoring apparatus and process. In the stroke monitoring process and apparatus described herein, the MRI (or X-ray or EMT) registered image is implemented as a priori information to provide the initial geometrical shape and dielectric properties of the stroke region. However, instead of the inefficient prior art gradient-based EMT processes that calculate variables at every pixel of the imaged region, the shape, dimensions, and dielectric properties of the stroke region are calculated as global parameters. Accordingly, the global values of these parameters are continuously updated at later monitoring times by a gradient-free optimization process, as described below.

[0061] For example, in the described embodiments the parameters are the two geometrical ones defining the shape and dimensions of the stroke region as the semi-major axes of two overlapping ellipses with fixed minor axes of 1 cm, as shown in FIG. 1A, i.e. a.sub.0 and b.sub.0 (with initial values of a.sub.0=1.5 cm and b.sub.0=2 cm), and a third parameter defining the relative permittivity i.e. ε.sub.r of the stroke region. Treating the shape and dielectric properties of the stroke as global parameters reduces the number of unknown quantities from, say, around 1000 pixel-based variables for the image in FIG. 1A, to only 3 global parameters. This reformulation enables gradient-free optimizations to be used for EMT as highly efficient techniques for parameter-optimization (since parameters are optimized globally, not locally at each pixel). Clearly, such gradient-free optimisations are unable to provide pixel-based images such as those of FIGS. 1D to 1F. However, and regarding the medical needs during CRP, rapid and accurate tracking of stroke parameters (e.g., a.sub.0, b.sub.0 and ε.sub.r) every two minutes is much more medically relevant than obtaining a coarse image of stroke on an hourly time-basis by gradient-based optimization.

[0062] Gradient-free optimization methods were developed to solve problems for which gradient-based optimizations were not applicable, in particular, when the function to be minimized during optimization is not differentiable or smooth.

[0063] Various gradient-free optimization methods have been developed for different electromagnetic and antenna applications, including the Nelder-Mead (“NM”), genetic algorithm (“GA”), and particle swarm optimization (“PSO”) methods. Where the computational time is the main concern in the optimization procedure apart from accuracy, the NM optimization method is usually the fastest gradient-free optimization method. In the context of stroke monitoring during CRP, where time is life, the inventors consider that NM optimization best suits the monitoring requirements (although other gradient-free optimization method may be used in other embodiments). As NM optimization has not been previously used for EMT applications, an overview of the NM optimization process for this particular application is described below. The general NM methodology is described in N. Pham, A. Malinowski and T. Bartczak, “Comparative study of derivative free optimization algorithms,” IEEE Trans. Industr. Inform., vol. 7, no. 4, pp. 592-600, November 2011 (“Pham”).

[0064] In the described embodiments, the described processes are executed by a stroke monitoring apparatus, as shown in FIG. 7. In use, the apparatus is in communication with a vector network analyser (VNA) or transceiver 701 that is in turn connected to an array of microwave antennas 705.

[0065] The array of microwave antennas 701 is arranged to receive the head 704 of a patient whose brain is to be imaged, as shown, so that each antenna of the array can be selectively energised to radiate electromagnetic waves or signals of microwave frequency into and through the subject's head to be scattered and the corresponding scattered signals detected by all of the antennas of the array, including the antenna that transmitted the corresponding signal.

[0066] As will be apparent to those skilled in the art, the vector network analyser (VNA) 701 energises the antennas as described above, and records the corresponding signals from the antennas as data (referred to herein as ‘scattering’ data) representing the amplitudes and phases of the scattered microwaves in a form that is known in the art as “scattering parameters” or “S-parameters”. The VNA 701 sends this data to the apparatus for processing to generate information on internal features of the imaged object (e.g., brain clots, bleeding sites, and other features). In the described embodiments, a VNA which has a large dynamic range of more than 700 dB and a noise floor below −700 dBm, can be used to activate the antennas to transmit electromagnetic signals across the frequency band of 0.5 to 4 GHz and receive the scattered signals from those antennas.

[0067] Although the apparatus of the described embodiments is in the form of a computer, this need not be the case in other embodiments. As shown in FIG. 7, the stroke monitoring apparatus of the described embodiments is a 64-bit Intel Architecture computer system, and the stroke monitoring processes executed by the stroke monitoring apparatus are implemented as programming instructions of one or more software modules 702 stored on non-volatile (e.g., hard disk or solid-state drive) storage 704 associated with the computer system. However, it will be apparent that at least parts of these processes could alternatively be implemented in one or more other forms, for example as configuration data of a field-programmable gate array (FPGA), or as one or more dedicated hardware components, such as application-specific integrated circuits (ASICs), or as any combination of such forms.

[0068] The stroke monitoring apparatus includes random access memory (RAM) 706, at least one processor 708, and external interfaces 710, 712, 713, 714, all interconnected by a bus 716. The external interfaces include a network interface connector (NIC) 712 which connects the stroke monitoring apparatus to a communications network such as the Internet 720, and universal serial bus (USB) interfaces 710, at least one of which may be connected to a keyboard 718 and a pointing device such as a mouse 719, and a display adapter 714, which may be connected to a display device such as an LCD panel display 722.

[0069] The stroke monitoring apparatus also includes an operating system 724 such as Linux or Microsoft Windows, and in some embodiments includes additional software modules 726 to 730, including web server software 726 such as Apache, available at http://www.apache.org, scripting language support 728 such as PHP, available at http://www.php.net, or Microsoft ASP, and structured query language (SQL) support 730 such as MySQL, available from http://www.mysql.com, which allows data to be stored in and retrieved from an SQL database 732.

[0070] Together, the web server 726, scripting language module 728, and SQL module 730 provide the stroke monitoring apparatus with the general ability to allow remote users with standard computing devices equipped with standard web browser software to access the stroke monitoring apparatus and in particular to monitor the progress of a stroke during the CRP.

[0071] A. NM Optimization for EMT: Initialization

[0072] For the EMT problems shown in FIGS. 1B and 1C (initiated by FIG. 1A), eight imaging dipole antennas are arranged in a circular configuration surrounding the head of a subject so that the 2D MRI model of the subject's head crosses the ports of the dipole antennas (because the S-matrix that is implemented in the optimization process is recorded at the dipole ports).

[0073] The antennas illuminate the head with a sinusoidal electromagnetic wave at 1.5 GHz, say five and a half hours after symptoms onset, when thrombolysis has already been applied. This timing example is taken from a stroke case described in Schellinger, where the subject's arrival time is around 3 hours after symptom onset, and the 2D single-slice stroke-MRI is prepared 0.75 hours after arrival. The antennas then record the corresponding scattered fields in the form of an S-matrix (the effects of the number of antennas on the accuracy and retrieval time are discussed below). This S-matrix is then implemented to update the three stroke parameters a.sub.0, b.sub.0, ε.sub.r whose initial values in FIG. 1A are X.sub.0=(a.sub.0=1.5 cm, b.sub.0=2 cm, ε.sub.r=39), which can be considered to define a point (or vertex as described below) in 3D space (as the number of parameters is 3).

[0074] These parameters were selected for the following reasons. Since the effect of any medication applied to a stroke is reflected by the change in the relative permittivity ε.sub.r of the stroke region of the subject's brain, it can be considered as a reliable parameter to monitor the stroke and converge to a value that matches the S-matrix at later times. As described above, since the retrieval accuracy of the conductivity σ of the stroke region is usually poor by comparison with the relative permittivity ε.sub.r, the conductivity is excluded from retrieval. Moreover, the geometrical or shape parameters a.sub.0, b.sub.0 are chosen to best reflect the geometrical change (extension or shrinkage) of the stroke region. For the smooth stroke shape shown in FIG. 1A, a.sub.0, b.sub.0 are the semi-major axes of two ellipses that are considered to most accurately represent the geometrical changes of the stroke region with respect to possible alternative geometrical parameters (more complicated geometries are discussed below). The inventors have found that these shape parameters also continuously converge to values that best fit the corresponding S-matrix. The 2D physical problem domain in FIG. 1A is different from the 3D mathematical space of parameters constructed in FIG. 2A.

[0075] In addition to X.sub.0, the variation range of the three parameters defines the available parameter space, assuming that the stroke region subject to hypoperfusion (FIG. 1B) can be extended across the entire left hemisphere of the brain, and that successful clot breakdown (FIG. 1C) can entirely remove the clot (i.e., shrink the stroke region to nothing). This large variation range corresponds to a parameters space of 0≤a.sub.0≤4 cm, 0≤b.sub.0≤7 cm, and 39≤ε.sub.r≤43.5, indicating that the optimisation process does not require very accurate parameter ranging in finding the optimum vertex, even if there is no solid prediction about the stroke response to medication. To optimize these 3 parameters, the NM optimization first constructs an equal-length simplex (generalized triangle) with 3+1 vertices within the 3D space of parameters, as shown in FIG. 2A. One of the vertices is, indeed, the initial guess (i.e., starting point) X.sub.0 provided by the initial MRI image. All other vertices with equal distance c are derived by adding the following vectors to this initial guess, as follows:

X.sub.1=X.sub.0+(p, q, q)

X.sub.2=X.sub.0+(q, p, q)

X.sub.3=X.sub.0+(q, q, p)   (1)

where the vector components are,

[00001]$\begin{array}{cc}q=\frac{c}{N\sqrt{2}}\left(\sqrt{N+1}-1\right)p=q+\frac{c}{\sqrt{2}}& \left(2\right)\end{array}$

and where N is the number of parameters (i.e., 3 in the described embodiments). Typically, c=1 to allow the process to search in a sufficiently large volume at the initial step. Small values for c normally require a long computational time and can mislead the process to find only a local minimum. Moreover, as the presence of noise can cause slightly different vertices to result in the same frequency response (S-matrix), locating the vertices far enough from one another (distanced with c≥1) makes the process robust with respect to noise at early iterations.

[0076] By constructing the simplex, the next step is to evaluate an objective function (i.e., the function to be minimized by the NM optimization) at every vertex X.sub.0, X.sub.1, X.sub.2, X.sub.3. In the EMT problem of FIG. 1B or 1C, this objective function is defined as the L.sub.2 norm mismatch between the measured (“meas.”) S-matrix and the retrieved (“retr.”) S-matrix, as follows:

[00002]$\begin{array}{cc}F\left(a_0,b_0,{\mathrm{.Math.}}_r\right)={.Math.S^{\mathrm{meas}.}-S^{\mathrm{retr}.}.Math.}^2={.Math.S^{\mathrm{meas}.}-k_0^2\underset{\Omega }{\int }G\left(a_0,b_0,{\mathrm{.Math.}}_r\right).Math.\left(\chi \left(a_0,b_0,{\mathrm{.Math.}}_r\right)E\left(a_0,b_0,{\mathrm{.Math.}}_r\right)\right)d\Omega .Math.}^2& \left(3\right)\end{array}$

where Ω denotes the imaged domain (in cylindrical coordinates â.sub.ρ, â.sub.φ, â.sub.z), χ is the contrast in the dielectric properties of the human head defined as

[00003]$\chi =\left(\frac{\mathrm{.Math.}}{{\mathrm{.Math.}}_0}-1\right)$

where

[00004]$\mathrm{.Math.}={\mathrm{.Math.}}_0{\mathrm{.Math.}}_r-j\frac{\sigma }{\omega }$

is the complex permittivity,

[00005]${\mathrm{.Math.}}_0=\frac{1}{36\pi }\times {10}^{-9}$

is the free-space permittivity, and ω=2π×1.5 GHz is the angular frequency, respectively. In the described example with eight antennas, the size of the S-matrix is therefore 8×8; the mismatch is thus the difference between the corresponding matrix elements of S.sup.meas. and S.sup.retr.. E is the total electric field across the imaged domain of FIG. 1B or 1C, derived as described in Chapter 9 of W. C. Chew, Waves and Fields in Inhomogeneous Media, IEEE Press, New York, 1995 (“Chew”) and A. Afsari, A. Abbosh, and Y. Rahmat-Samii, “Modified Born iterative method in medical electromagnetic tomography using magnetic field fluctuation contrast source operator,” IEEE Trans. Microw. Theory Techn., DOI: 10.1109/TMTT.2018.2876228, and given by:

[00006]$\begin{array}{cc}E\left(a_0,b_0,{\mathrm{.Math.}}_r\right)=E^{\mathrm{inc}}-k_0^2\underset{\Omega }{\int }G\left(a_0,b_0,{\mathrm{.Math.}}_r\right).Math.\left(\chi E\left(a_0,b_0,{\mathrm{.Math.}}_r\right)\right)d\Omega & \left(4\right)\end{array}$

where E.sup.inc is the incident electric field in absence of any object in Ω. Finally, G is the dyadic Green's function given in Chapter 1 of Chew, and k.sub.0=ω√{square root over (μ.sub.0ε.sub.0)} is the free-space wavenumber, wherein μ.sub.0=4π×10.sup.−7 is the free-space permeability.

[0077] After evaluating the objective function for all of the vertices, three of the vertices possess special importance in the NM technique, as shown in FIG. 2B. The vertex which gives the smallest value of the objective function (given by Equation (3) in the described embodiments) is referred to as the “best” vertex, and is denoted as X.sub.b. The vertex with highest L.sub.2 norm mismatch is referred to as the “worst” vertex, X.sub.w. In order to guide the process into a direction that always reduces the mismatch in Equation (3), the process determines an auxiliary point for which Equation (3) has the second highest value. This vertex is referred to as the “second worst” vertex, X.sub.sw. The average value of all of the vertices except X.sub.w, is then calculated as X.sub.a. The line segment between X.sub.a and X.sub.w (L.sub.XwXa) is always downward, and contains some useful points wherein Equation (3) may have a lower mismatch than X.sub.w during the optimization processes. As per this downward-direction search for the optimum point, the process is also referred to as the “downhill” optimization.

[0078] After the above initialization, the process performs at least two, and at most five “error-reduction” operations. In each of these steps, the old value of x.sub.w is removed from the computer memory (i.e., is not stored for the next iteration), and all the other vertices are rearranged to provide new values for x.sub.b, x.sub.w, x.sub.sw. Hereafter, the following operations are introduced:

[0079] Reflection: The first optimization step in the NM approach is to reflect the worst vertex x across L, with the same length, as follows:

X.sub.r=2X.sub.a−X.sub.w   (5)

as shown in FIG. 2c. This operation is to check whether moving in the L.sub.X.sub.w.sub.X.sub.a direction should be continued by the process or another direction will lead the process to the optimum vertex. If Equation (3) has a lower mismatch at X.sub.r with respect to X.sub.b, that is F (X.sub.r)<F (X.sub.b), before replacing

X.sub.b←X.sub.r

X←X.sub.b

X.sub.sw←X

X.sub.w←X.sub.sw   (6)

the process evaluates the chance of finding even a better vertex (where the programming convention A←B represents that the old value A is substituted by the new value B). To this end, the expansion operation is always performed by further moving in the same direction L.sub.X.sub.x.sub.X.sub.a.

[0080] Expansion: As per FIG. 2d, the process further moves along L.sub.X.sub.w.sub.X.sub.a with the same step-length i.e.

X.sub.e=2X.sub.r−X.sub.a   (7)

[0081] Then, Equation (3) is also evaluated at this expansion vertex. If its value is lower than X.sub.b (even if it is worse than X.sub.r) i.e. F(X.sub.e)<F(X.sub.b), the process replaces

X.sub.b←X.sub.e

X←X.sub.b

X.sub.sw←X

X.sub.w←X.sub.sw   (8)

and (iteratively) returns to the reflection step. The reason that the process does not immediately accept X.sub.b←X.sub.r, despite it being the best-found vertex among the other vertices, comes from the fact that this vertex is reserved by the process, as it lies inside the new simplex formed by X.sub.e. Hence, by performing the expansion, the neighborhood domain of X.sub.r is merely safeguarded as the subdomain wherein some other good or even better vertices may exist to minimize Equation (3). Nevertheless, if Equation (3) at X.sub.e is not lower than X.sub.b, then the substitutions of Equation (6) are performed and the process iteratively returns to the first operation (i.e., reflection).

[0082] Forward Contraction: Either Equation (6) or Equation (8) assumes F(X.sub.r)<F(X.sub.b). If this is not realized, but F(X.sub.sw)<F(X.sub.r)<F(X.sub.w), the process has excessively moved along the L.sub.X.sub.w.sub.X.sub.a direction, and better vertices may lie at distances closer than X.sub.r. Accordingly, the forward contraction along L.sub.X.sub.w.sub.X.sub.a is performed by returning back half a step-length (typically the balanced step-length) from X.sub.r toward X.sub.a as shown in FIG. 2E, i.e.

X.sub.fc=1.5X.sub.a−0.5X.sub.w   (9)

[0083] If F(X.sub.fc)<F(X.sub.r), a new simplex is formed on vertices X, X.sub.b, X.sub.sw, X.sub.fc by returning to the initialization step and rearranging these vertices from the worst to the best one.

[0084] Backward Contraction: If F(X.sub.w)<F(X.sub.r), then L.sub.X.sub.w.sub.X.sub.a may still contain some vertices that can improve the mismatch in Equation (3), and in the same way as forward contraction, a backward contraction is performed by moving backward half a step-size from X.sub.a toward X.sub.w, as shown in FIG. 2f, as follows:

X.sub.bc=0.5X.sub.a+0.5X.sub.w   (10)

[0085] If F(X.sub.bc)<F(X.sub.w), then a new simplex as shown in FIG. 2f is constructed on X, X.sub.b, X.sub.sw, X.sub.bc by returning to the initialization step and rearranging the vertices from the worst to the best.

[0086] Shrinking: If, nonetheless, none of the above conditions takes place, then the last step to find a better direction toward the optimum vertex is to shrink the simplex. To this end, only the best vertex X.sub.b is kept, and for the other vertices, the shrinking operation is performed as follows (for each ith vertex):

X.sub.i(new)=0.5X.sub.b+0.5X.sub.i(old)   (11)

[0087] Then, the process returns to the initialization step to rearrange the new vertices formed in the shrinking step shown in FIG. 2g. This iterative process of NM optimization is continued until F(X.sub.b) meets the truncation condition of the objective function, in the described embodiment being:

F(X.sub.b)<10.sup.−7   (12)

[0088] The value of X.sub.b that satisfies the truncation condition is stored as the final result. The truncation condition in Equation (12) is chosen to be very small, so as to ensure that the required accuracy in retrieving the parameters is satisfied. Larger values of truncation conditions do not lead to very accurate parameter retrieval. To demonstrate all these steps at once, FIG. 3 is a flowchart of the NM optimization process for EMT, and corresponding pseudocode is given below.

TABLE-US-00001 Steps Commands 1: Input: introduce parameters a.sub.0, b.sub.0, ε.sub.r Input: parameters ' range 0 ≤ a.sub.0 ≤ 4 cm,0 ≤ b.sub.0 ≤ 7 cm,39 ≤ ε.sub.r ≤ 43.5 Input: X.sub.0 Input: Measured S-parameters of Fig. 1B or 1C 2: Do: vertex construction: eq. (1) Do: vertex rearrangement: Fig. 3A Do: calculation of X.sub.a 3: if (12) is satisfied Output: X.sub.b else Do Reflection if F(X.sub.r ) < F(X.sub.b ) Do Expansion if F(X.sub.e ) < F(X.sub.b ) Perform (8) by returning to Vertex Rearrangement else Perform (6) by returning to Vertex Rearrangement end else if F (X.sub.w ) < F (X.sub.r ) Do Backward Contraction if F (X.sub.bc ) < F (X.sub.w ) Accept X.sub.bc and Return to Vertex Rearrangement else Shrink and Return to Vertex Rearrangement end else if F (X.sub.sw ) < F(X.sub.r ) < F(X.sub.w ) Do Forward Contraction if F(X.sub.fc ) < F(X.sub.r ) Accept X .sub.fc and Return to Vertex Rearrangement else Shrink and Return to Vertex Rearrangement end else Accept X.sub.r and Return to Vertex Rearrangement end end end end

[0089] I. NM Gradient-Free Optimization in Practice: 2D Retrieval

[0090] To efficiently monitor the different stroke behaviours known as hypoperfusion and shown in FIG. 1B or the clot-breakdown in FIG. 1C, the NM optimization process described above is utilized to optimize the three parameters a.sub.0, b.sub.0, and ε.sub.r. The reason for merely choosing the shape parameters a.sub.0 and b.sub.0 for optimization and excluding location parameters such as the coordinate of the center of the stroke in FIG. 1B or 1C comes from the fact that in either hypoperfusion or clot-breakdown process, the stroke location does not change. Instead, the stroke merely extends or shrinks around its center.

[0091] A. NM Optimization Performance

[0092] For hypoperfusion, the actual values of the parameters in the given example are a.sub.0=2.5 cm, b.sub.0=3 cm, ε.sub.r=39, and the retrieved ones, whose evolution over successive iterations are shown in FIG. 4A, are a.sub.0=2.68 cm, b.sub.0=3.32 cm, ε.sub.r=39.08. As FIG. 5A shows, the NM process variations gradually become stable when the simplex becomes smaller during iterations. For the clot-breakdown of FIG. 1C, the actual and retrieved values of the NM optimization process are a.sub.0=1 cm, b.sub.0=1 cm, ε.sub.r=40, and a.sub.0=0.997cm, b.sub.0=1.002 cm, ε.sub.r=40.07 (FIG. 4B), respectively. It is to be noted that the relative permittivity of the affected area, during clot-breakdown, is again approaching that of the brain, when being recovered. In both cases, although the simulation environment is highly noisy (the SNR of the S-matrix is set to 15 dB), the parameters are nevertheless accurately updated. To provide a quantitative analysis of this accuracy, a “retrieval error” is defined as:

[00007]$\begin{array}{cc}10\log \left(\frac{.Math.\underset{i=1}{\overset{3}{.Math.}}({\mathrm{parameter}}_i^{\mathrm{actual}}-{\mathrm{parameter}}_i^{\mathrm{retrieved}}.Math.}{.Math.\underset{i=1}{\overset{3}{.Math.}}{\mathrm{parameter}}_i^{\mathrm{actual}}.Math.}\right)& \left(13\right)\end{array}$

[0093] The retrieval error is low in each case: −18.86 dB for the hypoperfusion outcome, and −27.95 dB for the clot-breakdown outcome.

[0094] The critical point, however, is the retrieval time. As per FIG. 4C, the process terminates once the termination condition of Equation (12) is satisfied. In the described examples, the hypoperfusion and clot-breakdown retrievals require 94 and 69 iterations, respectively, and each iteration takes 1.25 seconds. Hence, within every two minute time-frame, the stroke is continuously monitored to evaluate the treatment process and make instant preparations for craniotomy if required (the surgical removal of a part of the subject's skull to directly access the affected area).

[0095] FIGS. 4A and 4B provide a graphical illustration of the performance of the NM process respectively over the first 70 and 40 iterations only (the variations of the retrieved parameters in further iterations being too small to be usefully plotted in the same graph). The greater geometrical change caused by hypoperfusion with respect to the initial topology (i.e. FIG. 1A) in comparison with the clot-breakdown process results in a higher retrieval error and a longer convergence time when compared to the clot-breakdown retrieval.

[0096] A. Influential Factors on Accuracy and Computational Time

[0097] Among the different factors affecting the accuracy (retrieval error) and computational time of EMT based on NM optimization, the major contributors are: the total number of imaging antennas, the SNR, the shape parameters a.sub.0, b.sub.0, and the minor difference in the dielectric properties of each subject's head tissues with respect to the database described in C. Gabriel, S. Gabriel and E. Corthout, “The dielectric properties of biological tissues: I. Literature survey”, Phys. Med. Biol., vol. 41, no. 1, pp. 2231-2249, 1996 (“Gabriel”) utilized to register the initial stroke-MRI image. To depict these influences, as the retrieval of hypoperfusion is more challenging (having a higher retrieval error), the effects of the number of antennas, SNR, and the accuracy level of Gabriel are described for this medical scenario.

[0098] FIG. 5A shows that increasing the number of antennas up to a specific level remarkably improves the accuracy at the expense of increasing the computation time per iteration. This improvement comes from the fact that having additional antennas is equal to having more information. For the same reason, the total number of iterations reduces from 94 to 78, as the space of possible solutions is now more confined, and the process requires fewer iterations to satisfy the truncation condition. Nevertheless, due to the increase in the size of the objective function, whose size corresponds to the size of the S-matrix, the computational time per iteration drastically increases. Indeed, having M imaging antennas results in an M×M S-matrix. In the literature, 8 has been suggested as the minimum number of antennas to provide a medical image of the head with acceptable accuracy.

[0099] The second factor to consider is the SNR. When the measured data are highly contaminated by noise, the vertices that are close to each other (having close parameter values in FIG. 2A) possess very similar frequency responses that are most likely dominated by noise. Thus, the accuracy of the process is degraded as it is no longer able to discriminate between close vertices. If the medical imaged domain is highly isolated with respect to noise, both the retrieval accuracy and the required number of iterations are improved as shown in FIG. 5B. This is because the measured S-matrix in Equation (3) contains more reliable information to retrieve the stroke response (hypoperfusion), and thus, the vertex that best minimizes the objective function in Equation (3) is closer to the actual one. The computational time per iteration remains unchanged i.e. 1.25 seconds, as this time is a function of the number of imaging antennas and the retrieval parameters.

[0100] The next factor to consider is the size-range of the stroke region that can be retrieved accurately. To this end, the cross-sectional factor a.sub.0×b.sub.0, as a rule of thumb, is introduced to represent the stroke size. FIG. 5c shows that the retrieval process becomes inaccurate for the example in FIG. 1B if a.sub.0×b.sub.0<0.64 cm.sup.2 or a.sub.0×b.sub.0>9 cm.sup.2. As described above, the diffraction effect and the presence of evanescent waves are noticeable in the UHF and S bands, and these effects become highly destructive when the size of the stroke region becomes less than a quarter-wavelength. As the wavelength at 1.5 GHz is 3.2 cm inside the clot (with ε.sub.r=39), the quarter-wavelength is 0.8 cm, resulting in 0.64 cm.sup.2 cross-sectional factor, below which the retrieval accuracy decreases substantially. On the other hand, for a very high cross-sectional factor e.g., a.sub.0×b.sub.0>9 cm.sup.2, the change in the stroke shape is not accurately trackable as this change unavoidably requires additional shape parameters (besides a.sub.0, b.sub.0) to more accurately model the irregularity and asymmetry of the stroke, especially when the stroke is deformed by the skull from one side, but is still extendable from the other side. From the computational time perspective, if only two shape parameters are employed for large-size stroke regions, the process requires more iterations with respect to the stroke in FIG. 1B to partially conform the initial shape of the stroke to the extended one, as shown in FIG. 5c. This results from the fact that adding each shape parameter with its corresponding variation range will more confine the space of possible solutions. Dealing with only two shape parameters does require more iterations to partially match the irregularity and asymmetry of extended stroke (this is demonstrated below). Nevertheless, it is noted that increasing the number of stroke parameters, as a solution to improve the retrieval error and reduce the total iterations, does not necessarily reduce the computational time, as adding each parameter to the problem adds another dimension to the space of parameters. Accordingly, the computational time per iteration increases exponentially.

[0101] The last factor is the effect of the level of agreement between the database in Gabriel and the dielectric properties of each individual human head. As the gradient-free optimization process merely focuses on retrieving the global stroke parameters, it is very important to ensure that the dielectric properties in Gabriel by which the stroke-MRI is registered into the EMT apparatus are highly accurate to avoid a noticeable retrieval deviation from the desired values. Fortunately, this agreement is typically very high, as the materials constructing the tissues in FIG. 1A are the same for different individuals, and the minor difference between these materials and the database in Gabriel mainly comes from the partial statistical randomness. To take this factor into consideration, FIG. 5d applies up to 5% randomdifference to each tissue in FIG. 1A with respect to their corresponding values reported in Gabriel, to demonstrate the robustness of the process with respect to minor differences in dielectric properties of different head tissues. As seen, both the computational speed and retrieval error slightly drop, while still falling within the acceptable requirements of the CRP.

[0102] By studying the effects of influential factors on the accuracy and computational time of the NM optimization process in medical EMT application, the process is applied below to a more complicated problem where a 3D MRI-derived stroke model, as per FIG. 4 of Leijsen, supplements the NM process as shown in FIG. 6. Four shape parameters are implemented in 3D to monitor the more complicated clot breakdown process.

[0103] I. NM Gradient-Free Optimization in Practice: 3D Retrieval

[0104] For realistic 3D EMT problems, the accuracy of the retrieved parameters can be further improved if either the total electric field within the domain, or the retrieved S -matrix, is simulated using a well-developed numerical method such as finite element modelling (FEM). This can be realized by directly solving the wave equation and considering the entire three-dimensional physical structure of the imaging antennas (see FIG. 6A), rather than using the point-source-based Equations (3) and (4) given above. Indeed, the reason for implementing the point-source-based objective function of Equation (3) given above is that the same objective function is used in prior art gradient-based optimization methods in EMT. Therefore, such a point-source approximation provides a better basis for comparing the two methodologies. Nevertheless, due to the additional complexity of the realistic antenna structures of the problem shown in FIG. 6, the commercial COMSOL software package is used for the simulations. The stroke region within an entire head model is determined as a complex shape defined by triangular mesh elements, as shown in FIG. 6B, and the stroke parameters are illustrated in FIG. 6C with initial values X′.sub.0=(a.sub.1=0.41 cm, b.sub.1=0.4 cm, c.sub.1=0.41 cm, d.sub.1=0.42 cm, ε.sub.r=39), these values being the distances from the central point of the stroke region to the four furthest points on the perimeter of the region.

[0105] When thrombolytic treatment is performed in time, the clot breakdown process starts as seen in FIG. 6d. This process, slightly or significantly, changes the stroke parameters to X.sub.1=(a.sub.1=0.4 cm, b.sub.1=0.395 cm, c.sub.1=0.39 cm, d.sub.1=0.41 cm, ε.sub.r=42).

[0106] As per the graph of FIG. 6E, the retrieved stroke region shape parameters are X.sub.1=(a.sub.1=0.401 cm, b.sub.1=0.395 cm, c.sub.1=0.392 cm, d.sub.1=0.411 cm) and the retrieved permittivity of the stroke region is ε.sub.r=41.43, in excellent agreement with their actual values. As FIG. 6D shows, the changes in the spatial dimensions of the stroke region due to clot breakdown are greatest for the shape parameters c.sub.1, d.sub.1, which is why they require additional iterations to approach their actual values. FIG. 6F is a convergence plot of the NM optimization process. Due to having sufficient shape parameters, fewer than 50 iterations are required to reach the truncation condition. In COMSOL Multiphysics installed on the same PC, due to the large size of the 3D problem domain, each iteration requires 15 seconds to complete. In total, the behaviour of the stroke can be monitored every 11 minutes (i.e., the described apparatus can generate a new group of plots like FIG. 6E every 11 minutes), which is a reasonable time for 3D monitoring during the CRP. Taking the results described in Pham as the ground truth, at least 10 times more computational time is required when the GA or PSO gradient-free methods are used instead of the NM method, e.g. about 110 minutes for this specific problem. This time approaches that of prior art gradient-based optimization methods, and normally exceeds the CRP timeframe. Therefore, the inventors consider that, among the gradient-free methods and general-purpose computer hardware available at the time of writing, only EMT systems based on the NM gradient-free optimizations can be efficiently utilized as 2D or 3D stroke monitoring tools during the CRP.

[0107] The EMT monitoring process and apparatus described herein and based on Nelder-Mead gradient-free optimization provide the ability to monitor the expansion or contraction of stroke during the CRP, and can therefore potentially be considered as a translational medical advance to increase the chance of survival from stroke. The results described herein demonstrate that the process is highly efficient to retrieve a 2D stroke response within every 2 minutes, or a 3D stroke response within every 11 minutes on a general-purpose computer platform, while other gradient-free approaches such as GA or PSO can generate the same outputs but at the expense of much longer computational times. The described process can be initiated by stroke-MRI data available from early diagnosis. Then, the shape and dielectric properties (the real part of permittivity) of the stroke region are defined as global parameters. Following this, the patient can be successively imaged by a portable EMT system as described herein using a small number of imaging antennas, and the S-matrix recorded by these antennas at each imaging step is utilized to update the global parameters and thus identify the expansion or contraction of the stroke region, in particular in response to one or more treatments. The described process and apparatus can thus improve the treatment process, and consequently, the chance of survival for victims of stroke.

[0108] Many modifications will be apparent to those skilled in the art without departing from the scope of the present invention.