|
| 1 | +\chapter{Inverse Solutions} |
| 2 | +\label{ch:inv} |
| 3 | + |
| 4 | +\section{Overview} |
| 5 | + |
| 6 | +The inverse problem of electrocardiography is to find suitable electrical source parameters on the heart that adequately describe the observed body surface potentials. Regularization is typically employed in solution methods to reduce the sensitivity of the problem to relatively small errors in the observed body surface potentials, thereby stabilizing it. As a result, solution methods may be supplied body surface potentials, a forward model, and method-specific regularization parameters as input. Specific use of each method is described below. As output, modules primarily produce the resulting solution with extra outputs depending on the specific method. |
| 7 | + |
| 8 | +\section{Module Descriptions for Inverse Solution Methods} |
| 9 | + |
| 10 | +%%%%%%%%%%%%%%%%%%%%%%% |
| 11 | +\subsection{Tikhonov Regularization} |
| 12 | + |
| 13 | +\begin{figure}[H] |
| 14 | +\begin{center} |
| 15 | +\includegraphics[width=\textwidth]{ECGToolkitGuide_figures/SolveInverseProblemWithTikhonov.png} |
| 16 | +\caption{Tikhonov Regularization Module} |
| 17 | +\label{tikhonov} |
| 18 | +\end{center} |
| 19 | +\end{figure} |
| 20 | + |
| 21 | +\vspace{5pt}\textit{This is a method of solving the potential-based inverse problem.}\vspace{5pt} |
| 22 | + |
| 23 | +Tikhonov regularization is a common method of stabilizing solutions to many inverse problems including inverse ECG. The module that implements Tikhonov regularization (closed-form solution) in SCIRun is called ``SolveInverseProblemWithTikhonov''. The module requires that a forward solution matrix, regularization matrix (such as a surface Laplacian approximation), and a vector (single-column matrix) of observed/measured body surface potentials are supplied as input. |
| 24 | + |
| 25 | +In the user interface (UI), one can explicitly specify a scalar regularization parameter that weights the influence of the regularization matrix on the solution. The module can also automatically select a regularization parameter using the L-curve method. |
| 26 | + |
| 27 | +Upon completion, the module will output the inverse solution. In addition, it will either output the specified or automatically-chosen regularization parameter as well as the regularized inverse matrix (a closed-form Tikhonov inverse operator). |
| 28 | + |
| 29 | + |
| 30 | +%%%%%%%%%%%%%%%%%%%%%%% |
| 31 | +\subsection{Tikhonov Singular Value Decomposition (SVD)} |
| 32 | + |
| 33 | +\begin{figure}[H] |
| 34 | +\begin{center} |
| 35 | +\includegraphics[width=\textwidth]{ECGToolkitGuide_figures/SolveInverseProblemWithTikhonovSVD.png} |
| 36 | +\caption{Tikhonov SVD Module} |
| 37 | +\label{tikhonovsvd} |
| 38 | +\end{center} |
| 39 | +\end{figure} |
| 40 | + |
| 41 | +\vspace{5pt}\textit{This is a method of solving the potential-based inverse problem.}\vspace{5pt} |
| 42 | + |
| 43 | +This module implements Tikhonov regularization (closed-form solution) using the singular value decomposition (SVD) in SCIRun and is called ``SolveInverseProblemWithTikhonovSVD''. The module requires that the SVD of a forward solution matrix (left/right singular vector matrices and singular value matrix), regularization matrix (such as a surface Laplacian approximation), and a vector (single-column matrix) of observed/measured body surface potentials are supplied as input. |
| 44 | + |
| 45 | +In the user interface (UI), one can explicitly specify a scalar regularization parameter that weights the influence of the regularization matrix on the solution. The module can also automatically select a regularization parameter using the L-curve method. |
| 46 | + |
| 47 | +Upon completion, the module will output the inverse solution. In addition, it will either output the specified or automatically-chosen regularization parameter as well as the regularized inverse matrix (a closed-form Tikhonov inverse operator computed using the provided SVD). |
| 48 | + |
| 49 | +%%%%%%%%%%%%%%%%%%%%%%% |
| 50 | +\subsection{Truncated SVD} |
| 51 | + |
| 52 | +\begin{figure}[H] |
| 53 | +\begin{center} |
| 54 | +\includegraphics[width=\textwidth]{ECGToolkitGuide_figures/SolveInverseProblemWithTSVD.png} |
| 55 | +\caption{Truncated SVD Module} |
| 56 | +\label{tsvd} |
| 57 | +\end{center} |
| 58 | +\end{figure} |
| 59 | + |
| 60 | +\vspace{5pt}\textit{This is a method of solving the potential-based inverse problem.}\vspace{5pt} |
| 61 | + |
| 62 | +This module, called ``SolveInverseProblemWithTSV'', computes an inverse solution by creating a pseudo-inverse operator and multiplying that with the measured body surface potentials. The module requires that the SVD of a forward solution matrix (left/right singular vector matrices and singular value matrix) and a vector of observed/measured body surface potentials are supplied as input. |
| 63 | + |
| 64 | +In the user interface (UI), one can explicitly specify a scalar integer regularization parameter to control the amount of regularization imposed on the problem. Assuming the singular values are sorted from largest to smallest (and therefore most significant to least significant), the regularization parameter controls how many significant singular values to consider (truncation degree) when calculating the pseudo-inverse. Alternatively, the module will use the L-curve method to automatically obtain a suitable truncation degree. |
| 65 | + |
| 66 | +The primary output of this module is a vector containing the inverse solution. In addition, it will output the regularization parameter and regularized inverse operator used to obtain the inverse solution from the measured body surface potentials. |
| 67 | + |
| 68 | +%%%%%%%%%%%%%%%%%%%%%%%%% |
| 69 | + |
| 70 | +\subsection{Matlab Interface} |
| 71 | + |
| 72 | +The \textit{InterfaceWithMatlab} module allows a SCIRun network to make use of Matlab for some of its calculations. In Figure \ref{matlabinterfacemodule}, the user interface for this module shows how SCIRun data types from the module's input pipes are converted to Matlab variables for use in a Matlab script. Upon completion of the Matlab script, Matlab variables are converted back to SCIRun data types to be sent to the module's output pipes. |
| 73 | + |
| 74 | +\begin{figure}[H] |
| 75 | +\begin{center} |
| 76 | +\includegraphics[width=0.8\textwidth]{ECGToolkitGuide_figures/InterfaceWithMatlab.png} |
| 77 | +\caption{The InterfaceWithMatlab Module and its corresponding user interface} |
| 78 | +\label{matlabinterfacemodule} |
| 79 | +\end{center} |
| 80 | +\end{figure} |
| 81 | + |
| 82 | +In the rest of this section, we will describe some of the inverse solution methods distributed with SCIRun that have been implemented in Matlab. |
| 83 | + |
| 84 | + |
| 85 | +\subsubsection{Isotropy Method} |
| 86 | + |
| 87 | +\vspace{5pt}\textit{This is a method of solving the potential-based inverse problem.}\vspace{5pt} |
| 88 | + |
| 89 | +\begin{verbatim} |
| 90 | +function X_reg=greensite(A,Y,trunc_deg) |
| 91 | +% Function to calculate the Greensite inverse solution |
| 92 | +% A - forward matrix |
| 93 | +% Y - data |
| 94 | +% trunc_deg - truncation degree |
| 95 | +% X_reg - inverse solution |
| 96 | +\end{verbatim} |
| 97 | + |
| 98 | +The isotropy method (a.k.a. Greensite method) takes as input a forward solution matrix, measured body surface potential data, and an optional integer regularization parameter (truncation degree). In this method, the truncation degree controls the number of significant singular values of the measured body surface potential data to be used when computing the inverse solution. Consequently, the truncation degree partitions the measured body surface potential data into independent ``signal'' and ``noise'' components, ignoring the ``noise'' component in its calculation of the inverse solution. The only output of the method is the inverse solution. |
| 99 | + |
| 100 | +\subsubsection{Gauss-Newton Method} |
| 101 | + |
| 102 | +\vspace{5pt}\textit{This is a method of solving the activation-based inverse problem.}\vspace{5pt} |
| 103 | + |
| 104 | +\begin{verbatim} |
| 105 | +function tau = ActGaussNewton(A,Y,L,tauinit,lambda,w,minstep) |
| 106 | +% Implements the Gauss-Newton algorithm for solving the activation-based |
| 107 | +% inverse problem of electrocardiography. |
| 108 | +% => minimizes the objective function ||Y-A*X||^2+lambda*||L*X||^2 where |
| 109 | +% X is parameterized by the C^1 polynomial approximation to a step function |
| 110 | +% as explained in "The Depolarization Sequence of the Human Heart Surface |
| 111 | +% Computed from Measured Body Surface Potentials" by Geertjan Huiskamp and |
| 112 | +% Adriaan van Oosterom. |
| 113 | +% |
| 114 | +% Input Variables: |
| 115 | +% A: Forward matrix |
| 116 | +% Y: Observations (columns index time from 1 to T=size(Y,2)) |
| 117 | +% L: Regularization matrix (typically a surface Laplacian approximation) |
| 118 | +% lambda: Regularization parameter |
| 119 | +% w: Width parameter in step function approximation |
| 120 | +% tauinit: Initial phase shifts for starting the algorithm |
| 121 | +% |
| 122 | +% Output Variables: |
| 123 | +% tau: Solution phase shifts of the step functions |
| 124 | +\end{verbatim} |
| 125 | + |
| 126 | +The activation-based inverse problem of electrocardiography is to solve for electrical activation times (i.e. depolarization times) on the heart given body surface potentials during the QRS complex. This results in a nonlinear least squares optimization problem that this method solves using the Gauss-Newton algorithm. Solutions are typically highly dependent on choice of initial guess. This method requires a forward solution matrix, body surface potentials, regularization matrix, initial guess, regularization parameter, activation waveform transition width, and convergence parameter to be specified as input. The body surface potentials should be provided as a $M \times T$ matrix, where $M$ is the number of leads from which potentials are recorded and $T$ is the number of time samples recorded during the QRS complex of the heart beat whose activation times are to be estimated. The regularization matrix and corresponding parameter influence the inverse solution as in Tikhonov regularization (see above). The initial guess is the set of activation times from which the Gauss-Newton algorithm starts pursuing a more suitable inverse solution. The transition width controls the number of time samples (not necessarily integer) taken by each source to transition from inactivated to activated in this method. The convergence parameter is the minimum norm of the step taken by the iterations within the Gauss-Newton algorithm before the method is deemed to have converged to a suitable solution. The only output of the method is the inverse solution (in this case: an array/vector of activation times). |
| 127 | + |
| 128 | +\subsubsection{Wavefront-Based Potential Reconstruction} |
| 129 | + |
| 130 | +\vspace{5pt}\textit{This is a method of solving the potential-based inverse problem.}\vspace{5pt} |
| 131 | + |
| 132 | +\begin{verbatim} |
| 133 | +function [x_WBPR_forward,x_WBPR_backward] = WBPR(A,y,heart,first_act,last_act) |
| 134 | +% Function to calculate the WBPR solutions |
| 135 | +% Inputs: |
| 136 | +% A: forward matrix |
| 137 | +% y: torso data |
| 138 | +% heart: heart geometry |
| 139 | +% first_act: first activated node |
| 140 | +% last_act: last activated node |
| 141 | +% |
| 142 | +% Outputs: |
| 143 | +% x_WBPR_forward: inverse solution using forward WBPR |
| 144 | +% x_WBPR_backward: inverse solution using backward WBPR |
| 145 | +\end{verbatim} |
| 146 | + |
| 147 | +The Wavefront-Based Potential Reconstruction (WBPR) method imposes prior knowledge about the spatial patterns of electric potentials on the heart during the QRS complex to reconstruct an inverse solution. This method requires a forward solution matrix, measured body surface potentials, a structure containing the triangles and nodes of the heart geometry mesh, and the first/last nodes to activate during the QRS complex of the heart beat. The body surface potentials should be provided as a $M \times T$ matrix, where $M$ is the number of leads from which potentials are recorded and $T$ is the number of time samples recorded during the QRS complex of the heart beat whose electric potentials are to be estimated. The heart geometry is a Matlab ``struct'' that must contain a matrix of node coordinates and a matrix of indices for triangles that connect the nodes on the surface of the heart. The first/last nodes to activate are specified as the indices of the nodes. |
| 148 | + |
| 149 | +As output, WBPR produces two inverse solutions. The procedure by which the inverse solutions are obtained may be run both forward and backward in time. Therefore, the ``forward WBPR'' solution is that obtained by applying the method starting from the earliest sample time and ending at the latest. On the other hand, the ``backward WBPR'' solution is obtained by starting from the latest sample time and traversing the samples in reverse until ending at the earliest. |
| 150 | + |
| 151 | +%%%%%%%%%%%%%%%%%%%%%%%%% |
| 152 | + |
| 153 | +\section{Example Simulations} |
| 154 | + |
| 155 | +\subsection{Tikhonov Regularization} |
| 156 | + |
| 157 | +\vspace{5pt}\textit{The SCIRun network for this example can be found at:\\``\{ECGTOOLKITDIRECTORY\}/potential-based-inverse/tikhonov-inversion.srn''\\in the SCIRun source code directory.}\vspace{5pt} |
| 158 | + |
| 159 | +\begin{figure}[H] |
| 160 | +\begin{center} |
| 161 | +\includegraphics[width=0.9\textwidth]{ECGToolkitGuide_figures/TikhonovNetwork.png} |
| 162 | +\caption{The SCIRun network for the Tikhonov inverse solution example.} |
| 163 | +\label{TikhonovNetworkExample} |
| 164 | +\end{center} |
| 165 | +\end{figure} |
| 166 | + |
| 167 | +This example shows how to use the ``SolveInverseProblemWithTikhonov'' module in a SCIRun network to solve a potential-based ECG inverse problem using Tikhonov regularization. Specifically, this example simulates ``measured'' body surface potentials by applying a forward solution and adding pseudorandom noise with a specified signal-to-noise ratio (SNR). The forward solution is obtained by multiplying a known vector of epicardial potential data by a pre-computed forward solution matrix. The simulated measurements and forward solution matrix are the inputs to the ``SolveInverseProblemWithTikhonov'' module. The regularization parameter of this method is controlled via the user interface (UI) in this example. The output from the module is the inverse solution. This example compares the inverse solution to the known epicardial potentials from which the measurements were simulated. A visualization of both can be seen in the UI of the ``ViewScene'' module. |
| 168 | + |
| 169 | +\subsection{Gauss-Newton Method} |
| 170 | + |
| 171 | +\vspace{5pt}\textit{The SCIRun network for this example can be found at:\\``\{ECGTOOLKITDIRECTORY\}/activation-based-inverse/actgaussnewton-inversion.srn''\\in the SCIRun source code directory.}\vspace{5pt} |
| 172 | + |
| 173 | +\begin{figure}[H] |
| 174 | +\begin{center} |
| 175 | +\includegraphics[width=0.9\textwidth]{ECGToolkitGuide_figures/actgaussnewtonnetwork.png} |
| 176 | +\caption{The SCIRun network for the activation-based Gauss-Newton inverse solution example.} |
| 177 | +\label{GaussNewtonNetworkExample} |
| 178 | +\end{center} |
| 179 | +\end{figure} |
| 180 | + |
| 181 | +This example shows how to use the ``InterfaceWithMatlab'' module in a SCIRun network to solve an activation-based ECG inverse problem using the Gauss-Newton method. This example uses the ``normal male'' dataset and geometries from ECGSIM. Specifically, this method compares measured body surface potentials with body surface potentials simulated from activation times. We search the set of possible activation times for ones that minimize the squared error residual between the measured and simulated body surface potentials. The Gauss-Newton method in this example is implemented as a Matlab function and is used in SCIRun via the ``InterfaceWithMatlab'' module. As a numerical optimization method, Gauss-Newton will look for a local minimizer of the squared error residual starting from a suitable initialization point. In this example, we initialize the method with the activation times reported for the dataset in ECGSIM, but perturbed by noise. |
| 182 | + |
0 commit comments