Mugichoko's blog

Mugichoko’s blog

プログラミングを中心としたメモ書き.

RGB-D SLAMを実装する #1

目標

KellerらのRGB-D SLAM[1]が実装したい!と思い立ったので実装していく,というモチベーションの日誌.今回が初稿.ちょっとずつ実装していく.今回は,論文読解とアルゴリズムの整理を行う.

アルゴリズム

論文を読んで,今考えているアルゴリズムは以下のイメージ.ただし,動物体の検出に関しては,今回,ひとまず無視する.

  1. メモリの確保
    • グローバルなIndexマップを持ったある点群(以下は,1点当たりの要素,SSBOとして確保?)
      • 点群の配列 (float×3)
        • GPUに確保するので,一気に大きな配列を用意
      • 信頼度 (float)
        • 現在の点と対応する点がマージされると+1
        • 10以上(これは経験的な値)になると不安定から安定へ
      • 半径 (float)
        • サーフェル描画用
      • タイムスタンプ (int)
    • 4W×4HのIndexマップ(int,WとHは画面サイズの幅と高さ)
      • 全点を投影して,ルックアップテーブル (LUT) として使う
    • W{\times}Hの
      • 頂点群マップ (float×3)
      • 奥行マップ (float×3)
      • 法線マップ (float×3)
  2. 奥行画像の処理
    1. 奥行マップをバイラテラルフィルタでスムージング
    2. 奥行マップを頂点マップに変換
    3. 頂点マップ画像から法線マップの作成
    4. 頂点マップと奥行マップは3段階の解像度のものを用意
      • Lv0:元の解像度,Lv1:0の半分,Lv2:1の半分
      • Mipmapで実装できる?
  3. 位置合わせ(Iterative Closest Point (ICP) アルゴリズム等)
    • 参考文献[2]
  4. Indexマップの生成
    • 点群の配列をIndexマップに単なる点群として投影して,各点のindexを書き込む f:id:Mugichoko:20170725103831j:plain
  5. Indexマップから1点抽出(以下は,そのルール)
    1. {\pm \sigma_{depth}}以上離れている点は候補から外す
    2. 保存された法線と現在の法線との角度が20度以上の点は候補から外す
    3. 残りの点の内,信頼度の最も高い点を選択
    4. もし,3まで満たす点が複数あれば,2の内最近傍の点を使う
  6. 選択した点の統合/追加
    • 現フレームの3次元点と法線を以下の式に基づいて統合し,信頼度マップとタイムスタンプマップも更新する f:id:Mugichoko:20170725103915j:plain
    • ただし,{ \alpha = e^{- \gamma^2  / 2 \sigma^2} }
    • {\gamma}は正規化画像座標系における,光軸中心からの距離,{\alpha = 0.6}
    • 式1は,半径が現在の半径と比べて {(1 + \delta_r) \bar{r}} 以下だった場合のみ行う
      • ただし,{\gamma_r = 1 / 2}
    • 仮に,どの候補も見つからなければ,新しい点として追加される
  7. 点の削除
    1. 長い間,unstableであった点は削除する
    2. 点が統合された場合,それよりも前にある点は全て削除する
      • Indexマップ内で探索
    3. Indexマップ内の周囲に類似した点があれば1つにする

参考文献

  1. M. Keller, D. Lefloch, M. Lambers, S. Izadi, T. Weyrich, and A. Kolb, “Real-Time 3D Reconstruction in Dynamic Scenes Using Point-Based Fusion,” Proc. Int. Conf. on 3D Vision, pp. 1 - 8, 2013.
  2. Kok-Lim Low, “Linear Least-Squares Optimization for Point-to-plane ICP Surface Registration,” Technical Report TR04-004, pp. 1 - 3, 2004.

GLSL #11: Bilateral Filter

目標

前回(下記参照)の予告通り,Bilateral Filterを実装する.

mugichoko.hatenablog.com

実装環境

参考サイト

Shadertoy BETA

上記,参考サイトはOpenGL ESなので,OpenGL用に一部,処理を変更する必要がある.更に,3チャンネルの画像を想定しているので,1チャンネルの奥行画像を想定した処理に変更する.

main.cpp

  • csBilateralFilterクラス Compute Shaderを使って,奥行マップにBilateral Filterをかけるクラス.
  • キー操作 「b」キーを押すと,Bilateral Filterがかかった後の法線マップが表示され,それ以外のキーを押すと,元の法線マップが表示される.
#include <glm/glm.hpp>
#include <glm/ext.hpp>
#include <glm/detail/setup.hpp>
#include <opencv2/opencv.hpp>
#include "OpenGLWrapper/window.h"
#include "OpenGLWrapper/shaderUtil.h"
#include "OpenGLWrapper/modelTex.h"
#include "OpenGLWrapper/rgbdTex.h"

// Bilateral filter
// Ref: https://www.shadertoy.com/view/4dfGDH
class csBilateralFilter : public gl::compShader
{
public:
    csBilateralFilter(
        int w,
        int h,
        const string &compShaderName
    ) : gl::compShader(w, h, 1, GL_R32F, GL_RED, GL_FLOAT, compShaderName)
    {
    }

