# Generating Multiple Diverse Hypotheses for Human 3D Pose Consistent with 2D
Joint Detections

Ehsan Jahangiri, Alan L. Yuille
Johns Hopkins University, Baltimore, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a method to generate multiple diverse and valid human pose
hypotheses in 3D all consistent with the 2D detection of joints in a monocular
RGB image. We use a novel generative model uniform (unbiased) in the space of
anatomically plausible 3D poses. Our model is compositional (produces a pose
by combining parts) and since it is restricted only by anatomical constraints
it can generalize to every plausible human 3D pose. Removing the model bias
intrinsically helps to generate more diverse 3D pose hypotheses. We argue that
generating multiple pose hypotheses is more reasonable than generating only a
single 3D pose based on the 2D joint detection given the depth ambiguity and
the uncertainty due to occlusion and imperfect 2D joint detection. We hope
that the idea of generating multiple consistent pose hypotheses can give rise
to a new line of future work that has not received much attention in the
literature. We used the Human3.6M dataset for empirical evaluation.

##  1 Introduction

Figure 1: The input monocular image is first passed through a CNN-based 2D
joint detector which outputs a set of heatmaps for soft localization of 2D
joints. The 2D detections are then passed to a 2D-to-3D pose estimator to
obtain an estimate of the 3D torso and the projection matrix. Using the
estimated 3D torso, the projection matrix, and the output of the 2D detector
we generate multiple diverse 3D pose hypotheses consistent with the output of
2D joint detector.

Estimating the 3D pose configurations of complex articulated objects such as
humans from monocular RGB images is a challenging problem. There are multiple
factors contributing to the difficulty of this critical problem in computer
vision: (1) multiple 3D poses can have similar 2D projections. This renders 3D
human pose reconstruction from its projected 2D joints an ill-posed problem;
(2) the human motion and pose space is highly nonlinear which makes pose
modeling difficult; (3) detecting precise location of 2D joints is challenging
due to the variation in pose and appearance, occlusion, and cluttered
background. Also, minor errors in the detection of 2D joints can have a large
effect on the reconstructed 3D pose. These factors favor a 3D pose estimation
system that takes into account the uncertainties and suggests multiple
possible 3D poses constrained only by reliable evidence. Often in the image,
there exist much more detailed information about the 3D pose of a human than
the 2D location of the joints (such as contextual information and difference
in shading/texture due to depth disparity). Hence, most of the possible 3D
poses consistent with the 2D joint locations can be rejected based on more
detailed image information (_e.g_. in an analysis-by-synthesis framework or by
investigating the image with some mid-level queries such as “Is the left hand
in front of torso?”) or by physical laws (_e.g_. gravity). We can also imagine
scenarios where the image does not contain enough information to rule out or
favor one 3D pose configuration over another especially in the presence of
occlusion. In this paper, we focus on generating multiple plausible and
diverse 3D pose hypotheses which while satisfying humans anatomical
constraints are still consistent with the output of the 2D joint detector.
Figure 1 illustrates an overview of our approach.

The space of valid human poses is a non-convex complicated space constrained
by the anatomical and anthropomorphic limits. A bone never bends beyond
certain angles with respect to its parent bone in the kinematic chain and its
normalized length, with respect to other bones, cannot be much shorter/longer
than standard values. This inspired Akhter and Black [1] to build a motion
capture dataset composed of 3D poses of flexible subjects such as gymnasts and
martial artists to study the joint angle limits. The statistics of 3D poses in
this motion capture dataset is different from the previously existing motion
capture datasets such as CMU [11], Human 3.6M [15], and HumanEva [28], because
of their intention to explore the joint angle limits rather than performing
and recognizing typical human actions. Figure 2 shows the t-SNE visualization
[36] of poses from Akhter&Black motion Capture Dataset (ABCD) versus H36M in
two dimensions. One can see that the “ABCD” dataset is more uniformly
distributed compared to the H36M dataset. We randomly selected 4 poses from
the dense and surrounding sparse areas in the H36M t-SNE map and have shown
the corresponding images. One can see that all of the four samples selected
from the dense areas correspond to standing poses whereas all of the four
samples selected from sparse areas correspond to sitting poses.

Training and testing a 3D model on a similarly biased dataset with excessive
repetition of some poses will result in reduced performance on novel or rarely
seen poses. As a simple demonstration, we learned a GMM 3D pose model [29]
from a uniformly sampled set of Human 3.6M poses (all 15 actions) and
evaluated the likelihood of 3D poses per action under this model. The average
likelihood per action (up to a scaling factor) was: Directions 0.63,
Discussion 0.74, Eating 0.56 , Greeting 0.63 , Phoning 0.28 , Posing 0.38 ,
Purchases 0.55 , Sitting 0.07 , Sitting Down 0.07 , Smoking 0.47 , Taking
Photo 0.23 , Waiting 0.33 , Walking 0.64 , Walking Dog 0.29 , and Walk
Together 0.25. According to the GMM model, the “Discussion” poses are on
average almost 10 times more likely than “Sitting” poses which is due to the
dataset and consequently the model bias. The EM algorithm used to learn the
GMM model attempts to maximize the likelihood of all samples which will lead
to a biased model if the training dataset is biased. Obviously, any solely
data-driven model learned from a biased dataset that does not cover the full
range of motion of human body can suffer from lack of generalization to novel
or rarely seen yet anatomically plausible poses.

We propose a novel generative model on human 3D poses uniform in the space of
physically valid poses (satisfying the constraints from [1]). Since our model
is constrained only by the anatomical limits of human body it does not suffer
from dataset bias which is intrinsically helpful to diversify pose hypotheses.
Note that the pose-conditioned anatomical constraints calculated in [1] was
originally used in a constrained optimization framework for single 3D pose
estimation and turning those constraints into a generative model to produce
uniform samples is not trivial. One of our main contributions is a pose-
conditioned generative model which has not been done previously. We generate
multiple anatomically-valid and diverse pose hypotheses consistent with the 2D
joint detections to investigate the importance of having multiple pose
hypotheses under depth and missing-joints (_e.g_. caused by occlusion)
ambiguities. In the recent years, we have witnessed impressive progress in
accurate 2D pose estimation of human in various pose and appearances which is
made possible thanks to deep neural networks and lots of annotated 2D images.
We take advantage of the recent advancement in human 2D pose estimation and
seed our multi-hypotheses pose generator by an off-the-shelf 3D pose
estimator. Namely, we use the “Stacked Hourglass” 2D joint detector [19] and
the 2D-to-3D pose estimators of Akhter&Black [1] and Zhou et al. [42] to
estimate the 3D torso and projection matrix. However, note that to our generic
approach does not rely on any specific 2D/3D pose estimator and can easily
adopt various 2D/3D pose estimators.

After briefly discussing some related works in subsection 1.1 we propose our
approach in section 2. Our experimental results based on multiple 3D pose
estimation baselines is given in section 3. We conclude in section 4.

(a)

(b)

Figure 2: (a): The t-SNE visualization of poses from the H36M (fist from left)
and ABCD (second from left). (b): The images corresponding to the random
selection of poses from the dense (top row in right) and sparse (bottom row in
right) area of the H36M t-SNE map confirm the dataset bias toward standing
poses compared to sitting poses.

###  1.1 Related Work

There are quite a few works in the human pose estimation literature that are
directly or indirectly related to our work. Reviewing the entire literature is
obviously beyond the scope of this paper. Several areas of research are
related to our work such as 2D joiont detection, 3D pose estimation, and
generative 3D pose modeling. Due to the advancements made by deep neural
networks, the most recent works on 2D joint detection are based on
convolutional neural networks (CNN) [35, 9, 34, 10, 40, 39, 38, 19, 6, 26]
compared to the traditional hand-crafted feature based methods [27, 41, 12].
On the other hand, most of the 3D pose estimation methods use sparse coding
based on an overcomplete dictionary of basis poses to represent a 3D pose and
fit the 3D pose projection to the 2D joint detections [24, 37, 1, 42, 43].
Some works [8, 25, 26] try to train a deep network to directly predict 3D
poses. However, purely discriminative approaches for 3D structure prediction
(such as [8]) are usually very sensitive to data manipulation. On the other
hand, it has been shown that the deep networks are very effective and more
robust at detecting 2D templates (compared to 3D structures) such as human 2D
body parts in images [19].

We use conditional sampling from our generative model to generate multiple
consistent pose hypotheses. A number of previous works [7, 30, 2, 4, 5] have
used sampling for human pose estimation. However, the sampling performed by
these works are for purposes different from our goal to generate multiple
diverse and valid pose hypotheses. For example, Amin et al. [2] use a mixture
of pictorial structures and perform inference in two stages where the first
stage reduces the search space for the second inference stage by generating
samples for the 2D location of each part.

Some more closely related works include [33, 22, 16, 20, 23, 31, 17, 32].
Sminchisescu and Triggs [33] search for multiple local minima of their fitting
cost function using a sampling mechanism based on forwards/backwards link
flipping to generate pose candidates. Pons-Moll et al. [22] use inverse
kinematics to sample the pose manifold restricted by the input video and IMU
sensor cues in a particle filter framework. Lee and Cohen [16] use proposal
maps to consolidate the evidence and generating 3D pose candidates during the
MCMC search where they model the measurement uncertainty of 2D position of
joints using a Gaussian distribution. Their MCMC approach suffers from high
computational cost. Park and Ramanan [20] generate non-overlapping diverse
pose hypotheses (only in 2D) from a part model. One interesting work is the
“Posebit” by Pons-Moll et al. [23] that can retrieve pose candidates from a
MoCap dataset of 3D poses given answers to some mid-level queries such as “Is
the right hand in front of torso?” using decision trees. This approach is
heavily dependent on the choice of MoCap dataset and cannot generalize to
unseen poses. Simo-Serra1 et al. [31] model the 2D and 3D poses jointly in a
Bayesian framework by integrating a generative model and discriminative 2D
part detectors based on HOGs. Lehrmann et al. [17] learn a generative model
from the H36M MoCap dataset whose graph structure (not a Kinematic chain) is
learned using the Chow-Liu algorithm. Simo-Serra et al. [32] propagate the
error in the estimation of 2D joint locations (modeled using Gaussian
distributions) into the weights of dictionary elements in a sparse coding
framework; then by sampling the weights, some 3D pose samples are generated
and sorted based on the SVM score on joint distance features. However, their
approach does not guarantee that the joint angle constraints are satisfied and
do not address the depth ambiguity. We impose “pose-conditioned” joint angle
and bone length constrains to ensure pose validity of samples from our
generative model which has not been done before. In addition, our unbiased
generative model restricted only by anatomical constrains helps in generating
more diverse 3D pose hypotheses.

