Automatic landmark detection on Trochanter Minor in x-ray images

(1)

Automatic landmark detection on

Trochanter minor in x-ray images

Examensarbete utfört i Bildkodning vid Tekniska Högskolan i Linköping

av

Per Holm

Reg nr: LiTH-ISY-EX–05/3701–SE Link¨oping 2005

(2)

(3)

Automatic landmark detection on

Trochanter minor in x-ray images

Examensarbete utfört i Bildkodning vid Tekniska Högskolan i Linköping

av

Per Holm

Reg nr: LiTH-ISY-EX–05/3701–SE

Supervisor: Robert Forchheimer Examiner: Robert Forchheimer Link¨oping 11th April 2005.

(4)

(5)

Avdelning, Institution

Division, Department

Institutionen för systemteknik

581 83 LINKÖPING

Datum

Date

2005-02-11

Språk

Language

Rapporttyp

Report category

ISBN

Svenska/Swedish

X Engelska/English

Licentiatavhandling

X Examensarbete

ISRN LITH-ISY-EX--05/3701--SE

C-uppsats

D-uppsats

Serietitel och serienummer

_{Title of series, numbering}

ISSN

Övrig rapport

____

URL för elektronisk version

http://www.ep.liu.se/exjobb/isy/2005/3701/

Titel

Title

Automatisk landmärkesdetektering på Trochanter Minor i röntgenbilder

Automatic landmark detection on Trochanter Minor in x-ray images

Författare

Author

Per Holm

Sammanfattning

Abstract

During pre-operative planning for hip replacement, the choice of prosthesis can be aided by

measurements in x-ray images of the hip. Some measurements can be done automatically but this

require robust and precise image processing algorithms which can detect anatomical features. The

Trochanter minor is an important landmark on the femoral shaft. In this thesis, three di.erent

image processing algorithms are explained and tested for automatic landmark detection on

Trochanter minor. The algorithms handled are Active Shape Models, Shortest path algorithm and

a segmentation technique based on cumulated cost maps. The results indicate that cumulated cost

maps are an e.ective tool for rough segmentation of the Trochanter Minor. A snake algorithm was

then applied which could .nd the edge of the Trochanter minor in all images used in the test. The

edge can be used to locate a curvature extremum which can be used as a landmark point.

Nyckelord

Keyword

(6)

Abstract

During pre-operative planning for hip replacement, the choice of prosthesis can be aided by measurements in x-ray images of the hip. Some measurements can be done automatically but this require robust and precise image processing algorithms which can detect anatomical features. The Trochanter minor is an important land-mark on the femoral shaft. In this thesis, three different image processing algo-rithms are explained and tested for automatic landmark detection on Trochanter minor. The algorithms handled are Active Shape Models, Shortest path algorithm and a segmentation technique based on cumulated cost maps. The results indi-cate that cumulated cost maps are an effective tool for rough segmentation of the Trochanter Minor. A snake algorithm was then applied which could find the edge of the Trochanter minor in all images used in the test. The edge can be used to locate a curvature extremum which can be used as a landmark point.

Keywords: trochanter minor, landmark, edge detection, image processing, seg-mentation

(7)

(8)

Acknowledgment

First of all I would like to thank my supervisor and examiner, Prof. Robert Forch-heimer for all his help and encouragement during the writing of the thesis. I would also like to thank Dr. Anders Rosholm at Sectra Pronosco for starting the project and giving me the opportunity to do this work, Prof. Bj¨orn Kruse for fruitful discussions in the beginning of the work, MSc. Thomas Nordman for his patient waiting to ﬁnally complete his degree by being my opponent. At last, thank you Anna for being so wonderful despite my absent mind and late nights in front of the computer.

(9)

(10)

Notation

Some common notations:

Symbols

x, X Boldface letters are used for vectors, matrices and sets. Φ Matrix of eigenvectors.

λ Eigenvalue

Abbreviations

TM Trochanter Minor ASM Active Shape Models

(11)

(12)

Introduction

Hip replacement has become a common operation in patients with osteoporosis. The femoral shaft is cut and a prosthesis is inserted. The choice of prosthesis can be aided by doing measurements in x-ray images of the hip. During the prosthesis planning, several anatomical landmark points are marked and used to calculate a prosthesis template. Some of these landmarks can be automatically found by image processing methods while some has to be manually marked by a physician. Manual point detection is prone to human deviation and inexactness, so automatic methods is preferred if possible. One important point that is hard to detect automatically is a robust landmark on the Trochanter Minor (TM). The TM is an insertion point for muscles in the thigh and is located below the hip joint. Fig 1.1 shows an x-ray image of a hip with TM marked with circles on both sides.

Figure 1.1. X-ray image of a hip with the Trochanter minors marked with circles.

The TM is an important landmark because it is the only vertical reference point along the femoral shaft in the region. It is used in the prosthesis planning to calculate the cut of the leg to ensure equal leg length. In this thesis, three image processing methods for automatic landmark detection on TM are implemented and evaluated.

(15)

2 Introduction

1.1 Objectives

The goal of this thesis is to develop an image processing algorithm for automatic localization of a robust point landmark on Trochanter minor in x-ray images of the femoral shaft. The algorithm has to be robust enough to handle the large biological variety in shape typical to medical images.

1.2 Method

Several image processing methods were studied and implemented in Matlab. Ex-periments were performed on a set of 20 x-ray images of the Trochanter Minor from diﬀerent persons.

(16)

Chapter 2

Image analysis

In this chapter the basics of image analysis are explained. Methods and concepts later used in this thesis are reviewed.

2.1 What is image analysis?

Images are represented in computers as arrays of numbers. To visualize an image, each number in the array is mapped to a certain gray level. When a person looks at the image it makes sense, we interpret the gray level dots as lines, circles, faces, houses etc, without even have to think about it. The gap between representation (dots on the screen) and description (line, face, house) is bridged in fractions of a second in the human brain. Image analysis research strives to develop methods that make it possible for the computer itself to bridge the gap between image representation and image description. In this thesis, the visual task of locating a certain landmark in similar types of images is to be automated by the computer.

When standing before the task of designing an image analysis algorithm, there is no straightforward way of choosing the best possible method for a specific problem. In the literature there are numerous different approaches to solving all sorts of image analysis tasks. Yet no method has been discovered that can come even close to the human visual system when it comes to retrieving contextual information from images presented to the computer. So far, image analysis algorithms can only be practically helpful in applications where the task is simple and the structure of the objects in the image is known beforehand. When these constraints apply, we can optimize an algorithm for a specific problem but it will only be applicable as long as the initial assumptions about the image are correct.

