グレイコードパターン投影

プロジェクタとカメラを利用したアクティブ3次元計測の手法はたくさんありますが，グレイコード(Gray code)パターンを投影する方法は実装が簡単で，精度をそこまで気にしない時にはよく使われる方法だと思います．私が所属する研究室でも，新しく所属する学生に対して，毎年のように話題にのぼります．

原理

グレイコードパターン投影の目的は，プロジェクタ（あるいはディスプレイ）とカメラのピクセルの対応関係を求めることです．つまり，カメラのあるピクセルに映っている点は，プロジェクタの何番目のピクセルなのかを求めることです．

最も簡単な方法は，1枚目の投影では x=0 のピクセルを白，その他を黒，2枚目の投影では x=1 のピクセルのみを白… のように投影して，撮影された画像で注目しているピクセルが何回目の投影で白く光るのかをチェックすることです．しかし，この方法では，XGA (1024×768) のプロジェクタのピクセルの対応関係を求めるのに，1024枚の投影が必要になってしまいます（y方向を考えるとさらに+768！）．

ここで，プロジェクタの座標を2進数（バイナリコード）に変換して考えます．x=0 は [0 0 0 0 0 0 0 0 0 0] で，x=1023 は [1 1 1 1 1 1 1 1 1 1] です．1枚目の投影では最上位ビットが 1 のピクセルは白， 0 のピクセルは黒になるようなパターンを投影し，2枚目では，上から2番目のビットが 1 のピクセルを白， 0 のピクセルは黒 … と最下位ビットまで繰り返します．このように撮影された画像群から，すべてにピクセルに対して，白いところは１，黒いところは０となるようにビット列を復元すると，そのピクセルのプロジェクタの x 座標が求まります．この方法にすることで，投影枚数は10枚で済みます．

バイナリコードに代わって，グレイコードを利用した方法が，グレイコードパターン投影です．なぜ，グレイコードを使うのかというと，観測ノイズに対して頑健にするためです．ビット列は画像から復元するため，何個目のビットに対しても同じような確率で間違った値を復元してしまう可能性があります（とくに白と黒の境目）．もし，上記のようなバイナリコードを利用していると，上位ビットを間違えたとき，求まる対応関係の誤差はかなり大きくなってしまいますが，ハミング距離が１のグレイコードであれば，どのビットを間違えても，誤差は常に１で収まります．

実装

グレイコードはハミング距離が１のビット列です．グレイコードはバイナリコードと，それを右にビットシフトしたものとの排他的論理和 (C言語で書くと，n^(n>>1)) で計算できます．Python では10進数 n のグレイコードはバイナリコードを通して次のように計算できます．

def int2bin(n):
    '''
    自然数 n を2進数（バイナリコード）に変換します．
    '''
    if n:
        bits = []
            while n:
            n,remainder = divmod(n, 2)
            bits.insert(0, remainder)
            return bits
    else:
         return [0]

def bin2gray(bits):
    '''
    バイナリコードをグレイコードに変換します．
    '''
    return bits[:1] + [i ^ ishift for i, ishift in zip(bits[:-1], bits[1:])]
 
bits = int2bin(n)
gray = bin2gray(bits)

これを用いて，投影パターンを生成します．カメラで撮影された画像から，投影されている色が白か黒かを簡単に判別するために，実装の工夫として，それぞれのパターンの反転画像も投影するようにしておきます．

import numpy as np
import cv2

def zero_pad_func(length):    
    '''
    ビット長をlength（上位ビットを0埋め）にするための関数を作ります．
    '''
    def zero_padding(code):
        '''
        ビット列 code の上位ビットを0で埋めて，ビット列をlengthにします．
        '''
        pad = length - len(code)
        if pad == 0:
            return code
        bits = [0] * pad
        bits.extend(code)
        return bits
    return zero_padding
    
