Abstract

Bone tissue engineering has emerged as a promising strategy for the treatment of osseous fractures, defects, and ultimately diseases caused by, for example, bone tumor resection, accident trauma, and congenital malformation. Additive fabrication of stem cell-seeded, osteoconductive porous scaffolds has been an effective method in clinical practice for the treatment of bone pathologies (such as osteoporosis, osteoarthritis, and rheumatic diseases). Porosity is known to be one of the main morphological characteristics of bone tissues, which affects the functional performance of an implanted bone scaffold. Hence, in situ detection and quantification of scaffold porosity implemented to ensure functional integrity prior to implantation/surgery is an unavoidable need. The objective of this research work is to introduce a robust, image-based method for identification and subsequently characterization of the surface porosity and dimensional accuracy of additively manufactured bone tissue scaffolds, with a focus on pneumatic micro-extrusion (PME) process. It was observed that the presented method would be capable of detecting complex individual pores based on a micrograph. Using the proposed method, not only were scaffold pores detected, but also scaffold porosity was characterized on the basis of various defined quality metrics/traits (such as the relative standard deviation of distance to the nearest pore). The proposed method was validated by contrasting its performance in “surface pore detection” against that of a standard method, tested on a complex benchmark in four different simulated lighting environments. Besides, the performance of the method in terms of “pore filling” was compared to that of a standard method, tested on a real PME-fabricated bone scaffold. It was observed that the proposed method had a better performance in pore filling, detection, and consolidation. Overall, the outcomes of this work pave the way for high-resolution fabrication of patient-specific, structurally complex, and porous bone scaffolds with easily validatable, functional, and medical properties for the treatment of bone pathologies.

1 Introduction

1.1 Background.

Osseous fractures and defects are significant contributors to disability in the world. Millions of patients suffer yearly from bone fractures and disorders, which are the most common causes of severe long-term pain as well as physical disabilities [1,2]. The use of advanced manufacturing methods—such as pneumatic micro-extrusion (PME), selective laser sintering, and particulate leaching [35]—for the fabrication of biocompatible scaffolds (capable of hosting bone morrow-derived stem cells) has emerged as a clinically effective as well as a cost-effective approach. Most of the current treatment methods in clinical practice are centered on the implantation of bone grafts and/or osteoconductive porous scaffolds, seeded with autologous stem cells for bone regeneration. Patient-specific treatment of bone fractures requires fabrication of mechanically robust, dimensionally accurate, and porous bone tissue scaffolds. Therefore, there is a need for characterization of scaffold porosity, morphology, and functional integrity during and after scaffold fabrication process.

The functional, mechanical, and medical properties of bone scaffolds are significantly influenced by porosity (as an intrinsic scaffold property). For example, porosity affects the rate of osteogenesis and ultimately fracture healing. It is expected that bone scaffolds are composed of interconnected pores, enhancing material transport. Too small pores prevent cell migration, limiting the diffusion of nutrients as well as removal of waste products [6]. In addition, closed pores obstruct material transport (e.g., growth factors), cell proliferation, and osteoconduction inside a scaffold. Too large pores, while improving cellular infiltration, lead to a decrease in the surface area needed for cell adhesion and compromise the mechanical properties of scaffolds as a result of void volume [7]. Also, sub-optimal pores adversely affect not only scaffold functional integrity but also osteoconduction as well as osseointegration. It has been shown that achieving mean pore sizes of 96–150 µm would facilitate optimal cell attachment, while large pores in the range of 300–800 µm would be required for successful bone growth in scaffolds [6]. It has been observed that optimal pore size depends on different cell types [8].

Detection and quantification of scaffold porosity aid in (i) a better understanding of cell activity/behavior and (ii) identification and prediction of the underlying phenomena influential on the functional performance (including failure), accuracy, and functional integrity of fabricated bone tissue scaffolds. Besides, if implemented in real-time during the fabrication process, assessment of scaffold porosity allows for property control and thus achieving bone scaffolds with optimal functional and medical properties.

The research problem, which is tackled in this work, is based on the fact that the characterization of bone scaffold porosity and morphology is a complex image-processing problem since bone scaffolds (in an acquired image) are inherently composed of pores that are not discretely and completely bound (known as gradual pores). None of the current image processing methods can effectively address this problem without extensive use of, for example, X-ray computed tomography (CT) imaging, which cannot be implemented in situ nor practically in a common environment. Consequently, development of a robust image-based method that allows for in situ characterization of scaffold morphology with gradual pores would be inevitable. In the absence of such a method, optimal fabrication as well as in situ characterization of bone scaffolds will remain trial-and-error-driven (not patient-specific) with no control over desired medical and functional properties.

1.2 Review of Literature.

There is a broad range of image-based methods, developed for the characterization of additively manufactured parts. In a research work by Baldwin et al. [9], a new method was developed for the identification of pores and partitioning the void space into a well-defined collection of individual pores from volume images. It was observed that the proposed method would be able to process three-dimensional (3D) images as well as reconstructed porous media. However, this method is not optimal for the processing of two-dimensional (2D) images (in contrast to volume images) and for real-time monitoring and characterization of manufacturing processes through simple digital images.

Similarly, e-Silva et al. [10] introduced a new image-based method with the aim to measure pore size distribution and to estimate distribution curves. Utilizing an automatic algorithm, the method was established based on the concept of the largest inscribing opening size. In addition, the method was presented as a software program including a learning algorithm (based on Support Vector Machine) and digital image processing, developed to identify fibers and pores. Working based on expected amplitude distribution (instead of edge detection), this method is similar to the method presented in this study. However, it may be sub-optimal in performance for discriminating pores as well as managing pores without discrete, fully enclosed boundaries. In addition, it requires very specific, even light conditions over an image in order to gain accurate results, as the amplitude distribution method may not be as robust as an edge detection analysis on a pore-by-pore basis. Therefore, it would not be effective enough for the complex non-defined geometry intrinsically present in PME-fabricated bone scaffolds.

Chukwuma et al. [11] developed an image-based method for the identification and characterization of pore systems, established on the basis of ultrathin sections as well as field-emission scanning electron microscopy. Capable of processing a large number of images, the proposed method facilitated the characterization and quantification of the porosity of rock formation. Similarly, in a research work by Chen et al. [12], a pore quantification method was presented—working based on scanning electron microscopy (SEM) and a thresholding algorithm—with the aim to identify organic pore characteristics, e.g., pore size distribution curves of marine shale. It was observed that the resolution of the SEM technique might be a hindering factor, adversely affecting the characterization of organic pores. This technique is not viable though for in situ manufacturing or analyzing the internals of a fabricated scaffold without damaging it in the process.

X-Ray computed tomography has been utilized as a non-destructive method for the detection of a wide range of defects. Based on X-ray computed tomography and three-dimensional (3D) image analysis, Kim et al. [13] presented a method for investigation of the pore structure of additively manufactured metal parts. It was observed that the accuracy of pore structure quantification depended, to a great extent, on the segmentation of X-ray computed tomography images. The proposed method allowed for quantification of not only porosity, but also pore size distribution, pore shape, and particle size distribution. This method uses a thresholding algorithm developed to filter out noise. However, the thresholding method is not viable for use with pores lacking discrete, complete boundaries, nor is it practical to readily use X-ray analysis for most lab environments, in comparison to digital images (acquired using, for example, charge-coupled device cameras). Lastly, it has no potential for real-time analysis of the quality of material deposition and ultimately fabricated scaffolds.

Du Plessis et al. [14] developed a similar method for additive manufacturing quality control with a focus on porosity analysis, based on X-ray computed tomography. The goal of this work was to obtain porosity information in a standard fashion even from samples fabricated using different additive manufacturing systems. It was demonstrated that the method was capable of quantifying complex pores, such as tree-like pores (that grow in the direction of laser fusion) as well as sub-surface pores (that form due to the lack of fusion). Similarly, Léonard et al. [15] utilized X-ray computed tomography to characterize pore size, volume fraction, and spatial distribution, which were subsequently correlated to selective electron beam melting process parameters. Besides, Villarraga et al. [16] demonstrated dimensional metrology (namely, part-to-model comparison and wall thickness analysis) as well as inspection of material quality (e.g., void detection) on the basis of X-ray computed tomography. It was observed that void detection would be inherently limited by the resolution of the CT scanning system and therefore pores with a diameter of approximately one voxel or less could not be detected in a CT image. All in all, these methods are relatively expensive, difficult to validate, limited in penetration depth [17], and not viable for real-time morphology analysis during PME fabrication process.

Aminzadeh et al. [18] proposed an image-based method for in situ monitoring of porosity in metal powder-bed additive manufacturing. The method was centered on producing images in situ from each layer that would allow for visualization of surface quality as well as formed porosity. In addition, an intelligent (statistical Bayesian-based) framework was developed for identification of defective regions and ultimately qualitative assessment of porosity. However, the formation of microstructures and defects in the PME process is not comparable and is, to some extent, more complex (due to, for example, non-Newtonian molten polymer flow deposition and solidification) than microstructure/defect formation in powder-bed fusion. Therefore, a more sophisticated technique is needed for collecting accurate porosity information for PME-fabricated bone scaffolds.

Infrared-assisted methods have also been utilized for manufacturing process monitoring. Hu and Kovacevic [19] used infrared image sensing for real-time sensing and control of powder bed fusion (PBF) additive manufacturing of functionally graded materials. The proposed method was capable of sensing powder delivery rate in real-time with a high sampling frequency. Geometrical accuracy was achieved by controlling heat input as well as molten pool size as two consequential process parameters. Besides, the authors developed a finite element model to investigate the thermal behavior (including temperature distribution, geometrical feature, and cooling rate) of the molten pool in a building single-bead wall. The finite element model, in addition, paves the way for further residual stress analysis and provides guidance for the selection of optimal process parameters. In a similar work by Krauss et al. [20], thermography (using an infrared camera) was utilized for part quality assurance in the PBF process via monitoring the temperature distribution of not only a complete layer but also its temporal evolution. The authors developed a method for detecting pores based on the temperature distribution at each layer. The proposed method also allowed for observing process parameter deviations (based on the total area of the heat-affected zone in the PBF process) and detecting deviations and drifts (quantified based on circularity and aspect ratio). It was observed that adverse physical phenomena, such as insufficient heat dissipation, result in sub-optimal layer formation. Craeghs et al. [21] similarly developed a data processing algorithm for the detection of (i) process failures and (ii) deformation due to thermal stresses as well as overheating at overhang zones in the PBF process via real-time melt pool monitoring. The developed algorithm allowed for mapping the measured melt pool data in space (by simultaneously sampling and logging of the position and melt pool data). The proposed method also led to a reduction in data interpretation. Overall, while thermography is an effective method for layer-wise pore detection, real-time process monitoring, and understanding underlying physical phenomena, it (i) requires complex laser setup, (ii) is most suitable for PBF-like processes, where there is significant radiation given off by, for example, melt pool, and (iii) does not provide quantitative Z-axis data [17].

Holzmond et al. [17] introduced a novel non-destructive method for in situ real-time detection of localized and global defects (such as blob of filament as well as low flow) of additively manufactured parts, based on the certify-as-you-build concept. The proposed method not only captures the geometry (and thus allows for monitoring the surface geometry) of a part during the fabrication process (using 3D digital image correlation) but also compares the captured geometry with the computer-aided design (CAD) model to detect deviation or error in real-time. The proposed method requires a speckle texture pattern and optimal lighting conditions.

