ディストーション係数K3

「instantNeRFで遊ぶ Advent Calendar 2022」の24日目です。

今日はinstatnNeRFでディストーション係数K3を使用する方法を検討しました。12日目の「RealityCaptureのregistrationデータを流用する検討」にあった方針の案2です。

業務外の個人的な検討活動です。

すみません　今回は狙い通りの効果が得られませんでした。

警告
環境や入力内容によっては現在使用できている状況が壊れることがあります。
あなたの責任と判断で試してください。私は責任を取りません。

テストデータ

sketchfabのNefertiti statueを使用させていただきました。ありがとうございます。

instantNeRF / RealityCapture ディストーション係数の比較

instantNeRFとRealityCaptureが使うレンズのディストーション係数の組み合わせが合ってないようです。

アプリケーション	ディストーション係数
instantNeRF COLMAP	K1, K2, p1, p2; OPENCV
RealityCapture	K1, K2, K3, p1, p2; Brown3+tangential K1, K2, K3, K4, p1, p2; Brown4+tangential

実際の被写体とは別の被写体（特徴点が広く分布している例えば地面など）を撮りCOLMAPでディストーション係数を解析し、RealityCaptureのPrior lens distortionでFixにする（そのときのK3はゼロにする）方法もあります。

今回は試しにinstantNeRFをK3まで読み込めるように改造しRealityCaptureだけでカメラパラメータの解析を完結する方法を検討します。

Blender ディストーションした画像の生成

compositingの中でcompositorとMovie clip editorを使いディストーションが有効になった状態でRenderできました。

Movie clip editorでディストーション係数を設定しました。

ここでK3のような高次の係数は画像周辺の影響が強いので特徴点が出やすくなる背景を使用しました。

RealityCaptureで検算

Blenderでシミュレーションした画像を入力したところBlenderで設定したディストーション係数と近い値が表示されました。

RealityCaptureのカメラパラメータのcsvファイルを変換

RealityCaptureからexportしたカメラの内部パラメータと外部パラメータからそれぞれcameras.txtとimages.txtを変換しました。
長いので折りたたんでいます。▶を押すと展開します。

RC2COLMAPtxt.py

# csvファイル名と画像サイズは修正してください。
import csv
import math
import os

import numpy as np
from scipy.spatial.transform import Rotation as Rot

filePath = "./"
convertOrder = "ZXY"
eulars = []
trans = []
imageFiles = []
rt = []
camParams = []
imgSize = [1920, 1080]
bottom = np.array([0.0, 0.0, 0.0, 1.0]).reshape([1, 4])
dummyPoint = "99, 99, -1\n"

def importRCFile(imageFilesOut, eularsOut, transOut, camParamOut, filePath, fileName):
    with open(os.path.join(filePath, fileName), "r") as f:
        for line in f:
            eul = []
            tra = []
            camP = []
            if (line[0] == "#") or (line[0] == "\n"):
                continue
            lin = line.split(',')
            # print(lin)
            imageFilesOut.append(lin.pop(0))
            # print(lin)
            i = 0
            for i in range(3):
                tra.append((lin.pop(0)).strip())
            i = 0
            for i in range(3):
                eu = np.float_(lin.pop(0).strip())
                eul.append(eu)
            eul2 = [-eul[0], eul[1], eul[2]]
            for i in range(9):
                cams = np.float_(lin.pop(0).strip())
                camP.append(cams)
            # eul.reverse()
            transOut.append(tra)
            eularsOut.append(eul2)
            camParamOut.append(camP)

def exportCamerasTxt(camParams):
    with open(os.path.join(filePath, "cameras.txt"), "w") as f:
        f.write("# Camera list with one line of data per camera:\n")
        f.write("#   CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]\n")
        f.write("# Number of cameras: 1\n")
        f.write("1 OPENCV {} {} ".format(imgSize[0], imgSize[1]))

        camP = camParams[0]

        fx = camP[0] / 36.0 * imgSize[0]
        f.write("{} {} ".format(fx, fx))

        px = camP[1] / 36.0 * imgSize[0] + imgSize[0]/2
        py = camP[2] / 36.0 * imgSize[0] + imgSize[1]/2
        f.write("{} {} ".format(px, py))

        # K1, K2
        f.write("{} {} ".format(camP[3], camP[4]))

        # p1, p2
        f.write("{} {} ".format(camP[7], camP[8]))

        # K3
        f.write("{} \n".format(camP[5]))

