优化 C++ 代码以匹配参考运行时间

Question

我的任务是优化一些 c++ 代码，我不擅长编码，但我做了一些尝试，所以原文是：

#include "stdafx.h"
#include "HistogramStretching.h"



void CHistogramStretching::HistogramStretching(BYTE** pImage, int nW, int nH)
{
    //find minimal value
    int nMin = pImage[0][0];
    for(int j = 0; j < nW; j++)
        for(int i = 0; i < nH; i++)
            if(pImage[i][j] < nMin)
                nMin = pImage[i][j];

    //find maximal value
    int nMax = pImage[0][0];
    for(int j = 0; j < nW; j++)
        for(int i = 0; i < nH; i++)
            if(pImage[i][j] > nMax)
                nMax = pImage[i][j];

    //stretches histogram
    for(int j = 0; j < nW; j++)
        for(int i = 0; i < nH; i++)
        {
            if(nMax != nMin)
            {
                float fScale = (nMax - nMin)/100.0;//calculates scale
                float fVal = (pImage[i][j] - nMin)/fScale;//scales pixel value
                int nVal = (int)(fVal + 0.5);//rounds floating point number to integer

                //checks BYTE range (must be 0-255)
                if(nVal < 0)
                    nVal = 0;
                if(nVal > 255)
                    nVal = 255;

                pImage[i][j] = nVal;
            }
            else
                pImage[i][j] = 0;//if all pixel values are the same, the image is changed to black
        }

}

我的版本是：

#include "stdafx.h"
#include "HistogramStretching.h"

void CHistogramStretching::HistogramStretching(BYTE** pImage, int nW, int nH)
{
    //find minimal value
    int nMin = pImage[0][0];
    int nMax = pImage[0][0];
    for (int j = 0; j < nW; j++) {
        for (int i = 0; i < nH; i++) {
            if (pImage[i][j] < nMin)
                nMin = pImage[i][j];
            if (pImage[i][j] > nMax)
                nMax = pImage[i][j];

        }
    }

    if (nMax != nMin) {
        float fScale = (nMax - nMin) / 100.0;//calculates scale
        fScale = 1 / fScale;
        //stretches histogram
        for (int j = 0; j < nW; j++)
            for (int i = 0; i < nH; i++)
            {


                float fVal = (pImage[i][j] - nMin) * fScale;//scales pixel value
                int nVal = (int)(fVal + 0.5);//rounds floating point number to integer

                //checks BYTE range (must be 0-255)
                if (nVal < 0)
                    nVal = 0;
                if (nVal > 255)
                    nVal = 255;

                pImage[i][j] = nVal;
            }
        //if all pixel values are the same, the image is changed to black

    }
    else {
        pImage[0][0] = 0;

    }
}

所以我将前两个循环合并为一个，但仍然是第一个 if 占用约 15% 的 CPU 时间，下一步是将 if 语句拉到循环之外并更改乘法除法，这里除法占用约 8% 的 CPU时间和浮点到 int 铸造大约需要 5%，但我认为我不能对铸造做太多事情。有了这个“更正”，我的代码仍然比参考代码慢 6-7 倍。我在同一台机器上测试这两个代码。你能指出我可以做得更好的事情吗？

Answer 1

我认为 tadman 给了你正确的答案。 替换

for (int j = 0; j < nW; j++) {
    for (int i = 0; i < nH; i++) {
        if (pImage[i][j] < nMin)
        ...
    }
}

与

for (int i = 0; i < nH; i++) {
    for (int j = 0; j < nW; j++) {
        if (pImage[i][j] < nMin)
        ...
    }
}

通过这种方式，您的数据访问变得缓存/内存对齐，这应该会更快。

Answer 2

所有现代编译器在完全优化编译时都可以很好地对其进行矢量化（MSVC 为 /O2，gcc 和 clang 为 -O3）。

这个想法是给编译器一些帮助，以便它可以看到代码实际上可以向量化：

1. 让内部循环在单个指针上运行，而不是在索引上运行，并且除了指向的值之外不访问任何内容。

2. 将缩放执行为整数运算 - 不要忘记四舍五入 :)

3. 尝试设置操作，以便不需要额外的范围检查，例如您对BYTE 的检查小于0。通过正确设置偏移量和比例，可以保证结果在所需范围内。

内部循环将展开，并将被矢量化为一次处理 4 个字节。我尝试了最近的 gcc、clang 和 MSVC 版本，它们为此生成了非常快的代码。

您正在做一些“奇怪”的事情，因为您故意将结果缩放到 0-99 范围。因此，您失去了数据的分辨率 - 您有一个完整的字节可供使用，那么为什么不将其缩放到 255？

但是，如果您想缩放到 100 个值，那也没关系。请注意，100(dec) = 0x64。我们可以使 outputSpan 灵活 - 它适用于任何值

因此：

/* Code Part 1 */
#include <cstdint>

constexpr uint32_t outputSpan = 100;

static constexpr uint32_t scale_16(uint8_t min, uint8_t max)
{
    return (outputSpan * 0x10000) / (1+max-min);
}
// scale factor in 16.16 fixed point unsigned format
// empty histogram produces scale = outputSpan
static_assert(scale_16(10, 10) == outputSpan * 0x10000, "Scale calculation is wrong");

