Using time information to make a robust feature extraction and parameter estimation from unipennate muscles in B-mode ultrasound sequences

(1)

Using time information to make a robust feature extraction and

parameter estimation from unipennate muscles in B-mode ultrasound

sequences

Nirusan Tharmanathan

Thesis submitted for the degree of

Master in Digital Signal Processing and Image Analysis 60 credits

Department of Informatics

Faculty of mathematics and natural sciences

UNIVERSITY OF OSLO

(2)

(3)

Using time information to make a robust feature extraction and

parameter estimation from unipennate muscles in B-mode

ultrasound sequences

Nirusan Tharmanathan

Using time information to make a robust feature extraction and parameter estimation from unipennate muscles in B-mode ultrasound sequences

http://www.duo.uio.no/

Printed: Reprosentralen, University of Oslo

(4)

(5)

Abstract

Ultrasound imaging for medical applications is considered a cost-effective and non- invasive method of studying the internal structures of the body. One application is to study muscles and muscle architecture. Muscle architecture refers to a muscle and the arrangement of structures essential to the muscle’s function. One aspect of studying muscle architecture is to study how fascicles move and the fascicles length and pennation angle relative to an aponeurosis. These measurements are used to give us information on force production. In this thesis, the aim was to create an algorithm for B-mode ultrasound sequences that could detect the shallow and the deep aponeurosis and create a robust reference fascicle that could be used to measure the changes in pennation angle and fascicle length over time. We did not manage to create a robust method for detecting the deep aponeurosis, but we were able to find and detect the shallow aponeurosis in each frame consistently. We also developed a method of using à priori knowledge to find the location of the aponeuroses faster than previously. The construction of the reference fascicle as proposed in this thesis seems to be a robust way of creating a reference fascicle especially in sequences were some of the frames have little relevant information.

(6)

(7)

Acknowledgements

This thesis has been done in collaboration with the Department of Informatics at the University of Oslo. The thesis started in August 2018 and finished in May 2019.

I would like to thank my supervisor Professor Anne H. Solberg. I want to especially thank her for her patience, understanding, knowledge, and advice throughout the whole process.

I also would like to thank my family, for all of their support, encouragement, and love throughout my life and especially now this past year. And I want to thank Prableen for always believing in me and supporting me.

(8)

(9)

Introduction

Ultrasound imaging for medical applications is considered a cost-effective and non- invasive method of studying the internal structures of the body. One application is to study muscles and muscle architecture. Muscle architecture refers to a muscle and the arrangement of structures essential to the muscle’s function. One aspect of studying muscle architecture is to study how fascicles move and their length and pennation angle relative to an aponeurosis. These measurements are used to give us an indication of force production. Being able to measure force production can be used in several applications where one want to follow a subject over time or as an extended method of assessing a subject.

Traditionally, one has to do measurements by hand on one image at a time. Manual measurements have a higher chance of inconsistent results, especially from one frame to another or between people performing measurements on the same image.

Constructing reference fascicles by hand is also time-consuming. There exist automated methods for measuring parameters for one image at a time, but there is a need for applying and robustifying available methods on sequences where we follow some movement over time.

The scope of this thesis This thesis aims to develop an automated algorithm that can estimate the pennation angle and fascicles for ultrasound sequences. The aim is also to use the information provided from one or more previous frames to make a more robust tracking of the shallow and deep aponeurosis. We also want to find a robust way of representing the reference fascicle. Achieving these aims is going to provide a method that eliminates user interaction, potentially providing a more objective and consistent method of estimation.

The structure of this thesis In chapter 2 we are going to provide an introduction to the different topics making up this thesis and discuss previous research. In chapter 3 we are going to look into the data material provided. In chapter 4 we will discuss how we enhance the features so that we better extract relevant information for this thesis. We will then move on to talk about the Radon-transform in chapter 5; this will end the more theoretical part of the thesis. We will then in chapter 6 use the theory to find and track aponeuroses and propose a method to construct a reference fascicle in

(12)

chapter 7. We will then end the thesis by evaluating the proposed method in chapter 8 and discussing the results, before concluding and discuss further work in chapter 9.

(13)

Chapter 2

Background

In this chapter, we want to introduce the different topics that make up this thesis. We also provide a brief introduction to ultrasound imaging and muscle architecture. We will also look into how to describe an ultrasound image of the muscle architecture.

Finally, we look into the previous work relevant to this thesis.

2.1 Ultrasound and ultrasound imaging

Ultrasound is any soundwave above 20 kHz[1, p.6]. Soundwaves are compressional waves that propagate through a medium, like air or water in a longitudinal direction.

The compressional wave can be periodic or pulsed, and through water or biological tissue, the compressional wave propagates around 1500 m/s [1, p.4].

2.1.1 Generating and receiving ultrasound

Ultrasound is generated using a transducer/ultrasound-probe. To make a transducer one typically uses piezo-electric crystals. Applying a voltage to the crystal will produce a pressure wave. The piezo-electric crystals also have the property that pressure applied to the crystals is going to generate a voltage. This dual trait gives the transducers the ability to generate pressure-waves as well as to register incoming pressure/sound waves as electric current.

When using a probe on a person, the probe is placed against the skin, applying a voltage to the crystals. Applying a voltage to the probe generates the ultrasoundwaves that propagates through the skin into the tissue. As the ultrasoundwaves propagates through the body, the waves are either going to be reflected, absorbed or scattered by the tissue. Different types of tissue are going reflect, scatter or absorb the oncoming acoustic wave differently; this is due to differences in acoustic impedance ([1], [2]) between two different types of tissues/mediums. In an ultrasound-imaging we are interested in the pressure waves reflected back to the skin and then picked up by the transducer. The waves that return due to reflection are then converted to an appropriate voltage level based on the pressure exerted on the crystals. The intensity/voltage

(14)

level and the time it takes for a sound wave to get back to the probe is what is used to make a B-mode ultrasound image. A strong response indicates that there is a highly reflective object while a weaker response indicates that something is less reflective. At the same time, the time the sound wave to get back also indicates the distance from the probe. The pressure waves also get attenuated as they move farther into the tissue; this attenuation means that waves that get reflected deeper in the tissue does not return to the probe. By using all this information, we can generate an image that shows the different area of the body where the audio waves have been reflected. How much of the soundwave that has been reflected is going to decide the contrast of the image relative to the surrounding area [1, p.8].

