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

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

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

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

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

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的實(shí)現(xiàn)函數(shù)，這是實(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;
  }
}

這里有個(gè)小技巧是，實(shí)現(xiàn)排序函數(shù)的時(shí)候因?yàn)槲覀冎皇菫榱饲笾兄?，所?strong>只需計(jì)算出前一半的有序元素即可，比如數(shù)組：

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

一般是全部排完得到：

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

中值就是34，但其實(shí)外部循環(huán)迭代只需要迭代到原來(lái)的一半（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

到這一步就可停止，因?yàn)橹兄狄呀?jīng)得到了，不過(guò)剛好這個(gè)例子是排到這一步就全部排好了而已。

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

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

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

因?yàn)橐话阄覀兲幚韴D像的數(shù)據(jù)像rgb類型的數(shù)據(jù)其起取值范圍是[0 ~ 255]，這時(shí)候其實(shí)完全不需要用float來(lái)存儲(chǔ)，用uint16_t類型就足夠了，中間計(jì)算也是全部用uint16_t替換，完整代碼：
https://github.com/Ldpe2G/ArmNeonOptimization/blob/master/ConstantTimeMedianFilter/src/normal_median_filter_uint16.cpp

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

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

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

這版優(yōu)化其實(shí)非常的暴力，就是既然每個(gè)窗口單次排序這樣子太慢，那么就利用并行計(jì)算一次同時(shí)計(jì)算8個(gè)窗口的排序結(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

從代碼上可以看到，因?yàn)橛玫氖莡int16_t類型的數(shù)據(jù)，所以可以一次處理8個(gè)窗口，相當(dāng)于把從左到右8個(gè)窗口內(nèi)的數(shù)據(jù)打包成C8的結(jié)構(gòu)，然后看下排序函數(shù)的改動(dòng)：

#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

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

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

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

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

ok，我們接著來(lái)看下這版的耗時(shí)：

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

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

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

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

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

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

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

之后來(lái)看下這一版的耗時(shí)：

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

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

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

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

A fast two-dimensional median filtering algorithm

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

審核編輯 黃昊宇