def generate_projection_images(size=(768,1024)):
    h, w = size
    y = range(h)
    x = range(w)
    ii = 0
    
    # y 座標のグレイコードを生成します．
    bins = map(int2bin, y)
    gray = map(bin2gray, bins)    
    # 最大ビット長を求めます．
    maxlen = max(map(len, gray))
    # 上位ビットを0埋めして，ビット長を揃えます．
    gray2 = map(zero_pad_func(maxlen), gray)
    gray2 = np.array(gray2).T
    
    # 投影パターンを生成していきます．
    for i in range(maxlen):
        img = np.tile(gray2[i], w).reshape((w, h)).T
        cv2.imwrite('graycode/'+str(ii)+'.png', img * 255)
        # 反転パターンも作ります．
        img = 1 - img
        cv2.imwrite('graycode/'+str(ii)+'n.png', img * 255)
        ii += 1

    # x 座標も同様に行います．
    bins = map(int2bin, x)
    gray = map(bin2gray, bins)    
    maxlen = max(map(len, gray))
    gray2 = map(zero_pad_func(maxlen), gray)
    gray2 = np.array(gray2).T    
    for i in range(maxlen):
        img = np.tile(gray2[i], h).reshape((h, w))
        cv2.imwrite('graycode/'+str(ii)+'.png', img * 255)
        img = 1 - img
        cv2.imwrite('graycode/'+str(ii)+'n.png', img * 255)
        ii += 1

複数の投影画像が生成されるので，それをプロジェクタから投影しながら，カメラで撮影します．

つぎに，撮影された画像群からビット列（グレイコード）を復元します．

def decode(x_bits=10, y_bits=10, shape=(768,1024)):
    '''
    撮影画像から，プロジェクタの座標を復元します．
    '''
    
    ii = 0
    gray = []
    for i in range(y_bits):
        # i 番目のビットの投影画像を読み込む．
        posi = cv2.imread('images/'+str(ii)+'.png', 0).flatten()
        nega = cv2.imread('images/'+str(ii)+'n.png', 0).flatten()
        # 各画素のビットを復元．（正負を比較）
        bits = [1 if p > n else 0 for p, n in zip(posi, nega)]
        gray.append(bits)
        ii += 1
    # 復元されたビット列（グレイコード）を自然数に変換
    gray = np.array(gray).T
    bins = map(gray2bin, gray)
    y = map(bin2int, bins)
    # 各画素におけるプロジェクタの y 座標
    y = np.array(y).reshape(shape)    
    # 画像で対応関係を見てみる
    cv2.imwrite('y.png', y * 255 / shape[0])
    
    # x 座標も同様に行う．
    gray = []
    for i in range(9, 18):
        posi = cv2.imread('images/'+str(ii)+'.png', 0).flatten()
        nega = cv2.imread('images/'+str(ii)+'n.png', 0).flatten()
        bits = [1 if p > n else 0 for p, n in zip(posi, nega)]
        gray.append(bits)
        ii += 1
    gray = np.array(gray).T
    bins = map(gray2bin, gray)
    x = map(bin2int, bins)
    # 各画素におけるプロジェクタの x 座標
    x = np.array(x).reshape(shape)
    cv2.imwrite('x.png', x * 255 / shape[1])

    return x, y

グレイコードから10進数の変換は次のように行います．

def bin2int(bits):
    '''
    2進数のビット列を自然数に変換します．
    '''
    i = 0
    for bit in bits:
        i = i * 2 + bit
    return i

def gray2bin(bits):
    '''
    グレイコードをバイナリコードに変換します．
    '''
    b = [bits[0]]
    for nextb in bits[1:]: b.append(b[-1] ^ nextb)
    return b

これにより，カメラのすべてのピクセルで，プロジェクタの座標が求まります．

利点と欠点

グレイコードパターン投影は，実装が簡単で，投影枚数も log(n) ，かつノイズに頑健です．しかしながら，求まる対応関係は整数値であるため，精度が荒いという問題があります．サブピクセルの精度で対応関係を求めたければ，位相シフト法（フェーズシフト・Phase shifting）等を検討しましょう． また，プロジェクタのピントが合っている場所しか対応点は取れません．とりあえず，簡単に対応関係を求めたい，そんなに精度は必要としていない，という場合には大いに役に立つ方法だと思います．