Heatmap Regression

Heatmap regression locates keypoints, such as body joints or facial landmarks, by training a network to output one image-like map per keypoint instead of the coordinates themselves. The training target for each map is a small bump, usually a 2D Gaussian, centred on the annotated position; at test time the keypoint is read off the map, classically at its maximum. Tompson et al. (2014) trained a human pose estimation network this way, and later papers date the heatmap representation to their work (Gu, Yang and Yao, 2021); it became the most common output representation of pose networks (Zhang et al., 2020), including the stacked hourglass network, HRNet and OpenPose. Its weak point is resolution: the map is smaller than the input, so the maximum lands on a coarse grid. Much of the method’s later history is about recovering sub-pixel coordinates from that grid (the quarter-cell shift, DARK, offset fields), making the readout differentiable (integral regression, also called soft-argmax), or replacing the 2D map with coordinate representations that avoid it (direct regression with a learned likelihood, coordinate classification).

The name is used loosely. Gu, Yang and Yao (2021) point out that heatmap-based methods are called “detection-based”, “indirect regression” or “heatmap regression” in different papers, and Sun et al. (2018) call the same methods detection based. This article uses heatmap regression for the representation: a dense map regressed against a target map and decoded into a coordinate.

Definition

For a keypoint at position $\mathbf{x}$ in an input image, heatmap regression trains a convolutional network to output a map $H$ on a grid of stride $\lambda$ (input pixels per map cell), whose target $G$ peaks at $\mathbf{x}/\lambda$. The network is trained by a per-cell loss between $H$ and $G$, and a decoder $D$ maps the predicted $H$ back to a coordinate, $\hat{\mathbf{x}} = \lambda\, D(H)$ plus any change of frame from the crop to the image. The three design choices are the target (encoding), the loss, and the decoder.

Intuition

A network built from convolutions produces outputs that move with the image: shift a wrist by eight pixels and the wrist’s activation moves by eight pixels, give or take the stride. A heatmap uses that property directly. Each cell only has to answer a local question, “is the keypoint here?”, and the answer for every cell is supervised, which gives far more training signal than two numbers per keypoint. Sun et al. (2018) give this dense, per-pixel supervision as a likely reason that heatmap learning is easier than regressing coordinates. Tompson, Goroshin et al. (2015) add that direct regression copes badly with multi-modal outputs, for example two plausible positions for a joint; a map can hold two bumps, while a regressed coordinate must commit to one.

The price is that the map is a grid. Pooling and strided convolutions shrink it to a fraction of the input, often a quarter, and the maximum of a grid can only be at a cell. A keypoint that falls between cells is reported at the nearest one, an error of up to half a cell in each direction. Since a cell is several input pixels wide, the error is several pixels, which strict metrics such as COCO’s average precision at high keypoint similarity, or PCKh at a tenth of the head size, can see.

Formal Definition

Gaussian targets

Tompson et al. (2014) trained their part detector with a mean squared error between the predicted map and a target that is a 2D Gaussian with a small variance centred at the joint. Tompson, Goroshin et al. (2015) wrote the objective out, summing over joints $j$ and map positions $(x, y)$,