##  2 The Proposed Method

Figure 3: “Stacked Hourglass” 2D joint detector [19] in the absence and
presence of occlusion. On the right-hand-side of each image are the
corresponding heatmaps for joints.

Since our approach is closely related to the joint-angle constraints used in
[1], we find it helpful for better readability to briefly review this work. To
represent the human 3D pose by its joints let $\mathbf{X}$ denote the matrix
corresponding to $P$ kinematic joints in the 3D space namely
$\mathbf{X}=[\mathbf{X}_{1}...\mathbf{X}_{P}]\in\mathcal{X}\subset{\rm
I\\!R}^{3\times P}$ where $\mathcal{X}$ denotes the space of valid human
poses. Akhter&Black [1] (similar to [24, 42]) assumed that all of the 2D
joints are observed and estimated a single 3D pose by solving the following
optimization problem:

$\displaystyle\min_{\mathbf{\omega},s,\mathbf{R}}\hskip
1.42262ptC_{r}+C_{p}+\beta C_{l},\vspace{-2mm}$ (1)

where, $C_{r}$ is a measure of fitness between the estimated 2D joints
$\hat{\mathbf{x}}\in{\rm I\\!R}^{2\times P}$ and the projection and
translation of estimated 3D pose
$\hat{\mathbf{X}}=[\hat{\mathbf{X}}_{1}...\hat{\mathbf{X}}_{P}]\in{\rm
I\\!R}^{3\times P}$ to the 2D image coordinate system in a weak perspective
camera model (orthographic projection) with scaling factor $s\in{\rm
I\\!R}^{+}$, rotation $\mathbf{R}\in SO(3)$, and translation
$\mathbf{t}\in{\rm I\\!R}^{2\times 1}$, defined as:

$\displaystyle
C_{r}=\sum_{i=1}^{P}\lVert\hat{\mathbf{x}}_{i}-s\mathbf{R}_{1:2}\hskip
1.42262pt\hat{\mathbf{X}}_{i}+\mathbf{t}\rVert_{2}^{2},\vspace{-2mm}$ (2)

where, $\mathbf{R}_{1:2}$ denotes the first two rows of the rotation matrix.
Note that if the origin of the 3D world coordinate system gets mapped to the
origin of the 2D image coordinate system then $\mathbf{t}=\mathbf{0}$; this is
usually implemented by centering the 2D and 3D poses. Authors used a sparse
representation of the 3D poses similar to [24] where the 3D pose is
represented by a sparse linear combination of bases selected using the
Orthogonal Matching Pursuit (OMP) algorithm [18] from an overcomplete
dictionary of pose atoms, namely
$\hat{\mathbf{X}}=\mathbf{\mu}+\sum_{i\in\mathcal{I^{*}}}\omega_{i}\mathbf{D}_{i}$,
where $\mathbf{\mu}$ is the mean pose obtained by averaging poses from the CMU
motion capture dataset [11] and $\mathcal{I^{*}}$ denotes the indices of
selected bases using OMP with weights $\mathbf{\omega}_{i}$. An overcomplete
dictionary of bases was built by concatenating PCA bases from poses of
different action classes in the CMU dataset after bone length normalization
and Procrustes aligned. The second term $C_{p}$ in equation (1) is equal to
zero if the estimated pose $\hat{\mathbf{X}}$ has valid joint angles for limbs
and infinity otherwise. According to the pose-conditioned constraints in [1] a
pose has valid joint angles if the upper arms/legs’ joint angles map to a 1 in
the corresponding occupancy matrix (learned from the ABCD dataset) and the
lower arms/legs satisfy two conditions that prevent these bones from bending
beyond feasible joint-angle limits (inequalities (4) and (2.1)). The term
$C_{l}$ in equation (1) penalizes the difference between the squares of the
estimated $i^{\text{th}}$ bone length $l_{i}$ and the normalized mean bone
length $\bar{l}_{i}$ _i.e_., $C_{l}=\sum_{i=1}^{N}|l_{i}^{2}-\bar{l}_{i}^{2}|$
(normalized mean bones calculated from the CMU dataset) with weight $\beta$.
Note that [1] does not introduce any generative pose model.

As we mentioned earlier, 3D pose estimation from 2D landmark points in
monocular RGB images is inherently an ill-posed problem because of losing the
depth information. There can be multiple valid 3D poses with similar 2D
projection even if all of the 2D joints are observed (see Figure 1). The
uncertainty and number of possible valid poses can further increase if some of
the joints are missing. The missing joints scenario is more realistic because
it happens when either these joints exist in the image but are not confidently
detected, due to occlusion and clutter, or do not exist within the borders of
the image _e.g_. when only the upper body is visible similar to images from
the FLIC dataset [27]. It is observed that thresholding the confidence score
obtained from some deep 2D joint detectors (_e.g_. [19, 21, 14]) can be
reasonably used as an indicator for the confident detection of a joint. Figure
3 shows the the output of “Stacked Hourglass” 2D joint detector [19] in the
absence and presence of a table occluder segmented out from the Pascal VOC
dataset [13] and pasted on the left hand of the human subject. On the right-
hand-side of each image is shown the heatmap for each joint. It can be seen
that the level of the two heatmaps corresponding to the left elbow and left
wrist drop after placing the table occluder on the left hand. Newell et al.
[19] used the heatmap mean as a confidence measure for detection and threshold
it at 0.002 to determine visibility of a joint. Obviously, invisibility of
some joints in the image can result in multiple hallucinations for the 2D/3D
locations of the joints. Let $S_{\text{o}}$ and $S_{\text{m}}$ denote the set
of observed and missing joints, respectively. We have $S_{\text{o}}\cap
S_{\text{m}}=\emptyset$ and $S_{\text{o}}\cup S_{\text{m}}=\\{1,2,...,P\\}$,
and let $\mathbf{\alpha}=\\{\alpha_{i}\\}_{i\in S_{\text{o}}}$ denote a set of
normalized joint scores from the 2D joint detectors such that
$\frac{1}{|S_{\text{o}}|}\sum_{i\in S_{\text{o}}}\alpha_{i}=1$. The missing
joints are detected by comparing the confidence score of 2D joint detector
with a threshold (0.002 in the case of using Hourglass). For the case of
missing joints, we modify the fitness measure to:

$\displaystyle C_{r}=\sum_{i\in
S_{\text{o}}}\alpha_{i}\lVert\hat{\mathbf{x}}_{i}-s\mathbf{R}_{1:2}\hskip
1.42262pt\hat{\mathbf{X}}_{i}+\mathbf{t}\rVert_{2}^{2}.\vspace{-2mm}$ (3)

The scores are normalized because they have to be in a comparable range with
respect to the $C_{l}$ term in equation (1) otherwise either $C_{r}$ is
suppressed/ignored in the case of very small confidence scores or the same
happens to $C_{l}$ in the case of very large scores. For example, if the mean
of heatmaps from the Hourglass joint detector are directly (without
normalization) used as scores the $C_{r}$ term will be drastically suppressed
since the heatmaps are full of close-to-zero values. Note that the
optimization problem in equation (1) with the updated $C_{r}$ term according
to equation (3) still outputs a full 3D pose even under missing joints
scenario because the 3D pose is constructed by a linear combination of full
body basis. However, there is no reason that the output 3D pose should have a
close to correct 2D projection due to the missing joint ambiguity added to the
depth ambiguity. Optimizing $C_{r}$ is a non-convex optimization problem over
the 3D pose and projection matrix. To obtain an estimate of the 3D torso and
projection matrix, we tried both iterating between optimizing over the
projection matrix and 3D pose used in [1] as well as the convex relaxation
method in [42] as will be presented in the experimental results section. Note
that the torso pose variations are much fewer than the full-body. The torso
plane is usually vertical and not as flexible as the full body. Hence, it is
much easier to robustly estimate its 3D pose and the corresponding camera
parameters.

To generate multiple diverse 3D pose hypotheses consistent with the output of
2D joint detector, we cluster samples from a conditional distribution given
the collected 2D evidence. For this purpose, we follow a rejection sampling
strategy. Before discussing conditional sampling in subsection 2.2 we describe
unconditional sampling as follows.

###  2.1 Unconditional Sampling

Given the rigidity of human torso compared to the limbs (hands/legs), the
joints corresponding to the torso including thorax, left/right hips, and
left/right shoulders can be represented using a small size dictionary after an
affine transformation/normalization. Given the torso, the upper arms/legs and
head are anatomically restricted to be within certain angular limits. The
plausible angular regions for the upper arms/legs and head can be represented
using an occupancy matrix [1]. This occupancy matrix is a binary matrix that
assigns 1 to a discretized azimuthal $\theta$ and polar $\phi$ angle if these
angles are anatomically plausible and 0 otherwise. These angular positions are
calculated in the local Cartesian coordinate system whose two axis are the
“backbone” vector and either the “right shoulder $\rightarrow$ left shoulder”
vector (for the upper arms and head) or the “right hip $\rightarrow$ left hip”
vector (for the upper hips). Hence, to generate samples for the upper
arms/legs and head we just need to take samples from the occupancy matrix at
places where the value is 1 and get the corresponding azimuthal and polar
angles. Given the azimuthal and polar angles of the head we just need to
travel in this direction for the length of the head; we do the same for the
length of upper arms and legs to reach the elbows and knees, respectively. The
normalized length of the bones is sampled from a Beta distribution with
limited range under the constraint that similar bones have similar length
_e.g_. both upper arms have the same length.

According to [1], the lower arm/leg bone $\mathbf{b}_{p_{1}\rightarrow
p_{2}}=\mathbf{X}_{p_{2}}-\mathbf{X}_{p_{1}}$, where $p_{2}$ and $p_{1}$
respectively correspond to either “wrist and elbow” or “ankle and knee” is at
a plausible angle if it satisfies two constraints. The first constraint is:

$\displaystyle\mathbf{b}^{\top}\mathbf{n}+d<0,$ (4)

