Systems, devices, and methods for tracking moving targets

Abstract

A system for tracking a moving target having up to six degrees of freedom and rapidly determining positions of the target, said system includes an easy to locate precision optical target fixed to the target. This system includes at least two cameras positioned so as to view the optical camera from different directions with each of the at least two cameras being adapted to record two dimensional images of the precision optical target defining precise target point. A computer processor is programmed to determine the target position of x, y and z and pitch, roll and yaw. In an embodiment, the system can be configured to utilize an iteration procedure whereby an approximate first-order solution is proposed and tested against the identified precise target points to determine residual errors which can be divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction.

Description

BACKGROUND

There are various modalities for performing medical imaging of patients. For example, magnetic resonance imaging (MRI) is a medical imaging technique used in radiology to visualize internal structures of the body in detail. An MRI scanner is a device in which the patient or a portion of the patient's body is positioned within a powerful magnet where a magnetic field is used to align the magnetization of some atomic nuclei (usually hydrogen nuclei—protons) and radio frequency magnetic fields are applied to systematically alter the alignment of this magnetization. This causes the nuclei to produce a rotating magnetic field detectable by the scanner and this information is recorded to construct an image of the scanned region of the body. These scans typically take several minutes (up to about 40 minutes in some scanners) and in prior art devices any significant movement can ruin the images and require the scan to be repeated.

U.S. Pat. No. 8,121,361 issued Feb. 21, 2012, entitled Motion Tracking System for Real Time Adaptive Imaging and Spectroscopy describes a system that adaptively compensates for subject motion. and the disclosure therein is hereby incorporated herein by reference.

SUMMARY

An accurate and reliable method of determining the dynamic position and orientation of a patient's head or other body portion during MRI scanning is a requirement in any attempt to compensate for subject motion during the procedure. Toward this end, disclosed herein are systems and methods that include practical optical head tracking capability using at least a first sensor, e.g., a first camera, and a second sensor, e.g., a second camera, such as a pair of cameras, for example ordinary CCD cameras, ordinary white light or LED illumination, and a marker target, such as a compact, inexpensive target which mounts comfortably and rigidly to the subject's skeletal frame The camera can be configured to detect any desired wavelength or range of wavelengths of energy, including one or more of the infrared, near-infrared, visible, or ultraviolet spectra for example. Some preferred embodiments can track head and other body motion having up to and including six degrees of freedom (sometimes referred to as 6-DOF).

In an embodiment, the system is configured for tracking a moving target having up to six degrees of freedom and rapidly determining positions of the target, said system includes an easy to locate precision optical target fixed to the target. The system can also include at least two cameras positioned so as to view the optical camera from different directions with each of the at least two cameras being adapted to record two dimensional images of the precision optical target defining a precise target point. In an embodiment, a computer processor is programmed to determine the target movement in Cartesian coordinates of x, y and z and pitch, roll and yaw utilizing an algorithm adapted to identify a set of precise target points on the precision optical target and the x, y and z displacement and the pitch, roll and yaw rotation of the precise target points based on optical images collected by the at least two cameras. The system can utilize an iteration procedure whereby an approximate first-order solution is proposed and tested against the identified precise target point projections on the cameras to determine residual errors which are then divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction. The system can be configured to repeat the above actions until residual error becomes smaller than desired accuracy. Using this process the system can be configured to determine the position of the target at rates of at least 100 times per second with translations accuracies of about or no more than about 0.1 mm and angle accuracies of about or no more than about 0.1 degrees. With repetition rates in the range of 100 times per second, the full 6-DOF movement determination can be performed for each repetition. In these embodiments the results of each movement determination is used for the initial first order solution during the next iteration.

The six degrees of freedom movements are over orthogonal directions x, y, and z and roll, pitch and yaw angles. Direction x is along the spinal axis. Direction y perpendicular to x is along the shoulder to shoulder direction and direction z is perpendicular to both x and y and in the floor-to-ceiling direction assuming the patient is lying on his back parallel to the floor. The roll angle is about the x-axis; the angle made by a shaking head “No”. The pitch angle is about the y-axis; the angle made by shaking head “Yes” and the Yaw angle is about the z-axis, the angle made by leaning head toward a shoulder.

In an embodiment, the desired accuracy is about 0.1 mm for each of the directions and about 0.1 degrees for each of the angles. Movements are measured relative to a pivot point in the patient's neck. In an embodiment the pivot point is located at the base of the patient's neck where the head swivels for nod turn and lean motions. The offset of the precision optical target from this pivot point position is Δy=0, Δx −4.5″, Δz=5.5″. The precision of these offsets is not critical since all motions of interest are relative motions. The six measurements are x, y, and z distances and roll, pitch and yaw angles. In some embodiments, the measurements are up-dated at a rate of about 100 solutions per second with a latency of about 10 milliseconds. The system can be configured to report to MRI systems the exact position or the approximate position of the head with accuracies of about or better than about 0.1 mm in distances and about 0.1 degree in angles.

One possible coordinate system for reporting 6-DOF motions to the MRI field compensation system is a Cartesian system aligned with the symmetry axis of the head coil. The head coil coordinate system is coincident with body coordinates in the nominal (“square”) head position. Target displacements and rotations can be reported to the coil field compensation system using this system of coordinates.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a drawing of a precision optical target, according to some embodiments of the invention.

FIG. 1B show the target of FIG. 1A attached to the top front teeth of a patient.

FIG. 1C show the size of one embodiment of the optical target relative to a U.S. penny.

FIGS. 2A and 2B show how two cameras together provide sensitivity needed to track motion, according to some embodiments of the invention.

FIGS. 3A and 3B show how a patient's head and two cameras are located in an MRI device, according to some embodiments of the invention.

FIGS. 4A and 4B show how Cartesian coordinates are used relative to a patient's head for the purpose of tracking motions, according to some embodiments of the invention.

FIG. 5 shows how three points of the precision optical target are imaged on the focal plane of each of the two cameras, according to some embodiments of the invention.

FIG. 6 shows the results on one camera image of a 15 degree yaw movement (about the z-axis), according to some embodiments of the invention.

FIGS. 7A and 7B shows how two cameras are able to monitor precisely a pitch movement (about the y-axis), according to some embodiments of the invention.

FIGS. 8A and 8B show how a roll movement (about the x-axis) is monitored, according to some embodiments of the invention.