Below will follow an overview of concepts, tools and methods later used in this thesis.

(17)

4 Image analysis

2.2 Edge detection

Everything that is important when retrieving information from images is associated with changes in the image gray level value from one place to another. Interesting areas always have a border of some kind which means that the gray level changes rapidly on the edge of the object. Therefore edge detection is a fundamental area within image processing.

A fact that is valid on all edge points is that they have high derivative in some direction. Points that fulfill this constraint can easily be located by convolving the image with a derivative kernel, such as the derivative of a Gaussian and applying a threshold. However, this step does not provide any knowledge of the relation between the detected edge pixels. In fact after edge filtering, all we get is just another image where the intensity of the pixels represents edge strength instead of brightness. To gain actual knowledge about the image we have to refine the process to be able to group edge pixels into meaningful structures, such as connected curves around object borders.

2.3 Orientation

All relevant edges in an image has a specific orientation, i.e. they can be horizontal, vertical, or any angle in between. Image orientation is an important property that is not explicitly given in an edge magnitude image. Before we can decide how orien-tation can be calculated, we must define a proper represenorien-tation. One could easily believe that the angle of the image gradient is sufficient but this representation has two problems that need to be handled:

• The angle of the gradient is 180 degrees diﬀerent depending on if the edge is

positive or negative. The orientation representation has to be invariant to if the edge is positive or negative.

• A numerical angle measure wraps at 360 degrees, thus edges oriented almost

the same but on diﬀerent sides of the wrap point will have very diﬀerent measures. We want the orientation representation to be continuous.

An orientation representation that handles these issues is the double angle vec-tor of the gradient response [2]. It is convenient to use a complex representation for the vector, which gives this expression for the orientation:

o(x, y) = (gx(x, y) + igy(x, y))2 (2.1)

Where o(x, y) is the orientation at point (x, y) and gx, gy are the x- and

y-components of the gradient response. The square doubles the angle of the vector. The norm of the vector is not important in means of orientation but the expression can be divided bydx2+ dy2 to make the norm equal to the edge magnitude. This gives a uniﬁed representation of edge magnitude and orientation. Fig 2.1 shows the behaviour of the orientation representation around a circle.

(18)

2.4 Utilizing ´a priori knowledge 5

Figure 2.1. Double angle representation for orientations around a circle. Note that

opposite sides of the circle has the same orientation despite oppositely directed edges.

With this representation, orientations can be compared in a continuous way. The diﬀerence between two orientations can be represented as a real angle, φ, and is calculated as:

φ =argo1 o2

(2.2)

Note that the argument of a fraction between two complex numbers is the diﬀerence in argument between the nominator and the denominator.

2.4 Utilizing ´

a priori knowledge

The main reason why humans are so good at interpreting images is that we have the ability to compare what we see with memories of things we have seen before. For example, we can say that we have a basic model in our memory of what a car looks like. When we see a new sort of car that we have never seen before, we can still recognize it as a car because our model generalizes to ﬁt the new car. When we have decided what model to use to describe what we see, we can analyze the object by for instance locating the wheels or the windows of the car by using prior knowledge from the model in our memory.

This model approach can also be implemented in computer vision in a variety of ways. In chapter 4 ´a priori information is collected from training examples and a model of the shape variation is built. In chapter 5 and 6, characteristics speciﬁc to all TM-images are used as prior information to the algorithms.

The opposite to model based methods is data driven methods. In data driven methods no model is used but the image is pre-processed to reveal features of interest, such as edges, lines or more complex patterns. The goal is to reﬁne the information in the pixels so that interesting structures easily can be distinguished.

(19)

6 Image analysis

2.5 Landmark detection

In this thesis the ﬁnal goal of the algorithm is to locate a robust landmark on the Trochanter minor. A landmark is a distinguishable point that is present on each example image. For instance, good landmarks in images of faces could be the corners of the eyes and the tip of the chin. Typically, good landmark points are boundary corners, T-junctions and boundary edge points with high curvature.

The procedure of locating a landmark is diﬀerent depending on the pre-processing of the data. In a model based approach the concept of landmarks is in fact an es-sential part of the method itself. The model is built up by manually annotating landmarks. When a data driven approach is used, the landmark detection is a post-processing step where a search is performed in the pre-processed image.

(20)

Chapter 3

Problem description

In this chapter the image processing problem to be solved is explained.

3.1 X-ray images of Trochanter minor

The x-ray images that are input to the algorithm are gray level images with 8 bits/pixel. The region of interest i.e. the TM is fairly well defined beforehand by information from the orthopedic software. The images can vary in a number of ways which puts demands on the algorithm to be invariant to the differences and still locate a common landmark. In fig 3.1 are shown four typical input images. Lots of false edges are present, emerging from details in the muscle region outside the bone. The gray level of the bone can not be used alone to distinguish the TM, sometimes the muscle region has intensity equal to parts of the TM.

3.2 A robust landmark

To find a robust landmark on the TM, we must first decide which point to choose. This is not an easy task in this case because of the large variation in shape and the lack of distinct points such as crossings and sharp corners. The gray level inside the T.M. does not have any significant features, the only reliable structure present in all examples is the edge towards the muscle region. If that edge could be detected, its curvature can be analyzed. A possible landmark could then be the upper curvature extremum of the edge. The upper bend is always sharper than the lower and thus a more precise landmark. Fig 3.2 shows the proposed landmark on the detected edge.

Once the edge is located and translated into a curve, (x(t), y(t)), there are nu-merous methods in the literature for calculating curvature and locating extremals. See for example [3]. However the detection of the correct edge showed up to be complicated enough for the scope of this thesis, so in the reminder only methods for correct edge detection is handled.

(21)

8 Problem description

Figure 3.1. Typical input images. The algorithm has to be invariant to the diﬀerence

in shape and appearance and still locate a common landmark.

Figure 3.2. The upper extremum in curvature could be used as a feasible landmark if

the edge can be located.

3.3 Method proposal

The first method implemented and tested was a landmark based method called Active Shape model. It was chosen because of its ability to use á priori information learnt from training examples. This is very favorable in the case of large variation in shape and appearance of the input data, typical to medical images. In [6], ASM is successfully used to locate finger bones in x-ray images of the hand. Active shape models are handled in chapter 4. In chapter 5 and 6 graph optimization techniques

(22)

3.3 Method proposal 9

are used to locate the edge towards the muscle region in the TM-images. These methods focus on the fact that the TM-border is an edge that reaches through the whole image. The optimization tools are used to ﬁnd optimal paths through the image with respect to edge properties.

(23)

(24)

Chapter 4