where $\mathbf{n}$ and $d$ are functions of the azimuthal $\theta$ and polar
$\phi$ angles of their parent bone namely the upper arm or leg (resulting in
pose-dependent joint angle limits) learned from the ABCD dataset. The above
inequality defines a separating plane, with normal vector $\mathbf{n}$ and
distance from origin $d$, that attempts to prevent the wrist and ankle from
bending in a direction that is anatomically impossible. Obviously, for a very
negative offset vector $d$ this constrain is always satisfied. Therefore,
during learning of $\mathbf{n}$ and $d$ the second norm of $d$ is minimized,
namely $\min_{\mathbf{n},d}\hskip 1.42262pt\lVert d\rVert_{2}\hskip
5.69054pt\text{s.t.}\hskip 5.69054pt\mathbf{B}^{\top}\mathbf{n}<-d\mathbf{1}$,
where $\mathbf{B}$ is a matrix built by column-wise concatenation of all
$\mathbf{b}$ instances in the ABCD dataset whose parents are at the same
$\theta$ and $\phi$ angular location. The second constraint to satisfy is that
the projection of normalized $\mathbf{b}$ (to unit length) onto the separating
plane using the orthonormal projection matrix
$\mathbf{T}=[\mathbf{T}_{1};\mathbf{T}_{2};\mathbf{T}_{3}]$, whose first row
$\mathbf{T}_{1}$ is along $\mathbf{n}$, has to fall inside a bounding box with
bounds $[bnd_{1},bnd_{2}]$ and $[bnd_{3},bnd_{4}]$, namely:

$\displaystyle
bnd_{1}\leq\mathbf{T}_{2}\mathbf{b}/\lVert\mathbf{b}\rVert_{2}\leq bnd_{2},$
$\displaystyle
bnd_{3}\leq\mathbf{T}_{3}\mathbf{b}/\lVert\mathbf{b}\rVert_{2}\leq bnd_{4},$
(5)

where, bounds $bnd_{1}$, $bnd_{2}$, $bnd_{3}$, and $bnd_{4}$ are also learned
from the ABCD dataset. To generate a sample $\mathbf{b}$ that satisfies the
above constraints, we first generate two random values
$u_{2}\in[bnd_{1},bnd_{2}]$ and $u_{3}\in[bnd_{3},bnd_{4}]$ and set
$u_{1}=(\max(1-u_{2}^{2}-u_{3}^{2},0))^{1/2}$. We then generate two candidates
$\mathbf{u}^{\pm}=(\pm u_{1},u_{2},u_{3})/\lVert(u_{1},u_{2},u_{3})\rVert_{2}$
from which only one can be on the valid side of the separating plane
satisfying inequality (4). To check, we first undo the projection and
normalization by $\mathbf{b}^{\pm}=l\hskip
1.42262pt\mathbf{T}^{-1}\mathbf{u}^{\pm}$, where $l$ is a sample from the bone
length distribution on $\mathbf{b}$. A sample “$\mathbf{b}$” is accepted only
if it satisfies inequality (4). Note that similar bones have the same length
therefore we sample their length only once for each pose. The prior model can
be written as below according to a Bayesian graph on the kinematic chain:

$\displaystyle
p(\mathbf{X})=p(\mathbf{X}_{i\in\text{torso}})p(\mathbf{X}_{\text{head}}|\mathbf{X}_{i\in\text{torso}})\times$
$\displaystyle p(\mathbf{X}_{i\in\text{ l/r
elbow}}|\mathbf{X}_{i\in\text{torso}})p(\mathbf{X}_{i\in\text{ l/r
wrist}}|\mathbf{X}_{i\in\text{ l/r
elbow}},\mathbf{X}_{i\in\text{torso}})\times$ $\displaystyle
p(\mathbf{X}_{i\in\text{ l/r
knee}}|\mathbf{X}_{i\in\text{torso}})p(\mathbf{X}_{i\in\text{ l/r
ankle}}|\mathbf{X}_{i\in\text{ l/r
knee}},\mathbf{X}_{i\in\text{torso}}),\vspace{-1.5mm}$ (6)

where $p(\mathbf{X}_{i\in\text{torso}})$ is the probability of selecting a
torso from the torso dictionary which we assumed is uniform. The torso joints
$\mathbf{X}_{i\in\text{torso}}$ are used to determine the local coordinate
system for the rest of the joints. We have removed torso joints in the
equations below for notational convenience. We have:

$\displaystyle
p(\mathbf{X}_{i})=\frac{1}{l_{\text{bone}}^{2}|\sin(\phi_{i})|}p(l_{\text{bone}})\hskip
1.42262ptp(\mathbf{\theta}_{i},\mathbf{\phi}_{i}),$ (7)

for $(i,\text{bone})$ being from (l/r knee, upper leg) , (head, neck + head
bone), or (l/r elbow, upper arm). The multiplier factor in (7), which is the
inverse of Jacobian determinant for a transformation from the Cartesian to
spherical coordinate system, is to ensure that the left side sums up to one if
$\int_{l}\int_{\theta}\int_{\phi}p(l)\hskip
1.42262ptp(\mathbf{\theta},\mathbf{\phi})d\phi\hskip 1.42262ptd\theta\hskip
1.42262ptdl=1$, since $dx\hskip 1.42262ptdy\hskip
1.42262ptdz=l^{2}|\sin(\phi)|\hskip 1.42262ptdl\hskip 1.42262ptd\theta\hskip
1.42262ptd\phi$. For lower limbs we have:

$\displaystyle p(\mathbf{X}_{i}|\mathbf{X}_{pa(i)})\propto
p(l_{\text{bone}})\mathbf{1}_{valid}(\mathbf{X}_{i},\mathbf{X}_{pa(i)})$ (8)

where $(i,pa(i),\text{bone})$ is from (l/r wrist, l/r elbow, forearm) or (l/r
ankle, l/r knee, lower leg), and
$\mathbf{1}_{valid}(\mathbf{X}_{i},\mathbf{X}_{pa(i)})$ is an indicator
function that nulls the probability of configurations whose angles does not
satisfy the constraints in inequalities (4) and (2.1) for
$\mathbf{b}=\mathbf{X}_{i}-\mathbf{X}_{pa(i)}$. Conditional sampling is
carried out by rejection sampling discussed in the next subsection.

(a) (b)

Figure 4: (a): The input image and the corresponding 3D pose. (b): Generation
of five diverse 3D pose hypotheses consistent with the 2D joint detections.

###  2.2 Conditional Sampling

We run a 2D joint detector on the input image $I$ and get an estimate of the
2D joint locations $\hat{\mathbf{x}}$ with confidence scores
$\mathbf{\alpha}$. Then, to obtain a reasonable estimate of torso
$\hat{\mathbf{X}}_{i\in\text{torso}}$ and camera parameters namely
$(\hat{\mathbf{R}},\hat{\mathbf{t}},\hat{s})$, we run a 2D-to-3D pose
estimator capable of handling missing joints (we modified [1] and [42] to
handle missing joints; see equation (3)). Note that we are not restricted to
any particular 2D/3D pose estimator and any 2D joint detector that estimates
2D joint locations $\hat{\mathbf{x}}$ and their confidence scores
$\mathbf{\alpha}$ and any 2D-to-3D pose estimator can be used in the initial
stage. We then assume that the estimated camera parameters and
$\hat{\mathbf{X}}_{i\in\text{torso}}$ are reasonably well estimated and keep
them fixed. Note that the human torso and its pose (usually vertical) does not
vary much compared to the whole body pose. We do not include the estimated
camera parameters and 3D torso in our formulation below for notational
convenience. From the Bayes rule we have:

$\displaystyle p(\mathbf{X}|\hat{\mathbf{x}},\mathbf{\alpha})\propto
p(\mathbf{X})\hskip
1.42262ptp(\hat{\mathbf{x}},\mathbf{\alpha}|\mathbf{X}).\vspace{-2mm}$ (9)

We define:

$\displaystyle
p(\hat{\mathbf{x}},\mathbf{\alpha}|\mathbf{X})\propto\prod_{i\in\text{ limb
}\cap S_{\text{o}}}\mathbf{1}(\lVert\hat{\mathbf{x}}_{i}-\hat{s}\hskip
1.42262pt\hat{\mathbf{R}}_{1:2}\hskip
1.42262pt\mathbf{X}_{i}+\hat{\mathbf{t}}\rVert_{2}<\tau_{i})\vspace{-2mm}$

where $\mathbf{1}(.)$ is the indicator function depending on the 2D distance
between detected joints and the projected 3D pose under an acceptance
threshold defined by $\tau_{i}=0.25\hskip 1.42262pt\hat{s}\hskip
1.42262pt\bar{l}_{\text{limb}}/\alpha_{i}$, where $\bar{l}_{\text{limb}}$ is
the mean limb length, $\hat{s}$ is the estimated scaling factor, $\alpha_{i}$
is the $i^{\text{th}}$ joint normalized confidence score, and the factor
$0.25$ was chosen empirically. The likelihood function defined above accepts
prior (unconditional) samples $\mathbf{X}^{(q)}\sim p(\mathbf{X})$ whose
projected joints to the image coordinate system are within a distance not
greater than thresholds $\tau_{i}$ from detected limb joints. The inverse
proportion of the threshold to the confidence $\alpha_{i}$ allows acceptance
in a larger area if the confidence score is smaller for the $i^{\text{th}}$
limb joint and therefore considering the 2D joint detection uncertainty. Note
that there is no indicator function in the likelihood function for the missing
limb joints which allows acceptance of all anatomically plausible samples for
limb joints from $S_{\text{m}}$. Note that even though torso pose estimation
is a much easier problem compared to the full body pose estimation, a poorly
estimated torso, _e.g_. due to occlusion, can adversely affect the quality of
conditional 3D pose samples.

###  2.3 Generating Diverse Hypotheses

The diversification is implemented in two stages: (I) we sampled the occupancy
matrix at 15 equidistant azimuth and 15 equidistant polar angles for the upper
limbs and accept the samples if the occupancy matrix had a 1 at these
locations. For the lower limbs, we sampled 5 equidistant points along each
$u_{2}$ and $u_{3}$ directions between $[bnd_{1},bnd_{2}]$ and
$[bnd_{3},bnd_{4}]$, respectively. (II) To generate fewer number of pose
hypothesis, we use the kmeans++ algorithm [3] to cluster the posterior samples
into a desired number of diverse clusters and take the nearest neighbor 3D
pose sample to each centroid as one hypothesis. Kmeans++ operates the same as
Kmeans clustering except that it uses a diverse initialization method to help
with diversification of final clusters. Note that we cannot take the centroids
as hypotheses since there is no guarantee that the mean of 3D poses is still a
valid 3D pose. Figure 4 shows five hypotheses given the output of Hourglass 2D
joint detector for the top-left image and detections shown by yellow points.
In Figure 4, the 2D detection of joints are shown by the black skeleton and
the diversified hypotheses that are consistent with the 2D detections are
shown by the blue skeletons. It can be seen that even though the 2D projection
of these pose hypotheses are very similar, they are quite different in 3D. To
generate the pose hypotheses in Figure 4, we estimated the 3D torso and
projection matrix using [1]. s

##  3 Experimental Results

