【16位RAW圖像處理五】任意位深16位圖像的中值模糊快速實(shí)現(xiàn)及其應(yīng)用。

　　在我博客里，也多次提到了中值模糊的優(yōu)化，比如以下兩篇文章：

　　【算法隨記三】小半徑中值模糊的急速實(shí)現(xiàn)（16MB圖7.5ms實(shí)現(xiàn)） + Photoshop中蒙塵和劃痕算法解讀。

　    任意半徑中值濾波（擴(kuò)展至百分比濾波器）O(1)時間復(fù)雜度算法的原理、實(shí)現(xiàn)及效果。

　　但是，這些都是這對8位圖像的優(yōu)化，也就是說圖像的色階最多只有256，如果把這個優(yōu)化算法直接擴(kuò)展到16位的RAW圖像，有以下幾個情況：

　　1、小半徑的（3*3 / 5*5）的快速實(shí)現(xiàn)，是完全可以借鑒8位的處理方式的，而且速度依舊是非常出色的。 

　　2、更大的半徑的，如果直接沿用8位的處理方式（上述第二篇參考文章），有以下幾個問題：

　?。?）對于10、12位的也許還可以，這種情況需要將直方圖分別分解為32及64個粗分及細(xì)分直方圖，然后進(jìn)行處理，這個時候相關(guān)的直方圖相加的計(jì)算量也會對應(yīng)的增加，內(nèi)部的中值的更新計(jì)算量也會相應(yīng)的增加，整體耗時要比8位的增加2倍及4倍以上。

　?。?）對于14位及16位的圖像，如果使用同樣的處理方式，細(xì)分及粗分直方圖的數(shù)量增加至128和256個，這個時候的增加的計(jì)算量非常客觀了，已經(jīng)不具備任何的性能優(yōu)勢，另外一個重要的問題是，相關(guān)的內(nèi)存分配可能會失敗，因?yàn)樵贠(1)時間復(fù)雜度的文章中，相關(guān)的內(nèi)存需要下如以下代碼所示：

 　　當(dāng)處理14位圖時，H_Fine需要分配128*128*Width*sizeof(unsigned short)字節(jié)大小的內(nèi)存，假定寬度為3072像素，這個尺寸在處理RAW數(shù)據(jù)時并不是很夸張的，那么H_Fine需要96MB的內(nèi)存，如果是16位，則增加到384MB，了解計(jì)算機(jī)內(nèi)存的都知道，分配這么大的連續(xù)內(nèi)存是很有可能失敗的。

　　當(dāng)然，內(nèi)存問題也許有其他的方案可以解決，但是核心的速度問題還是無法滿足實(shí)際的需求的。 

　　（3）如果是9位、11位、13位、15位這種圖像，那其實(shí)直接使用8位的方式，除了上述兩個問題外，還有就是粗分和細(xì)分直方圖如何取舍，也是個研究點(diǎn)，當(dāng)然，實(shí)際中我們很少遇到這樣的格式，而且即使遇到了，也可以把他們當(dāng)成向上一位的10/12/14/16格式處理，因?yàn)橹兄堤幚硭惴ㄊ遣粫黾有碌南袼刂担ㄔ谠瓐D中原先不存在的像素色階）的。

　　因此，我一直在尋找或者說構(gòu)思一個能夠快速處理unsigned short數(shù)據(jù)類型的通用的中值算法，不過這么多年也一直沒有找到答案。 

　　最近，來自安道爾（誰都知道怎么回事）微信朋友江向我咨詢8位的中值和自適應(yīng)中值能不能用到16位上，又讓我在這個問題上摸索了半個月，無意中在GIMP的中值模糊中找到了問題的答案。

　　我比較了幾個能處理16位圖像的中值，粗略的做了個耗時統(tǒng)計(jì)：

　　同樣大小的8位灰度圖，使用SSE優(yōu)化后的速度耗時約為0.13S。

　　后續(xù)補(bǔ)充： 后續(xù)用halcon測試同樣的數(shù)據(jù)，結(jié)果耗時40s，我真的感謝他，可能halcon根本不關(guān)系16位的圖，畢竟工業(yè)上基本都是8位的。 