FIG. 9 shows how x-axis translation (positive toward the top of the patient's head) is monitored on one camera, according to some embodiments of the invention.

FIGS. 10A and 10B shows the effect of y-axis translation (positive to the patient's right side) as monitored on the two cameras, according to some embodiments of the invention.

FIGS. 11A and 11B show the effect of z-axis translation (toward the ceiling), according to some embodiments of the invention.

FIGS. 12A and 12B show the effect of simultaneous pitch and x-axis and z-axis translation, according to some embodiments of the invention.

FIGS. 13A and 13B show the effect of simultaneous roll and y-axis and z-axis translation, according to some embodiments of the invention.

FIGS. 14B and 14C display features of an iteration technique utilized to precisely monitor head movement utilizing the camera images of the precision optical target, according to some embodiments of the invention.

FIG. 14D is a flow diagram of an iterative process for tracking movement.

FIG. 15 is a block diagram depicting an embodiment of a computer system configured to implement one or more embodiments of the methods, devices, and systems described herein.

FIGS. 16A and 16B show techniques for camera calibration, according to some embodiments of the invention.

DETAILED DESCRIPTION

Optical Target Fixed to an Anatomical Location, e.g., the Head

To overcome the challenge of the elastic nature of skin, in an embodiment, an optical tracking target can be coupled to the upper teeth of the patient. One accessible feature which is rigid to the skull is the upper teeth. Unlike the teeth on the lower jawbone, the upper teeth are rigidly affixed to the skull of the patient. In an embodiment, a compact and reliable optical tracking target can be attached to one or more of the upper teeth with a clip-on or other coupling device. Such attachment devices can be configured to be extremely comfortable. In an embodiment, a printed precision optical target is attached to the top front teeth of a patient.

An optical target can be configured to be easy to locate with a high degree of accuracy regardless of orientation in a sensor field of view. A circle or series of concentric circles or ellipses can be potentially advantageous in this regard. Furthermore, to accommodate the fastest composite 2D data processing methods, a number (at least 3) of centroid positions can be discernible at every instant in time. The target can be, in some embodiments, composed of three sets of concentric circles or ellipses located at the vertices of an equilateral triangle. Compactness is desired for practical reasons, but the minimum size and spacing of the targets is dictated to large extent by characteristics of the sensors and the available non-occluded optical lines of sight through the MRI field compensation coil. A tradeoff arises, for instance, between the minimum size of the target and the cost of the imaging cameras used to sense the head motion—the smaller the edge dimension of the target triangle, the more pixels required in the camera sensor, and the faster the readout and processing electronics required.

As a reasonable compromise, in some embodiments an equilateral triangle side length of 0.5 inches can be adopted. The printed target pattern includes a solid central elliptical dot of 1/16″ minor diameter at each triangle vertex, and each dot is surrounded by a first concentric ellipse of 3/16″ minor diameter and 1/32″ line width, and a second concentric ellipse of 5/16″ minor diameter and 1/32″ line width (ellipses scaled to look circular from camera nominal 45° look angle). In this embodiment, the entire target measures about 1 inch wide by about 0.93 inches high. Other dimensions are possible.

A camera viewing this target is able to determine the centroid of each ellipse on the target pattern using a simple brightness moment calculation, independent of orientation of the target. The target itself subtends only a small portion of the camera field of view, but is recognizable by its high contrast and lack of gray scale. In embodiments the computer processor is programmed to track each of the three sub-targets by enclosing each of the three sub-targets within a sub-pixel array of 48×48 pixels and to calculate centroids of each sub-target by dividing (a) the sum of the product of pixel darkness and pixel position by (b) the sum of the pixel darkness of all of the pixels in the 48×48 sub-pixel array. The processor is also programmed to move each of the 48×48 pixel arrays so that its target is always located fully within the sub-pixel array. With sufficient camera spatial and brightness resolution and target illumination and contrast, centroid positional accuracy of about 0.1 pixels in row and/or column or less is achievable using this target.

FIG. 1A is an enlarged view of the optical target and two of the three Cartesian axes. FIG. 1B shows a full-scale target (compared to a U.S. penny) affixed to NTI coupling device for placement on the upper teeth of the patient. (Lower right) Subject with optical target and night guard clipped onto front teeth.

Latency

Latency in the measurement of head motion using optical tracking techniques is comprised of the camera sensor integration and readout time, the target centroid determination time and the 6-DOF decomposition time. In order to reliably track head motions as fast as 2 cm/second and head rotations as fast as 10 degrees per second, a camera frame rate of about 100 Hz is desired, with electronic shuttering to freeze motion at rates up to 10 times this speed for sharp resolution of the optical target without blurring. A significant field of view is required to accommodate large motions, so fast camera readout without expensive mechanical tracking capabilities will require either a low pixel density or a camera with a larger focal plane but the ability to window a smaller region of interest for readout. Centroid and 6-DOF decomposition algorithms running in composite 2D, rather than full 3D space, and utilizing rapidly converging solution methods can be capable of returning solutions to the compensating head coil electronics at 100 solutions per second, with about 10 ms of latency. In some embodiments, the system can be configured to operate with a latency that enables it to update the scanner in between each image acquisition

Cameras

For a subject wearing or coupled with the optical head tracking target, the target size and subject rotation angles and translation position determine the physical location of the three target centroids precisely in three dimensions. With precise knowledge of these angles and the optical sensor (camera and lens) parameters—pixel pitch, lens focal length and radial distortion, camera location and orientation relative to nominal target position—the location of the target centroid projections on the focal plane sensor can be predicted to any level of accuracy even prior to measurement.

In principle, the inverse problem should be equally simple as long as the 3D position of the target centroids can be ascertained optically. Using two cameras, a stereo view of the centroid projections can be used to determine the 3D location in space of each target centroid, and the 6-DOF displacement vector can then be determined through a simple matrix inversion. In practice, however, this approach leads to expensive and complicated requirements on camera pixel density, pixel count, camera alignment and camera calibration.

An alternate unfolding approach dispenses with stereo ranging but uses separate 2D projections from two cameras without attempting to correlate absolute target positions on the two cameras. This approach eliminates the strict requirements on camera alignment and magnification matching characteristic of the stereo vision approach, and also relaxes the pixel density and count requirements needed to obtain the required positional accuracy (about 0.1 mm in translation and about 0.1 degrees in rotation) by about a factor of 20, resulting in significant savings in cost and processing speed.

Even for this 2D measurement approach some basic steps can be taken to calibrate camera parameters once the cameras are integrated with the head coil; these can be performed at the manufacturing facility. These include measuring the projected pixel location of a single reference point on both cameras, as well as the camera magnification factors for pixel displacement per degree of rotation in pitch, yaw and roll, and per mm of translation along x, y and z. However, as stated before, it is not necessary that the cameras be exactly aligned in space (e.g. perfectly normal) or that their magnifications (lens focal length and distance to reference point) be identical, as is easily verified by simulation.

Stereo Versus Composite 2D Vision Requirements

With a single camera viewing the target from 45 degrees off of vertical in the target plane, the camera sees very little centroid displacement when the target moves in the direction of the camera (e.g. upward vertical translation equal to horizontal translation in the camera direction, with no rotation). Assuming a 7 micron pixel pitch, a 25 mm lens, and a working distance of 14 inches, target displacement in the camera direction may be at least 0.6 mm before the target can be detected as a 0.1-pixel increase in target centroid separation. However, as shown in FIGS. 2A and 2B a second camera placed orthogonally, e.g. at −45 degrees relative to vertical in the same plane, is maximally sensitive to this same motion, seeing a full pixel displacement of each centroid for a diagonal translation of only 0.1 mm. The second camera eliminates the “blind spot” that a single camera has to motion along its optical axis. While certain embodiments described systems in which cameras are positioned orthogonally, cameras can also be placed at relative angles other than orthogonal with respect to vertical in the same plane, depending on the desired clinical result.

Camera Depth of Field

To accommodate head roll of +/−15 degrees plus the 0.85-inch target width at a working distance of 14 inches, the lens can be configured to provide sharp focus for distances between 13″ and 15.5″. At f/22, assuming a circle of confusion slightly smaller than a camera pixel (7 microns), a 25 mm focal-length lens provides this necessary depth of field a nominal 14-inch focus. At this working distance, the optical path can be folded with a turning mirror (FIG. 3) or otherwise configured to fit within the 70 cm diameter bore of the main MRI coil. A non-ferrous camera can be utilized in the MRI environment. In an embodiment, it can be cost effective to repackage a commercial camera for use in the strong magnetic field.

In some embodiments, one possible camera that can be utilized or modified for use with systems and methods as disclosed herein, is produced by Allied Vision Technologies and designated the Prosilica GE-680 Monochrome CCD Camera. This camera features a Kodak KAI-0340 ⅓″ 640×480 VGA focal plane sensor with 7.4 μm square pixels and a fast Gigabit Ethernet output delivering up to 205 frames per second at 12-bit pixel depth. An inexpensive possible lens for use is an Edmund Optics TechSpec 25 mm high-resolution fixed focal length lens.

For this camera and lens, at 14 inches from the target at 45° incidence, the 5/16″ diameter target circles project to ellipses on the camera, with the minor diameter of the largest ellipses at about 28 pixels and the major diameter at about 40 pixels. With sufficient S/N ratio (target illumination) and lens MTF (sharpness), this pattern should allow accurate centroiding to about 0.1 pixels in row and/or column or less. The entire projected target subtends about 128 H×168 V pixels, and allowing for head roll of +/−11.5 degrees, a camera with 640 horizontal pixels (pixel columns) can accommodate the entire field of interest without mechanical tracking provisions.

FIGS. 3A and 3B show a modified head coil with cameras mounted longitudinally and turning mirrors to accommodate a longer working distance than is possible with a straight optical path in the constrained space of the main MRI coil bore. In embodiment, the system is configured with two or more cameras with a direct view of the optical tracking targets without the use of mirrors.

Six Degree-of-Freedom Measurement and Reporting Algorithm

In some embodiments, the MRI Head Tracker takes real-time input from two 2D imaging sensors and analyzes these data to determine and report motions in six degrees of freedom with minimal latency. This task can be performed by detecting and measuring the three centroid positions on the target and utilizing those positions with a reporting algorithm to determine the position of the patient's head.

Six-Degree-of-Freedom Coordinate System

In an embodiment, the system is configured to use a coordinate system for reporting 6-DOF motions to the MRI field compensation system that is a Cartesian system aligned with the symmetry axis of the head coil as shown in FIGS. 4A and 4B. Head coil coordinate system shown in FIG. 4A is coincident with body coordinates in the nominal (“square”) head position as shown in FIG. 4B. The z direction is into and out of the plane of the drawing. Target displacements and rotations are reported to the coil field compensation system using this system of coordinates.

Coordinate definitions are adopted by the same conventions used in defining aircraft motion, except that the rotation directions are taken to be right-handed (positive for counter-clockwise rotation about the basis direction vectors):

x is the longitudinal (chin-to-crown) direction, with values increasing toward the top of the head

y is the transverse (left-to-right) direction, with increasing values toward the patient's right ear

z is the up-down direction, with increasing values toward the ceiling

$\psi ={\tan }^{-1}\left(\frac{\Delta y}{\Delta x}\right)$ is the yaw angle or right-handed rotation about the z-axis (head lean toward shoulder while facing forward, zero at normal “square” position, positive values for patient leaning toward patient's right shoulder)

$\theta ={\tan }^{-1}\left(\frac{\Delta x}{\Delta z}\right)$ is the pitch angle or right-handed rotation about the y-axis (nodding “yes,” zero at normal “square” position, positive values for patient looking “upward”)

$\phi ={\tan }^{-1}\left(\frac{-\Delta y}{\Delta z}\right)$ is the roll angle or right-handed rotation about the x-axis (shaking the head “no,” zero at normal “square” position, positive values for patient looking toward patient's left side).

The origin of coordinates and angle zero references are arbitrary, as only relative motions are reported, however two convenient reference origin positions exist: 1) at the center of the target in its normal (“square”) head position, and 2) at the base of the neck at the point where the head swivels for nod, turn and lean motions. The latter is adopted here (as shown in FIG. 2), simply for ease in orthogonalizing the set of principal observation parameters with common motion directions in the 6-DOF decomposition algorithm.

Target Displacement Equations

The full 6-DOF translation is composed of a 3-D displacement as well as a 3-axis rotation. To first order we assume that the skull moves as a rigid body about a single rotation point somewhere in the neck. From this point the translation becomes separable from the rotation, so this is chosen as the coordinate origin. The rotations are separated into roll, pitch and yaw as described above, and the translated position through rotation follows the Euler rotation matrix formulation as follows (using right-handed angle conventions). The x, y, and z displacement coordinates then follow the independent translations:

$\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \phi & -\sin \phi \\ 0& \sin \phi & \cos \phi \end{array}\right)\left(\begin{array}{ccc}\cos \theta & 0& \sin \theta \\ 0& 1& 0\\ -\sin \theta & 0& \cos \theta \end{array}\right)\left(\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right).$

Decomposing the six independent translations from the absolute and relative displacements of the measured target centroids is the subject of this effort. The 2D inverse problem is somewhat more difficult than the 3D problem, in that after the target centroid projections in focal plane row and column are determined, significant degeneracies remain in the unfolding matrices for each camera. Combining the data from both cameras removes these degeneracies through a series of interrelated, nonlinear equations. The fastest procedure for solving this inverse problem is obtained by the Newton-Raphson method or a variant thereof, whereby an approximate first-order solution is proposed and tested against the known (measured) centroid locations on the two camera focal planes. The residual error is divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction. The first-order solution is chosen by considering the features of the projected target pattern which are most strongly affected by a single rotation angle or displacement, and linearizing the inversion problem along these feature axes.

A 6-DOF motion simulation and decomposition algorithm was developed and tested to allow simulation of arbitrary motions and then verify the ability of a pair of orthogonal cameras to decompose centroid measurements at the 0.1-pixel level into distinct x, y, z, roll, pitch and yaw components at the requisite level of accuracy.

Six-Degree-of-Freedom Motion Determination Algorithm

General subject motion is a superposition of translation along x, y, and z as well as rotation about the x, y and z axes (designated roll, pitch and yaw respectively). Displacements along each of these degrees of freedom are not sensitive to coordinate system origin; however it is convenient (as explained above) for modeling purposes to place an origin near the region of the spine about which the head rotates and swivels, and a secondary reference point at the center of the optical tracking target in the nominal (“correct”) head position and orientation. This secondary reference is typically offset from the spinal origin by ˜10 cm in x and ˜10 cm in z.

The target shown in FIG. 1, as viewed from a single camera, appears as three sets of concentric ellipses with centroids projected onto three different positions (column, row) on the camera focal plane. The camera is centered along the (x=constant) plane of the target and aligned such that the central pixel row images this plane, at an angle of 45 degrees with respect to both the y and z axes and with the nominal target center projected to the central pixel column. Using a camera with 7.4 micron pixels and a 25 mm lens, positioned at a distance of 14.1 inches from the nominal target center, centroids from the vertices of an equilateral triangle target with sides of length 0.5 inches are projected onto the camera focal plane as shown in FIG. 5. This figure shows projected positions of target centroids for a target with sets of concentric circles arranged about the vertices of an equilateral triangle of side length 0.5 inches, using a camera focal length 25 mm, pixel pitch 7.4 microns and view angle 45 degrees (camera to right and above paper), aligned with the camera centerline. The 45-degree view angle results in the foreshortening of the equilateral triangle from left to right across the focal plane.

Yaw

Rotation about the z-axis is designated as yaw; a positive or “right handed” rotation about this axis (head leaning to subject's right shoulder) results in a counterclockwise rotation of the target. Because this rotation usually occurs about a point lower in the neck, it is typically accompanied by a translation to the subject's right side (camera left), as seen in FIG. 6. Projected positions of target centroids for same conditions as FIG. 5, but before and after inducing a yaw of 15 degrees.

The median of the centered target triangle (as shown at the right in FIG. 6) is aligned approximately with a single column of camera pixels, but is rotated out of this alignment (as shown at the left side of FIG. 6) by yaw. For the camera, lens, target distance and target size described above, a yaw of only 0.1 degrees results in a relative displacement of 0.13 pixel columns between the two ends of the median. Assuming that the centroid algorithm is able to determine position of the triangle vertices to 0.1 pixels in row and column, the yaw angle determination is measurable down to and accurate to about 0.1 degrees.

Pitch

Rotation about the y-axis is designated as pitch; a positive or “right-handed” rotation about this axis (head tipped back) results in motion of the target upward off the gantry (+z) and toward the top of the head (+x). For a single camera this projection is not easily distinguishable from a simultaneous target displacement in x and y (see FIG. 7), but for two cameras at opposite sides of the head the apparent y-displacement is in the opposite direction, removing this degeneracy. A second degeneracy with pitch rotation remains, for simultaneous target translation in +x and +z—this is discussed in more detail later—but the tilt of the target plane during a pitch rotation yields a small difference in the amount of translation of the base of the target triangle relative to its apex, thus resulting in a slight apparent rotation of the target triangle as shown in FIG. 7, which is not a characteristic of simple translation. This becomes in some embodiments the defining characteristic of pitch motion.

FIGS. 7A and 7B show the projected positions of target centroids for same conditions as FIG. 5, but before and after a target pitch of 8 degrees. Left is view from a camera at the left side and above the paper, right is view from a camera at the right side and above the paper. In each case motion away from the gantry (+z) makes the target appear more distant from the observer.

Roll

Rotation about the x-axis is designated as roll; a positive or “right-handed” rotation about this axis (head pointing toward subject's left side) results in a motion of the target toward the subject's left (−y). For a single camera this motion is not easily distinguishable from a displacement in y (see FIG. 8), but for two cameras the difference in position and in apparent foreshortening of the triangle is much more pronounced for rotation than for translation. This is because the roll moves the target plane closer to normal incidence with one camera sightline and further from normal incidence with the other camera sightline, at a rate which is much larger than that for a simple translation (or yaw). There is a significant degeneracy between roll and simultaneous +y and +z translation which is only resolved comparing the lengths of the triangle base as seen between the two cameras. A large difference in the base lengths is a characteristic of roll motions and not a characteristic of y+z translation, hence this is the distinguishing characteristic for roll.

As shown in FIGS. 8A and 8B the projected positions of target centroids for same conditions as for FIG. 5, but before and after target roll of 12 degrees. Left is view from a camera at the left side and above the paper, right is view from a camera at the right side and above the paper. The camera at the left side sees much wider triangle because target plane is closer to normal to this camera sightline. The camera at the left also sees much larger displacement of triangle center.

X-Axis Translation

Translation along the x-axis (positive toward top of head) results in a motion of the target along the vertical direction of the camera focal plane (see FIG. 9). Unlike for pitch rotation (which also involves a translation in z), the target does not move significantly between pixel columns, and rotation of the target triangle is minimal. This up-down camera translation without accompanying rotation is the distinguishing characteristic of x-axis translation. FIG. 9 shows the projected positions of target centroids for same conditions as for FIG. 5, but before and after target translation of 12 mm in x.

Y-Axis Translation

Translation along the y-axis (positive toward subject's right side) results in a motion of the target along the horizontal axis of the camera focal plane (see FIGS. 10A and 10B). Unlike for roll (which also involves a differential rotation of the target plane as seen by the left and right side cameras), the target's projected size, displacement and rotation varies only slightly between left and right camera views for y-axis translation; this is the distinguishing characteristic for y-displacement. FIGS. 10A and 10B show projected positions of target centroids for same conditions as FIG. 5, but before and after target translation of 15 mm along y-axis. Left is view from a camera at the left side and above the paper, right is view from a camera at the right side and above the paper. Unlike roll, target displacements and sizes are similar for two cameras viewing y-axis translation.

Z-Axis Translation

Translation along the z-axis (positive toward the ceiling) results in apparent motion of the target along the horizontal axis of the camera focal plane. Unlike for y translation, however, the direction of the horizontal displacement is opposite between the left-side and right-side cameras (see FIGS. 11A and 11B). This is the distinguishing characteristic for z-axis translation. FIGS. 11A and 11B show projected positions of target centroids for same conditions as for FIG. 5, but before and after target translation of 15 mm along the z-axis. Left is view from a camera at the left side and above the paper, right is view from a camera at the right side and above the paper. Unlike translation along y, apparent target displacement is in opposite direction in two camera views.

Non-Degenerate Target Motion Parameters

Pitch Versus (X+Z) Translation Degeneracy

Pitch is nearly degenerate with simultaneous x and z translation, except for a small tilt in the triangle vertical which results from the tilt of the target plane about the y axis. This tilt creates an apparent clockwise rotation of the triangle from the left-side view and an apparent counterclockwise rotation from the right side view, as shown in FIGS. 12A and 12B. These drawings show projected positions of target centroids for same conditions as FIG. 5, but before and after target pitch of 4 degrees and translations in x and z of −9.5 mm and +7.8 mm respectively. FIG. 12A is view from a camera at the left side and above the paper, FIG. 12B is view from a camera at the right side and above the paper. The camera at left sees triangle rotated clockwise, with upper vertices rotated away from the camera because of an increase in z relative to the lower vertex. The camera at the right sees triangle rotated counterclockwise for the same reason. For a pitch motion of 0.1 degrees accompanied by translations in x and z of −0.244 mm and 0.187 mm respectively, the triangle apex centroid does not move in either camera view. However, in this case, the left-side camera sees the triangle base displaced by 0.13 pixels to the right while the right-side camera sees the triangle base displaced by 0.13 pixels to the left. Assuming the centroiding routine can locate the vertices of the target triangle to an accuracy of 0.1 pixels, a pitch as small as 0.1 degrees is distinguishable from a simple translation by comparison of the vertical tilts.

Roll Versus (Y+Z) Translation Degeneracy

Roll is nearly degenerate with simultaneous y and z translation, except for larger camera-to-camera differences in apparent target size encountered with roll, resulting from tilt of the target's plane about the x-axis. A significant difference in the apparent length of the target triangle base is a reliable distinguishing characteristic of roll motion rather than simple translation. FIGS. 13A and 13B show projected positions of target centroids for the same conditions as in FIG. 5, but before and after target roll of 4 degrees and translations in y and z of 9.75 mm and 0.34 mm respectively. FIG. 13A is view from a camera at the left side and above the paper, FIG. 13B is view from a camera at the right side and above the paper. Camera at left sees triangle base shrink due to rotation about the x-axis away from camera normal, while camera at right sees triangle base grow due to rotation toward camera normal.

For a roll of 0.1 degrees and translations in y and z of −0.244 mm and 0.0002 mm respectively, the lower centroid is unchanged in both camera views. In this case, the left-side camera sees the target triangle base shrink by 0.15 pixels while the right-side camera sees the triangle base grow by 0.15 pixels. Assuming the centroiding routine can locate the target centroids to an accuracy of 0.1 pixels, shifts of 0.14 pixels should be discernible, so a pitch as small as 0.1 degrees is distinguishable from a simple translation by comparison of the length of the target triangle base.

Six-Degree-of-Freedom Motion Determination Algorithm Architecture

Complementary Projections Versus Stereo Imaging

The target size, rotation angles and translation vector determine the relative displacement of the three target centroids precisely in three dimensions. Precise knowledge of camera and lens parameters (e.g., pixel pitch, lens focal length and radial distortion, camera location and orientation relative to nominal target position), are then sufficient to predict the location of the target centroid projections to better than 0.1 pixels in row and column for each camera. In principle, the inverse problem should be equally simple; the stereo view of the centroid projections determine the 3D location in space of each target centroid, and the 6-DOF displacement vector can then be determined through a simple matrix inversion. In practice, however, this approach leads to expensive and complicated requirements on camera pixel density, pixel count, camera alignment and camera calibration. An alternate unfolding approach dispenses with stereo ranging and uses the two camera projections separately without strict requirements on precise matching of camera alignment and magnification, to determine the 6-DOF displacement vector to within 0.1 degrees in each rotation angle and 0.1 mm along each translation axis. This approach relaxes the pixel density and count requirements by about a factor of 20 relative to the stereo approach, resulting in significant savings in cost and processing speed.

Even for this 2D approach some basic measurements can be made to calibrate camera parameters once the cameras are integrated with the head coil; these can be easily performed at the manufacturing facility. These measurements include the projected pixel location of a single reference point on both cameras, as well as the camera magnification factors for pixel displacement per degree of rotation in pitch, yaw and roll, and per mm of translation along x, y and z. However, as stated before, it is not necessary that the cameras be exactly aligned in space (e.g. perfectly normal) or that their magnifications (lens focal length and distance to reference point) be identical, as has been easily verified by simulation.

Inversion Equations

The 2D inversion problem is somewhat more difficult than the 3D problem, in that after the target centroid projections in focal plane row and column are determined, significant degeneracies remain in the unfolding matrices for each camera. Combining the data from both cameras removes these degeneracies through a series of interrelated, nonlinear equations. The fastest procedure for solving this inverse problem is obtained by a variant of the Newton-Raphson method, whereby an approximate first-order solution is proposed and tested against the known (measured) centroid locations on the two camera focal planes. The residual error is divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction. The first-order solution is chosen by considering the features of the projected target pattern which are most strongly affected by a single rotation angle or displacement, and linearizing the inversion problem along these feature axes.

6-DOF Extraction Algorithm

The method for extracting the 6 degree of freedom displacement matrix from the observed target location on two focal plane cameras is described.

Step 1: Characterizing the Target Images

The optical target consists of elliptical targets shown in FIG. 1A that are drawn so as to appear as circular patterns when imaged at 45 degrees by the two cameras shown in FIG. 2A. The center of each of the three circular patterns define one of the three vertices of an equilateral triangle at the focal plane of each of the two cameras. A centroid calculation routine determines the positions of the centroids at each of the three vertices, on each of two independent cameras. These centroids are displayed on a computer monitor displaying the 640×480 pixels of each of the two cameras. FIG. 5 shows the three vertices being displayed on one of the cameras. These vertex positions are designated (X_i,j, Y_i,j), for vertex index i from 1 to 3, and camera index j from 1 to 2, resulting in twelve measured coordinates. From the twelve measured coordinates, and initialized values of these coordinates, six principal quantities are computed to characterize the two camera views of the equilateral triangle target:

a) Σ_HD—the sum of the horizontal displacements (in pixels) of the target center on camera 1 and camera 2; the formula used is

$\sum _{i=1}^3\sum _{j=1}^2\left(X_{i,j}-X_{0i,j}\right),$ where X_0i,j is the initial (zero displacement) horizontal camera coordinate of each centroid projection.

b) Δ_HD—the difference between the horizontal displacements (in pixels) of the target center for camera 1 and camera 2; the formula used is

$\sum _{i=1}^3\left(X_{i,1}-X_{0i,1}\right)-\left(X_{i,2}-X_{0i,2}\right).$

c) Σ_VD—the sum of the vertical displacements (in pixels) of the target center for camera 1 and camera 2; the formula used is

$\sum _{i=1}^3\sum _{j=1}^2\left(Y_{i,j}-Y_{0i,j}\right),$ where Y_0i,j is the initial (zero displacement) vertical camera coordinate of each centroid projection.

d) Δ_BL—the difference in the apparent base length of the target triangle (in pixels) for camera 1 and camera 2; the formula used is {√{square root over (X_3,1−X_1,1)²+(Y_3,1−Y_1,1)²)}{square root over (X_3,1−X_1,1)²+(Y_3,1−Y_1,1)²)}−√{square root over ((X_3,2−X_1,2)²+(Y_3,2−Y_1,2)²)}{square root over ((X_3,2−X_1,2)²+(Y_3,2−Y_1,2)²)}}.

