Reconstructing a 3D heart surface with stereo-endoscope by learning eigen-shapes

: An eﬃcient approach to dynamically reconstruct a region of interest (ROI) on a beating heart from stereo-endoscopic video is developed. A ROI is ﬁrst pre-reconstructed with a decoupled high-rank thin plate spline model. Eigen-shapes are learned from the pre-reconstructed data by using principal component analysis (PCA) to build a low-rank statistical deformable model for reconstructing subsequent frames. The linear transferability of PCA is proved, which allows fast eigen-shape learning. A general dynamic reconstruction framework is developed that formulates ROI reconstruction as an optimization problem of model parameters, and an eﬃcient second-order minimization algorithm is derived to iteratively solve it. The performance of the proposed method is ﬁnally validated on stereo-endoscopic videos recorded by da Vinci robots.


Introduction
In robotic-assisted off-pump heart surgery, heart beating considerably influences the accuracy of surgical operations, resulting in longer surgical duration and increased surgical risks.By measuring heart motion and actively synchronizing surgical instruments with this motion, a technique called active motion compensation, it is possible to provide a virtually stable operating environment to surgeons [1,2].Passive vision using stereo-endoscope, compared with other sensing techniques such as structured lighting [3,4] and ultrasound [5], is more appropriate for measuring heart motion in a Minimally Invasive Surgery (MIS) because no additional instrument ports are required, as indicated in [6].The 3D reconstruction with endoscope is also important for other advanced surgical techniques, e.g.augmented reality guidance [7], and multispectral [8] or multimodal [9] imaging.In addition, it is useful for offline applications as well, such as surgery simulation [10] and visual medical record [11], for learning, training or evaluation purposes.
However, in a highly dynamic MIS, it is very challenging to track and reconstruct Regions of Interest (ROIs) with complex soft-tissue deformations from real-time endoscopic videos.Model-based methods have been explored [12][13][14][15][16][17][18].With a parameterized deformable model, an initial ROI is warped as a template to match pixels at subsequent stereo-pair frames.Model complexity (rank) can be simply measured by the number of model parameters.Low-rank models, such as rigid and affine models, are robust and computationally efficient but difficult to deal with complex soft-tissue deformations, while high-rank models generally suffer from problems of parameter convergence and heavy computational burden, difficult to meet real-time and robustness requirements, as indicated in [16].
On the other hand, most existing models, typically based on the mesh [2,11,13,19] or spline [12,14,[16][17][18], are designed for general deformations with certain continuity and smoothness constraints, which do not take into account the statistics and specificities of ROIs.Statistical shape models have been explored to reconstruct static anatomical structures, such as bone surfaces from X-ray images [20] and employed as a priori information in a wide range of segmentation tasks from CT [21], MRI [22], 3D ultrasound [23], and optical coherence tomography [24,25].The training shapes in the literature are usually from CT databases or obtained manually by experts from pre-operative imaging.To the author's knowledge, statistical modeling in full 3D space for highly dynamic soft-tissues from stereo-endoscope is still open.It is almost impossible to establish a valid pre-operative shape database for beating heart.The shape similarity of anatomical structures from a given population, which has been used for statistical modeling in the literature, is unavailable for modeling beating heart due to high dynamics of cardiac soft-tissues.
In this work, we first derive a general 3D shape representation of ROI by decoupling a Thin Plate Spline (TPS) mapping into shape and position components.The high-rank version of the decoupled TPS is employed to pre-reconstruct a ROI for a number of frames to obtain ROI-specific shape data.Eigen-shapes are then learnt from these data, on which a low-rank 3D model is built and used to reconstruct subsequent frames.The linear transferability of Principal Component Analysis (PCA) is proved, which significantly reduces the computational burden of eigen-shape learning and thus makes real-time statistical modeling possible.A general model-based reconstruction framework with an iterative algorithm for estimating model parameter is presented.The feasibility and effectiveness of the proposed method is finally validated on stereo-endoscopic videos recorded by da Vinci (Intuitive Surgical, CA) robots.

Representing 3D shape with TPS
The TPS has been widely used in image alignment and shape matching, and recently used to handle soft-tissue deformations [14,17].We extend and decouple the classical TPS mapping into two independent components related to shape and position, respectively.
Let m ∈ R 2 denote pixel coordinates.The TPS of K Control Points (CPs) is written as where ϕ (m) is a row vector of elements {ϕ ( m − c k )} K k=1 with ϕ (r) = r 2 ln r and c k being the k-th CP, and {α, β} are parameter vectors associated with non-affine and affine transformations, respectively.All vectors are column vectors, unless otherwise specified.To obtain squareintegrable second derivatives, the following constraint is imposed.
where x denotes the augmented vectors of x by adding 1 as the last element.By substituting it into Eq.( 1), we can re-parameterize the non-affine part of the TPS as ϕ (m)α , where α is a K − 3 dimensional subvector of α (three elements of α are eliminated).Let p ∈ R 3 denote 3D coordinates in the real space.By separately mapping in the x, y, z directions, the TPS can be extended to R 2 → R 3 that maps m to p with the parameter vectors {α x , β x , α y , β y , α z , β z }, which contain 3K model parameters (scalars) in total.
Let R = {m n } N n=1 be a ROI.Assume that m o is the center of R, and its 3D counterpart p o can denote the position of R. The relative 3D coordinates of N pixels with respect to p o can represent the 3D shape of R that is position-independent.Let p = p − p o , where b(m) = ϕ (m) − ϕ (m o ), (m − m o ) T , and θ x/y/z is a parameter vector stacked by α x/y/z and the first two elements of β x/y/z .Accordingly, the 3D shape of the ROI is represented as where s is a N -dimensional shape vector stacked by {p n } N n=1 with N = 3N, and B is a N × K transformation matrix stacked by {B(m n )} N n=1 with K = 3(K − 1).The shape and position parameters {θ, p o } together determine a 3D surface in the real space, and the decoupled TPS model can be rewritten as p = B(m)θ + p o .

Low-rank shape modeling based on eigen-shapes
With enough CPs, the TPS model can accurately describe any 3D surface.Such general models, however, are inefficient for the ROI reconstruction of beating heart.Since heart motion is normally quasi-periodic with stable statistics [26], it is possible to learn a low-rank model for a specific ROI from its past shapes.For this, the PCA is applied to the shape data.
Assume that {s l } L l=1 are past shapes.A shape data matrix of row-wise zero mean is built: The PCA of S can be computed by using the Eigen Decomposition (ED): where Λ is a diagonal matrix of the eigenvalues {λ j } N j=1 sorted in descending order, and the columns of U are the corresponding eigenvectors.Each eigenvector identifies a specific shape in the 3D space, so we call it eigen-shape.
According to the Eckart-Young theorem, a data matrix can be optimally rebuilt on the first few eigenvectors.Assume that Ŝ is the rebuilt matrix on the first J eigenvectors.The rebuilding error defined with the Frobenius norm is Since the eigenvectors corresponding to the largest eigenvalues are reserved, information loss caused by Ŝ is minimal of all possible J-rank rebuilding.In fact, the rebuilt data may be closer to ground truth, as the original data are usually contaminated by measurement noise, which is more likely to rest on the discarded eigenvectors.
For quasi-periodic heartbeat, the shape vectors of a ROI at different moments are correlated and sparse-distributed in the high-dimensional shape vector space.Accordingly, the past shapes of the ROI can be recovered accurately on the principal eigen-shapes.With the assumption of stable statistics, the shapes of the ROI in the future can also be reconstructed accurately with where U J is the submatrix of U formed by the first J eigenvectors, and w is a new parameter vector of J dimension.Compared with the general shape model in Eq. ( 4), the new model in Eq. ( 8) is more efficient as statistics of the ROI are introduced.However, the PCA for a large matrix is computationally expensive in both time and memory.The complexity for analyzing S is O N L 2 by assuming L ≤ N .It is impractical for most real-time systems, as N and L are very large for a typical reconstruction task.Here, we develop a fast eigen-shape algorithm, which reduces the complexity to O LK 2 where K L.
Since rank (B) = K , for the shape matrix S there are at most K valid eigen-shapes that correspond to non-zero eigenvalues.Assume that Θ is the shape-parameter matrix of S, i.e. S = BΘ.Let U Θ Λ Θ U T Θ be the ED of ΘΘ T , then it can be easily proved that: Proposition 1: Λ Θ and BU Θ are the non-zero eigenvalues and eigenvectors of SS T , respectively, if and only if B is column-orthonormal.
The proposition shows that with an orthonormal transformation matrix (called transferring condition), the PCA can be transferred from a large data matrix to a small parameter matrix.Unfortunately, the shape transformation matrix B of the TPS is not orthonormal, which distorts the mapping from the eigen-parameters (U Θ ) to eigen-shapes (U).
To overcome this problem, we orthonormalize B with the QR decomposition B = B R, and re-parameterize Eq. ( 4) as s = B θ .Correspondingly, the parameter matrix Θ is built.By using Proposition 1, the eigen-shapes can be computed fast with B U Θ , and the eigenvalues are obtained from Λ Θ directly.The computational complexity is hence reduced to O(LK 2 ).
Though proposed for shape modeling, Proposition 1 is useful to other PCA tasks associated with linear systems.We can extend it from the point view of function composition: Proposition 2: Given two real vector spaces that are isomorphic, there is a linear function U X, mapping one vector set (matrix) X in U to another in V, where B U/V is a transformation matrix formed by an orthonormal basis of U/V.Let P : X → U X and S : X → Λ X be the functions of calculating eigenvectors and eigenvalues from a vector set, respectively.Then there are P • F = F • P and S • F = S for any vector set in U.
The proof is given in Appendix A. Proposition 2 shows that between any two isomorphic real vector spaces, we can always find an invertible linear transformation (bijection) that is commutative with the eigenvector operation and transparent to the eigenvalue operation, such that we can map the eigenvectors from one space to another and keep the eigenvalues unchanged.

Model-based tracking and reconstruction
The past shapes can be collected with high computational cost and even manual intervention by employing a high-rank TPS for pre-reconstruction.The shape transformation matrix of the TPS can be orthonormalized offline because it is constant for any ROI of given size.
The complete steps of the proposed reconstruction method are summarized as follows: 1) Given a ROI (R) in the template image, the shape transformation matrix of a K-CP TPS is constructed and orthonormalized with B = BR −1 .
2) The first L frames are pre-reconstructed with the orthonormalized TPS model, written as where B (m) is the submatrix of B associated with m.As a result, {θ l } L l=1 are saved.
3) The parameter matrix Θ is built and decomposed by where 5) The subsequent frames from L + 1 are reconstructed with the new model in Eq. (10).With a calibrated stereo-endoscope, the ROI reconstruction involved in the steps 2) and 5) can be converted into an optimization problem, i.e., finding the best model parameters that minimize the alignment error of ROI between the template image and the current stereo-pair images.
We discuss the problem, without loss of generality, based on a general R 2 → R 3 model: where ξ is the parameter vector, T (m) and d (m) are the transformation matrix and offset vector associated with m, respectively.Many existing models can be written in this form.For linear models, both T (m) and d (m) are linear with respect to m, such as the affine model [2], while for nonlinear models at least one of them is nonlinear, such as the TPS, B-Spline [12], Piecewise Bilinear Mapping (PBM) [13], as well as the new statistical model in Eq. ( 10).The TPS in Eq. ( 9), for instance, is obtained by letting T (m) = [B (m), I (3) ], d(m) = 0, and ξ be a combined vector of θ and p o , where I (3) denotes a 3-dimensional identity matrix.
As illustrated in Fig. 1, given a pixel m in the template image, its 3D coordinates p in the real space are determined by the model parameter ξ to be estimated by using Eq.(11).According to the pinhole camera model, the perspective projection of p onto the left/right camera plane of the stereo-endoscope can be represented by the homogeneous coordinates where C L/R is the 3×4 left/right camera matrix, which can be calibrated with Zhang's method [27] before surgery.The corresponding pixel coordinates can be calculated by where H maps homogeneous coordinates to pixel coordinates, with an ideal parameter ξ * .However, in most cases, the ideal parameter does not exist, involving the solution of an overdetermined nonlinear system of equations.With the sum of squared pixel differences.The reconstruction problem can be formulated as arg min where I L/R and I T are N-dimensional intensity vectors stacked by , respectively.Regularization items can be added in Eq. ( 16) for robust solution, as in [28].
The Efficient Second-order Minimization (ESM) technique [29] is used to solve Eq. ( 16), which has a wider convergence basin and a faster convergence rate than traditional minimization algorithms.In the ESM, the model parameter is updated iteratively at each frame, ξ ← ξ + ∆ξ, until ∆ξ is less than a preset threshold.The increment is calculated by where ∇ x Y denotes the Jacobian matrix of Y at x, and "+" denotes the matrix pseudo-inverse.
The computations of the Jacobians in Eq. ( 17) are given in Appendix B.
It should be noted that the proposed modeling method as well as the reconstruction algorithm, though derived from the TPS, apply equally to other deformable models.The shape representation in Eq. ( 4) and thus the statistical shape model in Eq. ( 8) can be derived from any R 2 → R 3 deformable model that can be written in the form of Eq. ( 11), including the B-spline and the PBM mentioned before.In fact, the R 2 → R B-spline based statistical shape model used for tissue segmentation in [23] is sub-optimal, because it is learnt simply from the parameter spaces without orthonormalizing the transformation matrix, which distorts the mapping from eigen-parameters to eigen-shapes, as proved in Proposition 1.