　　測試時，我注意觀察了下CPU的使用情況，可以確認(rèn)PS、GIMP、ImageJ等都使用多線程，CPU使用率都接近100%，而OpenCv實(shí)測只支持濾波器為3*3或者5*5的16位圖像的中值，后面我看了相關(guān)說明，確實(shí)有如下的講法：

　　不過可以確認(rèn)的是16位的中值也可以是做到非?？焖俚?。既然GIMP也能做到那么快，那我們當(dāng)然可以參考他的實(shí)現(xiàn)方式。 

　　在GIMP源代碼的有關(guān)的文件夾里搜索median，可以在gegl-master\operations\common\找到這個median-blur.c這個文件。

　　源代碼分析：

　　GIMP的源代碼其實(shí)對普通用戶來說是不太優(yōu)化的，因?yàn)檫@畢竟是一個大工程，我們也不期望直接去運(yùn)行或者調(diào)試這個軟件的代碼，而是從一個代碼里去大概得猜測他的實(shí)現(xiàn)思路和想法。

　　我已經(jīng)基本吃透了這個代碼，這里就對他做一些簡單易懂的說明：

　　整體來說，這個median-blur也是基于直方圖的操作，先統(tǒng)計(jì)第一個位置的直方圖，然后依據(jù)直方圖計(jì)算出中值，并記錄對應(yīng)的中值統(tǒng)計(jì)累加量，對于下一個像素，更新直方圖，更新完成后不是直接用所有直方圖的信息來計(jì)算中值，而是依據(jù)前一次中值的信息，在其附近搜索新的中值，因?yàn)閷τ趫D像來說，不管是8位還是16位的，相鄰的位置中值其實(shí)不會差異很大，因此在舊的中值處搜索新的中值，將是一個效率很高的過程。

　　我們關(guān)注下這個過程中的三個函數(shù)histogram_modify_vals、histogram_get_median、histogram_update，對于他們做一個簡單的解釋：

 1 static inline void histogram_modify_val (Histogram  *hist,   const gint32 *src, gint  diff,  gint  n_color_components, gboolean  has_alpha)
 2 {
 3   gint alpha = diff;
 4   gint c;
 5   if (has_alpha)
 6     alpha *= hist->alpha_values[src[n_color_components]];
 7   for (c = 0; c < n_color_components; c++)
 8     {
 9       HistogramComponent *comp = &hist->components[c];
10       gint                bin  = src[c];
11       comp->bins[bin] += alpha;
12       /* this is shorthand for:
13        *
14        *   if (bin <= comp->last_median)
15        *     comp->last_median_sum += alpha;
16        *
17        * but with a notable speed boost.
18        */
19       comp->last_median_sum += (bin <= comp->last_median) * alpha;
20     }
21   if (has_alpha)
22     {
23       HistogramComponent *comp = &hist->components[n_color_components];
24       gint                bin  = src[n_color_components];
25       comp->bins[bin] += diff;
26       comp->last_median_sum += (bin <= comp->last_median) * diff;
27     }
28   hist->count += alpha;
29 }

　　這個histogram_modify_val是個修改直方圖值的一個內(nèi)聯(lián)函數(shù)，代碼一堆，實(shí)際的意思呢，其實(shí)也很簡單，就是如果要把某個位置的像素信息（像素值為bin）添加到現(xiàn)有的直方圖中，則alpha設(shè)置為1，對應(yīng)Bins[bin]就增加1，如果這個時候bin值小于或等于之前的中值last_median，則需要將last_median_sum增加1，而last_median_sum實(shí)際上是保存了之前第一次超過或等于中值時所有元素的總數(shù)量。

　　如果是要將某個位置的像素信息從現(xiàn)有直方圖中刪除掉，則則alpha設(shè)置為-1，對應(yīng)Bins[bin]就減少1，如果這個時候bin值于小于或等于之前的中值last_median，則需要將last_median_sum減少1。這樣就動態(tài)的保持了直方圖信息和中值累加值的更新。

