在復現(xiàn) Side window 中值濾波的時候就在思考中值濾波能怎么優(yōu)化，直觀上看中值濾波好像沒什么可優(yōu)化的點，因為中值濾波需要涉及到排序，而且半徑越大，排序的耗時也越大。那么中值濾波能否進一步加速呢？或者像均值濾波一樣，可以不受濾波半徑的影響呢？
作者：梁德澎

最近在復現(xiàn) Side window 中值濾波的時候就在思考中值濾波能怎么優(yōu)化，直觀上看中值濾波好像沒什么可優(yōu)化的點，因為中值濾波需要涉及到排序，而且半徑越大，排序的耗時也越大。那么中值濾波能否進一步加速呢？或者像均值濾波一樣，可以不受濾波半徑的影響呢？

答案是能！這篇博客就是記錄了我是怎么去優(yōu)化中值濾波的實踐過程。而前面的3小節(jié)都是介紹我自己嘗試的優(yōu)化思路，最后一節(jié)才是講本文標題提到的常量階時間復雜度中值濾波的實現(xiàn)思路，想直接看其實現(xiàn)思路的讀者可以跳到最后一小節(jié)。

1、一般中值濾波的實現(xiàn)

一開始能想到的中值濾波最直觀的實現(xiàn)就是，把每個濾波窗口的內(nèi)的值放進一個數(shù)組里面，然后排序，排序結(jié)果的排中間的值就是濾波結(jié)果。下面給出中值濾波的一般實現(xiàn)的示例代碼（下面展示的所有代碼只是為了用于說明，不保證能運行，實際代碼以github上的代碼為準）：

median_filter(const float  *input,
              const int     radius,
              const int     height,
              const int     width,
              float        *output) {

  int out_idx = 0;
  for (int h = 0; h < height; ++h) {
    const int h_lower_bound = std::max(0, h - radius);
    const int h_upper_bound = std::min(height - 1, h + radius);
    const int h_interval = h_upper_bound - h_lower_bound + 1;

    for (int w = 0; w < width; ++w) {
      const int w_left_bound = std::max(0, w - radius);
      const int w_right_bound = std::min(width - 1, w + radius);
      const int arr_len = h_interval * (w_right_bound - w_left_bound + 1);

      int idx = 0;
      for (int i = h_lower_bound; i <= h_upper_bound; ++i) {
        const int h_idx = i * width;
        for (int j = w_left_bound; j <= w_right_bound; ++j) {
          m_cache[idx ++] = input[h_idx + j];
        }
      }

      sortArr(m_cache.data(), arr_len);
      output[out_idx ++] = m_cache[arr_len / 2];
    }
  }
}

排序函數(shù)sortArr的實現(xiàn)函數(shù)，這是實現(xiàn)的是選擇排序法：

static void sortArr(float *arr, int len) {
  const int middle = len / 2;
  for (int i = 0; i <= middle; ++i) {
    float min = arr[i];
    int min_idx = i;
    for (int j = i + 1; j < len; ++j) {
      if (min > arr[j]) {
        min_idx = j;
        min = arr[j];
      }
    }
    // swap idx i and min_idx
    float tmp = arr[min_idx];
    arr[min_idx] = arr[i];
    arr[i] = tmp;
  }
}

這里有個小技巧是，實現(xiàn)排序函數(shù)的時候因為我們只是為了求中值，所以只需計算出前一半的有序元素即可，比如數(shù)組：

132, 45, 8, 1, 9, 100, 34

一般是全部排完得到：

1, 8, 9, 34, 45, 100, 132

中值就是34，但其實外部循環(huán)迭代只需要迭代到原來的一半（7 / 2）= 3 就行了就可停止了，下面看下選擇排序中間每一步結(jié)果：

第0步，1和132交換：

132, 45, 8, 1, 9, 100, 34 -> 1, 45, 8, 132, 9, 100, 34

第1步，8和45交換：

1, 45, 8, 132, 9, 100, 34 -> 1, 8, 45, 132, 9, 100, 34

第2步，9和45交換：

1, 8, 45, 132,9, 100, 34 -> 1, 8, 9, 132, 45, 100, 34

第3步，34和132交換：

