inline float ssimWindowMean()

in cpp/testutils/src/comparison/Ssim.cpp [38:106]

• 48 lines of code
• 5 McCabe index (conditional complexity)

inline float ssimWindowMean(
    const image::pixel::Specification& spec,
    const std::vector<std::unique_ptr<image::Scanline>>& imgA,
    const std::vector<std::unique_ptr<image::Scanline>>& imgB,
    const int x,
    const int y) {
  // compute SSIM over all channels
  const auto numChannels = spec.numberOfComponents();

  // sumsA/B = sum of pixels in window by channel
  std::array<int, ssimMaxChannels> sumsA{};
  std::array<int, ssimMaxChannels> sumsB{};
  // sumSqsA/B = sum of pixels^2 in window by channel
  std::array<int, ssimMaxChannels> sumSqsA{};
  std::array<int, ssimMaxChannels> sumSqsB{};
  // sumProdAB = sum of pixel_A * pixel_B in window by channel
  std::array<int, ssimMaxChannels> sumProdAB{};

  // compute sumsA/B / sumSqsA/B / sumProdAB
  for (int i = y; i < y + ssimWindowSize; i += 1) {
    const auto* rowA = imgA[i].get();
    const auto* rowB = imgB[i].get();

    for (int j = x; j < x + ssimWindowSize; j += 1) {
      const auto* pixelA = rowA->dataAtPixel(j);
      const auto* pixelB = rowB->dataAtPixel(j);

      for (int c = 0; c < numChannels; c += 1) {
        const auto chanA = pixelA[c];
        const auto chanB = pixelB[c];
        sumsA[c] += chanA;
        sumsB[c] += chanB;
        sumSqsA[c] += chanA * chanA;
        sumSqsB[c] += chanB * chanB;
        sumProdAB[c] += chanA * chanB;
      }
    }
  }

  // the mean SSIM of all channels in the window
  float meanSsim = 0.0f;

  // compute channel SSIMs and fold into the mean SSIM
  for (int c = 0; c < numChannels; c += 1) {
    // compute channel means, variances, and covariance
    const float meanA = float(sumsA[c]) / ssimWindowArea;
    const float meanB = float(sumsB[c]) / ssimWindowArea;
    const float varA = float(sumSqsA[c]) / ssimWindowArea -
        float(sumsA[c] * sumsA[c]) / ssimWindowAreaSq;
    const float varB = float(sumSqsB[c]) / ssimWindowArea -
        float(sumsB[c] * sumsB[c]) / ssimWindowAreaSq;
    const float covAB = float(sumProdAB[c]) / ssimWindowArea -
        float(sumsA[c] * sumsB[c]) / ssimWindowAreaSq;

    // compute SSIM of the channel
    const float t1 = 2 * meanA * meanB + ssimC1;
    const float t2 = 2 * covAB + ssimC2;
    const float t3 = meanA * meanA + meanB * meanB + ssimC1;
    const float t4 = varA + varB + ssimC2;
    const float ssim = (t1 * t2) / (t3 * t4);

    // numerically stable mean
    // mu_n = (\sum_{i=1}^n x_i) / n
    //      = mu_{n-1} + (x_n + mu_{n-1}) / n
    meanSsim = meanSsim + (ssim - meanSsim) / (c + 1);
  }

  return meanSsim;
}