

ORIGINAL ARTICLE 

Year : 2019  Volume
: 9
 Issue : 1  Page : 4249 

Sensitivity uniformity ratio as a new index to optimize the scanning geometry for fluorescent molecular tomography
Anita Ebrahimpour^{1}, Seyed Salman Zakariaee^{2}, Marjaneh Hejazi^{1}
^{1} Department of Medical Physics and Biomedical Engineering, Faculty of Medicine, Tehran University of Medical Sciences, Tehran, Iran ^{2} Department of Medical Physics, Faculty of Paramedical Sciences, Ilam University of Medical Sciences, Ilam, Iran
Date of Web Publication  6Sep2019 
Correspondence Address: Dr. Marjaneh Hejazi Department of Medical Physics and Biomedical Engineering, Faculty of Medicine, Tehran University of Medical Sciences, Tehran Iran Dr. Seyed Salman Zakariaee Department of Medical Physics, Faculty of Paramedical Sciences, Ilam University of Medical Sciences, Ilam Iran
Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.JMSS_22_18
Background: Molecular fluorescence imaging is widely used as a noninvasive method to study the cellular and molecular mechanisms. In the optical imaging system, the sensitivity is the change of the intensity received by the detector for the changed optical characteristics (fluorescence) in each sample point. Sensitivity could be considered as a function of imaging geometry. A favor imaging system has a uniform and highsensitivity coefficient for each point of the sample. In this study, a new parameter was proposed which the optimal angle between the source and detector could be determined based on this parameter. Methods: For evaluation of the new method, a twodimensional mesh with a radius of 20 mm and 5133 nodes was built. In each reconstruction, 0.5mm fluorescence heterogeneity with a contrasttopurpose ratio of fluorescence yield of 10 was randomly added at different points of the sample. The source and the detector were simulated in different geometric conditions. The calculations were performed using the NIRFAST and MATLAB software. The relationship between mean squared error (MSE) and sensitivity uniformity ratio (SUR) was evaluated using the correlation coefficient. Finally, based on the new index, an optimal geometrical strategy was introduced. Results: There was a negative correlation coefficient (R = −0.78) with the inverse relationship between the SUR and MSE indices. The reconstructed images showed that the better image quality achieved using the optimal geometry for all scanning depths. For the conventional geometry, there is an artifact in the opposite side of the inhomogeneity at the shallow depths, which has been eliminated in the reconstructed images achieved using the optimal geometry. Conclusion: The SUR is a powerful computational tool which could be used to determine the optimal angles between the source and detector for each scanning depth. Keywords: Fluorescent molecular tomography, geometry optimization, NIRFAST, sensitivity matrix, sensitivity uniformity ratio
How to cite this article: Ebrahimpour A, Zakariaee SS, Hejazi M. Sensitivity uniformity ratio as a new index to optimize the scanning geometry for fluorescent molecular tomography. J Med Signals Sens 2019;9:429 
How to cite this URL: Ebrahimpour A, Zakariaee SS, Hejazi M. Sensitivity uniformity ratio as a new index to optimize the scanning geometry for fluorescent molecular tomography. J Med Signals Sens [serial online] 2019 [cited 2020 Feb 26];9:429. Available from: http://www.jmssjournal.net/text.asp?2019/9/1/42/253747 
Introduction   
Molecular fluorescence imaging is widely used as a noninvasive method to study the cellular and molecular mechanisms. The elastic light scattering in the biological tissues is one of the main causes of the information destruction which significantly reduced the quality of the reconstructed images. Therefore, different methods were used to compensate information loss at different steps of image acquisition,^{[1],[2],[3]} forward, and inverse stages.^{[4],[5]} The optimization of the sourcedetector geometry is one of the most practical and effective methods that have been used in recent studies.^{[6],[7]} Various methods, such as singularvalue analysis ^{[8],[9]} and orthogonality of the Jacobian matrix, were used to evaluate the sourcedetector geometries and to optimize sampling frequency, the field of view, etc.^{[10],[11]} However, the singularvalue decomposition and orthogonal methods could not assess different geometries. They did not change steadily with the changes of the sourcedetector geometry.^{[12]} In many studies, the uniformity of the sensitivity matrix has been used to evaluate different geometries. The sensitivity matrix, which depends on the imaging geometry and sampling strategy, is the sum of the rows of the Jacobian matrix. This matrix represents the sensitivity of the sample points for each positioning of source and detector elements. In pioneer studies, the uniformity of the sensitivity matrix was qualitatively evaluated by a simple profile.^{[13],[14],[15]} Then, the nonuniformity of the sensitivity matrix was quantitatively calculated using the Laplacian operator as a proper method for evaluating different sourcedetector geometries.^{[12],[16]} In recent years, new methods, such as CramerRao Lower Bound and Effective Independence, have been used to optimize sourcedetector geometry, which is relatively more sophisticated methods.^{[17],[18]} In this study, a new parameter was introduced which is applied to the sensitivity matrix. Both the mean sensitivity and the variation of the sensitivity are taken into account in determining the best geometry using the proposed parameter. It seems that the uniformity of the sensitivity matrix is not a comprehensive criterion to determine the best sourcedetector geometry. Achieving high sensitivity with maintaining the sensitivity uniformity is the essential goal in this method. A favor imaging system has the uniform highsensitivity coefficient for each point of the sample. The optimal angle between the source and detector could be determined based on the proposed parameter. Different depths of the sample have been separately studied, and the optimal scanning conditions were determined.
Materials and Methods   
Sensitivity uniformity ratio
In this study, sensitivity uniformity ratio (SUR) was introduced to determine the optimal geometry (arrangement of the source and detector around the sample) in the optical imaging system. The computational algorithm of the SUR index is summarized in [Figure 1].  Figure 1: The computational steps of the sensitivity uniformity ratio index
Click here to view 
For the sourcedetector geometry, the Jacobian matrix maps differential changes in the spatial distribution of fluorescent sources to differential changes in the fluorescent light captured at the surface [Figure 2]. The Jacobian elements were defined using a Green's function as the Eq. 1:^{[16],[19]}  Figure 2: The schematic mapping for each row of Jacobian and sensitivity matrices on the mesh space. The Jacobian matrix has the row numbers equivalent to the number of sourcedetector pairs. The columns number of the Jacobian matrix is equivalent to the number of points defined for the sample simulation (the number of mesh nodes). The number of each element indicates the sensitivity of the corresponding point (the change in the intensity received by the detector for the changed optical characteristics in each point) in the imaging system. If there are a source and one detector around the source, the imaging sensitivity could be seen in the form of a banana in the sample space, which is a row of the Jacobian matrix overlaid on the sample. The sensitivity matrix (s_{j}) is calculated based on the raw summation of the Jacobian matrix
Click here to view 
Where φi is the light intensity measured at the i^{th} projection, μ_{j} is the optical absorption at the j^{th} spatial location, η is the quantum efficiency of the fluorophore, and r_{si} and r_{di} are the spatial location of the source and detector, respectively. For the ith optical projection, G ^{x} and G ^{m} are the Green's functions for light propagation at the excitation and emission wavelengths, respectively.^{[16]}
The Jacobian matrix represents the sensitivity of the sample points (nodes) for all sourcedetector pairs. Each column of the Jacobian matrix shows the sensitivity of a point for all sourcedetector pairs. For each point, the overall sensitivity is determined by the sum of row elements as follows:^{[16]}
Where s_{j} is the sensitivity matrix, m is the number of sourcedetector pairs.
The normalized sensitivity matrix was calculated using the Eq. 3:
Where s_{j} (n) is the normalized sensitivity matrix.
The mean normalized sensitivity was calculated based on the normalized sensitivity matrix, as the following Eq. 4
Where is the mean normalized sensitivity and n is the voxel numbers.
The mean normalized sensitivity represents the overall sensitivity of the subject for a specific scanning geometry. The higher number means the greater sensitivity, and therefore, increasing of the detection capability.
In the applied imaging geometry, the sensitivity nonuniformity is very important and should be reduced as much as possible. The following equation is defined to calculate the mean sensitivity nonuniformity, ū:
Gradient operator (∇) is used to calculate the sensitivity changes in different nodes, which is based on the sensitivity difference between the adjacent nodes. A final vector with n − 1 elements is obtained (including the sensitivity of adjacent nodes), after determining the absolute magnitudes. The mean changes were achieved by dividing by n − 1. This parameter could represent the average changes of the sensitivity in the sample points. In a homogeneous subject with an ideal imaging system, it must be zero. However, it is not and should be minimized.
The ratio of the sensitivity to uniformity is calculated according to the following Eq. 6:
The SUR is defined to compare different imaging geometries. Increasing the SUR quantities means the increase of the sensitivity or decreasing the nonuniformity. Hence, its efficiency and performance have been firstly investigated.
Validation of the sensitivity uniformity ratio index
For the new method, the power of the geometric evaluation was evaluated based on its relationship with the quality of the reconstructed images. NIRFAST simulation software developed at Dartmouth college (USA), University of Birmingham (UK), and Kitware Inc. (USA) was used for evaluation of the new method.^{[20],[21]} This software is based on the finite element method (a numerical method for solving problems of engineering and mathematical physics which is available at www.nirfast.org). It is an open source optical modeling package developed in Dartmouth. A twodimensional mesh with radius of 20 mm (equivalent to the radius of the mouse body) and 5133 nodes was built as the sample [Figure 3]. The nods distance was set so that their mean sensitivity statistically indicates the sensitivity of the depth, while the reconstruction time was not much increased. In each reconstruction, 0.5mm fluorescence heterogeneity with a contrasttobackground ratio of fluorescence yield of 10 was randomly added to different points of the sample. The contrasttobackground ratio was set as a high number so that the inherent contrast of the tissue could be high enough and only effective geometric factors on the sensitivity (such as the sourcedetector arrangement) play a role in the quality of the reconstructed image. For the scanning geometries, heterogeneities were randomly placed in 36 different coordinates within the sample, and a reconstruction was made for each simulated scenario. After sample preparing, the source and the detector were placed in different geometric conditions [including different angles between the source and detector as shown in [Figure 3]. To create different geometric conditions, the angle between the source and detector was changed from 0° to 180° (with a step size of 10°). Then, the forward stage and reconstruction were performed in the simulation environment. For each reconstructed image related to a specific imaging geometry, the mean squared error (MSE) was calculated as follows:  Figure 3: (a) The defined mesh and its features. (b) Five examples of the sourcedetector arrangement around the sample. The source and detector were placed around the sample at different angles from 0° to 180° with a step size of 10°. The angle between the source and the detector in one projection is shown in the figures. The tomographic imaging was carried out using 36 projections
Click here to view 
Where n, and are the voxel number, the fluorescence absorption coefficient of the reconstructed image and the fluorescence absorption coefficient of the true image, respectively.
Then, the normalized sensitivity matrix and the SUR parameter were achieved. All computational procedures were performed using MATLAB software version 2008a, (The MathWorks TM, Natick, Massachusetts, USA). The relationship between MSE and SUR was evaluated using the correlation coefficient. A significant relationship between SUR and MSE could indicate the potential of SUR for the comparing and evaluation of the different imaging geometries.
Determining the optimal imaging geometry
After validation of the new parameter, it was used as a tool to determine the optimal angle between the source and detector for the different depths. The homogeneous mesh with the mentioned features was defined using the NIRFAST software [Figure 3]. Different imaging geometries were created around the mesh (in different angles of the sourcedetector from 0° to 180° with a step size of 10°). For each imaging geometry, the normalized sensitivity matrix was calculated using Eq. 3. The sensitivity matrix was mapped to the mesh space. For the mapped matrix in different imaging geometries, ten region of interests (ROIs) were selected at different depths to evaluate each depth separately. ROIs were in the donate shapes with r2r1 of 2 mm. For each scanning geometry, different depths have different sensitivity regarding quantity and uniformity. Therefore, the different depths must be separately investigated. Selected ROIs are shown in [Figure 4].  Figure 4: Defined region of interest at different depths of the sensitivity mapped on the mesh space. Each depth was separately evaluated using the sensitivity uniformity ratio, and the geometry with the highest sensitivity uniformity ratio is determined as the optimal geometry for the depth
Click here to view 
The SUR was calculated for each specific angle between the source and detector. The angle of the sourcedetector, which generates the maximum SUR, was determined as the optimal angle for the depth. There is an optimal imaging geometry for each depth. Imaging could be done using the optimal geometries for each depth. Based on this imaging strategy, the final sensitivity matrix would be achieved using Eq. 8.
Where θ_{i} (opt) is the optimized angle between the source and detector for the i^{th} depth. L is the number of separate depths considered to determine the optimal angle for the specific detectorsource geometry. S_{θi} (opt) is the sensitivity of the different angles between the source and detector.
Comparison of the optimal and conventional geometries
The quality of the reconstructed images was compared for the optimal and conventional geometries [Figure 3]. At each depth shown in [Figure 4] (0–2, 2–4, 4–6, 6–8, 8–10, 10–12, 12–14, 14–16, 16–18, and 18–20 mm), a heterogeneous media with a radius of 0.5 mm and a contrast ratio of 10 were placed. The defined subjects in the conventional and optimal imaging geometries were simulated using the NIRFAST software. For statistical evaluations, heterogeneities were randomly placed at six points, and simulations were repeated six times. MSE of the reconstructed images was compared for two geometries (Eq. 7).
Results   
For the reconstructed images, there was a negative correlation coefficient (R = −0.78) with the inverse relationship between the SUR and MSE indices. The comparison of the SUR and MSE indices is shown as a bin scatter plot in [Figure 5].  Figure 5: Bin scatter plot between the mean squared error of the reconstructed images and sensitivity uniformity ratio applied to the sensitivity matrix in different sourcedetector geometries. There was a high correlation coefficient (R = −0.78)
Click here to view 
The achieved correlation coefficient shows the inverse relationship between SUR and MSE, which confirms the benefits of SUR for geometric evaluation. The calculation of SURs at different angles between the source and detector shows that for each depth, the maximum SUR could be achieved at a particular angle between the source and detector [Figure 6]. For sample scanning with appropriate sensitivity, the optimal angles between the source and detector must be used. Therefore, new geometry was proposed.  Figure 6: The sensitivity uniformity ratio magnitudes for the different angles between the source and detector at ten separate depths
Click here to view 
For all studied depths, the sensitivity map of the proposed geometry shows the higher uniformity and sensitivity than that of the conventional geometry. These results are clearly seen in the histograms [Figure 7].  Figure 7: (a) and (b) The sensitivity map in the conventional (which sourcedetector is placed in the opposite direction [180°]) and proposed geometries, respectively. (c) Comparison of the sensitivity histograms for the conventional and optimal geometries. In the optimal geometry, the sensitivities are focused on higher numbers and have a limited distribution with lower diversity
Click here to view 
For all scanning depths, the comparison of the reconstructed images showed that the better image quality achieved using the optimal geometry [Figure 8]. Reconstructed images using both geometries are shown in [Figure 9]a. For the conventional geometry, there is an artifact in the opposite side of the inhomogeneity at the shallow depths. This artifact has been eliminated in the reconstructed images achieved using the optimal geometry.  Figure 8: Mean squared error of the reconstructed images obtained using the conventional and optimal geometries. The horizontal axis denotes the scanning depth. The sample space is divided into ten depths with a 2mm gap. Depth numbers increased from the center to the surface
Click here to view 
 Figure 9: (a) The first row of the images represents the true figures. The second row represents the fluorescence yield (mm^{− 1}) of the reconstructed images for the conventional geometry. The artifact could be seen at the shallow depths in opposite side of the inhomogeneity (artifact location is indicated using the arrow). The last row is the fluorescence yield (mm^{− 1}) of the reconstructed images for the optimal geometry, which the reference geometryrelated artifacts have been eliminated. Reconstructions have been made by placing inhomogeneities at ten separate depths, which each depth indicated using a number. (b) The sensitivity map on the mesh for the different angles between the source and detector. The numbers below the images denote angle between the source and detector (in degrees)
Click here to view 
Discussion   
In many studies, the sensitivity of an imaging procedure was evaluated based on the sensitivity uniformity.^{[16]} The differences of the sensitivity coefficients in different points of the sample [Figure 9]b could result in artifacts or misinterpretation of results. In addition, the sensitivity coefficients at different points of the sample must be high enough. More signals would be detected in highsensitivity area which increases the signaltonoise ratio. In the present study, a new parameter was proposed to evaluate the sourcedetector geometries. The SUR was defined based on the sensitivity matrix because this matrix is a function of the imaging geometry. The SUR represents the ratio of sensitivity to nonuniformity (Eq. 6). The sensitivity and uniformity are the effective parameters on the image quality. The optimal scanning arrangement is defined as the geometry which has a uniform sensitivity matrix with a higher average value. Sensitivity variations were calculated using a gradient operator, which derived from the differential of the discrete points. The numbers of discrete points are different at the various depths. Therefore, the mean of the differential of the discrete points was considered to reduce the error caused due to the inconstant number of the nodes. For the reconstructed images, there was a negative correlation coefficient with an inverse relationship between the proposed index and the MSE (R = −0.78), which confirms the accuracy of the parameter for geometric evaluation of the imaging system [Figure 5]. For the geometric conditions, the sensitivity is different for each depth of the sample [Figure 6]. It is not possible to select an overall geometry that has been optimal for all depths. Therefore, for each depth of the sample, the optimal angle between the source and detector was separately determined.
For each depth, the angle between the source and detector with maximum SUR was chosen. Therefore, the optimal angles related to the specific depths could be used for sample scanning. The sensitivity map of the new geometry has higher uniformity and sensitivity than that of the conventional geometry [Figure 7]. The comparison of the quality of the reconstructed images with the conventional and optimized geometries confirms this hypothesis [Figure 8] and [Figure 9]. It must be noted that the using of optimal geometry increases the imaging time by a factor of L (according to Eq. 8) which could be partially compensated by increasing the mechanical speed of the system. For each depth, the tomography with the optimal geometry must be repeated L times; unlike the conventional geometry with 360° data collection which the sourcedetector is placed in the opposite direction (180°). For using the optimal geometry, system upgrading (regarding speed and time control) is a critical issue that must be addressed. In this study, the SUR parameter has been used to optimize the angle between the source and detector. However, this parameter could be used to optimize other geometric parameters including the step numbers of the sampling, the numbers of sources and detectors, the angle of detector arrangement in the scanning geometry, and the consideration of geometric parameters in reflective imaging mode. The mentioned parameters are affecting the sensitivity matrix; hence, the SUR could be used to compare the various geometric conditions.
Conclusion   
The optimization of the geometric conditions of the molecular fluorescence imaging is one of the important issues which have been widely studied. Determining the optimal angle between the source and detector is one of the most widely used methods to increase the image quality. The SUR is a powerful computational tool which could be used to determine the optimal angles between the source and detector for each scanning depth. For the fluorescence heterogeneities placed at specified depths of the sample, their images were properly reconstructed using the optimal scanning geometry.
Acknowledgment
The authors would like to thank the research affair of Medicine faculty of Tehran University of Medical Sciences.
Financial support and sponsorship
This research has been financially supported by Tehran University of Medical Sciences.
Conflicts of interest
There are no conflicts of interest.
Biographies   
Anita Ebrahimpour is Ph.D. candidate of Medical Physics in Department of Medical Physics of Tehran University of Medical Sciences, Iran. Ebrahimpour received her BSc in Radiologic technology from Iran University of Medical Sciences. She received her Msc in Medical Physics from Tehran University of medical Sciences in 2016.
Email: anitaebrahimpor@gmail.com
Seyed Salman Zakariaee was born in 1987, Sanandaj, Iran. He finished his B.S. in solidstate physics from Urmia University in 2009. Then he followed his postgraduate study in Medical Physics and he was graduated in 2013 from Tabriz University of Medical Sciences. He entered the Ph.D. program of Tehran University of Medical Sciences in the same year and finished his Ph.D. in 2018. He has been working in Medical Physics department of Ilam University of Medical Sciences as assistant professor from 2018.
Email: salman_zakariaee@yahoo.com
Marjaneh Hejazi received her BSc in Radiation Technology from Iran University of Medical Sciences and MSc in Medical Physics from Tehran University of Medical Sciences (TUMS). She graduated on 2004 and continued her academic activity with PhD in Tehran University of Medical Sciences, Medical physics dep., and she is associate professor now.
Email: mhejazi@sina.tums.ac.ir
References   
1.  Boas D, Gaudette T, Arridge S. Simultaneous imaging and optode calibration with diffuse optical tomography. Opt Express 2001;8:26370. 
2.  Serdaroglu A, Yazici B, Kwon K, editors. Optimum Source Design for Detection of Heterogeneities in Diffuse Optical Imaging. Proc SPIE. Conference Paper, 2006, DOI: 10.1117/12.647209. 
3.  Stott JJ, Culver JP, Arridge SR, Boas DA. Optode positional calibration in diffuse optical tomography. Appl Opt 2003;42:315462. 
4.  Zhu Y, Jha AK, Wong DF, Rahmim A. Improved sparse reconstruction for fluorescence molecular tomography with poisson noise modeling. Microscopy Histopathology and Analytics. Optical Society of America; 2018. 
5.  Zhu Y, Jha AK, Wong DF, Rahmim A. Image reconstruction in fluorescence molecular tomography with sparsityinitialized maximumlikelihood expectation maximization. Biomed Opt Express 2018;9:310621. 
6.  Lasser T, Ntziachristos V. Optimization of 360 degrees projection fluorescence molecular tomography. Med Image Anal 2007;11:38999. 
7.  Pogue B, McBride T, Osterberg U, Paulsen K. Comparison of imaging geometries for diffuse optical tomography of tissue. Opt Express 1999;4:27086. 
8.  Graves EE, Culver JP, Ripoll J, Weissleder R, Ntziachristos V. Singularvalue analysis and optimization of experimental parameters in fluorescence molecular tomography. J Opt Soc Am A Opt Image Sci Vis 2004;21:23141. 
9.  Culver JP, Ntziachristos V, Holboke MJ, Yodh AG. Optimization of optode arrangements for diffuse optical tomography: A singularvalue analysis. Opt Lett 2001;26:7013. 
10.  Xing X, Lee K, Choe R, Yodh A, editors. Optimization of the Sourcedetector Geometry for Diffuse Optical Tomography. Biomedical Topical Meeting. Fort Lauderdale, Florida United States: Optical Society of America; 2006. 
11.  Zhang X, Badea C. Effects of sampling strategy on image quality in noncontact panoramic fluorescence diffuse optical tomography for small animal imaging. Opt Express 2009;17:512538. 
12.  Holt RW, Leblond F, Pogue BW, editors. Toward Ideal Imaging Geometry for Recovery Independence Fluorescence Molecular Tomography. Proc SPIE. Conference Paper, 2013, DOI:10.1117/12.2004870. 
13.  Kepshire DS, Davis SC, Dehghani H, Paulsen KD, Pogue BW. Subsurface diffuse optical tomography can localize absorber and fluorescent objects but recovered image sensitivity is nonlinear with depth. Appl Opt 2007;46:166978. 
14.  Correia T, Banga A, Everdell NL, Gibson AP, Hebden JC. A quantitative assessment of the depth sensitivity of an optical topography system using a solid dynamic tissuephantom. Phys Med Biol 2009;54:627786. 
15.  Dehghani H, White BR, Zeff BW, Tizzard A, Culver JP. Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography. Appl Opt 2009;48:D13743. 
16.  Holt RW, Leblond FL, Pogue BW. Methodology to optimize detector geometry in fluorescence tomography of tissue using the minimized curvature of the summed diffuse sensitivity projections. J Opt Soc Am A Opt Image Sci Vis 2013;30:16139. 
17.  Pera V, Brooks DH, Niedre M. On the use of the cramérrao lower bound for diffuse optical imaging system design. J Biomed Opt 2014;19:025002. 
18.  Sabir S, Kim K, Heo D, Cho S, editors. Adaptive Selection of Minimally Correlated Data for Optimization of SourceDetector Configuration in Diffuse Optical Tomography. Proc SPIE. Conference Paper, 2016, DOI: 10.1117/12.2211580. 
19.  Corlu A, Choe R, Durduran T, Rosen MA, Schweiger M, Arridge SR, et al. Threedimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans. Opt Express 2007;15:6696716. 
20.  Jermyn M, Ghadyani H, Mastanduno MA, Turner W, Davis SC, Dehghani H, et al. Fast segmentation and highquality threedimensional volume mesh creation from medical images for diffuse optical tomography. J Biomed Opt 2013;18:86007. 
21.  Dehghani H, Eames ME, Yalavarthy PK, Davis SC, Srinivasan S, Carpenter CM, et al. Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction. Commun Numer Methods Eng 2008;25:71132. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9]