Results
Taking shape statistics into account, the low-rank model can be efficient in computation and robust against divergence of model parameter, when incorporated into the model-based reconstruction framework.Four stereo-endoscopic videos captured by the da Vinci robots are used to validate the proposed method.Each video consists of 800 stereo-pair frames with a frame rate of 25 Hz.As shown in Fig. 2, the videos I and II record the motion of a phantom heart model [30,31], and the videos III and IV are in vivo beating hearts in real TECAB procedures [32].Color images are converted into 8-bit grayscale for efficient processing.

Sparsity and statistical stability of shape data
The rank of eigen-shape based model (J) is adjustable according to the efficiency and robustness requirements of reconstruction tasks in practice, as well as the sparsity of shape data itself, which can be indicated by the eigenvalues.A guideline for choosing J is given in our experiments.
To obtain shape data, an orthonormalized 9-CP TPS is employed to reconstruct a 120 × 120 ROI on each video for 600 frames.At each frame, the ESM is used to estimate the shape and position parameters iteratively.Specular reflections caused by glossy tissue surfaces are removed by using a simple thresholding method [14].The reconstruction flow is restarted manually when it is interrupted.We finally obtain four 24 × 600 shape-parameter matrices, which identify four hyper-high dimensional shape matrices (43200 × 600).
The eigenvalues and eigen-shapes are computed from shape-parameter matrices by using Proposition 1.As shown in Fig. 3(a), the eigenvalues on each video are sparse, indicating that there is obvious correlation between the shape vectors.The 24-dimensional shape subspace spanned by the 9-CP TPS is redundant for the ROI.Therefore, it is feasible to rebuild the shapes in a lower-dimensional space spanned by few eigen-shapes.
According to Eq. ( 7), the rebuilding error can be indicated more intuitively with the Signal-to-Noise Ratio (SNR) γ and the Root-Mean-Square Error (RMSE) σ p , where  To show the statistical stability of shape data, we compared the PCA results computed from the first 600 frame and from all 800 frames.Figure 4 shows the corresponding mean shapes and first eigen-shapes on each video.The differences between the corresponding shapes are trivial (of the order of 0.01 mm/pixel), showing that the statistics of heartbeats are stable and thus reconstructing future frames with the statistical model trained from past frames is feasible.