　　而從更新后的直方圖中獲取新的中值則需要使用histogram_get_median函數(shù)。

 1 static inline gfloat histogram_get_median (Histogram *hist,  gint component,  gdouble  percentile)
 2 {
 3   gint                count = hist->count;
 4   HistogramComponent *comp  = &hist->components[component];
 5   gint                i     = comp->last_median;
 6   gint                sum   = comp->last_median_sum;
 7   if (component == hist->n_color_components)
 8     count = hist->size;
 9   if (count == 0)　　 return 0.0f;
10   count = (gint) ceil (count * percentile);
11   count = MAX (count, 1);
12   if (sum < count)
13     {
14       while ((sum += comp->bins[++i]) < count);
15     }
16   else
17     {
18       while ((sum -= comp->bins[i--]) >= count);
19       sum += comp->bins[++i];
20     }
21   comp->last_median     = i;
22   comp->last_median_sum = sum;
23   return comp->bin_values[i];
24 }

　　不要過分的在意代碼里的一些不知所謂的判斷和變量啊，我們只關(guān)注下第12行到第20行，這里sum變量表示上一次統(tǒng)計(jì)直方圖時，第一次超過或等于中值時所有元素的總數(shù)量，如果sum小于我們設(shè)定的中值統(tǒng)計(jì)停止值，則需要增加中值，直到他再次大于或等于停止值，如果sum已經(jīng)大于了停止值，則需要減少中值，直到他第一次小于停止值，但是此時，我們需要將得到的臨時中值增加1，同時sum也要增加對應(yīng)的數(shù)量，以便讓新的中值滿足大于或等于停止值的要求。 

　　這個代碼里充分體現(xiàn)了++i和i++的靈活運(yùn)用，不過我感覺我還是不要用這種方式書寫代碼，寧愿分開寫，也要讓理解變得更為簡單。

　　下面是histogram_update的過程，這個在GIMP里支持矩形、圓形、菱形的中值，我只貼出矩形部分的代碼：

static inline void histogram_update (Histogram  *hist, const gint32  *src, gint  stride,GeglMedianBlurNeighborhood  neighborhood, gint  radius,const gint  *neighborhood_outline,Direction  dir)
{
  gint i;
  switch (neighborhood)
    {
    case GEGL_MEDIAN_BLUR_NEIGHBORHOOD_SQUARE:
      switch (dir)
        {
          case LEFT_TO_RIGHT:
            histogram_modify_vals (hist, src, stride, -radius - 1, -radius, -radius - 1, +radius, -1);
            histogram_modify_vals (hist, src, stride, +radius, -radius, +radius, +radius, +1);
            break;
          case RIGHT_TO_LEFT:
            histogram_modify_vals (hist, src, stride,+radius + 1, -radius, +radius + 1, +radius, -1);
            histogram_modify_vals (hist, src, stride,-radius, -radius,-radius, +radius, +1);
            break;
          case TOP_TO_BOTTOM:
            histogram_modify_vals (hist, src, stride,-radius, -radius - 1, +radius, -radius - 1,-1);
            histogram_modify_vals (hist, src, stride,-radius, +radius,+radius, +radius, +1);
            break;
        }
      break;
}

　　所謂的直方圖更新，意思是當(dāng)計(jì)算完一個位置的像素中值后，我們根據(jù)當(dāng)前像素已經(jīng)統(tǒng)計(jì)好的直方圖信息，去利用相關(guān)重疊信息來獲取下一個位置的新的直方圖，在GIMP這個的代碼里，采用了一個更新策略，即先從左到右（LEFT_TO_RIGHT）更新，到一行像素的最后一個位置時，再從上到向下（TOP_TO_BOTTOM）更新，到下一行時，則從右到左（RIGHT_TO_LEFT）更新，處理到下一行的第一個元素是，再次從上到下更新，然后接著又是從左到右，如下圖所示，如此往復(fù)循環(huán)。

 　　　在更為具體的GIMP代碼里，我們還注意到GIMP還有很有意思的convert_values_to_bins函數(shù)，這個東西是個有點(diǎn)意思的玩意，他實(shí)際上是干啥呢，說白了就是減少冗余的信息，正常來說一個16位的圖像，那么他可能所含有的色階數(shù)最多就是有65536中，但是實(shí)際上一副圖里真正都用到的色階呢很有可能是不到65536個的，那這個現(xiàn)象有什么意義呢，他的核心就在于可以減少統(tǒng)計(jì)直方圖獲取中值的時間，具體來說，convert_values_to_bins就是把原始的圖像信息做適當(dāng)壓縮，使得每個色階都有至少一個值存在于圖像中。我們舉個例子，一副只有10個像素的16位圖，他們的值分別為：

　　　　　1  　　　　100　　　　2000　　　　150　　　　100　　　　40000　　　　350　　　　1　　　　2000　　　　300

　　這個時候如果我們定義的直方圖為Histgram[65536]，那么只有稀疏幾個色階有對應(yīng)的直方圖信息，大部分都為0，這樣我們要統(tǒng)計(jì)中值還是要從0開始循環(huán)，一次掃過一堆為0的直方圖信息，直到某個值為止符合中值的條件才停止。但是如果在進(jìn)行直方圖統(tǒng)計(jì)前已經(jīng)把圖像信息進(jìn)行過統(tǒng)計(jì)和重新賦值，則有可能極大的改進(jìn)這個統(tǒng)計(jì)過程。

　　具體如下操作，首先統(tǒng)計(jì)只10個數(shù)據(jù)有幾個不同的值，明顯有 1、100、150、300、350、2000、40000等7個不同的值，然后把不同大小的值按從小到大排序，再把原始10個數(shù)據(jù)修改為這些排序后的不同的值的索引，則新的10個值變?yōu)椋?/span>

