はじめに
NTSC加重平均法という有名なグレースケール変換がある． このグレースケール変換は，ある画素の輝度値 $Y$ を，次の式（\ref{eq:ntscAverage}）のように求めるものだ． （$r, g, b$ は赤，緑，青それぞれの値）
\begin{equation} Y = 0.298912 r + 0.586611 g + 0.114478 b \label{eq:ntscAverage} \end{equation}
単純な平均をとったグレースケール変換
\begin{equation} Y = \dfrac{r + g + b}{3} \end{equation}
と比較すると，NTSC加重平均法によるグレースケール変換は人間の赤，緑，青の感じ方を考慮したものになっており，自然なグレースケール化を行うことができる．
あるとき，このNTSC加重平均法によるグレースケール変換を，次のように整数のみで行っているプログラムを見ることがあった．
\begin{equation} Y = \dfrac{77 r + 150 g + 29 b}{256} \label{eq:ntscAverageOptimized08} \end{equation}
浮動小数点演算を整数演算に変換する発想に感動したので，この記事では，どのような過程で式（\ref{eq:ntscAverage}）から式（\ref{eq:ntscAverageOptimized08}）を得るのかについて解説する．
導出
記述の簡略化のため，
\begin{equation} W_r = 0.298912, \:\: W_g = 0.586611, \:\: W_b = 0.114478 \end{equation}
とおく． まず，理想的なNTSC加重平均法による輝度値を求める関数を以下の式（\ref{eq:ntscIdealFunction}）ように定める．
\begin{equation} F(r, g, b) = W_r r + W_g g + W_b b \label{eq:ntscIdealFunction} \end{equation}
これを計算機で計算し，最終的に得られる輝度値（小数点を切り捨てたもの）を求める関数を以下の式（\ref{eq:ntscFlooredFunction}）ように定める．
\begin{equation} F'(r, g, b) = \lfloor W_r r + W_g g + W_b b \rfloor \label{eq:ntscFlooredFunction} \end{equation}
ここで，関数 $F (r, g, b)$ を以下のように変形し， $f_N(r, g, b)$ を $F (r, g, b)$ の近似値を求める関数として定義する．
\begin{eqnarray} F(r, g, b) & = & F(r, g, b) \times \dfrac{2^N}{2^N} \nonumber \\ & = & \dfrac{2^N (W_r r + W_g g + W_b b)}{2^N} \nonumber \\ & \approx & \dfrac{\lfloor 2^N W_r + 0.5 \rfloor r + \lfloor 2^N W_g + 0.5 \rfloor g + \lfloor 2^N W_b + 0.5 \rfloor b}{2^N} \nonumber \\ & \equiv & f_N(r, g, b) \label{eq:ntscAverageApproxDerivation} \end{eqnarray}
$2^N$ を掛けているのは，2の累乗の計算はビット演算として表現可能なためである． ビット演算は掛け算や割り算と比べて，比較的高速に計算可能であり，容易に高速化が見込める．
さて， $N = 8$ のとき，すなわち， $f_8(r, g, b)$ は以下のようになる．
\begin{equation} f_8 (r, g, b) = \dfrac{77r + 150g + 29b}{256} \end{equation}
誤差についての考察
前述の導出より，式（\ref{eq:ntscAverageOptimized08}）がNTSC加重平均法の近似式であることがわかった． しかし，近似計算に誤差はつきものである． ここでは，$N$ によって，特に $N = 8$ のとき，どの程度の誤差が生じるのかについて調べてみる．
まず，記述の簡略化のために，
\begin{equation} w_{rN} = 2^N \times W_r, \:\: w_{gN} = 2^N \times W_g, \:\: w_{bN} = 2^N \times W_b \end{equation}
\begin{equation} w'_{rN} = \lfloor w_{rN} + 0.5 \rfloor, \:\: w'_{gN} = \lfloor w_{gN} + 0.5 \rfloor, \:\: w'_{bN} = \lfloor w_{bN} + 0.5 \rfloor \end{equation}
とおく． 前述の導出式（\ref{eq:ntscAverageApproxDerivation}）より，
\begin{equation} F (r, g, b) = f_N (r, g, b) + \dfrac{(w_{rN} - w'_{rN}) r + (w_{gN} - w'_{gN}) g + (w_{bN} - w'_{bN}) b}{2^N} \end{equation}
となる．右辺の第二項を理想の値との誤差を求める関数として，
\begin{equation} \epsilon_N (r, g, b) = \dfrac{(w_{rN} - w'_{rN}) r + (w_{gN} - w'_{gN}) g + (w_{bN} - w'_{bN}) b}{2^N} \end{equation}
とおく． ここに $N = 8$ を代入すると，
\begin{equation} \epsilon_8 (r, g, b) = \dfrac{-0.478528 r + 0.172416 g + 0.306368 b}{2^8} \end{equation}
を得る． $r, g, b$ の各係数の符号に注目すると，
\begin{equation} \epsilon_8 (255, 0, 0) \leq \epsilon_8 (r, g, b) \leq \epsilon_8 (0, 255, 255) \end{equation}
すなわち，
\begin{equation} -0.479 \leq \epsilon_8 (r, g, b) \leq 0.479 \end{equation}
を得る（小数第四位で四捨五入している）． これは，やや誤差が大きいので，実際に計算機に計算させるときは， $N = 10$ として，
\begin{equation} \epsilon_{10} (0, 255, 0) = -0.0773 \leq \epsilon_{10} (r, g, b) \leq 0.0775 = \epsilon_{10} (255, 0, 255) \end{equation}
程度に誤差を抑えるのがよいだろう． もちろん， $N$ をより大きくすれば，誤差をさらに抑えることができる．そのときには，オーバーフローに気をつけなくてはならない．
ここで，r, g, bの各ビット深度が8bitであり，途中の計算結果を32bit整数として扱うと考える．
\begin{equation} W_r + W_g + W_b = 1 \end{equation}
かつ，
\begin{equation} 0 \leq r, g, b \leq 255 \end{equation}
であるので，
\begin{equation} 0 \leq W_r r + W_g g + W_b b < 2^8 \end{equation}
したがって，
\begin{equation} 0 \leq 2^{24} (W_r r + W_g g + W_b b) < 2^{32} \label{eq:ntscIdealShift} \end{equation}
となり， $N = 24$ がオーバーフローを起こさず，誤差を最も小さくできるビットシフト量である...と結論付けたいところであるが，実際には，$2^N W_r, \:\: 2^N W_g, \:\: 2^N W_b$ をそれぞれ四捨五入した係数 $w'_{rN}, \:\: w'_{gN}, \:\: w'_{bN}$ を用いるので，$N = 24$ のとき，式（\ref{eq:ntscIdealShift}）が $2^{32}$ で抑えられるとは限らない． したがって， $N = 23$ が誤差を最も小さくできるビットシフト量であり，そのときの誤差 $\epsilon_{23} (r, g, b)$ は
\begin{equation} -2.061 \times 10^{-5} \leq \epsilon_{23} (r, g, b) \leq 2.025 \times 10^{-6} \end{equation}
となる．
まとめ
この記事では，NTSC加重平均法によるグレースケール変換を，整数演算のみで行う近似計算の導出を行い，その誤差について述べた．
具体的には，以下のような計算を行っているプログラム
y = (unsigned char) (0.298912 * r + 0.586611 * g + 0.114478 * b)
を
y = (unsigned char) ((306 * r + 601 * g + 117 + b) >> 10)
のように変更すると，元の計算とほとんど相違なく，高速に計算を行うことが可能となる．
グレースケール変換はOpenCVにも実装されているので，自分で実装する機会はほとんどないかもしれないが，整数 -> 整数という計算の途中に浮動小数点が介入する場合の高速化のアイデアとして，示唆に富むものがあるだろう．