Bowoto et al. [22] developed a method for in situ layer-wise defect detection and examination of porosity in fused filament fabrication. The calculation of porosity was based on a thresholding process. It was demonstrated that the proposed method was able to detect the presence or absence of defects. However, a simple thresholding method is not sophisticated enough to be able to consistently and accurately identify gradually developing indentations in varying lighting conditions and varying shapes and sizes. This will be demonstrated in the study.

1.3 Significance.

The porosity of PME-fabricated bone scaffolds, captured by a micrograph, is intrinsically composed of a wide range of complex gradual pores characterized by intensity gradients and no discrete, sharp edges. In fact, it is not easy to identify the regions near an edge and thus to define regions corresponding to a pore from a single micrograph. Hence, “topographical or depth mapping” (proposed in this work) is a critical step toward identification of gradual pores and quantification of scaffold surface porosity. In the absence of such a topographical mapping-driven method, the detected pores of a bone scaffold may significantly overlap or be internally discontinuous/heterogeneous, resulting in inaccurate quantification of scaffold porosity (in terms of, for example, estimation of surface area per volume, pore size distribution, and material diffusion rate) as well as subsequent cell–scaffold interactions.

Topographical maps of images with intensity (represented as altitude) will be presented in this work as a critical component of the porosity filling algorithm (discussed in Sec. 2.4.4). This concept was originally presented in a 1990’s work [23] and was referred to as the “watershed method”. Delineated in Sec. 2.4, this concept is used extensively for analyzing the porosity of complex PME-fabricated scaffolds, as a key principle for identifying porosity in a wide range of lighting conditions, and on a wide range of pore geometries (most specifically on pores, which are gradually forming without a discretely defined boundary). This will be a key principle, contributing to the novelty of the proposed algorithm.

Another key principle of this work, which sets it as an improvement above the standard methodologies (presented in Sec. 1.2), is the robustness of the algorithm to different, complex, real-world lighting conditions without calibration. An image processing technique, which is often used to solve this problem, is Otsu’s method from a 1979 work [24]. This was used by Bowoto et al. [22]; however, it only works when there are clearly defined dark and light regions with a clear partition in intensities. For complex, gradually forming indentations in variable lighting conditions, there are often more than discrete groupings of intensities, which cause Otsu’s method to fail.

Another problem addressed in this work is the challenge of variable lighting conditions across a single image, to which a single threshold cannot be set. The common image processing technique, used to solve this, is known as “adaptive thresholding,” with a work from 2007 [25] introducing the major implementation of the technique. The problem with using adaptive thresholding is that it requires a predefined ratio of dark to light, and if that is not available, then it requires the use of Otsu’s method with the aforementioned challenges. All in all, the proposed work uses a new technique of (i) analyzing porosity on a pore-by-pore basis, (ii) filling pores (until overflow is detected) using the watershed technique, and (iii) analyzing the properties of the fill.

1.4 Scope of the Work.

The long-term goal of this research work is to fabricate biocompatible, biodegradable, and porous bone tissue scaffolds, designed to incorporate autologous human bone marrow mesenchymal stem cells, for the treatment of osseous fractures, defects, and structural deformities in clinical practice. In pursuit of this goal, the objective of the work is to introduce a robust, image-based method for identification and subsequently characterization of the surface porosity and ultimately the dimensional accuracy of additively manufactured bone tissue scaffolds, with a focus on pneumatic micro-extrusion (PME) process. The aim is (i) to detect complex surface pores and (ii) to characterize scaffold surface porosity using a range of developed quality metrics/traits based on top-layer or layer-by-layer micrographs acquired from PME-fabricated scaffolds.

While several image processing techniques have been presented in the literature (as discussed in Sec. 1.2), there are few works, which may be capable of detecting and characterizing the complex morphology and the intrinsic gradual pores (defined as pores that are not discretely, completely bound) of PME-fabricated bone scaffolds as exemplified in Fig. 1; note that due to the unique depth of the scaffold pores, there is an intensity gradient associated with each pore, leading to the formation of edges that are not clearly defined. The method implemented in this work is unique in that it is the only method known, which can detect and label gradually forming pores from a single camera image under a robust set of lighting conditions. In addition, a set of morphology/porosity metrics are introduced and integrated into the algorithm for quantitative characterization and analysis of scaffold surface porosity; it is shown that the porosity metrics also have a correlation with scaffold mechanical (and other functional) properties.

Fig. 1
(a)–(c) Porous, biocompatible, and biodegradable bone scaffolds, allowing for cell incorporation and diffuse proliferation. The scaffolds are composed of complex “gradual pores,” characterized by intensity gradients and no discretely defined boundaries [26,27].
Fig. 1
(a)–(c) Porous, biocompatible, and biodegradable bone scaffolds, allowing for cell incorporation and diffuse proliferation. The scaffolds are composed of complex “gradual pores,” characterized by intensity gradients and no discretely defined boundaries [26,27].
Close modal

The proposed method allows for a comparison between the porosity of a layer-by-layer fabricated scaffold and that of a reference model and thus verifying whether the newly fabricated scaffold meets porosity/morphology requirements or not. Such a capability, in addition, will pave the way for near real-time characterization of scaffold dimensional accuracy as well as process monitoring by determining (i) how well a printed scaffold matches its target CAD geometry and (ii) how deviation from the target geometry is likely to weaken the fabrication process (given the same material deposition rate). In this work, the material, proposed method (including model validation), as well as experimental designs, is presented in Sec. 2. The results of the work are discussed in Sec. 3. Finally, the conclusions and future work are delineated in Sec. 4.

2. Material and Methods

2.1 Pneumatic Micro-Extrusion Process.

PME is a high-resolution material-extrusion additive manufacturing method, which has been widely utilized for the fabrication of biological tissues, structures, and constructs. Having a layer resolution of 100 μm, the PME process allows for non-contact, multi-material deposition of a wide range of functional bio-inks for tissue engineering applications. However, the PME process is inherently complex (highly nonlinear), governed by complex multi-physics phenomena (e.g., phase-change, coalescence, receding-relaxation, and wetting equilibrium). The PME complexity, in addition, stems from the presence of a broad spectrum of process parameters—e.g., material viscosity, flow pressure, material deposition temperature, and solidification rate—as well as material–process interactions. Consequently, investigation of the effects of significant process parameters (and their interactions) in addition to the assessment of the dimensional and functional properties of fabricated bone scaffolds would be required for optimal fabrication of mechanically strong, dimensionally accurate, and patient-specific bone scaffolds.

The medium of transport and deposition in the PME process is a high-pressure flow of air, supplied by a compressor. As demonstrated in Fig. 2, the high-pressure inlet flow is injected into the deposition head’s cartridge, loaded with polymer material in powder form and heated to a temperature above the melting temperature of the loaded polymer; this leads to the formation of a molten polymer flow in the cartridge prior to deposition [27]. Subsequently, the molten polymer flow is deposited on a heated/cooled substrate via a converging micro-capillary nozzle. With the aid of a computational fluid dynamics (CFD) model, Yeow et al. [28] demonstrated that the transport of molten PCL through a micro-capillary nozzle (200 µm in diameter) under a flow pressure of 550 kPa would be a viscous flow, characterized with a Reynolds number (Re) of ≪ 1. This is unlike other direct-write additive manufacturing methods (e.g., aerosol jet printing) where material deposition is intrinsically turbulent [2931]. Filtered air is delivered to the chamber at a constant rate using a fan, which aids in maintaining a fixed level (rate) of polymer solidification (and thus layer adhesion) after deposition. In general, an optimal substrate temperature is critical for proper layer adhesion and thus accurate material deposition [32].

Fig. 2
(a) The main components of the PME process and (b) PCL deposition on a heated glass substrate (leading to bone scaffold formation). The filtered air is introduced to the chamber at a constant rate, which aids in obtaining a constant polymer solidification rate.
Fig. 2
(a) The main components of the PME process and (b) PCL deposition on a heated glass substrate (leading to bone scaffold formation). The filtered air is introduced to the chamber at a constant rate, which aids in obtaining a constant polymer solidification rate.
Close modal

2.2 Material.

Composed of PCL (Cellink, Boston, MA, USA), porous bone scaffolds were fabricated using the PME process. A polyester-based polymer, PCL is derived from caprolactone monomer using ring-opening polymerization. It is not only biocompatible and biodegradable, but also semi-crystalline and hydrophobic. As discussed in Sec. 1.4, the long-term goal of this work is to treat osseous fractures. Hence, fabricated bone scaffolds need to be biocompatible (to create a non-cytotoxic environment for cell growth and diffuse proliferation) as well as biodegradable (to facilitate scaffold-free, autologous bone remodeling and metabolism). Jansen et al. [33] observed that scaffold hydrophobicity, in addition, is a significant design factor that enhances the treatment of bone defects.

In this study, the PCL material (in powder form) was used as received. It has a molecular weight (Mn) and density of approximately 80,000 and 1.145 g/mL (at 25 °C), respectively. In addition, PCL has a melting temperature in the range of 59–64 °C as well as a glass transition temperature of −60 °C. Furthermore, PCL has a tensile strength and elasticity modulus of 16 MPa and 0.4 GPa, respectively [34]. The PCL material was stored in a freezer at −80 °C prior to consumption to prolong the shelf life of the material (approximately one year).

In the second portion of experiments, bone scaffolds were instead fabricated using a hydroxyapatite-polysaccharide starch material. This was to allow for more rapid material deposition and data collection, consisting of layer-by-layer image acquisition. The material was synthesized with 11.42% polysaccharide potato starch (Pure Supplements Co., Lindon, UT, USA), 8.43% hydroxyapatite (CAS number: 1306-06-5, molecular weight: 502.31, Sigma-Aldrich, St. Louis, MO, USA), and 4.00% titanium dioxide (Pure Supplements Co., Lindon, UT, USA) by mass (merely used to change the color of the synthesized materials from translucent to opaque), with the remainder of the weight consisting of water.

2.3 Experimental Design.

As demonstrated in Table 1, a single-factor experiment was designed and conducted with the aim to fabricate bone scaffolds with three levels of porosity, i.e., low, medium, and high (in order to assess the performance of the proposed pore detection method). The fabricated scaffolds were correspondingly based on three triply periodic minimal surface (TPMS) designs, i.e., Schwarz Diamond (Design 1), Schwarz Primitive (Design 2), and Neovius (Design 3). Note that TPMS surfaces repeat periodically in 3D and intrinsically have zero mean curvature [3537]; in other words, at each point, the sum of the principal curvatures is zero. Klemstine et al. investigated the mechanical properties of a wide range of TPMS designs [38]. Being of two intertwined congruent labyrinths, the Primitive design resembles an inflated tubular cubic lattice [39], while the Diamond design resembles the diamond bond structure (i.e., tubular and inflated in shape) consisting of two intertwined congruent labyrinths [39]. Both designs have a high surface-to-volume ratio and thus are suited for the fabrication of porous scaffolds. In addition, the Neovius design (as a cubical unit cell) is composed of outward necks that are extended toward the middle of each edge [35,39]. The Schwarz Primitive design, the Schwarz Diamond design, and the Neovius design are mathematically defined by Eq. (1)(3), respectively [40], as follows.
cosx+cosz+cosy=0
(1)
sinxsinysinz+sinxcosycosz+cosxsinycosz+cosxcosysinz=0
(2)
3(cosx+cosy+cosz)+4(cosxcosycosz)=0
(3)
Table 1

PME process parameters, defined in this study, for fabrication of porous bone tissue scaffolds

