Effect of subvoxel processes on non-destructive characterization of β-tricalcium phosphate bone graft substitutes

Acta Biomaterialia 7 (2011) 4045–4056

Contents lists available at ScienceDirect

Acta Biomaterialia journal homepage: www.elsevier.com/locate/actabiomat

Effect of subvoxel processes on non-destructive characterization of b-tricalcium phosphate bone graft substitutes M. Bashoor-Zadeh a, G. Baroud a,⇑, M. Bohner b a b

Laboratoire de Biomécanique, Département de Génie, Université de Sherbrooke, Sherbrooke, Canada QC J1K 2R1 RMS Foundation, 2544 Bettlach, Switzerland

a r t i c l e

i n f o

Article history: Received 27 January 2011 Received in revised form 27 May 2011 Accepted 7 July 2011 Available online 14 July 2011 Keywords: Subvoxelization Bone substitute Characterization Pore size Pore interconnection size

a b s t r a c t The geometric features of bone graft substitutes, such as the pore and pore interconnection sizes, are of paramount importance for their biological performance. Such features are generally characterized by micro-computed tomography (lCT). Unfortunately, the resolution of lCT is often too limited. The aim of this study was to look at the effect of lCT resolution on the geometric characterization of four different bone graft substitutes. An attempt was also made to improve the characterization of these materials by applying a subvoxelization algorithm. The results revealed that both approaches increased the accuracy of the geometric characterization. They also showed that the interconnection size in particular was affected. Comparing the results obtained from the scanned and numerical subvoxelization datasets revealed a minor difference of less than 2.5% for the porosity values. The difference for the pore sizes was up to 10%. Considerable differences of up to 35–50% were found for the interconnection sizes. The present study demonstrates how complex geometric characterization is and how important it is for biomaterial researchers to be aware of the impact of lCT resolution on the pore and pore interconnection sizes. Ó 2011 Published by Elsevier Ltd. on behalf of Acta Materialia Inc.

1. Introduction Geometric features of bone substitutes have a significant impact on the biological response of bone and on the healing process [1]. Generally, in porous bone substitutes the pores provide the space for bone formation and cell attachment, while pore interconnections act as a path for cell movement, nutrient transport and vascularization [2–4]. Clinical observations have shown that cells cannot penetrate narrow pores [2]. Hence, an effective bone substitute should have an interconnected porous structure with adequate pore interconnection sizes. Accurate characterization of the substitute microstructure has, therefore, received attention in designing effective scaffolds [1]. Characterization techniques reported in the literature include gas pycnometry, Archimedes method, mercury intrusion porosimetry, gas adsorption, optical microscopy and, more recently, microcomputed tomography (lCT) [1]. For example, Lin-Gibson et al. [5] compared different techniques (e.g. lCT, mercury porosimetry, and gravimetric analysis) to systematically study tissue engineering scaffolds and underlined the ability of lCT to evaluate the three-dimensional (3D) interior structure. Among different techniques they showed that lCT is the only technique allowing the

⇑ Corresponding author. Tel.: +1 819821 8000/61344; fax: +1 819 821 7163. E-mail address: [email protected] (G. Baroud).

determination of accurate qualitative and quantitative information on the internal structure of porous bone substitutes in three dimensions [6–9]. In addition, lCT is known to be a non-destructive method, as the samples remain intact and can be used for additional investigations [1,9]. The use of lCT has grown significantly for the non-invasive study of trabecular bone and porous bone substitutes [1,8,9]. The literature shows that the image resolution has a significant impact on the accuracy of microstructural analyses [10]. In particular, when the voxel size is relatively large compared with the structural features of interest, it becomes difficult to obtain accurate geometric information [11,12]. Therefore, there is increased interest in studying the effects of resolution and segmentation methods on the accuracy of microstructural analyses of micro-porous substitutes. In the following, research studies devoted to the relation between voxel size and geometric parameters are briefly reviewed. A number of studies have focused on the impact of voxel size on trabecular bone characterization using lCT and/or magnetic resonance (MR) images [11–17]. In some of these comparative studies the accuracy of the trabecular bone structural parameters obtained from MR images were compared with high resolution lCT images as a reference [11–13]. Specifically, Majumdar et al. [11] scanned trabecular bone specimens using a MR scanner at 156  156  300 lm voxel size and a high resolution X-ray topographic microscopy scanner at 18 lm isotropic voxel size. They showed that structural parameters such as trabecular thickness and volume fraction tended

1742-7061/$ - see front matter Ó 2011 Published by Elsevier Ltd. on behalf of Acta Materialia Inc. doi:10.1016/j.actbio.2011.07.014

4046

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

