Calculation of C-arm Orbit Information

Introduction and Motivation

Computer assisted surgical systems that rely on 2D fluoroscopy to perform 3D navigation often rely on multiple views to recover 3D position and orientation information of important objects, such as tools or bony anatomy. Unfortunately, most C-arms found in the operating theater are not encoded and the relative pose information needed to establish a multiple-view geometry is not available. Although this information may be recovered using external fiducial objects or by using the patient's anatomy as the fiducial object, these approaches come at the cost of additional surgical equipment, potentially obscured diagnostic fluoroscopic fields of view, or the requirement for some manual, intraoperative, annotation of anatomy in order properly initialize a 6-degrees-of-freedom (DoF), rigid, 2D/3D registration of each 2D image. Assuming multiple views are collected by only undergoing an orbital rotation with respect to each other, then it is often feasible to perform patient-as-fiducial recovery of the multiple-view geometry in a computationally efficient manner while also minimizing any manual interaction from the user. This may be accomplished by registering a single-view well, potentially leveraging some manual annotation of landmarks for initialization of a 6-DoF search, and then subsequently optimizing over a single orbital rotation angle (1-DoF) for each subsquent view.

This intraoperative strategy is employed extensively by Grupp, et al. 2020 in order to track the motion of osteotomy bone fragments. The orbit information computed and used in Grupp, et al. 2020 was calculated at the start of a 16 month period of experimentation and remained sufficiently well-calibrated over the entire duration to successfully execute the registration strategy. Moreover, multiple other unrelated experiments were conducted throughout this period, which involved affixing structures of non-trivial weight to the X-ray detector and relocation of the C-arm unit amongst different laboratory areas.

This example walks through the process of calculating the orbital rotation axis and center of orbital rotation for a C-arm using:

• a fiducial object,
• an existing CT scan of the fiducial object,
• and 2D fluoroscopic views of the fiducial object collected using an orbital rotation of the C-arm.

The general strategy is as follows:

1. Collect CT of a phantom object,
2. Roughly annotate point landmarks in 3D object,
3. Collect a batch of 2D fluoroscopy views of the object solely through the use of orbital rotation,
4. Perform a 2D/3D registration of the object to each view,
5. Recover the multiple-view geometry by using the coordinate frame of the fiducial object,
6. Using the relative pose information between views, estimate the center of orbital rotation and orbital rotation axis.

This pipeline does not make assumptions on the specific geometry of the fiducial object. Since general point-landmark features and image-intensities are used to perform registrations, this strategy should be robust to many different types of fiducial objects, making it ideal for use in a laboratory environment.

Running the Example

Acquiring the Data

Download the sawbones dataset here. The 10_sawbones.zip file contains a CT scan of several sawbones phantom objects and a collection 2D fluoroscopy images of one phantom. Extract the contents:

unzip 10_sawbones.zip

This will create a new sawbones directory. The CT scan data is located in sawbones/ct and the fluoroscopy data is located in sawbones/fluoro.

Although some manual annotation is required in order to fully complete the processing steps of this example, example annotations are also available for use:

• A region of interest (ROI) used to tightly crop the original CT volume to the appropriate phantom object is available here,
• A labeling labeling of the phantom objection from the cropped CT is available here,
• 3D landmarks of the phantom object are available here,
• 2D landmarks for fluoroscopy frames 39-69 are available here and may be extracted into the working directory:

unzip sawbones_2d_landmarks.zip

• Python script for batch fluoroscopy conversion, 2D/3D registration and world coordinate frame calculation is available here.

Convert the DICOM CT Series to 3D NIFTI Files

Conversion of the DICOM CT series is accomplished using the following command:

xreg-convert-dicom-vols --name-w-conv ./sawbones/ct .

This will retain the native spacing of each DICOM series - no resampling is performed. The --name-w-conv flag will incorporate the convolution kernel used for 3D reconstruction into the output file names. This will allow us to more easily choose the appropriate series. After the command completes, the following files should exist:

• saw_bone_00_FC18.nii.gz
• saw_bone_01_FC30.nii.gz

