RGB24の1677万7216色のうち、YUVが表せる色の数は？

※2012/03/09追記
　　　この記事での計算方法はあまり適切ではなかったかもしれません。
　　　追加記事を書きましたのでそちらもご覧ください。

　　　　　1677万7216色の画像をAvisynth 2.6でYV24にしてからRGB24に戻して色を数えてみた

　　　　　YUV(8bit～10bit)で表せる色数の再計算

---

　
「RGB24の1677万7216色のうち、
　8bit-depthのYUVで表せるのは、最大でも440万色程度」

この計算を行ったChikuzen氏のブログ記事に触発され、自分でもちょっと計算してみました。
まず最初にChikuzen氏の記事を読んだうえでご覧になってください。

　　　Ch's barn: ConvertToRGB　（Chikuzen氏のブログ記事）

　
Chikuzen氏の記事では8bit-depthの計算を行っていたので、
こちらではそれに加えて9bit-depthと10bit-depthについても計算をしてみました。

ただし、正直なところ、こんな計算でいいのか不安です。
記事の最後にソースコードも示していますので、間違ってるとこがあればご指摘ください。
プログラミングも初心者レベルです。

さて、まずは計算結果から。

　
 8bit-Rec601: 2955936
 8bit-Pc.601: 4262360
 8bit-Rec709: 3046424
 8bit-PC.709: 4400226
 9bit-Rec601: 15831400
 9bit-Pc.601: 16713229
 9bit-Rec709: 16149193
 9bit-PC.709: 16777216
10bit-Rec601: 16777216
10bit-PC.601: 16777216
10bit-Rec709: 16777216
10bit-PC.709: 16777216

　　■8bit-depthについては、おおまかな数値はChikuzenさんの計算結果と一致しています。
　　　　フルレンジでも最大440万色。実際にはTVレンジで扱うことが多いので実質300万色前後というところですね。

　　■9bit-depthについては、RGB24の全ての色を表現しうるのはPC.709だけという結果となりました。
　　　　とはいえ、Rec601でも1583万色を表現できますので実用上はほぼ問題ないレベルでしょうか。

　　■10bit-depthについては、どのパターンでもRGB24の全ての色を表現できるという結果となりました。

　
計算に使ったソースコードを以下に示します。表示に「SyntaxHighlighter 3.0」を使っています。
下の「expand source」の部分をクリックするとソースコードが展開表示されます。
ソース部分をダブルクリックすると全選択され、コピーできます。

/* yuvcolor_calc.cpp */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define CLIP_Y(Y) Y = (Y < 0) ? 0.0 : (Y > 1.0) ? 1.0 : Y
#define CLIP_C(C) C = (C < -0.5) ? -0.5 : (C > 0.5) ? 0.5 : C
#define SATURATE(X) X = (X < 0) ? 0 : (X > 255) ? 255 : X

/* RGB24の色数（256*256*256色） */
#define RGB_TABLE_SIZE 16777216

/* H.264のRound()がこんな感じなので */
int myRound(double x)
{
	return (x >= 0) ? (int)floor(x + 0.5) : (int)(-floor(abs(x) + 0.5));
}

