Electrocardiographic imaging (ECGI) is a diagnostic tool designed for the noninvasive detection of electrical activity in the heart. Mathematically, this imaging process can be modelled by an inverse problem for the coupled elliptic system with Cauchy data. For the bidomain model describing the electrical activity in the myocardium, we reconstruct the transmembrane potential on the heart surface from ECG recordings. The reconstruction process is split into two steps: firstly computing the electrical potential and current on heart surface from torso surface recordings and then recovering the transmembrane potential on the heart surface. The first step is essentially a Cauchy problem for elliptic equation which is well-known to be severely ill-posed. We realize this step in terms of boundary integral system by a quasi-regularization scheme to deal with the ill-posedness. The second step is implemented by an integral equation of the second kind with nontrivial null space. We remove the non-uniqueness and then provide a numerical scheme by the boundary element method to obtain the transmembrane potential on the heart surface. The regularized solution obtained from noisy measurement data is proven rigorously to be convergent to the exact solution as the measurement error in the torso surface tends to zero. For solving the regularizing system numerically, we establish quadrature formulas for boundary potentials for 3-dimensional bio-tissue in terms of the concept of solid angles. The numerical realizations for different configurations are finally presented to show the validity of the proposed scheme.