盒子濾波算是很基礎(chǔ)和經(jīng)典的函數(shù)，但是在PC上實(shí)現(xiàn)的話因?yàn)橛?a href="http://www.brongaenegriffin.com/tags/gpu/" target="_blank">GPU，借助其強(qiáng)大的算力所以可以很暴力的實(shí)現(xiàn)，每個(gè)thread計(jì)算以某點(diǎn)為中心給定半徑下的區(qū)域大小的和即可。那如果在移動(dòng)端cpu上如何寫高效的盒子濾波操作呢？
作者：梁德澎

最近一段時(shí)間做比較多移動(dòng)端開發(fā)相關(guān)的工作，感覺移動(dòng)端優(yōu)化相關(guān)的對(duì)我來說挺有趣的，以前都是在PC上寫代碼，寫代碼的時(shí)候?qū)τ诖a的性能沒有過多的思考和感覺。但是在移動(dòng)端上寫代碼明顯能察覺到一段代碼寫的好不好，對(duì)于在移動(dòng)端上運(yùn)行性能有很大的影響，尤其在一些比較老舊的機(jī)型上測(cè)試更能有感覺。

然后最近剛好在復(fù)現(xiàn)一篇論文，要在MXNet中實(shí)現(xiàn)類似盒子濾波（box filter）的操作子，其實(shí)就是步長為1的sum pooling，盒子濾波算是很基礎(chǔ)和經(jīng)典的函數(shù)，但是在PC上實(shí)現(xiàn)的話因?yàn)橛蠫PU，借助其強(qiáng)大的算力所以可以很暴力的實(shí)現(xiàn)，每個(gè)thread計(jì)算以某點(diǎn)為中心給定半徑下的區(qū)域大小的和即可。然后突發(fā)奇想想試試在移動(dòng)端cpu上試試如何寫高效的盒子濾波操作。

這篇文章就是把我的實(shí)踐過程記錄下來，首先給出最簡(jiǎn)單的實(shí)現(xiàn)然后如何一步步優(yōu)化，到最后給出一個(gè)性能優(yōu)化還不錯(cuò)的版本。由于我正式接觸移動(dòng)端優(yōu)化的時(shí)間不長，很多東西理解的不深，所以有哪里論述不正確的地方請(qǐng)讀者指出。

本文的代碼：
https://github.com/Ldpe2G/ArmNeonOptimization/tree/master/boxFilter

1.首先來看下Boxfilter最簡(jiǎn)單最直觀的實(shí)現(xiàn)

void BoxFilter::filter(float *input, int radius, int height, int width, float *output) {
  for (int h = 0; h < height; ++h) {
    int height_sift = h * width;
    for (int w = 0; w < width; ++w) {
      int start_h = std::max(0, h - radius);
      int end_h = std::min(height - 1, h + radius);
      int start_w = std::max(0, w - radius);
      int end_w = std::min(width - 1, w + radius);

      float tmp = 0;
      for (int sh = start_h; sh <= end_h; ++sh) {
        for (int sw = start_w; sw <= end_w; ++ sw) {
          tmp += input[sh * width + sw];
        }
      }
      output[height_sift + w] = tmp;
    }
  }
}

對(duì)每個(gè)點(diǎn)，計(jì)算給定半徑下區(qū)域的和，需要注意下邊界的處理。

其時(shí)間復(fù)雜度是 O( height x width x (radius x 2 + 1) x (radius x 2 + 1) )，

這個(gè)最簡(jiǎn)單的實(shí)現(xiàn)在輸入大小固定的情況下，半徑越大耗時(shí)越大，有很多重復(fù)計(jì)算的地方，相鄰元素在計(jì)算各自區(qū)域內(nèi)和的時(shí)候其實(shí)是有重疊的。然后第一個(gè)優(yōu)化的思路就是boxfilter的計(jì)算是行列可分離的，具體可參考[4]。

