System Configuration

The system consists of a consumer-grade 3D depth sensor (Xtion Pro Live, ASUSTeK Computer Inc. Taipei, Republic of China) and a laptop computer (Core-i7, 7200U-16 GB CF-SZ5, Panasonic Inc, Tokyo, Japan). The depth sensor has depth measurement range of 0.8–3.5 m; horizontal, vertical, and diagonal view angles of 58°, 45°, and 70°; depth image sizes of 320 × 240 pixels (QVGA) and 640 × 480 pixels (VGA); and space sampling resolution at horizontal and vertical distances of 800 mm, 2.8 mm (QVGA), and 1.4 mm (VGA), and diagonal distances of 2.2 mm (QVGA and VGA). This system can be immediately used following installation of the external 3D deep sensor, which can be set up within a few minutes (Fig. 1). In addition, the sensor needs only a single calibration.

Figure 1 Experimental device and algorithm for the detection of 3D asymmetries. The device consists of a 3D scanner and a laptop computer. A subject bends forward and the surface of the back is captured by the scanner. Full size image

Asymmetry Analysis Algorithm

An overview of our fully automated algorithm for detecting 3D asymmetry is shown in Fig. 1. Asymmetry analysis is performed and an asymmetry index I asym is calculated as follows:

Stage 1: Capture point clouds of a patient’s back

The surface of the patient’s back is scanned by the depth sensor so that the transverse lines of right and left of the waist are roughly aligned with the recommended lines displayed on the computer monitor. A 3D point cloud \({{\boldsymbol{P}}}^{0}(\,=\,\{{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{0}\})\) is captured that includes the back, breech, neck, and occipital surfaces and a floor, where \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{0}\) denotes the 3D coordinates of a point included in point cloud P0.

The original points in P0 from the sensor typically include a non-negligible amount of errors and outliers. To reduce them, a smoothing process with a 3D moving average filter24 is first applied to P0 to obtain the smoothed point cloud P1, as shown in Supplementary Fig. 1A.

Stage 2: Segment the point cloud

Since the original point cloud P1 still includes unnecessary points from the floor, P1 is first partitioned into several connected point regions using a smoothness-constraint-based region growing algorithm25. As shown in Supplementary Fig. 1A, a point cloud P2 that includes only the back, breech, neck, and occipital surfaces is extracted by taking the region closest to the scanner for which the sum of the side lengths of the region bounding box exceeds 650 mm.

Stage 3: Estimate a median sagittal plane

As shown in Supplementary Fig. 1B, principal component analysis is applied to P2, and the pose-normalized point cloud P3 is obtained from P2 where the approximate median sagittal plane is aligned along the x and z axes of P3.

Stage 4: Estimate the boundaries of the back surface

P3 still includes points on the breech, neck, and occipital surfaces—which are unnecessary for asymmetry analysis of the back surface and need to be removed. To this end, as shown in Supplementary Fig. 2A, P3 is first projected to its xy plane to generate a binary image \({{\boldsymbol{I}}}^{3}=\) \(\{{\boldsymbol{b}}({\boldsymbol{i}},{\boldsymbol{j}})\,|{\boldsymbol{i}}=\frac{{\boldsymbol{x}}-{{\boldsymbol{x}}}_{{\boldsymbol{\min }}}}{{\rm{\Delta }}{\boldsymbol{x}}}\in \{0,{{\boldsymbol{i}}}_{{\boldsymbol{\max }}}\},{\boldsymbol{j}}=\frac{{\boldsymbol{y}}-{{\boldsymbol{y}}}_{{\boldsymbol{\min }}}}{{\rm{\Delta }}{\boldsymbol{y}}}\in \{0,{{\boldsymbol{j}}}_{{\boldsymbol{\max }}}\},{\boldsymbol{b}}\in \{0,1\}\}\) where a black pixel (b(i, j) = 1) includes at least one projected point, and a white pixel (b(i, j) = 0) does not. Δx and Δy denote pixel resolutions of the x and y direction, and x min , x max , y min and y max denote the ranges of x, y coordinates of the points in P3. \(\lfloor x\rfloor \) denotes a floor function that gives the largest integer that does not exceed x. Then, the width of the black pixel region along the y direction at each discrete x position in I3 is evaluated as the discrete width function