FC18 indicates a soft-tissue reconstruction kernel and FC30 indicates a bone reconstruction kernel. The soft-tissue reconstruction file, saw_bone_00_FC18.nii.gz, will be used in this example. The reconstruction kernel type may also be inferred by reading the Series Description DICOM tag (0008|103e) using the xreg-report-dcm tool. The series with description text containing "Bone" indicates a bone tissue reconstruction, while the description text with "Body" indicates a soft-tissue reconstruction.

Cropping to the Phantom Object

When opening up the CT volume in 3D Slicer, it is easily seen that more than one sawbones object is present in the scan. This is demonstrated in the following coronal slice and volume rendering:

Manual adjustment of the windowing enhances the contrast between the objects and air:

The most superior (top) sawbones phantom object will be used in this example. Osteotomies were performed on this object prior to imaging and is the reason for the missing left acetabulum, which is actually located inferior to this object in the scan (most easily seen in the volume rendering above).

Next, annotate a ROI about the sawbones of interest so that it may be cropped into a sub-volume:

Save the ROI to cut_pelvis_roi.acsv.

A previously created ROI file is also available here. When loading this file, check the Show Options checkbox and select the ROI radial button.

The full CT volume may be cropped with the following command:

xreg-crop-vol saw_bone_00_FC18.nii.gz cut_pelvis_roi.acsv saw_bone_00_FC18_crop.nii.gz

The cropped volume may then be visualized in 3D Slicer:

The separate acetabulum fragment and the scanner bed are still visible in this cropped volume, which will most likely reduce the quality of the 2D/3D registrations. One possible solution, which is discussed in the following section, would be to segment the object we are concerned with.

Segmentation

Segmenting the larger sawbones pelvis fragment is possible with some basic operations using the Segmentation module in 3D Slicer.

First, a thresholding is performed, labeling all voxels with intensities greater than -819 HU as foreground:

This threshold value was found by manually adjusting the slider and interpreting the segmentation preview.

Next, the "islands" tools are used to interactively click on the acetabulum fragment group of labels and remove it:

The pelvis object labels are disconnected from the scanner bed labels by first performing an erosion using the "Shrink" option of the "Margin" tools with a size of 2 mm. The first screenshot below is just prior to erosion and the second screenshot is just afterwards:

Once again, the "islands" tool is used. This time, only the large pelvis group of labels is kept by using the "Keep Selected Island" tool:

Finally, correct pelvis labels that were removed during the erosion operation are restored using a dilation operation. This is accomplished using the "Grow" option of the "Margin" tools with a size of 2 mm:

That completes the labeling process. It is important to emphasize that the labeling does not need to be perfect, but rather "good enough" to enable quality 2D/3D registrations.

This segmentation is exported to a labelmap and saved to disk as saw_bone_00_FC18_crop-label.nii.gz.

A previously created labeling is available here.

For futher validation of the segmentation, a 3D mesh may be created and visualized using the following commands:

xreg-create-mesh saw_bone_00_FC18_crop-label.nii.gz pelvis_mesh.h5
xreg-show-mesh pelvis_mesh.h5 --mesh-color-bone

This should produce a view similar to the following:

Landmark Annotations

Corresponding sets of landmarks must be established in the 3D CT coordinate frame and each 2D view in order to initialize each 2D/3D registration. At least 4 corresponding landmarks must be annotated in each view to perform this initialization. The following anatomical landmarks will be used in this example:

• Left and right anterior superior iliac spine (ASIS-l, ASIS-r),
• Right anterior inferior iliac spine (AIIS-r),
• Right femoral head center (FH-r),
• Right greater sciatic notch (GSN-r),
• Inferior tip of sacrum (IS),
• Left and right inferior obturator foramen (IOF-l, IOF-r),
• Left and right medial obturator foramen (MOF-l, MOF-r),
• Left and right superior pubis symphysis (SPS-l SPS-r),
• Left and right inferior pubis symphysis (IPS-l, IPS-r).

The sawbones object has two metallic screws which connect each hemipelvis to the sacrum. The left and right screw heads (SH-l, SH-r) and screw tips (ST-l, ST-r) are used as landmarks.