Active Shape models

Active shape models is a method for locating shapes in images by utilizing prior knowledge learnt from training examples. A shape model is built up by annotating corresponding point landmarks in a set of images representing diﬀerent instances of the object. In the search procedure the algorithm compares instances of the model with image data and updates its parameters in order to minimize the error.

4.0.1 Shape representation

In active shape models the built-in knowledge consists of a statistical model of the objects shape variation. To be able to handle such a model we have to clarify what we mean by shape.

D.G. Kendall [5] states this deﬁnition

Deﬁnition 4.1 Shape is all the geometrical information that remains when

loca-tion, scale and rotational eﬀects are ﬁltered out from an object.

The location, rotation and scale of a shape is called the pose of the shape. To be able to compare and model the variation between a number of shapes the calculations has to be carried out in a pose invariant representation.

One way of describing a shape is to annotate a number of landmarks around the border. Fig 4.1 shows the same shape under diﬀerent poses.

For a shape with n landmarks, the x and y coordinates of the landmarks can be concatenated into a 2n-point vector as a form of shape representation.

x = [x₁, x₂, . . . , xn, y1, y2, . . . , yn] (4.1)

However, this representation also includes the pose of the shape. To obtain a true shape representation, location, rotation and scale has to be ﬁltered out. One way of doing this is carried out by aligning the shapes to a common coordinate reference with respect to location, rotation and scale. A commonly used method of aligning shapes is called Procrustes analysis and is described below. After aligning, shapes can be analyzed and compared in means only by the inﬂuence of shape.

(25)

12 Active Shape models

Figure 4.1. The same shape under diﬀerent poses.

Procrustes analysis aligns the shapes so that the sum of squares diﬀerence between each shape and the mean shape is minimized. This can be expressed as min D = |xi− ¯x|. For this measure to be useful, constraints have to be

applied on the mean during the alignment. Typically we want the mean shape to be centered at the origin and have unit scale. The alignment of a set of shapes can be done as follows:

1. Translate each example so that its center of gravity is at the origin.

2. Choose one example as an initial estimate of the mean shape and scale so that|¯x| = 1.

3. Record the ﬁrst estimate as ¯x₀to deﬁne the default reference frame. 4. Align all the shapes with the current estimate of the mean shape. 5. Re-estimate mean shape from aligned shapes.

6. Apply constraints on the current estimate of the mean shape by aligning it with ¯x₀ and scaling so that¯x| = 1.

7. If not converged, return to 4.

Convergence is declared if the estimate of the mean does not change signiﬁcantly after an iteration.

In step 1 the center of mass is calculated as:

(¯x, ¯y) = ₁ n n j=1 xj, 1 n n j=1 yj (4.2)

To perform step 2 we need to establish a shape metric. The most common metric used is the frobenius norm or 2-norm:

S(x) = n j=1 [(x_j− ¯x)2+ (y_j− ¯y)2] (4.3)

(26)

4.1 Shape modeling 13

In step 4 the actual alignment to the mean is carried out for one shape at a time. Alignment here means to minimize the difference between each shape and the mean shape by changing the pose parameters of the shape to be aligned. The operations allowed during the alignment with the mean will affect the final distribution of the shapes. For example, if the norm of each shape is fixed to 1 before aligning with the mean, all shapes will lie on a hypersphere, which will introduce unnecessary non-linearities. It is shown in [1] that the best approach is to allow scaling and rotation during the alignment and then project the result into the tangent space to ¯x by scaling with 1/x.¯x.

The alignment w.r.t. scaling and rotation is done by ﬁnding the transformation

T, that minimizes the expression|T(x) − ¯x|, where T is a transformation matrix

allowing scaling and rotation. The correct T is obtained by diﬀerentiating the expression for the error and equating to zero. The calculation will result in a linear system of equations. See Appendix B in [1].

When this is done, we have a representation which includes only the true shape of the objects. With this at hand we can start analyzing the shape variation and build a model to describe the group of shapes.

4.1 Shape modeling

Suppose now we have s sets of points xi, aligned into the pose-free representation.

These shape vectors form a distribution cloud in the shape space. The size of the cloud depends on the variation of shape between the annotated examples. We want to build a model of that distribution to be able to generate new examples similar to those in the training set. To be able to adjust the model to handle the variation within the set, we need to give it a number of parameters that allows for its degrees of freedom. The model can be expressed like x = M(b) where M is the model function, b is the model parameters and x is the resulting model instance, i.e. the shape vector x that the model has generated.

The original dimensionality of the shape space which reside all possible 2D-shapes with n points is 2n-4 (2n for the x- and y- coordinates minus the four pose dimensions). We want to reduce the dimensionality in our model to span only shapes which lie within the point cloud of our training set. A very eﬀective method for reducing the dimensionality of a dataset with minimal loss of information is PCA, Principal Component Analysis. PCA produces an orthogonal basis for the shape space where the base vectors correspond to the major axes of variance in the point cloud. The PCA components (base vectors) are the eigenvectors to the covariance matrix of the dataset. The approach is as follows.

(27)

1. Compute the mean of the data (The mean shape vector),

¯ x = 1 s s i=1 x_i (4.4)

2. Compute the covariance matrix of the dataset,

S = 1 s− 1 s i=1 (x_i− ¯x)(x_i− ¯x)T (4.5)

3. Compute the eigenvectors φ_iand the eigenvalues λ_i of the covariance matrix

S and sort them so that λ_i≥ λ_i+1.

If we take t, t < s of the eigenvectors corresponding to the largest eigenvalues and arrange them side by side into a n by t matrix Φ, we can approximate any shape

x in the training set by

x≈ ¯x + Φb (4.6)

where b is a t-dimensional vector given by

b = ΦT(x− ¯x) (4.7)

Algebraically this can be explained as approximating the shape vectors in the training set with their projection into the best possible subspace with t dimensions. The projection matrix in this sense is Φ while b is the coordinate vector in the new subspace basis spanned by [φ_i]t

1

We now have a model of our shapes in the training set with t parameters. The deformability of the model depends on the choice of t. How do we know how many parameters we need? The variance of the dataset in the direction of the ith

eigenvector is given by λi. As the components in the PCA basis are orthogonal,

the correlation between them is zero. This means that we can add components by summing the sorted eigenvalues until we reach a desired (e.g. 98%) part of the total variance.

One might wonder why the eigenvectors to the covariance matrix is the opti-mal subspace basis. The solution to this can be found by solving a minimization problem for the correlation of the data in the direction of the basis vectors, it will not be proven here.

(28)

4.2 Fitting a model to new points 15

1

2 b

