聞きかじりめも

主にC++やメディア処理技術などに関して気付いたことを書いていきます.ここが俺のメモ帳だ!

デルタ型3Dプリンタキット「Zonestar D810」を組み立てた(そのに)

 前回は到着まで書いたので,ここから組み立てに入ります(記事作成時点で動いてます).

f:id:Mzawa2:20160718223807j:plain

とは言っても殆どマニュアル通りなのでそんなに書くことはないです.中華3Dプリンタの場合はマニュアルが途中で終わっているとか部品表がないなんて事もあるらしいですが,D810ではそんな事はありませんでした.動画解説もあるのでイメージしやすいし,丁寧な方でしょう.ここではつまづいたところを雑記しておきます.

1.部品が足りない

・・・まあある程度予想はしてましたが,自分の場合テフロンチューブと細かなネジ類が足りませんでした.致命的な部品ではないので良しとしましょう.ちなみにテフロンチューブについて問い合わせたら,

Hello ,Dear 
Could you pls send more photos ? we had checked very carefully before shipment ,it should be included ,teflon tube is a white tube ,it looks the same color with the packaging ,pls double check and advise ,thanks ~ 

 と返ってきました.Zonestarは丁寧な業者ですし確かに入れたはずと言っているので,それを信じるならば輸入の際の開封確認作業時に失くされた可能性が高いです.テフロンチューブだけもう一回送ってよこせなんて時間かかってアホらしいので,近くのホームセンターでΦ4のエアチューブを購入して中にシリコンスプレーを吹き付けることで代用しました.ついでに足りないネジ類も購入し,先方にメッセージを送信.

Hello, seller. The picture is the rest of hardware setup D810 documentation slides. I cannot find teflon tube, but it's OK. Maybe it was lost while checking goods imported. I'll buy another tube.And I have finished hardware installation. Thanks a lot! 

f:id:Mzawa2:20160719032602j:plain

 これにて事なきを得ました.

2.3Dプリントパーツが割れる

D810はパーツの一部を3Dプリントで作ってありますが,3Dプリント品は積層方向への荷重に弱く,またプリント不良があるとそこから一気に割れます.自分の場合,キャリッジのアームをネジ止めする部分が割れちゃいました.(写真は割れてないやつ)

f:id:Mzawa2:20160719023637j:plain

精度に非常に関係する部分なので焦りましたが,瞬間接着剤で固めたらなんとかなりました.但し,プリント不良の部分がバネみたいになってしまったせいでキャリブレーションに影響が出てしまっています.このままでは問題なので,送ってもらったSTLファイルを早くプリントして交換する予定です.

3.平面出しができない

電装系も組みあがり,さてキャリブレーションしようと動かしたところ,ちょんちょんとビルドプレートをクリックする間隔が明らかに前後で異なり,ついには動かす段階でガツンとプレートにぶつかってキャリブレーションが失敗しました.明らかにホットエンドに対してプレートが斜めになっているのが原因ですが,プレートのネジ調整の範囲では直りません.

真の原因は支柱にありました.上部のコーナーパーツの奥までしっかりと挿入されていなかったため,そこでベルトの長さが変わってしまったようです.ただでさえキツキツな作りなのに上部コーナーパーツはちゃんと差さっているか見えないため,奥まで入り込んでいないことに気付きづらいです.もうちょっと隙間広くして入れやすくしてくれてもいいのに・・・

4.フィラメントが送れない

キャリブレーションが正常に動くようになり,プリントしないテスト用のGコードもうまく動いてくれたので,次はエクストルーダをチェックしたところ全く動きません.X軸モータとケーブル差し込み位置を入れ替えたところちゃんと動いたので,モータやケーブルが原因ではないようです.マルチメータで測ったところ,どうやらメイン基板の時点で出力信号が出ていないようです.ヤバい.壊れてるのか・・・?

youtu.be

これ以上AliExpress担当に技術的トラブルを報告するのも良くないので,テクニカルサポートエンジニアにメールを送りました.すると,

