Finding 3D Teeth Positions by Using 2D Uncalibrated Dental X-ray Images

(1)

MEE 10:90

Finding 3D Teeth Positions by Using 2D Uncalibrated Dental

X-ray Images

Bitra Sridhar Dandey Venkata Prasad

This thesis is presented as part of Degree of

Master of Science in Electrical Engineering with emphasis on Signal Processing

Blekinge Institute of Technology November 2010

Blekinge Institute of Technology School of Engineering (ING) Department of Signal Processing Supervisor: Siamak Khatibi Examiner: Siamak Khatibi

(2)

ii

(3)

iii

ABSTRACT

In Dental Radiology very often several radiographs (uncalibrated in position) are taken from the same person. The radiographs do not provide the depth details, and there is often requirement of three dimensional (3D) data to achieve better diagnosis by radiologist. The purpose of this project is a step forward to solve needs of dentists for evaluating the degree of severity of teeth cavities by 3D reconstruction implementing the uncalibrated radiographs. The 3D information retrieval from two radiographs can be achieved when the 3D position of radiographs are known and when we can find correspondent intensity matching points in the pair of radiographs of the same scene. However the intensity of a point in X-ray images is changed due to variation of proportional distance between X-ray source and detector. In the thesis we propose a novel approach to retrieve 3D information by retendering two radiographs from the same scene. For the retendering, the 3D plane of radiographs is moved to new positions where the intensity distortion due to the distance is minimized. By having the 3D retendered position of radiographs we are able to find accurately relative 3D positions of original radiographs to each other. The proposed methodology approach is tested on different data sets for visually and quantitatively verification and validation.

A data set of simulated radiographs with known 3D information was required to validate our work. Due to lack of such data set a novel methodology is used to simulate the dental X-ray radiographs to validate the 3D positions of teeth obtained from the 2D dental X-ray images by the proposed methodology approach. The simulation of these images is generated by a virtual source of X-rays and optical images of the jaw model which bone and tissue were segmented and used separately. The simulation result is an image, namely the virtual X-ray image. The bone and tissue segmented part of an optical image undergo a process defined by an appropriate response function depending to optical density conversion of absorbed dose to a certain material (bone or tissue). The variables of changing thickness of the material, tube voltage and material of the detector are optional in the simulator. This method of simulating X-ray images is also useful in diagnosis phases for comparison between original and synthesized one.

(4)

iv

List of Figures

1 Failure Maxilla first molar Intraoral . . . 1

2.1 Permanent Jaw A: Upper, B: Lower. . . 3

2.2 Anatomy of tooth. . . 3

2.3 Baby Jaw A: Upper, B: Lower. . . 4

2.4 Maxilla and Mandible of Skull. . . 4

3.1 Positioning of radiograph in teeth structure. . . 6

3.2 Positioning of Radiograph in maxilla. . . 7

3.3 A: Bitewing Radiograph B: Periapical Radiograph . . . 7

4.1 Perspective model . . . 8

4.2 Model of the 2D Projective Space . . . 10

4.3 Epipolar geometry . . . 13

4.4 Camera imaging system vs x-ray imaging system . . . 21

4.5 Radiography system geometry. . . 22

4.6 Vanishing points and Vanishing line of a Hexagonal calibration object. . . . . 23

4.7 Initial position of the calibration grid for vanishing lines method. . . 23

4.8 Geometry of Biplane X-ray Imaging system. . . 24

5.1 Effect of intensity on florescent screen. . . 26

5.2 X-ray imaging . . . 27

5.3 Linear attenuation coefficient v/s X-ray energy . . . 27

5.4 Characteristic Response curve . . . 28

5.5 Blurring due to imaging system. . . 29

5.6 A: Isometric display of point stimulus . . . . . . . 30

5.6 B: PSF rotationally symmetric. . . 30

5.7 Imaging of Clock by lens . . . 30

5.8 Spatial filtering in an optical system. . . 31

5.9 Frequency domain filtering in an optical system . . . 32

6.1 System Model . . . 33

6.2 Flow chat of system model . . . 34

6.3.1 Setup to take X-ray. . . 35

6.3 Periapical X-ray Radiograph of teeth . . . 36

6.4 Computed Epipolar Geometry . . . 36

6.5 3D position recovered images . . . 37

6.6 3D points for two X-ray images . . . 38

6.7 Process of taking X-ray radiograph to teeth . . . 38

6.8 Dental X-ray images of teeth. . . 39

6.9 3D position of centre of x-ray of teeth. . . 40

6.11 Bitewing X-ray radiograph . . . 40

6.12 Computed Epipolar geometry . . . 41

6.13 3D positions recovered Bitewing Radiograph . . . 41

6.14 Full mouth set of Bitewing X-rays . . . 41

6.15 3D positions of centre of x-rays of teeth . . . 41

(7)

vii

6.16 Models of dental lower jaw . . . 42

6.17 Experimental Setup to take optical images . . . 42

6.18 Sequence of images around the model. . . .43

6.19 3D Positions of Centres of Optical Images of teeth . . . 43

6.20 Flow diagram of x-ray simulation . . . .45

6.21 Segmented image of teeth for bone and tissue . . . 46

6.22 Bone and tissue 2D-PSF. . . 46

6.23 Magnitude of 2D-PSF Bone and tissue . . . 47

6.24 Simulated X-ray Radiograph . . . 47

6.25 Pair of simulated X-ray Radiographs. . . 47

6.26 epipolar geometry for simulated X-rays . . . 48

6.27 3D position recovered for simulated X-rays . . . 48

6.28 3D Positions of Simulated Pair of X-rays . . . 48

List of Tables

2.1 Effective dose from diagnostic X-ray examinations. . . 5

3.1 Comparison of different diagnostic system. . . .7

