Automatic pore size measurements from scanning electron microscopy images of porous scaffolds

Pore sizes and distribution are amongst the main morphological characteristics of porous scaffolds which indicate the suitability of scaffolds for many biological applications. Scaffolds usually have complex structures and are designed to have a specific range of pore sizes appropriate for target cells. Pore sizes are commonly estimated manually or based on semi-automatic techniques requiring high level of human intervention. Such methods are time consuming and subject to error, mainly due to lack of consistency in the process and subjective nature of the results following operator involvement. In this work, we present a novel image processing method for the measurement pore size distribution (the main morphological characteristics of scaffolds) independent from their complexity. We use thresholding, based on the histogram analysis, to segment pore areas from scaffold, followed by morphological filters to separate pores from each other. This algorithm provides robust detection and measurement of pore sizes and the distribution. The performance of the algorithm is assessed using standard calibration kit which is used for calibration of Scanning Electron Microscopy (SEM) imaging systems. The results showed consistent output with 1.3% average error as compared against their true size. The algorithm was applied to 3D Apatite-Wollastonite scaffolds manufactured using the Thermally Induced Phase Separation technique. The results were robust and consistent with visual evaluation of SEM images. The algorithm also provides the morphology of each pore and, subsequently, offering further comprehension of the influence of microstructures across a range of fields, such as tissue engineering processes.


Introduction
One of the greatest challenges posed by researchers is to fabricate scaffolds which mimic the highly complex structure of natural tissue. Pore size and pore distribution are key factors in scaffold manufacture, as they directly and indirectly influence other factors including its interconnectivity, surface area, and mechanical characteristics [1][2][3]. Pore size is calculated as the area of void space between solid materials, and can be calculated in metric units using high resolution SEM images. An increase in the number and size of pores is usually associated with a decrease in mechanical stability and compressive strength. It is, therefore, essential that both the pore sizes and their distribution are adequately balanced for the desired outcome and also to ensure adequate mechanical stability within the scaffolds [2,4].
Scaffolds are usually designed to have a diverse range of pore sizes which are required for their targeted application and cell type [5,6]. The range of targeted pore sizes and their distribution play a vital role in providing an adequate supply of nutrition for cells, space for cell infiltration, and the required internal surface area for cell adhesion and proliferation. In effect, the pore size is an essential parameter in the formation of new tissue as pore size directly impacts on the spatial distribution, location and mass transfer of cells, nutrients, oxygen and waste metabolic products in a scaffold [7][8][9].
The estimation of pore size is highly dependent upon the method employed most of which are highly dependent on manual intervention [10,11]. The measurements are commonly performed on a number of pores in a couple of images, based on visual judgement, and extrapolated to represent the distribution of pores to the whole area. The length 1 3 of a pore is commonly measured by marking the coordinates of one side of a pore against the coordinates of another side, and computing their distance [2,12]. Doi et al. [10] have reported measuring pore diameters by placing a fixed circle over a pore, asserting that the diameter of the pore is that of the circle, as long as two furthest edges of the pore touches the edges of the circle. This technique is suitable for circular pores but inevitably will quantify pore sizes larger than the reality [10]. An example of the result of applying this technique is shown on one of our AW2 scaffolds in Fig. 1. Considering human error and dependency of these techniques to human interpretation, such techniques are subject to error. This will be more critical if measurement on a large number of images is the aim, mainly due to the time required to measure hundreds of pores in each image, which potentially causes high level of error due to the subsequent fatigue. Bartoš et al. [11] focused on the comparison of two widely used methods for the characterisation of porosity: scanning electron microscopy (SEM) and micro-computed tomography (micro-CT). Bartoš et al. [11] used the biggest inner circle diameter of pores and area-equivalent circle diameter to represent the pore size. Such measurement tends to underestimate the size of pores. They manually measured pore sizes of 40 pores in each SEM image of a scaffold [11]. The pore size dimensions were measured by means of ImageJ software [13][14][15][16].
A limited number of automatic image analysis algorithms have been developed for measuring pore sizes [17][18][19]. Elhadidy et al. [19] used thresholding, followed by watershed algorithm to segment pore regions and delete and presented only the pores that fit a conical shape. They tested the algorithm on atomic microscopy images of scaffolds and evaluated against manually delineated pores [19]. Additional research has identified the use of the maximum distance between two sides of pores, measured in x-axis and y-axis, in SEM images and reported the distances in two directions as the representation of pore size [20,21]. Applying such measurement techniques can lead to different pore sizes being reported, as measurement is dependent on the orientation of scaffold versus camera axis, and therefore such measurement techniques can only be suitable for scaffolds with circular pores. Techniques based on template matching are also proposed for measuring pore size distribution. Kim et al. [18] proposed an entirely new approach: using cross correlations with kernels of different shapes and sizes. They used rectangular shapes and oval shapes of different sizes and orientation to fit the pores and represent pore size based on the best fit shapes [18]. However, with 3D scaffolds characterised by non-circular and irregular pores, no assumption can be made about their shapes and sizes and the pores cannot fit into any template shapes, therefore such a technique is not suitable.
Leal-Egaña et al. [17] suggested an algorithm based on image segmentation and skeletonisation of the transmission electron microscopy (TEM) images to measure pore radii. The technique obtains information about the radial deformation of the hydrogel matrix due to cell growth, and the consequent modification of the pore size distribution pattern surrounding the cells [17]. The technique needs further evaluation. Tomlins et al. [6] used scanning confocal microscopy to measure pore size, they studied measuring pores using different imaging techniques [6]. Stachewicz et al. [22] studied the effect of pore size and distribution on cell growth into electrospun fibre scaffolds. Sizes of pores and fibres, pores shapes and porosity, before and after cell culture, were verified by imaging with scanning electron microscopy (SEM) and combination of focus ion beam (FIB) and SEM to obtain 3D reconstructions of samples [22]. This paper presents a novel automated method for the detection and accurate measurement of pores. The first step of algorithm uses a thresholding technique based on histogram analysis to segment SEM images into pore and non-pore areas. A set of morphological filters is then applied to separate connected pores from each other to allow measurement of individual pore size. The size of each is then measured in pixel size and converted to metric unit for further analysis. This algorithm does not make any assumptions about the shape or size of the pores on the scaffold. ImageJ is the primary software used for this study and analysis is performed on SEM images of 3D scaffolds. ImageJ is a public domain image-processing tool developed by the National Institutes of Health (NIH) for simple image processing applications (ImageJ, 2018).

