3D Reconstruction of Wrist Bones from C-Arm Fluoroscopy Using Planar Markers

概要

diagnostics
Article

3D Reconstruction of Wrist Bones from C-Arm Fluoroscopy
Using Planar Markers
Pragyan Shrestha 1 , Chun Xie 2 , Hidehiko Shishido 2 , Yuichi Yoshii 3, *
1

2
3

*

and Itaru Kitahara 2, *

Graduate School of Science and Technology, Degree Programs in Systems and Information Engineering,
Doctoral Program in Empowerment Informatics, University of Tsukuba, Tsukuba 305-8577, Japan
Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan
Department of Orthopedic Surgery, Tokyo Medical University Ibaraki Medical Center, Ami 300-0395, Japan
Correspondence: yyoshii@tokyo-med.ac.jp (Y.Y.); kitahara@ccs.tsukuba.ac.jp (I.K.)

Abstract: In orthopedic surgeries, such as osteotomy and osteosynthesis, an intraoperative 3D
reconstruction of the bone would enable surgeons to quickly assess the fracture reduction procedure
with preoperative planning. Scanners equipped with such functionality are often more expensive
than a conventional C-arm fluoroscopy device. Moreover, a C-arm fluoroscopy device is commonly
available in many orthopedic facilities. Based on the widespread use of such equipment, this paper
proposes a method to reconstruct the 3D structure of bone with a conventional C-arm fluoroscopy
device. We focus on wrist bones as the target of reconstruction in this research as this will facilitate
a flexible imaging scheme. Planar markers are attached to the target object and are tracked in the
fluoroscopic image for C-arm pose estimation. The initial calibration of the device is conducted
using a checkerboard pattern. In general, reconstruction algorithms are sensitive to geometric
calibration errors. To assess the practicality of the method for reconstruction, a simulation study
demonstrating the effect of checkerboard thickness and spherical marker size on reconstruction
quality was conducted.
Keywords: three-dimensional data; tracking; computed tomography; fluoroscopy; preoperative plan;
distal radius fracture
Citation: Shrestha, P.; Xie, C.;
Shishido, H.; Yoshii, Y.; Kitahara, I.
3D Reconstruction of Wrist Bones
from C-Arm Fluoroscopy Using
Planar Markers. Diagnostics 2023, 13,
330. https://doi.org/10.3390/
diagnostics13020330
Academic Editor: Tomoki
Nakamura
Received: 8 December 2022
Revised: 1 January 2023
Accepted: 13 January 2023
Published: 16 January 2023

Copyright: © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

1. Introduction
Radiographic imaging is a crucial technology in modern healthcare systems and
medical diagnostics. Interventional radiologists use minimally invasive procedures guided
by various modalities of medical imaging. Orthopedic surgeries, such as osteosynthesis and
osteotomies, also benefit from intra-operative X-ray imaging [1,2]. Computed tomography
(CT) scans prior to surgery are used to plan and simulate the fracture reduction process,
while intra-operative X-ray images are used as guidance for positioning implants or cutting
bones during surgery [3–7]. However, an apparent drawback with such an image-based
navigation system is the loss of depth perception due to the transmission nature of X-ray
images. To compensate for the loss of 3D information, a line of research focuses on
registering preoperative 3D models to the intraoperative X-ray images [8]. While this
serves as visual guidance during surgery, it would be beneficial to also have a registered
3D model of the current subject undergoing the surgery. This need could be fulfilled by
mobile cone beam CT (CBCT) devices which have become more common in recent years [9].
These devices have a motorized iso-centric arm that can capture motion trajectory with
high precision and reconstruct the volume inside the region of interest. However, many
C-arm devices that are currently in practical use do not have tomographic capabilities since
they are not designed for such use cases. Moreover, the initial cost of installing a mobile
CBCT device is too high for rural sites. Therefore, a simple and easy-to-setup system for
3D reconstruction using a conventional C-arm fluoroscopy device would be beneficial in
such cases.

Diagnostics 2023, 13, 330. https://doi.org/10.3390/diagnostics13020330

https://www.mdpi.com/journal/diagnostics

Diagnostics 2023, 13, 330

2 of 12

In the field of computer vision, 3D reconstruction from multi-view images of a scene
has been researched extensively [10,11]. In structure from motion [11], the problems of
camera calibration and pose estimation are solved by decomposing the fundamental matrix
obtained by point-to-point correspondences in multi-view images, which is followed by the
triangulation step for initial reconstruction. Unfortunately, these methods cannot be applied
to X-ray images since point correspondences cannot be established due to overlapping
structures. In the field of CT reconstruction, techniques such as FDK- (Feldkamp, Davis,
Kress) [12] and ART-based (Algebraic Reconstruction Technique) [13] methods are stateof-the-art techniques for tomographic reconstruction. These techniques require accurate
geometry information beforehand. Therefore, accurate geometric calibration is required for
applying these techniques in devices that are not designed for the use of these techniques.
To circumvent the problem of accurate geometric calibration, some of the related works use
a data-driven approach for reconstruction. In [14–16], the authors reconstructed the shape
of a distal femur from multiple X-ray images using statistical shape models. These methods
first require a statistical shape model, which is obtained by applying a principal component
analysis to various CT models of similar structures. The obtained model is deformed so that
its simulated projection is consistent with the acquired real projection. In [17], the author
proposed a deep learning-based approach that uses a generative adversarial network to
reconstruct CT models from biplanar X-rays. On the other hand, other works have focused
either on C-arm calibration [18,19], reconstruction [20], or both given the initial geometry
estimates [21].
Closely related to our work, Abella et al. [22] proposed a low-cost solution for tomographic reconstruction using a 3D scanner mounted on the C-arm. They identified three
major issues that occur during tomographic reconstruction using a conventional C-arm
device. First, the acquisition trajectory may deviate from the circular path and the mechanical stress causes each projection to have different acquisition parameters (i.e., source
position, detector position, detector rotations). Second, these acquisition parameters are not
repeatable for consecutive acquisitions. Third, the problem of limited angle tomography
needs to be solved. To solve these issues, they used the calibration phantom developed
by Cho et al. [23] for the geometric calibration of the device. Furthermore, the geometric
errors originating from the non-repeatability of the C-arm trajectory were corrected using
the information from the 3D scanner. An algorithm incorporating surface information was
also introduced to account for the limited-angle tomography problem.
In this paper, an image-based framework for the tomographic reconstruction of wrist
bones is proposed. A simulation study to evaluate its practicality was conducted. The
proposed method assumes rotation of the target subject instead of the C-arm, which is
possible for wrist bones. This way, the internal parameters, such as detector-to-source
distance and detector rotations, remain constant between projections. For estimating the
internal parameters of the C-arm device, a single checkerboard-based calibration was used,
as demonstrated in [24]. A simple, easy-to-use, attachable calibration plane is proposed
to estimate the C-arm pose without relying on external sensors. An open-source library
“Tomographic Iterative GPU-based Reconstruction Toolbox” (TIGRE) [25] in MATLAB was
used for reconstructing the volume. Equations were derived to convert results from the
camera system to the tomography device system used in TIGRE.
2. Materials and Methods
Figure 1 shows the flow diagram of the proposed method. The C-arm calibration
step is required before the acquisition of the subject to estimate the internal parameters
such as the distance between the source to the detector and the piercing point (i.e., the
point of intersection of the detector plane and the ray originating from the source, which
is perpendicular to the detector). Pose estimation involves tracking the planar markers
and solving for extrinsic parameters in the camera setup. The camera parameters are then
converted into CBCT geometry for TIGRE and the simultaneous iterative reconstruction
technique (SIRT) algorithm is applied for 20 iterations to reconstruct the volume.