Hi Friend,Do you have heat the hotend before you load the filament?? In order to protect the extruder motor, the printer will prohibit rotate extruder motor when the hotend temperature less than 150 degree.
 な,なんだってー!(画像略)
モータの安全のため150℃以上になるまで動かないとかいうFAQにも書いてない仕様があるらしい.試しに180℃に設定して同じく"Change Filament"コマンドを実行するとちゃんと動きました.

 5.プリントしようとするとファームウェアが落ちる

上記と同時に判明した問題です.プリント時のオートレベリングまでは問題なく動くのですが,その直後に強制終了してファームウェアが再起動してしまう問題にぶち当たり,途方に暮れてしまいました.

youtu.be

見た感じ,ベルト用モータとエクストルーダ用モータとホットエンドを同時に起動した時に落ちるっぽいので,電力不足が原因に思えました.なぜかというと,ホットエンドやエクストルーダを使用しないテスト用Gコードはちゃんと動くからです.メインボード電源は12V5Aなので,試しに12V10Aのヒーテッドベッド用電源をメインボードに挿すという危険極まりない行為をしてみたところ,ちゃんとプリントできるようになりました.

でもこれをやっちゃうと,下のヒーテッドベッドが使えなくなっちゃうんだよね・・・どうしよう.

結果

youtu.be

f:id:Mzawa2:20160719032155j:plain

だいぶいい感じです.側面ががたついてますが,多分これは割れたキャリッジのせいでしょう.ついでに花びんも作ってみた.

youtu.be

次はキャリッジパーツを作っておこう・・・

デルタ型3Dプリンタキット「Zonestar D810」を組み立てた(そのいち)

前々からマイ3Dプリンタが欲しかった(自作したかった)のですが,最近やっとお金が貯まったので思い切ってキットを買っちゃいました.格安になるように自分で一から設計しようかとも思ったのですが,流石に本体からフィラメントまで諸々が揃っていながら2万円台ってのを見るとそっちのがトータルコスト安いよね,ってことで香港から直輸入しました.それでも本当にバラバラなところから作るので,ちゃんと自作感があります.中学英語が読み書きできるならオススメ.

実は今まで研究室の3Dプリンタ(Repricator2)を実効支配してたのでマイ3Dプリンタは不要だったのですが,最近そいつのホットエンドに樹脂がべったりこびりついて使えなくなっちゃいました.それも購入に踏み切ったひとつの理由です.

はじめてのAliExpress

買うなら絶対デルタ型!って思ってたので(だって超カッコいいし),デルタロボットタイプの3Dプリンタが格安で売っていることで有名なAliExpressでの輸入に踏み切りました.RepRapプロジェクトのKosselコピー品がズラリと並んでいて,しかもどれも安いのです.Kosselは廉価を目的としてるので当たり前ですけど.

私の購入条件

以下の条件で購入を検討しました.

  • 予算3万円未満(送料・関税込なので実質$250位が目安)
  • デルタロボット型(カッコいいだろう!?(ギャキィ))
  • ヒーテッドベッド付属(ABSを使いたい)
  • オートレベリング機能あり(デルタ型ならあるとめちゃくちゃ便利)
  • フィラメント・SDカード等備品一式が付属(楽だし)
  • 最小積層ピッチ0.1mm以下(最近のだったらこれくらい普通)
  • 液晶は4段もしくはグラフィック液晶(2段キャラクタは流石に見づらい)
  • ビルド範囲がなるべく広い(安いのはやたら小さかったりする)
  • 2色刷り可能タイプ(多分無理だけど可能であれば…)
  • 店の評価が高いこと(中国のショップなのでこれが重要)

ぶっちゃけこんな構成で3万未満とかできるのか?って思うような超我儘な構成ですが,なんとお眼鏡に適うものがちゃんとあるんですね.流石中国.しかも送料無料だし.

Buy Products Online from China Wholesalers at Aliexpress.com

