盒子濾波算是很基礎和經典的函數，但是在PC上實現的話因為有GPU，借助其強大的算力所以可以很暴力的實現，每個thread計算以某點為中心給定半徑下的區(qū)域大小的和即可。那如果在移動端cpu上如何寫高效的盒子濾波操作呢？
作者：梁德澎

最近一段時間做比較多移動端開發(fā)相關的工作，感覺移動端優(yōu)化相關的對我來說挺有趣的，以前都是在PC上寫代碼，寫代碼的時候對于代碼的性能沒有過多的思考和感覺。但是在移動端上寫代碼明顯能察覺到一段代碼寫的好不好，對于在移動端上運行性能有很大的影響，尤其在一些比較老舊的機型上測試更能有感覺。

然后最近剛好在復現一篇論文，要在MXNet中實現類似盒子濾波（box filter）的操作子，其實就是步長為1的sum pooling，盒子濾波算是很基礎和經典的函數，但是在PC上實現的話因為有GPU，借助其強大的算力所以可以很暴力的實現，每個thread計算以某點為中心給定半徑下的區(qū)域大小的和即可。然后突發(fā)奇想想試試在移動端cpu上試試如何寫高效的盒子濾波操作。

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

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

1.首先來看下Boxfilter最簡單最直觀的實現

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;
    }
  }
}

對每個點，計算給定半徑下區(qū)域的和，需要注意下邊界的處理。

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

這個最簡單的實現在輸入大小固定的情況下，半徑越大耗時越大，有很多重復計算的地方，相鄰元素在計算各自區(qū)域內和的時候其實是有重疊的。然后第一個優(yōu)化的思路就是boxfilter的計算是行列可分離的，具體可參考[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];
      }
    }
  }
}

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

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

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];
    }
  }
}

這一版時間復雜度大概是O(height x width x 4 )。不算邊界只看中間部分的計算就是一次加法和一次減法，行方向和列方向都一樣。這里行方向的部分很好理解，因為邊界部分需要特殊處理，比如開始部分只有加，結尾部分只有減法，所以計算分成了3部分。列方向計算的話按照常規(guī)思路，那就是按一列列來處理，可是我們知道數據一般是按照行來存儲的，這樣子跳行取數據，會造成很多次cache miss，這樣子性能肯定會受很大的影響，所以這里用了一個大小是width的向量colSum，來存儲每一列對應點的半徑區(qū)域內的和，然后遍歷的時候還是按照行來遍歷，如果一下子理解不了這個思路的話，可以想象如果width為1的情況，那么應該可以更好的理解。

然后我們來看下實驗結果，這三版boxfilter在輸入是2000x2000的情況下，在不同半徑下的運行耗時，測試手機是華為榮耀4C（CHM-TL00），每個函數運行10次取平均為其耗時：

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

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 ++;
    }
  }

上面的代碼是截取列方向中間計算部分來展示如何使用arm Neon Intrinsics函數，完整代碼可以看
https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/boxFilter/src/boxFilter.cpp#L143
行方向是沒辦法并行的，因為相鄰元素有依賴。而列方向上則可以，所以在列方向上做neon加速。
以上代碼其實挺好理解的，vld1q/_f32指令就是加載4個浮點數，然后vaddq/_f32，為把兩個float32x4/_t向量相加，相當于同時計算了4個輸出，然后再把結果用vst1q/_f32存回去對應的地址，然后所有參與運算的地址都是每次加4，具體可以參考官網文檔。

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

可以看到耗時又少了一點，但是收益已經不大了。然后還想嘗試進一步優(yōu)化把Intrinsics部分改寫成內聯匯編試試。

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

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

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

匯編指令其實和intrinsics函數有對應的具體可參考官方文檔。

然后我們來看下耗時：

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

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"

      // 感謝 @隨風漂 指出這里錯誤，用錯了寄存器，輸出結果是錯的
      // "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

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

然后我們來看下最終版的耗時：

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

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

審核編輯 黃昊宇