ParameterTypeLevel (unit)
Variables
Scaffold designDesign
  1. Schwarz—Diamond (Design 1)

  2. Schwarz—Primitive (Design 2)

  3. Neovius (Design 3)

Fixed Parameters
Scaffold dimensionsDesign15 × 15 × 15 (mm) (Cubic)
Layer height (thickness)Design200 (µm)
Infill patternDesignConcentric
Nozzle sizeMachine200 (µm)
Bed temperatureMachine10 (°C)
Print speedMachine2.5 (mm/s)
Deposition head temperatureMachine180 (°C)
Deposition flow pressureMachine300 (kPa)
PreFlow delayMachine900 (ms)
PostFlow delayMachine400 (ms)
ParameterTypeLevel (unit)
Variables
Scaffold designDesign
  1. Schwarz—Diamond (Design 1)

  2. Schwarz—Primitive (Design 2)

  3. Neovius (Design 3)

Fixed Parameters
Scaffold dimensionsDesign15 × 15 × 15 (mm) (Cubic)
Layer height (thickness)Design200 (µm)
Infill patternDesignConcentric
Nozzle sizeMachine200 (µm)
Bed temperatureMachine10 (°C)
Print speedMachine2.5 (mm/s)
Deposition head temperatureMachine180 (°C)
Deposition flow pressureMachine300 (kPa)
PreFlow delayMachine900 (ms)
PostFlow delayMachine400 (ms)

Listed in Table 1, the PME process parameters were set at their optimal values, based on the authors’ prior characterization studies [5,27,38,41]. PCL powder was loaded into the cartridge and maintained at 180 °C. To ensure steady-state and uniform material deposition, the loaded PCL was kept in the heated cartridge for 15 min prior to deposition. Supplied by an oilless, rust-free air compressor, the flow pressure was set to 300 kPa. Laminar molten PCL flow was deposited on a cooled glass surface (uniformly kept at 10 °C) with a print speed of 2.5 mm/s, using a 200-μm nozzle. The scaffolds were fabricated using a layer height of 200 μm, as illustrated in Fig. 3.

Fig. 3
Composed of PCL, bone scaffolds with different morphological characteristics, fabricated using the PME process
Fig. 3
Composed of PCL, bone scaffolds with different morphological characteristics, fabricated using the PME process
Close modal

The geometry of the fabricated scaffolds is cubic, having dimensions of 15 × 15 × 15 mm. The PME platform used in this study was BIO X (Cellink, Boston, MA, USA). In addition, DNA Studio (Cellink, Boston, MA, USA) was the slicer software of choice, used to convert the 3D models into a G-code and thus to create a tool-path for fabrication. Exemplified in Fig. 3, monochromatic images were captured off-line from each scaffold using a high-resolution (9.1 MP) charge-coupled device (CCD) camera (GS3-U3-91S6M-C, FLIR Systems, Inc., Richmond, BC, Canada). The imaging system is illuminated by a LED ring light having a brightness of 30,000–40,000 Lux and a color temperature of 6000 K (AmScope, Irvine, CA, USA). The acquired images were saved in raw8 data format and .tif file format. In the second part of the experiments (discussed in Sec. 3.4), the layer-by-layer images were captured in situ (after material deposition) using the process camera toolhead (Cellink, Boston, MA, USA) mounted on the PME system with the chamber lights on.

The compressive properties of the fabricated bone scaffolds were measured using a compression testing machine (MTI-10 K, Measurements Technology Inc., Marietta, GA, USA). As shown in Fig. 3, the fabricated scaffolds were cubic, having a gauge length of 10.5 mm. For the compression test, they were placed between two compression plates and subjected to a compression load applied with a rate of 0.08 mm/s. The elasticity modulus of each specimen was accurately calculated using a computer program developed in the matlab environment, processing the stress–strain data points obtained from the compression test.

2.4 Image-Based Scaffold Pore Characterization.

The following sections will describe: (i) methods to find the dimensional boundaries of a 3D-printed bone scaffolding sample, (ii) the process of edge detection, continuous edge identification, pore-highlight discrimination, and (iii) the methods of filling identified pores. Finally, there will be a discussion of the final porosity analysis and validation. Figure 4 illustrates a flowchart, detailing the main steps of the proposed method for characterization of porosity and analysis of dimensional accuracy of additively manufactured bone tissue scaffolds.

Fig. 4
A flowchart depicting the complete process of the proposed method, including noise filtering steps. Note the more robust edge-focused methodology, where edges from Canny edge detection are the building block, and thresholding is not done until the pore-filling block.
Fig. 4
A flowchart depicting the complete process of the proposed method, including noise filtering steps. Note the more robust edge-focused methodology, where edges from Canny edge detection are the building block, and thresholding is not done until the pore-filling block.
Close modal

2.4.1 Scaffold Bounding.

The aim of scaffold bounding is to find the initial dimensions of a scaffold. The photos of the samples are slightly rotated in one direction or the other. The sample photos contain noise in the form of material threads, leftover from the PME process. Therefore, a simple method was devised to effectively bound the samples for dimensional measurements.

Initially assuming that the image is correctly rotated, the objective is to find the correct positioning for a bounding box, bounding only the solid scaffold structure, discarding the slight wisps of material and protruding strings, as shown in Fig. 5.

Fig. 5
(a) A real monochromatic image of a fabricated porous bone scaffold, captured using a high-resolution CCD camera and (b) image thresholding, image rotation, as well as scaffold bounding
Fig. 5
(a) A real monochromatic image of a fabricated porous bone scaffold, captured using a high-resolution CCD camera and (b) image thresholding, image rotation, as well as scaffold bounding
Close modal

To start, an initial bounding box encloses all visible material within the image, but in order to discard noise, the box boundaries are incrementally shrunk, monitoring the amount of material discarded in each incremental step. Once the rate of change of bounded material within the bounding box exceeds a specific threshold parameter, collision with the solid structure of the sample is assumed, and the process terminates.

To find the optimal image rotation, the bounding process is performed multiple times at different angles, and through successive approximations, the angle which results in the bounding box with the least resulting area is used, as the structure is now correctly rotated.

2.4.2 Edge Tracing and Consolidation.

The initial processing step to each bone scaffolding printed sample is to utilize matlab’s Canny edge detection algorithm to find the initial edges of what could be pores, highlights (reflection), or other surface features, as illustrated in Fig. 6. The Canny edge detection has a thresholding parameter that defines how sensitive it is in defining edges along differences in pixel intensities. A threshold is manually selected to be the minimum needed to identify all major surface pores; as going beyond the minimum, the threshold introduces unnecessary noise from surface highlights. This step needs only to be done once given specific lighting conditions and is generally inconsequential to final measurements, as the following algorithm is quite robust in discarding edges not corresponding to surface pores. Before the Canny edge detection algorithm was run, a slight Gaussian blur was added to assist in reducing noise in the scaffold image. The usage of the Canny edge algorithm is as follows: (i) edges are grouped into continuous edges, (ii) the curvature of these edges is analyzed in comparison to image intensity, and (iii) the curvature is used to determine pores, highlights, and those features which need further analysis.

Fig. 6
(a) A PME-fabricated bone scaffold, composed of PCL and (b) after Canny edge detection with an appropriate threshold. Subsequently, the detected edges are color-coded after the consolidation of unique edges and further refined by pore highlight (reflection) removal.
Fig. 6
(a) A PME-fabricated bone scaffold, composed of PCL and (b) after Canny edge detection with an appropriate threshold. Subsequently, the detected edges are color-coded after the consolidation of unique edges and further refined by pore highlight (reflection) removal.
Close modal

Once an image is generated through Canny edge detection, the data representation is merely a set of discrete pixels representing edges of contrast. To begin, because a continuous edge often corresponds to a unique surface feature, a marching algorithm is used to chain continuous edges together by analyzing pixel-by-pixel contacts, resulting in a list of separate non-contacting edges. The marching algorithm is more sophisticated than needed for simply grouping contacting lines of pixels. As the algorithm traces the edge, the path trajectory is monitored in order to simultaneously record the curvature along the current point of the edge, which will become relevant in the next step. Details on the path trajectory and the curvature calculation are delineated in the Appendix (Sec.  A.1 and  A.2).

2.4.3 Initial Pore-Highlight Discrimination Through Gradient-Based Curvature Analysis.

Early on, it is important to filter out as many surface edges as possible which correspond to highlights (as opposed to pores) to minimize computation time and to minimize the noise passed through to the next stages of the algorithm.

From an initial consideration, if the curvature along a continuous path and divergence of the image intensity could be calculated, then it would be simple to discriminate between pores and highlight surface features. This could be done by stepping along the path, taking the dot product of a curvature vector (pointing inwards along the radius), multiplied by the gradient vector, followed by integrating the result. However, due to the non-continuous nature of the problem, with the intensity being represented by individual pixels on a grid, there is more nuance in the process. In Fig. 7, a 3D height mapping of image intensities is shown for further reference on how the intensity gradient corresponds to pores as indentations and highlights as tall peaks.

Fig. 7
A 3D heightmap comparison of a gradual forming pore on the surface of a scaffold and its corresponding image reference (left). The challenge of identifying the pore is the lack of a discrete border which defines where it begins and where it ends, as it is a gradually forming indentation.
Fig. 7
A 3D heightmap comparison of a gradual forming pore on the surface of a scaffold and its corresponding image reference (left). The challenge of identifying the pore is the lack of a discrete border which defines where it begins and where it ends, as it is a gradually forming indentation.
Close modal

The gradient of image intensity must first be calculated along every point of the edge being analyzed. Due to the discrete nature of a raster-pixelized image, calculating the gradient is done as an approximation. The aim is to calculate the general direction of intensity difference. It is calculated by taking the average gradient along the x-axis, two pixels ahead and two pixels behind the target pixel; this gives the x-component of the gradient at that pixel. In addition, the same is done to calculate the y-component of the gradient. A visual diagram is given in Fig. 8(a) of the method described, and the resulting gradients are shown in Fig. 8(b). Note that in both parts of the figure, the gradient vectors are normalized to more easily visualize the direction of the gradient.

Fig. 8
(a) A demonstration of the calculation method of gradients for each pixel and (b) the resulting normalized gradient vectors of the detected edge
Fig. 8
(a) A demonstration of the calculation method of gradients for each pixel and (b) the resulting normalized gradient vectors of the detected edge
Close modal

Using the calculated gradients and calculated curvatures of the path along each edge, a numerical value is calculated. Due to the discrete nature of raster images and the specific objective of merely discriminating between pores and highlights, divergence along the curve is calculated in a modified manner for this problem. To begin, the algorithm starts marching through individual image pixels along the detected edges. As the algorithm marches along the edges, pixel by pixel, of each surface feature, referring to Fig. 20 (in the  Appendix), the difference in the angle between the average direction vector and the newly selected direction is calculated at each step along an edge. The sign convention used is determined as follows: a delta difference vector is calculated between the new direction selected and the average direction vector. The result is dot-product multiplied by the intensity gradient vector at that point along the edge; if the result is positive, then the angle is positive; if the result is negative, the angle is negative. Thus, macroscopically, the angle difference in trajectory is positive if the trajectory is curving toward the gradient, as in a highlight.

By taking the average trajectory angle difference along each point of an edge, with the sign convention described above, a discretized version of divergence can be calculated along an edge, and because the image intensity gradient affects the sign convention only, it would be more accurate to refer to the calculation as curvature toward the intensity gradient from this point forwards. The result is shown in Fig. 9.

