to provide a comprehensive understanding of the invention , its specific illustrative embodiments are described below ; however , those of ordinary skill in the art will recognize that methods and systems may be modified within the scope of the invention as defined by the appended claims . methods and systems disclosed here use a device of surface ecg mapping , visualization techniques of computer ( ct ) or magneto - resonance ( mrt ) tomography , computing techniques , as well as mathematical algorithms of solution of the inverse problem of electrocardiography for non - invasive reconstructing electrograms at internal points of the chest and on the heart epicardial surface and for constructing isopotential and isochronous epicardial maps on a realistic three - dimensional ( 3 - d ) computer model of the heart . fig1 illustrates a general schematic view of the method . the method includes ( 1 ) a registration of 240 unipolar ecg on the chest surface , ( 2 ) an implementation of ct or mrt of the chest , ( 3 ) data processing of surface ecg mapping and of computer ( mrt ) tomography using computing techniques and ( 4 ) a representation of the obtained electrophysiological information with means of computer graphics . fig2 illustrates a schematic view of the methodology of surface ecg mapping . a mapping device comprises a digital multi - channel electrocardiograph ( 1 ) connected with a personal computer ( 2 ). the digital multi - channel electrocardiograph allows one to register ecg - signals in 12 standard leads and in up to 240 unipolar leads from the chest surface . fig3 illustrates a scheme of imposing electrodes . for surface ecg mapping , one - off chlorine - silver electrodes are used ( 1 ). electrodes are imposed in the form of 5 - 8 horizontal strips ( belts ) positioned at similar distances along the vertical . the upper strip is positioned at the level of sterno - cleidal articulation , the lower strip — at the level of lower edge of rib - arch . each strip includes from 16 to 30 electrodes placed at similar distances in circumference of the chest ( 2 ). when a roentgen computer tomography is used as a visualization methodology , one - off metal chlorine - silver electrodes are applied as they are well visualized in roentgen tomography images and give a minimum level of artifacts . when a magneto - resonance therapy is used as a visualization methodology , one - off graphite electrodes are applied as they show the similar properties for this tomography technique . fig4 depicts the main stages of computer processing of the information . the stage ( 1 ) is a real - time processing of ecg - signals in the course of multi - channel ecg registration from the chest surface . the stage ( 2 ) is a retrospective processing of ecg - signals . the stage ( 3 ) includes constructing voxel models of the chest , heart and its compartments on ct or mrt data . the stage ( 4 ) comprises constructing polygonal surfaces of the chest , heart and its compartments . the stage ( 5 ) includes an automatic determination of coordinates of registration electrodes on the chest surface according to ct or mrt data . at stage ( 6 ) a surface interpolation of values of surface mapping ecg - signals at each discrete moment and a construction of isopotential maps on the chest surface are performed . the stage ( 7 ) includes a computational reconstruction of the heart electric field potential at internal points of the chest and on the heart epicardial surface . at the last stage , reconstructing epicardial electrograms ( 8 ) and constructing epicardial isopotential , isochronous maps with using means of computer graphics ( 9 ) on a realistic computer model of the heart , and visualizing the dynamics of electrophysiological processes of the myocardium in animation mode ( propagation mapping ) ( 10 ) are performed , respectively . fig5 illustrates processing of ecg - signals in the course of real - time ecg mapping . ecg - signals registered are reflected in computer display . an operator controls the quality of an ecg - signal in each of the leads ; if necessary , a programmed suppression of power - line ( 1 ) muscle ( 2 ) noises and of isoline drift ( 3 ) is used . automatic control of the contact of an electrode with skin and correctness of imposing electrodes based on spectral and mutual - correlation analyses of ecg - signals are also performed . a result of stage ( 1 ) represents digitalized and filtered values of ecg - signals in 240 unipolar leads from the chest surface and in 12 standard leads with duration up to 3 minutes . an operator looks through ecg - signals registered and selects one or several cardiocycles ( 1 , 2 ) for posterior processing . further , a reduction of ecg to a unity isoline ( 3 , 4 ) is implemented : to this end , operator selects in one of ecgs such a time interval τ within which an ecg - signal coincides with an isoline ( as a rule , this interval belongs to the segment pq ). correction of ecg - signals is performed according to the formula : where u 0 ( t ) is the corrigiert ecg - signal , u ( t ) is an initial ecg - signal , u 0 is an average value of an initial ecg - signal at time interval τ . afterwards , operator selects a cardiocycle fragment , being under interest , for subsequent calculations . fig7 illustrates constructing a voxel model of the torso and heart in voxel graphics editor . on ct or mrt data of the chest and heart , a voxel rendering of anatomical structures of the chest is carried out . to this end , a shear - warp factorization of the viewing transformation algorithm is used , which belongs to a group of scanline - order volume rendering methods . the concept of the voxel rendering method applied here consists of three main steps ( philippe lacroute fast volume rendering using a shear - warp factorization of the viewing transformation .— ph . d . dissertation , technical report csl - tr - 95 - 678 , stanford university , 1995 ). at first step , volume data are transformed by a shear matrix in the corresponding object space , each parallel slice of volume data after transformations passing through a special filter for diminishing distortions . at second step , an intermediate 2d image within the same shear space is formed from a combined set of filtered and sheared slices by their direct - order superposition . at third step , the intermediate 2d image obtained is transferred in a normal image space with using a shear matrix , and further it again passes through a filter for formation of a final image . an operator with the help of instruments of voxel edition makes ready a voxel model of the torso , heart or one of its structures . fig8 illustrates constructing polygonal surfaces ( triangulation grids ) of the torso and heart on the basis on voxel models . based on the obtained voxel models , polygonal surfaces consisting of united planar triangles and volume tetrahedral finite - element meshes are automatically constructed . initial data represent a three - dimensional scalar field of densities in a voxel representation , i . e ., a three - dimensional right - angled grid , in whose nodes values of conditional densities of chest tissues are given . constructing triangulation grids of the torso and heart represents a construction of polygonal surfaces that by the nearest way repeat surfaces of aforesaid structures given by the definite level of density . filtrating initial voxel models for diminishing a casual noise level ; constructing a finite - element volume and surface grid on the basis of the exhaustion method , more known in english - written literature as an advancing front algorithm . advancing front algorithm is described in more detail in lo s . h . volume discretization into tetrahedra — ii . 3d triangulation by advancing front approach // computers and structures , pergamon , vol . 39 , no 5 , p . p . 501 - 511 , 1991 ; rassineux a . generation and optimization of tetrahedral meshes by advancing front technique // international journal for numerical methods in engineering , wiley , vol . 41 , p . p . 651 - 674 , 1998 ; gol &# 39 ; nik e r ., vdovichenko a . a ., uspekhov a . a . construction and application of a preprocessor of generation , quality control , and optimization of triangulation grids of contact systems // information technologies , 2004 , no . 4 , p . 2 - 10 [ in russian ]. at the next step , a specific electroconductivity coefficient of a biological tissue is determined for each node of a finite - element grid . firstly , a type of a biological tissue is determined based on hounsfield numbers in computer tomograms or values of a mr - signal in magneto - resonance tomograms . afterwards , a specific electroconductivity coefficient is ascribed to every type of a biological tissue on the basis of published data . an example of conformity between hounsfield numbers and specific electroconductivity values of chest tissues is given below ( hofer m . computer tomography teaching manual [ russian translation ]. moscow : meditsinskaya literatura , 2006 ; martirosov e . g ., nikolaev d . v ., rudnev s . g . technologies and methods for determination of human body composition . [ in russian ].— moscow : nauka , 2006 ). fig1 illustrates constructing isopotential maps on the torso surface . constructing isopotential maps is performed by surface interpolation of ecg - signal values at each discrete moment with using radial basis functions . the electric field potential on the chest surface , s , is represented in the form of a decomposition according to the system of radial basis functions ( rbf ): where u ( x ) is the electric field potential , f i ( x ) are radial basis functions , a i are indefinite coefficients . as rbf , functions of the following kind given at ecg - registration points are used f j ⁡ ( x ) = exp ⁡ ( -  x - x j  c 2 ) , where x is a random point on the body surface , x j are ecg - registration points , ∥ x − x j ∥ is the minimal length of a line belonging to the surface s and connecting points x and x j , c is an experimentally chosen coefficient that determines approximation properties of a function . coefficients a j are found from the condition for the minimum of functional j : j = 1 2 ⁢ ∑ i = 1 n ⁢ [ ( ∑ j = 1 n ⁢ a j ⁢ f j ⁡ ( x i ) + a 0 ) - u ⁡ ( x i ) ] 2 where u ( x i ) are values of the electric field potential at x i points of ecg - registration on the chest surface , n is a number of ecg - registration points . for finding coefficients a j , the corresponding system of linear algebraic equations with a matrix of n × n size is solved . the potential u ( x i ) is calculated in nodes of the torso triangulation surface x i according to the formula to calculate the potential at each point of the torso surface , a bilinear interpolation on values in vertices of a grid triangle , which this point belongs to , is applied . the claimed method includes a method for noninvasive reconstructing the heart electric field potential at internal points of the chest on measured values of the electric field potential on the chest surface by a numerical solution of the inverse problem of electrocardiography for a model of the chest with a variable electroconductivity coefficient with using the finite element method based on iteration algorithms . for realizing this method , the following model is used . let ωεr 3 be a part of the chest bounded by a sufficiently smooth border ∂ ω , which includes the torso surface being in contact with external medium γ b , cross - sections of the chest at the level of the diaphragm and clavicles γ t1 γ t2 , as well as the heart epicardial surface γ e . chest tissues in domain ω are assumed to have a variable continuous positive limited specific electroconductivity coefficient k ( x ), xεω ∪∂ ω . the heart electric field potential in domain ω is assumed to satisfy the laplace equation in an inhomogeneous medium : where x =( x 1 , x 2 , x 3 ) t εω ⊂ r 3 is a point in three - dimensional ( 3d ) space , ∇ ≡ ( ∂ ∂ x 1 , ∂ ∂ x 2 , ∂ ∂ x 3 ) at the part of border γ b of domain ω , the dirichlet condition ( electric field potential measured as a result of surface ecg mapping ) is assumed to be known u ( x )= u ( x ), xεγ b , uεl 2 ( γ b ) ( 2 ) the dirichlet condition contains a noise component as the result of experimental measurements : u ( x )= u 0 ( x )+ ξ ( x ), xεγ b , u 0 εc 28 , ξεγ 2 ( γ b ), ( 3 ) where u 0 ( x ) is the exact value of the potential on the chest surface , ξ ( x ) is an measurement error estimated as ∥ ξ ( x )∥ l 2 & gt ; δ . ∂ u ⁡ ( x ) ∂ n = p ⁡ ( x ) = 0 , x ∈ γ b , p ∈ l 2 ⁡ ( γ b ) , ( 4 ) is a potential derivative u ( x ) along an internal normal to the surface . solution of the inverse problem of electrocardiography consists in finding in the class of functions l 2 ( γ h ) a potential trace u ( y ) on the surface γ h that satisfies the laplace equation in domain ω ( 4 ) and the boundary conditions ( 5 )-( 7 ) at borders of regions . 1 . it is required to find a potential u ( x ) such as that : let us name this problem as a direct one in respect of the inverse boundary problem under study . 2 . it is required to find a potential u ( x ) such as that : let us name this problem as a conjugate problem in respect of the direct problem . let u ( y ) be a trace of solution of the direct problem ( 5 )-( 7 ) on the surface γ b . let us introduce an operator of the direct problem a that reflects the given on the surface γ h dirichlet condition v ( y ) into the trace of solution of the direct problem u ( x ) on the surface γ b , multiplied by an electroconductivity coefficient k ( x ), xεγ b at fixed and equal - to - zero neumann condition on γ b : then , solution of the inverse problem is reduced to solution of an operator equation regarding an unknown function v ( y ): this functional is positive and strongly convex , and its exact lower border equals to zero . therefore , a problem of solving the equation ( 12 ) and a variation problem of finding the function v , on which the present functional reaches its minimum , are equivalents : the claimed method includes algorithms of solution of the inverse problem of electrocardiography by numerical minimization of functional ( 13 ) based on methods of gradient optimization or iteration solution of the euler equation that is a necessary condition for the minimum of functional . algorithms indicated are iteration ones , at each iteration a solution of direct and conjugate problems by the boundary element method being performed . according to hadamard , the problem of minimization of functional ( 13 ) is ill - posed because of an incorrect statement of the inverse problem of electrocardiography . the claimed method involves algorithms of solution of the inverse problem of electrocardiography by numerical minimization of functional ( 13 ) with using regularization methods based on restricting a number of iterations and on the tikhonov method . 1 . finite - element discretization of computational domain is carried out : domain ω is split into tetrahedral elements and its borders γ h and γ h — into triangle elements , functions u ( x ), xεγ h , p ⁡ ( x ) ≡ ∂ u ⁡ ( x ) ∂ n , u ⁡ ( x ) = ∑ i = 1 n ⁢ ⁢ u i · φ i ⁡ ( x ) , x ∈ γ h , p ⁡ ( x ) = ∑ i = 1 n ⁢ p i · φ i ⁡ ( x ) , x ∈ γ h u ⁡ ( x ) = ∑ i = 1 m ⁢ u i · φ i ⁡ ( x ) , x ∈ γ h , p ⁡ ( x ) = ∑ i = 1 m ⁢ p i · φ i ⁡ ( x ) , x ∈ γ h , p ⁡ ( x ) ≡ ∂ u ⁡ ( x ) ∂ n , where u i , p i , p i , u i are values of functions u ( x ), p ( x ), u ( x ), p ( x ) in surface nodes of a finite - element grid , φ i ( x ) are linearly independent finite basis functions given in nodes of a finite - element grid . 2 . vector u ={ u 1 , u 2 , . . . , u n } is computed based on an iteration procedure of numerical minimization of functional ( 21 ). 3 . function u ( x ), xεγ h , which is a final solution of the problem , is found according to the formula : the method includes the following iteration methods for finding a vector u ={ u 1 , u 2 , . . . , u n }. see : ( gill f ., murray y , wright m practical optimization [ russian translation ]. moscow : mir , 1985 ). ⁢ s ( 0 ) = - j ′ ⁡ ( u ( i ) ) , ⁢ w ( i ) = j ⁡ ( u ( i ) ) 2 j ⁡ ( u ( i - 1 ) ) 2 s ( i ) = - j ′ ⁡ ( u ( i ) ) + ω ( i ) · s ( i - 1 ) , ⁢ τ ( i ) = argmin ⁡ ( j ⁡ ( u ( i ) + τ · s ( i ) ) ) , ⁢ u ( i + 1 ) = u ( i ) + τ · s ( i ) , where i = 1 , 2 . . . , n is the iteration number , u ( 0 ) ={ u i 0 , u 2 0 , . . . , u n 0 } is an initial approximation of vector u , u ( i ) ={ u 1 i , u 2 i , . . . u n i } is the next approximation of vector u . exit from an iteration procedure is performed according to the principle of the residual ( the morozov principle ): an iteration process is stopped as soon as the following condition is reached : when reaching the iteration number , divisible by m , it is assumed that : s ( im ) = s ( 0 ) where m is an integer parameter being chosen by an experimental way . 2 . quasi - newton methods , which involve a gradient descent method ( the cauchy method ) and the newton method but use iterative ways for computing the hessian inverse matrix τ = arg min [ j ( u ( i ) − τ ( i ) · a ( i ) · j ′( u i ))], u ( i + 1 ) = u ( i ) − τ ( i ) · a ( i ) · j ( u i ) ), δ g ( i ) = j ′( u i + 1 ) )− j ′( u ( i ) ), where a ( i ) is the next approximation of the hessian inverse matrix of functional j ( u ), e is a unit matrix , b ( i ) is a correcting matrix being computed by different techniques described lower . exit from an iterative procedure is performed according to the principal of the residual ( the morozov principle ): an iteration process is stopped as soon as the following condition is reached : when reaching the iteration number , divisible by m , it is assumed that : a ( im ) = a ( 0 ) where m is an integer parameter being chosen by an experimental way . in iterative procedures described , it is necessary to compute a functional j ( u i ) and its gradient j ′( u ( i ) ). the method includes calculations of the indicated objects by the following methods . 1 . 1 . the function v ( i ) ( x ), xεω is found by solving the following mixed boundary problem for the laplace equation in an inhomogeneous medium with using the boundary element method ( the direct problem ( 5 )-( 7 )): a solution trace v ( i ) ( x ) at the border γ b is found . j ( u ( i ) ) is calculated by numerical integration according to the formula : 2 . 1 . the function g ( i ) ( x ), xεω is found by solving the following mixed boundary problem for the laplace equation in an inhomogeneous medium with using the boundary element method ( the conjugate problem ( 8 )-( 10 )): 2 . 2 . a normal derivative of solution at the border γ h : xεγ h is calculated by numerical differentiation of the found solution g ( i ) ( x ). 2 . 3 . the obtained normal derivative is multiplied by a coefficient of electroconductivity k ( y ) with inverse sign on the surface γ h : 3 . 3 . values of the obtained function p ( i ) ( x ), xεγ h in surface nodes of a boundary - element grid are assumed to be values of vector j ′( u ( i ) ): j ′( u ( i ) )=−{ p 1 ( i ) , p 2 ( i ) , . . . , p n ( i ) }. the method involves a minimization of functional ( 13 ) based on numerical solution of the euler equation that is the necessary condition for the minimum of functional ( 13 ): solving the equation ( 11 ) is performed on the basis of the following iteration algorithm : where ū is a random initial approximation , i = 1 , 2 , . . . , n is the iteration number , τ is a parameter of an iterative method . exit from an iterative procedure is implemented according to the principle of the residual : u i + 1 ) ( y )= u ( i ) ( y )− γ ·( p ( i ) ( y )− p ( y )), calculations of the function p ( y )= a * u ( x ) are carried out by the following way . 1 . with using the finite element method , the following mixed boundary problem for the laplace equation in an inhomogeneous medium is solved ( the conjugate problem ( 8 )-( 10 )): 2 . by numerical differentiation of the found solution g ( y ), a normal derivative of solution at the border γ h : γ h : ∂ g ⁡ ( y ) ∂ n , 3 . the obtained normal derivative is multiplied by a coefficient of electroconductivity k ( y ) with inverse sign at the border γ h : calculation of functions p ( y )= a * v ( i ) ( x ) is carried out by the same way . 1 . with using the finite element method , the following mixed boundary problem for the laplace equation in an inhomogeneous medium is solved at each iteration ( the conjugate problem ( 8 )-( 10 )): 2 . by numerical differentiation of the obtained solution g ( y ), a normal derivative of solution is calculated at the border γ h : γ h : ∂ g ( i ) ⁡ ( y ) ∂ n , 3 . the obtained normal derivative is multiplied by a coefficient of electroconductivity k ( y ) with inverse sign at the border γ h : γ h : p ( i ) ⁡ ( y ) = - k ⁡ ( y ) · ∂ g ( i ) ⁡ ( y ) ∂ n , calculation of the function v ( i ) ( x )= a · u ( i ) ( y ) is performed by the following way . 1 . the function v ( i ) ( x ), xεω is found by solving the following mixed boundary problem for the laplace equation in an inhomogeneous medium with using the boundary element method ( the direct problem ( 5 )-( 7 )): 2 . a trace of the obtained solution at the border is multiplied by a coefficient of electroconductivity k ( x ) at the border γ b : v ( i ) ( x )= k ( x )· u ( x ), xεγ b . the method involves a minimization of functional ( 13 ) with using the tikhonov regularization on the basis of solving the corresponding euler equation : solving the equation ( 15 ) is implemented based on an iterative procedure : u ( i + 1 ) = u ( i ) − τ ·( a * ·( a · u ( i ) + α · u ( i ) − a *· u ). exit from an iterative procedure is performed when the following condition is reached : | u i + 1 − u i |& lt ; ε where ε is a small positive parameter depending on the machine accuracy . the choice of a regularization parameter α is carried out according to the principle of the residual : such α is chosen at which the following equality is fulfilled the most exactly : where u ( α ) is a parameter α - depending solution obtained as a result of implementing the iterative procedure . block - diagrams of algorithms are shown in fig1 , 12 , 13 . fig1 gives convergence diagrams of a quasi - newton iterative procedure of davidon - fletcher - powell ( 14a ) and of iteration solution of the euler equation ( 14b ). in calculations , a model of the torso and heart of a real patient was used . for modeling the standard electric field of the heart , a quadruple source to be placed in the geometric center of the heart was used . fig1 gives imposed on realistic models of the heart isopotential maps of the exact electric potential ( 14a ) calculated by the disclosed in the present invention algorithm with taking into account an electrical inhomogeneity of the chest ( 14b ) and by an algorithm based on a homogeneous model of the chest and disclosed in the patent - prototype of the present invention ( 14c ). fig1 shows examples of visualizing results of noninvasive electrophysiological study of the heart . 1 . constructing electrograms at interactively chosen points of the heart epicardial surface , endocardial surfaces of interventricular and interatrial septa , as well as at internal points of the chest on tomography cross - sections ( fig1 a ). 2 . constructing isopotential maps on tomography cross - sections of the chest ( fig1 b ). 3 . constructing isopotential and isochronous maps on the heart epicardial surface , endocardial surfaces of interventricular and interatrial septa ( fig1 c ). 4 . visualizing the dynamics of the myocardium excitation on the heart epicardial surface , endocardial surfaces of interventricular and interatrial septa in animation mode ( propagation maps ) ( fig1 d ). unipolar electrograms are constructed by interpolation of computed values of the heart electric field potential for all the moments of the cardiocycle at a given point . bipolar electrograms are constructed as the difference of electrograms in chosen node and at the point located in the vicinity to this node at a distance δl in the direction to i . parameters δl and 1 are interactively given . isopotential maps are constructed on the basis of bilinear interpolation of computed values of the heart electric field potential in nodes of a grid at given moment of the cardiocycle by a gradient painting method or constructing isopotential lines . for constructing isochronous maps two modes — manual and automatic — are provided . in manual mode at interactively chosen node of a grid an unipolar electrogram u ( t ), bipolar electrogram u b = u 1 ( t )− u 2 ( tt ), as well as a differential electrogram i . e ., a diagram of first derivative of an unipolar electrogram over time , are reconstructed . an operator in interactive mode marks in indicated diagrams a time - point τ corresponding to the start of the myocardium activation at a given point . the choice of corresponding mark of a time - point τ in automatic mode proceeds without operator &# 39 ; s interference . the time - point τ is determined as a maximum of a negative differential unipolar electrogram : isochronous maps are visualized on the basis of bilinear interpolation of τ values in grid nodes by the gradient painting method or constructing isochronous lines . the same data are represented in animation mode in the form of so - called excitation propagation maps . fig1 presents reconstructed by the described method epicardial isochronous maps of the extrasystole caused by an ectopic source in the region of excretory tract of the right ventricle . the mini - circle indicates a localization of the ablation electrode with the help of which a successful radio - frequency ablation of this ectopic source was implemented .
the invention relates to medicine , namely to cardiology , cardiovascular surgery , functional diagnosis and clinical electrophysiology of the heart . the invention consists in reconstructing electrograms , whose experimental registration requires an invasive access , by computational way on unipolar ecgs recorded at 80 and more points of the chest surface . based on reconstructed electrograms , isopotential , isochronous maps on realistic models of the heart are constructed , the dynamics of the myocardium excitation is reconstructed and electrophysiological processes in the cardiac muscle are diagnosed . application of the method allows one to improve the accuracy of non - invasive diagnosis of cardiac rhythm disturbances and other cardio - vascular diseases .