2.Boxfilter優(yōu)化第一版

void BoxFilter::fastFilter(float *input, int radius, int height, int width, float *output) {
  float *cachePtr = &(cache[0]);
  // sum horizonal
  for (int h = 0; h < height; ++h) {
    int sift = h * width;
    for (int w = 0; w < width; ++w) {
      int start_w = std::max(0, w - radius);
      int end_w = std::min(width - 1, w + radius);

      float tmp = 0;
      for (int sw = start_w; sw <= end_w; ++ sw) {
        tmp += input[sift + sw];
      }

      cachePtr[sift + w] = tmp;
    }
  }

  // sum vertical
  for (int h = 0; h < height; ++h) {
    int shift = h * width;
    int start_h = std::max(0, h - radius);
    int end_h = std::min(height - 1, h + radius);

    for (int sh = start_h; sh <= end_h; ++sh) {
      int out_shift = sh * width;
      for (int w = 0; w < width; ++w) {
        output[out_shift + w] += cachePtr[shift + w];
      }
    }
  }
}

所謂行列可分離就是，把行列分開計(jì)算，從代碼里可以看到，對(duì)每個(gè)元素，首先計(jì)算行方向上半徑內(nèi)的和，然后再計(jì)算列半徑內(nèi)的和，所以這時(shí)候的時(shí)間復(fù)雜度是O(height x width x (radius x 2 + 1) x 2)。

可以看到行列分離之后，時(shí)間復(fù)雜度減少了不少，尤其半徑越大減少的越多，但是還是有重復(fù)計(jì)算的地方。而且在固定輸入下時(shí)間復(fù)雜度還是會(huì)隨半徑的變大而變大。那么有沒有方法可以使得計(jì)算復(fù)雜度不受半徑的影響呢？優(yōu)化思路就是比如在算某一行每個(gè)點(diǎn)的半徑區(qū)域內(nèi)的和時(shí)，對(duì)于行開頭第一個(gè)點(diǎn)，首先計(jì)算其半徑內(nèi)和，然后對(duì)于接下來的點(diǎn)，不需要重新計(jì)算其半徑區(qū)域內(nèi)和，而是只需要把前一個(gè)元素半徑內(nèi)的和，按半徑窗口偏移之后減去舊的點(diǎn)和加上新加入的點(diǎn)即可。

3.Boxfilter優(yōu)化第二版

void BoxFilter::fastFilterV2(float *input, int radius, int height, int width, float *output) {
  float *cachePtr = &(cache[0]);
  // sum horizonal
  for (int h = 0; h < height; ++h) {
    int shift = h * width;

    float tmp = 0;
    for (int w = 0; w < radius; ++w) {
      tmp += input[shift + w];
    }

    for (int w = 0; w <= radius; ++w) {
      tmp += input[shift + w + radius];
      cachePtr[shift + w] = tmp;
    }

    int start = radius + 1;
    int end = width - 1 - radius;
    for (int w = start; w <= end; ++w) {
      tmp += input[shift + w + radius];
      tmp -= input[shift + w - radius - 1];
      cachePtr[shift + w] = tmp;
    }

    start = width - radius;
    for (int w = start; w < width; ++w) {
      tmp -= input[shift + w - radius - 1];
      cachePtr[shift + w] = tmp;
    }
  }

  float *colSumPtr = &(colSum[0]);
  for (int indexW = 0; indexW < width; ++indexW) {
    colSumPtr[indexW] = 0;
  } 
  // sum vertical
  for (int h = 0; h < radius; ++h) {
    int shift = h * width;
    for (int w = 0; w < width; ++w) {
      colSumPtr[w] += cachePtr[shift + w];
    }
  }

  for (int h = 0; h <= radius; ++h) {
    float *addPtr = cachePtr + (h + radius) * width;
    int shift = h * width;
    float *outPtr = output + shift; 
    for (int w = 0; w < width; ++w) {
      colSumPtr[w] += addPtr[w];
      outPtr[w] = colSumPtr[w];
    }
  }

  int start = radius + 1;
  int end = height - 1 - radius;
  for (int h = start; h <= end; ++h) {
    float *addPtr = cachePtr + (h + radius) * width;
    float *subPtr = cachePtr + (h - radius - 1) * width;
    int shift = h * width;
    float *outPtr = output + shift;
    for (int w = 0; w < width; ++w) {
      colSumPtr[w] += addPtr[w];
      colSumPtr[w] -= subPtr[w];
      outPtr[w] = colSumPtr[w];
    }
  }

  start = height - radius;
  for (int h = start; h < height; ++h) {
    float *subPtr = cachePtr + (h - radius - 1) * width;
    int shift = h * width;
    float *outPtr = output + shift; 
    for (int w = 0; w < width; ++w) {
      colSumPtr[w] -= subPtr[w];
      outPtr[w] = colSumPtr[w];
    }
  }
}