1, 8, 9, 132, 45, 100, 34 -> 1, 8, 9, 34, 45, 100, 132

到這一步就可停止，因為中值已經(jīng)得到了，不過剛好這個例子是排到這一步就全部排好了而已。

然后看下這個最普通的實現(xiàn)在手機上的耗時，測試機型是華為P30（麒麟980），下面所有實驗設置輸入分辨率都是512x512，濾波半徑大小從1到5，耗時跑30次取平均：

可以看到性能很一般，而且隨著半徑增加耗時也急劇增加。下面來看下第一版的優(yōu)化，首先可以優(yōu)化的點就是計算的數(shù)據(jù)類型。

2、第一版優(yōu)化，float數(shù)據(jù)類型改uint16_t

因為一般我們處理圖像的數(shù)據(jù)像rgb類型的數(shù)據(jù)其起取值范圍是[0 ~ 255]，這時候其實完全不需要用float來存儲，用uint16_t類型就足夠了，中間計算也是全部用uint16_t替換，完整代碼：
https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/ConstantTimeMedianFilter/src/normal_median_filter_uint16.cpp

這樣子簡單改一下數(shù)據(jù)類型之后，我們來看下其耗時：

可以看到就是簡單改下運算數(shù)據(jù)類型，其運行耗時就可以下降不少。

3，第二版優(yōu)化，簡單利用并行計算指令

這版優(yōu)化其實非常的暴力，就是既然每個窗口單次排序這樣子太慢，那么就利用并行計算一次同時計算8個窗口的排序結(jié)果，下面是示例代碼：

#if defined(USE_NEON_INTRINSIC) && defined(__ARM_NEON)
    int neon_arr_len = h_interval * (w_end - w_start + 1) * 8;
    for (int w = w_second_loop_start; w < remain_start; w += 8) {
      const int w_left_bound = std::max(0, w + w_start);
      const int w_right_bound = std::min(width - 1, w + w_end);

      int idx = 0;
      for (int i = h_lower_bound; i <= h_upper_bound; ++i) {
        const int h_idx = i * width;
        for (int j = w_left_bound; j <= w_right_bound; ++j) {
          for (int q = 0; q < 8; ++q) {
            m_cache[idx ++] = input[h_idx + j + q];
          }
        }
      }

      sortC4ArrNeon(m_cache.data(), neon_arr_len);
      for (int i = 0; i < 8; ++i) {
        m_out_buffer[out_idx ++] = m_cache[(neon_arr_len / 8 / 2) * 8 + i];
      }
    }
#endif

完整代碼見：
https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/ConstantTimeMedianFilter/src/normal_median_filter_uint16.cpp#L102

從代碼上可以看到，因為用的是uint16_t類型的數(shù)據(jù)，所以可以一次處理8個窗口，相當于把從左到右8個窗口內(nèi)的數(shù)據(jù)打包成C8的結(jié)構(gòu)，然后看下排序函數(shù)的改動：

#if defined(USE_NEON_INTRINSIC) && defined(__ARM_NEON)
static void sortC4ArrNeon(uint16_t *arr, int len) {
  const int actual_len = len / 8;
  const int middle = actual_len / 2;
  uint16_t *arr_ptr = arr;
  for (int i = 0; i <= middle; ++i) {
    uint16x8_t  min = vld1q_u16(arr_ptr);
    uint16x8_t   min_idx = vdupq_n_u16(i);

    uint16_t *inner_arr_ptr = arr_ptr + 8;
    for (int j = i + 1; j < actual_len; ++j) {
      uint16x8_t curr =  vld1q_u16(inner_arr_ptr);
      uint16x8_t   curr_idx = vdupq_n_u16(j);
      uint16x8_t  if_greater_than = vcgtq_u16(min, curr);
      min     = vbslq_u16(if_greater_than, curr, min);
      min_idx = vbslq_u16(if_greater_than, curr_idx, min_idx);
      inner_arr_ptr += 8;
    }
    // swap idx i and min_idx
    for (int q = 0; q < 8; ++q) {
      float tmp = arr[min_idx[q] * 8 + q];
      arr[min_idx[q] * 8 + q] = arr[i * 8 + q];
      arr[i * 8 + q] = tmp;
    }
    arr_ptr += 8;
  }
}
#endif // __ARM_NEON