2.1.2 B-mode ultrasound imaging

We can split ultrasound imaging into several modes. The mode we are interested in this thesis is B-mode [2, p.374]. B-mode stands for Brightness-mode, and it is a brightness modulated image where we observe the cross-section of the tissue in different depths.

Figure 2.1 shows an example of an B-mode ultrasound image. Using reflected pressure waves as a way of generating images has some challenges associated with it, the most common ones are resolution, framerate, and speckle.

Resolution

Resolution is dependent on the frequency of the emitted pulse. Higher frequency usually means better resolution. The challenge of using acoustic waves with a high frequency is that the waves are going to get attenuated more as they travel farther into the body. The way to remedy this is to use a pulsed wave with a lower frequency; this increases the range we can see into the body. However, using a wave with a lower frequency is going to affect the resolution of the image since lowering the frequency is going to affect how much of the wave is going to interact with the tissue. More delicate details and smaller objects in the image is not going to be visible.

Usually, it is easier to image tissue closer to the surface since one can use waves with higher a higher frequency without it getting "lost" in the tissue. Resolution is also dependent on framerate.

Frame rate

The frame rate (how often a new image is taken) is also an important consideration when dealing with ultrasound images. As discussed earlier an ultrasound image is generated by acoustic waves traveling through the body and hitting tissues that make the waves reflect back to the surface. This process takes time since the wave is propagating at a finite speed and has to travel from the probe, hit some tissue that is going to reflect it and get back. The travel time of sound-waves have implications when it comes to observing movement in tissue, the implication being reducing our ability to accurately and fluently observe movement.

(15)

Figure 2.1: A B-mode ultrasound image of the muscle architecture. We see the arrangement of different structures

Speckle

Speckle is a type of artifact in ultrasound image that comes as a result of us using waves and its reflections to image something through a medium. One often observes speckle as white dots throughout the ultrasound image. In ultrasound imaging speckle is caused by constructive and destructive interference in the tissue. The reason for this interference is that small irregularities on the tissue surface scatter the soundwave [2, p.278].

2.2 Muscle architecture and movement

For this thesis, we are interested in muscle architecture[3]. More specifically, we are studying how the muscle fibers act during flexion. As previously mentioned studying muscle fibers during movement can give us some insight when it comes to how the body generates force. Before we can do this, we have to establish the arrangement of the muscles in the human body; often referred to as muscle architecture.

(16)

Figure 2.2: Illustration of how the muscles are constructed [5].

2.2.1 Muscle architecture

The arrangement of skeletal muscle fibers also referred to as the skeletal muscle architecture is the primary determinant of muscle function when it comes to movement.

Understanding muscle architecture can give us insights when it comes to the physiological basis of force production.

The building blocks of a muscle

To able to talk about muscles and muscle architecture we have to make some definitions.

These definitions are not complete and are only intended to help our understanding when it comes to talking about and noticing features and structures we can see in ultrasound images.

Muscle mostly consists of proteins, the two most abundant ones being actin and myosin. The actin and myosin are assembled making up a sarcomere. When stacking these sarcomeres, one after the another we create a myofibrils also known as muscle fibril/muscle fibers. A bundle of muscle fibers makes up afascicle. The collection of these fascicles are what makes up what we call a skeletal muscle [4]. All this is visualized in Figure 2.2.

(17)

Figure 2.3: An illustration from [6] showing different orientations.

Orientation

How the skeletal muscles move are based on their arrangement and their angle relative to the axis of force generation. The angle between the muscle fiber and the force generating axis for mammalian muscles can vary from about 0%−30%, and most muscles fall into the category of being multipennate muscle, meaning that the fibers in the muscle are oriented in several angles relative to the axis of force generation. There are two other categories when it comes to muscle architecture, these being longitudinal and unipennate. Longitudinal muscles have an orientation parallel to the axis of force generation. Unipennate muscle’s orientation is at a fixed angle relative to the force generating axis[6]. Muscles that have a unipennate architecture is going to the focus in this thesis. An illustration of the different orientations is illustrated in Figure 2.3

2.2.2 Using B-mode ultrasound imaging to study muscle arcithecture We now have some more knowledge and tools to study and make sense of what we see in Figure 2.1. Figure 2.4 is the same as Figure 2.1, but some of the more distinctive parts of the image are colored.

Starting at the top of the image in Figure 2.4 we see a highly reflective area that we have colored in green this is the skin is and where the probe is placed. The further we look downward the deeper we are looking into the body. As we look further downwards we first meet the blue line, this is the upper/shallow-aponeurosis. After this we see some white lines (some of them are highlighted in yellow) going downwards at an angle, these lines are the fascicles. As the fascicles go downwards, it eventually meets

(18)

Figure 2.4: A B-mode ultrasound image where we have highlighted different areas that are of interest. The skin(green), the upper/shallow-aponeurosis(blue), fascicles(yellow) and finally the lower/deep aponeurosis(red).

(19)

the red line. This is the lower/deep-aponeurosis. During flexion, it is possible to see the fascicles move. Since we can see how the fascicles are moving, it might be possible to automatically track the movements and measure the changes in fascicle lengths.

2.3 Previous work

In [7] the authors proposed a way of estimating the fascicle orientation in B-mode ultrasound images. The orientation estimates were done on the fascicles located in the vastus lateralis muscle. The authors first enhanced the fascicles by using a vessel enhancing filter [8], in section 7.2.2 we discuss this method in detail. The authors then performed a Radon transform and a wavelet analysis on the filtered B-mode ultrasound images. How to perform a Radon transform and how to extract the dominant orientation will be discussed in chapter 5. To find the dominant orientation using wavelet analysis the images the authors manually cropped the image so that it only spanned from the superficial to the deep aponeurosis. The acquired orientation was compared to manually acquired measurements. The manually acquired measurements were done by measuring fascicles spanning from the shallow to the deep aponeurosis.

The given orientation is the orientation relative to the x-axis. To test the two methods the authors used synthetic and sixty digitized ultrasound images. The synthetic images are grids of parallel lines with a known orientation that had sinusoidal changes in intensity across the lines. There was also random noise added to the synthetic images.