Method | Directions | Discussion | Eating | Greeting | Phoning | Posing | Purchases | Sitting | SitDown
---|---|---|---|---|---|---|---|---|---
Ours (No KM++/[42]) | 63.12 | 55.91 | 58.11 | 64.48 | 68.69 | 61.27 | 55.57 | 86.06 | 117.57
Ours (k=20/[42]) | 77.08 | 71.15 | 75.39 | 79.01 | 84.68 | 74.90 | 72.37 | 102.17 | 131.46
Ours (k=5/[42]) | 82.86 | 77.52 | 81.60 | 85.20 | 90.93 | 80.46 | 78.75 | 109.27 | 138.71
Zhou et al. [42] | 80.51 | 74.56 | 73.95 | 85.43 | 88.96 | 82.02 | 76.21 | 107.43 | 146.47
|  |  |  |  |  |  |  |  |
Ours (k=5/[1]) | 105.14 | 100.28 | 107.75 | 106.88 | 111.44 | 105.74 | 101.18 | 124.87 | 147.48
Akhter&Black [1] | 133.80 | 128.03 | 124.47 | 133.47 | 133.93 | 136.63 | 128.30 | 133.61 | 162.01
|  |  |  |  |  |  |  |  |
Chen et al. [8] | 145.37 | 139.11 | 140.24 | 149.13 | 149.61 | 154.30 | 147.04 | 161.49 | 200.06
|  |  |  |  |  |  |  |  |

| Smoking | TakingPhoto | Waiting | Walking | WalkingDog | WalkTogether | Average |
---|---|---|---|---|---|---|---|---
Ours (No KM++/[42]) | 71.02 | 71.21 | 66.29 | 57.07 | 62.50 | 61.02 | 67.99 |
Ours (k=20/[42]) | 85.90 | 84.49 | 80.41 | 71.57 | 78.41 | 74.92 | 82.93 |
Ours (k=5/[42]) | 91.79 | 90.06 | 86.43 | 77.93 | 85.45 | 81.49 | 89.23 |
Zhou et al. [42] | 90.61 | 93.43 | 85.71 | 80.03 | 90.89 | 85.73 | 89.46 |
|  |  |  |  |  |  |  |
Ours (k=5/[1]) | 113.61 | 105.58 | 105.80 | 100.28 | 106.25 | 104.63 | 109.79 |
Akhter&Black [1] | 135.75 | 132.92 | 133.93 | 133.84 | 131.77 | 134.80 | 134.48 |
|  |  |  |  |  |  |  |
Chen et al. [8] | 152.37 | 159.18 | 152.67 | 148.20 | 156.10 | 147.71 | 153.51 |
|  |  |  |  |  |  |  |

Table 1: Quantitative comparison on the Human3.6M dataset evaluated in 3D by
mean per joint error (mm) for all actions and subjects whose ground-truth 3D
poses were provided.

We empirically evaluated the proposed “multi-pose hypotheses” approach on the
recently published Human3.6M dataset [15]. For evaluation, we used images from
all 4 cameras and all 15 actions associated with 7 subjects for whom ground-
truth 3D poses were provided namely subjects S1, S5, S6, S7, S8, S9, and S11.
The original videos (50 fps) were downsampled (in order to reduce the
correlation of consecutive frames) to built a dataset of 26385 images. For
further evaluation, we also built two rotation datasets by rotating H36M
images by $30$ and $60$ degrees. We evaluated the performance by the mean per
joint error (millimeter) in 3D by comparing the reconstructed pose hypotheses
against the ground truth. The error was calculated up to a similarity
transformation obtained by Procrustes alignment. The results are summarized in
Table 1 for various methods and actions. For a fair comparison, the limb
length of the reconstructed poses from all methods were scaled to match the
limb length of the ground-truth pose. The bone length matching obviously
lowers the mean joint errors but makes no difference in our comparisons. One
can see that the best (lowest Euclidean distance from the ground-truth pose)
out of only 5 generated hypotheses by using [1] as baseline for 3D torso and
projection matrix estimation is considerably better than the single 3D pose
output by [1] for all actions. We also used the 2D-to-3D pose estimator by
Zhou et al. [42] with convex-relaxation as baseline and observed considerable
improvement compared to [1] in both 3D pose and projection matrix estimation.
Using [42] as baseline to estimate the 3D torso and projection matrix we
generated multiple 3D pose hypotheses. Since the accuracy of [42] is already
high, the best out of 5 pose hypotheses cannot significantly lower the average
joint distance from the single 3D pose output by [42]. However, by increasing
the number of hypotheses we started to observe improvement. Table 1 also
includes the best hypothesis out of conditional samples from only the first
diversification stage _i.e_., by diversifying conditional samples and using no
kmeans++ clustering (shown by No KM++), using [42] as base. This achieves the
lowest joint error in comparison to other baselines. The pose hypotheses can
be generated very quickly ($<$ 2 seconds) in Matlab on an Intel i7-4790K
processor.

Figure 5: PCK curves for the H36M dataset (original), H36M rotated by 30 and
60 degrees respectively from left to right. The y-axis is the percentage of
correctly detected joints in 3D for a given distance threshold in millimeter
(x-axis).

Method | Directions | Discussion | Eating | Greeting | Phoning | Posing | Purchases | Sitting | SitDown
---|---|---|---|---|---|---|---|---|---
Ours (k=5/[1]) | 98.44 | 93.70 | 102.62 | 97.50 | 96.29 | 98.90 | 93.32 | 105.51 | 110.07
Akhter&Black [1] | 118.02 | 112.55 | 111.27 | 117.46 | 111.77 | 122.27 | 112.23 | 107.27 | 126.95
|  |  |  |  |  |  |  |  |
Ours (k=5/[1]) | 108.60 | 105.85 | 105.63 | 109.01 | 105.47 | 109.93 | 102.01 | 111.25 | 119.57
Akhter&Black [1] | 153.80 | 149.14 | 135.44 | 155.06 | 139.62 | 156.46 | 149.05 | 126.33 | 141.89
|  |  |  |  |  |  |  |  |
Ours (k=5/[1]) | 125.03 | 121.77 | 115.13 | 124.11 | 116.92 | 123.75 | 116.42 | 119.63 | 130.81
Akhter&Black [1] | 185.57 | 180.43 | 158.55 | 185.65 | 162.39 | 185.78 | 178.81 | 145.15 | 155.29
|  |  |  |  |  |  |  |  |

| Smoking | TakingPhoto | Waiting | Walking | WalkingDog | WalkTogether | Average | Average Diff.
---|---|---|---|---|---|---|---|---
Ours (k=5/[1]) | 97.53 | 97.63 | 99.43 | 90.23 | 97.27 | 95.21 | 98.24 |
Akhter&Black [1] | 113.22 | 120.61 | 119.97 | 115.81 | 116.60 | 115.62 | 116.11 | $\boxed{17.87}$
|  |  |  |  |  |  |  |
Ours (k=5/[1]) | 107.76 | 107.05 | 111.34 | 108.38 | 106.96 | 110.28 | 108.61 |
Akhter&Black [1] | 142.98 | 152.65 | 155.27 | 155.18 | 151.88 | 155.00 | 147.98 | $\boxed{39.37}$
|  |  |  |  |  |  |  |
Ours (k=5/[1]) | 120.60 | 118.38 | 127.13 | 125.89 | 121.61 | 127.62 | 122.32 |
Akhter&Black [1] | 165.47 | 177.44 | 186.20 | 189.66 | 183.01 | 186.25 | 175.04 | $\boxed{\textbf{52.72}}$
|  |  |  |  |  |  |  |

Table 2: Quantitative comparison on the Human3.6M dataset when 0 (top pair), 1
(middle pair), and 2 (bottom pair) limb joints are missing.

We also used Deep3D of Chen et al. [8] as another baseline. The Deep3D [8] is
a 3D pose estimator that directly regresses to the 3D joint locations directly
from a monocular RGB input image. Deep3D had the highest mean joint errors as
shown in Table 1. We also observed that the pre-trained Deep3D is very
sensitive to image rotation and usually outputs an anatomically implausible 3D
pose if the input image is rotated. But other 2D-to-3D pose estimation
baselines which decouple the projection matrix and the 3D pose are quite
robust to rotation of the input image. Figure 5 shows the Percentage of
Correct Keypoints (PCK) versus an acceptance distance threshold in millimeter
for various baselines and H36M dataset variations namely the original H36M and
30/60 degree rotations. One can see that the PCK of Deep3D drops drastically
by rotating the input image. This is partly due to insufficient number of
tilted samples in the training set (H36M plus synthetic images). One of the
main problems of purely discriminative approaches such as [8] is their extreme
sensitivity to data manipulation. On the other hand, humans can learn from a
few examples and still not suppress the rarely seen cases compared to the
frequently seen ones.

In a realistic scenario with occlusion, the location of some 2D joints cannot
be accurately detected. The added uncertainty caused by occlusion makes one
expect a larger average estimation error for the estimated 3D pose from a
single-output pose estimator compared to the best 3D pose hypothesis. To test
this, we ran experiments with different number of missing joints (0, 1 and 2)
selected randomly from the limb joints including l/r elbow, l/r wrist, l/r
knee, and l/r ankle. Table 2 shows the mean per joint errors for the 3D pose
estimated by the modified version of Akhter&Black [1] that can handle missing
joints compared to the best out of five hypotheses generated by our method
when 0, 1, and 2 limb joints are missing. In this test, we used the ground-
truth 2D location of the joints and randomly selected the missing joints. One
can see that by increasing the number of missing joints the performance gap
between the estimated 3D pose and the best 3D pose hypothesis increases. This
underscores the importance of having multiple hypothesis for more realistic
scenarios.

##  4 Conclusion

There usually exist multiple 3D poses consistent with the 2D location of
joints because of losing the depth information in monocular images. The
uncertainty in 3D pose estimation increases in the presence of occlusion and
imperfect 2D detection of joints. In this paper, we proposed a way to generate
multiple valid and diverse 3D pose hypotheses consistent with the 2D joint
detections. These pose hypotheses can be ranked later by more detailed
investigation of the image beyond the 2D joint locations or based on some
contextual information. To generate these pose hypotheses we used a novel
unbiased generative model that only enforces pose-conditioned anatomical
constraints on the joint-angle limits and limb length ratios. This was
motivated by the pose-conditioned joint limits from [1] after identifying bias
in typical MoCap datasets. Our compositional generative model uniformly spans
the full variability of human 3D pose which helps in generating more diverse
hypotheses. We performed empirical evaluation on the H36M dataset and achieved
lower mean joint errors for the best pose hypothesis compared to the estimated
pose by other recent baselines. The 3D pose output by the baseline methods
could also be included as one hypothesis but to investigate our hypothesis
generation approach we did not do so in the experimental results. Our
experiments show the importance of having multiple 3D pose hypotheses given
only the 2D location of joints especially when some of the joints are missing.
We hope our idea of generating multiple pose hypotheses inspire a new line of
future work in 3D pose estimation considering various ambiguity sources.