Image processing algorithm
An image processing algorithm was developed to measure the pore sizes from the SEM images. The proposed algorithm is presented in five main steps as shown in flowchart, Fig. 2, (1) pixel size unit conversion to metric units; (2) image segmentation; (3) morphological filter to split connected pores; (4) verification step; (5) extract statistics from the detected pores.

Pixel size unit conversion
The size of pixels in SEM images depends on the magnification selected during imaging. It is necessary to convert the pixel to metric unit (ie. µm) in order to obtain the final measurements in unified scale, independently from the scale factor of imaging. To convert the image sizes to micrometres (µm), we performed the following sequence of steps in ImageJ: • Click on Line Tool • Click across the scale bar on the image so the bar is covered by the line • Click Analyse > Set Scale in order to select the scale bar in pixels • Click on Known distance and input the value of the original scale bar and the units to the scale units (e.g. 100 and μm as in Fig. 3)

Image segmentation
Image segmentation was performed using an automatic global (histogram-derived) thresholding technique described in Gonzales and Woods [16,23]. The algorithm finds two distinct peaks in the histogram of the image representing the scaffold material (bright parts of the image) and pores (dark parts of the image) respectively and detects the lowest intensity between these two peaks. The segmentation step generates a binary image with two categories: scaffold material (AW2) and void spaces (pores), as shown in Fig. 3.
The assumption is based on the fact that the histogram of a scaffold image usually has two distinct peaks representing scaffold material (bright parts of image) and pores (dark parts of image). Figure 3 shows a sample SEM image (a), its histogram with position of the threshold (b), and the segmentation results (c). Figure 3 (b) shows the threshold is located close to the lowest point between the two peaks in the histogram. ImageJ plugin allows adjustment of the threshold to optimise the output, and this can be done within the verification step. The output of the segmentation step is a binary image with scaffolds having value 1 and pore areas value 0. In ImageJ, after loading the image, the segmentation can be performed by the following stages: • Click on Image > Adjust > Threshold. The segmented regions represent the scaffold pores.
Then it is necessary to change the segmented image to binary image by following these steps: • Click on Process > Binary > Make Binary. This conversion changes the intensity values 0 to 255 to black and white parts of image with only two values 0 and 1, respectively.