這一版時(shí)間復(fù)雜度大概是O(height x width x 4 )。不算邊界只看中間部分的計(jì)算就是一次加法和一次減法，行方向和列方向都一樣。這里行方向的部分很好理解，因?yàn)檫吔绮糠中枰厥馓幚?，比如開始部分只有加，結(jié)尾部分只有減法，所以計(jì)算分成了3部分。列方向計(jì)算的話按照常規(guī)思路，那就是按一列列來處理，可是我們知道數(shù)據(jù)一般是按照行來存儲(chǔ)的，這樣子跳行取數(shù)據(jù)，會(huì)造成很多次cache miss，這樣子性能肯定會(huì)受很大的影響，所以這里用了一個(gè)大小是width的向量colSum，來存儲(chǔ)每一列對(duì)應(yīng)點(diǎn)的半徑區(qū)域內(nèi)的和，然后遍歷的時(shí)候還是按照行來遍歷，如果一下子理解不了這個(gè)思路的話，可以想象如果width為1的情況，那么應(yīng)該可以更好的理解。

然后我們來看下實(shí)驗(yàn)結(jié)果，這三版boxfilter在輸入是2000x2000的情況下，在不同半徑下的運(yùn)行耗時(shí)，測(cè)試手機(jī)是華為榮耀4C（CHM-TL00），每個(gè)函數(shù)運(yùn)行10次取平均為其耗時(shí)：

可以看到第二版優(yōu)化的耗時(shí)在不同半徑下的表現(xiàn)都很穩(wěn)定，基本不受影響。然后接下來的優(yōu)化思路就是在確定了C++ 的代碼之后可以采用arm Neon Intrinsics來加速了，就是利用向量計(jì)算指令同時(shí)處理多個(gè)數(shù)據(jù)，把獨(dú)立的運(yùn)算同時(shí)做，比寫匯編要容易。

4.Boxfilter優(yōu)化第二版 Neon Intrinsics

int n = width >> 2;
  int re = width - (n << 2);

  int start = radius + 1;
  int end = height - 1 - radius;
  for (int h = start; h <= end; ++h) {
    float *addPtr = cachePtr + (h + radius) * width;
    float *subPtr = cachePtr + (h - radius - 1) * width;
    int shift = h * width;
    float *outPtr = output + shift; 
    int indexW = 0;
    float *tmpOutPtr = outPtr;
    float *tmpColSumPtr = colSumPtr;
    float *tmpAddPtr = addPtr;
    float *tmpSubPtr = subPtr;

    int nn = n;
    int remain = re;
#if __ARM_NEON
    for (; nn > 0; nn--) {
      float32x4_t _add = vld1q_f32(tmpAddPtr);
      float32x4_t _sub = vld1q_f32(tmpSubPtr);
      float32x4_t _colSum = vld1q_f32(tmpColSumPtr);

      float32x4_t _tmp = vaddq_f32(_colSum, _add);
      _tmp = vsubq_f32(_tmp, _sub);

      vst1q_f32(tmpColSumPtr, _tmp);
      vst1q_f32(tmpOutPtr, _tmp);

      tmpAddPtr += 4;
      tmpSubPtr += 4;
      tmpColSumPtr += 4;
      tmpOutPtr += 4;
    }
#endif // __ARM_NEON
    for (; remain > 0; --remain) {
      *tmpColSumPtr += *tmpAddPtr;
      *tmpColSumPtr -= *tmpSubPtr;
      *tmpOutPtr = *tmpColSumPtr;
      tmpAddPtr ++;
      tmpColSumPtr ++;
      tmpOutPtr ++;
      tmpSubPtr ++;
    }
  }

