Online Calibration of a Linear Micro Tomosynthesis Scanner

In a linear tomosynthesis scanner designed for imaging histologic samples of several centimeters size at 10 µm resolution, the mechanical instability of the scanning stage (±10 µm) exceeded the resolution of the image system, making it necessary to determine the trajectory of the stage for each scan to avoid blurring and artifacts in the images that would arise from the errors in the geometric information used in 3D reconstruction. We present a method for online calibration by attaching a layer of randomly dispersed micro glass beads or calcium particles to the bottom of the sample stage. The method was based on a parametric representation of the rigid body motion of the sample stage-marker layer assembly. The marker layer was easy to produce and proven effective in the calibration procedure.


Adaptation of Linear Tomosynthesis to Microscopy of Pathology Samples
Tomosynthesis is an imaging technique that involves taking a series of X-ray images at different projection angles and reconstructing those images to form a stack of crosssectional images at a range of depths [1]. It has the advantage of being able to resolve overlapping structures at different depths with relatively simple hardware. In the field of breast imaging, tomosynthesis has helped overcome challenges such as tissue superposition and false positive rates, thereby becoming widely accepted today [2,3]. Examples of other clinical applications besides mammography have also been published [4][5][6]. A particular type of linear tomosynthesis scan involves moving the sample down a straight line parallel to a flat panel detector [7][8][9][10]. It is related in concept to an early form of tomography called planigraphy [11]. In the field of non-destructive testing, it is sometimes called laminography, because it is particularly suited to scanning objects of flat shapes [7,8,10]. While this type of tomosynthesis was developed for materials testing and security screening, we were the first to realize that it was also advantageous for pathology samples that have flat shapes, such as standard paraffin embedding cassettes and tissue in Petri dishes. Therefore, we adapted it in the form of a linear micro tomosynthesis scanner for microscopy of tissue samples [12,13]. It provided 7 to 10 µm resolution with good soft-tissue contrast in a 15 min scan, which was an effective scouting tool to guide thin sectioning, staining and light microscopy [12,13]. The scanner is an X-ray cabinet system ( Figure 1). The X-ray hardware is stationary, whereas the sample is carried on a motorized stage and scanned across a wide cone beam. The scanning geometry allows flat samples to be laid close to the X-ray source, to increase the geometric magnification and beam intensity through the sample. These factors then help to improve resolution and signal-to-noise ratio. Image reconstruction in tomosynthesis requires knowing the projection matrix from the 3D coordinates in the sample space to the 2D coordinates in the projection image. In practice, the geometry of the X-ray hardware and the sample deviate from the ideal design, and thus a calibration procedure is used to determine the deviations.
practice, the geometry of the X-ray hardware and the sample deviate from the ideal design, and thus a calibration procedure is used to determine the deviations. Figure 1. The X-ray micro tomosynthesis scanner with the radiation enclosure. The outer dimensions of the enclosure are 45 × 52 × 66 cm (width × depth × height). The motorized sample stage moves nominally in the horizontal plane. The sample is scanned horizontally across the vertical cone beam in either x or y direction. Flat samples are usually scanned in a horizontal plane near the X-ray focal spot to maximize the photon flux density through the sample and the geometric magnification for greater contrast-to-noise ration and resolution.

What Is New about the Present Calibration Method
Several recent methods have been published for the calibration of clinical tomosynthesis imaging systems [14][15][16][17][18][19]. These are based on offline scans of calibration phantoms. What is different about the current scanner is that the mechanical instability of the motorized sample stage (±10 µm) exceeded the resolution of the imaging system, and additionally the instability fluctuated with the scan range and speed. It meant that the geometric deviations were not repeatable, and needed to be calibrated online for each scan.
As a brief review of the literature, there is a wealth of calibration methods for conebeam computed tomography (CT) and tomosynthesis imaging in diverse applications . They can broadly be separated into offline phantom-based calibration methods and online phantom-less methods [39]. Graetz [47] and Jiang et al. [17] provided representative literature reviews of offline methods. For offline methods, the geometry of the imaging system is determined with dedicated scans of calibration phantoms. Calibration phantoms contain distinct positional markers, such as radio-opaquebeads [15,18,34,35,41,42,45,49] or wires [44]. Some phantom-based methods do not require Figure 1. The X-ray micro tomosynthesis scanner with the radiation enclosure. The outer dimensions of the enclosure are 45 × 52 × 66 cm (width × depth × height). The motorized sample stage moves nominally in the horizontal plane. The sample is scanned horizontally across the vertical cone beam in either x or y direction. Flat samples are usually scanned in a horizontal plane near the X-ray focal spot to maximize the photon flux density through the sample and the geometric magnification for greater contrast-to-noise ration and resolution.