def exportImages(imageFiles, eulars, trans):
    with open(os.path.join(filePath, "images.txt"), "w") as f:
        f.write("# Image list with two lines of data per image:\n")
        f.write("#   IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME\n")
        f.write("#   POINTS2D[] as (X, Y, POINT3D_ID)\n")
        f.write("# Number of images: {}, mean observations per image: 999\n".format(len(imageFiles)))

        rot2 = Rot.from_euler("X", 180, degrees=True)

        lines = len(imageFiles)
        i = 0
        for (img, eular, tran) in zip(imageFiles, eulars, trans):    
            f.write("{} ".format(i))

            rot1 = Rot.from_euler(convertOrder, eular, degrees=True)
            rot1 = rot1 * rot2
            rot1 = rot1.inv()   # r^t = r^-1
            qua1 = rot1.as_quat()

            f.write("{} {} {} {} ".format(qua1[3], qua1[0], qua1[1], qua1[2]))

            mat1a = rot1.as_matrix()

            tranSrc = list(np.float_(tran))
            tranSrc = [n*-1 for n in tranSrc]

            tranDst = np.dot(mat1a, tranSrc)
            f.write("{} {} {} 1 ".format(tranDst[0], tranDst[1], tranDst[2]))

            f.write("{}\n".format(img))

            f.write(dummyPoint)

            i += 1

if __name__ == "__main__":

    TEXT_FOLDER = "."
    fileName = "221208.csv"
    sharpness = 100
    imageID = 1
    downScale = 1/1

    importRCFile(imageFiles, eulars, trans, camParams, TEXT_FOLDER, fileName)

    exportCamerasTxt(camParams)

    exportImages(imageFiles, eulars, trans)

回転を変換

カメラの回転角度の値は、RealityCaptureの場合は真下を向いた状態からの変化、COLMAPの場合は真上の向いた状態からの変化を示していました。つまりRealityCaptureのカメラの角度をそのままクォータニオンを計算すると正反対の方向になるようです。
したがって下記のように回転を後乗せしました。

      # カメラの初期設定の違いを後乗せする分
      rot2 = Rot.from_euler("X", 180, degrees=True)
...
            # RealityCaptureのカメラ向きrot1をCOLMAP用に変換する
            rot1 = Rot.from_euler(convertOrder, eular, degrees=True)
            rot1 = rot1 * rot2

掛け算の順番が違うと欲しい結果が変わるのをBlenderで確認しました。
上の段のようにrot2 x rot1とするとrot2の回転の前にrot1の回転が反映されました。
今回欲しかったのは下の段の方です。

使用したデータ

長いので折りたたんでいます。▶を押すと展開します。

RealityCaptureがexportしたカメラパラメータ

#name,x,y,alt,heading,pitch,roll,f,px,py,k1,k2,k3,k4,t1,t2
cam1.png,-2.056299069,-63.85277789,19.10340477,-177.1107172,89.9971589,-177.2701738,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0
cam2.png,-2.055653469,-63.85258886,6.75191013,-179.8419707,75.00129077,179.9992156,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0
cam3.png,-8.231942591,-63.83614473,6.751286222,-171.8410939,75.00086291,-179.9983833,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0
cam4.png,-8.231477635,-63.83610116,19.10307359,-171.8431785,89.99685868,180,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0
cam5.png,-8.230959122,-63.83774112,31.45513982,8.157316563,75.00681903,1.10E-03,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0
cam6.png,-2.056038613,-63.85450434,31.45553271,0.159211651,75.00650065,0,50,0,0,-9.89E-02,0.496456017,1.001879389,0,0,0

生成されたcameras.txt

# Camera list with one line of data per camera:
#   CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]
# Number of cameras: 1
1 OPENCV 1920 1080 2666.6666666666665 2666.6666666666665 960.0 540.0 -0.0989 0.496456017 0.0 0.0 1.001879389 

生成されたimages.txt

# Image list with two lines of data per image:
#   IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME
#   POINTS2D[] as (X, Y, POINT3D_ID)
# Number of images: 6, mean observations per image: 999
0 0.7071236067921164 0.7070885859532999 -0.0009830945298747243 0.000984813188646648 1.8785412097644945 19.106572100173107 63.857307095127375 1 cam1.png
99, 99, -1
1 -0.7933457229281179 -0.6087697789861539 0.0008449647200686603 -0.0010982437949960673 1.8787114967970358 23.04816480481967 59.93512866386523 1 cam2.png
99, 99, -1
2 0.7913393057973016 0.6072258011075182 -0.04329639367758816 0.056430061705555194 -0.9092362789428219 23.1774808335617 60.41833984860332 1 cam3.png
99, 99, -1
3 -0.7053354758729568 -0.7052968059444011 0.05028918707048552 -0.050291944319993453 -0.909053600881123 19.106602088533435 64.35715837520632 1 cam4.png
99, 99, -1
4 0.6072672777968686 0.7913083305078705 -0.05641957176419029 0.04329447202557145 -0.9090058656610509 13.734171030819649 70.30631518505511 1 cam5.png
99, 99, -1
5 0.6088058464540544 0.7933180389587617 -0.001102223908086942 0.0008458654995738078 1.87859408094634 13.863456825088504 69.82371378330606 1 cam6.png
99, 99, -1

このcameras.txtとimages.txtをinstant-ngpにあるcolmap_txtへコピーし、使用した画像ファイルをdataフォルダの所定のフォルダへ置き、colmap2nerf.pyを実行します。ここではCOLMAPを起動しないので、下記のように入力しました。

(ngp) D:\git\instant-ngp>python scripts/colmap2nerf.py --aabb_scale 4 --images data/221216 --out 221216K3.json

生成された221216K3.json (transform.json)

