Book
Book
Medical Robotics
Achim Schweikard
Institute for Robotics and Cognitive Systems
Universität zu Lübeck
Ratzeburger Allee 160
23538 Lübeck, Germany
schweikard@rob.uni-luebeck.de
Phone: +49 – (0)451 – 5005201
Fax: +49 – (0)451 – 5005202
Floris Ernst
Institute for Robotics and Cognitive Systems
Universität zu Lübeck
Ratzeburger Allee 160
23538 Lübeck, Germany
ernst@rob.uni-luebeck.de
Phone: +49 – (0)451 – 5005208
Fax: +49 – (0)451 – 5005202
Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1 Robots for Navigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Movement Replication . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.3 Robots for Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4 Rehabilitation and Prosthetics . . . . . . . . . . . . . . . . . . . . . . 24
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3 Robot Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.1 Three-joint robot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.2 Six-joint robot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.3 Inverse Solution for the Seven-Joint DLR-Kuka Robot . 89
3.4 Eight-Joint Robot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
vii
viii Contents
Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247
E Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 419
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435
Preface
1
2 Preface
In the appendix, we derive the geometric Jacobian matrix for the six-
joint elbow manipulator and for the DLR-Kuka seven-joint robot. We
also include solutions to selected exercises.
Acknowledgements
5
6 1 Introduction
Prismatic joint
Tool
Revolute joints
Fig. 1.1: Surgical saw mounted to a robotic arm. The arm itself consists
of revolute joints. The saw is mounted to a passive component (prismatic
joint). The surgeon moves the saw manually (sliding motion only). By
construction, the manual motion of the saw is restricted to a single axis.
[1]
1.1 Robots for Navigation 7
Remark 1.1
z0
y0
x0
In the figure, notice that the first two joint axes (revolute joints) in-
tersect in space. Furthermore, the axes of joints two, three and four
are parallel. We also see the axes x0 , y0 and z0 of a base coordinate
system.
Drilling direction
Fig. 1.3: Navigating a surgical drill. Femur figure from [2, Fig. 245].
the drilling direction). The target is on the left side. The tumor is well
visible in the CT image, but often not visible in an x-ray image. To
remove the tumor, we must guide the drill. Notice also that the drill
must not penetrate the joint surface. Here, the pre-operative 3D im-
ages must be matched to the intra-operative situation.
Figure 1.4 shows the head of a femur bone. A small gap between the
head and the rest of the bone is visible (see arrow). Epiphyseolysis is
a condition, in which a fracture along this gap results in a downward
slippage of the femur head (figure 1.4) [3].
To stabilize the femur head, screws are used (figure 1.5). However,
placing the screws is difficult. The reason is that only the surface of
the bone consists of hard material, the cortical bone. The interior of
the bone, called cancellous bone, or spongy bone, is much softer. The
screw tip must reach (but not penetrate) the cortical bone surface.
Figure 1.6 shows a cross-section (slice) of a CT-data set, side-by-side
with the corresponding MRI slice. The cortical bone is visible as a
ring (white/light grey in the CT image, and black in the MRI data).
The spongy bone is inside this ring. When placing a screw (as in fig-
ure 1.5), the tip of the screw should reach, but not penetrate the thin
cortical bone layer.
Figure 1.7 shows a similar application in spine surgery. Two or more
vertebrae must be stabilized. Screws through the pedicle hold a fixa-
1.1 Robots for Navigation 9
Fig. 1.6: Cortical bone and spongy bone in CT (left) and MRI (right).
In the CT scan, the cortical bone shows up white (black in MRI). The
spongy bone is inside.
1.1 Robots for Navigation 11
Body
Costal fovea
Fig. 1.7: Placement of pedicle screws. The insertion path for the screws
(red lines) must not touch the spinal canal. Source of drawing [2, Fig.
82].
Bony wedge
Cutting plane
Fig. 1.8: Corrective osteotomy for the femur bone. Source of the draw-
ing: [2, Fig. 245].
12 1 Introduction
J5
J3 J4
J2
J1
base coordinate system. Notice that infrared tracking not only outputs
the position of the pointer tip of the marker, but also the orientation of
the marker in space. This means we not only know the tip coordinates
(e.g. as x−, y−, z-coordinates) with respect to the camera coordinate
system, but we also obtain information about the angle and the point-
ing direction of the marker in space. To output such angular inform-
ation, standard tracking systems use matrices or quaternions. Typical
systems report the marker position to within 0.1 mm.
Figure 1.12 shows the principle of radiological navigation. A first
marker of the infrared tracking system is attached to the C-arm. The
second marker is attached to the saw. After making a skin incision
(opening the skin), a third marker is rigidly attached to the bone. Then
an x-ray image of the bone is taken with the C-arm. Given the position
of the C-arm in space, we can compute the exact spatial coordinates
of points on the bone visible in the C-arm image. These coordinates
are computed with respect to the (infrared) camera coordinate system.
We can now steer the robot to a predefined point on the bone surface,
or guide the robot to find a pre-planned angle with respect to the bone.
14 1 Introduction
Marker 1
Marker 2
Bone
Fig. 1.11: Infrared tracking system with two markers. One of the mark-
ers is attached to a bone.
If the bone moves during the procedure, we can calculate and correct
the displacement by subtraction.
When attaching a marker to the robot or the saw, we do not know the
spatial reference between the robot’s coordinate system, and the in-
ternal coordinate system of the infrared tracking camera. The process
of finding the spatial mapping between the two coordinate systems is
called hand-eye calibration. A similar problem occurs when calibrat-
ing C-arm images with respect to the infrared marker.
As noted above, tumors are often not visible in x-ray images, but well
visible in 3D CT images taken before the operation. The navigation
problem is now to align the bone in the x-ray image to the same bone
in the CT image. The alignment should be done such that the two im-
ages match, i.e. corresponding structures should be at the same posi-
tion. This process is called image registration. Not only x-ray images
are registered to CT images, but any pair of image modalities (CT,
MRI, ultrasound and many other modalities) can be considered in this
context.
1.1 Robots for Navigation 15
C-arm
Saw Bone
Fig. 1.12: Radiologic navigation. Markers for the infrared tracking sys-
tem are attached to the saw, the C-arm and the bone.
z J1
x J3
y J2 J4
J5 Needle
Frame base
Head-to-frame fixture
Fiducial markers
Target
Frame base
Localizer frame Fiducial markers
rods. Six of these rods are horizontal, and three of them are oblique
(see figure).
After taking a CT or an MR image (with the frame and localizers in
place), we see vertical cross sections of the head (figure 1.14, right).
Vertical cross sections (with the patient on the treatment table, as in
figure 1.14) are also called axial cross sections. The three viewing
directions for imaging are shown in figure 1.15.
In the cross-sections, we will see points (indicated as small circles in
figure 1.14, right), stemming from the fiducial rods. From the distance
between two adjacent points (one from an oblique rod and one from a
horizontal rod) we derive the z-coordinate of the image cross-section.
With similar methods, we can obtain x- and y-coordinates of the target
in the image.
In the next step of the procedure, we remove the localizer box with
the fiducial rods, and attach the jointed mechanism to the frame base.
As noted above, the jointed mechanism carries the instrument (e.g. a
biopsy needle or an electrode).
Having determined the coordinates of the target in the images, we
can compute the angle settings for the jointed mechanism and insert a
needle along a predefined path. Again, calibration methods are needed
to map the CT/MR coordinate system to the frame coordinate system.
Example 1.1
The joints J4 and J5 in figure 1.13 are revolute joints. The joint angles
of J4 and J5 determine the orientation of the needle. The needle tip
moves when we change the values of the angles. Let θ4 and θ5 be
the values of the joint angles for J4 and J5 . The current values for
θ4 and θ5 can be read from a scale imprinted onto the metal frame.
We can compute the position of the needle tip for given θ4 , θ5 . The
position is described with respect to the coordinate system shown in
the figure. We can also derive a closed form expression, stating the
coordinates of the needle tip as a function of θ4 , θ5 . Similarly, we can
include the prismatic joints into the computation. This derivation of
a closed formula is one example of a forward kinematic analysis of
a mechanism. Conversely, computing angle values from a given tip
18 1 Introduction
Coronal
Sagittal Axial
Coil
Markers
tient correlates to the tumor motion. Then we can use skin motion as
a surrogate signal, to track the internal target.
Linear
accelerator
Beam
collimator
Patient
couch
As noted above, robots are not only used for navigation. Using ap-
propriate sensors, it is not difficult to record all hand motions of a
surgeon during an intervention. One can then replicate hand motion
by a robotic arm. With this, a number of problems in microsurgery
can be addressed. An example is shown in the next figure. The motion
of the surgeon can be down-scaled. Assume the motion range of the
surgeon’s hand is 1 cm, and our robot replicates the same motion, but
down-scales it to a range of 1 mm. Thus, delicate interventions can be
performed with very high accuracy, see figure 1.18.
1.3 Robots for Imaging 21
Fig. 1.18: Replicating the surgeon’s hand motion with a robotic inter-
face mounted to the patient couch (da Vinci Surgical System, Intuitive
Surgical, Inc.)
Source
Robotic manipulator
Flat panel
detector
X-ray source
Exercises
Fig. 1.24: Joint angle θ4 for the neurosurgical frame in figure 1.13
Exercise 1.1
The joint angles (θ4 and θ5 ) for the neurosurgical frame in figure 1.13
are shown in figure 1.25 and figure 1.24.
Place a coordinate system at the centroid of the frame base, with axes
(x-,y-,z-) as shown in figure 1.25.
Assume the distance of the needle tip (figure 1.25) from the coordinate
origin is 100 mm.
a) For angle values taken from the set θ4 ∈ {−90, 0, 90} and θ5 ∈
{0, 90, 180}verify that the expression
cos(θ5 )
100 · − sin(θ4 ) sin(θ5 ) (1.1)
cos(θ4 ) sin(θ5 )
gives the coordinates of the needle tip in mm.
b) Given the coordinates
px
p = py
(1.2)
pz
of a point in space, compute the joint angles θ4 and θ5 for reaching
this point.
Hint: Assume kpk = 100. Start from the equation
px cos(θ5 )
py = 100 · − sin(θ4 ) sin(θ5 ) (1.3)
pz cos(θ4 ) sin(θ5 )
The first line of equation 1.3 gives
c) Discuss the case kpk 6= 100 with the prismatic joints in figure 1.13.
Exercise 1.2
Exercise 1.3
Detector
Figure 1.26 illustrates the input data for reconstruction. For each
angle, we have a projection image. The goal is to compute a gray
level for each voxel, i.e. values for the variables xi . Input data are the
pixel values measured at the detector. In the figure, the pixel value for
the topmost pixel is the sum of the values x1 , . . . , x4 in the first row.
In the case of several projection images, we obtain a linear system of
equations of the following form
Here, the values bi are the pixel values measured at the detector, the
coefficients ai j are constant, and each ai j is either 0 or 1. (We could
1.4 Rehabilitation and Prosthetics 29
use other values for ai j , for example a factor expressing the length of
the path of the ray through the corresponding voxel.) The resulting
linear equation system can be solved in many ways. The ART (algeb-
raic reconstruction technique) algorithm is particularly useful for this
application. To simplify the notation, we assume there are only two
equations and two variables.
a11 x1 + a12 x2 = b1
a21 x1 + a22 x2 = b2
(1.8)
x2 x2
a11 x1 + a12 x2 = b1
Start point
x(0) = [0, 0]T
x1 x1
a21 x1 + a22 x2 = b2 x(1)
Fig. 1.27: ART algorithm. Left: visualization of the two lines in equa-
tion 1.8. Right: Projecting the start point (0,0) onto one of the lines.
ai x(n) − bi
x(n+1) = x(n) − ai (1.9)
a2i
30 1 Introduction
We write the scalar product of ai and x(n) as ai x(n) instead of aTi x(n)
(see also the notation page in appendix E). Likewise, a2i denotes the
scalar product aTi ai . This iteration formula computes the (n + 1)-th
estimate x(n+1) from the n-th x(n) , i.e. it projects the point x(n) onto
the i-th line of the system, which is given by ai and bi . The line index
i can be chosen at random. The coefficients ai j form a matrix which
will be called A.
Instead of equation 1.8, we can thus write Ax = b. The ART algorithm
is now given by:
1. Start by setting x(0) = (0, 0)T .
2. Choose a random index i and compute x(n+1) from x(n) via equa-
tion 1.9
3. Compute the forward projection from the (n + 1)-th estimate, i.e.
b(n+1) = Ax(n+1)
4. Iterate steps 2 and 3 until the difference between new forward pro-
jection b(n+1) , computed in step 2, and values measured at the de-
tector is below tolerance
Notice that the process of reconstruction just described can be rever-
ted. This will give rise to interesting applications in registration and
navigation. We will discuss such applications in the chapter on regis-
tration.
a) Explain the difference between the two types of projections (pro-
jection in figure 1.26) and projection in figure 1.27 in the ART
algorithm.
b) Explain why reconstruction is needed to compute 2D slice images
(i.e. regular CT image stacks) from 2D x-ray images.
c) Explain why we cannot cancel some of the ai ’s in equation 1.9.
d) Prove equation 1.9.
e) Implement the ART algorithm for planar images. For simplicity,
first use synthetic images such as the binary image of a rectangle
shown in figure 1.28 (left) as input.
f) Compare different implementations of linear equation solvers,
such as the one provided in MATLAB (linsolve()) (or any
other linear equation system solver) to the implementation of the
ART algorithm for larger images.
References 31
0000000000000 000*000000000
0000000000000 000*000000000
0011111111100 001*111111100
0010000000100 0010*00000100
0010000000100 0010*00000100
0010000000100 0010*00000100
0011111111100 00111*1111100
0000000000000 00000*0000000
0000000000000 00000*0000000
References
2.1 Matrices
33
34 2 Describing Spatial Position and Orientation
zs
ys
z S xs
p = (px , py , pz )T
y
B x
Fig. 2.1: Vector p describes the position of the coordinate system S with
respect to the base coordinate system B.
100 px
0 1 0 py
M=
0 0 1
. (2.1)
pz
000 1
xs
ys
z S
zs
p
y
B x
Example 2.1
ys
S zs
z xs
B x
Consider the surgical saw in figure 1.1 (first figure of the introduction
chapter). Our goal is to move the saw with a robot. Thus, we must
command the positions and orientations of the saw. First we attach a
coordinate system S to the saw. In figure 2.3, the saw is simply shown
as a block. The origin of S is a vertex of this block. The axes xS , yS and
zS of S coincide with edges of the block.
According to our notation, we find that the following matrix B MS de-
scribes the position of the saw:
0 0 1 px
B
0 1 0 py
MS = −1 0 0 pz
(2.4)
0 00 1
The notation B MS indicates that the matrix gives the position and ori-
entation of S with respect to B.
Notice that here we have again rotated about the y-axis, but now in the
opposite direction. We see that we need a convention indicating the
direction of rotation. We will do this next, after introducing conven-
tions for rotations about fixed axes.
2.2 Angles 37
2.2 Angles
B
Y PRS = (px , py , pz , α, β , γ)T (2.6)
is called the Yaw-Pitch-Roll vector, or YPR-vector. However, there are
several other conventions (see the exercises at the end of this chapter).
We will soon see that the matrix representation is often preferable over
the YPR-representation.
The medical literature has its own naming conventions for rotations.
They are:
medicine engineering
flexion (+) / extension (-) yaw
varus (+) / valgus (-) pitch
rotation (+) / derotation (-) roll
38 2 Describing Spatial Position and Orientation
Suppose we replace the saw by a drill (see figure 2.4). The drill has
its own coordinate system S. At the tip of the drill, we place a new
coordinate system S0 . We have the position of S with respect to the
base B. Now we would like to compute the position of S0 with respect
to the base system B. Thus, suppose we have defined two matrices
B M and S M 0 .
S S
yS
S zS
z xS
yS0
S0 zS 0
y
xS0
B x
Here we have again used the convention that the matrix B MS describes
a transition from coordinate system B to coordinate system S. The
first of the two matrices (B MS ) has been computed above. The second
matrix is also easy to compute. From the figure, we see that it is simply
given by:
1 0 0 p0x
S
0 1 0 p0y
MS 0 =
0 0 1 p0z
(2.7)
000 1
Here p’ = (p0x , p0y , p0z )T denotes the origin of S0 , given in coordinates
of S.
2.2 Angles 39
S
BM
S
SM
S0
S0
BM
S0
Now the matrix we are looking for (drill tip with respect to base) is
simply given by the product of the two matrices already computed,
i.e.
B
MS0 = B MS S MS0 . (2.8)
This is illustrated in figure 2.5. It is easy to verify this equation with
the above interpretation of matrix columns.
Above we looked at relative transformations between coordinate sys-
tems. The following example discusses a different situation: a vector
is given in a coordinate system S, and our goal is to express its co-
ordinates in another coordinate system.
Example 2.2
Assume the vector q in figure 2.6 is given with respect to the drill
coordinate system S. q points along the z-direction of S, and
0
S
q = 0 .
(2.9)
1
40 2 Describing Spatial Position and Orientation
ys
S
q = zs
z xs
B x
B
q =B MS S q. (2.12)
2.3 Linkages
Example 2.3
Consider the robot shown in figure 2.7. It only has a single joint, and
a single link. But it does have a gripper, mounted to the end of the
single link.
This robot is not very interesting (it cannot even close its gripper) but
it illustrates how we apply coordinate systems. A null position of a
robot is a fixed reference position. Assume the null position of our
one-joint robot is the position shown. In this null position, the joint
angle is zero. If we change the joint angle, our robot will start to move.
42 2 Describing Spatial Position and Orientation
zB
yB
xB
Fig. 2.8: Base coordinate system for the robot in figure 2.7. The robot
has a single revolute joint. The axis of rotation of this joint is the z-axis
of the base system.
zG
yG
xG
Fig. 2.9: Gripper coordinate system for the robot in figure 2.7.
zS1
yS1
xS1
Fig. 2.10: Intermediate coordinate system S1 for the robot in figure 2.7.
Notice that the origin of S1 coincides with the origin of B. S1 rotates with
the robot, while B remains fixed.
The gripper system is fixed at a point in the center between the two
gripper jaws. This point is called the tool center point, and depends on
the tool we use, e.g. a saw, a drill, or a gripper.
Throughout this book, we will often use the term effector instead of
gripper. We can replace the gripper by a saw or a drill, or some other
tool. An effector can thus be any tool mounted to the end of a kin-
ematic chain.
The placement of the system S1 , as the robot moves, is shown in fig-
ure 2.11.
zB , zS1
yB
xS1
xB
θ
Fig. 2.11: The intermediate coordinate system S1 rotates with the robot,
(angle θ ). The two z-axes (of B and S1 ) coincide. The y-axis of S1 is not
shown.
44 2 Describing Spatial Position and Orientation
We now look at the scene from above: the placement of S1 within the
base system is shown in figure 2.12:
yB
yS1 xS1
θ sin θ
cos θ xB
We consider the vector xS1 in the figure. In terms of the base, this vec-
tor has the x-coordinate cos θ , and the y-coordinate sin θ . Likewise,
we see that the y-vector of S1 is given by
− sin θ
(2.13)
cos θ
Putting this into a 4 × 4-matrix as described above, we obtain the pos-
ition of the intermediate system S1 with respect to the base, depending
on the joint angle θ :
cos θ − sin θ 0 0
B
sin θ cos θ 0 0
MS1 =
0
(2.14)
0 1 0
0 0 01
We pause a minute to look at this matrix. Clearly, the z-vector should
not move as we rotate about it. And indeed it does not move, namely
the third column of the matrix simply shows the unchanged z-vector.
The first and second columns of the matrix show the motions of the x-
and y-vectors, depending on the angle θ .
This completes the first step of our computation. For the second step,
we find the transformation from S1 to the gripper. Since this trans-
formation is only a translation by the link length (called L), the matrix
is given by:
2.3 Linkages 45
100L
S1
0 1 0 0
MG =
0 0 1 0 (2.15)
0001
Thus the desired position of the gripper in terms of the base coordinate
system is given by the matrix product
cos θ − sin θ 0 L cos θ
B
sin θ cos θ 0 L sin θ
MS1 · S1 MG =
0
(2.16)
0 1 0
0 0 0 1
Note that this matrix contains the angle θ . It expresses the location
of the gripper coordinate system, depending on the joint angle of the
robot. We have completed the forward kinematic analysis in a simple
example. In general, the goal of forward analysis is to find the location
of the effector, given values for each joint angle.
1 0 0 0 cos θ 0 sin θ 0
0 cos θ − sin θ 0 0 1 0 0
R(x, θ ) =
0 sin θ cos θ
R(y, θ ) =
− sin θ
0 0 cos θ 0
0 0 0 1 0 0 0 1
46 2 Describing Spatial Position and Orientation
cos θ − sin θ 00
sin θ cos θ 0 0
R(z, θ ) =
0
(2.17)
0 1 0
0 0 01
Even more straightforward is the derivation of the elementary transla-
tion matrices:
1 0 0 px 100 0
0 1 0 0 0 1 0 py
T (px , 0, 0) =
0 0 1 0
T (0, p y , 0) =
0 0 1 0
000 1 000 1
100 0
0 1 0 0
T (0, 0, pz ) =
0 0 1 p z
(2.18)
000 1
Suppose now we have a robot with three joints. Such a robot is shown
in figure 2.13. We call this robot the elbow manipulator. Our goal is
to express the position and orientation of the gripper with respect to
the base coordinate system. We place a coordinate system into each
link of the robot. We obtain a chain of coordinate systems, leading
from the base system to the gripper system in several steps. Notice
that this robot has only three joints, and typical robots have six joints.
However, all computations for the three-joint robot can later be used
for more complex robots.
The three joint angles are called θ1 , θ2 , θ3 . The position shown in the
figure is the null position of the robot, meaning θ1 = θ2 = θ3 = 0.
As above, we define intermediate coordinate systems S0 = (x0 , y0 , z0 )
and Si = (xi , yi , zi ), i = 1, . . . , 3 (see figure 2.14).
Think of the first system B = S0 = (x0 , y0 , z0 ) as the base system. It
is attached to the table, and does not move, although it rests with its
origin at the upper tip of the first link. This may seem odd, but we do
2.4 Three-Joint Robot 47
Fig. 2.13: Three-joint robot (elbow manipulator). All three joints are
revolute joints
zG
yG
z0 z3
xG
y0 y3
d4
x0 x3
a2
z1 z2
x1 x2
y1 y2
Fig. 2.14: Joint coordinate systems for transforming between base sys-
tem S0 = (x0 , y0 , z0 ) via Si = (xi , yi , zi ) to the gripper coordinate system
SG = (xG , yG , zG ).
48 2 Describing Spatial Position and Orientation
this to simplify the overall matrix product, i.e. we thus do not need an
extra transformation leading from the table to the tip of the first link.
Our goal is to derive a matrix which transforms from the base co-
ordinate system B = S0 to the gripper system SG . Having defined the
intermediate matrices Si , we can partition this goal into a series of
small steps. The first step is the transformation from S0 to S1 . We will
do this next.
The rotation axis of the first joint is the axis z0 . To get from S0 to S1 ,
we must take into account any rotations about this first axis. As in the
case of our single joint robot in the above example, the matrix R(z, θ )
in equation 2.17 describes this rotation. But now we are not done,
since we have not reached S1 yet. To do so, we must rotate about the
x0 -axis by the fixed angle of -90 degrees (see figure 2.14). Overall, the
following matrix product gives the transformation from S0 to S1 .
cos θ1 − sin θ1 0 0 1 0 00
S0
sin θ1 cos θ1 0 0 0 0 1 0
MS1 = R(z, θ1 )R(x, −90) = 0
0 1 0 0 −1 0 0
0 0 01 0 0 01
(2.19)
Hence (simplifying the notation S0 MS1 to 0 M1 ):
cos θ1 0 − sin θ1 0
0
sin θ1 0 cos θ1 0
M1 =
0 −1
(2.20)
0 0
0 0 0 1
Similarly (denoting the constant link lengths by a2 and d4 ) we obtain:
cos θ2 − sin θ2 0 a2 cos θ2
1
sin θ2 cos θ2 0 a2 sin θ2
M2 = 0
(2.21)
0 1 0
0 0 0 1
cos θ3 0 sin θ3 0
2
sin θ3 0 − cos θ3 0
M3 = 0 1
(2.22)
0 0
0 0 0 1
2.5 Standardizing Kinematic Analysis 49
100 0
3
0 1 0 0
MG =
0 0 1 d 4 (2.23)
000 1
Multiplying, we have (abbreviating s1 for sin θ1 etc.) :
0
M3 =0 M1 · 1 M2 · 2 M3 (2.24)
c1 (c2 c3 − s2 s3 ) −s1 c1 (c2 s3 + s2 c3 ) a2 c1 c2
s1 (c2 c3 − s2 s3 ) c1 s1 (c2 s3 + s2 c3 ) a2 s1 c2
(2.25)
−s2 c3 − c2 s3 0 −s2 s3 + c2 c3 −a2 s2
0 0 0 1
The last equation can be simplified with the trigonometric addition
formulas. Therefore (setting s12 for sin(θ1 + θ2 ), etc.):
c1 c23 −s1 c1 s23 a2 c1 c2
0
s1 c23 c1 s1 s23 a2 s1 c2
M3 = −s23 0 c23 −a2 s2
(2.26)
0 0 0 1
For the final result, we must multiply with the constant transformation
3 M , given above, which transforms from S to the gripper coordinate
G 3
system:
c1 c23 −s1 c1 s23 a2 c1 c2 + d4 c1 s23
0
s1 c23 c1 s1 s23 a2 s1 c2 + d4 s1 s23
MG = −s23 0 c23 −a2 s2 + d4 c23
(2.27)
0 0 0 1
You may wish to check equation 2.27 for various settings of the angles
θ1 , θ2 , θ3 (see also the exercises at the end of this chapter).
Oi zi
Joint normal
zi−1
We noted above that the main observation in the original paper by De-
navit and Hartenberg [1] is that four elementary transformations will
always suffice for transforming from Si−1 to Si . Notice that we must
follow the rules when placing Si−1 and Si . The four specific trans-
formations to be given now will not suffice, if you did not follow the
rules. (To understand why, you should check that the rules amount to
‘removing’ two degrees of freedom.)
Denavit and Hartenberg show that the following elementary transla-
tions/rotations always suffice:
1. Translate along zi−1 by an amount di
2. Rotate about zi−1 by θi
3. Translate along xi by ai
4. Rotate about xi by αi
In our example, the values di , ai , αi are constant, while θi is the (vari-
able) joint angle. We will hence call θi the joint variable. We have
given names to the elementary transformations, and we have written
them as matrices in equations 2.17 and 2.18. Thus the transformations
can be written as follows:
T (0, 0, di )
R(z, θi )
T (ai , 0, 0)
R(x, αi ) (2.28)
52 2 Describing Spatial Position and Orientation
i−1
Mi = T (0, 0, di )R(z, θi )T (ai , 0, 0)R(x, αi ) (2.29)
transforms from Si−1 to Si .
With equations 2.17 and 2.18, we can write
100 0 cθi −sθi 00 1 0 0 ai 1 0 0 0
0 1 0 0 sθi cθi 0 0
i−1
Mi = 0 1 0 0 0 cαi −sαi 0
0 0 1 d i 0 0 1 0 0 0 1 0 0 sαi cαi 0
000 1 0 0 01 000 1 0 0 0 1
(2.30)
Finally, we obtain
cθi −cαi sθi sαi sθi ai cθi
i−1
sθi cαi cθi −sαi cθi ai sθi
Mi =
0 sαi
. (2.31)
cαi di
0 0 0 1
We see that the matrix describing a robot with six joints is a product
of six matrices as in equation 2.31. Although we have taken care to re-
duce the number of elementary transformations, the matrix product is
not simple. However, having set up the rules for the placement of co-
ordinate systems, we can arrive at a yet more compact description of a
robot. We simply set up a table containing the parameters αi , ai , di , θi .
Each line in the table corresponds to one joint, and we obtain a table
of the form in table 2.1.
i αi ai di θi
1 α1 a1 d1 θ1
2 α2 a2 d2 θ2
.. .. .. .. ..
. . . . .
It should be noted that a number of special cases can arise when ap-
plying the DH-rules. For example, the axes zi−1 and zi can intersect
or be parallel, or even coincide. In all such cases, the rules can be re-
laxed, and we obtain two or more equivalent ways for fixing the axis
xi or the origin Oi .
We saw that we can describe the position of the gripper, given angles
θ1 , θ2 .... We now wish to do the opposite. Namely, given a desired
position for the gripper, how do we choose the angles? In general, this
is what we need to program a robot. The next example discusses a
simple case.
Example 2.4
Fig. 2.16: Calculating the joint angle θ for reaching a given point p
sin θ py
= (2.32)
cos θ px
p p
From this, we could conclude tan θ = pxy , and θ = arctan pxy .
But this simple solution has a drawback. By definition, the function
arctan returns values in the range (−π/2, π/2). Thus, it can never
return a value larger than π/2. By consequence, our robot can never
reach the point p shown in the figure, since its null position points
along the x-axis. In fact, the entire half-plane to the left of the y-axis
is lost for our robot.
There is a standard way to repair this. We extend the function arctan.
The function atan2 is defined on the basis on the arctan-function. This
function atan2 takes two arguments x and y, hence the name. It is
defined by distinguishing several cases for the values of the arguments
x and y.
atan2(y, x) = (2.33)
0 for x = 0, y = 0
π/2 for x = 0, y > 0
3 · π/2 for x = 0, y < 0
arctan(y/x) for x > 0, y ≥ 0
2π + arctan(y/x) for x > 0, y < 0
π + arctan(y/x) for x < 0, y ≥ 0
π + arctan(y/x) for x < 0, y < 0
θ = atan2(py , px ). (2.34)
2.7 Quaternions 55
You should verify that our one-joint robot can reach all positions on a
full circle with this definition.
A matrix equation as in equation 2.27 is called the forward kinematic
solution of the robot. The result in equation 2.34 shows a simple
example of an inverse kinematic solution. Thus the inverse solution
states the required joint angles necessary for reaching a given position
and orientation.
We will discuss the forward and inverse kinematics for medical link-
ages in more detail in chapter 3.
2.7 Quaternions
a = (a1 , a2 , a3 , a4 )T (2.35)
We can add quaternions just like regular 4-vectors i.e.
a + b = (a1 + b1 , a2 + b2 , a3 + b3 , a4 + b4 )T (2.36)
We would like to define a multiplication for quaternions. One may
argue that we already have a multiplication, i.e. the standard scalar
product for two vectors. However, this is not what we need. We need
a multiplication for two quaternions such that the result is also a qua-
ternion, i.e. a 4-vector. The cross-product for 3-vectors could be a can-
56 2 Describing Spatial Position and Orientation
didate, but it is only defined for 3-vectors. Here we define for two
quaternions a, b:
a1 b1 − a2 b2 − a3 b3 − a4 b4
a1 b2 + a2 b1 + a3 b4 − a4 b3
ab =
a1 b3 − a2 b4 + a3 b1 + a4 b2
(2.37)
a1 b4 + a2 b3 − a3 b2 + a4 b1
This is certainly not easy to memorize. To simplify memorizing the
formula in equation 2.37, quaternions can be written in a slightly dif-
ferent, but equivalent way.
Instead of a = (a1 , a2 , a3 , a4 )T , we write
a = a1 + a2 i + a3 j + a4 k (2.38)
Note that i, j and k are symbols, not numbers, and we have a table
stating symbolic product rules for i, j and k. The symbolic rules are
ii = −1 j j = −1 kk = −1
ij = k jk = i ki = j (2.39)
ji = −k k j = −i ik = − j
Now we can multiply i, j and k just as they were numbers, and one can
easily verify that the product ab for quaternions defined in this way is
equivalent to this new definition.
The new rules are somewhat easier to memorize than the original
definition in equation 2.37, and they are reminiscent of the rules for
the positive/negative sense of rotation. For all purposes other than re-
membering the specific multiplication above, we will remain with our
first definition in equation 2.37. Hence, from now on, a quaternion is
simply a 4-vector.
Notice that in general ab 6= ba.
Just as for any 4-vector, we have a norm, namely
q
kak = a21 + a22 + a23 + a24 (2.40)
Unit quaternions are quaternions having norm 1. We represent rota-
tions in the following way. Given a 3-vector p = (px , py , pz )T . Set
2.7 Quaternions 57
q = apa−1 (2.42)
defines quaternion rotation, where a is the (unit) quaternion by which
we rotate, and p is the object of rotation. Thus the quaternion a takes
the role of a 3 × 3-rotation matrix. Notice that a quaternion is much
more compact than a 3 × 3-matrix, which is why quaternions are also
often used in graphics.
To understand exactly which rotation is represented by a given unit
quaternion, we observe that for a = (a1 , a2 , a3 , a4 )T , the component
a1 defines the angle, and the 3-vector (a2 , a3 , a4 )T defines the axis of
the rotation. Details are explained in the following example.
Example 2.5
Suppose we wish to rotate about the z-axis of the base system by the
angle θ = π. To select the axis of rotation, we set p = (0, 0, 1)T . Then
the quaternion
p T
a = cos π/2, sin π/2 (2.43)
kpk
58 2 Describing Spatial Position and Orientation
Exercises
The 3 × 3 matrix
nx ox ax
R = ny oy ay (2.44)
nz oz az
is called orthogonal if the product RRT is the 3 × 3 unit matrix. Here
RT denotes the transpose of R. The upper left 3 × 3 matrix in a ho-
mogeneous matrix must always be orthogonal, since the first three
columns in the homogeneous matrix represent vectors in an ortho-
gonal coordinate system.
a) Show that for an orthogonal matrix R as in equation 2.44, we have
RxRx = xx for any 3-vector x. Interpret this observation.
b) From the definition of an orthogonal matrix, we see that the inverse
of an orthogonal matrix is its transpose, i.e. RT = R−1 . Based on
this, derive the inverse of a homogeneous 4 × 4 matrix.
c) Show that RxRy = xy for any 3-vectors x, y, if R is orthogonal.
(z) twice. But notice that now our rotations do not refer to the axes
of the base system, but in each case we rotate about the new axes. In
this case (as opposed to the YPR angles in the previous exercise), we
multiply the matrices from left to right.
a) Given base and saw coordinate systems as in figure 2.3, derive the
Euler angles for the saw coordinate system.
Compare to the YPR angles for the same systems.
b) Derive and implement a conversion between 3×3 rotation matrices
and Euler angles. The derivation closely follows exercise 2.4, with
one difference: you will have to start from
nx ox ax
R(z, α)R(y, β )R(z, γ) = ny oy ay (2.47)
nz oz az
Notice here that the matrices appear in reverse order from what
you would expect! This corresponds to the fact that now we rotate
about the new axes each time.
Notice also that some text books define Euler angles with the angle
conventions z, x, z instead of z, y, z. In fact there are 12 different
ways to define the Euler angles, all equivalent.
a) Verify the rotation results in example 2.4 for the input vectors
(1, 0, 0)T and (0, 0, 1)T .
b) Derive a conversion from a quaternion to a 3 × 3 rotation matrix.
Hint: If u = (ux , uy , uz )T is a given axis, and θ is a given angle, we
set c = cos θ and s = sin θ .
Then the 3 × 3 matrix
c + u2x (1 − c) ux uy (1 − c) − uz s ux uz (1 − c) + uy s
R = uy ux (1 − c) + uz s c + u2y (1 − c) uy uz (1 − c) − ux s
uz ux (1 − c) − uy s uz uy (1 − c) + ux s c + u2z (1 − c)
(2.48)
describes a rotation about u by θ . Given a quaternion q, extract
the axis u and the angle θ , and insert into the definition for R.
62 2 Describing Spatial Position and Orientation
Derive the forward equation (equation 1.3) for the neurosurgical frame
in figure 1.13.
Summary
the joint angles for a given position and orientation. Inverse kinemat-
ics typically requires a forward analysis. Stated in simpler terms (with
the notation in figure 2.16), forward kinematics finds the position and
orientation of the effector, given the joint angles, while inverse kin-
ematics finds the joint angles given the target position and orientation.
For the kinematic analysis of a single-joint robot, we introduced three
coordinate systems, namely the base system, gripper system and an
intermediate system. The coordinate systems thus form a chain. We
assign a matrix to each pair of consecutive systems in this chain. Then
the product matrix gives the transformation from the base to the grip-
per. The analysis is very similar for robots with several joints.
The atan2-function is an extension of the arctan function. It takes two
arguments and implements the inverse kinematics for a robot with
a single revolute joint. More complex robots can be analyzed with
this function as well, again on the basis of a forward analysis. Fur-
thermore, the atan2-function is a general tool for solving kinematic
equations.
Quaternions are 4-vectors. To add two quaternions, we simply add
their components. The multiplication of two quaternions resembles
the cross-product for 3-vectors. Quaternion multiplication is not com-
mutative. The rotation defined by a quaternion is derived from equa-
tion 2.42 and equation 2.43. If a = [a1 , a2 , a3 , a4 ]T is a quaternion, then
the first element a1 describes the angle of rotation, while [a2 , a3 , a4 ]T
describe the axis. To convert from a (unit) quaternion to a 3x3-rotation
matrix, we must set up a rotation matrix given an angle and an axis
(exercise 2.7b). The reverse direction of the conversion consists of
extracting an angle and an axis from a rotation matrix (exercise 2.7c).
Notes
References
67
68 3 Robot Kinematics
c1 c23 −s1 c1 s23 a2 c1 c2 + d4 c1 s23 nx ox ax px
s1 c23 c1 s1 s23 a2 s1 c2 + d4 s1 s23 ny oy ay py
(3.2)
−s23 0 c23 −a2 s2 + d4 c23 = nz oz az pz
0 0 0 1 0 0 0 1
− s1 = ox (3.3)
From the second entry of the second matrix line, we obtain
c1 = oy (3.4)
Again, we use our function atan2, which we have already applied in
several examples and exercises in chapter 2.
Set
θ1 = atan2(−ox , oy ) (3.5)
Once we have solved for θ1 , the value for θ1 becomes a constant.
There is an easy way to move θ1 from the left side of the equation to
the right side. We multiply both sides of the matrix equation with the
matrix 0 M−11 (see equation 2.20). This ensures that all variables are on
the left side, and the matrix on the left side is simplified considerably.
To further simplify the expressions, we remove the constant matrix
3 M−1 from the left side of the equation.
G
We thus transform the equation from which we started (see equa-
tion 3.1)
0
M1 · 1 M2 · 2 M3 · 3 MG = T (3.6)
into
1
M2 · 2 M3 = 0 M−1 3 −1
1 T · MG . (3.7)
We have
3.1 Three-joint robot 69
100 0
3 −1
0 1 0 0
MG =
0 0 1 −d4 (3.8)
000 1
and
nx ox ax −ax d4 + px
ny oy ay −ay d4 + py
T · 3 M−1
G = n o a −a d + p
(3.9)
z z z z 4 z
0 0 0 1
and set
nx ox ax −ax d4 + px nx ox ax p0x
ny oy ay −ay d4 + py ny oy ay p0y
= T0 . (3.10)
nz oz az −az d4 + pz = nz oz az p0z
0 0 0 1 0 0 0 1
Now
0 0 0 1
a2 s2 = −p0z
a2 c2 = p0x c1 + p0y s1 (3.13)
70 3 Robot Kinematics
Thus
−p0z p0x c1 + p0y s1
θ2 = atan2 , , (3.14)
a2 a2
where
p0x = −ax d4 + px ,
p0y = −ay d4 + py ,
p0z = −az d4 + pz . (3.15)
cos(θ2 + θ3 ) = az
sin(θ2 + θ3 ) = −nz (3.16)
so that
θ3 = atan2(−nz , az ) − θ2 (3.17)
This completes our inverse analysis of the three-joint robot.
Before we look into the analysis of more complex robots, we briefly
summarize the overall procedure in the analysis of our three-joint ro-
bot for future reference: we first set up the forward matrix equation.
In this equation, we found a pair of equations of the form si = u, and
ci = v, with constants u, v. We solved for the i-th joint angle on the
basis of the atan2-function. Then we moved the matrix i−1 Mi to the
right-hand side of the equation, since it no longer contained any un-
solved variables. To find an initial pair of equations (of the form si = u,
and ci = v), we could be somewhat more systematic. We could gen-
erate several matrix equations by moving matrices to the right side,
while looking for appropriate pairs of equations. We will call this
simple technique equation-searching.
In the next section, we will analyze a six-joint robot. We will see that
this robot can be regarded as an assembly of two three-joint robots in
our analysis, i.e. we can decouple the analysis for joints 1-3 and joints
4-6.
3.2 Six-joint robot 71
Figure 3.1 shows a six-joint robot. This six-joint robot can be obtained
by assembling two three-joint robots. If we take a closer look at the
kinematic structure, we see that the base (first three joints) is identical
to the three-joint elbow manipulator discussed above. We mount a
spherical wrist (with three joints) to the base. Figure 3.2 shows the res-
ulting kinematic structure, an elbow manipulator with spherical wrist
(see also exercise 2.8).
Joint 6
Joint 5
Joint 4
Joint 2
Joint 3
Joint 1
Forward Analysis
For the forward analysis of the six-joint robot, we must assign co-
ordinate systems to the joints. We have already derived the placement
of the coordinate systems S0 , S1 , S2 , S3 in figure 2.14.
Figure 3.3 adds the remaining coordinate systems for the hand as-
sembly (systems S4 , S5 , S6 ).
We obtain the following DH-table for the six-joint robot.
72 3 Robot Kinematics
z0
y0
x0
i αi ai di θi
1 −90 0 0 θ1
2 0 a2 0 θ2
3 90 0 0 θ3
4 −90 0 d4 θ4
5 90 0 0 θ5
6 0 0 d6 θ6
Since the hand assembly is also a (small) three joint robot, its analysis
should be simple. But can we decouple the kinematic analysis of the
two robots thus glued together? It is remarkable that this is indeed
possible.
The main idea for the inverse analysis is the following.
We have the target coordinate system S6 at the end of the hand as-
sembly (figure 3.3). We look at the z-vector of S6 . Now revert the
direction of this vector, to point into the opposite direction. We follow
this opposite direction until we reach the intersection point between
3.2 Six-joint robot 73
z6
y6
x6
z5
q y5
x5
z4
x4
y4
z0
z3
y0
x0 y3
x3
Fig. 3.3: Coordinate systems for the hand assembly of the six-joint robot.
the axes of joints five and six. Call this intersection point q, the wrist
center. q is also the origin of the coordinate systems S4 and S5 (see fig-
ure 3.3). Taking a second look at the kinematics of our robot, we see
that the axes for the last three joints intersect in one point. Thus q will
not move, when we move any of the axes four, five or six. This means
that the position of q only depends on the joints one, two and three.
It is important to note that the inverse analysis takes a homogeneous
4 × 4-matrix as an input, which gives the position and orientation of
the last coordinate system (tool system). In this matrix, we can find
the point q in the following way: we do have the gripper’s z-vector in
column three of the matrix (referred to as vector (ax , ay , az )T ) in the
matrix notation from equation 2.3). It is a unit vector. A vector point-
ing into the opposite direction is easy to find. Simply take the negat-
ive of the gripper’s z-vector, i.e. −(ax , ay , az )T . But now the distance
74 3 Robot Kinematics
between the origin of the gripper coordinate system and the point q
is a kinematic constant, and we can take it from our DH-table! Thus,
multiply the gripper’s z-vector by d6 , and subtract that from the vector
in column four of the matrix (also called vector p). Overall,
q = p − d6 a. (3.18)
All elements on the right side of this equation are known, and we
can thus compute q from the input matrix. So now we can find the
angles θ1 , θ2 , θ3 from the point q alone, by analyzing a three-joint
robot. After having found the angles θ1 , θ2 , θ3 , we rewrite the forward
matrix equation
nx ox ax px
0
ny oy ay py
M1 · · · 5 M6 =
nz oz az pz (3.19)
0 0 0 1
to
nx ox ax px
3 M · 4 M · 5 M = 2 M−1 · 1 M−1 · 0 M−1 ny oy ay py
4 5 6
(3.20)
3 2 1 n o a pz
z z z
0 0 0 1
This again is the kinematic equation for a three-joint robot (since
everything on the right hand side is now constant), and we have
already solved variants of such equations before.
This is the basic idea for solving the six-joint robot. However, some
subtleties are easy to overlook here: When we analyzed the three-joint
robot above, we had a full gripper matrix to start from. Now we only
have a single point, namely q. This will cause trouble, and we shall
see that we must work through the three-joint solution again, to adjust
to the new situation. The new solution will be much more detailed,
and will address the case of multiple solutions.
Specifically, look at the so-called left-shoulder and right-shoulder
configurations shown in figure 3.4.
Figure 3.5 shows two distinct configurations of the robot, where again
the resulting tool position is the same. We will call the configurations
in figure 3.5 the Elbow-up and the Elbow-down configurations.
3.2 Six-joint robot 75
A
x
y z
B
x
y
z
C
y
z
x
D
x
y z
Fig. 3.6: Hand configuration. First note that for configurations A and
D, the position of the hand coordinate system is the same. But the joint
angles are not the same! To see why, we move through three steps, i.e.
from A to B, then from B to C, finally from C to D. For A to B: rotate
joint J4 by 180 degrees. B to C: joint J5 rotates. C to D: rotate joint J6 ,
again by 180 degrees. Hence joints four five and six have all taken new
values.
As noted above, the solution for the six-joint robot proceeds in two
steps. First we solve for the angles θ1 , θ2 , θ3 . Our sole input at this
stage will be the point q and the desired configuration. Notice that this
kinematic problem is slightly different from our above solution for the
three-joint robot (see example 3.1), where our input was a target 4x4-
matrix. Having found θ1 , θ2 , θ3 , we will then solve for the remaining
parameters via equation 3.20.
From equation 3.2 we see that q has the coordinates
qx a2 c1 c2 + d4 c1 s23
q = qy = a2 s1 c2 + d4 s1 s23 . (3.21)
qz −a2 s2 + d4 c23
But we have a second equation for q, given our target position as a
homogeneous 4 × 4 matrix T. Given T, we can take the two vectors
p and a from this matrix, and apply equation 3.18, i.e. q = p − d6 a.
Recall that d6 is a kinematic constant from the DH- table. Overall we
have:
a2 c1 c2 + d4 c1 s23
p − d6 a = a2 s1 c2 + d4 s1 s23 (3.22)
−a2 s2 + d4 c23
Our next goal will thus be to extract θ1 , θ2 , θ3 from equation 3.22.
Solving for θ1
Our solution for θ1 is very simple, and we can apply the same tech-
nique as for the one-joint robot in chapter 2. We look at the robot from
above (figure 3.7).
Figure 3.8 further illustrates the situation. Projecting the entire con-
figuration onto the plane of the robot’s base plate (x-y-plane), we see
that the angle θ1 can be taken from the position of the point q directly.
After projecting q onto the x-y-plane, we see that θ1 can be computed
as
θ1 = atan2(qy , qx ) (3.23)
Equation 3.23 refers to one possible solution for θ1 . We can obtain
a second solution by changing the configuration of the shoulder (see
78 3 Robot Kinematics
Fig. 3.7: Solving for angle θ1 . The origin of the base system (O) and q
are marked (black dots).
y0 y0
q q
θ1 θ1 + π
O x0 O x0
a) b)
Fig. 3.8: Solutions for the angle θ1 . The robot has several ways for reach-
ing the given point q. Figure a shows the so-called left-shoulder config-
uration. To reach the right shoulder configuration (figure b) we move
the first joint to the angle θ1 + π, and then flip the second joint.
θ1 = atan2(qy , qx ) + π. (3.24)
From the definition of the atan2-function, we see that
We will next solve for θ3 rather than for θ2 , since this will simplify
the solution path.
Solving for θ3
We take a step back and look at the robot. In figure 3.9 we see that
the distance between q and the origin of the base system should only
depend on the angle θ3 , but not on the joint angles θ1 and θ2 .
z0 x0 -z0 -plane
y0
x0
q
θ2 ||q|| θ3
Fig. 3.9: Solving for angle θ3 . The distance kqk between the origin of the
base system and the point q only depends on θ3 .
We take the square of the norm on both sides of equation 3.21, and
obtain:
Solving for θ2
x = a2 cos(θ2 ) + d4 sin(θ2 + θ3 ),
z = −a2 sin(θ2 ) + d4 cos(θ2 + θ3 ) (3.30)
y = a2 cos(θ2 ) + d4 sin(θ2 + θ3 ),
z = −a2 sin(θ2 ) + d4 cos(θ2 + θ3 ) (3.31)
with y = qy / sin(θ1 ).
Furthermore, we set
3.2 Six-joint robot 81
u = a2 + d4 s3 ,
v = d4 c3 . (3.32)
Then we have
x = u cos(θ2 ) + v sin(θ2 ),
z = −u sin(θ2 ) + v cos(θ2 ). (3.33)
You may wish to check equation 3.33 by inserting the definitions for
u, v above into equation 3.33. The result should lead back to equa-
tion 3.30 via the trigonometric addition formulas.
Now we set
p
r = u2 + v2 . (3.34)
For r 6= 0 we set
Then (again if r 6= 0)
u/r = sin(γ),
v/r = cos(γ) (3.37)
and thus
x = r sin(γ)c2 + r cos(γ)s2 ,
z = −r sin(γ)s2 + r cos(γ)c2 . (3.38)
x = r sin(γ + θ2 ),
z = r cos(γ + θ2 ). (3.39)
Now we have
x/r = sin(γ + θ2 ),
z/r = cos(γ + θ2 ) (3.40)
or
Remark 3.1
Remark 3.2
For both θ1 and θ3 we each have two solutions. This gives rise to a
total of four solutions, which correspond to the four configurations of
the robot obtained by choosing values ±1 for each of the two variables
shoulder and elbow.
3.2 Six-joint robot 83
u = a2 + d4 s3 (3.47)
and
v = d4 c3 . (3.48)
Then
3 −1
M6 = 0 M3 · T. (3.51)
Here again, T denotes our target matrix as in equation 3.1, i.e.
nx ox ax px
ny oy ay py
T=
nz oz az pz . (3.52)
0 0 0 1
Let M denote the matrix on the right hand side of equation 3.51.
Notice that we can now assume that M is constant, since the angles
θ1 , θ2 , θ3 are the only variables on the right hand side of equation 3.51,
and we have already solved for these variables.
We evaluate both sides of equation 3.51, by multiplying the matrices.
First, we must compute 3 M6 :
c4 c5 c6 − s4 s6 −c4 c5 s6 − s4 c6 c4 s5 d6 c4 s5
3
s4 c5 c6 + c4 s6 −s4 c5 s6 + c4 c6 s4 s5 d6 s4 s5
M6 =
−s5 c6
. (3.53)
s5 s6 c5 d4 + d6 c5
0 0 0 1
θ5 = atan2(hand · u, v) (3.61)
Here the parameter ’hand’ represents the two solutions of equa-
tion 3.58, and the two possible configurations of the hand (see fig-
ure 3.6). As for the values ’shoulder’ and ’elbow’, ’hand’ is an input
value, to be specified by the user.
We return to equation 3.51.
86 3 Robot Kinematics
We compare the elements (3,1) and (3,2) in this matrix equation, and
conclude that
c4 s5 = m13 ,
s4 s5 = m23 ; (3.62)
c4 = m13 /s5 ,
s4 = m23 /s5 . (3.63)
u4 = −ax s1 + ay c1 ,
v4 = ax c23 c1 + ay c23 s1 − az s23 . (3.67)
Then
3.2 Six-joint robot 87
Solving for θ6
In our matrix in equation 3.53, we see that the same method which led
to the solution for θ4 can also give us θ6 . Specifically, the entries m31
and m32 of the matrix give us:
and obtain:
Remark 3.3
We conclude that solutions for our six-joint robot are unique except
for cases with s5 = 0, as long as we provide input values for ’shoulder’,
’elbow’ and ’hand’. This completes the analysis of the six-joint robot.
3.3 Inverse Solution for the Seven-Joint DLR-Kuka Robot 89
i αi ai di θi
1 -90 0 d1 θ1
2 90 0 0 θ2
3 -90 0 d3 θ3
4 90 0 0 θ4
5 -90 0 d5 θ5
6 90 0 0 θ6
7 0 0 d7 θ7
Remark 3.4
due to the elbow motion (see also figure 3.12). Given a target matrix,
we call the curve corresponding to this target matrix the null-space of
the robot.
Fig. 3.12: The null-space of a seven-joint robot, for elbow-up and elbow-
down configurations.
J4
J3 J5
J2 J6 J7
J1 P = x0 -z0 -plane
J4
J3 J5 J6 J7
J2
J1
Remark 3.5
In practice the angular ranges for revolute joints are limited. Typical
values are ±170 degrees for the R-joints and ±120 degrees for the
P-joints. Taking into account these limitations, the range of the elbow
angle (null-space in figure 3.12) is limited as well. It turns out that
the null-space is fragmented for some positions of the effector, i.e.
the null-space consists of several disconnected pieces. This effect is
illustrated in figure 3.15.
Inverse Solution
We begin with an informal overview of the solution path. We solve
for the first joint variable in much the same way as in the case of the
six-joint robot. That is, we look at the robot from above, and pro-
ject the wrist center point q onto the x-y-plane. As in the case of the
six-joint robot, we could then simply find the first joint angle as in
equation 3.44, i.e.
500
0
500 500
0 0
-500 -500
y0
q
∆
θ1
O x0
Fig. 3.16: Slack value ±∆ for the first joint of the seven-joint robot. The
value of ∆ controls the placement of the elbow.
But this will be not quite the solution for the first joint. Rather, we
allow for some extra ‘slack’ for choosing our first angle value. This is
illustrated in figure 3.16. The slack is denoted by ∆ .
To illustrate the solution path, we compare the seven-joint robot to
the six-joint robot under the R-P-notation introduced above. For the
six-joint robot, we had an RPPRPR structure. To find the inverse solu-
tion, we grouped the structure into two parts, namely RPP − RPR. The
robot’s wrist has an RPR-structure.
If we have assigned a value for the first joint variable (including some
slack value ∆ ), the remainder of the seven-joint robot is an PRP −
RPR structure. This is a six-joint robot, but different from the one we
analyzed above. The wrist is the same, but the robot’s base is a PRP-
3.3 Inverse Solution for the Seven-Joint DLR-Kuka Robot 95
q = p − d7 a (3.78)
where
q = (qx , qy , qz )T (3.79)
is again our wrist center point.
Solving for θ1
Set
Solving for θ2 , θ3 , θ4
Let
1 5
M2 · ... · M6 (3.81)
[4]
96 3 Robot Kinematics
Solving for θ2 , θ3 , θ4
q0x = d5 (c2 c3 s4 + s2 c4 ) + s2 d3 ,
q0y = (s2 c3 s4 − c2 c4 ) − c2 d3 ,
q0z = s3 s4 d5 . (3.88)
k q’ k2 = q02 02 02
x + qy + qz . (3.89)
Evaluating the right hand side yields
k q’ k2 −d32 − d52
c4 = , (3.90)
2d3 d5
which implies
q
s4 = ± 1 − c24 , (3.91)
and finally
θ4 = atan2(s4 , c4 ). (3.92)
Having solved for θ4 , we can use the last line in equation 3.88 to solve
for θ3 . If s4 is not zero, we set
q0z
s3 = (3.93)
s4 d5
and
q
c3 = ± 1 − s23 (3.94)
so that
θ3 = atan2(s3 , c3 ). (3.95)
Otherwise (s4 = 0), we set s3 = 0.
98 3 Robot Kinematics
Finally, for the last angle (θ2 ), a procedure similar to the case of the
six-joint robot applies.
Set
u = c3 s4 d5 ,
v = c4 d5 + d3 . (3.96)
Then
p
r= u2 + v2
v u
γ = atan2 , (3.98)
r r
As for the six-joint robot, the case r = 0 can be excluded. We thus
have
v
= sin(γ),
r
u
= cos(γ). (3.99)
r
Therefore, we can rewrite the expressions for q0x , q0y as
q0x = r cos(γ + θ2 ),
q0y = r sin(γ + θ2 ), (3.101)
3.3 Inverse Solution for the Seven-Joint DLR-Kuka Robot 99
so that
Solving for θ5 , θ6 , θ7
4
M7 = (0 M4 )−1 · T. (3.104)
Evaluating the left side, we find
c5 c6 c7 − s5 s7 −c7 s5 − c5 c6 s7 c5 s6 d7 c5 s6
4
c6 c7 s5 + c5 s7 c5 c7 − c6 s5 s7 s5 s6 d7 s5 s6
M7 =
−c7 s6
. (3.105)
s6 s7 c6 d5 + d7 c6
0 0 0 1
This matrix is of course the same matrix which we obtained for the
wrist of the six-joint arm, the only difference being the names of the
angles.
The right side of equation 3.104 is constant, and we can evaluate it,
since we have already solved for θ1 , ..., θ3 .
Thus, we can define
m13 = 4 M7(1,3) ,
m23 = 4 M7(2,3) ,
m33 = 4 M7(3,3) ,
m31 = 4 M7(3,1) ,
m32 = 4 M7(3,2) . (3.106)
100 3 Robot Kinematics
q0z
s3 = (3.118)
s4 d5
and
q
c3 = 1 − s23 (3.119)
(Otherwise set θ3 = 0.)
Then
where
u = c3 s4 d5
v = d5 c4 + d3 (3.121)
For
−1
M =4 M7 = 0 M4 T, (3.122)
102 3 Robot Kinematics
3.5 C-arm
For the analysis of the six-joint robot and the seven-joint robot we ap-
plied a useful method: We partitioned the joint set into two parts (base
and wrist), and solved the two parts separately. Clearly, this strategy
is not always applicable. An example is the C-arm.
Originally, C-arms were not designed to be robots. They were meant
to be positioned by hand. However, newer C-arms have actuated
joints, and can be regarded as robots. In this section we will look
at geometric methods for analyzing linkages. Thus, we will not use
our standard tools (e.g. inverting DH-matrices), but we will compute
the joint angles directly from geometric constraints. This results in a
simple solution, with several practical advantages. Geometric meth-
ods can be regarded as an alternative to the algebraic methods presen-
ted before.
104 3 Robot Kinematics
O4 a5 O5
a4
O3 θ5
z1
S1 y1 θ4
x1
z0
S0
y0
x0
Fig. 3.18: The first two coordinate systems S0 , S1 and the origins O3 to
O5 for the forward kinematic analysis of a C-arm.
d3
x3
S3 /O3
S2 /O2 z2 θ4
y3 z3
x2
y2
x4
z4
S4 y4
O4 O5
z5
y5
S5 x5
Given the coordinate systems, we can set up the DH-table for the C-
arm (see table 3.2). Recall that the first line in the table describes the
transition from the coordinate system S0 to system S1 (figure 3.18).
i αi ai di θi
1 0 0 d1 0
2 -90 0 0 θ2
3 0 0 d3 -90
4 90 a4 0 θ4
5 90 a5 0 θ5 + 90
From the table, we derive the matrices for the forward transformation.
100 0 c2 0 −s2 0
0 1 0 0 1
0
M1 = M2 = s2 0 c2 0
0 0 1 d1 0 −1 0 0
000 1 0 0 0 1
0 10 0 c4 0 s4 a4 c4
−1 0 0 0 3
2
M3 = M4 = s4 0 −c4 a4 s4
0 0 1 d3 0 1 0 0
0 00 1 0 0 0 1
−s5 0 c5 −a5 s5
4
c5 0 s5 a5 c5
M5 =
0 1 0 0
0 0 0 1
We multiply the matrices just derived for the C-arm and obtain:
0
M5 [1] =
−c2 s4 s5 − s2 c5
−s2 s4 s5 + c2 c5
−c4 s5
0
3.5 C-arm 107
0
M5 [2] =
−c2 c4
−s2 c4
s4
0
0
M5 [3] =
c2 s4 c5 − s2 s5
s2 s4 c5 + c2 s5
c4 c5
0
0
M5 [4] =
−a5 c2 s4 s5 − a5 s2 c5 + a4 c2 s4 − s2 d3
−a5 s2 s4 s5 + a5 c2 c5 + a4 s2 s4 + c2 d3
−a5 c4 s5 + a4 c4 + d1
1
(3.126)
a5
a4
None of these ideas will work here. We cannot decouple rotation from
translation. Also, there are no pairs of equations of the form sin(θi ) =
a, and cos(θi ) = b in equation 3.126.
Furthermore, we do not have a full 4-by-4 matrix to start from. Rather,
our target is given as a direction from which to image, and a point on
the beam’s central axis.
Thus, we assume we only have a point p and a unit vector u with
px ux
p = py
and u = uy , (3.127)
pz uz
rather than the full 4 × 4 target matrix.
Here p is our target position of the point O5 , and u is the axis of the x-
ray beam. Recall that O5 is the origin of S5 . O4 is the center of rotation
of the C-arm’s C. The vector u is the z-vector of the last coordinate
frame. The inputs for the inverse analysis are shown in figure 3.21.
Our goal is to extract the angles and translation values from p and u
alone.
3.5 C-arm 109
z0 p
x0 y0
Fig. 3.21: Target position and orientation for the inverse analysis of the
C-arm.
v
v⊥u u
u⊥E
v||u u u||v v
For the inverse analysis, we start from point O5 . As noted above, our
target is not a homogeneous matrix, but rather the position of O5 and
a direction vector u pointing along this central axis. We assume u is a
unit vector. Also, by p we denoted the coordinate vector for the point
O5 .
We place a line g along the z2 -axis (see figure 3.19). Notice in the
figure, that the axis z2 will always coincide with axis z3 . Our line g
has a direction vector, which we will call v.
The second joint of the C-arm is a revolute joint with axis z1 and joint
angle θ2 . We will now assume that we have already solved for this
angle θ2 . Hence our line g and its direction vector v is fixed from now
on. From figures 3.18 and 3.19, we see that the coordinates of v can
be calculated explicitly, as v = (− sin θ2 , cos θ2 , 0)T .
Our strategy will be the following: we solve for the open variables
d1 , d3 , θ4 , θ5 , based on the assumption that we already have a value
for θ2 . Thus we will have equations for the values of d1 , d3 , θ4 , θ5 as
functions of θ2 .
Solving for θ4
Detector Fv
z0
θ4
u
v
Source
Solving for θ5
We again look at the two vectors u and v. The two vectors span a
plane, which is the mechanical plane of the C-arm’s C. The vector
u⊥v is obtained as u − (uv)v. Then θ5 is the angle between u and u⊥v ,
hence we set
u⊥v
θ5 = sgn(θ5 ) arccos u . (3.133)
||u⊥v ||
θ5
u⊥v
u
v
Solving for d1
pz = −a5 c4 s5 + a4 c4 + d1 . (3.134)
Thus, d1 = pz + a5 c4 s5 − a4 c4 .
3.5 C-arm 113
Solving for d3
The points O3 , O4 and O5 span a plane. This plane contains the vectors
u and v defined above. The situation is illustrated in figure 3.25.
z0
u
O4 a5
O5 = p
a4
v
θ4
O3
O0 y0
d3
Notice that the lengths of the offsets are a4 and a5 , while the vectors
u⊥v and v⊥u point into the directions of these offsets (figure 3.26).
We have
v||u = (uv)u,
v⊥u = v − v||u ,
v − (uv)u. (3.135)
O4 = O5 − ||vv⊥u || a5 (3.136)
⊥u
O3 = O4 − ||uu⊥v || a4 (3.137)
⊥v
114 3 Robot Kinematics
u
z0
θ5
u⊥v
v
O0 y0
z0
θ4
v
O0 y0
Fig. 3.26: The vectors u⊥v and v⊥u point into the directions of the offsets.
Left: rotation about angle θ5 . Right: rotation about angle θ4 .
Thus,
O5 = p, (3.139)
3.5 C-arm 115
d3 = ||O2 − O3 || (3.140)
and
O2 = (0, 0, d1 )T . (3.141)
Solving for θ2
3.5.3 Applications
Given the forward and inverse kinematic solution for the C-arm, we
can obtain several interesting applications of robotic C-arms.
1. image stitching
2. intraoperative 3D reconstruction and CT imaging
3. 4D imaging
4. cartesian control
In the first application, we paste several images to obtain a panorama
image, e.g. of a long bone or full leg (see figure 3.27). For the second
application, we move joint five, to obtain images from a large number
116 3 Robot Kinematics
Example 3.1
Laser beam
z0 u
v⊥u
a5
u⊥v a4
h
v
d3
d1
y
O0 xp 0
y p
0
Fig. 3.29: Vector chain for the linear equation system in equation 3.142.
From the figure, we see that a chain of vectors leads from the origin
to the point q = (x p , y p , 0)T . In this chain there are three variables:
d1 , d3 and λ (scale factor of the direction vector for h). Thus, if h
is in general position, then d1 and d3 will occur as solutions of the
following linear equation system:
0 vx ux xp
u⊥v v⊥u
d1 0 + d3 vy + a4 + a5 + λ uy = y p
||u⊥v || ||v⊥u ||
1 vz uz 0
(3.142)
Exercise 3.6 at the end of this chapter discusses a short program for
solving the four-axes problem for C-arms.
We saw that the methods for kinematic analysis not only apply for
revolute joints, but also for translational joints. The C-arm is one such
example, as it consists of several translational and revolute joints.
Many medical devices are based on a so-called center-of-arc kin-
ematic construction. This means that a rotation about the patient’s
head or body is the basic motion type built into the system. The ef-
fector is thus mounted to an arc in space. This arc is frequently moun-
ted to a base, which itself can be rotated. Examples for devices follow-
ing this general convention are: neurosurgical frames for stereotaxis,
angiography systems and radiation therapy machines.
Figure 3.30 shows an example of a particularly simple center-of-arc
construction. This robotic construction is used in functional neurosur-
gery. There are four translational axes (t1 ,t2 ,t3 ,t4 ) and two revolute
axes r1 , r2 . The two revolute axes intersect in space.
It turns out that the inverse analysis for the center-of-arc mechanism
in figure 3.30 is much simpler than the analysis of the C-arm and the
six-joint robot, and is left as an exercise.
120 3 Robot Kinematics
and check whether the robot can reach p under this orientation. The
last step can be implemented by a routine for inverse kinematics.
An alternative is to require a full sphere be reachable at p, and then
plot all points p in the work space, admissible under this measure. We
can then compare different constructions, and we will do this in the
following example.
Example 3.2
i αi ai di joint range
1 -90 0 340 -170...+170
2 90 0 0 -120...+120
3 -90 0 400 -170...+170
4 90 0 0 -120...+120
5 -90 0 400 -170...+170
6 90 0 0 -120...+120
7 0 0 111 -170...+170
Table 3.3: Link parameters (in mm) and joint ranges the seven-joint
robot in the example
Figures 3.31-3.32 show the results for the six-joint version compared
to the seven-joint version. As expected, the seven-joint robot does
have a larger range of points reachable in arbitrary orientation. But
what is the influence of the joint ranges? In other words, can we in-
122 3 Robot Kinematics
Fig. 3.31: Dexterity volume for the six-joint version of the seven-joint ro-
bot. Grid points marked in black are reachable in arbitrary orientation
(full sphere). The grid points p have a grid distance of 100 mm.
crease the joint ranges of the six-joint version in such a way that it
matches (or outperforms) the seven-joint version?
To answer this second question, we increase the joint ranges of the
six-joint version from
to
and recompute the dexterity map. The result is shown in figure 3.33.
Fig. 3.33: Dexterity volume for the six-joint robot, with joint ranges en-
hanced to −150◦ , . . . , +150◦ for the P-joints.
z0
x0
Fig. 3.34: Six-joint Dextor robot with unlimited joint ranges. Self-
collisions are mechanically excluded.
z0
x0
Exercises
Exercise 3.1
Verify the results in equations 3.5, 3.14 and 3.17 for the settings θ1 =
0, θ2 = 0, θ3 = 0, and θ1 = π/2, θ2 = 0, θ3 = −π/2, and θ1 = 0, θ2 =
−π/4, θ3 = −π/4.
Exercise 3.2
Show that the value r defined in equation 3.34 is not zero, unless a2 =
d4 and sin θ3 = 0.
Exercise 3.4
Implement the kinematic solution for the six-joint robot. Inputs: tar-
get matrix T and configuration parameters shoulder, elbow and hand.
Extract angles θ1 , ..., θ6 , and resubstitute these values into the forward
kinematic equation, giving a 4×4 matrix. Compare to the input matrix
T.
Exercise 3.5
Exercise 3.6
Summary
Notes
References
θ2
l1 l2
y
θ1 x p
This robot has two revolute joints. How do we move the joints, in or-
der to produce a motion of the point p along the x-axis? As a first
guess, assume we move the two joints with constant and equal velo-
cities, ignoring all necessary accelerations. The two links have equal
lengths (i.e. l1 = l2 ), and we hope to produce our linear motion in this
way. Unfortunately, a quick simulation will show that this will not
work. Instead, our tool will move along a circular arc. If we modify
131
132 4 Joint Velocities and Jacobi-Matrices
our input parameters, the arc will become bigger or smaller, or it will
be deformed. But, alas, it will be difficult to deform the arc into a line.
What now, if we not only need a line motion, but we also require
that the velocity of the tool remains constant throughout the motion?
Hence, even if we succeed to produce a line, we may not traverse
this line with a very even velocity. As an example, we again look at
figure 4.1. Assume, the tool is very near the furthest point reachable
by the linkage along the x-axis. As we move the joints to stretch out
the arm even further, the tool will move more and more slowly, until
the velocity will become zero, for the arm fully stretched.
We started the above discussion for the case of a line motion. How-
ever, there is an even simpler case for velocities. This most simple
case arises when we require to move our joints, but to keep the tool in
a fixed position. It will be discussed in the first section of this chapter.
4.1 C-arm
O4 a5 O5
a4
O3
z0
y0
Fig. 4.2: The point O4 is the center of rotation of the C-arm’s C. a4 and
a5 are the offsets. To obtain a CT image from C-arm x-ray images, we
rotate the C in the drawing plane.
O5
Fig. 4.3: Rotating the C-arm’s C, while keeping the target point O5 static
requires compensating motion with the two prismatic joints. The motion
of the two prismatic joints is visualized by the two arrows.
134 4 Joint Velocities and Jacobi-Matrices
within the interval [0, 2π]. It may seem odd to choose this interval for
a prismatic joint, but the reason will soon become clear.
We insert the functions of t thus defined into the forward kinematic
matrix, and obtain a matrix, representing tool orientation and tool po-
sition as a function of time.
Can we state explicit time functions for d1 (t) and d3 (t) which will
keep the point O5 static during orbital rotation of the C-arm’s C?
First, note that in practice our orbital rotation should have constant
velocity. Hence, the derivative θ50 (t) should be constant, or, to be more
specific, we set
θ5 (t) = t (4.1)
Here θ5 runs from 0 to 2π, and this is why we chose the same range
for t. A full rotation of 2π is unrealistic for mechanical reasons, but
we will ignore this for now.
z0
a5
θ a5 sinθ
a4
a5 cosθ
d1
y0
d3
Fig. 4.4: Deriving explicit joint velocity functions d1 (t) and d3 (t).
Figure 4.4 shows how to obtain explicit joint motion functions for
d1 (t) and d3 (t). We see that the displacement required for d3 amounts
to −a5 cos(θ ), where θ is the angle marked in the figure. Likewise,
we find that the displacement for d1 over time is a5 sin(θ ).
4.1 C-arm 135
Subtracting the kinematic offset a4 (see figure 4.2), we obtain the fol-
lowing motion functions
d1 (t) = a5 sin(t) − a4 + z p
d3 (t) = −a5 cos(t) + y p (4.2)
Here the constants y p and z p denote the target position for O5 . Having
found the motion functions for the joints, we can obtain the velocities
for the prismatic joints by taking derivatives.
By d˙1 we denote the joint velocity d10 (t). With this notation, we obtain
d˙1 = a5 cos(t)
d˙3 = a5 sin(t) (4.3)
In practice, explicit joint motion functions are not always needed. In-
stead, we can discretize the required motion. For example, to move
136 4 Joint Velocities and Jacobi-Matrices
Example 4.1
4.2 Jacobi-Matrices
side is given by the constant matrix in equation 4.7 and the right side
is a matrix containing the six angles θ1 , ..., θ6 as variables.
This matrix equation consists of 12 individual equations, since we
have 4 × 4 matrices, where we can ignore the last matrix line. Each
equation has a constant on the left side, and variables θi on the right
side.
We take partial derivatives of these 12 equations with respect to the 6
variables. This works as follows: assume the current position of our
robot is known and given both in terms of the joint variables and in
terms of the tool matrix. When taking the partial derivative of one
such equation with respect to, say, θ1 , we insert the constant values
for all other joint variables θ2 , ..., θ6 into the equation, and obtain an
equation of a single variable (here variable θ1 ). We can now take the
derivative of this single-variable equation. To see why we have the
joint angle values for all joints, at the current position, we look at the
following argument: We step through the sequence of points p1 , ..., pn ,
and we can assume that we have the joint angles corresponding to pi ,
when we compute the joint angles for pi+1 . This argument holds for
p1 , since here we can read the encoders of our joints.
For the six-joint robot, we obtain a matrix of the form
∂ m11 ∂ m11
∂ θ ... ∂ θ6
1
∂ m21 ∂ m21
∂ θ ... ∂ θ
J= 1 6 . (4.9)
..
.
∂ m34 ∂ m34
∂ θ1 ... ∂ θ6
J is called the Jacobi-matrix for our robot, or the Jacobian, for short.
We can now set up an equation, relating joint increments to path in-
crements:
dnx
dny dθ1
..
.. = J . . (4.10)
.
dθ6
d pz
Notice that the equation gives an approximation. Thus, we require the
increments to be small.
4.2 Jacobi-Matrices 139
(n o a) (4.12)
is the current (known) orientation.
We assume our angular increments dα, dβ , dγ refer to the x-, y- and
z-axes respectively. Hence, we can set up a matrix with these incre-
ments.
The matrix is obtained by multiplying the elementary rotations, i.e.
140 4 Joint Velocities and Jacobi-Matrices
sdα = dα (4.14)
for small angles. Furthermore, we set any products of two small angles
to zero.
Then, after multiplying, we obtain the approximation matrix D, ap-
proximating the matrix in equation 4.13, where
1 −dγ dβ
D = dγ 1 −dα . (4.15)
−dβ dα 1
D · (n o a) = (n + dn o + do a + da). (4.16)
We rewrite this equation:
0 −dγ dβ nx ox ax dnx dox dax
dγ 0 −dα ny oy ay = dny doy day . (4.18)
−dβ dα 0 nz oz az dnz doz daz
4.2 Jacobi-Matrices 141
All matrix entries on the left side of equation 4.18 are known. We
can thus compute the unknowns on the right hand side, simply by
multiplying the two matrices on the left side.
Remark 4.1
The matrix
0 −dγ dβ
S = dγ 0 −dα (4.19)
−dβ dα 0
is obtained as S = D − I, and has a special property, namely
S + ST = 0 (4.20)
Matrices with this property are called skew-symmetric. Notice that S
is not orthogonal, and does not represent a rotation.
After having found all constant values, equation 4.10 is a linear equa-
tion system, and we can solve for the variables dθ1 , ..., dθ6 .
Remark 4.2
When looking at the six-joint robot in figure 3.1 or even the eight-
joint robot in figure 3.17 we see that the most important element is
the recurrent two-link sub-assembly with parallel joint axes. We have
discussed this two-link sub-assembly in several examples throughout
the preceding chapters and also in the introduction of this chapter (fig-
ure 4.1).
Why is this element the most important element? We have already
observed that our six-joint robot consists of two three-joint sub-
assemblies. In the kinematic analysis, we separated these two sub-
assemblies.
We now look at the first of these sub-assemblies. It is the well-known
three-joint robot from chapter 2. Here, the base joint determines a
coarse pre-positioning in the work space. Recall that the kinematic
analysis for this joint was particularly simple, and we were able to
decouple it from the entire rest of the robot. Now the next two joints
(joints 2 and 3 for the six-joint robot) do the actual manipulation at the
target, while all remaining joints only provide hand orientation. Thus,
roughly speaking, joints 2 and 3 provide the main elements of the
robot, when it comes to handling and manipulation. We give a name
to the corresponding sub-assembly, and refer to it as the (planar) two-
link manipulator.
We compute the forward matrix for a planar two-link manipulator. We
obtain
c12 −s12 0 l2 c12 + l1 c1
0
s12 c12 0 l2 s12 + l1 s1
M2 =
0 0 1
.
(4.21)
0
0 0 0 1
Notice that here we refer to the two joints of our two-link manipulator
as joint 1 and joint 2, ignoring the fact that the two-link manipulator
may be part of a larger robot with more joints.
4.3 Jacobi-matrices and velocity functions 143
Our goal is to find explicit velocity functions for the two-link ma-
nipulator. We have already done this for the case of a C-arm, in the
introduction of this chapter, but this was a particularly simple case,
and we did not need any dedicated methods.
We consider the elements (1,4) and (2,4) of the matrix in equa-
tion 4.21. These two elements give the x- and y-positions of our tool.
Hence, we have
x = l2 c12 + l1 c1 ,
y = l2 s12 + l1 s1 . (4.22)
Here θ1 (t) and θ2 (t) are two separate functions, describing the mo-
tions of the two joints over time.
We differentiate these two functions with respect to the time parameter
t.
This gives
x0 (t) = −l1 sin(θ1 (t)) · θ10 (t) − l2 sin(θ1 (t) + θ2 (t)) · (θ10 (t) + θ20 (t)),
y0 (t) = l1 cos(θ1 (t)) · θ10 (t) + l2 cos(θ1 (t) + θ2 (t)) · (θ10 (t) + θ20 (t)).
(4.24)
We recall the notation from equation 4.3. With this notation (replace
θi0 (t) by θ̇i ) and the familiar short forms si = sin(θi ) we rewrite equa-
tion 4.24 to
x θ̇
p= , θ̇ = 1 . (4.26)
y θ̇2
−l1 s1 − l2 s12 − l2 s12
J= . (4.27)
l1 c1 + l2 c12 l2 c12
ṗ = Jθ̇ . (4.28)
You may wish to take a minute to multiply the product Jθ̇ and see that
this will indeed result in the right hand side of equation 4.25.
The matrix J is a 2-by-2 matrix, and has an interesting property: we
can invert it! Thus we can write
Now we can derive explicit velocity functions for the joints of our
two-link manipulator.
Suppose we wish to move the tool along the x-axis. Thus we set the
velocity along the x-axis to 1 and the velocity along the y-axis to 0.
Then
1
ṗ = . (4.33)
0
Multiplying, we obtain
c12
l1 s2
θ̇ = J ṗ =
−1
. (4.34)
− l2cs12 − lc112s2
This shows that the velocity will tend to infinity as our angle θ2 tends
to 0, i.e. as the arm stretches.
In equation 4.34 we have explicitly stated the non-linear joint velocity
functions which will move the tip along a straight line.
Remark 4.3
From equation 4.31 we see that we cannot invert the matrix A in equa-
tion 4.30 if ad − bc = 0. In this case the inverse A−1 does not exist,
and A is called singular. We apply this to equation 4.32. Here, J has
no inverse, whenever s2 = 0. This means that J is singular if θ2 = 0
or θ2 = π/2. Looking at the robot, we see that such a configuration
occurs, when the two links are aligned, i.e. the arm is stretched.
ω = θ̇ a (4.36)
Here θ̇ is a scalar function of time. Thus, the angular velocity ω is
simply a vector, the length of which gives the magnitude of the velo-
4.4 Geometric Jacobi-matrix 147
city, and the direction gives the axis of the rotation (here unit vector
a).
Now suppose we have a point p which is given in the coordinates of
the frame S, and rotates with S. We define the linear velocity of p as
v = ω × p. (4.37)
This definition may seem odd, and we visualize the effects of the latter
definition for the case of a small example.
θ1
z0
p1 v
x0 -y0 -plane
Example 4.2
Figure 4.5 shows a one-joint robot, which moves in the x0 -y0 -plane.
The robot rotates about the z0 -axis. The joint angle is θ1 . The position
of the end effector is given by the point p1 in the center of the grip-
per. Now we move the only joint. The angular velocity is the scaling
factor θ̇1 of the vector z0 . The linear velocity of p1 is v. This linear
velocity v is the cross-product of θ̇1 z0 and p1 . We see that an increase
of the link length of the robot, or an increase of θ̇1 will increase the
linear velocity. By contrast, an increase in link length will not increase
angular velocity.
Thus, our linear velocity is in fact a linear velocity resulting from an
angular velocity, and it is an instantaneous velocity. Notice that the
vector v in figure 4.5 changes direction as the robot moves.
In figure 4.6 we see a planar robot with several revolute joints, again
moving in the x0 -y0 -plane. Suppose now only a single joint of this
148 4 Joint Velocities and Jacobi-Matrices
pi−1 /zi−1
v = ω × (pn − pi−1 )
θi
y0 pn
z0 x0
Fig. 4.6: Linear velocity v of the point pn , resulting from the motion
of a single joint (here joint i, with axis zi−1 ). The axes z0 and zi−1 are
orthogonal to the drawing plane, and are indicated by circles with dots.
We should also note that the angular velocity, according to the above
definition, is a property of a frame, while the linear velocity is a prop-
erty of a point.
4.4 Geometric Jacobi-matrix 149
ω = ω0 + 0 R1 (θ )ω1 . (4.40)
Then ω is again a vector, given with respect to the base frame, so that
we can iterate the definition, for more (revolute) joints.
Here we have simply defined the composite angular velocity by a vec-
tor addition. We must now confirm that this definition is in agree-
ment with our intuitive understanding of a combined rotation of two
coordinate frames. This will be done in the following example. The
example will also illustrate the application of angular velocities.
Example 4.3
We are given a robot with four revolute joints. The four joints form a
kinematic chain, as in typical robots. Now suppose the four rotations
are the following:
1. Rotation about the z-axis of the base frame S0 (angle α)
2. Rotation about the x-axis of the new coordinate system S1 (angle
β)
3. Rotation about the x-axis of S2 (angle −β )
4. Rotation about the z-axis of S3 (angle −α)
150 4 Joint Velocities and Jacobi-Matrices
α(t) = t,
β (t) = t. (4.41)
Hence,
−α(t) = −t,
−β (t) = −t. (4.42)
Then we have
0 1
ω0 = 0 ,
ω1 = 0 ,
(4.43)
1 0
and
−1 0
ω2 = 0 , ω3 = 0 . (4.44)
0 −1
From chapter 3 we know that the combined rotation is simply given
by the matrix product
100
R(z, α(t))R(x, β (t))R(x, −β (t))R(z, −α(t))) = 0 1 0 (4.46)
001
i.e. a zero rotation. But we are not supposed to know this yet. Instead
we would like to apply the composition rule for angular velocities. We
should then obtain a total angular velocity of zero.
Recall that we must express all angular velocities in the common base
frame. Now ω0 is already given in base coordinates.
To express ω1 in base coordinates, we simply multiply with the matrix
4.4 Geometric Jacobi-matrix 151
0
R1 = R(z, α(t)). (4.47)
Thus, the angular velocity of S2 with respect to S1 , given in base co-
ordinates is
With the definitions of angular and linear velocity, we can specify our
goals with respect to the geometric Jacobian. Our goal is thus to find
two 3 × n matrices Jv and Jω , such that the linear velocity v and the
angular velocity ω of the end effector can be expressed in terms of the
joint velocity vector θ̇ as
v = Jv θ̇ ,
ω = Jω θ̇ . (4.52)
152 4 Joint Velocities and Jacobi-Matrices
The two parts Jv and Jω of the geometric Jacobian will then be as-
sembled into a single matrix
Jv
J= . (4.53)
Jω
We then have
v
= Jθ̇ . (4.54)
ω
This is a direct extension of equation 4.28. Notice again that in this
equation, both the angular and the linear velocity (of our effector) are
given in terms of the base frame.
Theorem 4.1:
Jv z0 × (pn − p0 ) z1 × (pn − p1 ) . . . zn−1 × (pn − pn−1 )
J= = .
Jω z0 z1 ... zn−1
(4.55)
Here again, pi denotes the position of the origin of the i-th DH-frame,
and can be read from the last column of the DH-matrix 0 Mi (see
chapter 3). Likewise, zi is short for
0
0
Ri 0 , (4.56)
1
where 0 Ri is the 3 × 3 rotational part of 0 Mi . Clearly, both pi and zi
can be computed with the forward analysis methods from chapter 3.
Notice that the theorem requires robots with revolute joints. In case
joint i is a prismatic joint, we replace the column
zi−1 × (pn − pi−1 )
(4.57)
zi−1
by the column
4.4 Geometric Jacobi-matrix 153
zi−1
. (4.58)
0
Before we prove the theorem, we will look at two examples.
Example 4.4
We begin with the one-joint robot from figure 4.5. Here, we have
c1
pn = p1 = s1 .
(4.59)
0
Also,
0
z0 = 0 . (4.60)
1
Thus
0 c1 −s1
z0 × (p1 − p0 ) = z0 × p1 = 0 × s1 = c1 . (4.61)
1 0 0
This gives Jv1 .
Moving joint 1 with a magnitude of θ̇1 will hence result in a linear
velocity of p1 of
−s1
θ̇1 c1 . (4.62)
0
We see that this is consistent with the linear velocity v in our defini-
tion, which we had looked at in example 4.2 (see equation 4.39).
For Jω1 , we have
0
Jω1 = z0 = 0 .
(4.63)
1
Thus the angular velocity (again joint 1 moving at magnitude θ̇1 ) is
154 4 Joint Velocities and Jacobi-Matrices
0
0 . (4.64)
θ̇1
Therefore the full geometric Jacobian is
−s1
c1
Jv 0
= 0 .
(4.65)
Jω
0
1
Very similarly, for the two-link manipulator (see figure 4.1), we obtain
−l1 s1 − l2 s12 −l2 s12
l1 c1 + l2 c12 +l2 c12
Jv 0 0
J= =
. (4.66)
Jω 0 0
0 0
1 1
We compare to the matrix derived in equation 4.27. Clearly, the first
two lines are the same. The third line in equation 4.66 has only zeros.
This is simply due to the fact that we now look at the Jacobi-matrix for
a full three-dimensional robot, and our planar two-link manipulator
can never leave the plane. Finally, the three lines at the bottom of
equation 4.66 describe angular velocity, and we see that we add the
angular velocity from the two joints.
Finally, we look at an extension of the two-link manipulator. We add
a third revolute joint at the end of the two-link manipulator. The entire
mechanism can only move in the plane. The length of the third link is
l3 .
In this case, we obtain the geometric Jacobi-matrix
4.4 Geometric Jacobi-matrix 155
−l1 s1 − l2 s12 − l3 s123 −l2 s12 − l3 s123 −l3 s123
l1 c1 + l2 c12 + l3 c123 l2 c12 + l3 c123 l3 c123
0 0 0
J= . (4.67)
0 0 0
0 0 0
1 1 1
The matrix in equation 4.55 consists of two lines. We will first con-
sider the second line.
We will see that it follows directly from the composition rule for an-
gular velocities.
We must show (compare equation 4.54) that
θ̇1
..
ω = Jω . . (4.68)
˙
θn
Given the definition of J in equation 4.55, we have
Jω = z0 z1 . . . zn−1 , (4.69)
and we must show that
v1 = ω1 × p, (4.77)
v2 = ω2 × p. (4.78)
(4.79)
v1 + v2 = ω1 × p + ω2 × p (4.80)
By linearity, the right hand side of the last equation reduces to
ω1 × p + ω2 × p = (ω1 + ω2 ) × p. (4.81)
We see that the vector v1 + v2 also represents a linear velocity, act-
ing on the same point p. But this new linear velocity stems from an
angular velocity ω, for which we have
ω1 + ω2 = ω. (4.82)
Similar to the composition rule for angular velocities, we thus obtain
a composition rule for linear velocities (stemming from angular velo-
cities), which holds whenever the vectors are defined within the same
base frame. But now the desired result in equation 4.76 follows from
the fact that all vectors zi and pi are given with respect to the base
frame.
(End of Proof)
Example 4.5
that the center point of the beam axis (point p) remains at the same
position. We saw that we must compensate with the two prismatic
joints.
The derivation in equation 4.3 was elementary, and we did not use
Jacobi-matrices.
We will now set up the geometric Jacobian for CT-imaging with the
C-arm. We can simplify the construction. As discussed in the intro-
duction of this chapter, only three joints are needed in this case. Thus,
we set the joint angles θ2 and θ4 to zero. Figure 4.7 shows the null
position and the coordinate systems.
x2 z3 y3
z2
a5
y2 p x3
a4
d3
z1
x1
y1
d1
z0
y0
x0
Fig. 4.7: Null position and coordinate systems for CT imaging with the
C-arm.
i αi ai di θi
1 −90 0 d1 0
2 90 a4 d3 −90
3 90 a5 0 θ5 + 90
From the DH-parameters, we obtain the matrices
1 0 0 0
0
0 0 1 0
M1 =
0 −1 0 d1 , (4.83)
0 0 0 1
0 0 −1 0
1
−1 0 0 −a4
M2 = 0 1 0 d3 ,
(4.84)
0 0 0 1
−s5 0 c5 −a5 s5
2
c5 0 s5 a5 c5
M3 = 0 1 0 0 .
(4.85)
0 0 0 1
Multiplying, we find the following:
0 0 −1 0
0
0 1 0 d3
M1 · 1 M2 =
1 0 0 a 4 + d 1 (4.86)
00 0 1
0 −1 0 0
0
c5 0 s5 a5 c5 + d3
M1 · 1 M2 · 2 M3 =
−s5 0 c5 −a5 s5 + a4 + d1 (4.87)
0 0 0 1
Hence,
0 0
p0 = 0 ,
p1 = 0
(4.88)
0 d1
and
160 4 Joint Velocities and Jacobi-Matrices
0 0
p2 = d3 , p3 = a5 c5 + d3 . (4.89)
a4 + d1 −a5 s5 + a4 + d1
Finally,
0 0
z0 = 0 , z1 = 1
(4.90)
1 0
and
−1 0
z2 = 0 , z3 = s5 . (4.91)
0 c5
The first two joints are prismatic. Thus
Jv = z0 z1 z2 × (p3 − p2 ) . (4.92)
Evaluating, we obtain
00 0
Jv = 0 1 −a5 s5 . (4.93)
1 0 −a5 c5
We have
ṗ = Jv θ̇ . (4.94)
Since we require p to remain fixed, we have
0
ṗ = 0 . (4.95)
0
Inserting this into equation 4.94, we obtain
0 00 0 d˙1
0 = 0 1 −a5 s5 d˙3 . (4.96)
0 1 0 −a5 c5 θ˙5
Furthermore, we move the orbital joint (joint 5) with constant velocity
during CT-imaging. Therefore, we can assume that
4.5 Singularities and Dexterity 161
d˙1 − a5 cos(t) = 0,
d˙3 − a5 sin(t) = 0 (4.98)
or that
d˙1 = a5 cos(t),
d˙3 = a5 sin(t). (4.99)
We have thus confirmed the result in equation 4.3 using the geomet-
ric Jacobian. The approach taken in equation 4.3 was elementary. We
see that the geometric Jacobian is a general tool for finding velocity
functions. We will see other applications of the geometric Jacobian
in later chapters. Specifically, the geometric Jacobian will be helpful
when analyzing forces and joint torques.
Example 4.6
Exercises
Insert the functions from equations 4.1 and 4.2 into the forward 4 × 4-
matrix for the C-arm, and verify that the target point O5 remains fixed
while moving t from 0 to 2π.
Write a short program for verifying the result in equation 4.34. Your
result should look like the image in figure 4.8.
0.03
0.02
0.01
−0.01
−0.02
−0.03
−0.04
−0.05
Fig. 4.8: Joint velocity (for joint 2) for moving the tool tip of a two-link
manipulator along a line segment, see equation 4.34.
Consider the two-link manipulator in figure 4.1. Assume the two link
lengths are equal to 1. For this robot, the position of the tool tip is
given by the equation
px c1 + c12
= . (4.101)
py s1 + s12
164 4 Joint Velocities and Jacobi-Matrices
dx = −dθ1 ,
dy = dθ1 + dθ2 . (4.105)
Solving this for dθ1 and dθ2 , we obtain
dθ1 = −dx,
dθ2 = dx + dy. (4.106)
Visualize this result by moving the robot tip to four positions in the
vicinity of the current tip position.
The C-arm has 5 joints. We fix two of these joints, namely joints 2
and 4, and set the joint angles for these two joints to zero. Then the re-
maining mechanism has 3 joints. For non-vertical imaging, we include
4.5 Singularities and Dexterity 165
l2
p
l1
y
Fig. 4.9: Two link manipulator for position θ1 = π/2 and θ2 = −π/2
the angulation axis, i.e. joint 4. Set up the upper part of the geometric
Jacobian matrix Jv for arbitrary non-vertical positions. Derive the ex-
plicit velocity functions for non-vertical CT-acquisition based on the
Jacobian. Compare to the result from section 4.1.
166 4 Joint Velocities and Jacobi-Matrices
p = (a1 , ..., a6 )T ,
q = (b1 , ..., b6 )T .
(4.107)
We look at the path of a single angle, e.g. the first one. Then we must
move θ = θ1 , i.e. the first angle, from a1 to b1 .
Furthermore, we can constrain the velocity of the joint during the mo-
tion. It is useful to require θ 0 (0) = 0 and θ 0 (1) = 0, i.e. the velocity
at the two end points of the motion should be zero. We thus have four
constraints on the trajectory (for the first joint).
Thus
θ (0) = a1 ,
θ (1) = b1 ,
θ 0 (0) = 0,
θ 0 (1) = 0.
(4.108)
4.5 Singularities and Dexterity 167
Summary
z0 × (pn − p0 ) z1 × (pn − p1 ) . . . zn−1 × (pn − pn−1 )
J= .
z0 z1 ... zn−1
(4.110)
With the Jacobi-matrix, we can express the angular and linear velo-
cities of the end effector as functions of the joint velocities. In appro-
priate cases, the Jacobi-matrix is invertible. By inverting the Jacobi-
matrix, we can do the opposite, namely derive explicit functions for
the joint motions needed to produce a given velocity at the effector.
The Jacobi-matrix is useful for many other applications. Thus, we can
characterize singular configurations of the robot. As a further applic-
168 4 Joint Velocities and Jacobi-Matrices
Notes
The text book by Spong, Hutchinson and Vidyasagar [1] gives a de-
tailed analytical derivation of the geometric Jacobi-matrix (including
prismatic joints), together with an analysis of singularities and ma-
nipulability for many robot constructions. This includes an analytical
proof for the composition rule for angular velocity, based on algebraic
properties of skew-symmetric matrices. A similar approach is taken in
the text book by John Craig [2]. Our description summarizes the res-
ults from [1] and [2].
References
Output: θ, Dx, Dy
169
170 5 Navigation and Registration
Fig. 5.2: Registration of a femur head with C-arm x-ray images (2D-3D
registration)
CT image stack
target
Fig. 5.4: Stereo x-ray imaging. Two fixed x-ray sources and two fixed
detectors.
Assume we have two clouds of points in space. Call the clouds A and
B. Our goal is to bring the two clouds into registration. This means
we must move cloud B in such a way that the two clouds match. We
allow the point cloud B to move under rotation and translation, but we
do not allow B to deform. Thus our goal is to compute a translation
vector t and a 3 × 3 orientation matrix R such that the set of points
Rb + t with points b from B will match A. We will also call B the
floating point set, and A the fixed point set.
For simplicity we will first assume that the number of points in both
clouds is the same, and that we know in advance which point b from
B should be matched to which point a from A. With other words,
we know the correspondence between the two sets. We thus assume
the points in A and B are ordered in such a way that the corres-
pondences are simply given by the indexing, i.e. A = (a1 , . . . , an ) and
B = (b1 , . . . , bn ), where ai should by matched to bi .
174 5 Navigation and Registration
For small angles, we can replace the expressions sin θ and cos θ (oc-
curing in R) by linear approximations: Set cos θ to 1 and sin θ to θ .
We thus replace R from equation 5.2 by
1 −θ
. (5.4)
θ 1
The function f is now a quadratic function in the three parameters tx ,ty
and θ . Notice that the coordinates of the points ai , bi are constants
here.
Taking the derivatives with respect to tx , then with respect to ty and
with respect to θ gives three linear expressions. We set these three lin-
ear expressions to zero and obtain a system of three linear equations
with three variables. The solution of this system gives an approxima-
tion for the desired parameters tx ,ty and θ .
In the three-dimensional case, we can use the same method with one
minor difference. Here the matrix R depends on three angles α, β , γ.
We can linearize this matrix in the following way. Replace R by
5.2 Points and Landmarks 175
1 −γ β
D = γ 1 −α . (5.5)
−β α 1
We see that the matrix D is the matrix from equation 4.15. Thus, recall
that D − I is a skew-symmetric matrix. The derivation of the matrix D
is not difficult, and is given in chapter 4.
We again obtain a linear system of equations after taking derivatives.
In this case we have a linear system of six equations in the six vari-
ables tx ,ty ,tz , α, β , γ.
Example 5.1
+ o
+
+ o
+ + o
+ o
+ +
+ + + + o o
+ + o o
+ + o o oo
+ ++ + ++
+ + + oo
+
+ + +
+ +
+ + + o o o
+ + o
Here the points are on the head surface, automatically extracted from
an MR data set. A second cloud of points is acquired later, and
our goal is to register the two point sets. In this case, an additional
procedure must iteratively compute the best set of correspondences.
Only heuristic methods are known for this problem, since the number
of possible pairwise correspondences for two sets of n points is n!.
5.3 Contour-based Registration 177
To this end we want to place the surface model in such a way that the
2D bone contours match optimally.
Contour-based methods need contours, and hence segmented im-
ages to start from. Let us first assume the entire problem was two-
dimensional. This 2D version of the problem does not have much
practical value. But the discussion can be simplified quite substan-
tially. Thus, we assume we have a 2D object and take images of this
object from two camera positions C and C0 (figure 5.7).
In the 2D case the projected object contours each consist of two points
and there are four tangents la , lb , la0 , lb0 to the object. We can determine
these four tangents, since C and C0 are known (figure 5.7).
We manually place the bone model in an approximate position over
the x-ray images. This gives us an initial estimate for the bone posi-
tion. Of course, this initial placement is not optimal and needs to be
refined. In many cases, CT images are taken in standard positions, and
x-ray images are also taken from standard angles, so that this estimate
for the pose of the model is already known, and the step of manually
placing the model is not needed.
la lb
la0
C0
lb0
After having obtained the initial estimate, we can extract four points
a, b, a0 , b0 on the bone model in the following way: We extract a
second set of four tangents, but now these tangents touch the bone
model (in the estimated pose).
5.3 Contour-based Registration 179
(a − u − λ (v − u))(v − u) = 0. (5.8)
The point on the line la which is closest to the point a thus has the
λ -value
180 5 Navigation and Registration
(a − u)(v − u)
λ0 = . (5.9)
(v − u)2
Thus, for λ = λ0 we obtain a closed formula for the quadratic distance
d(a, la )2 . We insert the values for λ0 into the above function f . This
way, we can resolve the distance expressions:
la la
a
Fig. 5.9: Small rotations may move tangent points by large amounts.
5.4 Intensity-based Registration 181
Example 5.2
1 1 3 3 3 3 3 3
1 1 3 3 3 3 3 3
2 2 3 3 1 1 2 2
2 2 3 3 1 1 2 2
Fig. 5.10: Abstract example images. 16 pixels, only three potential gray
values.
However, images from two different sources might not assign the
same gray levels to a given structure. As an example, look at the im-
ages from figure 1.6, showing two images of the same structure (femur
bone) imaged with CT (left) and MR (right).
Since CT is x-ray-based, the bone appears white in the image, whereas
it is black in the MR image. Then it would seem hard to match the
structures, based on the pixel gray levels.
CT MR
1 1 3 3 2 2 2 2
1 1 3 3 2 2 2 2
2 2 3 3 3 3 1 1
2 2 3 3 3 3 1 1
Look at the example in figure 5.11. Which should be the correct re-
gistration position of the images from figure 5.11? The answer is that
the MR image should again be rotated clockwise by 90 degrees. But
now, we must additionally map the pixel intensities according to the
following table:
5.4 Intensity-based Registration 183
CT MR
1 →3
2 →1
3 →2
pA,B (a, b)
I(A, B) = ∑ pA,B (a, b) log . (5.11)
a,b pA (a)pB (b)
The sum in this definition ranges over all outcomes a, b of the two
random variables A, B. It can be shown that I(A, B) ≥ 0 [1].
The connection between mutual information and independence be-
comes visible, if we observe the following: If A, B are independent,
then
and
pA,B (a, b)
log = log 1 = 0. (5.13)
pA (a)pB (b)
The next example illustrates mutual information for image pairs.
Example 5.3
Consider the two images A and B with only four ‘pixels’ each (fig-
ure 5.12). They have only two gray levels: 0, or black, and 1, or white.
A B
0 1 0 1
1 0 1 0
Fig. 5.12: Abstract example images. Four pixels, two gray values (0 and
1).
and of course
Thus, we can calculate the mutual information between the two im-
ages as
I(A, B) = 0.5 log 0.5/(0.5)2 + 0.5 log 0.5/(0.5)2 = 1. (5.18)
The next example shows that the two images do not have to use the
same colors.
Example 5.4
A B
0 1 1 0
1 0 0 1
Fig. 5.13: Gray level 0 is mapped to gray level 1 and vice versa. Here the
mutual information I(A, B) also equals one!
Now, by contrast, assume one outcome has probability 1 and the other
0, i.e.
pA (0) = 1 pA (1) = 0.
Example 5.5
Figure 5.14 shows a case in which the mutual information of two im-
ages is zero.
A B
0 1 1 1
1 0 1 1
Example 5.6
One final example is shown in figure 5.15. Three pixels are the same,
but not the fourth.
A B
0 1 0 1
1 0 1 1
Assume there are 256 gray levels in both images A and B, (0, . . . , 255).
The first step is to estimate the distribution pA . Set up a vector with
256 entries, one entry for each possible gray level. We count the num-
ber of voxels (in 2D: pixels) for each gray level value. This means:
we step through the image voxel by voxel. Whenever we meet a voxel
with gray level j, we increment the entry with index j in the vector by
1.
From this vector, we compute an approximation for pA : we normalize
the vector by dividing by the total number of voxels, so that the sum
of all entries becomes 1. We call this vector the gray level histogram
for pA . In the same way, we compute the histogram for pB .
After calculating the two histograms, the only thing missing is the
joint distribution of the two random variables. Therefore, consider a
table with 256 × 256 entries as a joint histogram with all entries ini-
tialized to 0. First, bring the two images into an approximate match-
ing position (i.e. by superimposing the images). For each voxel pair
5.4 Intensity-based Registration 189
with gray levels (a, b) increment the table entry with index (a, b) by
1. Compute the estimate of the joint distribution by normalizing the
table entries (divide by number of voxel pairs), so that the sum of all
table entries becomes 1. The computation is illustrated in figure 5.16.
a
CT
MR 0 255
0
b #(a,b)
255
Fig. 5.16: Table for computing the joint distribution pA,B . Place the num-
ber of voxels having gray level a in A and gray level b in B at the table
position with index a, b.
Let image A be the fixed image, image B be the floating image, i.e. B
moves over image A (rotation and translation), until the final registra-
tion position is reached.
We move the floating image under six transformation parameters:
three angles α, β and γ and three translations tx ,ty and tz .
Criterion: For each position of B, evaluate the mutual information via
equation 5.11. The images match optimally for the position where the
mutual information is maximal.
190 5 Navigation and Registration
One problem that arises is that the pixel (voxel) center in the floating
image may not always be exactly above the pixel (voxel) center of the
fixed image (Figure 5.17). Thus, in figure 5.17, we have assumed that
in one box-section of our CT image (fixed image) we have four voxels
with gray levels of 97, 98, 102 and 103, respectively. And we assume
a voxel of the MR, inside this box has gray level 13.
MR
97 98
CT 13
102 103
For the example of figure 5.17, the CT-voxel with gray level 98 is the
nearest neighbour of the MR-voxel with grey level 13. We thus incre-
ment the table entry for the index pair (98, 13) in the joint distribution.
5.4 Intensity-based Registration 191
Trilinear Interpolation
The search for the registration position starts with the initialization,
i.e. images are overlaid so that centers and axes match. The optimiz-
ation is then decomposed into a separate search along several axes.
Thus, we compute the mutual information I(A, B) for the current po-
sition (start position) according to equation 5.11. We then move the
floating image by a small amount along one of the six axes (three
translational and three revolute axes). For this new position, we re-
compute the mutual information I(A, B). This process is repeated to
maximize the mutual information. Notice that the mutual information
depends on the position of the floating image relative to the fixed im-
age. As a consequence, we must recompute the joint distribution for
each new position in this process. However, the distributions pA and
pB remain unchanged when moving the floating image. The maxim-
ization process can, for example, be guided by Levenberg-Marquardt
optimization [2].
192 5 Navigation and Registration
Binning
Suppose the grid points in G and G0 are landmarks in two images. Our
goal is to deform the image underlying G0 , while moving each land-
mark point in G0 to its target position in G. We find an initial regis-
tration position with landmark-based registration, as described above.
After this step, the positions of the landmarks in G and G0 may not
match exactly, and we deform G0 appropriately, to bring the points
Interpolated value
du dw dv
u w v
Original values
du dw
dv
u v w
Given an image point with coordinates (x, y)T , we then define a poly-
nomial
a0 + a1 x + a2 y + a3 xy + a4 x2 + a5 y2 + a6 x2 y + a7 xy2 + a8 x3 + a9 y3 .
(5.20)
5.5 Image Deformation 195
This polynomial (in the two variables x and y) describes the shift of
this image point in x-direction. Likewise,
b0 + b1 x + b2 y + b3 xy + b4 x2 + b5 y2 + b6 x2 y + b7 xy2 + b8 x3 + b9 y3
(5.21)
describes the shift in y-direction.
The parameters ai and bi in this polynomial must be adjusted in such
a way that the distances between the model points (xm , ym )T (i.e. grid
points in G) and the corresponding landmarks in G0 are minimized.
The coordinates of the landmarks in G0 are given by (xd , yd )T . Shifting
under the correction polynomials will give points
pi
zi
Fig. 5.21: Deforming a thin metal plate. The plate is deformed by a force
pushing downwards at point pi . The magnitude of the displacement at
pi is zi .
X-ray detector
Calibration
disk
Example 5.7
Example 5.8
Table 5.1 defines two-dimensional grid points and shift values. Fig-
ure 5.23 shows bilinear interpolation for the grid in table 5.1 and fig-
ure 5.24 shows interpolation with thin-plate splines on the same grid.
(0, 0)T 1
(1, 1)T 1
(1, 0)T 0.5
(0, 1)T 0
(2, 1)T 0
(0, 2)T 1
Table 5.1: grid points and shift values for figure 5.23 and figure 5.24.
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0
0.6 0.5
0.8 1
1.5
1 2
y x
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0
0.6 0.5
0.8 1
1.5
1 2
y x
∑x ∑y T (x, y)A(x + u, y + v)
C(u, v) = q q (5.28)
∑x ∑y [T (x, y)] ∑x ∑y [A(x + u, y + v)]2
2
5.6 Hand-eye Calibration 201
Here T is the template, and T (x, y) is the gray level at pixel (x, y) in
the template.
Maximizing the cross correlation by moving the floating template
gives the desired estimate for the registration position of the template.
In this way, we can find landmarks and correspondence pairs. But how
do we find good templates? Obviously, one would like to select inter-
esting windows rather than uninteresting ones. An uninteresting win-
dow has simply no contrast at all, i.e. all pixels have the same color
or gray level. Thus we choose windows having two very distinct gray
levels, assuming such windows will show delicate bony structures or
salient landmarks. To find such windows, we compute histograms for
each window in the floating image. A histogram simply plots the num-
ber of pixels having a given gray value over the x-axis with the gray
levels. The histogram is called bimodal, if we can distinguish two dis-
tinct peaks.
This matching process described here can also be extended to mutual
information matching.
Elastic registration in real-time is difficult and error-prone. In later
chapters we will describe tracking methods based on motion correla-
tion which allow for motion tracking in real-time.
The position and orientation data will then be returned to the host
computer.
However, you may argue that we now have the position of the pointer
with respect to the internal system of the camera. How does the robot
know where the internal system of the camera is? This type of prob-
lem is a calibration problem, and this is what we will address in this
section.
To illustrate hand-eye calibration in a more specific case, assume an
LED marker is attached to a robot gripper. The tracking system reads
the 3D position of the marker. We have attached the marker on the
robot in a rather arbitrary position, and we have placed the tracking
system on a floor stand in front of the robot. The situation is illustrated
in figure 5.25.
CT
M
ET
M
RT RT
E C
Two of the four matrices are known. Firstly, the position of the tool
flange/gripper coordinate system with respect to the robot base is
known. The robot is a kinematic chain, and matrices describe the po-
sition of the gripper or tool flange with respect to the robot’s base
coordinate system. Thus R TE is known.
Secondly, the tracking system computes the position of the marker,
so that the matrix C TM is known. (Most tracking systems output the
marker orientation as a quaternion, but we can convert this to a mat-
rix.)
However, the two other matrices are unknown. We typically do not
know where exactly the internal coordinate system of the tracking
system is located with respect to the robot’s base, since we just put
the tracking system on a floor stand, free to move in the room. Also
we attach a marker to the robot, but we are not sure how to determine
the transformation matrix from robot tool system to the marker. Now,
hand-eye-calibration is the process of computing the two unknown
matrices in this scenario.
We note that the transformation from robot base to marker via the
robot gripper must give the same result as the transformation from
robot base to marker via the tracking system. This observation can be
expressed in an equation:
R E
TE TM =R TC C
TM (5.29)
Now recall that two of the four matrices are known, namely R TE and
C T . Thus only two matrices will have to be computed from this
M
equation.
We rename the above matrices to obtain a simpler notation.
204 5 Navigation and Registration
AX = YB (5.30)
So the last two equations are the same, only we have changed the
names of the matrices. Now A and B are known, and X and Y are
unknown.
But notice that we can move the robot. After moving to a second pos-
ition, we obtain a second equation.
So assuming the first position gives
A1 · X = Y · B1 . (5.31)
The second position will give the equation
A2 · X = Y · B2 . (5.32)
We rearrange the two equations to have
A1 · X · B−1 −1
1 = A2 · X · B2 , (5.33)
or
A’ X = X B’, (5.34)
where
A’ = A−1
2 A1 , (5.35)
B’ = B−1
2 B1 . (5.36)
Ai X = Y Bi , i = 1, . . . , n. (5.37)
5.6 Hand-eye Calibration 205
Setting Ci = B−1
i , we have
Ai XCi = Y. (5.38)
The matrix X is a homogeneous 4x4-matrix, and hence there are
12 non-trivial entries in this matrix. Likewise, Y has 12 non-trivial
entries. Thus equation 5.38 gives rise to a linear system of equations
with 12 · n equations and 24 variables.
Since we assume n > 2, the resulting system has more equations
than variables. We have already encountered a method for solving
such over-determined systems in the context of image reconstruction.
The ART algorithm is a suitable method for this problem (see exer-
cise 1.2). After solving the linear system we obtain the values for the
variables x11 , . . . , x34 and y11 , . . . , y34 . This gives matrices describing
the transformations needed for the calibration.
RT
E
ET
Ultrasound camera M
Marker
RT
C
CT
M
Water tank
Exercises
Exercise 5.1
Verify the values I(A, B) for the examples in figures 5.13 to 5.15.
Exercise 5.2
Show that I(A, B) = H(A) if the two images A, B are the same, i.e.
A = B. Discuss this in the light of example 5.3 above. Find the value
for H(A) in case the image A has an even distribution of gray levels,
i.e. there are 256 gray levels, and all gray levels occur equally often in
the image.
Exercise 5.3
Show that
Exercise 5.4
Exercise 5.5
Exercise 5.6
Exercise 5.7
The χ 2 -function tests the hypothesis that two random variables are
independent. It can also measure the degree of dependency between A
and B. (Pronounce χ 2 as ’ki-squared’, as in ’kite’). Thus, we measure
the ‘distance’ between the observed joint distribution pA,B (a, b) and
the hypothetic distribution which would be observed if A and B were
independent:
Exercise 5.8
Set
p1 = 0 f (p1 ) = 0
p2 = 1 f (p2 ) = 1
p3 = 2 f (p3 ) = 0, (5.43)
208 5 Navigation and Registration
X-ray detector
z0
x0
S0
z1
x1
S1
X-ray source
Fig. 5.27: Projection geometry of a C-arm. The focal length f is the dis-
tance between the source and the detector.
Exercise 5.9
Exercise 5.10
b1 x<2
b2 y n
b3 y<3 y<2
a b4 y n y n
b3 b1 b4 b2
a) Find an example where this simplified version does not report the
correct nearest neighbor.
b) Describe and implement a way to the repair the problem from part
a).
Exercise 5.11
Summary
are not known. ICP alternates between a move step (as in landmark-
based registration) and a match step. The match step computes the
correspondences by a search for the nearest neighbor of each point.
Contour-based registration is a direct extension of ICP. We replace
the point-to-point distance in ICP by a point-to-line distance. MI re-
gistration allows for multi-modal image registration, by automatically
finding correspondences between pixel intensities (or voxel intensit-
ies). We start by placing image B on top of image A. Then move B, and
compute the value I(A, B) for each position. The position giving the
maximum is the registration position. The movement of B is guided by
Levenberg-Marquardt optimization. To compute I(A, B), we need the
histograms for pA , pB , and pA,B . Each histogram is a table. The table
entries correspond to the image intensities. For each table entry, we
count the number of pixels having the corresponding intensity. Nor-
malizing gives an approximation for pA , pB , and pA,B . Intensity-based
registration is rarely used in navigation, due to limitations in speed
and robustness. One typical application of MI-registration is CT-MRI
fusion. In radiation therapy, tumors must be delineated in the CT. If
the tumor boundary is not well visible, then the MRI can help improve
the delineation.
To apply medical imaging in navigation, distortions must be correc-
ted, and we must determine parameters describing the projection geo-
metry. For x-ray camera calibration, calibration objects with grid pat-
terns are placed into the beam path. The geometry of the grid pattern is
known from measurements. Based on this, we deform the image such
that the markers match the geometry of the grid. Thin-plate splines
and cubic splines reduce the deformation problem to a linear equation
system.
Hand-eye calibration aligns internal coordinate systems of cameras
and robots. The problem of hand-eye calibration leads to a matrix
equation of the form AX = YB, where A and B are known, and X and
Y are unknown. We can set up a number of such equations, by moving
the robot through a series of n different positions, giving n equations
Ai X = YBi , where i = 1, ..., n. The resulting system is overdetermined
for larger n. The ART algorithm provides a fast and practical solution
method.
212 5 Navigation and Registration
Notes
References
219
220 6 Treatment Planning
isthmus of the
femur neck
wedge angle
cutting plane
shaft
lateral corticalis
it hard to drill the channel. But now it is difficult to find the correct
angle for the drill, since the implant shaft must touch the lateral femur
shaft after removing the wedge.
Cut
Fig. 6.2: The implant must be placed before the cut is made.
6.1 Planning for Orthopedic Surgery 221
To plan the intervention, a desired new position of the femur head with
respect to the bone must be determined. Hence, three parameters are
specified in the planning phase: channel direction, cutting plane and
wedge plane.
The planning procedure consists of the following steps:
1. Take two x-ray images (anterior/posterior and axial)
2. Construct a simplified model of the femur bone before wedge re-
moval
3. Determine the desired pose, based on this model
4. Plan the implant position
Images for step one are taken during the intervention. Planning is
based on the x-ray images, hence no registration step is needed.
For step two, we extract the following five features from the x-ray im-
ages: femur shaft axis, center of the femoral neck isthmus (see also
figure 6.1), femoral head center point, femoral shaft radius, and os-
teotomy plane. In the navigation set-up with an optical tracking sys-
tem and the C-arm, we can obtain the spatial position of these features
in the following way. In the simplest case, the feature is a point. The
point is marked in both images on the computer screen. Cast a line in
space from the x-ray source point to the marked image point on the
detector. This is done for both images. The intersection point of the
two lines gives the spatial position of the point. In a similar way, the
spatial positions of the remaining features can be extracted.
The features thus extracted define a simplified bone model, and step
three is done as follows: We define a coordinate system for the femur
shaft. This coordinate system does not move with the shaft. All mo-
tions of the femur head and neck caused by repositioning are now de-
scribed with respect to this coordinate system. The coordinate system
is defined in the following way:
• The origin is the intersection point between the cutting plane and
the shaft axis (defined in step 2).
• The z-axis is the extension of the shaft axis and points into the
proximal direction (i.e. towards the patient’s head).
• The y-axis is a vector in the AP-direction in the anatomy. Note that
the AP-direction is defined implicitly by the user, when taking the
x-ray image in AP-direction under navigation. To obtain a valid
222 6 Treatment Planning
z
Femur
Head
Wedge
x
Fig. 6.3: Forward and inverse planning for intertrochanteric femur os-
teotomy. Frontal view onto the right femur. In this view, a valgus of 15
degrees is specified. Original posture of the femur (solid) and new pos-
ture (transparent).
With the features defined above, a model of the femur can be visual-
ized on the computer screen (figure 6.3).
Three rotational parameters and three translational parameters are set
by the surgeon. In an inverse planning step, we compute the result-
6.2 Planning for Radiosurgery 223
timal, since the healthy tissue along the beam path will receive the
same dose as the target. On the right side of the figure, two beams are
used. Generally, the right plan is preferrable over the left plan because
we can ensure that the dose in most of the healthy tissue remains be-
low a given threshold.
If we increase the number of beams in figure 6.4, we obtain a so-
called isocentric treatment (figure 6.5). All beams intersect, so that
the treated volume is a sphere.
Fig. 6.5: Isocentric treatment. Beam axes cross in a single point, the
treated volume is spherical.
6.2 Planning for Radiosurgery 225
T
.
.
.
.
.
In practice, the target region often is not spherical. Figure 6.6 shows
an elongated tumor region T and a critical structure B.
During treatment, the beam is placed in a series of fixed positions
(called treatment nodes). At each position the beam is activated for
a certain amount of time. Thus, treatment planning consists of two
phases:
For phase 1, we must choose the placement of the beams. The number
of beams can be large (i.e. 1000 beams). For phase 2, the duration of
activation of each beam must be computed. This is an optimization
problem. The computations for phase 1 and phase 2 interact. Thus it
is useful to iterate over both phases.
After we have computed a treatment plan, we must simulate the effect
of this plan in terms of dose delivered to the individual regions in
the target area. Such a simulation must then be visualized. Thus we
distinguish between forward and inverse planning:
226 6 Treatment Planning
Inverse planning. Given: only the regions (here T for target and C for
critical). Find appropriate beam directions and weights.
Fig. 6.7: Sphere packing for non-spherical target shape results in uneven
coverage (hot spots).
Consider the target volume in the figure 6.8. This shape is obtained
by sweeping a sphere. Now, instead of having a single isocenter, we
sweep the isocenter along a line segment. This gives rise to a series of
discretized isocenter points (figure 6.9b). Each of these new isocen-
ter points is then treated as if it was the center of a spherical target
volume. The result is an overlay of a series of spherical volumes. The
main difference between this scheme and sphere packing is that now
the spheres do overlap, while homogeneity is achieved by the sweep-
ing process.
Now the question remains of how to place beams for arbitrarily
shaped target volumes. A straight-forward extension of the above line-
sweeping placement scheme is illustrated in figure 6.10. Thus, the
outer boundary of the target volume is retracted, and a series of iso-
center points are placed on the retracted surface. Sweeping is then
applied to this retracted surface.
228 6 Treatment Planning
a)
b)
C
T
We extend each line segment to a line. The lines partition the plane
into so-called cells. Here, we define a cell as a maximally connec-
ted region containing no points on the given lines. (Thus, we consider
geometric cells, not anatomic cells.) All cells thus defined are convex.
Now we take a reference point in each cell. This is an arbitrary point
inside the cell. For each such point we can compute a label. This label
simply lists all regions which contain the cell. The label also includes
the beams containing the cell. Our definition of geometric cells can
be extended in two ways: first, the target is three-dimensional, and
second, our beams are typically cylinders. For both extensions the ne-
cessary modifications are straightforward.
The partitioning of space which has been induced by the beams and
the delineation gives rise to an equation system: We take the label of a
single cell c. Assume c is contained in the target, and also in three of
the beams, called s1 , s2 and s3 . Then we wish to bound the dose in the
target from below. Thus the dose should be above a given threshold
value t.
Our goal is to compute the duration of activation of each beam. In our
example (figure 6.11), we assign variables x1 , x2 , x3 for these unknown
durations of activation for the beams s1 , s2 and s3 .
Overall we have an equation of the form
x1 + x2 + x3 ≥ t. (6.1)
230 6 Treatment Planning
Here, the value t specifies the lower threshold for the dose in the target.
In the same way we can set a bound for the dose in critical regions,
i.e. we can state that the dose should not exceed a value t 0 .
Let n be the total number of beam directions. Assume the beams are
numbered in such a way that c is contained in the beams s1 , ..., sk . By
x1 , ..., xk we denote the weights of the beams s1 , ..., sk . Let c be a cell
in a critical region. The dose in c should not exceed t 0 . We obtain the
equation
x1 + ... + xk ≤ t 0 . (6.2)
After having computed all cells, this results in a system of inequalities
of the form (for n beams)
x1 + x2 ≤ 10 (6.5)
x1 ≥ 5
x2 ≥ 4
Obviously, these constraints are feasible, i.e. there are values for x1
and x2 satisfying these inequalities. The constraints are visualized in
figure 6.12.
x2
10
5 10 x1
Fig. 6.12: Feasible region for the linear constraints in equation 6.5
x1 + x2 ≤ 10
x1 ≥ 5
x2 ≥ 6
(6.6)
232 6 Treatment Planning
x2
10
5 10 x1
a11 x1 + · · · + a1n xn ≤ b1
a21 x1 + · · · + a2n xn ≤ b2
..
.
am1 x1 + · · · + amn xn ≤ bm
where m and n can be different. Notice that the ≤-sign used here is
not a restriction. Inequalities in linear programming can also take the
form
ai1 x1 + · · · + ain xn ≥ bi .
We have used only the form with a ≤-sign, since we could always
transform a greater-than inequality into a less-than inequality by mul-
tiplying the entire inequality by -1.
The inequality system determines a polytope in n-dimensional space
(n is the number of variables, i.e. x1 , .., xn ). Notice that n is not equal
to two or three, although the system describes the constraints for the
three-dimensional case. Rather, the dimension n is the number of vari-
ables in the system. A feasibility test (F) only checks whether this
polytope is not empty, i.e. if there is a solution or not. Two problems
are distinguished in linear programming:
and
Notice that the above system of inequalities in equation 6.6 is not the
standard form of linear programming in the literature. The standard
linear programming (LP) is defined in the following way: Given con-
straints in equation 6.6, maximize the expression
c1 x1 + ... + cn xn , (6.7)
where all xi ≥ 0. The expression in equation 6.7 is called the objective
function . Figure 6.14 shows the situation. Here the objective function
is 1x1 + 0x2 .
allowed
maximum
The gray region is the feasible region. Notice that it does not extend
into the negative quadrants, since we require all xi to be positive. The
values ci determine a vector c. We maximize in direction of this vector
c, while remaining within the feasible region.
Thus (LP) contains a maximization, which is more than (F) alone.
However, in standard form, (LP) imposes the restriction that all xi be
positive. Convince yourself that this restriction can be removed: sub-
stitute each occurrence of the variable xi in the given LP-problem by
the expression (yi − y0i ), where yi and y0i are new variables. After sub-
stitution, the new problem has the standard LP-form with all variables
yi , y0i ≥ 0.
234 6 Treatment Planning
The matrix form of the LP problem in equations 6.6 and 6.7 is ob-
tained by setting up a coefficient matrix A from the coefficients ai j ,
and two vectors b = (b1 , . . . , bm )T , c = (c1 , . . . , cn )T . An LP problem
can then be stated in the compact form
feasible
polytope
Fig. 6.15: Solving a system of linear equations does not necessarily yield
a feasible point
After extending the input linear system in this way, the system will
always be feasible. In the optimization step, we then minimize the
sum of the slack variables, instead of the sum of the variables xi .
As noted above, many of the grid points will have the same labels.
Obviously, it is useful to have a small number of inequalities, since
the computing time for solving the resulting system is often in the
range of several hours. To this end, we remove grid points having the
same label as a grid point found earlier. This is straightforward, but
we can do better.
Some grid points not having the same labels may still be redundant,
since the restrictions they impose are already subsumed in other con-
straints. To discuss the process of removing redundant labels, we re-
turn to the definition of geometric cells in the plane. Although the
geometric structures we use are not all lines (e.g. they are cylinders in
space and line segments in CT slices) the concepts can readily be visu-
alized in the plane. Consider oriented lines as in figure 6.16. One cell
(shown shaded) is marked. It is the intersection of several half-planes.
The half-planes are bounded by lines.
Fig. 6.16: Maximal cell. Lines are oriented, as indicated by arrows. The
number of lines having the shaded cell in their interior half-space attains
a local maximum.
6.2 Planning for Radiosurgery 237
Fig. 6.17: CT slice for sample case with polygon delineating the target
238 6 Treatment Planning
Fig. 6.18: Inverse planning input and output for sample case in fig-
ure 6.17. (b) Stack of tumor polygons in each positive slice. (c) 80%
isodose curve after computing the beam weights with linear program-
ming. (d) overlay of (b) and (c).
6.2.2.4 Collimators
100
80
60
Volume [%]
40
20
2769
Dose [cGy]
Fig. 6.21: Beam collimator. The beam is not cylindrical, but a cone. The
radiation source is at the tip of the cone.
Target
Movable
leaves
Fig. 6.24: Motion of organs in soft tissue. The target moves faster than a
critical structure in the vicinity of the target.
Exercises
Maximize
5x1 + x2
subject to
x1 + x2 ≤ 10
x1 ≥ 5
x2 ≥ 4.
Summary
tions for irradiation, the activation duration (called weight) for each
beam must be found. Linear programming provides an effective and
fast way for computing such weights.
In linear programming, a set of linear constraints is checked with re-
spect to feasibility. After computing a feasible point (i.e. a point sat-
isfying the constraints), we can move this point in such a way that an
objective function is maximized or minimized. This objective func-
tion is also a linear function. Assuming the dose from different beams
add up to a total dose, we obtain a linear system of inequalities. In this
system, the dose from each beam is represented as one variable, and
points covering the target region give rise to linear constraints in these
variables. We saw that linear programming directly extends methods
for solving linear equation systems, which we have frequently applied
in the previous chapters.
To reflect organ motion, we modify the linear program accordingly.
Based on four-dimensional imaging, the organ motion can be quanti-
fied. We can thus protect critical structures which would receive too
much dose if organ motion would bring this structure closer to the
target during respiration or pulsation. In the exercises, at the end of
this chapter, we look at duals of linear programs. Variables in the dual
linear program correspond to inequalities in the primal and vice versa.
In the next chapter, we will extend the mathematical tools used so far
to address the problem of tracking soft tissue motion.
Notes
References
ning for robotic radiosurgery.” Medical Physics, vol. 32, no. 12,
pp. 3786–3792, Dec 2005.
[5] G. Fichtinger, J. P. Fiene, C. W. Kennedy, G. Kronreif, I. Ior-
dachita, D. Y. Song, E. C. Burdette, and P. Kazanzides, “Robotic
assistance for ultrasound-guided prostate brachytherapy,” Med-
ical Image Analysis, vol. 12, no. 5, pp. 535 – 545, 2008, special
issue on the 10th international conference on medical imaging
and computer assisted intervention - {MICCAI} 2007.
[6] H. Gottschling, M. Roth, A. Schweikard, and R. Burgkart, “In-
traoperative, fluoroscopy-based planning for complex osteotom-
ies of the proximal femur,” International Journal of Medical Ro-
botics and Computer Assisted Surgery, vol. 1, no. 3, pp. 67–73,
2005.
[7] S. Nakajima, H. Atsumi, A. H. Bhalerao, F. A. Jolesz, R. Kikinis,
T. Yoshimine, T. M. Moriarty, and P. E. Stieg, “Computer-
assisted surgical planning for cerebrovascular neurosurgery,”
Neurosurgery, vol. 2 41, pp. 403–410, 1997.
[8] N. Wilson, K. Wang, R. W. Dutton, and C. Taylor, “A soft-
ware framework for creating patient specific geometric models
from medical imaging data for simulation based medical plan-
ning of vascular surgery,” in Medical Image Computing and
Computer-Assisted Intervention - MICCAI 2001, ser. Lecture
Notes in Computer Science, W. J. Niessen and M. A. Viergever,
Eds. Springer Berlin Heidelberg, 2001, vol. 2208, pp. 449–456.
[9] C. M. C. Tempany, E. A. Stewart, N. McDannold, B. J. Quade,
F. A. Jolesz, and K. Hynynen, “MR imaging–guided focused ul-
trasound surgery of uterine leiomyomas: A feasibility study,” Ra-
diology, vol. 226, no. 3, pp. 897–905, 2003, pMID: 12616023.
[10] A. Jordan, R. Scholz, P. Wust, H. Fähling, and R. Felix, “Mag-
netic fluid hyperthermia (MFH): Cancer treatment with AC mag-
netic field induced excitation of biocompatible superparamag-
netic nanoparticles,” Journal of Magnetism and Magnetic Ma-
terials, vol. 201, no. 1-3, pp. 413 – 419, 1999.
[11] S. Rossi, M. Di Stasi, E. Buscarini, P. Quaretti, F. Garbagnati,
L. Squassante, C. T. Paties, D. E. Silverman, and L. Buscarini,
“Percutaneous rf interstitial thermal ablation in the treatment of
References 249
251
252 7 Motion Correlation and Tracking
Figure 7.2 shows an example. In the figure, the motion of the skin
surface is illustrated by the two vertical arrows. Thus, the skin surface
(arc line in the figure) moves up and down during the respiratory cycle.
The target (gray circle) moves in left-right direction. At full exhalation
(left image), the target is on the left side. At full inhalation, the target
is on the right side.
Our main observation is the following: tracking the skin surface is
easier than tracking a tumor inside the body.
7.1 Motion Correlation 253
infrared data
(continuous)
Fig. 7.3: Correlation-based tracking: the bold curve stems from the real-
time sensor (here optical tracking). A series of points on the dashed
curve are known (here: x-ray images), but these points are in the past,
since we need segmentation time. The bold curve and the dashed curve
(here: ground truth and target position) are in correlation. Our goal is
to estimate the unknown dashed curve from this data.
By taking more than two positions of the target, we can refine this
correlation model. Thus curvilinear motion of the target can be ap-
proximated by a series of line segments, giving a model for this mo-
tion. The correlation model is continuously updated during treatment
to reflect changes in breathing pattern.
As in the case of registration, we can replace specific image mod-
alities but maintain the general method. In our description, optical
254 7 Motion Correlation and Tracking
tracking takes the role of a real-time sensor, and x-ray imaging com-
putes ground truth intermittently. Motion correlation then determines
the ground truth in real-time.
Example 7.1
Fig. 7.4: Data points from optical tracking and x-ray imaging. Both
types of sensor data points have time stamps. This allows for establishing
the correspondences between the data points.
(x-coordinate)
X-ray
optical
(x-coordinate)
(1, 9)
(2, 21)
(3, 33)
(4, 39)
(5, 50)
If we take a closer look, we see that the data points are roughly on a
line. We wish to find the line that best fits the data. This line does not
necessarily contain all the points, but it should match the data. Our
goal is thus to fit the line y = mx + b to the data, i.e. solve the system
b + 1m = 9
b + 2m = 21
b + 3m = 33
(7.1)
b + 4m = 39
b + 5m = 50.
Bw = y (7.4)
Since B is a 5x2-matrix, BT will be a 2x5-matrix, and we can rewrite
our system to:
BT Bw = BT y (7.5)
This equation is called the normal equation for the system in equa-
tion 7.4. BT B is a 2x2-matrix. Now this system can be solved by
Gaussian elimination. Notice that the transition from Bw = y to
BT Bw = BT y substantially reduces the matrix size. It can be shown
that solutions of the least-squares problem (equation 7.2) are solu-
tions of equation 7.5 and vice versa (see exercise 7.1 at the end of this
chapter).
258 7 Motion Correlation and Tracking
inhale
exhale
Fig. 7.6: Hystersis: inhale path of target differs from exhale path.
The difference can be quite large. A simple way to address this prob-
lem is to fit two lines to the data instead of just one. Here, we would
use the optical tracking data to determine a gradient in the data, i.e.
decide whether a given data point is an inhale or an exhale point. We
would then map to one of the two regression lines. It turns out that res-
piratory traces are indeed curves and not just lines, so obviously we
should attempt to fit curves to the data points. Polynomials of degree
two or higher can be fit to the points by least squares methods. How-
ever, it is useful to include more sensor information. Beyond gradient
information on the motion, speed and acceleration can be measured. A
breath temperature sensor distinguishes between respiratory phases. A
spirometer measures respiratory air flow. Similarly, data for pulsation
can be acquired in a variety of ways.
Support vector methods collect heterogeneous data to set up training
samples. Support vector regression computes a functional representa-
tion of correlations in the data. In the most simple form, the process
resembles line regression. However, very complex correlations can be
captured.
7.3 Support Vectors 259
output
input
w
x1 x4
x2 x3
w1 x1 + w2 x2 + d = 0. (7.6)
To visualize this representation of lines, notice that w = (w1 , w2 )T is
the normal vector of the line. But furthermore, the direction of this
normal vector w indicates an orientation of the line. Thus a point p =
7.3 Support Vectors 261
(p1 , p2 )T is above the line (in the direction into which w is pointing),
if
w1 p1 + w2 p2 + d > 0. (7.7)
And p will be below the line if w1 p1 + w2 p2 + d < 0. The representa-
tion of lines discussed here also has a close relationship to distances.
Namely, if w is a unit vector, then we know the distance of p from the
line. It is simply the magnitude of the value w1 p1 + w2 p2 + d.
Notice also that this representation of lines is not unique. If wx+d = 0
is a line, then (λ wx) + λ d = 0 is the same line. Thus, we can scale
our line equations. The scaling that turns the normal vector into a
unit vector (i.e. λ = 1/kwk) is often useful. One advantage of this
representation of lines is that it works in any dimension. Hence,
w1 x1 + w2 x2 + w3 x3 + d = 0 (7.8)
is a plane in space, with much the same distance properties we just
discussed. And the same holds for higher dimensions.
x11 w1 + x12 w2 + d ≥ 0
x21 w1 + x22 w2 + d ≥ 0
x31 w1 + x32 w2 + d ≤ 0
x41 w1 + x42 w2 + d ≤ 0.
(7.9)
y1 (wx1 + d) ≥ 0
y2 (wx2 + d) ≥ 0
y3 (wx3 + d) ≥ 0
y4 (wx4 + d) ≥ 0
(7.10)
The system in equation 7.10 is the same as the system in equation 7.9.
We have turned all inequalities into ≥-inequalities.
Now the linear program above will return a line such as the line shown
in the next figure (figure 7.9). Notice that two of the points in figure 7.9
are actually on the separating line, hence not strictly below or above
it. These are the points x2 and x4 . Indeed linear programming will
typically return this type of separating line. One can argue that this
type of separating line does not represent the data points very well.
x1
x4
x2
x3
Fig. 7.9: Separator line for the classification problem in equation 7.10
Assume the points x1 , x2 , x3 , and x4 are samples from a larger data set
with more (unknown) sample points. We would again like to classify
7.3 Support Vectors 263
the data into two subsets. It would then seem more reasonable to use
the dashed line in figure 7.10 as a separator line.
x1
x4
x2
x3
Fig. 7.10: Maximum margin separator line for the classification problem
in equation 7.10
The advantage of this dashed line is that the margin between the sep-
arating line and the data points is maximal. In the picture, the data
points can be thought of as repelling points which push the separating
line away by the maximum amount. Likewise, we may say the data
points ‘support’ the dashed separator line.
The basis of support vector theory is the observation that one can find
such a maximum margin separator plane with quadratic programming.
Specifically, we replace the above linear program for finding the sep-
arator plane by the program:
Minimize
w21 + w22
subject to
y1 (wx1 + d) ≥ 0
y2 (wx2 + d) ≥ 0
y3 (wx3 + d) ≥ 0
y4 (wx4 + d) ≥ 0.
(7.11)
264 7 Motion Correlation and Tracking
7.3.3 Regression
yi − (wxi + d) ≤ ε
(wxi + d) − yi ≤ ε.
The constraints mean that the distance between the value (wxi + d)
and the value yi should not be bigger than ε. Recall that the values yi
in the case of regression are the function values at the points xi . Thus
we no longer have yi = ±1. Instead the values yi are real numbers.
Notice also that there are two constraints per sample point xi . The
two constraints bound the distance from below and above. This means
that the regression line must not have a distance more than ε from all
samples (xi , yi ). This may not always be feasible, i.e. if ε is too small.
To allow for finding a regression line, even if the points are far apart,
we introduce slack variables. Slack variables allow for a violation of
7.3 Support Vectors 265
yi − (wxi + d) ≤ ε, (7.13)
we ask for a little less, namely
yi − (wxi + d) ≤ ε + ξ , (7.14)
where ξ is a variable (ε is a constant)!
ξ is treated just like any other variable in quadratic programming (or
linear programming), i.e. ξ ≥ 0. You should take a minute to see why
equation 7.14 is indeed a weaker constraint than equation 7.14.
We would like ξ to be as small as possible, since ξ is the amount by
which the original constraint is violated.
Thus we include the minimization of the ξ -values into the above quad-
ratic program:
Minimize
ww + ξ1 + ξ10 + ξ2 + ξ20 . . .
subject to
y1 − (wx1 + d) ≤ ε + ξ1
(wx1 + d) − y1 ≤ ε + ξ10
y2 − (wx2 + d) ≤ ε + ξ2
(wx2 + d) − y2 ≤ ε + ξ20
...
ξi , ξi ≥ 0.
0
(7.15)
7.3.3.1 Dualization
1 m m m
− ∑ (αi − αi )(α j − α j )xi x j − ε ∑ (αi + αi ) + ∑ yi (αi − αi0 )
0 0 0
2 i, j=1 i=1 i=1
subject to
m
∑ (αi − αi0) = 0
i=1
and
Remark 7.1
We mentioned the fact that the dual has several most surprising ad-
vantages over the primal. We will now set out to include kernel func-
tions into the dual system in equation 7.16. This will allow us to find
non-linear regression functions, which we need in our application.
But before we do that, one more thing remains:
Having solved the quadratic program in equation 7.16, e.g. with the
Matlab routine quadprog, we get values for the variables αi , αi0 . How
can we get back to the values for w and d from the values for αi , αi0 ?
The answer is: As part of the dualization process, we obtain equa-
tions (setting to zero derivatives with respect to dual variables). These
equations state that
m
w = ∑ (αi − αi0 )xi = 0. (7.17)
i=1
The value d can be computed in the following way:
m
d = y j − ∑ (αi − αi0 )xi x j − ε, (7.18)
i=1
for a j with α j 6= 0, or
268 7 Motion Correlation and Tracking
m
d = y j − ∑ (αi − αi0 )xi x j + ε, (7.19)
i=1
for a j with α 0j 6= 0.
Finally, we must evaluate the function value y for an unknown sample
x.
Here, of course, we have
f (x) = wx + d, (7.20)
but it is preferable to write
m
f (x) = ∑ (αi − αi0 )xi x + d, (7.21)
i=1
where we have simply inserted our formula for w derived above. The
reason for this rewriting will become clear after we have discussed
kernel functions. We will do this next.
7.3.4 Kernels
x1 x1 y1 y1
x1 x2 y1 y2 2 2 2 2
F(x)F(y) =
x2 x1 y2 y1 = x1 y1 + 2x1 x2 y1 y2 + x2 y2 (7.24)
x2 x2 y2 y2
1.5
0.5
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
0
4
2 4
0 3
2
−2 1
−4 0
subject to
m
∑ (αi − αi0) = 0
i=1
and
αi , αi0 ∈ [0,C]. (7.26)
However: We can also solve the following dual QP (version II):
Maximize
m m m
−1/2 ∑ (αi − αi0)(α j − α 0j )(xix j )2 − ε ∑ (αi + αi0) + ∑ yi(αi − αi0)
i, j=1 i=1 i=1
subject to
m
∑ (αi − αi0) = 0
i=1
and
αi , αi0 ε[0,C]. (7.27)
Why is it better to solve the dual QP in version II, rather than solving
it in version I above? The reason is: We have
Example 7.2
272 7 Motion Correlation and Tracking
Remark 7.2
An important point to note is that the last observation also holds for
the step of computing w and d, from the values αi and αi0 . You may
wish to check this in the formulas defining the values for w and d
(equation 7.17 and equation 7.19). (End of Remark 7.2)
Hence, equation 7.27 gives the final method for finding a regres-
sion function with support vector machines. Notice that the use of
kernels allows for finding non-linear regression functions as in fig-
ure 7.7. Equation 7.27 could be implemented in Matlab, e.g. by calling
the Matlab-function quadprog. However, some minor rearrangement
would be necessary to do this. Furthermore, the quadprog() module in
Matlab does not give very good results in this case. Only very small
data sets can be handled. A much better way to implement support
vector methods is to use specialized algorithms for solving the optim-
ization problem in equation 7.27. While this problem is indeed a quad-
ratic program, more efficient algorithms for solving this special type
of problem have been developed [1]. Such algorithms are provided in
dedicated libraries for support vector methods [2].
Example 7.3
288
286
284
282
280
278
141 142 143 144 145 146 147 148 149
288
286
284
282
280
278
141 142 143 144 145 146 147 148 149
Fig. 7.14: Support vector regression for motion correlation. Large gray
dots: training data. Small black dots: fluoroscopy data. Small gray dots:
points computed by support vector regression. Same data as in fig-
ure 7.13.
290
288
286
284
282
280
278
141 142 143 144 145 146 147 148 149
amplitude [mm]
286
284
282
280
278
0 20 40 60 80 100 120
time [s]
288
amplitude [mm]
286
284
282
280
278
0 20 40 60 80 100 120
time [s]
fiducial correlation output model points error
Fig. 7.16: Motion traces for the data in figure 7.13. Top row: polynomial
fitting. Bottom row: support vector regression.
Ultrasound probe
with pressure sensor
Ultrasound
probe
Phantom
260
250
240
cause inaccuracies. We see the rib cage in x-ray images, but the tu-
mor is not always visible. The rib cage moves during respiration. To
obviate the need for implanted markers, we can use two correlation
models, both of which are established before treatment:
• model 1: external marker motion to rib cage motion
• model 2: rib cage motion to target motion
We combine the two correlation models in a chain to infer the target
position in real-time (figure 7.20).
DRRs 4D CT
double correlation
Fig. 7.20: Double correlation. The motion of the rib cage is detected by
digitally reconstructed radiographs (DRRs). Before treatment the target
is marked in a series of CT scans, showing distinct respiratory states (4D
CT). DRRs are computed from a 4D CT data set. Live shots are matched
to the DRRs to detect the current respiratory state, and compute ground
truth location of the target.
Exercises
Exercise 7.1
For the example in equation 7.1, show that the least squares solution
obtained by direct differentiation of equation 7.2 with respect to b and
m is the same as the solution obtained from the normal equation.
Exercise 7.2
Looking at only two sample points x1 and x2 , explain why the quad-
ratic program
Minimize
w21 + w22
subject to
y1 (wx1 + d) ≥ 0
y2 (wx2 + d) ≥ 0,
(7.30)
Hints:
Represent the line by its parameters (w, d), where w is the normal
vector of the line. Then this representation of a line can be multiplied
by a factor λ , without changing the line. Thus (λ w, λ d) for non-zero
λ is the same line as (w, d). First show that it suffices to consider only
those lines having
y1 (wx1 + d) = 1,
y2 (wx2 + d) = 1.
(7.31)
280 7 Motion Correlation and Tracking
Now recall that for any line, we obtain the point-line distance by scal-
ing the normal vector w to unit length, so that the margin is given
by:
w d
x1 + (7.32)
||w|| ||w||
1
Hence the margin is simply ||w|| . To maximize the margin, we must
minimize the length of w.
Exercise 7.3
a) Write a program for the system in equation 7.11 with input points
as in figure 7.5, and visualize the resulting regression line. This
is an exceedingly simple support vector machine. Again, it will
not be a very effective support vector machine, but will serve to
illustrate the basics.
b) Show that the system in equation 7.16 is actually a quadratic pro-
gram. Thus, show that input parameters in the form H, f, A, b can
be set up for the system in equation 7.16. One seemingly simple
way to do this would be to set θi = αi − αi0 , and use the θi as vari-
ables, i.e. minimze 12 θ T Hθ + fT θ subject to Aθ ≤ b However, this
will fail, since we have the subexpression αi + αi0 in the objective
function. Hint: Use the 2m variables θ1 , . . . , θ2m by setting:
θi = αi − αi0 for i = 1, . . . , m
and
θi+m = αi0 for i = 1, . . . , m.
7.4 Double Correlation 281
Max f (x, y)
Subject to g(x, y) = c
The function f can be visualized via its level curves (figure 7.21).
The level curves are obtained as curves where f (x, y) = c for various
constants c. Likewise, we can visualize the constraint g(x, y) = c as
a curve. Looking again at the above optimization problem, two cases
are possible: the curves f (x, y) = c and g(x, y) = c can cross, or they
can meet tangentially. The key observation for Lagrange multipliers is
the following: f can attain a maximum under the constraint g(x, y) = c
only if the level curves meet as tangents.
The gradient vector
∂f ∂f T
∇f = , (7.35)
∂x ∂y
is orthogonal to the level curve. Thus, if f and g meet tangentially,
their gradient vectors must be an α-multiple of each other for a non-
zero α. Therefore, we have
282 7 Motion Correlation and Tracking
∇ f = α∇g. (7.36)
α is called Lagrange-multiplier. To summarize, we can state the so-
called Lagrange system for the above optimization problem in the fol-
lowing form:
∂f ∂g
=α
∂x ∂x
∂f ∂g
=α
∂y ∂y
g(x, y) = c
∇α Λ (x, y, α) = 0 (7.39)
is equivalent to g(x, y) = 0.
(Here we have included the constant c into the definition of the func-
tion g, so that we have replaced the original constraint g(x, y) = c
by the constraint g(x, y) = 0). The compact version of this system,
namely the function Λ has a name. It is called the Lagrangian of the
given optimization problem, or also the Lagrange function.
Now consider the optimization problem:
Minimize
f (x, y) = x2 + y2
subject to
g(x, y) = x + y = 2.
a) Find the Lagrange function for this system
b) Find the optimum with the help of the Lagrange function
Notice that the theory of Lagrange multipliers extends to the case of
inequality constraints (see e.g. [3, 4]).
Exercise 7.5
Max c1 x + c2 y
subject to
a11 x + a12 y = b1
a21 x + a22 y = b2 (7.40)
284 7 Motion Correlation and Tracking
Maximize
f (x1 , . . . , xn )
subject to
g1 (x1 , . . . , xn ) = 0
..
.
gm (x1 , . . . , xn ) = 0, (7.41)
we must use m Lagrange multipliers α1 , . . . , αm .
Here, the tangent condition used above (gradient vectors of f and g
aligned) is not as obvious. In fact, for multiple constraints, the tangent
condition requires that the gradient vector of f is a linear combination
of the gradient vectors of the g-functions.
Specifically, in this case, the Lagrange function is
m
Λ (x1 , . . . , xn , α1 , . . . , αm ) = f (x1 , . . . , xn ) − ∑ αk gk (x1 , . . . , xn ).
k=1
(7.42)
At the optimum, the partial derivatives of the Lagrange function with
respect to the variables xi and αi must vanish.
a) Set up the Lagrange function Λ for the optimization problem in
equation 7.40. (You will need two Lagrange multipliers α1 , α2 ,
for the two constraints.) Take derivatives of Λ with respect to the
Lagrange multipliers αi and the primal variables x, y, and set the
resulting expressions to zero. Inserting these equations back into
Λ , show that Λ = bα, and AT α = c. Compare to the dual of a
linear program defined in exercise 6.2.
7.4 Double Correlation 285
Minimize
ww + ξ1 + ξ10 + ξ2 + ξ20 . . .
subject to
y1 − (wx1 + d) = ε + ξ1
(wx1 + d) − y1 = ε + ξ10
y2 − (wx2 + d) = ε + ξ2
(wx2 + d) − y2 = ε + ξ20
...,
taking derivatives with respect to the slack variables ξi , ξi0 as well.
This is a slight modification of the primal quadratic program for
regression in equation 7.15. Compare your result to the dual quad-
ratic program for regression in equation 7.16.
Exercise 7.6
Summary
Notes
ory motion tracking with registration (see e.g. [17, 18]), or by com-
bining motion correlation with direct registration.
A library for support vectors is available from [2], see also [1, 19, 20].
Additional background on normal equations and least squares prob-
lems is given in [21]. The text book by C. Bishop [22] discusses learn-
ing methods and regression.
References
1 patient
robot
0.5
position
−0.5
−1
0 1 2 3 4 5 6
t
time [sec]
The figure shows an ideal situation, with a perfect sine wave. At time
point t shown in the figure, the position of the target and the robot
position are known. Time is discretized, i.e. at fixed time intervals (i.e.
every 50 milliseconds) we must send a motion command to the robot.
By the time we have computed the target position and the robot has
reached its new position, the target will have moved. This time lag
291
292 8 Motion Prediction
1 + patient
o
robot
0.5
position
−0.5
−1
0 1 2 3 4 5 6
t + δt
time [sec]
124
122
120
position
118
116
114
0 20 40 60 80 100
time [sec]
110
100
90
position
80
70
60
50
40 45 50 55 60 65
time [sec]
Fig. 8.4: Cardiac motion and respiratory motion, measured with ultra-
sound. Typically, the frequency of cardiac motion is four times larger
than the frequency of the respiratory signal alone. The amplitude of res-
piratory motion is larger than the amplitude of cardiac motion at this
target location.
−0.5
−1
0 1 2 3 4 5 6
time [sec]
e∗ (n + 1) = 2e(n) − e(n − 1) +
[2e(n) − e(n − 1) − 2e(n − 1) + e(n − 2)],
thus
e∗ (n + 1) = 4e(n) − 4e(n − 1) + e(n − 2). (8.9)
Overall, the MULIN algorithm of order 2 is now given by
p(n + 1) =
p(n) + 4[p(n) − r(n)] −
4[p(n − 1) − r(n − 1)] +
[p(n − 2) − r(n − 2)]. (8.10)
We see that MULIN of higher order looks more deeply into the past
of the signal, as it uses the older values p(n − 1), p(n − 2).
MULIN can be generalized to order k (see exercise 8.1 at the end of
this chapter).
When deriving MULIN, we have always assumed that the value of p
for the immediate next time step p(n + 1) must be predicted. Thus our
prediction horizon was a single time step. However, it may also be
necessary to predict several time steps into the future.
To look further into the signal’s future, we would predict values
p(n + s) (8.14)
we use the same method as above, but we also compute p∗ (n + s − 1)
and then set
Fig. 8.6: MULIN predictor for a motion signal. Thick solid line: patient
position. Dashed line: robot position, dotted line: predictor output for
MULIN order two.
Figure 8.6 shows an example for MULIN of order two. The dotted
line shows the predictor output. Notice that overshooting occurs at the
maxima and minima.
298 8 Motion Prediction
en = yn − yLMS
n (8.16)
yLMS
n+1 = w y
n n (8.17)
wn+1 = wn + µen yn (8.18)
error
w1 w2 w
Fig. 8.7: Gradient descent
yLMS
1 = · · · = yLMS
M = y0 (8.25)
with
en = yn − ynLMS
n (8.29)
nLMS
yn+1 = wn yn (8.30)
en
wn+1 = wn + µ (8.31)
yn
If you compare the definitions in equation 8.18 and equation 8.31, you
will notice that a multiplication is simply replaced by a division. How
can this be right? It may seem astonishing that we can simply replace
this, and still obtain a valid definition. And how can this division be a
normalization?
The riddle is solved if we look again at the multi-dimensional version
of LMS: Here, a multi-dimensional version is obtained by replacing
the definition in equation 8.27 by
µ yn − ynLMS
n un
wn+1 = wn + 2
, n ≥ M. (8.32)
kun k2
Thus, we simply divided by the square of the magnitude of our in-
put signal, and this was done to normalize! After normalization, it
becomes easier to choose the value µ, especially in cases where
the magnitude of our input changes rapidly. Of course, in the one-
dimensional-case, the division will cancel out, so that we obtain a
derivation for the one-dimensional nLMS method.
Note that it is necessary to add a constant α (typically 10−4 or smaller)
to the denominators of equations 8.31 and 8.32 to avoid division by
zero.
The next section will deal with a principle similar to LMS, applied to
the bands of a wavelet decomposition.
But the resemblance of the cycles can be quite faint, even if the patient
breathes normally.
To motivate the use of wavelets in our application, we again look at
the example in figure 8.6. Here we applied the MULIN predictor to
an aperiodic, but nonetheless very smooth signal. MULIN gives good
results in this case. However, if the signal is not smooth, problems will
arise. Often, anatomical motion data are far from being smooth (fig-
ure 8.8). If a sudden change in direction occurs near a peak point along
the curve (see arrow in figure 8.8), then the predictor will assume that
the peak has been reached, and initiate a downwards motion, thereby
causing errors.
626
625.5
position
625
624.5
20 25 30 35 40 45 50 55 60 65 70
time [sec]
c0,n = yn (8.33)
1
c j+1,n = c j,n−2 j + c j,n , j≥0 (8.34)
2
c0,n = yn , n = 0, . . . , m. (8.36)
This will give the first line of the c-table. I.e. assume m = 8, so we
have computed the first line, giving
Likewise we now compute c1,3 , c1,4 , . . . c1,8 from the values in the first
line of the table. Our table will look like this:
We now compute
1
c2,3 = (c1,1 + c1,3 ). (8.39)
2
Likewise we obtain c2,4 , . . . , c2,8 . Finally, we obtain the full c-table
(for J = 3):
Notice that our c-table is not rectangular, i.e. table entries c1,0 , c2,0 ,
c2,1 , c2,2 , c3,0 etc. do not exist.
Remember, we set
Each line in the table is called a band. For J = 2, we obtain three bands
(table 8.2).
0
2.5 3 3.5 4 4.5 5
Decomposition (W1, W2, W3, c)
6
W1
3
0
2.5 3 3.5 4 4.5
6
W2
position
0
2.5 3 3.5 4 4.5
6
W3
3
0
2.5 3 3.5 4 4.5
6
c
3
0
2.5 3 3.5 4 4.5
time [sec]
Fig. 8.9: Wavelet decomposition for an input signal (with J = 3). Top
row: input signal. Bottom rows: Wavelet bands W1 , W2 , W3 and c .
W1 ,W2 ,W3 and c (four bottom rows). Notice that the input signal is
noisy when compared to the c-band in the decomposition. Thus, a
standard application of wavelet decomposition is noise removal in sig-
nal processing. Definitions 8.33 to 8.35 are just one way of defining
a wavelet decomposition. There are many other ways. The specific
decomposition introduced above is called the à trous wavelet decom-
position. Other types of wavelets are discussed in the exercises at the
end of this chapter, and this will also explain why we chose the à trous
wavelets for our application.
Equation 8.45 is the key to the predictor, to be defined next. The signal
is first decomposed into several bands. Then we predict the next values
for the individual bands separately.
We define a predictor in the following way:
8.3 Wavelet-based LMS Prediction 307
J
yn+1 = wJ+1 cJ,n + vJ+1 cJ,n−2J + ∑ w jW j,n + v jW j,n−2 j (8.45)
j=1
(8.46)
and
W̃1,n−1 W̃2,n−1 . . . W̃J−1,n−1 W̃J,n−1 c̃J,n−1
.. .. .. .. .. ..
B= . . . . . . (8.50)
W̃1,n−M W̃2,n−M . . . W̃J−1,n−M W̃J,n−M c̃J,n−M
Note that the entries of the matrix B are themselves vectors, namely
1-by-2 vectors. But this does not change the fact that the following
equation is simply a linear system of equations!
308 8 Motion Prediction
w̃1 yn
B ... = ...
(8.51)
w̃J+1 yn−M+1
This system has the 2J + 2 variables w1 , . . . , wJ+1 and v1 , . . . , vJ+1 and
a number M of equations. If M is larger than 2J + 2, then the system
is over-determined.
Equation 8.51 is the training equation for w1 , . . . , wJ+1 and v1 , . . . , vJ+1 .
Each matrix line in B corresponds to one training sample, and we thus
have M training samples. The training samples stem from the signal
history y0 , y1 , . . . , yn . As we go down the matrix lines, we step back-
wards in the signal history by one sample for each matrix line.
It remains to discuss how the system in equation 8.51 should be
solved. Naturally, one could use Gaussian elimination. However, over-
determined systems of this type can be solved with advantage by nor-
mal equations, so here again the method from the previous chapter can
be used. This will fit a regression line in the least-mean-square sense
to the data points represented by the matrix lines. Obviously, we can
also use support vectors for this regression. There are several ways to
do this, and we will discuss this possibility below.
But before we do this, we will look at two issues related to the im-
plementation and practical improvements of the wLMS algorithm.
Remember that the original version of the LMS algorithm (without
wavelets) used a smoothing parameter µ. Indeed we have used such a
parameter for MULIN as well.
To emphasize the fact that we had set up the B-matrix from the de-
composition data with the indices n − 1, . . . , n − M, we will add these
indices to the notation. Hence, write
B(n−1,...,n−M) (8.52)
for the above B-matrix in equation 8.51.
Likewise we extend the notation for the weight-vector, giving
(n−1,...,n−M)
w̃1 yn
B(n−1,...,n−M) ... = ... .
(8.53)
w̃J+1 yn−M+1
8.3 Wavelet-based LMS Prediction 309
(new)
w̃1
..
. =
w̃J+1
(n−2,...,n−M−1) (n−1,...,n−M)
w̃1 w̃1
(1 − µ) ... + µ ...
.
w̃J+1 w̃J+1
Instead of this latter definition, with just two weight variables w and
v per value of j, we use a series of such variables (for each j), now
called
(0) (a )
wj ,...,wj j .
Thus we obtain
(0) (1)
yn+1 = wJ+1 cJ,n + wJ+1 cJ,n−2J ·1 + ... +
(a
J+1 −1)
wJ+1 cJ,n−2J ·(aJ+1 −1) +
J
∑ wJ+1W j,n−2 j ·0 + w j
(0) (1)
W j,n−2 j ·1 + ... +
j=1
(a j −1)
wj W j,n−2 j ·(a j −1) . (8.56)
the training data after a new sample has been added. It may also be
necessary to “un-learn” single samples.
The RMS error thus gives an average over the squared distances
between the signal itself and the value returned by the predictor. N
is the total number of time steps.
This simple measure provides an absolute error, but does not allow for
much comparison. Thus, normalization is needed. In order to normal-
ize the RMS error, we divide by the error which would have occurred
if we had not applied any prediction at all. Stated in simpler terms,
we compute the RMS-distance between the patient curve and the pre-
dicted curve (figure 8.6), then the RMS-distance between the patient
curve and the delayed curve, then divide the two values. This yields
the relative RMS error. With yδ denoting the delayed curve:
314 8 Motion Prediction
s
1 N 2
RMS(y, yδ ) = ∑ yi − yδ ,i
N i=1
(8.60)
RMS(y, yPRED )
RMSrel = (8.61)
RMS(y, yδ )
In this way, the improvement obtained by prediction can be expressed
as a percentage of reduction in terms of RMS error. Thus, if success-
ful, prediction reduces RMS error by a certain percentage.
Relative RMS error is not the only measure for evaluating motion
prediction for robots. Maximum absolute error can be an important
measure as well. For technical reasons, it is often useful to limit the
so-called jitter of the predictor output signal. The jitter measures dis-
tances between two consecutive (discrete) points along the signal:
1 1 N PRED
J(y PRED
)= ∑ yi − yPRED
N ∆t i=0 i+1 (8.62)
1
Here, the factor ∆t normalizes with respect to the spacing between the
samples. Thus the jitter value J measures the average distance (per
second) traveled by the robot following the predictor output.
Notice that the measures just described may be in conflict with one
another. Respiration patterns vary considerably among patients. For
patients with large lung tumors, the patterns can become very irregu-
lar. Unfortunately, for lung tumors, the respiratory motion excursion
is particularly large. As noted above, no single algorithm works best
for all patients. For further empirical evaluation, and for the imple-
mentation of a fast-lane method, histograms are very useful. Given a
data base of test cases, a histogram plots the number of data base cases
(x-axis) versus the percentage gains (relative RMS) obtained with the
individual algorithm (y-axis). An example of such a histogram for a
number of predictors is shown in figure 8.10. The figure shows both,
RMS error and jitter. Here six algorithms are shown. The RLS method
is a variant of the LMS method, whereas EKF denotes the extended
Kalman filter [2] [6], applied to motion prediction.
8.5 Fast-Lane Methods and Performance Measures 315
1
0.9 0.8487
0.7768
0.8
0.7796
0.7
0.6
content
0.5789
0.5
0.4704
0.4
0.3783
0.3
0.2
0.1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
relative RMS
1
MULIN
0.9 nLMS
RLS
0.8 wLMS
EKF
0.7
SVRpred
0.6
content
0.5
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3 3.5 4
relative jitter
Table 8.3 states results of relative RMS errors for a large database of
respiratory traces from patient treatments, containing data from 304
full treatment sessions [7].
Table 8.3: Results of relative RMS errors obtained with different pre-
diction methods
It turns out that wLMS gives the best results in terms of RMS for over
60 % of the data base cases. MULIN gives best results on 28 % of
316 8 Motion Prediction
Fig. 8.11: Grid search in parameter space for the wLMS algorithm
8.5 Fast-Lane Methods and Performance Measures 317
Exercises
and
y1 − y2 y −y y −y
d1 = √ m.
√ , d2 = 3√ 4 , ..., dm/2 = m−1 (8.65)
2 2 2
Let T (y) = (a1 , ..., am/2 , d1 , ..., dm/2 ) denote the transformed signal.
The transformation of the signal y thus obtained is called the Haar
transform. The signal T (y) has the same length m as the input signal
y.
a) Show that the transformation thus defined is invertible, i.e. there is
a mapping I(y), such that I(T (y)) = y.
b) Show that T preserves the energy of the input signal, i.e. the scalar
product yy is equal to T (y)T (y).
c) Define and implement a ’denoising’ scheme based on the Haar
transform. Your denoiser should work as follows: Transform the
signal under T , set to zero all values in T (y) below a fixed
threshold γ, then transform the thresholded signal back under I.
d) Define a synthetic signal, add random noise, and derive a threshold
γ such that 99.9 percent of the energy of the original signal is pre-
served.
The Haar transform T was defined in the previous exercise. The defin-
ition of Haar wavelets is based on the Haar transform.
a) Iterate the construction leading to the Haar transform T , giving
a multi-level transform. Hint: for level 2, assume m is divisible
by 4. Construct the second level transform T (2) by transforming
only the values a1 , ..., am/2 , but leaving the values d1 , ..., dm/2 at
the end of the signal T (y) unchanged. This yields a signal T (2) y =
(a(2) , d(2) , d(1) ). Here, the lengths of the strings a(2) and d(2) are
m/4, while d(1) has length m/2 .
b) Define a component-wise multiplication for signals, i.e. yz =
(y1 z1 , ..., ym zm ). Show that the transform T can be represented by
this component-wise multiplication with so-called Haar wavelets
and Haar scaling functions. Hint: The Haar wavelets and Haar
8.5 Fast-Lane Methods and Performance Measures 319
The Daubechies wavelets are defined in much the same way as the
Haar wavelets. However, four adjacent values from the original signal
are used to compute two values ai and di . Specifically, we have
√ √ √ √
1+ 3 3+ 3 3− 3 1− 3
a1 = y1 √ + y2 √ + y3 √ + y4 √ . (8.66)
4 2 4 2 4 2 4 2
The values defined in equations 8.33 to 8.35 resemble the definitions
in equations 8.64 and 8.65. In both cases, an average and a difference
of two values is used to define the new values. Thus we could rename
the à trous method and call it the ‘à trous Haar’ method. In fact the
so-called Haar-function is the basis of both definitions: The definition
in equations 8.33 to 8.35 and the definition in equations 8.64 and 8.65.
(The Haar function is defined as the step function with values -1 in the
interval [0, 1/2), +1 in the interval [1/2, 1), and 0 elsewhere.)
Draw a graph of a Daubechies wavelet, and discuss why Daubechies
wavelets (although highly effective in signal processing) would seem
less adequate for prediction than the à trous scheme in equations 8.33
to 8.35.
320 8 Motion Prediction
Summary
Motion prediction is used to overcome the time lag required for data
and image processing, as well as robot motion time, and collimator
motion time. Thus the command sequence sent to the robot will not
refer to the current position of the anatomical target, but to its future
position. The future position of the target can be extrapolated in time,
if the motion is stable and periodic. We obtain a straightforward, yet
very practical method simply by adding the error from previous time
steps. Classical methods from signal processing attempt to capture
unknown characteristics of an incoming signal, and can be applied to
signal prediction and feed-forward loops. In practice, the convergence
radius of the signal processing methods for our application is rather
small, and better results are obtained with more advanced methods,
such as wavelet decomposition and support vector regression. Wave-
let decomposition filters the input signal into several bands. For each
band, a least squares method estimates parameter values for describ-
ing the curve corresponding to this band. Support vector regression
can be applied in several ways to the problem of prediction. Firstly,
we can simply set up samples (t, y(t)), or samples containing sev-
eral points of the form (t, y(t)) with equidistant time spacing. But
secondly, we can apply support vector regression to the bands in a
wavelet decomposition. Finally, fast online methods are particularly
suitable for the problem, since typically only a single new sample is
added at each time step. Measurements for evaluating predictors can
themselves be used in new predictors, which gives rise to a fast-lane
scheme.
Notes
Kalman filtering and RLS were developed in the same context. [2] is
a text book on filtering methods from signal processing. [15] gives a
survey of time series prediction.
The combination of correlation and prediction for moving anatomic
targets was suggested in 2000 [9]. It is in routine clinical use since
2002 [16]. The MULIN method was specifically developed for respir-
ation tracking [17]. Newer methods for respiration prediction include
methods based on kernel adaptive filtering [18].
The wavelet-based predictor is based on the à trous wavelet decom-
position [19]. The term à trous refers to the multi-level computation
series described in the chapter. The authors in [20] have turned this
general decomposer into a predictor, and this was applied to respira-
tion prediction in [21]. [8] is a text book with an introduction to dis-
crete wavelet decomposition.
As we have seen, support vector methods require several external
parameters to be chosen beforehand. One can estimate and modify
the values of the parameters during the application, via fast-lane meth-
ods (see also [5]). Probabilistic methods from machine learning, also
based on kernels, can also be applied to estimate the parameters
[22, 23].
References
325
326 9 Motion Replication
Operation site
Handle
Replicating robot
Tremor data
122
120
118
116
x−position
114
112
110
108
106
−2 0 2 4 6 8 10 12 14 16 18
time [sec]
amplitude
298
296
−5 0 5 10 15 20
Decomposition
2
W1
0
−2
−5 0 5 10 15 20
2
W2
0
−2
−5 0 5 10 15 20
300
c
298
296
−5 0 5 10 15 20
time [sec]
Fig. 9.3: Wavelet decomposition of tremor data, original signal in the top
row, recorded with infrared tracking.
9.2 Forces
measured force, and after this comply with the acting force by mov-
ing to a new position (figure 9.4). During this process, we must com-
pensate for gravity from a payload, since gravity could be mistaken
by the robot as a force exerted by the operator, pushing downwards.
One example for the need for force feedback is suturing in cardiac
surgery. If the robot pulls too hard while tying a knot, tissue can be
damaged. On the other hand, the knots should be tightened as much
as possible. It is then difficult to measure excessive forces.
Force sensor
Effector
Gravity
shoulder
Suppose our goal is to apply a specific force at the tool tip of the
replicator robot. What are the torques needed a the individual joints
of the robot?
Clearly, in practice, we would need sensors in this application. Non-
etheless, it is useful to understand the connection between forces ap-
plied at the tip of the robot, and torques acting on the joints. We again
consider the planar two-link manipulator (with revolute joints), in fig-
ure 4.1 and we again assume equal link lengths, i.e. l1 = l2 = 1. Before
we begin, we must define the term joint torque.
f
y l
θ1 x
nitude of the force is given by the length of f, and the direction of the
force is the direction of f.
Here, f is a 3D vector with coordinates
fx
f = fy .
(9.1)
fz
Since our stick rotates in the drawing plane we have fz = 0.
Let r be a vector of length l, pointing along our stick. Then applying
a force f to the tip of the stick will result in a torque at the joint. We
define the joint torque t by
t = r × f. (9.2)
Hence, t is simply defined as the cross-product of r and f. Since r and
f are in the drawing plane, t will be orthogonal to the drawing plane.
Hence, t coincides with the axis of rotation. The length τ of t is the
magnitude of our torque, i.e.
τ = ||t||. (9.3)
Our next goal is to analyze the two-link manipulator with respect to
forces and joint torques.
y y r1 r2
x p x p
Place vectors r1 and r2 along the two links, as shown in figure 9.7.
Then,
cos(θ1 )
r1 = , (9.4)
sin(θ1 )
9.2 Forces 333
and
cos(θ1 + θ2 )
r2 = . (9.5)
sin(θ1 + θ2 )
We apply a force to the tip (point p). Following our definition of joint
torque, the torque at joint 2 is then simply given by
t2 = r2 × f. (9.6)
The tip p of our robot is obtained simply by adding the two vectors r1
and r2 .
We now think of the vector r1 + r2 as a single stick (in the sense of
our above definition of torque). This stick connects the origin and the
point p. Thus the joint torque at joint 1, caused by f, is simply
t1 = (r1 + r2 ) × f. (9.7)
This last equation again follows directly from our definition of joint
torques.
Now both vectors t1 and t2 are orthogonal to our robot plane (i.e. the
x-y-plane), hence they are of the form
0
t = 0 . (9.8)
τ
Again τ denotes the magnitude of our torque. We consider the mag-
nitudes of the two torques (acting at joints 1 and 2) separately, and we
will call these magnitudes τ1 and τ2 . Inserting our definitions for r1
and r2 , we see that
τ = A · B f, (9.12)
where
−l1 s1 − l2 s12 l1 c1 + l2 c12
A= . (9.13)
−l2 s12 l2 c12
We see that the joint torques are directly proportional to the arm
lengths. Thus, typically, for a six-joint robot, the torques for joints
2 and 3 are the dominant torques.
Remark 9.1
Remark 9.2
B c12 −s12 E
f= · f. (9.14)
s12 c12
Inserting this into equation 9.12, we see that
c12 −s12 E
τ = A· · f. (9.15)
s12 c12
We now evaluate the product
c12 −s12
A· (9.16)
s12 c12
and obtain
c12 −s12 l1 s2 l2 + l1 c2
A· = . (9.17)
s12 c12 0 l2
z0
r
f
x0 -y0 -plane
What if f is not planar? In this case, t is also aligned with the joint axis
z = z0 , and its length τ is the length of the projection of r × f onto the
joint axis.
τ = (r × f)z
f
z0
Thus,
τ = (r × f)z (9.19)
and
In the following theorem, we will do the same for the general case of
n revolute joints and three-dimensional forces.
Theorem 9.1
Proof:
pi-1 pn
zi-1
r = (pn - pi-1)
f
y0
z0 p0 x0
Remark 9.3
Jv = z0 × (pn − p0 ) z1 × (pn − p1 ) . . . zn−1 × (pn − pn−1 ) , (9.28)
and
9.3 Joint Torques and Jacobi-Matrices 339
Jω = z0 z1 . . . zn−1 . (9.29)
Thus
z0 × (pn − p0 )
z1 × (pn − p1 )
T
Jv = .. . (9.30)
.
zn−1 × (pn − pn−1 )
Comparing this to equation 9.21 we see that indeed
τ1 z0 × (pn − p0 )
.. .. T
. = . · f = Jv f. (9.31)
τn zn−1 × (pn − pn−1 )
We can thus express the relationship between forces at the effector and
torques at the joints by the transpose of the upper part of the geometric
Jacobian. The following remark shows that the lower part Jω of the
geometric Jacobian allows for computing joint torques from torques
acting upon the effector.
Remark 9.4
We can not only apply a force to the effector, but we can also apply
torque to the effector. Thus, we can induce joint torque in two different
ways: (1) A force is applied to the effector. (2) A torque acts on the
effector.
In case (1) we push the effector with a finger tip. In case (2) we grasp
and twist the effector. Both cases result in joint torque.
Let m = (mx , my , mz )T be a torque, applied to the effector. Then the
induced joint torque is
t = (mz)z. (9.32)
The length of t is again the magnitude τ of the torque
340 9 Motion Replication
τ = mz. (9.33)
As a result of this definition, we obtain t as an orthogonal projection
of the effector torque vector m onto z = z0 (see figure 9.11).
m = (mx , my , mz )T
z = z0
y0
x0
Example 9.1
Table 9.1: Maximum joint torques (in Nm) for a force vector f of 50 N
acting at the effector. Positions of the effector ranging over a point grid
in the work space, and directions of f defined as grid points on a hemi-
sphere.
342 9 Motion Replication
Table 9.1 lists the maximum joint torques at each joint. We see that the
torques at the last joint (joint six for the six-joint robot and joint seven
for the seven-joint robot) are always zero. This is due to the orientation
of the last joint, and to the fact that we do not apply torque to the
effector. To allow for comparision, we adjusted the DH-parameters
in such a way that the total arm length is the same for both robots
(1251mm). The DH-parameters are the same as in table 3.3.
In practice, gear ratios should be as small as possible. For typical
small-size and light-weight motors, the maximum motor torques range
between 0.5 Nm and 1.2 Nm. Thus, from table 9.1 we see that a gear
ratio of 1:60 suffices for forces of up to 50 N, depending on the max-
imum motor torque.
Exercises
Exercise 9.1
Extend the static force equation in equation 9.10 for the planar two-
link manipulator to a planar three-link manipulator. Compare to the
force equation obtained by transposing the Jacobian.
Exercise 9.2
Summary
Notes
References
10.1 Radiosurgery
347
348 10 Applications of Surgical Robotics
Step 3. The relative position of the tumor with respect to the frame
coordinates is calculated from the tomography.
Step 4. The frame base remains attached to the patient’s head, and is
additionally attached to the irradiation device.
Step 5. With the double attachment from step 4, we find the exact
spatial reference. The patient couch is repositioned such that
the target is centered at the isocenter of the treatment beams
(typically beams of photon radiation).
Step 6. Treatment
There are two types of radiosurgical systems:
• static beam systems (GammaKnife)
• moving beam systems (LINAC systems, figure 10.1)
gantry rotation
gantry
patient table
table rotation
For the first type (Gamma Knife), beams are placed in a helmet, and
do not move during treatment. Thus the frame is attached to the pa-
10.1 Radiosurgery 349
tient’s head, and also to the treament couch. We thus have a static spa-
tial relationship between target and the radiation source(s). By con-
trast, LINAC systems move the radiation beam along circular arcs in
space, while the target is at the center of rotation. By moving the pa-
tient couch, different angles can be reached. Initially, radiosurgery was
limited to the head, due to organ motion in the chest or abdomen.
Robotic radiosurgery with a six-joint robot was introduced in 1995 to
obviate the need for attaching a head frame, and to obtain more kin-
ematic flexibility for beam motion paths [2]. Inverse treatment plan-
ning with linear programming (see chapter 6) computes beam direc-
tions and beam weights [3]. Planning for multi-leaf collimators is dif-
ferent when compared to planning for cylinder collimators, but the
general principles are similar: beam directions for the multi-leaf col-
limator are chosen in the same way as for cylindrical collimators. It is
then necessary to partition the cross-section of each beam into smal-
ler subsections to protect critical tissue in the vicinity of the target. As
for the case of beam selection, this is done by a randomized search.
Random geometric sub-patterns for beam cross-sections are gener-
ated, and a large set of ‘beam-lets’ thus chosen is submitted to linear
programming for weight computation. In the weighting step, a large
number of sub-patterns will receive weight zero. Such patterns are
removed, so that a re-sampling scheme is obtained for inverse plan-
ning [4]. The registration procedure relies on stereo x-ray imaging
and digitally reconstructed radiographs (DRRs) (see chapter 5). Thus,
no frame is needed. An on-line position correction is performed after
periodic x-ray imaging.
Respiratory motion compensation (based on motion correlation) was
then added to exisiting robotic radiosurgery systems [5]. This proced-
ure is called stereotaxic body radiosurgery (SBRS). Fast-lane methods
for motion prediction improve the accuracy of motion tracking in the
chest and abdomen. For correlation-based tracking, 1-5 gold markers
are implanted near the target. Marker-less treatment for lung tumors
uses the following procedure for image registration: Digitally recon-
structed radiographs are matched to live images (stereo x-ray image
pairs). The x-ray images of the target region are partitioned into small
patches. Typically, a 5x5 array of patches is used. Registration match
coefficients based on cross-correlation are evaluated for each pair of
350 10 Applications of Surgical Robotics
Cutting
plane
Implant Implant
cavity
distal femur
cutting
planes
bone to
remove for implant
placing
implant
femur
free-form
unicondylar
implant
Skin surface
Superior
Vena Cava
Sup. Pulmonary
Veins
Cut lines
Inf. Pulmonary
Veins
Inferior Vena Cava
10.5 Neurosurgery
ate robotic instrument was available, the following method for robotic
tumor removal with minimal user interaction would become realistic:
Via a small skin incision or opening, we approach the target. The in-
strument and the endoscopic camera are mounted to a snake-robot.
Having reached the target, we step through the cell array in the vi-
cinity of the target, where small regions are classified by imaging.
Tumor cells are removed step by step, while the instrument plans a
path around healthy cells. The control mode for this scheme is more
autonomous than the above scheme and is based on real-time endo-
scopic imaging.
Summary
Notes
Vinci system. This positioning system was later used in the develop-
ment of the Zeus system [22], in the context of a telesurgery paradigm
similar to the da Vinci system. The CyberKnife system was originally
developed to provide frameless radiosurgery in the head [23]. It in-
cluded an inverse planning system based on linear programming, and
first patients were treated with the system in 1994. It was later ex-
tended to incorporate full body radiosurgery and robotic motion com-
pensation [24]. Recently, other radiosurgical systems were equipped
with respiration tracking based on correlation and prediction, so that
robotic methods are now available for conventional radiation therapy
as well [25].
References
363
364 11 Rehabilitation, Neuroprosthetics and Brain-Machine Interfaces
Belt
Exosceleton
Cart
A B C D E F
G H I J K L
M N O P Q R
S T U V W X
Y Z 1 2 3 4
5 6 7 8 9
The so-called speller matrix (figure 11.2) further illustrates the use,
but also the limitiations of EEG for BMIs [8]. The matrix contains the
letters of the alphabet and the numbers 0, . . . , 9 in a quadratic table.
There are six rows and columns. The matrix is equipped with a small
light source for each letter and number. With the light sources, we can
highlight either an entire row of the matrix, or a column, as illustrated
in the figure. Assume a subject (or patient) wears a cap with EEG
electrodes. Our goal is to establish communication between the sub-
ject and a second person in another room, solely via the EEG signals.
We do this letter by letter.
Suppose now the subject wishes to transmit the letter ‘P’. We highlight
the rows and columns of the matrix in a random pattern. Thus, for
example, row 1 flashes for several seconds. After this, column 3 would
flash. The target letter ‘P’ is neither contained in row 1 nor in column
3. However, when we activate column 4 (containing ‘P’), we obtain
a detectable reaction in the EEG. More specifically, mismatch and
match conditions in visual tasks are typical data extractable from EEG
signals.
Similar EEG-based interfaces have been developed for wheel-chairs
and avatars in virtual reality environments. Here, EEG analysis de-
tects directional intention (left-right, forward-reverse) [9]. It should
11.2 Brain-Machine Interfaces 367
be noted that EEG analysis for this application is difficult, and long
training phases are needed for each individual patient. For example,
the EEG signal will look very different if the patient simply closes
the eyes. To address such difficulties, an EEG must be combined with
a system for artifact correction, e.g. an electro-oculogram (EOG) for
recording eye movement.
The inherent limitations of EEG raise the requirements for the control
intelligence of the robotic exoskeleton. One further shortcoming of
EEG in this context is the lack of an upstream sensory feedback path.
EMG
Nerves
normally Chest Electrodes
connecting muscles
the spine to
the hand
Robotic
arm
Missing arm
Finally, we can also listen to the local signal processing in the brain
via imaging. The first such method is fMRI. Here, the patient or sub-
ject must be inside an MRI scanner, i.e. cannot walk around. A second
method is optical coherence tomography (OCT). It is capable of de-
tecting local morphology in tissue to a depth of 1-2 mm with very high
resolution. We conjecture [17] that neuronal activity is associated with
small morphologic changes, i.e. neurons move or contract when they
fire. OCT allows for several extensions, which can be of help here.
Firstly, it is possible to combine OCT with microscopy, giving a tech-
nique called OCM (optical coherence microscopy). Secondly, OCT
and OCM devices can be miniaturized and placed at the tip of a thin
needle, giving a technique called endoscopic OCM.
Electrode
Pulse
Generator
Tissue damage in the brain (after injury, disease, or stroke) can cause
cognitive deficits and movement disorders. For example, Parkinson’s
disease can be treated with deep brain stimulation (DBS), if medic-
ation does not or no longer help (figure 11.4). With DBS, electrodes
are implanted into the deep brain, targeting specific nuclei respons-
ible for the symptoms. A pulse generator is implanted in the chest and
connected to the electrode. The electrode thus stimulates the affected
11.3 Steerable Needles 371
By taking more than two such needles, we can produce more complex
curvature patterns.
The neurosurgeon controls the position of the needle tip via acoustic
signals. The firing patterns of individual regions of the brain are char-
acteristic. This can be detected acoustically during needle insertion,
by reading from the electrodes at the tip of the needle and sending
the signals to a speaker. If we notice a slight deviation from the path,
we would have to retract the needle pair and insert it again. This is a
disadvantage of the above steering method with hollow needles.
This difficulty can be addressed by an alternative steering scheme (fig-
ure 11.6). Here we have a single highly flexible needle with a bevel tip.
The bevel tip causes the needle to move into the tissue along curved
paths. By rotating the needle appropriately during the inserting pro-
cess, we can steer the tip [18].
As noted above, MEAs can record signals from the deep brain. To
apply MEAs, we need to place the electrodes into the appropriate loc-
ation in the deep brain. Similarly, in the future, endosopic imaging
(e.g. OCT or OCM imaging) may be capable of recording such signals
from individual neurons. The goal of OCM in this case is to read the
entire firing pattern of three-dimensional clusters of neurons, which
would not be possible with MEAs [17].
Finally, it is possible to give sensory and haptic feedback from an
artificial hand back to the brain [19]. To this end the activity of the
peripheral nerve of a patient’s arm is recorded with flexible micro
electrode arrays. This enables the patient not only to send control sig-
nals to a prosthetic hand, but also to use the afferent nerve structures
to obtain shape information via the artificial hand.
11.3 Steerable Needles 373
Exercises
Outer needle
l1 z0
x0
Inner needle l2 = π
Exercise 11.1
References
377
378 A Image Modalities and Sensors
arm images. For older C-arms, the images must be calibrated, i.e. we
must place a calibration phantom into the beam path. The geometry of
the phantom is assumed to be known and we can thus determine the
projection geometry and distortion features.
MRI
Four-dimensional Imaging
In most cases, CT and MR are used in single image mode. This means
that a single image stack is taken at a fixed point in time. Breathing
artifacts can be reduced with respiratory gating. An alternative is four-
dimensional CT or MR imaging. A four dimensional respiratory CT
data set consists of 10 image stacks, where each stack corresponds
to one single time point in the respiratory cycle [2]. Applications of
four-dimensional imaging (e.g. in radiosurgery and cardiology) are
discussed in chapter 7.
Ultrasound
Infrared Tracking
line cameras
Magnetic Tracking
The physical principle for magnetic tracking is very different from in-
frared tracking. The position of small active coils is computed from
readings of a magnetic sensor array. The advantage of magnetic track-
ing over infrared tracking is the following: the marker (coil) does not
have to be visible in a camera image. Thus magnetic tracking is suit-
able for navigating laparoscopic instruments, endoscopes and biopsy
needle tips. However, the coil must be connected to the host unit by
a cable. In surgical applications, the cable must be sterile. Newer sys-
tems use small cable-less electromagnetic transponders, which can be
implanted into the target region before the operation, and will then
serve as landmarks during the operation.
Fig. A.2: Motor cortex mapping with fMRI (left), and robotic TMS
(right).
fMRI
References
b) For
nx ox ax px
ny oy ay py
M=
nz oz az
(B.1)
pz
0 0 0 1
the inverse is given by
nx ny nz −np
ox oy oz −op
M−1 =
ax ay az −ap
(B.2)
0 0 0 1
385
386 B Selected Exercise Solutions
There are several equivalent ways to find the YPR-angles given a ro-
tation matrix. One is the following:
Start from equation
nx ox ax
R = R(z, γ)R(y, β )R(x, α) = ny oy ay (B.3)
nz oz az
Then insert the elementary rotations (see equation 2.17) giving
cγ cβ cγ sβ sα − sγ cα cγ sβ cα + sγ sα nx ox ax
sγ cβ sγ sβ sα + cγ cα sγ sβ cα − cγ sα = ny oy ay (B.4)
−sβ cβ sα cβ cα nz oz az
We thus have
sβ = −nz (B.5)
We could solve this equation for β using the arcsin-function, i.e. we
could simply set β = arcsin −nz . However, arcsin only returns values
in the range [−π/2, π/2]. There may be solutions outside this range.
We could define a function asin2 similar to the atan2-function (see
equation 2.34). Instead we prefer to use atan2 since we have already
defined it. Since sin2 β + cos2 β = 1, we have
q
cβ = ± (1 − n2z ) (B.6)
Given equations B.5 and B.6, we are in a position to apply the function
atan2, i.e.
q
β1,2 = atan2(−nz , ± (1 − n2z )) (B.7)
After having solved for β , we can insert our value into equation B.4,
and obtain values for α and γ.
From the upper two elements of the first matrix column in equa-
tion B.4 we obtain the value for γ.
ny nx
γ = atan2 , (B.8)
cβ cβ
B Selected Exercise Solutions 387
The last two elements of the last matrix line in equation B.4 give α.
oz az
α = atan2 , (B.9)
cβ cβ
Notice that in equations B.8 and B.9 we assume cβ 6= 0. We must
handle the case cβ = 0 separately. In this case, we have two possibil-
ities:
1. β = π/2 and sβ = 1,
2. β = −π/2 and sβ = −1.
For the first possibility, the matrix equation B.4 simplifies to
0 sα−γ cα−γ nx ox ax
0 cα−γ −sα−γ = ny oy ay (B.10)
−1 0 0 nz oz az
Here we have used the trigonometric addition formulas to simplify
expressions (for example: cγ sα − sγ cα is written as sα−γ ).
Likewise, for the second possibility we obtain:
0 −sα+γ −cα+γ nx ox ax
0 cα+γ −sα+γ = ny oy ay (B.11)
1 0 0 nz oz az
In both cases the solutions for α and γ can again be found with the
atan2-function from the matrix equations B.10 and B.11. Notice that
in both cases the solutions for γ and α will not be unique. We can only
obtain expressions for α + γ or α − γ. The interpretation of this fact
is the following: at sβ = 0, we have a singularity, and the two angles
γ and α can compensate each other.
c) Consider
c + u2x (1 − c) ux uy (1 − c) − uz s ux uz (1 − c) + uy s
uy ux (1 − c) + uz s c + u2y (1 − c) uy uz (1 − c) − ux s =
uz ux (1 − c) − uy s uz uy (1 − c) + ux s c + u2z (1 − c)
nx ox ax
ny oy ay
nz oz az
(B.12)
oz − ay = 2sux
ax − nz = 2suy
ny − ax = 2suz
θ 0 (0) = 0 (B.15)
θ 0 (1) = 0 (B.16)
θ (0) = a1 (B.17)
B Selected Exercise Solutions 389
θ (1) = b1 (B.18)
and
c0 = a1 (B.21)
Combining equation B.15 and equation B.20 we have
c1 = 0 (B.22)
Hence
b1 = a1 + c2 + c3 (B.24)
Combining equation B.16 and equation B.20
c0 = a1
c1 = 0
c2 = 3(b1 − a1 )
c3 = −2(b1 − a1 )
(B.26)
0 ≤ pA (a) ≤ 1. (B.28)
Thus for each a, we have
Looking at only two sample points x1 and x2 , explain why the quad-
ratic program
Minimize
w21 + w22
Subject to
y1 (wx1 + d) ≥ 0
y2 (wx2 + d) ≥ 0
(B.31)
y1 (wx1 + d) = 1
y2 (wx2 + d) = 1
(B.32)
y1 (wx1 + d) = a
where a > 0, and
y2 (wx2 + d) = b
with b > 0 as well. Otherwise our line would cross one of the points
x1 , x2 , and then it can certainly not be a maximum margin separator
line. But now we can also assume a = b. Otherwise, if a were not
equal to b for a given line we could simply move this line slightly and
obtain a line with a yet bigger margin. But now, since a = b, we can
scale with a factor λ to have the above form equation B.32. But now
recall that for any line, we obtain the point-line distance by scaling the
normal vector w to unit length, so that the margin is given by:
w d
x1 + (B.33)
kwk kwk
To see why this is indeed the margin, recall that the distance of the
origin from a line with equation wx + d = 0 is d, if the normal vector
w has unit length, i.e. kwk = 1. Likewise distances of the points x1 , x2
from the given line can be calculated. But equation B.33 is the same
as
1
(wx1 + d) (B.34)
kwk
And by equation B.32 we have wx1 + d = 1. Hence the margin is
1
simply kwk . To maximize the margin, we must minimize the length of
w.
392 B Selected Exercise Solutions
−1 0 −1 0 (B.36)
and the first entry of the vector b is 0.
Now define the constraint matrix A and the vector b according to
B Selected Exercise Solutions 393
A b meaning
−1 0 −1 0 0 α1 ≥ 0
0 −1 0 −1 0 α2 ≥ 0
1 0 1 0 C α1 ≤ C
0 1 0 1 C α2 ≤ C
0 0 −1 0 0 α10 ≥ 0 (B.37)
0 0 0 −1 0 α20 ≥ 0
0 0 1 0 C α10 ≤ C
0 0 0 1 C α20 ≤ C
1 1 0 0 0 (α1 − α1 ) + (α2 − α20 ) ≤ 0
0
Fig. B.1: Regression with dual quadratic programming, result for two
points x1 = 1, x2 = 2, y1 = 1, y2 = 2
Finally, we define f as
ε − y1
ε − y2
f=
2ε
(B.38)
2ε
394 B Selected Exercise Solutions
Now call θ = quad prog(H, f, A, b). Your result should look like
figure B.1.
f (x, y) = x2 + y2
Subject to
g(x, y) = x + y = 2
a) Find the Lagrange function for this system.
b) Find the optimum with the help of the Lagrange function.
We have L(x, y, α) = f (x, y) − αg(x, y) = x2 + y2 + α(x + y − 2).
Taking 5L = 0, we have
x+y = 2
2x = α
2y = α
This is a linear system with the solution x = 1, y = 1, α = 2.
c1 x + c2 y (B.39)
subject to
a11 x + a12 y = b1
a21 x + a22 y = b2 (B.40)
B Selected Exercise Solutions 395
∂Λ
= a11 x + a12 y − b1 = 0 (B.42)
∂ α1
∂Λ
= a21 x + a22 y − b2 = 0 (B.43)
∂ α2
∂Λ
= c1 − α1 a11 − α2 a21 = 0 (B.44)
∂x
∂Λ
= c2 − α1 a12 − α2 a22 = 0 (B.45)
∂y
From the last two equations, we directly obtain
Λ = c1 x + c2 y (B.48)
From equations B.46 and B.47 we have
hence Λ = αb.
z0
l1
q
x0
l2
θ = atan2(py , px ) (B.57)
From first two lines in equation B.56 we obtain two cases:
If sin2 (θ ) > ε, we set
l1 = − sin(l2 ) − pz (B.60)
398 B Selected Exercise Solutions
References
i αi ai di θi
1 -90 0 0 θ1
2 0 a2 0 θ2
3 90 0 0 θ3
4 -90 0 d4 θ4
5 90 0 0 θ5
6 0 0 d6 θ6
Table C.1: DH-table for the six-joint elbow manipulator with spherical
wrist
cos (θi ) − sin (θi ) cos (αi ) sin (θi ) sin (αi ) cos (θi ) ai
i−1 M =
sin (θ i ) cos (θ i ) cos (α i ) − cos (θ i ) sin (α i ) sin (θ i ) ai
i
0 sin (α ) cos (α ) d
i i i
0 0 0 1
399
400 C Geometric Jacobian for the Six-Joint Elbow Manipulator
DH-matrices:
c1 0 −s1 0
0M =
s1 0 c1 0
1
0 −1 0 0
0 0 0 1
c2 −s2 0 c2 a2
1M
s2 c2 0 s2 a2
2
=
0 0 1 0
0 0 0 1
c3 0 s3 0
2M
s3 0 −c3 0
3 =
0 1 0 0
0 0 0 1
c4 0 −s4 0
3M
s4 0 c4 0
4
=
0 −1 0 d
4
0 0 0 1
c5 0 s5 0
4M
s5 0 −c5 0
5 =
0 1 0 0
0 0 0 1
c6 −s6 0 0
5M =
s6 c6 0 0
6
0 0 1d
6
0 0 0 1
C Geometric Jacobian for the Six-Joint Elbow Manipulator 401
Forward kinematics:
c1 0 −s1 0
0M =
s1 0 c1 0
1
0 −1 0 0
0 0 0 1
c1 c2 −c1 s2 −s1 c1 c2 a2
0M
s1 c2 −s1 s2 c1 s1 c2 a2
2
=
−s2 −c2 0 −s2 a2
0 0 0 1
c1 (−s2 s3 + c2 c3 ) −s1 c1 (c2 s3 + s2 c3 ) c1 c2 a2
0M
s1 (−s2 s3 + c2 c3 ) c1 s1 (c2 s3 + s2 c3 ) s1 c2 a2
3 =
−s c − c s
2 3 2 3 0 −s2 s3 + c2 c3 −s2 a2
0 0 0 1
m11 m12 m13 m14
0M
m21 m22 m23 m24
4
=
m m m m
31 32 33 34
0 0 0 1
where
m11 = −c1 c4 s2 s3 + c1 c4 c2 c3 − s1 s4
m12 = −c1 (c2 s3 + s2 c3 )
m13 = c1 s4 s2 s3 − c1 s4 c2 c3 − s1 c4
m14 = c1 (d4 c2 s3 + d4 s2 c3 + c2 a2 )
402 C Geometric Jacobian for the Six-Joint Elbow Manipulator
m21 = −s1 c4 s2 s3 + s1 c4 c2 c3 + c1 s4
m22 = −s1 (c2 s3 + s2 c3 )
m23 = s1 s4 s2 s3 − s1 s4 c2 c3 + c1 c4
m24 = s1 (d4 c2 s3 + d4 s2 c3 + c2 a2 )
m11 m12 m13 m14
0M
m21 m22 m23 m24
5
=
m31 m32 m33 m34
0 0 0 1
where
m11 = −c5 c1 c4 s2 s3 + c5 c1 c4 c2 c3 − c5 s1 s4 − c1 s5 c2 s3 − c1 s5 s2 c3
m12 = c1 s4 s2 s3 − c1 s4 c2 c3 − s1 c4
m13 = −s5 c1 c4 s2 s3 + s5 c1 c4 c2 c3 − s5 s1 s4 + c1 c5 c2 s3 + c1 c5 s2 c3
m14 = c1 (d4 c2 s3 + d4 s2 c3 + c2 a2 )
m21 = −c5 s1 c4 s2 s3 + c5 s1 c4 c2 c3 + c5 c1 s4 − s1 s5 c2 s3 − s1 s5 s2 c3
m22 = s1 s4 s2 s3 − s1 s4 c2 c3 + c1 c4
m23 = −s5 s1 c4 s2 s3 + s5 s1 c4 c2 c3 + s5 c1 s4 + s1 c5 c2 s3 + s1 c5 s2 c3
m24 = s1 (d4 c2 s3 + d4 s2 c3 + c2 a2 )
C Geometric Jacobian for the Six-Joint Elbow Manipulator 403
m31 = −c4 c5 c2 s3 − c4 c5 s2 c3 + s5 s2 s3 − s5 c2 c3
m32 = (c2 s3 + s2 c3 ) s4
m33 = −c4 s5 c2 s3 − c4 s5 s2 c3 − c5 s2 s3 + c5 c2 c3
m34 = −d4 s2 s3 + d4 c2 c3 − s2 a2
m11 m12 m13 m14
0M =
m21 m22 m23 m 24
6
m m m m
31 32 33 34
0 0 0 1
where
m11 = −c6 c5 c1 c4 s2 s3 + c6 c5 c1 c4 c2 c3 − c6 c5 s1 s4 − c6 c1 s5 c2 s3 −
c6 c1 s5 s2 c3 + s6 c1 s4 s2 s3 − s6 c1 s4 c2 c3 − s6 s1 c4
m12 = s6 c5 c1 c4 s2 s3 − s6 c5 c1 c4 c2 c3 + s6 c5 s1 s4 + s6 c1 s5 c2 s3 +
s6 c1 s5 s2 c3 + c6 c1 s4 s2 s3 − c6 c1 s4 c2 c3 − c6 s1 c4
m13 = −s5 c1 c4 s2 s3 + s5 c1 c4 c2 c3 − s5 s1 s4 + c1 c5 c2 s3 + c1 c5 s2 c3
m14 = −d6 s5 c1 c4 s2 s3 + d6 s5 c1 c4 c2 c3 − d6 s5 s1 s4 + d6 c1 c5 c2 s3 +
d6 c1 c5 s2 c3 + c1 d4 c2 s3 + c1 d4 s2 c3 + c1 c2 a2
m21 = −c6 c5 s1 c4 s2 s3 + c6 c5 s1 c4 c2 c3 + c6 c5 c1 s4 − c6 s1 s5 c2 s3 −
c6 s1 s5 s2 c3 + s6 s1 s4 s2 s3 − s6 s1 s4 c2 c3 + s6 c1 c4
m22 = s6 c5 s1 c4 s2 s3 − s6 c5 s1 c4 c2 c3 − s6 c5 c1 s4 + s6 s1 s5 c2 s3 +
s6 s1 s5 s2 c3 + c6 s1 s4 s2 s3 − c6 s1 s4 c2 c3 + c6 c1 c4
m23 = −s5 s1 c4 s2 s3 + s5 s1 c4 c2 c3 + s5 c1 s4 + s1 c5 c2 s3 + s1 c5 s2 c3
m24 = −d6 s5 s1 c4 s2 s3 + d6 s5 s1 c4 c2 c3 + d6 s5 c1 s4 + d6 s1 c5 c2 s3 +
d6 s1 c5 s2 c3 + s1 d4 c2 s3 + s1 d4 s2 c3 + s1 c2 a2
404 C Geometric Jacobian for the Six-Joint Elbow Manipulator
m31 = −c6 c4 c5 c2 s3 − c6 c4 c5 s2 c3 + c6 s5 s2 s3 − c6 s5 c2 c3 + s4 s6 c2 s3 +
s4 s6 s2 c3
m32 = s6 c4 c5 c2 s3 + s6 c4 c5 s2 c3 − s6 s5 s2 s3 + s6 s5 c2 c3 + s4 c6 c2 s3 +
s4 c6 s2 c3
m33 = −c4 s5 c2 s3 − c4 s5 s2 c3 − c5 s2 s3 + c5 c2 c3
m34 = −d6 c4 s5 c2 s3 − d6 c4 s5 s2 c3 − d6 c5 s2 s3 + d6 c5 c2 c3 − d4 s2 s3 +
d4 c2 c3 − s2 a2
Geometric Jacobian
" #
Jν
J=
Jω
m11 m12 . . . m16
m21 m22 . . . m26
J = .. .. . . . ..
. . .
m61 m62 . . . m66
where
C Geometric Jacobian for the Six-Joint Elbow Manipulator 405
m11 = d6 s5 s1 c4 s2 s3 − d6 s5 s1 c4 c2 c3 − d6 s5 c1 s4 − d6 s1 c5 c2 s3 −
d6 s1 c5 s2 c3 − s1 d4 c2 s3 − s1 d4 s2 c3 − s1 c2 a2
m12 = −c1 (d6 c4 s5 c2 s3 + d6 c4 s5 s2 c3 + d6 c5 s2 s3 − d6 c5 c2 c3 +
d4 s2 s3 − d4 c2 c3 + s2 a2 )
m13 = −c1 (d6 c4 s5 c2 s3 + d6 c4 s5 s2 c3 + d6 c5 s2 s3 − d6 c5 c2 c3 +
d4 s2 s3 − d4 c2 c3 )
m14 = s5 (c1 s4 s2 s3 − c1 s4 c2 c3 − s1 c4 ) d6
m15 = −d6 (c5 c1 c4 s2 s3 − c5 c1 c4 c2 c3 + c5 s1 s4 + c1 s5 c2 s3 + c1 s5 s2 c3 )
m16 = 0
m21 = −d6 s5 c1 c4 s2 s3 + d6 s5 c1 c4 c2 c3 − d6 s5 s1 s4 + d6 c1 c5 c2 s3 +
d6 c1 c5 s2 c3 + c1 d4 c2 s3 + c1 d4 s2 c3 + c1 c2 a2
m22 = −s1 (d6 c4 s5 c2 s3 + d6 c4 s5 s2 c3 + d6 c5 s2 s3 − d6 c5 c2 c3 +
d4 s2 s3 − d4 c2 c3 + s2 a2 )
m23 = −s1 (d6 c4 s5 c2 s3 + d6 c4 s5 s2 c3 + d6 c5 s2 s3 − d6 c5 c2 c3 +
d4 s2 s3 − d4 c2 c3 )
m24 = s5 (s1 s4 s2 s3 − s1 s4 c2 c3 + c1 c4 ) d6
m25 = −d6 (s1 s5 c2 s3 + s1 s5 s2 c3 + c5 s1 c4 s2 s3 − c5 s1 c4 c2 c3 − c5 c1 s4 )
m26 = 0
406 C Geometric Jacobian for the Six-Joint Elbow Manipulator
m31 = 0
m32 = d6 s5 c4 s2 s3 − d6 s5 c4 c2 c3 − d6 c5 c2 s3 − d6 c5 s2 c3 − d4 c2 s3 −
d4 s2 c3 − c2 a2
m33 = d6 s5 c4 s2 s3 − d6 s5 c4 c2 c3 − d6 c5 c2 s3 − d6 c5 s2 c3 −
d4 c2 s3 − d4 s2 c3
m34 = s5 s4 (c2 s3 + s2 c3 ) d6
m35 = d6 (−c4 c5 c2 s3 − c4 c5 s2 c3 + s5 s2 s3 − s5 c2 c3 )
m36 = 0
m41 = 0
m42 = −s1
m43 = −s1
m44 = c1 (c2 s3 + s2 c3 )
m45 = c1 s4 s2 s3 − c1 s4 c2 c3 − s1 c4
m46 = −s5 c1 c4 s2 s3 + s5 c1 c4 c2 c3 − s5 s1 s4 + c1 c5 c2 s3 + c1 c5 s2 c3
m51 = 0
m52 = c1
m53 = c1
m54 = s1 (c2 s3 + s2 c3 )
m55 = s1 s4 s2 s3 − s1 s4 c2 c3 + c1 c4
m56 = −s5 s1 c4 s2 s3 + s5 s1 c4 c2 c3 + s5 c1 s4 + s1 c5 c2 s3 + s1 c5 s2 c3
C Geometric Jacobian for the Six-Joint Elbow Manipulator 407
m61 = 1
m62 = 0
m63 = 0
m64 = −s2 s3 + c2 c3
m65 = (c2 s3 + s2 c3 ) s4
m66 = −c4 s5 c2 s3 − c4 s5 s2 c3 − c5 s2 s3 + c5 c2 c3
Appendix D
Geometric Jacobian for the Seven-Joint
DLR-Kuka Robot
DH-table:
i αi ai di θi
1 -90 0 d1 θ1
2 90 0 0 θ2
3 -90 0 d3 θ3
4 90 0 0 θ4
5 -90 0 d5 θ5
6 90 0 0 θ6
7 0 0 d7 θ7
DH-matrices:
c1 0 −s1 0
0M =
s1 0 c1 0
1
0 −1 0 d
1
0 0 0 1
409
410 D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot
c2 0 s2 0
1M =
s2 0 −c2 0
2
0 1 0 0
0 0 0 1
c3 0 −s3 0
2M
s3 0 c3 0
3
=
0 −1 0 d3
0 0 0 1
c4 0 s4 0
3M
s4 0 −c4 0
4 =
0 1 0 0
0 0 0 1
c5 0 −s5 0
4M
s5 0 c5 0
5
=
0 −1 0 d
5
0 0 0 1
c6 0 s6 0
5M
s6 0 −c6 0
6 =
0 1 0 0
0 0 0 1
c7 −s7 0 0
6M =
s7 c7 0 0
7
0 0 1d
7
0 0 0 1
D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot 411
Forward kinematics
c1 0 −s1 0
0M
s1 0 c1 0
1
=
0 −1 0 d1
0 0 0 1
c1 c2 −s1 c1 s2 0
0M
s1 c2 c1 s1 s2 0
2 =
−s 0 c d
2 2 1
0 0 0 1
c1 c2 c3 − s1 s3 −c1 s2 −c1 c2 s3 − s1 c3 c1 s2 d3
0M
s1 c2 c3 + c1 s3 −s1 s2 −s1 c2 s3 + c1 c3 s1 s2 d3
3
=
−s2 c3 −c2 s2 s3 c2 d3 + d1
0 0 0 1
m11 m12 m13 m14
0M
m21 m22 m23 m24
4
=
m31 m32 m33 m34
0 0 0 1
where
m11 = c4 c1 c2 c3 − c4 s1 s3 − c1 s2 s4
m12 = −c1 c2 s3 − s1 c3
m13 = s4 c1 c2 c3 − s4 s1 s3 + c1 s2 c4
m14 = c1 s2 d3
412 D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot
m21 = c4 s1 c2 c3 + c4 c1 s3 − s1 s2 s4
m22 = −s1 c2 s3 + c1 c3
m23 = s4 s1 c2 c3 + s4 c1 s3 + s1 s2 c4
m24 = s1 s2 d3
m31 = −s2 c3 c4 − c2 s4
m32 = s2 s3
m33 = −s2 c3 s4 + c2 c4
m34 = c2 d3 + d1
m11 m12 m13 m14
0M
m21 m22 m23 m24
5 =
m m m m
31 32 33 34
0 0 0 1
where
m11 = c5 c4 c1 c2 c3 − c5 c4 s1 s3 − c5 c1 s2 s4 − s5 c1 c2 s3 − s5 s1 c3
m12 = −s4 c1 c2 c3 + s4 s1 s3 − c1 s2 c4
m13 = −s5 c4 c1 c2 c3 + s5 c4 s1 s3 + s5 c1 s2 s4 − c5 c1 c2 s3 − c5 s1 c3
m14 = d5 s4 c1 c2 c3 − d5 s4 s1 s3 + d5 c1 s2 c4 + c1 s2 d3
m21 = c5 c4 s1 c2 c3 + c5 c4 c1 s3 − c5 s1 s2 s4 − s5 s1 c2 s3 + s5 c1 c3
m22 = −s4 s1 c2 c3 − s4 c1 s3 − s1 s2 c4
m23 = −s5 c4 s1 c2 c3 − s5 c4 c1 s3 + s5 s1 s2 s4 − c5 s1 c2 s3 + c5 c1 c3
m24 = d5 s4 s1 c2 c3 + d5 s4 c1 s3 + d5 s1 s2 c4 + s1 s2 d3
D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot 413
m31 = −c5 s2 c3 c4 − c5 c2 s4 + s2 s3 s5
m32 = s2 c3 s4 − c2 c4
m33 = s5 s2 c3 c4 + s5 c2 s4 + s2 s3 c5
m34 = −d5 s2 c3 s4 + d5 c2 c4 + c2 d3 + d1
m11 m12 m13 m14
0M
m21 m22 m23 m24
6 =
m m m m
31 32 33 34
0 0 0 1
where
m11 = c6 c5 c4 c1 c2 c3 − c6 c5 c4 s1 s3 − c6 c5 c1 s2 s4 − c6 s5 c1 c2 s3 −
c6 s5 s1 c3 − s6 s4 c1 c2 c3 + s6 s4 s1 s3 − s6 c1 s2 c4
m12 = −s5 c4 c1 c2 c3 + s5 c4 s1 s3 + s5 c1 s2 s4 − c5 c1 c2 s3 − c5 s1 c3
m13 = s6 c5 c4 c1 c2 c3 − s6 c5 c4 s1 s3 − s6 c5 c1 s2 s4 − s6 s5 c1 c2 s3 −
s6 s5 s1 c3 + c6 s4 c1 c2 c3 − c6 s4 s1 s3 + c6 c1 s2 c4
m14 = d5 s4 c1 c2 c3 − d5 s4 s1 s3 + d5 c1 s2 c4 + c1 s2 d3
m21 = c6 c5 c4 s1 c2 c3 + c6 c5 c4 c1 s3 − c6 c5 s1 s2 s4 − c6 s5 s1 c2 s3 +
c6 s5 c1 c3 − s6 s4 s1 c2 c3 − s6 s4 c1 s3 − s6 s1 s2 c4
m22 = −s5 c4 s1 c2 c3 − s5 c4 c1 s3 + s5 s1 s2 s4 − c5 s1 c2 s3 + c5 c1 c3
m23 = s6 c5 c4 s1 c2 c3 + s6 c5 c4 c1 s3 − s6 c5 s1 s2 s4 − s6 s5 s1 c2 s3 +
s6 s5 c1 c3 + c6 s4 s1 c2 c3 + c6 s4 c1 s3 + c6 s1 s2 c4
m24 = d5 s4 s1 c2 c3 + d5 s4 c1 s3 + d5 s1 s2 c4 + s1 s2 d3
414 D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot
m31 = −c6 c5 s2 c3 c4 − c6 c5 c2 s4 + c6 s2 s3 s5 + s6 s2 c3 s4 − s6 c2 c4
m32 = s5 s2 c3 c4 + s5 c2 s4 + s2 s3 c5
m33 = −s6 c5 s2 c3 c4 − s6 c5 c2 s4 + s6 s2 s3 s5 − c6 s2 c3 s4 + c6 c2 c4
m34 = −d5 s2 c3 s4 + d5 c2 c4 + c2 d3 + d1
m11 m12 m13 m14
0M =
m 21 m22 m23 m24
7
m m m m
31 32 33 34
0 0 0 1
where
m11 = c7 c6 c5 c4 c1 c2 c3 − c7 c6 c5 c4 s1 s3 − c7 c6 c5 c1 s2 s4 − c7 c6 s5 c1 c2 s3 −
c7 c6 s5 s1 c3 − c7 s6 s4 c1 c2 c3 + c7 s6 s4 s1 s3 − c7 s6 c1 s2 c4 −
s7 s5 c4 c1 c2 c3 + s7 s5 c4 s1 s3 + s7 s5 c1 s2 s4 − s7 c5 c1 c2 s3 − s7 c5 s1 c3
m12 = −s7 c6 c5 c4 c1 c2 c3 + s7 c6 c5 c4 s1 s3 + s7 c6 c5 c1 s2 s4 + s7 c6 s5 c1 c2 s3 +
s7 c6 s5 s1 c3 + s7 s6 s4 c1 c2 c3 − s7 s6 s4 s1 s3 + s7 s6 c1 s2 c4 −
c7 s5 c4 c1 c2 c3 + c7 s5 c4 s1 s3 + c7 s5 c1 s2 s4 − c7 c5 c1 c2 s3 − c7 c5 s1 c3
m13 = s6 c5 c4 c1 c2 c3 − s6 c5 c4 s1 s3 − s6 c5 c1 s2 s4 − s6 s5 c1 c2 s3 − s6 s5 s1 c3 +
c6 s4 c1 c2 c3 − c6 s4 s1 s3 + c6 c1 s2 c4
m14 = d7 s6 c5 c4 c1 c2 c3 − d7 s6 c5 c4 s1 s3 − d7 s6 c5 c1 s2 s4 − d7 s6 s5 c1 c2 s3 −
d7 s6 s5 s1 c3 + d7 c6 s4 c1 c2 c3 − d7 c6 s4 s1 s3 + d7 c6 c1 s2 c4 +
d5 s4 c1 c2 c3 − d5 s4 s1 s3 + d5 c1 s2 c4 + c1 s2 d3
D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot 415
m21 = c7 c6 c5 c4 s1 c2 c3 + c7 c6 c5 c4 c1 s3 − c7 c6 c5 s1 s2 s4 − c7 c6 s5 s1 c2 s3 +
c7 c6 s5 c1 c3 − c7 s6 s4 s1 c2 c3 − c7 s6 s4 c1 s3 − c7 s6 s1 s2 c4 −
s7 s5 c4 s1 c2 c3 − s7 s5 c4 c1 s3 + s7 s5 s1 s2 s4 − s7 c5 s1 c2 s3 + s7 c5 c1 c3
m22 = −s7 c6 c5 c4 s1 c2 c3 − s7 c6 c5 c4 c1 s3 + s7 c6 c5 s1 s2 s4 + s7 c6 s5 s1 c2 s3 −
s7 c6 s5 c1 c3 + s7 s6 s4 s1 c2 c3 + s7 s6 s4 c1 s3 + s7 s6 s1 s2 c4 −
c7 s5 c4 s1 c2 c3 − c7 s5 c4 c1 s3 + c7 s5 s1 s2 s4 − c7 c5 s1 c2 s3 + c7 c5 c1 c3
m23 = s6 c5 c4 s1 c2 c3 + s6 c5 c4 c1 s3 − s6 c5 s1 s2 s4 − s6 s5 s1 c2 s3 + s6 s5 c1 c3 +
c6 s4 s1 c2 c3 + c6 s4 c1 s3 + c6 s1 s2 c4
m24 = d7 s6 c5 c4 s1 c2 c3 + d7 s6 c5 c4 c1 s3 − d7 s6 c5 s1 s2 s4 − d7 s6 s5 s1 c2 s3 +
d7 s6 s5 c1 c3 + d7 c6 s4 s1 c2 c3 + d7 c6 s4 c1 s3 + d7 c6 s1 s2 c4 +
d5 s4 s1 c2 c3 + d5 s4 c1 s3 + d5 s1 s2 c4 + s1 s2 d3
m31 = −c7 c6 c5 s2 c3 c4 − c7 c6 c5 c2 s4 + c7 c6 s2 s3 s5 + c7 s6 s2 c3 s4 −
c7 s6 c2 c4 + s7 s5 s2 c3 c4 + s7 s5 c2 s4 + s7 s2 s3 c5
m32 = s7 c6 c5 s2 c3 c4 + s7 c6 c5 c2 s4 − s7 c6 s2 s3 s5 − s7 s6 s2 c3 s4 +
s7 s6 c2 c4 + c7 s5 s2 c3 c4 + c7 s5 c2 s4 + c7 s2 s3 c5
m33 = −s6 c5 s2 c3 c4 − s6 c5 c2 s4 + s6 s2 s3 s5 − c6 s2 c3 s4 + c6 c2 c4
m34 = −d7 s6 c5 s2 c3 c4 − d7 s6 c5 c2 s4 + d7 s6 s2 s3 s5 − d7 c6 s2 c3 s4 +
d7 c6 c2 c4 − d5 s2 c3 s4 + d5 c2 c4 + c2 d3 + d1
416 D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot
Geometric Jacobian
" #
Jν
J=
Jω
m11 m12 . . . m16
m21 m22 . . . m26
J = .. .. . . ..
. . . .
m61 m62 . . . m66
where
m11 = −d7 s6 c5 c4 s1 c2 c3 − d7 s6 c5 c4 c1 s3 + d7 s6 c5 s1 s2 s4 +
d7 s6 s5 s1 c2 s3 − d7 s6 s5 c1 c3 − d7 c6 s4 s1 c2 c3 − d7 c6 s4 c1 s3 −
d7 c6 s1 s2 c4 − d5 s4 s1 c2 c3 − d5 s4 c1 s3 − d5 s1 s2 c4 − s1 s2 d3
m12 = −c1 (d7 s6 c5 s2 c3 c4 + d7 s6 c5 c2 s4 − d7 s6 s2 s3 s5 + d7 c6 s2 c3 s4 −
d7 c6 c2 c4 + d5 s2 c3 s4 − d5 c2 c4 − c2 d3 )
m13 = −s1 d7 s6 c5 c3 c4 + s1 d7 s6 s3 s5 − s1 d7 c6 c3 s4 − s1 d5 c3 s4 −
c2 d7 s6 c5 c4 c1 s3 − c2 d7 s6 s5 c1 c3 − c2 d7 c6 s4 c1 s3 − c2 d5 s4 c1 s3
m14 = −c1 c3 d7 s6 c5 c2 s4 − c1 d7 s6 c5 s2 c4 + s1 s3 d7 s6 c5 s4 +
c1 c3 d7 c6 c2 c4 − c1 d7 c6 s2 s4 − s1 s3 d7 c6 c4 + c1 c3 d5 c2 c4 −
c1 d5 s2 s4 − s1 s3 d5 c4
m15 = −s6 (c5 c1 c2 s3 − s5 c1 s2 s4 − s5 c4 s1 s3 + c5 s1 c3 + s5 c4 c1 c2 c3 ) d7
m16 = −d7 (c6 s5 c1 c2 s3 + c6 s5 s1 c3 + s6 s4 c1 c2 c3 − c6 c5 c4 c1 c2 c3 +
c6 c5 c1 s2 s4 + c6 c5 c4 s1 s3 + s6 c1 s2 c4 − s6 s4 s1 s3 )
m17 = 0
D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot 417
m21 = d7 s6 c5 c4 c1 c2 c3 − d7 s6 c5 c4 s1 s3 − d7 s6 c5 c1 s2 s4 − d7 s6 s5 c1 c2 s3 −
d7 s6 s5 s1 c3 + d7 c6 s4 c1 c2 c3 − d7 c6 s4 s1 s3 + d7 c6 c1 s2 c4 +
d5 s4 c1 c2 c3 − d5 s4 s1 s3 + d5 c1 s2 c4 + c1 s2 d3
m22 = −s1 (d7 s6 c5 s2 c3 c4 + d7 s6 c5 c2 s4 − d7 s6 s2 s3 s5 + d7 c6 s2 c3 s4 −
d7 c6 c2 c4 + d5 s2 c3 s4 − d5 c2 c4 − c2 d3 )
m23 = −c2 d7 s6 c5 c4 s1 s3 − c2 d7 s6 s5 s1 c3 − c2 d7 c6 s4 s1 s3 − c2 d5 s4 s1 s3 +
d7 s6 c5 c4 c1 c3 − d7 s6 s5 c1 s3 + d7 c6 s4 c1 c3 + d5 s4 c1 c3
m24 = −s2 d7 c6 s4 s1 + s3 d7 c6 c1 c4 + s1 c3 d7 c6 c2 c4 − s2 d5 s4 s1 + s3 d5 c1 c4 +
s1 c3 d5 c2 c4 − s2 d7 s6 c5 c4 s1 − s3 d7 s6 c5 c1 s4 − s1 c3 d7 s6 c5 c2 s4
m25 = −s6 (−c5 c1 c3 + s5 c4 s1 c2 c3 + c5 s1 c2 s3 − s5 s1 s2 s4 + s5 c4 c1 s3 )d7
m26 = d7 (c6 c5 c4 s1 c2 c3 + c6 c5 c4 c1 s3 − c6 c5 s1 s2 s4 − c6 s5 s1 c2 s3 +
c6 s5 c1 c3 − s6 s4 s1 c2 c3 − s6 s4 c1 s3 − s6 s1 s2 c4 )
m27 = 0
m31 = 0
m32 = −d7 s6 c5 c2 c3 c4 + d7 s6 c5 s2 s4 + d7 s6 c2 s3 s5 − d7 c6 c2 c3 s4 −
d7 c6 s2 c4 − d5 c2 c3 s4 − d5 s2 c4 − s2 d3
m33 = s2 (d7 s6 c5 c4 s3 + d7 s6 s5 c3 + d7 c6 s4 s3 + d5 s4 s3 )
m34 = −c2 d5 s4 − c3 d7 c6 s2 c4 + c3 d7 s6 c5 s2 s4 − c2 d7 c6 s4 −
c2 d7 s6 c5 c4 − c3 d5 s2 c4
m35 = s6 (s5 s2 c3 c4 + s5 c2 s4 + s2 s3 c5 ) d7
m36 = −d7 (−c6 s2 s3 s5 + c6 c5 c2 s4 − s6 s2 c3 s4 + c6 c5 s2 c3 c4 + s6 c2 c4 )
m37 = 0
418 D Geometric Jacobian for the Seven-Joint DLR-Kuka Robot
m41 = 0
m42 = −s1
m43 = c1 s2
m44 = −c1 c2 s3 − s1 c3
m45 = s4 c1 c2 c3 − s4 s1 s3 + c1 s2 c4
m46 = −s5 c4 c1 c2 c3 + s5 c4 s1 s3 + s5 c1 s2 s4 − c5 c1 c2 s3 − c5 s1 c3
m47 = s6 c5 c4 c1 c2 c3 − s6 c5 c4 s1 s3 − s6 c5 c1 s2 s4 − s6 s5 c1 c2 s3 −
s6 s5 s1 c3 + c6 s4 c1 c2 c3 − c6 s4 s1 s3 + c6 c1 s2 c4
m51 = 0
m52 = c1
m53 = s1 s2
m54 = −s1 c2 s3 + c1 c3
m55 = s4 s1 c2 c3 + s4 c1 s3 + s1 s2 c4
m56 = −s5 c4 s1 c2 c3 − s5 c4 c1 s3 + s5 s1 s2 s4 − c5 s1 c2 s3 + c5 c1 c3
m57 = s6 c5 c4 s1 c2 c3 + s6 c5 c4 c1 s3 − s6 c5 s1 s2 s4 − s6 s5 s1 c2 s3 +
s6 s5 c1 c3 + c6 s4 s1 c2 c3 + c6 s4 c1 s3 + c6 s1 s2 c4
m61 = 1
m62 = 0
m63 = c2
m64 = s2 s3
m65 = −s2 c3 s4 + c2 c4
m66 = s5 s2 c3 c4 + s5 c2 s4 + s2 s3 c5
m67 = −s6 c5 s2 c3 c4 − s6 c5 c2 s4 + s6 s2 s3 s5 − c6 s2 c3 s4 + c6 c2 c4
Appendix E
Notation
Kinematics
B, S0 base coordinate system
Si DH-coordinate system i
pi origin of Si (given as vector in base coordinates)
x0 , y0 , z0 coordinate axes of system S0
xi , yi , zi coordinate axes of system Si
ai , αi , θi , di DH-parameters
0M DH-matrix transforming fom S0 to S1
1
i−1 M DH-matrix transforming fom Si−1 to Si
i
li length of a link in a kinematic chain
419
420 E Notation
t joint torque
f force vector
m moment vector
τ magnitude of a joint torque
v linear velocity
ω angular velocity
Jv linear velocity part of a geometric Jacobian
Jω angular velocity part of a geometric Jacobian
Machine Learning
421
422 List of Figures
2.1 DH-table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.2 DH-table for the three-joint robot . . . . . . . . . . . . . . . . . . 53
5.1 grid points and shift values for Figure 5.23 and
Figure 5.24. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
433
434 List of Tables
Index
435
436 INDEX