Quantifying motion artifacts in 4D-CT using edge strength analysis

(1)

Free University/ University of Amsterdam

Bachelor Thesis

Physics & Astronomy (joint degree)

Size: 15 EC, conducted between March 30 and June 25, 2020 Faculteit der Natuurwetenschappen, Wiskunde en Informatica (FNWI)

Quantifying motion artifacts in 4D-CT using edge

strength analysis

N.J.J. Deen

Supervisor: Dr. ir. G.J. Streekstra Daily Supervisor: ir. J.G.M. Oonk Examiner: Prof. Dr. ir. G. J. Strijkers

Biomedical Engineering and Physics

Amsterdam UMC

(2)

Abstract

Early diagnosis of wrist joint motion pathologies is considered beneficial for the health of patients, since conditions such as arthritis can be treated in an earlier phase. 4D-CT imaging, which is a sequence of 3D-4D-CT scans in time, provides the ability to closely inspect wrist motion making it easier to detect motion related pathologies which remain hidden in regular CT imaging.

For accurate wrist motion analysis, segmentation and registration of the carpal bones is required, which are time consuming processes. Imaging is subject to motion arti-facts, in particular motion blur, which are known to affect the liability of registration results. Localization of the cortical layer is often part of the segmentation and regis-tration process. In this study, a method to quantify motion blur is introduced, that uses both first- and second order Gaussian derivative analysis to calculate the edge-and ridge strength of the cortical layer of carpal bone, using data of a moving wrist phantom. Both methods are evaluated through visual assessment of 3D edge/ridge strength visualizations. Additionally, statistical analysis was performed to check for a correlation between edge/ridge strength behaviour and wrist motion speed. Finally, the effect of Gaussian smoothing on the performance of the model was investigated for three smoothing levels.

It was found that at a low smoothing level, both edge- and ridge strength showed a significant, negative correlation with wrist motion speed. Medium smoothing only showed a significant correlation for the ridge strength method and the highest level of smoothing led to no correlations. This indicated that low smoothing levels provided the most significant results. Visual assessment concluded that at low smoothing invalid edges and noise were affecting the reliability of edge- and ridge strength cal-culations. Medium smoothing was considered the best option with a sharp reduction in noise and invalid edges, while distortion of the edges representing the cortical layer of the carpal bones was limited. The highest smoothing level led to a distortion of the edge of the cortical layer of carpal bones and was not advised. Results following from Statistical analysis and visual assessment were contradicting. It was therefore concluded that low smoothing should be applied for best performance of edge/ridge strength analysis, but that these results might not be entirely representative for the edge/ridge strength of the cortical layer due to the presence of invalid edges in this setup at low smoothing. Furthermore, it was found that the ridge strength model was more sensitive to a change in smoothing than the edge strength model. Both models were more sensitive to a change in smoothing around a sharp edge than around a smooth edge, indicating more distortion around sharp edges as more smoothing is applied. Automatic scale selection is recommended as an extension to this model, since it will improve the ability to suppress invalid edges and noise, while maintaining the natural structure of the carpal bone edges.

(3)

Populair Wetenschappelijke Samenvatting Voor Scholieren

Aandoeningen in het polsgewricht, zoals arthritis, worden vaak pas vastgesteld wan-neer het al te laat is. Dit komt doordat dit soort aandoeningen in de beginfase vaak alleen nog waar te nemen zijn door de beweging van de polsgewrichten onder de loep te nemen. Bij aanhoudende klachten in de polsgewrichten wordt vaak een CT-scan gemaakt. Door middel van een CT-scan kan een 3D reconstructie verkregen worden van de botten. Hieronder is een voorbeeld te zien van zo’n 3D reconstructie. Echter, een CT-scan is een momentopname en dus is het niet mogelijk hierop bewegingen waar te nemen. Daarom wordt er nu onderzoek gedaan in de ontwikkeling van 4D-CT-scan technieken. Een 4D-CT scan is een reeks 4D-CT-scans in de tijd waarmee een 3D filmpje gemaakt kan worden van de botten en bewegingen dus w´el kunnen worden waargenomen en geanalyseerd. Wanneer tijdens het maken van een 4D-CT scan te snel wordt bewogen zullen de beelden vervagen, waardoor het lastiger wordt een 3D reconstructie te vormen. Daardoor wordt het analyseren van de bewegingen bemoeil-ijkt of zelfs onmogelijk. Aangezien het analyseren van de bewegingen veel tijd kost, zou je graag van te voren willen weten of je met de gemaakte scans een nauwkeurige weergave van de beweging kunt cre¨eren. De betrouwbaarheid van deze resultaten hangt nauw samen met de hoeveelheid vervaging in de beelden. Daarom is het doel van dit onderzoek geweest om deze vervaging te kwantificeren. Dit is gedaan door een algoritme te bouwen die de beelden van 4D-CT scans verwerkt en dit algoritme vervolgens te testen en evalueren.