6.1 Data sets used for Experiment. . . 33

6.2 Angle measurements between individual optical images . . . .44

6.3 Angle measurements between individual X-ray radiographs . . . .44

(8)

1

Chapter 1 Introduction

Dental X-rays are pictures of the teeth, bones, and soft tissues around them, helps to find problems with the teeth, mouth, and jaw. X-ray pictures can show cavities, hidden dental structures (such as wisdom teeth), and bone loss that cannot be seen during a visual examination.

They are very useful in detecting the early stages of decay between teeth. These X-rays use a small amount of radiation and have some common forms. The panoramic X-ray is seldom taken, but shows a broad view of the jaws, teeth, sinuses, nasal area, and temporomandibular joints.

These X-rays do not find cavities. These X-rays do show problems such as impacted teeth, bone abnormalities, cysts, solid growths (tumours), infections, and fractures [53].

More common are the Bitewing and the Per-apical X-rays. These methods are important as they are widely used even in the small clinics for finding out the problems mentioned above.

Dentist manipulates the indicator cone behind the teeth where area of diagnosis is required. The indicator cone is operated from outside the position and orientation of the film adjusted inside the mouth to get exact projection. The bitewing shows the upper and lower back teeth and how the teeth touch each other in a single view. These X-rays are used to check for decay between the teeth and how well the upper and lower teeth line up. They also show bone loss when severe gum disease or a dental infection is present. The planes of the detector and the cone are aligned parallel in Bitewing X-rays. This arrangement makes Bitewing X-rays give exact view of the internal structure of the teeth. On the other hand, Intra oral radiography gives the perspective vision of the teeth and so is the difficulty of the doctor to find the individual differences of patients. In [24, 25, 27, 28] report about the problem of raising demand in the requirement of skilled students, technicians and dentists. For obtaining good quality X-ray radiograph. It also reports about the failures of the dental students as shown in the Figure 1(a) and fig 1(b) for a maxilla first molar Intra oral radiograph.

Figure 1: Failure Maxilla first molar Intraoral radiographs Ref. [29].

Most dentists still use X-ray plates to capture the images, but more and more are going over to digital sensors. Less radiation is needed to make an image with digital radiography, than with standard dental X-rays, and images can be recalled from digital storage easily. In general cases more than one X-rays is usually taken, for a full mouth up to 21 radiographs are needed, there is enough replicate information to build a three dimensional image of the mouth. Recently, 3D structure recovery through self-calibration of the camera has been actively researched.

(9)

2

Traditional calibration algorithm requires known 3D coordinates of the control points where as self-calibration only requires the corresponding points of images, thus it has more flexibility in real application. In general, self-calibration algorithm results in the nonlinear optimization problem using constraints from the intrinsic parameters of the camera. Thus, it requires initial value for the nonlinear minimization. Traditional approaches get the initial values assuming they have the same intrinsic parameters while they are dealing with the situation where the intrinsic parameters of the camera may change. In this paper, we propose new initialization method using the minimum 2 images.

The proposed 3D reconstruction method is validated on simulated and real data and is shown to perform robustly and accurately in the presence of noise. A new methodology for simulation of radiographs is introduced which is one of the two important folds of this project.

Once a three dimensional model of the entire mouth, or area of interest in the mouth, has been built up, the model can be shown on an auto stereoscopic display. An auto stereoscopic display produces depth perception in the viewer even though the image is displayed on a flat device.

Such a display would allow the dentist to visualize the entire mouth in three dimensions, without the need for glasses or other equipment. The extra dimension in the image means that the model on the display matches very closely to the spatial reality of the patient‟s mouth. This close fit with reality makes misdiagnosis less likely, and gives the dentist freedom to explore the hidden areas of the mouth as easily as the surface of the teeth.

1.1 Problem Statement and Thesis outline

 Given a set of radiographs of the dental structure - is it possible to transform 2D x-ray images to 3D positions using projective model of x-ray system?

 Simulation of X-Ray images - is it possible to get the optical transfer function of the dental x-ray system?

In this thesis, the proposed approach is uses the epipolar constraint that relates corresponding point in both X-ray images to directly estimate the rectifying transformations. Some methods have been proposed in uncalibrated rectification of the image. However, the constraints that the thesis is structured as follows: This chapter printed a brief introduction to the problem. The main purpose of this chapter is to put forth some questions that will be answered in the thesis. Chapter 2 explains the structure of the tooth and its anatomy, importance of the role of radiographs in dental disease detection and effective radiation in dental imaging. Chapter 3 shows the importance of bitewing and per apical imaging, and its positioning. Chapter 4 presents a brief, basic concept of projective space transformation, followed by a description of epipolar geometry and it is been concluded by methodologies for X-ray calibration. Chapter 5 deals with the simulation of X-ray Radiographs from optical images. Chapter 6 concludes with the experimental results and discussion of the work done.

(10)

3

Chapter 2 Anatomy of Dental teeth

2.1 Structure of teeth

Teeth are a hard bonelike structure in the jaws of humans which are used for biting and chewing. An adult human normally has 32 teeth, 16 in each upper jaw and lower jaw has on each side two incisors, one canine, two premolars and, at the back three molars as shown in Figure 2.1. Projecting from the gum is the crown of the tooth in Figure 2.2. The part embedded in the gum and reaching into a socket in the jawbone is known as the root. The body of the tooth is made up of a hard pale yellowish bonelike substance called dentine. Inside this there is pulp cavity, which contains blood vessels and nerves.

Figure 2.1: Permanent Jaw A: Upper, B: Lower.

Covering the crown of the tooth is a layer of enamel of varying thickness. It is made up almost entirely of apatite crystals with calcium phosphate filling. Calcium phosphate is also mainly responsible for the hardness of dentine. The enamel layer crystals are elongated and are all

