def _endpoint_heatmap(endpoints_zyx, shape, sigma):
    d, h, w = (int(s) for s in shape)
    heat = np.zeros((d, h, w), dtype=np.float32)

    rad = max(1, int(np.ceil(3 * sigma)))
    two_s2 = 2.0 * sigma * sigma

    for z0, y0, x0 in np.asarray(endpoints_zyx, float):
        zc, yc, xc = int(round(z0)), int(round(y0)), int(round(x0)))

        zl, zr = max(0, zc - rad), min(d, zc + rad + 1)
        yl, yr = max(0, yc - rad), min(h, yc + rad + 1)
        xl, xr = max(0, xc - rad), min(w, xc + rad + 1)

        zz, yy, xx = np.mgrid[
            zl:zr, yl:yr, xl:xr
        ].astype(np.float32)

        g = np.exp(
            -((zz-z0)**2 + (yy-y0)**2 + (xx-x0)**2) / two_s2
        )

        np.maximum(
            heat[zl:zr, yl:yr, xl:xr],
            g,
            out=heat[zl:zr, yl:yr, xl:xr]
        )

    return heat