Reconstruction experiments
In the reconstruction experiments, the first L = 600 frames in each video, serving as training frames, are pre-reconstructed with the 9-CP TPS, while the subsequent 200 frames, serving as Fig. 4. Mean shapes and 1st eigen-shapes calculated from all 800 frames vs. from the first 600 frames.For each video, the first row shows the mean shapes, the second row the 1st eigen-shapes; the left shapes are from all 800 frames, the right shapes from the first 600 frame.
test frames, are reconstructed with the eigen-shape based Statistical Deformable Models (SDMs).

Qualitative analyses
Since the ground truth is unknown, reconstructed ROI is usually evaluated qualitatively by visually inspecting its concordance with heart beating and by analyzing its motion trajectories.
Figure 5 shows some reconstructed frames in each video.The SDMs can successfully handle soft-tissue deformations on beating hearts.The 3D shapes reconstructed at various frames are consistent with those perceived by human vision, and the center of the ROI is located accurately.
For comparison, the corresponding frames reconstructed by a 4-CP TPS on the video I are given in Fig. 6.Even though the 4-CP TPS has 12 parameters more than those (8 eigen-shape weights + 3 position parameters) of the SDM learnt for the video I, the 3D surfaces represented by the 4-CP TPS appear to be over-smooth with many 3D details missed out, comparing with those in Fig. 5(a), and therefore the advantages of the SDM are demonstrated.
The motion trajectories of the ROIs on the phantom videos (I and II) and their Power Spectral Densities (PSD) are shown in Fig. 7.The trajectories of p o along x, y, z axes exhibit highly regular periodicity, and their power concentrates at 1.5 Hz and its second (3.0 Hz) and third (4.5 Hz) harmonics, which is consistent with the given rate of phantom heart beating.