to be increasingly overestimated with lower resolution. In a study by Last et al. [12] the architectural parameters of calcaneus specimens were measured using MR images of 66 lm voxel size and high resolution lCT images of 10 lm voxel size. The results showed that the parameters quantified from the MR images were overestimated compared with the corresponding lCT data. Another approach to determine the impact of voxel size on trabecular bone characterization has been to acquire high resolution images by scanning the samples and then to subsequently produce low resolution images by artificial image degradation methods [14–17]. The artificial image degradation methods combined a group of neighboring voxels into one voxel with an intensity equal to the mean intensity of the voxel group [14–17]. For example, Tabor [14] took binary images at 30 lm resolution and artificially degraded them up to a pixel size of 150 lm. In a study by Kothari et al. [15] the original scanned images at 40 lm resolution were degraded in all directions to create low resolution images with isotropic and non-isotropic voxel sizes and then the dependency of bone structural parameters on slice thickness was investigated. Consistent with other investigations [11–13], this study showed that trabecular thickness increased with increasing slice thickness. Muller et al. [16] used trabecular bone images scanned at 14 lm and then artificially increase the voxel size by factors ranging from 2 to 20. Similarly to Kothari [15], they showed that a lower resolution led to an overestimation of trabecular thickness and trabecular separation. The above mentioned studies demonstrated that the resolution significantly affected the results for the structural parameters of bone. Specifically, when using lCT at higher resolution the partial volume effect (i.e. a voxel consists of a mixture of different tissue values) is decreased and more accurate structural and geometric information can be obtained. In addition to the use of fuzzy-based segmentation, two approaches can be used to improve the characterization of scaffold microstructures: (i) scanning the sample directly at higher resolution, leading to higher scanning costs; (ii) decreasing the voxel size artificially using subvoxelization algorithms, in which the high resolution images are subsequently generated by dividing voxels into subvoxels. There are several algorithms in the literature to artificially decrease the voxel size. For example, the linear interpolation method calculates the material volume fraction (MVF) of a voxel located between two adjacent voxels as the average MVF of the two voxels [10]. In addition, the Bayesian subvoxel approach classifies the subvoxels based on a Gibbs prior distribution. In this method the output is a binarized image, in other words the subvoxels are assigned as object or non-object [18]. Hwang and Wehrli [10] also developed a subvoxelization algorithm in which the intensity of each subvoxel is calculated based on local neighboring intensities with due attention to conservation principles. The output of this algorithm consists of fuzzy images because the subvoxel processing assigns partial fractions to each subvoxel [10]. They also showed that this method is superior to the interpolation method. This subvoxelization algorithm was developed by Hwang and Wehrli [10] for the study of trabecular bone using MR images. In this study the algorithm is adapted to study the structure of bone graft substitutes using lCT images. Previous studies have mainly investigated the impact of voxel size on the quantification of trabecular bone structure. To our knowledge no study has shown the relation between voxel size and the geometric parameters of bone substitutes, particularly the pore and pore interconnection sizes. Therefore, due to the significant effects of pore and pore interconnection size of a substitute on bone formation and cellular resorption, this article has focused on examining and improving the characterization of these pore-related parameters. In a previous study one resolution was used to perform an extensive characterization of bone graft substitutes with four different

pore sizes [19]. Here the aim is to use different resolutions and apply the subvoxelization algorithm of Hwang and Wehrli [10] to analyze their impact on the geometric features of these four bone graft substitutes. 2. Materials and methods 2.1. Scaffold fabrication and preparation Four classes of porous b-tricalcium phosphate (b-TCP) bone substitutes were fabricated by the calcium phosphate emulsion method [20,21]. Using this method, calcium phosphate cement pastes were produced by mixing the cement components and then oil droplets were dispersed in the resulting paste to form pores and interconnections. Variations in the emulsifier concentration led to different macropore sizes, while the volume fraction was constant. The resulting paste was then molded, after hardening. The samples were cleaned, dried and sintered at 1250 °C. The samples were then machined to produce cylinders 8 mm in diameter and 13 mm in height. One sample was randomly selected from each group of calcium phosphate bone substitutes with variable geometric features and analyzed at different resolutions to examine the effect of spacial resolution on the evaluation of microstructural parameters. 2.2. Image acquisition Four cylindrical samples of b-TCP bone substitutes with different geometric features were scanned with a SKYSCAN1172 lCT scanner (Desktop X-ray Microtomography, Aartselaar, Belgium). Each sample was scanned at three different resolution levels to generate the 7.5, 15 and 30 lm isotropic voxel sizes. For all configurations the X-ray tube voltage and current were fixed at 55 kV and 181 lA, respectively. A copper plus aluminum filter was used to absorb the soft X-rays and, consequently, reduce the beam hardening effects. A frame average of four and a rotation step of 0.4° were set within an angular range of 360°, thus a total of 900 radiographic images were acquired for each sample at different resolution levels. The total scanning time varied from approximately 30 to 100 min, depending on the resolution level. As the resolution increased the acquisition time increased accordingly. After acquiring the images, NRecon software v. 15.1.5 (Skyscan, Aartselaar, Belgium) was used to reconstruct three-dimensional (3D) structures from two-dimensional (2D) radiographic images and to generate 2D cross-sectional images. 2.3. Subvoxelization process The low resolution images (30 lm) were generated using the subvoxel processing algorithm established by Hwang and Wehrli [10]. The basic assumptions used to define this algorithm are: (a) subvoxels can be assigned to higher volume fractions; (b) continuity of the material phase [10,22]. According to this algorithm each reference voxel was divided into eight subvoxels and the gray value of each subvoxel was determined by adjacent neighboring values and strict conservation principles. In this algorithm the gray value allotted to a newly formed subvoxel is determined in two steps. First, a gray value for the index subvoxel is calculated based on the gray values of all voxels adjacent to the subvoxel. Second, the subvoxel gray value is adjusted relative to all newly formed subvoxels while taking into account the gray value of the original voxel in an attempt to conserve the initial gray value and the distribution of the scanned structure. Further details of the algorithm are given in Hwang and Wehrli [10]. Based on this process, images of 15 lm resolution were artificially generated and then analyzed to determine the pore and pore

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

interconnection sizes. The results obtained from artificial 15 lm images were compared with those obtained from 15 lm scanned images. 2.4. Thresholding The lCT images consist of an assembly of gray pixels. The intensity of the gray color is inversely proportional to the material radiodensity. To differentiate between material and pores two threshold values are defined. Voxels with a gray level value below the first threshold value are considered to belong to pores. Voxels with a gray level value above the second threshold value are considered to belong to the material. Voxels with intermediate gray level values consist of a mixture of material and pore (see also below). The threshold values were selected using the gray level histogram of each image. Their exact value was determined by visual inspection and by comparing the calculated porosity with the measured porosity of the bone substitute. Since the shape of the histograms changed between samples and resolution levels, determination of a general rule to obtain the threshold values proved to be difficult. When the image resolution was high enough the histograms became bimodal. The two peaks corresponded to the calcium phosphate material and the void space. The valley between the two peaks is known as the fuzzy zone, which corresponds to the intensities of voxels in the interface between the material and the pores. Threshold values around the minimum value of the fuzzy zone were selected to provide appropriate porosity values. The threshold values in unimodal histograms (obtained with too low a resolution) were selected by looking at the sphericity of the pores in the bone substitute images and by comparing their calculated and measured porosities.