    void csBilateralFilter::execute(
        GLuint inTexID,
        float sigma = 10.0f,
        float bSigma = 0.1f
    )
    {
        glUseProgram(prog);
        {
            GLuint locSigma = glGetUniformLocation(prog, "sigma");
            glUniform1f(locSigma, sigma);
            GLuint locBSigma = glGetUniformLocation(prog, "bSigma");
            glUniform1f(locBSigma, bSigma);

            // Input texture
            glBindImageTexture(0, inTexID, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);
            // Output texture
            glBindImageTexture(1, texID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R32F);
            glDispatchCompute(width / 32, height / 32, depth);
        }
        glUseProgram(0);
    }
};


// Vertex map calculation(前回と同じであるため割愛)
...

// Normal map calculation(前回と同じであるため割愛)
...

class MyWindow : gl::window
{
private:
    csBilateralFilter *myCSBilateralFilter;
    csCalcVertMap *myCSCalcVertMap;
    csCalcNormMap *myCSCalcNormMap;
    gl::rgbdTex *myRGBDTex;
    gl::modelTex *myModelTex;

public:
    MyWindow(
        int w,
        int h,
        const string &wndName,
        const glm::mat3 &K,    // Intrinsic parameters
        const float factor,   // See: http://vision.in.tum.de/data/datasets/rgbd-dataset/file_formats#intrinsic_camera_calibration_of_the_kinect
        const string &colorImgName,
        const string &depthImgName,
        const string &compBilateralFilterShaderName,
        const string &compVertMapShaderName,
        const string &compNormMapShaderName
    ) : window(w, h, wndName)
    {
        // Read images
        cv::Mat imgC = cv::imread(colorImgName);        // Color image
        cv::Mat imgD = cv::imread(depthImgName, -1);   // Depth image
        imgD.convertTo(imgD, CV_32F);
        imgD /= factor;

        // Upload the RGB-D image to your GPU
        myRGBDTex = new gl::rgbdTex(imgC.cols, imgC.rows);
        myRGBDTex->updateRGBDTex(imgC.data, imgD.data);
        
        // Create computer shaders
        myCSBilateralFilter = new csBilateralFilter(imgC.cols, imgC.rows, compBilateralFilterShaderName);
        myCSCalcVertMap = new csCalcVertMap(imgC.cols, imgC.rows, K, compVertMapShaderName);
        myCSCalcNormMap = new csCalcNormMap(imgC.cols, imgC.rows, compNormMapShaderName);

        // Create a texture polygon to display the results
        myModelTex = new gl::modelTex(imgC.cols, imgC.rows, "../shader/modelTex.vert", "../shader/modelTex.frag");
    }
    ~MyWindow()
    {
        if (myRGBDTex)
        {
            delete myRGBDTex;
            myRGBDTex = NULL;
        }
        if (myCSBilateralFilter)
        {
            delete myCSBilateralFilter;
            myCSBilateralFilter = NULL;
        }
        if (myCSCalcVertMap)
        {
            delete myCSCalcVertMap;
            myCSCalcVertMap = NULL;
        }
        if (myCSCalcNormMap)
        {
            delete myCSCalcNormMap;
            myCSCalcNormMap = NULL;
        }
        if (myModelTex)
        {
            delete myModelTex;
            myModelTex = NULL;
        }
    }

    void render()
    {
        const glm::mat4 rotMat180X = glm::rotate(glm::pi<float>(), glm::vec3(1.0f, 0.0f, 0.0f));
        myModelTex->setModelMat(rotMat180X);

        // Rendering loop
        while (!glfwWindowShouldClose(wnd))
        {
            // Clear color and depth buffers
            glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);


            if (keyID == 'B')
            {
                // Bilateral filter
                myCSBilateralFilter->execute(myRGBDTex->getDepthTexID(), 2.0f, 0.05f);
                // Calculate a vertex map from a depth map
                myCSCalcVertMap->execute(myCSBilateralFilter->getTexID());
            }
            else
            {
                // Calculate a vertex map from a depth map
                myCSCalcVertMap->execute(myRGBDTex->getDepthTexID());
            }
            // Calculate a normal map from a depth map
            myCSCalcNormMap->execute(myCSCalcVertMap->getTexID());

            // Show the results
            myModelTex->render(myCSCalcNormMap->getTexID());

            // Swap front and back buffers
            glfwSwapBuffers(wnd);
            // Poll for and process events
            glfwPollEvents();
        }
    }
};

int main()
{
    glm::mat3 K(525.0f, 0.0f, 0.0f, 0.0f, 525.0f, 0.0f, 319.5f, 239.5f, 1.0f);

    MyWindow wnd(
        640, 480,
        "GLFW 3 Window",
        K, 5000.0f,
        "../data/1311868164.363181.png",
        "../data/1311868164.338541.png",
        "../shader/bilateralFilter.comp",
        "../shader/calcVertMap.comp",
        "../shader/calcNormMap.comp"
    );
    wnd.render();

    return 0;
}

bilateralFilter.comp

現状,カーネルサイズ「mSize」は固定.

// Ref: https://www.shadertoy.com/view/4dfGDH
#version 430

layout (binding = 0, r32f) readonly uniform image2D dataIn;        // Input depth map
layout (binding = 1, r32f) writeonly uniform image2D dataOut;  // Output depth map

uniform float sigma;
uniform float bSigma;
const int mSize = 7;

layout (local_size_x = 32, local_size_y = 32) in;

float normPdf(float x, float sigma)
{
    return 0.39894 * exp(-0.5 * x * x / (sigma * sigma)) / sigma;
}