Fig. 9
A demonstration of the process for discriminating pores from highlights with positive curvatures (+ curv.) as well as negative curvatures (- curv.); please note that intensity is indicative of absolute curvature intensity
Fig. 9
A demonstration of the process for discriminating pores from highlights with positive curvatures (+ curv.) as well as negative curvatures (- curv.); please note that intensity is indicative of absolute curvature intensity
Close modal

At this point, while most edges have been classified as a pore or a highlight, some edges have values with an insignificant leaning to either end of the range of values. These are usually in the form of straight-line edges due to a surface imperfection on the bone scaffold. To filter these out, a simple post-processing step is done which checks the ratio of width to height of unclassified features, easily discarding straight edges before continuing. Details on discarding uncertain edges via analysis of minimizing dimensions are given in the Appendix (Sec.  A.3).

2.4.4 Pore Filling Algorithm by Thresholding Until Collision.

The purpose of this step was to analyze the surface porosity and quality of a fabricated bone scaffold; note that the edge-analysis described earlier is simply a precursor step. It can be observed from Fig. 9 that the defining feature of a pore appears darker than that of its direct surroundings. However, thresholding the image, selecting pixels below a specific brightness, and calling them a pore are not sufficient, as different regions of the image require different threshold limits, and it is not trivial to define such. Therefore, the first purpose of defining unique edges corresponding to pores was to be able to go through on a pore-by-pore basis.

The major challenge of this work was that there was no simple way to define regions, which do and do not correspond to a pore, as the majority of the pores are gradual, as in, they are not completely bounded; instead, they tend more to “spill in” being bound by a “C” shaped boundary, as shown in Fig. 7. As a result, it can be observed in Fig. 10 (in the 3D height map) that it is similar to an indentation in a surface, and if it were to be filled with water, similar to a puddle, it could be filled until it spilled over the edge. Therefore, that was the direct approach taken in the following steps described below, in order to identify the regions near an edge that correspond to a pore. For further information regading 3D height map techniques, referred to as the “watershed technique”, please see Ref. [23].

Fig. 10
(a)–(f) The process of filling up a pore in the presence of noises (such as highlights)
Fig. 10
(a)–(f) The process of filling up a pore in the presence of noises (such as highlights)
Close modal

Given all edges which have the possibility to bound pores, the pore-filling algorithm functions as follows: (1) The pixels drawing the edge are extracted into an image, and the same region is extracted from the original grayscale picture of the printed bone scaffolding sample. (2) Darker pixels below an initial threshold from the greyscale image are extracted. That thresholded pore is checked for “complete intersection” with the corresponding edge. (3) The threshold is iteratively increased, with lighter pixels forming the “puddle,” until the intersection with the edge exceeds a certain amount, and the process terminates when the filling has spilled over.

Note that for strangely shaped pores, a robust intersection algorithm was developed named “complete intersection,” as detailed in the Appendix (Sec.  A.4). In addition, there is a final major noise filtering step (following the pore filling), which checks whether pores are filled outside-in (instead of inside-out, in which they are labeled as noise and omitted). The noise filtering step also has been delineated in the Appendix (Sec.  A.5).

2.4.5 Analysis Methods of Produced Surface Porosity Mapping.

At the end of the process, a mapping is produced with color-coded, distinct pores in an array of images which are compiled together and consolidated; note that if two fillings overlap, they are consolidated. Using this image array, many metrics of surface quality can be determined. There were two routes of image analysis taken in the pursuit of quantifying quality using this array of pore images. One route that has taken of quantifying quality consists of calculating consistency metrics of the pore image array, such as variation in pore areas and pore distances, while the other, more comprehensive route, compares the produced pore image array directly to the bone scaffold model geometry.

Given data on each individual pore on the surface of the bone scaffold as shaded regions of a grid, the quality of the bone scaffold can be quantified to a significant extent without even referring to the target geometry of the surface. The reference-independent analysis, to be described, makes assumptions that the surface pores are regularly spaced in a consistent pattern and that all pores on the surface of the scaffold are expected to be at the same size.

Thus, two metrics are produced from the reference-independent analysis step: the relative standard deviation of pore areas (rsd_area), and the relative standard deviation of nearest-neighbor distance (rsd_distance). More specifically for the latter, every pore on the surface is iterated over and the distance to the closest pore is added to a list. The relative standard deviation is then calculated over these distances. The weighted average of both the “rsd_area” and “rsd_distance” metrics is then taken and referred to as the “precision error” in the rest of this work. The ratio is 1:2, respectively. A perfect printed scaffold would have zero standard deviation in both metrics, assuming the model geometry met the original assumption stated earlier.

The second route of analysis involves a comparison to the target model (CAD) geometry. This is achieved by taking a top-down projection of the target model geometry to extract a shaded cross-section. The cross-section image is then aligned with the actual surface pore-shaded image data, and the overlap can then be investigated. To ensure that the two images are properly aligned, their alignment is iteratively optimized until an alignment is found, which minimizes the area not overlapping between the two images (as illustrated later in Fig. 17).

The advantage of this reference-dependent analysis route is that it can determine between insufficient material deposition or excess material deposition with a quantitative amount. If the printed pores are too large compared to the target geometry, meaning that insufficient material was deposited, then the data are representing under-extrusion; if the printed pores are too small compared to the target geometry, then the data are representing over-extrusion. This extrusion error can then be quantified by a single metric, called “percent over/under-extrusion” with over-extrusion being positive and under-extrusion being negative. An additional metric is calculated, called “percent deviation from reference,” which takes the absolute value of all deviation from the target model. The outcome of this metric is that while a model with half the pores over-extruded and half the pores under-extruded, may net zero overall error, the “percent deviation from reference” metric would quantify the absolute error.

3 Results and Discussion

3.1 Validation of the Proposed Method.

During the development of the work, in order to validate the method developed, a test sample with known parameters was designed and 3D rendered with a similar opaque white material with protrusions and pores behind the pores, as illustrated in Fig. 11. The sample was designed with a surface porosity of 38.22%, and a dimensional ratio (width to height) of 1:1.1. The proposed method detected and labeled every unique pore on the surface of the simulated sample (shown in Fig. 11). In addition, it calculated a surface porosity of 38.12. The measured dimensional ratio through the initial sample bounding was 1:1.099 (width to height).

Fig. 11
Initial validation of the proposed method using varied simulated geometries. Colored regions correspond to detected pores. The algorithm appears to identify every surface-level pore in the rendering, based on inspection of the colored regions.
Fig. 11
Initial validation of the proposed method using varied simulated geometries. Colored regions correspond to detected pores. The algorithm appears to identify every surface-level pore in the rendering, based on inspection of the colored regions.
Close modal

However, this simple model is not complex enough to validate the method proposed in this work. In Fig. 11, the pores are clearly defined, and the model has no complex noise and geometry beyond the few edges placed behind the pores. In addition, this model does not contain any pores, which gradually develop and are vaguely defined. One of the claims of this work is that in regards to quantifying the quality of a printed scaffold with a single photo, the proposed method is robust to different lighting environments and materials, functioning with a single set of parameters and solving for the ones that remain. In order to test this hypothesis, a more thorough validation is performed, as follows.

A more robust validation approach was taken by directly implementing a complex geometry into a proper, production rendering environment. Using Blender’s Cycles, a benchmark was produced (in a ray-tracing rendering environment) by simulating a 3-meter-by-3-meter laboratory room, with standard white walls. The design selected for rendering was based on a triply periodic minimal surface (i.e., Schwarz Gyroid [38]). It was placed in the center of the room on a black floor to be tested under various lighting conditions for validation, as illustrated in Fig. 12. The camera in the simulation was a low field-of-view camera placed close to the scaffold to be similar to the camera module mounted on the BIO X 3D-printing system, which was used to photograph many of the fabricated bone scaffolds in this work. In the rendering environment, to match a complex, reflective, real-world material, a physically based material was developed to model a white, glossy, reflective material similar to PCL’s visual appearance. For validation of material invariance, a brown, rough physically based material was generated as well; this is perhaps representing a wood-filled material or similar; however, it was selected to be in contrast to the glossy white material, to be used in the comparison.

Fig. 12
A more sophisticated approach to rendering a realistic benchmark. The Schwarz gyroid model design was loaded into a physically based rendering software (Blender Cycles), and a simulated material was composed to match a glossy white PCL polymer. The complexity can be observed in the rendered output (right), containing only gradually forming pores, and highlights from reflected light.
Fig. 12
A more sophisticated approach to rendering a realistic benchmark. The Schwarz gyroid model design was loaded into a physically based rendering software (Blender Cycles), and a simulated material was composed to match a glossy white PCL polymer. The complexity can be observed in the rendered output (right), containing only gradually forming pores, and highlights from reflected light.
Close modal

There were three different distinct lighting environments that were selected with the intention of (a) creating variation in the visual appearance of the scaffold design, and (b) modeling possible different lighting environments that may be found in a standard lab or workspace setting. The first lighting environment tested in the simulated environment was the presence of a floor lamp with a white 50-W bulb, positioned close to a meter away from the scaffold, shown in the first row of Fig. 13. This offset light source causes a pronounced cast shadow in the deeper gradual indents of the pores, causing a perception of closer to nine pores on the surface. The next lighting environment is a large white fluorescent light panel in the ceiling of the room slightly offset, commonly found in laboratories, with the rendering shown in the second row of Fig. 13. And the third lighting environment simulated was the presence of a bright flashlight around half a meter directly above the sample. This also illuminates the gradual indents of the scaffold, but harsh reflected highlights can be seen produced as a result; this can be observed in the third row of Fig. 13. The final row of Fig. 13 displays the brown material variation in the 50-W lamp environment for comparison to glossy white polymeric material.

Fig. 13
Depicted in the figure in the first column are four different simulated lighting environments, with (a) a 50-W floor lamp, (b) a fluorescent ceiling lighting, (c) a flashlight shining 50 cm above, and (d) the 50 W floor lamp, but with a clay material. The second column depicts pore-filling image arrays produced by the proposed method. Each consolidated pore is marked with a “+” at the centroid. The color hue depicts relative size, with the redder hues representing larger pores and the bluer hues representing smaller pores. The final column depicts a comparison in pore filling to the standard method used by Bowoto et al. [22]. (Color version online.)
Fig. 13
Depicted in the figure in the first column are four different simulated lighting environments, with (a) a 50-W floor lamp, (b) a fluorescent ceiling lighting, (c) a flashlight shining 50 cm above, and (d) the 50 W floor lamp, but with a clay material. The second column depicts pore-filling image arrays produced by the proposed method. Each consolidated pore is marked with a “+” at the centroid. The color hue depicts relative size, with the redder hues representing larger pores and the bluer hues representing smaller pores. The final column depicts a comparison in pore filling to the standard method used by Bowoto et al. [22]. (Color version online.)
Close modal

The second column of Fig. 13 displays the visual pore, highlighting the results of the method described in this work. The color represents the size of an individual pore which was identified, with a marker representing the centroid of the pore; note that redder colors represent relatively larger pores and bluer colors represent relatively smaller pores.

The third column of Fig. 13 displays the efficacy of a standard methodology for identifying regions of porosity in recent scientific literature and industry. Specifically, the standard method used for a benchmark comparison is the method recently used by Bowoto et al. [22] for in situ monitoring of internal porosity using an optical camera for layer-by-layer material deposition. The method used functions on contrast thresholding, where the image is converted to grayscale, and a specific threshold is selected, where any pixel darker than a certain value is marked as porosity. The white-marked regions (in the third column) can be observed in comparison to the individually highlighted regions (in the second column). Note that the highlighted regions of the standard method of the third column all appear the same color and without markers, as the method used was for a comparison of pore-filling techniques.