上面的代碼是截取列方向中間計(jì)算部分來展示如何使用arm Neon Intrinsics函數(shù)，完整代碼可以看
https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/boxFilter/src/boxFilter.cpp#L143
行方向是沒辦法并行的，因?yàn)橄噜徳赜幸蕾?。而列方向上則可以，所以在列方向上做neon加速。
以上代碼其實(shí)挺好理解的，vld1q/_f32指令就是加載4個(gè)浮點(diǎn)數(shù)，然后vaddq/_f32，為把兩個(gè)float32x4/_t向量相加，相當(dāng)于同時(shí)計(jì)算了4個(gè)輸出，然后再把結(jié)果用vst1q/_f32存回去對(duì)應(yīng)的地址，然后所有參與運(yùn)算的地址都是每次加4，具體可以參考官網(wǎng)文檔。

然后來看下這版優(yōu)化的耗時(shí)如何：

可以看到耗時(shí)又少了一點(diǎn)，但是收益已經(jīng)不大了。然后還想嘗試進(jìn)一步優(yōu)化把Intrinsics部分改寫成內(nèi)聯(lián)匯編試試。

5.Boxfilter優(yōu)化第二版 Neon Assembly

int n = width >> 2;
  int re = width - (n << 2);

  int start = radius + 1;
  int end = height - 1 - radius;
  for (int h = start; h <= end; ++h) {
    float *addPtr = cachePtr + (h + radius) * width;
    float *subPtr = cachePtr + (h - radius - 1) * width;
    int shift = h * width;
    float *outPtr = output + shift; 
    int indexW = 0;
    float *tmpOutPtr = outPtr;
    float *tmpColSumPtr = colSumPtr;
    float *tmpAddPtr = addPtr;
    float *tmpSubPtr = subPtr;

    int nn = n;
    int remain = re;
#if __ARM_NEON
    asm volatile(
      "0:                       /n"
      "vld1.s32 {d0-d1}, [%0]!  /n"
      "vld1.s32 {d2-d3}, [%1]!  /n"
      "vld1.s32 {d4-d5}, [%2]   /n"
      "vadd.f32 q4, q0, q2      /n"
      "vsub.f32 q3, q4, q1      /n"
      "vst1.s32 {d6-d7}, [%3]!  /n"
      "vst1.s32 {d6-d7}, [%2]!  /n"
      "subs %4, #1              /n"
      "bne  0b                  /n"
      : "=r"(tmpAddPtr), //
      "=r"(tmpSubPtr),
      "=r"(tmpColSumPtr),
      "=r"(tmpOutPtr),
      "=r"(nn)
      : "0"(tmpAddPtr),
      "1"(tmpSubPtr),
      "2"(tmpColSumPtr),
      "3"(tmpOutPtr),
      "4"(nn)
      : "cc", "memory", "q0", "q1", "q2", "q3", "q4"
    );

#endif // __ARM_NEON
    for (; remain > 0; --remain) {
      *tmpColSumPtr += *tmpAddPtr;
      *tmpColSumPtr -= *tmpSubPtr;
      *tmpOutPtr = *tmpColSumPtr;
      tmpAddPtr ++;
      tmpColSumPtr ++;
      tmpOutPtr ++;
      tmpSubPtr ++;
    }
  }

