Accuracy assessment of Tri-plane B-mode ultrasound for non-invasive 3D kinematic analysis of knee joints
© Masum et al.; licensee BioMed Central Ltd. 2014
Received: 7 February 2014
Accepted: 1 August 2014
Published: 26 August 2014
Currently the clinical standard for measuring the motion of the bones in knee joints with sufficient precision involves implanting tantalum beads into the bones. These beads appear as high intensity features in radiographs and can be used for precise kinematic measurements. This procedure imposes a strong coupling between accuracy and invasiveness. In this paper, a tri-plane B-mode ultrasound (US) based non-invasive approach is proposed for use in kinematic analysis of knee joints in 3D space.
The 3D analysis is performed using image processing procedures on the 2D US slices. The novelty of the proposed procedure and its applicability to the unconstrained 3D kinematic analysis of knee joints is outlined. An error analysis for establishing the method’s feasibility is included for different artificial compositions of a knee joint phantom. Some in-vivo and in-vitro scans are presented to demonstrate that US scans reveal enough anatomical details, which further supports the experimental setup used using knee bone phantoms.
The error between the displacements measured by the registration of the US image slices and the true displacements of the respective slices measured using the precision mechanical stages on the experimental apparatus is evaluated for translation and rotation in two simulated environments. The mean and standard deviation of errors are shown in tabular form. This method provides an average measurement precision of less than 0.1 mm and 0.1 degrees, respectively.
In this paper, we have presented a novel non-invasive approach to measuring the motion of the bones in a knee using tri-plane B-mode ultrasound and image registration. In our study, the image registration method determines the position of bony landmarks relative to a B-mode ultrasound sensor array with sub-pixel accuracy. The advantages of our proposed system over previous techniques are that it is non-invasive, does not require the use of ionizing radiation and can be used conveniently if miniaturized.
KeywordsMotion analysis Ultrasound Image registration Knee joint diagnostics
Kinematic analysis allows the motion of individual bones in the knee joint to be measured [1, 2]. These motion measurements provide significant insights into normal and abnormal joint trajectories which can lead to improved artificial joint component design. They also allow improved diagnosis for ligament injuries and facilitate accurate therapeutic strategies to be formulated. Knee joint kinematics is complex in a sense that the natural motions of flexion, extension and internal, external rotation are coupled [3, 4]. During flexion, the two main bones in a knee joint, do not bend like a door hinge, instead they rotate about several specific axes. This type of motion makes it crucial that the joint kinematics be evaluated in 3D space.
Currently, the clinical standard for 3D modelling of joint kinematics, implant performance and implant bearing wear is Roentgen Stereo Analysis (RSA) [5–8]. This approach involves implanting tantalum beads or markers in the bones and then taking x-rays projected through the joint in two imaging planes to measure the 3D kinematics. The reported precision (standard deviation of error) for locating the beads using this technique is 0.2° to 0.8° for rotation and 0.1 to 0.5 mm for translation [8, 9]. While the precision of this technique is considered to be acceptable, the main disadvantages of the method are a) the method can only be applied for post-operative assessment of artificial knee joints due to the requirement of implanting tantalum beads into the bone during knee replacement surgery, b) the procedure requires the patient to be exposed to ionizing radiation during the capture of the dual x-ray images. This exposure has an associated elevated risk of the patient contracting cancer and c) the x-ray images are captured while the patient is lying in a horizontal position. Hence this procedure does not allow the 3D motion analysis of the knee joint while the patient performs normal everyday activities.
Another technique is based on intra-cortical pins and bone marker [3–5]. This technique is accurate, but still highly invasive. An alternate solution for performing the required 3D motion analysis is achieved by merging 2D x-ray videos with 3D CT scans [10, 11]. It provides almost the same accuracy as RSA . While this technique does not require tantalum beads to be implanted it does still have the following disadvantages: a) ionizing radiation exposure to the patient is required and b) the range of movement of the patient is limited to within the small field of view of the x-ray equipment.
The most non-invasive technique for kinematic analysis of joints is the optical tracking of skin mounted markers . Commonly used products in this category are Vicon, Qualysis, Polaris, Opto-trak (Certus) and Optical Tracking Systems (Northern Digital Inc., Waterloo, ON, Canada). Although optical tracking provides sub-millimeter accuracy for kinematic analysis using skin mounted markers and is very fast for tracking multiple targets, the main disadvantage of this technique is that the markers on the skin can move independently of the underlying bone due to the flexibility of the soft muscle tissue between the bone and the skin-mounted marker. This independent movement naturally leads to soft tissue artifacts (STA).
In summary, conventional motion analysis systems utilize optoelectronic, ultrasonic or video based systems to track markers attached to the skin. These systems are non-invasive, easy-to-operate and widely used in motion analysis, computer graphics, animations and have also been proposed for computer assisted surgery [10, 14–18]. However, studies have shown that they do not provide the precision necessary which is generally 2-4 mm  for many clinical applications due to the large relative movement of skin and soft tissues with respect to the underlying bones during dynamic activities [20, 21]. A non-invasive ultrasound and image registration based approach will be described in the following sections to mitigate the error due to soft-tissue artifact in motion analysis of human knee joints.
Measuring bone position
In our proposed framework, a novel H-shaped ultrasound array geometry is used to capture all six motion parameters for 3D rigid body movement estimation between bone and sensor. The three probes were arranged in such a way as to capture three 2D US image slices in a novel H shaped orientation. For US image acquisition, the three US arrays are mounted firmly to prevent any relative motion among the image slices as well as calibrated to ensure they preserve the H shape orientation. In this study, US slices are acquired with a frame rate of 30 frames per second (fps) at 24 MHz with a resolution of 0.132 mm. Five out of the six motion parameters required to describe the 3D rigid body motion can be extracted from any two orthogonal US slices but the remaining out-of-plane rotation can only be measured by using twoparallel US slices.
The registration algorithm
Image registration is the process of spatially aligning one or more images to another image, generally called the reference image. It is the process of determining point to point correspondences between images of the same entity. Image registration is a versatiletechnique used in computer vision, video surveillance, automatic target recognition, medical imaging, satellite imaging and some military applications . The goal of image registration in our case is to find the geometric transformation between target and reference images in order to find the relative movement between the US sensor and the underlying bone.
In the registration algorithm, a gradient descent optimization approach proposed by Lucas and Kanade  is used with a new similarity metric called the sum of conditional variance (SCV) to efficiently determine the rigid body transform parameters required to align the bone surface between two corresponding images, and this is the first article to use SCV for US-US image registration. Originally SCV was developed by our group and proposed for multi-modal medical image registration  but recently it has been used for mono-modal registration techniques where non-linear illumination change was an issue. In  the authors showed that its performance was superior when compared to the mutual information (MI) and cross-cumulative residual entropy (CCRE) similarity measures. The motivation of using SCV in our case comes from its resistance to view dependent non-linear intensity variations in US images as well as its ability to handle a large convergence radius with computational ease compared to the other popular similarity measures (e.g. sum-of-squared difference (SSD), sum-of-absolute difference (SAD) etc.).
The SCV can be calculated as follows: suppose the current US image of the bone is I(x i ,y i ) to be registered with the initial US image . The coordinates (x i ,y i ) and denote the locations of pixels in I and R respectively for i=1…N, where N is the total number of pixels in each image. To construct the joint histogram, the images I i and R i were quantized from the original 256 possible values to 64 possible values. Hence, the joint histogram has 64×64 bins corresponding to the 64 possible quantized values in I and R.
Any spatial misalignment between I and R will mean that each pixel value R i will correspond to a range of I i values. However, typically this range of values is approximately equally distributed above and below the value for the registered version of I. The reason for this approximately equal distribution is explained by observing that any spatial misalignment will cause pixel values to change in opposite directions for regions of the image that exhibit opposite spatial gradients. Hence, if a spatial shift of the image causes a value of I i to increase at positions with a positive spatial gradient, then this same shift will cause the same value of I i to decrease in areas with negative spatial gradient and vice-versa. For the approximately equal distribution to occur, it is assumed that there is an approximately equal number of positive and negative spatial gradients in the images which is typically true for the application addressed in this paper.
Since this measure is the sum of the squared difference between values of I i and their conditional mean, this equation actually represents a sum of the conditional variances (SCV) over all values of R i . Conveniently, since the equation for S now has the same form as the equation for the sum-of-squared difference, any standard optimization technique (e.g. Gauss-Newton, Levenberg Marquardt) can be used to find the geometric transform required to register and I. In our experiments we minimized the SCV measure using a gradient descent optimization approach proposed by Lucas and Kanade .
In the first validation experiment, an artificial model of a femur and tibia were placed in a water-filled container as shown in Figure 3a. The use of a water tank is common practice for validating ultrasound related biomedical experiments [26–28]. The water in the tank acted as a coupling medium and simulated the behaviour of muscle tissue. The frequency of the ultrasound signals used to capture the images was 24 MHz and the resolution of the images captured using these probes was 0.132 mm.
The probes were then translated by ±4 mm in steps of 0.5 mm in the x, y and z directions away from the initial positions and rotated by ±3° in steps of 1° around the x, y and z axes. At each position, the US probes were used to capture 2D B-mode scans of the surface of the model bone. These images were then registered to US images captured at the initial position of the probes. The amount of translation or rotation required to register the images was taken as an estimate of the movement of the probe relative to the surface of the bone. In-plane translations (along the x, y and z axes) and two in plane rotations (around the x and y axes) were determined from the two orthogonal US scans captured at the positions of AA’ and CC’ as well as the last out-of-plane rotation (about the z axis) was determined using parallel image slices at AA’ and BB’, according to Figure 2. The estimates of the in-plane translations and rotations came directly from the output of the 2D-2D image registration. From the known center of rotation, the angular displacement of the two parallel probes was calculated and averaged to produce the estimate of the rotation around the z axis.
In the second validation experiment, we deliberately hide the internal bone structure by using an opaque balloon filled with a low concentration (10%, w/w at 20°C) collagen-water solution to mimic the properties of muscle tissue. This setup is shown in Figure 5. Figure 5a shows the experimental setup while Figure 5b shows the balloon phantom which covers the bone and is filled with the collagen-water solution. A typical scan from each of the three sensor arrays is also shown in Figure 5c. This kind of experimental setup can be justified by the fact that, in a real scenario we are not able to see which part of the femur and tibia are currently being scanned. As a result, there is no experimental bias introduced by placing the US probes at the best possible location for the subsequent registration process.
The low concentration collagen-water solution has been shown to closely resemble the acoustic properties of real muscle tissue [29, 30]. The speed of sound in this medium is almost 1540 m/s and collagen, the fibrous protein, is the major constituent of mammalian flesh and connective tissues [31–33]. To ensure successful transmission and reception of ultrasonic pulses through the balloon, ultrasound transmission gel was placed on the transducer and balloon interface. Then the H shaped US transducer assembly was moved carefully on the balloon surface to capture the translated and rotated images of the bone surface in a similar manner to the first water tank experiment. These images were then processed to extract the 6 3D rigid body transform parameters. The experiment as described in the following section, which was carried out on a human subject using accelerometer was approved by the human research ethics advisory panel (HREA) of UNSW Canberra and the approval number is A-14-31.
Results and discussion
Mean, standard deviation (SD) of error & their averages
Translation & rotation error due to different noise levels
In this study, the images were acquired statically at strategic anatomical positions during the experiments. However, in a dynamic activity, there may be significant amount of motion-blur induced in the acquired images, which could reduce the accuracy of this technique. Another considerable fact is the range of translation and rotation analysis performed in this study and this should be consistent with the error analysis due to the soft tissue artefact (STA), which is mainly induced by inertial effects, skin deformation, sliding and deformation caused by muscle contractions in of the knee joints. In , authors reported that maximum translational and rotational error occurred during dynamic activities are 11 mm and 10°, respectively, however, in this study we considered the dynamic ranges of ±4 mm and ±3°, respectively, which may require sufficiently high frame rates to capture motion, which might occur between two subsequent frames. Another issue of the proposed system may be the possible interference between two sensors attached on the knee joint to capture the relative motion between femur and tibia. However, careful selection of the anatomical positions above the joint and controlling the US beam width and directing the beam into the region of interests can virtually reduce the risk of possible interference between these sensors. In  the authors demonstrated a prototype for a low profile US probe with three transducer arrays each containing 96 piezoceramic elements. This probe is called low profile because of its small footprint and light-weight compared to the conventional US probes. Their prototype system contained integrated multiplexing electronics to reduce the size and weight of the probe and provided a viewing angle of almost 41°. This probe was specially designed for trans-rectal imaging for diagnostic assessment of prostate cancer. It is feasible however, that a similar design could be used to develop a low profile H shaped multi-plane US probe with embedded electronic switching circuitry which could be attached to the skin above the ends of the femur or tibia as shown in Figure 1(a). In this study, results show that the shape of the bone’s surface at the positions which were scanned by the US probe has enough detail to allow an accurate measurement of the relative position of the US probe from one scan to another, however, in some anatomical positions, bone details can be insufficient to be used in motion analysis. Also, some measurement uncertainty is introduced due to view dependent intensity variation from one US scan to another, which could further reduce the accuracy of this system. Moreover, in this method, some curvature of the bone contours is a prerequisite for successful tracking. As a result, this proposed technique can not only be used for the motion analysis of knee joints, but can be used for other human joints such as elbow, hand and hip joints provided that anatomical details in the US image slices can feasibly be captured. However, studies have shown that the necessary precision is generally 2-4 mm  for many clinical applications due to the large relative movement of skin and soft tissues with respect to the underlying bones during dynamic activities [20, 21]. And the precision of the proposed system compares favourably with the current clinical standard of RSA which has a reported precision of 0.2° to 0.8° for rotation and 0.1 to 0.5 mm for translation [8, 9].
In this paper we have presented a novel non-invasive approach to measuring motion of the bones in a knee using tri-plane B-mode ultrasound and image registration. However, in practice, we may require to use an optical tracking system to transfer the local bone kinematic information into a global coordinate system and also build a low profileH-shaped US sensor which can be mounted on the patient’s skin using surgical tape. The experimental results reported in this paper prove that sub-millimeter precision is achievable. The use of this type of sensor would allow a system to be developed which could measure the kinematics of knee joints without restricting the patient’s movement to be confined to the field of view of imaging equipment. Such a system would provide accurate measurements while being non-invasive and would not require the use of ionising radiation.
This work has been supported by School of Engineering & IT, UNSW Canberra, under Research Publication Fellowship.
- Kitagawa A, Tsumura N, Chin T, Gamada K, Banks SA, Kurosaka M: In vivo comparison of knee kinematics before and after high-flexion posterior cruciate-retaining total knee arthroplasty. J Arthroplasty 2010,25(6):964–969. 10.1016/j.arth.2009.07.008View ArticleGoogle Scholar
- Smith PN, Refshauge KM, Scarvell JM: Development of the concepts of knee kinematics. Arch Phys Med Rehabil 2003,84(12):1895–1902. 10.1016/S0003-9993(03)00281-8View ArticleGoogle Scholar
- Giphart JE, Zirker CA, Myers CA, Pennington WW, LaPrade RF: Accuracy of a contour-based biplane fluoroscopy technique for tracking knee joint kinematics of different speeds. J Biomec 2012,45(16):2935–2938. 10.1016/j.jbiomech.2012.08.045View ArticleGoogle Scholar
- Lafortune MA, Cavanagh PR, Sommer HJ III, Kalenak A: Three-dimensional kinematics of the human knee during walking. J Biomech 1992,25(4):347–357. 10.1016/0021-9290(92)90254-XView ArticleGoogle Scholar
- Leardini A, Chiari L, Croce UD, Cappozzo A: Human movement analysis using stereophotogrammetry: part 3. soft tissue artifact assessment and compensation. Gait Posture 2005,21(2):212–225. 10.1016/j.gaitpost.2004.05.002View ArticleGoogle Scholar
- You B-M, Siy P, Anderst W, Tashman S: In vivo measurement of 3-d skeletal kinematics from sequences of biplane radiographs: application to knee kinematics. Med Imaging IEEE Trans 2001,20(6):514–525. 10.1109/42.929617View ArticleGoogle Scholar
- Selvik G: Roentgen stereophotogrammetric analysis. Acta Radiol 1990,31(2):113–126. 10.3109/02841859009177472View ArticleGoogle Scholar
- Kärrholm J, Selvik G, Elmqvist L-G, Hansson LI: Active knee motion after cruciate ligament rupture. Acta Orthopaedica 1988,59(2):158–164.Google Scholar
- Trozzi C, Kaptein BL, Garling EH, Shelyakova T, Russo A, Bragonzoni L, Martelli S: Precision assessment of model-based rsa for a total knee prosthesis in a biplanar set-up. Knee 2008,15(5):396–402. 10.1016/j.knee.2008.05.001View ArticleGoogle Scholar
- Fregly BJ, Rahman HA, Banks SA: Theoretical accuracy of model-based shape matching for measuring natural knee kinematics with single-plane fluoroscopy. J Biomech Eng 2005,127(4):692–699. 10.1115/1.1933949View ArticleGoogle Scholar
- Scarvell JM, Pickering MR, Smith PN: New registration algorithm for determining 3d knee kinematics using ct and single-plane fluoroscopy with improved out-of-plane translation accuracy. J Orthop Res 2010,28(3):334–340.Google Scholar
- Zwingmann J, Konrad G, Kotter E, Südkamp N, Oberst M: Computer-navigated iliosacral screw insertion reduces malposition rate and radiation exposure. Clin Orthop Relat Res 2009,467(7):1833–1838. 10.1007/s11999-008-0632-6View ArticleGoogle Scholar
- Tranberg R, Saari T, Zügner R, Kärrholm J: Simultaneous measurements of knee motion using an optical tracking system and radiostereometric analysis (rsa). Acta Orthop 2011,82(2):171–176. 10.3109/17453674.2011.570675View ArticleGoogle Scholar
- Miranda DL, Rainbow MJ, Crisco JJ, Fleming BC: Kinematic differences between optical motion capture and biplanar videoradiography during a jumpcut maneuver. J Biomech 2013,46(3):567–573. 10.1016/j.jbiomech.2012.09.023View ArticleGoogle Scholar
- Banks SA, Hodge WA: Accurate measurement of three-dimensional knee replacement kinematics using single-plane fluoroscopy. Biomed Eng IEEE Trans 1996,43(6):638–649. 10.1109/10.495283View ArticleGoogle Scholar
- Zuffi S, Leardini A, Catani F, Fantozzi S, Cappello A: A model-based method for the reconstruction of total knee replacement kinematics. Med Imaging IEEE Trans 1999,18(10):981–991. 10.1109/42.811310View ArticleGoogle Scholar
- Tashman S, Anderst W: In vivo measurement of dynamic joint motion using high speed biplane radiography and ct: application to canine acl deficiency. J Biomech Eng 2003, 125: 238–245. 10.1115/1.1559896View ArticleGoogle Scholar
- Mozes A, Chang T-C, Arata L, Zhao W: Three-dimensional a-mode ultrasound calibration and registration for robotic orthopaedic knee surgery. Int J Med Robot Comput Assist Surg 2010,6(1):91–101.Google Scholar
- Phillips R: (ii) the accuracy of surgical navigation for orthopaedic surgery. Curr Orthop 2007,21(3):180–192. 10.1016/j.cuor.2007.06.006View ArticleGoogle Scholar
- Benoit DL, Ramsey DK, Lamontagne M, Xu L, Wretenberg P, Renström P: Effect of skin movement artifact on knee kinematics during gait and cutting motions measured in vivo. Gait Amp Posture 2006,24(2):152–164. 10.1016/j.gaitpost.2005.04.012View ArticleGoogle Scholar
- Reinschmidt C, van den Bogert AJ, Nigg BM, Lundberg A, Murphy N: Effect of skin movement on the analysis of skeletal knee joint motion during running. J Biomech 1997,30(7):729–732. 10.1016/S0021-9290(97)00001-8View ArticleGoogle Scholar
- Goshtasby A: 2-D and 3-D Image Registration for Medical, Remote Sensing, and Industrial Applications. John Wiley & Sons; 2005.Google Scholar
- Baker S, Matthews I: Lucas-kanade 20 years on: a unifying framework. Int J Comput Vis 2004,56(3):221–255.View ArticleGoogle Scholar
- Pickering MR, Muhit AA, Scarvell JM, Smith PN: A new multi-modal similarity measure for fast gradient-based 2d-3d image registration. Engineering in Medicine and Biology Society, 2009. EMBC 2009. Annual International Conference of the IEEE 2009, 5821–5824.View ArticleGoogle Scholar
- Richa R, Sznitman R, Taylor R, Hager G: Visual tracking using the sum of conditional variance. Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ International Conference On 2011, 2953–2958.View ArticleGoogle Scholar
- Chen W, Lou EHM, Le LH: Using ultrasound imaging to identify landmarks in vertebra models to assess spinal deformity. Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE 2011, 8495–8498.View ArticleGoogle Scholar
- Mitri FG, Greenleaf JF, Fatemi M: Chirp imaging vibro-acoustography for removing the ultrasound standing wave artifact. Med Imaging IEEE Trans 2005,24(10):1249–1255.View ArticleGoogle Scholar
- Ashfaq M, Brunke SS, Dahl JJ, Ermert H, Hansen C, Insana MF: An ultrasound research interface for a clinical system. Ultrason Ferroelectr Freq Control IEEE Trans 2006,53(10):1759–1771.View ArticleGoogle Scholar
- Parker NG, Povey MJW: Ultrasonic study of the gelation of gelatin: phase diagram, hysteresis and kinetics. Food Hydrocolloids 2012,26(1):99–107. 10.1016/j.foodhyd.2011.04.016View ArticleGoogle Scholar
- Frizzell LA, Gindorf JD: Measurement of ultrasonic velocity in several biological tissues. Ultrasound Med Amp Biol 1981,7(4):385–387. 10.1016/0301-5629(81)90049-1View ArticleGoogle Scholar
- Pillen S, Tak RO, Zwarts MJ, Lammens MMY, Verrijp KN, Arts IMP, van der Laak JA, Hoogerbrugge PM, van Engelen BGM, Verrips A: Skeletal muscle ultrasound: correlation between fibrous tissue and echo intensity. Ultrasound Med Biol 2009,35(3):443–446. 10.1016/j.ultrasmedbio.2008.09.016View ArticleGoogle Scholar
- Zaidman CM, Holland MR, Hughes MS: Quantitative ultrasound of skeletal muscle: Reliable measurements of calibrated muscle backscatter from different ultrasound systems. Ultrasound Med Biol 2012,38(9):1618–1625. 10.1016/j.ultrasmedbio.2012.04.020View ArticleGoogle Scholar
- Saxon GI, Tracy JM: The Molecular Structure of Collagen. 1965.Google Scholar
- UMI: Philips IE33 Ultrasound Machine An Imaging Revelation. . http://@www.umiultrasound.com/[http://www.umiultrasound.com/ultrasound-machine/philips/ie33]
- Uddin MS, Tahtali M, Pickering MR: Complex wavelet based speckle reduction using multiple ultrasound images. Proc. SPIE, Sixth International Conference on Digital Image Processing (ICDIP 2014) 2014, 91591I–91591I-6.Google Scholar
- Lee M-S, Yen C-L, Ueng S-K: Speckle reduction with edges preservation for ultrasound images: using function spaces approach. Image Process IET 2012,6(7):813–821. 10.1049/iet-ipr.2011.0366MathSciNetView ArticleGoogle Scholar
- Liu G, Zeng X, Tian F, Li Z, Chaibou K: Speckle reduction by adaptive window anisotropic diffusion. Signal Process 2009,89(11):2233–2243. 10.1016/j.sigpro.2009.04.042View ArticleGoogle Scholar
- Uddin MS, Tahtali M, Lambert A, Pickering MR: Speckle reduction for ultrasound images using nonlinear multi-scale complex wavelet diffusion. Signal and Image Processing Applications (ICSIPA), 2013 IEEE International Conference On 2013, 31–36.View ArticleGoogle Scholar
- Garling EH, Kaptein BL, Mertens B, Barendregt W, Veeger HEJ, Nelissen RGHH, Valstar ER: Soft-tissue artefact assessment during step-up using fluoroscopy and skin-mounted markers. J Biomech 2007,40(Supplement 1):18–24.View ArticleGoogle Scholar
- Felix N, Voisin D, Clade O, Dufait R: Biplane ultrasound arrays with integrated multiplexing solution for enhanced diagnostic accuracy in endorectal and transvaginal imaging. Ultrasonics Symposium, 2005 IEEE 2005, 2251–2254.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.