To quantitatively demonstrate the robustness of the proposed method, Table 2 contains quantitative data scores for each of the four simulated environments in Fig. 13. The objective was to have similar quality metrics irrespective of lighting conditions. Any standard deviation demonstrates the quantified error due to variable lighting conditions of the algorithm.

Table 2

Measurements of comparison for the identical scaffold model, rendered in different lighting conditions, to test the robustness of the model

EnvironmentRSD Area (%)RSD Dist. (%)Precision Error (%)
(a) 50 W Floor Lamp31.6318.3122.75
(b) Fluorescent Ceiling Lighting43.5011.6422.26
(c) Flashlight (50 cm above)81.3923.5542.83
(d) 50-W Floor Lamp, Clay Material9.161.514.06
Standard (Std.) Deviation30.219.5115.84
EnvironmentRSD Area (%)RSD Dist. (%)Precision Error (%)
(a) 50 W Floor Lamp31.6318.3122.75
(b) Fluorescent Ceiling Lighting43.5011.6422.26
(c) Flashlight (50 cm above)81.3923.5542.83
(d) 50-W Floor Lamp, Clay Material9.161.514.06
Standard (Std.) Deviation30.219.5115.84

Note: The “RSD Area” represents the percent relative standard deviation (RSD) of pore areas, with “RSD Distance” representing the percent RSD of pore pattern spacing distance, and with the “Precision Error” being the weighted average of the two.

3.2 Comparing the Presented Approach to the Standard Method

3.2.1 Comparing the Accuracy of Pore Filling.

In comparing the accuracy of the pore-filling algorithms on the benchmarks, it is important to establish the objective of the pore-filling algorithm, and the challenge with gradually forming pores. The challenge with identifying the porosity of a printed scaffold through only optical photographs is that without the use of depth mapping, the only information that can be gathered on the scaffold is through reflected light. If the scaffold has discrete, sharp edges, then identifying the surface porosity is simple; however, with gradually forming smooth edges, as shown in Fig. 9, it is not trivial to define a point along the dropping curvature where a pore is defined. Furthermore, identifying the same point along the curvature in dynamic lighting conditions is a complex challenge and beyond the scope and objectives of this work. The objective was not to identify the percent of the surface, which was marked as belonging to a pore, but to define quality metrics that measure the consistency of the porosity, which can be approximated in most lighting conditions (as discussed in Sec. 2.4.5).

In comparison to the “gold standard” technique for identifying porosity, as used in works such as Bowoto et al. [22]; in 2020, thresholding techniques are the dominant technique due to their simplicity and efficacy. Testing the standard method on the rendered images of the scaffolds produced quality results across lighting conditions after being processed through the wisp cropping module defined in this work. The standard thresholding techniques function well with a well-selected threshold with even lighting conditions, as depicted in the third column of Fig. 13. In comparison, the proposed method performed very similarly, except for the presence of some introduced noise due to the addition of all the processing steps.

The feature-by-feature analysis for filling the pores until collision with the boundaries has the benefit of functioning automatically across lighting conditions, while traditional thresholding techniques are highly sensitive to the selected threshold for images without high contrast. Figures 14(a) and 14(b) depict a real-world print in PCL of the Schwarz gyroid scaffold used in the renderings. Due to challenges in printing PCL, as discussed in Ref. [38], the scaffold does not appear very similar to the rendering. The final three images, i.e., Figs. 14(d)14(f), demonstrate the sensitivity of small threshold changes to the resulting porosity labeling. The standard solution is to use automatic thresholding selection algorithms, such as Otsu’s algorithm [24], used by Bowoto et al. [22]. However, Otsu’s algorithm only works on bimodal intensity distributions, where there are clear dark regions and clear light regions in the image. Bowoto et al. [22] had less challenge with this by using only diffusive materials with a light directly horizontal to the surface to create sharp shadows; also the pores that their work dealt with were sharply defined hemispherical crater defects. Otsu’s algorithm fails to function on a multimodal image as shown in Fig. 14(b), as in real images, it is rare to see the dark base plate or deep pores that are seen in the rendering. Due to compiled errors, it is more common to see complex material pattern below the surface pores. Other algorithms, such as adaptive thresholding in Ref. [25] used for optical character recognition, failed to function in testing on the large gradual pores. In complex real-world scenarios, the proposed method functioned far better in pore filling, as shown in Fig. 14(c), with seldom adjustment to program parameters.

Fig. 14
(a) A fabricated TPMS (Schwarz Gyroid-based) PCL bone scaffold, (b) the internal structure of a bone scaffold after cropping and grayscaling, (c) pore identification using the proposed method. Pore identification using a standard thresholding method (i.e., Otsu’s method [24]) at a threshold level of, (d) 0.50, (e) 0.70, and (f) 0.85 (implying that Otsu’s method unlike the proposed method fails to function on the bone scaffold image and in general on multimodal images).
Fig. 14
(a) A fabricated TPMS (Schwarz Gyroid-based) PCL bone scaffold, (b) the internal structure of a bone scaffold after cropping and grayscaling, (c) pore identification using the proposed method. Pore identification using a standard thresholding method (i.e., Otsu’s method [24]) at a threshold level of, (d) 0.50, (e) 0.70, and (f) 0.85 (implying that Otsu’s method unlike the proposed method fails to function on the bone scaffold image and in general on multimodal images).
Close modal

In addition, a quantitative analysis was performed to further assess the performance of the proposed method. It was observed that the proposed method had an accuracy of 60.86% (defined as total accurate overlap with the real pore area minus total falsely labeled area divided by total real pore area), while the standard method (i.e., based on the global contrasting method) had an accuracy of 7.74%, −5.53%, and −534.26% for the manually selected threshold of 0.50, 0.70, and 0.85, respectively, as shown in Fig. 14.

3.2.2 Validation of Analysis Methods.

In comparing data analysis methods to the standard method, the most common metrics are measurements such as the average area of porosity, standard deviation of pore size, centroid position identification, etc. An important difference of the proposed method is the preprocessing step to the data analysis, where pores that are not contacting are consolidated by estimating if different regions correspond to the same pore. This is in contrast to the standard “connected-component labeling” techniques used to simply group all contacting pixels. This method can be seen to be working effectively in the renderings in Figs. 13(a-ii) and 13(d-ii) where even groups of pixels which are not contacting are successfully grouped into, corresponding to the same pore; however, the model does introduce noise in very specific conditions, such as Fig. 13(c-ii), where it can be observed that portions of the bottom three pores were all consolidated into one pore (as shown in red for its size). This module of the proposed method is very effective in improving the consistency of results in messy, real-world scaffolds; however, the pore in Fig. 13(c-ii) was consolidated due to more subtle shadowed connections between the pores in the rendering, shown in Fig. 13(c-i).

In reviewing the data quantification metrics in Table 2, comparing the quality scores of the different scaffolds rendered in different lighting conditions, it can be seen that the issues related to using a rendering for a benchmark and the difficulty of the benchmark used, caused minor errors in the pore filling and consolidation algorithm, which led to the introduction of major noise into the quality results. Referring to Fig. 13, first, when observing the noise present in Fig. 13(a-ii), the algorithm decided that the hole pairs belonged to the same pore, as they were present in the same surface indentation, in eight out of the nine indentations, with breaking the ninth pore up into two smaller pores. This resulted in a high RSD of Areas, and RSD of Distances (Sec. 2.4.5 for details on the quality metrics). In the next environment, in Fig. 13(b-ii), the algorithm opted to break up each pore except for one, creating the same issue. In Fig. 13(c-ii), the pores were all separate, except for one additional ridge line pore as well as the major error of consolidating the bottom three pores. In addition, only in the last environment with the diffuse material, did the algorithm maintain consistency in consolidation, as shown in Fig. 13(d-ii). As a result, over the four environments, the fourth environment gave the correct answer, with the RSD of Distance being close to zero at 1.51%. The precision error, i.e., the weighted average of both metrics with a bias toward RSD of Distance had a final standard deviation of 15.84 due to the noise in variation of widely varying lighting conditions in a rendering environment. Note that it is challenging to achieve substantial lighting varying to the extent demonstrated in a lab environment, as there is always a combination of lights illuminated in the room, with more surfaces to diffuse, as compared to the simulation environment where the light sources shown were the only sources present. Interestingly, the pore fillings appear near identical in Figs. 13(a-ii) and 13(d-ii) implying that the failure to consolidate one of the pores correctly in Fig. 13(a-ii) is an outlier. Analysis of outliers will be demonstrated in Sec. 3.4.

3.3 Comparing Quality Estimates of the Proposed Method to Compression Testing Metrics.

In the authors’ prior work [38], multiple scaffold designs were selected for printing and analysis of mechanical properties. With the authors printing multiple samples of each design, the compression modulus of each of the samples was assessed in a compression testing apparatus. For this work, before the samples were destroyed (compressed), surface images were collected for three of the scaffold designs (as discussed in Sec. 2.3), demonstrated in Fig. 15(a)15(c).

Fig. 15
PME-fabricated bone scaffolds with a high, medium, and low level of pore quality, including (a) Diamond geometry, (b) Schoen I-WP geometry, and (c) Neovius geometry; (d)–(f) the production of the porosity surface mappings with redder colors representing relatively larger individual pores, and bluer colors representing relatively smaller pores (Color version online.)
Fig. 15
PME-fabricated bone scaffolds with a high, medium, and low level of pore quality, including (a) Diamond geometry, (b) Schoen I-WP geometry, and (c) Neovius geometry; (d)–(f) the production of the porosity surface mappings with redder colors representing relatively larger individual pores, and bluer colors representing relatively smaller pores (Color version online.)
Close modal

As demonstrated in Fig. 15, there was a discernable variation in manufacturing quality for the three scaffold designs; however, in Fig. 15(d)15(f), the proposed method was able to function well on the inconsistent scaffolds. The average quality metrics for the three scaffold designs are displayed in Figs. 16(a) and 16(b); it can be observed that under real-world conditions, the standard deviation for metrics on the scaffold sets was lower than that of the renderings.

Fig. 16
The data in (a) and (b) were collected from a single photo (per scaffold) of the top layer of a PME-fabricated scaffold in general lighting. Graph (a) depicts the average percent relative standard deviation (RSD) of pore areas. Graph (b) depicts the average % RSD of pore pattern spacing distances. Graph (c) depicts the average compression moduli found in a compression testing machine. The error bars represent the standard deviation in measurements across the three sets of six samples for each of the metrics.
Fig. 16
The data in (a) and (b) were collected from a single photo (per scaffold) of the top layer of a PME-fabricated scaffold in general lighting. Graph (a) depicts the average percent relative standard deviation (RSD) of pore areas. Graph (b) depicts the average % RSD of pore pattern spacing distances. Graph (c) depicts the average compression moduli found in a compression testing machine. The error bars represent the standard deviation in measurements across the three sets of six samples for each of the metrics.
Close modal

Referring to Fig. 16, using the proposed method on a single photograph in general lighting conditions, as the average quality scores of surface porosity improved (the error decreased), the average compression modulus of the PCL printed scaffolds increased. Note that the corresponding error bars represent the standard deviation for the quality metrics. There is substantial variation between the quality scores of the printed scaffolds. While there is actual variation between fabricated scaffolds, as can be seen in the standard deviation of compression modulus tests, there is significant noise in the proposed method’s final quality measurements, which is likely the source of the substantially higher standard deviation in the quality measurements.