Four metallic beads that were inserted into the ilium surface will also be used. Starting from most superior and in clockwise order (when viewing at an AP perspective in front of the patient), they are numbered 1, 2, 3, 4.

The large number of landmarks is necessary in order to have a sufficient number of visible landmarks in each 2D view.

The following images show example 3D annotations overlaid on a 3D mesh created from the segmentation and a volume rendering thresholded to best show the metallic objects:

Save these landmarks as ct_lands.fcsv. A previous set of 3D annotations is available here.

The 2D landmarks for each fluoroscopy image may be annotated by opening each fluoroscopy DICOM file (located in sawbones/fluoro/dicom) in 3D Slicer and picking the visible landmarks. The desired DICOM files may be loaded into 3D Slicer by either using the "Data" button or by dragging and dropping the files into the 3D Slicer window. It is not necessary to use the separate DICOM module in 3D Slicer.

This example will assume that landmarks are annotated for projections 39-69. However, users are free to annotate different sets of images and update the appropriate commands as needed.

Each projection's annotations are assumed to be saved in their own FCSV file within a sawbones_2d_landmarks subdirectory with filename equal to the DICOM file name with a .fcsv suffix. For example, sawbones_2d_landmarks/39_1.fcsv would store the landmarks for the projection image: sawbones/fluoro/39_1.

Example annotations for projections 40, 42, 43, 57 and 69 are shown in the following images:

Previous 2D annotations saved in a single .zip file are available here.

Conversion of Fluoroscopy DICOM

Prior to registering each 2D fluoroscopy view, each 2D DICOM file needs to be converted to an HDF5 projection data file that will store the pixel data, appropriate geometry metadata data and locations of the 2D landmarks.

The following command will perform this conversion for projection 57:

xreg-convert-dicom-radiograph --no-bayview-check sawbones/fluoro/dicom/57_1 57_1_pd.h5 sawbones_2d_landmarks/57_1.fcsv

This converts the 57_1 DICOM file and the 57_1.fcsv landmarks file into a single 57_1_pd.h5 file.

The --no-bayview-check indicates that the DICOM file should not be checked for any metadata that would indicate that the DICOM file was created by the Siemens CIOS Fusion C-arm from the Johns Hopkins Bayview Lab. An extrinsic camera model is populated using a previously completed calibration process when this flag is absent and the DICOM is determined to have been created at the Bayview lab. We want to pass this flag since this example uses data from the Bayview lab and our current goal is to demonstrate this extrinsic cacluation from scratch.

The geometry of the X-ray imaging scene may be viewed with the following command:

xreg-draw-xray-scene -i -l 57_1_pd.h5

Which should produce a view similar to the following:

The -i flag indicates that the fluoroscopy image should be rendered on the physical location of the X-ray detector. The -l flag indicates that 2D landmarks should be drawn as small spheres on the X-ray detector. The large green sphere in view indicates the location of the X-ray source and the line connecting the X-ray source and principal point on the X-ray detector is dneoted by the thin green line.

Projection 39 may be converted using the following command:

xreg-convert-dicom-radiograph --no-bayview-check sawbones/fluoro/dicom/39_1 39_1_pd.h5 sawbones_2d_landmarks/39_1.fcsv

Just as with projection 57, the X-ray scene may be visualized with this command:

xreg-draw-xray-scene -i -l 39_1_pd.h5

A view similar to the following should be produced:

Unlike the view for projection 57, the above image includes a rendering of the coordinate system axes, which may be toggled on and off by the user by pressing the a key.

Note that the X-axis runs parallel with increasing column direction, the Y-axis runs parallel with increasing row direction and the Z-axis is orthogonal to the detector plane and increasing from the detector towards the X-ray source. The X-ray source is located at the origin, (0,0,0) and the principal point on the detector is located in the center of the fluoroscopy image.

This conversion process, which needs to be performed on each projection, is best executed as a script in order to avoid repetitive typing and introducing possible mistakes. The python script available here will perform this conversion in addition to the 2D/3D registrations and world frame calculation described in subsequent sections.

2D/3D Registrations to Each 2D Fluoroscopy View