φ

Figure 4.2. PCA applied to a 2 dimensional dataset. The points can be well

ap-proximated with their projection onto the largest eigenvector Φ1 and modelled with the

parameter b along that vector.

4.2 Fitting a model to new points

Suppose now we are presented a new set of points in an image, a new shape, and want to ﬁt our model to those points. The model x = ¯x + Φb does only model

shape, not pose, it means that apart from b, we must also ﬁnd the pose parameters translation (Xt, Yt), rotation θ and scale s, that makes the model instance as similar

to the new points as possible. This can be expressed as a similarity transform of the shape model:

x = T_X_t_,Y_t_,s,θ(¯x + Φb) (4.8) Where the function TXt,Yt,s,θ performs a translation by (Xt, Yt), a rotation by

θ and a scaling by the factor s of the shape model instance ¯x + Φb. Thus, the

optimization of the model ﬁt depends on 4 + t parameters.

For one single pointx_y, the similarity transform TXt,Yt,s,θ can be expressed as

TXt,Yt,s,θ x y = Xt Y_t + s cos θ s sin θ −s sin θ s cos θ x y (4.9)

To ﬁnd the best possible shape and pose parameters to ﬁt to new image points

Y is equal to minimizing the sum of squares diﬀerence between the corresponding

points, i.e. minimizing the expression

|Y − TXt,Yt,s,θ(¯x + Φb)|2 (4.10)

If we want the model x = ¯x + Φb to generate examples similar to those in the

training set, we must choose b from a distribution learnt from the training set. The cloud of all shapes, x, in the training set corresponds to a distribution cloud of

b-parameters through their projection into the model space. We want our choices

(29)

set is diﬃcult to model but if we assume that the bi’s are independently gaussian,

hard limits for each bi can be set to i.e.

|bi| ≤ 3

λi (4.11)

An iterative approach for solving the minimization 4.10, is, according to Cootes [1], as follows:

1. Initialize the shape parameters,b, to zero. 2. Generate the model instance x = ¯x + Φb

3. Find the pose parameters (Xt, Yt, s, θ) which best map x to Y (See Appendix

B in [1]).

4. Invert the pose parameters and use to project Y into the model coordinate frame:

y = T_X−1_t_,Y_t_,s,θ(Y) (4.12) 5. Project y into the tangent plane to ¯x by scaling by 1/(y.¯x)

6. Update the model parameters to match to y

b = ΦT(y− ¯x) (4.13) 7. Apply constraints on b. (Truncate b so that 4.11 holds)

8. If not converged, return to step 2.

Convergence is declared if the model parameters does not change signiﬁcantly after an iteration, it usually just take a few iterations.

The iterative approach is needed because of the truncation of b. If b is truncated during an iteration, the pose parameters might have to be changed for an optimal ﬁt.

Note that the new shape, Y, must have equally many points which must corre-spond to the same landmarks as the shapes in the training set. This in fact is one of the drawbacks with the method. Classes of shapes where the shape variation is such that distinct landmarks are hard to distinguish, are also hard to model.

4.3 Image interpretation

We have now seen how to build a model of a certain class of shapes, and how to ﬁt this model to a new set of points. In this section we will show how to search in the image for candidate points and iteratively ﬁt the model to the new points.

Given an input image with an object we want to locate, the goal of the algorithm is to ﬁnd the best pose and shape parameters that ﬁt the model instance to the image object. This is an optimization problem which can be seen as a minimization

(30)

4.3 Image interpretation 17 Sampled profile Model point Image object Sample Intensity Profile sample

Figure 4.3. The search is done by sampling along the normal to the model points, and

look for structures similar to those in the training set.

of the difference between the model and the underlying image. So far, the shape model only includes geometrical information about the object and no connection is made to the pixel values, the appearance of the object. To determine if a suggested model instance is a good fit and to propose a direction of search, such a connection has to be made. Objects in images are most likely limited by an edge or line. This means that when the algorithm has reached its goal, the shape model points will be located at distinguishable image points, deviating from the neighbourhood in the direction normal to the object. We can utilize this in the search by looking for structures around each model point similar to those in the training set. In active shape models the training images are sampled along the normal to the annotated shape in each point. A mean sample is calculated for each point. The search is done by sampling a longer piece in the same way and update the model point location with the best fit along the sample. Given a rough starting approximation, the model instance can be iteratively fit to the image using this approach:

1. Initialize a model instance roughly in the image.

2. Sample the image along the normal to each model point and ﬁnd the best match for the mean sample from the same point in the training images. Record the new proposed shape.

3. Update the model parameters (Xt, Yt, s, θ), b to best ﬁt the new points.

(Ac-cording to 4.2)

4. Repeat until convergence, which is declared when an iteration does not change the model parameters signiﬁcantly.

(31)

When creating the normal samples in the images, the derivative instead of the gray level is sampled to avoid dependence on global gray level in the images. This approach is used:

1. Sample the image k pixels on either side of the model point i into a 2k + 1 point vector g_i.

2. Compute the derivative (point to point diﬀerence) of gi.

3. Normalize the sample by dividing through with the sum of absolute element values,

1

ij|gij|

g_i→ g_i (4.14)

This is repeated for all points of all shapes in the training set and a mean sample can be computed for every model point.

4.4 Implementation

The active shape model was implemented in Matlab and applied to a set of TM-images. The ASM code was divided into three sub programs:

• Annotation program used to collect the training data.

– Input: Set of images

– Output: Set of raw shapes and sampled proﬁles

• Shape statistics generator

– Input: Set of raw shapes

– Output: ¯x (Mean shape), Φ (Eigenvectors), θ (Eigenvalues)

• Search algorithm

– Input: ¯x, Φ, θ, Start shape, Sampled proﬁles, Input image – Output: Resulting model instance, hopefully around image object

(32)

4.4 Implementation 19

4.4.1 Annotation

Figure 4.4. Annotated examples in the training set.

5 10 15 20 25 30

2 4 6 8

Figure 4.5. Mean of sampled proﬁles over all images. Each sample is one column, there

are 32 columns, one for each model point. Note that column 16 to 22, from the right edge, is oriented in the opposite direction.

Fig 4.4 shows three annotated TM-images. The program was made to handle non continuous shapes, i.e the shape can be divided into several sub shapes, all in the same covariance matrix. This makes it possible to include information from all structures specific to the examples, even if they are not along the same border. If only the outer edge was used in the shape model, the possibility of the algorithm to converge at a wrong edge was increased. The only implementational issue when allowing non continuous shapes is that the angle for the normal sample must be specially handled at all ends (perpendicular to the end segment). The lack of distinct points in the examples makes it difficult to manually annotate robust point landmarks that correspond between the examples. The annotation was made by fitting a spline polynomial through the user annotated points in each sub shape and place the actual model points equally spaced along the spline.