e) Σ_MT—the sum of the apparent median tilt of the target triangle (offset in horizontal pixels between center-of-base and apex) for camera 1 and camera 2; the formula used is

$\sum _{j=1}^2\left\{\left(X_{2,j}-\frac{X_{3,j}+X_{1,j}}{2}\right)-\left(X_{02,j}-\frac{X_{03,j}+X_{01,j}}{2}\right)\right\}.$

f) Δ_MT—the difference between the apparent median tilt of the target triangle (in pixels) for camera 1 and camera 2; the formula used is

$\left\{\left(X_{2,1}-\frac{X_{3,1}+X_{1,1}}{2}\right)-\left(X_{02,1}-\frac{X_{03,1}+X_{01,1}}{2}\right)\right\}-\left\{\left(X_{2,2}-\frac{X_{3,2}+X_{1,2}}{2}\right)-\left(X_{02,2}-\frac{X_{03,2}+X_{01,2}}{2}\right)\right\}$ Step 2: Characterizing Global Variation in Principal Quantities with 6-DOF Motions

Partial derivatives relative to subject displacements and rotations (φ, θ, ψ, Δx, Δy, Δz), of the principal quantities described above, about the initial (non-displaced) position, are computed numerically. Here:

Roll φ is right-handed rotation about the x-axis