Figure 1: Voorbeeld van een 3D-reconstructie van de botten in het polsgewricht, ook wel de carpale botjes genoemd. Deze reconstructies worden verkregen door segmentatie van de botten uit een CT-scan. Adapted from [1]

(4)

1 Introduction

The use of 4D-CT-imaging, a sequence of 3D-CT scans in time, is considered to be a promising method for diagnosing joint motion pathology. Although this method is hardly being used clinically yet, it could contribute to detecting pathologies in an earlier stage. In previous wrist motion studies for example, it was found that Osteoarthritis and carpal motion instability are highly correlated [2], and that early detection and treatment of carpal motion instability could limit the consequences for the patient’s health in the future [3]. Dynamic wrist imaging using 4D-CT is described as a promising method to detect these instabilities in their early phase, because of the ability to detect pathology that only emerges during wrist movement that could not or would be harder to detect with methods using static CT- or vide-ofluoroscopy scans [4].

For accurate carpal motion analysis, segmentation and registration of the bones is required, where segmentation aims to distinguish relevant structures from an image for visualization and processing and by registration these structures can be mapped on itself to align scans with each other [5]. Registration can be applied to align scans of different types, such as in CT-MR imaging [6], or to align scans in time so that the motion of image structures can be monitored. This last application of registration is used in carpal motion analysis. Dobbe et al. have introduced a semi-automatic seg-mentation and registration method to analyze 3D carpal bone motion and evaluated the accuracy of this method, by measuring the translation and rotation registration errors for multiple constant rotation speeds of a wrist phantom [7]. In this study the effect of motion artifacts on the bone registration error is found through varying the rotation speed of the phantom. Registration is a time consuming process, and motion blur affects the accuracy of this process significantly [7]. In other studies on motion artifact in 4D-CT it was found that an increased motion speed was correlated to a higher degree of observed motion blur [8]. Furthermore, streaking- and shading artifacts were observed, besides blurring [9].

Ideally, in clinical environments, one would want to know if segmentation and reg-istration will give realistic and usable results before proceeding with such tedious processes. For example, if the quality of the images resulting from the scan seems in-sufficient to be used for accurate bone registration, a new scan can be made instantly with different patient instructions or scanner settings. To determine if segmentation and registration is feasible, a method for assessment of image quality could be devel-oped. Using this method, the registration- and segmentation process for inadequate scans could be avoided and time would be saved. Besides that, it would not be nec-essary for the patient to return to the clinic for a second scan.

Edge strength analysis has been studied in the past for acquiring a measure for image degradation. In Zhang et al. , the edge strength is used as a measure for the probability of a pixel belonging to the edge of an object. The goal in this study was to mimic the (subjective) human visual system, which is the way how humans would rate the quality of an image based on distortion of information. Distortion is

(6)

usually better visible at the edges, which is defined as a set of pixels belonging to the boundary of an object [10]. In gray-scale images, the pixel intensity changes sharply at these boundaries, so that the edge can be recognised as a border between areas of different intensities. Edge strength analysis has also been proven to be effective in registration of combined scans such as CT-MR brain images [6]. In multiple studies, bones are segmented and registered by locating the cortical layer (outer edge) [7] [11] [12] [13], which can be identified in a CT-scan as a pointset shaped as a thin layer of bright points, enclosed by darker points representing the soft tissue surrounding the bone. On the basis of this definition, the cortical layer of bone could be identified both as a ridge and an edge, in which a ridge is defined as a thin line of pixels with a different intensity than the pixels surrounding it. Image distortions, that appear in case of motion-induced artifacts, are mostly present at the edges. This could imply a decrease in edge strength, since it will be harder to track down the pixels that belong to edges due to these distortions. This could have consequences for the quality of bone registration and -segmentation.

The aim of this study is to elaborate a method that quantifies motion blur in 4D-CT scans, by using edge strength analysis on the cortical layer of carpal bone data. Since the cortical layer of bone can be regarded as both an edge and a ridge, an edge detection and ridge detection method will be used and evaluated for this situation, using first- and second order Gaussian derivative analysis respectively. These meth-ods could be used as an extra tool alongside visual assessment to quantify motion blur faster and more systematically.

2 Theory

2.1 Edge strength