$$E_1 = \frac{1}{N} \sum_{j=1}^{N} \sum_{x,y} \left( H'_j(x, y) - H_j(x, y) \right)^2,$$

with predicted maps $H'_j$, target maps $H_j$, and a constant variance of about 1.5 pixels. For a keypoint at $\boldsymbol{\mu} = \mathbf{x}/\lambda$ in map coordinates, the target used by most later work is the unnormalized Gaussian

$$G(\mathbf{p}) = \exp\left( -\frac{\lVert \mathbf{p} - \boldsymbol{\mu} \rVert^2}{2\sigma^2} \right),$$

with peak value 1 and spread $\sigma$ in map cells. The released code of the stacked hourglass network uses $\sigma = 1$ on $64 \times 64$ maps for $256 \times 256$ inputs (pose-hg-train, opts.lua); HRNet’s released COCO configuration uses $\sigma = 2$ on $64 \times 48$ maps for $256 \times 192$ inputs, draws the Gaussian only within $3\sigma$ of its centre, and zeroes the training weight of a joint whose Gaussian falls entirely outside the map (HRNet code, JointsDataset.py and the w32_256x192 configuration). Its $384 \times 288$ configuration, also at stride 4, uses $\sigma = 3$ on $96 \times 72$ maps, so the spread grows with the crop rather than staying fixed in cells (w48_384x288 configuration). DARK’s paper writes the target with the normalizing factor $1/(2\pi\sigma^2)$ (Zhang et al., 2020, Eq. 14), while its released code, like HRNet’s, draws a unit-height Gaussian (DarkPose code, JointsDataset.py).

Biased and unbiased encoding

Where exactly the Gaussian is centred matters. HRNet’s code computes the centre as int(x / stride + 0.5), rounding the keypoint to the nearest cell before drawing the target (HRNet code, JointsDataset.py). Zhang et al. (2020) call such targets biased: the training label itself carries a quantization error of up to half a cell, so even a perfect network would learn to put its peak in the wrong place. Their remedy is to centre the Gaussian at the unrounded position $\mathbf{x}/\lambda$, so the target’s true peak can lie between cells. Huang et al. (2020), in their study of unbiased data processing (UDP), analysed the same rounding as part of the “standard encoding-decoding”, alongside a second source of bias in the coordinate transforms used for flip testing.

Losses

The default loss is the per-cell squared error above, averaged over cells and joints. In HRNet’s released code it is half the mean squared error per joint, with each joint’s predicted and target maps multiplied by its training weight, so that unlabelled joints contribute nothing (HRNet code, loss.py); SimpleBaseline used the same mean squared error against Gaussian targets (Xiao, Wu and Wei, 2018). Variants change either the per-cell penalty or which errors count:

• Online hard keypoint mining. In the cascaded pyramid network, the second stage keeps, for each person, only the largest $M$ of the $N$ per-keypoint losses. On COCO minival, $M = 8$ of 17 gave the best result, 69.4 AP against 68.6 with all keypoints (Chen et al., 2018, Table 6). HRNet’s code includes the same loss with a default of eight keypoints (loss.py).
• Classification losses. Sun et al. (2018) list two alternatives to the Gaussian target with mean squared error: a one-hot target with a cross-entropy over a softmax of all cells, as in Mask R-CNN’s keypoint head, and per-cell binary cross-entropy against a disk of positives within 15 pixels of the joint. Papandreou et al. (2017) used binary classification against disks of radius $R = 25$ pixels, with offsets for precision (see Offset Fields below).
• Adaptive Wing loss. For facial landmarks, Wang, Bo and Li (2019) argued that mean squared error is insensitive to small errors and lets the many background cells dominate, producing blurry, low peaks. Their loss behaves like a logarithmic, Wing-like penalty with a large influence of small errors on cells where the target is near 1, and like mean squared error where the target is near 0, by letting the exponent depend on the target value $y$:

$$\mathrm{AWing}(y, \hat{y}) = \begin{cases} \omega \ln\left( 1 + \left| \dfrac{y - \hat{y}}{\epsilon} \right|^{\alpha - y} \right) & \text{if } |y - \hat{y}| < \theta, \\ A |y - \hat{y}| - C & \text{otherwise,} \end{cases}$$

with $A$ and $C$ chosen to make the two pieces join smoothly, and a weighted loss map that up-weights foreground and nearby difficult background cells.

Decoders

The classical decoder is the argmax, $D(H) = \arg\max_{\mathbf{p}} H(\mathbf{p})$. The refinements below keep the argmax as a starting point and estimate a sub-cell correction from the neighbouring cells or from extra network outputs. Integral regression replaces the argmax by an expectation; it has its own section below.

Resolution and Quantization Error

If keypoints fall uniformly within cells and the argmax always lands on the nearest cell, the error along each axis is uniform on $[-\tfrac12, \tfrac12]$ cells. Huang et al. (2020) derive the consequences per axis: an expected absolute error of $1/4$ cell with variance $1/48$ for the plain argmax, falling to $1/8$ cell with variance $1/192$ for the quarter-cell shift described below. In input pixels the errors are multiplied by the stride (the variances by its square), so halving the heatmap resolution doubles the error.

Tompson, Goroshin et al. (2015) measured the effect on FLIC with networks that used 4×, 8× or 16× pooling. More pooling lowered the detection rate at small distance thresholds, and with 8× pooling the coarse model’s horizontal error, in the pixel coordinates of its $360 \times 240$ input, was approximately uniform between −4 and +4 pixels, as quantization predicts; the face’s error had a standard deviation of 3.13 pixels (Fig. 11a). A second network that refined each joint from features cropped around the coarse maximum reduced it: for their full model, Table 1 lists 1.46 pixels for the face with 8× pooling and 1.09 with 4× pooling, against 0.65 pixels for the spread of their own repeated human annotations (13 annotators on 10 images). The histogram of the cascaded model in Fig. 11b, which the text presents after the 8× coarse histogram without naming its pooling, is labelled 1.09 for the face and 2.43 for the shoulder, the table’s 4× values.

The standard ways to fight quantization trade one cost for another. Larger maps cost computation and memory: SimpleBaseline appends three transposed convolutions to its ResNet (Xiao, Wu and Wei, 2018), which recover the resolution from 1/32 to 1/4 of the input (Li et al., 2022), and HRNet keeps a high-resolution branch throughout the network. Larger inputs scale every map. In the integral regression paper’s MPII ablation, halving the heatmap from $64 \times 64$ to $32 \times 32$ on $256 \times 256$ inputs cut the heatmap baseline’s PCKh@0.1 from 28.0 to 17.2 but saved only 15% of the computation (Sun et al., 2018, Table 2). DARK’s results show the same dependence on input size from the other side: with HRNet-W32 on COCO val2017, DARK added 3.8 AP at $128 \times 96$ input, 1.2 AP at $256 \times 192$ and 0.8 AP at $384 \times 288$ (Zhang et al., 2020, Table 4). Sub-pixel decoding matters most when maps are coarse.

Sub-Pixel Decoding

The quarter-cell shift

The oldest refinement moves the argmax a quarter of a cell toward the higher of its two neighbours, separately along each axis. The released training code of the stacked hourglass network applies it as a “very simple post-processing step to improve performance at tight PCK thresholds” (pose-hg-train, pose.lua), and HRNet’s code does the same, adding $0.25 \cdot \operatorname{sign}(H(p_x + 1, p_y) - H(p_x - 1, p_y))$ to the $x$ coordinate and likewise for $y$ (HRNet code, inference.py). The rule is a guess at which half of the cell holds the true peak; when the guess is right, the error in that half-cell interval is at most a quarter of a cell. DARK’s paper describes the standard shift differently, as a step of 0.25 pixels from the maximum toward the second-highest activation (Zhang et al., 2020, Eq. 1); the released hourglass and HRNet codes, both among the models DARK was compared with, use the per-axis sign rule.

The shift has a large effect on reported numbers. Zhang et al. (2020) measured 61.2 AP without it and 66.9 with it for HRNet-W32 at $128 \times 96$ input on COCO val2017; the human pose estimation article gives the full comparison with DARK.

DARK: Taylor expansion of the log-heatmap

Distribution-aware coordinate representation of keypoints (DARK) treats the predicted map as a sampled Gaussian (Zhang et al., 2020). The logarithm of a Gaussian, $P(\mathbf{x}) = \ln G(\mathbf{x})$, is a quadratic with its maximum at $\boldsymbol{\mu}$, so a second-order Taylor expansion of $P$ at the argmax $\mathbf{m}$ locates the true peak exactly:

$$\boldsymbol{\mu} = \mathbf{m} - \left( D''(\mathbf{m}) \right)^{-1} D'(\mathbf{m}),$$

where $D'$ and $D''$ are the gradient and Hessian of $P$, estimated by finite differences of the log-heatmap. This is one Newton step on the log-map. Because predicted maps are rarely clean Gaussians and often show several small peaks near the maximum, DARK first smooths the map with a Gaussian kernel of the training spread (“distribution modulation”), then rescales it. The paper rescales by min–max normalization to the original maximum (Eq. 11); the released code smooths with a $11 \times 11$ OpenCV Gaussian kernel by default and rescales only by the ratio of the maxima, then takes the logarithm after clipping at $10^{-10}$ (DarkPose code, inference.py and default.py). Modulation alone added 0.3 AP in the paper’s ablation (Table 2). The paper’s other component is the unbiased encoding above; with both, HRNet-W32 at $128 \times 96$ gained 3.8 AP over the standard pipeline (Table 3), and the gains held across stacked hourglass, SimpleBaseline and HRNet models at that input size, from 2.4 to 4.4 AP (Table 5).

Offset fields

Papandreou et al. (2017) made the sub-cell correction a learned output. Their network (G-RMI) predicts, for each keypoint $k$, a heatmap $h_k$ that classifies whether each position lies within a disk of radius $R$ of the keypoint, and an offset field $F_k(\mathbf{x}_i) = \mathbf{l}_k - \mathbf{x}_i$ from each position to the keypoint. The two are fused by voting,

$$f_k(\mathbf{x}_i) = \sum_j \frac{1}{\pi R^2}\, G\!\left( \mathbf{x}_j + F_k(\mathbf{x}_j) - \mathbf{x}_i \right) h_k(\mathbf{x}_j),$$

with $G$ the bilinear interpolation kernel: every position casts a vote at the location its offset points to, weighted by its heatmap probability, a form of Hough voting. Training combines logistic losses on the heatmaps with a Huber loss on offsets inside the disks, and the predictions, made at stride 8, are bilinearly upsampled to the $353 \times 257$ crop. UDP adopted this combined classification and regression encoding to remove the expected decoding error of the standard method: disk-shaped heatmaps and offset maps, decoded by taking the maximum of a Gaussian-smoothed heatmap and adding the smoothed offsets there (Huang et al., 2020). In their ablation on COCO val2017 with HRNet-W32 at $256 \times 192$, replacing the standard encoding and decoding raised AP from 75.6 to 76.6 with flip testing and the usual one-pixel shift of the flipped maps (Table 3, methods F and J). The same ablation shows how much the flip-testing details matter: without that one-pixel shift, flip testing gave 73.3 AP, lower than the 74.5 AP of no flip testing at all.

Integral Regression

The argmax is not differentiable, so heatmap networks are trained on maps rather than on the coordinates that are evaluated; Tompson, Goroshin et al. (2015) noted that a more direct objective would measure the argmax, and regressed to target maps because they could not. Sun et al. (2018) replaced “taking-maximum” by “taking-expectation”. The predicted map is normalized with a softmax over all positions,

$$\tilde{H}_k(\mathbf{p}) = \frac{e^{H_k(\mathbf{p})}}{\sum_{\mathbf{q} \in \Omega} e^{H_k(\mathbf{q})}},$$

and the joint is its expected position,

$$\mathbf{J}_k = \sum_{\mathbf{p} \in \Omega} \mathbf{p}\, \tilde{H}_k(\mathbf{p}).$$

The operation is differentiable, adds almost no computation, produces continuous coordinates and can be trained with a loss on the coordinates; Sun et al. found an L1 loss consistently better than L2 for this. The operation itself was known as soft-argmax, and Sun et al. name two contemporary uses for human pose, by Luvizon et al. and by Nibali et al. (the latter’s differentiable spatial to numerical transform, DSNT); their own contribution was an extensive set of ablations, including 3D pose with volumetric heatmaps, where the expectation is taken over depth as well.

In their MPII ablation on a validation split of about 3,000 images with ResNet-50, the integral version of the Gaussian-target model raised PCKh@0.5 only from 86.8 to 87.3 but PCKh@0.1 from 17.2 to 29.3 (Table 1), and it degraded much less when the heatmap or input was made smaller (Table 2). The table’s default setting, which reproduces the Table 1 numbers, is the $32 \times 32$ heatmap, although the text describes the head as upsampling to $64 \times 64$ by default. On COCO test-dev with ResNet-101 at $256 \times 256$, the integral model scored 67.8 AP against 66.3 for the same network decoded by argmax (Table 6).

Expectation has its own failure mode. The softmax gives every cell some probability, including cells where the network outputs nothing, and those cells pull the expectation toward the centre of the map. Gu, Yang and Yao (2021) derived this bias: it vanishes only when the peak is at the map’s centre and grows with the distance from it and with a smaller softmax temperature parameter $\beta$ in $e^{\beta H}$. With $\beta = 10$ on a $64 \times 48$ map, they computed biases of 0.2, 0.7 and 1.2 pixels for joints at 0.2, 0.5 and 0.8 of the maximum distance from the centre (Fig. 2 caption). A network can learn to push its output further out to compensate, which they suggest makes training harder and slower: in their COCO experiments the heatmap-trained model was close to its final accuracy after about 40 epochs, while the integral model needed about 90 to match it (Fig. 4). They proposed a closed-form correction of the expectation, assuming that the mass belonging to the peak lies within a region symmetric about it, and a training schedule that adds the heatmap loss to the corrected integral loss for the first epochs and drops it afterwards.

Alternatives to 2D Heatmaps

Direct regression

DeepPose regressed all joint coordinates, normalized by the person’s box, with an L2 loss and a cascade of refinement stages (Toshev and Szegedy, 2014). Tompson et al. (2014) attributed its inaccuracy in the high-precision region to the difficulty of learning that mapping directly, and Sun et al. (2018) found their regression baseline, a ResNet with average pooling and a fully connected output, below both the heatmap and integral models at PCKh@0.5 on MPII, although above the heatmap model at PCKh@0.1 (Table 1).

Residual log-likelihood estimation (RLE) revisited regression as maximum likelihood estimation (Li et al., 2021). A regression loss assumes an output density: an L2 loss is the negative log-likelihood of a Gaussian with fixed variance, an L1 loss of a Laplace distribution with fixed scale. RLE has the network predict a location $\hat{\mu}$ and a scale $\hat{\sigma}$ per joint and lets a normalizing flow (RealNVP) learn the shape of the error density, as a residual on top of a simple base density, during training only. The flow is not run at test time, so inference costs the same as plain regression. On COCO val2017, a ResNet-50 regression model rose from 58.1 AP with an L1 loss to 70.5 AP with RLE, and to 71.3 AP when the backbone was first trained with a heatmap loss (Table 1), which the authors compare with 71.0 AP for SimpleBaseline, listed as SimplePose (Table 2). The text’s claimed gain over integral regression, 7.5 AP, matches the 70.5 AP row; against the starred 71.3 AP row printed beside the integral model’s 63.0 AP in Table 2, the gap is 8.3. On COCO test-dev, HRNet-W48 with RLE reached 75.7 AP, against 75.5 for HRNet-W48 decoded from heatmaps as listed in the same table (Table 3).

Coordinate classification: SimCC

SimCC turns each coordinate into a classification over bins (Li et al., 2022). For an input of width $W$ and height $H$ and a splitting factor $k \geq 1$, the horizontal and vertical coordinates are each classified into $W \cdot k$ and $H \cdot k$ bins by one linear layer on the flattened keypoint features, trained with the Kullback–Leibler divergence against smoothed labels. Since a bin can be smaller than an input pixel, the quantization error is below a pixel without a large 2D map, and the decoder needs no refinement. On COCO val2017, SimCC helped most at small inputs: for HRNet-W48 at $64 \times 64$ it scored 59.7 AP against 48.5 for the heatmap model with the quarter shift and 36.9 without it, while at $256 \times 192$ the margin was 75.9 against 75.1 (Table 1). With SimpleBaseline-ResNet-50, SimCC could drop the three transposed convolutions altogether and still match or exceed the heatmap model, 70.8 against 70.4 AP at $256 \times 192$ (Table 4). Larger $k$ was not always better: performance rose and then fell as $k$ went from 1 to 4, with $k = 2$ or 3 recommended (Fig. 4). Gaussian label smoothing gave the best of the variants tried, 73.4 AP against 71.3 without smoothing, with SimpleBaseline-ResNet-50 at $384 \times 288$ (Table 5). On MPII validation with HRNet-W32 at $256 \times 256$, SimCC’s PCKh@0.5 was slightly below the heatmap model’s (89.6 against 90.3), while its PCKh@0.1 was higher (37.8 against 33.1; Table 7).

Comparisons across representations

Comparisons between these families are hard to read because networks, input sizes, detectors and training schedules differ. Some papers hold the backbone fixed: DARK’s authors, with HRNet-W32 at $128 \times 96$ on COCO val2017, reported 57.6 AP for DSNT, 68.0 for integral regression and 70.7 for DARK (Zhang et al., 2020, Table 6); Gu, Yang and Yao (2021), with a common SimpleBaseline network, found that heatmap and integral models reached a similar final accuracy, with integral regression having lower endpoint error on people with few visible joints, heavy occlusion or small boxes (Fig. 4 and Table 1). In our reading, these results agree that the gap between heatmaps and coordinates is narrowest at high resolution and widest at low resolution, where quantization dominates.

Properties

• Translation equivariance. A fully convolutional heatmap head keeps the spatial structure of the features; Gu, Yang and Yao (2021) give this as the usual reason that heatmap methods outperform regression from pooled features.
• Multiple hypotheses. A map can hold several peaks, so the same representation serves bottom-up methods that detect every person’s joints at once, as in OpenPose, and lets a spatial model or a later stage choose between candidates.
• Resolution-bound accuracy. Argmax decoding has an error floor proportional to the stride; refinement, offsets, expectation or bins are needed to go below a cell.
• Spread $\sigma$ as a trade-off. A wider Gaussian gives more cells a non-zero target and more gradient, but a flatter top, where noise moves the maximum more easily (see the example). The released configurations above set $\sigma$ between 1 and 3 cells.
• Calibration. The peak height is commonly used as a keypoint confidence, for example by averaging over keypoints to rescore people (Papandreou et al., 2017), and RLE derives its confidence from the predicted scale $\hat{\sigma}$ instead.

Examples

Decoding error at several heatmap resolutions

The example measures how far each decoder lands from a known keypoint. It places keypoints uniformly at random, at sub-pixel positions, in the central part of a $256 \times 256$ input crop and draws, for each, the map a heatmap network would ideally output: a unit-height Gaussian centred at the unrounded position, on maps of stride 2, 4, 8 and 16 (from $128 \times 128$ down to $16 \times 16$ cells), with or without additive Gaussian noise of standard deviation 0.05 (5% of the peak). No network is trained, so the error comes only from decoding. Five decoders are compared:

• argmax;
• the quarter-cell shift, following HRNet’s code;
• DARK, smoothing with the training $\sigma$ as the paper specifies (the code’s default $11 \times 11$ OpenCV kernel corresponds to $\sigma = 2$) and otherwise following DarkPose’s code (rescaling by the maxima, logarithm, one Newton step), with SciPy’s Gaussian filter in place of OpenCV’s;
• integral regression applied to the heatmap itself (negative values clipped, then normalized to sum to 1);
• integral regression of an output trained for it, modelled as logits whose softmax is the Gaussian, plus the same noise.

Two choices that the sources leave open are printed both ways. The first is whether $\sigma$ is fixed in map cells (2 cells at every stride) or in input pixels (8 pixels, which is 2 cells at stride 4). The second is what a network’s output looks like before the expectation is taken, the last two decoders above.

import numpy as np
from scipy.ndimage import gaussian_filter

INPUT = 256                     # input crop size in pixels


def gaussian(shape, mu, sigma):
    """Unit-height Gaussian centred at mu = (x, y) in heatmap cells (unrounded)."""
    yy, xx = np.mgrid[0:shape[0], 0:shape[1]]
    return np.exp(-((xx - mu[0]) ** 2 + (yy - mu[1]) ** 2) / (2 * sigma ** 2))


def argmax(h):
    py, px = np.unravel_index(np.argmax(h), h.shape)
    return np.array([px, py], float)


def quarter(h):
    """HRNet code: shift the argmax 0.25 cell toward the higher neighbour on each axis."""
    p = argmax(h)
    px, py = int(p[0]), int(p[1])
    H, W = h.shape
    if 1 < px < W - 1 and 1 < py < H - 1:
        p += 0.25 * np.sign([h[py, px + 1] - h[py, px - 1], h[py + 1, px] - h[py - 1, px]])
    return p


def dark(h, sigma):
    """DarkPose code: smooth, rescale to the original max, log, one Taylor (Newton) step."""
    p = argmax(h)
    px, py = int(p[0]), int(p[1])
    H, W = h.shape
    g = gaussian_filter(h, sigma, mode="constant")
    g = np.log(np.maximum(g * h.max() / g.max(), 1e-10))
    if 1 < px < W - 2 and 1 < py < H - 2:
        d = 0.5 * np.array([g[py, px + 1] - g[py, px - 1], g[py + 1, px] - g[py - 1, px]])
        dxx = 0.25 * (g[py, px + 2] - 2 * g[py, px] + g[py, px - 2])
        dyy = 0.25 * (g[py + 2, px] - 2 * g[py, px] + g[py - 2, px])
        dxy = 0.25 * (g[py + 1, px + 1] - g[py - 1, px + 1] - g[py + 1, px - 1] + g[py - 1, px - 1])
        Hs = np.array([[dxx, dxy], [dxy, dyy]])
        if np.linalg.det(Hs) != 0:
            p = p - np.linalg.solve(Hs, d)
    return p


def expectation(p):
    """Integral regression: expected coordinate under a normalised map p."""
    yy, xx = np.mgrid[0:p.shape[0], 0:p.shape[1]]
    return np.array([(p * xx).sum(), (p * yy).sum()])


def soft_argmax_heatmap(h):
    """Integral decoding of a heatmap-trained output: clip negatives, normalise."""
    p = np.maximum(h, 0)
    return expectation(p / p.sum())


def soft_argmax_logits(z):
    """Integral decoding of logits: softmax over all cells, then expectation."""
    p = np.exp(z - z.max())
    return expectation(p / p.sum())


DECODERS = ["argmax", "quarter", "DARK", "int-heatmap", "int-logits"]
STRIDES = [2, 4, 8, 16]
N = 2000


def trial(rng, stride, sigma, noise):
    n = INPUT // stride
    err = {d: [] for d in DECODERS}
    for _ in range(N):
        kp = rng.uniform(64, 192, 2)                     # true keypoint, input pixels
        mu = kp / stride                                 # heatmap coordinates
        g = gaussian((n, n), mu, sigma)
        h = g + noise * rng.standard_normal((n, n))      # heatmap-trained output
        z = np.log(np.maximum(g, 1e-300)) + noise * rng.standard_normal((n, n))  # logits
        est = {"argmax": argmax(h), "quarter": quarter(h), "DARK": dark(h, sigma),
               "int-heatmap": soft_argmax_heatmap(h), "int-logits": soft_argmax_logits(z)}
        for d in DECODERS:
            err[d].append(np.linalg.norm(est[d] * stride - kp))
    return {d: np.mean(e) for d, e in err.items()}


rng = np.random.default_rng(0)
print("mean distance of a uniform point in a unit square to its centre:",
      f"{np.mean(np.linalg.norm(rng.uniform(-0.5, 0.5, (10**6, 2)), axis=1)):.4f}")

for label, sig in [("sigma = 2 heatmap cells at every stride", lambda s: 2.0),
                   ("sigma = 8 input pixels (2 cells at stride 4)", lambda s: 8.0 / s)]:
    for noise in (0.0, 0.05):
        print(f"\n{label}, noise std {noise}: mean error in input pixels")
        print("stride  map    sigma  " + "".join(f"{d:>12s}" for d in DECODERS))
        for s in STRIDES:
            r = trial(rng, s, sig(s), noise)
            print(f"{s:6d}  {INPUT // s:3d}^2  {sig(s):5.2f}  "
                  + "".join(f"{r[d]:12.2f}" for d in DECODERS))

mean distance of a uniform point in a unit square to its centre: 0.3826

sigma = 2 heatmap cells at every stride, noise std 0.0: mean error in input pixels
stride  map    sigma        argmax     quarter        DARK int-heatmap  int-logits
     2  128^2   2.00          0.77        0.38        0.00        0.00        0.00
     4   64^2   2.00          1.52        0.76        0.00        0.00        0.00
     8   32^2   2.00          3.06        1.54        0.00        0.00        0.00
    16   16^2   2.00          6.17        3.06        0.11        0.73        0.73

sigma = 2 heatmap cells at every stride, noise std 0.05: mean error in input pixels
stride  map    sigma        argmax     quarter        DARK int-heatmap  int-logits
     2  128^2   2.00          1.01        0.66        0.10       45.76        0.03
     4   64^2   2.00          1.98        1.31        0.20       37.48        0.05
     8   32^2   2.00          3.98        2.62        0.40       23.47        0.10
    16   16^2   2.00          8.04        5.32        0.81        9.68        0.79

sigma = 8 input pixels (2 cells at stride 4), noise std 0.0: mean error in input pixels
stride  map    sigma        argmax     quarter        DARK int-heatmap  int-logits
     2  128^2   4.00          0.77        0.38        0.00        0.00        0.00
     4   64^2   2.00          1.52        0.77        0.00        0.00        0.00
     8   32^2   1.00          3.06        1.53        0.00        0.00        0.00
    16   16^2   0.50          6.03        3.08        0.00        0.35        0.35

sigma = 8 input pixels (2 cells at stride 4), noise std 0.05: mean error in input pixels
stride  map    sigma        argmax     quarter        DARK int-heatmap  int-logits
     2  128^2   4.00          1.76        1.46        0.10       38.17        0.02
     4   64^2   2.00          2.03        1.32        0.20       37.94        0.05
     8   32^2   1.00          3.18        1.68        0.41       38.42        0.10
    16   16^2   0.50          6.11        3.24       27.70       39.56        0.40

The script, run with NumPy 2.3.5 and SciPy 1.16.3, takes about 8 seconds. It was run as first written; nothing was changed after the first run. Each cell averages 2,000 keypoints from one random stream (seed 0).

On clean maps the argmax error is 0.383 cells, the mean distance from a random point in a cell to the cell’s centre printed on the first line, and so it doubles with every doubling of the stride: from under a pixel at stride 2 to about six pixels at stride 16. The quarter shift halves it at every stride, consistent with Huang et al.’s per-axis analysis. DARK and both integral decoders are exact on a clean Gaussian, as expected from the mathematics: the logarithm of a sampled Gaussian is exactly quadratic, and the expectation of a symmetric bump is its centre. The exception is the $16 \times 16$ map, where keypoints lie at cell coordinates 4 to 12, as close as three to four cells to the edge; with $\sigma = 2$ cells the map’s edge cuts off part of the bump, which pulls the expectation inward (0.73 pixels), and DARK’s zero-padded smoothing has a smaller effect of the same kind (0.11 pixels). In our reading this border effect also applies to real crops, where joints near the edge of the box are common.

Noise separates the decoders. With $\sigma$ fixed at 2 cells, argmax and the quarter shift get worse, because noise can move the maximum to a neighbouring cell of nearly equal height; DARK’s error stays near 0.05 cells, growing in pixels only with the stride. Integral regression applied to the heatmap is useless here: every background cell contributes its positive noise to the expectation, and the result is pulled tens of pixels toward the centre of the map. It is worse on finer maps, which have more background cells, the reverse of every other decoder. This is an exaggerated form of the bias that Gu, Yang and Yao (2021) derived for the softmax. When the output is modelled instead as logits that the softmax turns into the Gaussian, so that background cells carry almost no probability, the expectation is the most accurate decoder in every noisy row, though at stride 16, where the border effect adds to its error, only marginally ahead of DARK. Which picture describes a trained integral network depends on its training, which is why the two are printed side by side; in our reading, the gap between them shows why integral models are trained through the softmax rather than applied to a heatmap trained with mean squared error.

Fixing $\sigma$ in input pixels instead changes two things. At stride 2 the Gaussian is 4 cells wide and its top is flat, so noise moves the argmax further (1.76 pixels against 1.01 with $\sigma = 2$ cells). At stride 16 it is half a cell wide: argmax and the quarter shift are barely affected, but DARK, which reads the map two cells away from the maximum for its second differences, breaks down under noise, presumably because a Gaussian of $\sigma = 0.5$ has fallen to $e^{-8}$ of its peak two cells away, far below the noise. The released configurations cited above keep $\sigma$ between one and three cells, away from both extremes. The example leaves out everything a trained network adds, such as misplaced or double peaks and systematic offsets, so its numbers isolate the error of decoding and are not predictions of benchmark accuracy.

Common Misconceptions

• “Heatmaps are probability distributions.” Maps trained with mean squared error against unit-height Gaussians are not normalized and can be negative; treating them as densities, as in the example’s clipped integral decoder, biases the result. Integral regression normalizes with a softmax, and its output is a distribution only in that sense.
• “Higher resolution is the only fix for quantization.” At a fixed network and resolution, decoding alone moved COCO AP by several points in DARK’s ablation, and replacing the encoding and decoding moved it by about one point in UDP’s (Zhang et al., 2020; Huang et al., 2020).
• “Integral regression has no quantization error, so it has no bias.” It avoids the grid but introduces a bias toward the map’s centre through the softmax (Gu, Yang and Yao, 2021), and truncation at the map border has the same effect.
• “Small preprocessing details do not matter.” UDP’s ablation shows a 2.3 AP difference from a one-pixel shift of flipped heatmaps that was not reported in the papers whose code used it (Huang et al., 2020); comparisons across papers should hold such details fixed.
• “Coordinate regression is inherently worse than heatmaps.” With a learned error density (RLE) or bins finer than the input pixel (SimCC), coordinate outputs matched or exceeded heatmap models of the same backbone in their authors’ COCO experiments (Li et al., 2021; Li et al., 2022).

Where It Is Used

• Human pose estimation. Top-down networks such as the stacked hourglass network, SimpleBaseline and HRNet output one heatmap per joint for each detected person; bottom-up methods such as OpenPose detect all joints in an image as heatmap peaks and group them, as described in human pose estimation. Benchmarks on MPII Human Pose and COCO report the effects of decoding on PCKh and keypoint AP.
• 3D pose. Integral regression extends to volumetric heatmaps whose expectation gives 3D joint positions (Sun et al., 2018).
• Facial landmarks. Heatmap regression is one of the main approaches to face alignment, with losses designed for it such as the Adaptive Wing loss (Wang, Bo and Li, 2019).
• Object detection by keypoints. CornerNet detects objects as pairs of corners, predicted as peaks in per-class corner heatmaps; because mapping heatmap locations back to the image loses precision, which matters most for small boxes, it also predicts offsets to adjust each corner (Law and Deng, 2018), the same pattern as G-RMI’s heatmaps and offsets.