Quantitative evaluation
The SDMs are evaluated quantitatively and compared with well-established deformable models based on the ESM solution in terms of robustness, accuracy and speed.The reference models  Robustness: The ROI reconstruction with stereo-endoscope is very challenging.The processing chain for continuous image sequence is often interrupted by dynamic effects in the MIS.As recursive structure is usually adopted for efficiency, in which the final estimate at the current frame serves as an initial estimate for the next frame, reconstruction algorithms normally cannot recover by themselves when interrupted and must be restarted.The robustness of reconstruction algorithms can be measured by the number of frames processed continuously.
All 200 test frames on the phantom videos can be continuously processed by all investigated models, due to ideal imaging conditions and less dynamic disturbance.On the in vivo videos, the SDMs and the QST can take all 200 test frames, performing much better than the 9-CP TPS, 9-CP PBM and 9-CP B-Spline models which take averagely 58, 47 and 42 frames, respectively.It verifies that the high-rank models with too many parameters are susceptible to dynamic effects in the MIS, making the ESM prone to divergence, while the low-rank models with fewer parameters are easier to find optimal solutions in the noisy environments.The SDMs also perform better than the 5-CP TPS (taking averagely 152 frames) on the in vivo videos, though the ranks of these models are identical or similar.The result shows that for a specific ROI, the solution-space designed online by statistical modelling is more efficient than the fixed solution-space spanned by a general deformable model of identical or similar rank.
Accuracy: For evaluating accuracy, natural landmarks in the ROIs are manually tracked in both left and right image sequences.The 3D coordinates of the landmarks obtained by the model-based reconstruction are projected back onto the left and right images, and compared with the manually tracked results.To ensure the comparability of the results, interrupted algorithms are restarted manually.The joint pixel error is calculated at each frame where m L/R denotes the pixel coordinates projected onto the left/right image and ḿL/R denotes those manually tracked.The Mean Error (ME) and Standard Deviation (SD) of over 200 test frames on each video are shown in Table 1.
On the phantom videos, good results are obtained by all models, as in the robustness test.The high-rank (9-CP) models have a slight advantage over the low-rank models (4-CP TPS, SDMs and QST).However, the SDM of 8/6 eigen-shapes on the video I/II performs better than the 4-CP TPS and QST, though have similar ranks.On the in vivo videos, dynamic effects in the MIS complicate the reconstruction task.The best results are obtained by the SDMs, showing that the proposed method is efficient for handling complicated deformations on in vivo hearts.Although the SDMs are trained from the shape data obtained by the 9-CP TPS, they are still able to achieve more accurate results than the latter in in vivo settings, thanks to the intrinsic advantage that the eigen-shapes extracted by the SDMs are insensitive to small noises in the training data.
Speed: The computational cost of reconstruction algorithms depends on the rank of the employed model.The most time-consuming operation in the ESM iteration is the pseudo-inverse calculation of the Jacobian matrix in Eq. ( 17).For a K-rank model on a N-size ROI, the computational complexity for one frame is O(µNK 2 ), where µ is the number of the ESM iterations.Typically, on the in vivo videos, the average µ of the 12 eigen-shape SDM is 6.8, close to that of the QST and about half those of the 9-CP models, confirming that models with fewer parameters are easier to converge.For the same ROI, the computational cost of a typical 12 eigen-shape SDM is nearly 15% that of a 9-CP model.
To evaluate practical speed, the average CPU time is computed based on our experimental platform (Matlab R2012a, Intel Xeon E3-1231 CPU and 8GB RAM).The speed of the 12 eigen-shape SDM is 2.9 s/frame, much faster than 9.2 s/frame of the 9-CP TPS.By means of statistical modeling, real-time reconstruction is feasible if Graphics Processor Units (GPU) and more efficient programming languages (e.g.C/C++) are adopted.
Finally, the SDMs built with different SNR guidelines (γ > 30 dB and γ > 10 dB) are tested.For γ > 30 dB, the first 12/10/19/16 eigen-shapes are reserved on the video I/II/III/IV.There is not significant improvement in the reconstruction accuracy though more eigen-shapes are employed, while the robustness and speed are decreased obviously.Their performances are between those of the SDM (γ > 20 dB) and the 9-CP TPS.For γ > 10 dB, only the first 4/2/5/5 eigen-shapes are reserved on the video I/II/III/IV, which results in faster and more robust reconstructions.However, the reconstructed surfaces are inaccurate due to too few shape degrees of freedom.The tracking errors are even larger than those of the QST in Table 1.
Overall, high-rank models can render tissue surfaces more accurately in theory; however, in practice they often confuse parameter estimation algorithms due to too many degrees of freedom.For general deformable models based on spline or mesh, increasing CPs will result in a quadratic growth of computational burden and cause model parameters difficult to convergence, while reducing CPs may lead to over-smoothing and loss of 3D details.It is a dilemma for these general models to choose an appropriate rank, since they were designed for handling general surfaces without taking into account the statistics and specificities of the ROI.The proposed method  The Jacobian cannot be computed directly because the ideal parameter ξ * is unknown.An effective approximation is derived by using the template image and the current parameter ξ.Given a parameter ξ, the mapping M ξ : m → m Y is normally invertible, i.e., there is m = M −1 ξ m Y .By using Eq. ( 15), the n-th row of ∇ ξ * I Y can be written as and similarly decomposed into can be approximated by letting ξ * = ξ.In Eq. ( 28), ∂T(m)/∂η (with η denoting u or v for convenience) and ∇ m d(m) are constant and can be pre-computed for a given ROI.Take the statistical model in Eq. ( 10) for example:

