Pictorial Structures

A pictorial structure represents an object as a collection of parts, each with its own model of local appearance, connected in pairs by spring-like costs that penalize departures from the parts’ expected relative positions. Matching the model to an image means choosing a location for every part at once, trading how well each part fits the image against how much the springs are stretched. The idea goes back to Fischler and Elschlager (1973). Felzenszwalb and Huttenlocher (2005) made it practical: when the connections form a tree and the springs have a quadratic form, the globally best match can be found exactly, in time roughly linear in the number of candidate locations per part, by dynamic programming with generalized distance transforms. Their framework also gave the model a probabilistic reading, a way to sample matches and a way to learn the parameters from labelled examples. The family includes the constellation models of the Caltech and Oxford groups, the classical body models of human pose estimation, the star-shaped deformable part models used for detection, Yang and Ramanan’s flexible mixtures of parts, and several deep pose networks that kept a part graph and its inference on top of convolutional features. The terms part-based model, deformable part model and, for the probabilistic version, tree-structured Markov random field overlap with it.

Definition

A pictorial structure model is an undirected graph $G = (V, E)$ whose vertices $v_1, \dots, v_n$ are the parts of an object and whose edges join the pairs of parts that are constrained relative to each other. A configuration $L = (l_1, \dots, l_n)$ gives a location $l_i$ to every part. The location can be just an image position, or it can include the part’s orientation, scale or foreshortening. A match function $m_i(l_i)$ measures how badly part $i$ fits the image at $l_i$, and a deformation function $d_{ij}(l_i, l_j)$ measures how far the pair $(v_i, v_j)$ is from its ideal relative placement. The best match is the configuration of lowest total cost. Felzenszwalb and Huttenlocher stress that the framework does not fix how appearance is modelled or what form the connections take, and that matching makes no early decision about where any single part is: all parts are placed together.

Intuition

Picture the parts as rigid templates laid on the image and joined by springs. Each template is pulled toward the places where the image looks like it, and the springs pull the templates toward a plausible arrangement. A forearm on its own is hard to find, since any pair of parallel edges looks like one, but a forearm attached to an elbow attached to an upper arm attached to a torso is far less ambiguous. Felzenszwalb and Huttenlocher make this point about the human body: simple part appearance models together with the relations between parts give enough context to find the body as a whole, even where generic parts such as a lower leg would be hard to find alone.

The difficulty is search. With $h$ candidate locations per part, there are $h^n$ configurations; with a million locations per part, even a pair of parts has $10^{12}$ joint placements. The two restrictions that make pictorial structures practical attack this from two sides. A tree of connections lets the search be split into a chain of small problems, one per edge, each costing $O(h^2)$. A spring with a quadratic cost lets each of those problems be solved in $O(h)$, because finding the best child location for every parent location turns out to be a distance transform of the child’s cost map.

Formal Definition

The energy

Felzenszwalb and Huttenlocher write the matching problem of Fischler and Elschlager as

$$L^* = \arg\min_{L} \left( \sum_{i=1}^{n} m_i(l_i) + \sum_{(v_i, v_j) \in E} d_{ij}(l_i, l_j) \right).$$

The deformation costs usually depend only on the relative position of the two parts, which makes the model invariant to certain global transformations, such as translating the whole object. For arbitrary graphs and arbitrary costs this minimization is NP-hard; the authors note that it has the same form as the Markov random field energies of low-level vision, such as stereo and image restoration, but that the structure of the graph and the space of solutions differ substantially.

The statistical reading

Felzenszwalb and Huttenlocher recast the same energy as a posterior. With model parameters $\theta = (u, E, c)$, where $u$ holds the appearance parameters of every part and $c$ the connection parameters of every edge, the likelihood of the image is approximated by a product over parts, $p(I \mid L, \theta) \propto \prod_i p(I \mid l_i, u_i)$, and the prior over configurations is a tree-structured Markov random field with no preference for absolute position, $p(L \mid \theta) = \prod_{(v_i, v_j) \in E} p(l_i, l_j \mid c_{ij})$. Taking the negative logarithm of the posterior gives the energy above with

$$m_i(l_i) = -\log p(I \mid l_i, u_i), \qquad d_{ij}(l_i, l_j) = -\log p(l_i, l_j \mid c_{ij}),$$

so the best match is the maximum a posteriori (MAP) configuration. Two other questions now have answers: several good matches can be drawn by sampling from the posterior, and all parameters can be fitted to labelled examples by maximum likelihood. The product likelihood is a good approximation only when parts do not overlap; the authors point out that it counts the same image evidence twice when two parts lie on top of each other.

The two restrictions

The efficient algorithms require that $G$ be a tree and that each deformation cost be a Mahalanobis distance between transformed locations,

$$d_{ij}(l_i, l_j) = \big(T_{ij}(l_i) - T_{ji}(l_j)\big)^\top M_{ij}^{-1} \big(T_{ij}(l_i) - T_{ji}(l_j)\big),$$

with $M_{ij}$ diagonal and the transformed locations lying on a grid. In the statistical reading this is a Gaussian prior on the displacement between transformed locations, $p(l_i, l_j \mid c_{ij}) \propto \mathcal{N}\big(T_{ij}(l_i) - T_{ji}(l_j);\, 0, D_{ij}\big)$, with $D_{ij} = M_{ij}/2$ (the similarity and distance measures article defines the Mahalanobis distance). The transformations do a lot of work. For a face model whose part positions are related by a Gaussian with full covariance $\Sigma_{ij} = U_{ij} D_{ij} U_{ij}^\top$, the paper sets $T_{ij}(l_i) = U_{ij}^\top(l_i - s_{ij})$ and $T_{ji}(l_j) = U_{ij}^\top l_j$, which rotates the displacement into axes where the covariance is diagonal. For a limb model, the transformation maps each part to the position of the joint it shares with its neighbour, so that the spring connects the two ends of a joint rather than the two part centres.