2色刷りタイプもありましたが,そっちは予算オーバーだったので除外.(最初は$200以下じゃんやったー!と思ってたのですが,実はその値段で買えるのはフィラメントだけで,本体を買うにはもう$250必要だったというオチでした.)

電装系が頭の方にあるので安定性がちょっと気になりますが,黒ボディに赤パーツという厨二心をくすぐるマシンだったのでまあいいやって感じで購入ボタンポチーッ!

注文から到着まで

幸いにしてZonestarは対応が丁寧で,平易な英語でスムーズなやり取りができました.

実はヒーテッドベッドが付属するのはOption 3を選んだときだけだったのですが,注文時点でOption 3が選択できなかったのとOptionの意味がよくわかっていなかったため,Option 2を注文していました.そこで「Option 2買っちゃったけどヒーテッドベッドが欲しい.買える?」と先方に伝えた所,「いいよ!じゃあ追加の$31をここから買ってね!」と柔軟に対応してくれました.これはありがたい.

FedExで送ってもらいましたが,通関処理のところで「コマーシャルインボイスの記入が不十分です」とか言われて止まってしまいました.これはコマーシャルインボイス(送付状のようなもの)の記載が不十分のため通関できないということです.私の場合,宛先に自分の名字が書かれていなかったのでNGを食らったみたいです.FedExカスタマーサービスに電話口で修正することで事なきを得ました.

そして到着.注文から1週間で追跡番号の連絡,それから4営業日で届きました.海外からなのに意外と早い気がする.

開封

f:id:Mzawa2:20160717141909j:plain

キターーーーーーーーーーー!!来た来た来たーーーーーーーーーーーッ!

否が応でもテンション上がる瞬間ですよね!!

f:id:Mzawa2:20160717142246j:plain

これがキットの中身全部.中国輸入だと壊れてるとか足りないとかザラらしいので,注意深く見ていきたいと思います.(そのにへ続く?)

OpenCVの内部パラメータでOpenGLの透視投影行列を作成(そのに)

以前にもOpenCVの内部パラメータでOpenGLの透視投影行列(フラスタム行列)を作成するコードを掲載しましたが,あれが数学的に何をやっているのかちょこっと解説します.元はパワポですが,SlideShareにアップロードするまでもない量と内容なので画像でスライドを掲載します.

OpenCVOpenGLの座標系

f:id:Mzawa2:20160612200757p:plain

OpenCVのデフォルトの座標系は左図,OpenGLの座標系は右図の通りです.OpenGLの方では暗黙的に,主点(光軸中心)が原点でZ軸と光軸が完全に一致しています.ですので,OpenCVOpenGLを渡り歩くときには,この図に従って適切に座標変換する必要があります.OpenCVのカメラモデルでは画素の剪断変形はありませんが,画素自体のアスペクト比が考慮されているため,fx=fyとは限りません.

透視投影行列

f:id:Mzawa2:20160612201114p:plain

さて,透視投影行列とは,空間の3次元位置を正規化デバイス座標系([-1,1]の3次元空間)にうまく圧縮する行列です.透視投影では視体積が錐体の頭を切り取った形をしているので,このような複雑な行列で変換する必要があります.なお,Z軸は出力写像が入力の反比例の関係になるように定められています.

OpenCV内部パラメータ表示

f:id:Mzawa2:20160612201704p:plain

導出と考え方は画像の通りです.znのあたりは透視投影行列の導出方法と合わせて読むと理解できると思います.

以上より,OpenCVのパラメータからOpenGLの透視投影行列を導くことができました.

FlyCapture2 SDKがVS2015で動かなくなった件

生存報告も兼ねて.

もう半年以上放置していたのは,しゅーろんとかいうクソイベとDの新生活が思っていた以上に忙しかったためです.

PointGray社のFlyCapture2 SDKにバグがあったのでその報告と対策.

何が起こったのか

【以前】 Windows 8.1(64bit)環境で,VisualStudio2013で開発していたときには何も問題なくビルドができていた.