The next steps consists of using the previous 3D annotations and newly converted 2D projection data files to register each 2D view.

The following command will register projection 57:

xreg-hip-surg-pelvis-single-view-regi-2d-3d -s saw_bone_00_FC18_crop-label.nii.gz saw_bone_00_FC18_crop.nii.gz ct_lands.fcsv 57_1_pd.h5 regi_57.h5 regi_debug_57.h5

The -s saw_bone_00_FC18_crop-label.nii.gz indicates that the 3D segmentation of the phantom object should be used to mask the 3D CT volume. regi_57.h5 and regi_debug_57.h5 are output files storing the registration rigid transformation and debug information saved during the registration process, respectively.

The 3D geometry of the X-ray scene along with the registered phantom object may be displayed with the following command:

xreg-draw-xray-scene -i --mesh-color-bone 57_1_pd.h5 pelvis_mesh.h5 regi_57.h5

Providing the pelvis_mesh.h5 argument will include the phantom object mesh in the rendering and the regi_57.h5 argument uses the registration transform to update the pose of the phantom object.

A view similar to the following should be obtained:

The following command will create a pair of movies (edges_57.mp4 and mov_57.mp4) that help visualize the phantom object's pose during each optimization iteration during the registration:

xreg-regi2d3d-replay --flip-rows --flip-cols --proj-ds 0.5 --video-fps 10 regi_debug_57.h5 edges_57 mov_57 --boundary-edge-thresh-hu -800

Each frame of mov_57.mp4 is a digitally reconstructed radiograph (DRR) of the phantom object created using the current estimate of the phantom's pose. A video similar to the following should be produced:

Each frame of edges_57.mp4 contains a base image of the measured fluoroscopic view (e.g. pixels of DICOM 57_1) with green edge features overlaid in green. These edge features are calculated using the phantom object's current pose (boundary edges of the projected object and Canny edges computed from the DRR). The edges_57.mp4 video should appear similar to the following:

The following command will register projection 39:

xreg-hip-surg-pelvis-single-view-regi-2d-3d -s saw_bone_00_FC18_crop-label.nii.gz saw_bone_00_FC18_crop.nii.gz ct_lands.fcsv 39_1_pd.h5 regi_39.h5 regi_debug_39.h5

Just as with projection 57, the X-ray scene may be plotted using the registered pose for projection 39:

xreg-draw-xray-scene -i --mesh-color-bone 39_1_pd.h5 pelvis_mesh.h5 regi_39.h5

This should produce a view similar to the following:

The debug movies (mov_39.mp4 and edges_39.mp4) are created using the following command:

xreg-regi2d3d-replay --flip-rows --flip-cols --proj-ds 0.5 --video-fps 10 regi_debug_39.h5 edges_39 mov_39 --boundary-edge-thresh-hu -800

Videos similar to the following should be produced:

Registration needs to be performed on each projection. Just as with the DICOM conversion process, this is best executed as a script. An example python script (the same script identified previously) is available here.

Updating Extrinsic Coordinate Frame

Now that each view is registered to the phantom object, the extrinsic coordinate frames for each view may be updated using the registration transformations. This yields relative pose information between views and will eventually be used when calculating the orbital rotation axis and center of rotation.

First, all of the projection files are concatenated into a single projection data file using the following command:

xreg-proj-cat 39_1_pd.h5 40_1_pd.h5 41_1_pd.h5 42_1_pd.h5 43_1_pd.h5 44_1_pd.h5 45_1_pd.h5 46_1_pd.h5 47_1_pd.h5 48_1_pd.h5 49_1_pd.h5 50_1_pd.h5 51_1_pd.h5 52_1_pd.h5 53_1_pd.h5 54_1_pd.h5 55_1_pd.h5 56_1_pd.h5 57_1_pd.h5 58_1_pd.h5 59_1_pd.h5 60_1_pd.h5 61_1_pd.h5 62_1_pd.h5 63_1_pd.h5 64_1_pd.h5 65_1_pd.h5 66_1_pd.h5 67_1_pd.h5 68_1_pd.h5 69_1_pd.h5 all_pd.h5

