はじめに
DeNA 20 新卒 Advent Calendar 2019の記事です．
来年4月からDeNAでWebエンジニアになるのですが，普段は大学院で光についての研究をしています．
雑に研究紹介をするとレーザーの形状を変化させてぶっ放したり，光を使って物体計測をしたりしています（ちなみに，僕のQiitaアカウント画像をSLMというものに表示してレーザーを当てると，レーザーが100本に分岐します）．
本記事ですが，せっかく光の研究をしているので今回はPythonを使って簡単に実装できる「単一画素カメラ」というものを紹介したいと思います（行列計算ができれば何でも実装できるのでMATLABやOctaveでも楽々実装できます)
単一画素カメラ(シングルピクセルカメラ)
普通のカメラ
単一画素カメラの原理の話をする前に，まずは通常のカメラの仕組みからお話しします．
めっちゃザックリ言うと，通常のカメラはCCDと呼ばれるディテクター(光の強さを検出できるもの)がいっぱい並んでいるものを使って撮影を行っています．
例えば，100×100pixelの画像が撮影できるカメラなら，100×100個のディテクターが並んでいるわけです．
そして，対象から反射してくる光を各ディテクターで検出して，検出した光の強さを画素値にして画像にしているのです．
単一画素カメラの原理
では単一画素カメラはどうかというと，このディテクターがたったの一個で画像が撮れてしいます．
もうちょっと言うと
「一つのディテクターを使って複数回計測した後，行列計算で画像を再構成できるカメラ（手法）」
です（一つのディテクターで撮影できるから単一画素という名前なんです）．
これだけ聞いても実際のイメージが湧かない方が大半だと思います．
大丈夫です．単一画素カメラの仕組み自体はマジでめっちゃ簡単です．ぶっちゃけ僕が記事のネタに選んだ理由も解説が簡単そうだったから．．．
単一画素カメラで実際に撮影をしている様子は下の図のようになります．
ここで，Sceneが撮影対象です．レンズは光をディテクターに集めるために使っています．
上図のように，物体から来た光はマスクを通った後，すべて一つのディテクターに集められ検出されます．
ここで急に登場したマスクパターン．これが重要です．こいつのおかげで単一画素カメラくんは一つのディテクターだけで撮影ができるのです．
「物体から来た光がマスクパターンを通った後の光強度を検出する」(マスクと撮影対象の順番はどちらが先でも構いません)．
自分であらかじめ用意しておいたマスクパターンを次々と変更し，上記の作業を繰り返すことで単一画素カメラの撮影は行われます．
ここで，i回目の光強度検出の様子は
マスクパターン(符号化行列)をWi, 撮影対象をX, 検出値をYiとし，Wiが以下のように(1,N)に変形でき，Xは(N,1)に変形できるとすると
のようになります．ここで、Nはマスクの画素数であり，このNはそのまま撮影対象の画素数になります．
上図のように，i回目の計測結果YiはWiとXの内積
$$ Wi・x = Yi $$
と表現できます．
この計測をM回行った時の結果は以下の式のように，行列の積$(M,N)×(N,1)=(M,1)$として表すことができます(通常 N = Mの時，対象を完全に再構成できます)．
W,Yの一行はそれぞれ，各計測回で使用したマスク，計測結果です．
今，X,W,Yは
- X : 観測対象・・・今撮りたいもの．どんなものかは分からない
- W : 符号化行列・・・自分で用意したものだから分かっている
- Y : 検出値・・・検出した値だから分かっている
という状態です．
そしてYは
$$ WX=Y $$
とWとXの行列の掛け算で表すことができます．
ここから，Xの様子を知るにはどうすればよいでしょうか？？
そう，簡単ですね．両辺にWの逆行列$W^{-1}$を掛けてやればいいだけです．そうすると
$$ X = W^{-1}Y $$
となりXを得ることができました．
これが単一画素カメラの撮影原理です．
メリット・デメリット
わざわざこんな作業をして単一画素カメラを使うメリットの一つとして，めっちゃ弱い光でも撮影ができるという点があります．
なぜめっちゃ弱い光でも撮影ができるのかというと，単一画素カメラは一つのディテクターに光を集めて検出するからです．
たとえば $10$ の光が撮影対象からきているとして，10×10pixelのカメラでは各画素ごとにディテクターを分けているため，各ディテクターで検出できる光量は$10/(10×10)$となります．
しかし，単一画素カメラは光を一つのディテクターに集めて検出するため光量は10のままです（実際はマスクを通るので4とか5とかに落ちますが）．
そのため，通常のカメラでは撮影できないような光量でも撮影を行うことができるのです．
デメリットは一枚の画像を取るのに複数回の計測が必要なことですね．なので「ある一瞬」を捉える用途であったり，動画の撮影には向いていません．
Pythonでのシミュレーション
それではいよいよPyhonで単一画素カメラの流れを実装してみます．
まぁ実装と言っても，原理で述べたように，行列計算をしているだけなんで，その計算をPythonにさせるだけです．
まずは使用するライブラリをインポートします．行列を扱うためにnumpy，画像を扱うためにopenCVを使います．
import numpy as np
import cv2
その後，撮影対象とする画像を読み込みます．また画像の縦横サイズを取得しておき，そこからマスクの画素数Nを計算しておきます
X = cv2.imread("img/target.jpg",0)
h,w = X.shape
N = h*w
また，画像を表示する処理をshow関数として用意しておきます.
def show(img):
img = (img/np.max(img))*255
img = img.astype("uint8")
cv2.imshow("img", img)
cv2.waitKey(0)
読み込んだ画像を表示してみましょう
show(X)
今回は16×16pixelの十字の画像を用います．
次に，マスクを用意します．今回は画素数Nのランダムな白黒マスクでM回計測すると想定します（M=Nとします）．
M = N
W = np.random.rand(M, N)
W = np.where(W>0.5, 1, 0)
show(W)
print(W.shape)
画像が表示され, printの結果256,256が表示されたでしょうか?
これは256画素のマスクで256回計測するという意味です．
1回目の計測で使われるマスクも確認しておきましょう
show(W[0].reshape(h,w))
では次に計測結果であるYを取得していきましょう.
Wから順次マスクを取り出して，撮影対象との内積をとっていき，それをYに追加していきます
Y = np.dot(W,X.reshape(N,1))で計算すれば一発ですが，今回は実際の計測と同様一回一回値をとっていきます．
この時，Xのshapeを(h,w)から(N,1)にして置くのを忘れずに．
Y = []
for mask in W[:M]:
# i回目の計測
yi = np.dot(mask.reshape(1,N), X.reshape(N,1))
Y.append(yi)
"""
同じ意味
Y = np.dot(W,X.reshape(N,1))
"""
Y = np.array(Y).reshape(M,1)
print(Y.shape)
(256, 1)と表示されたでしょうか？ これは256回分の計測結果があるよという意味です．これで
$$WX = Y$$
ができ，無事計測が終わった状態です．
次にWとYを使ってXを再構成してみましょう．
$$ X = W^{-1}Y$$
Xを得るためには上式のようにWの逆行列である$W^{-1}$が必要でした．
そこでnumpyのpinvを使って逆行列を計算します
InvW = np.linalg.pinv(W[:M])
print(InvW.shape)
ではこのInvWとYを使ってXを取得，表示します．Xは(N,1)の状態で再構成されるため，画像として表示する場合は(h,w)を使ってreshapeする必要があることに注意してください.
rec = np.dot(InvW, Y)
print(rec.shape)
show(rec.reshape(h,w))
無事に元の画像が取得できていれば成功です．
以上が単一画素カメラの一連の流れになります．めちゃくちゃ簡単なコードで実装できることが分かっていただけたと思います．
実際の実験でも上記と同様に，Wに使用したマスク，Yにディテクターの検出値を入れてやれば撮影対象Xを再構成できます．
ちなみに，ここでもし計測回数Mを減らすとうまく画像を取得できないと思います．
これは通常計測回数M=Nである必要があるためです．
しかし圧縮センシングという技術を使うとこの計測回数Mを大幅に減らすことができます．通常はこの技術を組み合わせることによって計測回数を減らします．
圧縮センシングも解説すると長くなるので今回は触れませんが，興味のある方はぜひ実装してみてください．
こちらの解説なんかは分かりやすいです．
http://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/Presentation/lecturenote20160822.pdf
以上で僕の担当記事を終わります．
引き続きDeNA 20 新卒 Advent Calendar 2019をよろしくお願いします．
それでは :D