In Zhang et al. edge strength is defined as a measure for the probability of a pixel belonging to the edge of an image object. Furthermore, it is stated that an edge always shows discontinuity in intensity in a certain direction (moving across the edge), while showing continuity orthogonal to that direction (moving along the edge). As an extreme example, an image that would show discontinuity in both directions would be an image with only noise and an image that would show continuity in both directions would be an image with constant intensity. These factors combined were used to detect and measure the strength of salient edges, by defining the edge strength of a pixel as the absolute difference between the directional image derivative in direction u and the directional image derivative in direction v, or S = |DuI − DvI|

[10], where u and v are 2D unit vectors, u ⊥ v and

DuI = ∇I · u =

∂I ∂xu1+

∂I

∂yu2 (1)

(7)

2.2 Directional Gaussian derivatives

For calculating the directional first order derivatives of the images Gaussian kernels can be used. Gaussian kernels are commonly used in image processing operations to apply smoothing. By convolving an image I with the directional derivative of a Gaussian kernel, such that

∂ ∂i G(σ) ∗ I = ∂ ∂iG(σ) ∗ I, [i = x, y], G(σ) = 1 2πσ2e −x2+y2 2σ2 , (2)

a smooth derivative image is obtained. The Gaussian kernel hereby acts as a noise suppressor in the calculation of the derivatives. This method turned out to be a pow-erful tool in calculation of derivative images, providing meaningful derivative values [14] [15]. Directional derivatives are obtained by combining equations 1 and 2 for directions u1 = [1 0] 0 ⊥ v1 = [0 1] 0 and u2 = √1₂[1 1] 0 ⊥ v2 = √1₂[−1 1] 0 . A visual representation of these directions can be found in the Appendix. With these direc-tional derivatives, two edge strengths (S1 and S2) can be calculated at each pixel,

where S1 represents the edge strength in the horizontal/vertical direction and S2 the

edge strength in diagonal directions. Finally, the total edge strength magnitude in that pixel will then be defined as Edgestrength =p(S1)2+ (S2)2

For second order derivative analysis of images, the eigenvalues of the Hessian matrix can be used. The values of the Hessian matrix, H, must be calculated for each pixel, where H =Ixx Ixy Iyx Iyy , Imn= ∂

∂m∂nG(σ) ∗ I, [mn = xx, xy, yx, yy]. (3) Values of the Hessian are also calculated using Gaussian derivatives of the image, Imn, to suppress unwanted noise. It is expected that this method is more sensitive

to noise, since second order derivative filters are used [16]. The absolute value of the largest eigenvalue of the Hessian corresponds to the eigenvector moving accros the ridge, which is the direction where the second order derivative takes on its largest absolute value[17]. The absolute value of the smallest eigenvalue corresponds to the eigenvector moving along the ridge, which is the direction where the second order derivative takes on its smallest absolute value [14]. Similar to edge analysis, the in-tensity profile must be discontinuous across the ridge and continuous along the ridge. The difference of the absolute largest eigenvalue with the absolute smallest eigenvalue is defined as a measure for Ridge strength in this study: Ridgestrength = |λ1| − |λ2|,

with λ1 the largest and λ2 the smallest eigenvalue of the Hessian Matrix.

The smoothing operator is denoted by σ. By adjusting this parameter, the desired amount of smoothing can be applied to suppress noise, prior to derivative operations. A higher σ results in more smoothing. However, if too much smoothing is applied valuable image information might be erased along with the unwanted noise. It can lead to poor localization of the edges as well, since increased smoothing can cause shape distortions [18]. Gaussian-smoothing an image also increases the variance of

(8)

the pixel values, since neighbouring pixels contribute to the value of the center pixel of the convolution kernel. This leads to a spacing out of the edges in the Gaussian smoothed image. According to definitions above, this will cause a decrease in edge-and ridge strength at edges in general. Especially at sharp edges this effect is ex-pected to be large. Hence, the value of σ can have a significant influence on the edge/ridge strength results.

3 Methods

A subset of the data from the research of Dobbe et al. into carpal kinematics was used, since the data needed to be cropped for proper isolation of carpal bones and this was not achievable for the entire dataset. Subsequently, edge- and ridge strength cal-culations were performed using this data. These calcal-culations were done three times with different amounts of smoothing, by adjusting the σ of the Gaussian filters. More explanation about these processes will follow below. Finally, the correlation between motion blur and the change in edge/ridge strength was investigated. Multiple linear regressions were executed for this, each time with edge- or ridge strength as the de-pendent variable and the rotation speed of the phantom as the indede-pendent variable. Since 3 different amounts of smoothing were applied for both edge- and ridge strength calculations separately, 6 linear regressions were performed. How these edge/ridge strengths were calculated will be elaborated below. Since σ was expected to have influence on the edge/ridge strength, two regressions with the edge/ridge strength as dependent variable and σ as independent variable were also executed. All linear regressions were performed using R (version 3.5.2. The R Foundation for Statistical Computing, Vienna, Austria). With the linear models that were obtained a test for correlation between the variables of interest could be performed, with the null hy-pothesis being that there was no correlation and rejection of the null when p<0.05