## References

  * [1] I. Akhter and M. J. Black.  Pose-conditioned joint angle limits for 3D human pose reconstruction.  In CVPR, pages 1446–1455, June 2015.
  * [2] S. Amin, M. Andriluka, M. Rohrbach, and B. Schiele.  Multi-view pictorial structures for 3d human pose estimation.  In British Machine Vision Conference (BMVC), September 2013.
  * [3] D. Arthur and S. Vassilvitskii.  K-means++: The advantages of careful seeding.  In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’07, pages 1027–1035, Philadelphia, PA, USA, 2007\. Society for Industrial and Applied Mathematics.
  * [4] V. Belagiannis, S. Amin, M. Andriluka, B. Schiele, N. Navab, and S. Ilic.  3d pictorial structures for multiple human pose estimation.  In CVPR, 2014.
  * [5] V. Belagiannis, S. Amin, M. Andriluka, B. Schiele, N. Navab, and S. Ilic.  3d pictorial structures revisited: Multiple human pose estimation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(10):1929–1942, 2016.
  * [6] A. Bulat and G. Tzimiropoulos.  Human pose estimation via convolutional part heatmap regression.  In ECCV, 2016.
  * [7] M. Burenius, J. Sullivan, and S. Carlsson.  3d pictorial structures for multiple view articulated pose estimation.  In CVPR, pages 3618–3625, 2013.
  * [8] W. Chen, H. Wang, Y. Li, H. Su, Z. Wang, C. Tu, D. Lischinski, D. Cohen-Or, and B. Chen.  Synthesizing training images for boosting human 3d pose estimation.  In 3D Vision (3DV), 2016.
  * [9] X. Chen and A. L. Yuille.  Articulated pose estimation by a graphical model with image dependent pairwise relations.  In Advances in Neural Information Processing Systems, pages 1736–1744, 2014.
  * [10] X. Chu, W. Ouyang, H. Li, and X. Wang.  Structured feature learning for pose estimation.  In CVPR, 2016.
  * [11] P. Doe.  Cmu human motion capture database. availabel online at:, 2003.
  * [12] M. Eichner, M. Marin-Jimenez, A. Zisserman, and V. Ferrari.  2d articulated human pose estimation and retrieval in (almost) unconstrained still images.  International Journal of Computer Vision, 99:190–214, 2012.
  * [13] M. Everingham, S. M. A. Eslami, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman.  The pascal visual object classes challenge: A retrospective.  International Journal of Computer Vision, 111(1):98–136, Jan. 2015\.
  * [14] E. Insafutdinov, L. Pishchulin, B. Andres, M. Andriluka, and B. Schiele.  Deepercut: A deeper, stronger, and faster multi-person pose estimation model.  In ECCV, May 2016.
  * [15] C. Ionescu, D. Papava, V. Olaru, and C. Sminchisescu.  Human3.6m: Large scale datasets and predictive methods for 3d human sensing in natural environments.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(7):1325–1339, jul 2014.
  * [16] M. W. Lee and I. Cohen.  Proposal maps driven mcmc for estimating human body pose in static images.  In CVPR, 2004.
  * [17] A. M. Lehrmann, P. V. Gehler, and S. Nowozin.  A non-parametric bayesian network prior of human pose.  In CVPR, pages 1281–1288, 2013.
  * [18] S. G. Mallat and Z. Zhang.  Matching pursuits with time-frequency dictionaries.  IEEE Transactions on Signal Processing, pages 3397–3415, Dec. 1993\.
  * [19] A. Newell, K. Yang, and J. Deng.  Stacked hourglass networks for human pose estimation.  In ECCV, May 2016.
  * [20] D. Park and D. Ramanan.  N-best maximal decoders for part models.  In ICCV, 2011.
  * [21] L. Pishchulin, E. Insafutdinov, S. Tang, B. Andres, M. Andriluka, P. Gehler, and B. Schiele.  Deepcut: Joint subset partition and labeling for multi person pose estimation.  In CVPR, June 2016.
  * [22] G. Pons-Moll, A. Baak, J. Gall, L. Leal-Taixé, M. Müller, H.-P. Seidel, and B. Rosenhahn.  Outdoor human motion capture using inverse kinematics and von mises-fisher sampling.  In ICCV, 2011.
  * [23] G. Pons-Moll, D. J. Fleet, and B. Rosenhahn.  Posebits for monocular human pose estimation.  In CVPR, pages 2345–2352, June 2014.
  * [24] V. Ramakrishna, T. Kanade, and Y. Sheikh.  Reconstructing 3d human pose from 2d image landmarks.  In ECCV, 2012.
  * [25] G. Rogez and C. Schmid.  Mocap-guided data augmentation for 3d pose estimation in the wild.  2016\.
  * [26] G. Rogez, P. Weinzaepfel, and C. Schmid.  Lcr-net: Localization-classification-regression for human pose.  In CVPR, 2017.
  * [27] B. Sapp and B. Taskar.  Modec: Multimodal decomposable models for human pose estimation.  In CVPR, pages 3674–3681, 2013.
  * [28] L. Sigal, A. Balan, and M. J. Black.  Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion.  International Journal of Computer Vision, 87:4–27, 2010.
  * [29] L. Sigal, S. Bhatia, S. Roth, M. J. Black, and M. Isard.  Tracking loose-limbed people.  In CVPR, June 2004.
  * [30] L. Sigal, M. Isard, H. Haussecker, and M. J. Black.  Loose-limbed people: Estimating 3D human pose and motion using non-parametric belief propagation.  International Journal of Computer Vision, 98(1):15–48, May 2011\.
  * [31] E. Simo-Serra, A. Quattoni, C. Torras, and F. Moreno-Noguer.  A joint model for 2d and 3d pose estimation from a single image.  In CVPR, pages 3634–3641, 2013.
  * [32] E. Simo-Serra, A. Ramisa, G. Alenyà, C. Torras, and F. Moreno-Noguer.  Single image 3d human pose estimation from noisy observations.  In CVPR, 2012.
  * [33] C. Sminchisescu and B. Triggs.  Kinematic jump processes for monocular 3d human tracking.  In CVPR, 2003.
  * [34] J. J. Tompson, A. Jain, Y. LeCun, and C. Bregler.  Joint training of a convolutional network and a graphical model for human pose estimation.  In Advances in neural information processing systems, pages 1799–1807, 2014.
  * [35] A. Toshev and C. Szegedy.  Deeppose: Human pose estimation via deep neural networks.  In CVPR, pages 1653–1660, 2014.
  * [36] L. van der Maaten and G. E. Hinton.  Visualizing high-dimensional data using t-sne.  Journal of Machine Learning Research, 9:2579–2605, 2008.
  * [37] C. Wang, Y. Wang, Z. Lin, A. L. Yuille, and W. Gao.  Robust estimation of 3d human poses from a single image.  In CVPR, 2014.
  * [38] S.-E. Wei, V. Ramakrishna, T. Kanade, and Y. Sheikh.  Convolutional pose machines.  In CVPR, June 2016.
  * [39] H. L. X. W. Xiao Chu, Wanli Ouyang.  Crf-cnn: Modeling structured information in human pose estimation.  In NIPS, 2016.
  * [40] W. Yang, W. Ouyang, H. Li, and X. Wang.  End-to-end learning of deformable mixture of parts and deep convolutional neural networks for human pose estimation.  In CVPR, 2016.
  * [41] Y. Yang and D. Ramanan.  Articulated pose estimation with flexible mixtures-of-parts.  In CVPR, 2011.
  * [42] X. Zhou, S. Leonardos, X. Hu, and K. Daniilidis.  3d shape estimation from 2d landmarks: A convex relaxation approach.  In CVPR, pages 4447–4455, June 2015.
  * [43] X. Zhou, M. Zhu, S. Leonardos, K. G. Derpanis, and K. Daniilidis.  Sparseness meets deepness: 3d human pose estimation from monocular video.  In CVPR, June 2016.


††thanks: Corresponding author

# Multi-modal transducer-waveguide construct coupled to a medical needle

Yohann Le Bourlout  Medical Ultrasonics Laboratory (MEDUSA), Department of
Neuroscience and Biomedical Engineering, Aalto University, Rakentajanaukio 2,
Espoo, 02150, Finland  Gösta Ehnholm  Medical Ultrasonics Laboratory (MEDUSA),
Department of Neuroscience and Biomedical Engineering, Aalto University,
Rakentajanaukio 2, Espoo, 02150, Finland  Heikki J. Nieminen
<EMAIL_ADDRESS>Medical Ultrasonics Laboratory (MEDUSA), Department
of Neuroscience and Biomedical Engineering, Aalto University, Rakentajanaukio
2, Espoo, 02150, Finland

(September 3, 2024)

††preprint: Yohann Le Bourlout, JASA

## Abstract

Annually, more than 16 billion medical needles are consumed worldwide.
However, the functions of the medical needle are still limited to cutting and
delivering or drawing material through the needle to a target site. Ultrasound
combined with hypodermic needle could potentially add value to many medical
applications such as pain reduction, adding precision, deflection reduction in
tissues and even improve tissue collection. In this study we introduce a
waveguide construct enabling an efficient way to convert a longitudinal wave
mode to flexural mode and to couple the converted wave mode to a conventional
medial needle, while maintaining high electric-to-acoustic power efficiency.
The structural optimization of the waveguide was realized _in silico_ using
the finite element method followed by prototyping the construct and
experimental characterization of the prototypes. The experiments at 30 kHz
demonstrate flexural tip displacement up to 200 µm, at low electrical power
consumption (under 5 W), with up to 69% of efficiency. The high electric-to-
acoustic efficiency and small size of the transducer would facilitate design
of medical needle and biopsy devices, potentially enabling portability,
batterization and high patient safety with low electric powers.

##  I Introduction