完整版代碼：https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/boxFilter/src/boxFilter.cpp#L331

這里我只對(duì)列計(jì)算中間部分做了改寫，neon匯編下面的"cc"，"memory"之后跟的寄存器，是為了告訴編譯器（主要是q開頭的，q和d是一樣的，q表示128位向量寄存器（16個(gè)），d表示64位（32個(gè)），q0 =（d0 + d1）），這些寄存器會(huì)在匯編內(nèi)被用到，然后編譯器在進(jìn)入這段代碼之前，要緩存這些寄存器的內(nèi)容，然后在離開這段匯編之后恢復(fù)原來的值。一定要記得寫上用了哪些向量寄存器。

簡(jiǎn)單解釋一下，指令的意思，"vld1.s32 {d0-d1}, [%0]! /n"，相當(dāng)?shù)扔趶膖mpAddPtr這個(gè)地址連續(xù)讀取4個(gè)浮點(diǎn)數(shù)到{d0-d1}也就是q0寄存器，浮點(diǎn)數(shù)每個(gè)32位，乘以四就是128位。最后的感嘆號(hào)表示，這個(gè)指令完成之后tmpAddPtr地址加4的意思，沒有就是不變。"vadd.f32 q4, q0, q2 /n" 就是把 q0和q2相加的結(jié)果放到q4，"vsub.f32 q3, q4, q1 /n" 就是把q4減去q1的結(jié)果放到q3，和上面的intrinsics指令對(duì)應(yīng)。
然后vst1.s32就是把寄存器的內(nèi)容存到tmpOutPtr和tmpColSumPtr地址指向的內(nèi)存。
最后的subs指令和bne相當(dāng)于for循環(huán)的功能，最后對(duì)nn減一然后bne判斷是否為0， 不為0則繼續(xù)循環(huán)跳到開頭0標(biāo)記出繼續(xù)執(zhí)行。

匯編指令其實(shí)和intrinsics函數(shù)有對(duì)應(yīng)的具體可參考官方文檔。

然后我們來看下耗時(shí)：

什么鬼，竟然還慢了，一定是我使用的方式不對(duì)。去查了下資料，看到這篇博客里面提到，指令vld和vst都是需要消耗兩個(gè)時(shí)鐘周期，其他指令基本都是一個(gè)時(shí)鐘周期，但是卻不意味著一個(gè)時(shí)鐘周期之后能立刻得到結(jié)果。那么看下來 vsub.f32 指令依賴 vadd.f32 的結(jié)果，所以白白浪費(fèi)了不少時(shí)鐘周期。而且現(xiàn)代的處理器支持雙發(fā)射流水線，也就意味著CPU可以同時(shí)拾取兩條數(shù)據(jù)無關(guān)指令，那么能否利用這點(diǎn)來更進(jìn)一步加速呢。

6.Boxfilter優(yōu)化第二版 Neon Assembly 第二版