$${{\boldsymbol{w}}}_{{\boldsymbol{d}}}({\boldsymbol{x}})=\sum _{{\boldsymbol{j}}\in \{0,{{\boldsymbol{j}}}_{{\boldsymbol{\max }}}\}}\,{\boldsymbol{b}}(\lfloor ({\boldsymbol{x}}-{{\boldsymbol{x}}}_{{\boldsymbol{\min }}})/{\rm{\Delta }}{\boldsymbol{x}}\rfloor ,{\boldsymbol{j}}).$$

Finally, a high-order polynomial w y (x) is fit to w d (x) between x min and x max . We adopted the 10th-order polynomial as w y (x), which gives appropriate results. Next, by taking some local minima of polynomial w y (x), the boundaries between the back and neck and the back and breech can be identified, and a rectangular region I4 is extracted from I3 that excludes the breech, neck, and occipital surfaces and is suitable for asymmetry analysis (Supplementary Fig. 2B). A subset of points of P3 contained only in I4 are extracted as a point cloud for asymmetry analysis of P4, as shown in Supplementary Fig. 2C.

Stage 5: Generate a reflected point cloud

To carry out the asymmetry analysis, a reflected point cloud P4r is first generated by taking a mirror projection of P4 with respect to the approximated median sagittal (xz) plane (Supplementary Fig. 2C). Since the points close to the outer boundaries of P4 and P4r may contain larger measurement errors than those in the center, width-restricted point clouds P5 and P5r are extracted whose ranges in y direction is limited to double-sided 90% range of those of P4 and P4r. Width-restricted point cloud P5 and its reflected version P5r are then used to find the iterative closest point (ICP)-based best fit.

Stage 6: Find the best fit between the two point clouds

The point-to-plane ICP algorithm26 is carried out to obtain an optimum position and orientation \(\langle R,t\rangle \) for P5r best-fitted to P5. As shown in Supplementary Fig. 3A, in ICP, the optimum R, t is derived as a solution of a minimization problem defined in Eq. (1):

$$\mathop{\mathrm{minimize}\,}\limits_{\langle {\boldsymbol{R}},{\boldsymbol{t}}\rangle }\frac{1}{|\{({{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5},{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}})\}|}\sum _{{{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5}\in {{\boldsymbol{P}}}^{5},\,{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}}\in {{\boldsymbol{P}}}^{5{\boldsymbol{r}}}}{|{{\boldsymbol{n}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5}\cdot ({\boldsymbol{R}}{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}}+{\boldsymbol{t}}-{{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5})|}^{2},$$ (1)

where \(({{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5},{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}})\) gives a valid pair of the closest points on P5 and P5r, R is a 3 × 3 rotation matrix, and t represents a 3D translation vector for transforming reflected point cloud P5r. \(|\{({{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5},\,{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}})\}|\) denotes the number of valid pairs of the closest points. A normal vector \({{\boldsymbol{n}}}_{{\boldsymbol{i}}}^{5}\) and a radius of curvature \({{\boldsymbol{\rho }}}_{{\boldsymbol{i}}}^{5}\)at every point \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5}\) in P5, and \({{\boldsymbol{n}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}}\) and \({{\boldsymbol{\rho }}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}}\)at \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}}\) in P5r are precomputed by principal component analysis using a set of neighboring points \({\boldsymbol{N}}({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5},{{\boldsymbol{r}}}^{5})\) and \({\boldsymbol{N}}({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}},{{\boldsymbol{r}}}^{5})\) respectively, where N(p, r) is a set of points included in a small sphere centered at p with a radius r. r5 = 25 mm has been shown to give adequate results.