Pitch θ is right-handed rotation about the y-axis

Yaw ψ s right-handed rotation about the z-axis

Δx is toe-to-head direction

Δy is left-to-right direction

Δz is down-to-up direction

Starting from an initial target position in 3-D world space, defined as (φ, θ, ψ, Δx, Δy, Δz)=(0, 0, 0, 0, 0, 0), the initial target vertex world coordinates (x_0i, y_0i, z_0i) are determined for vertex index i=1 to 3, based on the geometric size and shape of the target triangle and definition of a convenient coordinate origin.

Local partial derivatives of each of the principal quantities, with respect to each of the 6 degrees of freedom (roll, pitch, yaw, dx, dy, dz), are performed numerically by evaluating changes in these quantities for small increments in each degree of freedom. Changes in the target vertex positions for specified motions along the six degrees of freedom are computed using the Euler rotation matrices and translation vector:

$\begin{array}{cc}\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \phi & -\sin \phi \\ 0& \sin \phi & \cos \phi \end{array}\right)\left(\begin{array}{ccc}\cos \theta & 0& \sin \theta \\ 0& 1& 0\\ -\sin \theta & 0& \cos \theta \end{array}\right)\left(\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{x}_{0i}\\ y_{0i}\\ z_{0i}\end{array}\right)+\left(\begin{array}{c}\Delta x\\ \Delta y\\ \Delta z\end{array}\right)& \left[1\right]\end{array}$

Subsequently, the camera projections of these new target vertex positions are determined using a geometric projection calculation. Given accurate knowledge of camera positions, pixel size and lens focal length, the horizontal and vertical pixel numbers on each camera focal plane (camera index j equal to 1 or 2) that these new 3-D positions in space should project onto is as follows:

$\begin{array}{cc}{X}_{i,j}=X_{0,j}+\left(\frac{f.l.}{s_{\mathrm{pix}}}\right)\sin {\alpha }_{i,j}\cos {\beta }_{i,j},Y_{i,j}=Y_{0,j}+\left(\frac{f.l.}{s_{\mathrm{pix}}}\right)\sin {\alpha }_{i,j}\sin {\beta }_{i,j}& \left[2\right]\end{array}$

Here X_i,j and Y_i,j are the horizontal and vertical pixel numbers for translated target vertex i projected onto the camera j sensor, X_0,j and Y_0,j are the horizontal and vertical number of the pixel column and row intersected by the optical axis of that camera (typically at or very near the camera center), f.l. and s_pix are the lens focal length and camera pixel pitch, and the angles α_i,j and β_i,j are the polar and azimuth angles locating target vertex i, relative to the camera j focal axis. These angles are calculated from the vertex world coordinates as follows:

$\begin{array}{cc}{\alpha }_{i,j}={\sin }^{-1}\left(\frac{\sqrt{{\left(x_{\perp ,j}-x_i\right)}^2+{\left(y_{\perp ,j}-y_i\right)}^2+{\left(z_{\perp ,j}-z_i\right)}^2}}{\sqrt{{\left(x_{\mathrm{cj}}-x_i\right)}^2+{\left(y_{\mathrm{cj}}-y_i\right)}^2+{\left(z_{\mathrm{cj}}-z_i\right)}^2}}\right),& \left[3\right]\\ {\beta }_{i,j}=\pm {\cos }^{-1}\left(\frac{\left(x_{\perp ,j}-x_i\right)\left(y_{\mathrm{cj}}-y_0\right)-\left(y_{\perp ,j}-y_i\right)\left(x_{\mathrm{cj}}-x_0\right)}{\begin{array}{c}\sqrt{{\left(x_{\perp ,j}-x_i\right)}^2+{\left(y_{\perp ,j}-y_i\right)}^2+{\left(z_{\perp ,j}-z_i\right)}^2}\\ \sqrt{{\left(x_{\mathrm{cj}}-x_0\right)}^2+{\left(y_{\mathrm{cj}}-y_0\right)}^2}\end{array}}\right),& \left[4\right]\end{array}$ where the point (x_⊥,j, y_⊥,j, z_⊥,j) is the point of intersection between the camera optical axis and the plane perpendicular to the optical axis which includes the translated target vertex (x_i, y_i, z_i):x_⊥,j=x₀+κ(x_cj−x₀);y_⊥,j=y₀+κ(y_cj−y₀);z_⊥,j=z₀+κ(z_cj−z₀),  [5]with (x_cj, y_cj, z_cj) defining the 3-D position of camera j, (x₀, y₀, z₀) defining the nominal boresight position of both cameras at the un-displaced target center and the constant κ based on geometric projection and given by:

$\begin{array}{cc}\kappa =\left\{\frac{\left(x_{\mathrm{cj}}-x_0\right)\left(x_i-x_0\right)+\left(y_{\mathrm{cj}}-y_0\right)\left(y_i-y_0\right)+\left(z_{\mathrm{cj}}-z_0\right)\left(z_i-z_0\right)}{{\left(x_{\mathrm{cj}}-x_0\right)}^2+{\left(y_{\mathrm{cj}}-y_0\right)}^2+{\left(z_{\mathrm{cj}}-z_0\right)}^2}\right\}.& \left[6\right]\end{array}$

In equation [4], the inverse cosine function is taken to range from 0 to π, and the appropriate sign for β_i,j is given by:

$\mathrm{sign}\left[{\beta }_{i,j}\right]=\mathrm{sign}\left[\left(z_{\mathrm{cj}}-z_i\right)\left\{\left(x_{\mathrm{cj}}-x_0\right)\left(x_{\perp ,j}-x_i\right)+\left(y_{\mathrm{cj}}-y_0\right)\left(y_{\perp ,j}-y_i\right)\right\}-\left(z_{\perp ,j}-z_i\right)\left\{{\left(x_{\mathrm{cj}}-x_0\right)}^2+{\left(y_{\mathrm{cj}}-y_0\right)}^2\right\}\right]$

During this determination of the camera projection of the 3-D target vertices, a compensation function may be applied for large values of the polar angle α_i,j to account for barrel distortion in the lens, based on prior lens calibration measurements. The geometric value for α_i,j is first computed based on equation [3] and then adjusted for lens distortion by way of a pre-determined look-up table or measured fit function, and this new compensated value for is then used in the calculation of X_i,j and Y_i,j through equation [2].

To numerically evaluate the partial derivatives of the principal quantities about the initialized target position, the un-displaced 3-D target vertex coordinates (x_0i, y_0i, z_0i) are first projected to camera coordinates using equations [2] through [6] above, and initial values are computed for each of the principal quantities described in Step 1 (most should be zero or near-zero at the starting position). Then small increments of roll, pitch, yaw, x-, y- and z-axis displacements are introduced one at a time; for each increment the new world coordinates and the new camera projections of the target vertices are computed and the principal quantities are re-calculated. The change in each principal quantity is divided by the small angular or displacement increment to determine the partial derivative.

For instance, to determine the partial derivatives with respect to roll, the displacement vector (φ, θ, Δx, ix, Δy, Δz)=(δφ, 0, 0, 0, 0, 0) is introduced to the general displacement equation [1] to determine the translated target vertex positions (x_i, y_i, z_i). The conversion to camera coordinates (X_i,j, Y_i,j) is then performed using equations [2] through [6], and the principal quantities are calculated as outlined in Step 1. The difference between each principal quantity and the corresponding value of that quantity for the un-displaced calculation is divided by the small increment in roll, to give the partial derivative of each quantity with respect to roll. To determine partial derivatives with respect to pitch, the displacement vector (φ, θ, Δx, Δy, Δz) (0, δθ, 0, 0, 0, 0) is used to initiate the calculations, and so on for all six degrees of freedom.

Each of these six repetitions produces one column of the global partial derivative matrix:

$\left(\begin{array}{cccccc}\frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial z}\end{array}\right){|}_{\left(0,0,0,0,0,0\right)}.$ Step 3: Determining First-Order Displacement Vector

A first-order approximation to the displacement matrix is determined by multiplying the matrix of measured principal quantities, as determined in Step 1, by the inverse of the partial derivative matrix computed in Step 2:

