Im

Documentation/ECGTool_inverse.tex

@@ -254,6 +254,82 @@ \subsection{Truncated SVD Method (TSVD)}
 % \subsection{Total Variation}
 
 %     \subsubsection{Options and Modes of Operation}
+
+
+
+\subsection{Method of Fundamental Solutions}
+
+Method of Fundamental Solutions (MFS) is a method to approximate
+solutions to partial differential equations, such as Laplace's
+equation used in bioelectricity problems.  MFS is similar to BEM in
+that it is formulated to take into account boundary conditions and
+solutions, but not volumetric ones.  However, MFS is considered
+meshless, even though points representing the heart and torso regions
+are needed, because the points do not need tessellation, as do BEM and
+FEM. In MFS the points representing these regions must be outside the
+computational domain, \ie{} in the myocardium and outside the torso,
+but as close to the boundary as possible.  Also, unlike BEM and FEM
+methods, computing a forward matrix is not necessary.  Instead, the
+MFS computes the coefficients representing pericardial potential
+values based on analytical evaluation of the potential distribution
+and the boundary condition \cite{JDT:Mat77,JDT:Wan2006}.  These
+coefficients are are usually solved for with Tikhonov regularization
+\cite{JDT:Wan2006,JDT:Joh2018}.  In this toolkit, we implemented MFS
+as described by Wang and Rudy \cite{JDT:Wan2006} with the methods to
+choose the regularization parameters described by Johnston
+\cite{JDT:Joh2018}.  It is implemented in a python library
+(`PythonLibrary/MFS\_inverse/mfs\_inverse.py') that can be called from
+SCIRun using the InterfaceWithPython module, as shown in
+`MSF\_inverse\_Cage.srn5' and `MSF\_inverse\_python.srn5'.
+
+\noindent{\bf Inputs:}
+    \begin{enumerate}
+        \item Heart boundary ($H\in\Re^{N,3}$)
+        \item Torso boundary ($T\in\Re^{M,3}$)
+		\item Heart recording points ($P\in\Re^{F,3}$)
+		\item Measured Potentials ($Y\in\Re^{K,T}$)
+		\item Torso Electrodes ($C$ list of $K$ indices)
+        \item options
+    \end{enumerate}
+    {\bf Outputs:}
+     \begin{enumerate}
+        \item Inverse Solution ($X\in\Re^{F,T}$)
+        \item Regularization Parameter ($\lambda$)
+        \item Regularized Curve 
+    \end{enumerate}
+
+
+\subsubsection{Options and Modes of Operation}
+
+Options for the algorithm are the scale method (along the normals or
+by scaling the point cloud), scale factor, regularization parameter
+choosing method \cite{JDT:Joh2018}, lambda range, gamma (for some
+techniques).  The implementation has an option to plot the
+regularization function.
+
+
+
+
+
+
+
+
+\subsection{Isotropy Method}
+
+    The Isotropy Method is an approach that considers temporal correlations
+    to imrove the conditioning of the inverse taken. In SCIRun we have
+    implemented this approach with a combination of modules and
+    \href{http://scirundocwiki.sci.utah.edu/SCIRunDocs/index.php/CIBC:Documentation:SCIRun:Reference:BioPSE:SolveInverseProblemWithTikhonov}{{\tt
+    SolveInverseProblemWithTikhonov}}. Our implementation of the Isotropy
+    Method can be found in the example network
+    ``potential-based-inverse/greensite-inverse-IsotropyMethod.srn5'',
+    shown in \autoref{fig:IsotropyMethod}.
+
+    The network is composed of the standard sections of a simulation
+    network to loading data, visualizing and synthesizing the ECG
+    measurements. The specific blocks that correspond to the Isotropy
+    Method are the ``Temporal Decomposition'', ``solving'' and ``Temporal
+    Reconstruction''. Here we describe this blocks in more detail: 
 
 
 
@@ -262,26 +338,41 @@ \subsection{Truncated SVD Method (TSVD)}
 %     \subsubsection{Options and Modes of Operation}
 
 
-\subsection{Isotropy Method}
 