To avoid to taking an improper (noncorresponding) point pair \(({{\boldsymbol{p}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5},{{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}})\) in evaluating Eq. (1), if the angle between the normal vectors at the closest point pairs \(({{\boldsymbol{n}}}_{{\boldsymbol{c}}({\boldsymbol{i}})}^{5},{{\boldsymbol{n}}}_{{\boldsymbol{i}}}^{5{\boldsymbol{r}}})\) exceeds the threshold (50°), such a point pair is treated as invalid and is ignored in minimizing Eq. (1).

Unfortunately, \(\langle R,t\rangle \) derived from Eq. (1) may change slightly depending on the initial orientation of P5r relative to P5, and this change may induce a non-negligible error in the asymmetry analysis. To avoid this, we perform small rotational perturbations from −5° to 5° at intervals of 1.25° about the x and y axes to the initial orientation of \({{\boldsymbol{P}}}_{{\boldsymbol{r}}}^{5}\) and then conduct ICP on Eq. (1) several times. In each optimization process from the initial orientation, we iterate the process 50 times. Finally, the best fitted \(\langle R,t\rangle \), which gives the smallest value of Eq. (1), is adopted as the optimum rotation and translation \(\langle {R}^{\ast },{t}^{\ast }\rangle \) for P5r, as shown in Supplementary Fig. 3B. By applying \(\langle {R}^{\ast },{t}^{\ast }\rangle \) to the original reflected point cloud P4r, the optimal reflected point cloud P6r—which is best fitted to the original point cloud P4—is finally generated (Supplementary Fig. 3C).

Stage 7: Roughly extract deviations for colormap rendering

The difference in the position of P6r relative to P4 is roughly evaluated, and the difference distribution is rendered as a color map as Supplementary Fig. 4. In the evaluation, a triangle mesh surface M6r is generated from point cloud P6r. Then, distance \({{\boldsymbol{d}}}_{{\boldsymbol{i}}}^{4}\) is evaluated from a point \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{4}\) in P4 to the point on M6r closest to \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{4}\) along the normal vector \({{\boldsymbol{n}}}_{{\boldsymbol{i}}}^{4}\) at \({{\boldsymbol{p}}}_{{\boldsymbol{i}}}^{4}\). If \({{\boldsymbol{d}}}_{{\boldsymbol{i}}}^{4}\) exceeds the threshold (100 mm), it is ignored in the evaluation. Finally, the distribution of the distances \(\{{{\boldsymbol{d}}}_{{\boldsymbol{i}}}^{4}\}\) is rendered as a colormap on P4 representing how much the geometry of the surface of the subject’s back deviates from its ideally symmetric shape, and where the deviations are most prominent.

Stage 8: Precisely extract deviations and estimate the asymmetry index

Finally, the asymmetry index is evaluated, which precisely indicates the extent to which the geometry of the surface of the patient’s back deviates from perfect symmety. For this purpose, the centroid and principal axes of the union of point clouds P6r and P4 are calculated (Supplementary Fig. 5A), and the points clouds P4 and P6r are transformed to P7 and P7r such that the principal axes are aligned to x, y, and z axes where the xy plane is aligned to the frontal plane and the x axis is oriented to the craniocaudal axis. Then, as shown in Supplementary Fig. 5B, two point clouds P7 and P7r are regularly resampled at two-dimensional rectangular grid points, each of which (l, m) is located on the xy plane and whose grid spacings are ∆x and ∆y in the x and y directions of the new axes. As a result, two resampled point clouds,

$$\begin{array}{c}{{\boldsymbol{P}}}^{{\bf{8}}}=\{{{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}}=({{\boldsymbol{p}}}_{{\boldsymbol{x}},\,{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}},\,{{\boldsymbol{p}}}_{{\boldsymbol{y}},{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}},{{\boldsymbol{p}}}_{{\boldsymbol{z}},{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}})|{\boldsymbol{l}}=\frac{{\boldsymbol{x}}}{{\rm{\Delta }}{\boldsymbol{x}}},{\boldsymbol{m}}=\frac{{\boldsymbol{y}}}{{\rm{\Delta }}{\boldsymbol{y}}}\}\,\,{\rm{and}}\\ {{\boldsymbol{P}}}^{{\bf{8}}{\boldsymbol{r}}}=\{{{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}{\boldsymbol{r}}}=({{\boldsymbol{p}}}_{{\boldsymbol{x}},{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}{\boldsymbol{r}}},{{\boldsymbol{p}}}_{{\boldsymbol{y}},{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}{\boldsymbol{r}}},{{\boldsymbol{p}}}_{{\boldsymbol{z}},{\boldsymbol{l}},{\boldsymbol{m}}}^{{\bf{8}}{\boldsymbol{r}}})|{\boldsymbol{l}}=\frac{{\boldsymbol{x}}}{{\rm{\Delta }}{\boldsymbol{x}}},{\boldsymbol{m}}=\frac{{\boldsymbol{y}}}{{\rm{\Delta }}{\boldsymbol{y}}}\},\end{array}$$

