Ultrasound Image Based Human Gallbladder 3D Modelling along with Volume and Stress Level Assessment

Purpose Three-dimensional (3D) gallbladder (GB) geometrical models are essential to GB motor function evaluation and GB wall biomechanical property identification by employing finite element analysis (FEA) in GB disease diagnosis with ultrasound systems. Methods for establishing such 3D geometrical models based on static two-dimensional (2D) ultrasound images scanned along the long-axis/sagittal and short-axis/transverse cross-sections in routine GB disease diagnosis at the beginning of emptying phase have not been documented in the literature so far. Methods Based on two custom MATLAB codes composed, two images were segmented manually to secure two sets of the scattered points for the long- and short-axis GB cross-section edges; and the points were best fitted with a piecewise cubic spline function, and the short-axis cross-section edges were lofted along the long-axis to yield a 3D geometrical model, then GB volume of the model was figured out. The model was read into SolidWorks for real surface generation and involved in ABAQUS for FEA. Results 3D geometrical models of seven typical GB samples were established. Their GB volumes are with 15.5% and − 4.4% mean errors in comparison with those estimated with the ellipsoid model and sum-of-cylinders method but can be correlated to the latter very well. The maximum first principal in-plane stress in the 3D models is higher than in the ellipsoid model by a factor of 1.76. Conclusions A numerical method was put forward here to create 3D GB geometrical models and can be applied to GB disease diagnosis and GB shape analysis with principal component method potentially in the future.

A few meta-analyses on the existing outcome after LC based on CCK-cholescintigraphy or ultrasound scans for ABP patients were performed [24][25][26][27], and it was indicated that the use of EF to select patient with suspected ABP is not supported, even though the symptom is more significantly relived than nonoperative therapy; consequently, randomized clinical trials are on demand [24,25,28,29]. To standardize the trials and make the outcomes of different trials comparable, a protocol was recommended based on 0.02 µg/kg sincalide over 60 min with normal EF specified ≥ 38% [30]. A pilot randomized controlled trial of LC was launched in [31] based on that protocol for 30 ABP patients.

Ultrasound Images of Human GB
A series of two ultrasound static 2D images of human GB scanned in the long-axis/sagittal and short-axis/transverse directions has been obtained from a teaching hospital for seven patients. The images at the beginning of emptying phase are shown in Fig. 1 to reflect a variety of GB shapes. These images were taken when the patients were routinely inspected for GB diseases with an ordinary clinical ultrasound system. The hospital did not provide us with system parameter settings and patient information when the images were scanned. Such information seems to be unimportant here because the method for generating GB 3D geometrical models was concerned only in the paper.