void main(void)
{
    ivec2 u = ivec2(gl_GlobalInvocationID.xy);
    float c = imageLoad(dataIn, u).r;

    const int kSize = (mSize - 1) / 2;
    float kernel[mSize];
    float finalColor = 0.0;

    // Create the 1-D kernel
    float Z = 0.0;
    for (int i = 0; i <= kSize; ++i)
    {
        kernel[kSize + i] = kernel[kSize - i] = normPdf(float(i), sigma);
    }

    float cc;
    float factor;
    float bZ = 1.0 / normPdf(0.0, bSigma);
    // Read out the texels
    for (int x = -kSize; x <= kSize; ++x)
    {
        for (int y = -kSize; y <= kSize; ++y)
        {
            cc = imageLoad(dataIn, u + ivec2(x, y)).r;
            factor = normPdf(cc - c, bSigma) * bZ * kernel[kSize + y] * kernel[kSize + x];
            Z += factor;
            finalColor += factor * cc;
        }
    }

    imageStore(dataOut, u, vec4(finalColor / Z, 0.0, 0.0, 1.0));
}

結果

以下に,Bilateral Filterをかける前後の法線マップを示す,尚,かける前の画像は前回の結果である.この通り,目標達成だ!Bilateral Filterをかけた結果は,法線マップの見た目が滑らかになる.

現状,Bilateral Filterのカーネルサイズが固定だが,将来的には可変にしたい.これは,シェーダのmain関数内に固定長の配列が宣言されていることが原因なのだが,何かいい方法はないものだろうか…?

f:id:Mugichoko:20170725001059j:plain f:id:Mugichoko:20170725003714j:plain

GLSL #10: RGB-D画像からの頂点・法線マップ生成

目標

前回(下記参照)は,Compute Shaderを使って,3DCGを描画した際に生成されるZバッファの線形化を行うプログラムを作成した.今回は,RGB-D画像を用いて,頂点・法線マップを生成するプログラムを作成する.尚,RGB-D画像は,Zバッファの様に非線形化は行われていないため,線形化する必要はありません.

mugichoko.hatenablog.com

実装環境

処理内容

流れ

  1. RGB-D画像読み込み
    • TUM RGB-D Dataset等から適当なものをダウンロードしておく
    • もしくは,自らKinectやXtionから取得しておく/実時間で取得する
  2. Compute Shaderを使って頂点マップを生成
    • 頂点マップとは,RGB-D画像の各画素に対応するfloat型3チャンネルの画像
  3. Compute Shaderを使って法線マップを生成
    • 法線マップとは,RGB-D画像の各画素に対応するfloat型3チャンネルの画像
  4. 生成した頂点・法線マップを描画

頂点マップの計算

参考文献1の3章によると,以下の式によって,{i}フレーム目の奥行マップ{{\mathcal D}_i}から{i}フレーム目の頂点マップ{{\mathcal V}_i}が生成できる.

{
{\boldsymbol v}_i ({\boldsymbol u}) = {\mathcal D}_i({\boldsymbol u}) {\boldsymbol K}^{-1}_i ({\boldsymbol u}^{\rm T}, 1)^{\rm T} \in {\mathbb R}^3
}

この時,{{\boldsymbol u} = (x, y)^{\rm T}}は画像上の2次元位置,{{\boldsymbol v}_i} = (X, Y, Z)^{\rm T}は求めたい頂点マップ{{\mathcal V}_i}{{\boldsymbol u}}に対応する頂点,{{\mathcal D}_i({\boldsymbol u})}{{\boldsymbol u}}に対応する奥行値,{{\boldsymbol K}}は事前の校正により既知である内部パラメータ行列を表す.

{\in {\mathbb R}^3}」の部分は,出力される{{\boldsymbol v}_i ({\boldsymbol u})}が実数の3次元ベクトルであることを表しているだけである.尚,ここでは,1フレーム(1枚の画像)のみ入力することを考えて,{i}は無視する.以上より,シェーダで記述する式は以下のようになる.

{
{\boldsymbol v} ({\boldsymbol u}) = {\mathcal D}({\boldsymbol u}) {\boldsymbol K}^{-1}\ ({\boldsymbol u}^{\rm T}, 1)^{\rm T}
}

法線マップの計算

ここでも参考文献1の3章に倣って,Central Differenceを取ることで,頂点マップ{{\mathcal V}}から法線マップ{{\mathcal N}}を計算する.具体的には,法線マップ{{\mathcal N}}の要素{{\boldsymbol n}}以下の式より求める.

{
{\boldsymbol v}_x = {\mathcal V}(x+1,y) - {\mathcal V}(x-1,y)
}

{
{\boldsymbol v}_y = {\mathcal V}(x,y+1) - {\mathcal V}(x,y-1)
}

{
{\boldsymbol n}({\boldsymbol u}) = {\rm norm} \left( {\boldsymbol v}_x \times {\boldsymbol v}_y \right)
}

尚,{{\rm norm}({\boldsymbol x})}は,ベクトル{{\boldsymbol x}}を正規化する関数である.

サンプルプログラム

前回までのプログラムをそのまま利用し,以下の様に書き換えます.

main.cpp

  • csCalcVertMapクラス Compute Shaderを使って,奥行マップから頂点マップを生成するためのクラス.
  • csCalcNormMapクラス Compute Shaderを使って,頂点マップから法線マップを生成するためのクラス.
  • キー操作 「2」「3」キーを押すと,それぞれ頂点マップ,法線マップが表示される.それ以外を押すとRGB画像が表示される.