$\left(\begin{array}{c}{\phi }_0\\ {\theta }_0\\ {\psi }_0\\ {\left(\Delta x\right)}_0\\ {\left(\Delta y\right)}_0\\ {\left(\Delta z\right)}_0\end{array}\right)={\left(\begin{array}{cccccc}\frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial x}& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial y}& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial z}\end{array}\right)}^{-1}\left(\begin{array}{c}{\Sigma }_{\mathrm{HD}}\\ {\Delta }_{\mathrm{HD}}\\ {\Sigma }_{\mathrm{VD}}\\ {\Delta }_{\mathrm{BL}}\\ {\Sigma }_{\mathrm{MT}}\\ {\Delta }_{\mathrm{MT}}\end{array}\right).$ Step 4: Characterizing Local Variation in Principal Quantities with 6-DOF Motions

First order values for (φ, θ, ψ, Δx, Δx, Δy, Δz) determined in Step 3 are entered into the translation equation [1] to determine the corresponding translated 3-D target position (x_i, y_i, z_i) for each of the three target vertices. These world coordinates are projected to camera coordinates (X_i,j, Y_i,j) using equations [2] through [6], and the principal quantities are re-calculated. These six quantities are compared against the measured values of these quantities determined in Step 1, to create a residual error matrix:(σ_ΣHD,σ_ΔHD,σ_ΣVD,σ_ΔBL,σ_ΣMT,σ_ΔMT)

Local partial derivatives of the principal quantities are calculated by introducing small increments in roll, pitch, yaw, x-, y- and z-axis displacements one at a time as before, but this time the increments are relative to the first-order displacement vector. For each increment, the new world coordinates and the new camera projections of the target vertices are re-computed and the principal quantities are re-calculated. The change in each principal quantity is divided by the small angular or displacement increment to determine a local partial derivative. For instance, to calculate partial derivatives with respect to roll, the first-order displacement vector {φ₀, θ₀, ψ₀, (Δx)₀, (Δy)₀, (Δz)₀} is replaced by {φ₀+δφ, θ₀, ψ₀, (Δx)₀, (Δy)₀, (Δz)₀} and resulting changes to each of the principal quantities is divided by δφ to determine the local derivative with respect to roll. This is repeated for each of the six degrees of freedom.

Each of these six repetitions produces one column of the new local partial derivative matrix:

$\left(\begin{array}{cccccc}\frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial x}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial y}& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \chi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \gamma }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \chi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \gamma }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial z}\end{array}\right){|}_{\left\{{\phi }_0,{\theta }_0,{\psi }_0,{\left(\Delta \chi \right)}_0,{\left(\Delta \gamma \right)}_0,{\left(\Delta z\right)}_0\right\}}$ Step 5: Determining Coarse Correction to First-Order Displacement Vector

A coarse correction is computed to improve the first-order displacement vector and reduce residual error, by multiplying the residual error matrix determined in Step 4 by the inverse of the local partial derivative matrix, also determined in Step 4:

$\left(\begin{array}{c}\Delta \phi \\ \Delta \theta \\ \Delta \psi \\ \Delta \left(\Delta x\right)\\ \Delta \left(\Delta y\right)\\ \Delta \left(\Delta z\right)\end{array}\right)={\left(\begin{array}{cccccc}\frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{HD}}}{\partial \gamma }\\ \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \chi }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial \gamma }& \frac{\partial {\Delta }_{\mathrm{HD}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{VD}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \chi }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial \gamma }& \frac{\partial {\Delta }_{\mathrm{BL}}}{\partial z}\\ \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \chi }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial \gamma }& \frac{\partial {\Sigma }_{\mathrm{MT}}}{\partial z}\\ \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \phi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \theta }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \psi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \chi }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial \gamma }& \frac{\partial {\Delta }_{\mathrm{MT}}}{\partial z}\end{array}\right)}^{-1}\left(\begin{array}{c}{\sigma }_{{\Sigma }_{\mathrm{HD}}}\\ {\sigma }_{{\Delta }_{\mathrm{HD}}}\\ {\sigma }_{{\Sigma }_{\mathrm{\gamma D}}}\\ {\sigma }_{{\Delta }_{\mathrm{BL}}}\\ {\sigma }_{{\Sigma }_{\mathrm{MT}}}\\ {\sigma }_{{\Delta }_{\mathrm{MT}}}\end{array}\right)$

The first-order displacement vector is incremented by the coarse correction matrix to create a better approximation to the displacement vector:{φ₀+Δφ,θ₀+Δθ,ψ₀+Δψ,(Δx)₀+Δ(Δx),(Δy)₀+Δ(Δy),(Δz)₀+Δ(Δz)}.Step 6: Performing Fine Correction to Determine Final 6DOF Displacement Vector

Steps 4 and 5 are repeated, starting with the coarse-corrected displacement vector, to determine a final fine correction to the displacement vector. After this iteration, the resultant fine correction increments are added to the coarse-corrected vector to create the final 6-DOF displacement vector. Empirical results from a general simulation indicate that this fine correction is sufficient in all cases to reduce residual errors to well below the stated 0.1-degree, 0.1-mm tolerances.

Algorithm Numerical Simulation to Verify Absolute Convergence

As cameras, targets and rotation stages are being procured and assembled, the 6 DOF decomposition algorithm can be coded and tested for a test set of rotations. It is clear that the routine will converge for small translations and rotations, but it can be potentially advantageous to determine whether there are limitations on its convergence for extreme displacements in all six degrees of freedom. To this end, we imagine an extreme target displacement, calculate the 3D position of the displaced target, calculate the centroid positions that will be seen on each of the two cameras, and run the decomposition algorithm to determine speed of convergence.

In some embodiments, to demonstrate absolute convergence of the iterative 6DOF unfolding algorithm, the simulation is started with a test set of very large rotations and displacements, as listed in Table 1 below.