-    The Isotropy Method is an approach that considers temporal correlations to imrove the conditioning of the inverse taken.
-    In SCIRun we have implemented this approach with a combination of modules and \href{http://scirundocwiki.sci.utah.edu/SCIRunDocs/index.php/CIBC:Documentation:SCIRun:Reference:BioPSE:SolveInverseProblemWithTikhonov}{{\tt SolveInverseProblemWithTikhonov}}.
-    Our implementation of the Isotropy Method can be found in the example network ``potential-based-inverse/greensite-inverse-IsotropyMethod.srn5'', shown in \autoref{fig:IsotropyMethod}.
 
-    The network is composed of the standard sections of a simulation network to loading data, visualizing and synthesizing the ECG measurements.
-    The specific blocks that correspond to the Isotropy Method are the ``Temporal Decomposition'', ``solving'' and ``Temporal Reconstruction''. Here we describe this blocks in more detail:
     \begin{enumerate}
-        \item The {\bf temporal decomposition} block computes the singular value decomposition of the input data, truncates the right singular vectors (corresponding to time) and projects the input data into this lower dimensional space.
-        \item The {\bf solving } block consists on a Tikhonov module that solves the inverse problem for the input data projected onto the low dimensional temporal space.
-        \item The {\bf temporal reconstruction} block, uses the truncated right singular vectors from (1) to reconstruct the full temporal behavior of the inverse solutions obtained in the Tikhonov solver.
+		\item The {\bf temporal decomposition} block computes the
+		singular value decomposition of the input data, truncates the
+		right singular vectors (corresponding to time) and projects
+		the input data into this lower dimensional space.  
+		\item The {\bf solving } block consists on a Tikhonov module
+		that solves the inverse problem for the input data projected
+		onto the low dimensional temporal space.
+		\item The {\bf temporal reconstruction} block, uses the
+		truncated right singular vectors from (1) to reconstruct the
+		full temporal behavior of the inverse solutions obtained in
+		the Tikhonov solver.
     \end{enumerate}
 
     \subsubsection{Options and Modes of Operation}
 
-    As implemented, the Isotropy Method does not have many options for operation. The most important parameter that needs to be determined by the user is the truncation point of the right singular vectors, which can be accessed in the ``SelectSubMatrix'' module within the ``temporal decomposition'' block.
-    The GUI for this module, shown in \autoref{fig:IsotropyMethod_gui}, allows to select a submatrix from an input matrix. In the case of the isotropy method, the users are only interested in the ``end'' parameter from the column range selector (lower right entry) since it determined where to truncate the right singular vectors.
-
-    This implementation of the Isotropy Method has extra parameters that determine the operation of the Tikhonov inverse. These parameters are specific to the \href{http://scirundocwiki.sci.utah.edu/SCIRunDocs/index.php/CIBC:Documentation:SCIRun:Reference:BioPSE:SolveInverseProblemWithTikhonov}{{\tt SolveInverseProblemWithTikhonov}} module and we refer the user to \autoref{sec:inverse:tikhonov} for more details.
+    As implemented, the Isotropy Method does not have many options for
+    operation. The most important parameter that needs to be determined by
+    the user is the truncation point of the right singular vectors, which
+    can be accessed in the ``SelectSubMatrix'' module within the ``temporal
+    decomposition'' block. The GUI for this module, shown in
+    \autoref{fig:IsotropyMethod_gui}, allows to select a submatrix from an
+    input matrix. In the case of the isotropy method, the users are only
+    interested in the ``end'' parameter from the column range selector
+    (lower right entry) since it determined where to truncate the right
+    singular vectors.
+
+    This implementation of the Isotropy Method has extra parameters that
+    determine the operation of the Tikhonov inverse. These parameters are
+    specific to the
+    \href{http://scirundocwiki.sci.utah.edu/SCIRunDocs/index.php/CIBC:Documentation:SCIRun:Reference:BioPSE:SolveInverseProblemWithTikhonov}{{\tt
+    SolveInverseProblemWithTikhonov}} module and we refer the user to
+    \autoref{sec:inverse:tikhonov} for more details.
 
    \begin{figure}
        \begin{center}

Documentation/ECGTool_math.tex

@@ -300,36 +300,34 @@ \subsection{BEM in the Forward/Inverse Toolkit}
 that assumption, the surfaces of those subdomains become a sufficient
 domain upon which to solve the problem for the entire domain.
 
-Briefly, one of the Green's Theorems from vector
-calculus is applied to an integrated form of Laplace's equation to transform the
-differential problem into a Fredholm integral problem \cite{RSM:Bar77}. The surfaces are
-each subdivided (tessellated) into a collection of small surface (or
-boundary) elements. Then (two-dimensional) basis functions (again usually
-low-order polynomials) are used to
-approximate the quantities of interest between the nodes of the resulting
-surface meshes. Given this discretization, after
-manipulation of the resulting integral equation,
-the integrals required can
+
+Briefly, one of the Green's Theorems from vector calculus is applied to an
+integrated form of Laplace's equation to transform the differential problem
+into a Fredholm integral problem \cite{RSM:Bar77}. The surfaces are each
+subdivided (tessellated) into a collection of small surface (or boundary)
+elements. Then (two-dimensional) basis functions (again usually low-order
+polynomials) are used to approximate the quantities of interest between the
+nodes of the resulting surface meshes. Given this discretization, after
+manipulation of the resulting integral equation, the integrals required can
 be computed through a series of numerical integrations over the mesh
-elements.
-In the BEM method, these integrals involve as unknowns the
+elements. In the BEM method, these integrals involve as unknowns the
 potential and its gradient. The integration involves the computation of the
 distance between each node within the surface and to all others surfaces.
 In complicated geometries, and in all cases when the node is integrated
-against the points on its ``own''
-surface, there are numerical difficulties computing these integrals. In those cases
-there are a number of sophisticated solutions which have been proposed in the literature (and
-some of them are adopted in the SCIRun implementation).
-The result of all these integrals is a transfer matrix, which again we will
-denote $\mathbf{A}$, relating the source potentials or currents to the
-unknown measurement potentials. In the BEM case, because of the all-to-all
-nature of the integrations required, this matrix will be dense, not
-sparse. On the other hand, the size of the equation will be directly
-determined by the number of measurements and sources rather than the
-number of nodes in the entire domain. (We note that there is an alternative
-formulation of the BEM method which retains the potentials at all nodes on
-all surfaces, and which can be reduced to the transfer matrix described
-here, but we omit the details as usual.)
+against the points on its ``own'' surface, there are numerical difficulties
+computing these integrals. In those cases there are a number of
+sophisticated solutions which have been proposed in the literature (and
+some of them are adopted in the SCIRun implementation). The result of all
+these integrals is a transfer matrix, which again we will denote
+$\mathbf{A}$, relating the source potentials or currents to the unknown
+measurement potentials. In the BEM case, because of the all-to-all nature
+of the integrations required, this matrix will be dense, not sparse. On the
+other hand, the size of the equation will be directly determined by the
+number of measurements and sources rather than the number of nodes in the
+entire domain. (We note that there is an alternative formulation of the BEM
+method which retains the potentials at all nodes on all surfaces, and which
+can be reduced to the transfer matrix described here, but we omit the
+details as usual.)
 
 
 \section{Solutions to the Inverse Problem in the Forward/Inverse Toolkit}
@@ -341,16 +339,17 @@ \section{Solutions to the Inverse Problem in the Forward/Inverse Toolkit}
 $\mathbf{A}$, calculated by any appropriate method, including either FEM or
 BEM.
 
-In the activation-based case, the source model is that the unknowns are an
-activation surface. (That is, activation times as a function of position on
-the heart surface, which we denote as $\tau(x)$, where $x$ indicates
-position on the heart surface.)
-The assumptions required for the activation-based model imply that the temporal waveform of the potential
-(and current) at each source node has a fixed form, the same at all
-locations on the surface.
-This is assumed to be either a step function or a smoothed version of a step function (using
-piecewise polynomials or inverse trigonometric functions). We denote this
-function as $u(t)$. Thus the relevant forward equation can be written as
+In the activation-based case, the source model is that the unknowns
+are an activation surface. (That is, activation times as a function of
+position on the heart surface, which we denote as $\tau(x)$, where $x$
+indicates position on the heart surface.) The assumptions required for
+the activation-based model imply that the temporal waveform of the
+potential (and current) at each source node has a fixed form, the same
+at all locations on the surface. This is assumed to be either a step
+function or a smoothed version of a step function (using piecewise
+polynomials or inverse trigonometric functions). We denote this
+function as $u(t)$. Thus the relevant forward equation can be written
+as
 
 \begin{equation} y(p,t) = \int_{x} \mathbf{A}_{p,x}u(t-\tau(x))\,dx \label{eq:act}
 \end{equation}
@@ -430,9 +429,14 @@ \subsubsection{Standard Tikhonov regularization}
 \end{eqnarray}
 \end{center}
 
-\noindent where $\lambda$ is the regularization parameter, which is a user defined scalar value. The matrix $\mathbf{P}$ represents the \textit{a priori} knowledge of the measurements. The matrix $\mathbf{L}$ describes the property of the solution $x$ to be constrained.
-Conceptually, $\lambda$ trades off between the misfit between predicted and measured data (the first term in the equation) and the \textit{a priori} constraint.
-An approximate solution $\hat{x}$ of \autoref{tik_problem} is given for the
+
+\noindent where $\lambda$ is the regularization parameter, which is a user
+defined scalar value. The matrix $\mathbf{P}$ represents the \textit{a
+priori} knowledge of the measurements. The matrix $\mathbf{L}$ describes
+the property of the solution $x$ to be constrained. Conceptually, $\lambda$
+trades off between the misfit between predicted and measured data (the
+first term in the equation) and the \textit{a priori} constraint. An
+approximate solution $\hat{x}$ of \autoref{tik_problem} is given for the
 overdetermined case ($n > m$) as follows:
 \begin{center}
 \begin{eqnarray}

Documentation/ECGTool_overview.tex

@@ -135,7 +135,7 @@ \section{Toolkit overview and capabilities}
 \caption{Algorithms currently included in the CIBC ECG Forward/Inverse
 toolkit \newline
 $^\dag$Activation-based forward solution is not yet
-unavailable\newline
+available\newline
 $^{\&}$Based on a Gauss-Newton optimization \newline
 $^{*}$ Matlab implementation
 $^\ddag$ Python implementation

Documentation/toolkit.bib

@@ -10,6 +10,20 @@ @INPROCEEDINGS{JDT:Gul97
     month="Oct"
 }
 
+@article{JDT:Joh2018,
+title = "Accuracy of electrocardiographic imaging using the method of fundamental solutions",
+journal = "Computers in Biology and Medicine",
+volume = "102",
+pages = "433 - 448",
+year = "2018",
+issn = "0010-4825",
+doi = "https://doi.org/10.1016/j.compbiomed.2018.09.016",
+url = "http://www.sciencedirect.com/science/article/pii/S0010482518302786",
+author = "Peter R. Johnston",
+keywords = "Inverse problems, Electrocardiology, Tikhonov regularisation, Method of fundamental solutions, Electrophysiology",
+abstract = "Solving the inverse problem of electrocardiology via the Method of Fundamental Solutions has been proposed previously. The advantage of this approach is that it is a meshless method, so it is far easier to implement numerically than many other approaches. However, determining the heart surface potential distribution is still an ill-posed problem and thus requires some form of Tikhonov regularisation to obtain the required distributions. In this study, several methods for determining an “optimal” regularisation parameter are compared in the context of solving the inverse problem of electrocardiology via the Method of Fundamental Solutions. It is found that the Robust Generalised Cross-Validation method most often yields epicardial potential distributions with the least relative error when compared to the input distribution. The study also compares the inverse solutions obtained with the Method of Fundamental Solutions with those obtained in a previous study using the boundary element method. It is found that choosing the best solution methodology and regularisation parameter determination method depends on the particular scenario being considered."
+}
+
 @Article{RSM:Gul98,
   author =       "R.M. Gulrajani",
   title =        "The forward and inverse problems of
@@ -41,6 +55,19 @@ @Article{RSM:Bar77
                  development of the problem., ECG055",
 }
 
+@article{JDT:Mat77,
+author = "Mathon, R. and Johnston, R.",
+title = "The Approximate Solution of Elliptic Boundary-Value Problems by 
+Fundamental Solutions",
+journal = j-SIAM-J-NUM,
+volume = "14",
+number = "4",
+pages = "638-650",
+year = "1977",
+doi = "10.1137/0714043",
+URL = "https://doi.org/10.1137/0714043",
+}
+
 
 @Article{RSM:Oos2004,
   author =       "A. van Oosterom and T.F. Oostendorp",
@@ -73,3 +100,22 @@ @Article{RSM:Oos2004
                  www.ecgsim.org.",
   bibdate =      "Wed Aug 11 06:37:42 2004",
 }
+
+@article{JDT:Wan2006,
+	Author = "Wang, Yong and Rudy, Yoram",
+	Doi = "10.1007/s10439-006-9131-7",
+	Journal = j-ABE,
+	Keywords = "Action Potentials/*physiology; *Algorithms; Animals; Body 
+    Surface Potential Mapping/*methods; Computer Simulation; Diagnosis,
+    Computer-Assisted/*methods; Electrocardiography/methods; Electromagnetic
+    Fields; Heart Conduction System/*physiology; Humans; *Models,
+    Cardiovascular",
+	Month = "08",
+	Number = "8",
+	Pages = "1272--1288",
+	Title = "Application of the method of fundamental solutions to 
+    potential-based inverse electrocardiography",
+	Url = "https://www.ncbi.nlm.nih.gov/pubmed/16807788",
+	Volume = "34",
+	Year = "2006",
+}