其實代碼上看主體框架改動不大，還是采用選擇排序法，不過如何利用neon intrinsic并行計算指令，同時對8個窗口內(nèi)的數(shù)據(jù)進行排序呢？借助 vcgtq 和 vbslq 這兩個指令就可以做到。

vcgtq 表示將第一個參數(shù)內(nèi)的數(shù)組元素與第二個參數(shù)對應元素比較，如果第一個數(shù)組的元素，大于等于對應第二個數(shù)組的對應元素，則結(jié)果對應位置會置為1，否則為0。

vbslq 指令有三個輸入，第一個輸入可以看做是判斷條件，如果第一個輸入的元素位置是1則結(jié)果的對應的位置就取第二個輸入的對應位置，否則從第三個輸入對應位置取值。其實這和mxnet的where操作子很像。

然后一次循環(huán)迭代完了之后，min_idx 數(shù)組就包含了這8個窗口當前迭代的各自最小值的位置。

ok，我們接著來看下這版的耗時：

可以看到用了neon加速之后，耗時減少了很多，大概是3~4倍的提速。

4，第三版優(yōu)化，算法上的改進

經(jīng)過前面的鋪墊，終于到了本文的重點部分。如何讓中值濾波的耗時不受濾波半徑的影響，其實本質(zhì)來說就是改變一下計算濾波窗口內(nèi)中值的思路，不再采用排序，而是采用統(tǒng)計直方圖的方式，因為一般圖像數(shù)據(jù)rgb取值范圍就是[0~255]，那么求一個窗口內(nèi)的的中值完全可以采統(tǒng)計這個窗口內(nèi)的長度是256的直方圖，然后中值就是從左到右遍歷直方圖，累加直方圖內(nèi)每個bin內(nèi)的值，當求和結(jié)果大于等于窗口內(nèi)元素個數(shù)的一半，那么這個位置的索引值就是這個窗口的中值。

不過這也不能解決濾波半徑增大的影響，那么如何去除半徑的影響呢，本文開頭提到的這篇“Median Filtering in Constant Time ”文章里面引入了列直方圖的方法，也就是除了統(tǒng)計濾波窗口的直方圖，還對于圖像的每一列，都初始化一個長度是256的直方圖，所以濾波圖像太寬的話需要的內(nèi)存消耗也會更多。

然后不考慮邊界部分，對于中間部分的濾波窗口，其直方圖不需要重新統(tǒng)計，只需要減去移出窗口的列直方圖，然后加上新進來的列直方圖即可，然后再計算中值，這三步加起來時間復雜度不會超過O(256*3)，不受濾波半徑影響，所以在行方向上是常量階時間復雜度。

然后列方向就是同樣的，列直方圖在往下一行移動的時候也是采用同樣方法更新，減去上一行和加上下一行的值，然后這樣子列方向上也不受濾波半徑影響了。

論文里采用的計算方式，當從左到右濾波的時候，第一次用到列直方圖的時候才去更新列直方圖，而我在實現(xiàn)的時候是移動到新的一行從頭開始濾波之前，首先更新所有的列直方圖，然后再計算每個濾波窗口的中值。而且我在申請直方圖緩存的時候是所有直方圖都放在同一段緩存內(nèi)。

之后來看下這一版的耗時：

可以看到耗時很穩(wěn)，基本不受濾波半徑影響，不過由于需要涉及到直方圖的計算，在濾波窗口比較小的時候，這個算法相對于直接算是沒有優(yōu)勢的，但是當濾波窗口大于等于3的時候，其優(yōu)勢就開始體現(xiàn)了，而且越大越有優(yōu)勢。

論文里還提到了其他的優(yōu)化思路，比如下面這篇文章的對直方圖分級，不過我目前還沒看懂怎么做，這個先挖個坑吧，等以后有機會再深挖:)。

A coarse-to-fine algorithm for fast median filtering of image data with a huge number of levels

?
還有在計算中值的時候其實是不需要每次都從頭開始遍歷直方圖來計算中值的，下面這篇論文介紹了一個計算技巧可以減少計算中值的時間，有興趣的讀者可以看下:

A fast two-dimensional median filtering algorithm

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

審核編輯 黃昊宇