The real images are taken from the distal part of the left vastus lateralis while the subject was cycling a stationary ergometer. The sixty ultrasound images from this sequence were then manually digitized by ten different researchers. Each researcher digitized the sequence twice. The authors concluded that the Radon-transform gave orientation estimates that most closely matched the manually measured orientations on actual ultrasound images. It is worth noting that this paper extracts the dominant orientation in the whole region determined by the user. This method will therefore not consider local variations in orientation. Local orientation is essential to consider if one wants to estimated the orientation of the fascicles close to the aponeurosis as these tend to curve more than the rest of the fascicle.

In [9], the authors devolped a tracking method based on the Lucas-Kanade optical flow algorithm [10]. The algorithm was also extended using affine optical flow. Developing this tracking method makes it possible to track how the fascicle length changes during movement. The authors collected images of the fascicles from the medial gastrocnemius muscle, with a sampling rate of 80Hz. The authors defined the length of muscle fascicle as the straight line distance between the superficial and the deep aponeurosis. The first frame of every sequence was repositioned so that a fascicle was visible from the the superficial to the deep aponeurosis. After positioning the fascicle in the middle, the authors started the tracking. The results of the automated tracking were then compared with manual measurements done on the same sequence.

The results where deemed promising with a CMC¹of 0.90±0.09 compared with the manual measurements. As mentioned by the authors, some assumptions had to be

1The coeffienct of multiple correlation(CMC) produces a scalar R that serves as a measure of how strong the association between independent variables and one dependent variable. The R scalar can be any value from 0 to 1. R = 1 being a strong linear association, and R = 0 being no linear association.

(20)

made to make the method more robust. The first assumption being that the muscles during movement are only going to change in a way that an affine transformation which considers translation, rotation, dilation, shear, and skew is going to find the same area in the next frame. The method also assumes that the fascicles are straight between the shallow and the deep aponeurosis.

The authors in [11] proposed an automatic method of identifying the pennation angle and fascicle length of the gastrocnemius muscle. The method was also able to identify the superficial and the deep aponeurosis automatically and to estimate the length of a fascicle situated between the superficial and the deep aponeurosis. The authors performed testing on both synthetic and real data. The authors wanted to enhance the line-structures associated with the fascicle in the ultrasound images. The enhancement was done using a set of wavelets. After filtering the images, they were Hough-transformed. The Hough transformed images were used by the authors to find the dominant orientation of a fascicle. The orientation was found by finding the area in the Hough transformed image that had the highest variance. To find the dominant orientation of the aponeurosis the authors used the normalized Radon transform to transform the image. As with the Hough transform, the dominant orientation was found by finding the orientation that had the highest variance. Now having the measurements, the authors calculated the fascicle length. To verify the results, the same sequences were also manually measured and compared using Pearson product-moment correlations². The correlation between the proposed method and the manually measured results whereby the authors described as highly correlated with a PCC ofr =0.89±0.05. In situations where the previous method in [9] failed to measure the whole fascicles length the proposed managed to estimate more of the fascicles length. The automatic method did fail to estimate some sequences as the method lost its targets in some frames. Losing the target was due to high speckle noise and inhomogenous deformation across the if interest between consecutive images.

This error could not be corrected making it necessary to examine the fascicle length estimation to see if the results were valid. As with the other papers, the authors assumed that the fascicles are straight between the shallow and the deep aponeurosis.

In the thesis [12] the author worked on B-mode ultrasound images of the vastus lateralis muscle. The author’s main aim was to enhance the structures in the ultrasound images and reduce noise. The author also provided a method of estimating the pennation angle and fascicle length in a single frame. The method also includes a way of detecting the region of interest, this being the region between the shallow and the deep aponeurosis. To enhance the structures in the images, the author used Knutsson tensor filters. On the structure enhanced images Radon-transform was used to find the shallow and the deep aponeurosis. The Radon-transform was also used to find the angle of the fascicles so that one could construct a reference-fascicle.

This reference fascicle is a fascicle that is generated based on the measurements done in the region of interest considering the orientation of fascicles in a specific depth.

This method produces a single fascicle that depicts an "average" fascicle is going to behave. The thesis then concludes with a method of calculating the fascicle length and the pennation angle. The data set used in the thesis is the same as used in as the

2The Pearson product-moment correlations (PCC) is a measure of the linear correlation between two variables. The PCC produces a value r that ranges from -1 to 1. r = 1 is total linear correlation, r = 0 is no linear correlation and -1 is total negative correlation.

(21)

two first training sequences in this thesis (chapter 3). The author did not have any manual measurements to compare the results with so a visual inspection was done instead. When detecting the superficial and the deep aponeuroses the method worked well when the aponeuroses were straight and clearly defined. When the aponeuroses had a substantial curve to it or were not visible, the method struggled to find it and highlight it appropriately. The fascicle measurements had similar results. The author also experienced that measurements in the edges were not so reliable as the other measurements. The algorithm also had difficulties when the aponeuroses had a shape other than a straight line or a concave curve. When measuring the pennation angle the results were dependent on how accurate the locating of shallow aponeuroses was, this naturally also affected the measurement of the fascicle length.

In the thesis [13], the author proposed an automatic method that identifies the shallow and deep aponeuroses and marks them. Between the shallow and deep aponeurosis where the fascicles are visible, the algorithm places points are randomly distributed points. Around these points, a normalized Radon-transform is performed, and by finding the area on the radon domain with the most variance and determining that area as the dominant orientation. Thes measurements were then used to construct a reference fascicles that was used to extract a potential fascicle length and the and pennation angle. The data set used in the thesis is the same as used in this thesis (chapter 3). The author improved the fascicle measurements and the aponeurosis detection proposed in [12]. The algorithm sometimes struggles to locate the aponeuroses through the sequence but when it fits the results are good. The results when measuring fascicles were mostly good, with pennation angles and manual measurement being somewhat close. This thesis is most closely related to the work done in this author, so through this thesis, we will refer and explain the approach suggested in [13] closer.

2.4 Summary

We have now looked into the different topics making up this thesis. We have also been looking into previous work. When it comes to measuring the fascicle length and pennation angles, there are some different approaches. To start every implementation has a pre-processing stage to enhance the structures in the ultrasound images. From the literature, we see that there are three primary methods of enhancing fascicle-structure.