What Is New about the Present Calibration Method
Several recent methods have been published for the calibration of clinical tomosynthesis imaging systems [14][15][16][17][18][19]. These are based on offline scans of calibration phantoms. What is different about the current scanner is that the mechanical instability of the motorized sample stage (±10 µm) exceeded the resolution of the imaging system, and additionally the instability fluctuated with the scan range and speed. It meant that the geometric deviations were not repeatable, and needed to be calibrated online for each scan.
As a brief review of the literature, there is a wealth of calibration methods for conebeam computed tomography (CT) and tomosynthesis imaging in diverse applications . They can broadly be separated into offline phantom-based calibration methods and online phantom-less methods [39]. Graetz [47] and Jiang et al. [17] provided representative literature reviews of offline methods. For offline methods, the geometry of the imaging system is determined with dedicated scans of calibration phantoms. Calibration phantoms contain distinct positional markers, such as radio-opaquebeads [15,18,34,35,41,42,45,49] or wires [44]. Some phantom-based methods do not require knowledge of the positions of the markers in the phantom [23,25,35,36,45,47]. The second category of methods are online phantom-less calibration. These are based on the data from the sample scan itself, and do not require separate calibration scans or calibration phantoms [21,[29][30][31][32]39,40]. To our knowledge, there were no published methods for online calibration of tomosynthesis scanners. Existing online calibration methods are phantom-less methods for CT scans of the full 180 • or 360 • projection angles. Unlike CT scans, tomosynthesis collects data from a limited range of projection angles and has inherent artifacts associated with the limited range. Therefore, it can be unreliable to measure geometric errors from artifacts or blurring of the sample itself, unless the sample contains an abundance of edges or points. Most of our pathology samples were soft tissue specimens that lack such sharp features.
From the above considerations, we created a hybrid method of online phantom-based calibration for tomosynthesis scans. The phantom is a dispersed layer of point makers (glass micro beads) placed at a sufficient distance below the sample platform, such that any tomosynthesis artifacts from the markers are confined in depth and do not interfere with the sample images. Our calibration algorithm has similarities to that of Stevens et al. for a circular-trajectory tomosynthesis scanner [23]. We operate on unfiltered back-projection images of the marker layer at and near focus. From the distribution of local deviations, we directly calculate the geometric parameters in a parametrized version of the sample stage movement. The parameters were then fed into the tomosynthesis reconstruction to improve the quality of the sample reconstruction.
The goal of the calibration is to reduce image blurring and preserve resolution. Since errors in the assumed movement of the sample stage lead directly to misalignment among the back-projected images in the image reconstruction process, which would lead directly to blurring when these back-projections were summed together to give the reconstructed image, calibration is verified based on the level of misalignment among all the backprojected images. The closer is the calibrated movement to the true movement, the less misalignment remains. The calibration is sufficient when the misalignment falls below the system resolution.