(33)

4.4.2 Shape statistics

Figure 4.6. Shape statistics for training set.

The shape statistics part aligns the raw shapes and calculates a mean shape, eigenvector basis and eigenvalues, according to 4.1 and 4.2. When the calculation is done, a result window (Fig 4.6) is shown with the following information:

• Mean shape with ellipses in each model point representing one standard

de-viation in x- and y-direction.

• Color coded covariance matrix. Bright means high variance. The values in the

diagonal is the actual point variances whereas the other elements represent covariance between diﬀerent points. The upper half represent the x-values and the lower part the y-values.

• Mean shape and ﬁrst mode of variance. This is the major axis shape variance

in the model,± one standard deviation.

• Distribution of eigenvalues. The eigenvalues to the covariance matrix sorted

in decreasing order. For 98% of the variance to be covered, 8 eigenvectors had to be included in the model basis.

(34)

4.4 Implementation 21

4.4.3 Search algorithm

Figure 4.7. (a) Badly initialized start points. (b) Circles represent new proposed position

for each landmark based on a search along the normal. (c) Proposed shape projected into the model space (d) Model after 5 iterations (e) Model after 60 iterations

The input to the search algorithm is the image to be searched, the output from the shape statistics program, a starting shape and the mean sampled proﬁles from the annotation program. The algorithm implements 4.3.

In fig 4.7 the progress of the algorithm can be seen. The start shape in (a) is quite badly initialized to show the effect of the projection into the model space in (c). Note how the obviously mislocated points are moved into a more decent position. After 5 iteration in (d), the shape is still wrong in the lower region because the inner edge is found instead of the correct outer edge. In this case it took 60 iterations for the model to converge to the result in (d). The shape is not perfectly fit along the TM. When this happens in an image within the training set, which should be known to the model, it is either because of the reduction of basis vectors to include only 95% of the variance, or because the model points probably have not found the same points as was annotated during the training. The latter is quite probable in the TM-images because of the lack of horizontal features. The model can easily slide up or down and still locate strong (but wrong) edges for all model points. The lack of a vertical reference points also makes it hard to initialize the model so that the model points are somewhat vertically correct. A horizontal edge somewhere in the region would have helped to solve this problem.

(35)

4.5 Results

Figure 4.8. Result of ASM-search in six diﬀerent images.

A search was performed on images both within the training set and on new images. In 4.8 are shown the result of six searches. The search was initialized with the mean shape and the mean pose from the training set. As can be seen the algorithm is far from perfect. The model points have found the wrong image points in several examples.

(36)

4.6 Discussion 23

4.6 Discussion

Before stating that the method does not work, there are several things that need to be accounted for:

• The size of the training set is not as large as it should be, this limits the

variability of the ﬁnal model.

• The choice of model points has not been optimized. There could be other

ways of placing the model landmarks that give better results.

• The implemented version of the ASM-method is quite basic. There are lots of

reﬁnements in the literature that could be applied to the algorithm, includ-ing better modellinclud-ing of the b-parameter distribution, mahalanobis distance measure on the sampled proﬁles etc.

• The initialization of the model in the image could be done better by utilizing

(37)

(38)

Chapter 5

Shortest path algorithm

When segmenting the Trochanter minor the most important feature is the edge towards the muscle region. This edge is smooth and connected from the top of the image to the bottom. In this chapter a shortest path algorithm is used to ﬁnd the optimal path between the top and the bottom of the image with respect to edge strength. This is a true data driven method. The result is a direct function of the pixels in the image.

5.1 Shortest path

A graph in the context of discrete mathematics is a set of points M and a set of lines connecting them. The points are often called nodes and the connecting lines are called edges, or arcs if they have a given direction. Let (i, j) be an arc that exits node i and enters node j. In this chapter we will look at a special type of graph, where all arcs (i, j) are assigned a number, or cost, w(i, j).

B

D

E

F

G

C

2 5 4 9 1 8 2 6 7

A

Figure 5.1. An acyclic digraph. The numbers on the arcs represent the cost of moving

through the arc. We want to go from node A to node G as cheap as possible.

If we deﬁne this number to be the cost of moving between the nodes connected

(39)

26 Shortest path algorithm

to the arc, we can state the problem of ﬁnding the cheapest (optimal) path between two nodes in the graph. A path is a set of arcs {(x1, x2), (x2, x3), ..., (xn−1, xn)}

such that all xi are nodes in the graph and no arc is crossed more than once.

Let P (i, j) be the optimal path between node i and node j. Optimal means that

(k,l)∈P (i,j)w(k, l) is the lowest possible sum of costs over all possible combinations

of arcs between node i and j. Algorithms solving this kind of problems are often referred to as shortest path algorithms.

When using this framework on images to ﬁnd an optimal edge, we will let the pixels of the image represent arcs. The cost of each pixel will be proportional to the inverse of it’s edge magnitude. In this way the optimal path will be the one passing through as strong edge pixels as possible the whole way.

In ﬁg 5.1 a simple graph is shown. The graph in Fig 5.1 holds a special property that we will use when describing an image. There are no paths that start and end in the same node, i.e no loops. Graphs with this property are called acyclic digraphs. In particular for acyclic digraphs, the nodes can be labeled so that this deﬁnition holds.

Deﬁnition 5.1 A graph is an acyclic digraph ⇐⇒ (i, j) exists ⇒ i < j

In other words, the nodes can be ordered so that all arcs will point to a node with higher order number. It can easily be veriﬁed that this is true for the graph in ﬁg 5.1 where the nodes are ordered from A to G.

Before we present the general algorithm of ﬁnding the shortest path in an acyclic digraph we need a fundamental theorem in optimization, Bellman’s principle of

optimality.

Theorem 5.1 If the optimal path P (i, k) between node i and node k passes through

node j, then P (i, k) = P (i, j) + P (j, k).

The theorem states that an optimal path can be split up in two optimal sub paths. The consequences of this is that we don’t have to do an exhaustive search to ﬁnd the optimal path but can divide the problem into smaller sub-problems. The principle of optimality together with the property of acyclic digraphs makes it possible to solve the shortest path problem in O(N ), where N is the number of nodes. This is very fast compared to general algorithms for more complex graphs, e.g. Dijkstras algorithm is a popular shortest path algorithm with time complexity

O(N2).