The model geometry used in Ref. [38] was designed for strength given a certain global density. If the printed scaffold deviates from the target geometry, it is expected that the strength of the scaffold will suffer, and the low and medium-quality scaffolds in Fig. 15, which appear completely different from their target geometry patterns demonstrate such. Since each of the scaffold geometries meets the previously defined requirements for valid quality metrics (as defined in Sec. 2.4.5), any error measured by the proposed method would result in a weaker printed scaffold, assuming no measurement noise was present.

The data collected were only from that of a single photo directed at the top of the printed scaffold, which disregards the internal structure of the scaffold which may affect the scaffold strength, but as shown in Fig. 16, there is a correspondence to just the surface quality and the compression modulus. This is perhaps a result of the PME process being sensitive to internal defects, where since new layers are printed on previous layers, and depend on the quality of previous layers, a defect significantly affects quality for many layers after, especially with a biomaterial such as PCL, which takes a long time to cool. Another explanation for this phenomenon could be a result of the compression strength being dependent on the weakest load-bearing elements of the scaffold. Besides, the weakest portion of the scaffold is the lowest quality portion, which is often the last layer, as implied from Fig. 18 (discussed in Sec. 3.4).

In the authors’ prior work [38], prints with identical process parameters and printer instructions would vary between prints due to variations in environmental factors and material factors; as a result, it is unknown if printed scaffolds met the same mechanical strength as previously printed and tested scaffolds, without testing the newly printed scaffold and destroying it. The proposed method demonstrates a rough analysis of comparing the material properties of a porous printed part to those of a printed part that was previously tested through a single photo, in general lighting conditions. The issue with this approach is the noise of surface porosity mappings, which can be improved in future works and mitigated with further results (as demonstrated in Sec. 3.4 below), but more so, the issue with the proposed method in estimating quality for mechanical strength standards is the lack of information regarding the internal structure of the printed parts, which will be remedied in further results below.

3.4 Layer-by-Layer Analysis and an Extension of the Method With Target Model Awareness.

The previous results were using analysis techniques general to any scaffold design which met the requirements (defined in Sec. 2.4.5), with the analysis simply checking the consistency of porosity sizing and consistency of pattern spacing. However, a target model aware methodology opens the door for far more analysis to get more consistent results. It even allows for noise filtering, which drastically improves the consistency of results.

By taking the target scaffold geometry and using CAD software to take a cross section at a specific layer, the exact target surface porosity at that layer can then be stored as a binary image. Performing this for every printed layer height, a library of binary images is generated for comparing the printed scaffold to the target geometry at any given height.

The first and most significant advantage of this approach is that ‘after optimally matching the generated reference image to the photograph of the scaffold, the reference image is used for substantial noise reduction in discarding minor pores or indentations which are likely not significant to the overall quality of the measurements. A depiction of this more robust approach is depicted in Fig. 17.

Fig. 17
The new workflow for analyzing pore geometry given a target STL image. Depicted above is layer 11 from Fig. 19, where a quality improvement was experienced midway through the print. Image (a) depicts a cropped photo of a hydroxyapatite-polysaccharide-based PME scaffold. Image (b) depicts a cross-section image of the target geometry. Image (c) depicts the pores fitted optimally to the target geometry, with purple pixels depicting overflow (OF), and green pixels depicting underflow (UF). Image (d) depicts the final noise-filtered results with centroid positions labeled with a marker, and color depicting relative size, with redder hues depicting larger pores, and bluer hues depicting smaller pores (Color version online.)
Fig. 17
The new workflow for analyzing pore geometry given a target STL image. Depicted above is layer 11 from Fig. 19, where a quality improvement was experienced midway through the print. Image (a) depicts a cropped photo of a hydroxyapatite-polysaccharide-based PME scaffold. Image (b) depicts a cross-section image of the target geometry. Image (c) depicts the pores fitted optimally to the target geometry, with purple pixels depicting overflow (OF), and green pixels depicting underflow (UF). Image (d) depicts the final noise-filtered results with centroid positions labeled with a marker, and color depicting relative size, with redder hues depicting larger pores, and bluer hues depicting smaller pores (Color version online.)
Close modal

To demonstrate the efficacy of this approach, utilizing the 3D-printing system's camera module, a photo was automatically taken off the 3D-printed material after each layer, with the analysis program receiving the photos and building a quality metric graph in real-time. Note that for layer-by-layer experiments, the material was changed to a hydroxyapatite-polysaccharide-based bone scaffolding material which prints more quickly and tends to deviate from target geometry more quickly; this was done for more rapid data collection.

Figure 18 depicts the quality metrics measured of a sample bone scaffold, which had a gradual quality degradation as the material deposition progressed. Observing the first 20 layers, the RSD-Area and RSD-Distance (pattern spacing consistency) metrics are holding very stable, degrading gradually (as demonstrated layer-wise in Fig. 19). This consistency of the “precision error” metric exhibits the efficacy of the benefits of the model aware noise reduction, with the only major errors being on layer 21, when a pore was unrecognized, and on layer 24, when two pores were wrongly consolidated. The error spikes that follow are a result of actual errors in the fabricated scaffold. These errors can easily be resolved in future work, which takes a model-aware approach into consideration when initially performing the pore filling, as opposed to just at a final step where target geometry is only used for noise reduction.

Fig. 18
The layer-by-layer quality metrics produced in real-time by the proposed method on a PME-fabricated hydroxyapatite-polysaccharide scaffold (depicted in Fig. 19). The “RSD Area” represents the percent relative standard deviation (RSD) of pore areas, with “RSD Distance” representing the percent RSD of pore pattern spacing distance, and with the “Precision Error” being the weighted average of the two. “Over-Under Flow %” represents the percent overflow or underflow when compared to target geometry, and “% Deviation STL” represents the absolute value of deviation, when compared to target geometry. The “Overall Error” is the most important metric, being a weighted average of all the previous metrics, with more details in Sec. 2.4.5. Note the improvement in quality as the scaffold walls begin to thin more and the pores begin to open at layer 11.
Fig. 18
The layer-by-layer quality metrics produced in real-time by the proposed method on a PME-fabricated hydroxyapatite-polysaccharide scaffold (depicted in Fig. 19). The “RSD Area” represents the percent relative standard deviation (RSD) of pore areas, with “RSD Distance” representing the percent RSD of pore pattern spacing distance, and with the “Precision Error” being the weighted average of the two. “Over-Under Flow %” represents the percent overflow or underflow when compared to target geometry, and “% Deviation STL” represents the absolute value of deviation, when compared to target geometry. The “Overall Error” is the most important metric, being a weighted average of all the previous metrics, with more details in Sec. 2.4.5. Note the improvement in quality as the scaffold walls begin to thin more and the pores begin to open at layer 11.
Close modal
Fig. 19
Deviation of the PME-fabricated hydroxyapatite-polysaccharide bone scaffold over time. Note the improvement in quality as the scaffold walls begin to thin more and the pores begin to open at layer 11. Also, note the complete closing of a pore by layer 29 (as a major defect) corresponding to the spiking of the error metrics (plotted in Fig. 18).
Fig. 19
Deviation of the PME-fabricated hydroxyapatite-polysaccharide bone scaffold over time. Note the improvement in quality as the scaffold walls begin to thin more and the pores begin to open at layer 11. Also, note the complete closing of a pore by layer 29 (as a major defect) corresponding to the spiking of the error metrics (plotted in Fig. 18).
Close modal

Beyond noise reduction and improving pore filling, using a model aware approach, far more data metrics can be calculated by comparing to target geometry, with the metrics previously defined in Sec. 2.4.5. As shown in Fig. 18, the metrics “% deviation from STL” and “over-underflow %” are calculated through a model aware comparison, and the two metrics depict much more consistent results. Because in the PME printing process, the quality of a surface layer is dependent on the quality of the layer on which it is deposited, it was expected that a curve depicting a quantified measure of print quality would not have a sudden discontinuous jump in improvement (while a sudden jump in error is possible in the case of print errors). While the “precision error” metrics are seen to have sudden improvements in quality, such as from layer 21 to 22, or from layer 25 to 26, which simply corresponds to noise caused by the distortions at the end of the print, the model-aware metrics maintain an impressive consistency and ability to filter noise, even under extreme distortions found at the end of the print, shown in Fig. 19.

The resulting final metric, titled “overall error”, is a weighted average of the “precision error” metrics and the model-aware metric, utilizing all the previously collected data, and it represents a final calculated numerical quantification of quality for a singular layer of the print of the bone scaffold. Due to the wide array of data used and the improved metrics, the “overall error” score is consistent and robust to noise. Earlier results in this work depicted a correlation between the compressive strength of fabricated scaffolds and their single-image noisy “precision error” metrics, without any reference to the target geometry. Hence, it is expected that by utilizing the layer-by-layer curve for “overall error,” a much stronger correlation will be derived in future works.

4 Conclusions and Future Work

In this study, a novel image-based method was introduced for the detection of complex gradual pores and subsequently quantification of the porosity of PME-fabricated bone tissue scaffolds. The presented method is unique in the sense that it is capable of (i) detecting individual pores of varying geometry from a single camera image under a robust set of lighting conditions using a handcrafted image processing algorithm, and (ii) conducting a range of porosity analyses with or without knowledge from the target geometry of a scaffold. The proposed method allows for both labeling of scaffold pores and quantification of porosity on a layer-by-layer basis, used for the non-destructive assessment of the functional performance of fabricated bone tissue scaffolds.

It was observed that, while the algorithm was able to robustly adjust pore-filling parameters in different lighting conditions, the noise from varying lighting conditions had an impact on the quality metrics (when lacking the model-aware noise rejection analysis). In addition, in comparison to the standard approach used in recent literature, tested on the ideal scaffold rendering, after manually adjusting the threshold values, the standard approach performed comparably in filling pores to the pore filling of the proposed method. However, in real-world testing (based on PME-fabricated PCL scaffolds), the proposed method performed substantially better in pore filling and consolidation than the standard method.

To validate the applicability of the proposed method, a relationship was demonstrated between the quality scores from a single top-layer photo of fabricated scaffolds and their corresponding compressive strength, using model-unaware metrics, which led to substantial noise present in the results. However, the noise was substantially reduced by referencing binary image cross sections of the target geometry along with new quality metrics (utilizing the target geometry cross sections). Overall, by combining the model-aware and general quality metrics, the noise was near entirely eliminated.

With the noise-free layer-by-layer quality assessment, it is expected that a much stronger relationship can be derived, comparing the compressive strength of the fabricated bone scaffolds to the quantified quality error curves/trends (addressed as part of the authors’ future work). Also, the proposed method allows for a real-time comparison between layer-by-layer fabricated scaffold quality and that of a reference scaffold (meeting a target strength/stiffness) and thus is capable of verifying whether newly fabricated scaffolds meet quality standards or not. The results of the proposed method will be of significant value in a clinical environment (as well as in other fields highly sensitive to print quality, e.g., aerospace), as there is a need to validate the strength of fabricated scaffolds (having a complex internal structure) prior to implantation into a patient. It is worth mentioning that based on the model-aware approach and enhanced noise reduction, the resulting layer-by-layer pore fillings can be used to create a rough, approximate 3D model of the scaffolds for analysis of other parameters, such as estimating overall volumetric porosity. Last but not least, the proposed method is of great value in characterizing and monitoring changes/drifts in the quality of fabricated scaffolds for optimization of the PME process.

Acknowledgment