#include <glm/glm.hpp>
#include <glm/ext.hpp>
#include <glm/detail/setup.hpp>
#include <opencv2/opencv.hpp>
#include "OpenGLWrapper/window.h"
#include "OpenGLWrapper/shaderUtil.h"
#include "OpenGLWrapper/modelTex.h"
#include "OpenGLWrapper/rgbdTex.h"

// Vertex map calculation
class csCalcVertMap : public gl::compShader
{
private:
    const glm::mat3 K; // Intrinsic parameters

public:
    csCalcVertMap(
        int w,
        int h,
        const glm::mat3 &K,
        const string &compShaderName
    ) : gl::compShader(w, h, 1, GL_RGBA32F, GL_RGB, GL_FLOAT, compShaderName), K(K)
    {

        glUseProgram(prog);
        {
            // Update the uniform
            GLuint locInvK = glGetUniformLocation(prog, "invK");
            glUniformMatrix3fv(locInvK, 1, GL_FALSE, glm::value_ptr(glm::inverse(K)));
        }
        glUseProgram(0);
    }

    void csCalcVertMap::execute(
        GLuint inTexID
    )
    {
        glUseProgram(prog);
        {
            // Input texture
            glBindImageTexture(0, inTexID, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);
            // Output texture
            glBindImageTexture(1, texID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);
            glDispatchCompute(width / 32, height / 32, depth);
        }
        glUseProgram(0);
    }
};

// Normal map calculation
class csCalcNormMap : public gl::compShader
{
public:
    csCalcNormMap(
        int w,
        int h,
        const string &compShaderName
    ) : gl::compShader(w, h, 1, GL_RGBA32F, GL_RGB, GL_FLOAT, compShaderName)
    {
    }

    void csCalcNormMap::execute(
        GLuint inTexID
    )
    {
        glUseProgram(prog);
        {
            // Input texture
            glBindImageTexture(0, inTexID, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RGBA32F);
            // Output texture
            glBindImageTexture(1, texID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);
            glDispatchCompute(width / 32, height / 32, depth);
        }
        glUseProgram(0);
    }
};

class MyWindow : gl::window
{
private:
    csCalcVertMap *myCSCalcVertMap;
    csCalcNormMap *myCSCalcNormMap;
    gl::rgbdTex *myRGBDTex;
    gl::modelTex *myModelTex;

public:
    MyWindow(
        int w,
        int h,
        const string &wndName,
        const glm::mat3 &K,    // Intrinsic parameters
        const float factor,   // See: http://vision.in.tum.de/data/datasets/rgbd-dataset/file_formats#intrinsic_camera_calibration_of_the_kinect
        const string &colorImgName,
        const string &depthImgName,
        const string &compVertMapShaderName,
        const string &compNormMapShaderName
    ) : window(w, h, wndName)
    {
        // Read images
        cv::Mat imgC = cv::imread(colorImgName);        // Color image
        cv::Mat imgD = cv::imread(depthImgName, -1);   // Depth image
        imgD.convertTo(imgD, CV_32F);
        imgD /= factor;

        // Upload the RGB-D image to your GPU
        myRGBDTex = new gl::rgbdTex(imgC.cols, imgC.rows);
        myRGBDTex->updateRGBDTex(imgC.data, imgD.data);
        
        // Create computer shaders
        myCSCalcVertMap = new csCalcVertMap(imgC.cols, imgC.rows, K, compVertMapShaderName);
        myCSCalcNormMap = new csCalcNormMap(imgC.cols, imgC.rows, compNormMapShaderName);

        // Create a texture polygon to display the results
        myModelTex = new gl::modelTex(imgC.cols, imgC.rows, "../shader/modelTex.vert", "../shader/modelTex.frag");
    }
    ~MyWindow()
    {
        if (myRGBDTex)
        {
            delete myRGBDTex;
            myRGBDTex = NULL;
        }
        if (myCSCalcVertMap)
        {
            delete myCSCalcVertMap;
            myCSCalcVertMap = NULL;
        }
        if (myCSCalcNormMap)
        {
            delete myCSCalcNormMap;
            myCSCalcNormMap = NULL;
        }
        if (myModelTex)
        {
            delete myModelTex;
            myModelTex = NULL;
        }
    }

    void render()
    {
        const glm::mat4 rotMat180X = glm::rotate(glm::pi<float>(), glm::vec3(1.0f, 0.0f, 0.0f));
        myModelTex->setModelMat(rotMat180X);

        // Rendering loop
        while (!glfwWindowShouldClose(wnd))
        {
            // Clear color and depth buffers
            glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

            // Calculate a vertex map from a depth map
            myCSCalcVertMap->execute(myRGBDTex->getDepthTexID());

            // Calculate a normal map from a depth map
            myCSCalcNormMap->execute(myCSCalcVertMap->getTexID());


            // Show the results
            switch (keyID) {
            case '2':
                // Vertex map
                myModelTex->render(myCSCalcVertMap->getTexID());
                break;
            case '3':
                // Normal map
                myModelTex->render(myCSCalcNormMap->getTexID());
                break;
            default:
                // Color image
                myModelTex->render(myRGBDTex->getColorTexID());
                break;
            }

            // Swap front and back buffers
            glfwSwapBuffers(wnd);
            // Poll for and process events
            glfwPollEvents();
        }
    }
};