The first pre-processing method is to use wavelet as in [7], the second one is to us vessel enhancement method proposed by Frangi et al. as done in [7] and [13] . The final method is to use Knutson-tensor filters as done in [12]. In this thesis, we have decided to use the vessel enhancement technique proposed by Frangi et al. [8]. The reason for this is that this is the method most widely-used in different situations, and based on the results it seems to be able to enhance the structures that we are interested in so that the estimates are within an acceptable range. When it comes to observing how the length of fascicles changes during flexion the Frangi-filter approach makes a better alternative. There is generally a difference between the approaches when it comes to how the data is acquired. For real-time-ultrasound sequences where the data/images are more or less directly received from the ultrasound scanner, using an optical flow algorithm in conjunction with the Radon-transform to estimate the

(22)

orientation is preferred. For non-realtime settings, there is a tendency to mostly rely on the Radon-transform to constructing a reference-fascicles that represent the underlying data to measure changes in the fascicle length changes. This is done by using the relationship between the shallow and deep aponeuroses to extrapolate how to the fascicles length changes in conjunction with the apparent pennation angle. For this thesis, we determined that the best approach is to continue/ base this thesis of the work done in [13] and [12] as the datasets are more or less the same and the results are promising.

(23)

Chapter 3

Data material

In this thesis, we had three ultrasound sequences available. The Norwegian School of sports science provided the ultrasound sequences. The sequences were collected using a Philips HD11 XE scanner operating at 12Mhz, resulting in a sequencs with a framerate of 43Hz. The sequences show how fascicles located in the vastus lateralis moves during movement. Two of the sequences are used to train the algorithm while the last sequence is used as a test sequence to asses how robust the proposed implementation is.

3.1 Data format

The captured sequences are handled according to the Digital Imaging and Communica- tions in Medicine (DICOM) protocol [14]. The DICOM-protocol is an all-encompassing protocol that handles data transfer, storage and how to display the images. The data transfers happen from any modality, as in this thesis an ultrasound scanner to a storage device. In the storage device, the data as DICOM-formated file, after archiving the files they can be accessed by different workstations. In this thesis, we are working with images saved as tagged image file format (tif) images.

3.2 The sequences

Generally, we observe that the visible fascicle’s movement varies through the sequence.

In some parts of the sequence, the fascicles are stationary, while in other parts of the sequence the fascicles move slowly. We also have a final case where the movement is so fast that the fascicles "jumps" from one frame to another. We also observe that different parts of the fascicles in the image move at different speeds due to flexion.

In some parts of the sequence, the fascicles are visible while in some cases it is not the case. During fast movements, the fascicles also tend to "smear" as well to "jump"

between frames. A final thing to notice is that during the sequence, the fascicles can disappear from the frame due to movement in space. Either by moving out of frame or disappearing by changing the depth of its location. The disappearing fascicles can

(24)

be a result of things happening during acquisition. During acquisition movement by either the subject or the operator can change the ultrasound probe’s position relative to the part of the muscle architecture we want to study.

3.2.1 Training sequence 1

Training sequence 1 consists of 391 images; two images from this is sequence are shown in Figure 3.1. This sequence has mostly visible fascicles throughout the whole sequence except in the top right corner. The shallow aponeurosis is always visible and mostly have an orientation of 0^◦. The deep aponeurosis this sequence is less visible.

The deep aponeurosis has a slight curve to it, and during flexion, it moves towards the shallow aponeurosis and then down again to its starting position.

Training sequence 2 consists of 226 images and is shown in Figure 3.2. The first thing that one might notice is that three distinct horizontal lines are the aponeurosis In this thesis we are most interested in the region closest to the probe, so we are interested in the top-most bright aponeurosis and the aponeurosis below that. Here the shallow aponeurosis is the less bright aponeurosis in the sequence, while the shallow aponeurosis is the most bright aponeurosis one throughout the sequence. The fascicles are in some parts of the sequence not shown as a continuous line, but rather as discontinuous line structures curving in a general direction. We also observe more of a smearing effect where regions only consist of a bright texture.

3.2.3 Training sequence 3 and the test sequence

The validation sequence consists of 676 images and is by far the most complex sequence when it comes to testing the robustness of the method. This sequence was provided late in the process of developing this algorithm; thus it seemed fitting that we could use this sequence for testing. We start by looking at the different aponeurosis in the image. We see that both the shallow and the deep aponeurosis will not necessarily have an orientation around 0^◦. The two brightest objects in the ultrasound image are the shallow-aponeurosis and the deep aponeurosis closest to the skin. We also see that throughout the sequences we have some bright lines in the bottom part of the ultrasound sequence. The fascicles throughout the sequences are mostly only partially visible; this makes them more similar to the more difficult sections in training sequence 2. Since this sequence is so different from the two other training sequences, we decided to use the first 200 to make a third training sequence. The last 476 images in the sequence are only used to test the algorithm. We can see an image from this sequence in Figure 3.3.

(25)

(a) Here we see a frame where the fascicle are clearly defined. The shallow aponeurosis is also clearly defined, the deep aponeurosis is less so.

(b) In this frame the fascicles are less defined, but the shallow aponeurosis is clearly visible.

Figure 3.1: Here we see two images from training sequence 1, showing the variation of visible features within one sequence.

(26)

(a) Here we see an image where large parts of the fascicles are visible. We also see a clearly defined deep aponeurosis. The shallow aponeurosis is less bright.

(b) Here we see a more challenging part of the sequence. We observe some "smearing" due to fast movement.

Figure 3.2: Here we see two images from training sequence 2, showing the variation of visible features within one sequence.

(27)

Figure 3.3: A frame from the training set part of the third sequence. Here we see that the fascicles are not that clearly visible, and that the deep aponeurosis is at an angle.

3.3 Summary

In this chapter, we have looked at the different sequences we will be using in this thesis. The first two sequences will be frequently used to modify and develop the algorithm throughout the whole process. We want to avoid using the final training sequence as much as possible, before testing our algorithm on the last 400 frames from the third sequence.

(28)

(29)

Chapter 4

Feature enhancement

We want to enhance the line structures associated with the fascicles in our ultrasound- sequence. Enhancing the fascicles is going to aid us when estimating the orientation of the fascicles. One thing that makes it more challenging to enhance lines in ultrasound- images is speckle noise. During an ultrasound sequence, the gray-scale intensities in an area are going to vary. The variation of gray-scale intensity is problematic since finding lines and edges in an image is traditionally done by finding the gradients of the image, where an area of large change/gradient is an indication of an edge. Using simple gradient-filters for high-pass filtering as Sobel is not going to enhance the structures sufficiently as discussed in [12].