{
  "camera_angle_x": 0.6911111611634243,
  "camera_angle_y": 0.39959649248062956,
  "fl_x": 2666.6666666666665,
  "fl_y": 2666.6666666666665,
  "k1": -0.0989,
  "k2": 0.496456017,
  "k3": 1.001879389, 
  "p1": 0.000001,
  "p2": 0.000001,
  "cx": 960.0,
  "cy": 540.0,
  "w": 1920.0,
  "h": 1080.0,
  "black_transparent":true,
  "aabb_scale": 4,
  "frames": [
    {
      "file_path": "./data\\221216/cam1.png",
      "sharpness": 241.32418836442383,
      "transform_matrix": [
        [
          -0.17188851282110787,
          -4.14955160892982e-06,
          0.9066038351385488,
          3.8238249300878944
        ],
        [
          0.9845229052962888,
          1.6522707932588106e-06,
          0.1716716601672445,
          0.7134249882238327
        ],
        [
          -2.2103154515103417e-06,
          0.9220806280830399,
          3.801322236269476e-06,
          1.7370481813517463e-05
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    },
    {
      "file_path": "./data\\221216/cam2.png",
      "sharpness": 221.6554947487402,
      "transform_matrix": [
        [
          -0.17189910816637333,
          0.2345777272492189,
          0.8757283604150276,
          3.8238551169875983
        ],
        [
          0.9845208988933357,
          0.04443013516687471,
          0.16583447377329139,
          0.7134857486412896
        ],
        [
          -7.655465396356803e-06,
          0.8906796707270588,
          -0.23858417550265393,
          -1.0513006329621049
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    },
    {
      "file_path": "./data\\221216/cam3.png",
      "sharpness": 239.74123384368445,
      "transform_matrix": [
        [
          -0.04407481613502864,
          0.23849301962508987,
          0.8903115549726209,
          3.9219164616568802
        ],
        [
          0.9988324567829607,
          0.00853396538394311,
          0.031896781531901466,
          0.15219028039479165
        ],
        [
          9.271752804071631e-06,
          0.8906779225368785,
          -0.23859070167291532,
          -1.0513512328321077
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    },
    {
      "file_path": "./data\\221216/cam4.png",
      "sharpness": 238.46831398479796,
      "transform_matrix": [
        [
          -0.04408322883149452,
          1.0736996157079683e-06,
          0.9217010598234293,
          3.9218489382485333
        ],
        [
          0.9988321554148737,
          -2.714800923986508e-06,
          0.03302779529821646,
          0.15222665375307728
        ],
        [
          2.5376968199949358e-06,
          0.9220806280895403,
          -9.527685110624582e-07,
          -8.314883668141512e-06
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    },
    {
      "file_path": "./data\\221216/cam5.png",
      "sharpness": 247.37153101838447,
      "transform_matrix": [
        [
          -0.044092357623335895,
          -0.23849570270312323,
          0.8903099676721334,
          3.9219213677000724
        ],
        [
          0.9988318282503761,
          -0.008546564291621338,
          0.03191308623705689,
          0.1522949963428471
        ],
        [
          -2.0425493510426027e-06,
          0.8906770559307186,
          0.23859393693001565,
          1.0513583465315044
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    },
    {
      "file_path": "./data\\221216/cam6.png",
      "sharpness": 238.84922712626152,
      "transform_matrix": [
        [
          -0.17189238863861847,
          -0.2345945665909449,
          0.8757251685312653,
          3.823908748242755
        ],
        [
          0.9845221713657459,
          -0.04441854148044574,
          0.1658300250931058,
          0.7134712823353254
        ],
        [
          -4.388140606329475e-06,
          0.8906757635632908,
          0.23859876128564045,
          1.0513892849489217
        ],
        [
          0.0,
          0.0,
          0.0,
          1.0
        ]
      ]
    }
  ]
}

通常の手順の通り、画像ファイルとtransform.jsonをdata以下の所定のフォルダに置きました。

common_device.cuh

カメラのディストーション係数のK2でプロジェクトを検索したところ、ディストーション計算らしい箇所を見つけましたので、radialにK3 * r2 * r2 * r2を追加しました。

そのほかnerf_loader.cuのread_camera_distortionやtestbed_nerf.cuのset_calera_intrinsics, python_api.cuなど、プロジェクトを検索して見つかったカメラのひずみ係数のK2の付近でK3の記述を追加しビルドして実行しました。

testbed.exe

実行はできましたが浮遊物が増加してしまいました。視錐台に沿うように浮遊物が見えたので、像高が大きい部分の誤差が余計に発生してしまったのかもしれません。目論見通りになりませんでした。

ディストーション係数をCOLMAPであらかじめ解析する必要がありそうです。

備考　openCVのキャリブレーションでK3=0にする

openCVのキャリブレーションを使いK2まで出力する方法を試しました。
こちらのスクリプトを部分的に改造しました。cv2.calibrateCameraのflagsを設定するため入力するパラメータの個数に合わせNoneを増やしました。

ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1], None, None, None, None, cv2.CALIB_FIX_K3)