Morphological filtering
This step allows separation of pores that are connected to each other in order to detect and measure them individually. The separation is performed by applying morphological dilation followed by erosion, two fundamental image processing operations in mathematical morphology [23]. A structuring element is used to perform the operation. The dilation of image I, using structuring element B, is defined by operator ⊕ [23].
Image O is the output image. If B is a flat structuring element, the dilation adds one pixel to the boundaries of objects in the binary image resulting in expanding the scaffold parts making the pores smaller in size. Repeating dilation will result in disconnecting pores from each other at the smallest junctions where pores are joined. To expand the pores toward the original size, erosion is applied. Erosion removes pixels on object boundaries and results in enlargement of pore areas. The erosion of image I using structuring element B is defined by operator ⊖ .
If B is the same structuring element used for dilation, the size of pore expands/recovers by the same amount. However, small connections between the pores which were connected in the dilation step will not be reconnected, which results in separation of pores so they can be individually evaluated in the next steps.
For these morphological operations, a flat structure element of small size is scanned over the image and every point in the centre of structure element is replaced with one of the pixel values (the lowest or highest pixel value) under the structure element depending on the operation (erosion or dilation). Figure 4 shows a segmented SEM image following twice dilating the image using a structure element of 3 × 3 pixels and twice eroding the image using the same structuring element. In ImageJ, the dilation and erosion can be performed using: • For Dilation; click on Process > Binary > Dilate • For Erosion; click on Process > Binary > Erode

Verification
A verification step is included in the image processing algorithm to ensure that the threshold selected in the earlier step leads to a correct segmentation and separation of pores after morphological filtering. This is performed by analysing if the boundaries of the segmented regions are aligned and if the pores which are inter-connected (connected with small gaps) in the original image are separated from each other. If the result is not satisfying any of the two criteria, the image segmentation and morphological filter steps can be repeated with new threshold using a simple rule. Particularly, if pore sizes are smaller than expected, it is necessary to move the threshold to lower values (left cursor of the middle image of Fig. 3), while if too large, it is necessary to move the threshold to higher values (right cursor of the middle image of Fig. 3). If pores are not separated, it is necessary to apply a greater number of dilations and the same number of erosions. The repeat of the process may be needed in rare occasions but was not required in this study.

Extract statistics from the detected pores
The areas in black shown in Fig. 4b, c, d represents the void spaces (pores area) within the scaffold material. Using Analyse Particles tool of ImageJ, the size of each pore can be measured and reported with an index representing the pore number in the image. The steps to extract and export pore sizes using ImageJ are as below: • Click on Analyse Particles > Show-Outlines (to see outlines) • Click on Analyse the particles (to extract the area of each segmented region) The area of each pore with its index is reported and can be simply saved in an Excel sheet in metric units (µm 2 in our application).

Validation technique
Validation of the image processing algorithm was necessary to establish the performance of the algorithm. The validation is performed quantitively by comparing the results of the algorithm on a known input and measuring the error. The error is measured as the percentage of incorrectly detected pixels in a pore versus the total number of pixels in the pore. The robustness of the algorithm can be established by evaluating the standard deviation of error over large varieties of known areas/pixels [23].
A reliable reference for evaluating the results of the algorithm is to use calibration standard kit of the SEM imaging instrument, and comparing the results of the algorithm against the size of standard calibration circles on the kit. This process was repeated five times to generate the standard deviations of the error. The SEM images of a standard calibration kit (Fig. 5) were analysed by the algorithm to measure their sizes, where the error of algorithm (%) can be measured using Eq. (3).
The real size of objects i are provided by the manufacturer.
In brief, PLA was dissolved in Dioxane (99.8% purity, Sigma-Aldrich) to obtain a 5% w/v solution at 85 °C for 3 h using a hot plate with magnetic stirrer (IKA C-MAG HS7). Then, 50 g of AW2 powder (20-50 µm fractions) were added to the solution at 40 °C, mixing for 2 h. The solution was poured into ice-cold glass moulds and then frozen (− 20 °C) overnight. The moulds were repeatedly washed with cold (− 20 °C) 80% Ethanol (200 proof, Sigma-Aldrich) over four consecutive days. The scaffolds were removed from the moulds and then freeze-dried for 48-h (ALPHA 1-2 LD plus; Martin Christ Freeze Dryers, UK). Finally, the freezedried scaffolds underwent specific heat treatment cycles in the furnace (SNOL10/1300LHM01 Eurotherm 3216), consisting of a heating step to 779 °C for 1 h at a heat rate of 5 °C per minute, followed by a heat rate of 10 °C per minute until the required sintering temperature was reached, then held at that temperature for 1 h. The sintering temperature ranged from 1230 °C to 1260 °C in 5 °C increments.