In the previous work section (section 2.3) we saw that there were mainly three types of filters used to enhance the structures in an image. The first one is made using wavelets [11], the second one being Knutson-tensor-filters used in [12], and finally multiscale vessel enhancement filtering in as used in [7] and [13]. In this thesis, we decided to use multiscale vessel enhancement filtering [8]. This thesis is closely related to the work done in [13] we, therefore, chose that enhancement technique. The method is also more widely used on different types of images making it potentially more robust.

In this chapter, we want to give an introduction to the Frangi-filter and how it enhances line structures.

4.1 A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A.

Viergever, “Multiscale vessel enhancement filtering”

The method is described in [8]. Originally the method was developed for enhancing blood-vessels in magnetic resonance angiography images. This method of enhancing tubular structures has then since also been used on ultrasound-images([7], [12] and [13]). The transition from MR to ultrasound has been possible since blood vessels and fascicles share some of the same features as they are long-tubular structures.

The proposed algorithm by Frangi et. al uses Gaussian second derivative filtershxx,

(30)

h_yyandh_xy:

h_xx(x,y) = ¹ 2πσ⁴

x² σ²−1exp

−(x²+y²) 2σ²

(4.1) hyy(x,y) = ¹

2πσ⁴ y² σ²−1exp

−(x²+_y²) 2σ²

(4.2) h_xy(x,y) = ^xy

2πσ⁶exp

−(x²+y²) 2σ²

(4.3) Where σ is the size of standard deviation for each kernel. By using the second derivative filters, we find the true edge and not the region in the image with the steepest gradient as one would find using a first derivative filter. The filters are then convolved with an image Ito find the second-order partial derivatives of the image (I_xx,I_yyandI_xy):

Ixx= ^∂

2I

∂x² = I∗hxx (4.4)

I_xy = ^∂

2I

∂x∂y = I∗h_xy (4.5)

I_yy = ^∂

2I

∂y² = I∗h_yy (4.6)

The three second-order partial derivatives make the necessary elements for constructing a Hessian matrix. Having the Hessian matrix of the current image makes it possible to extract the principal-directions from the eigenvalues associated with the Hessian matrix.

H_i,j =^hÎ^xx⁽î,j⁾Î^xy⁽î,j⁾

Ixy(i,j)Iyy(i,j)

i

(4.7) When looking for tubular structures the relationship between the first and second eigenvalue(λ₁,λ₂) are assumed to beλ₂≈0 whileλ₁is larger thanλ₂.

After finding a tubular structure with the eigenvalues that have the relationship described, we useλ₁andλ₂to define the two relationshipsRandS.Ris a measure of if there are any tubular structures present. The largerR, the more "tubular" the area is.

Sis the norm of the Hessian matrix; a high S indicates that the area has high contrast.

Ris given as:

R= ^λ¹

λ₂ (4.8)

whileSis given as:

S= q

λ²₁+λ²₂ (4.9)

Based on these relationships Frangi et. al proposed a "vesselness" measure:

V₀=