TABLE 1
Example of a set of extreme angular rotations and
linear translations of an imaginary patient for
purposes of testing algorithm convergence.
Head Yaw (Lean) Psi (deg toward patient's right shoulder) 8.0000
Head Pitch (Nod) Theta (deg relative to level; pos is toward −15.0000
top of head)
Head Roll (Shake) Phi (deg relative to square; pos toward 12.0000
patient's left side)
Head shift dx (mm toward top of head) −9.0000
Head shift dy (mm to patient's right) 3.0000
Head shift dz (mm away from table) 7.0000

The simulation begins by determining the locations of the displaced centroids that will be seen by each camera, allowing for some degree of mispointing and misalignment of each camera. The original (nominal) target location is rotated and displaced by the Euler rotation formalism presented in Section 2.5.2.2, to determine the three displaced target centroid locations in three-dimensional space. Next these “world coordinates” are translated to 2-D “camera coordinates” for each of the two cameras independently, as described in the same Section.

Assuming the target is imaged into these camera coordinates, but that the operator has no prior knowledge of the displacement matrix giving rise to this target position, we use the algorithm as described in Section 2.5.2 from end to end to recreate the displacement matrix. By the end of Step 3 (Section 2.5.2.3), the algorithm returns an initial estimate of the 6DOF displacement vector, as shown in Table 2 below.

TABLE 2
First estimate of 6DOF displacement based
on method described in Section 2.5.2.
First Approximation Yaw (degrees) 4.4313
First Approximation Pitch (degrees) −19.4474
First Approximation Roll (degrees) 8.8784
First Approximation X displacement (mm) −6.4257
First Approximation Y displacement (mm) −2.5639
First Approximation Z displacement (mm) −5.9428

As expected, residual errors at this stage are atypically large, due to the extreme magnitudes of the translations and rotations chosen for this simulation along and about each axis; this situation creates a good test for absolute convergence of the Newton Raphson algorithm methodology. Assuming this estimate to be correct, the algorithm in Step 4 (Section 2.5.2.4) again calculates the displaced position of the target, the resulting centroid positions seen by each camera, and the principal quantities (vertical tip sum and difference, base length difference, vertical displacement sum, and horizontal displacement sum and difference) which would result, for comparison with the actual observed values. The residual errors, in pixels, and the local derivatives of each of the principal values for small changes (pixels per 0.1 degrees) in yaw, pitch, and roll, and for small changes (pixels per 0.1 mm) in dx, dy and dz are calculated as described in Section 2.5.2.4, and tabulated as shown in Table 3 below.

TABLE 3
Residual Error (in pixels) and local derivatives with respect to Yaw, Pitch,
Roll (pixels per 0.1 deg), x-displacement, y-displacement, and z-displacement
(pixels per 0.1 mm), of the principal quantities Vertical Tip Sum, Vertical
Tip Difference, Base Length Difference, Vertical Displacement Sum, Horizontal
Displacement Sum, and Horizontal Displacement Difference.
							Residual
	∂/∂Y	∂/∂P	∂/∂R	∂/∂x	∂/∂y	∂/∂z	Error
VT1 + VT2	0.2575	0.0383	−0.0994	0.0021	0.0045	0.0021	−3.6558
VT1 − VT2	0.0657	−0.2756	−0.0131	0.0006	0.0018	0.0274	6.8709
BL1 − BL2	−0.3223	0.0277	0.4988	0.0109	−0.0702	0.0106	−2.9918
VDT + VD2	−0.3118	5.8134	0.0350	1.8813	0.0112	−0.2223	−168.5591
HD1 + HD2	−2.5875	−0.1680	3.8651	0.0117	−1.3090	−0.0124	58.1859
HD1 − HD2	−0.5823	1.4452	0.7697	−0.0140	−0.1114	−1.4280	120.7937

The matrix of derivatives at the left of Table 3 is inverted and multiplied by the residual error vector at the right, to yield first-order corrections to the initial estimate of the displacement vector, as described in Section 2.5.2.5, and as shown at the left of Table 4 below. These are added to the initial estimates, to produce the more refined estimate of the 6 DOF displacement vector, shown at the right of Table 4.

TABLE 4
(left) First-Order Corrections to Initial Estimates of Yaw,
Pitch, Roll, dx, dy and dz, obtained by inverting the matrix
of derivatives at left of Table 3 above and multiplying this
inverse matrix by the residual error vector at right of Table
3. These corrections are added to initial 6DOF motion estimates
to produce improved estimates at right above.
Yaw Adjustment	3.8632	First Newton Iteration Yaw	8.2945
(deg)		(deg)
Pitch Adjustment	4.5672	First Newton Iteration Pitch	−14.8803
(deg)		(deg)
Roll Adjustment	3.5642	First Newton Iteration Roll	12.4426
(deg)		(deg)
dx Adjustment	−3.0846	First Newton Iteration Delta	−9.5103
(mm)		X (mm)
dy Adjustment	6.5969	First Newton Iteration Delta	4.0329
(mm)		Y (mm)
dz Adjustment	12.9426	First Newton Iteration Delta	6.9998
(mm)		Z (mm)

This process is repeated for a second and final time as described in Section 2.5.2.6, assuming again that the (now refined) 6 DOF displacement vector is accurate, and calculating first the 3D target centroid positions and then the locations of the target centroids as projected onto each of the two camera focal planes. Again the six principal quantities are computed and compared with the actual observations to produce a vector of residual errors. Again the local derivatives are computed, this time at the location of the first-order displacement vector. The results are tabulated as shown in Table 5 below.

TABLE 5
First-Order Residual Error (in pixels) and new local derivatives with respect
to Yaw, Pitch, Roll (pixels per 0.1 deg), x-displacement, y-displacement, and
z-displacement (pixels per 0.1 mm), of the principal quantities Vertical Tip
Sum, Vertical Tip Difference, Base Length Difference, Vertical Displacement
Sum, Horizontal Displacement Sum, and Horizontal Displacement Difference.
							Residual
	∂/∂Y	∂/∂P	∂/∂R	∂/∂x	∂/∂y	∂/∂z	Error
VT1 + VT2	0.2498	0.0545	−0.0785	0.0020	0.0028	0.0007	0.4715
VT1 − VT2	0.0682	−0.2935	0.0223	−0.0012	−0.0034	0.0242	−0.0827
BL1 − BL2	−0.3146	0.0536	0.4966	0.0171	−0.0723	0.0094	0.5096
VD1 + VD2	−0.5927	5.7797	0.0405	1.9353	0.0084	−0.1911	−4.3941
HD1 + HD2	−2.5462	−0.3237	3.7395	0.0074	−1.3067	−0.0135	−4.8578
HD1 − HD2	−0.6876	1.7791	0.7547	−0.0177	−0.0884	−1.4784	2.5723

The matrix of derivatives at the left of Table 5 is inverted and multiplied by the residual error vector at the right, to yield final corrections to the first-order estimate of the displacement vector, as shown at the left of Table 6 below. These corrections are added to the first-order estimates, to produce the final second-order estimate of the 6 DOF displacement vector, shown at the right of Table 6.

TABLE 6
(left) Second-Order Corrections to First-Order Estimates of Yaw, Pitch,
Roll, dx, dy and dz, obtained by inverting the matrix of derivatives
at left of Table 5 above and multiplying this inverse matrix by the
residual error vector at right Table 5. These corrections are added
to first-order correction obtained by the same method, to produce final
values for each of the 6 DOF motions used in the simulation.
	Yaw Adjustment	−0.2947	Final Yaw	7.9999
	(deg)		(deg)
	Pitch Adjustment	−0.1210	Final Pitch	−15.0013
	(deg)		(deg)
	Roll Adjustment	−0.4448	Final Roll	11.9978
	(deg)		(deg)
	dx Adjustment	0.5114	Final Delta X	−8.9989
	(mm)		(mm)
	dy Adjustment	−1.0377	Final Delta Y	2.9952
	(mm)		(mm)
	dz Adjustment	−0.0058	Final Delta Z	6.9941
	(mm)		(mm)

Even for the extreme rotations and displacements used in this model, the algorithm is shown to converge to within 0.003 degrees and 0.006 mm in only two iterations. Given the number of floating-point operations needed to perform the initial estimate and two successive iterations of the Newton method, the algorithm can produce a solution on a typical laptop computer in less than 5 milliseconds.

Quaternion Representation

The head coil ICD specifies the rotation vector in terms of the quaternion, for which (still using right-handed Euler angle rotation conventions):

$q=\left[\begin{array}{c}{q}_r\\ q_x\\ q_y\\ q_z\end{array}\right]=\left[\begin{array}{c}\cos \left(\phi /2\right)\cos \left(\theta /2\right)\cos \left(\psi /2\right)-\sin \left(\phi /2\right)\sin \left(\theta /2\right)\sin \left(\psi /2\right)\\ -\sin \left(\phi /2\right)\cos \left(\theta /2\right)\cos \left(\psi /2\right)-\cos \left(\phi /2\right)\sin \left(\theta /2\right)\sin \left(\psi /2\right)\\ -\cos \left(\phi /2\right)\sin \left(\theta /2\right)\cos \left(\psi /2\right)+\sin \left(\phi /2\right)\cos \left(\theta /2\right)\sin \left(\psi /2\right)\\ -\cos \left(\phi /2\right)\cos \left(\theta /2\right)\sin \left(\psi /2\right)-\sin \left(\phi /2\right)\sin \left(\theta /2\right)\cos \left(\psi /2\right)\end{array}\right]$ The translation vector is unchanged from the form calculated here.Centroid Determination Algorithm

The centroid location on the focal plane is given by:

$\begin{array}{cc}{x}_c=\frac{\sum _{\mathrm{ij}}{x}_{\mathrm{ij}}{I}_{\mathrm{ij}}}{\sum _{\mathrm{ij}}{I}_{\mathrm{ij}}},& y_c=\frac{\sum _{\mathrm{ij}}{y}_{\mathrm{ij}}{I}_{\mathrm{ij}}}{\sum _{\mathrm{ij}}{I}_{\mathrm{ij}}}.\end{array}$

This calculation is performed for three subregions on the target as shown in FIG. 14C (the dashed lines do not appear on the real target), inverting the image such that large count numbers correspond to black (near 4095, for the 12-bit monochrome camera readout) and small count numbers for white (near 0). With a minimal amount of sophistication, the routine can detect the pattern of circles and approximately locate these subregions automatically. In some embodiments, the routine can be initialized with a key click to identify the approximate position of each centroid at startup. Subsequently, the three regions of interest for each new frame will be centered at the centroid locations from the previous frame, plus and minus 48 pixel rows and plus and minus 48 pixel columns. Regions of interest around each of the three target circles which can be integrated to determine target centroids.

Centroid Determination

In some embodiments, a test target can be printed and mounted in the view field of a monochrome camera at an angle of approximately 45 degrees. At this angle the elliptical target projected to an approximately round target on the camera focal plane. FIG. 16 shows the camera focused at full-scale printed target oriented at 45 degrees at a distance of 14.1 inches. Camera field of view is roughly the size of the rectangle in the center of the camera calibration target mounted next to the target.

The calculated target centroid is displayed as a red dot at the center of the LabView image in FIG. 17, and displayed as a floating point (x,y) pair to the right of the image. At illumination levels above about 20% of full scale, the measured centroid location does not fluctuate above the 0.1-pixel level in row or column; for lower intensity levels, statistical fluctuations exceed this threshold. It is noted, however, that for the black-on-white printed target, uniformity of illumination can be potentially important—if the target is illuminated significantly more strongly from the left or right side, for instance, the moment calculation could add bias in the horizontal direction and would shift the centroid outside of the specified error threshold. This effect could in some cases put an undesirable cost constraint on the illumination approach, so an intensity thresholding algorithm is first implemented, by which the target histogram is clipped near the lower extrema for the bright and dark region intensities, eliminating the undesirable effect. In some embodiments, a Camera Control screen view can allow control of camera frame rate and readout resolution, showing manually-selected region of interest. Full camera field of view is approximately represented by a black region on the screen. The centroid can be displayed as a red dot at the center of the circular target, and camera x-y coordinates are displayed as floating point numbers to 2-decimal precision to the right of the display.

EXAMPLE 1

Camera Calibration

As with any camera lens, the lens used for the head tracker could have some level of distortion as a function of distance from imaging axis. Azimuthal distortion should be negligible, but radial distortion can be measured after lens installation and fit to a polynomial curve to allow rapid compensation of centroid positions near the edges of the camera field of view. The 6DOF unfolding algorithm can be constructed to accommodate typical levels of radial distortion as a second-order compensation during the application of the Newton Raphson iteration method.

Radial distortion can be determined using a printed reference target with concentric circles of diameter ⅓″, ⅔″, 1″, and so on up to a maximum diameter of 4 inches, as shown in FIGS. 16A and 16B. The approximate FOV of the camera and 25 mm lens at a working distance of 14.08 inches is 2″×2.67″, as indicated by the inner rectangle printed on the target. The camera is mounted 14.24 inches from the target such that the inner rectangle is visible at the edges of the camera FOV, and the target is centered in this field. A single image frame is captured and the intersections of the circles and radial lines are identified and precisely located through local centroid calculations on the image. The polar angles of the world coordinate system are compared against the polar angles recorded on the camera to determine the radial distortion. FIG. 16A is the Camera Calibration Target and FIG. 16B is the off-axis radial distortion of the 25 mm fixed-focal length camera lens, measured by comparing the diameters of circles recorded on the camera focal plane.

In one embodiment, the measured radial distortion measured for the TechSpec High Resolution Fixed Focus 25 mm lens follows camera polar angle θ_c=(1+0.0053144θ−0.0016804θ²+0.0002483θ³−0.0000138θ⁴)θ, with laboratory polar angle θ in degrees. At the extreme corner of the viewing field, where θ˜6.75°, camera aberration results in a radial growth in camera angle of about 0.7% relative to true angle, or about 2.8 pixels in radius.

Full 6-DOF Tracking

The full 6-DOF tracking algorithm was coded in LabView with the Graphical User Interface (GUI). The upper left side of the GUI screen gives centroid information for target circles in the current frame, and the lower left side gives the same information for the prior frame. For each, one nested target circle from the set of three is displayed in negative (white on black) along with a histogram of its pixel brightness within a 48-by-48 pixel box centered on the centroid location of the previous frame. This histogram is split into two sections to display (at left) the peak from background pixels at one end of the brightness scale, and (at right) the peak from the pixels of the target itself, at the other end of the brightness scale. A long continuum of pixels in between represents pixels at dark-light boundaries in the target frame. From analysis of the two histograms, the target field is clipped at the lower-brightness shoulder on the bright side, and the upper brightness shoulder on the dark side, to create a binary target field that is not sensitive to variations in illumination across the target. Although displayed in real time for only one target circle, all three target circles are processed in this way.

Next to the target histograms, the x-y camera centroid locations are displayed to two-decimal precision for each of the three nested circle targets; again at the upper half of the screen for the current data and at the lower half of the screen for the prior frame.

The right side of the screen displays the processed 6-DOF data, after analysis using the approach described in Section 2.5. An analog meter-style display shows the acquisition and processing time per frame, which is limited at its low end to the camera frame integration and readout time of about 8 milliseconds. Using a single iteration of the Newton-Raphson routine described in Section 2.5, the algorithm runs during the integration period for the successive frame, so the processing time is approximately 8 milliseconds, corresponding to a 120 Hz camera readout rate. The 6-DOF data can be displayed in either analog or digital format, but the digital format can be read to precision of 0.01 mm and 0.01 degree for comparison with the 0.1 mm, 0.1 degree accuracy requirements.

Laboratory Mechanical Layout for Head Tracking Simulation

The laboratory setup was designed to mimic head rotation and displacement using a six-degree-of-freedom optical rotation mount. This mount included three ganged translation stages along the x-, y-, and z-axes of the optical table, and three ganged rotation stages corresponding to yaw, roll and pitch respectively. The two monochrome cameras and turning mirrors were mounted in the appropriate geometry for use with an existing 12-channel head coil. The two monochrome cameras are in foreground, mounted at ±45° relative to horizontal to accommodate rotation by the turning mirrors. The turning mirrors are mounted 10 inches behind cameras (slightly obscured by the cameras in the picture). The target is partially visible in the reflection of each mirror. The 6-DOF rotation stage is at center in foreground, with the y-axis stage at bottom, x-axis stage next, and z-axis stage above that, followed by the yaw rotation stage, the roll stage, and finally the pitch stage with target at the top (the pitch rotation handle is obscured by the stage). A near-IR illumination LED is at the center in background; light from this stage is within the camera spectral range, but hardly visible to the human eye.

X-Axis Translation

The second translation stage from the bottom in the 6-DOF displacement assembly controls x-axis displacement (aligned with the patient's spine). The x-axis translation stage control knob is turned four full rotations (corresponding to −2.54 mm), and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. Results are: the displacement determined by the unfolding algorithm is −2.56 mm in x, less than 0.1 mm in y and z, and less than 0.1° in roll, pitch and yaw. The target displacement by dx=−2.54 mm, with zoom on lower right display section of GUI showed calculated dx=−2.56 mm, dy=0.08 mm, dz=0.02 mm, dφ=0.05°, dθ=−0.03°, and dψ=−0.01°.

Y-Axis Translation

The bottom translation stage in the 6-DOF displacement assembly controls y-axis displacement (patient's left-to-right). The y-axis translation stage control knob is turned four full rotations (corresponding to −2.54 mm), and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. This resulted in a target displacement by dy=−2.54 mm, with zoom on lower right display section of GUI showing dx=0.00 mm, dy=−2.47 mm, dz=−0.01 mm, dφ=0.64°, dθ=−0.04°, and dψ=−0.03°.

Z-Axis Translation

The top translation stage in the 6-DOF displacement assembly controls z-axis displacement (patient's down to up, with the patient lying on his back). The z-axis translation stage control knob is turned four full rotations (corresponding to −2.54 cm), and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. The displacement determined by the unfolding algorithm was −2.54 mm in z, less than 0.1 mm in x and y, and less than 0.1° in roll, pitch and yaw. The results were a target displacement by dz=−2.54 mm, with zoom on lower right display section of GUI showing dx=0.01 mm, dy=−0.01 mm, dz=−2.59 mm, dφ=−0.02°, dθ=−0.06° and dψ=0.01°.

Yaw Rotation

The bottom rotation stage in the 6-DOF displacement assembly controls yaw rotation (patient's left shoulder-to-right shoulder lean direction). The yaw rotation stage control knob is turned by +4° degrees (heading 315° to heading 311° on stage, corresponds to movement toward right shoulder), and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. The displacement determined by the unfolding algorithm is less than 0.1 mm in dx, dy and dz, 0.1° in roll and less than 0.1° in pitch, and 3.94° in yaw. The results were a target rotation by dψ=+4.00°, with zoom on lower right display section of GUI showing dx=0.07 mm, dy=−0.05 mm, dz=0.02 mm, dφ=0.10°, dθ=−0.01°, and dψ=3.94°.

Roll Rotation

The middle rotation stage in the 6-DOF displacement assembly controls roll rotation (patient's right shoulder-to-left shoulder “head shaking” direction). The roll goniometer control knob is turned by +5° degrees, and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. The displacement determined by the unfolding algorithm is less than 0.1 mm in dx, and dz, 1.78 mm in dy, 4.97° in roll and less than 0.1° in pitch and yaw. Displacement in y is expected due to the fact that the center of rotation for the Thorlabs GNL18 goniometer stage is 44.5 mm above the mount surface, while the target is only 22 mm above the stage. For the resulting −20.5 mm lever arm, the y-displacement due to a 5° roll rotation is −(−20.5 mm)*sin(5°)=+1.79 mm, in good agreement with the measured data.

The results were a target rotation by dφ=+5.00°, with zoom on lower right display section of GUI showing dx=0.07 mm, dy=1.78 mm, dz=−0.01 mm, dφ=4.97°, dθ=−0.03°, and dψ=0.08°.

Pitch Rotation

The top rotation stage in the 6-DOF displacement assembly controls pitch rotation (patient's “nodding” direction). The pitch goniometer control knob is turned by +5° degrees, and the absolute position change is calculated from the resulting motion of the centroid camera coordinates. The calculated pitch is 4.95°, with less than 0.1° in yaw. The center of rotation for the Thorlabs GNL10 goniometer stage is 25.4 mm above the mount surface, while the target is only 6.4 mm above the stage. For the resulting −19 mm lever arm, the x-displacement due to a 5° rotation is −19 mm*sin(5°=−1.66 mm, the y-displacement is 0.00 mm, and the z-displacement is −19 mm*[1−cos(5°)]=0.07 mm. These displacements are all within 0.1 mm of measured data.

The results were a target pitch rotation by dθ=+5.00°, with zoom on lower right display section of GUI showing dx=−1.63 mm, dy=0.09 mm, dz=0.17 mm, dφ=0.21°, dθ=4.95°, and dψ=−0.07°.

VARIATIONS

Specific embodiments have been described in detail above with emphasis on medical application and in particular MRI examination of a patient's head. However, the teachings of the present invention can be utilized for other MRI examinations of other body parts where movements of up to six degrees of freedom are possible. In addition medical procedures involving imaging devices other than MRI equipment (e.g., CT, PET, ultrasound, plain radiography, and others) may benefit from the teaching of the present invention. The teachings of the present invention may be useful in many non-medical applications where tracking of a target having several degrees of freedom are possible. Some of these applications could be military applications. Furthermore, while particular algorithms are disclosed, variations, combinations, and subcombinations are also possible.

Computing System

In some embodiments, the computer clients and/or servers described above take the form of a computing system 1500 illustrated in FIG. 15, which is a block diagram of one embodiment of a computing system that is in communication with one or more computing systems 1520 and/or one or more data sources 1522 via one or more networks 1518. The computing system 1500 may be used to implement one or more of the systems and methods described herein. In addition, in one embodiment, the computing system 1500 may be configured to apply one or more of the methods and systems described herein. While FIG. 15 illustrates an embodiment of a computing system 1500, it is recognized that the functionality provided for in the components and modules of computing system 1500 may be combined into fewer components and modules or further separated into additional components and modules.

Motion Correction Control Systems

In an embodiment, the system 700 comprises a motion correction control system module 1514 that carries out the functions described herein with reference to motion correction mechanism, including any one of the motion correction methods described above. The motion correction control system module 1514 may be executed on the computing system 1500 by a central processing unit 1504 discussed further below.

In general, the word “module,” as used herein, refers to logic embodied in hardware or firmware, or to a collection of software instructions, possibly having entry and exit points, written in a programming language, such as, for example, COBOL, CICS, Java, Lua, C or C++ or Objective C. A software module may be compiled and linked into an executable program, installed in a dynamic link library, or may be written in an interpreted programming language such as, for example, BASIC, Perl, or Python. It will be appreciated that software modules may be callable from other modules or from themselves, and/or may be invoked in response to detected events or interrupts. Software instructions may be embedded in firmware, such as an EPROM. It will be further appreciated that hardware modules may be comprised of connected logic units, such as gates and flip-flops, and/or may be comprised of programmable units, such as programmable gate arrays or processors. The modules described herein are preferably implemented as software modules, but may be represented in hardware or firmware. Generally, the modules described herein refer to logical modules that may be combined with other modules or divided into sub-modules despite their physical organization or storage.

Computing System Components

In an embodiment, the computing system 1500 also comprises a workstation or other computing devices suitable for controlling and/or communicating with large databases, performing transaction processing, and generating reports from large databases. The computing system 1500 also comprises a central processing unit (“CPU”) 1504, which may comprise a conventional microprocessor. The computing system 1500 further comprises a memory 1508, such as random access memory (“RAM”) for temporary storage of information and/or a read only memory (“ROM”) for permanent storage of information, and a mass storage device 1502, such as a hard drive, diskette, or optical media storage device. Typically, the modules of the computing system 1500 are connected to the computer using a standards based bus system. In different embodiments, the standards based bus system could be Peripheral Component Interconnect (PCI), Microchannel, SCSI, Industrial Standard Architecture (ISA) and Extended ISA (EISA) architectures, for example.

The computing system 1500 comprises one or more commonly available input/output (I/O) devices and interfaces 1512, such as a keyboard, mouse, touchpad, and printer. In one embodiment, the I/O devices and interfaces 1512 comprise one or more display devices, such as a monitor, that allows the visual presentation of data to a user. More particularly, a display device provides for the presentation of GUIs, application software data, and multimedia presentations, for example. In the embodiment of FIG. 15, the I/O devices and interfaces 1512 also provide a communications interface to various external devices. The computing system 1500 may also comprise one or more multimedia devices 1506, such as speakers, video cards, graphics accelerators, and microphones, for example.

Computing System Device/Operating System

The computing system 1500 may run on a variety of computing devices, such as, for example, a mobile device or a server or a desktop or a workstation, a Windows server, an Structure Query Language server, a Unix server, a personal computer, a mainframe computer, a laptop computer, a cell phone, a personal digital assistant, a kiosk, an audio player, a smartphone, a tablet computing device, and so forth. The computing system 1500 is generally controlled and coordinated by operating system software, such as iOS, z/OS, Windows 95, Windows 98, Windows NT, Windows 2000, Windows XP, Windows Vista, Windows 7, Linux, BSD, SunOS, Solaris, or other compatible operating systems. In Macintosh systems, the operating system may be any available operating system, such as MAC OS X. In other embodiments, the computing system 1500 may be controlled by a proprietary operating system. Conventional operating systems control and schedule computer processes for execution, perform memory management, provide file system, networking, and I/O services, and provide a user interface, such as a graphical user interface (“GUI”), among other things.

Network

In the embodiment of FIG. 15, the computing system 1500 is coupled to a network 1518, such as a LAN, WAN, or the Internet, for example, via a wired, wireless, or combination of wired and wireless, communication link 1516. The network 1518 communicates with various computing devices and/or other electronic devices via wired or wireless communication links. In the embodiment of FIG. 15, the network 1518 is communicating with one or more computing systems 1520 and/or one or more data sources 1522.

Access to the motion correction control system module 1514 of the computer system 1500 by computing systems 1520 and/or by data sources 1522 may be through a web-enabled user access point such as the computing systems' 1520 or data source's 1522 personal computer, cellular phone, laptop, or other device capable of connecting to the network 1518. Such a device may have a browser module is implemented as a module that uses text, graphics, audio, video, and other media to present data and to allow interaction with data via the network 1518.

The browser module may be implemented as a combination of an all points addressable display such as a cathode-ray tube (CRT), a liquid crystal display (LCD), a plasma display, touch screen display or other types and/or combinations of displays. In addition, the browser module may be implemented to communicate with input devices 1512 and may also comprise software with the appropriate interfaces which allow a user to access data through the use of stylized screen elements such as, for example, menus, windows, dialog boxes, toolbars, and controls (for example, radio buttons, check boxes, sliding scales, and so forth). Furthermore, the browser module may communicate with a set of input and output devices to receive signals from the user.

The input device(s) may comprise a keyboard, roller ball, pen and stylus, mouse, trackball, voice recognition system, or pre-designated switches or buttons. The output device(s) may comprise a speaker, a display screen, a printer, or a voice synthesizer. In addition a touch screen may act as a hybrid input/output device. In another embodiment, a user may interact with the system more directly such as through a system terminal connected to the score generator without communications over the Internet, a WAN, or LAN, or similar network.

In some embodiments, the system 1500 may comprise a physical or logical connection established between a remote microprocessor and a mainframe host computer for the express purpose of uploading, downloading, or viewing interactive data and databases on-line in real time. The remote microprocessor may be operated by an entity operating the computer system 1500, including the client server systems or the main server system, an/or may be operated by one or more of the data sources 1522 and/or one or more of the computing systems. In some embodiments, terminal emulation software may be used on the microprocessor for participating in the micro-mainframe link.

In some embodiments, computing systems 1520 that are internal to an entity operating the computer system 1500 may access the motion correction control system module 1514 internally as an application or process run by the CPU 1504.

User Access Point

In an embodiment, the computing system 1500 comprises a computing system, a smartphone, a tablet computing device, a mobile device, a personal computer, a laptop computer, a portable computing device, a server, a computer workstation, a local area network of individual computers, an interactive kiosk, a personal digital assistant, an interactive wireless communications device, a handheld computer, an embedded computing device, or the like.

Other Systems

In addition to the systems that are illustrated in FIG. 15, the network 1518 may communicate with other data sources or other computing devices. The computing system 1500 may also comprise one or more internal and/or external data sources. In some embodiments, one or more of the data repositories and the data sources may be implemented using a relational database, such as DB2, Sybase, Oracle, CodeBase and Microsoft® SQL Server as well as other types of databases such as, for example, a signal database, object-oriented database, and/or a record-based database.

Conditional language, such as, among others, “can,” “could,” “might,” or “may,” unless specifically stated otherwise, or otherwise understood within the context as used, is generally intended to convey that certain embodiments include, while other embodiments do not include, certain features, elements and/or steps. Thus, such conditional language is not generally intended to imply that features, elements and/or steps are in any way required for one or more embodiments or that one or more embodiments necessarily include logic for deciding, with or without user input or prompting, whether these features, elements and/or steps are included or are to be performed in any particular embodiment. The headings used herein are for the convenience of the reader only and are not meant to limit the scope of the inventions or claims.

Although this invention has been disclosed in the context of certain preferred embodiments and examples, it will be understood by those skilled in the art that the present invention extends beyond the specifically disclosed embodiments to other alternative embodiments and/or uses of the invention and obvious modifications and equivalents thereof. Additionally, the skilled artisan will recognize that any of the above-described methods can be carried out using any appropriate apparatus. Further, the disclosure herein of any particular feature, aspect, method, property, characteristic, quality, attribute, element, or the like in connection with an embodiment can be used in all other embodiments set forth herein. For all of the embodiments described herein the steps of the methods need not be performed sequentially. Thus, it is intended that the scope of the present invention herein disclosed should not be limited by the particular disclosed embodiments described above.

Claims

1. A system for tracking a moving target having up to six degrees of freedom and rapidly determining positions of the moving target, said system comprising: A) A precision optical target fixed to the moving target,B) at least two cameras positioned so as to view the optical target from different directions with each of the at least two cameras being adapted to record two dimensional images of the precision optical target defining a precise target point,C) a computer processor programmed to determine a target position in Cartesian coordinates of x, y, and z and pitch, roll and yaw utilizing an algorithm adapted to: 1) identify the precise target point on the precision optical target and the x, y and z displacement of the precise target point based on optical images collected by the at least two cameras,2) utilize an iteration procedure whereby an approximate first-order solution is proposed and tested against the identified precise target points to determine residual errors which are then divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction,3) repeat step 2 until residual error becomes smaller than desired accuracy, and4) utilize the results of steps 1), 2) and 3) to determine the position of the target at rates of at least 100 times per second.