【今回】 同じ環境で,VisualStudio2015にアップグレードしたら以前のコードに大量のエラー. なんかFlyCapture2::Cameraとか重要なやつに赤い波線が…

クラス宣言の外側に無効な型指定子があります.

namespace"FlyCapture2"にメンバー"Camera"がありません.

ぶっちゃけありえな~い!

原因

FlyCapture2::Camera等の型指定子FLYCAPTURE2_APIの定義は FlyCapture2Platform.h内のマクロで決められているが,こいつがVS2015以降に仕事しなくなったのだ.

// 前略
#if defined(WIN32) || defined(WIN64)
// Windows 32-bit and 64-bit
// 中略
#elif defined(MAC_OSX)

// Mac OSX

#else
// Linux and all others
// 後略

元々,VisualStudioの定義済みマクロは_WIN32_WIN64など,_が前に付いているのが普通であり, 今までWIN32とかWIN64が使えていたほうが不思議だったらしい. FlyCapture2 SDKは,この謎のマクロに依存した書き方をしていたため,VS2015以降の仕様変更(?)で動かなくなったのだろう.

参考:

定義済みマクロ

visual studio - C++ MSVC14.0 での WIN32マクロについて - スタック・オーバーフロー

対策

FlyCapture2.hをインクルードする前にWIN32をdefineしてしまえばいい.たとえば次の様に.

#ifdef _WIN32
#define WIN32
#endif
#ifdef _WIN64
#define WIN64
#endif

#include <FlyCapture2.h>
#include <FlyCapture2GUI.h>

以上.

この件については公式サポートに問い合わせメールを出しておきました.

cv::waitKey()の処理時間を計測してみた

ニコ動っぽいタイトルになっちゃったけど嫌いな人は許して(笑)

なんか遅くね?

PointGray社のFlea3は最大120fpsまで出る素晴らしい高速カメラですが,OpenGLを使わずOpenCVだけと連携させたところ,どうしても64fps以上出ないという問題にぶち当たりました.これではこのカメラを使う意味がないのでどうにか80fps以上は出せるようにしたい.そこで処理時間を計測してみたところ,どうやらcv::waitKey()で相当時間をかけているようなので,実際に計測してみました.

計測用コード

画像取得の部分に使ってるのは自作クラスですが,これについては以前の記事を参照のこと.

#include "FlyCap2CVWrapper.h"

using namespace FlyCapture2;

int main(void)
{
    FlyCap2CVWrapper cam;
    int count = 0;
    double time100 = 0;
    int waitparam = 1;
    // capture loop
    char key = 0;
    while (key != 'q')
    {

        // Get the image
        cv::Mat image = cam.readImage();
        flip(image, image, 1);
        cv::imshow("image", image);

        double f = 1000.0 / cv::getTickFrequency();
        int64 time = cv::getTickCount();

        key = cv::waitKey(waitparam);

        if (count == 1000)
        {// TickCountの変化を[ms]単位で表示 1000回平均
            std::cout << "param = " << waitparam << ", time = " << time100 / (count+1) << " [ms]" << std::endl;
            time100 = 0;
            count = 0;
            waitparam++;
        }
        else
        {
            count++;
            time100 += (cv::getTickCount() - time)*f;
        }
    }

    return 0;
}

1000回ループしたときの平均時間を表示したらcv::waitKey()の引数を1ずつ上昇させていくコードです. 風の噂によるとcv::getTickCount()は1msくらいの精度はあるみたいなので使っています.本当かどうか良く分かりませんが...

ちなみに計測環境は,Windows8.1 64bit,Intel Xeon X5675,12GB RAMです.

結果

f:id:Mzawa2:20151228175054j:plain

これはひどい

なんと,cv::waitKey()はせいぜい15ms程度の精度しかないことが分かりました.

ちなみにカメラ入力を描画するまでの処理は1.2msでした.

原因

これはcv::waitKey()の実装を見ればわかります.sources\modules\highgui\src\window_w32.cppの1900行目にあります.