This research work was made possible by the NASA Established Program to Stimulate Competitive Research (EPSCoR)—Rapid Response Research (R3) (Grant No. 80NSSC22M0249). Dr. Salary would like to sincerely thank the WV Space Grant Consortium.

In addition, Dr. Salary would like to sincerely acknowledge the Marshall University Research Corporation (MURC) for supporting this research work via the John Marshall Scholar Faculty Award. Also, Dr. Salary would like to thank the College of Engineering & Computer Sciences (CECS) at Marshall University as well as the Cabell Huntington Hospital & Medical Center (Huntington, WV, USA).

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix

A.1 Scoring Viable Next Pixels While Following an Edge

In order to select the next pixel in the chain, all possible pixels must be scored in case there are multiple options available, as shown in Fig. 20. The objective is to favor pixels, which are closest in distance and whose direction is most similar to the present average direction. First, each viable pixel’s distance from the current pixel is calculated, with directly adjacent pixels having a distance of 1 and diagonal pixels having a distance of approximately 1.41. Subsequently, to calculate the similarity of trajectory, the running average of the direction vectors of the previous few pixels in the chain is calculated, with the direction vectors being shown in Fig. 20. The average direction is then compared to all the new possible directions through the use of the dot product. Finally, the score of each pixel is calculated by taking the negative of the distance and adding it to one-fourth of the dot product, with one-fourth being the factor in order to reduce the impact of direction so that minimal distance is always favored.

Fig. 20
A selected pixel in a chain with two possible options for the next pixel in the chain. The final arrow corresponds to the average direction of the previous few steps.
Fig. 20
A selected pixel in a chain with two possible options for the next pixel in the chain. The final arrow corresponds to the average direction of the previous few steps.
Close modal

Technically, if no scoring were performed, the edge labeling algorithm would still be fully capable of finding all connected edges. However, the purpose behind such scoring will become clear in the following discussion of edge curvature analysis. The average of the previous eight pixels' direction differences is used to calculate direction. Its selection was roughly through trial and error according to the following criteria: while the current direction should be as tangential as possible to the current path, due to the pixelized discrete nature, a small pixel history is used to smoothen out the edges. Small particles on the surface are blurred and less likely to be detected, while larger features (such as pores) become smoother, resulting in more continuous edge detection.

A.2 Arbitrary Starting Points and Handling Skipped Portions of Edges

Up until this point, the method described has specified what takes place given a starting pixel, and nothing has been mentioned as to what happens to a skipped off-shoot of an edge due to a branching point, as schematically shown in Fig. 20. At the beginning of the process, after edge detection, it was mentioned that xy coordinates of every pixel from the image generated by the Canny-algorithm were placed in a list. Simply, the first pixel used to begin the edge tracing is the first pixel from that list, and then after that, the edge trace terminates due to a dead end; the next starting pixel selected is the next unlabeled pixel in that list. The result is that no matter what edge is skipped due to a branching point or due to starting in the middle of an edge, the starting point of the next pixel is located somewhere on the unlabeled edge; hence, inevitably all edges are traced, regardless of the arbitrary starting pixels.

A.3 Discarding Uncertain Edges via Analysis of Minimizing Dimensions

Some edges have uncertain curvature averaging near zero due to some pores having “S” shaped edges, where for some of the edge, the curvature is bounding a pore, while for the other portion of the edge, the curvature is bounding a highlight or surface protrusion. Also, it is hard to discriminate between straight lines with little curvature and edges with varying curvature along the path. To solve this problem in the next stage of the algorithm, in order to minimize the amount of non-pore edges passed on to the next step, the minimal bounding dimensions of the edge are considered.

An edge with little curvature, a straight line, would have a minimal bounding box with a much more extreme ratio of width to height. As in, a straight line may be minimally bounded by a box with a dimensional ratio of 7:1, while an “S” shaped edge with the same curvature could only be bounded by a box with a ratio of 3:1. The challenge is that the edges are not oriented in such a way that a bounding box would be the minimal bounding box; therefore, before all edges can continue, any edge with a curvature of ±a threshold are further inspected.

The method of inspection is to extract the suspect edge into an image, and simply through a method of successive approximations; the edge is rotated and bounded, repeatedly, in the direction of the greatest decrease in the area of the bounding box, until a minimum acceptable tolerance is achieved. This method is highly effective in discriminating straight lines as compared to “S” shaped edges. For the algorithm, the minimum acceptable ratio used is 1:4, found by observing the “S” shaped edge with the smallest ratio, which was approximately 0.27. Note that the first few steps in the marching algorithm along an edge are weighted very little in calculating the curvature, as the average trajectory along those points is still highly variable, due to the discrete nature of the image and method of calculating a trajectory.

A.4 Method of Calculating Complete Intersection

Initially, in the development of this algorithm, the method of calculating the intersection between the filling (the thresholded image) and the edge was simply to check the direct overlap between them; however, because a single parameter controls the percent intersection limit on all pores, the challenge was that some edges suffered substantially far more intersection early in the filling process due to protruding edges, which spiraled into the region that should be filled. This is usually caused by pores within pores, as the Canny edge-detection algorithm tends to connect them into a single edge.

A method was developed to treat edges with and without internally branching portions the same, by requiring that a “complete intersection” occurs before being included. Essentially, instead of counting the number of edge pixels, which are overlapped by the filling, the edge image is broken up into rows and columns, where instead, the algorithm counts the total number of rows and columns experiencing all of their corresponding pixels in the intersection.

For example, in a rectangular edge, there would be two pixels contained in each row and two pixels in each column, except at the perimeters. If the rectangle has dimensions of 20-by-30 pixels, then there would be 50 rows and columns with the potential of intersection. The only way for any row or column to be counted as completely intersected is if all the pixels within them are overlapped. Therefore, with the actual overflow parameter used being around 50%, the algorithm will continue to increase the filling threshold until 50% of the rows and columns of the rectangular edge experience a complete intersection, requiring both parallel sides of the rectangle to be intersecting. The outcome of this method of determining a complete intersection is that any stray edges spiraling into the rectangle example above have no impact on the outcome, as the intersection must take place end to end along a row or column, regardless of an internal edge.

A.5 Final Discrimination of Edges by Analysis of Fill

Lastly, before the porosity fillings are all compiled and passed on for analysis, there is a final quality control test to ensure that no edges corresponding to straight-edge surface features nor highlights can contaminate the data. Shown in Fig. 21 is an example of a pair of two edges that were labeled with no significant curvature due to a break (discontinuity) in the continuous Canny-detected edge right before the U-turn which would have skewed its curvature calculation to that of a highlight. In addition, the substantial bend resulted in passing the minimal dimension ratio test. Fortunately, there is an obvious feature of highlights, when imagining the filling process on the 3D intensity heightmap of Fig. 10, which is that they fill outside in.

Fig. 21
(a) A real image acquired from the surface of a bone scaffold, showing a highlight, (b) inaccurate highlight edge detection due to edge discontinuity, and (c) analysis of fill together with discrimination of edges, leading to accurate detection of the highlight’s edges
Fig. 21
(a) A real image acquired from the surface of a bone scaffold, showing a highlight, (b) inaccurate highlight edge detection due to edge discontinuity, and (c) analysis of fill together with discrimination of edges, leading to accurate detection of the highlight’s edges
Close modal

A final analysis is performed after the filling process on all fillings; if the percent collision with the outer boundary of the region corresponding to the edge exceeds a specific parameter (around 55%), the filling is flagged and discarded. This method of only accepting pores is the most effective barrier to preventing noise from affecting the final analysis. In fact, it is effective enough to accomplish most of the task of discarding highlights and straight edges on its own, with the previous checks only assisting to reduce computation time and to reduce noise passed on. Note that the fillings which intersect with the image perimeter are discarded, for they are most often just a result of indentations along the perimeter of the top surface of the pore, as there are no pores that could be on the perimeter of the sample.

References