2. The system as in claim 1 wherein the moving target is a human head.

3. The system as in claim 1 and configured to interface as a component of an MRI device.

4. The system as in claim 1 wherein the iteration procedure is a variant of the Newton-Raphson method.

5. The system as in claim 2 wherein movements are measured relative to a pivot point in the patient's neck.

6. The system as in claim 2 wherein measurements are updated at a rate of at least 100 solutions per second with a latency of less than 10 milliseconds.

7. The system as in claim 2 wherein measurements are updated at a rate of at least 200 solutions per second with a latency of less than 10 milliseconds.

8. The system as in claim 3 wherein the system is adapted to report to MRI systems the exact position of the head with accuracies better than 0.1 mm in distances and 0.1 degree in angles.

9. The system as in claim 1 wherein the precision optical target comprises at least three concentric sub-targets.

10. The system as in claim 2 wherein the precision optical target is fixed to at least one of the patient's upper teeth.

11. The system as in claim 9 wherein the three concentric sub-targets are each concentric ellipses.

12. The system as in claim 9 wherein the computer processor is programmed to calculate centroids of each sub-target by dividing (a) a sum of a product of pixel intensity and pixel position by (b) a sum of pixel intensity in a sub-pixel array.

13. The system as in claim 11 wherein the pixel array is about 48×48 pixels.

