A 3D Sparse Motion Field Filtering for Quantitative Analysis of Fascial Layers Mobility Based on 3D Ultrasound Scans G. Turini, S. Condino, A. Stecco, V. Ferrari, M. Ferrari, and M. Gesi
Abstract— In the last few years, there has been an increasing interest in the role of deep fascia mobility in musculoskeletal dynamics and chronic pain mechanisms. In a previous paper we presented an innovative semiautomatic approach to evaluate the 3D motion of the fascia using ultrasound (US) imaging, generating a sparse deformation vector field. This paper presents an improvement of our original method, focusing on the filtering of the sparse vector field and its validation. Moreover, in order to evaluate the performance of the algorithm, a method is proposed to generate synthetic deformation vector fields, including: expansion, rotation, horizontal shear, and oblique shear components. Preliminary tests on the final synthetic deformation vector fields showed promising results. Further experiments are required in order to optimize the tuning of the algorithm.
contraction. The proposed approach, which attempts to contribute to the knowledge of the myofascial system and of deep fascia role in musculoskeletal dynamics and dysfunctions, relies on the acquisition of two 3D US datasets: at rest (Rest Volume, Volrest) and during a voluntary muscular contraction (Contraction Volume, Volcont) [6]. More particularly, our method consists of two phases: the 3D US dataset analysis and generation of a displacement vector field using a block matching technique (Phase 1), and the validation of the resulting displacement vector field filtered for outliers removal (Phase 2) [6]. The objectives of this research study are the following:
to propose a method for the validation of Phase 2 using a configurable synthetic vector field representing a generic deformation (resulting from the combination of different basic deformations, and integrating custom random noise);
to enhance Phase 2 adding a new filter to select the best feature match couple;
to show the preliminary results of the validation of the algorithm, and an example of a comparison between the original and the current versions.
I. INTRODUCTION The deep fascia is a continuous fibrous membrane that protects, enfolds and divides several anatomical parts: muscles, nerves, vessels, ligaments, and joints just to name a few. This protective membrane is organized into different layers depending on the anatomical region, and it is adaptable to the dynamics of the underlying muscles thanks to its sliding layers [1] [2]. In previous years, several studies have demonstrated that the deep fascia is of utmost importance for movement perception and coordination, muscular interaction and force transmission, tensional stress and etiology of pain [1]. Recently, evidence has been found to establish the significance of fascial layers mobility for optimal muscular activity and joint movement [1], and it has been hypothesized that any issue in the smooth sliding of fascial layers could result in improper coordinated movements [3]. In the past few years, researchers have investigated the deep fascia exploiting ultrasound (US) imaging [4], making possible to observe pathological changes of the fascial layers, and to evaluate the deep fascia physiology [5]. In a previous study [6] we have defined a novel approach to study in-vivo fascial mobility with the semiautomatic generation of 3D motion vector fields describing the displacement of salient fascial features during a muscular
G. Turini is with the Computer Science Department, Kettering University, Flint, MI 48504 USA (e-mail: [email protected]). S. Condino is with the EndoCAS Center, Department of Translational Research on New Technologies in Medicine and Surgery, University of Pisa, Ospedale di Cisanello, Via Paradisa 2, 56124 Pisa, Italy (corresponding author; phone: +39-050-99-5689; fax: +39-050-99-5689; e-mail: [email protected]). A. Stecco is with the Department of Internal Medicine, Sport Medicine Unit, University of Padua, Padua, Italy (e-mail: [email protected]).
978-1-4244-9270-1/15/$31.00 ©2015 IEEE
II. ALGORITHM DESCRIPTION This section describes our approach to analyze the deep fascia using volumetric US scans. In the proposed method, the motion of each fascial layer is studied individually by analyzing the displacement of points which correspond to clearly distinguishable anatomical features (e.g. fascia/intermuscular septum intersections). A. Phase 1: Generation of 3D Motion Field from US Scans This phase is performed through the following sequence of steps: Image Filtering to improve the visualization of fascial layers [7, 8], Regions of Interest Selection to select the target fascial layer [9, 10], Feature Detection and Extraction to identify clearly distinguishable anatomical features [11], and Feature Matching [12, 13].
G. Turini, V. Ferrari, M. Ferrari, and M. Gesi are with the EndoCAS Center, Department of Translational Research on New Technologies in Medicine and Surgery, University of Pisa, Ospedale di Cisanello, Via Paradisa 2, 56124 Pisa, Italy (e-mail: [email protected], [email protected], [email protected]). V. Ferrari is with the Information Engineering Department, University of Pisa, Pisa, Italy.
775
The most important step in this phase is the Feature Matching, which aims to find corresponding pairs of feature points in two different datasets: the first feature is located in a Reference Volume, whereas the second feature is located in a Target Volume (Figure 1). In our approach, this procedure is performed twice: the first time to find for each feature point in Volrest (Reference Volume), its corresponding feature point in Volcont (Target Volume); the second time vice versa, to find for each feature point in Volcont (Reference Volume), its corresponding feature point in Volrest (Target Volume) (Figure 1). This procedure creates two sets of feature matches: RC Matches (feature matches with Volrest as Reference Volume and Volcont as Target Volume), and CR Matches (feature matches with Volcont as Reference Volume and Volrest as Target Volume). Each RC and CR feature match is then represented as a 3D displacement vector describing the motion of the relative anatomical feature between the contraction and rest phases. The double execution of the Feature Matching has two advantages: an increased number of potential feature matches (default RC Matches are extended adding CR Matches) by using salient features in both volumes to generate a more dense motion vector field; in addition, it improves the validation of feature matches exploiting the fact that RC Matches and CR Matches have been generated independently (i.e. a strong similarity between a RC match and a CR match is a significant evidence of the reliability of both feature matches).
B. Phase 2: Validation and Filtering of 3D Motion Field The sparse motion vector field resulting from Phase 1 is then filtered by Phase 2, which consists in an automatic validation process designed to select only reliable feature matches, removing uncertain matches. This phase is split into two steps: Match Certification, and Match Recovery; and the final result consists in a classification of all the feature matches into three distinct categories: Certified Matches (highly reliable), Compatible Matches (reliable), and Conflicting Matches (unreliable) (Figure 2). The first stage of Phase 2 is the Match Certification (Figure 2a), and performs the identification of Certified Matches. For each RC Match RCm = {Rm, Cm} we search the closest CR Match CRn = {Cn, Rn}, where: Rm and Rn are 3D feature points in Volrest and Cm and Cn are 3D feature points in Volcont. The search for CRn is based on the evaluation of two distance vectors: Rnm = Rn – Rm, and Cnm = Cn – Cm. Then, the pair RCm/CRn undergoes a Distance Check (1) and a Motion Compatibility Check (2) (Figure 1 and Figure 2). ( ‖Rnm ‖ < Th1 ) ∧ ( ‖Cnm ‖ < Th1 )
(
‖ Dm - Dn ‖ ‖ Dm ‖
< Th2 ) ∧ (
‖ Dm - Dn ‖ ‖ Dn ‖
< Th2 ) ∧ (Amn < Th3 )
In (1) and (2): Dm = Rm – Cm, Dn = Rn – Cn, Amn is the angle between Dm and Dn, and Th1, Th2 and Th3 are threshold
Figure 1: An example of a generic RCm/CRn pair including the representation of all the vectors used in Phase 2. (A) and (B) show two cross sections of the Reference Volume (Volrest in this case) and Target Volume (Volcont in this case), illustrating the two feature points Rm and Cm of a RC feature match RCm. (D), (E), and (F) illustrate the reference systems used in the Rest Volume Vol rest and in the Contraction Volume Volcont, and how these are unified to compute the candidate displacement vectors Dm and Dn.
776
Figure 2: Diagram of Phase 2 illustrating the algorithm to validate feature matches: (A) step 2.1 Match Certification is used to select highly reliable matches (Certified Matches), whereas (B) step 2.2 aims to exploit the Certified Matches in order to recover previously discarded feature matches, improving in this way the density of the final motion vector field.
values. Equation (1) ensures that RCm and CRn are located in the same anatomical region, whereas (2) ensures that RCm and CRn are similar in terms of their normalized modulus and direction. If these two feature matches don’t satisfy one of these checks (Distance Check and Motion Compatibility Check) they are temporarily skipped, and will be analyzed again in the second stage of Phase 2. Otherwise, if both checks are passed, the couple RCm/CRn is inserted in a temporary set containing all the couples RCm/CRk...r that have passed both checks. Then, we use a cost function fcost (3-9) to evaluate each couple in this set in order to select the best CR match CRn. Finally we classify both RCm and CRn as Certified Matches, and repeat the process for the next RC match.
fcost ( RCm / CRn ) = ∑6i=1 vi wi ,
v1 =
v2 =
| ‖ Rmn ‖ - Th1 |
| ‖ Cmn ‖ - Th1 |
Th1 Th1
v3 = max ( 0, 1 - |
v4 =
v5 =
|
‖ Rmn ‖ ‖ Cmn ‖
‖ Dm - Dn ‖ - Th2 | ‖ Dm ‖
Th2 |
∑6i=1 wi =1
‖ Dm - Dn ‖ - Th2 | ‖ Dn ‖
Th2
| )
v6 =
| Amn - Th3 | Th3
The cost function fcost in (3) evaluates the quality of a feature match couple by computing a weighted arithmetic mean of six different properties: the distance (Rmn) between the feature points in Volrest (Rm, and Rn) normalized on a threshold value (4), the distance (Cmn) between the feature points in Volcont (Cm, and Cn) normalized on a threshold value (5), the ratio between the previous two distances (Rmn, and Cmn) (6), normalized ratios of the displacement vectors Dm and Dn (7) (8), and angle between RCm and CRn normalized on a threshold value (9). This improved version of the algorithm includes new stages (in bright red in Figure 2) in order to select the best CR Match to pair with the input RC Match RCm. In the original version of the algorithm we considered only the closest CR Match to RCm. Now, instead, we collect a group of CR Matches close to RCm, and then we choose the best CR Match CRn evaluating the quality of all the matches in the group using the cost function fcost (3). Match Recovery is the second stage of Phase 2, and starts selecting all the feature matches which have not been certified: from now on called Non-Certified Match (NCM). Then, for each NCMi (with its relative Ri and Ci), the closest Certified Match CMj (with its relative Rj and Cj) is searched evaluating the distance vectors: Rij = Ri – Rj, and Cij = Ci – Cj. NCMi is
777
then classified as Compatible Match if it satisfies the Motion Compatibility Check as rephrased in (3).
(
‖ Di - Dj ‖ ‖ Di ‖
< Th4 ) ∧ (
‖ Di - Dj ‖ ‖ Dj ‖
< Th4 ) ∧ (Aij < Th5 )
In (11): Di = Ri – Ci, Dj = Rj – Cj, Aij is the angle between Di and Dj, and Th4 and Th5 are threshold values. Otherwise, if (10) is not satisfied, NCMi is classified as Conflicting Match and discarded. At the end, reliable matches (Certified Matches and Compatible Matches) are used for fascial motion estimation. III. TESTING METHODOLOGIES The main goal or our tests is the evaluation of the performance of our algorithm on a configurable deformation vector field, enabling a per-vector analysis of the final results. To do this, we generated a custom synthetic sparse deformation vector field which simulates the deformation of a fascial layer, and then we processed this vector field using our algorithm in order to obtain a validated set of deformation vectors discarding outliers. A. Deformation Modeling and Sparse Vector Field Generation To generate a generic 3D deformation vector field, we computed a deformation vector for each voxel of an input volumetric dataset. So, we multiplied the 3D position vi of each voxel i in the volumetric dataset by a deformation matrix L, obtaining the 3D deformation vector di applied on the voxel i (11).
vi L = di
In (11), di is the displacement vector of the feature located in voxel i, and it simulates a feature match as if generated by Phase 1. In this case, we approximated a generic three– dimensional deformation by a linear 3D vector field described by the matrix L [14, 15]: L11 L12 L13 1 0 0 L = [L21 L22 L23 ] = E [ 0 1 0 ] + L31 L32 L33 0 0 1 1 0 0 1 0 0 1 0 + R [ 0 0 -1 ] + S1 [ 0 1 0 ] + S2 [ 0 0 0 1 0 0 0 -1 0 1 where
E= S1 =
L11 + L22 + L33 3 L11 + L22 - L33 3
R= S2 =
0 1] 0
fascial layer. Now, to simulate the output of Phase 1, two different random selections of displacement vectors are performed to create the two sets of feature matches: RC Matches and CR Matches, the input required by Phase 2 (Figure 2). Then, the resulting sparse vector field (including both RC and CR matches) is affected by random noise in order to test Phase 2 robustness. B. Experimental Setup The generation of both the dense and sparse vector fields, and the entire algorithm (both Phase 1 and Phase 2) were performed in Matlab® version 7.14, release 2012a (MathWorks, Natick, Massachusetts, USA). The size of the input volumetric datasets used to test the algorithm are 5×370×100 voxels (Sx×Sy×Sz) and the density of the simulated sparse deformation field is 0.05% (for a total of 92 vectors per field). These values were chosen to obtain a deformation field similar to those generated analyzing real fascial volumes as described in our previous work [6]. The thresholds Th1, Th2, Th3, Th4, and Th5 of Phase 2 were set to 2.89 voxels, 0.4, 20°, 1, and 40°, respectively. These values are those preliminarily set in [6] to optimize the algorithm performance on the acquired datasets. The weights used in the cost function fcost (3) were set as follows: w1=w2=w4=w5=0.1; w3= w6=0.3. Phase 2 was tested processing 15 different deformations. The scalar deformation coefficients E, R, S1, and S2 were set to 0 or 0.1. All the possible combinations of the elementary components (expansion, rotation, horizontal shear, and oblique shear) were tested as illustrated in Figure 3. During the sparse vector field generation, we applied a random error to each displacement vector adding a 3D vector with x, y, and z coordinates ranging respectively from 0 to 2% of Sx, Sy, and Sz. For each of the 15 deformations tested, we processed 10
L11 - L23 + L32 3 L11 + L23 + L32
3
In (12), L is decomposed as a sum of four matrices each representing an elementary deformation component (13): expansion, rotation, horizontal shear, and oblique shear (Figure 3). Each deformation component (matrix) is controlled by a single scalar coefficient: E, R, S1, and S2 respectively [14, 15]. Equation (12) describes how our synthetic deformation field is generated, that is: the 3D position of (the centroid of) each voxel v is multiplied by matrix L, and the result is the 3D deformation vector d affecting that voxel, as described in (11). The vector field obtained so far is a single dense deformation field that represents the deformation affecting a 778
Figure 3: An example of the four different elementary deformations used to compose the generic deformation matrix L: expansion (A), rotation (B), horizontal shear (C), and oblique shear (D).
different sparse 3D vector fields, and the results were summarized in Table 1 as: the ratio between the modulus of the error vector and the modulus of the relative original displacement vector (without noise) expressed as a percentage (Mod), and the angle between the error vector and the relative original displacement vector (Ang). Then, we performed the same tests increasing the maximum error percentage of 1% each time: from 2% to 8%. The diagram in Figure 5 illustrates how the error affects the classification in terms of the percentages of Certified, Compatible, and Conflicting Matches. IV. PRELIMINARY RESULTS As expected, values reported in Table 1 show that the error (Mod and Ang) affecting Certified and Compatible Matches is always lower that the error of Conflicting Matches. The only exception is relative to a pure rotational deformation (row 4, gray cells), for which we have an error in the angle slightly greater for Compatible Matches. Figure 4 shows an example of a section of a 3D deformation vector field as resulting after Phase 2, and illustrates: Certified, Compatible (both RC and CR) and Conflicting Matches with the relative error vector. In the same figure, black circles show a pair of feature matches classified by the current improved version of Phase 2 and the original version of our algorithm: highlighting that Phase 2 is now able to classify as Certified Matches two displacement vectors previously classified as Compatible and Conflicting Matches. ERROR OF CLASSIFIED MATCHES AS DEFORMATION CHANGES
Deformation Type
Error of Classified Matches Certified
Compatible
Conflicting
E
R
S1
S2
Mod (%)
Ang (°)
Mod (%)
Ang (°)
Mod (%)
Ang (°)
0
0
0
0.1
16.2
41.5
26.8
64.4
45.4
65.2
0
0
0.1
0
25.9
42.6
32.8
53.4
37.1
70.8
0
0
0.1
0.1
14.8
45.4
18.8
56.6
26.3
62.2
0
0.1
0
0
18.6
44.5
25.1
58.7
41.6
57.1
0
0.1
0
0.1
11.7
50.9
14.2
57.8
125.9
58.5
0
0.1
0.1
0
17.9
44.6
34.4
55.8
79.9
59.3
0
0.1
0.1
0.1
10.0
50.8
11.9
53.3
40.4
60.4
0.1
0
0
0
23.0
39.4
31.4
50.3
38.2
68.0
0.1
0
0
0.1
12.2
51.5
19.7
55.1
87.2
66.4
0.1
0
0.1
0
20.9
46.8
25.3
47.8
98.2
80.4
0.1
0
0.1
0.1
10.6
55.3
18.1
54.1
46.9
64.3
0.1
0.1
0
0
13.8
43.5
19.4
56.3
30.0
59.9
0.1
0.1
0
0.1
15.4
51.5
13.2
55.1
50.3
64.4
0.1
0.1
0.1
0
11.3
52.6
15.7
55.4
40.7
59.3
0.1
0.1
0.1
0.1
10.7
49.6
14.1
58.9
86.8
59.9
Analyzing the classification of feature matches as the maximum error increases, it is clear (as illustrated in Figure 5) that the percentage of Conflicting Matches increases accordingly. This happens because while increasing the error, we are at the same time decreasing the percentage of displacement vectors congruent with the original deformation vector field. V. DISCUSSION In conclusion, this paper presents an improvement of the method proposed in [6] to study in-vivo the mobility of the deep fascia using a semiautomatic approach. In particular, this study is focused on the final filtering of the 3D sparse vector field representing the deformation of a fascial layer, proposing an enhanced algorithm and a method for its validation. Our algorithm has been tested on different synthetic deformation vector fields and as, preliminarily demonstrated by the results obtained, it works for deformations including expansion, rotation, horizontal shear, and oblique shear components. All our tests were performed using a fixed configuration of the algorithm settings. MATCH CLASSIFICATION AS NOISE CHANGES
Mean Match Type (%)
TABLE 1
Figure 4: A section of a 3D deformation vector field as resulting after Phase 2. The black circles show the comparison between the current improved version of Phase 2 (top) and the original version of the same stage (bottom).
100% 80%
Conflicting Matches
60% 40%
Compatible Matches
20% 0% 2
3
4
5
6
7
8
Certified Matches
Max Noise Amplitude (%) Figure 5: Differences in the classification of feature matches as the maximum noise applied increases from 2% to 8%. Each value illustrated represents the average percentage over 15 different deformations each tested using 10 different sparse vector fields.
779
In the future, we plan to extensively test the proposed method, varying the configuration of the deformation and the applied error in order to optimally tune the parameters of our algorithm. ACKNOWLEDGMENT The authors would like to thank Eng. Simone Parrini for his essential contribution throughout the initial research project that has led to this study. REFERENCES [1] C. Stecco, V. Macchi, A. Porzionato, F. Duparc and R. De Caro, "The fascia: the forgotten structure," Italian Journal of Anatomy and Embryology, vol. 116, no. 3, pp. 127-138, 2011. [2] G. Natale, S. Condino, P. Soldani, F. Fornai, M. Mattioli Belmonte and M. Gesi, "Natale et. al.'s response to Stecco's fascial nomenclature editorial," Journal of Bodywork and Movement Therapies, vol. 18, no. 4, pp. 588-590, 2014. [3] C. Stecco and W. I. Hammer, Functional atlas of the human fascial system, Elsevier Health Sciences, 2014. [4] H. M. Langevin, D. Stevens-Tuttle, J. R. Fox, G. J. Badger, N. A. Bouffard, M. H. Krag, J. Wu and S. M. Henry, "Ultrasound evidence of altered lumbar connective tissue structure in human subjects with chronic low back pain," BMC Musculoskeletal Disorders, vol. 10, no. 1, p. 151, 2009. [5] A. Stecco, A. Meneghini, R. Stern, C. Stecco and M. Imamura, "Ultrasonography in myofascial neck pain: randomized clinical trial for diagnosis and follow-up," Surgical and Radiologic Anatomy, vol. 36, no. 3, pp. 243-253, 2014. [6] S. Condino, G. Turini, S. Parrini, A. Stecco, F. Busoni, V. Ferrari, M. Ferrari and M. Gesi, "A semiautomatic method for in vivo three-dimensional quantitative analysis of fascial layers mobility based on 3D ultrasound scans," International Journal of Computer Assisted Radiology and Surgery (IJCARS), pp. 1-15, 2015. [7] A. F. Frangi, W. J. Niessen, K. L. Vincken and M. A. Viergever, "Multiscale vessel enhancement filtering," in Medical Image Computing and Computer-Assisted Intervention (MICCAI), 1998. [8] M. Rana and J. M. Wakeling, "In-vivo determination of 3D muscle architecture of human muscle using free hand ultrasound," Journal of Biomechanics, vol. 44, no. 11, pp. 2129-2135, 2011. [9] P. A. Yushkevich, J. Piven, H. C. Hazlett, R. G. Smith, S. Ho, J. C. Gee and G. Gerig, "User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability," Neuroimage, vol. 31, no. 3, pp. 1116-1128, 2006. [10] V. Ferrari, C. Cappelli, G. Megali and A. Pietrabissa, "An anatomy driven approach for generation of 3D 780
models from multi-phase CT images," International Journal of Computer Assisted Radiology and Surgery (IJCARS), vol. 3, no. 1, pp. 272-273, 2008. [11] J. Shi and C. Tomasi, "Good features to track," in Proceedings of IEEE Computer Vision and Pattern Recognition (CVPR), 1994. [12] O. Commowick, N. Wiest-Daesslé and S. Prima, "Block-matching strategies for rigid registration of multimodal medical images," in Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI), 2012. [13] A. Gastounioti, S. Golemati, J. S. Stoitsis and K. S. Nikita, "Carotid artery wall motion analysis from Bmode ultrasound using adaptive block matching: in silico evaluation and in vivo application," Physics in Medicine and Biology, vol. 58, no. 24, p. 8647, 2013. [14] D. Zoccolan and V. Torre, "Using optical flow to characterize sensory-motor interactions in a segment of the medicinal leech," The Journal of neuroscience, vol. 22, no. 6, pp. 2283-2298, 2002. [15] A. Giachetti and V. Torre, "The use of optical flow for the analysis of non-rigid motions," International Journal of Computer Vision, vol. 18, no. 3, pp. 255279, 1996.
Abstract— In the last few years, there has been an increasing interest in the role of deep fascia mobility in musculoskeletal dynamics and chronic pain mechanisms. In a previous paper we presented an innovative semiautomatic approach to evaluate the 3D motion of the fascia using ultrasound (US) imaging, generating a sparse deformation vector field. This paper presents an improvement of our original method, focusing on the filtering of the sparse vector field and its validation. Moreover, in order to evaluate the performance of the algorithm, a method is proposed to generate synthetic deformation vector fields, including: expansion, rotation, horizontal shear, and oblique shear components. Preliminary tests on the final synthetic deformation vector fields showed promising results. Further experiments are required in order to optimize the tuning of the algorithm.
contraction. The proposed approach, which attempts to contribute to the knowledge of the myofascial system and of deep fascia role in musculoskeletal dynamics and dysfunctions, relies on the acquisition of two 3D US datasets: at rest (Rest Volume, Volrest) and during a voluntary muscular contraction (Contraction Volume, Volcont) [6]. More particularly, our method consists of two phases: the 3D US dataset analysis and generation of a displacement vector field using a block matching technique (Phase 1), and the validation of the resulting displacement vector field filtered for outliers removal (Phase 2) [6]. The objectives of this research study are the following:
to propose a method for the validation of Phase 2 using a configurable synthetic vector field representing a generic deformation (resulting from the combination of different basic deformations, and integrating custom random noise);
to enhance Phase 2 adding a new filter to select the best feature match couple;
to show the preliminary results of the validation of the algorithm, and an example of a comparison between the original and the current versions.
I. INTRODUCTION The deep fascia is a continuous fibrous membrane that protects, enfolds and divides several anatomical parts: muscles, nerves, vessels, ligaments, and joints just to name a few. This protective membrane is organized into different layers depending on the anatomical region, and it is adaptable to the dynamics of the underlying muscles thanks to its sliding layers [1] [2]. In previous years, several studies have demonstrated that the deep fascia is of utmost importance for movement perception and coordination, muscular interaction and force transmission, tensional stress and etiology of pain [1]. Recently, evidence has been found to establish the significance of fascial layers mobility for optimal muscular activity and joint movement [1], and it has been hypothesized that any issue in the smooth sliding of fascial layers could result in improper coordinated movements [3]. In the past few years, researchers have investigated the deep fascia exploiting ultrasound (US) imaging [4], making possible to observe pathological changes of the fascial layers, and to evaluate the deep fascia physiology [5]. In a previous study [6] we have defined a novel approach to study in-vivo fascial mobility with the semiautomatic generation of 3D motion vector fields describing the displacement of salient fascial features during a muscular
G. Turini is with the Computer Science Department, Kettering University, Flint, MI 48504 USA (e-mail: [email protected]). S. Condino is with the EndoCAS Center, Department of Translational Research on New Technologies in Medicine and Surgery, University of Pisa, Ospedale di Cisanello, Via Paradisa 2, 56124 Pisa, Italy (corresponding author; phone: +39-050-99-5689; fax: +39-050-99-5689; e-mail: [email protected]). A. Stecco is with the Department of Internal Medicine, Sport Medicine Unit, University of Padua, Padua, Italy (e-mail: [email protected]).
978-1-4244-9270-1/15/$31.00 ©2015 IEEE
II. ALGORITHM DESCRIPTION This section describes our approach to analyze the deep fascia using volumetric US scans. In the proposed method, the motion of each fascial layer is studied individually by analyzing the displacement of points which correspond to clearly distinguishable anatomical features (e.g. fascia/intermuscular septum intersections). A. Phase 1: Generation of 3D Motion Field from US Scans This phase is performed through the following sequence of steps: Image Filtering to improve the visualization of fascial layers [7, 8], Regions of Interest Selection to select the target fascial layer [9, 10], Feature Detection and Extraction to identify clearly distinguishable anatomical features [11], and Feature Matching [12, 13].
G. Turini, V. Ferrari, M. Ferrari, and M. Gesi are with the EndoCAS Center, Department of Translational Research on New Technologies in Medicine and Surgery, University of Pisa, Ospedale di Cisanello, Via Paradisa 2, 56124 Pisa, Italy (e-mail: [email protected], [email protected], [email protected]). V. Ferrari is with the Information Engineering Department, University of Pisa, Pisa, Italy.
775
The most important step in this phase is the Feature Matching, which aims to find corresponding pairs of feature points in two different datasets: the first feature is located in a Reference Volume, whereas the second feature is located in a Target Volume (Figure 1). In our approach, this procedure is performed twice: the first time to find for each feature point in Volrest (Reference Volume), its corresponding feature point in Volcont (Target Volume); the second time vice versa, to find for each feature point in Volcont (Reference Volume), its corresponding feature point in Volrest (Target Volume) (Figure 1). This procedure creates two sets of feature matches: RC Matches (feature matches with Volrest as Reference Volume and Volcont as Target Volume), and CR Matches (feature matches with Volcont as Reference Volume and Volrest as Target Volume). Each RC and CR feature match is then represented as a 3D displacement vector describing the motion of the relative anatomical feature between the contraction and rest phases. The double execution of the Feature Matching has two advantages: an increased number of potential feature matches (default RC Matches are extended adding CR Matches) by using salient features in both volumes to generate a more dense motion vector field; in addition, it improves the validation of feature matches exploiting the fact that RC Matches and CR Matches have been generated independently (i.e. a strong similarity between a RC match and a CR match is a significant evidence of the reliability of both feature matches).
B. Phase 2: Validation and Filtering of 3D Motion Field The sparse motion vector field resulting from Phase 1 is then filtered by Phase 2, which consists in an automatic validation process designed to select only reliable feature matches, removing uncertain matches. This phase is split into two steps: Match Certification, and Match Recovery; and the final result consists in a classification of all the feature matches into three distinct categories: Certified Matches (highly reliable), Compatible Matches (reliable), and Conflicting Matches (unreliable) (Figure 2). The first stage of Phase 2 is the Match Certification (Figure 2a), and performs the identification of Certified Matches. For each RC Match RCm = {Rm, Cm} we search the closest CR Match CRn = {Cn, Rn}, where: Rm and Rn are 3D feature points in Volrest and Cm and Cn are 3D feature points in Volcont. The search for CRn is based on the evaluation of two distance vectors: Rnm = Rn – Rm, and Cnm = Cn – Cm. Then, the pair RCm/CRn undergoes a Distance Check (1) and a Motion Compatibility Check (2) (Figure 1 and Figure 2). ( ‖Rnm ‖ < Th1 ) ∧ ( ‖Cnm ‖ < Th1 )
(
‖ Dm - Dn ‖ ‖ Dm ‖
< Th2 ) ∧ (
‖ Dm - Dn ‖ ‖ Dn ‖
< Th2 ) ∧ (Amn < Th3 )
In (1) and (2): Dm = Rm – Cm, Dn = Rn – Cn, Amn is the angle between Dm and Dn, and Th1, Th2 and Th3 are threshold
Figure 1: An example of a generic RCm/CRn pair including the representation of all the vectors used in Phase 2. (A) and (B) show two cross sections of the Reference Volume (Volrest in this case) and Target Volume (Volcont in this case), illustrating the two feature points Rm and Cm of a RC feature match RCm. (D), (E), and (F) illustrate the reference systems used in the Rest Volume Vol rest and in the Contraction Volume Volcont, and how these are unified to compute the candidate displacement vectors Dm and Dn.
776
Figure 2: Diagram of Phase 2 illustrating the algorithm to validate feature matches: (A) step 2.1 Match Certification is used to select highly reliable matches (Certified Matches), whereas (B) step 2.2 aims to exploit the Certified Matches in order to recover previously discarded feature matches, improving in this way the density of the final motion vector field.
values. Equation (1) ensures that RCm and CRn are located in the same anatomical region, whereas (2) ensures that RCm and CRn are similar in terms of their normalized modulus and direction. If these two feature matches don’t satisfy one of these checks (Distance Check and Motion Compatibility Check) they are temporarily skipped, and will be analyzed again in the second stage of Phase 2. Otherwise, if both checks are passed, the couple RCm/CRn is inserted in a temporary set containing all the couples RCm/CRk...r that have passed both checks. Then, we use a cost function fcost (3-9) to evaluate each couple in this set in order to select the best CR match CRn. Finally we classify both RCm and CRn as Certified Matches, and repeat the process for the next RC match.
fcost ( RCm / CRn ) = ∑6i=1 vi wi ,
v1 =
v2 =
| ‖ Rmn ‖ - Th1 |
| ‖ Cmn ‖ - Th1 |
Th1 Th1
v3 = max ( 0, 1 - |
v4 =
v5 =
|
‖ Rmn ‖ ‖ Cmn ‖
‖ Dm - Dn ‖ - Th2 | ‖ Dm ‖
Th2 |
∑6i=1 wi =1
‖ Dm - Dn ‖ - Th2 | ‖ Dn ‖
Th2
| )
v6 =
| Amn - Th3 | Th3
The cost function fcost in (3) evaluates the quality of a feature match couple by computing a weighted arithmetic mean of six different properties: the distance (Rmn) between the feature points in Volrest (Rm, and Rn) normalized on a threshold value (4), the distance (Cmn) between the feature points in Volcont (Cm, and Cn) normalized on a threshold value (5), the ratio between the previous two distances (Rmn, and Cmn) (6), normalized ratios of the displacement vectors Dm and Dn (7) (8), and angle between RCm and CRn normalized on a threshold value (9). This improved version of the algorithm includes new stages (in bright red in Figure 2) in order to select the best CR Match to pair with the input RC Match RCm. In the original version of the algorithm we considered only the closest CR Match to RCm. Now, instead, we collect a group of CR Matches close to RCm, and then we choose the best CR Match CRn evaluating the quality of all the matches in the group using the cost function fcost (3). Match Recovery is the second stage of Phase 2, and starts selecting all the feature matches which have not been certified: from now on called Non-Certified Match (NCM). Then, for each NCMi (with its relative Ri and Ci), the closest Certified Match CMj (with its relative Rj and Cj) is searched evaluating the distance vectors: Rij = Ri – Rj, and Cij = Ci – Cj. NCMi is
777
then classified as Compatible Match if it satisfies the Motion Compatibility Check as rephrased in (3).
(
‖ Di - Dj ‖ ‖ Di ‖
< Th4 ) ∧ (
‖ Di - Dj ‖ ‖ Dj ‖
< Th4 ) ∧ (Aij < Th5 )
In (11): Di = Ri – Ci, Dj = Rj – Cj, Aij is the angle between Di and Dj, and Th4 and Th5 are threshold values. Otherwise, if (10) is not satisfied, NCMi is classified as Conflicting Match and discarded. At the end, reliable matches (Certified Matches and Compatible Matches) are used for fascial motion estimation. III. TESTING METHODOLOGIES The main goal or our tests is the evaluation of the performance of our algorithm on a configurable deformation vector field, enabling a per-vector analysis of the final results. To do this, we generated a custom synthetic sparse deformation vector field which simulates the deformation of a fascial layer, and then we processed this vector field using our algorithm in order to obtain a validated set of deformation vectors discarding outliers. A. Deformation Modeling and Sparse Vector Field Generation To generate a generic 3D deformation vector field, we computed a deformation vector for each voxel of an input volumetric dataset. So, we multiplied the 3D position vi of each voxel i in the volumetric dataset by a deformation matrix L, obtaining the 3D deformation vector di applied on the voxel i (11).
vi L = di
In (11), di is the displacement vector of the feature located in voxel i, and it simulates a feature match as if generated by Phase 1. In this case, we approximated a generic three– dimensional deformation by a linear 3D vector field described by the matrix L [14, 15]: L11 L12 L13 1 0 0 L = [L21 L22 L23 ] = E [ 0 1 0 ] + L31 L32 L33 0 0 1 1 0 0 1 0 0 1 0 + R [ 0 0 -1 ] + S1 [ 0 1 0 ] + S2 [ 0 0 0 1 0 0 0 -1 0 1 where
E= S1 =
L11 + L22 + L33 3 L11 + L22 - L33 3
R= S2 =
0 1] 0
fascial layer. Now, to simulate the output of Phase 1, two different random selections of displacement vectors are performed to create the two sets of feature matches: RC Matches and CR Matches, the input required by Phase 2 (Figure 2). Then, the resulting sparse vector field (including both RC and CR matches) is affected by random noise in order to test Phase 2 robustness. B. Experimental Setup The generation of both the dense and sparse vector fields, and the entire algorithm (both Phase 1 and Phase 2) were performed in Matlab® version 7.14, release 2012a (MathWorks, Natick, Massachusetts, USA). The size of the input volumetric datasets used to test the algorithm are 5×370×100 voxels (Sx×Sy×Sz) and the density of the simulated sparse deformation field is 0.05% (for a total of 92 vectors per field). These values were chosen to obtain a deformation field similar to those generated analyzing real fascial volumes as described in our previous work [6]. The thresholds Th1, Th2, Th3, Th4, and Th5 of Phase 2 were set to 2.89 voxels, 0.4, 20°, 1, and 40°, respectively. These values are those preliminarily set in [6] to optimize the algorithm performance on the acquired datasets. The weights used in the cost function fcost (3) were set as follows: w1=w2=w4=w5=0.1; w3= w6=0.3. Phase 2 was tested processing 15 different deformations. The scalar deformation coefficients E, R, S1, and S2 were set to 0 or 0.1. All the possible combinations of the elementary components (expansion, rotation, horizontal shear, and oblique shear) were tested as illustrated in Figure 3. During the sparse vector field generation, we applied a random error to each displacement vector adding a 3D vector with x, y, and z coordinates ranging respectively from 0 to 2% of Sx, Sy, and Sz. For each of the 15 deformations tested, we processed 10
L11 - L23 + L32 3 L11 + L23 + L32
3
In (12), L is decomposed as a sum of four matrices each representing an elementary deformation component (13): expansion, rotation, horizontal shear, and oblique shear (Figure 3). Each deformation component (matrix) is controlled by a single scalar coefficient: E, R, S1, and S2 respectively [14, 15]. Equation (12) describes how our synthetic deformation field is generated, that is: the 3D position of (the centroid of) each voxel v is multiplied by matrix L, and the result is the 3D deformation vector d affecting that voxel, as described in (11). The vector field obtained so far is a single dense deformation field that represents the deformation affecting a 778
Figure 3: An example of the four different elementary deformations used to compose the generic deformation matrix L: expansion (A), rotation (B), horizontal shear (C), and oblique shear (D).
different sparse 3D vector fields, and the results were summarized in Table 1 as: the ratio between the modulus of the error vector and the modulus of the relative original displacement vector (without noise) expressed as a percentage (Mod), and the angle between the error vector and the relative original displacement vector (Ang). Then, we performed the same tests increasing the maximum error percentage of 1% each time: from 2% to 8%. The diagram in Figure 5 illustrates how the error affects the classification in terms of the percentages of Certified, Compatible, and Conflicting Matches. IV. PRELIMINARY RESULTS As expected, values reported in Table 1 show that the error (Mod and Ang) affecting Certified and Compatible Matches is always lower that the error of Conflicting Matches. The only exception is relative to a pure rotational deformation (row 4, gray cells), for which we have an error in the angle slightly greater for Compatible Matches. Figure 4 shows an example of a section of a 3D deformation vector field as resulting after Phase 2, and illustrates: Certified, Compatible (both RC and CR) and Conflicting Matches with the relative error vector. In the same figure, black circles show a pair of feature matches classified by the current improved version of Phase 2 and the original version of our algorithm: highlighting that Phase 2 is now able to classify as Certified Matches two displacement vectors previously classified as Compatible and Conflicting Matches. ERROR OF CLASSIFIED MATCHES AS DEFORMATION CHANGES
Deformation Type
Error of Classified Matches Certified
Compatible
Conflicting
E
R
S1
S2
Mod (%)
Ang (°)
Mod (%)
Ang (°)
Mod (%)
Ang (°)
0
0
0
0.1
16.2
41.5
26.8
64.4
45.4
65.2
0
0
0.1
0
25.9
42.6
32.8
53.4
37.1
70.8
0
0
0.1
0.1
14.8
45.4
18.8
56.6
26.3
62.2
0
0.1
0
0
18.6
44.5
25.1
58.7
41.6
57.1
0
0.1
0
0.1
11.7
50.9
14.2
57.8
125.9
58.5
0
0.1
0.1
0
17.9
44.6
34.4
55.8
79.9
59.3
0
0.1
0.1
0.1
10.0
50.8
11.9
53.3
40.4
60.4
0.1
0
0
0
23.0
39.4
31.4
50.3
38.2
68.0
0.1
0
0
0.1
12.2
51.5
19.7
55.1
87.2
66.4
0.1
0
0.1
0
20.9
46.8
25.3
47.8
98.2
80.4
0.1
0
0.1
0.1
10.6
55.3
18.1
54.1
46.9
64.3
0.1
0.1
0
0
13.8
43.5
19.4
56.3
30.0
59.9
0.1
0.1
0
0.1
15.4
51.5
13.2
55.1
50.3
64.4
0.1
0.1
0.1
0
11.3
52.6
15.7
55.4
40.7
59.3
0.1
0.1
0.1
0.1
10.7
49.6
14.1
58.9
86.8
59.9
Analyzing the classification of feature matches as the maximum error increases, it is clear (as illustrated in Figure 5) that the percentage of Conflicting Matches increases accordingly. This happens because while increasing the error, we are at the same time decreasing the percentage of displacement vectors congruent with the original deformation vector field. V. DISCUSSION In conclusion, this paper presents an improvement of the method proposed in [6] to study in-vivo the mobility of the deep fascia using a semiautomatic approach. In particular, this study is focused on the final filtering of the 3D sparse vector field representing the deformation of a fascial layer, proposing an enhanced algorithm and a method for its validation. Our algorithm has been tested on different synthetic deformation vector fields and as, preliminarily demonstrated by the results obtained, it works for deformations including expansion, rotation, horizontal shear, and oblique shear components. All our tests were performed using a fixed configuration of the algorithm settings. MATCH CLASSIFICATION AS NOISE CHANGES
Mean Match Type (%)
TABLE 1
Figure 4: A section of a 3D deformation vector field as resulting after Phase 2. The black circles show the comparison between the current improved version of Phase 2 (top) and the original version of the same stage (bottom).
100% 80%
Conflicting Matches
60% 40%
Compatible Matches
20% 0% 2
3
4
5
6
7
8
Certified Matches
Max Noise Amplitude (%) Figure 5: Differences in the classification of feature matches as the maximum noise applied increases from 2% to 8%. Each value illustrated represents the average percentage over 15 different deformations each tested using 10 different sparse vector fields.
779
In the future, we plan to extensively test the proposed method, varying the configuration of the deformation and the applied error in order to optimally tune the parameters of our algorithm. ACKNOWLEDGMENT The authors would like to thank Eng. Simone Parrini for his essential contribution throughout the initial research project that has led to this study. REFERENCES [1] C. Stecco, V. Macchi, A. Porzionato, F. Duparc and R. De Caro, "The fascia: the forgotten structure," Italian Journal of Anatomy and Embryology, vol. 116, no. 3, pp. 127-138, 2011. [2] G. Natale, S. Condino, P. Soldani, F. Fornai, M. Mattioli Belmonte and M. Gesi, "Natale et. al.'s response to Stecco's fascial nomenclature editorial," Journal of Bodywork and Movement Therapies, vol. 18, no. 4, pp. 588-590, 2014. [3] C. Stecco and W. I. Hammer, Functional atlas of the human fascial system, Elsevier Health Sciences, 2014. [4] H. M. Langevin, D. Stevens-Tuttle, J. R. Fox, G. J. Badger, N. A. Bouffard, M. H. Krag, J. Wu and S. M. Henry, "Ultrasound evidence of altered lumbar connective tissue structure in human subjects with chronic low back pain," BMC Musculoskeletal Disorders, vol. 10, no. 1, p. 151, 2009. [5] A. Stecco, A. Meneghini, R. Stern, C. Stecco and M. Imamura, "Ultrasonography in myofascial neck pain: randomized clinical trial for diagnosis and follow-up," Surgical and Radiologic Anatomy, vol. 36, no. 3, pp. 243-253, 2014. [6] S. Condino, G. Turini, S. Parrini, A. Stecco, F. Busoni, V. Ferrari, M. Ferrari and M. Gesi, "A semiautomatic method for in vivo three-dimensional quantitative analysis of fascial layers mobility based on 3D ultrasound scans," International Journal of Computer Assisted Radiology and Surgery (IJCARS), pp. 1-15, 2015. [7] A. F. Frangi, W. J. Niessen, K. L. Vincken and M. A. Viergever, "Multiscale vessel enhancement filtering," in Medical Image Computing and Computer-Assisted Intervention (MICCAI), 1998. [8] M. Rana and J. M. Wakeling, "In-vivo determination of 3D muscle architecture of human muscle using free hand ultrasound," Journal of Biomechanics, vol. 44, no. 11, pp. 2129-2135, 2011. [9] P. A. Yushkevich, J. Piven, H. C. Hazlett, R. G. Smith, S. Ho, J. C. Gee and G. Gerig, "User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability," Neuroimage, vol. 31, no. 3, pp. 1116-1128, 2006. [10] V. Ferrari, C. Cappelli, G. Megali and A. Pietrabissa, "An anatomy driven approach for generation of 3D 780
models from multi-phase CT images," International Journal of Computer Assisted Radiology and Surgery (IJCARS), vol. 3, no. 1, pp. 272-273, 2008. [11] J. Shi and C. Tomasi, "Good features to track," in Proceedings of IEEE Computer Vision and Pattern Recognition (CVPR), 1994. [12] O. Commowick, N. Wiest-Daesslé and S. Prima, "Block-matching strategies for rigid registration of multimodal medical images," in Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI), 2012. [13] A. Gastounioti, S. Golemati, J. S. Stoitsis and K. S. Nikita, "Carotid artery wall motion analysis from Bmode ultrasound using adaptive block matching: in silico evaluation and in vivo application," Physics in Medicine and Biology, vol. 58, no. 24, p. 8647, 2013. [14] D. Zoccolan and V. Torre, "Using optical flow to characterize sensory-motor interactions in a segment of the medicinal leech," The Journal of neuroscience, vol. 22, no. 6, pp. 2283-2298, 2002. [15] A. Giachetti and V. Torre, "The use of optical flow for the analysis of non-rigid motions," International Journal of Computer Vision, vol. 18, no. 3, pp. 255279, 1996.