Dynamic programming on a tree

Choose any part $v_r$ as the root. For a part $v_j$ with parent $v_i$ and children $C_j$, define the cost of the best placement of $v_j$ and its whole subtree, as a function of the parent’s location:

$$B_j(l_i) = \min_{l_j} \left( m_j(l_j) + d_{ij}(l_i, l_j) + \sum_{v_c \in C_j} B_c(l_j) \right).$$

For a leaf the sum is empty. The functions are computed from the deepest parts up to the root, recording with each $B_j$ the minimizing $l_j$ for every $l_i$. The root’s best location is

$$l_r^* = \arg\min_{l_r} \left( m_r(l_r) + \sum_{v_c \in C_r} B_c(l_r) \right),$$

and the other parts follow by tracing the recorded choices back down the tree. Felzenszwalb and Huttenlocher describe this as a generalization of the Viterbi algorithm, and note the link to belief propagation. The choice of root does not change the result. Each $B_j$ evaluated naively compares every child location with every parent location, so the whole algorithm takes $O(h^2 n)$ time, which they judged too slow when a part has hundreds of thousands or millions of possible locations.

Generalized distance transforms

The classical distance transform of a set of grid points $B$ gives, at every grid location, the distance to the nearest point of $B$. Writing it as $\mathcal{D}_B(x) = \min_y \big(\rho(x, y) + 1_B(y)\big)$, with an indicator that is 0 on $B$ and $\infty$ elsewhere, and replacing the indicator by an arbitrary function $f$ gives the generalized distance transform

$$\mathcal{D}_f(x) = \min_{y} \big(\rho(x, y) + f(y)\big),$$

