1w LI 1 M 'I :II. ll‘l‘HbHH 2m --1 w 1%: l mco—n I V LIBRARY Michigan State University This is to certify that the thesis entitled A TWO-STAGE METHOD FOR ESTIMATING RELATIVE PRESSURE FIELDS FROM NOISY VELOCITY DATA presented by Christopher David Bolin has been accepted towards fulfillment of the requirements for the MS. degree in Mechanical Engineering 4an - Major Professor’s Signature cST/i%IO-‘I Date MSU is an Affirmative Action/Equal Opportunity Employer PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5/08 K:IProj/Acc&Pres/ClRC/DaIeDue.indd A TWO-STAGE METHOD FOR ESTIMATING RELATIVE PRESSURE FIELDS FROM NOISY VELOCITY DATA By Christopher David Bolin A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IVIechanical Engineering 2009 ABSTRACT A TWO-STAGE METHOD FOR ESTIMATING RELATIVE PRESSURE FIELDS FROM NOISY VELOCITY DATA By Christopher David Bolin The determination of intravascular pressure fields is important to the character- ization of cardiovascular pathology. A two-stage method that solves the problem of estimating the relative pressure field from noisy velocity fields on an irregular do- main with limited spatial resolution and includes a filter for the experimental noise is presented here. For the pressure calculation, the pressure Poisson equation is solved by embedding the irregular flow domain into a regular domain. To lessen the prop- agation of the noise inherent to the velocity measurements, three filters — a median filter and two physics-based filters — are evaluated and compared to three filters from the literature using a pair of two—dimensional mathematical phantoms. The most accurate pressure field results from a filter that applies in a least-squares sense three constraints simultaneously: consistency between measured and filtered velocity fields, divergence-free and removes noise in the Laplacian of the velocity field. This filter leads to a 5-fold gain in accuracy for the estimated relative pressure field compared to without noise filtering, using spatial resolutions that are consistent with phase— contrast magnetic resonance imaging (MRI) of the carotid artery on a clinical MRI scanner in the more complex of the two phantoms. To my friends and family. iii ACKNOWLEDGMENTS I wish to acknowledge those with out whom this work would not have been possi- ble. First, I would like to thank my committee members Professor Seungik Baek and Professor David Zhu for aiding me in the preparation of this thesis. I also wish to thank all of the course instructors I have had through my undergraduate and grad- uate studies for providing me with necessary the background knowledge to complete this work. Last, but certainly not least, I would especially like to thank my research advisor, Professor L. Guy Raguin for introducing me to this topic area and guiding me throughout this process. I would also like to express my appreciation to those who provided non-academic support. I am grateful to my fellow graduate students who proved to be excellent distractions whenever they were necessary, and sometimes when they were wholly unnecessary. Thank you to my friends from my undergraduate years for entertaining me with your constant stream of emails about life in the ‘real world’ and the occasional visit. Finally, I would like to thank my parents and family for giving me ample opportunities to perform manual labor, feeding me good food whenever I came home and the constant encouragement. iv “If we knew what it was we were doing, it would not be called research, would it?" - Albert Einstein TABLE OF CONTENTS LIST OF TABLES ............................... LIST OF FIGURES .............................. LIST OF SYMBOLS .............................. LIST OF ABBREVIATIONS ........................ 1 Introduction ................................. 1.1 Motivation ................................. 1.2 Methods of Pressure Estimation ..................... 1.3 Minimally Invasive Blood Flow Assessment ............... 1.4 Thesis Layout ............................... 2 Methodology ................................. 2.1 Pressure Estimation ................... ' ........ 2.1.1 Theory ............................... 2.1.2 Numerical Implementation .................... 2.2 Proposed Noise Filters .......................... 2.2.1 Filter Description ......................... 2.2.2 Filter Implementation ...................... 2.3 Mathematical Phantoms & Validation .................. 2.3.1 Mathematical Phantoms ..................... 2.3.2 Resolution and Embedding Experiments ............ 2.3.3 Filter Parameter Estimation and Noise Experiments ...... 3 Results & Discussion ............................ 3.1 Influence of Spatial Resolution and Embedding Strategy ....... 3.2 Calibration and Performance of Proposed Noise Filters ........ 3.3 Discussion ................................. 4 Conclusions .................................. 4.1 Insight ................................... 4.2 Future Work ................................ APPENDICES ................................. A Finite Difference Schemes ......................... vi 16 16 17 19 22 23 28 29 3O 34 39 57 64 66 B Fourier Poisson Equation Solver ..................... 74 B.1 2D Neumann Problem .......................... 74 8.2 2D Dirichlet Problem ........................... 80 C Mathematical Phantom .......................... 84 Cl Couette Flow Between Rotating Cylinders ............... 84 C2 Poiseuille Channel Flow ......................... 86 D Experimental Parameters ......................... 88 E Tabulated Experimental Results ..................... 91 BIBLIOGRAPHY ............................... 95 vii B.1 D1 D2 D3 El E2 E3 E4 E.5 LIST OF TABLES The weights for the perturbation required to solve the 2D Neumann Poisson equation. ............................. Parameters used in the Couette flow mathematical phantom in the resolution, embedding and initial noise experiments ........... Parameters used in the Couette flow mathematical phantom noise ex- periments with more realistic viscous terms. .............. Parameters used in the Poiseuille flow mathematical phantom noise experiments ................................. The normalized relative RMS error, 5, results of the isotropic and anisotropic resolution tests (see Figures 3.1 and 3.2). ......... The normalized relative RMS error, E (‘76), results of the embedding tests (see Figure 3.3) ............................ The filter parameters that yield optimal performance at SNR = 5. These are the parameters used in all future experiments with noisy velocity fields (see Figures 3.4 and 3.5) .................. The average normalized relative RMS error and standard deviation, E i a (‘70), in the estimated relative pressure field, for the initial Cou- ette flow parameters, a resolution of n = 0.05 on a 20 x 20 flow do- main .0p and an embedding level P = 25%, after applying the three filters proposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering and three filters from the lit- erature (PSDF [1], (d/d:1.‘)2 [2], Visc. Diss. [3]) over the range of SNR tested (see Figure 3.6). .......................... The average performance factor F and standard deviation, F :l: a, for all filters tested using the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 x 20 flow domain .01: and an embedding level P = 25% across the range of SNR tested (see Figure 3.8) ..... viii 79 89 89 90 91 92 92 93 93 E6 E7 The average normalized relative RMS error and standard deviation, 5 :l: a (%), in the estimated relative pressure field, for the second Cou- ette flow parameters, a resolution of r} = 0.05 on a 20 x 20 flow do— main .QF and an embedding level P = 25%, after applying the three filters proposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering (see Figure 3.10). ...... The average normalized relative RMS error and standard deviation, E :t o (%), in the estimated relative pressure field, for the Poiseuille flow parameters, a resolution of 7) = 0.05 on a 20 x 20 flow domain .01: and an embedding level P = 25%, after applying the three filters pro- posed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering and two filters from the literature ((d/dx)2 [2], Visc. Diss. [3]) over the range of SNR tested (see Fig- ure 3.11). ................................. ix 94 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 3.1 LIST OF FIGURES The algorithm for solving the PPE on an irregular domain as described in [4]. The Poisson equation is solved in step 4 using the direct method outlined in Figure 2.2 and detailed in Appendix B.1. ......... The algorithm for the direct solver of Poisson equations based on Fourier transforms as described in [5] that completes step 4 in Fig- ure 2.1. Greater detail can be found in Appendix B.1 .......... The algorithm for projecting the noisy velocity field into the space of divergence-free fields as proposed in [1] .................. A schematic of the Couette flow system used in the 2D validation experiments with the important geometric features identified ...... Example sections of the analytical solution for the Couette flow math- ematical phantom: (a) the velocity field, V A; (b) the magnitude of the pressure field relative to its spatial average, PA, and normalized by its RMS, FARMS, on 0F. .......................... A schematic of the Poiseuille flow system used in the 2D validation experiments with the important geometric features identified ...... Example sections of the analytical solution for the Poiseuille flow math- ematical phantom: (a) the velocity field at a arbitrary axial location, V A; (b) the magnitude of the pressure field relative to its spatial av- erage, PA, and normalized by its RMS, FARMS, on .012 ......... Illustrations useful in visualizing the embedding tests, (a) the defini- tion of the level of embedding F and an schematic of the fluid do- main .QF embedded in the center of the extended computational do- main .Q, (b) 0p embedded centered in one direction and off-center in the other and (c) .01: embedded in a corner of (2. ........... Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of the the relative pixel size, 77 = Ax/ng. ................... 20 22 29 31 33 34 35 38 44 3.2 3.3 3.4 3.5 3.6 3.7 Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of Xa : a/nconv ................................ Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of the embedding level F with a 20 x 20, which corresponds to r) = 0.05, How domain 0p. ................................ Plot of the average normalized relative RMS error, E, of the estimated relative pressure field for the initial Couette flow experimental parame- ters, a resolution of r) = 0.05 on a 20 x 20 flow domain 0p, an embedding level P = 25% and SNR = 5, as a function of: (a) a’l for Filter 1; (b) 0’2 and 7’2 for Filter 2, the minimum is circled. ............ Contour plots of the average normalized relative RMS error, E, of the estimated relative pressure field for the initial Couette flow experimen— tal parameters, a resolution of 17 = 0.05 on a 20 x 20 flow domain .01: and an embedding level P = 25% and SNR = 5, as a function of: (a) as and 73 for the filter with filtering term proposed in [2]; (b) 0:1 and 7:1 for filter with the filtering term proposed in [3]. The minimum is circled in both figures. .......................... Plot comparing the average normalized relative RMS error, E, in the es— timated relative pressure field, for the initial Couette flow parameters, a resolution of r) = 0.05 on a 20 x 20 flow domain {21: and an embedding level P = 25%, after applying the three filters proposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimation with the field estimated with no filtering and three filters from the literature (PSDF [1], (d/dar)2 [2], Vise. Diss. [3]) over the range of SNR tested. Sample images of the magnitude of the normalized relative error dis— tribution, |e(:r, y)|, for the estimated relative pressure fields obtained using the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 x 20 flow domain 0p, an embedding level P = 25% and SNR = 5 (a) no filter; (b) a median filter with radius 1 pixel; (c) Filter 1; (d) Filter 2; (e) a filter using the term of [2]; and (f) a filter using the term of [3] on the noisy velocity field. The normalized RMS error, E, is specified for each case. ................ xi 45 46 47 48 49 3.8 3.9 3.10 3.11 B.1 C.1 Plot comparing the average performance factor F to SNR for all fil- ters tested using the initial Couette flow experimental parameters, a resolution of n = 0.05 on a 20 X 20 flow domain {21: and an embedding level P = 25%. Note that the median filter, the PSDF method and the filter using the viscous dissipation term of [3] all become ineffective by SN R = 15 .................................. Sample velocity vector fields of the initial Couette flow experimental parameters, a resolution of n = 0.05 on a 20 x 20 flow domain 01:, an embedding level P = 25% and SNR = 5 illustrating the effect of (a) no filter (Vexp); (b) a median filter with radius 1 pixel; (c) Filter 1; (d) Filter 2; (e) the filter using the term of [2], see (2.11); and (f) the filter using the term of [3], see (2.12), on the noisy velocity field. The subscript f indicates the filtered field. The performance factor, F, is specified for each case. .......................... Plot of the average normalized relative RMS error, E, in the estimated relative pressure field versus SNR for the more realistic Couette flow parameters, resolution 77 = 0.05 on a 20 X 20 flow domain .01: and an embedding level of F = 25%, comparing the results of the three filters 52 53 proposed here to the result with no filter over the range of SNR tested. 54 Plot of the average normalized relative RMS error, E, in the estimated relative pressure field versus SNR for the Poiseuille flow parameters, resolution 77 = 0.05 on a 20 x 20 flow domain 9F and an embedding level of F = 25%, comparing the results of the three filters proposed here to the result with no filter and two filters from the literature (d/ d:1:)2 [2], Visc. Diss. [3]) over the range of SNR = 5 — 20. ............ A schematic illustrating the computational domain .0 and defining the names of the boundaries and indices used in the direct solver of the Poisson equation in two—dimensions. .................. A schematic illustrating the components of v9 for the case where or = in the Cartesian coordinate system ................... xii 55 85 LIST OF SYMBOLS Roman Letters TI 7'2 Ra External magnetic field strength (T) Arterial diameter (n1) Performance factor of velocity filters (-) Height of the channel of the Poiseuille flow mathematical phantom (m) Index in the x—direction in real space (-) Index in the y-direction in real space (-) Step-size change for line-search algorithm (-) Index in Fourier space corresponding to z' (-) Maximum index in the x-direction (-) Maximum index in the y—direction (—) Number of points across the computational domain (-) Number of points across fluid domain (-) Pressure field (N m—2) Estimated relative pressure field (N m‘2) True relative pressure field (N m‘2) Flow rate per unit depth (m2 s—l) Radial coordinate (m) Radius of inner cylinder of Couette flow mathematical phantom (n1) Radius of outer cylinder of Couette flow mathematical phantom (In) Radius of a blood vessel (m) Step index for algorithms (-) xiii Ur ’Ug Time (8) Average arterial velocity (m s—l) Velocity component in the x-direction (111 8—1) Velocity component in the y-direction (m s-l) Radial velocity component (m s‘l) Azimuthal velocity component (m S—l) Spatial coordinate (In) Spatial coordinate (m) Reynolds number (-) Womersley number (-) Average Reynolds number of Couette mathematical phantom (-) Average Reynolds number of arterial flow (-) Bold symbols and vectors V b , 0" ~11 Will i} The filter of [1] was proposed as a method of filtering PC-MRI meeasurements n == 7’4 {HV . Vexpll}[ respectively. based on the notion that the measured velocity field must be divergence free. The velocity field is projected onto the space of divergence—free fields using a discrete projector P that has the form PV = V —— GL‘IDV (2.17) 27 where G is the discrete gradient operator, L is the discrete Laplacian operator and D is the discrete divergence operator. L—IDV is found by solving the Poisson equation V2L_1DV = V-V on a 2.1 L’lDV = 0 on an ( 8) The projection generates a divergence free velocity field, but it does not eliminate all of the noise. It also does not. guarantee that the filtered velocity field will still be physically relevant. 2.2.2 Filter Implementation The filters tested in this work are implemented in manners of varying complexity. For example, the medf 111:2 function in MatLab is used to perform the median filtering whereas each cost function and the PSDF filter require their own special functions. The complexity of the special functions affects the time required to filter the ve- locity field. Whenever possible MatLab’s built-in functions are used to reduce the complexity of the algorithms. - Filter 1, Filter 2, and the filters based on [2] and [3] considerably more compu- tational resources than the median filter. The divergence operator in f1 and f2 and the Laplacian operator in f2 are applied using matrix algebra in custom functions written for the 2D validation. This is done to make use of the speed of MatLab at performing these operations. Derivatives are numerically approximated in the same method as described in Section 2.1.2 for the approximation of bp and are detailed in Appendix A. Care is taken on 80p, as identified by a mask, not to include points from outside 0F. The matrices containing the operators are generated once and are inputs to the functions that evaluate f1 and f2. These matrices are very large and sparse. In 2D, if .0 of size M x N points, then the coefficient matrices are of size 2(M N ) x 2(M N ) This can lead to storage issues for large computational domains. The evaluation of f3 and f4 is done element by element due to the complexity of their 28 additional filtering terms. All four cost functions f1, f2, f3 and f4 are minimized us- ing the unconstrained nonlinear minimization algorithm, fminunc, in MatLab. This algorithm uses a simple cubic approximation method to minimize the two functions. The PSDF filter is conviently implemented using the the same Poisson solver used in the relative pressure field estimation. The algorithm for filtering the velocity field is in Figure 2.3. This system has Dirichlet boundary conditions which require sine transforms and inverse transforms and do not require the extra steps necessary for the solution of a Neumann problem. This version of the Poisson solver is described in detail in Appendix B.2. step 1 Compute the divergence of the velocity field DV step 2 Solve the Poisson equation using the algorithm described in 2.1.2 yielding L“1DV step 3 Take the gradient of the solution to the Poisson equation re— sulting in GL—IDV step 4 Subtract the result of step 3 from the original velocity field V yielding the filtered velocity field V — GL‘IDV Figure 2.3. The algorithm for projecting the noisy velocity field into the space of divergence—free fields as proposed in [1]. 2.3 Mathematical Phantoms & Validation Several components of the the proposed method require validation and optimization as a sort of ‘proof of concept’ to assess which approach is best. This section is divided into three parts. First, a mathematical phantom of suitable complexity to test the algorithms is described. The second part of the section outlines a testing procedure for determining the effect of resolution and embedding on the estimation of the pressure 29 field. In [4] the proposed algorithm for estimating the relative pressure field was tested for robustness to noise at a higher resolution than that expected for PC-MRI data of the carotid artery. No results were reported on the effect of resolution or the size of .0 relative to 0p. In the final part the procedures for optimizing the filter parameters described in Section 2.2.1, the tests used to determine the method’s robustness to noise at a realistic resolution for PC-MRI measurements of the carotid artery, determining the effect of the proposed filters and testing the scaling of the filter parameters proposed in Section 2.2.2. 2.3.1 Mathematical Phantoms Examination of the usefulness of the relative pressure estimation method proposed in this work, like all numerical methods, requires comparing results to an example with known properties. This can be done by an experiment on a physical system performed with additional measurements provided by a trusted second method as in [13, 14, 17, 20], experiments on a physical phantom with well understood behavior like in [12], or numerical experiments using either commercial CFD software [11], or a mathematical phantom of sufficient complexity [4]. Experiments with physical phantoms and PC- MRI are expensive and noise free tests cannot be guaranteed. Comparing a numerical solution to CFD in a complex geometry allows for noise-free tests, but questions about the accuracy of the CF D model can still be raised. Instead, a mathematical phantom of a suitably complex system with a known analytical solution is chosen for use here because of the flexibility it allows in the testing and the simplicity of its implementation. For the validation tests, a mathematical phantom of Couette flow between rotating cylinders is used. Certain comprimises must be made when choosing a mathematical phantom. Couette flow between rotating cylinders provides a nontrivial system with an analytical solution for both the pressure and velocity fields. The cylindrical shape 30 Figure 2.4. A schematic of the Couette flow system used in the 2D validation exper- iments with the important geometric features identified. of this system also allows for identification of potential issues with the calculation of bp near boundaries that are similar to those anticipated in vessel geometry. The two boundaries in the system, on the inner and outer cylinder walls, also test the effect of multiple, moving embedded boundaries on the PPE solver. This is expected in viva for data collected above the internal/ external bifurcation in the carotid artery. However, one of the characteristics of this flow field is that the viscous terms of (2.3) are identically zero. This is a comprimise that must be made to ensure that the majority of the components of the algorithms are tested sufficiently. The derivation of the flow field from (2.3) is detailed in Appendix C1 The final results are presented in this section. The Cartesian velocity field (a, v) in polar coordinates (r, 6) is given by: u(r,0) = — (Ar + E) sin0 B 7‘ (2.19) v(r, 0) 2 (Ar + 7) c050 31 where r := V452 + 312, 6 := arctan (y/z‘) and wgrg — 0217‘? 2_ 2 7'2 ’1 A: B = rfrg (tal — (.02) 2 2 7’2‘7‘1 r1 and r2 are the radii of the inner and outer cylinders, respectively, and <01 and 022 are their respective angular velocities, see Figure 2.4 for a schematic of the system. Though the flow is 1D in polar coordinates, the transformation to Cartesian coor- dinates generates a flow with non-zero gradients in both directions for at least one component of velocity at every point. Figure 2.5(a) is a vector plot of half of the analytical velocity field, V A, which illustrates its complexity. An analytical solution for the relative pressure field is similarly obtained: 2 2 2 P(r)=p A r +2ABlnr—B—+C (2.20) 2 2r2 where C is an arbitrary integration constant, such that P is the relative pressure field as defined earlier. An example of half of the analytical pressure field relative to its spatial average, PA, and normalized by its root-mean-square, PARA-is, on .01: (RMS) is found in Figure 2.5(b). A second simple phantom is necessary to determine if any bias has been introduced by the combination of filter and mathematical phantom. Poiseuille channel flow is a simple system with a 1D analytical velocity field and linear pressure gradient, see Figure 2.6 for a schematic. The Cartesian velocity field is ad) = 13% (y2 -— by) (2.21) where u is the velocity component along the :r-direction, QC is the volume flow rate per unit depth and h is the height of the channel. Figure 2.7(a) is a vector plot of the analytical velocity field at an arbitrary :r-location, V A- The analytical pressure field is 12QCH 13(1)) = h3 (:1: — 1‘0) (2.22) 32 (a) Analytical velocity field [PA/FARMS] (b) Normalized analytical pressure field Figure 2.5. Example sections of the analytical solution for the Couette flow mathe— matical phantom: (a) the velocity field, VA; (b) the magnitude of the pressure field relative to its spatial average, PA, and normalized by its RMS, FARMS, on .017. where 2:0 is an arbitrary datum. An example of the analytical pressure field relative to its spatial average, PA, and normalized by its root-mean-square, PA,RMSi on 01; (RMS) is found in Figure 2.7(b). In this mathematical phantom the viscous terms of (2.3) are not identically zero, see Appendix C.2 for a detailed derivation. The mathematical phantoms described here are implemented using custom func- tions in MatLab. The physical domain in which they are contained is fixed, but the resolution in all directions is variable to facilitate the resolution experiments described in Section 2.3.2. These functions also generate the mask used in the algorithms of Section 2.1.2 and Section 2.2.2 for identifying points on 01:. 33 “y//////[/// ' h - --][—-F--P--v-_- -" 0 _ //////// 7x7 Figure 2.6. A schematic of the Poiseuille flow system used in the 2D validation experiments with the important geometric features identified. 2.3.2 Resolution and Embedding Experiments As stated above, no testing of the algorithm in Section 2.1.2 for robustness to changes in resolution or the level of embedding have been, to the best of the author’s knowl- edge, published. In [4] the algorithm was tested on a phantom similar to the 2D Couette flow between rotating cylinders used here with a grid size of 256 x 256. The carotid artery has a diameter of ~ 10 mm and PC-MRI velocity data is expected to have an in—plane voxel size of ~ 0.5 mm x 0.5 mm. Thus, approximately twenty voxels are expected across the vessel. This low resolution could be enhanced using interpolation to increase the number of data points on 01:, but upsampling like this with noisy data is problematic. The effect resolution has on the algorithm is investi- gated here to determine if upsampling can be avoided. The resolution of the velocity data cannot be guaranteed to be isotropic, so tests on the effect of anisotropic resolu- tion in both directions are undertaken here to explore its influence on the estimated pressure field. The size of 0 can also be an issue from a computational point of view. 34 [[11111]] (a) Analytical velocity field I PA/PA,RMS l (b) Normalized analytical pressure field Figure 2.7. Example sections of the analytical solution for the Poiseuille flow math- ematical phantom: (a) the velocity field at a arbitrary axial location, VA; (b) the magnitude of the pressure field relative to its spatial average, PA, and normalized by its RIVIS, PA,RMSv on 0F- For maximum efficiency .0 should be as small as possible. An investigation into the relative size of 0 compared to 0p is performed to determine a reasonable relative size for .0 or if an optimum size exists. In 3D, a single contiguous domain is desired to eliminate the difficulties of point registration. It is expected that the carotid artery will not remain centered in each plane of this domain. Therefore. the effect of the placement of .017 Within .0 is also investigated. The behavior of the pressure estimating algorithm with respect to resolution is required before any of the other experiments can be performed. If the estimation can 35 be performed accurately in noise-free conditions, the rest of the tests should be con— ducted at realistic PC-MRI resolutions. However, if the algorithm does not output results of sufficient accuracy at realistic resolutions without noise, then additional steps, like the upsampling mentioned above, will be required for the following ex- periments. Only the image matrix size and not the size of the voxels was reported in [4], so determining the resolution in that study is impossible. For these resolution experiments, the 2D mathematical phantom described in Section 2.3.1 is used. The parameters used to generate the mathematical velocity and pressure fields for these experiments are in Appendix D. The size of .01: and the size of .0 relative to .01: are both fixed. To vary the resolution, the number of points in the computational mesh is increased or decreased. In these experiments the resolution is isotropic. The size of the mesh is varied from 2048 X 2048 to 16 X 16 which corresponds to relative pixel size, 77 = Amt/2m, of 7 X 10—4 and 0.09, respectively. PC-MRI data of the carotid artery, with 7) = 0.05, falls within this range. The metric (E) used to determine the accuracy of the algorithms is RMS (PA) a = (2.23) where Pnum is defined as before and PA is the analytical pressure field minus its spatial average (i.e., PA is an absolute pressure field). At 77 = 0.05, E is less than 1% (see Section 3.1). This is deemed sufficient for the rest. of the experiments to be conducted at this resolution without upsampling. It is possible that the measured velocity field will have an anisotropic resolution and that this will affect the estimated pressure field. The Poisson equation solver used here relies on a FF T to decouple the two spatial directions and was derived to be robust to anisotropic resolution over a regular domain. However, the additional steps required for solutions on an irregular domain have the potential to favor isotropic resolution for more accurate solutions. Because the FF T is applied in one direction, the effect anisotrOpic resolution on the irregular solution can be dependent on which 36 direction has lower resolution. Two tests are carried out to determine these effects. In the first, the relative pixel size in the direction of the FFT, "FFT = A3'313‘FT/(2T2li is held constant at 77mm, = 0.0056 while 7750! = A1350] / 2r2, the relative pixel size in the solution direction, is allowed to vary from 0.0004 to 0.09. The second test holds 77501 = 77mm, and allows TIFFT to vary from 0.0004 to 0.09. The fixed relative pixel size was chosen such that it fell in the converged region as determined by the initial isotropic resolution tests. The size of .0 relative to .017 in both directions is fixed. Determining the effect of embedding on the algorithm is important not only from the standpoint of the accuracy of the solution, but also from a computational stand- point. The size of .0 relative to 01: is referred to as level of embedding (F) and is defined as NQ F = —— — 2.24 Na ( > where N g is the number of voxels across .0 and N 9F is the number of voxels across the diameter of the outer cylinder in the 2D mathematical phantom, see Figure 2.8(a) for an illustration. This experiment is performed with N 9F = 20 and 77 = 0.05. Ex- periments are conducted with .0]: located in three different regions of .0 to simulate the expected meandering of 0p in an in viva data set. In the first experiment 0F is embedded centered in .0 as shown in Figure 2.8(a) while in the other two experiments 0F is embedded centered in one direction but off center in the other, see Figure 2.8(b), and in the corner of 0, see Figure 2.8(c). For all three experiments the level of em- bedding is allowed to vary from 0% to 250%. Based on the results of this experiment (see Section 3.1) a level of embedding of 25% is deemed suitable for the rest of the 2D experiments. 37 (C) (b) Figure 2.8. Illustrations useful in visual- izing the embedding tests, (a) the defini- tion of the level of embedding F and an schematic of the fluid domain 0p em- bedded in the center of the extended computational domain .0, (b) .01: em- bedded centered in one direction and off- center in the other and (c) 01: embedded in a comer of .0. 38 2.3.3 Filter Parameter Estimation and Noise Experiments Once the determination of how the method is affected by resolution has been made, the level of embedding and the location of .01; within .0, tests to determine the usefulness of filtering the velocity field before estimating the relative pressure field can begin. This is done in two steps. In the first step the same set of experimental parameters as were used in the tests of Section 2.3.2 to optimize scaled versions filter parameters of Section 2.2.2. Once the scaled parameters are found, the filters are tested using the same experimental parameters over a range of SNR to determine their effectiveness. Based on this testing, three filters are chosen to be representatve of the filtering options in terms of speed of implementation and accuracy. These representatives are put through another set of noise tests using different experimental parameters to test the scaling of the filter parameters. The SNR for experiments with a noisy velocity field is defined as RMS (VA) SW = ‘ RMS (Vexp — V A) (2.25) where RMS is defined as before, Vexp is the noisy experimental velocity and V A is the analytical velocity field. Noise is added to the V A with the randn function in MatLab which generates random numbers with a Gaussian distribution and the level is adjusted based on the desired SNR. The two proposed physics-based filters and two of those described in the literature require the tuning of the parameters (11’, 012’ and 72’, 03’ and 73’, and 04’ and 74’. This is done using the worst conditions considered here (SNR=5) and optimizing the parameters such that the pressure solution resulting from the filtered velocity fields has the smallest E. These parameters are optimized empirically to determine the order of magnitude for the penalty terms in the cost functions. It is expected that these terms should be small compared to the term representing (2.2) so this level of optimization is deemed sufficient because it is unlikely that a set of parameters can 39 be found that will be optimal for all noisy velocity fields. The optimized values of the parameters, 01 ’ = 0.1, 02' = 0.01 and 7’2 = 0.53, are used in the application of the two proposed physics-based filters Filter 1 and Filter 2, see Section 2.2.2. The optimized parameters for the filters from the literature, 03 = 0.01 and '73 = 0.1 and of, = 0.0001 and ya = 0.1, are used in the application of those filters. The first set of tests with the noisy velocity fields are used to determine the effect of decreasing SNR on the algorithm and the effectiveness of the filters using the initial filter parameters. It was shown in [4] that the algorithm was robust to noise in 2D experiments using a mathematical phantom at a higher resolution than that of interest here. White Gaussian noise is added to the analytical velocity field at SNRs of 5, 10, 15, 20, 25 and 30 and the resulting noisy velocity field is input either directly into the algorithm of Section 2.1.2 or filtered before being input into the algorithm. A PC-MRI data set is expected to have a SNR between 5 and 10. However, with advances in sequence design or technology it is possible for this to increase, so the higher SNRs are tested to verify the filters do not enhance the proportion of noise in these conditions. Several noisy velocity fields are used to test each of the filters and the relative pressure estimation without a filter and the results are averaged. A method of determining the quality of the filtered velocity field is of use as a way to directly compare the effect of the filters on the velocity field rather than on their effect on the resulting estimated relative pressure field. Such a method was described in [1] to evaluate the velocity filter pr0posed in that work. The performance factor F of a filter is defined as F ___ “VA " Vfll “VA — VeXPll (2.26) where Vf is the filtered velocity field. F is the ratio of the error in the filtered velocity field to the error in the noisy velocity field. If F < 1, then the filter is effective in moving the the velocity field towards the true field, whereas F > 1 indicates that the filter is no longer effective and is moving the velocity field away from the true 40 field. This metric is used in the initial noise experiments as an additional method of differentiating between the proposed filters and evaluating their performance as compared to those filters from the literature. A second set of experiments with noisy velocity data are performed using only the median filter and the two proposed physics-based filters proposed here. Based on the initial set of experiments, the three proposed filters illustrate the range of filtering options adequately. This second set of experiments are carried out using a different set of flow parameters to determine the appropriatness of the proposed scaling of the filter parameters. The second set of experimental parameters attempt to match the scale of the velocities that occur within the carotid artery. Setting experimental parameters that attempt to more closely match those in arterial flow is somewhat difficult with the Couette flow mathematical phantom used here. There are no direct connections in terms of length scales between this phantom and Carotid flow which makes matching the Reynold’s number (Re) complex. Additionally, arterial flow is normally described, as stated in 2.1.1, with W0 not Re. As time dependence is not being considered at this stage, an attempt is made to mathematically approximate a flow with a more accurately scaled velocity field than that in the initial set of tests by matching an average Re rather than Womersley numbers for the two flows. To more closely match the behavior of blood, its average density and viscosity are used. The average Re in the Couette flow phantom defined as 0217‘1 + (.0’27’2 2 (2.27) where the subscript C indicates the Couette flow phantom, is approximately matched to an average reynolds number in carotid artery flow defined as D _ Rea = BLUE (2.28) p where the subscript 0 indicates the approximate conditions in the carotid artery. This is not a direct scaling between the two flows. The mathematical phantom parameters 41 used to match these numbers are in Appendix D. The matching of these two numbers allows for only a very rough approximation of the velocity field which is adequate for testing the scaling of the filter parameters. These two sets of experiments represent only one possible flow geometry. As dis- cussed in Section 2.3.1, certain compromises are made in the selection of the Coutette flow mathematical phantom. In particular the viscous terms are identically zero in the analytical description of the phantom. The constraint in Filter 2 that minimizes the Laplacian of the noisy velocity field is attempting accomplish this same requirement. This indicates a potential bias in the results of this work. To determine the impor- tance of this bias, a set of simple experiments were undertaken using the Poiseuille flow phantom. The flow domain is embedded in the center of an extended domain as in the Couette noise experiments, noise is added to the analytical velocity field at SNR = 5 — 20 and the filter parameters optimized for the Couette flow mathemat- ical phantom are used. The median filter, Filter 1, Filter 2, (d / dc)2 Filter and the Visc. Diss. Filter are tested and compared to the relative pressure estimation with no filtering. The mathematical phantom parameters used are in Appendix D. 42 CHAPTER 3 Results & Discussion The results of the 2D validation of both stages of the proposed methodology are presented in this chapter. The results of these experiments serve as a ‘proof of con- cept’ for the method. Without them it is unreasonable to assume the method would work in more realistic conditions or with velocity fields derived from real PC—MRI measurements. Starting with the resolution and embedding experiments described in Section 2.3.2, each additional experiment relies on the results of those before it. For this reason, this section is divided into three parts. The first two parts mirror the presentation of the experimental methodology in Section 2.3. Tabulated versions of the graphical resluts presented in these two sections are in Appendix E. The final section is dedicated to a discussion of the results. 3.1 Influence of Spatial Resolution and Embed- ding Strategy The examination of the influence of isotropic spatial resolution, represented here by the relative voxel size n, on the normalized relative RMS error, E, reveals that the relative pressure field estimation is second-order accurate, as shown in Figure 3.1. 43 0 10E ’3 e10".- Ito E -2 10 r f l —3 10 _3 _2 -1 10 10 10 77 Figure 3.1. Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of the the relative pixel size, 77 = Ax/2’l‘2. More importantly, Ex 1% for the resolution expected for MRI studies of the carotid artery, 7] = 0.05. Thus, it is not necessary to artificially increase the resolution of the velocity field with interpolation in the remaining experiments. Anisotropic resolution has a different effect depending on which direction is coarser in resolution. Above X301 = 10 and XFFT = 10, E differs between the two directions, see Figures 3.2(a) and (b), respectively. For X301 > 10 E is larger compared to the opposite, XFFT > 10, maximum E: 2.6% and maximum E: 1.3%, respectively. Below X30] 2 10 and XFFT = 10, the solver behaves identically to decreasing values of X30] and XFFT with a clear minimum E at x501 = XFFT = 1. The effect of embedding 0p in increasingly larger 0 is shown in Figure 3.3. If 0p is embedded in the center of .0, E decays monotonically towards an asymptotic value of 0.8% with increasing F. In addition, Figure 3.3 indicates that there is an error 44 E (%) 10_ 10 10 Xa Figure 3.2. Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of X0 = a/nconv. minimum (Ex 0.9%) at F of approximately 45% when 0F is located in the corner of .0 and at about 110% for the case when 0;: is located off to one side of .0. All in all, the effects of the embedding are relatively negligible (E< 1%) provided the embedding level is greater than 20%. A F of 25% is used for all other experiments presented here. 45 1.5 — Center - - -Off Center - 1.4 ----- Corner 1.3 ~, ' 0 50 100 150 200 250 F (%) Figure 3.3. Plot of the normalized relative RMS error, E, of the estimated relative pressure field for Couette flow with no noise added as a function of the embedding level P with a 20 x 20, which corresponds to 77 = 0.05, flow domain 01:. 3.2 Calibration and Performance of Proposed Noise Filters The Lagrange multipliers used for Filter 1 (01), Filter 2 (02 and 72) and the modified form of Filter 2 using the filtering terms of [2] (03 and '73) and [3] (a4 and ’74), are calibrated by optimizing the results for the worst noise conditions considered here (SNR = 5). The normalized relative RMS error for the estimated relative pressure field, E, is plotted against (1,1,2 and 7'2 in Figure 3.4(a,b). The optimal performance is obtained for a’l = 0.1 and (1'2 = 0.01,’7’2 20.53, which are the conditions used henceforth. Similar plots for the filters using the terms of [2], see (2.11), and [3], see (2.12), are shown in Figure 3.5(a,b). The optimal performances of these filters is 46 (a) Filter 1 calibration (b) Filter 2 calibration Figure 3.4. Plot of the average normalized relative RMS error, E, of the estimated rel- ative pressure field for the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain 01:, an embedding level 1" = 25% and SNR = 5, as a function of: (a) (1’1 for Filter 1; (b) at? and 7’2 for Filter 2, the minimum is circled. 47 -3 -2 —l 0 1 log (1’3 (a) (d/d.z)'2 calibration logy’4 —4 —3 .5 —3 log (1’4 (b) Viscous dissipation calibration Figure 3.5. Contour plots of the average normalized relative RMS error, E, of the estimated relative pressure field for the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain .01: and an embedding level P = 25% and SNR = 5, as a function of: (a) a]; and 73 for the filter with filtering term proposed in [2]; (b) “£1 and 72 for filter with the filtering term proposed in [3]. The minimum is circled in both figures. 48 obtained using as = 0.01,)f3 = 0.1 and d1, = 0.0001“1 = 0.1 for (2.11) and (2.12), respectively. 25 i l“. No Filter K‘s. ' - - Median Filter 20— ’ '- - -- Filter 1 X Filter 2 X PSDF “ .. .. .. 2 15~ \‘i (d/dx) $5 ‘ ‘ ‘ "- - ‘ Visc. Diss. v 0‘ ~ . ~ ~ N) 10 Ire \ \'\ \\\ h ----------------- J E E 's \\ SKETIE'W: - - -. .. .. .. “a? . - _ _LQ.‘ ‘ n. ‘ Ni! V luv-n un- - 1' "-7; :1 O i 1 5 10 15 20 25 30 Figure 3.6. Plot comparing the average normalized relative RMS error, E, in the estimated relative pressure field, for the initial Couette flow parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain .01: and an embedding level P = 25%, after applying the three filters pr0posed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimation with the field estimated with no filtering and three filters from the literature (PSDF [1], (d/dax)2 [2], Visc. Diss. [3]) over the range of SN R tested. The normalized relaitve RMS error, E, obtained before and after filtering is plot- ted against SNR in Figure 3.6. The tests were carried out using the initial set of experimental parameters with resolution of n = 0.05, on a 20 X 20 flow domain 0p embedded in a 25 X 25 computational domain .0. The embedding level P = 25% is deemed sufficient in View of the results presented in Section 3.1. For all filters, E decreases with increasing SN R, as expected. Moreover, the two physics-based filters proposed here outperform the median filter and the PSDF filter of [1] over the full range of SNR values tested. Filter 2 outperforms all filters tested throughout the 49 range of SNR tested. Both the median filter and PSDF filter become ineffective at low SNR, SNR = 10 and SN R = 7.5, respectively. The other filters all remain effec- tive through the SNR tested. For Filter 1, E decreases from an average of 14.3% at SNR = 5 to an average of 2.5% at SNR = 30. For the filter of [2], E goes from 10.4% to 2.8% for the same SNR range. For the filter of [3], E goes from 7.5% to 4.3%. For Filter 2, E goes from 5.3% to 2%. The dependence of the estimated relative pressure field, Pnum(:r, y), on the filtering of the noisy experimental velocity data with SNR = 5 is illustrated in Figure 3.7. Pixelated images of the magnitude of the normalized error distribution, e(a:, y), are plotted when no filtering is applied in Figure 3.7(a) and then for the median filter in Figure 3.7(b), Filter 1 in Figure 3.7(c), Filter 2 in Figure 3.7(d) and filters using the terms of of [2] in Figure 3.7(e) and [3] in Figure 3.7(f). The results of the PSDF filter are not included because the reduction in error is indistinguishable. The improved performance of Filter 2 is demonstrated by the clearly lower magnitude of 5(22, 3;) and the significantly lower relative RMS error (e) as compared to the results for the other filters (by factor of 1.35 w.r.t. the filtering term of [3], 1.84 w.r.t. the filtering term of [2], 2.7 w.r.t. Filter 1 and factor of 3.64 w.r.t. the median filter). The performance factor F is plotted against SNR in Figure 3.8 for all filters tested here with the initial set of Couette flow experimental parameters. Without filtering F is expected to be 1, F < 1 indicates that the filtered velocity field is closer to the true velocity field and F > 1 means the filtered velocity field is further from the true velocity field than the original noisy field. The two physics-based filters proposed here, Filter 1 and Filter 2, both have F values less than 1 for the entire range of SNR tested. The median filter and the PSDF filter both begin to push the filtered velocity field away from the true velocity field after an SNR of 5. The filter using the term of [2], (d/das)2, remains effective through the range of SNR tested, but F approaches one at SNR = 30. The filter using the term of [3] remains effective 50 (a) No filter: E: 23.6% (1)) Median filter: E: 19.3% (c) Filter 1: E: 14.3% (d) Filter 2: E: 5.3% (e) (d/dx)2: E: 9.8% (f) Visc. Diss.: E: 7.2% Figure 3.7. Sample images of the magnitude of the normalized relative error distri- bution, s(.v, y)|, for the estimated relative pressure fields obtained using the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain 01:, an embedding level P = 25% and SNR = 5 (a) no filter; (b) a median filter with radius 1 pixel; (c) Filter 1; (d) Filter 2; (e) a filter using the term of [2]; and (f) a filter using the term of [3] on the noisy velocity field. The normalized RMS error, E, is specified for each case. 51 4.5 - - -Med. Filter ’ , ’ ’ 4” -----Filter1 a’ ’ 35 —Filter2 I” ‘ PSDF ,’ 3- ~ .. “(d/dX)2 x’ --~ Visc.Diss. ” 2.5” I’ 4 LL. ”’ 9‘ "d... 2)- ” .. ,. a “I J I” " E E 1.5- ,’ .. ’ ’ ’ r 1'” , ’ I ~~~~~ - IIIII -.-:v;fld'n;-=-fl-T._f:.:,:r_7— ----------- q ass-"'7 0 l 1 1 l 5 10 15 20 25 30 SNR Figure 3.8. Plot comparing the average performance factor F to SN R for all filters tested using the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain 0F and an embedding level P = 25%. Note that the median filter, the PSDF method and the filter using the viscous dissipation term of [3] all become ineffective by SNR = 15. through an SNR of 15. The effect of filtering on the noisy experimental velocity data with SNR = 5 is illustrated ill Figure 3.9. Vector plots of the analytical velocity field over half of 0p are superposed with the velocity field resulting with the application of no filtering in Figure 3.9(a), median filtering in Figure 3.9(b), Filter 1 in Figure 3.9(c), Filter 2 in Figure 3.9(d), a filter using the term of [2] in Figure 3.9(e) and a filter using the term of [3] in Figure 3.9(f). The results of the PSDF filter are not included because the reduction in error is indistinguishable. The F value for the respective filtered velocity fields are indicated beneath the individual vector plots. The efficacy of Filter 2 at producing a velocity field that is closer to the true velocity field is evidenced by the difference in F as compared to the other filters (a. factor of 4.17 smaller w.r.t. Filter 1, 52 (a) No filter: F = 1 \ “9V ifl ..L‘.‘R\ ////I -..~\\\\ \ / "5 ¢:::--x\\\ «\\ it"" ,//j(,,. H” ill (c) Filter 1: F = 0.75 . I?“ c x “I“? Vf If If“:\:\"‘\\\ ( fifth/I ‘\ \\ J§ZZ¢::::-.EE:§§§\ I , K K \ ‘ \. ix ‘ l I ? ’V . i l "/v’t""" . \ eff/H” it “H" M (e) (d/dz)2: F = 0.64 Figure 3.9. Sample velocity vector fields of the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 X 20 flow domain 0F, an embedding level P = 25% and SNR = 5 illustrating the effect of (a) no filter (Vexp); (b) a median filter with radius 1 pixel; (c) Filter 1; (d) Filter 2; (e) the filter using the term of [2], see (2.11); and (f) the filter using the term of [3], see (2.12), on the noisy velocity field. The subscript f indicates the filtered field. The performance factor, F, is specified for each case. 53 ,. ,. w KR V I“ s5.\‘\ ”“9 '_' ‘ E ii‘\_ ”f/r/’ _‘ \\\\V\ ,x’fii/Vfl’iutuufi f”//” \\\‘ oft/11' ’ i ‘1'], I ill will .iil (f) Visc. Diss.: F = 0.62 a factor of 4.9 w.r.t. the term of [3], a factor of 5 w.r.t. the term of [2] and a factor of 5.22 smaller w.r.t. the median filter). 40 . ' ' fl 1., ---- No Filter 35 - - - Median Filter ‘ _ "\' """ Filter 1 30 "i — Filter 2 25 l' \i,\ s r, \ "0 ”is \. 15 I ".‘s \\\‘ ‘1 \_ ~ ~ ~ “-\_ 10* ‘.~,~'~" EK‘N 5\ -'-.-'-"'~'~.-,-.___:P O 1 4 l l 5 10 15 20 25 30 SNR Figure 3.10. Plot of the average normalized relative RMS error, E, in the estimated relative pressure field versus SNR for the more realistic Couette flow parameters, resolution 77 = 0.05 on a 20 x 20 flow domain 01: and an embedding level of F = 25%, comparing the results of the three filters proposed here to the result with no filter over the range of SNR tested. The second set of noise experiments using the optimized scaled filter parameters on the Couette flow mathematical flow phantom with a more realistically scaled velocity field show that the scaling of the filter parameters is appropriate, see Figure 3.10 for a plot of the RMS error E versus SNR. In this case the PPE solver does amplify the noise to a greater extent, E: 38%, than in the previous set of experiments, E= 24.3%. The filters remain effective in reducing the propagation of noise into the estimated relative pressure field reducing E to 21.9% using the median filter, 22.9% using Filter 1, and 6.1% using Filter 2. In these experiments the median filter out performs Filter 1 by 1% and continues to be effective until SNR = 20. Filter 1 becomes more effective 54 than the median filter as the SN R increases and remains effective through the range of SNR tested with E: 4% at SNR = 30 compared to 6.7% for the median filter and 6% without a filter. Filter 2 decreases E to 2% at SNR = 30. 1000 . x ”W No Filter 900 b X - - -Median Filter 800: \‘x ’ ,, ''''' Filter 1 ~-~ ->=.,.-\—‘-..-.--' —Fllter2 700. \\\ \‘\ - .. -(d/dx)2 2;; 600: \\ s‘ _- ~ - Visc. Diss. IT: '\ “ ____________ 500 i ". \\ 400 ~ 300 - 200 ~ 100 5 Figure 3.11. Plot of the average normalized relative RMS error, E, in the estimated relative pressure field versus SN R for the Poiseuille flow parameters, resolution 7) = 0.05 on a 20 X 20 flow domain 01: and an embedding level of F = 25%, comparing the results of the three filters proposed here to the result with no filter and two filters from the literature (d/das)2 [2], Visc. Diss. [3]) over the range of SNR = 5 — 20. The set of noise experiments using the optimized filter parameters on the Poiseuille flow mathematical flow phantom demonstrate that the solver has difficulty with noise in such simple flows, see Figure 3.11 for a plot of the RMS error E versus SNR. In this case the PPE solver does amplify the noise to a greater extent, E: 924%, than in the previous experiments with Couette flow, E: 24.3 — 38%. The median filter is entirely ineffective. The filter with the additional term of [3] based on viscous dissipation also reduces the error at SNR = 5 ((E= 765%), but this filter also becomes ineffective 55 before SNR = 10. Filter 1, Filter 2 and the filter with the additional term of [2] that minimizes the square of the first spatial derivatives of the velocity field remain effective through the range of SN R tested. For this flow, the filtering term of [2] is the most effective reducing E to 131% at SNR = 5 whereas Filter 1 reduces it. to 573% and Filter 2 reduces it. to 483%. 56 3.3 Discussion The investigation of the influence of spatial resolution and embedding level P on the estimation of the relative pressure field reveals that the PPE solution technique is capable of generating accurate relative pressure fields at spatial resolutions that are consistent with PC-MRI of small blood vessels. Based on the examination of effect of anisotropic resolution, the Poisson solver combined with embedding is affected by changes in the coarseness of the resolution and which direction is coarser. Greater coarseness in the direction of the FF T does not affect the accuracy of the estimated pressure field as much as coarseness in the solving direction when the resolutions differ by a factor of 10 or more. In situations such as these the Poisson solver can be adjusted to maximize the accuracy by applying the FF T in the direction of coarser resolution. It is clearly beneficial to use isotropic resolution as this maximizes the accuracy of the estimated pressure field. This fact is more important in three—dimensional velocity fields measured by MRI where it is more common to have anisotropic resolution due to the multi—slice acquisition of the spatial information (preferred to lengthier 3D data acquisition). The accuracy of the solution is affected by the value for 1",\especially when .017 is not embedded at the center of .0, and there exists optimal embedding levels that provide more accurate solutions. From a practical aspect however, these optimal embedding levels may require an excessively large computational domain 0, which then requires longer computational times and / or more computational power. Normalized RMS errors of the estimated relative pressure field on the order of 1% or less (5‘ S 1%) are achievable for realistic in-plane resolutions and embedding levels P > 20% in the absence of noise in the velocity input. For PC-MRI, realistic RMS noise levels are about 10% to 20%, corresponding to SNR = 5 to 10. In the initial set of noise tests, see Figure 3.6, the curves for 5‘ vs. SNR for no filter and the median filter cross at SNR = 10, which reveals that the median filter actually amplifies the propagation of the experimental noise into 57 the estimated relative pressure field for SNR > 10. Similarly, the PSDF filter of [1] begins to amplify the propagation of noise for SNR > 7.5. By implementing two physics-based filters on the experimental velocity field with a normalized RMS noise of 20%, the normalized RMS error in the estimated relative pressure field is reduced from E = 23.6% to 14.3% by Filter 1 and down to only 5% by Filter 2, which is quite encouraging. Implementing the filter of [2], see (2.11), reduces the error to 10.4%, and implementing the filter of [3], see (2.12), reduces the error to 7.5%. Filter 1 requires less computational time than Filter 2: 30 min compared to 75 to 90 min on a Dell Optiplex 755 with a 2.99 GHz processor, 1.96 GB of RAM, cache speed of 3 GHz, and cache size of 4 MB (Dell Computers, Round Rock, TX). Thus, the significant gain in accuracy provided by Filter 2 vs. Filter 1 is inversely proportional to their relative computational time, which is an expected trade—off. The filters using the terms of [2] and [3] are implemented in the same manner as Filter 2 and have the same computational time. Because Filter 2 as proposed here is more generates a more accurate relative pressure field estimation than either of these two filters with the same computational time, it is clear that they offer no advantage over Filter 2. The comparison of the F values of the filtered velocity fields is also revealing. Over the range of SNR tested the two physics-based filters outperformed all other filters tested. In some cases the effect of the filter on the velocity field can be inferred from Figure 3.6. For instance, it is apparent that both the projection onto the space of divergence free fields, PSDF, proposed in [1] and the median filter become ineffective in moving the noisy velocity field towards the true velocity field because their relative pressure field 5 curves cross the curve for estimation with no filter. From Figure 3.8 the F values for these two filters start at an average of F = 0.95 and F = 0.98, respectively, for SNR = 5 and quickly become greater than one. The results of the PSDF method do not agree with those published at similar SNR, F = 0.88 in [4] compared to an average of F = 0.95 here for SNR = 5. However, those experiments 58 were performed with a larger computational domain (32 x 32). The PSDF is fast (0.01 sec), produces a velocity field closer to the true velocity field than the median filter and is implemented using a similar Poisson solver to that used for the pressure estimation, but the resulting relative pressure field estimation is worse than that of a median filtered velocity field. The other two filters proposed in the literature have similar, though less severe, effects on the velocity field. Based on Figure 3.6 both of the physics-based filters proposed here, Filter 1 and Filter 2, should generate velocity fields that are nearer the true field because their E curves stay near or below the expected error for the range of SNR tested. Again the results in Figure 3.8 agree with this expectation. The performance of Filter 1 decreases from an average of F = 0.75 to F = 0.69 over the range of SN R tested while the performance of Filter 2 increases from an average of F = 0.18 to F = 0.25. It is more difficult to determine the effectiveness of the filters using the terms of [2], (d/d.r)2, and [3], Visc. Diss. based solely either on E or F. Over the lower portion of the range of SNR tested, both of these filters generate velocity fields that better estimate the relative pressure field than the noisy velocity field as compared to PSDF, median filtering and Filter 1 according to Figure 3.6. The F value at SNR = 5 for these filters are on average F = 0.56 and F = 0.46 for the filtering terms of [2] and [3], respectively, which agrees with the E interpretation. According to the F values of these two filters, they start becoming less effective than Filter 1 at SN R z 10 whereas according to E they become less effective between SN R = 15 and 20. The F value of the filter using the term of [2] remains less than one only through SN R = 25, but its E curve remains below that for no filtering throughout the range of SNR tested. The filter using the term of [3] becomes ineffective at SNR x 12.5 according to its F value which is before the SNR that its E curve crosses that of Filter 1 and near the SNR where it crosses the curve of the filter with the term of [2]. At SNR = 25 its E curve crosses the curve of estimation without a filter and its F = 2. 59 The second set of noise tests demonstrate that the method and the optimized scaled filter parameters do transfer to more realistic conditions. While the RMS error in the estimated relative pressure field increases by 60%, from 23.6% to 38%, between the initial experiments and these experiments when no filter is applied to the velocity field at SNR = 5, it only increases 20% from 5% to 6%, when Filter 2 is applied. This indicates that the proposed scaling is appropriate. However, the results of applying Filter 1 seem to indicate the opposite. At SNR = 5 the median filter actually outperforms Filter 1 and the E increases from 14.3% in the initial set of tests to 22.9% in the second set. This increase of z 60% mirrors that of the case with no filter, but it does not mirror the increase in the value of (11 which increases from an average of 0.65 to 89. More likely the additional filtering term of minimizing the Laplacian of the velocity field is an effective way of combating the propagation of noise from the velocity field to the estimated relative pressure field. The approximate optimization of a’l used in Filter 1 is also not as precise as that for 015 and 7’2 used in Filter 2. Fine-tunning the value of az'l may increase its accuracy. The median filter remains effective over a greater range of SNR in these experiments decreasing E through SNR = 20 as opposed to SNR = 10 in the initial set of tests. Both Filter 1 and Filter 2 remain effective through the range of SN R tested as they did in the initial set of experiments. The times to execute Filter 1 and Filter 2 in these experiments remains the same as it was in the initial experiments. The final set of experiments with a mathematical phantom, Poiseuille flow, demon- strate that the proposed filters reduce the error in the estimated relative pressure field when compared to the estimation without filtering in a second flow system. The mathematical simplicity of the Poiseuille flow belies the difficulty encountered by the proposed method in estimating the relative pressure field from the noisy ve- locity field. There is a drastic increase of E in the estimated relative pressure fields with and without first filtering the velocity field. Estimating a linear gradient using 60 Fourier methods is difficult because of the practical limitation to use a finite rather than infinite representation. The solver is capable of accurately estimating the rel- ative pressure field with no noise in the velocity field at these resolution conditions (E = 0.2%). The additions of noise to the velocity field, however, allows the pressure estimation to converge to a periodic field which is preferred by this type of solver. It is worth noting that the filters, with the exception of the filter using the term of [3], follow the same trend of accuracy improvement as was seen in the Couette mathe- matical phantom experiments. The median filter is effective initially, but becomes ineffective before SNR = 10. Filter 1, Filter 2 and the filter using the term of [2] all reduce the error through the range of SNR tested. Unlike before, the simpler of the two filters proposed here, Filter 1, is either as or more effective at higher SNR than Filter 2. Also unlike the Couette experiments, the filter using the term of [2] is the best performer throughout the range of SNR tested. Interestingly, the filter using the term of [3] increases E compared to the case with no filter at higher SNR which is a behavior not seen in the previous experiments. Because the error was so large across all cases in using this phantom, it is difficult to tell if any bias was introduced by the combination of filter and mathematical phantom. The fact that Filter 2 continues to be one of the better performers indicates that it has the potential to be a useful filter in cases where the flow physics are less well understood. These ‘proof of concept’ experiments to determine if this type of two—stage method is worth investigating further are encouraging. PC-MRI resolution conditions do not greatly hamper the PPE solver’s ability to estimate an accurate relative pressure field in zero noise conditions. The level of embedding, F, and the position of 0p within .0 do factor greatly into the accuracy of the solution, and the level of embedding required to generate a solution that is within 1% relative RMS error of the true relative pressure field does not require .0 to be excessively larger than 0F. Based on the two sets of noise testing, the scaling of the filter parameters is appropriate. Additionally, the 61 Filter 2 proposed here outperforms all of the other filters tested in the more complex mathematical phantom. Filter 1 is faster than Filter 2 and the filters using the terms of [2] and [3], but it produces a less accurate velocity field than Filter 2 and struggles in conditions a with more realisitically scaled velocity field. The results of this 2D validation show that the two—stage method in general, and using either one of the two physics-based filters proposed here in particular, is a reasonable approact to estimating the relative pressure field from noisy velocity fields. 62 CHAPTER 4 Conclusions Engineering is a discipline whose goal is to find the best path from problem to solution. In many cases this is an iterative process where the ideal path is found through a series of small improvements. Here the problem being considered directly is estimating the relative pressure field from noisy velocity measurements, but the underlying goal is to provide a diagnostic tool. The work presented here represents one iteration toward the best solution for both of these aspirations. The previous iterations include the minimally invasive blood flow assessments reviewed in Section 1.3 and the various algorithms for estimating the relative pressure field reviewed in Section 1.2. This step adds a noise filter between the velocity measurements and the pressure estimation to improve the resulting relative pressure field. The results of this work show that the inclusion of either one of two proposed physics-based filters as well as two filters from the literature improve the accuracy of the estimated relative pressure field. In this chapter the insight into the problem gained in this iteration of estimating relative pressure fields and the future of the proposed method are presented. 63 4. 1 Insight It is important to understand how an algorithm will behave in the presence of various complicating factors before it is used to interpret real-world measurements. The mathematical phantom experiments were useful in elucidating the complexities of the problem of estimating the relative pressure field on an irregular domain. The effect of spatial resolution on the algorithm is the first issue that many numerical methods are required to overcome. The results of the presented experiments have demonstrated that the pressure estimation algorithm is adequately capable of resolving relative pressure fields accurately at the relatively coarse resolution achievable in clinical PC- MRI. However, it is recommended that isotropic resolution be used when possible. Embedding the irregular domain as done here is useful for simplifying the application of the boundary conditions, but its effect on the accuracy of the solution has not been published. The results presented here show that using a computational domain that is 25% larger than the major dimension of the irregular physical domain produces suitable results in noise-free conditions. These results provide the necessary level of knowledge of the behavior of the algorithm to then proceed with tests using noisy velocity fields. Though the resolution and embedding tests discussed in Section 3.1 demonstrate that the pressure algorithm is accurate at the low resolution expected for PC-MRI, the noise tests are equally demonstrative of the inaccuracies that are introduced by noisy velocity measurements. The first stage of the proposed method, namely removing noise in the velocity field using one of three proposed filters or filters described in the literature, proves useful in reducing the negative effects of noisy data. Of the tested filters, the filter that simultaneously applies in a least-squares sense the continuity equation for an incompressible fluid, a constraint on the variation between the filtered velocity field and the noisy one and an additional constraint that limits the size of the Laplacian of the velocity field, Filter 2, performs the best. It succeeds in reducing 64 the error in the estimated pressure field from an average of 23.6% with no filter to an average of 5.3% with filtering at a realistic resolution in the worst conditions tested. This non-negligible improvement was also demonstrated in a second set of experiments with different flow parameters. The trend was also seen in a set of experiments using a second mathematical phantom. By comparing these results to those of Filter 1, which does not apply the additional term and improves the pressure estimation to an average error of 14.6%, it is evident that removing noise in the Laplacian is useful. Other second terms such as those proposed in [2] and [3] when applied using scaling factors like those proposed here are also effective in mitigating the pollution of the estimated relative pressure field with noise from the measured velocity field. However, accuracy is not the only factor to be considered when comparing the filters, the time required to complete the filtering and estimation of the pressure field is also of importance. At this point, it is uncertain how accurate the pressure field needs to be to fulfill the eventual goal of providing a useful diagnostic technique. For this reason, it is not prudent to ignore the capabilities of Filter 1 which is faster to implement than filters with additional terms. Using a physics-based approach in the design of the filters is an obvious improve— ment over ad hoc methods like median filtering. An understanding of the principles of fluid mechanics that govern blood flow in vessels like the carotid artery is necessary for these approaches to be workable. Minimizing the continuity equation for a noisy velocity field results in a velocity field that satisfies the conservation of mass, but there is no guarantee that the resulting field is still relevant to the physical system. Even though they are polluted by noise, the measured velocity fields still represent a system that satisfies all governing equations. When removing noise from this field the filtered field should not stray far from this starting point. From these physical insights the cost function of Filter 1, see Section 2.2.1, was written. Filter 2 adds an additional constraint that is designed to allow for a more accurate characterization of the effect of viscosity. Though it is not directly based on a governing equation, the minimization of the Laplacian of the velocity field is a method of removing noise from the flow field that directly affects the estimation of the viscous terms in the NSE. It is reasonable to assume that the velocity field does not contain discontinuities or large spikes at this level, so this additional constraint seeks to eliminate them. Other terms can also be applied if a priori knowledge of the flow field or experience dictates physical constraints separate from those tested here. The difficulty encountered by the method in estimating the relative pressure field from a noisy velocity field with no filterng in the Poiseuille flow experiments underscores the importance of filtering stage. The scaling of the filter parameters proposed here appears appropriate and is a useful way to avoid reoptimizing the parameters for each new application. Horn these physical principles, as well as the results of the resolution and embedding experiments, this two-stage methodology is proposed for estimating the relative pressure field from a velocity field at conditions that are realistic for clinical PC-MRI experiments. 4.2 Future Work This work endeavors towards providing the first step of a noninvasive technique for the diagnosis of the physiological relevance of arteriosclerotic stenoses in the carotid artery. To that end, a two—stage method is proposed to estimate the relative pressure field from a noisy velocity field. While the velocity field by itself does provide enough information to calculate such potentially important parameters as wall shear stress and oscillatory shear index, the relative pressure field is also necessary before any serious investigation into the dynamics of the vessel wall can be investigated. It is these dynamics that are expected to provide information about the danger associated with a particular stenosis. The extension of the two—stage method proposed here to the evaluation of the dynamics of the vessel wall is an obvious next step. 66 Before that step can be taken some refinement to the details of the method are desired. The time required to implement either of the two physics-based filters pro- posed here is one area of concern. Different minimization algorithms, rewriting the functions in a lower-level computing language like C or Fortran and faster computers are all possible solutions to this problem. Another area with potential for improving the method concerns the parameters used in the implementation of the filters. To this point the parameters are only roughly optimized. It is unlikely that there exists a set of exactly optimized parameters that will provide ideal results independent of the SNR. However, there is the possibility that further optimization of these param- eters could lead to better results than those presented here especially in the case of Poiseuille flow. It is possible that some sort of experiment with a physical phantom will be necessary to determine the best set of parameters to be used. This potential for improvement does not preclude the discussion, in general terms, of the future use of the proposed method that follows. Experiments with 2D mathematical flow phantoms are presented here, but blood flow in the carotid artery and the dynamics of the vessel wall are 3D and time-varying. It is difficult to construct a mathematical phantom with an analytical solution for both the velocity field and pressure field that is rich in all three spatial dimensions and as well as the temporal dimension. Thus, testing of the method proposed here is difficult to implement in higher dimensions. However, the time-varying velocity field is not a factor in the design of either of the two proposed physics-based filters and its effects are only present in the acceleration term in the N SE. Additionally, the pressure field is calculated relative to its own spatial average not the time average. Therefore, once the time derivative of the velocity field is estimated, the 3D pressure field that represents each time step can be calculated as a stand-alone system. Once the method is expanded and validated in 3D, tests on a simple time-varying flow, such as the Womersley flow, can be used to determine the effect of time the derivative on 67 the pressure estimation which allows for the focus to be shifted to the determination of the characteristics of the vessel wall. Interaction between the wall and the blood flowing through a vessel will have an impact on the estimation of the pressure field. The experiments carried out in this work did use flow systems with moving boundaries, but not systems with changing volumes. While it is reasonable to assume that the volume of a segment of a blood vessel will be unchanged between cardiac cycles, it is not reasonable to assume that it will remain unchanged during a cycle. If the properties of the vessel wall are to be understood, some coupling between the methodology presented here and the appropriate tissue mechanics models will be required. This is the next major step in developing the minimally invasive tool that is the underlying desire of this work. 68 APPENDIX A Finite Difference Schemes Finite difference schemes use Taylor series expansions to approximate the derivative of a function at a point. In this work there are no mixed derivatives so this appendix focuses on the derivation of schemes for derivatives in a single direction. The Taylor series expansion of a function f (.13) about the point f (.27 + A33) is OO . ’ _ (A2)" din) f f(.i: + An) _ 112:0: n! F”) 1, (A1) If f is evaluated on a discrete grid with a fixed spacing of As: between points, then the index notation i :l: n, where n = 0, 1, 2, . . . ,oo, can be used to expand the first few terms of the summation of (A1) at the point (i + 1) _d_f (Am)2 (12_f , (AI)3d_:f +A( _)Ilr“d_:f A similar expansion for the point i — 1 yields may (A. >2 .12) (M3 d3f (A. >4 d‘f 5 f._1)=f.-— .d_—,.+ 2 8,7,1.— 6 d—,,,l.+ ,4 .1 —,I.+ 0(Arv) (A3) By subtracting (A.3) from (A2) and rearranging terms an approximation of the first derivative at point i is possible using the value of f at the locations i + 1 and i — 1 fi|'_f(i+l)— f(i—l) dzz— 2A2: +O(A:L‘2) (A4) 69 Similarly, by adding (A2) to (A3) an approximation of the second derivative of f at point i is possible using the value of f at i + 1, i and i — 1. (Rffl_flHU-Zfl+fWJ) E5l, _ (Ax)? + 0(Ax2) (A5) Both (A4) and (A5) are second-order accurate central-difference schemes. The second-order accuracy stems from the expected magnitude of the first term of the Taylor series that has been which is of the order of the size of the step in :7: squared. First-order accurate schemes can also be derived from (A2) and (A3) df __ f(i+1) _ fi fi._fi-flau 8;]; — T + 0(ACII) (A.7) These are a forward scheme, (A6), and a backward scheme, (A.7). Schemes like these are of use on the boundaries of a discrete domain. If second-order accurate forward and backward schemes are desired or second derivatives are required, then more points are necessary. The expansions for points (i + 2) and (i + 3) are d_f +A<2 )42‘47 (242)3d3f (242)4d4f 5 f(i)=+2 fi+2A$ d—- 2 Egg—l2 __6—003—3" TEI'I+O(A‘T)(A°8) and df +(3A222)2 d2 f (_.__3A6i)3 d3 f (_-__3A.i:)4 d4f f(i+3)—fi+3A-r _li+ 7jli+ 0130135) (A-Q) (1.13 ——d.i:3l4+ 24 ___d—ffll’fl’ respectively. A second-order accurate forward scheme for a first derivative is found using (A2) and (A8) df _ _f(i+2) + 4f(i+1) — 3f,- _ . _ , 2 d1: 2A2: + 0(A.r ) (A.10) These points are also required to find a first-order accurate forward scheme for the second derivative d_2fl_ f(i+2) _2f(i+1) +fz' drr2 z — (Ax)2 + 0(Ai) (All) 70 The point approximated in (A9) is necessary to find the second-order forward scheme of the derivative in (A.11) d2)” |_ f(i+3) + 4f(i+2) — 5fi+1 + 2ft 2) (11:2 z — (Air)2 + 0(Aa: (A.12) A similar approach is used to find backward equivalents of (A.10), (All) and (A12). The expansions about points (i — 2) and (i — 3) are (1 2A d2 2A 3d3 d4 f( -2)=f.-— 2—Axd—f-+ ( 2") —dei-( 6‘” -d—.fl.-+ (,f) d——f|io:5+(Ar)(A-13) and d_f (A 0432);, (3A )3331 (3A 3)4d4f, 5 as): f.— 3A 0,, + 2 33,).— 6 731.324 FI.+0(A3)(A.14) respectively. The backward scheme equivalents of (A.10), (All) and (A12) are df f(i_2) — 4f(i_1) + 3ft Elf = 2m + 0(Ax2) (A15) (12 fi_2 _Zfi— +fz' ,3. = < > at v ..( .3) (A...) and (12 —fi_ +4f27_ — 5fi— +2fi 35:93 = ( 3) ((3)2 ( 1) +O(ALE2) (A.17) For those special cases where the proximity of 0.01: prevents the use of any of the above approximations, a central difference approximation that allows for variations in A57: is required. To start, substitute (A.4) into (A.2) f(i+1) f(i—1)+(A$)2 (12f|_+ (A$)3§3_f 2A3: 2 d932’+ 6 (13:3 f(i+l) = f; 'l' Ail} [i + 0(Al‘4) (A.18) Rearranging yields a new expression for f(i+1) in terms of f, and f(,-_1) 2d2f (A303 d3f _ . __ . 4 dzr2l’ 3 (1333], + 0(Aa: ) (A.19) f(i+1) = sz' — f(i—1)+(A$) For the proposes of this derivation it is convenient to once again use (:rinAr) instead of (i :i: n). Let us now assume that the grid spacing in the direction of (:1: + Ax) is 71 different from the spacing in the direction of (.‘L’ — A17), denoted as A$+ and Ar., respectively, and both of these spacings are different than the regular spacing on the rest of the grid, Ax. Rewritting (A2) using this new notation, solving for the second derivative and substituting (A.4) yields {IQ—f 2 f(1'+AI) _ f(:r.-A.r.) )] (A20) (1x2 (A$+)2 [f(1+A:r+) f1? l‘+ ( 2A1: A new expression for f(1‘+ A1.) in terms of f(I+A$+) and f($_A:,I) is found by substi- tuting (A20) into (A.19) 2(11—(A)?n+(AYfIIIm(AA—1mm} 1+3— f (1+Ar) = (A21) Following a similar procedure, a new expression for f(J:—A:r) in terms of f(:r—A:I:_) and f(;c+A1‘) iS 2((1- (A1 A) I IIIA-oII) _ 1+3Ai l‘_ (A22) By substituting (A22) into (A21) and vice versa new expressions for f(1:+A;r) and . . A A f(:r—A:1:) in terms of the ratlos R+ = Zfi and R_ = mg: f(:r.+Aa:+) and f(1,_A$_) are obtained R f(1‘.+A;r) = 2{ [1-(R+)2+ 1:12.1(1— 133)] far. + (3+)2 f(;r+A:r+) +R+_1(R—) )2f(x— Ax- l} ()1+R+)1[ (R+—1)(R—-1)i}_l (A23) 1+R_ (R++1)(R_+1) and Substituting (A23) and (A24) into (A.4) and (A5) yields approximations for first and second derivatives in areas where the schemes that assume a fixed grid are not possible. When these schemes are applied in any of the algorithms in this work it is assumed that. Ami = 0.5Aa: and that the velocity at that point is zero. This allows for a rudimentary approximation of the derivatives in virtually all regions of 0p. 73 APPENDIX B Fourier Poisson Equation Solver The solution of the Poisson equation in step 4 of Section 2.1.2 is performed using the solver based on fast Fourier transforms (FF T) presented in [5]. The solver is direct and more efficient than the Gaussian method. The two different versions of the solvers used in this work will be worked through in the following order: 2D problem with Neumann boundary conditions followed by 2D problem with Dirichlet boundary conditions. B.1 2D Neumann Problem Following the general algorithm laid out in Section 2.1.2, for a 2D computational domain .0 the first step is to incorporate the boundary conditions b - ft on 8.0 into V . b, see Figure 8.1 for an illustration as well as the definition of the boundary names. On the respective boundaries the boundary conditions are (VPfi 0,]- = (bT)0,j (VP'fiMIj = (bit-)Mj I ’ ‘ B.l (VP-72),“) = (bi/hp ( ) (VP'fihw = (by)i,N where bx is the :c-cornponent of the discrete form of the b and by is its y-component. Using the second-order finite difference scheme (AS) in the center of {2, the discretized 74 j=0 east j: i=0 9 th south y nor x i =M west Figure B.1. A schematic illustrating the computational domain .0 and defining the names of the boundaries and indices used in the direct solver of the Poisson equation in two-dimensions. version of V2P is P(z'-+1)Ij - 2PM + P(-2:-1)Ij PAW) — 2PM + Paw—1) (Ave)2 (Ay)2 = (V ' b)i,j = 9m (B?) By using the same discretization on the east boundary then for points i = 0 and j: 1,2,...,(N—1) P14 — 2130,,- + P—IIj P0I(j+1) ‘ 21003“ + POM—1) (Ax)2 (Ay)2 = (V ' b)0,j (33) However, the point P_1,j does not exist in this computational domain. To eliminate this point from (33), the Neumann boundary condition are applied as prescribed in (31) discretized using (A.4) Plaj — P_13j _ Substituting (B.4) into (33) yields the expression 2(DU—2100i P0(j+1)—2P0,j+P0(j—1) 2 i 1 7 ’ = V.b .+_ b , = , 85 (A113)2 (Ay)2 ( )0,] ASII( 1‘)0,] 90,] ( ) 75 Similarly for the west boundary for points i = AI and j = 1,2,. . . , (N —- 1) -2PM,j+2P(M—1),j + PM,(j+1)‘2PMIJ‘+PMI(j—1) (Arc)2 (Ag)2 2 . = (V ' b)1t[,j — A; (bI)1\[,j = 9M,j (36) the south boundary for points i = 1,2,. . . , (M — 1) and j = 0 (Av)2 (Ay)2 P(i+1),0 — 213210 + P(i—1) 0 2pm - 2PM) 2 = (V ' bhyo + E (bx/>130 = 9230 (8'7) and the north boundary for points i = 1,2,. . . , (M — 1) and j = N P(-z'+1),N '2Pz‘IN+P(2i—1),N + ‘2Pz'IN+2PiI(N—1) (A32)2 (Ag)2 2 = (V ' WAN — A; (bl/by = 92',N (38) In the corners, boundary conditions apply in both directions, but the same method is used to incorporate the boundary conditions as on the sides of the domain. The southeast corner (i = j = 0) of .0 has contributions from the south. and the east boundary conditions 2PM) — 21301) 2P0,1 — 213030 (Ass)2 (Ay)2 2 2 = (V ' b)0,0 + $011203 + E (by)0,0 = 90,0 (B9) The remaining corners proceed in the same fashion. At i = 0 and j = N 2P1,N - 2P0,N + ‘2P0,N + 2P0,(N—1) (Arc)2 (Ag)2 2 2 = (V ° b)0,N + E (brlow — A—y (by)0,N = 90,N- (310) At i = .M and j = 0 -2PM,0 + 2P(M—1),0 + 2PMII — 2PMIO (Ax)2 (Ag)? 2 2 = - — — b — b , = 8.11 (V b)M,O A,“ I)M,0 + Ag ( y)M,o 9M,0 ( ) 76 Ati=Mandj=N —2PM,N + 2P(M—1),N ‘2PMIN + ZPAIIW-l) (Ax)2 (Ag)2 2 2 =(V'b)M,N — A—$(bx)M,N—A—y(by)m,gv=9M,N (312) The next step in the solver is to use a FFT to decouple (8.5). This requires a cosine transform of the form 2 M 7r I ll . gkfij = IT! 9in cos “CM (B.13) i=0 II where the indicates that the first and last. terms are multiplied by one half. The inverse of (B.13) is 9233': :Z—[OI’ g) j cos ik— (8.14) The transform of (8.14) can also be applied to the discretized pressure field A! A . 7r P'iJ = Z HPkJ cos “CM (B.15) k=0 Substituting (B.13) and (B.15) into (B2) yields M 1 ,, A 7r Artfg Pk] [ +1k—-2coszkfi +cos(i—1)kM] + M P I I 2P -+ P. -_ :0” PMj+1) m2 M] 1) cosik— = :0 (Ay) M :Z—O” 9;, j cos ik— (B.16) The first term on the RHS of (B.16) can be simplified using trigonometric identities cos(A + B) = cosAcosB — sinAsinB and 2sin2A :1— cos 2A, to k—7r—— — — —=— — — . cos(i +1) AI 2cosikAI+cos(i 1)kM 45in2 k2AICOSikM (B17) 77 Substituting (B.17) into (316) and combining the sums on the LHS yields Ag 2 - . 21::I/ y)2 [3% 010+” —(2 "l” dig—sin 2193—11) kaj + Pk,(j—1)] COS lkfi ' II~ .1 7T = z} ng- coszkfi (B.18) Now apply the orthogonality of the basis function cos iTls [to yield the tridiagonal system * 2 ,I 2 7r * * 2- Pk,(j+1) — (2 + 4X s1n km) Pk,j + Pk,(j——l) = A3] 913,} (8.19) where X = Ay/ A13. On the south (j = 0) and north. (j = N) boundaries (B.13) and (8.15) to (8.7) and (38) respectively are applied yielding AI 1 ,, A ‘ , 77 (Ax)2 lg) PM) [cos(i + ”161111 — 2cos ikfi + cos(i— INC—jg] A A A P. — P I! k,1 LO +2]; ( Ay2 ) COSZk-M- A, 71 _ ”I - _ — E: 9“) cos zklM (B20) and 1 AI n * , I, _. _ _ _ _ _ (Ar)? [2%) Pk,N [LOb(l+1)k7r\I 2cosilc7r M +cos(i 1)k Ml Pk ,N + Pk ,1(N~ ) 7r II - . +2 2]: (_ Ay2 cos MEI— — —Z” 9!: N COS 211— (8.21) Following the same simplification as for (319) the boundaries become 2Pk’1 — (2 + 4x2 sin2 k2—71l1) Pkg = Ay2§h0 (B22) forj = 0 and —(2 + 4x2 sin 21:27:11) PM, + 2Pk,(N_,) = AyQQkVN (B23) 78 for j = N. The result of the transform is (N + 1) tridiagonal systems of (M + 1) equations which can be solved in a variety of standard ways. For the problem with Neumann boundary conditions, however, an additional step is required before the system can be solved. In the case of the Neumann problem, a least-squares solution to the Poisson equation is sought. To complete the solution the problem the discrete form of Green’s theorem is applied to this system [W N M N Z (hr/(LN) ‘ bum») + Z (bum) — bro») = Z Z 914 (B24) i=0 i=0 i=0 j=0 Noisy velocity data and embedding mean that this system may not satisfy (B24) initially and, therefore, that a solution does not exist. The RHS of (B24) is per- turbed slightly so that the resulting new system gm- does obey Green’s theorem. The solution to this system will be a least-squares solution to the original system 9M if the perturbation is small compared to 913]" The perturbed system is defined as M N i=0 j=ocicjgiJ Al N 21:0 23:0 Cicj gig = 9233‘ - (B25) Table 31. The weights for the perturbation required to solve the 2D Neumann Poisson equation. 6: 0.5 1 05 j 0 1,2, .,N—1 N c] 0.5 1 05 where c,- and cj are weights whose value depends on i, j, see Table B.1. The Fourier coefficients are solved for using the (M + 1) tridiagonal systems of (N + 1) equations defined in (B.19), (B22), (B23) and gm- yielding PM. As mentioned before one of 79 the tridiagonal systems will be singular. The output of this system can be set to an arbitrary value, in this work it is set to zero. The final step is to transform back into physical space using (B.15) to obtain PM. B2 2D Dirichlet Problem The solver for the 2D Dirichlet problem starts with the same discrete version of VQP as the Neumann problem, but as the boundary conditions are different, so is their incorporation. For this problem (B2) is the same, but the boundary conditions are P0,j = (Peast)j PALj : (pwest)j (B26) Pi,N = (pnorth)i Pi,0 : (psouth)z' where Peasta pwesta Pnorth and Psouth are discrete functions that define the value of PM on their respective boundaries. At 2' = 1 and j = 1,2,. . . , (N— 2), (B2) becomes P2,j — 2P1,j + P0.j P1,j+1" 2P1,j+ P1,(j—1) . = V - b - B27 (Ax)2 my)2 ( )1” f ) Applying the appropriate boundary condition from (B26) and rearranging P2j - 2P1 j P1(j+1) — 2P1,j + P1 (j—l) 1 ___’—__’- + ’ ’ = V . b . _ — . . = ' (A513)2 (Amg ( )1,] (A$)2 (Feast)] 91,] (B28) The other boundaries follow in the same manner. For i = (M — 1) and j = 1,2,...,(N—2) ‘2P(M—1),j + P 86 Substituting (Cll), integrating and solving (C12) for C yields C = 332% = 111%; ((3.13) The analytical velocity field is now My) = 31% (y2 - by) (C14) The analytical pressure field is found by integrating (C13) w.r.t. .1; 10(1) = -— ”QC" (.1- — 230) ((3.15) h3 where .270 is an arbitrary datum. 87 APPENDIX D Experimental Parameters The experiments determining the robustness of the PPE solution method to reso- lution, level of embedding and location of {21: in .0 used a set of parameters in the Couette flow mathematical phantom that were not intended to match physiological conditions, see Table D1 Instead it was intended to be a purely mathematical sys- tem. As such, the parameters were chosen for convenience. Using the values from Table DI the Rec of this flow is 2.10 x 106 which is far beyond the transition to tur- bulent flow in this system. Therefore, the tests performed using this set of parameters serve as a ‘proof of concept’ for the method. In the second set of noise experiments an attempt was made to have a more realistically scaled velocity field. This was done by matching Rec and Rea as defined in Section 2.3.3 as closely as possible. Using the conditions from the calculation of the W0 from Section 2.1.1 (p = 1040kg m’3, u = 0.004Ns, Da = 0.01m) and assuming an average velocity of U0 = 0.1 m s‘l, ReA % 260. The parameters in Table D2 closely approximate these Viscous conditions leading to Reg = 235. For the set of Poiseuille flow of experiments, similar values to those chosen for the second set of Couette flow experiments. The volumetric flow rate QC was chosen to approximate the average flow rate in the carotid artery and normalized by the average 88 radius of the artery to yield the desired units and velocity scales. This is based on the average velocity and approximate cross-sectional area. The parameters in Table D3 were used to generate the analytical velocity field and relative pressure field. Table D1 Parameters used in the Couette flow mathematical phantom in the reso- lution, embedding and initial noise experiments. 53min (m) —1 firmax (In) 1 ymin (m) "1 ymax (fl) 1 r1 (m) 0.25 7‘2 (m) 0.79 wl (rad/s) 7r w2 (rad/s) 27r u (N s m-2) 0.656 x 10-3 p (kg m—3) 992.2 7] (-) 0.05 Table D2. Parameters used in the Couette flow mathematical phantom noise exper- iments with more realistic viscous terms. mm (m) —6.2 x 10—3 :rmax (m) 6.2 x 10-3 ymin (In) -6.2 X 10_3 ymax (m) 6.2 x 10—3 7‘1 (m) 2 X 10—3 3 (m) 5 x 10—3 w; (rad/s) 3671' (.02 (rad/s) 187r p (N s m’2) 4 X 10—3 p (kg 111-3) 1040 7) () 0.05 89 Table D3. Parameters used in the Poiseuille flow mathematical phantom noise ex- periments. mm (m) —6.2 x 10-3 xmax (m) 6.2 x 10-3 ymin (m) -6.2 x 10—3 ymax (m) 6.2 x 10—3 h. (m) 10—3 QC (m3 is) 1.57 x 10—3 p (N s m‘2) 4 x 10—3 p (kg 111—3) 1040 7] (—) 0.05 90 APPENDIX E Tabulated Experimental Results This appendix contains the tabulated results from Sections 3.1 and 3.2. The tables are labeled with their corresponding figures. Table E1 The normalized relative RMS error, 5, results of the isotropic and anisotropic resolution tests (see Figures 3.1 and 3.2). 7) 5" (%) xa a = 71301 5 (%) a = 77301 F? (%) 0.0007 0.0032 0.1250 0.0217 0.0217 0.0014 0.0028 0.2500 0.0246 0.0246 0.0056 0.0087 0.5000 0.0222 0.0222 0.0112 0.0243 1.0000 0.0183 0.0183 0.0223 0.0966 2.0000 0.0430 0.0430 0.0255 0.1278 4.0000 0.1485 0.1485 0.0298 0.1930 5.3333 0.1520 0.1520 0.0357 0.3249 6.4000 0.2139 0.2139 0.0446 0.3943 7.1111 0.2563 0.2563 0.0595 0.5899 8.0000 0.3650 0.3650 0.0893 1.8457 9.1429 0.5880 0.4227 10.6667 0.6379 0.6379 12.8000 1.2395 0.7797 16.0000 2.6315 1.3225 91 Table E2. The normalized relative RMS error, E (%), results of the embedding tests (see Figure 3.3). F (%) 0.0000 8.6957 17.391 26.807 34.783 43.478 60.699 69.565 Center 1.4908 1.1014 1.0029 0.9595 0.9223 0.8923 0.8714 0.8575 Off Center 1.4320 1.0835 1.0024 0.9875 0.9833 0.9755 0.9639 0.9515 Corner 1.2316 1.0178 0.9584 0.9269 0.9187 0.9260 0.9412 0.9581 F(%) 78.261 86.957 95.652 104.35 113.04 121.74 130.44 139.13 Center 0.8487 0.8430 0.8390 0.8360 0.8339 0.8322 0.8308 0.8296 Off Center 0.9407 0.9329 0.9286 0.9275 0.9288 0.9323 0.9369 0.0.9423 Corner 0.9733 0.9851 0.9933 0.0.9983 1.007 1.0014 1.0009 0.9996 F(%) 147.83 156.52 165.22 173.91 182.61 191.30 200.00 208.70 Center 0.8286 0.8277 0.8268 0.8261 0.8255 0.8249 0.8243 0.8239 Off Center 09481 0.9536 0.9590 0.9640 0.9684 0.9723 0.9756 0.9784 Corner 0.9979 0.9960 0.9940 0.9920 0.9900 0.9882 0.9865 0.9849 F(%) 217.39 226.09 234.78 243.48 252.17 260.87 269.57 278.26 Center 0.8234 0.8229 0.8225 0.8222 0.8219 0.8916 0.8213 0.8211 Off Center 0.9808 0.9828 0.9844 0.9858 0.9869 0.9877 0.9885 0.9890 Corner 0.9834 0.9821 0.9808 0.9797 0.9786 0.9777 0.9768 0.9761 Table E3. The filter parameters that yield optimal performance at SN R = 5. These are the parameters used in all future experiments with noisy velocity fields (see Fig- ures 3.4 and 3.5). I I I 01 0‘2 72 I 0‘3 73 I 0‘4 7f; Parameter 0.1 0.01 0.53 0.01 0.1 0.0001 0.1 92 Table E.4. E: i 0 (%), in the estimated relative pressure field, for the initial Couette flow pa- rameters, a resolution of 77 = 0.05 on a 20 x 20 flow domain 0F and an embedding level P = 25%, after applying the three filters preposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering and three filters from the literature (PSDF [1], (d/d:r)2 [2], Visc. Diss. [3]) over the range of SNR tested (see Figure 3.6). The average normalized relative RMS error and standard deviation, SNR 5 10 15 No Filter 23.6 i 4.10 11.9 :1: 1.88 8.09 :t 1.48 Median Filter 19.3 i 1.81 12.1 :1: 0.91 10.1 :1: 0.52 Filter 1 14.3 i 2.54 7.65 :t 1.52 5.0834 i 1.06 Filter 2 5.30 :1: 1.96 2.64 :1: 1.28 2.29 a: 0.45 PSDF 22.4 a: 4.38 13.2 i: 1.63 11.11 1.13 d/dr2 10.5 :t 1.82 6.48 i 0.71 4.23 a: 0.60 Visc. Diss. 7.56 i 1.34 5.30 i 0.73 5.14 :t 0.56 SNR 20 25 30 No Filter 5.89 :l: 1.04 4.78 :t 0.81 3.98 :t 0.71 Median Filter 10.1 :1: 0.40 10.2 :1: 0.36 10.2 :1: 0.34 Filter 1 3.63 :l: 0.60 3.03 :5 0.61 2.55 i 0.55 Filter 2 2.94 :1: 0.23 1.95 :t 0.46 2.00 i: 0.26 PSDF 10.32 :1: 0.70 9.85 :1: 0.58 9.58 i 0.46 d/dr2 3.72 :1: 0.52 3.16 i 0.43 2.88 i 0.35 Vise. Diss. 4.78 :t 0.51 4.55 :1: 0.36 4.30 i 0.49 Table E5. The average performance factor F and standard deviation, F :l: 0, for all filters tested using the initial Couette flow experimental parameters, a resolution of 77 = 0.05 on a 20 x 20 flow domain 0F and an embedding level 1" = 25% across the range of SNR tested (see Figure 3.8). SNR 5 10 15 Median Filter 0.95 :l: 0.03 1.60 i 0.05 2.22 :1: 0.22 Filter 1 0.75 :1: 0.04 0.72 a: 0.05 0.74 d: 0.06 Filter 2 0.18 a: 0.03 0.22 :1: 0.04 0.20 i 0.02 PSDF 0.98 i 0.05 1.34 i 0.07 1.82 d: 0.08 d/dr2 0.56 i 0.04 0.74 i 0.05 0.80 d: 0.06 Visc. Diss. 0.46 :t 0.03 0.82 :l: 0.07 1.34 :L- 0.09 SNR 20 25 30 Median Filter 2.99 i 0.14 3.79 i 0.21 4.41 :t 0.34 Filter 1 0.68 :t 0.06 0.76 :t 0.06 0.69 i 0.05 Filter 2 0.19 i 0.03 0.23 1’: 0.04 0.25 :t 0.02 PSDF 2.34 a: 0.15 2.91 i 0.16 3.36 a: 0.19 gar? 0.87 a: 0.02 0.99 i: 0.03 1.12 :t 0.06 Vise. Diss. 1.73 :1: 0.16 2.01 5: 0.22 2.31 :1: 0.16 93 Table E6. The average normalized relative RMS error and standard deviation, E: :l: a (%), in the estimated relative pressure field, for the second Couette flow pa- rameters, a resolution of 77 = 0.05 on a 20 x 20 flow domain {21: and an embedding level P = 25%, after applying the three filters proposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering (see Figure 3.10). SNR 5 10 15 No Filter 37.9 :1: 8.57 18.4 :1: 4.05 12.2 :t 2.66 Median Filter 21.9 i 2.76 13.1 :1: 1.49 9.99 :t 1.12 Filter 1 22.9 :1: 3.48 10.6 d: 1.47 7.01 i 1.24 Filter 2 5.85 :l: 0.60 3.37 :l: 0.97 3.08 :l: 0.36 SNR 20 25 30 No Filter 9.11 :t 1.99 7.29 :t 1.58 6.08 :l: 1.32 Median Filter 8.36 :L- 0.85 7.35 :l: 0.67 6.72 :l: 0.54 Filter 1 5.97 :l: 1.33 4.16 :l: 0.69 4.06 :l: 0.56 Filter 2 2.41 i 0.51 2.05 d: 0.21 2.09 i 0.16 Table 13.7. The average normalized relative RMS error and standard deviation, 5 i 0 (%), in the estimated relative pressure field, for the Poiseuille flow parame— ters, a resolution of n = 0.05 on a 20 x 20 flow domain 0F and an embedding level P = 25%, after applying the three filters proposed here (Median Filter, Filter 1, Filter 2) to the relative pressure field estimated with no filtering and two filters from the literature ((d/d:l:)2 [2], Visc. Diss. [3]) over the range of SNR tested (see Figure 3.11). SNR 5 10 15 20 No Filter 923.4 :1: 78.2 456.1 i 15.8 332.9 :1: 16.3 2399 d: 33.7 Median Filter 968.4 i 51.4 576.7 :1: 27.2 526.4 d: 25.6 488.1 :1: 25.7 Filter 1 572.5 :t 61.1 234.1 a: 73.6 139.9 :t 9.27 145.7 :t 21.4 Filter 2 483.5 :1: 25.8 217.1 :1: 28.1 166.6 i 18.1 113.3 i 6.93 d/dr2 157.5 :1: 14.8 130.5 :t 5.81 120.9 :t 6.94 119.4 d: 19.0 Vise. Diss. 764.7 :t 22.9 742.3 i 85.7 793.2 i 58.6 727.1 i 31.9 94 [ll [2] l3] l4] [5] [6l M [8] l9] BIBLIOGRAPHY S. M. Song, S. Napel, G. H. Glover, and N. J. Pele. Noise reduction in three- dimensional phase contrast MR velocity measurements. J. Magn. Reson. Imag- ing, 3:587—596, 1993. N. Fatouraee and A. A. Amini. Regularization of flow streamlines in multislice phase-contrast MR imaging. IEEE Trans. Med. Imag., 22:699—709, 2003. L. G. Raguin, A. K. Kodali, D. V. Rovas, and J. G. Georgiadis. An inverse problem for reduced-encoding MRI velocimetry in potential flow. Conf. Proc. IEEE Engineering in Medicine and Biology Society, 2:1100—3, 2004. S. M. Song, R. M. Leahy, D. P. Boyd, B. H. Brundage, and S. N apel. Determining cardiac velocity fields and intraventricular pressure distribution from a sequence of ultrafast CT cardiac images. IEEE Trans. Med. Imag., 13:386—397, 1994. P. N. Swarztrauber. Fast Poisson solvers. In G. H. Goulb, editor, Studies in Na- merical Analysis, volume 24, pages 319—70. Mathematical Association of Amer- ica, Washington DC, 1984. J. P. Mohr, L. R. Caplan, J. W. Melski, R. J. Goldstein, G. W. Duncan, J. P Kistler, M. S. Pessin, and H. L. Bleich. The Harvard cooperative stroke registry: A prospecitve registry. Neurology, 28:754—762, 1978. J. Bogousslavsky, G. Van Melle, and F. Regli. The Lausanne stroke registry: Analysis of 1000 consecutive patients with first stroke. Stroke, 19:445—453, 1988. C. W. White, C. B. Wright, D. B. Doty, L. F. Hiratza, C. L. Eastham, D. G. Harrison, and M. L. Marcus. Does visual interpretation of the coronary arteri- ogram predict the physiologic importance of a coronary stenosis? N. Engl. J. Med, 310:819—824, 1984. W. Grossman. Cardiac Catheterization, Angiography, and Intervention, chapter Current Practice Standards, pages 9—10. Willimas & Wilkins, Baltimore, 5th edition, 1996. 95 [10] [11] [12] [131 [14] Ml [16] [17] [18] [191 1201 T. Ebbers, L. Wigstriim, A. F. Bolger, B. Wranne, and M. Karlsson. Noninvasive measurement of time-varying three-dimensional relative pressure fields within the human heart. J. Biomech. Eng-T. ASME, 124:288—293, 2002. A. Nasiraei-Moghaddam, G. Behrens, N. Fatouraee, R. Agarwal, E. T. Choi, and A. A. Amini. Factors affecting the accuracy of pressure measurements in vascular stenosis from phase-contrast MRI. Magn. Resonan. Med, 52:300—309, 2004. G.-Z. Yang, P. J. Kilner, N. B. Wood, S. R. Underwood, and D. N. Firmin. Computation of flow pressure fields from magnetic. resonance velocity mapping. Magn. Resonan. Med, 36:520—526, 1996. J. M. Tyszka, D. H. Laidlaw, J. W. Asa, and J. M. Silverman. Three-dimensional, time—resolved (4D) relative pressure mapping using magnetic resonance imaging. J. Magn. Reson. Imaging, 12:321—329, 2000. D. P. Lum, K. M. Johnson, R. K. Paul, A. S. Turk, D. W. Consigny, J. R. Grinde, C. A. Mistretta, and T. M. Grist. Transstenotic pressure gradients: Measurement in swine — retrospectively ECG-gated 3D phase-contrast MR angiography versus endovascular pressure-sensing guidewires. Radiology, 245:751—760, 2007. L. Hatle, A. Brubakk, A. T‘romsdal, and B. Angelsen. Noninvasive assessment of pressure—drop in mitral-stenosis by Doppler ultrasound. Br. Heart J., 40: 131—140, 1978. S. N. Urchuk and D. B. Plewes. MR measurement of time-dependent blood pressure variations. J. Magn. Reson. Imaging, 51621—627 , 1995. S. N. Urchuk, S. E. Fremes, and D. B. Plewes. In vivo validation of MR pulse pressure measurement in an aortic flow model: Preliminary results. Magn. Res- onan. Med, 38:215—223, 1997. T. Ebbers, L. Wigstrém, A. F. Bolger, J. Engvall, and M. Karlsson. Estimation of relative cardiovascular pressures using time-resolved three-dimensional phase contrast MRI. Magn. Reson. Med, 45:872—879, 2001. R. B. Thompson and E. R. McVeigh. Fast measurement of intracardiac pres— sure differences with 2D breath—hold phase-contrast MRI. Magn. Reson. Med, 49:1056—1066, 2003. J.-P. Tasu, E. Mousseaux, A. Delouche, C. Oddou, O. Jolivet, and J. Bittoun. Estimation of pressure gradients in pulsatile flow from magnetic resonance ac- celeration measurements. Magn. Resonan. Med, 44:66—72, 2000. 96 [21] F. Buyens, O. Jolivet, A. De Cesare, J. Bittoun, A. Herment, J.-P. Tasu, and E. Mousseaux. Calculation of left ventricle relative pressure distribution in MRI using acceleration data. Magn. Resonan. Med, 53:877—884, 2005. [22] N. B. Wood. Aspects of fluid dynamics applied to the larger arteries. J. Theor. Biol, 1992137-161, 1999. [23] F. P. Glor, J. J. M. Westenberg, J. Vierendeels, M. Danilouchkine, and P. Ver- donck. Validation of the coupling of magnetic resonance imaging velocity mea- surements with computational fluid dynamics in a U bend. International Society for Artificial Organs, 26:622-635, 2002. [24] N. B. Wood, S. J. Weston, P. J. Kilner, A. D. Gosman, and D. N. Firmin. Combined MR imaging and CF D simulation of flow in the human descending aorta. J. Magn. Reson. Imaging, 13:699—713, 2001. [25] P. Yim, K. DeMarco, M. A. Castro, and J. Cebral. Plaque Imaging: Pirel to Molecular Level, chapter Characterization of Shear Stress on the Wall of the Carotid Artery Using Magnetic Resonance Imaging and Computational Fluid Dynamics, pages 412—442. IOS Press, 2005. [26] Q. Long, X. Y. Xu, B. Ariff, S. A. Thorn, A. D. Hughes, and A. V. Stanton. Re- construction of blood flow patterns in a human carotid bifurcation: A combined CFD and MRI study. J. Magn. Reson. Imaging, 11:299—311, 2000. [27] J. R. Cebral, P. J. Yim, R. L6hner, O. Soto, and P. L. Choyke. Blood flow modeling in carotid arteries with computational fluid dynamics and MR imaging. Acad. Radial, 9:1286—1299, 2002. [28] P. Papathanasopoulou, S. Zhao, U. Kéhler, M. B. Robertson, Q. Long, P. Hoskins, X. Y. Xu, and 1. Marshall. MRI measurement of time-resolved wall shear stress vectors in a carotid bifurcation model, and comparison with CFD predictions. J. Magn. Reson. Imaging, 17:153—162, 2003. [29] I. Marshall, S. Zhao, P. Papathanasopoulou, P. Hoskins, and X. Y. Xu. MRI and CFD studies of pulsatile flow in healthy and stenosed carotid bifurcation models. Journal of Biomechanics, 37:679-687, 2004. [30] A. Leuprecht, S. Kozerke, P. Boesiger, and K. Perktold. Blood flow in the human ascending aorta: a combined MRI and CF D study. Journal of Engineering Mathematics, 47:387—404, 2003. [31] S. D. Shpilfoygel, R. A. Close, D. J. Valentino, and G. R. Duckwiler. X-ray videodensitometric methods for blood flow and velocity measurement: A critical review of literature. Med. Phys, 27:2008—2023, 2000. 97 [32] [33] [341 [36] [37] [38] [39] [40] [41] [42] S. P. Huang, R. J. Decker, K. G. Goodrich, D. J. Parker, J. B. Muhlestein, D. D. Blatter, and D. L. Parker. Velocity measurement based on bolus tracking with the aid of three-dimensional reconstruction from digital subtraction angiography. Med. Phys, 24:677—686, 1997. N. Hangiandreou, J. Folts, W. Peppler, and C. Mistretta. Coronary blood flow measurement using an angiographic first pass distribution techique: A feasibility study. Med. Phys, 18:947—954, 1991. S. Molloi, Y. J. Qian, and A. Ersahin. Absolute volumetric blood flow measure- ments using dual-energy digital subtraction angiography. Med. Phys, 20:85—91, 1993. S. Molloi, G. Bednarz, J. Tang, Y. Zhou, and T. Mathur. Absolute volumetric coronary blood flow measurement with digital subtraction angiography. Inter- national Journal of Cardiac Imaging, 14:137—145, 1998. S. E. Nissen, J. L. Elion, D. C. Booth, J. Evans, and A. N. DeMaria. Value and limitations of computer analysis of digital subtraction angiography in the assessment of corornary flow reserve. Circulation, 73:562—571, 1986. P. Spiller, F. K. Schmiel, B. Politz, M. Block, U. Fermor, W. Hackbarth, J. Jehle, R. Korfer, and H. Pannek. Measurement of systolic and diastolic flow rates in the coronary artery system by x—ray densitometry. Circulation, 68:337—347, 1983. R. W. Barnes. Noninvasive Diagnostic Techniques in Vascular Disease, chap- ter 4: Continuous-wave Doppler Ultrasound, pages 19—24. The C. V. Mosby Company, 1985. K. W. Beach and D. E. Strandness. Noninvasive Diagnositc Techniques in Vas- cular Disease, chapter 5: Pulsed Doppler Ultrasound for Blood Velocity Mea- surements, pages 25—32. The C. V. Mosby Company, St. Lous, 1985. K. Fusejima. Noninvasive measurement of coronary-artery blood-flow using com- bined two-dimensional and Doppler echocardiography. J. AM. Coll. Cardiol, 10:1024—1031, 1987. H. Mast, J. P. Mohr, J. L. P. Thompson, A. Osipov, S. H. Trocio, S. Mayer, and W. L. Young. Transcranial Doppler ultrasonography in cerbral arteriove- nous malformations: Diagnostic sensitivity and association of flow velocity with spontaneous hemorrhage and focal neurological deficit. Stroke, 26:1024—1027, 1995. T. Hozumi, K. Yoshida, T. Aksaka, Y. Asami, Y. Ogata, T. Takagi, S. Kaji, T. Kawamoto, Y. Ueda, and S. Morioka. Noninvasive assessment of coronary 98 flow velocity and coronary flow velocity reserve in the left anterior descending coronary artery by Doppler echocardiography: Comparison with invasive tech- nique. J. Am. Coll. Cardiol, 32:1251—1259, 1998. [43] T. Hozumi, K. Yoshida, Y. Ogata, T. Akasaka, Y. Asami, T. Takagi, and S. Morioka. Noninvasive assessment of significant left anterior descending coro- nary artery stenosis by coronary flow velocity reserve with transthoracic color Doppler echocardiography. Circulation, 97:1557—1562, 1998. [44] G. N. Hounsfield. Computed medical imaging. Med. Phys, 7:283—290, 1980. [45] L. Axel. Cerebral blood-flow determination by rapid-sequence computed tomog- raphy — A theoretical-analysis. Radiology, 137:679—686, 1980. [46] N. A. Mullani and K. L. Gould. First-pass measurements of regional blood flow with external devices. J. Nucl. Med, 24:577—581, 1983. [47] W. Jaschke, R. G. Gould, P. A. Assimakopoulos, and M. J. Lipton. Flow mea- surements with a high-speed computed tomography scanner. Med. Phys, 14:238— 243, 1987. [48] P. F. Ludman, M. Darby, N. Tomlinson, P. A. Poole-Wilson, and S. Rees. Cardiac flow measurement by ultrafast CT - validation of continuous and pulsatile flow. J. Comput. Assist. Tomogr., 16:795—803, 1992. [49] R. J. Herfkens, L. Axel, M. J. Lipton, S. Napel, W. Berninger, and R. Redington. Measurement of cardiac-output by computed transmission tomography. Invest. Radiol, 17:550—553, 1982. [50] J. S. Garrett, P. Lanzer, W. Jaschke, E. Botvinick, R. Sievers, C. B. Higgins, and M. J. Lipton. Measurement of cardiac output by cine computed tomography. Am. J. Cardiol, 56:657—661, 1985. [01] T. Kaminaga, H. Naito, M. Takamiya, and T. Nishimura. Quantitative eval— uation of mitral regurgitation with ultrafast CT. J. Comput. Assist. Tomogr., 18:239—242, 1994. [52] U. J. Schoepf, R. Bruening, H. Konschitzky, C. R. Becker, A. Knez, J. Weber, O. Muehling, P. Herzog, A. Huber, R. Haberl, and M. F. Reiser. Pulmonary embolism: Comprehensive diagnosis by using electron-beam CT for detection of emboli and assessment of pulmonary blood flow. Radiology, 217:693—700, 2000. [53] P. F. Ludman, A. J. Coats, P. A. Poole—Wilson, and R. S. Rees. Measurement accuracy of cardiac output in humans: indicator dilution techique vs geometric analysis by ultrafast computed tomography. J. Am. Coll. Cardiol, 21:1482—1489, 1993. 99 [54] P. A. Doriot, P. A. Dorsaz, L. Dorsaz, and W. J. Rutishauser. Is the indica- tor dilution theory really the adequate base of many blood flow measurement techniques? Med. Phys, 24:1889—1898, 1997. [55] A. H. Mahnken, M. Grasruck, B. Schmidt, R. W. Gunther, and J. E. Wildberger. Flat panel computed tomography for non-invasive flow measurement: Initial results in In Vitro studies. Eur Radiol, 18:747—752, 2008. [56] S. M. Song and R. M. Leahy. Computation of 3-D velocity fields from 3-D cine CT images of a human heart. IEEE Trans. Med. Imaging, 10:295—306, 1991. [57] F. Bloch, W. W. Hansen, and M. Packard. Nuclear induction. Phys. Rev, 69:127, 1946. [58] E. M. Purcell, H. C. Torrey, and R. V. Pound. Resonance absorption by nuclear magnetic moments in a solid. Phys. Rev, 69:37—38, 1946. [59] P. C. Lauterbur. Image formation by induced local interactions: examples em- ploying nuclear magnetic-resonance. Nature, 242:190—191, 1973. [60] P. Mansfield. Multi-planar image formation using NMR spin echoes. J. Phys. C. Solid State Phys, 102L55—L58, 1977. [61] Z. P. Liang and P. C. Lauterbur. Principles of Magnetic Resonance Imaging: A Signal Processing Perspective. IEEE Press, 2000. [62] C. J. Elkins and M. T. Alley. Magnetic resonance velocimetry: applications of magnetic resonance imaging in the measurement of fluid motion. Exp. Fluids, 43:823—858, 2007. [63] L. Axel and L. Dougherty. MR imaging of motion with spatial modulation of magnetization. Radiology, 171:841—845, 1989. [64] T. J. Mosher and M. B. Smith. A DANTE tagging sequence for the evaluation of translational sample motion. Magn. Reson. Med, 15:334-339, 1990. [65] K. W. Moser, E. C. Kutter, J. G. Georgiadis, R. O. Buckius, H. D. Morris, and J. R. Torczynski. Velocity measurement of flow through a step stenosis using Magnetic Resonance Imaging. Exp. Fluids, 29:438—447, 2000. [66] K. W. Moser, J. G. Georgiadis, and R. O. Buckius. On the use of optical flow methods with spin-tagging magnetic resonance imaging. Annals of Biomedical Engineering, 2929—17, 2001. 100 [67] H. G. Bogren and M. H. Buonocore. 4D magnetic resonance velocity mapping of blood flow patterns in the aorta in young vs. elderly normal subjects. J. Magn. Reson. Imaging, 10:861—869, 1999. [68] M. Markl, F. P. Chan, M. T. Alley, K. L. Wedding, M. T. Draney, C. J. Elkins, D. W. Parker, R. Wicker, C. A. Taylor, R. J. Herfkens, and N. J. Pelc. Time-resolved three-dimensional phase-contrast MRI. J. Magn. Reson. Imaging, 17:499—506, 2003. [69] N. J. Pelc, M. A. Bernstein, A. Sllimakawa, and G. H. Glover. Encoding strategies for three-direction phase-contrast MR imaging of flow. J. Magn. Reson. Imaging, 1:405—413, 1991. [70] C. Tang, D. D. Blatter, and D. L. Parker. Accuracy of phase-contrast flow measurements in the presence of partial-volume effects. J. Magn. Reson. Imagin, 3:377—385, 1993. [71] D. N. Firmin, G. L. Nayler, P. J. Kilner, and D. B. Longmore. The application of phase-shifts in NMR for flow measurement. Magn. Reson. Med, 142230241, 1990. [72] N. J. Pelc, M. T. Alley, J. Listerud, and S. W. Atlas. Magnetic Resonance Imag- ing of the Brain and Spine, chapter Fundamentals of flow and hemodynamics, pages 101—126. Lippincott Willimas & Wilkins, 2002. [73] C. A. Taylor and M. T. Draney. Experimental and computational methods in cardiovascular fluid mechanics. Annu. Rev. Fluid Mech., 36:197—231, 2004. [74] M. F. Walker, S. P. Souza, and C. L. Dumoulin. Quantitative flow measurement in phase-contrast MR angiography. J. Comput. Assist. Tomogr., 12:304—313, 1988. [75] D. Meier, S. Maier, and P. Bésiger. Quantitative flow measurement on phantoms and on blood vessels with MR. Magn. Resonan. Med, 5:25-34, 1988. [76] H. Zhang, S. S. Halliburton, J. R. Moore, 0. P. Simonetti, P. R. Schvartzman, R. D. White, and G. P. Chatzimavroudis. Accurate quantification of steady and pulsatile flow with segmented k-space magnetic resonance velocimetry. Earp. Fluids, 33:458—463, 2002. [77] M. B. Robertson, U. Kéhler, P. R. Hoskins, and 1. Marshall. Quantitative analysis of PC MRI velocity maps: pulsatile flow in cylindrical vessels. Magn. Resonan. Imaging, 19:685—695, 2001. 101 [78] A. E. Ensley, A. Ramuzat, T. M. Healy, G. P. Chatzimavroudis, C. Lucas, S. Sharma, R. Pettigrew, and A. P. Yoganathan. Fluid mechanic assessment of the total cavopulmonary connection using magnetic resonance phase velocity mapping and digital particle image velocimetry. Annals of Biomedical Engineer- ing, 28:1172—1183, 2000. [79] G. A. Varaprasathan, P. A. Araoz, C. B. Higgins, and G. P. Reddy. Quantifi- cation of flow dynamics in congenital heart disease: Applications of velocity- encoded cine MR imaging. Radiographics, 22:895—905, 2002. [80] R. Botnar, G. Rappitsch, M. B. Scheidegger, D. Liepsch, K. Perktold, and P. Boe- siger. Hemodynamics in the carotid artery bifurcation: a comparison between numerical simulations and in vitro MRI measurements. J. Biomech., 33:137—144, 2000. [81] F. M. White. Fluid Mechanics. McGraw — Hill, Inc., 1994. [82] M. C. Potter and J. F. Foss. Fluid Mechanics. Great Lakes Press, Inc, Okemos, MI, 1982. [83] R. Fahraeus and T. Lindqvist. The viscosity of blood in narrow capillary tubes. Am. J. Physiol, 96:562—568, 1931. [84] J. R. Womersley. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J. Physiol, 127:553—563, 1955. [85] J. C. Tannehill, D. A. Anderson, and R. H. Pletcher. Computational Fluid Mehcanics and Heat Transfer. Taylor & Francis, Washington DC, 1984. [86] B. L. Buzbee, G. H. Goulb, and C. W. Nielsen. On direct methods for solving Poisson’s equations. SIAM J. Numer. Anal, 7:627—656, 1970. [87] C-C. J. Kuo and B. C. Levy. Discretization and solution of elliptic PDEs — a digital signal processing approach. Proc. IEEE, 78:1808—1842, 1990. [88] B. L. Buzbee, F. W. Dorr, J. A. George, and G. H. Goulb. The direct solution of the discrete Poisson equation on irregular regions. SIAM J. Numer. Anal, 8:722—736, 1971. [89] A. ’I‘ura, A. Sarti, T. Gaens, and C. Lamberti. Regularization of blood motion fields by modified Navier-Stokes equations. Medical Engineering and Physics, 21:27—36, 1999. [90] L. G. Raguin and J. G. Georgiadis. Kinematics of the stationary helical vortex mode in Taylor—Couette—Poiseuille flow. J. Fluid Mech., 516:125—154, 2004. 102 [91] J. C. Russ. The Image Processing Handbook. CRC Press, Ann Arbor, 1992. [92] B. K. P. Horn and B. G. Schunck. Determining optical-flow. Artificial Intelli- gence, 17:185—203, 1981. 103 NIVERSITY LI I] Ill [1111])“ 3062 8687 mellllllllllllll ll 3 12 9 3 0