Scanner Hardware Geometry and Scan Procedure
Referring to Figure 1, the scanner consists of a stationary micro-focus X-ray tube and a stationary flat panel detector. In between the two, a motorized sample stage travels in a plane that is ideally parallel to the image plane of the detector, and in directions that are ideally parallel to the rows or columns of the pixel matrix of the detector. The geometric magnification factor of the system is generally between 10× and 16× depending on the sample thickness [13]. The X-ray tube had a focal spot of 5 µm at the operating condition of 30 kV/160 µA. The source-to-detector distance was 126 mm, and the source-to-sample distance was typically between 8 and 20 mm depending on the sample size. The flat panel X-ray detector had a pixel size of 74.8 µm and matrix size of 3072 × 3840 [13].
During a sample scan, the stage moves across the X-ray cone-beam at a constant speed in a straight line. The detector acquires projection images at a constant frame rate throughout the scan. As the sample traverses the cone beam, each point in the sample experiences X-rays of continually changing directions, which is equivalent to acquiring projection images from a range of angles. The maximum tomosynthesis angle is thus the angle subtended by the cone beam. A typical setting of scan range and speed for samples in standard tissue-embedding cassettes was 24.5 mm at 0.0272 mm/s, with a 15 min scan time.
Weighted filtered back-projection was used to reconstruct z-stacks of cross-sectional images at user specified ranges of depth [13]. A parametrized version of the sample stage movement was part of the input for the image reconstruction, which is described below.

Theoretical Basis of the Calibration Method
The basic assumption is that the sample and the sample stage move as a rigid body during the scan. Referring to Figure 2, geometric calibration involves the transformation between two coordinate systems. The first is the stationary scanner coordinate system (x, y, z), with its origin at the X-ray focal spot, the Z axis pointing perpendicular to the detector, and the X and Y axes parallel to the rows and columns of the detector image matrix. The second is the sample coordinate system (x s , y s , z s ), which travels with the sample stage. It is defined to coincide with the stationary scanner coordinate system at the midpoint of the scan. The ideal scan is a translational movement along the X axis. Generally, small deviations from the ideal movement have six degrees of freedom. These are functions of the travel distance l of the stage, including the 3D positional errors d x (l), d y (l) and d z (l), and small rotations of the sample stage represented by the three Euler angles α(l), β(l) and γ(l). By definition of the coordinate systems, the deviations are zero at the midpoint of the scan. We define l = 0 at the midpoint. To the leading order in the deviations, the general form of the transformation matrix from the sample to the scanner coordinates is

Theoretical Basis of the Calibration Method
The basic assumption is that the sample and the sample stage move as a rigid body during the scan. Referring to Figure 2, geometric calibration involves the transformation between two coordinate systems. The first is the stationary scanner coordinate system (x, y, z), with its origin at the x-ray focal spot, the Z axis pointing perpendicular to the detector, and the X and Y axes parallel to the rows and columns of the detector image matrix. The second is the sample coordinate system (xs, ys, zs), which travels with the sample stage. It is defined to coincide with the stationary scanner coordinate system at the midpoint of the scan. The ideal scan is a translational movement along the X axis. Generally, small deviations from the ideal movement have six degrees of freedom. These are functions of the travel distance l of the stage, including the 3D positional errors dx(l), dy(l) and dz(l), and small rotations of the sample stage represented by the three Euler angles α(l), β(l) and γ(l). By definition of the coordinate systems, the deviations are zero at the midpoint of the scan. We define l = 0 at the midpoint. To the leading order in the deviations, the general form of the transformation matrix from the sample to the scanner coordinates is The definition of the sample coordinate system (Xs, Ys, Zs) and the stationary scanner coordinate system (X, Y, Z). The sample coordinate system is attached rigidly to the sample stage and moves with it during the scan. The two coordinate systems coincide with each other at the mid time point of the scan. The glass microbead layer serving as fiducial markers for online geometric calibration is rigidly attached to the sample stage, typically at 1 cm below the sample. The scan track may not be perfectly aligned with the X-axis and the sample stage may rotate slightly during the scan. An exaggerated version of such deviations is illustrated.