Suppose now we have a start node with order number 1 in an acyclic digraph and that the nodes are ordered so that deﬁnition 5.1 holds. We want to ﬁnd the optimal path P (1, i) from the start node to any other node i in the graph. Let c(i) be the cost of moving from the start node to node i along the optimal path, that is,

c(i) =

(k,l)∈P (1,i)

w(k, l) (5.1)

Assume we know c(i) for all i ≤ k. This is in fact all we need to extract the optimal path P (1, k). It is done by starting in node k and backtrace in the graph.

(40)

5.1 Shortest path 27

In node k, the only connecting arcs that can start taking us back to node 1 are the ones that lead backwards in the ordering, {(i, k) : i < k} (acyclic digraph property). As we have assumed to know c(i) for i < k we can get the optimal arc (i, k) in P (1, k) by choosing

i = argmin{c(i) + w(i, k) : (i, k) exists} (5.2)

By repeating this procedure until we have reached node 1, the optimal path P (1, k) can be extracted from the known total costs c(i).

We shall now prove that we can generate the total costs c(i) for all nodes in the graph by induction.

• The total cost c(1) of the start node is trivially 0.

• Assume we know c(i) for all i ≤ k, then we also know the optimal paths P (1, i)

as we have seen above. Now look at c(k + 1) and P (1, k + 1). According to the principle of optimality,

P (1, k + 1) = P (1, i) + P (i, k + 1) (5.3) for all nodes i≤ k belonging to P (1, k+1). Particularly it is true for the node connecting to k + 1 so that P (i, k) simply becomes one of the arcs entering node k. This means that we can reduce the expression for c(k + 1) to

c(k + 1) = min

{i:(i,k) exists}(c(i) + w(i, k + 1)) (5.4)

Note that i is only taken over all lower order neighbours to k so c(k + 1) is calculated as a min of only a few sums of node total costs and arc costs. This is what makes the algorithm as fast as O(N ).

In ﬁg 5.1 the calculation of the total costs would be: 1. c(A) = 0 2. c(B) = min (c(A) + 2) = 2 3. c(C) = min (c(A) + 4) = 4 4. c(D) = min (c(A) + 5) = 5 5. c(E) = min{(c(B) + 9), (c(C) + 1)} = 5 6. c(F ) = min{(c(C) + 8), (c(D) + 2)} = 7 7. c(G) = min{(c(E) + 6), (c(F ) + 7), (c(F ) + 1)} = 11

In ﬁg 5.2 the total costs are written above each node, and the optimal path to node G is marked with ﬁlled arrows.

(41)

28 Shortest path algorithm F G C 2 5 4 9 1 8 2 6 7 2 4 5 5 7 11 E A B D

Figure 5.2. The calculated total costs is written above each node and the optimal path

(ﬁlled arcs) can be backtraced from nodeG.

5.2 Edge detection using graph searching

We will now see how to put an image into the graph framework. In fig 5.3 a fictive edge image are described as an acyclic digraph with its optimal costs and path extracted in the right figure. The pixel values are inverted so that high edge magnitude gives a low cost. The cost values are written in each node, this means that the cost are valid for all arcs entering the node, this is the situation when edge magnitude is the only property used when calculating the cost.

In the example each pixel is allowed to connect to any of its three right side neighbours. This can be expanded vertically to allow for higher curvature in the paths. 6 3 7 0 3 3 1 4 2 8 6 3 8 9 0 1 5 2 2 3 8 6 5 6 5 9 0 1 1 1 11 3 0 1 24 3 8 1 0 6 2 1 1

Figure 5.3. (a) A model of an edge image in a graph framework. Low cost numbers

represent strong edges.

In ﬁgure 5.4 the algorithm is applied to two TM-images. The costs are cal-culated from the top of the image to the bottom. The connectiveness has been expanded to 7 pixels to allow for higher curvature. In the left example the path is fairly well found except for the ”short cut” in the lower edge junction. As the cost of making a strong turn across an edge is the same as going straight along the edge, the path may choose to take this kind of unwanted short cuts. In the right image the wrong edge is found because it is stronger. Clearly we have to build more knowledge into the algorithm to make it useful to ﬁnd the correct edge.

(42)

5.3 Applying constraints 29

Figure 5.4. Result of optimal path algorithm applied to two TM edge images.

5.3 Applying constraints

We have seen that a simple optimal path is in fact very sub optimal for the use of ﬁnding the outer edge on the TM. The inner edge is sometimes much stronger and the path tend to switch between two strong edges despite the non edge area in between. To handle this, more knowledge has to be built into the cost function. In particular we want to:

• Favour outer edge

• Suppress unwanted short cuts between two edges

In order to prevent unwanted short cuts through non edge areas, we need to punish such moves in the cost map. This can be done by taking into account not only the edge magnitude but also the image orientation. An arc having the same direction as the local image orientation should be cheaper than an arc crossing an edge.

To perform an orientation dependent cost calculation for the arcs we need to compare the orientation of the arcs with the orientation of the image. This can be done using the complex double angle representation of orientation. First we need to deﬁne the orientation of the arcs in this representation. If the arcs are allowed to connect to pixels within ±D steps vertically, there are 2D + 1 orientations to calculate. These are the same anywhere in the image so they only have to be calculated once.

The arcs (i, j) can be seen as vectors from pixel i to pixel j. j is always

one column ahead of i. If describing the arc-vectors as complex numbers, the 2D + 1 possible vectors are 1 + ni, −D ≤ n ≤ D, where the imaginary part is the vertical translation in pixels. We want to compare the direction of the arcs with the orientation of the image. However image orientation is deﬁned as the double angle vector of the vector pointing in the direction of the gradient, that is, orthogonal to the direction along the edge. To make the direction along the edge cheap, the

(43)

30 Shortest path algorithm

arc-vectors need to be rotated 90 degrees before comparing with the orientation. Multiplication with i generates the desired rotation, yielding this expression for the orientation of the arc-vectors:

oarc(n) = (i(1 + in))2= (i− n)2 (5.5)

Note that the square takes us to the double angle representation.

There are two pixels connected to each arc, so the mean orientation seems like a good value to compare with. Let o_pixel(i) be the complex double angle orientation of the pixel with node number i. The measure of diﬀerence between the orientation of arc (i, j) and pixels i and j can thus be calculated as

o_{dif f} =argopixel(i) + opixel(j)

o_arc(n)

(5.6)

Here o_{dif f} ∈ [0, π] is the absolute value of the diﬀerence in argument between the mean pixel orientation in the nominator and the arc orientation in the denominator.

The cost value for each arc can now be calculated as:

w(i, j) = F (e(j), odif f) (5.7)

where F is some increasing function. It was found that choosing F to that