The hypodermic needle is one of the most commonly used instruments in
healthcare, 16 billion needles are consumed annually Hauri2004TheSettings ;
Hutin2003UseEstimates . They are frequently used as standard procedure for
tissue sampling Amedee2001Fine-NeedleBiopsy ; Renshaw2007ComparisonBiopsy  or
injection of fluids Cardone2002JointInjection ; Wynaden2006BestSetting ;
Aiello2004EvolvingInjections ; Guarda-Nardini2008ArthrocentesisTechnique ;
Lenz2017AAdministration . However, only limited development has taken place in
their design for several decadesKirkup1998SurgicalInstruments. ;
Ball2004LocalNeedles . This has caused complications in some applications,
_e.g._ for targeting a specific area for anesthesia or targeting tumor in
biopsy Pritzker2019NeedleCare  (up to 34% failure rate in breast biopsy using
Fine Needle Aspiration (FNA) Pisano1998RateStudy  and up to 21% in Thyroid
Gharib1993Fine-NeedleBiopsies ). Moreover, approximatively 10% of the
population fears needlesNir2003FEARASSOCIATIONS , which can lead to avoiding
accessing healthcareHamilton1995NeedleDiagnosis .

During the operation of the needle, the operator can translate the needle
back-and-forth in a piston movement Kreula1990EffectBiopsy. , it can be bent
to guide the needle and the needle can be moved by rotating it around its
center axis by the operator’s hand Han2013StudyBiopsy . A more modern way to
introduce needle movement is to provide mechanical actuation of the needle to
produce acoustic waves, essentially providing micron-scale motions within the
needle at a rate much greater than possible to produce manually by a human
operator. Actuation of the needle tip by ultrasound has recently been
demonstrated to be useful for many medical purposes: to (i) decrease the
penetration force Liao2013ReducedStudy ; Liao2014PerformanceMimics ;
Gui2014AutomaticDevice ; Yang2004MicroneedleActuation  for an easier and more
accurate procedure, (ii) causing less pain by reducing the penetration
resistance Lehtinen1979PenetrationNeedles ;
Egekvist1999RegionalUltrasonography , and (iii) provide improved guidance of
the biopsy needle under ultrasound imaging Kuang2016ModellingBiopsy ;
Sadiq2014EnhancedNeedle ; Jiang2020LocalizationImaging .

We have recently developed an ultrasound-enhanced fine needle aspiration
biopsy (USe-FNAB) method, which allows introducing flexural ultrasound waves
to the tip of a medical needle to improve the tissue yield in biopsy
Perra2021UltrasonicYield . In USeFNAB, a piezoelectric actuator is employed to
produce longitudinal waves, which are converted into flexural waves
essentially producing mechanical movement of the needle tip. The system
previously studied was relatively large in size including a standalone
function generator and a power amplifier. However, clinically a small and a
fully portable device could be in some circumstances more preferable,
requiring high electrical-to-power-conversion efficiency.

In this study, we aimed to develop a new acoustic waveguide construct with the
intent to optimize the conversion efficiency from electrical to acoustic power
followed by validation of the technical performance.

##  II Material and methods

Using the Finite Element (FE) method, we chose to first model the Langevin
transducer that was selected as the sound source. This was followed by
optimization of the geometry of the full system comprising a piezoelectric
transducer, a waveguide and a biopsy needle. Based on the simulations, we then
built custom-made prototypes, followed by characterizing their performance and
comparing the performance with simulations.

We studied two waveguide constructs: (i) a system with a transducer center
axis placed perpendicular to the needle (Figure 2) and (ii) a system with a
transducer aligned at an offset and parallel to the needle (Figure 3), the
needle being a hypodermic needle (21G, length 120 mm, model: 466564/3, 100
Sterican, B Braun, Melsungen, Germany). The purpose of the waveguide was to
convert the longitudinal wave mode in the transducer to a flexural wave mode
in the needle, as well as to match the mechanical impedance of the transducer
to the mechanical impedance of the needle. The mass distribution, shape and
positioning of the structure are critical to achieve optimum performance of
the system. Therefore, to understand the wave propagation and optimize the
system, an acoustic simulation was conducted _in silico_ employing a FE method
(COMSOL Multiphysics 5.5., Inc., Burlington, Massachusetts, United States). In
the following, the steps realized during this optimization are elaborated.

###  II.1 Numerical simulations

For optimization of the arrangement, FE method was used (COMSOL Multiphysics
5.5). Electrostatics module was applied to the piezo elements using the
stress-charge constitutive relation, associating the electrical field applied
to the piezo elements with structural phenomena. Elastodynamics module (”Solid
Mechanics” module) was used to solve the equation of motion. To avoid infinite
vibration amplitude when the system is in resonance, a damping ration of the
structural component was set to 0.01. The different parts of the model were
meshed using the swept technique for the piezo stack and the free tetrahedral
technique for the rest of the geometry. The mesh size was implemented to
include at least 20 elements per wavelength. The simulations were conducted in
two steps. First, a stationary study was performed to calculate the
deformation of the structure, when pretension was applied to the piezo stack
as elaborated below. The second step was a frequency domain perturbation
study, allowing to simulate the small oscillations of the system.

The parameters of the different materials used in the models are given in the
table 1.

###  II.2 Simulation of the Langevin transducer

The modelled Langevin sandwich transducer was composed of twenty-two ring-type
piezo elements (PIC151; dimensions = ⌀out $\times$ ⌀in $\times$ thickness = 10
$\times$ 5 $\times$ 0.5 mm) stacked and fastened by a bolt. Each element was
arranged so that the direction of the electric fields through the piezo were
opposite to each other in adjacent elements. Electrically, the elements were
all connected in parallel: In this way the required driving amplitude voltage
was decreased by a factor of 22 compared to using a single element having a
thickness of the full stack. At each end of the stack, an extra inactivated
piezo element (PIC151, dimensions = ⌀out $\times$ ⌀in $\times$ thickness = 10
$\times$ 5 $\times$ 1 mm) was added to electrically isolate the transducer
from the rest of the system for safety. To simulate the glue used by the
manufacturer to hold the stack together, two layers of epoxy were added
(epoxy, dimensions = ⌀ $\times$ H = 10 $\times$ 0.2 mm). Two masses of brass
were added on each side of the piezo stack connected with a bolt, detailed
below, to shift the resonance frequency of the piezo stack to the targeted
resonance of the entire construct (approx. 30kHz). The front mass of the
transducer was cylindrical (brass, dimensions = H $\times$ ⌀ = 14 $\times$ 10
mm). The back mass was shaped to facilitate the fixation of the stack inside
the enclosure. An M4 screw (brass, H $\times$ ⌀ = 35 $\times$ 4 mm) was
connected to a cylindrical front mass (brass, H $\times$ ⌀= 3 $\times$ 10 mm)
and to a cube-shaped back mass (brass, H $\times$ W $\times$ L = 10 $\times$
10 $\times$ 10 mm) connected to a rectangular cuboid (brass, H $\times$ W
$\times$ L = 5 $\times$ 10 $\times$ 12 mm). The M4 screw, cylindrical, cube
and the rectangular cuboid were all machined from one piece (_c.f._ Figure
1.A).
To fasten all the elements mentioned above to create a stack, the M4 bolt
linked the back mass and the top mass in order to apply a pretension force
(3.1 kN) over the stack. The bolt connects both ends through the top mass via
a thread allowing to expose 3 mm of the bolt, which was used to connect the
Langevin transducer to the waveguide having a threaded hole. Finally, after a
pre-study calculating the bolt pretension effect over the system, a frequency
domain with perturbation simulation from 20 to 40 kHz was conducted to find
the resonance frequency of the system.

###  II.3 Simulations of the waveguides

Connecting the needle directly to the transducer without impedance matching
would lead to compromised transmission of acoustic power from the transducer
to the needle. The transducer stack is vibrating with a longitudinal mode with
very small amplitude, and has mechanically a high impedance. On the other
hand, the needle should vibrate with a large amplitude approx. perpendicularly
to its center axis and is quite flexible, _i.e._ it has a low mechanical
impedance. Therefore, mechanical impedance matching is needed.

To avoid reflection of power back to the transducer within the waveguide, a
mechanical impedance matching is required. This structure should also convert
the longitudinal mode from the transducer to a flexural mode of the needle
tip, while matching gradually the transducer mechanical impedance to the
needle. We investigated two waveguide geometries: a linear waveguide and a
right-angled waveguide.

The linear waveguide was made of EOS Stainless steel 316L and included a
rectangular resonator in the proximal end (dimensions = L $\times$ W $\times$
H = variable length (8:0.1:100) $\times$ 12 $\times$ 3 mm) (Figure 2). From
here the waveguide extrudes as a tapered section with the width (W) and the
height (H) decreasing in width exponentially along the length towards the
distal end (L = variable length (5:0.1:50)mm, dimensions starting point = W
$\times$ H = 12 $\times$ 3 mm, dimensions ending point = W $\times$ H = 0.9
$\times$ 0.5 mm). The tapering increases the amplitude of the acoustic signal
by mechanical impedance matching: The constant flow of power will propagate as
a wave with increasing amplitude into more flexible regions at smaller cross
sections of the waveguide minimizing reflections. The waveguide also included
a groove along the length of the waveguide to provide a coupling area for the
needle cannula (dimensions = ⌀in $\times$ ⌀out: = 0.514 $\times$ 0.8 mm).

The second right-angled waveguide design included a linear portion with a 90°
angle implemented in the beginning of the tapered portion (Figure 3). The
rationale for this geometry is to improve the ergonomics by aligning the
center axes of the transducer and the needle to be parallel.
The length dimensions of the waveguides were simulated with the resonance
frequency obtained from the transducer simulation.

###  II.4 Needle position

After optimization of the waveguide, the proximal part of the needle was
included into the simulation to understand its effect on the transducer and
find an optimal position. The distance from the tip of the waveguide to the
needle tip was varied from 8 mm to 60 mm at 100 µm increments. The simulated
needle follows the specification of a standard 21-gauge (21G) hypodermic
needle (dimensions = ⌀in $\times$ ⌀out = 0.514 $\times$ 0.8 mm) made of
stainless steel 316L. The needle tip, based on standard hypodermic needle tip,
is composed of two bevels which have opening angles of 17° and 11°,
respectively, observed as a side projection to _y_ _z_ plane (Figure 2.C).

###  II.5 Impedance characteristics of the experimental system and simulation

After building the experimental systems, the agreement between the
experimental system and the simulation was verified in terms of impedance
characteristics as follows: The Langevin transducer was coupled to the
waveguide, which was a 3D printed EOS stainless steel 316L object (3D Formtech
Oy, Jyvaskyla, Finland). The needle 21G x 120 mm (model: 466564/3, 100
STERICAN, B Braun, Melsungen, Germany), placed in the groove, was soldered
onto the last 20 mm of the waveguide’s distal end. The electrical impedance of
the arrangement, was obtained from COMSOL using frequency domain with
perturbation simulation (20-40 kHz with steps of 50 Hz). The impedance
measurement was conducted using a Digilent Analog Discovery 2 combined with
the WaveForms software (Digilent, Inc., Henley Court Pullman, 262 WA, United
States).