Numerical Method
In the images shown in Fig. 1, there are two views only: one is in the long-axis plane where one longitudinal GB cross-section is presented; the other is in the short-axis/ transverse plane where one cross-section with the maximum cross-sectional area is specified. The two images in these two views must be utilized to generate 3D GB geometrical models. There are three steps in the method, as demonstrated in Fig. 2.
In the first step, GB ultrasound images are segmented manually in MATLAB with a custom program to decide scattered points of GB wall edges. Then, two edges are generated by fitting the points with the piecewise cubic spline function in MATLAB. The closed short-axis cross-section is lofted and the GB volume is calculated. Finally, the data files of short-axis cross-sections and guide lines are made for SolidWorks.
In the second step, GB short-axis cross-sections and guide line profiles are established in SolidWorks, and a surface is lofted with these profiles, the surface is enclosed with a small sphere cap in the apex of GB fundus. Geometrical model information exchange files in the initial graphics exchange specification (IGES) and Parasolid formats are prepared for ABAQUS or ADINA to perform FEA under a specific bile pressure and boundary condition.
In the third step, a geometrical model information exchange file is read into ABAQUS or ADINA. A mesh is created and a FEA is conducted by means of homogenous isotropic thin shell solid mechanics model with a certain thickness and a specific internal bile pressure load as well as a boundary condition. Figure 3 illustrates the scattered points of GB 36 wall edges, fitted short-axis and long-axis cross-sections, surface of the GB wall, GB surface presentation and a quadrilateral mesh for FEA in ABAQUS.
In the second step, GB cross-section curve fitting, lofting and volume calculation are the essential in 3D GB geometrical model generation. The details of the MATLAB code for accomplishing these functions are as follows: (1) Determine image scales based on marks and number figures in the images in Fig. 1; (2) Read data files of two edges scattered points are made in the previous segmentation; (3) Calculate the geometrical centre coordinates of the short-axis/transverse cross-sectional curve (x c , y c ), determine the polar-coordinate system (θ, r), and perform the curve fitting with the MATLAB spline function r = csape(θ, r,'periodic'), which is cubic spline interpolation with end conditions, the coordinates of the edge are (x c + r cosθ, y c + r sinθ); (4) Locate the middle point of GB neck (y fneck , z fneck ) and specify it as the reference point of the long-axis/sagittal cross-section edge scattered points, determine the other reference point ( y bmax , z bmax ) at which the maximum length is achieved from the middle point (y fneck , z fneck ), rotate the cross-section edge points in an angle of − tan −1 y bmax − y fneck ∕ z bmax − z fneck to allow the cross-section to be laid horizontally in the new coordinates ( y rot , z rot ), compute the arc length s of the scattered points from the first one to the last one in the neck, define the MATLAB spline function z = csape(s , z rot , q ), y = csape(s , z rot , q ), and q = 1 1 + z rot2 − z rot1 3 6 is the end condition parameter, z rot1 and z rot2 are the horizontal coordinates of the first and second scattered points after rotation, then the rotated scattered points are interpolated, the reference point ( y bmax , z bmax ) is searched again, and the angle of rotation is updated, and the points are rotated accordingly, the MATLAB spline function is redefined, and a subsequent interpolation is launched once more; such a cycle is repeated for five times until the sagittal cross-section is horizontal exactly; (5) Interpolate more dense scattered points with the functions defined in (4), divide the whole long-axis crosssection edge into two parts: i.e. upper and lower parts, determine the location where the maximum height in the long-axis cross-section occurs, rescale the existing  short-axis cross-section height and insert this crosssection into the long-axis cross-section orthogonally; (6) Specify the short-axis cross-section at the GB neck and the short-axis cross-section near the GB fundus apex are a circle, and then loft the rest short-axis crosssections from three existing short-axis cross-sections based on the ratio of a local short-axis cross-section height to the maximum short-axis cross-section height; (7) Calculate the volume of GB by summarizing the volumes of each truncated cone formed by two neighbouring short-axis cross-sections, this method is slightly more complex than the sum-of-cylinders method [51,52,58]; (8) Write the Cartesian coordinates ( x surf , y surf , z surf ) of all short-axis cross-section edges (lofted plus existing) into separate data files for SolidWorks, and the Cartesian coordinates of the long-axis cross-section edges in the horizontal and vertical planes are also written into data files for SolidWorks as guide lines.

GB Shape
GB wall surfaces of the rest six GB samples, which were generated by means of the method mentioned above, are illustrated in Fig. 4. It is shown that GB 1 and 17 present a pear shape, while GB 9, 30 and 36 in Fig. 3 exhibit a slender structure, but GB 19 and 37 are in an odd shape.