int main(void)
{
	const struct {
		char *name;		// 色空間の名称
		int bitDepth;	// ビット深度
		int y_min;		// 輝度Yの最小値
		int y_max;		// 輝度Yの最大値
		int uv_min;		// 色差の最小値
		int uv_max;		// 色差の最大値
		double r_cr;	// YUV→RGB変換の係数１
		double g_cb;	// YUV→RGB変換の係数２
		double g_cr;	// YUV→RGB変換の係数３
		double b_cb;	// YUV→RGB変換の係数４
	} matrix[] = {
		{" 8bit-Rec601",  8, 16,  235, 16,  240, 1.402 , -0.344 , -0.714 , 1.772 },
		{" 8bit-Pc.601",  8,  0,  255,  0,  255, 1.402 , -0.344 , -0.714 , 1.772 },
		{" 8bit-Rec709",  8, 16,  235, 16,  240, 1.5748, -0.1873, -0.4681, 1.8556},
		{" 8bit-PC.709",  8,  0,  255,  0,  255, 1.5748, -0.1873, -0.4681, 1.8556},
		{" 9bit-Rec601",  9, 32,  470, 32,  480, 1.402 , -0.344 , -0.714 , 1.772 },
		{" 9bit-Pc.601",  9,  0,  511,  0,  511, 1.402 , -0.344 , -0.714 , 1.772 },
		{" 9bit-Rec709",  9, 32,  470, 32,  480, 1.5748, -0.1873, -0.4681, 1.8556},
		{" 9bit-PC.709",  9,  0,  511,  0,  511, 1.5748, -0.1873, -0.4681, 1.8556},
		{"10bit-Rec601", 10, 64,  940, 64,  960, 1.402 , -0.344 , -0.714 , 1.772 },
		{"10bit-Pc.601", 10,  0, 1023,  0, 1023, 1.402 , -0.344 , -0.714 , 1.772 },
		{"10bit-Rec709", 10, 64,  940, 64,  960, 1.5748, -0.1873, -0.4681, 1.8556},
		{"10bit-PC.709", 10,  0, 1023,  0, 1023, 1.5748, -0.1873, -0.4681, 1.8556},
		{NULL}
	};

	bool *rgbTable = (bool *)malloc(RGB_TABLE_SIZE*sizeof(bool));
	if(!rgbTable){
		fprintf(stderr,"malloc(rgbTable) failed\n");
		return -1;
	}

	double Y_a=0,Cb_a=0,Cr_a=0;
	int r=0,g=0,b=0;
	int i=0,count=0,num=0;
	int rgbNumber;

	for(int m=0; m<=11; m++){

		printf("%s: ", matrix[m].name);

		// rgbTableをクリア(falseにする)
		for(i=0;i < RGB_TABLE_SIZE;i++){
			rgbTable[i]=false;
		}

		count=0;
		for (int Y = matrix[m].y_min; Y <= matrix[m].y_max; Y++){
			for (int Cb = matrix[m].uv_min; Cb <= matrix[m].uv_max; Cb++){
				for (int Cr = matrix[m].uv_min; Cr <= matrix[m].uv_max; Cr++) {

					// まずはデジタルYCbCr→アナログYCbCr変換
					// (Y_a:0.0～1.0、Cb_a,Cr_a:-0.5～0.5)
					// "アナログYCbCr"という表現はよろしくないけどスルーの方向で。
					// bitDepthをnとすると、H.264で定義されてる量子化式はだいたいこんな感じなので、
					// ここから逆算してY_a、Cb_a、Cr_aを求める。
					//　●TVレンジ
					//　　　Y = ( 1 << (n-8) ) * (219 * Y_a + 16)
					//　　　Cb = ( 1 << (n-8) ) * (224 * Cb_a + 128)
					//　　　Cr = ( 1 << (n-8) ) * (224 * Cr_a + 128)
					//　●フルレンジ
					//　　　Y = ( (1 << n) - 1) * Y_a
					//　　　Cb = ( (1 << n) - 1) * Cb_a + ( 1 << (n-1) )
					//　　　Cr = ( (1 << n) - 1) * Cr_a + ( 1 << (n-1) )

					// minとmaxを用意しとけばTVレンジでもフルレンジでも
					// １つの式でいけるので、そのように変形してみた。
					Y_a  = (Y -  matrix[m].y_min) / (double)(matrix[m].y_max - matrix[m].y_min);
					Cb_a = (Cb - ( 1 << (matrix[m].bitDepth - 1) ) ) / (double)(matrix[m].uv_max - matrix[m].uv_min);
					Cr_a = (Cr - ( 1 << (matrix[m].bitDepth - 1) ) ) / (double)(matrix[m].uv_max - matrix[m].uv_min);
					CLIP_Y(Y_a);
					CLIP_C(Cb_a);
					CLIP_C(Cr_a);

					// 続いてアナログYCbCr→デジタルRGB変換
					// 基本的には、まるも氏の
					//　　http://www.marumo.ne.jp/db2002_5.htm#15
					// の記事を参照。
					// ただしこの記事は2002年のものであり、
					// BT.709の係数はBT.709-1のものになっている。
					// H.264ではBT.709-2以降の係数が使われているので、
					// ここではBT.709-2の係数で計算している。
					// YUV→RGBのアナログ変換式のみ載せておく。
					// (R,G,B,Y:0.0～1.0、U,V:-0.5～0.5)
					//　●BT.601
					//　　　R = Y     　　　　　＋ 1.402 × V 
					//　　　G = Y － 0.344 × U － 0.714 × V 
					//　　　B = Y ＋ 1.772 × U 
					//　●BT.709-2
					//　　　R = Y     　　　　　 ＋ 1.5748 × V
					//　　　G = Y － 0.1873 × U － 0.4681 × V
					//　　　B = Y ＋ 1.8556 × U

					r = myRound(255.0 * (Y_a                         + matrix[m].r_cr * Cr_a));
					g = myRound(255.0 * (Y_a + matrix[m].g_cb * Cb_a + matrix[m].g_cr * Cr_a));
					b = myRound(255.0 * (Y_a + matrix[m].b_cb * Cb_a                        ));
					SATURATE(r);
					SATURATE(g);
					SATURATE(b);

					// 再現できた色はtrueにする
					rgbNumber = ((r << 16) | (g << 8) | b);
					rgbTable[rgbNumber] = true;
				}
			}
		}

		// trueになっている色の数をカウントして出力
		num = 0;
		for (i = 0; i < RGB_TABLE_SIZE; i++){
			if (rgbTable[i]){
				num++;
			}
		}
		printf("%d\n", num);
	}
	free(rgbTable);
	return 0;
}