4047

background voxel [24]. The fuzzy distance is computed based on both the membership values and the Euclidean distance between the voxels [19,25]. The FDT algorithm of Darabi et al. [25] was used for this purpose. Once the FDT map had been determined the next step was the determination of geometric parameters, such as the pore and pore interconnection sizes. Generally, the pore size was defined by the diameter of the maximum sphere that could be completely inscribed in a pore and the pore interconnection size was defined as the diameter of the maximum sphere in an opening between two pores [19] (Fig. 1). The pore size and pore interconnection size were then obtained by applying the max–min detection algorithm to the FDT map [19]. The algorithm finds the local maxima and saddle points (or voxels) of the FDT map. According to this algorithm, the maxima voxels have the largest FDT value among their neighbors and the saddle voxels have a maximum FDT value in one direction and a minimum value in the other direction. The FDT values at a local maxima and at a saddle correspond to the radius of a pore and a pore interconnection, respectively. The average pore size is then calculated using two definitions, the arithmetic average value (numberbased, NB) and the volume weighted average value (volume-based, VB). The average interconnection size is calculated by arithmetically averaging the FDT value of saddle voxels. When verifying the 2D images produced by high resolution scans it was noticed that decreasing the voxel size led to the detection of local irregularities in the pore walls (boundary effects). The local surface irregularities affected the FDT values and, conse-

2.5. Geometric analysis and verification The method used to determine the pore and pore interconnection sizes has been described extensively in an earlier article [19]. Briefly, as a first step the scanned images were fuzzified using a sigmoidal function with two threshold values for the calcium material and void space. Specifically, the two thresholds were selected from histograms and by visual inspection of the interface between the calcium material and the void space. Applying a sigmoidal function, each voxel x of set X was assigned a value between 0 and 1, termed the ‘‘membership value’’ l(x),

ðx; lðxÞÞjx 2 X; l : X ! ½0; 1:

ð1Þ

The membership value of a particular voxel determines the fraction of its volume belonging to the object. Therefore, the volume of the investigated object is obtained by summing the membership values of all voxels present in the volume of interest (VOI). Since in this study the object of interest was the void space, the summation of membership values corresponded with the volume of the void space. Hence, the porosity was calculated as the ratio of the volume of void space to the total number of voxels in the VOI [23]. The voxel membership of the object was also used to compute the fuzzy distance transform (FDT) map. For this purpose an algorithm was used to calculate the FDT map based on the fuzzy distance from the background. The fuzzy distance is first defined as the length of a path between two points on a fuzzy subset. Since there are infinitely many paths between two points the fuzzy distance is finally defined as the length of the shortest path between these two points. Hence, the FDT map is presented as a set of images in which the FDT value at each voxel is given as the shortest fuzzy distance between that voxel (i.e. the object voxel) and the

Fig. 1. (a) Fuzzified image of an individual void space. (b) 2D fuzzy distance transform map. The pore size and interconnection size are defined as the diameter of circles in 2D (or spheres in 3D) centered at local maxima (white points) and saddle voxels (black points), respectively. The FDT value is in pixel units (or voxel units in 3D).

4048

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

quently, resulted in finding additional saddle voxels as representative of new interconnections. The extra saddle voxels have a size in the range of local neighboring maxima (Fig. 2). Therefore, additional conditions, the so-called removal conditions, were applied to eliminate the effect of such saddle voxels and, consequently, eliminate boundary effects. Specifically, saddle voxels whose corresponding spheres overlapped a maxima sphere and whose overlap volume is larger than a specific fraction of their sphere volume (e.g. 70% of their sphere volume) were removed and the average interconnection size was recalculated. Overlap volumes of 60, 70, 80 and 90% were examined. 3. Results Four samples of different geometric features were analyzed. The samples were denominated with the letter A for the smallest pore size sample to D for the largest pore size sample. Fig. 3 shows representative 2D cross-sectional images of samples scanned at 30, 15 and 7.5 lm resolution. Visual comparison of the 2D images showed that on decreasing the voxel size pores became sharper and presented more spherical walls, particularly for samples with small pores (e.g. samples A and B). The gray level histograms of samples obtained from images scanned at different resolutions are shown in Fig. 4. The histograms

of samples A and B show that the scanning resolution had a considerable impact on distinguishing the two phases (pore and material). At low resolution the histograms were unimodal (Fig. 4, samples A and B at 30 lm resolution). As a result it was difficult to set two threshold values enabling a distinction between pore and material. At high resolution the histograms were bimodal, even for the smallest pore size (Fig. 4, sample A at 7.5 lm resolution). The subvoxelization process had little effect on the gray level histogram. In other words, the intensity histograms of 30 lm scanned images and 15 lm artificially generated images were similar. As a result, the same threshold values were used in both cases (Table 1). The threshold values and geometric parameters for substitutes of different voxel sizes and prepared by different methods are summarized in Table 1. Generally, for all samples the NB average pore size decreased as the voxel size decreased (Table 1). A relatively small difference was observed between the results calculated from the 15 lm scan and the results calculated from the subvoxelized 30 lm scan (<10%) (Table 2). Comparison with the VB average pore size values calculated at 15 lm resolution exhibited a large difference between the scan and artificial subvoxel process results (<18%) (Table 2). Fig. 5 shows the pore size distribution of all samples scanned at different resolutions, with and without subvoxelization. For the