Data Structure and Processing

All data originates from the department of Biomedical Engineering in the Amsterdam UMC (location AMC). Image acquisition was performed with a Siemens Somatom Force scanner (120 kV, 75 mAs, slice thickness 0.6 mm, gantry rotation time 0.25s, with a voxel size of 0.43 X 0.43 X 0.30 mm [7]). Each 4D-CT scan consisted of 32 time frames, which were all 3D-CT scans that were reconstructed out of 191 DICOM images (512 X 512). In this study 3 4D-CT scans of a wrist phantom were used with phantom rotation speeds 0.01RPS, 0.2RPS and 0.4RPS. A cross section of the wrist phantom, with additional info about the setup and axis of rotation, is shown in figure 2. The data was first cropped in imageJ (National Institutes of Health, Laboratory for Optical and Computational Instrumentation, University of Wisconsin), using the ’freehand selection’ option. Cropping was done such that the Capitate could be iso-lated. This particular carpal bone was chosen because it was expected to be the bone that is most subjected to motion blur. It was positioned the furthest from the axis of rotation and therefore was moving at the highest velocity in this setup. Isolation of

(9)

the Capitate could only be achieved for a subset of the data, due to the orientation of the phantom. After cropping, the files were exported as TIF-files and loaded into programming software, were they were converted to 512 x 512 arrays of type double for further processing.

Figure 2: Cross section of the wrist phantom with cadaveric carpal bones that was used to acquire the 4D-CT data. The top bone in this image is the Capitate, the lower left bone is the Scaphoid and the lower right bone is the lunate. The phantom was indirectly attached to a stepper motor to give it a constant angular velocity. In this setup, the axis of rotation was placed from left to right through the center of the Lunate, such as the red line in this figure. Adapted from [7].

MATLAB code

MATLAB R2018b (Mathworks, Natick, MA, USA) has been used for programming. edge- and ridge strengths were calculated using custom code. Reading images, con-volving images with Gaussian kernels and displaying 3D edge/ridge strength repre-sentations were done using the Image Processing Toolbox (version 10.3, Mathworks, 2018) and MATLAB’s standard built-in functions.

The sum of all pixel edge/ridge strengths in a slice resulted in the 2D edge- or ridge strength of that slice. Since the data was cropped first, all edge/ridge strength contribution was expected to be mostly from the cortical layer of the Capitate. The contribution of invalid edges and ridges, originating from motion artifacts caused by the phantom surrounding the Capitate, was believed to be small as a result from the cropping step preceding the edge/ridge strength calculations. Overall, these were expected to have a much smaller edge/ridge strength than that of the cortical layer pixels, and thus not have a significant contribution to the sums. Next, the edge/ridge strength for the 3D reconstructions was determined by averaging the 2D edge/ridge strengths of the slices.

For the 3D edge/ridge strength visualisations, a colored heatmap with intensity scale was defined. This intensity scale included a blue color representing weak-, a green and yellow color representing stronger-, or even an orange/red color representing the strongest edges/ridges. This intensity scale was universally defined, so that the colors represent the same edge/ridge strength intensities in every visualisation. This was

(10)

done by calculating a baseline for the maximum intensity. For this baseline the mean of the maximum intensities of the 0.01rps 4D-CT was used, since these scans were expected to provide the highest edge/ridge strengths due to the motion blur being nil. A baseline was created for both the edge- and ridge strength method. After setting the baseline for the maximum intensity, a linear intensity scale was created to divide the intensities over 256 colors. A universal opacity map was also included, to make the lower intensities transparent and the higher intensities opaque. This was done to improve the visibility of the carpal bones in the representations. After all, lower intensities were expected not to provide a large contribution to the total edge/ridge strengths.

Edge- and Ridge strength calculations were performed for σ = 1, 2 and 3. Gaus-sian convolution kernels of 20 x 20 pixels were used. By consulting the 3D edge strength visualizations, an assessment could be acquired for which σ this method performed best. By adjusting σ a sweet spot was searched, with just the right amount of smoothing to get rid of unwanted information (such as invalid edges and noise), without altering the true nature bone edges severely to provide meaningfull edge/ridge strength data.