CV_IMPL int
cvWaitKey( int delay )
{
    int time0 = GetTickCount();

    for(;;)
    {
        CvWindow* window;
        MSG message;
        int is_processed = 0;

        if( (delay > 0 && abs((int)(GetTickCount() - time0)) >= delay) || hg_windows == 0 )
            return -1;

        if( delay <= 0 )
            GetMessage(&message, 0, 0, 0);
        else if( PeekMessage(&message, 0, 0, 0, PM_REMOVE) == FALSE )
        {
            Sleep(1);
            continue;
        }

        for( window = hg_windows; window != 0 && is_processed == 0; window = window->next )
        {
            if( window->hwnd == message.hwnd || window->frame == message.hwnd )
            {
                is_processed = 1;
                switch(message.message)
                {
                case WM_DESTROY:
                case WM_CHAR:
                    DispatchMessage(&message);
                    return (int)message.wParam;

                case WM_SYSKEYDOWN:
                    if( message.wParam == VK_F10 )
                    {
                        is_processed = 1;
                        return (int)(message.wParam << 16);
                    }
                    break;

                case WM_KEYDOWN:
                    TranslateMessage(&message);
                    if( (message.wParam >= VK_F1 && message.wParam <= VK_F24)       ||
                        message.wParam == VK_HOME   || message.wParam == VK_END     ||
                        message.wParam == VK_UP     || message.wParam == VK_DOWN    ||
                        message.wParam == VK_LEFT   || message.wParam == VK_RIGHT   ||
                        message.wParam == VK_INSERT || message.wParam == VK_DELETE  ||
                        message.wParam == VK_PRIOR  || message.wParam == VK_NEXT )
                    {
                        DispatchMessage(&message);
                        is_processed = 1;
                        return (int)(message.wParam << 16);
                    }

                    // Intercept Ctrl+C for copy to clipboard
                    if ('C' == message.wParam && (::GetKeyState(VK_CONTROL)>>15))
                        ::SendMessage(message.hwnd, WM_COPY, 0, 0);

                    // Intercept Ctrl+S for "save as" dialog
                    if ('S' == message.wParam && (::GetKeyState(VK_CONTROL)>>15))
                        showSaveDialog(window);

                default:
                    DispatchMessage(&message);
                    is_processed = 1;
                    break;
                }
            }
        }

        if( !is_processed )
        {
            TranslateMessage(&message);
            DispatchMessage(&message);
        }
    }
}

あの精度が低いことで有名なGetTickCount()を内部で使っているんですね.メディア処理プログラムでは御法度レベルなんですがこれ...Intelさんなんで未だにリアルタイム画像処理ライブラリでこんなの使ってんの?

ただし,cv::waitKey()は内部的にOpenGLでいうglutSwapBuffers()に相当する処理も兼ねてるらしいので,HighGUIを使う限りこいつからは逃れられません.もっとFPSを上げたい場合は画像処理だけOpenCVにやらせてあとは別の描画ライブラリを使うようにしましょう.

とはいえこれはwindows環境での話なので,別のOSならまた結果は違うかもしれませんね.Linuxだともっと早いかも.

ARToolKitとOpenCVの歪み補正の違い

本記事ではARToolKitPlusを対象としますが,その本家であるARToolKitでも通用する話ですので,タイトルはARToolKitにしました.

ARToolKitPlusでの歪み補正

ARToolKitPlusのレンズ歪みの式は次の通り.(ARToolKitの公式チュートリアルより引用)

f:id:Mzawa2:20151220221100p:plain f:id:Mzawa2:20151220220658p:plain

このように歪みベクトルとして,歪みの中心座標(x0, y0),歪み係数f,スケールファクターs(センサの物理的大きさの逆数に対応?)を持ちます.基本的に半径方向の歪みのみを考慮しています.

OpenCVでの歪み補正

一方OpenCVでは,レンズ歪みの式を次のように規定しています.(OpenCV2.2 C++リファレンスより引用)

f:id:Mzawa2:20151220221505p:plain