Fig. 2. (a) Representative illustration of a 2D fuzzified image (left) and a fuzzy distance transform (FDT) map (right) of individual void spaces extracted from 30 lm resolution images. The local maxima are shown as white points in the FDT map. (b) Corresponding illustration of a 2D fuzzified image (left) and a FDT map (right) after applying subvoxel processing (15 lm resolution). Subvoxelization resulted in the appearance of surface irregularities (red circle), generating new local saddle voxels (representative of interconnections) in the FDT map (black point). The FDT value is in pixel units (or voxel units in 3D).

4049

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

Voxel size

Sample A

Sample B

Sample C

Sample D

30 µm

15 µm Artificial (subvoxel)

15 µm Scan

7.5 µm

Fig. 3. 2D representative lCT slices of samples A–D scanned at 30, 15 and 7.5 lm resolution.

two smallest pore sizes (samples A and B) the pore size distribution of the high resolution images demonstrated more detail, in that two peaks corresponding to different populations of the size inside the structure became apparent. The pore size distribution of sample A (7.5 lm) (Fig. 5) reveals two populations of pore size around 60 and 120 lm. In addition, the pore size distribution of sample B (Fig. 5b) displays two populations of pores with two equal peaks at 15 and 7.5 lm resolution. Decreasing the voxel size from 15 to 7.5 lm moved the first peak, representative of smaller pores, from 120 to 90 lm. Also, as demonstrated in Fig. 5c and d, the main pore size distribution peaks of samples C and D remained constant. The VB average pore size values accordingly show various trends as the voxel size decreased (Table 1). A change in voxel size also affected the average interconnection size. Decreasing the voxel size decreased the interconnection size in samples A and B, but increased the average interconnection size in sample D (Table 1). For sample C the interconnection size obtained from the scan and artificial datasets showed opposite trends, increasing in the case of the artificial dataset and decreasing in the case of the scan dataset. In addition, for the scan dataset the average interconnection size dropped to 47 lm at 15 lm resolution and then increased to 66 lm at 7.5 lm resolution (Table 1). Interestingly, the interconnection size distributions of all samples (Fig. 6) showed that most interconnections were well below 60 lm. Fig. 6 also shows that the main peaks shifted to the left with decreasing voxel size and that this shift decreased the average interconnection size. The interconnection size distribution of sample D (Fig. 6d) also shows that the main peaks moved to slightly

smaller values with subvoxelization (scan and artificial). This implies that the interconnection size also decreased. In contrast to this observation, Table 1 shows that the average interconnection values increased for sample D. Additional conditions, the so-called removal conditions, were applied to eliminate boundary effects appearing in high resolution images. Accordingly, saddle voxels, representative of interconnections, whose respective spheres overlapped a maximum sphere by more than X% (X = 60, 70, 80, and 90%) of their sphere volume were removed (Fig. 2) and were not considered in determining the average pore size. Then the average interconnection size was recalculated. Table 2 summarizes the corrected average interconnection size for X = 60%. Fig. 7 displays the modified average interconnection size versus various overlapping volumes. The results show that the average interconnection size of samples A and B were less influenced by boundary effect than samples C and D. For all overlapping volumes (60–90%) the average interconnection size of sample C decreased as the voxel size decreased (Fig. 7c). In addition, Fig. 7d shows that the average interconnection size of sample D (15 lm) (scan and artificial) decreased considerably after applying the removal conditions. Fig. 8 shows the evolution of the pore interconnection size distribution of sample D as a function of scanning resolution in more detail. The results show that the distributions of interconnections larger than 300 lm were the ones mainly affected by removing the overlapping spheres of the saddle voxels. The removal process was based on the removal conditions which were applied to exclude the saddle voxels whose spheres overlapped a maxima sphere by more than 60, 70, 80 or 90% of the saddle sphere volume.

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

15 µ m

5

x 10

100

x 10

1.0

7.5 µ m

6

40

0.5

20

0.0

0 200

40

0.2

20

0.0

0

250

0

50

Gray scale 1.8 1.5

60

0.9 0.6

40

0.3

20

0.0

0 0

50

100

150

200

x 10

1.0

40

0.5

20

0.0

0 200

Histogram counts .

60

150

0.6

20

0.0

0 50

x 10

60

1.5

40

1.0 0.5

20

0.0

0 100

150

200

200

250

200

250

6

x 10

100 Histogram Cumulative

80 60

4.5 40

3.0

20

1.5 0.0

0 0

50

100

150

200

250

1.2 x 10

80

1.0

1.0

40

0.5

20

0.0

0 150

200

7

100

60

100

150

Gray scale

Histogram Cumulative

50

100

6.0

250

1.5

5.0 Histogram counts .

Histogram counts .

Sample D

80

2.0

50

150

6

0

100 Histogram Cumulative

0

0 50

100 Histogram Cumulative

0.8

60 40

0.4

20

0.2 0.0

0 0

250

80

0.6

50

Gray scale

5

2.5

100

2.0

250

Cumulative (%)

x 10

7.5

0.3

Gray scale

3.0

9.0

40

2.5

Cumulative (%)

Histogram counts .

Sample C

80

1.5

100

100 80 60

0

100 Histogram Cumulative

50

20

1.0

Gray scale

5

0

40

2.0

Gray scale

0.9

250

x 10

2.0

60

3.0

0

Histogram Cumulative

1.2

80

4.0

250

6

Gray scale

2.5

200

Histogram counts .

100 80

Histogram counts .

Histogram counts .

Sample B

Histogram Cumulative

Cumulative (%)

5

1.2

150

5.0

Gray scale

x 10

1.5

100

100 Histogram Cumulative

0.0

Histogram counts .

150

0.4

Cumulative (%)

100

60

Cumulative (%)

50

0.6

6

6.0

100

150

200

250

Gray scale

6

x 10

100 Histogram Cumulative

4.0

80

3.0

60

2.0

40

1.0

20

Cumulative (%)

0

80

Histogram counts .

1.0