are generated from P7 and P7r, respectively. We adopted 3 mm as ∆x and ∆y based on the space sampling resolution of our depth sensor. Usually, resampled points \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) and \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8{\boldsymbol{r}}}\) do not exactly locate on points in P7 and P7r. Therefore, as shown in Supplementary Fig. 5C, we first find a point \({{\boldsymbol{p}}}_{{\boldsymbol{j}}}^{7}\)in P7 whose projection on the xy plane is closest to \(({{\boldsymbol{p}}}_{{\boldsymbol{x}},{\boldsymbol{l}},{\boldsymbol{m}}}^{8},{{\boldsymbol{p}}}_{{\boldsymbol{y}},{\boldsymbol{l}},{\boldsymbol{m}}}^{8},0)\), and then fit a small spherical surface S, which passes through \({{\boldsymbol{p}}}_{{\boldsymbol{j}}}^{7}\) and has a normal vector \({{\boldsymbol{n}}}_{{\boldsymbol{j}}}^{7}\) and radius of curvature \({{\boldsymbol{\rho }}}_{{\boldsymbol{j}}}^{7}\) at \({{\boldsymbol{p}}}_{{\boldsymbol{j}}}^{7}\), which was computed in Step 6. Finally, we find a resampled point \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) intersecting a vertical line passing \(({{\boldsymbol{p}}}_{{\boldsymbol{x}},{\boldsymbol{l}},{\boldsymbol{m}}}^{8},{{\boldsymbol{p}}}_{{\boldsymbol{y}},{\boldsymbol{l}},{\boldsymbol{m}}}^{8},0)\) and surface S to build P8. The other resampled point \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8{\boldsymbol{r}}}\) in P8r is calculated similarly.

An example of P7 and P8 is shown in Supplementary Fig. 6. If both resampled points \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) and \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8{\boldsymbol{r}}}\) lie inside the effective range of P8 and P8r at a grid point (l, m), and these points can be resampled correctly, then grid point (l, m) is treated as effective, and is included in effective grid point set G E .