arranged with their ends toward the surface of the enamel.

Around the root of the tooth, enamel is replaced by cementum (the cement), another bonelike material which fixes the tooth firmly in the socket of the jaw. However, between the bone of the jaw and the cementum layer there is a layer of tissue, called the periodontal membrane, which is in contact with the tissues of the gums and the pulp cavity. Incisor and canine teeth have a single root, premolars have a double root and molars have three branches to the root.

Figure 2.2: Anatomy of tooth.

(11)

4

A child up to about the age of six has only 20 teeth. These are the milk teeth or deciduous teeth which are gradually replaced by the permanent teeth after the age of about six years. There are no molars in the milk set, but the teeth corresponding to the premolars are known as the milk-molars and they perform the same grinding function. The first milk teeth to erupt break through the gums are usually the two central incisors in the lower jaw. These normally appear between eight and ten months after birth. Then, in a fairly regular order which takes 18 months to 2 years, a further 18 teeth push their way through the gums. By the time, a child is 2 or 2½ years old it should have a full set of 20 teeth as shown in Figure 2.3. Although quite small, the deciduous teeth are very strong and well able to deal with the hardest foods that the child may come across.

Figure 2.3: Baby Jaw A: Upper, B:Lower

The teeth are formed with the substance of the jawbones, the maxilla above and the mandible below in Figure 2.4.

Figure 2.4: Maxilla and Mandible of Skull.

2.2 Role of radiographs in dental disease detection

A considerable amount of the information needed to assess periodontal diseases that can be obtained from the clinical examination. Gingivitis is detected by noting characteristic soft tissue changes, and many key features of periodontitis are also exclusively detectable clinically such as bleeding on probing, probing pocket depths, and loss of attachment, mobility and suppuration.

Much of the information is required to diagnosis and generates a treatment plan to the clinician by the end of the clinical examination. Some features relevant to the periodontal diseases may be detectable both by clinical and using radiographic techniques. Such features include supragingival and sub gingival calculus, other plaque retention factors and furcating defects.

(12)

5

However, neither method of investigation will detect every occurrence. Additionally radiographs provide information about the bone levels and pattern of bone loss that cannot be gained through routine clinical examination. Periapical pathology, including combined lesions of both endodontic and periodontal origin, may be suspected clinically but can only be visualised radio graphically on appropriate radiographs [31].

2.3 Effective radiation in dental imaging

The use of x rays for diagnostic purposes in dentistry encounters some risk with the exposure to the patient. The possible risk that such exposure entails and the methods used to affect exposure and reduce dose. This information in Table 2.1 provides the necessary background for explaining to concerned patients the benefits and possible hazards involved with the use of x-ray examinations, Dental x-ray are aimed in a tight beam at a small spot on the face. The only structures that receive the full dose of x-radiation are the tissues in the direct line of fire. The rest of the body receives only the radiation that is scattered off of the structures in the line of fire. Much less radiation scatters from an object in an x-ray beam than from an object in a beam of ordinary light due to the difference in the nature of the respective radiation sources.

Furthermore, the tissues in which dental x-rays are aimed are much less prone to injury from x- radiation than are tissues in other parts of the body, such as the intestinal lining or reproductive organs and other constantly reproducing tissues. The newest unit of measurement, the milisievert was designed to take this factor into account [55].

Examination Resolution Effective dose(uSv)

Full mouth periapical High 171

Panaromic Moderate 26

Cone Beam Computed Tomography

High 599

Computed Tomography High 860

Table 2.1: Effective dose from diagnostic X-ray examinations

(13)

6

Chapter 3 Dental imaging

3.1 Importance of Bitewing and Periapical imaging

In a bitewing imaging, all three elements, the teeth, the detector, and the x-ray beam are optimized to give the most undistorted shadows possible. The film and teeth are parallel, and the beam is aimed directly at 90 degree angle. Thus bitewing detector affords the most accurate representation of the true shape of the teeth and associated structures such as decay, fillings, shape of nerves and bone levels.

A periapical imaging is shot from an angle in which the three elements are not necessarily aligned parallel. Some distortion is introduced on purpose to be sure that the shadow of the entire tooth or teeth falls on the detector. This is done because in many instances, the space available in the mouth or the curvature of the roof of the mouth will not permit parallel placement of the detector, as shown in Figure 3.1.

3.2 Proper positioning of radiograph

Positioning of the Patient to radiograph, the maxillary arc and the patient‟s head should be positioned upright with the sagittal plane vertical and the occlusal plane horizontal. When the mandibular teeth are to be radiographed, the head is tilted back slightly to compensate for the changed occlusal plane when the mouth is opened. Detector Placement and the projections described for the paralleling technique may also be used for the bisecting-angle technique.

Figure 3.1: Positioning of radiograph in teeth structure.

Molars have an advantage of the entire length of the teeth can be projected on the detector. The detector resting on the palate with its midline centred at the midline of the arch, along the parallel axis of the maxillary central incisors.

(14)

7

Figure 3.2 Positioning of Radiograph in maxilla.

Projection of the central ray passes through the contact point of the central incisors and perpendicular to the plane of the detector and roots of the teeth. Because the axial inclination of the maxillary incisors is about 15 to 20 degrees, the vertical angulations of the tube should be at the same positive angle. The tube should have „0⁰‟ horizontal angulations, as shown in Figure 3.2. Like this to all the teeth detectors are positioned to take radiographs [55].

A) B)

Figure 3.3: A) Bitewing Radiographs B) Periapical Radiographs Comparison of MRI and CT scan methods with Intra Oral diagnosis:

MRI & CT SCAN Intra Oral Radiography

High Cost around 3000sek Low Cost – 300sek