int main()
{
    glm::mat3 K(525.0f, 0.0f, 0.0f, 0.0f, 525.0f, 0.0f, 319.5f, 239.5f, 1.0f);

    MyWindow wnd(
        640, 480,
        "GLFW 3 Window",
        K, 5000.0f,
        "../data/1311868164.363181.png",
        "../data/1311868164.338541.png",
        "../shader/calcVertMap.comp",
        "../shader/calcNormMap.comp"
    );
    wnd.render();

    return 0;
}

calcVertMap.comp

csCalcVertMapクラスで利用されるシェーダプログラム.頂点マップの計算をシェーダプログラムとして,具体的に書きだした.

#version 430

layout (binding = 0, r32f) readonly uniform image2D dataIn;            // Depth map
layout (binding = 1, rgba32f) writeonly uniform image2D dataOut;   // Vertex map

uniform mat3 invK;      // Inverse of the intrinsic matrix

layout (local_size_x = 32, local_size_y = 32) in;

void main(void)
{
    ivec2 u = ivec2(gl_GlobalInvocationID.xy);
    float z = imageLoad(dataIn, u).r;

    vec3 v = z * invK * vec3(u, 1.0);

    imageStore(dataOut, u, vec4(v, 1.0));
}

calcNormMap.comp

csCalcNormMapクラスで利用されるシェーダプログラム.法線マップの計算をシェーダプログラムとして,具体的に書きだした.

// REF: http://answers.opencv.org/question/82453/calculate-surface-normals-from-depth-image-using-neighboring-pixels-cross-product/
#version 430

layout (binding = 0, rgba32f) readonly uniform image2D dataIn;     // Depth map
layout (binding = 1, rgba32f) writeonly uniform image2D dataOut;   // Normal map

layout (local_size_x = 32, local_size_y = 32) in;

void main(void)
{
    vec3 xyz[4] = {
        imageLoad(dataIn, ivec2(gl_GlobalInvocationID.xy) - ivec2(1, 0)).rgb,
        imageLoad(dataIn, ivec2(gl_GlobalInvocationID.xy) + ivec2(1, 0)).rgb,
        imageLoad(dataIn, ivec2(gl_GlobalInvocationID.xy) - ivec2(0, 1)).rgb,
        imageLoad(dataIn, ivec2(gl_GlobalInvocationID.xy) + ivec2(0, 1)).rgb
    };

    vec3 vecX = xyz[1] - xyz[0];
    vec3 vecY = xyz[3] - xyz[2];

    imageStore(
        dataOut,
        ivec2(gl_GlobalInvocationID.xy),
        vec4(normalize(cross(vecX, vecY)), 1.0)
        );
}

結果

以下に示す,RGB画像,奥行マップ,頂点マップ,法線マップにある通り,目標達成!尚,奥行マップを表示させるようにはしていないので,これは入力画像をそのままここに載せている.

正直なところ,TUM RGB-D Datasetでのカメラ座標系の向き,OpenGLのカメラ座標系との関係が分かっていないので,画像だ出ただけ,と言えばそれだけ.座標軸が想定通りか,今後,検証する必要がある.また,法線マップがガタガタなので,次回はBilateral Filterを実装してみようと思う.

f:id:Mugichoko:20170724215954j:plain f:id:Mugichoko:20170724220027p:plain f:id:Mugichoko:20170724220044j:plain f:id:Mugichoko:20170725001059j:plain

参考文献

  1. M. Keller, D. Lefloch, M. Lambers, S. Izadi, T. Weyrich, and A. Kolb, “Real-Time 3D Reconstruction in Dynamic Scenes Using Point-Based Fusion,” Proc. Int. Conf. on 3D Vision, pp. 1 - 8, 2013.

オーストリア留学に際してのビザ申請

訪問研究員 (Visiting Researcher) として,オーストリア (Austria) のグラーツ (Graz) に1年間の留学が決まった.6カ月以上の滞在になるのでビザ申請が必要だ.オーストリア留学がマイナーなのか,研究員としての留学がマイナーなのかは不明だが,特段のまとまった情報が見つからなかったのでここにまとめる.

私の基本情報(状況)

  • 居住地:関東
  • 国内での立場:日本学術振興会特別研究員PD
    • 研究期間の半分まで(つまり,1.5年まで)の海外留学が許されている
  • 留学先:TUGraz
    • 留学先の教授に直接話しかけたところ,私の研究に興味を持って下さり,留学したい旨を伝え,了承頂いた
  • 留学期間:2017年8月~2018年7月の1年間(予定)
  • 留学先での立場:訪問研究員 (Visiting Researcher)
  • 備考:妻(無職)を留学先に連れて行く

ビザD取得(私の場合)

必要なビザの種類

上記のような状況の私は,日本国内でビザD取得(6カ月以内の滞在)の後,オーストリア国内で在留許可申請(6カ月以上の滞在)が必要らしい.こういった情報を入手した経緯や具体的な方法をまとめた.

必要な書類

まず,必要な書類だが,オーストリア大使館のWebサイトにある「必要提出書類(添付書類)」を参照のこと.

ビザDを取得した後の感想だが,結構色々な理由で突っ返され,何度も大使館を訪問することになる.どういうところで引っかかるかは,必要書類と以下のビザD取得までの流れとを照らし合わせて確認されることをお勧めしたい.