GB Bile Volume
The bile volumes of seven GB samples have been calculated based on the ellipsoid model with three major axis lengths measured in each 3D model dimensions from the image and the 3D model itself, respectively. The ellipsoid model [35] and its original 3D model for the image of GB 36 are illustrated in Fig. 5a, b as an example. The GB volume was estimated by means of the sum-of-cylinders method [51] and presented in the figure, too. The ellipsoid model volume is calculated by the expression V el = d 1 d 2 d 3 6 [52,58], the formula of the sum-of-cylinders method for GB volume can be found in [51,52,58]. The GB volume of 3D model is due to the method in Sect. 2.2. A comparison of GB volume is made between three methods in Fig. 5c and Table 1.
The GB volumes from 3D models, V 3D , are consistently larger than those from the simple ellipsoid model, V el , across all the GB samples with an error in the range of (0.8-29.2) % and the mean of 20.1%, which is basically close to a similar value of 15.0% in [58]. The GB volumes from the tedious sum-of-cylinders method, V sc , are basically smaller than those from 3D models with an error in (− 8.9 to 0.4)% and the mean of − 4.4%. Even though there is a noticed difference in GB volume across three methods, the GB volumes from 3D models can correlate to the volumes predicted with the sum-of-cylinders and ellipsoid models very well, respectively, as shown in Fig. 5d. This fact suggests that the method for generating 3D GB geometrical model from 2D ultrasound images proposed in the article is meaningful and reliable. Since GB volume from the ellipsoid model or sum-of-cylinders method is always higher or smaller than that from the 3D model, the ellipsoid method still can be used clinically due to its simplicity. Specially, one can use the correlations in Fig. 5d to predict a nearly true GB volume based on either V el or V sc .

FEA
Seven 3D GB geometrical models were put into ABAQUS 6.11 standard to carry out a linear FEA with an aim to identify whether the geometrical models can be applicable in FEA and what the stress pattern and level are different from those based on the analytical solution of an ellipsoid GB geometrical model in [35]. The material property constants and the model parameter setting in FEA are presented in Table 2.
The Young's modulus is 500 kPa approximately based on the experimental data summarized in [68], and the Poisson's ratio 0.49 is chosen for incompressible GB wall. The GB wall thickness was assumed to be 2.5 mm uniform [35]. Thin shell elements used in the FEA are included by choosing type S3 (3-node triangular general-purpose shell) and S4 (4-node quadrilateral general-purpose shell) elements. The pressure load was applied on the GB inside surface, and the load magnitudes are patient-specific and referred to [35]. The boundary condition was imposed in the GB neck inlet by fixing six degrees of freedom. The convergence tolerance is the default value of 1 × 10 −5 in displacements in ABAQUS. The independence of mesh size and effects of Young's modulus are demonstrated in Appendix. The results achieved with mesh1 in Table 5 are shown the following sections.
The ratio of the maximum first principal in-plane stress in a 3D GB geometrical model to the counterpart in its ellipsoid model, i.e. S = 1,3D 1,e , is shown in Fig. 6 to quantity the difference in stress level between two geometrical models. The first principal in-plane stress in the ellipsoid model was calculated based on the method in [35]. Clearly, under the same loading condition and material property constants, the stress level in the 3D geometrical model is higher than that in the ellipsoid model across the samples. The mean ratio is 1.76 with 0.11 standard deviation.  The contours of the first principal in-plane stress across seven 3D GB geometrical models are illustrated in Fig. 7. Further, the location for the maximum stress is coincident with the location for the largest strain, as shown in Fig. 8.
Since the curvature of the cross-sectional edge changes significantly in the short-axis/transverse plane, see Fig. 3b, the stress exhibits a bumpy variation, particularly with compression (negative) stress. In the contour of the first principal stress given by the analytical method based on the ellipsoid model in [35] demonstrates a bumpy characteristic as well, but in a different pattern and without compression stress, see Fig. 9 for GB 36.
These facts suggest that the GB ellipsoid model provides us with a simple clinical biomechanical model, however, its volume and stress level estimated can be lower by a certain factor than those based on 3D GB geometrical models from ultrasound images. Note that the contour of the first principal in-plane stress from the more realistic 3D GB models can differ from the pattern produced based on the ellipsoid model and analytical method.

