src/gpu/graphite/sparse_strips/MSAA_LUT.h

/*
 * Copyright 2026 Google LLC
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */
#ifndef skgpu_graphite_sparse_strips_MSAA_LUT_DEFINED
#define skgpu_graphite_sparse_strips_MSAA_LUT_DEFINED

#include "src/gpu/graphite/sparse_strips/SparseStripsConfig.h"
#include "src/gpu/graphite/sparse_strips/SparseStripsTypes.h"

#include <array>
#include <cmath>
#include <cstdint>
#include <type_traits>

namespace skgpu::graphite {

/*
 * Naively, calculating MSAA coverage in software is O(n) per each sample location. Instead we use a
 * variation of a lookup technique introduced by Laine and Karras, adapted to 2D by Li et al., and
 * refined by the Vello project:
 * https://dl.acm.org/doi/abs/10.1111/j.1467-8659.2010.01728.x,
 * https://dl.acm.org/doi/10.1145/2980179.2982434)
 * https://github.com/linebender/vello/blob/main/vello_shaders/shader/fine.wgsl
 * https://github.com/linebender/vello/pull/64
 *
 * * To precompute the MSAA coverage of a line segment slicing through a 1x1 pixel, in Li et al. the
 *   line is parameterized by its normal vector (angle) and perpendicular distance from the center.
 *   While effective, this approach has drawbacks: angle and distance do not map linearly to the
 *   Cartesian pixel grid, leading to higher quantization errors (especially near corners), and
 *   computing the normal requires expensive square root operations (e.g., std::hypot) during
 *   rendering.
 * * Instead, we use a continuous parameterization introduced by the Vello project. We
 *   map the infinite range of standard slopes (m) and y-intercepts (b) into two new variables, s
 *   (normalized slope) and t (normalized translation), strictly bounded between 0.0 and 1.0.
 *
 * ----------------------------------
 * Squashing the Slope (m -> s)
 * ----------------------------------
 * * We squash the infinite positive slope m in [0, infinity) using the curve:
 *   s = 1 / (m + 1)
 * * This maps the extremes while avoiding division-by-zero risks:
 *   A perfectly vertical line (m = infinity) maps to s = 0.0
 *   A 45-degree line (m = 1) maps to s = 0.5
 *   A perfectly horizontal line (m = 0) maps to s = 1.0
 * * To use this in the half-plane equation (m*x - y + b >= 0), we solve for m:
 *   m = (1 - s) / s
 *   Substituting and multiplying by s gives:
 *   x(1 - s) - y*s + b*s >= 0
 *
 * ----------------------------------
 * Squashing the Translation (b -> t)
 * ----------------------------------
 * * To ensure our spatial offset stays between 0.0 and 1.0, we scale the y-intercept against our
 *   new slope parameter:
 *   t = (b + m) / (m + 1)
 * * We can rewrite t entirely in terms of s:
 *   t = (b + m) * s  =>  t = b*s + m*s
 * * Since m*s = 1 - s, we get:
 *   t = b*s + 1 - s
 * * Solving for b*s gives:
 *   b*s = s + t - 1
 *
 * ----------------------------------
 * Final Equation
 * ----------------------------------
 * * Substitute our new definition of b*s back into the half-plane equation:
 *   x(1 - s) - y*s + (s + t - 1) >= 0
 * * To group this cleanly, we add zero (+ t*s - t*s):
 *   x(1 - s) - y*s + s + t - 1 + t*s - t*s >= 0
 * * Rearranging to group by (1 - s) and s:
 *   x(1 - s) - (1 - s - t + t*s) - y*s + t*s >= 0
 * * Factoring the middle group into (1 - t)(1 - s), and factoring -s from the end:
 *   x(1 - s) - (1 - t)(1 - s) - (y - t)s >= 0
 * * Factoring out (1 - s) yields the final parameterization:
 *   (x - (1 - t))(1 - s) - (y - t)s >= 0
 *
 * ----------------------------------
 * LUT Generation
 * ----------------------------------
 * * For each subsample location (x,y), we evaluate the final half-plane equation at discretized
 *   points, according to the resolution of the LUT (configured via SparseStripConfig::kLUTMaskWidth
 *   and SparseStripConfig::kLUTMaskHeight). Empirical testing indicates that 64x64 is sufficient
 *   for almost all rendering scenarios.
 * * Since our mathematical derivation assumes a positive slope (m >= 0), we partition the LUT into
 *   two halves. The bottom half (v >= SparseStripConfig::kLUTMaskHeight / 2) stores masks for
 *   positive slopes. The top half stores masks for negative slopes.
 * * Because s and t are bounded by [0.0, 1.0], we quantize them into the discrete 2D grid. The grid
 *   has `SparseStripConfig::kLUTMaskWidth` columns representing the translation t, and
 *   `SparseStripConfig::kLUTMaskHeight / 2` rows representing the slope s.
 * * For negative slopes, we reuse the exact same mathematical equation but geometrically flip the
 *   Y-axis of our sub-pixel sample points (y = 1.0 - y).
 * * To minimize maximum quantization error, we extract the continuous s and t values from the exact
 *   center of each grid cell (e.g., `(u + 0.5) / SparseStripConfig::kLUTMaskWidthF`).
 * * For each cell, we evaluate the half-plane equation against all N sub-pixel sample coordinates
 *   (SparseStripConfig::kNumSubSamples). If the result is >= 0.0, the sample is covered, and we set
 *   the corresponding bit in our integer mask using a bitwise OR (`mask |= (1 << k)`).
 *
 * * The actual positions of the subsample points use the D3D11 standard multisample pattern
 *   (provided by kMsaaPattern in SparseStripsTypes.h):
 *   https://learn.microsoft.com/en-us/windows/win32/api/d3d11/ne-d3d11-d3d11_standard_multisample_quality_levels
 *
 * * For more on sampling patterns:
 *   https://web.cs.wpi.edu/~emmanuel/courses/cs563/S10/talks/wk3_p1_wadii_sampling_techniques.pdf
 */

inline SparseStripConfig::LUTArray GenerateMSAALUT() {
    constexpr float scale = 1.0f / static_cast<float>(SparseStripConfig::kNumSubSamples);
    std::array<float, SparseStripConfig::kNumSubSamples> subX;
    std::array<float, SparseStripConfig::kNumSubSamples> subY;
    const auto& pattern = kMsaaPattern<SparseStripConfig::SubSampleType>;

    for (size_t k = 0; k < SparseStripConfig::kNumSubSamples; ++k) {
        subX[k] = (static_cast<float>(pattern[k]) + 0.5f) * scale;
        subY[k] = (static_cast<float>(k) + 0.5f) * scale;
    }

    SparseStripConfig::LUTArray lut;
    constexpr int32_t halfHeight = SparseStripConfig::kLUTMaskHeight / 2;
    constexpr float halfHeightF = static_cast<float>(halfHeight);

    for (size_t i = 0; i < SparseStripConfig::kLUTSize; ++i) {
        int32_t u = static_cast<int32_t>(i % SparseStripConfig::kLUTMaskWidth);
        int32_t v = static_cast<int32_t>(i / SparseStripConfig::kLUTMaskWidth);

        bool isPos = v >= halfHeight;

        // Extract continuous parameters from the center of the grid cells
        float t = (static_cast<float>(u) + 0.5f) / SparseStripConfig::kLUTMaskWidthF;
        float s = (static_cast<float>(v % halfHeight) + 0.5f) / halfHeightF;

        SparseStripConfig::SubSampleType mask = 0;
        for (size_t k = 0; k < SparseStripConfig::kNumSubSamples; ++k) {
            float x = subX[k];
            float y = subY[k];

            if (!isPos) {
                y = 1.0f - y;
            }

            float val = (x - (1.0f - t)) * (1.0f - s) - (y - t) * s;

            if (val >= 0.0f) {
                mask |= (static_cast<SparseStripConfig::SubSampleType>(1) << k);
            }
        }
        lut[i] = mask;
    }
    return lut;
}

}  // namespace skgpu::graphite

#endif  // skgpu_graphite_sparse_strips_MSAA_LUT_DEFINED