Histogram Cumulative

0.8

Cumulative (%)

60

Histogram counts .

1.5

80

Cumulative (%)

Histogram counts .

Sample A

Histogram Cumulative

x 10

7.0

100

2.0

Cumulative (%)

30 µ m

Cumulative (%)

Voxel size

Cumulative (%)

4050

0

0.0 0

50

Gray scale

100

150

200

250

Gray scale

Fig. 4. Intensity histograms of lCT images of samples A–D at various resolutions.

Table 1 Geometrical properties of bone substitutes presented as a function of resolution and method. Sample

Resolution

Method

Sample A

30 15 15 7.5 30 15 15 7.5 30 15 15 7.5 30 15 15 7.5

Scan Artificial Scan Scan Scan Artificial Scan Scan Scan Artificial Scan Scan Scan Artificial Scan Scan

Sample B

Sample C

Sample D

(subvoxel)

(subvoxel)

(subvoxel)

(subvoxel)

Threshold values

Porosity

Number-based pore size (lm)

Volume-based pore size (lm)

Interconnection size (lm)

Interconnection size (lm) after correction (X = 60%)a

55–80 55–80 40–90 40–80 50–100 50–100 45–85 45–80 55–100 55–100 50–90 40–75 50–110 50–110 30–80

50 52 54 55 52 53 54 54 52 54 54 55 53 53 54

135 ± 42 118 ± 37 107 ± 25 94 ± 34 193 ± 39 170 ± 45 181 ± 53 165 ± 64 365 ± 94 330 ± 108 353 ± 109 337 ± 125 909 ± 226 890 ± 236 899 ± 221

175 ± 52 152 ± 47 129 ± 46 138 ± 67 217 ± 44 204 ± 48 220 ± 48 223 ± 52 435 ± 119 428 ± 134 438 ± 122 439 ± 124 1030 ± 173 1019 ± 176 1015 ± 172

63 ± 29 54 ± 32 40 ± 19 22 ± 13 63 ± 28 52 ± 33 36 ± 23 30 ± 29 62 ± 41 66 ± 70 47 ± 53 66 ± 89 99 ± 106 234 ± 284 154 ± 217

51 ± 19 36 ± 19 37 ± 17 22 ± 13 59 ± 25 45 ± 27 33 ± 20 21 ± 19 58 ± 35 45 ± 40 35 ± 30 28 ± 33 85 ± 83 85 ± 100 66 ± 82

Due to the similarity of the average pore size and pore size distributions of sample D at 30 and 15 lm resolution this sample was not scanned and analyzed at 7.5 lm. a The average interconnection size after applying removal conditions which excluded the saddle voxels whose spheres overlapped a maximum sphere by more than 60% of their volume.

In other words, most of the large interconnections were generated by boundary effects or surface irregularities. The interconnection size distributions of sample D at 30 lm resolution (Fig. 8c) were not affected by applying the additional condition to remove overlapping spheres. Fig. 9 presents the differences between the pore size of each sample in terms of percentage error for the 7.5 and 15 lm scans.

Small pore size samples with higher voxel/pore size ratios revealed larger differences. 4. Discussion The objective of this study was to improve characterization of the pore and pore interconnection size of bone substitutes using lCT.

4051

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056 Table 2 A comparison of the voxel size dependency of structural parameters obtained from the scan and artificial datasets. Sample Sample A Sample B Sample C Sample D Difference

Method 15 micron-Scan 15 micron-Artificial (subvoxel) 15 micron-Scan 15 micron-Artificial (subvoxel) 15 micron-Scan 15 micron-Artificial (subvoxel) 15 micron-Scan 15 micron-Artificial (subvoxel) between actual scan results and artificial subvoxel process results

Porosity

Number-based pore size

Volume-based pore size

Interconnection size

N N N N N N N N <2.5%

. . . . . . . . <10%

. . N . N . . . <18%

. . . . . N N N 35–52%

Arrows (N, .) indicate the direction of change in the parameters as voxel sized shifts.

Fig. 5. Number-based pore size distributions of (a) sample A, (b) sample B, (c) sample C and (d) sample D derived from max–min operations and FDT values at different scan and artificial resolutions.

Pore and pore interconnection size are reported to have a significant impact on vascularization and new bone formation. In such porous materials pores provide a location for new bone cells and interconnections provide the path for cell transport [3,4]. In this study we focused on the accurate characterization of the bone substitute microstructural parameters. The aim was to improve characterization by increasing the voxel resolution using both numerical methods and high resolution scans. Previous studies on trabecular bone [10,17,26] have demonstrated that the structural parameters are highly influ-

enced by the voxel size. In particular, it was shown that higher resolution scans resulted in smaller and well-defined features. Since b-TCP bone substitutes have a porous structure, the structural parameters of the substitutes could also be affected by the resolution of the scanned images. To this end two strategies were adopted. In the first strategy the scanning resolution was improved. In the second strategy a subvoxelization algorithm [10] was applied to low resolution (30 lm) images of four types of bone substitutes to artificially generate 15 lm resolution images.

4052

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

Fig. 6. Interconnection size distributions of (a) sample A, (b) sample B, (c) sample C and (d) sample D derived from max–min operations and FDT values at different scan and artificial resolutions.

Analysis of the artificially generated high resolution images showed unexpected results with respect to the interconnection size, in that the average interconnection size increased for samples C and D compared with the original low resolution images. For example, the interconnection size of sample D increased from 99 ± 106 lm at 30 lm resolution to 154 ± 217 lm at 15 lm resolution (Table 1). To clarify this, the samples were scanned at higher resolutions (15 and 7.5 lm) and then the artificial and scanned datasets were compared. In this study the average pore size was calculated based on two definitions, VB and NB averages. Both VB and NB averages provide complementary information. The NB average is the arithmetic mean, whereas the VB average takes into consideration the volume of each pore and, thus, the smaller pores are less relevant when calculating the mean. Generally, due to the inability to detect smaller pores in low resolution images [17], the NB average pore size increased with increasing voxel size. A larger NB average pore size can also be related to a loss in the detection of smaller material parts, which results in a larger void space in the structure [17,26,27]. Furthermore, the NB average pore size showed a similar trend in relation to voxel size for all samples. However, the VB average pore size showed inconsistent trends among the samples (Table 1). This