To gain insight into how structures of edges/ridges are distorted as more smoothing is applied on a more local scale, it was investigated how sensitive the ridge- and edge strength models are to a change in σ around both sharp and smooth edges. Edge-and ridge strength analysis was performed on an image with a simulated sharp edge and an image with a simulated smooth edge, both of dimension 256 x 256. This anal-ysis was done for the values σ = 1, 2 and 3. Next, it was noted for both edges how much the total edge/ridge strength sum would change for σ = 2 and σ = 3, relative to the edge/ridge strength sum for σ = 1. The sharp- and smooth edge image were both generated in MATLAB using custom code and can be found in the Appendix.

(11)

4 Results

From the 4D-CT scan in which the phantom rotated at 0.01rps, nine time frames were usable for cropping and further analysis. In the 4D-CT in which the phantom rotated at 0.2rps and 0.4 rps, six and seven time frames were usable, respectively.

The average edge/ridge strength distributions of the DICOM images in the 3D recon-structions were calculated for each time frame. The medians of these mean edge/ridge strength distributions through time can be found in table 1 for each rotation speed and σ = 1, 2 and 3. Boxplots of the full distributions can be found in figures 5, 6 and 7 for σ = 1, 2 and 3 respectively.

0.01 rps 0.2 rps 0.4 rps σ = 1 Ridge 60713 59847 46848 Edge 62909 63443 48041 σ = 2 Ridge 8123 8178 7123 Edge 19324 21007 17470 σ = 3 Ridge 3275 3424 3095 Edge 12661 14419 12216

Table 1: Median values of the average slice edge/ridge strength distributions

By using the 3D edge/ridge strength visualizations, it was determined which value for σ could be used best judging by visual assessment. It was found that σ = 2 gave the best results, filtering away most of the noise and unwanted edges/ridges in the phantom without severely altering the carpal bone structures. See for example the edge- and ridge strength visualizations of the 0.01RPS 4D-CT at timestep 27 with σ = 1, 2 and 3 in figures 3 and 4. All edge- and ridge strength visualizations (each timestep and rotation speed) can be found in the Appendix.

For the method that calculated the edge strength using first order directional Gaus-sian derivatives, a significant negative correlation was found between the average slice edge strength and rotation speed of the phantom (p < 0.01) for σ =1. Calcu-lations performed with σ = 2 and σ = 3 did not provide any significant correlation between the edge strength and rotation speed of the phantom (p >> 0.05). For the method that calculated the ridge strength, using eigenvalues of the Hessian matrix, a significant negative correlation between the average slice ridge strength and rota-tion speed of the phantom was found for σ = 1 (p < 0.01) and σ = 2 (p < 0.05). Calculations performed with σ=3 did not provide any significant correlation between ridge strength and rotation speed (p >> 0.05). A significant, negative correlation

(12)

was also found between both edge- and ridge strength and σ, with (p < 0.01) in both regressions. Regression results follow from a Wald t-test. The R output of all regressions, as well as the scatter plots of edge/ridge strengths against σ, can be found in the Appendix.

By calculating the edge- and ridge strength sums at simulated edges, it was in-vestigated how sensitive both methods were to a change in σ around a sharp and smooth edge. For the edge strength method, a change in the edge strength sum of -3,86 % and -7,61 % was found for σ = 2 and σ = 3 around a sharp edge. Around a smooth edge the changes in edge strength were -0,35 % and -3,13 %, for σ = 2 and σ = 3, respectively. For the ridge strength method, a change in the ridge strength sum of -53,63 % and -70,42 % was found for σ = 2 and σ = 3 around a sharp edge. Around a smooth edge these values were -21,31 % and -36,07 % for σ = 2 and σ = 3, respectively. All changes in edge/ridge strength were measured relative to the values for the edge/ridge strength found for σ = 1.

Figure 3: Edge strength visualizations of the 0.01RPS 4D-CT at t27 with σ = 1 (left), σ = 2 (center) and σ = 3 (right)

Figure 4: Ridge strength visualizations of the 0.01RPS 4D-CT at t27 with σ = 1 (left), σ = 2 (center) and σ = 3 (right)

(13)

(a) Average slice edge strength for σ = 1 (b) Average slice ridge strength for σ = 1

Figure 5: The average slice edge- and ridge strengths of the cropped 3D reconstruction of the Capitate, distributed over time and presented for three different rotation speeds.

(14)

5 Discussion