ARToolKitPlusの方とは式が全く異なる上に,OpenCVの方では半径方向だけでなく円周方向の歪みも考慮されています.ちなみに歪み係数ベクトルdist_coeffs[]の並びは(k1, k2, p1, p2, k3(, k4, k5, k6))です.k4以降はデフォルトでは考慮されません.

ここで言いたいのは,「OpenCVdist_coeffs[]ARToolKitdist_factor[]にそのまま入れると破綻する」という事です.OpenCVで得たレンズ歪みパラメータはARToolKit(Plus)上ではそのまま使えないのです.(ただし内部パラメータは使える)

解決策

ではどうするか.簡単です.ARToolKitPlusに渡す前にOpenCV上で歪み補正してやればいいのです.つまり,

  1. OpenCVの内部パラメータをもとに,歪みベクトルが全て0のARToolKitPlus用のパラメータを用意する.
  2. ARToolKitPlus初期化の際に外部ファイルを使用せず,上記を渡す.
  3. ARToolKitPlusの歪み補正を無効化する.
  4. 画像を得る.
  5. OpenCVの機能で歪み補正する.
  6. ARToolKitPlusに渡す

という方針でやります.以下が必要な部分だけ取り出したコードです.

class OpenCVCamera : ARToolKitPlus::Camera
{
public:
    /**
   * Takes the OpenCV camera matrix and distortion coefficients, and generates
   * ARToolKitPlus compatible Camera.
   */
    static Camera * fromOpenCV(const cv::Mat& cameraMatrix, const cv::Mat& distCoeffs, cv::Size size)
    {
        Camera *cam = new ARToolKitPlus::CameraImpl;

        cam->xsize = size.width;
        cam->ysize = size.height;

        for (int i = 0; i < 3; i++)
        for (int j = 0; j < 4; j++)
            cam->mat[i][j] = 0;

        float fx = (float)cameraMatrix.at<double>(0, 0);
        float fy = (float)cameraMatrix.at<double>(1, 1);
        float cx = (float)cameraMatrix.at<double>(0, 2);
        float cy = (float)cameraMatrix.at<double>(1, 2);

        cam->mat[0][0] = fx;
        cam->mat[1][1] = fy;
        cam->mat[0][2] = cx;
        cam->mat[1][2] = cy;
        cam->mat[2][2] = 1.0;

        // OpenCVの歪み補正とARToolKitの歪み補正は全く別の計算式を使っているため,単純に値が使えない
        // ここではARToolKitの歪みベクトルをゼロとし,OpenCV側で補正してやることにする
        for (int i = 0; i < 4; i++) {
            cam->dist_factor[i] = 0;
        }
        return cam;
    }
}

int main()
{
    // ARToolKitPlusの初期化
    ARToolKitPlus::Camera *param = OpenCVCamera::fromOpenCV(cameraMatrix, distCoeffs, cameraSize);

    ARToolKitPlus::Logger *logger = nullptr;
    tracker = new ARToolKitPlus::TrackerSingleMarkerImpl<6, 6, 6, 1, 10>(cameraSize.width, cameraSize.height);
    tracker->init(NULL, 0.1f, 5000.0f);   // ファイルは使用しない
    tracker->setCamera(param);
    tracker->activateAutoThreshold(true);
    tracker->setNumAutoThresholdRetries(5);
    tracker->setBorderWidth(0.125f);            // BCH boader width = 12.5%
    tracker->setPatternWidth(60.0f);            // marker physical width = 60.0mm
    tracker->setPixelFormat(ARToolKitPlus::PIXEL_FORMAT_BGR);        // With OpenCV
    tracker->setUndistortionMode(ARToolKitPlus::UNDIST_NONE);        // UndistortionはOpenCV側で行う
    tracker->setMarkerMode(ARToolKitPlus::MARKER_ID_BCH);
    tracker->setPoseEstimator(ARToolKitPlus::POSE_ESTIMATOR_RPP);

    // Undistort Map
    initUndistortRectifyMap(
        cameraMatrix, distCoeffs,
        Mat(), cameraMatrix, cameraSize, CV_32FC1,
        mapC1, mapC2);

    // メインループ
    while (1)                            
    {
        colorImg = flycap.readImage();     // 適当なカメラ画像読込関数
        Mat temp;
        remap(colorImg, temp, mapC1, mapC2, INTER_LINEAR);
        ARToolKitPlus::ARMarkerInfo *markers;
        int markerID = tracker->calc(temp.data, -1, true, &markers);
        float conf = (float)tracker->getConfidence();      // 信頼度
        if (markerID == 4)
        {
                    // 適当なマーカー位置検出処理
                }
                // 適当な描画処理
    }
}