SEM imaging protocol
Scanning Electron Microscope (SEM) (Hitachi TM3030) was used for high resolution imaging of the scaffolds. Each sample was initially cleaned with pressurised air to remove any floating elements. Aluminium stubs were coated in double-sided Carbon tape and each scaffold was then gently loaded onto the SEM set up. A digital photo was taken of the loaded sample stage prior to insertion into the SEM device. Five areas of view were imaged from each scaffold (horizontally and vertically) at different positions. The SEM images were captured at 15 kV and at 800 magnification factors (saved as JPG or TIF file format).

Evaluation of SEM image analysis algorithm using standard template
The error was calculated for all measured circles in the SEM images of the calibration kit and presented as percentage of error using formula presented in Eq. (3). The calibration template consisted of six circles of varying sizes: 100, 200, 300, 500, 600 and 1000 µm diameter, all were used in this experiment. The circles were measured using the algorithm and then compared against the real area of each circle as supplied by the manufacturer. Figure 5 shows the original SEM image of the standard calibration sample and Fig. 6a shows a sample result using the algorithm. The results showed the percentage error of the detection algorithm against the real dimensions on average is-1.36% with standard deviation (STD) of 2.29%. The size of the circles detected by the algorithm was on average 1.36% smaller than their real size. This error in measuring the size is partially due to the digitisation effect and the low resolution of the SEM imaging system used, which was 4.4 µm per pixel. For the smaller circle with diameter 100 µm, the error due to the resolution can be up to 4.4%, however this margin of error is smaller for larger circles.
Any imaging system can have an error of one pixel at the boundary of the objects in the image, consequently this will result in an error in the measured size of circles from the real size (in our case, smaller than their true size). Such errors can also result in slight change in morphology of the pore structure being measured. To explain this, consider the outlines of two calibration circles (3) and (4) as shown in. Figure 6b. The circle number 3 has a symmetrical boundary as opposed to circle number 4 which has an asymmetrical boundary. Circle 3 has symmetrical outline (Fig. 6b) and  is measured more accurately than circle number 4 which has asymmetrical outline. This is due to the effect of pixelation on SEM images which results in an error of about one-pixel thickness in the algorithm output.

In-house 3D porous AW2 scaffolds manufactured by TIPS
Photographic images of 3D AW2 porous scaffolds, fabricated by TIPS, are shown in Fig. 7. The figure shows a set of 6 scaffolds after heat treatment on platinum foil, and an internal image of a scaffold broken on a fracture line. Figure 8 shows the SEM images of AW2 scaffolds fabricated at eight different sintering temperatures (1230 °C up to 1260 °C in 5 °C increments). There is a clear increase in the necking between AW2 particles as the temperature increases.
The SEM images were processed using the algorithm described in the method section and pore sizes were measured. Automatic thresholding was used and twice dilation followed by twice erosion was applied on all SEM images. Figure 9 shows the results of different steps of the algorithm, red arrow shows that the algorithms separated two connected pores allowing measurement of the size of each pore independently. In all cases, the quality of segmented images was suitable, as evaluated at the verification step, and there was no need to adjust or reapply the thresholding or morphology filters.
The size (in micrometres squared) and the number of pores for different scaffolds are presented as a histogram in Fig. 10a, by considering 5 different pore size ranges. It was observed that our image-processing algorithm detected very small pores, smaller than 0.5 μm 2 , to very large pores, larger than 1000 μm 2 . The figure shows that the highest number of large pores are present at 1245 °C. Processing time for one single image typically takes less than a minute.  . 10 a The distribution of pores in scaffolds fabricated using different sintering temperatures ( °C). Frequency of pores for different pore size ranges for a unit area of 300 µm × 300 µm. Brackets () is used for inclusive values and [] is used for exclusive values. b Aver-age pore area for large pores (pores larger than 5 µm 2 ) for ten different sintering temperatures from 1200 to 1260 °C. The average pore size is measured for a unit area. Error bars represent the standard error of the mean