hirax.net::Keywords::「OpenCV」のブログ

■「手持ちスマホ撮影動画からの超巨大開口レンズ撮影」に挑戦してみよう!?　

　かつて、スマホに搭載されているカメラのレンズはとても小さく、綺麗なボケとは無縁の存在でした。しかし、今や最新のスマホには特殊処理によるボケ生成機能などが備えられています。カメラレンズの光学開口径が小さくとも、たとえば２眼カメラなどを備えて距離情報を取得して、距離情報などからボケを人工的に合成するといった仕組みです。

　そんな最新スマホを持たずとも、大レンズのボケ味を手に入れるために、「手持ちスマホ撮影動画からの超巨大開口レンズ撮影」に挑戦してみました。スマホ動画から巨大開口レンズ撮影の手順はとても簡単、まずは目の前の風景にスマホを向けて・なるべくスマホが平面上を動くように意識しながら（スマホを動かしつつ）動画撮影します。そして、動画の各コマから画像ファイル群を生成し、それぞれの画像が撮影された位置や方向にもとづいて撮影された光情報を加算合成する、というものになります。

　細かい手順は、スマホ撮影動画（の展開画像をもとに） Bundler: Structure from Motion (SfM) から出力された刻々のカメラ位置・方向や特徴点情報ファイル(bundler.out)を読み込み、それらのカメラ情報にもとづいて、刻々の撮影画像をレンズ開口面に沿った（同じ方向を向く平行カメラが存在していた場合の）光線画像を位置・角度ズレを踏まえて重ねることで、任意のピント位置に焦点を合わせた超巨大開口レンズ撮影画像を生成する…というものです。

　試しに、iPhoneを約1.5m×1.0mの範囲で動かしつつ動画撮影し、つまり、レンズ直径約1.5ｍに相当する範囲で動かしつつ動画撮影し、その画像群から開口合成により超巨大カメラの撮影画像を作り出してみた結果が右上の画像です。

　右上画像を眺めてみても、良好なボケ味どころか、全くピントが合っていない画像にしか見えません。直径が1mを超える開口を持つカメラレンズとなると焦点深度もとても浅くなるのでピントがなかなか合わない…というわけでなく、手持ち撮影動画からのカメラ位置・方向精度が低いせいか、単一カメラに平行合成した後のズレが大きいようです。

　ちなみに、試しに各画像を（撮影方向による傾きを補正しつつ）位置毎に並べてみると、下の画像のようになります。動画撮影からのカメラ位置推定精度が果たして不十分なのかどうか、次は撮影カメラ位置を精度良く知る事ができる撮影治具でも作り、また再挑戦してみたいと思います。

　上記処理のコード手順、Python/OpenCVで書いたコード処理手順は、bundler.outからカメラ位置・方向・焦点距離や歪みパラメータを読み込み、cv2.initUndistortRectifyMapにカメラ情報を渡して、各撮影画像の向き補正用のホモグラフィーマップを作成してremapで変換した後に、各撮影画像を加算合成するという手順です。

　Bundlerの出力ファイルを読み込んでライトフィールド合成を行うOpenCV/Pythonコード、まだまだ間違い含まれているような気もしますが、とりあえずここに貼り付けておくことにします。

import cv2
import numpy as np
from matplotlib import pyplot as plt
from PIL import Image
import math
%matplotlib inline