以上,参考まで.

OpenCVで得たカメラ内部パラメータをOpenGLの射影行列に変換

なぜかネットに正しい情報が少ないのでメモ書き

これまでのあらすじ

ARアプリケーションはARToolKitを使うのが最も楽ですが,残念ながらARToolKitGLUTに大きく依存しており,GLFW+GLEW+GLMを使ってモダンOpenGLで開発している自分にとってこれは非常にありがたくないです.特にGLSLを使う関係上,どうしてもプロジェクション行列・モデルビュー行列を自分で導かなくてはならなくなりました.モデルビュー行列はまだいいとして,プロジェクション行列は今まで適当にしか学んでなかったため,gluPerspective()arglCameraFrustumRH()以外の方法を知らず,OpenCVのカメラキャリブレーションで得た内部パラメータをOpenGLでどう使ったものか途方に暮れていました.arglCameraFrustumRH()の中身もなんかよくわかんないし...

上手くいったコード

//  OpenCVカメラパラメータからOpenGL(GLM)プロジェクション行列を得る関数
void cameraFrustumRH(Mat camMat, Size camSz, glm::mat4 &projMat, double znear, double zfar)
{
    // Load camera parameters
    double fx = camMat.at<double>(0, 0);
    double fy = camMat.at<double>(1, 1);
    double cx = camMat.at<double>(0, 2);
    double cy = camMat.at<double>(1, 2);
    double w = camSz.width, h = camSz.height;

    // 参考:https://strawlab.org/2011/11/05/augmented-reality-with-OpenGL
    // With window_coords=="y_down", we have:
    // [2 * K00 / width,   -2 * K01 / width,   (width - 2 * K02 + 2 * x0) / width,     0]
    // [0,                 2 * K11 / height,   (-height + 2 * K12 + 2 * y0) / height,  0]
    // [0,                 0,                  (-zfar - znear) / (zfar - znear),       -2 * zfar*znear / (zfar - znear)]
    // [0,                 0,                  -1,                                     0]

    
    glm::mat4 projection(
        2.0 * fx / w,      0,                     0,                                     0,
        0,                 2.0 * fy / h,          0,                                     0,
        1.0 - 2.0 * cx / w,   - 1.0 + 2.0 * cy / h, -(zfar + znear) / (zfar - znear),       -1.0,
        0,                 0,                     -2.0 * zfar * znear / (zfar - znear),  0);
    projMat = projection;
}

いろんな実装例を試しましたが,これがうまくいっています(引用元がどこだったか忘れました...).数学的にどうなってるのかまだ良く分かってないのでこれから勉強します.projectionの初期化時に転置されているように見えますが,これはGLMライブラリ(ひいてはOpenGL)の特性によるもので,実際には転置されていません.

(追記:2015/12/04)引用元はMicrosoft RoomAlive Toolkitのソースコードです.元のコードはDirectXなので左手座標系になっていますが,OpenGLの右手座標系にするためにZ軸を全て-1倍したものが上記のソースコードになります.

(修正:2015/12/21) 上手くいってるように思ってましたが,よく見ると間違ってることに気付いたので修正.参考資料は→ https://strawlab.org/2011/11/05/augmented-reality-with-OpenGL

(追記:2016/06/16) 上の透視投影行列の導き方をアップしました.