As shown in Supplementary Fig. 7, a new normal vector \({{\boldsymbol{n}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) at \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) is estimated using the principal component analysis of the neighboring point set \({\boldsymbol{N}}({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8},{{\boldsymbol{r}}}^{8})\). Plane S′, which passes through \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) and has a normal \({{\boldsymbol{n}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\), is generated, and distance d l, m between plane S′ and \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8{\boldsymbol{r}}}\) is evaluated as a deviation at this resampled point. r8 = 25 mm gives adequate results. However, even if grid point (l, m) is included in effective grid point set G E , if (l, m) lies close to the boundary of P8, it is possible that resampled points \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8{\boldsymbol{r}}}\) or \({{\boldsymbol{p}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}^{8}\) still suffer from relatively large measurement error due to the high incident angle between the optical axis of the depth sensor and the surface normal. Therefore, we must evaluate deviation d l, m only at a valid grid point (l, m), which is included in G E but does not lie close to the boundary of P8. To do so, as shown in Supplementary Fig. 8, we adopt grid point (l, m)as valid only if it lies in a valid range located around the median sagittal plane such as M min (l) + 0.5(1 − α)W(l) ≤ m ≤ M max (l) − 0.5(1 − α)W(l), where M min (l) and M max (l) are minimum and maximum values of the y coordinate of the grid points at l included in G E , and W(l) = M max (l) − M min (l). The parameter α ∈ [0, 1] controls to what extent we adopt the grid points in G E as valid ones for the final asymmetry analysis. Based on this criterion, if a grid point (l, m) is classified as valid, then the point is included in valid grid point set G v .

Finally, the distribution of the sampled deviation between the point clouds of the surface of the patient’s back P8 and its best-fitted and reflected point cloud P8r at valid grid points is determined as \({\boldsymbol{D}}=\) \(\{{{\boldsymbol{d}}}_{{\boldsymbol{l}},{\boldsymbol{m}}}|{\boldsymbol{l}}=\lfloor \frac{{\boldsymbol{x}}}{{\rm{\Delta }}{\boldsymbol{x}}}\rfloor ,{\boldsymbol{m}}=\lfloor \frac{{\boldsymbol{y}}}{{\rm{\Delta }}{\boldsymbol{y}}}\rfloor ,({\boldsymbol{l}},{\boldsymbol{m}})\in {{\boldsymbol{G}}}_{{\boldsymbol{v}}}\}\). Based on distribution D, the asymmetry index I asym is evaluated as an averaged deviation over the grid points included in valid grid point set G v using Eq. (2), which precisely represents the degree of asymmetry quantitatively:

$${{\boldsymbol{I}}}_{{\boldsymbol{asym}}}=\frac{1}{|{{\boldsymbol{G}}}_{{\boldsymbol{v}}}|}\,\sum _{({\boldsymbol{i}},{\boldsymbol{j}})\in {{\boldsymbol{G}}}_{{\boldsymbol{v}}}}\,{{\boldsymbol{d}}}_{{\boldsymbol{i}},{\boldsymbol{j}}},$$ (2)

where |G v | denotes the number of valid grid points. If the shape of the back surface at G v is exactly symmetric with respect to the median sagittal plane, index I asym is zero; the more the back’s shape deviates from perfect symmetry, the greater I asym becomes.

Capturing Posture Images

For this fully automated algorithm for 3D asymmetry detection, the subject is asked to bend forward, a pose in which protrusions on the back are most easily observed and which is extensively used for the traditional scoliosis examination known as the Adams bending test15,16. To prevent cases in which the patient’s posture diverges significantly from the recommended pose, such as cases in which the subject’s back is not visible in the photographed data, the system implements the following countermeasures: (1) The positions in which subjects should place their feet are marked. In addition, a box that indicates the recommended positioning is displayed on the screen during photography. The subject can be positioned based on this frame and visually guided to approach the recommended position in the x (left–right) and y (craniocaudal) directions, and attain the proper z rotation angle. (2) The angular difference between the approximate plane of the 3D point cloud in the detection region of the photographed subject, and the direction of photography of the 3D sensor is calculated during photography. This value can be checked to assess whether it is within the threshold value; via analysis of the 3D data, the rotation angles along the x direction (<±7.5°), and the angular rotation along the y direction (<±15°) can be made to approach the recommended state.

Human Subjects

Informed consent was obtained from all the parents/legal guardians of paediatric patients and all work was performed in accordance with the guidelines of the institutional internal review board of the participating institution. The experimental protocols were approved by the institutional review board of Hokkaido University Hospital, and the registered number is 016–0236. The human subjects had been referred to our hospital on suspicion of IS and were asked to bend forward, as for the Adam’s forward bending test15,16, followed by scanning with the sensor. The inclusion criteria were as follows: (i) aged 7 to 18 years old; (ii) referred for confirmed diagnosis based on x-ray; (iii) willing and able to provide written informed consent. We explained the procedure, indications, and preparation, and participants indicated their understanding and signed the corresponding consent forms. We received the subjects’ medical records to comparing the accuracy of the developed technique with the current standard of care. The medical records included age, sex, body mass index (BMI), and Cobb angle measured from radiographs. Seventy-six 3D point clouds and Cobb angles were obtained for this study. The average age was 13.8 years (range, 7 to 18 years). The mean BMI and Cobb angle was 18.6 (range, 13.4 to 27.0) and 21.2° (range, 0° to 64°), respectively. When subjects with a Cobb angle of 0° were excluded, the mean Cobb angle of single thoracic curve (n = 33), double thoracic and thoracolumbar/lumbar curve (n = 19), and single thoracolumbar/lumbar (n = 17) was 22.7° (range, 4° to 64°), 29.1° (range, 18° to 43°) and 25.2° (range, 12° to 41°), and 18.2° (range, 8° to 35°), respectively. Our internal studies of intra- and inter-rater reliability have shown excellent kappa statistics27 for the Cobb angle measurements (0.95–0.99).

Repeatability Analysis and the Effects of Trunk Rotation

To assess the independent repeatability of and effects of trunk rotation on this system isolated from human subjects, phantom models were obtained from plaster wrap castings made when IS patients were fitted for hard braces to treat single thoracic curves, double thoracic and thoracolumbar/lumbar curves, and single thoracolumbar/lumbar curves, respectively. The phantom models were scanned 10 times from the front, and at 5° clockwise and counterclockwise rotations.

Intra-observer repeatability

To assess intra-observer repeatability, thirty subjects who agreed to the repeated test were scanned in the same position. The test-retest evaluation was performed in two acquisitions for the same subject, with repositioning between the acquisitions.

Effects of smoothing operations on the point clouds

We have assessed the precision of the method by comparing the asymmetry index between the point cloud captured from a phantom model’s surface by a high-definition non-contact 3D digitizer (Vivid 910, Konica Minolta Inc., Tokyo, Japan) and that by our consumer-grade sensor. The high-definition non-contact 3D digitizer has depth measurement range of 0.6–2.5 m; a horizontal and vertical scan range of 823 mm and 618 mm at 2.5 m depth, respectively. The measurement accuracy in the horizontal, vertical, and depth directions are ±0.38 mm, ±0.31 mm, and ±0.20 mm to the reference plane at the middle measurement mode. Around 98,000 dense points on the dummy surface were captured by the digitizer in a scan. Therefore, the accuracy of the non-contact 3D digitizer is much better than that of the 3D depth sensor, and the point clouds measured by the digitizer can be treated as reference data in the comparison. The phantom model was obtained from plaster casting when a patient was fitted for a hard brace (single thoracic curve, T5-T11,31°). We also evaluated the difference in the asymmetry index before and after our approximation, and smoothing operations for the point clouds measured both from the 3D digitizer and the consumer-grade sensor. The phantom model was scanned three times from the front.

Statistical analysis

Pearson’s correlation coefficient analyses were used to assess relationships between the asymmetry index and the Cobb angle. A Fisher r-to-z transformation was used to test for the difference between two correlation coefficients. Receiver operating characteristic (ROC) analyses were used to estimate the best cut-off values for the asymmetry index to predict Cobb angles greater than 15°, 20°, or 25°. The performance of the asymmetry indices were evaluated using the area under the ROC curve (AUC). AUCs were categorized as follows: no discrimination (AUC = 0.50), acceptable discrimination (0.7 ≤ AUC < 0.8), excellent discrimination (0.8 ≤ AUC < 0.9), and outstanding discrimination (AUC ≥ 0.9)28. To determine the cut-off values, we used the Youden index, which is calculated as sensitivity + specificity −1 for each potential cut-off point; the optimal cut-off point is the tool score with the highest value29. Based on the cut-off values, the sensitivities, specificities, positive predictive values, negative predictive values, accuracy, positive likelihood ratios, and negative likelihood ratios were determined. One-way ANOVA was used to compare differences across positions. Repeatability was determined in terms of the coefficient of variation (CV), which is a measure of variability relative to the mean and is calculated by dividing the standard deviations of intrasubject results by the mean and multiplying by 100 to present the result as a percentage30. A CV less than 10% was regarded as very good, 10% < CV ≤ 20% as good, 20% < CV ≤ 30% as fair/moderate, and CV > 30% as poor31. The intraclass correlation coefficient (ICC) was calculated to evaluate intra-observer repeatability. ICC values less than 0.5 are indicative of poor reliability, values between 0.5 and 0.75 indicate moderate reliability, values between 0.75 and 0.9 indicate good reliability, and values greater than 0.90 indicate excellent reliability32. Data analyses, except for ICC, were performed using JMP statistical software for Windows (version 12; SAS, Inc., Cary, NC, USA). The ICC was analysed using SPSS software for Windows (version 23; IBM, Inc., Chicago, IL, USA). P < 0.05 was considered statistically significant.