w(i, j) = e(j)(1 + αoβ_{dif f}H(odif f − φ)) (5.8)

produced good results. H is the Heaviside function, this means that the cost resulting from orientation difference will only exist when the difference between arc and image orientation is larger than φ. The reason why this is needed is that in the worst case, the smallest possible difference in angle between two consecutive arcs isπ₄. It is thus impossible for the path to always follow the exact orientation of the image in each step. If small differences in orientation would be expensive, edges not oriented exactly in any of the 2N + 1 arc directions would also be expensive. This is because the path has to go in small zig-zag steps to follow the edge. A suitable value for φ was shown to be a bit smaller than π₂ (π₄ angle difference due to double angle representation). The value _2.5π was used.

The constants α and β are used to tune the impact of the orientation difference. High values for β will increase the cost more the greater the difference in orientation. The values α = 10 and β = 2 was found to work well. If too high values was used, the path tend to be piecewise straight and not follow the edge smoothly (zig-zag moves gets too expensive). Too low values allows the path to take shortcuts across edges. The reason why e(j) is multiplied with the orientation dependent part is that the impact of the orientation difference becomes amplified if the edge magnitude is strong.

The result of the orientation dependence can be seen in fig 5.5. It is clear that the path better follows the edge with less tendency to jump between different edges. However we still need to make sure the correct edge is found which will not be the case if the wrong one is stronger. In fig 5.6 is an example of a problematic image where the wrong edge is much stronger. No robust method of solving this problem has yet been found.

(44)

5.4 Discussion 31

Figure 5.5. Path extracted from cost map in two images. In the left image in each pair,

no orientation dependence is used and in the right image orientation is accounted for as described above. Clearly the path in the right images follows the edge more smoothly.

Figure 5.6. Problematic image where the wrong edge is much stronger than the correct

edge.

5.4 Discussion

We have seen that the result of the optimal path algorithm can be improved by including the image orientation in the cost calculation. However it is still uncertain which edge is found, this problem needs to be solved before the shortest path algorithm can be used to locate the correct TM-border. Some experiments were performed with a weighting function applied to the edge image to amplify edges to the left. This method worked poorly where the edges were close and the diﬀerence in impact of the weighting function became too small.

(45)

(46)

Chapter 6

Segmentation with cost

function

In this chapter a method is presented that uses the cost map from the shortest path algorithm to segment parts of an edge image that consist of connected edges across the whole image. This segmentation technique is used to eliminate spurious edges outside the interesting bone edge. When the segmentation is done, other edge ﬁnding methods can be roughly initialized in a robust way and work without disturbance from false edges. The segmentation technique is a data driven step. After this, a model based snake algorithm is applied to ﬁne tune the detection.

6.1 Two way cost map

The shortest path algorithm in chapter 5 only used a cost map generated in one direction. The value of the pixels in that map reﬂected the cost to move along the optimal path back to the edge where the cost map calculation started.

Here we will utilize information from maps generated in both directions to get a map where each pixel represent the cost of moving through the pixel from one side to the other. In this way only pixels that are members of a good path through the whole image has a low cost. Spurious edges emerging from the bone edge and deviating into the muscle region will be covered with high cost in the map because they probably are not part of a complete edge path through the whole image.

The algorithm consist of the following steps: 1. Create an edge image.

2. Apply a threshold so that only 10% of the strongest edge pixels are left. 3. Set the cost of the edge pixels reversely proportional to the edge strength and

the cost of the removed pixels to inﬁnity.

4. Calculate a cost map in both directions, from top to bottom and vice versa.

(47)

34 Segmentation with cost function

5. Add the two cost maps.

6. Create a binary mask representing pixels in the map where the cost is ﬁnite, this will be all pixels that are part of a connected path of pixels through the whole image.

The steps of the algorithm are explained in ﬁg 6.1.

Figure 6.1. (a) An artiﬁcial segmented edge image with several possible connected paths

through the image as well as edges not part of a complete path. Light pixels are more expensive than dark, while the cost of the background pixels (white) is inﬁnity. (b) Cost map calculated from top down. Dark means low cost. (c) Cost map calculated from bottom up. (d) Sum of both cost maps. Only edges part of a complete path through the image has ﬁnite cost.

This method relies on the fact that an edge segment that ends in the middle of the image will always have infinite cost in one of the two cost maps (see (b) and (c) in fig 6.1). The sum of the maps thus becomes infinite in every point that not can be reached from both sides without crossing a non edge area.

Another interesting property of the sum of the two cost maps is that the sum is constant in every pixel of the optimal path. Once both cost maps has been calculated the optimal path can be found by searching for the minimum cost in every row of the image. The values in the two way cost map represent the cost of the cheapest possible path going through each pixel. It can thus be seen as something like a global edge map, which can be analysed like a normal edge image with the diﬀerence that every pixel reﬂects more than just local edge strength.

Note that the thresholding of the edge image is a crucial step for the result. By relaxing the threshold, weaker edges will have the possibility to be part of a complete path, whereas a high threshold increases the risk of gaps in a real edge, making it non complete through the image and thus inﬁnitely expensive.

The threshold was set to 90% of the cumulated histogram:

T hres n=1 H(n) = 0.9 N n=1 H(n) (6.1)

(48)

6.2 Segmentation results on TM 35

Where H(n) is the histogram of the edge image. Also two ad-hoc steps were added to the algorithm to improve the suppression of unwanted edges in the TM-images. Before the cost map calculation, the strongest edge point in both the top and bottom of the image was located, then the rest of the row except 10 pixels around the point was set to zero. This forces all paths to start and end in the 10 pixel window around the strongest edge in both sides of the image. To reduce the amount of non useful edges, all edges with a negative x-component was set to zero before thresholding, as the interesting edge is always in the same direction.

6.2 Segmentation results on TM

The algorithm was tested on the set of 20 TM-images. Overall, the results look promising. Fig 6.2 shows eight of the resulting two way cost maps.

Figure 6.2. Resulting two way cost map for eight TM-images. No spurious edges outside

the bone are left. Note the diﬀerent costs for the two possible main paths. Dark means low cost.

Within the set of 20 TM-images, no image caused the algorithm to fail segment-ing the real TM-edges from the rest of the image. This indicates that the method could be used as a robust way of initializing a more precise edge ﬁnder at the left border of the segmented region.

(49)

6.3 Snake algorithm

Edge magnitude

Center of mass Sample

Image object External forces

Internal forces

Figure 6.3. The snake algorithm uses a rubber band-like model with internal and