Part of this study involved evaluating what the effect of σ was on the quality of the results. In this quality review noise reduction was desired, while distortion of the edges must be avoided as much as possible. Through visual assessment of edge/ridge strength visualizations, it was found that a value of σ = 2 provided the best re-sults, out of the values for σ that have been considered. When calculations were performed with σ = 1, too much noise and edges/ridges were found in the phantom. These were not of interest and would be undesired in calculating the average slice edge/ridge strengths, since these calculations are performed to acquire information about distortions in the cortical layer. On the other hand, when σ = 3 was used for calculations, too much smoothing occurred and the edges of the Capitate were visually deformed. However, this kind of evaluation remains subjective, which should be kept in mind.

From the correlation test of the average edge/ridge strength sum against σ, it can be concluded that σ has a significant influence on the edge/ridge strength results, that expresses itself in a negative correlation. A possible explanation might be found in the spacing out effect of the edges, due to the behaviour of the Gaussian smoothing convolution kernel, and the effect of this phenomenon on the edge/ridge strength calculations. On the other hand, an explanation for this strong correlation could be searched for in the edge/ridge strengths of invalid edges that are found in the phan-tom. More and stronger invalid edges are detected at a low smoothing level than at higher smoothing levels. This would lead to a significant contribution in the aver-age edge/ridge strength of the slices at a low smoothing level, causing the negative correlation between edge/ridge strength and σ. More locally, both models showed a higher sensitivity to a change in σ around sharp edges, relative to smooth edges,

(15)

and thus implies a higher distortion of edge structures around sharp edges as more smoothing is applied. The cortical layer of the Capitate also has the properties of a sharp edge, making it more likely that increased amounts of smoothing would lead to undesired distortions in this region. This was also observed in visual assessment at σ = 3. On the other hand, from visual assessment followed that a higher degree of smoothing proved to be meaningful in noise reduction and removing invalid edges. So increasing the amount of smoothing leads to both positive and negative effects.

The conclusion of the virtual assessment is also accompanied by some contradic-tions, keeping in mind that calculations with σ = 1 led to the most significant results in the correlation tests between the average slice ridge/edge strength and rotation speed of the phantom. Thus, considering these results, one would say that σ = 1 would have to be used for best performance of the model. In that case, it would be assumed that the decline in edge/ridge strength at higher rotation speeds is caused by distortion of the edges of the Capitate and that this decline can not be monitored at higher smoothing, since smoothing itself also causes distortion of edges. Supported by implications that increased smoothing will lead to distortions in the edges rep-resenting the cortical layer, applying the lowest smoothing level would be advised. However, results following from these calculations may not be entirely representative for the edge/ridge strength of the edge representing the cortical layer, because of the presence of invalid edges and their contribution to the average edge/ridge strength sums.

At this point it is hard to say which of the two methods might be used best. When considering the correlation tests between the average slice edge/ridge strength and rotation speed, both edge- and ridge strength analysis found a significant correlation for σ = 1, while only in ridge strength analysis another significant correlation was found for σ = 2. According to these findings in combination with conclusions from visual assessment, ridge strength analysis seems to provide more significant results. Because of the suspicion of the presence of invalid edges, no conclusions can be drawn about which model is the better performing one yet. Therefore, more experimenting is advised on how to improve and expand both methods.

It is advised to perform edge/ridge strength analysis on more (carpal) bone data. By acquiring more edge/ridge strength data, quantification of motion blur through edge/ridge strength analysis will be more comprehensive and trustworthier, such that eventually could be determined at which edge/ridge strengths a scan can still be used for bone registration.

A possible favourable extension to this model would be automatic scale selection. Automatic scale selection has also been used in other imaging studies [18][19][20], in which the amount of smoothing was automatically changed for each pixel, based on local image structures. If the amount of smoothing can be regulated throughout the image, little smoothing can be applied at sharp edges that belong to the corti-cal bone layer, while more smoothing is then applied to the areas surrounding the

(16)

bone. In this way, little distortion will appear in the cortical layer of the bone, since less smoothing is applied there. Therefore, it will maintain its natural structure as much as possible. Besides that, invalid edge structures in the area surrounding the bone will be distorted. Thus, their contribution to the edge/ridge strength calcula-tions will be less- or not significant. Furthermore, some visualizacalcula-tions showed gaps or holes in the edge/ridge representations of the Capitate, meaning that the edge curves were not closed as they should be. Automatic scale selection would result in more closed curves, since the amount of smoothing would be able to vary along the edge as well [18]. In summary, isolation of the Capitate will improve, while maintaining the natural structure of the edge. The model could therefore provide more reliable and consistent results, using an automatic scale selection method.

6 Conclusion