　　　　0　　　　　　1　　　　　5　　　  　　2　　　　　1　　　　　　6　　　　　4　　　　0　　　　　5　　　　3　　　

　　這個時候直方圖的定義只需要敢為Histgram[7]就可以了，統(tǒng)計(jì)直方圖獲取中值的次數(shù)將大大減少。當(dāng)然這個時候獲得中值的值不是真正的像素值，而是一個索引，而根據(jù)這個索引結(jié)合上述排序后的數(shù)據(jù)，則就可以獲取真正的像素了。 

　　GIMP中為了上述效果，使用一個sort_input_values函數(shù)，對所有像素帶索引信息進(jìn)行排序，然后在獲取新的索引值，個人覺得這是個思路，但是其實(shí)針對圖像數(shù)據(jù)來說，可以不用這么復(fù)雜，完全可以根據(jù)直方圖信息來，我后面修改為如下簡單易理解的代碼？

    for (int Y = 0; Y < Width * Height; Y++)
    {
        Histgram[Src[Y]]++;
    }
    //    統(tǒng)計(jì)不同的色階數(shù)量,注意這里要搜索MaxV這個值
    for (int Y = 0; Y <= MaxV; Y++)
    {
        if (Histgram[Y] != 0)    BinAmount++;
    }
    //    分配合適的大小
    BinValue = (unsigned short*)malloc(BinAmount * sizeof(unsigned short));
    BinAmount = 0;
    for (int Y = 0; Y <= MaxV; Y++)
    {
        if (Histgram[Y] != 0)
        {
            Table[Y] = BinAmount;            //    通過索引Y(即像素實(shí)際值)能找到對應(yīng)的Bin
            BinValue[BinAmount] = Y;        //    通過Bin也能得到對應(yīng)的像素值
            BinAmount++;
        }
    }
    //    把Expand里的數(shù)據(jù)隱射到更小的范圍里，Expand是中間數(shù)據(jù)，可以隨意更改
    for (int Y = 0; Y < ExpandW * ExpandH; Y++)
    {
        Expand[Y] = Table[Expand[Y]];
    }

　　這個效率就比GIMP那個高很多了。