1.
Giannoudis
,
P. V.
,
Faour
,
O.
,
Goff
,
T.
,
Kanakaris
,
N.
, and
Dimitriou
,
R.
,
2011
, “
Masquelet Technique for the Treatment of Bone Defects: Tips-Tricks and Future Directions
,”
Injury
,
42
(
6
), pp.
591
598
.
2.
Pedrero
,
S. G.
,
Llamas-Sillero
,
P.
, and
Serrano-López
,
J.
,
2021
, “
A Multidisciplinary Journey Towards Bone Tissue Engineering
,”
Materials
,
14
(
17
), p.
4896
.
3.
Sun
,
Z.
,
Wu
,
F.
,
Gao
,
H.
,
Cui
,
K.
,
Xian
,
M.
,
Zhong
,
J.
,
Tian
,
Y.
,
Fan
,
S.
, and
Wu
,
G.
,
2020
, “
A Dexamethasone-Eluting Porous Scaffold for Bone Regeneration Fabricated by Selective Laser Sintering
,”
ACS Appl. Bio Mater.
,
3
(
12
), pp.
8739
8747
.
4.
Thadavirul
,
N.
,
Pavasant
,
P.
, and
Supaphol
,
P.
,
2014
, “
Development of Polycaprolactone Porous Scaffolds by Combining Solvent Casting, Particulate Leaching, and Polymer Leaching Techniques for Bone Tissue Engineering
,”
J. Biomed. Mater. Res. Part A
,
102
(
10
), pp.
3379
3392
.
5.
Yu
,
M.
,
Yeow
,
Y. J.
,
Lawrence
,
L.
,
Claudio
,
P. P.
,
Day
,
J. B.
, and
Salary
,
R.
,
2021
, “
Characterization of the Functional Properties of PCL Bone Scaffolds Fabricated Using Pneumatic Microextrusion
,”
ASME J. Micro- Nano-Manuf.
,
9
(
3
), p.
030905
.
6.
Murphy
,
C. M.
, and
O’Brien
,
F. J.
,
2010
, “
Understanding the Effect of Mean Pore Size on Cell Activity in Collagen-Glycosaminoglycan Scaffolds
,”
Cell Adhes. Migr.
,
4
(
3
), pp.
377
381
.
7.
Karageorgiou
,
V.
, and
Kaplan
,
D.
,
2005
, “
Porosity of 3D Biomaterial Scaffolds and Osteogenesis
,”
Biomaterials
,
26
(
27
), pp.
5474
5491
.
8.
O’Brien
,
F. J.
,
Harley
,
B.
,
Yannas
,
I. V.
, and
Gibson
,
L. J.
,
2005
, “
The Effect of Pore Size on Cell Adhesion in Collagen-GAG Scaffolds
,”
Biomaterials
,
26
(
4
), pp.
433
441
.
9.
Baldwin
,
C. A.
,
Sederman
,
A. J.
,
Mantle
,
M. D.
,
Alexander
,
P.
, and
Gladden
,
L. F.
,
1996
, “
Determination and Characterization of the Structure of a Pore Space From 3D Volume Images
,”
J. Colloid Interface Sci.
,
181
(
1
), pp.
79
92
.
10.
e Silva
,
R.
,
Negri
,
R.
, and
de Mattos Vidal
,
D.
,
2019
, “
A New Image-Based Technique for Measuring Pore Size Distribution of Nonwoven Geotextiles
,”
Geosynth. Int.
,
26
(
3
), pp.
261
272
.
11.
Chukwuma
,
K.
,
Bordy
,
E. M.
, and
Coetzer
,
A.
,
2018
, “
Evolution of Porosity and Pore Geometry in the Permian Whitehill Formation of South Africa–A FE-SEM Image Analysis Study
,”
Mar. Pet. Geol.
,
91
, pp.
262
278
.
12.
Chen
,
Z.
,
Song
,
Y.
,
Jiang
,
Z.
,
Liu
,
S.
,
Li
,
Z.
,
Shi
,
D.
,
Yang
,
W.
,
Yang
,
Y.
,
Song
,
J.
, and
Gao
,
F.
,
2019
, “
Identification of Organic Matter Components and Organic Pore Characteristics of Marine Shale: A Case Study of Wufeng-Longmaxi Shale in Southern Sichuan Basin, China
,”
Mar. Pet. Geol.
,
109
, pp.
56
69
.
13.
Kim
,
F.
,
Moylan
,
S.
,
Garboczi
,
E.
, and
Slotwinski
,
J.
,
2017
, “
Investigation of Pore Structure in Cobalt Chrome Additively Manufactured Parts Using X-Ray Computed Tomography and Three-Dimensional Image Analysis
,”
Addit. Manuf.
,
17
, pp.
23
38
.
14.
Du Plessis
,
A.
,
Sperling
,
P.
,
Beerlink
,
A.
,
Tshabalala
,
L.
,
Hoosain
,
S.
,
Mathe
,
N.
, and
Le Roux
,
S. G.
,
2018
, “
Standard Method for MicroCT-Based Additive Manufacturing Quality Control 1: Porosity Analysis
,”
MethodsX
,
5
, pp.
1102
1110
.
15.
Léonard
,
F.
,
Tammas-Williams
,
S.
, and
Todd
,
I.
,
2016
, “
CT for Additive Manufacturing Process Characterisation: Assessment of Melt Strategies on Defect Population
,”
Proceedings of the 6th Conference on Industrial Computed Tomography (ICT 2016)
,
Wels, Austria
,
Feb. 9–12
, pp.
1
8
.
16.
Villarraga
,
H.
,
Lee
,
C.
,
Corbett
,
T.
,
Tarbutton
,
J. A.
, and
Smith
,
S. T.
,
2015
, “
Assessing Additive Manufacturing Processes with X-Ray CT Metrology
,”
Proceedings of the ASPE Spring Topical Meeting: Achieving Precision Tolerances in Additive Manufacturing
,
North Carolina State University
,
Raleigh
,
Apr. 26–29
,
ASPE
, pp.
116
121
.
17.
Holzmond
,
O.
, and
Li
,
X.
,
2017
, “
In Situ Real Time Defect Detection of 3D Printed Parts
,”
Addit. Manuf.
,
17
, pp.
135
142
.
18.
Aminzadeh
,
M.
,
2016
, “
A Machine Vision System for In-Situ Quality Inspection in Metal Powder-Bed Additive Manufacturing
,”
Doctoral Dissertation
,
Advisor: Dr. Thomas Kurfess, Georgia Institute of Technology
,
Atlanta, GA
. http://hdl.handle.net/1853/56291
19.
Hu
,
D.
, and
Kovacevic
,
R.
,
2003
, “
Sensing, Modeling and Control for Laser-Based Additive Manufacturing
,”
Int. J. Mach. Tools Manuf.
,
43
(
1
), pp.
51
60
.
20.
Krauss
,
H.
,
Eschey
,
C.
, and
Zaeh
,
M.
,
2012
, “
Thermography for Monitoring the Selective Laser Melting Process
,”
Proceedings of the 2012 International Solid Freeform Fabrication Symposium
, The University of Texas at Austin, TX
,
Aug. 22
,
University of Texas Libraries
, pp.
999
1014
.
21.
Craeghs
,
T.
,
Clijsters
,
S.
,
Kruth
,
J.-P.
,
Bechmann
,
F.
, and
Ebert
,
M.-C.
,
2012
, “
Detection of Process Failures in Layerwise Laser Melting With Optical Process Monitoring
,”
Phys. Procedia
,
39
, pp.
753
759
.
22.
Bowoto
,
O. K.
,
Oladapo
,
B. I.
,
Zahedi
,
S.
,
Omigbodun
,
F. T.
, and
Emenuvwe
,
O. P.
,
2020
, “
Analytical Modelling of In Situ Layer-Wise Defect Detection in 3D-Printed Parts: Additive Manufacturing
,”
Int. J. Adv. Manuf. Technol.
,
111
(
7
), pp.
2311
2321
.
23.
Beucher
,
S.
,
1979
, “
Use of Watersheds in Contour Detection
,”
Proceedings of the International Workshop on Image Processing: Real Time Edge and Motion Detection/Estimation
,
Rennes, France
,
Sept. 17–21
,
CCETT
, pp.
1
12
.
24.
Otsu
,
N.
,
1979
, “
“A Threshold Selection Method From Gray-Level Histograms
,”
IEEE Trans. Syst. Man Cybern.
,
9
(
1
), pp.
62
66
.
25.
Bradley
,
D.
, and
Roth
,
G.
,
2007
, “
Adaptive Thresholding Using the Integral Image
,”
J. Graph. Tools
,
12
(
2
), pp.
13
21
.
26.
Yu
,
M.
,
Yeow
,
Y. J.
,
Lawrence
,
L.
,
Claudio
,
P. P.
,
Day
,
J. B.
, and
Salary
,
R. R.
,
2020
, “
Investigation of the Effects of Design and Process Parameters on the Mechanical Properties of Biodegradable Bone Scaffolds, Fabricated Using Pneumatic Microextrusion Process
,”
ASME 2020 International Mechanical Engineering Congress and Exposition (IMECE2020)
,
Portland, OR
,
Nov. 16–19
,
American Society of Mechanical Engineers
, Paper No. 24252.
27.
Zhao
,
D.
,
Yu
,
M.
,
Lawrence
,
L.
,
Claudio
,
P. P.
,
Day
,
J. B.
, and
Salary
,
R. R.
,
2020
, “
Investigation of the Influence of Consequential Design Parameters on the Mechanical Performance of Biodegradable Bone Scaffolds, Fabricated Using Pneumatic Micro-extrusion Additive Manufacturing Process
,”
ASME 2020 International Manufacturing Science and Engineering Conference (MSEC 2020)
,
Cincinnati, OH
,
June 22–26,
American Society of Mechanical Engineers,
Paper No. 8512.
28.
Yeow
,
Y. J.
,
Yu
,
M.
,
Day
,
J. B.
, and
Salary
,
R. R.
,
2020
, “
A Computational Fluid Dynamics (CFD) Study of Material Flow in Pneumatic MicroExtrusion (PME) Additive Manufacturing Process
,”
ASME 2020 International Mechanical Engineering Congress and Exposition (IMECE2020)
,
Portland, OR
,
Nov. 16–19
,
American Society of Mechanical Engineers
, Paper No. 24325.
29.
Salary
,
R.
,
Lombardi
,
J. P.
,
Weerawarne
,
D. L.
,
Rao
,
P. K.
, and
Poliks
,
M. D.
,
2018
, “
A Computational Fluid Dynamics (CFD) Study of Material Transport and Deposition in Aerosol Jet Printing (AJP) Process
,”
ASME 2018 International Mechanical Engineering Congress & Exposition (IMECE 2018)
,
Pittsburgh, PA
,
Nov. 9–15
,
American Society of Mechanical Engineers
, Paper No. 87647.
30.
Salary
,
R.
,
Lombardi
,
J. P.
,
Weerawarne
,
D. L.
,
Rao
,
P.
, and
Poliks
,
M. D.
,
2021
, “
A Computational Fluid Dynamics Investigation of Pneumatic Atomization, Aerosol Transport, and Deposition in Aerosol Jet Printing Process
,”
ASME J. Micro- Nano-Manuf.
,
9
(
1
), p.
010903
.
31.
Salary
,
R.
,
Lombardi
,
J. P.
,
Tootooni
,
M. S.
,
Donovan
,
R.
,
Rao
,
P. K.
,
Borgesen
,
P.
, and
Poliks
,
M. D.
,
2017
, “
Computational Fluid Dynamics Modeling and Online Monitoring of Aerosol Jet Printing Process
,”
ASME J. Manuf. Sci. Eng.
,
139
(
2
), p.
021015
.
32.
Spoerk
,
M.
,
Gonzalez-Gutierrez
,
J.
,
Sapkota
,
J.
,
Schuschnigg
,
S.
, and
Holzer
,
C.
,
2018
, “
Effect of the Printing Bed Temperature on the Adhesion of Parts Produced by Fused Filament Fabrication
,”
Plast., Rubber Compos.
,
47
(
1
), pp.
17
24
.
33.
Jansen
,
E. J.
,
Sladek
,
R. E.
,
Bahar
,
H.
,
Yaffe
,
A.
,
Gijbels
,
M. J.
,
Kuijer
,
R.
,
Bulstra
,
S. K.
,
Guldemond
,
N. A.
,
Binderman
,
I.
, and
Koole
,
L. H.
,
2005
, “
Hydrophobicity as a Design Criterion for Polymer Scaffolds in Bone Tissue Engineering
,”
Biomaterials
,
26
(
21
), pp.
4423
4431
.
34.
Jenkins
,
M.
, and
Stamboulis
,
A.
,
2012
,
Durability and Reliability of Medical Polymers
,
Woodhead Publishing
,
Philadelphia, PA
.
35.
Brakke
,
K.
,
2020
,
Triply Periodic Minimal Surfaces
,
Department of Mathematics and Computer Science, Susquehanna University
,
Selinsgrove, PA
.
36.
Schoen
,
A. H.
,
1970
, “
Infinite Periodic Minimal Surfaces Without Self-Intersections, National Aeronautics and Space Administration (NASA)
,” Technical Note (NASA TN D-5541),
Washington, DC
.
37.
Piper
,
S.
,
2018
,
Maths Models: Triply Periodic Minimal Surface Structures Mega Pack
,
MyMiniFactory
,
London, UK
.
38.
Klemstine
,
C.
,
Abdelgaber
,
Y.
,
Lawrence
,
L.
,
Day
,
J. B.
,
Claudio
,
P. P.
, and
Salary
,
R. R.
,
2021
, “
Characterization of the Compressive Properties of Triply Periodic Minimal Surface PCL Scaffolds for Bone Tissue Engineering
,”
ASME International Mechanical Engineering Congress & Exposition (IMECE 2021)
,
Virtual Conference
,
Nov. 1–5
,
American Society of Mechanical Engineers
, Paper No. IMECE2021-72125.
39.
Hoffman
,
D.
,
Hoffman
,
J.
,
Weber
,
M.
,
Trazet
,
M.
,
Wohlgemuth
,
M.
,
Boix
,
E.
,
Callahan
,
M.
,
Thayer
,
E.
, and
Wei
,
F.
,
2019
,
The Scientific Graphics Project—Triply Periodic Level Surfaces
,
Mathematical Sciences Research Institute
,
Berkeley, CA
.
40.
Chaffins
,
A.
,
Yu
,
M.
,
Claudio
,
P. P.
,
Day
,
J. B.
, and
Salary
,
R. R.
,
2021
, “
Investigation of the Functional Properties of Additively-Fabricated Triply Periodic Minimal Surface-Based Bone Scaffolds for the Treatment of Osseous Fractures
,”
ASME 2021 International Manufacturing Science and Engineering Conference (MSEC 2021), Virtual Conference (Due to the COVID-19 Pandemic), Hosted by the University of Cincinnati
,
Cincinnati, OH
,
June 21–25
,
American Society of Mechanical Engineers
, Paper No. 2004.
41.
Salary
,
R.
,
2022
,
Advanced Manufacturing for Bone Tissue Engineering and Regenerative Medicine
,
IntechOpen
,
London, UK
, p.
304
.