Figure 2.
The definition of the sample coordinate system (X s , Y s , Z s ) and the stationary scanner coordinate system (X, Y, Z). The sample coordinate system is attached rigidly to the sample stage and moves with it during the scan. The two coordinate systems coincide with each other at the mid time point of the scan. The glass microbead layer serving as fiducial markers for online geometric calibration is rigidly attached to the sample stage, typically at 1 cm below the sample. The scan track may not be perfectly aligned with the X-axis and the sample stage may rotate slightly during the scan. An exaggerated version of such deviations is illustrated.
The deviations can be expanded in Taylor series of the travel distance l. Based on experimental data explained below, only the linear terms in l were significant relative to the resolution of our image system. Thus, the deviations are simplified to the form It is further simplified by dropping the term δ x l in the X direction, which is a scaling of the reconstructed images and does not affect image resolution. The δ y and δ z represent lateral deviations of the scan direction from the ideal X direction, and the ω's represent the roll, pitch and yaw of the sample stage.
The projection matrix from the 3D sample coordinates to the 2D coordinates of the detector is given by where SID is the perpendicular distance from the focal spot to the image plane in the detector. Substituting Equation (2) into (1), and the results into (3), keeping only first order terms of the deviations and l, the 2D coordinates in the image plane are given by the expression Equation (4) shows that for any fixed point in the sample, its projected 2D trajectory on the image plane may deviate from the ideal horizontal direction. The in-plane deviation angle varies with the 3D coordinates of the point in the sample coordinate system, which is given by the expression: Equation (5) is the theoretical basis for the present calibration method. The next section explains how Equation (5) was observed experimentally for each sample scan, and how it was used to determine the parameters of the sample stage movement.

Fabrication of a Layer of Dispersed Markers
We used three protocols to create a dispersed layer of markers to be attached to the sample stage. One was dispersing glass microbeads of 20 µm diameter on the surface of a 1 mm thick acrylic plate by suspending both the microbeads and plate in water. The microbeads covered a square of 10 cm in size on the plate. The plate was rigidly attached to the sample stage at a level of 10 mm below the bottom surface. A second method was gluing a sheet of printer paper to the surface of an acrylic plate, whereby the cellulose fibers of the paper were the markers. The third method was dispersing a layer of hydroxyapatite powder on the acrylic plate and covering it with a 10 cm wide adhesive tape. Both the glass microbead and hydroxyapatite powder procedures produced a satisfactory marker layer, although the hydroxyapatite procedure was simpler. The printer paper had overly strong X-ray contrasts that made calculations unreliable.