class camera:
    def __init__(self):
        self.f = 1.0
        self.k1 = 0.0
        self.k2 = 0.0
        self.R = [[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]
        self.T = [0.0,0.0,0.0]

def readBundlerOut( filePath ):
    f = open(filePath, 'r')
    list = f.readlines()
    f.close()

    numberOfCameras = int((list[1].split())[0]) 
    cameras = []
    for i in range(numberOfCameras):
        aCamera = camera()
        fk1k2 = [float(j) for j in list[5*i+2].split()]
        aCamera.f  = fk1k2[0]
        aCamera.k1 = fk1k2[1]
        aCamera.k2 = fk1k2[2]
        rot = []
        for j in range(1,4):
            rot.append( [float(k) for k in list[ 5*i+2+j ].split()] )
        aCamera.R = rot
        aCamera.T = [float(k) for k in list[ 5*i+2+4 ].split()]
        cameras.append(aCamera)
    return cameras

def readImageList( listPath, imageDirPath ):
    f = open(listPath, 'r')
    list = f.readlines()
    list = [ i.rstrip() for i in list ]
    f.close()
    list = [imageDirPath+fileName for fileName in list]
    return list

class lightField:
    
    def __init__(self):
        self.w = 2000
        self.h = 2000
    
    def loadImageListAndMakeLightField( self, 
                    imagePathList, cameraList, workList, scaleA ):
        self.cimg  = np.zeros((self.h, self.w,3), dtype=np.uint8)
        sum = 1.0
        for i in workList:
            img = cv2.imread( imagePathList[i], cv2.IMREAD_COLOR )
            h, w = img.shape[:2]
            imageHeight = img.shape[0]
            imageWidth =  img.shape[1]
            focalLength = cameraList[i].f
            principalPointX = 0.500000
            principalPointY = 0.500000
            distCoef = np.array([ 0.0, 0.0, 0.0, 0.0, 0.0 ])
            cameraMatrix = np.array([ 
                [focalLength,   
                 0.0,          
                 imageWidth  * principalPointX], 
                [0.0,           
                 focalLength,  
                 imageHeight * principalPointY], 
                [0.0,           0.0,          1.0] 
                           ])
            newCameraMatrix, roi = cv2.getOptimalNewCameraMatrix(
                cameraMatrix,
                distCoef,
                (img.shape[1], img.shape[0]),1,
                (img.shape[1], img.shape[0]) )
            rotMatrix = np.array( cameraList[i].R )
            map = cv2.initUndistortRectifyMap( 
                newCameraMatrix,
                distCoef,
                rotMatrix,
                newCameraMatrix, 
                (img.shape[1], img.shape[0]),
                cv2.CV_32FC1)
            
            undistortedAndRotatedImg = cv2.remap( img, 
                        map[0], map[1],
                        cv2.INTER_LINEAR )
            
            scale = 1.0
            pt3 = np.array(cameraList[i].T) - np.array(cameraList[0].T)
            x =  ( self.w/2.0 - pt3[0] * scale * scaleA ) 
            y =  ( self.h/2.0 - pt3[1] * scale * scaleA ) 
            pts1 = np.float32( [[0,      0],
                                [w, 0],
                                [w, h],
                                [0, h]])
            pts2 = np.float32(  [[x,           y],
                                 [x + w*scale, y],
                                 [x + w*scale, y + h*scale],
                                 [x,           y + h*scale]] )
            M = cv2.getPerspectiveTransform( pts1, pts2 )
            img2 = cv2.warpPerspective( 
                undistortedAndRotatedImg, M, (self.w, self.h) )
            sum = sum + 1.0
            self.cimg = cv2.addWeighted(
                 self.cimg, (sum-1) / sum, 
                 img2,      (1.0) / sum, 0)
    
    def showImage(self):
        plt.figure( figsize=(14,14) )
        plt.imshow( np.array(self.cimg2) ) 
        plt.autoscale( False )

■Ricoh Theta S で24時間のHDR（多段露出）撮影をして遊んでみる!?　

　Ricoh Theta Sを制御するPython 3.xコードを書いた （「Ricoh Theta S で2017年旧暦1月満月の夜を撮ってみる」）ので、Theta Sを載せた三脚を高く立て、24時間を少し超えるくらいの時間にわたり、11段階の露出時間での撮影を交互に続けてみました。ちなみに、撮影に使ったPython コードは下記になります。

import osc
import theta
import time
thetas = RicohThetaS()
thetas.setCaptureMode( 'image' )
thetas.setOption("exposureProgram", 1 )
thetas.setOption("ISO", "100")
thetas.setOption("whiteBalance", "daylight")
shutterSpeeds1 = [0.00015625, 0.0005,
    0.002,0.008,0.03333333,0.125,0.5,
    2,8,30,60]
for i in range(400):
    for shutterSpeed in shutterSpeeds1:
      thetas.setOption("shutterSpeed", shutterSpeed)  
      time.sleep(1)
      response = thetas.takePicture()
      thetas.waitForProcessing(response['id'])
      time.sleep(1)
      thetas.getLatestImage()
      time.sleep(1)
thetas.closeSession()

　こうして、 1/6400秒～60秒までの11段階にわたる露出時間で撮影した画像群を、下記のPythonコードでHDR(ハイダイナミックレンジ)合成します。

import cv2
import numpy as np
iniID = 11006 # fileName
for i in range(400):
    files=["hoge"+"R00"+str(iniID+j+i*11)+".JPG" for j in range(11)]
    imgs = [cv2.imread(file) for file in files]
    mertens = cv2.createMergeMertens()
    resImg = mertens.process(imgs)
    resImg = np.clip( resImg*255, 0, 255 ).astype('uint8') 
    cv2.imwrite("hoge"+
                "out_"+str(1000+i)+".jpg", resImg)

　そして、11段階の多段露出画像を合成した連番画像ファイルを動画ファイルとして結合することで、たとえば、下に貼り付けたようなLittle Planet風のタイムラプス（微速度）動画として眺めることができます。撮影を行ったのは、都会に近い薄曇りの空の下ですが…それでも恒星が浮かぶ宇宙空間を回転する地球が進んでいるさまが見えるような気がしてきます。

■現実のリアルタイム風景と、そこから眺めた好きな時間・瞬間を合成してVR的に眺めてみよう!?　

　先日、周囲全方向を撮影することができるRicoh Theta Sを使い、24時間にわたり11段階の多段露出撮影画像からリトル・プラネット（とても小さな星を上空から撮影したような）風動画を作ってみました（Ricoh Theta S で24時間のHDR（多段露出）撮影をして遊んでみる!?）。けれど、そんな風に景色を離れたところから第三者的・客観的に眺めるのではなくて、その場に立って一人称的に景色を眺めたくなりました。そこで、合成した高画質4K動画をYoutubeに（３６０度映像だというタグを付けて）アップロードしてみました。

　こうすることで、スマホからYoutubeアプリケーションでアクセスし、VR用レンズセットを取り付けて眺めれば、あたかもこの場所に立ち・数百倍近く速く時間が流れる中で、夜空や太陽の影を眺めることができます。

　そして一番面白いのは、２４時間の撮影を行った場所に立って、目の前にあるリアルタイムな実際の景色を眺めつつ、「東の空に昇る太陽や夜星や夕焼けや」…２４時間の中にある美しい瞬間や一番好きな時間を、スマホに覗いて自由自在にVR合成しながら眺めてみると…とても不思議で面白い気持ちになってきます。

　今日は、２４時間にわたる、周囲全ての４πステラジアンの全天周風景を、11段階多段露出のハイダイナミックレンジ撮影をして、その映像を4K映像に変換して眺めてみました。

■草間彌生デザインの「水玉模様のモジモジ君ウェア」で人体の表面形状を推定してみよう！　

　数ヶ月前に買ったつもりのZOZOスーツ、服に取り付けられたセンサ群とスマホの間でさまざまなデータが交換され、自分の体を知ることができるという「面白さ」に惹かれて…はや数ヶ月、「大幅な性能向上」したものが届くという連絡内容を見ると、、そこには「コレジャナイ感が、超大盛りラーメン店のようにテンコ盛りされた、草間彌生デザインの「水玉模様のモジモジ君ウェア」でした。マーカー付けた衣服を使って（スマホによる）画像計測による採寸を行うデザインに変わっていた…というわけです。

　ネットには、新バージョンのZOZOスーツ試着写真もtwitterにはアップロードされ始めています。そこで、そんなマーカ模様が付けられた服を着用した画像からの人体形状推定をしてみることにしました。「ネットにアップされた画像から」というわけで、多視点撮影からの3次元推定を行うZOZOスーツ正式バージョンとは違い、一枚画像からの体表面形状の推定です。

　書いてみたのは、下に貼り付けたような、十数行ばかりのPython/Jupyterコードです。OpenCV でテクスチャ検出をして、その（円形マーカーの）大きさや形状の歪みを使って、模様部分の体の表面形状情報を推定する…というわけです。Pythonコードを実行した結果は、たとえば下の4 図のようになります。向かって左から、元画像・高さ（凹凸）画像・向き（180度の反対側を区別できない）画像・法線画像です。この処理例は、水玉模様の水着画像で処理を行ったものですが、新型ZOZOスーツ画像でも同様のことが行えます。

　この画像を見れば、水着の模様を介して胸やお腹の3次元形状が見えてくることがわかります。全身に模様が付けられたZOZOスーツを着用した自撮り写真がアップされたら…色んな体形状を可視化できそうです。

import numpy as np
import cv2
from matplotlib import pyplot as plt
from math import sin, cos

img = cv2.imread('zozo_polka202.jpg',2)
cimg = cv2.cvtColor(
    cv2.imread('zozo_polka202.jpg'),
    cv2.COLOR_RGB2BGR)
img = cv2.adaptiveThreshold(img, 128, 
    cv2.ADAPTIVE_THRESH_MEAN_C, 
    cv2.THRESH_BINARY_INV, 15, 5)
img, contours, hierarchy = cv2.findContours(
  img, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
if len(contours) > 0:    
    for i, contour in enumerate(contours):
        if (contour.shape)[0]>4:
            ellipse = cv2.fitEllipse(contour)
            aMax = max(ellipse[1][0],ellipse[1][1])
            aDiff = abs(ellipse[1][0]-ellipse[1][1])
            if aMax < 20 and aMax > 11 and aDiff < 12:
                xc = sin(ellipse[2]*np.pi/180)*255
                yc = cos(ellipse[2]*np.pi/180)*255
                cx=int(ellipse[0][0])
                cy=int(ellipse[0][1])
                cv2.ellipse(cimg,ellipse,
                            (xc,yc,255-xc),-1)
                #cv2.circle(cimg, (cx, cy), int((aMax-11)*2.0), 
                            (xc,yc,255-xc), -1, -1)
                #cv2.ellipse(cimg,ellipse,(int((aMax-11)*50.0),
                           int((aMax-11)*50.0),int((aMax-11)*50.0)),-1)
                #cv2.ellipse(cimg,ellipse,(int((aMax-11)*50.0),
                             yc+xc,int((aMax-11)*50.0)),-1)

plt.imshow(np.array(cimg))