ex-ternal forces. The inex-ternal forces strives to straighten the snake while the exex-ternal forces work like gravity towards the center of edge magnitude mass along samples normal to the snake. The edge image is sampled along the normal to the model points and the center of mass is calculated.

A snake algorithm was implemented and initialized at the left border of the segmented region. The snake algorithm is a popular edge finding technique, it includes ideas from mechanics to build a model of a flexible band drawn to edges. The snake is built up by a set of points. Every point is affected by both internal and external forces. The internal forces tries to straighten the snake while the external forces work like gravity towards high edge magnitude. There are many different implementations of the snake algorithm in the literature, but here we use a center of mass calculation on edge magnitude sampled normal to each model point as external attractor. Another popular method is to use the image gradient as a force field affecting the model points. The snake points are moved iteratively under the influence of the forces. The position update caused by the internal force is calculated as a fraction of the distance a point has to move to end up on a straight line between its neighbour points. The position update from the external force is a fraction of the distance along the normal to the center of edge magnitude mass. Figure 6.3 shows the snake model and the forces acting on it.

(50)

6.4 Result of snake algorithm 37

6.4 Result of snake algorithm

The snake algorithm was run on the whole set of images. After tuning of the parameters it managed to find the edge of the TM in all images in the test set. The snake was initialized at the left border as shown in fig 6.4. The final result from eight images is shown in fig 6.5.

Figure 6.4. Initialization of snake at left border of segmented region.

(51)

6.5 Discussion

The two way cost map segmentation together with the snake algorithm has given the best edge detection results to far. The parameters of the snake could be tuned so that the correct edge was found in all 20 images, however more test data is needed to ensure its robustness. The result is quite sensitive to the snake parameters. If the snake is made too stiﬀ, the curve tends to cross the edge and move inside the TM at sharp corners. On the other hand, if the snake is too ﬂexible it may fail to move into the maximal edge if the segmented region (and the initialization) reaches out far from the true edge in some spots.

There might be other methods than the snake that are better suited for the ﬁnal detection of the edge. The ASM-method could be a candidate, but there is still the problem of where to place corresponding landmarks in the model.

(52)

Chapter 7

Conclusion and further work

Among the methods tested in this thesis the two way cost map segmentation com-bined with a snake algorithm produced the best results for detecting the TM-border.

An indication of the robustness of a method is given by the number of param-eters that have to be tuned to make it work for the whole input set, as well as the size of the valid parameter interval. Few parameters and a large interval in which they give successful results indicates a robust algorithm. The two way cost map segmentation requires only one parameter, the threshold for the edge image. However, the snake algorithm has lots of parameters that needed quite a lot of tuning before valid results was produced. This indicates that other methods could be a better choice. For example the canny edge detector could be used as a step in the ﬁnal edge detection. It would further reduce the information in the segmented area and leave only local maximum edge pixels, which would limit the complexity of the ﬁnal step.

In the ASM-method the only parameters are the number of eigenmodes in the model and the length of the normal samples. Instead the method requires a lot from the model built by the user in means of good placement of landmark points and size of the training set. Further development of landmark placing techniques has to done before this method can be used.

No experiments have yet been done to ﬁnd the actual landmark along the found edge. Once the correct edge detection has been done, the location of the curvature extremum proposed in chapter 3 seems like a smaller obstacle. One problem that might arise is if the second derivative of the edge curve is rather constant around the maximal curvature. This would make a distinct extremum diﬃcult to distinguish. The position of the curvature extremum will also be dependent on the scale of smoothing used by the second derivative kernels.

(53)

(54)

Bibliography

[1] T. F. Cootes, C. J. Taylor. Statistical Models of appearance for computer vi-sion. Wolfson Image Analysis Unit, Imaging Science and Biomedical Engineering, University of Manchester. Online technical report, 2000.

[2] G¨osta H. Granlund, Hans Knutsson. Signal Processing for Computer Vision. Kluwer Academic Publishers, 1995.

[3] F. Leymarie and M. D. Levine. Curvature Morphology. Technical Report, TR-CIM-88-26, McGill University, Dec. 1988.

[4] Ronald L. Rardin. Optimization in Operations Research. Prentice Hall, 1998. [5] Mikkel B. Stegmann, David D. Gomez. A Brief Introduction to Statistical Shape

Analysis. Informatics and Mathematical Modelling, Technical University of Den-mark. Online technical report.

[6] H. H. Thodberg, A. Rosholm. Application of the Active Shape Model in a commercial medical device for bone densiometry. Online technical report.

(55)

Automatic landmark detection on Trochanter Minor in x-ray images

Automatic landmark detection on

Trochanter minor in x-ray images

Automatic landmark detection on

Trochanter minor in x-ray images

Avdelning, Institution

Division, Department

Institutionen för systemteknik

581 83 LINKÖPING

Datum

Date

2005-02-11

Språk

Language

Rapporttyp

Report category

ISBN

Svenska/Swedish

X Engelska/English

Licentiatavhandling

X Examensarbete

ISRN LITH-ISY-EX--05/3701--SE

C-uppsats

D-uppsats

Serietitel och serienummer

Title of series, numbering

ISSN

Övrig rapport

____

URL för elektronisk version

http://www.ep.liu.se/exjobb/isy/2005/3701/

Titel

Title

Automatisk landmärkesdetektering på Trochanter Minor i röntgenbilder

Automatic landmark detection on Trochanter Minor in x-ray images

Författare

Author

Per Holm

Sammanfattning

Abstract

During pre-operative planning for hip replacement, the choice of prosthesis can be aided by

measurements in x-ray images of the hip. Some measurements can be done automatically but this

require robust and precise image processing algorithms which can detect anatomical features. The

Trochanter minor is an important landmark on the femoral shaft. In this thesis, three di.erent

image processing algorithms are explained and tested for automatic landmark detection on

Trochanter minor. The algorithms handled are Active Shape Models, Shortest path algorithm and

a segmentation technique based on cumulated cost maps. The results indicate that cumulated cost

maps are an e.ective tool for rough segmentation of the Trochanter Minor. A snake algorithm was

then applied which could .nd the edge of the Trochanter minor in all images used in the test. The

edge can be used to locate a curvature extremum which can be used as a landmark point.

Nyckelord

Keyword

Abstract

Acknowledgment

Notation

Symbols

Abbreviations

Contents

Chapter 1

Introduction

1.1

Objectives

1.2

Method

Chapter 2

Image analysis

2.1

What is image analysis?

2.2

Edge detection

2.3

Orientation

2.4

Utilizing ´

a priori knowledge

2.5

Landmark detection

Chapter 3

Problem description

3.1

_{Title of series, numbering}