###  II.6 Quantifying the needle tip displacement

To observe the needle tip during sonication we used a high-speed camera
(model: Phantom V1612, Vision Research, Wayne, NJ, United States) in
conjunction with a macro lens (model: Canon MP-E 65 mm f / 2.8 1-5x Macro
Photo, Canon Inc., Ota, Tokyo, Japan). The camera was placed in front of an
acrylic chamber (external dimensions = L $\times$ W $\times$ H = 21 $\times$
14 $\times$ 15 mm, wall thickness = 5 mm) filled with air or deionized water
degassed to 6 mg/L at ambient temperature (22.3°C). The needle was immersed to
a depth of 30 mm from the water surface. The system was driven by a custom-
made amplifier coupled with a custom-made step-up transformer (11:20 turn
ratio) at 30.4 kHz, for 10 cycles with an incoming instantaneous power of 5 W.
High speed camera videos were taken in order to acquire the needle tip
displacement during sonication (sample rate = 304 000 frames per second,
exposure time = 0.8µs, resolution 128 x 128 pixels, angular magnification =
5x). The videos were then analyzed in MATLAB r2020b using an algorithm capable
of tracking the needle position using cross-correlation frequency technique
based on Perra et al. Perra2021UltrasonicYield  work.

###  II.7 Calorimetry

To define the total acoustic power emitted from the needle, a calorimetric
assessment was done. The calorimeter (material= thermal wall insulator, wall
thickness= 3 cm) was filled with 25 mL of deionized water (at room
temperature: 22°C) mixed with a magnetic stirrer (C-MAG HS series C-MAG HS 4
model, IKA, Staufen, Germany). The needle tip was then immersed to a depth of
30 mm from the water surface. The system was driven at 30.4 kHz using a burst
signal with 55% activation time (used for a better control at the needle tip)
during 10 seconds. The forward and reflected powers were measured and time-
averaged over a burst period: The time-averaged consumed power is equal to
their difference. The water temperature change was recorded using a PT 3000
thermometer probe for 10 seconds before and 10 seconds after the sonication.
After averaging the temperature data, the difference of the two temperatures
gave the heat deposed into a known volume of distilled water for a certain
incident energy. The analysis of the data was implemented using a custom-made
MATLAB r2020b software.

##  III Results

###  III.1 Langevin transducer

The first analysis aimed to define the resonance frequency of the transducer.
The computational results indicate a peak vibration amplitude around 27.4 kHz
as shown in Figure 1.

Figure 1: A. Cross-section of the Langevin transducer construct with twenty-
two rings arranged so that their electric fields direction through the piezo
were opposite with adjacents piezo elements. At each end of the piezo stack,
an extra non-activated piezo element was placed to isolate the transducer from
the rest of the system for safety. At the extremities of the transducer, two
masses made of brass were included. An M4 bolt connecting the two masses were
applied to prestress the piezo stack. B. Absolute displacement measured from a
frequency domain simulation using COMSOL multiphysics 5.5 in _y_ direction
with 1 W driving power.

###  III.2 Optimization of the resonator, linear waveguide and needle length

Figure 2: A. Displacement in _y_ and _z_ direction dependent of the resonator
length variation. B. Displacement in _y_ and _z_ direction dependent of the
amplifier length variation. C. Displacement in _y_ and _z_ direction dependent
of the needle length variation

The purpose of this simulation was to find a geometry of the waveguide, which
would provide a great displacement in _y_ direction while avoiding ‘spurious’
displacement in other directions due to different modes (As shown in the row A
of the normalized _z_ displacement graph around 82 mm). In Figure 2.A, there
is a clear recurrence, every 23 mm of the greater values which corresponds to
approx. half a wavelength. The selected dimension was 23 mm. The same applies
in Figure 2. B and C, every half a wavelength, an optimum in _y_ displacement
appearing. 39.5 mm was selected as the length of the waveguide. The
displacement is greatly amplified by the tapered shape compared to the simple
resonator. The needle displacement is also, greater than the amplifier
construct. This is due to the tapered proximal end of the needle which behaves
as an amplifier. The distance distal part of the needle to distal part of the
waveguide was 58 mm.

###  III.3 Optimization of the resonator, right-angled waveguide and needle
length

Figure 3: A. Displacement in _y_ and _z_ direction dependent of the resonator
length variation. B. Displacement in _y_ and _z_ direction for the right-
angled design dependent of the amplifier length variation. C. Displacement in
_y_ and _z_ direction for the right-angled design dependent of the needle
length variation.

The results show in the Figure 3, using the right-angled waveguide suggest a
similar functioning as the linear waveguide. For the optimization of the
amplifier, a resonator of half a wavelength was selected (23 mm). After the
simulation of the amplifier, a length of 39.5 mm was selected for the
compromise between length of the waveguide (long enough to reduce the
reflection of the tapering) and the useful needle length (greater the
waveguide length is, smaller is the proximal size of the needle). The proximal
length of the needle selected for this model was 58 mm.

###  III.4 Validation model - system

Figure 4: A. Comparison between the simulated impedance of the transducer
construct and the measured impedance of the custom-designed transducer. B.
Comparison of the whole system impedance including the transducer, linear
waveguide and needle simulation and experimental device with linear waveguide.
C. Comparison of the whole system with right-angled waveguide and needle
simulation and experimental device with right-angled waveguide.

To compare the difference between the simulations and our prototypes, an
impedance measurement with bare stack transducer and the whole systems were
compared with the respective simulated impedances in frequency range 20 - 40
kHz (Figure 4). The simulation curve and the measured peaks exhibit similar
behavior. The main impedance peak from the transducer is at 30.3 kHz, while
the linear waveguide system seems to be at 29.7 kHz and the right-angled one
appears to be at 29.1 kHz, which remains in the desired frequency range.

###  III.5 Experiment arrangement

Figure 5: Pictures of the two concepts. A. Shows the top and side view of the
device with the linear waveguide, while B. shows the device with right-angled
waveguide. Here the pictures exhibit the more compact size of the right-angled
waveguide.

Figure 5 shows the different prototypes. Figure 6 compares the waveguide
structures and shows that both of the designs are able to create flexural
waves at the needle tip (Figure 6.A-F). Both systems are behaving similarly
(Figure 6.C-F), a needle tip displacement up to 200 µm and 160µm, in air and
water, respectively. Cavitations were observable with both waveguides from the
high speed camera footage in water. Nevertheless, it seems that a slightly
greater electrical-to-acoustical power -efficiency is obtained with the linear
waveguide (Figure 6.G,H): 69% efficiency was recorded with the linear
waveguide compared to an efficiency of 62% with the right-angled waveguide.

Figure 6: Frequency domain simulation (COMSOL multiphysics 5.5) of the linear
(A) and right-angled (B) waveguides (an exaggerated deformation has been used
for better visualization), their experimental displacement of the needle in
air (C,D) or water (E,F) (30.4kHz, 10 cycles, instantaneous electrical power
5W, in blue background) recorded with a high speed camera (304000 fps) and
their needle tip heating deposition in water measurement (C,F) respectively.

##  IV Discussions

The results of the _in silico_ structural optimization demonstrated the
displacements at the needle tip were significantly dependent on the dimensions
of the waveguide design and needle-to-waveguide coupling point as expected,
greatest displacements obtained when the system was in resonance. In the
simulations, an increase in displacement was observed at the end of each
stage, _e.g._ transducer, waveguide resonator or waveguide tapering, to the
end of next system stage, i.e. waveguide resonator, waveguide tapering or
needle, respectively. Importantly, the experimentally demonstrated power-
efficiency was good, 69% with the linear waveguide, 62% with the right-angled
waveguide. Several factors may contribute to the power-efficiency, as outlined
in the following.

The extrusion of the linear and rectangular segment of the resonator extends
symmetrically from the transducer in opposite directions. Symmetricity
minimizes the generation of bending modes in the transducer at the coupling
point to the resonator, which would be more likely with an asymmetric
resonator. Bending waves in the transducer would yield in tension and
compression in different parts of the piezo stack simultaneously, compromising
the performance, which is now eliminated. Furthermore, the exponential tapered
shape of the waveguide gives an amplification of the signal by having
approximately the same energy in the smaller cross-section, while minimizing
the reflection from the boundaries towards the source. The geometry,
therefore, provides geometric amplification at the same time with successful
impedance matching from higher mechanical impedance of the transducer to the
lower mechanical impedance of the needle. The needle bevel provides further an
additional state of displacement amplification Perra2021UltrasonicYield .

It is notable that the studied waveguides converted the longitudinal waves
emitted by the transducer to a flexural displacement of the needle tip. The
experiments with needle in air provided large displacement amplitudes up to
200 µm and 160 µm in water. The generation of cavitation was observed using
both the linear and right-angled waveguides. The practical relevance of the
potential to generate large displacements is selecting the relevant magnitude
needed for a specific application in question. _E.g._ gentle actuation of
tissue by the needle or adjacent micro-bubbles could be used to improve tissue
yield in a biopsy application, whereas generation of large amplitudes and
violent cavitation clouds could be useful in an intervention application such
as histotripsy of a tumor.

Several reports Liao2013ReducedStudy ; Liao2014PerformanceMimics  have shown
that the actuation of the needle tip by ultrasound might be useful for many
medical purposes such as pain reduction or improving spatial accuracy of
needle insertion and avoiding the tissue deformation. Some devices have met
these expectations in terms of needle tip actuation, but none have been small
enough to be adopted by doctors. An initial objective of the project was to
create an ultrasonic actuated needle with high efficiency system. This could
permit batterization and miniaturization if developed to its full potential.
The instantaneous electrical input power into the transducer was below 5 W and
the displacements shown here are sufficient in the biopsy application
Perra2021UltrasonicYield , which could be feasible to provide in a
miniaturized fashion.

The system could be further improved with a resonance frequency tracking,
which was not implemented in this study. The loading of the needle affects the
resonance frequency. Increasing the load of the needle leads to a decrease of
the resistive power and an increase of the reactive power, which would
decrease the efficiency of the system. The frequency tracking for finding the
resonance frequency could thereby enable minimizing reactance and thereby
maintaining maximal needle tip action despite the mechanical environment of
the needle. The mechanical environment would change _e.g._ with the needle
penetration depth in tissue or with different mechanical properties of the
tissue. In our experiment, the effect of the water on the frequency change was
not investigated, but stiffer media could have a greater effect on the
frequency, leading to a non-optimal system.

##  V Conclusion