int start = radius + 1;
  int end = height - 1 - radius;
  for (int h = start; h <= end; ++h) {
    float *addPtr = cachePtr + (h + radius) * width;
    float *subPtr = cachePtr + (h - radius - 1) * width;
    int shift = h * width;
    float *outPtr = output + shift; 
    int indexW = 0;
    float *tmpOutPtr = outPtr;
    float *tmpColSumPtr = colSumPtr;
    float *tmpAddPtr = addPtr;
    float *tmpSubPtr = subPtr;

    int nn = width >> 3;
    int remain = width - (nn << 3);
#if __ARM_NEON
    asm volatile(
      "0:                       /n"
      "pld      [%0, #256]      /n"
      "vld1.s32 {d0-d3}, [%0]!  /n"
      "pld      [%2, #256]      /n"
      "vld1.s32 {d8-d11}, [%2]  /n"

      "vadd.f32 q6, q0, q4      /n"

      "pld      [%1, #256]      /n"
      "vld1.s32 {d4-d7}, [%1]!  /n"

      "vadd.f32 q7, q1, q5      /n"

      "vsub.f32 q6, q6, q2      /n"

      "vsub.f32 q7, q7, q3      /n"

      "vst1.s32 {d12-d15}, [%3]!  /n"

      // 感謝 @隨風(fēng)漂 指出這里錯(cuò)誤，用錯(cuò)了寄存器，輸出結(jié)果是錯(cuò)的
      // "vst1.s32 {d16-d19}, [%2]!  /n" 

      "vst1.s32 {d12-d15}, [%2]!  /n"

      "subs %4, #1              /n"
      "bne  0b                  /n"
      : "=r"(tmpAddPtr), //
      "=r"(tmpSubPtr),
      "=r"(tmpColSumPtr),
      "=r"(tmpOutPtr),
      "=r"(nn)
      : "0"(tmpAddPtr),
      "1"(tmpSubPtr),
      "2"(tmpColSumPtr),
      "3"(tmpOutPtr),
      "4"(nn)
      : "cc", "memory", "q0", "q1", "q2", "q3", "q4", "q5", "q6", "q7", "q8", "q9"
    );

#endif // __ARM_NEON
    for (; remain > 0; --remain) {
      *tmpColSumPtr += *tmpAddPtr;
      *tmpColSumPtr -= *tmpSubPtr;
      *tmpOutPtr = *tmpColSumPtr;
      tmpAddPtr ++;
      tmpColSumPtr ++;
      tmpOutPtr ++;
      tmpSubPtr ++;
    }
  }

完整版代碼：https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/boxFilter/src/boxFilter.cpp#L527

可以看到這里的改進(jìn)思路就是，把兩條 vadd.f32 指令放一起，然后跟兩條vsub.f32，然后把加載 vsub.f32 要用到部分?jǐn)?shù)據(jù)指令 vld1.s32 放到兩個(gè) vadd.f32之間，同時(shí) vld1.s32 指令之前加上 pld 指令。這個(gè)指令為什么能加速我問了下做移動(dòng)端優(yōu)化的同事，pld把數(shù)據(jù)從內(nèi)存加載到cache然后下一條指令把數(shù)據(jù)從 cache加載到寄存器，如果不用pld，數(shù)據(jù)若不在cache中，那么就是需要直接從內(nèi)存加載到寄存器，這里會(huì)比前者慢很多。

然后我們來看下最終版的耗時(shí)：

看表格最終版的耗時(shí)比起最原始的實(shí)現(xiàn)至少可以加速6~7倍，肯定是還有更好的優(yōu)化方式，比如如果能對(duì)輸入做量化把float類型數(shù)據(jù)轉(zhuǎn)成8bit整型，那么就可以在單位時(shí)間處理更多數(shù)據(jù)，當(dāng)然量化到8bit上計(jì)算溢出的風(fēng)險(xiǎn)也會(huì)增大許多。有時(shí)候煉丹煉久了，學(xué)習(xí)下優(yōu)化也挺好玩的，感覺可以很好的鍛煉下思維和代碼能力，現(xiàn)在深度學(xué)習(xí)在移動(dòng)端應(yīng)用越來越廣泛，訓(xùn)出來的模型如果部署到移動(dòng)端之后運(yùn)行的效率很低那么也是白費(fèi)功夫。所以感覺對(duì)移動(dòng)端優(yōu)化有一定的了解對(duì)于如何設(shè)計(jì)對(duì)移動(dòng)端更友好的模型還是有幫助的。

更多AI移動(dòng)端優(yōu)化的請(qǐng)關(guān)注專欄嵌入式AI以及知乎（@梁德澎）。

審核編輯 黃昊宇