( 0 i fλ₂ >0, exp

−^R²

2β²

(1−exp

−^S²

2c²

)otherwise (4.10)

(31)

whereβandcare threshold values that control the sensitivity of the measuresRandS.

The valuesβandcare chosen by the user as different types of images require different thresholds.

The final step in the algorithm is to determine the filter that gives the best "vesselness"- response for a range ofσs and the given thresholdsβandc.

V₀=max

σ V₀(σ) (4.11)

The result should be a filtered image that has the best "vesselness" given the two thresholds and the range ofσ-values.

4.1.1 Looking at a filtered image

When using the algorithm proposed by Frangi et al., three parameters have to be set. The setting of these parameters has been done differently from author to author implementing the algorithm. We ended up using the parameters proposed in [7] and in [13] since the parameters seemed to serve them well in both applications.

In [7] and [13] σwas set to range from 1 to 3 with an interval of 0.5. The threshold values were set toβ=0.5 andc=0.5. In Figure 4.2 we can see the results of filtering the image in Figure 4.1.

After filtering the image, we see that the lines associated with the fascicles and aponeuroses are enhanced. In some areas, especially in the region below the deep aponeurosis, we can see some artifacts in the form of "wavy" lines. These are not related to the fascicles or aponeuroses. These artifacts are therefore important to consider when trying to observe lines in the image.

4.2 Summary

In this chapter been looking at how to pre-process the image to enhance the structures we want to measure. We chose the vessel enhancement technique proposed by Frangi et al. [8]. For the rest of the thesis we setσto range from 1 to 3 with an interval of 0.5, and the threshold values are set toβ=0.5 andc=0.5.

(32)

Figure 4.1: An arbitrarily selected image from a B-mode ultrasound sequence

(33)

Figure 4.2: The results of using the parameters proposed by [7] and [13]

(34)

(35)

Chapter 5

The Radon Transform applied on 2D grayscale images

Formulated initially to study manifolds the Radon-transform [15] can also be used to study line segments and their orientation in a grayscale-image [7]. Studying the orientation of the images is done by transforming the image via the Radon-transform to a space referred to as the Radon-space. The Radon transform integrates over lines in the image, and this forms an intensity image where the axes now areρandθ, i.e., the polar coordinates. In this thesis, we are only going to look at how to transform discrete 2D grayscale images.

In this chapter, we will also look at how to use the Radon transform to estimate the dominant orientation in the image. Since there are different methods of doing this and some things one has to consider we briefly discuss this as well.

5.1 The Radon transform

The Radon-transform makes use of the relationship that a point in a Cartesian coordinate(x,y) can be described as a polar coordinate (ρ,θ), and thus makes it possible to describe a line asxcos(θ) +ysin(θ) =ρ.

From [16, p.390] we see how to transform an discrete gray-scale image to the Radon- space. We can look at an image as a 2D-matrix f with the sizeM×N, where every value/pixel in f is found by the points x and y. To transform an image f to the Radon-spaceg:

g(_ρ,θ) =

M−1 x

∑

=0

N−1 y

∑

=0

f(x,y)δ(xcos(θ) +ysin(θ)−ρ) _(5.1) Whereθis the angular-coordinate/the orientation of an accumulated response, andρ is the perpendicular distance to the accumulated response from the origin. Since we are performing a discrete Radon transform, theδrefers to the Kronecker delta which is defined asδ(x) =1, whenx=0.

(36)

Figure 5.1: This figure illustrates how the Radon transform is used to transform an image. The illustration is from [18]

We can look at the Radon transform of an image as the relationship between the number of accumulated responses for a giveng(ρ,θ), where an impulse response is made each time there is a transformation of an apparent line. In Figure 5.1 we see how the Radon transform projects and accumulates values asθandρchange relative to the Cartesian coordinates. It also worth noting that the origin of the Radon-transform is set to be in the center of the image(f(^M/2,^N/2)). The Radon transform accumulates the image intensities along the projected line from x⁰ toρ. It is this projected line’s distance from the origin we are referring to when talking aboutρ. When transforming the image, the orientationθ is set to span from 0^◦ to 180^◦. Defining the span ofθas such is possible since the radial line accumulates values in both directions at the same time. We have opted to use MATLAB’s implementation of the Radon-transform [17].

The reason for this was done to save time both from an implementation perspective, but also from a runtime and accuracy perspective.

5.1.1 Interpreting the orientation

From how we formulate and use Radon-transform, there are two different ways of describing the orientation of lines. The first way of describing the orientation is when the origin is in the center of the image as in Figure 5.1. Having the origin in the center means that when we are observing horizontal lines, they are going to be shown to have an orientation of 90^◦ in the Radon-domain. The second way of describing orientation is by placing the origin in the top-left corner of the image. In this thesis, we will mostly use the latter way of describing the orientation when we later construct the reference fascicle and compare our results with the fascicles in the image. This way of describing the orientation might make it more intuitive. For the same horizontal line, the orientation will now be 0^◦. Throughout the thesis, we will explicitly mention which orientation we are using.

(37)

5.2 Normalized Radon transform

We want to make the Radon-transform more robust. We have therefore decided to use the normalized Radon transform. The reason for this is two-fold, firstly the images consist of lines segments with different lengths and intensities. We are mostly interested in the orientation of the lines, so by normalizing we are focusing more on that aspect. Secondly, we want to use and compare the Radon-transforms performed in different areas within a region of interest. The only way of comparing different Radon-transform are if they are normalized similarly. We are later going to compare Radon-transforms that happens at different points in time, and this also requires a normalized Radon-transform.

We normalize the Radon-transform as such:

g(ρ,θ) = ^∑

M−1

x=0 ∑^Ny=⁻0¹ f(x,y)δ(xcosθ+ysinθ−ρ)

∑^Mx=⁻0¹∑y^N=⁻0¹δ(xcosθ+ysinθ−ρ) ^(5.2) The Radon transform is normalized to a Radon-transformed image of size(M,N) where every point/pixel in the image is a delta response.

5.3 Finding the dominant orientation in the Radon domain

By finding the dominant orientation in the Radon space, we are finding the most common orientation for all the line segments in the actual image/region as illustrated in Figure 5.2. In [7] found the dominant orientation by finding theθ-value that had the highest variance or the kurtosis, corresponded with the dominant orientation. For synthetic data, Rana et al. concluded that the variance seemed to give the mos accurate results while on actual ultrasound images finding theθwith the highest kurtosis was the most accurate.

In [12] and in [13] theθ-value that had the highest variance was used to determine the dominant orientation.

5.3.1 An example of estimating the dominant orientation

We want to verify that we can find the dominant orientation within a specific region on the image. We start with estimating the orientation of an aponeurosis as seen in Figure 5.3. The reason for starting with the aponeurosis in Figure 5.3 is that we visually can observe that should be around 90^◦or an orientation around 0^◦ if we imagine that the origin is in the top-left corner of the image. To enhance the image we also performed a Frangi filtering with the values from subsection 7.2.2. We can see the results before and after the filtering in Figure 5.4.

Transforming the region of interest

We decided to specify theθto span from 0 to 180 with a spacing of 0.1. Since we are only interested in how the Radon-transform transforms line segments to the Radon

(38)

Figure 5.2: An illustration from [7]. This figure shows how one can use the variance to find the dominant orientation in a Radon transformed image.

Figure 5.3: Image of a B-mode ultrasound image containing an aponeurosis we want to estimate the orientation.

(39)

(a) The region of interest before Frangi-filtration

(b) The region of interest after one performs vessel enhancement as proposed in [8]. We have altered the contrast for better visualization.

Figure 5.4: Region of interest before and after vessel enhancement

Figure 5.5: The Radon transform of Figure 5.4b

domain, we decided to use the regular Radon-transform and not the normalized radon transform.

Figure 5.5 shows how transforming the line segment shown in Figure 5.4b is going to look in the Radon-domain. When looking at the figure, we can see that the most accumulated values occur when the distanceρis around 0 from the image center. The θ-value that looks to have the most accumulated values are around 90^◦±10^◦.

In Figure 5.6 we can see the orientation estimate done by using the variance and the kurtosis. Using the variance to estimate the orientation we get an estimated orientation of 88.8^◦while using the kurtosis gives us an estimate of 88.6^◦. If we place the origin at the top-left corner the orientation is−1.2^◦and−1.4^◦. The estimated orientations are close to what we expect to get from a mostly horizontal line. The discrepancy is most likely due to the curving we see in the rightmost part in the image in Figure 5.4b and the other artifacts in the upper and bottom half of the image.

(40)

Figure 5.6: Here we see where theθvalue with the most variance and theθvalue with the highest kurtosis are within the Radon space.

5.3.2 Finding the dominant orientation in an area where a fascicle is present

We want to see if there is a difference when using the variance versus the kurtosis when we want to estimate the dominant orientation in an area where a fascicle is present if there is a difference, how large the difference is.

Before we can do that we have to consider how large the area around a fascicle should be to be able to estimate its orientation adequately. Choosing a too large window is going to make it sensitive to other structures in the image while making the window too small is going to make it the estimate more sensitive to noise.

In [12] the size of the local area/window was 70×70 and in [13] the size of the local area/window was 30×70. We will also be using a window size of 30×30 as a third window to test if that is going to contribute to a different result. In the test, we selected a single point in an ultrasound-image Figure 5.7 and then performed a Radon-transform on that point and its surrounding area. Before performing the Radon transform, we use the Frangi-filter, and we varied the size of the surrounding area we wanted to transform.

Radon transforming a point and its surrounding area

Figure 5.8 shows the result of Radon-transforming the point and the surrounding area in Figure 5.7 with different window sizes.

We see that using different window sizes gives us different results. When the size of the window is 70×70 (Figure 5.8a) and we place the origin to in the top left corner, the orientation is−16.9^◦when using variance as an estimator and−17.1^◦ when using kurtosis as an estimator. For a 30×70 window (Figure 5.8b ), the variance gives us

−_16.1^◦, while the kurtosis gives us−_16.3^◦_{. For a 30}×30 window (Figure 5.8c) using the variance gives us−16.0^◦, while the kurtosis gives us−16.8^◦.

Since different window sizes give us different results so deciding which is the best is going to be of interest. The reason for such a difference between the tree windows sizes might be as we discussed earlier that different window size is going to be affected differently by neighboring structures and noise. Transforming a larger area means that we also consider the orientation of a larger area. Intuitively we will, therefore, assume that a smaller area should give to most accurate results. However, this can only be true if we have enough data to make a proper estimate. To determine which of the

(41)

(a) A point selected for Radon transformation

(b) A point selected for Radon transformation. The same image as 5.7a but after one performs vessel enhancement as proposed in [8]. The contrast has been altered for better visualization.

Figure 5.7: An ultrasound image containing the point and the surrounding area we are interested in, before and after Frangi-filtering.

(42)

(a) The point in Figure 5.7 Radon transformed with a 71×71 window around the point. The orientation with maximum variance is whenθ=73.1^◦or −16.9^◦when placing the origin is at the top left corner. The orientation with maximum kurtosis is whenθ=72.9^◦or −17.1^◦when placing the origin is at the top left corner.

(b) The point in Figure 5.7 Radon transformed with a 31×71 window around the point. The orientation with maximum variance is whenθ=73.9^◦or −16.1^◦when placing the origin is at the top left corner. The orientation with maximum kurtosis is whenθ=73.7^◦or −16.3^◦when placing the origin is at the top left corner.

(c) The point in Figure 5.7 Radon transformed with a 31×31 window around the point. The orientation with maximum variance is whenθ=74.0^◦or −16.0^◦when placing the origin is at the top left corner. The orientation with maximum kurtosis is whenθ=73.5^◦or −16.8^◦when placing the origin is at the top left corner.

Figure 5.8: Region of interest in the Radon domain using different window sizes for the surrounding area to estimate the local orientation.

(43)

Figure 5.9: The normalized Radon transform of the point in Figure 5.7 and the 30×30 area around that point. The orientation with maximum kurtosis is whenθ =73.8^◦. In the bottom right corner we see an artifact that did not appear in the regular Radon- transform

estimates are the most accurate, we wanted to measure the orientation of the fascicle that we have select to Radon transform around.

From the image, we select three points in a way that they form a right-angled triangle.

Using the distances between these points and basic trigonometry we can make a rough estimate the angle in the top left corner of the triangle that we have created, we calculated the value to be approximately around 73^◦. The difference between using the variance and kurtosis on a single area is not enough to make a firm conclusion, but from our small experiment, the kurtosis seems to be more consistent. When using different window sizes, it seems that there is not a big difference between using a 30×_{30 or a 30}×_{70 window.}

Since Rana et al. in [7] concluded with that finding the dominant orientation by finding theθ-value with the highest kurtosis on real ultrasound images is the most accurate approach, we are going to do this in this thesis as well. When it comes to selecting a window size, we chose to use a 30×30 window around a point. The reason for this is that we see that there is not a big difference between using 30×30 and a 30×70 window and smaller window size is going to represent the local area better. One thing that a smaller window is more susceptible to is noise. Later in this thesis, we will talk about steps to make the measurements more robust for noise.

From now on we will use a 30×30 window when performing a Radon transform on an area within the region of interest. In the Radon space, we will use theθ-value with the highest kurtosis to determine the dominant orientation.

Performing a normalized Radon-transform

We want to use the normalized Radon transform as it helps us distinguish between the long lines and short lines, as well as making it possible to compare the results with other local Radon Transforms. Now that we know that we are going to use a 30×30 window we should expect the same measurements using the normalized Radon-transform. However when we look at the Radon-image Figure 5.9 and estimate the orientation using kurtosis we see that dominant orientation now isθ =73.8^◦, this is same as −16.2^◦if we set the origin to be in the top left corner. In the Radon-image Figure 5.9 we see a curved-line artifact appearing in the bottom of the image. This curved line did not appear in the non-normalized Radon transform Figure 5.8c.

Taking a closer look at the artifact in the figure we discovered that the normalization

(44)

Figure 5.10: The modified normalized Radon transform of the point in Figure 5.7 and the 30×30 area around that point. The orientation with maximum kurtosis is θ=73.5^◦/−16.3^◦if the origin is in the top left corner

matrix did not properly normalize the values farthest away from the center. To mediate this error, we masked matrix belonging to the Radon transform and the normalization matrix in such a way that the edges did not get included in the estimation. We illustrate the difference in the following code:

%Normalized Radon t r a n s f o r m t h e t a = 0 : 0 . 1 : 1 8 0

[ r , xp ] = radon ( img , t h e t a )

n o r m a l i z a t i o n = ones ( s i z e ( img ) ) ;

n o r m a l i z a t i o n = radon ( n o r m a l i z a t i o n , t h e t a ) ; idx = l o g i c a l ( n o r m a l i z a t i o n ) ;

r_n = r ;

r_n ( idx ) = r_n ( idx ) . / n o r m a l i z a t i o n ( idx ) ;

%Modified Normalized Radon Transform t h e t a = 0 : 0 . 1 : 1 8 0

[ r , xp ] = radon ( img , t h e t a )

n o r m a l i z a t i o n = ones ( s i z e ( img ) ) ;

n o r m a l i z a t i o n = radon ( n o r m a l i z a t i o n , t h e t a ) ; idx = l o g i c a l ( n o r m a l i z a t i o n ) ;

idx ( 5 : end−5 , : ) = f a l s e ; r_n = r ;

r_n ( idx ) = r_n ( idx ) . / n o r m a l i z a t i o n ( idx ) ;

In Figure 5.10 we can see the difference when modifying the normalized Radon- transform. We see that performing this extra masking helps in removing the artifacts.

An the measured orientation is nowθ=73.5^◦or −16.3^◦ if the origin is in the top left corner. The result is now identical to the results achieved using the standard Radon transform.

Limitingθ

Since the orientation of the fascicles usually stays between 30-85 degrees [13] we can limit the Radon-transform and still expect an angle estimates Figure 5.11. This will

(45)

Figure 5.11: The modified normalized Radon transform of the point in Figure 5.7 and the 30×30 area around that point. θ is limited to only span from 30^◦ to 85^◦. The orientation with maximum kurtosis is whenθ =73.5^◦

also make the transform more robust as noisy parts of the image can give ambiguous results.

5.4 Summary

We have now in this chapter talked about the Radon transform and the normalized Radon transform. From now on when performing a normalized Radon transform, we will always mean the normalized Radon transform that we modified. We will also be limiting the Radon-transform only to transform lines that have an orientation between 30-85 degrees, as this reduces computational time and makes the method more robust to a potential influence by the aponeuroses.

When doing some experiments where we tested different proposed methods for estimating the dominant orientation in a Radon transformed image, the different methods did not produce any significantly different results. However, from now on we will use the kurtosis to estimate the dominant orientation in an area. The area around a point we are going to transform is always going to be 30×30.

(46)

(47)

Chapter 6

Aponeuroses detection

We ended chapter 5 by looking at how the Radon-transform represents features in the image space in Radon-space. The author in [13] used this representation to detect and mark where the aponeurosis lies within the image. Because of this, we decided to use the method specific to the aponeurosis detection. In this chapter, we will summarize and discuss the method proposed in [13]. We also will look at possible improvements to the introduced method.

6.1 Baseline algorithm for detecting aponeurosis in single frames

The algorithm for detecting aponeurosis used in [13] is based on the proposed method by [19], [20] and [12]. The algorithm is designed to detect two different aponeuroses in a single frame. We will start this section with a brief introduce the algorithm proposed in [13].

The algorithm starts by finding the approximate location of any potential aponeuroses.

The apporoximate location is found by finding the dominant angle in neach dept of the image and locate sudden changes (in gradient). The area with largest gradient is an indication of the presence of an aponeurosis. This is an assumptiom made since the orientation of the fascicles should be mostly homogenous througout the dept realtive to the orientation of the aponeuroses. For each half of the image the algorithm then looks for a clear peak which indicates the presensce of an aponeurosis. If there is on clear peak in one of the halves that peak selected as the location of the aponeurosis. If There is no clear peak in one of the halves the alorithm performs a new Radon-transform to find the missing one. When Radon transforming the imageθis limited to span from 30^◦to 119^◦, this is based on the assumption that the aponeuroses are not going to have an orientation much more or less than that. If the Radon-space associated with one half has more than one peak the algorithm has to find the peak most likely associated with the aponeurosis. Finding the desired peak is done by thresholding. The peak that is larger than a given threshold is then deemed to belong to the aponeurosis.

Now that the algorithm has approximate knowledge of where the location of the

(48)

two aponeuroses is. The algorithm divides each approximate location into 20 pieces.

Each piece is Radon-transformed. The algorithm then finds the center of mass of each Radon-transformed piece and converts the location back to a cartesian coordinate system.

We now have 20 line pieces that make up one single aponeurosis. Eeach of these line pieces are interpolated together using a least squares spline curve.

6.1.1 Thoughts about the algorithm

When creating the algorithm, the author [13] made some assumptions; some of them will most likely make the algorithm not work on different types of sequences. We have to test how large the extent of these assumptions is going to affect the different sequences.

6.2 Using the baseling algorithm to detect aponeuroses for each frame in a sequence

As mentioned earlier the algorithm proposed in [13] is intended for detecting aponeuroses in single frames. We however want to thest the algorithm on sequences as it might gives us some isnight and help us understand the complexity. We ran the algorithm by using the first two training sequences and the first 100 frames from the third traing sequence.

We notice that the detection accuracy varies from one frame to another, and there are some parts in the proposed method that makes aponeurosis detection not robust enough for every type of frame. The lack of robustness and the assumption that one aponeurosis lies within one half of the image looks to affect the algorithm.

If we look at some images from test sequence one (Figure 6.1), we see that the shallow aponeurosis is easier to find than the deep aponeurosis. The line segments on the deep aponeurosis struggle to match. This mismatch is due to other features in the image have a more substantial accumulated value that to the aponeurosis we want to find. On the shallow aponeurosis, we see that the algorithm more easily can find that particular aponeurosis; however, in some frames, there is some mismatch there also. Since the assumption is that the shallow aponeurosis is somewhere in the upper half of the image the default settings is top place an aponeurosis when the algorithm it looks like it finds any region with the most accumulated response.

For test sequence two (Figure 6.2) we see good detections through the whole sequence, hower since the assumption is that the shallow aponeurosis is in the upper half of

(49)

(a) An example of a good detection

(b) An example of a detection where the deep aponeurosis is incorrectly identified.

Figure 6.1: Image of a frame after performing a aponeurosis detection on test sequence 1.

Using time information to make a robust feature extraction and parameter estimation from unipennate muscles in B-mode ultrasound sequences

Using time information to make a robust feature extraction and

parameter estimation from unipennate muscles in B-mode ultrasound

sequences

Nirusan Tharmanathan

Thesis submitted for the degree of

Master in Digital Signal Processing and Image Analysis 60 credits

Department of Informatics

Faculty of mathematics and natural sciences

UNIVERSITY OF OSLO

Using time information to make a robust feature extraction and

parameter estimation from unipennate muscles in B-mode

ultrasound sequences

Nirusan Tharmanathan

Abstract

Acknowledgements

Contents

Chapter 1

Introduction

Chapter 2

Background

2.1 Ultrasound and ultrasound imaging

2.2 Muscle architecture and movement

2.3 Previous work

2.4 Summary

Chapter 3

Data material

3.1 Data format

3.2 The sequences

3.3 Summary

Chapter 4

Feature enhancement

4.1 A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A.

Viergever, “Multiscale vessel enhancement filtering”

4.2 Summary

Chapter 5

The Radon Transform applied on 2D grayscale images

5.1 The Radon transform

∑

∑

5.2 Normalized Radon transform

5.3 Finding the dominant orientation in the Radon domain

5.4 Summary

Chapter 6

Aponeuroses detection

6.1 Baseline algorithm for detecting aponeurosis in single frames

6.2 Using the baseling algorithm to detect aponeuroses for each frame in a sequence