Ultrasound in medicine is a well known technology and is widely used for
different kind of therapeutic procedures. However, mechanical ultrasound in
healthcare is not widespread. Since the needle is a common tool in healthcare
but with clear limitations, this paper showed a novel way of improving
needles. By actuating the needle with a flexural wave, we might be able to
decrease the penetration force needed during the insertion of the needle or
improve the tissue collection during biopsies. Due to the efficiency of the
system, it suggests that the design could lead to a batterization and,
eventually, enhancing the miniaturization of the system.

##  VI Acknowledgment

We thank Mr. Emanuele Perra, M.Sc., Mr. Matti Mikkola, M.Sc., Mr. Saif Bunni,
M.Sc., Mr. Jouni Rantanen, M.Sc. and Dr. Nick Hayward, MBBS, Ph.D. for
constructive discussions related to the topic. Business Finland and Academy of
Finland are acknowledged for financial support (Business Finland: grant
5607/31/2018; Academy of Finland: grants 314286, 311586 and 335799). We thank
also Alto University, the Department of Neuroscience and Biomedical
Engineering and MEDUSA laboratory for the facilities and equipment.

##  VII Conflict of interest

Heikki Nieminen has stock ownership in Swan Cytologics Inc., Toronto, ON,
Canada and is an inventor within the patent application WO2018000102A1. Yohann
Le Bourlout, Gösta Ehnholm and Heikki Nieminen are inventors within the patent
application PCT/FI2020/050345.

##  VIII Contributions

All authors contributed to the design of the study, writing or reviewing the
manuscript and have approved it. Yohann Le Bourlout have produced the data and
conducted the data analysis. Gösta Ehnholm have assembled the transducer
construct and created the electronics for the system.

##  IX Data availability

The data is available upon request.

##  X Code availability

The codes are available upon request.

## References

  * (1) A. M. Hauri, G. L. Armstrong, and Y. J. Hutin, “The global burden of disease attributable to contaminated injections given in health care settings,” International Journal of STD and AIDS 15(1), 7–16 (2004) doi: 10.1258/095646204322637182.
  * (2) Y. J. Hutin, A. M. Hauri, and G. L. Armstrong, “Use of injections in healthcare settings worldwide, 2000: literature review and regional estimates,” BMJ 327(7423), 1075 (2003) doi: 10.1136/BMJ.327.7423.1075.
  * (3) R. G. Amedee and N. R. Dhurandhar, “Fine-Needle Aspiration Biopsy,” The Laryngoscope 111(9), 1551–1557 (2001) doi: 10.1097/00005537-200109000-00011.
  * (4) A. A. Renshaw and N. Pinnar, “Comparison of Thyroid Fine-Needle Aspiration and Core Needle Biopsy,” American Journal of Clinical Pathology 128(3), 370–374 (2007) doi: 10.1309/07TL3V58337TXHMC.
  * (5) D. A. Cardone and A. F. Tallia, “Joint and Soft Tissue Injection,” American Family Physician 66(2), 283 (2002).
  * (6) D. Wynaden, I. Landsborough, S. McGowan, Z. Baigmohamad, M. Finn, and D. Pennebaker, “Best practice guidelines for the administration of intramuscular injections in the mental health setting,” International Journal of Mental Health Nursing 15(3), 195–200 (2006) doi: 10.1111/J.1447-0349.2006.00423.X.
  * (7) L. P. Aiello, A. J. Brucker, S. Chang, E. T. Cunningham, D. J. D’Amico, H. W. Flynn, L. R. Grillone, S. Hutcherson, J. M. Liebmann, T. P. O’Brien, I. U. Scott, R. F. Spaide, C. Ta, and M. T. Trese, “Evolving guidelines for intravitreous injections,” Retina (Philadelphia, Pa.) 24(5 Suppl) (2004) doi: 10.1097/00006982-200410001-00002.
  * (8) L. Guarda-Nardini, D. Manfredini, and G. Ferronato, “Arthrocentesis of the temporomandibular joint: a proposal for a single-needle technique,” Oral Surgery, Oral Medicine, Oral Pathology, Oral Radiology, and Endodontology 106(4), 483–486 (2008) doi: 10.1016/J.TRIPLEO.2007.12.006.
  * (9) J. R. Lenz, D. D. Degnan, J. B. Hertig, and J. G. Stevenson, “A Review of Best Practices for Intravenous Push Medication Administration,” Journal of Infusion Nursing 40(6), 354–358 (2017) doi: 10.1097/NAN.0000000000000247.
  * (10) J. Kirkup, “Surgical history. The history and evolution of surgical instruments. VIII. Catheters, hollow needles and other tubular instruments.,” Annals of The Royal College of Surgeons of England 80(2), 81 (1998).
  * (11) C. Ball and R. Westhorpe, “Local anaesthesia—the continuing evolution of spinal needles,” Anaesthesia and Intensive Care 32(1), 3 (2004).
  * (12) K. P. Pritzker and H. J. Nieminen, “Needle Biopsy Adequacy in the Era of Precision Medicine and Value-Based Health Care,” Archives of Pathology & Laboratory Medicine 143(11), 1399–1415 (2019) doi: 10.5858/ARPA.2018-0463-RA.
  * (13) E. D. Pisano, L. L. Fajardo, J. Tsimikas, N. Sneige, W. J. Frable, C. A. Gatsonis, W. P. Evans, I. Tocino, and B. J. McNeil, “Rate of insufficient samples for fine-needle aspiration for nonpalpable breast lesions in a multicenter clinical trial: The Radiologic Diagnostic Oncology Group 5 Study,” Cancer: Interdisciplinary International Journal of the American Cancer Society 82(4), 679–688 (1998).
  * (14) H. Gharib, J. R. Goellner, and D. A. Johnson, “Fine-Needle Aspiration Cytology of the Thyroid: A 12-Year Experience With 11,000 Biopsies,” Clinics in Laboratory Medicine 13(3), 699–709 (1993) doi: 10.1016/S0272-2712(18)30434-7.
  * (15) Y. Nir, A. Paz, E. Sabo, and I. Potasman, “FEAR OF INJECTIONS IN YOUNG ADULTS: PREVALENCE AND ASSOCIATIONS,” The American Journal of Tropical Medicine and Hygiene 68(3), 341–344 (2003) doi: 10.4269/AJTMH.2003.68.341.
  * (16) J. G. Hamilton, “Needle phobia: a neglected diagnosis,” Journal of Family Practice 41(2), 169–182 (1995).
  * (17) J. Kreula, P. Virkkunen, and S. Bondestam, “Effect of suction on specimen size in fine-needle aspiration biopsy.,” Investigative Radiology 25(11), 1175–1181 (1990) doi: 10.1097/00004424-199011000-00005.
  * (18) P. Han and K. Ehmann, “Study of the effect of cannula rotation on tissue cutting for needle biopsy,” Medical Engineering & Physics 35(11), 1584–1590 (2013) doi: 10.1016/J.MEDENGPHY.2013.05.001.
  * (19) X. Liao, M. Sadiq, G. Corner, S. Cochran, and Z. Huang, “Reduced penetration force through ultrasound activation of a standard needle: An experimental and computational study,” IEEE International Ultrasonics Symposium, IUS 1436–1439 (2013) doi: 10.1109/ULTSYM.2013.0364.
  * (20) X. Liao, M. Sadiq, Y. Kuang, G. Corner, S. Cochran, and Z. Huang, “Performance optimization of ultrasonic needle actuating device for insertion operation into tissue mimics,” IEEE International Ultrasonics Symposium, IUS 823–826 (2014) doi: 10.1109/ULTSYM.2014.0202.
  * (21) L. Gui, R. Xiao, X. Liao, Y. Kuang, M. Sadiq, Z. Chunnian, S. Cochran, and Z. Huang, “Automatic frequency tracking system for needle actuating device,” IEEE International Ultrasonics Symposium, IUS 815–818 (2014) doi: 10.1109/ULTSYM.2014.0200.
  * (22) M. Yang and J. D. Zahn, “Microneedle Insertion Force Reduction Using Vibratory Actuation,” Biomedical Microdevices 2004 6:3 6(3), 177–182 (2004) doi: 10.1023/B:BMMD.0000042046.07678.2E.
  * (23) R. Lehtinen and E. Oksala, “Penetration of disposable needles,” International Journal of Oral Surgery 8(2), 145–148 (1979) doi: 10.1016/S0300-9785(79)80010-1.
  * (24) H. Egekvist, P. Bjerring, and L. Arendt-Nielsen, “Regional variations in pain to controlled mechanical skin traumas from automatic needle insertions and relations to ultrasonography,” Skin Research and Technology 5(4), 247–254 (1999) doi: 10.1111/J.1600-0846.1999.TB00137.X.
  * (25) Y. Kuang, A. Hilgers, M. Sadiq, S. Cochran, G. Corner, and Z. Huang, “Modelling and characterisation of a ultrasound-actuated needle for improved visibility in ultrasound-guided regional anaesthesia and tissue biopsy,” Ultrasonics 69, 38–46 (2016) doi: 10.1016/J.ULTRAS.2016.02.018.
  * (26) M. R. Sadiq, S. Cochran, X. Liao, and Z. Huang, “Enhanced US-guided needle intervention through ultrasound actuation of a standard needle,” IEEE International Ultrasonics Symposium, IUS 827–830 (2014) doi: 10.1109/ULTSYM.2014.0203.
  * (27) T. Jiang, G. McLeod, Z. Huang, X. Zhu, Y. Jiao, X. Li, Z. Shen, and Y. Cui, “Localization Accuracy of Ultrasound-Actuated Needle with Color Doppler Imaging,” Diagnostics 2020, Vol. 10, Page 1020 10(12), 1020 (2020) doi: 10.3390/DIAGNOSTICS10121020.
  * (28) E. Perra, E. Lampsijärvi, G. Barreto, M. Arif, T. Puranen, E. Hæggström, K. P. Pritzker, and H. J. Nieminen, “Ultrasonic actuation of a fine-needle improves biopsy yield,” Scientific Reports 2021 11:1 11(1), 1–15 (2021) doi: 10.1038/s41598-021-87303-x.

##  XI Annexes

Table 1: Material parameters used for the FE simulations Material | Density [kg/m3] | Young’s modulus [GPA] | Poisson’s ratio | Isotropy
---|---|---|---|---
Brass | 8730 | 82 | 0.31 | isotropic
Epoxy | 1250 | 3.5 | 0.3 | isotropic
PIC151 | 7760 | Manufacturer data sheet | Manufacturer data sheet | anisotropic
EOS Stainless steel 316L | 8000 | 193 | 0.265 | isotropic
Stainless steel 316L | 8070 | 205 | 0.275 | isotropic