Measurement of the Geometric Parameters
The overall geometric calibration procedure is a linear chain of 3 steps, consisting of: (A) measuring the in-plane tilt angle of the scan trajectories of the markers at different locations in the marker layer. This was carried out at a rectilinear grid of locations covering the entire field-of-view of the scan; (B) fitting the measurements of step A to the parametric model of Equation (5), from which the geometric parameters of the scan movement were extracted.; (C) applying the geometric parameters obtained in step B to image reconstruction of the sample. The following is a detailed description of the steps.
To illustrate the quantity that is being measured in the first step, a z-stack of images centered at the depth of the marker layer was reconstructed assuming no geometric deviations. The z-stack covered a z range of 1 to 2 mm at 5 µm increment with an in-plane pixel size of 7 µm. In slightly off-focus images of the marker layer, the markers defocused into short line segments (Figure 3). The line segments were straight but tilted from the horizontal direction. The tilt represents an in-plane deviation of the projected trajectory of the bead from the X axis. The straightness of the line segments meant that the deviations of the scan movement were approximately linear with respect to the scan position. This was the basis for the linear approximation expressed by Equation (2). The in-plane deviation angle varied with the location of the beads, which is illustrated in Figure 3. This was anticipated by Equation (5). locations in the marker layer. This was carried out at a rectilinear grid of locations covering the entire field-of-view of the scan; (B) fitting the measurements of step A to the parametric model of Equation (5), from which the geometric parameters of the scan movement were extracted.; (C) applying the geometric parameters obtained in step B to image reconstruction of the sample. The following is a detailed description of the steps.
To illustrate the quantity that is being measured in the first step, a z-stack of images centered at the depth of the marker layer was reconstructed assuming no geometric deviations. The z-stack covered a z range of 1 to 2 mm at 5 µm increment with an in-plane pixel size of 7 µm. In slightly off-focus images of the marker layer, the markers defocused into short line segments (Figure 3). The line segments were straight but tilted from the horizontal direction. The tilt represents an in-plane deviation of the projected trajectory of the bead from the X axis. The straightness of the line segments meant that the deviations of the scan movement were approximately linear with respect to the scan position. This was the basis for the linear approximation expressed by Equation (2). The in-plane deviation angle varied with the location of the beads, which is illustrated in Figure 3. This was anticipated by Equation (5). The first step of the calibration procedure was to obtain the in-plane deviation angle δp of the marker trajectories as a function of the location of the markers. For this purpose, the entire marker layer was divided into a rectangular grid of 60 to 200 grid points. At each grid point, a small z-stack of 400 images covering a volume of 1.4 × 1.4 mm in plane and 2 mm depth was reconstructed around the z depth of the marker layer. The reconstruction was repeated while assuming a range of in-plane tilt angles of the scan trajectory, resulting in multiple small z-stacks. The in-plane deviation angle δp at each grid point was determined as the scan trajectory angle that maximized the sharpness of the markers in The first step of the calibration procedure was to obtain the in-plane deviation angle δ p of the marker trajectories as a function of the location of the markers. For this purpose, the entire marker layer was divided into a rectangular grid of 60 to 200 grid points. At each grid point, a small z-stack of 400 images covering a volume of 1.4 × 1.4 mm in plane and 2 mm depth was reconstructed around the z depth of the marker layer. The reconstruction was repeated while assuming a range of in-plane tilt angles of the scan trajectory, resulting in multiple small z-stacks. The in-plane deviation angle δ p at each grid point was determined as the scan trajectory angle that maximized the sharpness of the markers in the local z-stack. For each z-stack, marker sharpness was measured by the standard deviation of the image pixel values in the slice where the markers come into focus. That slice was determined as the one with the highest standard deviation value. At each grid point, among the multiple z-stacks of different scan trajectory angles, the scan trajectory angle giving the maximal marker sharpness was the sought-after deviation angle δ p for that grid point. This process was repeated for all grid points, and collectively provided the function δ p (x s , y s , z s ) within the marker layer. Examples of the function is shown in the surface plots of Figure 4. the local z-stack. For each z-stack, marker sharpness was measured by the standard deviation of the image pixel values in the slice where the markers come into focus. That slice was determined as the one with the highest standard deviation value. At each grid point, among the multiple z-stacks of different scan trajectory angles, the scan trajectory angle giving the maximal marker sharpness was the sought-after deviation angle δp for that grid point. This process was repeated for all grid points, and collectively provided the function δp(xs, ys, zs) within the marker layer. Examples of the function is shown in the surface plots of Figure 4. In the next step of the calibration, the measured function δp(xs, ys, zs) from the marker layer was fitted to Equation (5) with polynomial fitting. The marker layer had a single depth zs, which reduced Equation (5) to a second-order function of the in-plane coordinates xs and ys. The fitted coefficients of the 6 terms in Equation (5) then provided all the geometric parameters of the sample stage movement, including δy, δz, ωx, ωy, ωz. Lastly, these parameters were applied in the reconstruction of the sample images.

Assessing the Effect of Calibration by the Misalignment of Back-Projected Images
Since the goal of the calibration procedure is to reduce image blurring and improve image resolution, the effect of calibration is evaluated by the range of misalignment among the back-projected images, which is the direct cause of image blurring when the back-projected images are summed to produce the reconstructed image. Experimentally the misalignment is seen as linear drifts of the back-projected marker positions (Figure 3). The horizontal drift is equivalent to a slight change in the reconstructed depth of the markers and does not cause image blurring (see Equation (4)). The vertical drift causes image blurring. The range of the vertical drift, or vertical misalignment, is measured as ( , , ) = ( ( ( = 0), ( = 0)), ( ( = ), ( = ))), In the next step of the calibration, the measured function δ p (x s , y s , z s ) from the marker layer was fitted to Equation (5) with polynomial fitting. The marker layer had a single depth z s , which reduced Equation (5) to a second-order function of the in-plane coordinates x s and y s . The fitted coefficients of the 6 terms in Equation (5) then provided all the geometric parameters of the sample stage movement, including δ y , δ z , ω x , ω y , ω z . Lastly, these parameters were applied in the reconstruction of the sample images.