　　另外，仔細(xì)看GIMP的代碼，在其process函數(shù)里，還增加了分開處理的部分，核心部分如下所示：

 if (! data->quantize &&
      (roi->width > MAX_CHUNK_WIDTH || roi->height > MAX_CHUNK_HEIGHT))
    {
      gint n_x = (roi->width  + MAX_CHUNK_WIDTH  - 1) / MAX_CHUNK_WIDTH;
      gint n_y = (roi->height + MAX_CHUNK_HEIGHT - 1) / MAX_CHUNK_HEIGHT;
      gint x;
      gint y;
      for (y = 0; y < n_y; y++)
        {
          for (x = 0; x < n_x; x++)
            {
              GeglRectangle chunk;
              chunk.x      = roi->x + roi->width  * x       / n_x;
              chunk.y      = roi->y + roi->height * y       / n_y;
              chunk.width  = roi->x + roi->width  * (x + 1) / n_x - chunk.x;
              chunk.height = roi->y + roi->height * (y + 1) / n_y - chunk.y;

              if (! process (operation, input, output, &chunk, level))
                return FALSE;
            }
        }
      return TRUE;
    }

　　其中MAX_CHUNK_HEIGHT和MAX_CHUNK_WIDTH 定義為128。

　　這里也是很有意思的部分，我理解他至少有幾重意義：

　　1、分塊后更加適合多線程處理了。原始的工作方式，從左到右，從上到下，從右到左在從上到下更新直方圖，這個過程是前后依賴的，是不能直接使用多線程的，而分塊后，每一個塊之間可以做到獨(dú)立處理，當(dāng)然就可以線程并行了，代價時整體的計(jì)算量其實(shí)是增加的，但是耗時會變少。

　　2、前面說了GIMP需要通過排序來壓縮數(shù)據(jù)，減少直方圖的總量，但是如果不分快，一個整個圖，實(shí)際上由于數(shù)據(jù)量大，實(shí)際使用過的色階數(shù)還是比較多的，但是分成小塊之后，這個數(shù)據(jù)就有可能極大的減少，特別有一些RAW圖像常有的背景區(qū)域，基本上就幾十個色階，這種對于速度提升來說是很有幫助的。

　　3、另外一個問題就是，排序是個耗時的工作，如果對整幅圖的數(shù)據(jù)進(jìn)行排序，那么這個意義就很小了，但是我們知道一個事實(shí)，24*24次 128*128的排序，要比單次3072*3072數(shù)量的排序快很多的，雖然實(shí)際上分塊后的排序?qū)挾群透叨壬线€要增加半徑值，但是也還是比整體排序快很多，所以分塊默認(rèn)還帶來了這個好處。

　　所以這個代碼也是一環(huán)扣一環(huán)，當(dāng)然，如果用我上面那種處理方式，就不存在排序的事情了。

　　最后在說一點(diǎn)，就是前面那種直方圖從左到右，從上到下，從右到左的更新方式還有一個好處，就是直方圖的清零工作只需要做一次，而以前我寫的此類算法一般都是在每一行第一個點(diǎn)位置處清零，然后直接從左到右進(jìn)行更新計(jì)算，這種寫法對于不分塊整體處理來說可能影響不大，但是對于分塊的算法來說，如果還是這種做法，需要 n_x  * width清零，其中n_x 表示水平分塊的數(shù)量，比如3072*3072的，則需要24*3072清零，如果直方圖大小平均為10000個元素，這個操作也是有點(diǎn)占用時間的。所以這些好處都是潛在的需要炸取的。

　　進(jìn)過這一系列的操作，加上我自己的一些其他的優(yōu)化，目前，我能在相同配置的機(jī)器上做到和PS差不多或者更強(qiáng)的速度，比如同樣的測試圖，我做到了多線程版本85ms的速度（16位的）。

　　基于中值的實(shí)現(xiàn)呢，可以輔助實(shí)現(xiàn)一些其他的功能，比如基于中值的銳化，基于中值的去噪等等。因此，也是非常有意義的一項(xiàng)工作。

　　關(guān)于16位RAW圖像，本人開發(fā)了一個簡易的增強(qiáng)和處理程序，可在 https://files.cnblogs.com/files/Imageshop/Optimization_Demo_16.rar下載測試。

　　 如果想時刻關(guān)注本人的最新文章，也可關(guān)注公眾號：