unclear trend can be attributed to the definition of VB averaging, in which the pore size is weighted using the volume of the respective sphere in order to calculate the average. It is evident that a small pore makes a smaller contribution to calculating the VB average pore size. Therefore, depending on the size and number of small pores detected in high resolution images, the VB average pore size showed various trends when the voxel size was reduced. Furthermore, since higher resolution images detected finer features and since such features were mainly present at the pore interconnections [17], the interconnection size decreased with a higher scanning resolution, as observed for samples A and B (Table 1 and Fig. 6). Surprisingly, an opposite trend was observed for sample D (15 lm) and sample C (7.5 lm). Decreasing the voxel size helped in detecting new features or local irregularities on the pore walls, and these new features markedly influenced the calculation of interconnection size, in particular for samples C and D (Table 1). It is speculated that the increase in interconnection size of samples C and D is caused by local irregularities on the pore walls that appear in the high resolution images. It is clear that detecting more detail on the surface provides sharper and more irregular material surfaces [26]. These irregularities on the surface affect the FDT values and change the FDT map slightly.

4053

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

70

70

60 50 40 30 20 SampleA_30scan SampleA_15scan

10 0

60%

70%

SampleA_15artificial SampleA_7.5scan 80%

90%

Interconnection size (micro meter)

(b) 80

Interconnection size (micro metre)

(a) 80

60 50 40 30 20 10 0

all intercon.

60%

SampleB_30scan

SampleB_15artificial

SampleB_15scan

SampleB_7.5scan

70%

80 70

all intercon.

(d) 250

60 50 40 30 20 SampleC_30scan SampleC_15scan

10

SampleC_15artificial SampleC_7.5scan

Interconnection size (micro meter)

(c) Interconnection size (micro meter)

Overlapping volume (%)

80% 90% Overlapping volume (%)

200

150

100

50

SampleD_30scan

SampleD_15artificial

SampleD_15scan 0

0 60%

70%

80%

90%

all intercon.

Overlapping volume (%)

60%

70%

80%

90%

all intercon.

Overlapping volume (%)

Fig. 7. Average interconnection size versus overlapping volume for (a) sample A, (b) sample B, (c) sample C and (d) sample D. At each point interconnections with an overlapping volume larger than X% (X = 60–90%) were removed and the average interconnection size was calculated.

According to the max–min algorithm the saddle voxels were locally extracted based on the FDT values (Fig. 1). Changes in the FDT map affected local numerical detection and resulted in finding new saddle voxels that represented new interconnections. As shown in Fig. 2, the sizes of these new interconnections are in the size range of the adjacent pores, therefore leading to an increase in the interconnection size. To confirm this, removal conditions were applied to eliminate the boundary effects and the distributions were carefully analyzed. Under these added conditions saddle voxels, representative of interconnections, a considerable volume of which overlaps with the maximum sphere (Fig. 2), were removed. As illustrated in Fig. 7c and d, removing these saddle voxels, for example in samples C (7.5 lm) and D (15 lm, scan and artificial dataset), led to a considerable decrease in the average interconnection size, with nearly stable mean values being obtained (Fig. 7). These results confirm that the removal of overlapping saddle voxels significantly decreased the mean interconnection size and produced results that are in agreement with the literature [17,26]. The same trends were not observed for samples A and B (7.5 lm), probably because a much higher resolution would be required for pore boundaries to become less fuzzy such that these surface irregularities become apparent on the pore wall of these samples. The readers are advised that these effects do not influence the pore size characterization because of the unique methods applied to find interconnections. Analysis of interconnection size using high resolution images demonstrated that for samples A and B the average interconnection size dropped from 63 to 21 and 29 lm, respectively. According to existing paradigms a minimum interconnection size of 50 lm is generally required for bone in-growth [3,28,29]. Accordingly, our results for samples A and B may imply that these samples do not

allow cell migration inside the structure. However, the in vivo outcomes of implanting these samples in bone defects of sheep [30] revealed positive results in terms of good osteoconduction and bone formation in these samples. The discrepancy between the biological results and interconnection size can be interpreted by following mechanisms: first, the properties of the scaffold material and, second, the effects of micropores on bone formation. Due to the bioresorbability of the scaffold material (i.e. b-TCP) the sizes of the pores and interconnections are modified in the in vivo environment over time. This results in enlargement of the small interconnections and the creation of new interconnections between pores. It is also in agreement with the description of Lu et al. [3], in which they mentioned ‘‘in resorbable material the density of pore and interconnection is more important than their size because the sizes are enlarged by degradation’’. In addition, a number of recent studies have shown that micropores (typically defined as pore sizes less than 10 lm) in CaP-based scaffolds can support bone in-growth [31–33]. Levengood et al. [31] observed the presence of bone, osteoid and osteogenic cells in micropores of biphasic calcium phosphate scaffolds. They showed that the migration of osteo-progenitor cells inside micropores and subsequent bone formation led to the in-growth of bone throughout CaP-based scaffolds [31]. Consequently, the biological results for samples A and B can be related to a combination of the above mechanisms, including bioresorbability of the CaP components of the scaffold material and the development of bone in-growth into micro-sized porous space. There is a third potential mechanism to explain these biological results. Specifically, the interconnection size distributions of samples A and B (Fig. 6a and b) show that 5% of all interconnections were >50 lm, and this could be sufficient to provide an

4054

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

60

3.5

Histogram count (%)

50