Fig. 1 .
Fig. 1.Illustration of model-based ROI reconstruction Let I T denote the intensity function of the template image and I L/R be that of the current left/right image.For the given ROI, there should be

Figure 3 (
Figure 3(b)  gives the SNR and RMSE curves with respect to J.In our study, we choose the least J that ensures γ > 20 dB.Following this guideline, the first 8/6/12/10 eigen-shapes are reserved for statistical modeling on the video I/II/III/IV, corresponding to the RMSE of 0.048/0.072/0.070/0.094mm/pixel, much lower than the measurement error (about 0.1 mm) of the stereo-endoscope image system itself and thus ignorable.

Fig. 3 .
Fig. 3. PCA of shape data on videos I-IV

Fig. 7 .
Fig. 7. Motion trajectory analysis of ROI (p o ) tracked by the statistical deformable model.

B. 2 .
where: a) ∇ m Y n I Y is the gradient (1 × 2) of the image I Y at m Y n .b) ∇ p Y n H is the derivative of function H at p Y n .From Eq. (14), one derives ∇ x H = p p Y ≡ [C Y ] 3 and ∇ ξ p| m n ≡ T(m n ) are constant for a given ROI.Computation of ∇ ξ * I L/R where: a) ∇ m n I T is the gradient of I T at m n .b) ∇ p p Y and ∇ ξ p| m n are constant, independent of ξ * , the same as in Eq.(23).c)∇ m Y n M −1 ξ * is the inverse of the Jacobian ∇ m n M ξ * = ∇ p Y n | ξ * H ∇ p p Y ∇ m n p| ξ * ,(27)where p Y n | ξ * , the ideal homogeneous coordinates determined by ξ * , can be approximated by p Y n , and ∇ m n p| ξ * , the derivative (3 × 2) of the deformable model at m n with ξ * , for m = [u v] T , ∇ m n p| ξ * = ∂T(m) ∂u u n ξ * ∂T(m) ∂v u n ξ * + ∇ m n d(m), to Eq. (3), ∂B(m)/∂η is formed diagonally with and C 2 be two complementary submatrices of the left matrix in Eq. (2), where C 1 is 3 × 3 invertible, corresponding to the three eliminated elements of α, and C 2 corresponds to α .Then ϕ (m) = ϕ (m) of the k-th radial basis kernel in ϕ (m) can be pre-computed by using∂ϕ ( m − c k ) ∂u ∂ϕ ( m − c k ) ∂v = (2 ln m − c k + 1) (m − c k ) T (33)