Diagnostics 2023, 13, x FOR PEER REVIEW

Diagnostics 2023, 13, 330

3 of 13

solving for extrinsic parameters in the camera setup. The camera parameters are then con3 of 12
verted into CBCT geometry for TIGRE and the simultaneous iterative reconstruction technique (SIRT) algorithm is applied for 20 iterations to reconstruct the volume.

Figure 1. Overview of the calibration and reconstruction pipeline. The C-arm is calibrated using
Figure
1. captures
Overview
of left)
the calibration
and reconstruction
pipeline.
The C-arm
is calibrated
using
multiple
(top
of the calibration
checkerboard to
obtain intrinsic
parameters
(top-center).
multiple captures (top left) of the calibration checkerboard to obtain intrinsic parameters (top-cenA planar marker board is attached to the subject before X-ray acquisition (bottom left). Markers in
ter). A planar marker board is attached to the subject before X-ray acquisition (bottom left). Markers
the acquired multi-view images are tracked for pose estimation (bottom center), followed by the
in the acquired multi-view images are tracked for pose estimation (bottom center), followed by the
reconstructionalgorithm
algorithm(bottom
(bottomright).
right).
reconstruction

2.1. C-Arm Calibration
2.1. C-Arm Calibration
A pinhole camera model was assumed for the image formation process in a C-arm
A pinhole
modelcoordinate
was assumed
the image
formation
in a C-arm
device.
A pointcamera
in the device
was for
projected
into the
detectorprocess
plane according
to
device.
A
point
in
the
device
coordinate
was
projected
into
the
detector
plane
according
Equation (1).
to Equation (1).
x = KR[ I | − t]X
(1)

(1)
x =orientation
KR I |−t X of the X-ray camera and t is its translation
where K is the intrinsic matrix, R is the
vector,Kasisseen
from the matrix,
world coordinate.
The intrinsic
matrix
was modeled
where
the intrinsic
R is the orientation
of the
X-rayKcamera
and t isas
itsfollows.
trans
lation vector, as seen from the world coordinate.
Theintrinsic matrix K was modeled as
f x 0 cx
follows.
K =  0 f y cy 
(2)
𝑓 00 0𝑐 1
K= 0 𝑓 𝑐
(2)
where f x and f y are the focal lengths in0respective
0 1 pixel units, and c x and cy are the pixel
coordinates of the principal point (i.e., piercing point). We did not model distortions for
where 𝑓 and 𝑓 are the focal lengths in respective pixel units, and 𝑐 and 𝑐 are the
the sake of simplicity, but they can be factored in easily using the polynomial distortion
pixel coordinates of the principal point (i.e., piercing point). We did not model distortions
models. We adopted the calibration method from [24]. The design of the checkerboard
for the sake of simplicity, but they can be factored in easily using the polynomial distortion
pattern is shown in Figure 2. It was a 4 by 5 squares checkerboard containing 12 identifiable
models. We adopted the calibration method from [24]. The design of the checkerboard
feature points. The feature points were localized and identified in the multi-view projection
pattern is shown in Figure 2. It was a 4 by 5 squares checkerboard containing 12 identifiimages of the checkerboard. Zhang’s camera calibration algorithm was applied to obtain
able feature points. The feature points were localized and identified in the multi-view
the intrinsic parameters.
projection images of the checkerboard. Zhang’s camera calibration algorithm was applied
to
obtain
the intrinsic parameters.
2.2.
Pose Estimation
The planar markers were designed with five spherical markers that formed a parallelogram, as shown in Figure 3 (left). The reason for adopting such a pattern was two-fold.
First, the 10 mm gap between the vertices P1 and P3 in the vertical axis ensured that the
markers remained visible at the areas where the plane normal was close to perpendicular to
the camera view direction. Second, the introduction of marker P5 helped to identify the line
P1-P5-P2 in the projection image by using the fact that perspectivity preserves collinearity.
Marker tracking required localizing the marker points by applying the circular Hough
transform to the edge image. Then, for each permutation of the markers, the collinearity
condition was checked to identify the P1-P5-P2 line segment. After the identification of the
line segment, the markers were assigned according to their y-position in the image, i.e., the

Diagnostics 2023, 13, 330

4 of 12

Diagnostics 2023, 13, x FOR PEER REVIEW

4 of 13

marker with the lowest y-value was assigned as P1. The remaining two markers were also
assigned according to their y-value.

Figure 2. The layout of the checkerboard. The detectable image feature points are marked in red
circles.

2.2. Pose Estimation
The planar markers were designed with five spherical markers that formed a parallelogram, as shown in Figure 3 (left). The reason for adopting such a pattern was two-fold.
First, the 10 mm gap between the vertices P1 and P3 in the vertical axis ensured that the
markers remained visible at the areas where the plane normal was close to perpendicular
to the camera view direction. Second, the introduction of marker P5 helped to identify the
line P1-P5-P2 in the projection image by using the fact that perspectivity preserves collinearity.2. The layout of the checkerboard. The detectable image feature points are marked in red
Figure

Figure 2. The layout of the checkerboard. The detectable image feature points are marked in red circles.
circles.

2.2. Pose Estimation
The planar markers were designed with five spherical markers that formed a parallelogram, as shown in Figure 3 (left). The reason for adopting such a pattern was two-fold.
First, the 10 mm gap between the vertices P1 and P3 in the vertical axis ensured that the
markers remained visible at the areas where the plane normal was close to perpendicular
to the camera view direction. Second, the introduction of marker P5 helped to identify the
line P1-P5-P2 in the projection image by using the fact that perspectivity preserves collinearity.

Figure 3.
3. The
Thelayout
layoutofofthe
theplanar
planar
markers
(left)
and
coordinate
system
used
in pose
estimaFigure
markers
(left)
and
thethe
coordinate
system
used
in pose
estimation
tion (right).
(right).

We continued
withrequired
decomposing
the homography
to obtain
the the
camera
pose
Marker
tracking
localizing
the markermatrix
pointsHby
applying
circular
for every
view. Antoestimate
of image.
cameraThen,
rotation
could
as the
follows:
Hough
transform
the edge
forand
eachtranslation
permutation
of be
thederived
markers,
collinearity condition was checked to identify the P1-P5-P2 line segment. After the identifiH = [h1 h2 h3 ]
(3)
cation of the line segment, the markers were assigned according to their y-position in the
image, i.e., the marker with the lowestr1y-value
= λK−1was
h1 assigned as P1. The remaining two
(4)
markers were also assigned according to their y-value.
−1
= λK
h2
(5)
We continued with decomposingr2the
homography
matrix H to obtain the camera
pose for every view. An estimate of camera
and translation could be derived (6)
as
r3 =and
rrotation
r
1×
Figure 3. The layout of the planar markers (left)
the2coordinate system used in pose estimation
follows:
(right).
t = λK−1 h3
(7)
(3)
H = h h h1
Marker tracking required localizing
points by applying the circular
λ = the−marker
1
||for
K each
h|| permutation of the markers, the colHough transform to the edge image.rThen,
(4)
= 𝜆K h
linearity
condition
was checked
to identify
the P1-P5-P2the
line
segment.
After
the by
identifiThe
rotation
matrix
was
obtained
by
approximating
best
rotation
defined
matrix


cation
segment,
the markers
according
their y-position
in the
Further
refinement
of rthe=were
external
wasto
achieved
by minimizing
r1 rof
r3 .line
(5)
2 the
𝜆K assigned
h parameters
image,
i.e.,
the
marker
with
the
lowest
y-value
was
assigned
as
P1.
The
remaining
two
the reprojection error. The obtained rotation matrix and translation vectors were defined
for the coordinate
system shown
in Figure
3 (right).
markers
were also assigned
according
to their
y-value.
We continued with decomposing the homography matrix H to obtain the camera
2.3. Reconstruction
pose
for every view. An estimate of camera rotation and translation could be derived as
The parameters estimated in the camera system needed to be converted to those
follows:
used in TIGRE. The parameters involved are summarized in Table 1. The conversion was
(3)
H in
= Figure
h h h4. During pose estimation, projection images
derived from the geometry observed
with mean reprojection errors larger than a given threshold (outliner data) were excluded
(4)
r = 𝜆K h

r = 𝜆K h

(5)

Diagnostics 2023, 13, x FOR PEER REVIEW
Diagnostics 2023, 13, 330

5 of 13
5 of 12

(6)
r = r image
r quality. For reconstruction of the volume,
because they degrade the reconstruction
the simultaneous iterative reconstruction technique (SIRT) in TIGRE was used with 20
(7)
t = 𝜆K h
iterations for all experiments.


T
R = x y z
(8)
1
𝜆=
t = −h‖
RT ∗t
(9)
‖K
The rotation matrix was obtained by
approximating
the best rotation defined by mad=
(10)
(t·z)z − t
trix r r r . Further refinement of the external parameters was achieved by minimizDSO = ||t − d ||
(11)
ing the reprojection error. The obtained rotation matrix and translation vectors were deDSD
= f x3 ∗(right).
α
(12)
fined for the coordinate system shown in
Figure


Rsource = −z x y
(13)
2.3. Reconstruction
where
is the camera
orientation
is the translation
vector,
d is the
offset
TheRparameters
estimated
in thematrix,
camerat system
needed to be
converted
to image
those used
vector,
α
is
the
size
of
the
detector
per
pixel
and
R
is
the
rotation
matrix
for
computing
source
in TIGRE. The parameters involved are summarized in Table 1. The conversion was deEulerfrom
angles.
The detector
rotation
offsets
are assumed
to be zeroprojection
in this conversion
rived
the geometry
observed
in and
Figure
4. During
pose estimation,
images
equations,
even
if
is
present
in
actual
system.
However,
it
does
not
not
affectexcluded
the reconwith mean reprojection errors larger than a given threshold (outliner data) were
struction
algorithm
the parameters
obtained
the above calculation
encodes
because
they
degradebecause
the reconstruction
image
quality.from
For reconstruction
of the volume,
such information (i.e., equivalent to redefining the coordinate system such that detector
the simultaneous iterative reconstruction technique (SIRT) in TIGRE was used with 20
rotation and offsets became zero in each view independently).
iterations for all experiments.
Table1.1.List
Listofofparameter
parameternames
namesused
usedininTIGRE
TIGREalong
alongwith
withtheir
theirdescriptions.
descriptions.
Table
Parameter Name
Parameter Name

DSD
DSD
DSO
DSO
offOrigin offOrigin
offDetector offDetector
angles
angles
nVoxel
nVoxel
dVoxel
dVoxel
sVoxel
nDetector sVoxel
sDetector nDetector
sDetector

Description
Description
Distance
between
the X-ray
source
and detector
planeplane
Distance
between
the X-ray
source
and detector
Distance
between
the X-ray
source
and world
Distance
between
the X-ray
source
and world
originorigin
applied
the volume
in world
OffsetOffset
applied
to thetovolume
in world
originorigin
Offset
applied
to
the
detector
Offset applied to the detector
The angle
of X-ray
source
in ZYZ
convention
The angle
of X-ray
source
from from
worldworld
originorigin
in ZYZ
convention
Number of voxels in the volume
Number of voxels in the volume
Size of one voxel in the world coordinate
Size of one voxel in the world coordinate
Total size occupied by the volume in the world coordinate
Total size occupied
by the
volume
the
world coordinate
Number
of pixels
inin
the
detector
panel
Number
of
pixels
in
the
detector
panel
The total size of the detector panel in the
world coordinate
The total size of the detector panel in the world coordinate

Figure 4. Relation between TIGRE geometry parameters and camera system setup.
Figure 4. Relation between TIGRE geometry parameters and camera system setup.

3. Experiments and Results
3.1. Simulation Setup

(8)
R = x y z
A simulation environment was set up to investigate the effect of different calibration
board thicknesses and spherical marker
(9)
∗ t on image feature points as well as image
t = −Rsizes
quality. TIGRE for MATLAB was used for simulating the X-ray generation, as well as

Diagnostics 2023, 13, 330

6 of 12

reconstruction algorithms. We used the internal parameters of the SIEMENS CIOS Select
tool for simulating the C-arm. The common configurations used for simulated X-ray
generation are shown in Table 2.
Table 2. Values of some common parameters used in the simulation.
Parameter Name

Value

DSD
DSO
offOrigin
offDetector
Angles
nDetector
sDetector

780 mm
390 mm
(0, 0, 0)
(0, 0)
180 samples between 0 and 180 degrees
(1024 px, 1024 px)
210 mm
Calibration Board
nVoxel [vx]
(500, 500, 500)
dVoxel [mm/vx]
(0.25, 0.25, 0.25)
sVoxel [mm]
(100, 100, 100)
CT Volume for Simulated X-ray
nVoxel [vx]
(604, 604, 644)
dVoxel [mm/vx]
(0.25, 0.25, 0.25)
sVoxel [mm]
(151, 151, 161)
Reconstruction Volume
nVoxel [vx]
(300, 300, 300)
dVoxel [mm/vx]
(0.5, 0.5, 0.5)
sVoxel [mm]
(150, 150, 150)

A 3D model of the checkerboard with a resolution of 0.25 mm/vx and a size of 100 mm
7 of 13
in all three dimensions was built with Blender and converted into voxel data with the
interior filled with a constant value. Figure 5 shows the 3D model along with an example
of an X-ray image in a particular pose. A CT scan of the wrist phantom with a resolution of
of
0.25
mm/vxininall
allthree
threedimensions
dimensionsand
andaasize
size of
of 151
151 mm
mm for width and depth,
0.25
mm/vx
depth, measurmeasuring
ing
161
mm
in
height,
was
used
for
simulating
X-rays.
Figure
6
shows
the
volume
render161 mm in height, was used for simulating X-rays. Figure 6 shows the
rendering
ing
of the
model
alongwith
withan
anexample
example X-ray
X-ray image
The
planar
of the
CTCT
model
along
image at
ataa0-degree
0-degreerotation.
rotation.
The
planar
markers
ofof
thethe
phantom
byby
editing
thethe
voxel
values.
markerswere
wereattached
attachedtotothe
theposterior
posteriorregion
region
phantom
editing
voxel
values.
During
outliers
(i.e.,
projection
images
Duringpose
poseestimation,
estimation,the
thethreshold
thresholdvalue
valuefor
forrejecting
rejecting
outliers
(i.e.,
projection
images
and
geometries)
based
on
the
reprojection
error
was
set
to
50
pixels.
This
threshold
waswas
and geometries) based on the reprojection error was set to 50 pixels. This threshold
selected
manually
by
inspecting
the
distribution
of
the
reprojection
error
so
that
at
least
selected manually by inspecting the distribution of the reprojection error so that at least
80%
sparse
view-related
artifacts
80% of
of the
the original
originalimages
imageswere
wereretained,
retained,while
whilepreventing
preventing
sparse
view-related
artifacts
from
dominating
the
reconstruction
error.
Evaluation
of
the
reconstructed
volume
waswas
from dominating the reconstruction error. Evaluation of the reconstructed volume
performed
using
the
result
of
the
SIRT
reconstructed
volume
with
ground
truth
acquisiperformed using the result of the SIRT reconstructed volume with ground truth acquisition
tion
geometry
for calibration
board
and
marker
size.Figure
Figure77shows
showsthe
the volume
volume rendering
geometry
for calibration
board
and
marker
size.
rendering of
of
SIRT
reconstructed
volume
with
groundtruth
truthgeometry
geometrywith
witha a1 1mm
mmboard
boardsize
sizeand
thethe
SIRT
reconstructed
volume
with
ground
and
2 mm
marker
2 mm
marker
size.size.

Diagnostics 2023, 13, x FOR PEER REVIEW

Figure5.5.AA3D
3Dview
viewofofcalibration
calibrationboard
boardininBlender
Blender
(left).
projection
board
(right).
Figure
(left).
AA
projection
of of
thethe
board
(right).

Diagnostics 2023, 13, 330

7 of 12

Figure 5. A 3D view of calibration board in Blender (left). A projection of the board (right).

Diagnostics 2023, 13, x FOR PEER REVIEW

8 of 13

Figure 6. Volume
thethe
CTCT
scan
of wrist
phantom
(left).(left).
A projection
of the volume
(right).
Figure
Volumerendering
renderingofof
scan
of wrist
phantom
A projection
of the volume
(right).

Figure7.7.Volume
Volumerendering
rendering
reconstructed
volume
using
ground
truth
projection
geFigure
ofof
thethe
reconstructed
volume
using
ground
truth
projection
and and
geomometries
(left).
(topmeta
left: carpals,
meta carpals,
topcarpals,
right: carpals,
bottom
right: bottom
carpals,
etries
(left).
FourFour
fixedfixed
slicesslices
(top left:
top right:
bottom right:
carpals,
left:
radius
and
ulna)and
of the
volume
bottom
left:
radius
ulna)
of the(right).
volume (right).

3.2.Effect
EffectofofBoard
BoardThicknesses
Thicknessesand
andMarker
MarkerSizes
Sizeson
onReconstructed
ReconstructedImage
ImageQuality
Quality
3.2.
Theproposed
proposedmethod
methodfor
forcalibration
calibrationuses
usesaacheckerboard
checkerboardpattern
patternto
toidentify
identifyfeature
feature
The
points.
However,
in
practice,
radio-opaque
lead
plates
used
for
black
squares
are
quite
thick.
points. However, in practice, radio-opaque lead plates used for black squares are quite
This
will
create
ambiguities
in
the
detected
feature
points.
An
example
of
such
ambiguities
thick. This will create ambiguities in the detected feature points. An example of such amis shown in Figure 8. For checkerboard patterns, the world feature points were in the corner
biguities is shown in Figure 8. For checkerboard patterns, the world feature points were
of the square, halfway through the board’s thickness. For spherical planar markers, the
in the corner of the square, halfway through the board’s thickness. For spherical planar
world feature points were considered to be the centers of each sphere. However, due to
markers, the world feature points were considered to be the centers of each sphere. Howperspectivity, its projection in the image plane deviated from the detected center of the
ever, due to perspectivity, its projection in the image plane deviated from the detected
ellipse. Due to these ambiguities in the image points, the resulting reconstruction quality
center of the ellipse. Due to these ambiguities in the image points, the resulting reconmay degrade depending on the chosen sizes of these materials. Therefore, we evaluated
struction quality may degrade depending on the chosen sizes of these materials. Therethe proposed method with four calibration board sizes ranging from 1 mm to 5 mm in
fore, we evaluated the proposed method with four calibration board sizes ranging from 1
thickness and five spherical marker sizes ranging from 2 mm to 6 mm in diameter, which
mm to 5 mm in thickness and five spherical marker sizes ranging from 2 mm to 6 mm in
are the commonly available sizes. A 1 mm spherical marker was eliminated due to the
diameter, which are the commonly available sizes. A 1 mm spherical marker was elimiresolution of the X-ray simulating volume.
natedTable
due to3 the
resolution
of the X-ray
simulating
shows
the structural
similarity
index volume.
measure (SSIM) compared to ground
truth reconstruction for each case. The best case of 0.9441 was achieved with 3 mm board
thickness and a marker size of 2 mm. The worst case of 0.8982 was obtained with 5 mm
board thickness and a marker size of 6 mm. Detection errors for planar markers and the
feature points in the calibration board are shown in Figure 9. Although a total of 180
projection images were captured around a 180 degrees angle of rotation, projection images
where the planar markers were projected in a lateral view were prone to pose estimation
errors, as well as wrong marker identification. These outlier images were rejected using

Diagnostics 2023, 13, 330

The proposed method for calibration uses a checkerboard pattern to identify feature
points. However, in practice, radio-opaque lead plates used for black squares are quite
thick. This will create ambiguities in the detected feature points. An example of such ambiguities is shown in Figure 8. For checkerboard patterns, the world feature points were
8 of 12
in the corner of the square, halfway through the board’s thickness. For spherical planar
markers, the world feature points were considered to be the centers of each sphere. However, due to perspectivity, its projection in the image plane deviated from the detected
center
the ellipse.
to these error.
ambiguities
image
points, the
resulting images
reconsimpleoffiltering
of theDue
reprojection
Figure in
10 the
shows
the number
of projection
struction
quality
may
degrade
depending
on
the
chosen
sizes
of
these
materials.
Therethat were used (i.e., after the filtering) for reconstruction for each marker size and board
fore,
we evaluated
the proposed
method
with fourand
calibration
sizes
from 1
thickness.
An example
of the volume
rendering
volume board
slices of
theranging
reconstruction
mm
thickness
and five spherical
sizes ranging
mm to
6 mm in
withtoa 53 mm
mm in
thickness
calibration
board andmarker
2 mm planar
markerfrom
(best2case)
is shown
in
diameter,
are the
commonly
sizes. A calibration
1 mm spherical
was planar
elimiFigure 11.which
A similar
example
with aavailable
5 mm thickness
boardmarker
and 6 mm
nated
due
to thecase)
resolution
of the
X-ray simulating
volume.
marker
(worst
is shown
in Figure
12.

Diagnostics 2023, 13, x FOR PEER REVIEW

9 of 13

where the planar markers were projected in a lateral view were prone to pose estimation
errors, as well as wrong marker identification. These outlier images were rejected using
simple filtering of the reprojection error. Figure 10 shows the number of projection images
that were used (i.e., after the filtering) for reconstruction for each marker size and board
thickness. An example of the volume rendering and volume slices of the reconstruction
Figure8.8.Close-up
Close-upimage
imageof
ofthe
theprojection
projectionof
ofthe
thecalibration
calibrationcheckerboard
checkerboard (left).Close-up
Close-upimage
imageof
of
Figure
with a 3 mm
thickness calibration
board and 2 mm planar(left).
marker (best case)
is shown in
theprojection
projectionof
ofspherical
sphericalmarker
marker(right).
(right).The
Thered
redcross
crossisisthe
theground
groundtruth
truthimage
imagepoint;
point;the
theblue
blue
the
Figure 11. A similar example with a 5 mm thickness calibration board and 6 mm planar
circle
image
point.
circleisisthe
thedetected
detected
image
point.
marker
(worst
case) is shown in Figure 12.
TableTable
3. SSIM
of the reconstructed
volume
with ground
truth
reconstructed
volume
for each configuration
3 shows
structural
similarity
index
measure
compared
to ground
Table the
3. SSIM
of the reconstructed
volume
with
ground(SSIM)
truth reconstructed
volume
for each conpair.
For
each
marker
size,
the
board
thickness
that
resulted
in
the
best
SSIM
are
marked
bold.
figuration
eachThe
marker
size,
the of
board
thickness
resultedwith
in thein
SSIM
are marked
truth reconstruction
for pair.
eachFor
case.
best
case
0.9441
was that
achieved
3best
mm
board
bold.
thickness
and ainmarker
size
of 2 mm.
The worst
case of 0.8982
5 mm
Board Thickness/Marker
Size
2 mm
3 mm
4 mm was obtained
5 mm with
6 mm
board thickness and
a marker
size of 6 Size
mm. Detection
for planar
and the6 mm
2 mm errors
3 mm
4 mmmarkers
5 mm
Board
Thickness/Marker
1 mm
0.9435
0.9205
0.9235
0.9120
0.9114
1 mm
0.9435
0.9205
0.9235
0.9114
feature points2in
the calibration
board
are shown
in Figure0.9236
9. Although
a total 0.9120
of 180
promm
0.9437
0.9206
0.9122
0.9111
2
mm
0.9437
0.9206
0.9236
0.9122
0.9111
jection images3 were
captured
around
a
180
degrees
angle
of
rotation,
projection
images
mm
0.9441
0.9206
0.9237
0.9125
0.9107
3 mm
0.9426
4 mm
0.9352
5 mm

4 mm
5 mm

0.9441

0.9225
0.9426
0.9218
0.9352

Marker detection error

Error [px]

Error [px]

1.5
1
0.5
0
3mm

4mm

5mm

0.9237
0.9116
0.9236
0.9107
0.9221

0.9125
0.9107
0.9043 0.9043
0.9116
0.8982 0.8982
0.9107

Calibration interest point detection error

2

2mm

0.9206

0.9236
0.9225
0.9221
0.9218

3
2.5
2
1.5
1
0.5
0

6mm

1mm

Marker Size

2mm

3mm

4mm

5mm

Board Thickness

Figure 9. MarkerFigure
detection
errordetection
according
to according
marker size
(left). Calibration
interest point
detection
9. Marker
error
to marker
size (left). Calibration
interest
point detection
error
according
to board
thickness (right).
error according to
board
thickness
(right).

Number of images

Effect of marker size on inlier projection images
175
170
165
160
155
150
145
140

Board
1mm
2mm
3mm
4mm
2mm

3mm

4mm

Marker Size

5mm

6mm

5mm

0

0
2mm

3mm

4mm

5mm

6mm

1mm

Marker Size
Diagnostics 2023, 13, 330

2mm

3mm

4mm

5mm

Board Thickness
9 of 12

Figure 9. Marker detection error according to marker size (left). Calibration interest point detection
error according to board thickness (right).

Number of images

Effect of marker size on inlier projection images
175
170
165
160
155
150
145
140

Board
1mm
2mm
3mm
4mm
2mm

3mm

4mm

5mm

6mm

5mm

Marker Size
Figure 10.
10. The
The number
number of
of inlier
inlier projection
projection images
images as
as the
the size
size of
ofthe
theplanar
planarmarker
marker
and
calibration
Diagnostics 2023, 13, x FOR PEER REVIEW
10 and
of
13calibration
Figure
Diagnostics 2023, 13, x FOR PEER REVIEW
10 of 13
boardthickness
thicknesschanges.
changes.
board

Figure11.11.
Volume
rendering
ofreconstructed
the reconstructed
volume
3 mm calibration
Figure
Volume
rendering
of the
volume
using ausing
3 mm a
calibration
board andboard
2 mm and 2 mm
Figure 11.
Volume
rendering
theleft:
reconstructed
volume
using
a 3 mm
calibration
and
2 mm
marker
(left).
FourFour
fixed
slicesof
(top
metacarpals,
top right:
carpals,
bottom
right:board
carpals,
bottom
marker
(left).
fixed
slices
(top
left: metacarpals,
top right:
carpals,
bottom
right: carpals,
marker
(left).
Four
fixed
slices
(top
left:
metacarpals,
top
right:
carpals,
bottom
right:
carpals,
bottom
left: radius and ulna) of the volume (right).
bottom
left:
radius
and
ulna)
of
the
volume
(right).
left: radius and ulna) of the volume (right).

Figure 12. Volume rendering of the reconstructed volume using a 5 mm calibration board and 6 mm
Figure12.
12.Volume
Volume
rendering
ofleft:
themetacarpals,
reconstructed
acalibration
5 mm
calibration
board
Figure
rendering
the
reconstructed
volume
using
ausing
5 mm
board
and
6 mm and 6 mm
marker
(left).
Four fixed
slicesof
(top
top volume
right:
carpals,
bottom
right:
carpals,
bottom
marker
(left).
fixed
slices
(top left:
metacarpals,
top right: carpals,
bottom
right: carpals,
bottom
marker
(left).
Four
slices
(top
left: metacarpals,
top right:
carpals,
bottom
right: carpals,
left:
radius
andFour
ulna)
of fixed
the
volume
(right).
left: radius and ulna) of the volume (right).

bottom left: radius and ulna) of the volume (right).
4. Discussion
4. Discussion
From Table 3, we can observe that the reconstructed image quality decreased as the
From Table 3, we can observe that the reconstructed image quality decreased as the
marker size increased. Although the image quality degraded when comparing reconstrucmarker size increased. Although the image quality degraded when comparing reconstruction results of a sample with 1 mm board thickness with one with a 5 mm board thickness,
tion results of a sample with 1 mm board thickness with one with a 5 mm board thickness,
the effect was not as significant as those obtained when changing the marker size. This
the effect was not as significant as those obtained when changing the marker size. This
suggests that the geometrical calibration of internal parameters obtained with checkersuggests that the geometrical calibration of internal parameters obtained with checkerboard calibration can be used for image reconstruction. As opposed to the calibration

Diagnostics 2023, 13, 330

10 of 12

4. Discussion
From Table 3, we can observe that the reconstructed image quality decreased as the
marker size increased. Although the image quality degraded when comparing reconstruction results of a sample with 1 mm board thickness with one with a 5 mm board thickness,
the effect was not as significant as those obtained when changing the marker size. This
suggests that the geometrical calibration of internal parameters obtained with checkerboard
calibration can be used for image reconstruction. As opposed to the calibration phantom
in [23], the checkerboard pattern was easier to prepare, and off-the-shelf software is readily
available for this kind of planar calibration.
On the other hand, the steady drop in the number of used images shown in Figure 10
was due to the larger detection error in larger marker sizes amplifying the reprojection
error, as can be observed from Figure 11, to above the threshold of 50 pixels for outlier
images. The effect of this variation in the number of images used can be observed in
the reconstructed image in Figures 11 and 12. The slice images in both cases had streak
artifacts that were a result of insufficient projection angles. The slice image in Figure 12,
however, contained metal artifacts resulting from a larger marker size as well. This suggests
that the size of the spherical marker in our proposed method was crucial in obtaining a good
reconstruction quality. From Table 3, it can be seen that approximately a 2% increase in SSIM
could be confirmed when compared to using a 3 mm marker size. Thus, it can be concluded
that the proposed method leads to better results when a 2 mm marker size is selected.
Since the proposed method cannot process projection images in which the projections
of the markers are close to being co-linear, there is an inherent limitation to the image
quality which can be seen in Figure 11, which shows the reconstructed image with the
configuration that led to the best SSIM. Additionally, our method assumes constant intrinsic
parameters throughout the acquisition. This limits its applicability to structures that cannot
be rotated easily such as pelvis.
Due to the recent advances in preoperative 3D programs for various procedures in
orthopedic surgery [3,26,27], it is necessary to establish a method for comparing 3D models
of a preoperative plan to 3D models during surgery. The method developed here has
the potential to facilitate the comparison of intraoperative 3D models with 3D images of
preoperative planning. In previous studies, the reduction accuracy of 3D preoperative
planning was moderate. This is because it is difficult to visualize the 3D model of the
reduction shape during surgery. Thus, this method will be useful in the clinical field with
comparisons to a preoperative 3D model. The advantage of this method is that there is no
need to install a new system in the operating room. Furthermore, the proposed method
does not require the installation of external devices, which makes it easier to adapt to many
existing systems.
In conclusion, this paper proposed an image-based reconstruction framework using
checkerboard calibration and rotation of the target subject with attached plane markers.
It was found that the image quality of the reconstruction depended heavily on the size
of the planar markers used. Smaller markers could be detected with low image point
error, resulting in an increased number of reliable projection images and geometries that
were used for reconstruction. The effect of calibration board thickness is not as significant as that of the planar marker size, but it should be kept below 4 mm for optimal
results. Further improvements can be achieved by installing another set of planar markers
attached perpendicular to the current one and by solving the camera pose estimation
problem with Perspective-n-Point [28], so that lateral views can also be incorporated into
the reconstruction pipeline. This is possible due to the cylindrical nature of wrist bones.
Author Contributions: P.S.: research design, writing the program, and writing of the manuscript,
C.X.: interpretation of results, and writing of the manuscript, H.S.: interpretation of results, and writing of the manuscript, I.K.: research design, interpretation of results, and writing of the manuscript,
Y.Y.: research design, interpretation of results, and writing of the manuscript. All authors have read
and agreed to the published version of the manuscript.

Diagnostics 2023, 13, 330

11 of 12

Funding: This research was supported by grants from the JSPS KAKENHI (Grant Number 19K09582),
AMED (A324TS), and the National Mutual Insurance Federation of Agricultural Cooperatives. These funders were not involved in data collection, data analysis, or the preparation or editing of the manuscript.
Institutional Review Board Statement: This research protocol was approved by our Institutional
Review Board (T2019-0178, 25 January 2020).
Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets analyzed during the present research are available from
the corresponding author upon reasonable request.
Conflicts of Interest: No benefits in any form have been received or will be received from a commercial party related directly or indirectly to the subject of this article.

References
1.
2.
3.
4.
5.

6.
7.

8.
9.
10.
11.
12.
13.
14.
15.
16.
17.

18.

19.
20.
21.

Keil, H.; Trapp, O. Fluoroscopic Imaging: New Advances. Injury 2022, 53 (Suppl. 3), S8–S15. [CrossRef]
Merloz, P.; Troccaz, J.; Vouaillat, H.; Vasile, C.; Tonetti, J.; Eid, A.; Plaweski, S. Fluoroscopy-Based Navigation System in Spine
Surgery. Proc. Inst. Mech. Eng. H 2007, 221, 813–820. [CrossRef] [PubMed]
Yoshii, Y.; Kusakabe, T.; Akita, K.; Tung, W.L.; Ishii, T. Reproducibility of Three Dimensional Digital Preoperative Planning for the
Osteosynthesis of Distal Radius Fractures. J. Orthop. Res. 2017, 35, 2646–2651. [CrossRef]
Harris, E.; Atthakomol, P.; Korteweg, J.; Li, G.; Hosseini, A.; Kachooei, A.; Chen, N.C. Non-Invasive 3D Fluoroscopic/CT Imaging
of the Wrist. Orthop. J. Harv. Med. Sch. 2021, 22, 40–48.
Wylie, J.D.; Ross, J.A.; Erickson, J.A.; Anderson, M.B.; Peters, C.L. Operative Fluoroscopic Correction Is Reliable and Correlates
with Postoperative Radiographic Correction in Periacetabular Osteotomy. Clin. Orthop. Relat. Res. 2017, 475, 1100–1106.
[CrossRef] [PubMed]
Reichert, J.C.; Hofer, A.; Matziolis, G.; Wassilew, G.I. Intraoperative Fluoroscopy Allows the Reliable Assessment of Deformity
Correction during Periacetabular Osteotomy. J. Clin. Med. Res. 2022, 11, 4817. [CrossRef]
Selles, C.A.; Beerekamp, M.S.H.; Leenhouts, P.A.; Segers, M.J.M.; Goslings, J.C.; Schep, N.W.L.; EF3X Study Group. The
Value of Intraoperative 3-Dimensional Fluoroscopy in the Treatment of Distal Radius Fractures: A Randomized Clinical Trial.
J. Hand Surg. Am. 2020, 45, 189–195. [CrossRef]
Markelj, P.; Tomaževiˇc, D.; Likar, B.; Pernuš, F. A Review of 3D/2D Registration Methods for Image-Guided Interventions.
Med. Image Anal. 2012, 16, 642–661. [CrossRef]
Posadzy, M.; Desimpel, J.; Vanhoenacker, F. Cone Beam CT of the Musculoskeletal System: Clinical Applications. Insights Imaging
2018, 9, 35–45. [CrossRef]
Schonberger, J.L.; Frahm, J.-M. Structure-from-Motion Revisited. In Proceedings of the IEEE Conference on Computer Vision and
Pattern Recognition, Las Vegas, NV, USA, 27–30 June 2016; pp. 4104–4113.
Cheung, G.K.; Baker, S.; Kanade, T. Shape-From-Silhouette Across Time Part I: Theory and Algorithms. Int. J. Comput. Vis. 2005,
62, 221–247. [CrossRef]
Feldkamp, L.A.; Davis, L.C.; Kress, J.W. Practical Cone-Beam Algorithm. J. Opt. Soc. Am. A JOSAA 1984, 1, 612–619. [CrossRef]
Andersen, A.H.; Kak, A.C. Simultaneous Algebraic Reconstruction Technique (SART): A Superior Implementation of the ART
Algorithm. Ultrason. Imaging 1984, 6, 81–94. [CrossRef]
Baka, N.; Kaptein, B.L.; de Bruijne, M.; van Walsum, T.; Giphart, J.E.; Niessen, W.J.; Lelieveldt, B.P.F. 2D-3D Shape Reconstruction
of the Distal Femur from Stereo X-ray Imaging Using Statistical Shape Models. Med. Image Anal. 2011, 15, 840–850. [CrossRef]
Zheng, G. Personalized X-ray Reconstruction of the Proximal Femur via Intensity-Based Non-Rigid 2D-3D Registration.
Med. Image Comput. Comput. Assist. Interv. 2011, 14, 598–606. [PubMed]
Yu, W.; Chu, C.; Tannast, M.; Zheng, G. Fully Automatic Reconstruction of Personalized 3D Volumes of the Proximal Femur from
2D X-ray Images. Int. J. Comput. Assist. Radiol. Surg. 2016, 11, 1673–1685. [CrossRef]
Ying, X.; Guo, H.; Ma, K.; Wu, J.; Weng, Z.; Zheng, Y. X2CT-GAN: Reconstructing CT from Biplanar X-rays with Generative
Adversarial Networks. In Proceedings of the 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR),
Long Beach, CA, USA, 15–20 June 2019; pp. 10611–10620. [CrossRef]
Vachon, É.; Miró, J.; Duong, L. Online C-Arm Calibration Using a Marked Guide Wire for 3D Reconstruction of Pulmonary
Arteries. In Proceedings of the Medical Imaging 2017: Image-Guided Procedures, Robotic Interventions, and Modeling, Orlando, FL, USA,
11–16 February 2017; Webster, R.J., Fei, B., Eds.; SPIE: Bellingham, WA, USA, 2017.
Albiol, F.; Corbi, A.; Albiol, A. Geometrical Calibration of X-ray Imaging With RGB Cameras for 3D Reconstruction. IEEE Trans.
Med. Imaging 2016, 35, 1952–1961. [CrossRef]
Kalmykova, M.; Poyda, A.; Ilyin, V. An Approach to Point-to-Point Reconstruction of 3D Structure of Coronary Arteries from 2D
X-ray Angiography, Based on Epipolar Constraints. Procedia Comput. Sci. 2018, 136, 380–389. [CrossRef]
Pande, K.; Donatelli, J.J.; Parkinson, D.Y.; Yan, H.; Sethian, J.A. Joint Iterative Reconstruction and 3D Rigid Alignment for X-ray
Tomography. Opt. Express 2022, 30, 8898–8916. [CrossRef]