3

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

40 30

Histogram count (%)

(a)

20 10

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

2.5 2 1.5 1 0.5

0 0

200

400

600

800

1000

1200

0 300 400 500 600 700 800 900 1000 1100 1200 1300 1400

1400

Interconnection size(micro meter)

Interconnection size (micro meter)

(b) 60

3.5 3

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

40 30

Histogram count (%)

Histogram count (%)

50

20

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

2.5 2 1.5 1 0.5

10

0 300 400 500 600 700 800 900 1000 1100 1200 1300 1400

0 0

200

400

600

800

1000

1200

1400

Interconnection size (micro meter)

Interconnection size (micro meter) 60

3.5

Histogram count (%)

50

3

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

40 30

Histogram count (%)

(c)

20

OV = 60% OV = 70% OV = 80% OV = 90% all interconnections

2.5 2 1.5 1

10 0.5

0 0

200

400

600

800

1000

1200

1400

0 200

300

400

500

600

700

Interconnection size (micro meter) Interconnection size (micro meter) Fig. 8. The impact of boundary effects on the interconnection size distributions of (a) sample D (15 lm), (b) sample D (15 lm) artificial and (c) sample D (30 lm) scans. Removing the overlapping interconnections significantly affected the distribution of interconnections larger than 300 lm. OV, overlapping interconnection.

interconnected network for bone cells to migrate into the structure of these resorbable bone substitutes. Since the pore size and interconnection size of sample D at 30 and 15 lm resolution exhibited similar distributions and the average pore size and porosity were almost the same, we did not analyze sample D at 7.5 lm. In summary, an attempt was made to assess the potential of a subvoxelization algorithm to more accurately describe the pore architecture of a bone substitute. A comparison between parameters calculated from scanned and artificial datasets at 15 lm resolution indicates that the change in geometric features due to an increase in resolution (through scanning or through subvoxelization) follows similar trends. However, the values obtained by each method (scan and artificial) were not identical. The presented error between the scan and artificial results (Table 2) can be ascribed to

the inability of the subvoxelization process to determine details which are undetectable in the original low resolution images.

5. Conclusion Material scientists increasingly use the lCT technique to characterize the pore architecture of bone graft substitutes with the aim of better understanding the complex process of bone repair. Indeed, existing paradigms interrelate the process of bone repair with the porous structure of bone graft substitutes. Scientists therefore seek the most accurate approaches to characterizing such structures. Considering the difficulties of scanning large samples at high resolution, it is fortunate that the scanning accuracy can be improved by using a subvoxelization approach, as shown

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

4055

50

Error relative to high resolution (%)

E1(30-15) 40

E2(15-7.5) E3(30-7.5)

30

20

10

0 SampleA(0.22)

SampleB(0.16)

SampleC(0.08)

SampleD(0.03)

Sample (voxel size/Pore size) Fig. 9. Percentage errors relative to high resolution datasets for four different bone substitute samples. Voxel size/pore size ratios calculated based on a voxel size of 30 lm and the corresponding computed pore size. E1 (30–15), the percentage error between the pore size analysis at 30 and 15 lm; E2 (15–7.5), the percentage error between the pore size analysis at 15 and 7.5 lm; E3 (30–7.5), the percentage error between the pore size analysis at 30 and 7.5 lm.

here. For that purpose subvoxelization was applied to four classes of b-TCP bone substitutes of constant porosity but of four different pore diameters ranging from 0.1 to 1.0 mm. To evaluate the numerical subvoxelization higher resolution scans were also obtained. The results demonstrate that the changes in porosity, pore size, and pore interconnection size calculated by the subvoxelization approach in most cases follow those of higher resolution scans. Also, using both subvoxelization and higher scan resolutions revealed distinguished pore distribution features that were not previously apparent. However, the extent of these changes was quite different due to undesired boundary effects. In the present study a method has been proposed to reduce these boundary effects.

Appendix A. Figures with essential colour discrimination Certain figures in this article, particularly Figures 1 and 2, are difficult to interpret in black and white. The full colour images can be found in the on-line version, at doi:10.1016/j.actbio.2011.07.014. References [1] Ho ST, Hutmacher DW. A comparison of micro CT with other techniques used in the characterization of scaffolds. Biomaterials 2006;27(8):1362–76. [2] Jones AC, Arns CH, Hutmacher DW, Milthorpe BK. The correlation of pore morphology, interconnectivity and physical properties of 3D ceramic scaffolds with bone ingrowth. Biomaterials 2009;30(7):1440–51. [3] Lu JX, Flautre B, Anselme K, Hardouin P, Gallur A, Descamps M, et al. Role of interconnections in porous bioceramics on bone recolonization in vitro and in vivo. J Mater Sci Mater Med 1999;10(2):111–20. [4] Mastrogiacomo M, Scaglione S, Martinetti R, Dolcini L, Beltrame F, Cancedda R, et al. Role of scaffold internal structure on in vivo bone formation in macroporous calcium phosphate bioceramics. Biomaterials 2006;27(17): 3230–7. [5] Lin-Gibson S, Cooper JA, Landis FA, Cicerone MT. Systematic investigation of porogen size and content on scaffold morphometric parameters and properties. Biomacromolecules 2007;8(5):1511–8. [6] Moore MJ, Jabbari E, Ritman EL, Lu L, Currier BL, Windebank AJ, et al. Quantitative analysis of interconnectivity of porous biodegradable scaffolds with micro-computed tomography. J Biomed Mater Res 2004;71(2):258–67. [7] Dorsey SM, Lin-Gibson S, Simon CG. X-ray microcomputed tomography for the measurement of cell adhesion and proliferation in polymer scaffolds. Biomaterials 2009;30(16):2967–74. [8] Ruegsegger P, Kollar B, Muller R. A microtomographic system for the nondestructive evaluation of bone architecture. Calcif Tiss Int 1996;58:24–5. [9] Van Lenthe GH, Hagenmueller H, Bohner M, Hollister SJ, Meinel L, Mueller R. Nondestructive micro-computed tomography for biological imaging and quantification of scaffold-bone interaction in vivo. Biomaterials 2007;28(15): 2479–90.