OpenCVの内部パラメータでOpenGLの透視投影行列を作成(そのに) - 聞きかじりめも

使い方

ちなみに,上のソースは次のように使います.

使用APIOpenCV, OpenGL3.3(with GLSL), GLFW, GLEW, GLM, ARToolKitPlusです.

// プロジェクション行列
    glm::mat4 Projection;
    cameraFrustumRH(cameraMatrix, cameraSize, Projection, 0.1, 5000);
// カメラビュー行列
// 光軸方向がZ軸正方向を向くカメラで,Y軸負方向がカメラの上ベクトルとする
// プロジェクション行列はそうなるように作っている
    glm::mat4 View = glm::mat4(1.0)
        * glm::lookAt(
        glm::vec3(0, 0, 0), // カメラの原点
        glm::vec3(0, 0, 1), // 見ている点
        glm::vec3(0, -1, 0)  // カメラの上方向
        );
// モデル行列
// ARマーカーの位置姿勢は,ARToolKitのarTransMat(),ARToolKitPlusのTrackerImpl::getModelViewMatrix()などで手に入れたやつを使う
    glm::mat4 markerTransMat = glm::make_mat4(tracker->getModelViewMatrix()); // ARマーカー位置姿勢
    glm::mat4 Model = glm::mat4(1.0)
        * markerTransMat; 
// モデルビュー行列,プロジェクション行列
// これをGLSLシェーダ―に転送する
    glm::mat4 MV = View * Model;
    glm::mat4 MVP = Projection * MV;

なお,カメラビュー行列をY軸正方向を上とした場合,つまり

 glm::mat4 View = glm::mat4(1.0)
        * glm::lookAt(
        glm::vec3(0, 0, 0), // カメラの原点
        glm::vec3(0, 0, 1), // 見ている点
        glm::vec3(0, 1, 0)  // カメラの上方向
        );

の場合は,cameraFrustumRH()の中身は次のようになります.

//  OpenCVカメラパラメータからOpenGL(GLM)プロジェクション行列を得る関数
void cameraFrustumRH(Mat camMat, Size camSz, glm::mat4 &projMat, double znear, double zfar)
{
    // Load camera parameters
    double fx = camMat.at<double>(0, 0);
    double fy = camMat.at<double>(1, 1);
    double cx = camMat.at<double>(0, 2);
    double cy = camMat.at<double>(1, 2);
    double w = camSz.width, h = camSz.height;

    // 参考:https://strawlab.org/2011/11/05/augmented-reality-with-OpenGL
    // With window_coords=="y_down", we have:
    // [2 * K00 / width,   -2 * K01 / width,   (width - 2 * K02 + 2 * x0) / width,     0]
    // [0,                 2 * K11 / height,   (-height + 2 * K12 + 2 * y0) / height,  0]
    // [0,                 0,                  (-zfar - znear) / (zfar - znear),       -2 * zfar*znear / (zfar - znear)]
    // [0,                 0,                  -1,                                     0]

    
    glm::mat4 projection(
        -2.0 * fx / w,     0,                     0,                                     0,
        0,                 -2.0 * fy / h,         0,                                     0,
        1.0 - 2.0 * cx / w,   - 1.0 + 2.0 * cy / h, -(zfar + znear) / (zfar - znear),       -1.0,
        0,                 0,                     -2.0 * zfar * znear / (zfar - znear),  0);
    projMat = projection;
}

気が向いたら実行結果のスクショを取りたいと思います.(今すぐ用意できる実行結果は企業秘密的にちょっとまずいファイルなので...)

(2015/12/22追記)この行列を使ってARした結果がこんな感じです.

f:id:Mzawa2:20151222163754p:plain:w120 f:id:Mzawa2:20151222163941p:plain:w120 f:id:Mzawa2:20151222163801p:plain:w120 f:id:Mzawa2:20151222163812p:plain:w120

(2016/06/12追記)この行列の導出過程をアップしました.

OpenCVの内部パラメータでOpenGLの透視投影行列を作成(そのに) - 聞きかじりめも