取得までの流れ

どうやら大阪と名古屋での申請に関してはビザ申請専用のWebサイトがあるようだが,なぜか東京版は存在しない. そこで,基本的な情報はオーストリア大使館 ビザインフォメーションで得る.しかし,抜けがあるようなので,私の体験をケーススタディとして以下にまとめる.

  1. オーストリア大使館 ビザインフォメーションを確認
  2. 問い合わせフォームで私の基本情報(上記)を伝え,必要なビザや申請手続きの詳細を問い合わせ
    • ビザDが必要だということが判明.「詳細はWebサイト(上記)をご確認下さい.」とのこと
    • 正直,大使館のWebサイトが分かりにくかったので,聞けば分かる精神で問い合わせたという経緯だったため,結局,在京オーストリア大使館を直接訪問することに(地図
  3. オーストリア大使館を訪問(1回目)
    • 留学先からのInvitation Letter留学先での居住証明書が必要ということ,6カ月以上の滞在には,オーストリア国内での在留許可申請が必要であることが分かった
  4. 大使館で必要と言われた書類を方々から取得
    • TUGrazのWelcome Centerに問い合わせて,滞在先証明書(サイン or スタンプ付)を取得
    • 留学先の教授に問い合わせて,Invitation Letter(サイン or スタンプ付)を取得
  5. 4の書類を在京オーストリア大使館にメールで送付し,書類に不備がないか確認を依頼
    • 「メールでの確認は行っていない.」と断られ,直接,大使館を訪問することに
  6. オーストリア大使館を訪問(2回目)するも,下記の理由で突き返される
    • 申請書の不備(これはその場で指摘を受けながら直した)
    • 雇用証明(英語可)がない(前回はInvitation Letterだけでいいと言っていたのに…)
    • 「充分な資金の証明、通帳(6ヶ月分の残高証明)」の不備
      • 新規に通帳を作ったので,過去数カ月分しかなく「6ヶ月分の残高証明」という要件に引っかかる
      • 前回は,記された金額が十分だと思われる場合はそれでよい,と言っていたのに…
  7. 学振に雇用証明に相当する書類(採用証明書様式10を提出)を取り寄せる
  8. オーストリア大使館を訪問(3回目)し,前回に不備のあった書類を提出
    • 過去の通帳を用意して,コピーと共に持って行ったが「最近の入出金の記録がない」という不備を指摘される(前回はその突っ込みなかったんですけど…)
    • 隣で「再三,電話で写真の規格について確認して,その時に言われた通りのものを持ってきたのに,ここに来てまた別のことを言うじゃないですか!」と怒っておられるのを「そうだそうだ!」と心の中で思う我々夫婦
  9. 通帳の不備だけだったので,後日,大使館に郵送
  10. 数日後にビザ取得
    • 妻が取りに行ってくれたので,私は委任状を書いた
    • 委任状の書き方は外務省のものを参考にした

在京オーストリア大使館の場所

時間的感覚

仕事などの都合もあったが,ビザDの申請から取得までおおよそ2カ月強かかった.

  • メールでの最初の問い合わせ(2月4日)
  • 第1回訪問(4月26日)← 申請開始
  • 第2回訪問(6月21日)
  • 第3回訪問(7月5日)
  • ビザ取得(7月11日)

オーストリア大使館のスタンス

感想に近いが,事前に把握しておくとイライラも少ない.

  • 申請作業のための訪問は要予約
  • 電話でよく応対してくれる
    • ただし,窓口での説明と結構異なる
  • 「基本的にWebページを確認してくれ.そこに書いてある.」と言われる
  • 情報は小出し(これが一番困る)
  • 基本的にビザ申請に関してのみ扱っている.
    • 在留許可申請に関しては現地が判断する,とのこと
    • 大使館から下手なことは言えない,そりゃそうか

ドロネー三角形分割:三角形に所属する各点のIDを取得

目標

前回ランダムに生成した2次元点を入力として,ドロネー三角形分割を行い,結果を描画した. 今回は,そこで用いたFade2Dを利用して,各三角形にどの2次元点が使われたかを特定する方法を実装する.

mugichoko.hatenablog.com

実装環境

プログラム

決め手は,「Point2::setCustomIndex()」によって各点にあらかじめIDを割り振っておき,ドロネー三角形分割の後に「Point2::getCustomIndex()」で,各三角形のIDを取得するということ.

#include <iostream>
using namespace std;
#include <random>
#include <Fade_2D.h>
namespace fade = GEOM_FADE2D;

// Constants
const int PTS_NUM(100);
const int PTS_RANGE(600);
const string WND_NAME("Delaunay Fade 2D");


int main()
{
    // ----- Generate 2D points
    // Mersenne twister (random value generator)
    mt19937 mt;
    uniform_int_distribution<int> distribution(0, ::PTS_RANGE - 1);
    
    vector<fade::Point2> pts;
    for (int i = 0; i < ::PTS_NUM; ++i)
    {
        int x = distribution(mt);
        int y = distribution(mt);
        fade::Point2 p(x, y);
        p.setCustomIndex(i);
        pts.push_back(p);
    }

    // ----- Initialize Fade 2D
    fade::Fade_2D dt((unsigned int)pts.size());
    vector<fade::Point2 *> vDtVertPtr;
    dt.insert(pts, vDtVertPtr);

    // ----- Draw Delaunay triangles
    cv::Mat img = cv::Mat::zeros(::PTS_RANGE, ::PTS_RANGE, CV_8UC3);

    vector<fade::Triangle2 *> vDtTri;
    dt.getTrianglePointers(vDtTri);
    for (auto tri = vDtTri.begin(); tri != vDtTri.end(); ++tri)
    {
        fade::Point2 *pt[3] = { (*tri)->getCorner(0), (*tri)->getCorner(1), (*tri)->getCorner(2) };
        cout << "Point ID: " << pt[0]->getCustomIndex() << ", " << pt[1]->getCustomIndex() << ", " << pt[2]->getCustomIndex() << ", " << endl;
    }

    return 0;
}

実行結果

コマンドプロンプトに各三角形に含まれる点のIDが表示された.

f:id:Mugichoko:20170606120102j:plain

Python 3.6 + OpenCV 3.2 on Windows 10 64-bit

目標

PythonOpenCVを使ったプログラミングを,Visual Studioでできるようにする. 今回は,Visual StudioPythonはインストール済みとして話を進めます. 代わりに,私が経験した,PythonOpenCVのインストール時のトラブルについて紹介します.

実装環境

遭遇した問題

PythonOpenCVのインストールの仕方を参考にしてインストールしてみた. 例えば,以下のリンクが参考になる.

prpr.hatenablog.jp

この際,正しくインストールできたかどうかを確かめるために,以下のコマンドをコマンドプロンプトから実行しろ,という記事によく出くわす.

C:\Users\XXXX> python
>>> import cv2
>>> print(cv.2__version__)

しかし,この2行目のコードを実行すると…

Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ImportError: DLL load failed: The specified module could not be found.

というエラーメッセージが表示される. どうやらDLLが見つからないようだ. この問題に関して,以下に列挙する,ネット上の種々の解決方法を試してみたが,全く効果がなかった.

解決策

解決策1:conda-forge

以下をコマンドプロンプトにコピペしてエンターキーで実行. conda-forgeを使ったPythonOpenCVのインストール用コマンドだ. しばらくしてインストールが終わる.

※この方法だと,ContribなしのOpenCVがインストールされると思われる.

conda install -c conda-forge opencv=3.2.0

参考

www.codesofinterest.com

解決策2:pipと.whil

ここにある,自分の環境にあった.whlファイルをダウンロードし,そのファイルを置いたパスで,コマンドプロンプトから以下を実行する.

※この方法だと,ContribありのOpenCVがインストールされる.

pip install opencv_python-3.1.0+contrib_opencl-cp35-cp35m-win_amd64.whl

私の場合,実装環境にある通りなので「opencv_python-3.2.0+contrib-cp36-cp36m-win_amd64.whl」を選択した.

参考

qiita.com

解決策3:Python3.5へのダウングレード(未確認)

下記のリンクにある回答によると,Python 3.5にダウングレードすれば,以下のコマンドでインストールできるようだ.

conda install -c https://conda.binstar.org/menpo opencv3

stackoverflow.com

テストコード

import numpy as np
import cv2

img = cv2.imread('me.jpg')
cv2.imshow('me', img)
cv2.waitKey()

実行結果

f:id:Mugichoko:20170531023511j:plain

GLSL #9: Zバッファの線形化

目標

前回(下記参照)までの実装を少し変更して,Compute Shaderを使って,3DCGを描画した際に生成されるZバッファの線形化を行うプログラムを作成する.

mugichoko.hatenablog.com

Zバッファの線形化とは?

Zバッファは,実際の奥行値から非線形変換が行われて0~1の大きさに正規化されており,実際の空間内の距離とは異なる.そこで,以下の式に基づいてZバッファを線形化(本来の奥行値)に戻す. 尚, Z'は出力, ZはZバッファの値, nは3DCGを描画する際に用いたNear平面までの距離, fはFar平面まで距離. 詳細は下記を参照のこと.

 {Z' = (2~n~f) / (f + n - (2~Z - 1) ~ (f - n))}

learnopengl.com

実装環境

レンダリングの流れ

  1. パス1:四角形ポリゴンに向いたカメラが奥行方向に0.5~5.5の範囲で移動する
  2. パス2:Compute Shaderに1のZバッファを渡し,上記の式を適用する
  3. パス3:パス1の結果はFBOに書き込まれているので,バックバッファにコピーする
  4. パス4:Compute Shaderで書き込んだテクスチャのデータをメインメモリに読み込んで,その中心の画素の値をコマンドプロンプトに表示する

※つまり,コマンドプロンプトには0.5~5.5の値が表示されるはずである!

サンプルプログラム

前回までのプログラムをそのまま利用し,以下の様に書き換えます.ただし,ここにあるバグを修正した.

main.cpp

//
// main.cpp
//
#include <glm/glm.hpp>
#include <glm/ext.hpp>
#include "OpenGLWrapper/window.h"
#include "OpenGLWrapper/shaderUtil.h"
#include "OpenGLWrapper/modelTex.h"
#include "OpenGLWrapper/fbo.h"

class imCompShader : public gl::compShader
{
public:
    imCompShader(
        int w,
        int h,
        const string &compShaderName,
        GLint internalFormat = GL_RGBA8,
        GLenum format = GL_RGBA,
        GLenum type = GL_UNSIGNED_BYTE
    ) : gl::compShader(w, h, 1, internalFormat, format, type, compShaderName)
    {
    }

    void imCompShader::execute(
        GLuint inTexID,
        float near,
        float far
    )
    {
        glUseProgram(prog);
        {
            // Uniforms
            glUniform1f(glGetUniformLocation(prog, "zNear"), near);
            glUniform1f(glGetUniformLocation(prog, "zFar"), far);

            // Input texture
            glBindImageTexture(0, inTexID, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);
            // Output texture
            glBindImageTexture(1, texID, 0, GL_FALSE, 0, GL_WRITE_ONLY, internalFormat);

            glDispatchCompute(width / 32, height / 32, depth);
        }
        glUseProgram(0);
    }
};

class MyWindow : gl::window
{
private:
    gl::modelTex *myTex;
    gl::fbo *myFBO;
    imCompShader *myZBuffComp;

    const float near;
    const float far;

    void path1()
    {
        glEnable(GL_DEPTH_TEST);

        glBindFramebuffer(GL_DRAW_FRAMEBUFFER, myFBO->getFBOID());
        {
            // Clear color and depth buffers
            glClearColor(0.3f, 0.5f, 0.8f, 0.0f);
            glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

            // Draw
            glm::vec3 eye(0.0f, 0.0f, 0.5f + 5.0f * abs(sin(glfwGetTime())));
            glm::vec3 center(0.0f);
            glm::vec3 up(0.0f, 1.0f, 0.0f);
            glm::mat4 V = glm::lookAt(eye, center, up);
            glm::mat4 P = glm::perspective(3.14f / 4.0f, (float)width / height, near, far);
            myTex->setViewMat(V);
            myTex->setProjMat(P);
            myTex->render();
        }
        glBindFramebuffer(GL_DRAW_FRAMEBUFFER, 0);
    }
    void path2()
    {
        myZBuffComp->execute(myFBO->getDepthTexID(), near, far);
    }
    void path3()
    {
        // Framebuffer to back buffer copy
        glBindFramebuffer(GL_READ_FRAMEBUFFER, myFBO->getFBOID());
        glBindFramebuffer(GL_DRAW_FRAMEBUFFER, 0);
        glBlitFramebuffer(
            0, 0, width, height, 0, 0, width, height,
            GL_COLOR_BUFFER_BIT, GL_NEAREST
        );
        glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
    }
    void path4()
    {
        float *pixels = new float[width * height];

        glBindTexture(GL_TEXTURE_2D, myZBuffComp->getTexID());
        glGetTexImage(GL_TEXTURE_2D, 0, myZBuffComp->getFormat(), myZBuffComp->getType(), pixels);

        cout << pixels[width * (height + 1) / 2] << endl;
        delete[] pixels;
    }

public:
    MyWindow(
        int w,
        int h,
        float near,
        float far,
        const string &name,
        const string &texFileName,
        const string &vertShaderName,
        const string &fragShaderName,
        const string &compShaderName
    ) : window(w, h, name), near(near), far(far)
    {
        myTex = new gl::modelTex(texFileName, vertShaderName, fragShaderName);
        myFBO = new gl::fbo(w, h, GL_RGBA, GL_RGBA, GL_UNSIGNED_BYTE);
        myZBuffComp = new imCompShader(w, h, compShaderName, GL_R32F, GL_RED, GL_FLOAT);
    }
    ~MyWindow()
    {
        if (myTex) delete myTex;
        if (myFBO) delete myFBO;
        if (myZBuffComp) delete myZBuffComp;
    }

    void render()
    {
        // Rendering loop
        while (!glfwWindowShouldClose(wnd))
        {
            path1();
            path2();
            path3();
            path4();
            
            // Swap front and back buffers
            glfwSwapBuffers(wnd);
            // Poll for and process events
            glfwPollEvents();
        }
    }
};

int main()
{
    MyWindow wnd(
        512, 512, 0.1f, 100.0f,
        "GLFW 3 Window",
        "../data/me.jpg",
        "../shader/modelTex.vert",
        "../shader/modelTex.frag",
        "../shader/zbuff.comp"
    );
    wnd.render();

    return 0;
}

zbuff.comp

//
// zbuff.comp
//
#version 430

layout (binding = 0, r32f) readonly uniform image2D dataIn;
layout (binding = 1, r32f) writeonly uniform image2D dataOut;

layout (local_size_x = 32, local_size_y = 32) in;

uniform float zNear;
uniform float zFar;

void main(void)
{
    ivec2 pos = ivec2(gl_GlobalInvocationID.xy);
    vec4 val = imageLoad(dataIn, pos);

    float z = (2.0 * zNear * zFar) / (zFar + zNear - (val.r * 2.0 - 1.0) * (zFar - zNear));

    imageStore(
        dataOut,
        pos,
        vec4(z, 0.0, 0.0, 0.0)
        );
}

結果

カメラが0.5の位置に近づいた時. f:id:Mugichoko:20170524172755j:plain

カメラが5.5の位置に近づいた時. f:id:Mugichoko:20170524172750j:plain

尚,0.5,5.5といった固定値を設定した場合にも,それぞれ0.5,5.5とコマンドプロンプトに表示された. ちなみに,線形化をしない場合,つまり「zbuff.comp」内で「float z = val.r;」とする場合,0.9辺りを前後した.