Diagnostics 2023, 13, 330

22.
23.
24.
25.
26.
27.
28.

12 of 12

Abella, M.; de Molina, C.; Ballesteros, N.; García-Santos, A.; Martínez, Á.; García, I.; Desco, M. Enabling Tomography with
Low-Cost C-Arm Systems. PLoS ONE 2018, 13, e0203817. [CrossRef] [PubMed]
Cho, Y.; Moseley, D.J.; Siewerdsen, J.H.; Jaffray, D.A. Accurate Technique for Complete Geometric Calibration of Cone-Beam
Computed Tomography Systems. Med. Phys. 2005, 32, 968–983. [CrossRef] [PubMed]
Zhang, Z. A Flexible New Technique for Camera Calibration. IEEE Trans. Pattern Anal. Mach. Intell. 2000, 22, 1330–1334.
[CrossRef]
Biguri, A.; Dosanjh, M.; Hancock, S.; Soleimani, M. TIGRE: A MATLAB-GPU Toolbox for CBCT Image Reconstruction. Biomed.
Phys. Eng. Express 2016, 2, 055010. [CrossRef]
Yoshii, Y.; Teramura, S.; Oyama, K.; Ogawa, T.; Hara, Y.; Ishii, T. Development of three-dimensional preoperative planning system
for the osteosynthesis of distal humerus fractures. BioMed Eng. Online 2020, 19, 56. [CrossRef]
Yoshii, Y.; Ogawa, T.; Hara, Y.; Totoki, Y.; Ishii, T. An image fusion system for corrective osteotomy of distal radius malunion.
BioMed Eng. Online 2021, 20, 66. [CrossRef]
Lu, X.X. A Review of Solutions for Perspective-n-Point Problem in Camera Pose Estimation. J. Phys. Conf. Ser. 2018, 1087, 052009.
[CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.