Discussions
In the paper, a numerical method was proposed to establish human 3D GB geometrical models from static ultrasound 2D images in the long-axis/sagittal and short-axis/transverse cross-sections in routine biliary disease diagnosis at the beginning of emptying phase in hospital. The extracted 3D models were adopted to estimate GB volume in MAT-LAB, stress level and distribution in ABAQUS and the corresponding results were compared with those due to the existing ellipsoid model and analytical method [35]. The idea and work are original, and have not been indicated in the literature, and can be meaningful in biomedical engineering and ABP diagnosis.
The cubic spline function in MATLAB was chosen to fit the scattered points in the long-axis/sagittal and shortaxis/transverse cross-sections for the edges. The function can represent the complex edge shape of human GB wall. In [69], the Fourier series with five terms was utilized to fit the scattered points in the short-axis/transverse cross-section. Initially, this method was employed to perform the fitting task. However, the fitted curve is unable to pass through the scattered points picked up from the image simultaneously, as seen Fig. 10 for GB 36. However, the fitted curve generated by the cubic spline function passes through the points and can preserve the edge geometrical feature. Thus, the cubic spline function in MATLAB was favoured in the paper.
In Figs. 3 and 4, seven human GB samples exhibit a variety of shapes. Fortunately, these GB shapes can be observed in GB 3D volumes rendered computed tomography (CT) cholangiographic images shown in Fig. 11 [48]. Among these shapes, GB 1 has a similar shape to the GB in Fig. 11a to some degree, so does GB 17 to the GB in Fig. 11b, GB30 to the GB in Fig. 11c, GB 36 to the GB in Fig. 11d, GB 19 to the GB in Fig. 11e and GB 37 to the GB in Fig. 11f. This information suggests that the method proposed in the paper is sensitive to GB shape.
Honestly, the paper is subject to a few limitations. Firstly, the method was not validated directly by using a balloon Fig. 9 Contour of the first principal in-plane stress in GB36 due to ellipsoid model and analytical method in [35]  filled by liquid or 3D ultrasound data or magnetic resonance imaging (MRI) images of GB, because these kinds of data are lack in open source databases presently. However, an indirect validation has been made in Sect. 3.2 in terms of GB bile volume. The mean errors in GB bile volume calculated by the ellipsoid model and the sum-of-cylinders method are 20.1% and − 4.4% against the 3D GB geometrical model, respectively, in the paper. In [58], the corresponding mean errors of the ellipsoid model and the sum-of-cylinders method are 15.0% and 2.7% against the GB bile volume determined by 3D ultrasound system. Thus, the bile volume based on the 3D GB geometrical model in the paper seems to be reasonable and its accuracy in GB volume estimations can be considered at least equivalent to those of the ellipsoid model and the sum-of-cylinders method. Secondly, human GB wall can exhibit both anisotropic and nonlinear behavior in biomechanical property shown by in vivo computational biomechanics approach [36,70,71] and by in vitro uniaxial tensile test [72,73]. However, in the paper, GB wall property has been considered isotropic and linear, even though the geometrical nonlinearity was involved when the biomechanical model being set up in ABAQUS. During the emptying phase of human GB, there are both active and passive tensions. The passive tension and anisotropic property can be modelled by using the numerical approach and fibre structure presented in [71]. The method for modelling the active tension developed in human GB in the emptying phase has been unavailable at all so far.
Thirdly, in a linear numerical GB biomechanical analysis, the reference configuration can be the same as the current configuration. However, the reference configuration should be different from the current configuration because the strain can be quite large in a nonlinear numerical GB biomechanical analysis. Accordingly, a reference configuration must be known in advance in the nonlinear analysis. How to determine the reference configuration for a GB based on in vivo ultrasound images remains unknown currently. Nonetheless, active tension modelling and reference configuration determination will be two key issues in nonlinear human GB biomechanics during the emptying phase.
Based on GB volume predicted, the method proposed herein is equivalent to the ellipsoid model or sum-of-cylinders method in accuracy based on two-view ultrasound images. However, the method in the paper can produce a more realistic GB shape in terms of the comaprison with the existing GB volumes rendered from CT images. Consequently, the stress and strain levels and patterns in terms of such a shape may more likely close to the true values and patterns than those based on the existing GB ellipsoid Fig. 11 GB 3D volumes rendered computed tomography (CT) cholangiographic images at the beginning of emptying phase, a-g GB pictures are after [48] ▸ model. These 3D GB shapes potentially can be applied to GB shape analysis by using principal component analysis (PCA) method to correlate the shape features to GB pain score in the future. Note that the method is preliminary and further validations performed by radiologists or other experts are on demand.