14. A method for tracking a moving target having up to six degrees of freedom and rapidly determining positions of the target comprising the steps of: A) attaching a precision optical target fixed to the moving target,B) positioning at least two cameras positioned so as to view the optical target from different directions with each of the at least two cameras being adapted to record two dimensional images of the precision optical target defining a precise target point,C) programming a computer processor to determine the target position in Cartesian coordinates of x, y and z and pitch, roll and yaw utilizing an algorithm adapted to 1) identify a precise target point on the precision optical target and the x, y and z displacement of the precise target point based on optical images collected by the at least two cameras,2) utilize an iteration procedure whereby an approximate first-order solution is proposed and tested against the identified precise target points to determine residual errors which are then divided by the local derivatives with respect to each component of rotation and translation, to determine an iterative correction,3) repeat step 2 until residual error becomes smaller than desired accuracy, and4) utilize the results of steps 1), 2) and 3) to determine the position of the target at rates of at least 100 times per second.

15. The method as in claim 14 wherein the target is a human head.

16. The method as in claim 15 wherein tracking results are utilized as an input to an MRI device so as to adjust its magnetic fields to compensate for movements of the human head.

17. The method as in claim 14 wherein the iteration procedure is a variant of the Newton-Raphson method.

18. The method as in claim 15 wherein movements are measured relative to a pivot point in the patient's neck.

19. The method as in claim 14 wherein measurements are up-dated at a rate of at least 100 solutions per second with a latency of less than 10 milliseconds.

20. The method as in claim 14 wherein measurements are up-dated at a rate of at least 200 solutions per second with a latency of less than 10 milliseconds.

21. The method as in claim 16 wherein the system is adapted to report to MRI systems the exact position of the head with accuracies better than 0.1 mm in distances and 0.1 degree in angles.

22. The method as in claim 14 wherein the precision optical target comprises at least three concentric sub-targets.

23. The method as in claim 22 wherein the three concentric sub-targets are each concentric ellipses.

24. The method as in claim 22 wherein the computer processor is programmed to calculate centroids of each sub-target by dividing (a) a sum of a product of pixel intensity and pixel position by (b) a sum of pixel intensity in a sub-pixel array.

25. The method as in claim 24 wherein the pixel array is about 48×48 pixels.

26. The system as in claim 4 wherein the computer processor is programmed with algorithms adapted to perform the iteration procedure in steps comprising the following steps: I) Locate the positions of the three target centroids on each of the two Cameras and calculate six principal quantities: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT;II) Make an initial guess of (0, 0, 0, 0, 0, 0) for the subject displacements and rotations (φ, θ, ψ, Δx, Δy, Δz) leading to those centroid positions;III) Enter the guess values for (φ, θ, ψ, Δx, Δy, Δz) into a translation matrix and calculate the corresponding translated 3-D target position (xi, yi, zi) for each of the three target centroids;IV) Calculate the position (horizontal and vertical pixel number) of the projection of each of the three target centroids calculated in Step III on each camera system;V) Calculate six principal quantities using Step IV data from the two cameras: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT;VI) Compare these six calculated principal quantities with measured values from Step I, List the differences as a matrix of errors in the guessed/calculated quantities relative to the measured quantities: (σΣHD,σΔHD,σΣVD,σΔBL,σΣMT,σΔMT);VII) Determine the local partial derivatives of: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT by repeating Steps III-V six times with small added displacements in each degree of freedom, one at a time;VIII) Determine a coarse correction matrix by dividing the error matrix by the derivative matrix, to improve the initial guess and reduce the residual error, Create a better displacement matrix: [φ+Δφ,θ+Δθ,ψ+Δψ,Δx+Δ(Δx),Δy+Δ(Δy),Δz+Δ(Δz)];IX) Repeat steps III through VIII for a second and final iteration, starting with the coarse-corrected displacement matrix computed in Step VIII as the input for Step III; andX) After each successive camera frame increment, repeat Steps I-IX but use the result of previous Step IX for the guess in Step II.

27. The method as in claim 14 wherein the computer processor is programmed with algorithms adapted to perform the iteration procedure in steps comprising the following steps: I) Locate the positions of the three target centroids on each of the two Cameras and calculate six principal quantities: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT;II) Make an initial guess of (0, 0, 0, 0, 0, 0) for the subject displacements and rotations (φ, θ, ψ, Δx, Δy, Δz) leading to those centroid positions;III) Enter the guess values for (φ, θ, Δx, Δx, Δy, Δz) into a translation matrix and calculate the corresponding translated 3-D target position (xi, yi, zi) for each of the three target centroids;IV) Calculate the position (horizontal and vertical pixel number) of the projection of each of the three target centroids calculated in Step III on each camera system;V) Calculate six principal quantities using Step IV data from the two cameras: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT;VI) Compare these six calculated principal quantities with measured values from Step I, List the differences as a matrix of errors in the guessed/calculated quantities relative to the measured quantities: (σΣHD,σΔHD,σΣVD,σΔBL,σΣMT,σΔMT);VII) Determine the local partial derivatives of: ΣHD ΔHD ΣVD ΔBL ΣMT ΔMT by repeating Steps III-V six times with small added displacements in each degree of freedom, one at a time;VIII) Determine a coarse correction matrix by dividing the error matrix by the derivative matrix, to improve the initial guess and reduce the residual error, Create a better displacement matrix: [φ+Δφ,θ+Δθ,ψ+Δψ,Δx+Δ(Δx),Δy+Δ(Δy),Δz+Δ(Δz)];IX) Repeat steps III through VIII for a second and final iteration, starting with the coarse-corrected displacement matrix computed in Step VIII as the input for Step III; and