[10] Hwang SN, Wehrli FW. Subvoxel processing: a method for reducing partial volume blurring with application to in vivo MR images of trabecular bone. Magn Reson Med 2002;47(5):948–57. [11] Majumdar S, Newitt D, Mathur A, Osman D, Gies A, Chiu E, et al. Magnetic resonance imaging of trabecular bone structure in the distal radius: relationship with X-ray tomographic microscopy and biomechanics. Osteoporos Int 1996;6(5):376–85. [12] Last D, Peyrin F, Guillot G. Accuracy of 3D MR microscopy for trabecular bone assessment: a comparative study on calcaneus samples using 3D synchrotron radiation microtomography. Magma 2005;18(1):26–34. [13] Pothuaud L, Laib A, Levitz P, Benhamou CL, Majumdar S. Three-dimensionalline skeleton graph analysis of high-resolution magnetic resonance images: a validation study from 34-micron-resolution microcomputed tomography. J Bone Miner Res 2002;17(10):1883–95. [14] Tabor Z. Analysis of the influence of image resolution on the discriminating power of trabecular bone architectural parameters. Bone 2004;34(1):170–9. [15] Kothari M, Keaveny TM, Lin JC, Newitt DC, Genant HK, Majumdar S. Impact of spatial resolution on the prediction of trabecular architecture parameters. Bone 1998;22(5):437–43. [16] Muller R, Koller B, Hildebrand T, Laib A, Gianolini S, Ruegsegger P. Resolution dependency of microstructural properties of cancellous bone based on threedimensional mu-tomography. Technol Health Care 1996;4(1):113–9. [17] Cooper D, Turinsky A, Sensen Ch, Hallgrimsson B. Effect of voxel size on 3D micro-CT analysis of cortical bone porosity. Calcified Tissue Inc 2007;80(3):211–9. [18] Wu Z, Chung H-W, Wehrli FW. A Bayesian approach to subvoxel tissue classification in NMR microscopic images of trabecular bone. Magn Reson Med 1994;31(3):302–8. [19] Bashoor-Zadeh M, Baroud G, Bohner M. Geometric analysis of porous bone substitutes using micro-computed tomography and fuzzy distance transform. Acta Biomater 2010;6(3):864–75. [20] Bohner M, van Lenthe GH, Gruenenfelder S, Hirsiger W, Evison R, Mueller R. Synthesis and characterization of porous b-tricalcium phosphate blocks. Biomaterials 2005;26(31):6099–105. [21] Bohner M. Calcium phosphate emulsions: possible applications. Key Eng Mater 2001;192–195:765–8. [22] Wehrli FW, Saha PK, Gomberg BR, Song HK. Noninvasive assessment of bone architecture by magnetic resonance micro-imaging-based virtual bone biopsy. Proc IEEE 2003;91(10):1520–42. [23] Sladoje N, Nystrom I, Saha PK. Measurements of digitized objects with fuzzy borders in 2D and 3D. Image Vis Comput 2005;23(2):123–32. [24] Saha PK, Wehrli FW, Gomberg BR. Fuzzy distance transform: theory, algorithms and applications. Comput Vis Image Understanding 2002;86(3):171–90. [25] Darabi A, Chandelier F, Baroud G. Thickness analysis and reconstruction of trabecular bone and bone substitute microstructure based on fuzzy distance map using both ridge and thinning skeletonization. Can J Elect Compt Eng 2009;34(1/2):57–62. [26] Peyrin F, Salome M, Cloetens P, Laval-Jeantet AM, Ritman E, Ruegsegger P. Micro-CT examinations of trabecular bone samples at different resolutions: 14, 7 and 2 micron level. Technol Health Care 1998;6(5/6):391–401. [27] Kim DG, Christopherson GT, Dong XN, Fyhrie DP, Yeni YN. The effect of microcomputed tomography scanning and reconstruction voxel size on the accuracy of stereological measurements in human cancellous bone. Bone 2004;35(6):1375–82.

4056

M. Bashoor-Zadeh et al. / Acta Biomaterialia 7 (2011) 4045–4056

[28] Chang BS, Lee C-K, Hong KS, Youn HJ, Ryu HS, Chung SS, et al. Osteoconduction at porous hydroxyapatite with various pore configurations. Biomaterials 2000;21:1291–8. [29] Otsuki B, Takemoto M, Fujibayashi S, Neo M, Kokubo T, Nakamura T. Pore throat size and connectivity determine bone and tissue ingrowth into porous implants: three-dimensional micro-CT based structural analyses of porous bioactive titanium implants. Biomaterials 2006;27(35):5892–900. [30] von Doernberg MC, von Rechenberg B, Bohner M, Gruenenfelder S, van Lenthe GH, Mueller R, et al. In vivo behavior of calcium phosphate scaffolds with four different pore sizes. Biomaterials 2006;27(30):5186–98.

[31] Lan Levengood SK, Polak SJ, Wheeler MB, Maki AJ, Clark SG, Jamison RD, et al. Multiscale osteointegration as a new paradigm for the design of calcium phosphate scaffolds for bone regeneration. Biomaterials 2010;31(13):3552–63. [32] Habibovic P, Sees TM, van den Doel MA, van Blitterswijk CA, de Groot K. Osteoinduction by biomaterials–physicochemical and structural influences. J Biomed Mater Res A 2006;77(4):747–62. [33] Hing KA, Annaz B, Saeed S, Revell PA, Buckland T. Microporosity enhances bioactivity of synthetic bone graft substitutes. J Mater Sci Mater Med 2005;16(5):467–75.