Titre: Forward and Inverse Modelling of Magnetic Induction Tomography Title: (MIT) for Biomedical Application Auteur: Amir Ahmad Roohi Noozadi Author: Date: 2017 Type: Mémoire ou thèse / Dissertation or Thesis Roohi Noozadi, A. A. (2017). Forward and Inverse Modelling of Magnetic Induction Référence: Tomography (MIT) for Biomedical Application [Ph.D. thesis, École Polytechnique Citation: de Montréal]. PolyPublie. https://publications.polymtl.ca/2536/ Document en libre accès dans PolyPublie Open Access document in PolyPublie URL de PolyPublie: https://publications.polymtl.ca/2536/ PolyPublie URL: Directeurs de recherche: Arthur Yelon, & David Ménard Advisors: Programme: Génie physique Program: Ce fichier a été téléchargé à partir de PolyPublie, le dépôt institutionnel de Polytechnique Montréal This file has been downloaded from PolyPublie, the institutional repository of Polytechnique Montréal https://publications.polymtl.ca (cid:19) (cid:19) UNIVERSITE DE MONTREAL FORWARD AND INVERSE MODELLING OF MAGNETIC INDUCTION TOMOGRAPHY (MIT) FOR BIOMEDICAL APPLICATION AMIR AHMAD ROOHI NOOZADI (cid:19) (cid:19) DEPARTEMENT DE GENIE PHYSIQUE (cid:19) (cid:19) ECOLE POLYTECHNIQUE DE MONTREAL (cid:18) (cid:19) (cid:19) THESE PRESENTEE EN VUE DE L’OBTENTION ^ DU DIPLOME DE PHILOSOPHI(cid:29) DOCTOR (cid:19) (GENIE PHYSIQUE) MARS 2017 ⃝c Amir Ahmad Roohi Noozadi, 2017. (cid:19) (cid:19) UNIVERSITE DE MONTREAL (cid:19) (cid:19) ECOLE POLYTECHNIQUE DE MONTREAL Cette th(cid:18)ese intitul(cid:19)ee : FORWARD AND INVERSE MODELLING OF MAGNETIC INDUCTION TOMOGRAPHY (MIT) FOR BIOMEDICAL APPLICATION pr(cid:19)esent(cid:19)ee par : ROOHI NOOZADI Amir Ahmad en vue de l’obtention du dipl^ome de : Philosophi(cid:26) Doctor a (cid:19)et(cid:19)e du^ment accept(cid:19)ee par le jury d’examen constitu(cid:19)e de : M. FRANCOEUR S(cid:19)ebastien, Ph. D., pr(cid:19)esident M. YELON Arthur, Ph. D., membre et directeur de recherche (cid:19) M. MENARD David, Ph. D., membre et codirecteur de recherche M. LEBLOND Fr(cid:19)ed(cid:19)eric, Ph. D., membre M. ADLER Andy, P. Eng., membre externe iii DEDICATION To my Mom, Dad and My dearest Ely iv ACKNOWLEDGMENT I would like to thank my advisers, Arthur Yelon and David M(cid:19)enard, for their consistent sup- port and inspiration. It was a pleasure working with these two brilliant scientists. I would also like to thank Herv(cid:19)e Gagnon for fruitful discussions, technical consultant and helpful suggestions. Herv(cid:19)e expertise had a noticeable effect on the progress rate of this re- search. I thank the members of the jury S(cid:19)ebastien Francoeur, Fr(cid:19)ed(cid:19)eric Leblond, and Andy Adler for evaluating this thesis and for their constructive comments. I would like to thank my former supervisors Reza Jafari and Shahram Jalalzadeh for their supports and fruitful discussions since I started studying physics till the day this thesis was submitted. I thank all my colleagues and friends, especially Nima Nateghi, Saeed Bohloul, Mehran Yazdizadeh, Antoine Morin, Nicolas Teyssedo, Sahar Fazeli and Saman Choubak, for their support and helpful discussions. And (cid:12)nally I thanks my wonderful Mom, Dad, and my sister Ely for their continuous support and endless love. v RE(cid:19)SUME(cid:19) Cette th(cid:18)ese d(cid:19)eveloppe un outil de simulation destin(cid:19)e au design d’instruments de tomographie par induction magn(cid:19)etique (MIT pour Magnetic Induction Tomography). Ce simulateur per- metd’investiguerlapossibilit(cid:19)ed’utiliserlam(cid:19)ethoded’imageriepartomographieparinduction magn(cid:19)etique a(cid:12)n de faire la conception d’un dispositif sans-contact capable de d(cid:19)etecter une h(cid:19)emorragie (cid:18)a l’int(cid:19)erieur du cr^ane humain. La m(cid:19)ethode sp(cid:19)eci(cid:12)que de calcul num(cid:19)erique utili- s(cid:19)ee pour la simulation du dispositif de m^eme que le calcul de la sensibilit(cid:19)e avec la m(cid:19)ethode directe (introduite dans cette th(cid:18)ese) et la r(cid:19)esolution du probl(cid:18)eme inverse, qui reconstruit une carte de la conductivit(cid:19)e (cid:18)a partir des r(cid:19)esultats de simulation, sont optimis(cid:19)es a(cid:12)n de simuler un dispositif op(cid:19)erant (cid:18)a 50 kHz. Ce dispositif est capable de d(cid:19)etecter le changement de la conductivit(cid:19)e dans une gamme se rapprochant de celle des tissus biologiques. Lefonctionnementdebasedelatomographieparinductionmagn(cid:19)etiquereposesurlesmesures des propri(cid:19)et(cid:19)es(cid:19)electromagn(cid:19)etiques dites passives telles que la conductivit(cid:19)e. L’utilisation d’une tellem(cid:19)ethodepourd(cid:19)etecterlesh(cid:19)emorragiesc(cid:19)er(cid:19)ebralessejusti(cid:12)eparlefaitquelaconductivit(cid:19)e du sang est plus(cid:19)elev(cid:19)ee que la conductivit(cid:19)e des autres tissus constituant le cerveau. Une autre application potentielle de cette m(cid:19)ethode est le suivi en temps r(cid:19)eel, de mani(cid:18)ere non invasive, de l’alt(cid:19)eration des tissus qui peut s’observer (cid:18)a partir d’un changement de conductivit(cid:19)e, par exemple les probl(cid:18)emes respiratoires, la gu(cid:19)erison de plaies ainsi que les processus isch(cid:19)emiques. Les bobines d’induction dans le dispositif de tomographie par induction magn(cid:19)etique pro- duisent un champ magn(cid:19)etique primaire dans la r(cid:19)egion d’int(cid:19)er^et (ROI pour Region of Interest) et ce champ alternatif induit des courants alternatifs (courants de Foucault) dans les r(cid:19)egions conductrices. Ces courants induits produisent (cid:18)a leur tour un champ magn(cid:19)etique secondaire dans la r(cid:19)egion d’int(cid:19)er^et. Ce champ magn(cid:19)etique secondaire produit un champ aux bobines de r(cid:19)eception qui est utilis(cid:19)e pour reconstruire la distribution de la conductivit(cid:19)e dans la r(cid:19)egion d’int(cid:19)er^et. Un d(cid:19)e(cid:12) important concernant l’(cid:19)etat de l’art de ces dispositifs est la d(cid:19)etection du signal se- condaire en pr(cid:19)esence du signal primaire, qui est plusieurs ordres de grandeur plus fort. Le dispositif pr(cid:19)esent(cid:19)e utilise une g(cid:19)eom(cid:19)etrie pour l’induction et la d(cid:19)etection sp(cid:19)ecialement con(cid:24)cue pour op(cid:19)erer dans les basses fr(cid:19)equences avec un ratio signal sur bruit acceptable de m^eme qu’une con(cid:12)guration du d(cid:19)etecteur qui est moins sensible au champ primaire. La m(cid:19)ethode directe du calcul de la matrice de sensibilit(cid:19)e introduite dans cette th(cid:18)ese nous fournit une m(cid:19)ethode num(cid:19)erique robuste pour la reconstruction d’images, ce qui r(cid:19)esulte en une qualit(cid:19)e d’image sup(cid:19)erieure par rapport aux autres m(cid:19)ethodes propos(cid:19)ees dans la litt(cid:19)erature. Les con(cid:12)gurations du syst(cid:18)eme que nous pr(cid:19)esentons impliquent une forme cylindrique avec vi 6 bobines d’excitation concentriques qui sont plac(cid:19)ees (cid:18)a diff(cid:19)erentes hauteurs sur la surface externe. Ces bobines produisent un signal primaire dont la majorit(cid:19)e des lignes de champ sont parall(cid:18)eles (cid:18)a l’axe principal du cylindre ou(cid:18) est plac(cid:19)e l’objet d’int(cid:19)er^et. Ce champ induit alors des courants de Foucault (courants de conduction et courants de d(cid:19)eplacement) dans les r(cid:19)e- gions conductrices, ce qui g(cid:19)en(cid:18)ere un champ magn(cid:19)etique secondaire. Ce champ secondaire est d(cid:19)etect(cid:19)e par 80 bobines de r(cid:19)eception (16 rang(cid:19)ees de d(cid:19)etecteurs plac(cid:19)es (cid:18)a 5 hauteurs diff(cid:19)erentes) qui sont plac(cid:19)ees sur les parois du cylindre de mani(cid:18)ere perpendiculaire au plan des bobines d’excitation. La simulation de ce dispositif n(cid:19)ecessite un mod(cid:18)ele math(cid:19)ematique comprenant le probl(cid:18)eme direct et le probl(cid:18)eme inverse. Le probl(cid:18)eme direct est la simulation math(cid:19)ematique du syst(cid:18)eme, qui r(cid:19)esout les (cid:19)equations de Maxwell dans la r(cid:19)egion d’int(cid:19)er^et avec les conditions fronti(cid:18)ere ad(cid:19)e- quates. Le probl(cid:18)eme impliquant les courants de Foucault est r(cid:19)esolu par l’utilisation d’un mod(cid:18)ele connu sous le nom d’(cid:19)equations de Maxwell compl(cid:18)etes (traduction libre de full Maxwell’s equations). Ce probl(cid:18)eme est r(cid:19)esolu par la m(cid:19)ethode des (cid:19)el(cid:19)ements (cid:12)nis ((cid:19)el(cid:19)ements en p(cid:19)eriph(cid:19)erie, traduction libre de edge elements) en consid(cid:19)erant les fonctions de Whitney de premier ordre. Le mod(cid:18)ele est compar(cid:19)e (cid:18)a des solutions analytiques connues pour des g(cid:19)eom(cid:19)etries simples. Le probl(cid:18)eme direct nous permet d’obtenir l’amplitude des champs magn(cid:19)etiques et (cid:19)electriques qui sont d(cid:19)etectables, la distribution et l’ordre de grandeur des courants de Foucault ainsi que le changement de phase du champ magn(cid:19)etique d(cid:19)etect(cid:19)e qui est caus(cid:19)e par la partie conductrice de la r(cid:19)egion d’int(cid:19)er^et. (cid:18) A partir du probl(cid:18)eme direct, la matrice de sensibilit(cid:19)e peut ^etre extraite, ce qui est n(cid:19)eces- saire a(cid:12)n de r(cid:19)esoudre le probl(cid:18)eme inverse. Les m(cid:19)ethodes disponibles pour extraire la matrice de sensibilit(cid:19)e, par exemple le th(cid:19)eor(cid:18)eme de r(cid:19)eciprocit(cid:19)e de Geselowitz, sont g(cid:19)en(cid:19)eralement des formulations ind(cid:19)ependantes du probl(cid:18)eme direct consid(cid:19)er(cid:19)e. Cependant, pour notre mod(cid:18)ele, la sensibilit(cid:19)e est calcul(cid:19)ee directement par la diff(cid:19)erentiation num(cid:19)erique de l’(cid:19)equation d’Helmholtz par rapport aux propri(cid:19)et(cid:19)es (cid:19)electriques. Cette m(cid:19)ethode pour calculer la matrice de sensibilit(cid:19)e (la m(cid:19)ethode directe) est introduite dans cette th(cid:18)ese et les r(cid:19)esultats de la reconstruction pour cette m(cid:19)ethode sont compar(cid:19)es avec ceux provenant du th(cid:19)eor(cid:18)eme de r(cid:19)eciprocit(cid:19)e appliqu(cid:19)e aux probl(cid:18)emes de tomographie par induction magn(cid:19)etique. La matrice de sensibilit(cid:19)e peut r(cid:19)ev(cid:19)eler (cid:18)a quel point la sensibilit(cid:19)e d’une partie de la r(cid:19)egion d’int(cid:19)er^et peut ^etre affect(cid:19)ee par un change- ment quelconque de courant ou de conductivit(cid:19)e (de m^eme que pour d’autres propri(cid:19)et(cid:19)es). De plus, extraire la matrice de sensibilit(cid:19)e permet de calculer la direction de sensibilit(cid:19)e maximale (cid:18)a un changement de conductivit(cid:19)e dans une certaine partie de la r(cid:19)egion d’int(cid:19)er^et, ce qui per- met de d(cid:19)eterminer le meilleur arrangement et emplacement des d(cid:19)etecteurs. L’extraction de la matrice de sensibilit(cid:19)e par la m(cid:19)ethode directe, (cid:18)a partir des outils permettant de r(cid:19)esoudre vii ce mod(cid:18)ele, n’est pas une t^ache simple. Ainsi, le mod(cid:18)ele direct doit ^etre d(cid:19)evelopp(cid:19)e dans son enti(cid:18)eret(cid:19)e dans MATLAB. Le probl(cid:18)eme inverse (la description de la structure interne du syst(cid:18)eme fournie par des don- n(cid:19)ees indirectes) est ensuite r(cid:19)esolu (cid:18)a l’aide de la matrice de sensibilit(cid:19)e et des m(cid:19)ethodes dispo- nibles dans la litt(cid:19)erature (ces m(cid:19)ethodes sont courantes pour la tomographie par imp(cid:19)edance (cid:19)electrique). Pour le probl(cid:18)eme inverse, r(cid:19)esolu par la m(cid:19)ethode du maximum a posteriori, les approches lin(cid:19)eaires et non lin(cid:19)eaires ont(cid:19)et(cid:19)e utilis(cid:19)ees a(cid:12)n de reconstruire la conductivit(cid:19)e. L’ap- proche lin(cid:19)eaire est prometteuse pour l’imagerie diff(cid:19)erentielle a(cid:12)n de con(cid:12)rmer l’apparition d’une l(cid:19)esion si nous disposons au pr(cid:19)ealable d’informations concernant la conductivit(cid:19)e nor- male de la t^ete et du cra^ne. Cependant, l’approche lin(cid:19)eaire am(cid:18)ene (cid:18)a des valeurs erron(cid:19)ees de la conductivit(cid:19)e. L’approche non lin(cid:19)eaire a (cid:19)egalement (cid:19)et(cid:19)e impl(cid:19)ement(cid:19)ee a(cid:12)n de permettre le calcul de la valeur absolue de la conductivit(cid:19)e sans aucune information pr(cid:19)ealable. Toutefois, l’approche non lin(cid:19)eaire est cou^teuse au point de vue des ressources computationnelles et donc faire la tomographie par induction magn(cid:19)etique en utilisant notre m(cid:19)ethode n’est pas envisa- geable pour calculer la valeur absolue de la conductivit(cid:19)e. La con(cid:12)guration des senseurs adopt(cid:19)ee pour notre dispositif, dont le positionnement des cap- teurs est perpendiculaire aux bobines d’excitation, am(cid:18)ene un avantage au point de vue de la soustraction du champ primaire, ce qui est un d(cid:19)e(cid:12) constant dans l’op(cid:19)eration des autres dispositifs de ce type. Les signaux ont (cid:19)et(cid:19)e simul(cid:19)es pour diff(cid:19)erentes hauteurs de disque dans le cylindre et diff(cid:19)erentes locations de la l(cid:19)esion (cid:18)a l’int(cid:19)erieur du disque. La distribution de la conductivit(cid:19)e est ensuite reconstruite (cid:18)a partir du changement de voltage induit dans les sen- seurs. Les r(cid:19)esultats du probl(cid:18)eme de reconstruction concernant le mod(cid:18)ele du disque dans le dispositif propos(cid:19)e, en consid(cid:19)erant que le niveau de bruit est constant et vaut 1% du signal maximal d(cid:19)etect(cid:19)e, montrent un ratio de signal sur bruit de 40 dB. Les capteurs magn(cid:19)etiques utilis(cid:19)es pour construire ce dispositif devraient ^etre capables de d(cid:19)etecter des champs de 10 pT avec ce m^eme ratio signal sur bruit. viii ABSTRACT This thesis develops a simulation package for the design of Magnetic Induction Tomogra- phy (MIT) instruments and exploit the simulator to investigate the possibility of using the magnetic induction tomography imaging method to design a non-contact device capable of detecting a blood hemorrhage inside the skull. The speci(cid:12)c numerical method (full Maxwell’s equation) used for simulation of the device, followed by calculation of the sensitivity with the direct method (introduced in this dissertation) and a regularized inverse solver which re- construct the conductivity map from the simulation outputs (magnetic (cid:12)eld), are optimized to simulate a device operating at 50 kHz. This device is capable of detecting the change in conductivity in ranges close to biological tissues. MIT operates based on the measurement of passive electromagnetic properties such as con- ductivity. The rationale behind using this method for detecting cerebral stroke is based on the fact that the conductivity of the blood is larger than that of the other tissues in the head. Other potential medical applications for this device are real-time, non-invasive monitoring of tissue alterations which are re(cid:13)ected in the change of the conductivity, e.g. ventilation disorders, wound healing and ischemic processes. The inductive coils in the MIT device produce a primary magnetic (cid:12)eld in the region of interest (ROI) and this alternating magnetic (cid:12)eld induces alternating (eddy) currents in the conductive regions. These eddy currents, in turn, generate a secondary magnetic (cid:12)eld in the ROI. This secondary magnetic (cid:12)eld generates a (cid:12)eld at the receivers, which is used to reconstruct the conductivity distribution of the ROI. An important challenge in the state of art devices is the detection of the secondary signal in the presence of the primary signal, which is orders of magnitude stronger. The device uses a geometry in induction and detection designed to operate in low frequencies with ac- ceptable SNR and detector con(cid:12)guration that is least sensitive to the primary (cid:12)eld. The direct method of calculating the sensitivity matrix introduced in this thesis provides us with a robust numerical method for image reconstruction which results in superior image quality compared to other proposed methods in the state of the art. The proposed con(cid:12)gurations involves a cylindrical shape device with 6 concentric excitation coils which are located at different heights on the outer surface of the cylinder. These coils produce a primary magnetic (cid:12)eld with majority of the (cid:12)eld lines parallel to the main axis of the cylinder, where we position the object of interest. This primary (cid:12)eld induces eddy currents (conduction and displacement current) in the conductive regions of the ROI, which generate a secondary magnetic (cid:12)eld. The secondary magnetic (cid:12)eld is detected in 80 sensors ix (16 arrays of sensors at 5 different heights), which are located on the cylinder walls perpen- dicular to the plane of excitation coils. The simulation of this device requires a mathematical model consisting of a forward and an Inverse problem. The forward problem is the mathematical simulation of the device, which solves Maxwell’s equations in the ROI with the appropriate boundary conditions. The problem is known as the eddy current problem is solved using a model known as full Maxwell’s equations. The eddy current problem is solved using a (cid:12)nite element method (edge basis) considering the (cid:12)rst order Whitney functions. The model is bench-marked with avail- able analytic solutions for known geometries. The forward problem provides us with the magnitude and direction of the detectable magnetic and electric (cid:12)eld, eddy current size and distribution and also the phase change in the detected magnetic (cid:12)eld, which is caused by the conductive area in the ROI. From the forward problem, the sensitivity matrix is extracted that is used in solving the reconstruction problem. Methods available for extracting the sensitivity matrix, like using the Geselowitz reciprocity theorem, are mainly independent formulations regardless of the method used in solving forward problem. However, in our model, the sensitivity is directly calculated by the rigorous numerical differentiation of the governing equation (Helmholtz) in the forward problem, with respect to the electrical properties. This method of calculating the sensitivity matrix (the direct method) is introduced in this thesis and the results of re- construction for this method are compared with results from the reciprocity theorem. The sensitivity matrix shows how sensitive is the secondary (cid:12)eld of an area in ROI to any change in current or conductivity of other areas (as well as other properties such as permittivity). Furthermore, extracting the sensitivity matrix could allows one to systematically investigate the most sensitive direction to the conductivity change in certain area, at any position in the region of interest to (cid:12)nd the the best arrangement and alignment for the sensors. Extracting the sensitivity matrix with the direct method, from available packages which are capable of solving the forward model, is not a straightforward task. Therefore, the forward model is fully developed in MATLAB. The inverse problem (the description of the internal structure of a system given by indirect data) is then solved using the sensitivity matrix and available methods known for solving inverse problems in the literature (These matrix inversion methods are common with Elec- trical Impedance Tomography). In the inverse problem, which is solved using the maximum a posterior method, both the linear and non-linear approach were taken to reconstruct con- ductivity. The linear approach is promising for differential imaging performed on a head modelled as a disk with different conductivity layers to con(cid:12)rm the appearance of the lesion, if we have prior information about the conductivity of the background (the head and the