The goal of this study was to quantify motion blur in the cortical layer of carpal bones in cropped 4D-CT wrist phantom data, by applying edge- and ridge strength analysis methods and using custom MATLAB code. Both methods were executed for multiple levels of smoothing and evaluated through statistical analysis and visual assessment. Trough visual assessment it was found that increased smoothing caused a decline in invalid edges found in the phantom. Since the goal was to focus on the analysis of the cortical layer, increased smoothing was advised. At increased smoothing an increase in distortion of the cortical layer was also observed, altering its natural structure. Simulation results also showed a sharp relative decline in edge/ridge strength around sharp edges as smoothing was increased, revealing that increased smoothing most likely affects results of analysis on the cortical layer. Statistical analysis was per-formed to check for correlation between edge/ridge strength and phantom motion. Results were significant for both methods when little smoothing was applied, indi-cating that low smoothing levels would have to be used for the best performance. This contradicted with conclusions that followed from visual assessment. Since ap-plying little smoothing indicated better performance in quantifying motion blur, this was eventually advised. However, edge/ridge strength results may not be completely representative for the edge/ridge strength of the cortical layer in this case, since at little smoothing invalid edges also had a contribution to the edge/ridge strengths. Recommendations are given in the Discussion.

7 Acknowledgements

Special thanks to ir. J.G.M. Oonk and Dr. ir. G.J. Streekstra for their guidance and support in this thesis process.

(17)

References

[1] Xin Chen, Jim Graham, Charles Hutchinson, and Lindsay Muir. Automatic inference and measurement of 3d carpal bone kinematics from single view fluo-roscopic sequences. IEEE transactions on medical imaging, 32(2):317–328, 2012. [2] DM Kushner, EM Braunstein, KA Buckwalter, K Krohn, and HA White. Carpal instability in rheumatoid arthritis. Canadian Association of Radiologists journal, 44(4):291–295, 1993.

[3] R Schmitt, S Froehner, G Coblenz, and G Christopoulos. Carpal instability. European radiology, 16(10):2161–2178, 2006.

[4] Yoon Seong Choi, Young Han Lee, Sungjun Kim, Hee Woo Cho, Ho-Taek Song, and Jin-Suck Suh. Four-dimensional real-time cine images of wrist joint kinemat-ics using dual source ct with minimal time increment scanning. Yonsei medical journal, 54(4):1026–1032, 2013.

[5] No´emie Debroux and Carole Le Guyader. A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and non-local shape descriptors. SIAM Journal on Imaging Sciences, 11(2):957–990, 2018.

[6] JB Antoine Maintz, Petra A van den Elsen, and Max A Viergever. Comparison of edge-based and ridge-based registration of ct and mr brain images. Medical image analysis, 1(2):151–161, 1996.

[7] Johannes GG Dobbe, Marieke GA de Roo, Jim C Visschers, Simon D Strackee, and Geert J Streekstra. Evaluation of a quantitative method for carpal motion analysis using clinical 3-d and 4-d ct protocols. IEEE transactions on medical imaging, 38(4):1048–1057, 2018.

[8] Mitsuhiro Nakamura, Yuichiro Narita, Akira Sawada, Kiyotomo Matsugi, Man-abu Nakata, Yukinori Matsuo, Takashi Mizowaki, and Masahiro Hiraoka. Im-pact of motion velocity on four-dimensional target volumes: a phantom study. Medical physics, 36(5):1610–1617, 2009.

[9] Julia F Barrett and Nicholas Keat. Artifacts in ct: recognition and avoidance. Radiographics, 24(6):1679–1691, 2004.

[10] Xuande Zhang, Xiangchu Feng, Weiwei Wang, and Wufeng Xue. Edge strength similarity for image quality assessment. IEEE Signal processing let-ters, 20(4):319–322, 2013.

[11] G Elisabeta Marai, David H Laidlaw, and Joseph J Crisco. Robust automated registration of wrist bones using tissue classification and distance fields. 2004. [12] GJM Tuijthof, L Beimers, R Jonges, ER Valstar, and L Blankevoort. Accuracy

of a ct-based bone contour registration method to measure relative bone motions in the hindfoot. Journal of biomechanics, 42(6):686–691, 2009.

[13] Lijkele Beimers, Gabri¨elle Josephine Maria Tuijthof, Leendert Blankevoort, Remmet Jonges, Mario Maas, and C Niek van Dijk. In-vivo range of motion of the subtalar joint using computed tomography. Journal of biomechanics, 41(7):1390–1397, 2008.

(18)

[14] Geert J Streekstra, Rein Van Den Boomgaard, and Arnold W. M. Smeulders. Scale dependency of image derivatives for feature measurement in curvilinear structures. International Journal of Computer Vision, 42(3):177–189, 2001. [15] Luc MJ Florack, Bart M ter Haar Romeny, Jan J Koenderink, and Max A