static constexpr uint8_t scale_pixel(uint8_t const pixel, uint8_t min, uint32_t const scale)
{
    uint32_t px  = (pixel - min) * scale;
    // result in 16.16 fixed point format
    return (px + 0x8080u) >> 16;
    // round to an integer value
}

我们使用定点数（而不是浮点数）。比例为 16.16 格式，因此整数部分为 16 位，小数部分为 16 位，例如0x1234.5678。值 1.0(dec) 将是 0x1.0000。

像素缩放只是将像素乘以比例，四舍五入，然后返回截断的整数部分。

四舍五入是“有趣的”。您会认为将 0.5(dec) = 0x0.8 添加到结果中即可对其进行舍入。事实并非如此。该值需要比这大一点，0x0.808 可以完成这项工作。它对值进行预偏置，以使精确值周围的误差范围为零均值。在所有情况下，误差最多为 ±0.5 - 因此，四舍五入为整数的结果不会失去准确性。

我们使用scale_16和scale_pixel函数来实现担架：

/* Code Part 2 */
void stretchHistogram(uint8_t **pImage, int const nW, int const nH)
{
    uint8_t nMin = 255, nMax = 0;

    for (uint8_t **row = pImage, **rowEnd = pImage + nH; row != rowEnd; ++row)
        for (const uint8_t *p = *row, *pEnd = p + nW; p != pEnd; ++p)
        {
            auto const px = *p;
            if (px < nMin) nMin = px;
            if (px > nMax) nMax = px;
        }

    auto const scale = scale_16(nMin, nMax);

    for (uint8_t **row = pImage, **rowEnd = pImage + nH; row != rowEnd; ++row)
        for (uint8_t *p = *row, *pEnd = p + nW; p != pEnd; ++p)
            *p = scale_pixel(*p, nMin, scale);
}

这也可以在没有 FPU 的架构上生成不错的代码，例如无 FPU 的 ARM 和 AVR。

我们还可以进行一些手动检查。假设min = 0x10、max = 0xEF 和pixel = 0x32。请记住，scale 是 16.16 格式：

scale = 0x64.0000 / (1 + max - min)
      = 0x64.0000 / (1 + 0xEF - 0x10)
      = 0x64.0000 / (1 + 0xDF)
      = 0x64.0000 / 0xE0

长除法：

       0x  .7249
       0x64.0000 / 0xE0
       ---------
         64.0
       - 62.0
       ------
          2.00
       -  1.C0
       -------
           .400
       -   .380
       --------
           . 800
       -   . 7E0
       ---------
           .  20

所以，我们有scale = 0x0.7249。它小于 1 (0x1.0)，也小于 1/2 (0x0.8)，因为我们将 224 个值映射到 100 个值 - 不到一半。

现在

px = (pixel - min) * scale 
   = (0x32 - 0x10) * 0x0.7249
   =  0x22 * 0x0.7249

长乘法：

   0x 0.7249
*  0x  .0022
------------
       .E492
+     E.492
------------
   0x F.2DB2

因此，px = 0xF.2DB2 ≈ 0xF。我们必须将其四舍五入为整数：

return = (px + 0x0.8080u) >> 16
       = (0xF.2DB2 + 0x0.8080) >> 16
       =  0xF.AE32 >> 16
       ≈  0xF 

让我们检查一下十进制：

100 / (max-min+1) * (pixel-min) =
                                = 100 / (239 - 16 + 1) * (50 - 16)
                                = 100 / 224 * 34
                                = 100 * 34 / 224
                                = 3400 / 224
                                ≈ 15.17 
                                ≈ 15
                                ≈ 0xF

这是一个测试用例，可确保min、max 和输入像素值的所有组合都没有舍入偏差，并且误差限制为 [-0.5, 0.5]。只需将其附加到上面的代码中，它就会编译并运行并产生以下输出：

-0.5 0.5 1

要缩放到outputSpan = 256 值（而不是 100），它会输出：

-0.498039 0.498039 0.996078

/* Code Part 3 */

#include <cassert>
#include <cmath>
#include <iostream>

int main()
{
    double errMin = 0, errMax = 0;

    for (uint16_t min = 0; min <= 255; ++min)
        for (uint16_t max = min; max <= 255; ++max)
            for (uint16_t val = min; val <= max; ++val)
        {
            uint8_t const nMin = min, nMax = max;
            uint8_t const span = nMax - nMin;
            uint8_t const val_src = val;
            uint8_t p_val = val_src;
            uint8_t *const p = &p_val;
            assert(nMin <= nMax);
            assert(val >= nMin && val <= nMax);

            auto const scale = scale_16(nMin, nMax);
            *p = scale_pixel(*p, nMin, scale);

            auto pValTarget = (val_src - nMin) * 256.0/(1.0+span);
            auto error = pValTarget - *p;
            if (error < errMin) errMin = error;
            if (error > errMax) errMax = error;
        }

    std::cout << '\n' << errMin << ' ' << errMax << ' ' << errMax-errMin << std::endl;
    assert((errMax-errMin) <= 1.0); // constrain the error
    assert(std::abs(errMax+errMin) == 0.0); // constrain the error average
}