which finds, for each $x$, a location $y$ that is both close to $x$ and cheap under $f$. With the Mahalanobis form of $d_{ij}$, each message is exactly such a transform evaluated at the transformed parent location, $B_j(l_i) = \mathcal{D}_f\big(T_{ij}(l_i)\big)$, where $f(y) = m_j\big(T_{ji}^{-1}(y)\big) + \sum_{c} B_c\big(T_{ji}^{-1}(y)\big)$ on the range of $T_{ji}$ and $\infty$ elsewhere. The parent’s transformed location need not fall on the grid, and the paper takes the value at the nearest grid point, an approximation it considers small because the transform changes slowly. With a diagonal $M_{ij}$ the whole matching runs in $O(h' n)$ time, where $h'$ is the number of grid locations in the transformed space. A diagonal matrix also makes the distance a sum over coordinates, so a multi-dimensional transform can be computed as a sequence of 1D transforms along each axis, the decomposition described in the 2012 paper below.

For the squared distance in one dimension, $\mathcal{D}_f$ is the lower envelope of the parabolas $(x - q)^2 + f(q)$ rooted at the grid points. Because they all have the same shape and arrive in order, the envelope can be built in one left-to-right pass in $O(n)$ time; the deformable part models article derives the intersection formula and gives the algorithm in the max form used by DPMs. Felzenszwalb and Huttenlocher describe this lower-envelope method in “Distance Transforms of Sampled Functions” (Theory of Computing, 2012), which also frames the transform as a minimum convolution, the operation of grayscale morphology, and lists Viterbi decoding and max-product belief propagation among its uses. The 2005 pictorial structures paper itself computes the transform with a method due to Karzanov, which it cites through Rucklidge’s description.

Sampling and marginals

Replacing minimization by summation turns the recursion into the sum-product algorithm on the tree, which computes the posterior marginal of the root,

$$S_j(l_i) \propto \sum_{l_j} \left( p(I \mid l_j, u_j)\, p(l_i, l_j \mid c_{ij}) \prod_{v_c \in C_j} S_c(l_j) \right),$$

after which a root location is drawn from its marginal and each child is drawn given its parent. The relation between the two algorithms, the paper notes, is that of the Viterbi algorithm to the forward-backward algorithm for hidden Markov models. With the Gaussian connection model, each $S_j$ is a Gaussian convolution in the transformed space, which the authors approximate in linear time with cascaded uniform filters. Once the $S$ functions are computed, each further sample costs little.

Learning

Given images with labelled configurations, the maximum-likelihood estimate splits into independent problems. Each part’s appearance parameters are estimated from that part’s examples alone, and each candidate connection’s parameters from the pairs of locations of its two parts. The tree itself is learned too: scoring each possible connection by the likelihood of the examples under its fitted joint distribution, the best tree is the minimum spanning tree of the complete graph with edge weights $-\log q(v_i, v_j)$, in the manner of Chow and Liu’s tree approximation, which Kruskal’s algorithm finds in $O(n^2 \log n)$ time.

History and Variants

Fischler and Elschlager

Fischler and Elschlager’s 1973 paper appeared in IEEE Transactions on Computers. It could not be read for this article, and what follows is how Felzenszwalb and Huttenlocher (2005) describe it. In their account the 1973 paper introduced the representation, posed matching as the minimization of a combined cost of part fit and spring deformation, and considered only the problem of finding the best match. They add that most earlier work used manually constructed models, and that previous matching methods relied on heuristics or local search that need a good initialization and do not guarantee the optimum. The three shortcomings they list, a hard minimization, many parameters, and the need for more than one good match, set the agenda for their paper.

Felzenszwalb and Huttenlocher

Their CVPR 2000 paper, “Efficient Matching of Pictorial Structures”, gave an exact matching algorithm under a different restriction on the connections, which, as the 2005 paper says, did not allow efficient sampling from the posterior. The 2005 journal paper in the International Journal of Computer Vision, read here in the authors’ copy, used the Mahalanobis form and added sampling and learning, and demonstrated two models.

• Faces. Five parts (the eyes, the nose and the corners of the mouth), located by $(x, y)$ position, with appearance modelled by a Gaussian with diagonal covariance over a 27-dimensional vector of Gaussian derivative filter responses at three scales. Training used twenty labelled images from the Yale face database, augmented by small scalings and rotations, and the learned tree was not specified by hand. MAP matching took under a second per image; the match stayed correct with up to two of the five parts occluded and failed with three. Multiple faces were found by taking the best root location, excluding its neighbourhood and repeating.
• People. Ten rectangular parts (torso, head, two per arm and two per leg) located by position, orientation and foreshortening, on binary images obtained by background subtraction. Each part’s likelihood counts foreground pixels inside the rectangle and in a border around it. Joint springs use Gaussian terms for the joint positions and the foreshortening difference and a von Mises term for the relative angle, which the paper turns into the required Gaussian form by representing each angle as a unit vector. The location grid was $70 \times 70 \times 10 \times 32$, over 1.5 million locations per part, with the 32 angles placed on an $11 \times 11$ grid that makes $h'$ about four times $h$. Because the product likelihood over-counts evidence when parts overlap, the MAP estimate could be wrong although the true configuration had a high posterior; the authors therefore drew 200 samples from the posterior with the likelihood smoothed by raising it to the power $1/T$, with $T = 10$, and picked the one with the best Chamfer distance to the binary image, at about one minute per image.

Constellation models

A parallel line of work at Caltech and Oxford modelled object categories as constellations of parts over sparse features. Felzenszwalb and Huttenlocher (2005) describe the constellation model of Burl, Weber and Perona (1998) as constraining all parts relative to a central coordinate system with a joint Gaussian rather than through pairwise connections, which suits rigid categories but, in their view, cannot represent articulated objects, and as matched with heuristics that need not find the optimum.

Fergus, Perona and Zisserman (CVPR 2003) extended the model of Burl et al. and Weber et al. and learned it without supervision from cluttered images. An entropy-based interest operator first detects up to about 30 regions with their scales. A hypothesis $h$ assigns one detected feature, or none, to each of $P$ parts (typically 3 to 7), and the model scores the image by a likelihood ratio against a background model, factored into appearance (a Gaussian per part with diagonal covariance), shape (a joint Gaussian with full covariance over the part positions after scale normalization), relative scale, and occlusion with feature counts (a table over the $2^P$ occlusion patterns and a Poisson term). Because the shape density couples all parts, the number of hypotheses is $O(N^P)$ for $N$ features, and the authors used A* and other search methods to keep learning and recognition tractable: learning a six- or seven-part model by expectation-maximization from 400 images with 20 to 30 features each took 24 to 36 hours, and recognition 2 to 3 seconds per image. Their task was image classification, deciding whether a category is present; with six parts they reported ROC equal error rates of 92.5% on motorbikes, 96.4% on faces and 90.2% on airplanes, and a recall-precision equal error rate of 88.5% on side views of cars (Figure 3).

The contrast with Felzenszwalb and Huttenlocher’s model is instructive. A constellation model evaluates a few dozen detected features with a fully connected shape model; a tree-structured pictorial structure evaluates every location for every part with a sparse shape model. The first makes a hard decision about where features are before matching, which Felzenszwalb and Huttenlocher’s framework is designed to avoid; the second gives up correlations that the tree does not contain.

Star models and deformable part models

A star, with one central part connected to all others, is the simplest tree. Felzenszwalb and Huttenlocher already name it as a useful structure, and deformable part models adopt it for object detection: a root filter on HOG features, part filters at twice the resolution with quadratic deformation costs, mixtures of components, and a latent SVM that learns everything from bounding boxes. Because the parts interact only through the root, each part’s distance transform is computed once and every root location is scored by a sum of shifted arrays, which makes a DPM a sliding window detector whose window score includes the best part placement. That article covers the model, its training and its relation to convolutional networks.

Pictorial structures revisited

Andriluka, Roth and Schiele (CVPR 2009) argued that a single pictorial structure model could serve pedestrian detection, upper-body pose estimation and full-body pose estimation, if its components were chosen well. They kept Felzenszwalb and Huttenlocher’s body model, with parts located by $(x, y, \theta, s)$ and Gaussian relations expressed in the coordinate system of the joint between two parts, and replaced the appearance model: each part was an AdaBoost classifier on densely sampled shape context descriptors (12 log-polar location bins by 8 orientation bins), trained with bootstrapped hard negatives, whose normalized margin, floored at $\epsilon_0 = 10^{-4}$, was read as the likelihood. Inference used sum-product belief propagation to compute each part’s marginal posterior, with the Gaussian convolutions done in the joint’s coordinate system as in the 2005 paper, on a grid with steps of 1 pixel, 15 degrees and 0.1 in scale. The models had 10 parts for the full body, 6 for the upper body and 8 for pedestrians (no arms, but feet).

On the Iterative Image Parsing dataset of Ramanan, with only 100 training images, their model localized 55.2% of parts against 27.2% for Ramanan’s method, under the part-detection criterion of Ferrari et al. (Table 2); plugging Ramanan’s edge templates into their kinematic model gave 37.5%. On Buffy they reported 73.5% with a front/back view prior (71.3% with the generic prior) against 57.9% for the progressive search-space reduction of Ferrari et al. (Table 1). For pedestrian detection on TUD-Pedestrians, a kinematic tree prior was slightly better than a star prior, and both beat the partISM detector by 4 and 5 points of equal error rate. They attributed their results to the combination of strong discriminatively trained part detectors with a flexible kinematic tree prior.

Flexible mixtures of parts

Yang and Ramanan (CVPR 2011; TPAMI 2013) dropped rotated and foreshortened limb templates. Parts sit at joints, with locations $p_i$ given by pixel position alone, and each part takes one of $T$ discrete types $t_i$ with its own HOG template and its own spring to its parent. The scoring function, given in the human pose estimation article, adds type biases and type co-occurrence terms $b_{ij}^{t_i, t_j}$ to the appearance and spring terms. Inference is the same tree recursion with types added. A part $i$ collects the messages of its children and passes a message to its parent $j$ (here $m$ denotes a message, not Felzenszwalb and Huttenlocher’s match cost):

$$\operatorname{score}_i(t_i, p_i) = b_i^{t_i} + w_i^{t_i} \cdot \phi(I, p_i) + \sum_{k \in \operatorname{kids}(i)} m_k(t_i, p_i),$$ $$m_i(t_j, p_j) = \max_{t_i} \left( b_{ij}^{t_i, t_j} + \max_{p_i} \Big( \operatorname{score}_i(t_i, p_i) + w_{ij}^{t_i, t_j} \cdot \psi(p_i - p_j) \Big) \right),$$

with $\psi(p_i - p_j) = [dx \;\, dx^2 \;\, dy \;\, dy^2]^\top$. Done naively, each message costs $O(L^2 T^2)$ for $L$ locations; the inner maximization is a distance transform for each pair of types, giving $O(L T^2)$ per part, and the authors’ articulated variant, whose springs depend only on the child’s type, needs $O(L T)$. Root scores serve directly as detections after thresholding and non-maximum suppression. The paper names special cases: $T = 1$ is a standard pictorial structure, and allowing only globally consistent types turns the model into a mixture of star models, as in DPMs.

Types are not annotated. The authors derive them by clustering each part’s position relative to its parent with k-means, so that a type corresponds roughly to an orientation, and train all templates, springs and biases together as a structural SVM whose negative constraints come from every root placement in images without people. In the TPAMI version, read in the authors’ copy, diagnostics on the Parse dataset with a 26-part model separate the ingredients (Tables 1 and 2, PCK at $\alpha = 0.1$): the full model with joint training gave 72.9% PCK, training the part templates independently 56.8%, and a star-structured model 56.7%, while switching off the latent updating of types during training gave 73.4%, which the text reads as latent updating not being helpful (the table caption says not crucial). The conference version describes the larger model as having 27 parts, the journal version 26. A caption in the journal version (Fig. 6) gives the number of global type combinations for 14 parts with 4 types each as about $2 \times 10^7$, although $4^{14}$ is about $2.7 \times 10^8$. The model took about a second per image; its benchmark results and the PCK and APK measures it introduced are covered in the human pose estimation article.

Beyond trees

Yang and Ramanan’s journal paper states the main weakness of trees plainly: given a torso, the two legs are placed independently and often respond to the same image region. This is double counting, the same over-counting that Felzenszwalb and Huttenlocher handled by sampling. They list the remedies tried in the literature, each at a cost: loopy models that need approximate inference (sampling, loopy belief propagation, iterative approximations), branch and bound with tree-based bounds, and stronger pose priors that risk fitting one dataset’s statistics. Their own position was that simple trees, trained jointly and discriminatively with their part models, were fairly effective.

Deep models with structured inference

As convolutional networks began to outperform part-based models, which Tompson et al. (2014) state of the recent work they cite, several pose networks kept a graph over parts and its inference. Tompson, Jain, LeCun and Bregler (NIPS 2014) put an “MRF-like” spatial model on top of a convolutional part detector that outputs a heatmap per joint. Their graph is fully connected rather than a tree, and each pairwise term is a learned convolution kernel, the conditional distribution of one joint’s position given another’s. The marginal for joint $A$ is

$$\bar{p}_A = \frac{1}{Z} \prod_{v \in V} \left( p_{A \mid v} * p_v + b_{v \to A} \right),$$

which the authors liken to a single round of sum-product belief propagation, with no guarantee of reaching a global optimum because the graph has loops. For training they replaced it with a log-space version using SoftPlus and ReLU stages, which they say is no longer exactly an MRF, used kernels as large as $128 \times 128$ on $90 \times 60$ heatmaps (computed with FFTs), and trained the detector and spatial model jointly by back-propagation. In their Figure 8b (wrist, FLIC) the spatial model added 8 to 12% to the detection rate at large distance thresholds, and joint training another 4 to 5%.

Chen and Yuille (NIPS 2014) kept a tree and exact dynamic programming but made the springs depend on the image. A deep network looks at the patch around each candidate part and predicts both whether the part is there and which of several discrete spatial relations it has with each neighbour; the logarithms of those probabilities become the unary terms and image-dependent pairwise terms, added to standard quadratic deformation costs and weighted by a structural SVM. With one relation type per edge the model reduces to a pictorial structure with learned unaries, and the authors show how Yang and Ramanan’s mixtures are a further special case. Inference with distance transforms costs $O(T^2 L K)$ for $K$ parts once the network has run.

Yang, Ouyang, Li and Wang (CVPR 2016) went one step further and trained the network and a deformable mixture of parts end to end. Their message passing layers perform one iteration of max-sum message passing each, with messages flowing along every edge at once, accelerated with distance transforms to $O(L T^2)$ per message; three such layers were enough in practice, and the same layers serve tree-structured or loopy graphs. For a tree, a single pass from the leaves to the root is the exact dynamic programme. In the same period, Girshick, Iandola, Darrell and Malik (CVPR 2015) showed that a DPM unrolls into a convolutional network with a distance transform pooling layer, as described in the deformable part models article. In the stacked hourglass network and HRNet, by contrast, spatial context is learned by the convolutional stages themselves, with no explicit part graph.

Properties

• Exact inference on trees. For a tree with the Mahalanobis springs, dynamic programming finds the global optimum over the discretized locations, with no initialization, in $O(h' n)$ time; the result does not depend on the choice of root. Felzenszwalb and Huttenlocher note that this is essentially the cost of matching each part separately.
• Cost of the alternatives. General graphs make exact minimization NP-hard. Constellation models with a joint Gaussian over all parts face $O(N^P)$ hypotheses and rely on search heuristics.
• Appearance is modular. Learning separates the appearance model of each part from the spatial model, so any part detector with a likelihood or a score can be used: Gaussian derivative responses, foreground counts, boosted shape context, HOG templates or network outputs.
• Double counting. The product likelihood and the tree prior do not stop two parts with similar appearance, such as the two legs, from claiming the same image region. Samples, non-tree terms, type co-occurrence or stronger priors are the usual remedies.
• Multiple instances. The root scores give the best configuration for every root location, so several objects can be read off by taking several maxima, as with faces in the 2005 paper and with detections followed by non-maximum suppression in DPMs and Yang and Ramanan’s model.
• Discretization. Locations live on a grid, and the parent’s transformed location is rounded to the nearest grid point of the child’s transform. The accuracy of a match is limited by the grid, as with heatmap decoding.

Examples

Tree inference with distance transforms on a stick figure

The example matches a 14-part body model to a synthetic $64 \times 64$ image of a stick figure. Joints are drawn as small Gaussian blobs (standard deviation 1.2 pixels) and the head as a disk of radius 3, with faint limb lines, 12 clutter blobs that look exactly like joints, and Gaussian noise of standard deviation 0.15. The figure’s joints are a fixed template plus random offsets of standard deviation 1.5 pixels. Parts are located by pixel position only, as in Yang and Ramanan’s model. Every joint part uses the same appearance score, a Gaussian-smoothed copy of the image scaled so that an ideal joint blob scores 1, and the head uses a coarser smoothing scaled for the disk; by appearance alone, a left wrist cannot be told from a knee or from clutter. The tree is rooted at the neck, and each edge has the quadratic spring $0.02\,\big((dx - \mu_x)^2 + (dy - \mu_y)^2\big)$, with $(\mu_x, \mu_y)$ the template’s offset from parent to child.

The code compares three ways of finding the best configuration. On a three-part chain on an $8 \times 8$ grid, it enumerates all $64^3$ configurations, to check the dynamic programme itself. On the full model, it computes every message either by brute force, maximizing over all child locations for every parent location, or with separable 1D generalized distance transforms, written in max form and evaluated at shifted query positions (the parent location plus the rest offset, which may fall outside the grid). Speed is measured by counting operations, not by timing: brute force counts parent–child pairs, and the distance transform counts parabola intersections plus envelope readouts. Image sizes from 24 to 64 pixels show how the counts grow. The last part compares the configuration found by the tree with the independent maximum of each part’s appearance score, on seed 0 and then over seeds 0 to 19. NumPy 2.3.5 and SciPy 1.16.3; the script takes about 10 seconds.

import numpy as np
from scipy.ndimage import gaussian_filter

# ---- the model: 14 joint parts, a kinematic tree rooted at the neck ----------
NAMES = ["neck", "head", "r_sho", "r_elb", "r_wri", "l_sho", "l_elb", "l_wri",
         "r_hip", "r_kne", "r_ank", "l_hip", "l_kne", "l_ank"]
PARENT = [-1, 0, 0, 2, 3, 0, 5, 6, 0, 8, 9, 0, 11, 12]
TEMPLATE = np.array([[32, 15], [32, 8], [25, 17], [21, 26], [19, 35], [39, 17],
                     [43, 26], [45, 35], [28, 36], [27, 47], [26, 58], [36, 36],
                     [37, 47], [38, 58]], float)          # (x, y) on a 64 x 64 canvas
A = 0.02                                              # spring stiffness per px^2
N_CLUTTER, NOISE, SIG_BLOB, SIG_HEAD = 12, 0.15, 1.2, 2.5


def blob(S, x, y, s):
    yy, xx = np.mgrid[0:S, 0:S]
    return np.exp(-((xx - x) ** 2 + (yy - y) ** 2) / (2 * s * s))


def disk(S, x, y, r):
    yy, xx = np.mgrid[0:S, 0:S]
    return gaussian_filter(((xx - x) ** 2 + (yy - y) ** 2 <= r * r).astype(float), 0.7)


def render(S, seed):
    """Synthetic image of a stick figure: joints are blobs, the head a disk,
    limbs faint lines, plus clutter blobs and noise. Returns image and truth."""
    rng = np.random.default_rng(seed)
    truth = TEMPLATE * S / 64 + rng.normal(0, 1.5, TEMPLATE.shape)
    img = np.zeros((S, S))
    for i, (x, y) in enumerate(truth):
        img += disk(S, x, y, 3) if NAMES[i] == "head" else blob(S, x, y, SIG_BLOB)
        if PARENT[i] >= 0:                            # faint limb line to the parent
            for t in np.linspace(0, 1, 12):
                px, py = truth[PARENT[i]] + t * (truth[i] - truth[PARENT[i]])
                img += 0.3 * blob(S, px, py, 0.6) / 12 * 3
    for x, y in rng.uniform(4, S - 4, (N_CLUTTER, 2)):
        img += blob(S, x, y, SIG_BLOB)
    img += rng.normal(0, NOISE, img.shape)
    return img, truth


def unaries(img):
    S = img.shape[0]
    c = S // 2
    jn = gaussian_filter(blob(S, c, c, SIG_BLOB), SIG_BLOB).max()   # ideal responses
    hn = gaussian_filter(disk(S, c, c, 3), SIG_HEAD).max()
    U_joint = gaussian_filter(img, SIG_BLOB) / jn
    U_head = gaussian_filter(img, SIG_HEAD) / hn
    return [U_head if n == "head" else U_joint for n in NAMES]


# ---- 1-D generalized distance transform (max form), with operation counts ------
def dt1d(f, a, queries, ops):
    """D[k] = max_p f[p] - a*(p - queries[k])**2 for sorted queries (may lie
    off the grid). Upper envelope of n parabolas = lower envelope for -f
    (Felzenszwalb & Huttenlocher's algorithm). Returns D and the maximizing p."""
    n = len(f)
    g = -f
    h = g + a * np.arange(n) ** 2
    v = np.zeros(n, int)
    z = np.empty(n + 1)
    z[0], z[1] = -np.inf, np.inf
    k = 0
    for q in range(1, n):
        s = (h[q] - h[v[k]]) / (2 * a * (q - v[k]))
        ops[0] += 1
        while s <= z[k]:
            k -= 1
            s = (h[q] - h[v[k]]) / (2 * a * (q - v[k]))
            ops[0] += 1
        k += 1
        v[k], z[k], z[k + 1] = q, s, np.inf
    D, arg = np.empty(len(queries)), np.empty(len(queries), int)
    k = 0
    for j, qy in enumerate(queries):
        while z[k + 1] < qy:
            k += 1
        p = v[k]
        D[j], arg[j] = f[p] - a * (p - qy) ** 2, p
        ops[0] += 1
    return D, arg


def message_dt(f, a, mu, ops):
    """m[yp, xp] = max over child (yc, xc) of f - a*((yc-yp-muy)^2 + (xc-xp-mux)^2),
    and the maximizing child location, by two separable 1-D passes."""
    H, W = f.shape
    T, ay = np.empty((H, W)), np.empty((H, W), int)
    for x in range(W):
        T[:, x], ay[:, x] = dt1d(f[:, x], a, np.arange(H) + mu[1], ops)
    M, ax = np.empty((H, W)), np.empty((H, W), int)
    for y in range(H):
        M[y], ax[y] = dt1d(T[y], a, np.arange(W) + mu[0], ops)
    by = np.take_along_axis(ay, ax, axis=1)
    return M, by * W + ax                       # flat index of the best child


def message_brute(f, a, mu, ops):
    H, W = f.shape
    yy, xx = np.divmod(np.arange(H * W), W)
    fl = f.ravel()
    M, arg = np.empty(H * W), np.empty(H * W, int)
    for s in range(0, H * W, 512):
        py, px = yy[s:s + 512, None] + mu[1], xx[s:s + 512, None] + mu[0]
        sc = fl[None, :] - a * ((yy[None, :] - py) ** 2 + (xx[None, :] - px) ** 2)
        arg[s:s + 512] = sc.argmax(1)
        M[s:s + 512] = sc[np.arange(len(sc)), arg[s:s + 512]]
    ops[0] += (H * W) ** 2
    return M.reshape(H, W), arg.reshape(H, W)


def infer(U, parent, mu, a, message):
    """Max-sum dynamic programming on a tree (leaves to root), then backtracking."""
    n, (H, W) = len(U), U[0].shape
    children = [[c for c in range(n) if parent[c] == i] for i in range(n)]
    order, score, best, ops = [0], [None] * n, [None] * n, [0]
    for i in order:
        order += children[i]
    for i in reversed(order):                    # children before parents
        score[i] = U[i] + sum(m for m, _ in (best[c] for c in children[i]))
        if parent[i] >= 0:
            best[i] = message(score[i], a, mu[i], ops)
    loc = [None] * n
    loc[0] = int(np.argmax(score[0]))
    for i in order[1:]:
        loc[i] = int(best[i][1].ravel()[loc[parent[i]]])
    return score[0].max(), np.array([[l % W, l // W] for l in loc]), ops[0]


def mu_for(S):
    t = TEMPLATE * S / 64
    return [None] + [t[c] - t[PARENT[c]] for c in range(1, len(NAMES))]


# ---- 1. dynamic programming against exhaustive enumeration (3-part chain) -----
rng = np.random.default_rng(1)
Ut = [gaussian_filter(rng.normal(0, 1, (8, 8)), 0.8) for _ in range(3)]
mut = [None, np.array([2.0, 1.0]), np.array([-1.0, 2.0])]
yy, xx = np.divmod(np.arange(64), 8)
def cost(m):                                      # cost[parent loc, child loc]
    return 0.3 * ((xx[None] - xx[:, None] - m[0]) ** 2 + (yy[None] - yy[:, None] - m[1]) ** 2)
E = (Ut[0].ravel()[:, None, None] + Ut[1].ravel()[None, :, None] + Ut[2].ravel()[None, None, :]
     - cost(mut[1])[:, :, None] - cost(mut[2])[None, :, :])
p0, p1, p2 = np.unravel_index(E.argmax(), E.shape)
s_dp, loc_dp, _ = infer(Ut, [-1, 0, 1], mut, 0.3, message_dt)
print(f"chain, 8x8 grid: exhaustive {E.max():.6f} over {E.size} configurations, "
      f"DP+DT {s_dp:.6f}, same configuration: "
      f"{[p0, p1, p2] == [l[1] * 8 + l[0] for l in loc_dp]}")

# ---- 2. distance transforms against brute-force messages, at several sizes ----
print(f"\n{'size':>4} {'h':>5} {'brute-force ops':>16} {'DT ops':>9} {'ratio':>6} "
      f"{'max |score diff|':>17} same")
for S in (24, 32, 48, 64):
    img, truth = render(S, 0)
    U = unaries(img)
    sb, lb, ob = infer(U, PARENT, mu_for(S), A, message_brute)
    sd, ld, od = infer(U, PARENT, mu_for(S), A, message_dt)
    print(f"{S:>4} {S * S:>5} {ob:>16,} {od:>9,} {ob / od:>6.0f} {abs(sb - sd):>17.2e} "
          f"{np.array_equal(lb, ld)}")

# ---- 3. localization on the 64 x 64 image: independent parts against the tree --
img, truth = render(64, 0)
U = unaries(img)
_, loc, _ = infer(U, PARENT, mu_for(64), A, message_dt)
indep = np.array([[np.argmax(u) % 64, np.argmax(u) // 64] for u in U])
e_tree = np.hypot(*(loc - truth).T)
e_ind = np.hypot(*(indep - truth).T)
print("\nseed 0, error in px  " + " ".join(f"{n:>5}" for n in NAMES))
print("independent argmax   " + " ".join(f"{e:5.1f}" for e in e_ind))
print("tree MAP             " + " ".join(f"{e:5.1f}" for e in e_tree))
print(f"independent argmax uses {len(np.unique(indep, axis=0))} distinct pixels for 14 parts")

def score_of(U, L, mu):                           # score of a given configuration
    return (sum(U[i][L[i][1], L[i][0]] for i in range(len(U)))
            - sum(A * np.sum((L[c] - L[PARENT[c]] - mu[c]) ** 2) for c in range(1, len(U))))

counts = []
print("\nseed  within 3 px  MAP - truth score  parts within 2 px of another part")
for seed in range(20):
    img, truth = render(64, seed)
    U = unaries(img)
    s_map, loc, _ = infer(U, PARENT, mu_for(64), A, message_dt)
    indep = np.array([[np.argmax(u) % 64, np.argmax(u) // 64] for u in U])
    counts.append(((np.hypot(*(indep - truth).T) <= 3).sum(),
                   (np.hypot(*(loc - truth).T) <= 3).sum()))
    if counts[-1][1] < 14:                       # explain the failures
        gap = s_map - score_of(U, np.clip(np.round(truth), 0, 63).astype(int), mu_for(64))
        d = np.hypot(*(loc[:, None] - loc[None]).transpose(2, 0, 1))
        shared = [NAMES[i] for i in range(len(NAMES)) if (d[i] <= 2).sum() > 1]
        print(f"{seed:>4} {counts[-1][1]:>12} {gap:>18.2f}  {len(shared)} {' '.join(shared)}")
c = np.array(counts)
print(f"\nseeds 0-19, joints within 3 px (of 14): independent {c[:, 0].min()}-{c[:, 0].max()} "
      f"(mean {c[:, 0].mean():.1f}), tree {c[:, 1].min()}-{c[:, 1].max()} "
      f"(mean {c[:, 1].mean():.1f}); tree better on {(c[:, 1] > c[:, 0]).sum()} of 20")

chain, 8x8 grid: exhaustive 1.025536 over 262144 configurations, DP+DT 1.025536, same configuration: True

size     h  brute-force ops    DT ops  ratio  max |score diff| same
  24   576        4,313,088    37,455    115          0.00e+00 True
  32  1024       13,631,488    66,395    205          0.00e+00 True
  48  2304       69,009,408   146,028    473          0.00e+00 True
  64  4096      218,103,808   255,754    853          1.78e-15 True

seed 0, error in px   neck  head r_sho r_elb r_wri l_sho l_elb l_wri r_hip r_kne r_ank l_hip l_kne l_ank
independent argmax     8.1  14.4  12.6  14.1  21.2   5.4   5.1  13.9  16.1  27.9  39.4  14.6  25.3  34.6
tree MAP               0.3   0.2   0.5   1.0   0.1   0.1   1.5   0.2   0.5   0.8   0.2   0.5   0.4   0.4
independent argmax uses 2 distinct pixels for 14 parts

seed  within 3 px  MAP - truth score  parts within 2 px of another part
   1           11               1.51  2 r_ank l_ank
   3           13               1.77  0 
   8           13               0.78  0 
   9           10               2.34  7 r_wri r_hip r_kne r_ank l_hip l_kne l_ank
  13            5               5.71  8 neck head r_sho r_elb r_wri l_sho l_elb l_wri
  14           13               0.95  0 
  15            1               0.49  4 neck head r_sho l_sho
  17           13               1.21  0 
  18           13               1.55  0 
  19            7               5.24  8 neck head r_sho r_elb r_wri l_sho l_elb l_wri

seeds 0-19, joints within 3 px (of 14): independent 0-2 (mean 1.2), tree 1-14 (mean 11.9); tree better on 19 of 20

The first two blocks check correctness. The dynamic programme reaches the exhaustive optimum on the small chain, and at every size the distance transform reproduces the brute-force messages: the same configuration, with scores equal up to floating-point rounding. The counts show the difference in growth. Brute force costs exactly $h^2$ per edge, $h^2 (n - 1)$ in all, while the distance transform needs about five operations per grid point per edge (by our count, two 1D passes, each with somewhat more than one intersection and exactly one readout per point), so the ratio grows in proportion to $h$, from about a hundred at $24 \times 24$ to over eight hundred at $64 \times 64$.

The localization block shows why the springs matter. Taken independently, the fourteen parts’ best locations fall on just two pixels, and on seed 0 no part is within 3 pixels of its true position. The tree places all fourteen parts within 1.5 pixels on that seed, and over the twenty seeds it got between 1 and 14 joints within 3 pixels, against 0 to 2 for independent maxima, and was better on 19 of the 20.

The failure table is the most instructive part. On every failing seed the MAP configuration scores higher than the ground-truth configuration (rounded to the grid), so the search did not fail: it found what the model prefers. Half the failures, including the three worst, involve parts piled close together: on seeds 13 and 19 the neck, the head and all six arm joints each lie within 2 pixels of another part, on seed 9 the six leg joints and the right wrist do, on seed 15 the neck, the head and both shoulders, and on seed 1 the two ankles lie within 2 pixels of each other. Seeds 3, 8, 14, 17 and 18 miss a single joint without any overlap. These are the double-counting failures that Yang and Ramanan describe for tree models, produced here by a model with identical appearance for all joints, which is a deliberately hard setting. Real models avoid some of them with part-specific appearance, as in Yang and Ramanan’s typed templates, and Felzenszwalb and Huttenlocher’s person model handled over-counting by sampling and re-scoring.

The script was changed twice after its first run, both times only in what it prints. The first run printed the correctness checks, the counts, the seed 0 errors and the summary over seeds. To find out why some seeds failed, two measurements were then added: the number of distinct pixels used by the independent maxima on seed 0 and, for each seed where the tree missed a joint, the score difference from the ground-truth configuration and the number of parts within 2 pixels of another part; a second edit printed those parts by name instead of only counting them. The image, model, parameters and inference are unchanged, and the lines printed by the first run are identical. The plan written before the first run expected that “only the head is likely right” under independent maxima and that “the tree is expected to be better on every seed”; neither held, as the output shows.

Common Misconceptions

• “Pictorial structures means an articulated body model.” The framework is generic. Felzenszwalb and Huttenlocher’s demonstrations included a face model with parts located by position only, and DPMs use the same machinery for cars and bicycles.
• “Dynamic programming over a tree is approximate.” On a tree it is exact over the discretized location space, for any root. Approximation enters with loopy graphs, with the rounding of transformed parent locations to the grid, and in later networks that run a fixed number of message passing rounds.
• “The distance transform approximates the maximization.” For quadratic springs it computes the same maximum as brute force, as the example’s identical configurations show; it is faster because of the structure of the cost, not because it skips candidates.
• “A higher-scoring configuration is a better pose.” Exact inference returns the model’s preferred configuration. When parts share appearance and the tree has no term preventing overlap, that preference can be wrong, as the failing seeds of the example show; Felzenszwalb and Huttenlocher’s sampling-and-rescoring step exists for this reason.
• “Constellation models are pictorial structures with fewer parts.” They differ in structure and search: a joint Gaussian over all parts evaluated on a sparse set of detected features, against a tree evaluated densely at every location.

Where It Is Used

Yang and Ramanan (2011) called pictorial structures with articulated parts the dominant approach to human pose estimation; the line runs from Felzenszwalb and Huttenlocher’s silhouettes through Andriluka, Roth and Schiele’s boosted part detectors to Yang and Ramanan’s mixtures of HOG parts, evaluated on small datasets such as Parse and Buffy and later on MPII Human Pose. In detection, the star-structured deformable part models carried the idea to PASCAL VOC, and Felzenszwalb and Huttenlocher note that models of local patches and their spatial relations had been popular for faces, as in their own model of eyes, nose and mouth corners (a subject of face detection). The generalized distance transform itself has wider uses, from Euclidean distance transforms of binary images to Viterbi decoding with large state spaces and belief propagation for stereo and optical flow, as Felzenszwalb and Huttenlocher’s 2012 paper lists.

In deep pose estimation, explicit part graphs appeared as spatial models on top of convolutional heatmaps, as image-dependent pairwise terms and as message passing layers trained end to end, before networks such as the stacked hourglass network and HRNet learned spatial context with convolutional stages instead.