Assessing the Effect of Calibration by the Misalignment of Back-Projected Images
Since the goal of the calibration procedure is to reduce image blurring and improve image resolution, the effect of calibration is evaluated by the range of misalignment among the back-projected images, which is the direct cause of image blurring when the backprojected images are summed to produce the reconstructed image. Experimentally the misalignment is seen as linear drifts of the back-projected marker positions (Figure 3). The horizontal drift is equivalent to a slight change in the reconstructed depth of the markers and does not cause image blurring (see Equation (4)). The vertical drift causes image blurring. The range of the vertical drift, or vertical misalignment, is measured as M y (x s , y s , z s ) = z s SID Disp y I x p (l = 0), y p (l = 0) , I x p (l = L), y p (l = L) , (6) where in addition to the variables previously defined in Section 2.2, L is the total travel distance of the sample stage, x p and y p are the projected coordinates in the image plane based on the assumed movement of the sample stage through Equations (1)- (3), and Disp y is the y component of the displacement between the two images. This measurement was obtained for each small window in the field of view as defined in Section 2.4.

Results
As described in the Methods section, the key intermediate step of the calibration process is to determine the in-plane deviation angles of the marker trajectories from the prescribed scan direction, as a function of the location of the markers. Examples of the results for three different scan settings are shown Figure 4A-C. These figures are surface plots of the measured in-plane deviation angle as a function of the x s and y s coordinates of the markers. Variability of these plots illustrates the instability of the geometric parameters. Table 1 summarizes the estimated geometric parameters for the three scan settings and a repeat scan of the first setting. The RMS of the residual of the deviation angles after subtraction of the model fit of Equation (5) was less than 7.5% of the measured values, indicating that the model approximates the actual distribution. The residual deviation angles of the marker trajectories as a function of the marker location are also plotted in Figure 4D,E. Also summarized in Table 1 are the RMS of the range of misalignment of the marker positions among the back-projected images over the course of the scan, before and after applying the geometric parameters. The measurement is defined in Section 2.5 above. In all three scan settings, the range of misalignment was less than 7.2% of the un-calibrated levels. The post-calibration misalignment was below the system resolution of 7 µm in all settings. The effect of the calibration procedure on the reconstructed images is demonstrated in two samples. The first is a plexiglass plate coated with a dispersed layer of hydroxyapatite particles, and the second a paraffin-embedded mouse heart sample containing calcified atherosclerotic plaques in the aorta. The results are summarized in Figures 5 and 6, respectively. In both cases, we observed an increase in the resolution of the reconstructed images after the calibration was applied.

Discussion
Since the tomosynthesis scanner described in this paper was designed to image tissue samples at 10 µm resolution, it was necessary to determine the precise movement and rotation of the sample stage each time a scan was performed, such that image reconstruction was corrected for minor mechanical instabilities. These instabilities were found to fluctuate from scan to scan (Table 1) and were significant enough to cause artifacts and blurriness in the resulting images if not corrected. Although the improvement can be subtle as in the example of Figure 6, the calibration procedure was needed to reach the true image resolution provided by the system hardware.
Although the calibration method used a layer of point markers attached to the bottom of the sample stage, no particular pattern of the markers was required so far as they were dispersed throughout the imaging field of view. In practice, a layer of dispersed hydroxyapatite powder on a thin acrylic plate was relatively simple to make and effective for the procedure. A necessary condition is that the marker layer is rigidly attached to the sample stage, such that the two move as a rigid body. Another key point of the method is that there is adequate vertical separation between the sample and the marker layer to avoid mutual interference. For tissue samples in standard embedding cassettes, a 10 mm separation was found to be sufficient. The limitation of this method is the linear parametric model of the sample stage translation and rotation. It limits the utility of this method to situations where the scan movement is not significantly non-linear with respect to the distance of travel.