Viergever. Scale and the differential structure of images. Image and vision computing, 10(6):376–388, 1992.

[16] Alasdair McAndrew. An introduction to digital image processing with matlab notes for scm2511 image processing. School of Computer Science and Mathe-matics, Victoria University of Technology, 264(1):236, 2004.

[17] Carsten Steger. An unbiased detector of curvilinear structures. IEEE Transac-tions on pattern analysis and machine intelligence, 20(2):113–125, 1998.

[18] Tony Lindeberg. Edge detection and ridge detection with automatic scale selec-tion. International journal of computer vision, 30(2):117–156, 1998.

[19] Ersin Bayram, Christopher L Wyatt, and Yaorong Ge. Automatic scale selection for medical image segmentation. In Medical Imaging 2001: Image Processing, volume 4322, pages 1399–1410. International Society for Optics and Photonics, 2001.

[20] Dorin Comaniciu and Visvanathan Ramesh. Systems and methods for automatic scale selection in real-time imaging, April 11 2006. US Patent 7,027,643.

(19)

8 Appendix

8.1 Simulated edges

Figure 8: simulated smooth edge, with f (x, y) = −x2+ 255 for −√255<x<√255 and f (x, y) = 0 otherwise.

(20)

8.2 Directional derivative vectors

(21)

8.3 Average slice edge/ridge strength distributions

(22)

(23)

(24)

Figure 14: Average slice ridge strength for σ = 3

8.4 3D edge/ridge strength visualisations

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

8.5 Correlation tests: R output

Dependent variable: Edge strength (σ = 1)

rps −35,694.080∗∗∗

(9,688.206) Constant 64,556.280∗∗∗

(3,175.063)

Note: ∗p<0.1; ∗∗p<0.05; ∗∗∗p<0.01

Table 2: Correlation test between edge strength and the rotation speed of the phantom with σ = 1. A significant correlation was found (p <0.01).

Dependent variable: Edge Strength (σ = 2)

rps −2,923.997

(2,804.406) Constant 19,830.150∗∗∗

(880.288)

Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 3: Correlation test between edge strength and the rotation speed of the phantom with σ = 2. No significant correlation was found.

(36)

Dependent variable: Edge Strength (σ = 3) rps −39.062 (1,818.081) Constant 13,106.310∗∗∗ (587.705) Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 4: Correlation test between edge strength and the rotation speed of the phantom with σ = 2. No significant correlation was found.

Dependent variable: Ridge Strength (σ = 1)

rps −36,055.980∗∗∗

(8,101.544) Constant 62,483.240∗∗∗

(2,649.211)

Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 5: Correlation test between Ridge strength and the rotation speed of the phantom with σ = 1. A significant correlation was found (p <0.01).

(37)

Dependent variable: Ridge Strength (σ = 2) rps −2,274.625∗∗ (967.011) Constant 8,292.871∗∗∗ (301.271) Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 6: Correlation test between Ridge strength and the rotation speed of the phantom with σ=2. A significant correlation was found (p <0.05).

Dependent variable: Ridge Strength (σ = 3)

rps −570.241

(386.721) Constant 3,375.216∗∗∗

(120.126)

Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 7: Correlation test between Ridge strength and the rotation speed of the phantom with σ=3. No significant correlation was found.

(38)

Dependent variable: Edge strength Sigma −22,248.620∗∗∗ (1,399.617) Constant 74,482.100∗∗∗ (3,538.140) Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 8: Correlation test between the Edge strength and σ. A significant correlation was found (p <0.01)

Dependent variable: Ridge strength

Sigma −26,094.150∗∗∗

(1,497.696) Constant 74,376.850∗∗∗

(3,695.019)

Note: ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01

Table 9: Correlation test between the Ridge Strength and σ. a significant correlation was found (p <0.01)

(39)

8.6 Scatterplots Ridge/Edge strength against Sigma

Figure 58: Scatter plot of the Edge strength against σ

Quantifying motion artifacts in 4D-CT using edge strength analysis

Free University/ University of Amsterdam

Bachelor Thesis

Quantifying motion artifacts in 4D-CT using edge

strength analysis

N.J.J. Deen

Biomedical Engineering and Physics

Amsterdam UMC

Contents

1

Introduction

2

Theory

2.1

Edge strength

2.2

Directional Gaussian derivatives

3

Methods

4

Results

5

Discussion

6

Conclusion

7

Acknowledgements

References

8

Appendix

8.1

Simulated edges

8.2

Directional derivative vectors

8.3

Average slice edge/ridge strength distributions

8.4

3D edge/ridge strength visualisations

8.5

Correlation tests: R output

8.6

Scatterplots Ridge/Edge strength against Sigma