概要

誤差拡散で画像を２値化すると、画像が明るく見える。
これはガンマが考慮されていないために発生する現象である。

近年、2値表示タイプのOLEDの普及(笑)により、画像を綺麗に表示したい場面がしばしば発生している。そのため、綺麗に見れるようにする。
(主な普及先は、RaspberryPiやら、BeagleBoneやら、Arduinoやら、NanoPiである)

主に↓のように、画像表示するときに必要となる。

や

既存の問題点

PythonではPILライブラリを使用することで、簡単に誤差拡散法での２値画像を生成することができる。

例：image = Image.open(cover_path).convert('1')

しかし、これには２つほど問題点がある。

RGB→Grayscale化が、Bt.601系で行われる
誤差拡散による２値化が、ガンマ値考慮されない

最近のディスプレイはsRGB(Bt.709)系であるため、Bt.601の係数でグレースケール化すると微妙な違和感が生じる。
誤差拡散は、ガンマ値考慮されないため、明るく出てしまう。

対策内容

floatは基本遅いので、すべて固定小数点の整数演算で実施する。

(1) RGB → Grayscale変換に、Bt.709の係数を使用

Bt.709での変換は下記の通り。

$ Luminance = R * 0.2126 + G * 0.7152 + B * 0.0722 $

これに65536掛けて固定小数点で扱う。
つまり、こうなる。

$ Luminance = (R * 13933 + G * 46871 + B * 4732) / 65536 $

ちなみに、Bt.601系は下記の通り。
OpenCVやPILのGrayscale化はこっちの式。おそらくJPEGのYUV―RGB変換がそうだから、こうなのであろう。

$ Luminance = R * 0.299 + G * 0.587 + B * 0.114 $

(2) 誤差拡散は、リニア空間で実施

リニア空間に戻してから誤差拡散の計算してあげる。
リニア空間に戻す際に使用するガンマ値はsRGBへの近似値2.2を使用し、リニア空間値は固定少数点20bit表現。
つまり、ガンマ値変換入れて、255が1048575になり、1が5にマップされます。
そして、1を1にマップすることが最低条件になるため、Gamma=2.2の場合 18bitは必要となるわけです。

※sRGBのガンマカーブは、機器実装の簡略化のため、暗部部が直線表現で下記のように規格化されています。

www.color.org/srgb.pdf

If R, G, B are less than or equal to 0.04045
  RL = R/12.92
  GL = G/12.92
  BL = B/12.02

If R, G, B are greater than 0.04045
  RL = ((R + 0.055)/1.055)2.4
  GL = ((R + 0.055)/1.055)2.4
  BL = ((R + 0.055)/1.055)2.4

なので、$ (1 / 255) / 12.92 * 4096 = 1.2432 $ となるので、ひとまず12bitあれば足りるわけである。
ディスプレイとかの 12bitLUT!! などを謳っているものは、実は最低条件なわけです。
※そして上記規格書の式、間違いが3つある。

(3) 誤差拡散は、FloydSteinbergを使用

FloydSteinberg

-	-	*	7/16	-
-	3/16	5/16	1/16	-

Python版コード

PILのimageを入れれば、PILのimageを返す。

import math
from PIL import Image

def ImageHalftoning_FloydSteinberg(image):

    shift   = 20;

    cx, cy = image.size;

    temp    = Image.new('I', (cx, cy));
    result  = Image.new('L', (cx, cy));

    tmp     = temp.load();
    dst     = result.load();

    # Setup Gamma tablw
    gamma   = [0]*256;
    for i in range(256):
        gamma[i]    = int( math.pow( i / 255.0, 2.2 ) * ((1 << shift) - 1) );

    # Convert to initial value
    if image.mode == 'L':
        src     = image.load();
        for y in range(cy):
            for x in range(cx):
                tmp[(x,y)]  =  gamma[ src[(x,y)] ];

    elif image.mode == 'RGB':
        src     = image.load();
        for y in range(cy):
            for x in range(cx):
                R,G,B   = src[(x,y)];
                Y       = (R * 13933 + G * 46871 + B * 4732) >> 16; # Bt.709
                tmp[(x,y)]  =  gamma[ Y ];

    elif image.mode == 'RGBA':
        src     = image.load();
        for y in range(cy):
            for x in range(cx):
                R,G,B,A = src[(x,y)];
                Y       = (R * 13933 + G * 46871 + B * 4732) >> 16; # Bt.709
                tmp[(x,y)]  =  gamma[ Y ];

    else:
        raise ValueError('Image.mode is not supported.')    

    # Error diffuse
    for y in range(cy):
        for x in range(cx):
            c   = tmp[(x,y)];
            e   = c if c < (1 << shift) else (c - ((1 << shift) - 1));

            dst[(x,y)]  = 0 if c < (1 << shift) else 255;

            # FloydSteinberg
            #   -       *       7/16
            #   3/16    5/16    1/16
            if  (x+1) < cx :
                tmp[(x+1,y)]    += e * 7 / 16;

            if (y+1) < cy :
                if 0 <= (x-1) :
                    tmp[(x-1,y+1)]  += e * 3 / 16;

                tmp[(x,y+1)]        += e * 5 / 16;

                if (x+1) < cx :
                    tmp[(x+1,y+1)]  += e * 1 / 16;

    return  result;

C++版のコード(RGB→L変換なし)は、下記参照。
NanoPi-NEOでOpenCV CapしてOLED表示 - Qiita

最後に

本Pythonコードは、下記への適用のために作成したものでした。
コードはgitにも登録済みです。

NanoPi-NEOとMPDとOLEDで音楽再生サーバ - Qiita