High Radiation 15 – 70mREM Low Radiation 1.5mREM

3D data – More information 2D data -- Less information Table 3.1: Comparison of different Diagnostic Systems

(15)

8

Chapter 4 Introduction to optical geometry for 3D vision

4.1 Basic concepts

“Problem of vision” is yet to be solved. A camera transforms a 3D scene point into an image point. Here the process is done by transformation of spaces. There are different geometric spaces which have an associated set of transformations with specific invariants. The behaviour of real cameras is ideally modelled by a perspective projection (i.e. pin hole) but how good this model has two aspects, which are extrinsic and intrinsic. In general a method which does not require knowledge of principle point is based on a fundamental result in projective geometry.

Euclidean space comes under the group of similarity transforms where the rigid displacements scale is uniformly changed. This space cannot tell you about variation of real world and reduced model images. But in this we can calculate angles and rotations. In this chapter, an introduction to homogenous coordinates system, the simple modelling of camera, the relation between two views and rectification of pairs of images are explained.

4.1.1 Phenomena of “at infinity”

The extension of finite geometry to handle “at infinity” phenomenon was initiated by painters. They used this to produce persuasive 3D depth. During their work, considerable use of vanishing points is done and many geometrical constructions are derived. The examples are parallel lines at 3D space try to converge in the images and stars are fixed as gone walk along lines etc.

4.1.2 Perspective model

This model is a pinhole camera where the image formed is by a simple optical device through a small hole. Hence pinhole is nothing but the centre of projection (C_p). The line passing through C_p and intersecting image plane forms a pixel in the image. This model can be formulated by an equation for mapping a point in space into an image and the Cp is on back of image plane. Finally this ends up with matrix multiplications where the matrices are in Euclidean coordinates. To make it even simple, homogenous coordinates are developed. So, theory about this is first presented and then to projective space.

Figure 4.1 Perspective model

(16)

9

4.1.3 Homogenous coordinates

An incoming light ray passes through many 3Dpoints and projects to an image point. So here for any 3D point on the ray has some image point. Suppose, (0, 0, 0) is origin of the camera and (X, Y, T) is a 3D point. Now,