The first batch of positional arguments indicate the projection data files which should be concatenated and the final positional argument indicates the output (destination) projection file. Projections are concatenated in the same order as their positions on the command line. Assuming identical directory contents as created in this example, wildcards may also be used to achieve the same result with the following command: xreg-proj-cat *_pd.h5 all_pd.h5. The full command, without wildcards, is also executed as part of the python example script. This concatenation operation does not update any relative pose information.

Relative pose information may be incorporated by updating each view's extrinsic parameters with the corresponding registration transformation. This process is also referred to as using a fiducial object to establish a world coordinate frame and may be performed using the following command:

xreg-proj-make-fiducial-world all_pd.h5 regi_39.h5 regi_40.h5 regi_41.h5 regi_42.h5 regi_43.h5 regi_44.h5 regi_45.h5 regi_46.h5 regi_47.h5 regi_48.h5 regi_49.h5 regi_50.h5 regi_51.h5 regi_52.h5 regi_53.h5 regi_54.h5 regi_55.h5 regi_56.h5 regi_57.h5 regi_58.h5 regi_59.h5 regi_60.h5 regi_61.h5 regi_62.h5 regi_63.h5 regi_64.h5 regi_65.h5 regi_66.h5 regi_67.h5 regi_68.h5 regi_69.h5

The wildcard command, xreg-proj-make-fiducial-world all_pd.h5 regi_*.h5, will not work properly since the working directory also contains regi_debug_*.h5 files. However, if the debug files are moved to a separate directory the wildcard command will execute properly. The full command, without wildcards, is also executed as part of the python example script.

The relative pose information may be verified by visually inspecting the X-ray scene with the following command:

xreg-draw-xray-scene -i --mesh-color-bone --draw-origin all_pd.h5 pelvis_mesh.h5 -

The final positional argument of - indicates the phantom object should be plotted using the identity pose. By passing --draw-origin, a yellow sphere will be drawn at the origin of the extrinsic coordinate frame.

Views similar to the following should be obtained:

Note that the coordinate axes relative to the phantom object match the (left/posterior/superior) LPS coordiantes of the CT scan exactly and the origin matches the origin established by the DICOM CT series. Although this new extrinsic frame may be useful for various multiple-view geometry tasks, it does not yet explicitly convey information regarding orbital rotation.

Executing the following command will calculate the orbital rotation axis, center of rotation and update the extrinsic transformations in all_pd.h5 accordingly:

xreg-proj-est-orbital-rot -u all_pd.h5

The -u flag indicates that the extrinsics of the projection data file should be updated. Omitting the -u flag may be useful to orbital information from the program's standard output, but updating the source projection data is not desired.

The updated X-ray scene may viewed with the following command:

xreg-draw-xray-scene -i --draw-origin all_pd.h5

Note that the phantom surface mesh is no longer included as a separate transformation would be needed to properly align it with the multiple-view geometry.

Views similar to the following should be produced:

The origin (yellow sphere) is now drawn at the center of orbital rotation and the X-axis is properly aligned so that rotation about it maps to the same rotation of the C-arm orbit. This is exactly what we set out to obtain!

Passing -d -v to xreg-proj-est-orbital-rot will enable verbose printing and a debug visualization. The verbose output will include the estimated amount of orbital rotation between consecutive views, an example from this data is shown below:

     0 -->  1 :  4.7 (deg.)
     1 -->  2 :  3.9 (deg.)
     2 -->  3 :  5.9 (deg.)
     3 -->  4 :  5.1 (deg.)
...
    27 --> 28 :  5.0 (deg.)
    28 --> 29 :  5.1 (deg.)
    29 --> 30 :  5.2 (deg.)

The debug visualization plots the X-ray source positions along with a corresponding best-fit circle. A yellow sphere is used to mark the center of the circle, which will be used as the origin (center of rotation) of the updated coordinate frame. Vectors of the updated coordinate axes are also plotted incident at the circle center. The debug visualization should appear similar to the following:

References

• Grupp, Robert B., et al. "Pose estimation of periacetabular osteotomy fragments with intraoperative X-ray navigation." IEEE Transactions on Biomedical Engineering 67.2 (2020): 441-452.