Conclusions
A numerical approach was proposed to establish 3D GB geometrical models from existing static 2D long-axis/sagittal and short-axis/transverse cross-sectional ultrasound images scanned in routine GB disease diagnosis practice in hospital at the beginning of emptying phase. GB volumes of the geometrical models were extracted and compared with those estimated with well-known ellipsoid model and sumof-cylinders method. The models were read into SolidWorks for surface generation and then involved in ABAQUS for FEA analysis based on homogenous isotropic linear material and thin shell mechanical model. The results were discussed against those of the ellipsoid model. The approach developed in the article can create a satisfactory 3D geometrical model from the ultrasound images of human GB. The GB volumes from 3D geometrical models are different from those estimated by using the ellipsoid model and sumof-cylinders method with 15.5% and − 4.4% mean errors, respectively, but can correlate positively to the latter with a correlation coefficient as large as 0.99. The first principal stress level in the 3D geometrical models can be higher than that in the ellipsoid model by a factor of 1.0-2.63 with the mean of 1.76. The bumpy stress pattern in the 3D geometrical models is significantly deviated from that in the ellipsoid model. The method proposed potentially can be applicable in clinical diagnosis of GB disease and GB shape analysis in terms of principal component method in the future. The forthcoming work includes generating 3D GB geometrical models in the whole emptying phase, determining GB reference configuration, identifying anisotropic nonlinear material properties of GB wall in passive state, and modelling active behavior of GB wall in the emptying phase based on the generated 3D models.

Compliance with ethical standards
Conflict of interest All authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix Independence of Mesh Size and Effects of Young's Modulus
Quadrilateral and triangle elements were utilized to discretize GB walls in the paper. Table 3 illustrates the element shape selection criteria limits in ABAQUS for those two kinds of elements. For triangle elements the shape factor of an element is defined as the ratio of the element area to the optimal element area, which is the area of an equilateral triangle with the same circumradius as the element.
To identify the independence of mesh size in FEA of seven GBs, three meshes, namely mesh1, mesh2 and mesh3 were created with different mesh sizes, i.e. the edge lengths used in the discretization of element. The mesh quality parameters and number of elements in these meshes are summarized in Table 4. The mesh quality parameters include mean minimum, maximum and worst face corner angles, mean and worst aspect ratios and shape factor (for triangle elements only). In these meshes, each mesh quality parameter varies in a consistent range from mesh1 to mesh3 and is not beyond the limits specified in Table 3 at all, suggesting a better mesh quality.
Based on these meshes, the independence of mesh size was examined with mesh1, mesh2 and mesh3 in terms of the material property constants and loading conditions presented in Table 2. The maximum first principal stresses were extracted and are tabulated in Table 5. When the mesh is changed from mesh1 to mesh2 and mesh3, the errors in the stress vary in the range of − 3.05% and + 2.76%, depending on GB samples themselves. This fact suggests that the stress with mesh1 can be considered mesh size-independent. Thus, the stress and strain at mesh1 are used in the paper.
The stress level may be influenced by Young's modulus. Based on mesh1 and the loading conditions in Table 2, a series of FEA was performed when the Young's modulus varied − 20% (− 100 kPa) and + 20% (+ 100 kPa), respectively from the Young's modulus E = 500 kPa. It is shown in Table 6 that the stress level rises with increasing Young's modulus or declines with decreasing Young's modulus. A