X, Y, T) = ( X, Y, (4.1) Also lies in the same ray.

(X, Y, T) (X, Y, T)=( X, Y, T) (4.2) This can be extended to any dimensions say (X, Y, Z, T) is 3D point in Cartesian coordinates and so,

XC = [X, Y, Z] ^T is mapped into homogenous form as

Xh= [wx wy wz w]^T

Here Xh is more than one and at any time it is one point for any value of „W‟. Mapping backwards can be done as,

(4.3) Though transformation between homogenous to Cartesian are irreversible, the space they represent are projective and Euclidean space respectively.

4.1.4 Line representation and duality

The homogeneous representations of point in homogeneous coordinates have a good connotation which links points and lines. An equation of the form, Ax + by + C = 0, in Cartesian coordinates represented a 2D line. Where as in homogeneous coordinates it is represented as Ax + By + Cz = 0. So, now [X, Y, Z]^T is similar to [A, B, C]^Twhich are representing a point and line respectively. This concludes that, points are lines represented as triplets may be inter changeable.

4.1.5 Ideal Points

The algebra of homogeneous coordinates will help to calculate the intersection of parallel lines, planes and hyper planes. Unlike in the Cartesian coordinates, the inter section between two parallel lines can be solved as, [see Figure 4.2]

A1x+B1y+C1z=0 (4.4) A₂x+B₂y+C₂z=0 (4.5)

(17)

10

Subtracting the two equations after dividing the first equation and second equation by B₁ and B₂ respectively and then applying the condition , we get

(C₂-C₁)Z = 0 (4.6) Now C₁ C_2.Since two different lines are used in the process. Finally „Z‟ is equal to zero.

Therefore in homogenous form, the inter section of two parallel lines is defined as Xh= [x y 0] T

(4.7) and for 3D

Xh= [x y z 0]^T

Thus, the points which have last co-ordinate with zero are called the points at infinity. These points are called ideal points.

Xw

Z Y

Xc

X

Figure 4.2: Model of the 2D Projective Space.

4.1.6 Projective space transformations

Algebraic representation of the transformations is the relevant aspect of homogenous coordinates. Rigid transformations are the one which deal with Cartesian coordinates and values of the angles changes. Algebraic representation involves matrix multiplication and addition which contain scale and translation.

As an example a 2D point transform from x1 to x2 is given as.

(4.8) Where „ ‟ is the rotation angle, S= [sx sy]^T for scale and T= [tx ty]^T for translation axis.

Projective Plane

Euclidean Plane Z=1

Ideal Plane Z=0

0

(18)

11 Its general form is

where „R‟ is an orthogonal matrix.

Affine transformation can take more general type of form by replacing „R‟ with „A‟ and also not orthogonal conditionally, as

(4.9) A more detail of affine transformation is explained in further topics. But the main point to be considered here is this transformation doesn‟t preserve the value of angles, but only parallel lines. Also, Euclidean geometry defines similarities related to principles and theorems and affine geometry explains affine transformation principles and theorems.

Form the above discussion we can learn a flow, which states that affine transformation has a special case of similarity transformation and affine is a special case of homography. Here cosmographies are the transformations in projective space. It preserves collinearly and cross ratios, which are defined by homogenous coordinates. It‟s generalised form as

X₂=HX1 (4.10) where H is a 3X3 matrix.

So, homographies are also included in rigid and affine transformation. It can be shown in general

form for rigid transformation as, X2= X1 (4.11)

Where

and for affine transformation as,

X2= X1

At Z=1 plane is explained by zero in last row which according [2] is a Euclidean plane and constraint to Euclidean points.

4.1.7 The affine group

An affine transformation preserves the points at infinity, and can be broken down as a non-singular linear transformation followed by a translation. A point x on the affine plane is transformed as

where a and t are the linear transformation and translation matrices defined as

(19)

12

(4.12)

The affine transformation therefore can be represented by a non-singular matrix

(4.13) In the projective space defined previously, the line at infinity for the 2-D space contains points with final co-ordinate zero. The affine space can be considered to be embedded in the projective space, corresponding to a map of all points with nonzero final co-ordinate. For every point that is not on the line at infinity in homogenous representation [x, y, z] , its affine coordinates may be computed as

(4.14) This group of transformations preserves parallel lines, the ratio of lengths of parallel line segments and the ratio of areas [4].

4.1.8 Perspective camera model analysis:

In the perspective camera model the projective space algebra are used to define space point and map to an image plane. From the previous discussion, a 3D point to a plane is transformed by special homography matrices which are homogenous coordinates known as projection. The Figure 4.1 used for introducing perspective projection is again considered here for further details. Let (x, y) be the 2D coordinates of P and (X, Y, Z) the 3D coordinates of P. A direct application of Thales thermo [4] shows that:

(4.15) Assume f=1, or other value which relates for scaling of the image. Now, complete calibration of the model is introduced which homogenous coordinates as

(4.16)

In Figure 4.1 origin of real image is not the principle point and different scaling of axis of image plane is required. A matrix „K‟ need to be defined for transformation and 3D coordinates miscoincidence must be rectified by multiply a matrix „M‟ which is a motion in Euclidian space.

(4.17)

(20)

13

where M is [R_3X3 T_1X3] with 6 degrees of freedom of a rigid motion, R defines a rotation matrix and T a translation. Finally the independent camera parameters of camera position and upper triangular matrix are K.

K = (4.17)

Sx and Sy are stands for axis scale factors, is the non orthogonality between the axis ( ) and the coordinates of principal point. An important point to consider is, camera model which is defined by 11 camera parameters: as angles parameters in the rotation matrix, tx, ty, tz in translation matrix and in K matrix. So, are intrinsic parameters and are extrinsic parameters. Finally explicit parameters representation of projection matrix can be written as,

P =

(4.18)

4.1.9 Projective stereoscopic viewing

Stereoscopic viewing in projective space will help to move in the way to make it happen with calibration. This means that it may work with completely uncalibrated cameras. This will reduce our work to calibrate the cameras where it is unable to do.

This may lead to provision of enough mathematical back ground for any type stereovision and multiple image algebra, like self calibration techniques. Next section of this chapter will introduce brief information about self calibration techniques and its basics.

4.2 Epipolar geometry

Figure 4.3: Epipolar geometry

The Figure 4.3 explains the epipolar geometry. The point „M‟ is imaged in two views at the positions „X‟ and X^‟. These points are in the image planes formed by the cameras with camera matrices O and O

´

. As „M‟ varies, the set of all planes defined by the centres of the cameras (O and O^‟) and „M‟ forms a family of planes with baseline OO

´

along ee

´.

The correspondence between two points is defined by epipolar plane XOO

´

. The epipolar plane I

´

will correspond to a point X, which is the intersection of epipolar plane and the other image

(21)

14

plane. This is symmetrically same for point X

´

. The expression for the back projected ray for the first camera P form x and containing the center of first camera O is given as,

M ( ) = P⁺X + O (4.19) where P⁺ is pseudo inverse of „P‟.

The epipolar line forming is

l

´

= (P

´

O) x (P

´

P⁺x)

The epipolar is P

´

O as shown in the Figure (4.3) and it is defined by the skew symmetric matrix by the cross product. So the epipolar line in the second image of the point in the first image is given by,

l

´

= [e

´

]_x(P

´

P⁺x) = F x (4.20) Here,

F = [e

´

] x (P^lP⁺) (4.21) Where [e

´

]x skew symmetric matrix which introduces the rank 2 constraint on fundamental matrix.

The relation between the 2D family of point and 1D family of lines is formed by a map induced by F. So, finally

x

´

(F x) = 0 (4.22) This equation explains that x

´

is required to prove the equation for an epipolar line of x in the image I

´

. Equation (4.22) is nothing but epipolar constraint, which is satisfied by the corresponding points in two images.

If [t]_x represents the cross product corresponding to translation vector, then the coplanar principle can be satisfied by

x

´

([t]x Rx) = 0, (4.23) Where the essential matrix is

E = [t]x R. (4.24) This expression will take the x

´

to the image in the first plane. For uncalibrated case K and K

´

are used which are calibration matrices to transform the image points in pixels with point in camera coordinate system. Comparing equations (4.23) and (4.24) will give,

(K^-l x

´

)^T E (K^-l x) = 0 (4.25)

(22)

15

Finally, the relation between essential and fundamental matrix is given as

F = K^-TEK-l

(4.26) The intrinsic parameters and extrinsic parameters of the camera are both described by the fundamental matrix, which is advantageous. This reveals that it only depends on matched points in both the images. If there is a pure translation which is horizontal parallel, then there is no change in internal parameters occurs and positions are differed by only translations.

It is given by

P=K [ P^l= K [I|-T] (4.27) The equation (4.21) can be written as

F = [e

´

] x K K-l

(4.28) (4.29) The equation (4.29) is fundamental matrix for pure translation and in matrix form it becomes.

(4.30)

When we apply projective transform, the epipole is moved to a point at infinity horizontally this process of estimation is called “projective rectification”.

After transformation the image points are given are given by,

= H x and (4.31) The relation between , and F becomes

(4.32) or

( ) = 0 (4.34) Finally

F = H

´

H (4.35) The estimation of fundamental matrix F is needed to get the epipolar geometry. To find F, initially we need point correspondences between images. If m

´

is a point in the second image with respect to the point m in the second image, then it must lie in epipolar line l

´

=FM and so m

´

^t^{. l}

´

=0 [12, eq. 5.1]. Then the constraint on epipolar geometry is therefore,

m

´

^tF m = 0 (4.36)

(23)

16

The equation (4.36) shows that F is estimated from 8 point matches in two images of same scene. F has seven degrees of freedom as defined before. The constraint det (F) = 0 is implied by rank 2 condition of 9 linear conditions modulo of matrix F (algebraically). So, F is calculated by 7 matching points and a rank constraint and three possible solutions.

In [32] this procedure of automatic epipolar geometry algorithm is explained for good out come and with some care. For 8 matched points (minimum), it may have noise. Then the least squares approach gives good solution as:

(4.37) The equation (4.37) calculates the eigen vector and this vector is for the smallest Eigen value in symmetric positive semi definite normal matrix A^tA. All these numerical techniques are explained in [44]. The rank constraint is not imposed on it and so another computation needed to be added to project the F on to rank 2 subspaces. This is done by SVD of F, SVD decomposes F into

F=QDR (4.38) where D is diagonal and Q, R is orthogonal. When D 0 the result is the desired one.

The method is unstable for application on images natively. Stability can be achieved when the pixels are normalized from [0, 512] to [-1, 1] prior to calculation. This provides a stable and accurate solution of F. But there is also a need of rejecting false correspondences in the input data as in [16].

Summary:

The estimation epipolar geometry is:

1. Extract points from images

2. Find out point correspondences as in initial set

3. Estimate the fundamental matrix according to above algorithm 4. Least median square approach is used to estimate the results 5. And make it more robust.

(24)

17

4.2.1 RANdom Sampling Consensus (RANSAC)

RANSAC estimates the fundamental matrix form a set of point matches in two images of same scene. The problem is that the matches obtained are not accurate. Point matches that satisfy the estimated epipolar constraints are inliers and without are outliers. Outliers are resulted from the inaccurate matching of points from the set of images. The estimation of error function eliminates outliers by minimization process using threshold, which maximizes the number of potential inliers.

4.2.2 RANSAC Algorithm for non-linear minimisation

The epipolar geometry of two images gives a set of matched points and the correspondences between the images gives 8 point matches. The non-linear iterative version is preferred for its robustness, and in combination with the RANSAC method, can be used to eliminate potentially false matches. The idea is to randomly choose n samples from a set of N point matches and find the corresponding fundamental matrix F by iterative minimization of a distance function (4.39) or (4.40).

(4.39)

(4.40) A threshold is applied to all of the N point matches to eliminate possible outliers using error function. This process is repeated to get the best fundamental matrix with the most number of inliers. Probability of selected set of point matches p of an outlier to all of n selected matches as inliers is given by . The minimum samples to be tried for a set of selected point matches consisting of inliers with probability P is given as,

(4.41) An important factor in using RANSAC is to estimate the fundamental matrix with the size of the sample set used in the iterations as in [5] [13] [3]. Once most of inliers are found, the algorithm can be terminated. The difference between inliers and outliers are differentiated with threshold epipolar distance as (4.42) and (4.43).

(4.42)

(4.43)

(25)

18

Computing F from (4.39) or (4.40) Using (4.43), determine the number of inliers corresponding to this F in the N - set of points. If Z values for the fundamental matrix are obtained, the set of fundamental matrices obtained, by choose the one that minimizes (4.39) or (4.40) with the most number of inliers. An adaptation of this method of estimating the epipolar geometry has been used to develop a rectification scheme that minimizes distortion in the images.

4.3 Rectifying transforms for uncalibrated

To make the epipolar geometry simple, the steps we go through is a modified form of epipolar constraint. As shown in (4.35). But to make the reconstruction look more close to real image, it is required much that the constraints are imposed which relates to physical properties of the image. The process starts directly with pre-computed matching points to rectify the pair of images of same view. Disparate positions in space related by arbitrary rotation and translation are constrained to lie on skewed lines in the image. These constraints on simplified geometry after rectification process are still needed to resemble the rectified images as original images.

The epipolar constraint is given as

X

´

^TF X = 0 (4.44) and then

X ^T F ^T X ´= 0 (4.45) This constraint reduces epipolar geometry to a canonical form and represents pure translation.

The fundamental matrix for rectifying transforms are given as

F = H

´

^T H (4.46) The transformations which satisfy (4.46) are not unique and so exact choice of the transform is not made. This may lead noticeable distortion in the rectified images.

4.3.1 Approach to uncalibrated rectification

The uncalibrated rectification is done in three folds. First, a parameterization of the transforms, second factorisation of fundamental matrix and followed by defining a cost function to estimate epipolar geometry for a detailed geometric condition and third constraints on the transforms at the same time to reduce visual distortion.

The factorisation and parameterisation is explained simultaneously, below as in [15] the transformation „H‟ in projective method can be represented as product of two matrices,

H = (H^-l)^-l= (QhRh)^-l= RHQH (4.47) This is true as Q is an upper triangular matrix and R is an orthogonal rotational matrix.

(26)

19

R is represented about x-, y- and z- axis and so R=R_xR_yR_z and R_z is orthogonal matrix so finally the homography matrix is represented as [17],

(4.48)

Replace in (4.48) with R_H for first matrix, Q_H for second matrix, R_y is third matrix and R_X for fourth matrix.

This look like segmenting the homography matrix which is same as

H = Hs Ha Hp (4.49) Here s - similarity, a - affine and p - projective component of H.

Finally the H is described with parameters as

(4.50) Where „ ‟ is rotation in image plane. It can be chosen to make the second coordinate of the epipole to zero for epipolar lines parallel to x-axis. Scaling and skewing are represented by QH.

Projective rectification for unrectified images by (4.44) will give fundamental matrix with independent of first rows. This is due to criteria of and so the estimation results from (2) as, [5]

(4.51) When we attribute only of the tilting and panning to one of image and the later to other image of the same seen, there is chance in reduction of degrees of freedom and so, the (4.49) reduces to

(4.51) Finally can be used to build , which is given by (4.51)

(4.52)

4.3.2 Transformation constraints

Considering which describes the fundamental matrices in terms of rectifying homographies reveal that constraints on the fundamental matrix estimation are also applied to estimate the rectifying transforms. Also, there are conditions to be considered to find the rectifying homographies minimized with constraints justified has physical significance. Firstly, the geometric distortion is reduced by using (4.42) and (4.43) extended to include epipolar lines to get symmetric distance criterion as [16] in terms of (6.10)

(4.53)

(27)

20

whenever there is uncertainty of matched points, the total cost (4.53) is inversely proportional to its variance [16]. Finally minimization over F to minimize the cost as, [16]

(4.54) Second one is visual distortion reduction by orthogonality and aspect ratio. The constraints also should consider the inclusion of orthogonality and relative sizes of the image by looking over it as in [17][18].

The orthogonality in the image is the perpendicular condition between the epipolar line and to pixel columns. To reduce this, error pixel information of the rectified image is used as in [17] and also aspect ratios are used to obtain inimitable cost function. The minimization problem for the estimation of rectifying transformations are given as in [5]

(4.55) Here (4.56)

By ortohgonality principle of

Also, x, y and xH, yH are axis defined as the difference vectors in the affine image plane of original and rectified images respectively.

4.3.3 Affine correction:

The estimation of transformations is more related to projectivities and there is a chance for skewing error for horizontally non uniform scaling. Therefore orthogonality and aspect ratio is re written with Ha included. As [9],

(Ha xH)^T(Ha yH) = 0 (4.57) and (4.58) where H_a isparameterized affine transformation given by,

(4.59)

And as in [5]

(28)

21

4.4 Methodologies for X-ray Calibration

To estimate the parameter for transforming a 3D space point to a 2D image point via a camera, the important requirement is to calibrate the camera. Once it is done we can transform between world and image plane.

4.4.1 Camera imaging

The Figure 4.4 explains the analogy between x ray imaging system and camera imaging system. The only difference is relative position of object, image plane and lens or image intensifier.

Figure 4.4: Camera imaging system and x-ray imaging system

4.4.2 Existing calibration methodologies of x ray source

X-ray calibration setup designed for NDT applications which is already existed is considered. The setup is a cone beam x-ray radiography system which has an x-ray source manipulator of the inspected specimen and x-ray detector. The system setup is analogous to our system but only the inspected specimen is replaced by patient and the detector used here is for NDT applications which have 200-400 Kev X-ray. The detector in medical applications is replaced with 20-180 ev image intensifiers. The schematic of the setup is shown in Figure (4.4).

In this method the X-rays are detected by the image intensifier which is amplified and depicted on a phosphor screen image intensifier which converts the x-rays projection to a visible radioscopic image.

(29)

22

4.4.3 Geometric model

The object coordinate system is similar to the one explained to the system in previous chapters. The notation is used as in [28] to differentiate between the projective geometry. In this method the origin O is the centre of rotation of the object.

There are two methods for x ray calibration of which reviewed first is the calibration of tomography and radiography systems and second is about calibration of 3D vision system. The tomography system calibrations will classify its geometric properties into two parts as intrinsic and extrinsic. In general the intrinsic parameters are constant and related to x-ray detection system geometry. For the later, there is a need of regular estimation as it relates to variations of the object under inspection. Alignment and nonlinear minimisation of projection matrices are used as [19] [20]. High precision, fast estimation and special calibration object are requirements in these methods. The camera calibration methods follow a two stage approach by geometric constraints. Precision is good for a compromise on speed of measurement. This is shown in [23].In [9] estimation by nonlinear minimization for only the extrinsic parameters and usage of [22] method to state with good initial estimate and fast convergence to solution is applied.

The geometry of the radiography system is shown in the Figure 4.5 which includes all geometric configurations explained. It has source, detector and object.

Figure 4.5: Radiography system geometry

The point M on the object is (x, y, z) in Figure 4.5. M is related to x-ray coordinate system by (4.61) where , , are coordinates of F in object coordinates . Also the three successive rotations of the object coordinate systems are R (δ), R(ξ) and R(ψ).

R= R (δ), R (ξ), R (ψ) (4.62)

(30)

23

4.4.4 Estimation of intrinsic and extrinsic parameters.

The parameters ( ) are intrinsic parameters. P and q are calculated by the detector screen plane in radiography system and as in [21] where image plane and detector plane is not same [tomography systems]. The extrinsic parameters are calculated by nonlinear minimization method where cost function is

(4.63) where i is a point in N calibration points.

Pixels coordinates are given by

(4.64) (4.65)

In [5], a detailed explanation and its performance with simulations are presented. The minimization is started with a initial guess by vanishing lines methods as in [4] as shown in figure. It uses the grid of sphere whose lines of joining are not parallel in image plane by vanishing line perspective geometry.

Figure 4.6: Vanishing points and Vanishing line of a Hexagonal calibration object [9]

Figure 4.7: Initial position of the calibration grid for vanishing lines method [9]

(31)

24

Means of vanishing point of all grid points which are center of spheres of the calibration object are calculated. The vanishing lines are calculated by least square method with respect to detector coordinates system „Gd‟, the method uses Euler given by . So, it becomes analogous as world to camera coordinates transformation with object to x ray source coordinate system described by (1)

4.5 Self Calibration of Biplane X-ray imaging system:

To show a good visualisation of structure of interest in mouth of a patient, the direction of camera is specified by two angles as right-left and angulations [30]. The right – left angulations are chosen by the clinician according to the morphology of total jaw. The different coordinate system is considered as shown in the Figure 4.7.

Figure 4.8: Geometry of Biplane X-ray Imaging system

The positions of the X-ray system at S1, S2 and R1 = (S1, x1, y1, z1), R2 = (S2, x2, y2, z2) are the positions of the two X-ray system. The linear transformations are applied for simplicity to model the perspective view of 3D object point on left and right image.

Now the system calibration matrices are calculated as in [30].

= M_i for i = 1, 2 (4.67)

where w is a scaling factor ; u, v are image coordinates and (X, Y, Z) is the 3D point.

Explicit representation geometrically,

M_i = P(d_i , s_p, u_si, v_si ) D_i for i = 1,2 (4.68)

(32)

25

If D_i is the displacement from object to source, s_p known pixel size, linear transformation in homogenous coordinates of perspective image of 3D point and R, Ri are object and source coordinates respectively. Now the matrix representation is given as [39],

=

(4.69)

If the movement of cone-beam is adjusted by the doctor, then angulations are α along x- axis, β rotation along y-axis and translation along z-axis which are measurable here. Then finally is given by,

= ( ) (4.70) Equation (4.68) and (4.70) will give,

(4.71) The mean square error of the system between the analytical coordinates in the object and the matched points calculated by the doctor from two images of same teeth is minimized. The matched points are found by RANSAC algorithm. The Mathematical expression is given by [30],

The equation (4.72) d [ ] gives Euclidean distance. Also, u,v are projection coordinates obtained from (4.67) and (4.71); n is n^th reference point of P_n; i is i^th plane. The iterations are continued with correction of at every time and this will continue until a predefined minimum error is achieved. The total 3D reconstruction process is clearly framed in [30] which include robust estimation and accuracy even in presence of noise.

Summary: In this Chapter the calibration methodologies of x-ray source is introduced. Also, self calibration of the X-ray imaging system is explained. This methodology reveals the fact of its complexity and uses linear transformation for simplicity. In chapter (6) in contrast to these methodologies, a novel simple and popular methodology is explained which re-renders bi-plane images and transformation or X-ray recovery matrix (XRM) is found to calculate 3D information like angle position of centre of images for 3D reconstruction.

(33)

26

Chapter 5 Simulation of X-ray Radiographs from optical images 5.1 Introduction of X-ray radiography

The real image is the one which is real, physical and it may be photography or a radiograph of X-ray. These images are helpful for scientific measurement.

X-ray when targeted to a sensitive material like fluorescent material will emit light which are of varying light intensities. These are similar to television screen. The schematic of this method is shown in this Figure (5.1).

Emission of intensity high for high intensity incident radiation

X-ray

Low emitted light intensity.

Florescent layer

Figure 5.1: Effect of intensity on florescent screen

The fluorescent layer doesn‟t need any light source to see. We can see it using naked eye.

5.2 Image Characteristics

The image formed by fluorescent screen has following characteristics.

1. Noise 2. Contrast 3. Sharpness 4. Resolution

5.2.1 Contrast

In order to identify the structure in a radiograph it must be of different optical density when compared to adjacent ones. Contrast is feature of the images which is the variation of intensity distribution of energy, phase, relaxation characteristics etc, the total outline is the contrast is the nothing but finding the variations in density or luminance if the differences are great we can say the contrast is high and so the vice versa.

(34)

27

To explain this contrast mathematically, let us take an example of situation of projection radiography. The area is exposed to uniform beam of x rays photons of (no). The area exposed is a slab of tissue with linear attenuation coefficient „u‟ and thickness x.

There is also variation towards the right which is small increases of „Z‟. The illustration is shown dogmatically in the figure (5.2).

Figure 5.2: X-ray imaging Ref. [10]

The increased thickness areas and remaining is exited by x ray fluency (photon/m²) of amount

„B‟ and „A‟ respectively. Now the contrast is measured as Now with the thickness

„x‟ we have attenuation which is equal to for

(5.1) (5.2) So, C (5.3) which explains that the variation in the surface cause contrast and increases overall contrast.

Also, varying both and z contrast varies proportionally. At lower energies changes „ ‟ increases by factor of „3‟ and so „C‟ also increases by „3‟. This low x ray energy makes mammography with rich soft tissue contrast. The relation between the subject contrast and x-ray energy is shown in the Figure (5.3).

Figure 5.3: Linear attenuation coefficient increases as the X-ray energy decreases. Ref. [10]

Z x

No No

(35)

28

5.2.2 Detector contrast:

For small variation in the input energy there will be considerable big changes in signal or contrast of image. So the detectors characteristics play an important role in producing the final contrast. A graph between the detected energy and output signal will till about detector contrast

„C_d‟. The curve explains the characteristics of it there are mainly two of them which are linear and logarithmic this shown in Figure (5.4).

Figure 5.4: Characteristic curve A: logarithmic response in radiographs. B: linear in flat imaging system. Ref. [10]

Most of the systems have linear characteristics; the slope(r) measurement of the curve will explain about the characteristics and is given as.

(5.4) Here A and B are adjacent regions. For radiographic contrast is given by RC=ODA-ODB, since the analog film is the output.

5.2.3 Contrast to noise ratio (CNR)

The CNR is a meaning full and frequently used measurement to justify the quality of the image. For digital image acquired by detector systems, a common processing technique is subtraction of a constant number „K‟.

Suppose let us take there are two regions which have average values of A and B. So, after processing it is

(A-K)(B-K) (5.5) As per the equation (5.5),

(5.6)

Finding 3D Teeth Positions by Using 2D Uncalibrated Dental X-ray Images