RealSenseをOpenCVで使う
今回紹介するクラスによって，RealSenseでVR or ARを実現するために必要なRealSenseの機能を簡単に使えるようになります．具体的には，
- RGBカメラ・デプスカメラ画像の取得
- デプスマップをRGBカメラから見た画像（とその逆）の取得
- デプスカメラを原点とする座標系（=RealSenseのワールド座標系）でのXYZ頂点座標の取得
- カメラ内部・外部パラメータ（デプスカメラから見たRGBカメラの位置姿勢）を取得
- 自動露出，自動ホワイトバランスの有効化・無効化
が可能となります．RGBカメラとデプスカメラは空間的に位置がずれているので，ARに使うためには外部パラメータが不可欠です．また，色に関する処理をする場合は自動露出がかえって邪魔になることがあるので，こちらも切れるようにしました．データはすべてOpenCVのcv::Matもしくはstd::vector<cv::Point3d>で取得できます．なお，カメラ内部パラメータはOpenCVで使える形式にしています．
サンプル
こんな感じに使います．
#include "RealSenseCVWrapper.h" int main(void) { RealSenseCVWrapper rscv(640, 480); rscv.getCalibrationStatus(); std::cout << "calibration data: " << "\n camera matrix = \n" << rscv.cameraMatrix << "\n camera distortion params = \n" << rscv.distCoeffs << "\n transform matrix of RGB cam = \n" << rscv.transform << std::endl; while (1) { cv::Mat src; // RealSenseのストリームを開始 rscv.queryFrames(); //rscv.getColorBuffer(src); // デプス画像にマップしたカラー画像を取得 rscv.getMappedColorBuffer(src); // 頂点座標の勾配を推定して陰影表示 rscv.getXYZBuffer(); cv::Mat XYZImage(rscv.bufferSize, CV_32FC1); for (int i = 1; i < rscv.bufferSize.height-1; i++) { for (int j = 1; j < rscv.bufferSize.width-1; j++) { // 隣4隅への勾配ベクトルを算出 float intensity = 0.0f; cv::Point3f v0, v1, v2, v3, v4; v0 = rscv.xyzBuffer[i*rscv.bufferSize.width + j]; if (v0.z != 0) { v1 = rscv.xyzBuffer[(i - 1)*rscv.bufferSize.width + (j - 1)] - v0; v2 = rscv.xyzBuffer[(i - 1)*rscv.bufferSize.width + (j + 1)] - v0; v3 = rscv.xyzBuffer[(i + 1)*rscv.bufferSize.width + (j + 1)] - v0; v4 = rscv.xyzBuffer[(i + 1)*rscv.bufferSize.width + (j - 1)] - v0; // それぞれの勾配から法線ベクトルを4種類算出して平均 // v1 v2 // v0 // v4 v3 // 法線方向は手前側が正 cv::Point3f n0, n1, n2, n3, n4; n1 = v2.cross(v1); n2 = v3.cross(v2); n3 = v4.cross(v3); n4 = v1.cross(v4); n0 = (n1 + n2 + n3 + n4) / cv::norm(n1 + n2 + n3 + n4); cv::Point3f light(3.f, 3.f, 0.f); intensity = n0.ddot(light/cv::norm(light)); } else { intensity = 0.f; } XYZImage.at<float>(i, j) = intensity; } } // RealSenseのストリームを終了 rscv.releaseFrames(); cv::imshow("ts", src); cv::imshow("xyz", XYZImage); int c = cv::waitKey(10); if (c == 27) break; if (c == 'a') { rscv.useAutoAdjust(false); rscv.setExposure(-8); rscv.setWhiteBalance(6000); } if (c == 'A') { rscv.useAutoAdjust(true); } } return 0; }
RealSenseCVWrapperクラスを初期化するには引数を指定する必要があります．あんまりC++の知識がないのでよくない実装な気もしますが，とりあえず動いているのでよしとしましょう．また，RealSenseCVWrapper::getXxxxBuffer()を利用する前と後をqueryFrames()とreleaseFrames()で挟む必要があります．デストラクタで終了処理をしているので，終了時については特に気にする必要はないと思います．
結果は次の画像のようになります．デプスカメラから見たカラー画像も，コマンド一つで取得可能です．また，XYZ座標値を可視化するために，法線を算出して簡単な拡散反射モデルで表示しています．
ソースコード
今回のクラスのソースコードです．自己責任でご利用ください．
RealSenseCVWrapper.h
#pragma once #include <opencv2\opencv.hpp> #include <RealSense\SenseManager.h> #include <RealSense\SampleReader.h> #pragma comment(lib, "libpxc.lib") class RealSenseCVWrapper { public: // constructor // @param // width : width of picture size // height: height of picture size RealSenseCVWrapper(); RealSenseCVWrapper(int width, int height); ~RealSenseCVWrapper(); void safeRelease(); // query RealSense frames // you should execute queryStream() function before get buffers // and releaseFrames(); function after all by every frame // @return // true = frame is succcessfully queried bool queryFrames(); void releaseFrames(); // get camera buffers void getColorBuffer(); void getDepthBuffer(); void getMappedDepthBuffer(); void getMappedColorBuffer(); void getXYZBuffer(); void getColorBuffer(cv::Mat &color); void getDepthBuffer(cv::Mat &depth); void getMappedDepthBuffer(cv::Mat &mappedDepth); void getMappedColorBuffer(cv::Mat &mappedColor); void getXYZBuffer(std::vector<cv::Point3f> xyz); // adjust camera settings // @param // use_aa: flag to use auto exposure and auto white balance void useAutoAdjust(bool use_aa); // set RGB camera exposure // @param // value: exposure time = 2 ^ value [s] // if you want 30fps you should set under 33 ~ 2 ^ -5 [ms] void setExposure(int32_t value); // set RGB camera white balance // @param // value: white point = value [K] // (ex: D65 standard light source = about 6500[K]) void setWhiteBalance(int32_t value); // get calibration status // - camera matrix // - distortion parameters // - transform matrix (RGB camera origin to depth camera origin) void getCalibrationStatus(); // OpenCV image buffer cv::Size bufferSize; // buffer size cv::Mat colorBuffer; // RGB camera buffer (BGR, 8UC3) cv::Mat depthBuffer; // depth camera buffer (Gray, 32FC1) cv::Mat colorBufferMapped; // RGB camera buffer mapped to depth camera (BGR, 8UC3) cv::Mat depthBufferMapped; // depth camera buffer mapped to RGB camera (Gray, 32FC1) std::vector<cv::Point3f> xyzBuffer; // XYZ point cloud buffer from depth camera (XYZ) // RGB camera calibration data writtern by OpenCV camera model cv::Mat cameraMatrix; // inclueds fx,fy,cx,cy cv::Mat distCoeffs; // k1,k2,p1,p2,k3 cv::Mat transform; // 4x4 coordinate transformation from RGB camera origin to the world (=depth) system origin protected: Intel::RealSense::SenseManager *rsm; Intel::RealSense::Sample *sample; };
RealSenseCVWrapper.cpp
(修正1/23)getColorBuffer()の問題を修正
#include "RealSenseCVWrapper.h" using namespace Intel::RealSense; using namespace cv; using namespace std; RealSenseCVWrapper::RealSenseCVWrapper() { } RealSenseCVWrapper::RealSenseCVWrapper(int w = 640, int h = 480) { rsm = SenseManager::CreateInstance(); if (!rsm) { cout << "Unable to create SenseManager" << endl; } // RealSense settings rsm->EnableStream(Capture::STREAM_TYPE_COLOR, w, h); rsm->EnableStream(Capture::STREAM_TYPE_DEPTH, w, h); if (rsm->Init() < Status::STATUS_NO_ERROR) { cout << "Unable to initialize SenseManager" << endl; } bufferSize = Size(w, h); } RealSenseCVWrapper::~RealSenseCVWrapper() { safeRelease(); } void RealSenseCVWrapper::safeRelease() { if (rsm) { //rsm->Release(); rsm->Close(); } } bool RealSenseCVWrapper::queryFrames() { if (rsm->AcquireFrame(true) < Status::STATUS_NO_ERROR) { return false; } sample = rsm->QuerySample(); } void RealSenseCVWrapper::releaseFrames() { rsm->ReleaseFrame(); } void RealSenseCVWrapper::getColorBuffer() { // Acquire access to image data Image *img_c = sample->color; Image::ImageData data_c; img_c->AcquireAccess(Image::ACCESS_READ_WRITE, Image::PIXEL_FORMAT_BGR, &data_c); // create OpenCV Mat from Image::ImageInfo Image::ImageInfo cinfo = img_c->QueryInfo(); colorBuffer = Mat(cinfo.height, cinfo.width, CV_8UC3); // copy data colorBuffer.data = data_c.planes[0]; colorBuffer = colorBuffer.clone(); // release access img_c->ReleaseAccess(&data_c); //img_c->Release(); } void RealSenseCVWrapper::getDepthBuffer() { Image *img_d = sample->depth; Image::ImageData data_d; img_d->AcquireAccess(Image::ACCESS_READ_WRITE, Image::PIXEL_FORMAT_DEPTH_F32, &data_d); Image::ImageInfo dinfo = img_d->QueryInfo(); depthBuffer = Mat(dinfo.height, dinfo.width, CV_32FC1); depthBuffer.data = data_d.planes[0]; depthBuffer = depthBuffer.clone(); img_d->ReleaseAccess(&data_d); //img_d->Release(); } void RealSenseCVWrapper::getMappedDepthBuffer() { // create projection stream to acquire mapped depth image Projection *projection = rsm->QueryCaptureManager()->QueryDevice()->CreateProjection(); Image *depth_mapped = projection->CreateDepthImageMappedToColor(sample->depth, sample->color); // acquire access to depth data Image::ImageData ddata_mapped; depth_mapped->AcquireAccess(Image::ACCESS_READ, Image::PIXEL_FORMAT_DEPTH_F32, &ddata_mapped); // copy to cv::Mat Image::ImageInfo dinfo = depth_mapped->QueryInfo(); depthBufferMapped = Mat(dinfo.height, dinfo.width, CV_32FC1); depthBufferMapped.data = ddata_mapped.planes[0]; depthBufferMapped = depthBufferMapped.clone(); // release access depth_mapped->ReleaseAccess(&ddata_mapped); depth_mapped->Release(); projection->Release(); } void RealSenseCVWrapper::getMappedColorBuffer() { // create projection stream to acquire mapped depth image Projection *projection = rsm->QueryCaptureManager()->QueryDevice()->CreateProjection(); Image *color_mapped = projection->CreateColorImageMappedToDepth(sample->depth, sample->color); // acquire access to depth data Image::ImageData cdata_mapped; color_mapped->AcquireAccess(Image::ACCESS_READ, Image::PIXEL_FORMAT_BGR, &cdata_mapped); // copy to cv::Mat Image::ImageInfo cinfo = color_mapped->QueryInfo(); colorBufferMapped = Mat(cinfo.height, cinfo.width, CV_8UC3); colorBufferMapped.data = cdata_mapped.planes[0]; colorBufferMapped = colorBufferMapped.clone(); // release access color_mapped->ReleaseAccess(&cdata_mapped); color_mapped->Release(); projection->Release(); } void RealSenseCVWrapper::getXYZBuffer() { Projection *projection = rsm->QueryCaptureManager()->QueryDevice()->CreateProjection(); std::vector<Point3DF32> vertices; vertices.resize(bufferSize.width * bufferSize.height); projection->QueryVertices(sample->depth, &vertices[0]); xyzBuffer.clear(); for (int i = 0; i < bufferSize.width*bufferSize.height; i++) { Point3f p; p.x = vertices[i].x; p.y = vertices[i].y; p.z = vertices[i].z; xyzBuffer.push_back(p); } projection->Release(); } void RealSenseCVWrapper::getColorBuffer(cv::Mat & color) { getColorBuffer(); color = colorBuffer.clone(); } void RealSenseCVWrapper::getDepthBuffer(cv::Mat & depth) { getDepthBuffer(); depth = depthBuffer.clone(); } void RealSenseCVWrapper::getMappedDepthBuffer(cv::Mat & mappedDepth) { getMappedDepthBuffer(); mappedDepth = depthBufferMapped.clone(); } void RealSenseCVWrapper::getMappedColorBuffer(cv::Mat & mappedColor) { getMappedColorBuffer(); mappedColor = colorBufferMapped.clone(); } void RealSenseCVWrapper::getXYZBuffer(std::vector<cv::Point3f> xyz) { getXYZBuffer(); xyz = xyzBuffer; } void RealSenseCVWrapper::useAutoAdjust(bool use_aa) { Device *device = rsm->QueryCaptureManager()->QueryDevice(); // get current values if you use manual mode device->SetColorAutoExposure(use_aa); device->SetColorAutoWhiteBalance(use_aa); } void RealSenseCVWrapper::setExposure(int32_t value) { rsm->QueryCaptureManager()->QueryDevice()->SetColorExposure(value); } void RealSenseCVWrapper::setWhiteBalance(int32_t value) { rsm->QueryCaptureManager()->QueryDevice()->SetColorWhiteBalance(value); } void RealSenseCVWrapper::getCalibrationStatus() { // aquire calibration data of RGB camera stream Projection * projection = rsm->QueryCaptureManager()->QueryDevice()->CreateProjection(); Calibration::StreamCalibration calib; Calibration::StreamTransform trans; projection->QueryCalibration() ->QueryStreamProjectionParameters(StreamType::STREAM_TYPE_COLOR, &calib, &trans); // copy RGB camera calibration data to OpenCV cameraMatrix = Mat::eye(3, 3, CV_32FC1); cameraMatrix.at<float>(0, 0) = calib.focalLength.x; cameraMatrix.at<float>(1, 1) = calib.focalLength.y; cameraMatrix.at<float>(0, 2) = calib.principalPoint.x; cameraMatrix.at<float>(1, 2) = calib.principalPoint.y; distCoeffs = Mat::zeros(5, 1, CV_32FC1); distCoeffs.at<float>(0) = calib.radialDistortion[0]; distCoeffs.at<float>(1) = calib.radialDistortion[1]; distCoeffs.at<float>(2) = calib.tangentialDistortion[0]; distCoeffs.at<float>(3) = calib.tangentialDistortion[1]; distCoeffs.at<float>(4) = calib.radialDistortion[2]; transform = Mat::eye(4, 4, CV_32FC1); for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { transform.at<float>(i, j) = trans.rotation[i][j]; // todo: rotationの行と列の順を確認 } transform.at<float>(i, 3) = trans.translation[i]; } projection->Release(); }
cv::Matにおけるclone()とcopyTo()の挙動の違い
※本記事はOpenCV Advent Calender 2016向けに書きました。 http://qiita.com/advent-calendar/2016/opencv
OpenCVをちょっと知ってる人はこのタイトルで「何言ってんだこいつ」ってなったかも知れませんが，一見同じに見える（※後述）この2種類のdeep-copyメソッドは，実はある特定のコードで挙動が異なります．
その特定のコードとはROI処理です．StackOverflowでも「ROIへのコピーはcopyTo()を使えばできるよ！」と書いてありましたが，なぜ（自分が）普段から使うclone()を使わないのかについて言及する人がいなかったので，今回の検証に繋がりました．
ROIに画像を貼り付けられない!?
例えば画像処理で顔検出して，そこにアニメアイコンを貼り付けるような自動クソコラ生成器を作ろうとする場合，ROIを設定してそこに既存の画像を貼り付ける，なんてことをすると思います．しかし，それを素直にcv::Mat::clone()メソッドで実装しようとするとなぜか貼り付けられない場合があります．
以下のコードは，次のような画像を作成しようとしているコードです．
// ROI Rect roi(100, 100, 100, 100); Mat img_roi = frame_c(roi); // 元画像のROIを生成 // clone Mat img_roi_clone = Mat(100, 100, CV_8UC3, Scalar(0, 255, 0)); // 貼り付ける画像を生成 img_roi = img_roi_clone.clone(); // 貼り付ける画像をcloneしてからROIに貼り付け imshow("Color ROI clone", frame_c); // なぜか貼り付けできてない
しかし，実際にこれで出てくる画像は・・・
このように，なぜか貼り付けできません．通常，ROIへの任意の画像処理は元画像frame_cに即座に反映されるため，このようにimg_roiに画像データをコピーすればROIの部分だけ塗りつぶされるはずなのですが，そうはならないようです．
そこでcv::Mat::clone()ではなくcv::Mat::copyTo()を使うと，目的の画像のようになります．
// ROI Rect roi(100, 100, 100, 100); Mat img_roi = frame_c(roi); // 元画像のROIを生成 // copyto Mat img_roi_copy = Mat(100, 100, CV_8UC3, Scalar(0,255,0)); img_roi_copy.copyTo(img_roi); // 貼り付ける画像をROIにコピー imshow("Color ROI copy", frame_c); // こうすると貼り付けできる
つまり，copyTo()でならROIに画像をコピーできるけど，clone()では画像をROIにコピーできないという仕様になっていることが確認できます．
※ちなみに，浅いコピー(shallow copy)img_roi = img_roi_shallow;は当然ながらうまくいきません．
やったね！今度からcopyTo()を使えばすべてうまくいくんだね！
更なる罠
…って簡単にはいかないのが困りものです．実は，一度「うまくいかないコピー」でROIをコピーしようとしてしまうと，その後いくらcopyTo()で正しく書こうが一切反映されなくなります．
// ROI
Rect roi(100, 100, 100, 100);
Mat img_roi = frame_c(roi);
// clone
Mat img_roi_clone = Mat(100, 100, CV_8UC3, Scalar(0, 0, 255));
img_roi = img_roi_clone.clone();
imshow("Color ROI clone", frame_c); // これが反映されないのは知ってる
// copyto after clone
Mat img_roi_copy_clone = Mat(100, 100, CV_8UC3, Scalar(255, 0, 0));
img_roi_copy_clone.copyTo(img_roi);
imshow("Color ROI copy after clone", frame_c); // ナンデ！コピーできないナンデ！？
このソースコードによれば，clone()を使う方では赤い矩形が，copyTo()を使う方では先ほど描いた赤の矩形が青で塗りつぶされるため青い矩形がそれぞれ描かれるはずです．ただし我々は既にclone()が上手くいかないことを知っているので，clone()の時点ではframe_cに何も描かれず，copyTo()の方で青い矩形が初めて描かれるだろうと期待します．しかし，実際には
このように，どちらの画像にも全く矩形が描かれなくなります．勿論，clone()ではなく=を先に書いたとしてもうまく描かれなくなります．
clone()の仕様
clone()もcopyTo()も，行列cv::Matの「深いコピー(deep copy)」に分類されます．じゃあclone()とcopyTo()の何が違うんじゃい！って思ったプログラマの皆さんなら，OpenCVのソースコードを覗き込みたくなると思います．実際，cv::Mat::clone()がどんな実装なのか見てみましょう．
inline Mat Mat::clone() const { Mat m; copyTo(m); return m; }
......
...
同じやんけ！！！
どうやら内部的にcopyTo()を呼び出しているようです（リファレンスにもそう書いてありましたが，挙動が違うことを発見したので信じられなくなって確認しました）．しかしよく考えてみてください．clone()ではわざわざ自分自身をcopyTo()した画像を作って返しています．そういう意味で本当に自身のクローンを生成しているわけです．
原因の考察
ROIは通常の画像と同等に扱えるように見えるけれども，実際には元画像とデータを共有しているため，やはり普通のcv::Matとはちょっと違います．これらのことから，clone()では画像のクローンですべての情報が上書きされてしまい，元画像との参照関係が崩れてしまうのだろうと推察できます．浅いコピーの方でうまくいかない理由は，浅いコピーが実際にはデータポインタを上書きする（＝元の参照関係を外す）処理であるため，これもまた元画像への参照が失われてしまうことが原因と思われます．そう考えると，一旦関係が崩れてしまったROIにいくら正しい方法で上書きしても，元画像には影響しないことが容易に察せられます．
結論
今回は，ROI処理においてclone()，copyTo()およびoperator=(cv::Mat)の挙動に違いがあることを確かめました．結論としては，ROIが元画像とデータ共有している関係上，copyTo()以外を使うとその参照関係が崩れてしまい，それ以降の変更が反映されなくなるのが原因と推察できました．皆さんもROIを扱う時にはcopyTo()を使っていきましょう．
cv::Mat画像のサブピクセル画素値を取得
久々にOpenCVのTips
画像処理をしていると，テクセルサンプリングのように浮動小数点座標の画素値を取り出したい場合が時々出てきます．そういう時って普通，整数座標の画素値を使ってバイリニア補間するんですけど，なぜか OpenCVでは「特定の1点のサブピクセル画素を取得する」という処理を標準で用意していません（ちなみにcv::Mat::at<>()メソッドの座標値に浮動小数点数を入れてもコンパイルエラーは起きませんが，実際には暗黙的にintに変換された座標値からサンプルされます）．サブピクセル精度でコーナー座標を検出するcv::findCornerSubPix()みたいに画像からサブピクセル座標を割り出す処理はあるんですが，今回やりたいのはその逆なんです．
ただし，サブピクセル精度で画像を切り出す処理ならあり，cv::getRectSubPix()でできます．今回のTipsのアイディアは，これで得たいサブピクセル座標を画像中心として3x3画像で切り取り中央座標値をサンプルする，というものです．意外と気づかない手法だと思うので，書き残すことにしました．
// グレースケール画像mの要素をサブピクセル座標pから取得する // 通常のMat::at()メソッドでは四捨五入された座標から取得されてしまうのでcv::getRectSubPix()を利用する // 画像範囲外は境界線と同じとする // サポートする型はCV_32F or CV_8U double sampleSubPix(cv::Mat m, cv::Point2d p) { cv::Mat r; cv::getRectSubPix(m, cv::Size(3, 3), p, r); return (double)r.at<float>(1, 1); }
上のコードはグレースケール画像のみ対応ですが，同じ方法でカラー画像もできます．ただし，実はcv::getRectSubPix()はCV_8UもしくはCV_32Fにしか対応していない関数なので，要素がdoubleのCV_64F画像はサポートしていません．その場合は仕方ないので一旦castしましょう．
透視投影変換にガンマをかけて嘘パースを実現
はじめに
今回の手法で次のような絵が撮れます．
本記事ではフリーのMMD描画アプリを利用して嘘パースの実験をしていますが，シェーダーレベルでアプリを改造する上，原理のみ紹介するので，CGデザイナよりもむしろエフェクトを開発するプログラマ向けの記事となります．予め御了承下さいませ．
嘘パースについて
皆さん，嘘パース（オーバーパース，漫画パース）というのはご存知でしょうか？
例えば，サンライズアニメでよく見かける，例の武器を画面いっぱいに突き出して構える迫力ある構図（サンライズパース）をCGアニメで再現したいとします．しかし，実寸で作られたCGに同じポーズをさせて普通に透視投影変換でレンダリングすると，あんなに大きく武器が映りませんのでちょっと迫力に欠けてしまいます．なぜそんなことが起こるかというと，サンライズパースでは正確な透視投影よりも手前の物体がより大きく映るようにあえて強調した投影に従って描いているからです．このように，通常の透視投影変換ではありえないくらい近い物体が大きく写っているような投影を嘘パースと呼びます．アニメや漫画では迫力があった構図なのにCGにすると物足りない，という現象の原因は嘘パースだったりします．
現在，CGアニメで嘘パースを実現する方法は，その構図専用の歪んだCGモデルを用意する というのが一般的です．でも，それってかなり面倒くさいですよね．だって同じキャラクターなのに超カッコいい構図の度に新しいCGモデルを用意しなきゃいけないんですから．他には，自動で人物CGモデルを局所的に肥大化させる手法（嘘パースプラグイン），画角の異なるカメラを複数用意してシームレスに映像を合成して実現する手法（マンガパース）などが知られています．
ここでは，より原理的な部分に着目して嘘パースを実現する手法を考案しました．それが，透視投影変換にガンマをかける という発想です．ディスプレイガンマのようにZ軸方向にガンマをかけることで，人間の知覚に近い歪み方をとっても簡単に再現できるのではないか？と思って試してみました．これは新しい投影変換を開発する必要がありますので，ちょっと考え方を示していきます．
原理
通常の透視投影
デザイン系の人には遠近法と言った方がわかりやすいかもしれないですが，透視投影とはすべての光線がある唯一の視点（原点）に集中するものとして考案された投影法です．
透視投影をCGで実現するには，レンダリング時に透視投影行列という同次変換行列（下式）をかませます．これによって，透視投影の視錐体がクリッピング空間に圧縮されるわけです．
透視投影にガンマをかける
透視投影では，物体の見た目の大きさは視点からの距離に反比例するという性質があります．3DCGプログラマにとっては非常に当たり前のことに思えますが，それはなぜでしょう？理由は，透視投影では光線が原点に直線的に集まると仮定しているからです．もちろん，この仮定は物理的には間違ってません．
対して嘘パースでは，視点の物体が通常の透視投影よりも極端に大きくなるという性質があります．透視投影では直線 に従って投影面
での見た目の大きさを決定していましたが，私はここでガンマ変換をかけることにしました．つまり，人が奥行きのあるシーンを観察したときに，スティーブンスの冪乘則に従って奥行知覚が実際の奥行きよりも強調される と仮定しました（これが事実かどうかは私は知りません）．この時の投影は次のようになります．
このようにすると，同じ大きさの物体でも遠くにあれば投影面ではより小さく映り，逆に近い物体はより大きく映ることになります．具体的には，同じ物体の見た目の大きさはの逆数に比例します．
ガンマ透視投影変換
この考えに基づいて投影行列を再構成すると，次のようになります．ここではガンマ透視投影変換と呼ぶことにします．とすると元の透視投影行列と一致するはずです．プログラム的には，モデルビュー行列をかけた後に一旦Z軸を編集する必要があるのでちょっと面倒ですが，バーテックスシェーダーで難なく実装できるでしょう．
実験
今回はMMD on WebGLのソースコードを改変して，初音ミクの標準モデル（あにまさ式）で試してみました．この場を借りて作者様に感謝します．
結果
まずは通常の透視投影変換でレンダリングした結果がこちらです．
これとほぼ同じ構図，ポーズで でガンマ透視投影変換でレンダリングした結果が次のようになります．
狙い通り，前方に突き出した手がかなり強調される嘘パースが実現できました．
今回の目的とは異なりますが，逆に とガンマを1より小さくした場合どうなるかというと...
こうなります．奥の方の手が透視投影よりも大きく映ります．偶然かもしれませんが，平行投影に似ています．
注意点
上の結果はかなりうまくいってるようですが，実はガンマ透視投影では構図をうまく調節しないとすぐにひどい結果になります．
このように，画面の外側は激しく歪みやすいです．しかし改善策はあります．ガンマ透視投影行列の光軸中心（行列の3列目の上２個）をずらしてやることで，ある程度問題を回避できるようです．
上２つはほぼ同じ方向からほぼ同じポーズをレンダリングしたもので，上は光軸中心を中央に固定しており，下では光軸中心を顔付近に持ってきています．このように，光軸を中央からずらして調節することで，ほぼ違和感のない嘘パースが様々な構図で実現できるようです．
経験則ですが，初音ミクのような人体モデルの場合はガンマによって頭部が歪むと強く違和感を覚えるので，頭部だけは歪まないように構図を調節すると良いでしょう．
倒立振子制作日記（そのに）
PWMによる電流制御
前回の記事でも取り上げましたが，この記事によるとDCモータで出力制御系を構成するときは電圧入力ではなく電流入力にすると理論的扱いがしやすいとのことですので，私もこの記事と同様に電流入力系を構成しようと思います．（注：以下，超基本的な所から説明します）
ちょっと復習
DCモータにはT-I特性（トルクvs電流）とT-N特性（トルクvs回転数）というものがあり，トルクと電流は比例し，負荷トルクが増すと直線的に回転数が減少する，という特性があります（ちなみに電圧を変えるとT-N特性が平行移動します）．つまり電流入力さえできればその時の負荷トルクさえ分かれば回転数（速度）が制御できることになります．
で，どーすんの？
さて一方で，例えDAコンバータを積んでいたとしてもマイコンにできるのはGPIOの電圧出力を変調することのみです．更に言うとArduinoUNOにはDACもないのでPWM出力をアナログ出力とみなすしかありません．モータの自己インダクタンスとか慣性とか非線形摩擦とかを考慮して，Hブリッジへの入力電圧だけでモータ入力電流値を予測するのは困難です．
そこで出てくるのが電流フィードバックです．これはモータへの入力電流を何らかの方法で計測し，その測定値をPWM出力にフィードバックして目標電流値を実現する，というものです．今回使用したモータドライバL298Nにはセンサ用PIN（と称するHブリッジの下流）がありますので，抵抗を挟んでGNDに接続し，抵抗間の電圧をアナログ入力に読ませることにしました．以前の記事では電流計測ICを用意していましたが，使い方が悪いのかセンサ自体の性能かわかりませんがノイズでかすぎて使い物にならなかったので抵抗にしました．
64種類 各20個（合計1280個）金属膜抵抗キット 1/4W（0.25W） 許容差±1%
- 出版社/メーカー: Ecowsera
- メディア: エレクトロニクス
- この商品を含むブログを見る
この抵抗セットがめちゃくちゃ便利でした．そんなに大量に使う予定がないなら，アキバでバラ100本入り100円とかで買うより良いかも．
コーディングする前に
電流センサ抵抗の選定
電流センサ抵抗の値はどうすればいいのか．これは中学生でも知っているオームの法則を利用します．下に今回の電流計測回路の簡単な回路図を示します．
Vss
|
[M] Motor (実際はHブリッジで接続)
|
|---- V (Analog Input)
|
I >
↓ > R
>
|
|
---
///
マイコンのアナログ入力はハイインピーダンスであるため，Analog Inputへは殆ど電流が流れません．そのためモータ入力電流は抵抗電流 \(I\) で近似でき，オームの法則 \(V=RI\) を利用してモータ入力電流をアナログ入力から得ることができます．従って，Analog Inputの許容最大電圧とMotor入力として想定する最大入力電流が分かれば \(R\) の値を見積もれます．
ArduinoUno(互換機)のアナログ入力の許容最大電圧はマイコン電源電圧と同じ \(V_{cc}=5[V]\) であり，モータドライバL298Nの定格最大電流は \(I_{max} = 2[A]\) であるので，
\(R = V_{cc} / I_{max} = 5[V]/2[A] = 2.5[\Omega]\)
と求められます．しかし現実には電池の供給能力の関係上2[A]も出せませんし，消費電力を考えると抵抗は小さい方がモータ入力に割く電流が大きくとれ，制御に有利そうです．そこで抵抗値をもっと下げて \(R=1[\Omega]\) にしてみるとどうでしょう．すると，定格最大電流は変わらないので想定最大電圧は
\(V_{max} = RI_{max} = 1[\Omega] \times 2[A] = 2[V]\)
と，先程の5[V]よりだいぶ低くなってしまいました．マイコンのアナログ入力の最大値は変わらず5[V]ですので，これはセンサ入力の分解能が低下することを意味します．計算は省略しますが，逆に抵抗値を増大させてセンサの分解能を上げようとする場合，そもそもの制御量であるモータ入力電流の最大値が減少してしまいます（或いは，最悪の場合マイコンのアナログ入力回路を破壊する程の電圧を許容することになります）．ということでここはバランスよく決めなくてはなりません．今回は手元に2.2[Ω]の抵抗があったので，最初の計算に近いということもありこいつを使います．
アナログ入力から電流値への変換
電子工作やPCに疎い人は分からないかもしれないので一応説明しておくと，実はマイコンのアナログ入力から直接返ってくる値は物理的な電圧値（実数）ではなく，それをADコンバータ(ADC)が読み取ってADCの精度(bit数)に応じて予め決められた最大値との比を整数値で返す，という動作をします．例えば，5[V]駆動の8bitADCに2[V]を入力すると，
\(256 \times {2[V] \over 5[V]} = 102.4 \simeq 102\)
という値が返ってきます．従ってこの逆算式
\(V = V_{max} \times {A \over A_{max}}\)
を通して物理値である電圧が得られますね（AはADCの出力）．更にセンサ値（今回は電流値）に変換するには電圧とセンシングしたい値との関係式をデータシートから読み取って実装すればいい訳です．今回はその関係式とはオームの法則に他ならないので， \(I=V/R\) として先程選定した抵抗値をぶち込めば晴れて電流値を読み取ることができます．
フィードバック制御
今回は速度型PI制御器を構成しました．同じPID制御でも位置型と速度型があるようで，位置型は制御量uそのものをPIDにより求めるもの，速度型は制御量uの増分をPIDにより求めるものです．今回の場合は電圧で直接電流を制御するので，以前の電流入力を保持する構造が何もありません．従って位置型にすると目標値に達した瞬間に電圧が0になり電流が急低下するといった現象が起こると考えられるので，速度型を採用する方が良さそうです．電流もいわば速度だし．
というわけで，制御量の数式は以下のようになります．
\(
u_t = u_{t-1} + \Delta u, \)
\(\Delta u = K_pe_t + K_i\sum_t{e_t},\)
\(\space e_t = I_\theta -I_t
\)
フィードバックゲインは限界感度法っぽい感じで適当に決めました．
制御量のPWM出力への変換
制御量をどうやってPWM出力に変換するのかも悩みどころです．制御量uは電流から生成したので，duty比が最大になった時に想定した最大電流が流れると仮定して，次の変換式に掛けました．
\(duty \space cycle = 255 \times {u_t/I_{max}}\)
255というのはArduinoのPWM出力が8bitだからです．単位がごちゃごちゃだし，あとあと考えてみると変な仮定な気がしますが，電流を電圧にしたところでゲインが変わるだけだし，これで動いたのでとりあえず良しとしましょう．
Arduinoスケッチ
CurrentFeedback.h
#ifndef CURRENTFEEDBACK_H_ #define CURRENTFEEDBACK_H_ #include "Arduino.h" #define LOOP_MAX 20 // アナログ入力平均化ループの周期 #define ANALOG_MAX 1024 // アナログ入力は10bit #define DUTY_MAX 255 // デューティ比は8bit // 定数 // current_max = analog_voltage / current_resist = 2.27A const float analog_voltage = 5.0; // アナログ最大電圧[V](analog = 1024) const float current_resist = 2.2; // モータ電流抵抗[ohm] const float current_max = 2.27; // モータ印加電流最大値[A] class CurrentFeedback{ protected: // Pin int analog_pin; // 電流センサ int input1_pin; // 正回転 int input2_pin; // 逆回転 int enable_pin; // PWM出力 // 内部変数 float current_target; // モータ印加電流の目標値[A] float Kp_current; // 電流比例フィードバックゲイン float Ki_current; // 電流積分フィードバックゲイン float error_integral; // 誤差積分値 float pwm_duty; // PWM出力のデューティ比（-255.0 ~ 255.0） public: // コンストラクタ CurrentFeedback(){ } CurrentFeedback(int analog, int input1, int input2, int enable, float feedback_gain_p = 1.0, float feedback_gain_i = 1.0){ begin(analog, input1, input2, enable, feedback_gain_p, feedback_gain_i); } // デストラクタ ~CurrentFeedback(){ } // 初期化関数(PIN, フィードバックゲインの登録) void begin(int analog, int input1, int input2, int enable, float feedback_gain_p = 1.0, float feedback_gain_i = 1.0){ // 引数の保存 analog_pin = analog; input1_pin = input1; input2_pin = input2; enable_pin = enable; Kp_current = feedback_gain_p; Ki_current = feedback_gain_i; // 値の初期化 current_target = 0; pwm_duty = 0; error_integral = 0; } // 目標値の登録 void setTargetCurrent(float target_current){ current_target = target_current; } // 現在の測定値を読みだす float getCurrent(){ // ADCの平均出力値を算出 float analog_mean = 0; for(int i=0;i<LOOP_MAX;i++){ int a = analogRead(analog_pin); analog_mean += (float)a/LOOP_MAX; } // ADC出力値から電流値への変換 float current_sense = analog_mean / ANALOG_MAX * analog_voltage / current_resist; current_sense = (digitalRead(input1_pin) == HIGH)? current_sense : -current_sense; // 印加電流の符号を決定 return current_sense; } // 登録した目標値に向けてフィードバック制御 void output(){ // sensor float current_sense = getCurrent(); // モータ印加電流の実測値[A] 符号は回転方向 // controller // フィードバックコントローラ 速度型PI制御 float current_error = current_target - current_sense; // 誤差 error_integral += current_error; float control_input = Kp_current * current_error + Ki_current * error_integral; // PI制御 pwm_duty += control_input / current_max * DUTY_MAX; // current -> duty cycle(-256 ~ 255) if(pwm_duty > DUTY_MAX) pwm_duty = DUTY_MAX; else if(pwm_duty < -DUTY_MAX) pwm_duty = -DUTY_MAX; int pwm_dst = (int)abs(pwm_duty); // actuator // モーター回転方向 if(pwm_duty > 0) { digitalWrite(input1_pin, HIGH); digitalWrite(input2_pin, LOW); } else if(pwm_duty < 0) { digitalWrite(input1_pin, LOW); digitalWrite(input2_pin, HIGH); } else{ digitalWrite(input1_pin, LOW); digitalWrite(input2_pin, LOW); } analogWrite(enable_pin, pwm_dst); } // 緊急停止 void quickStop(){ current_target = 0; pwm_duty = 0; digitalWrite(input1_pin, HIGH); digitalWrite(input2_pin, HIGH); analogWrite(enable_pin, 255); } }; #endif // CURRENTFEEDBACK_H_
こんな感じでクラス化してみました．Arduino IDEでクラスを書くのは初めてだし，ArduinoでC++記法がどこまで許されてるのか分からなかったのでちょっと大変でしたが，なんとかうまくいきました．やっぱC言語よりクラスが書けるC++の方が慣れると便利ですね．
動作テスト用スケッチ
#include "CurrentFeedback.h" #define MOTOR1_INPUT1_PIN 4 // Input 1 #define MOTOR1_INPUT2_PIN 5 // Input 2 #define MOTOR1_PWM_PIN 6 // Enable 1 #define MOTOR1_CURRENT_PIN 1 // Analog 1 CurrentFeedback current1; float current_target = 0; void setup() { // put your setup code here, to run once: Serial.begin(19200); Serial.println("hello"); current1.begin(MOTOR1_CURRENT_PIN, MOTOR1_INPUT1_PIN, MOTOR1_INPUT2_PIN, MOTOR1_PWM_PIN, 1.2, 0.4); } void loop() { // put your main code here, to run repeatedly: // シリアル入力（目標値の決定） char c = Serial.read(); if(c == '1'){ if(current_target > -current_max){ current_target -= 0.05; } else{ current_target = -current_max; } Serial.println(current_target, DEC); } if(c == '2'){ if(current_target < current_max){ current_target += 0.05; } else { current_target = current_max; } Serial.println(current_target, DEC); } if(c == '3'){ current_target = 0; Serial.println(current_target, DEC); } if(c == '4'){ current_target = 0; current1.quickStop(); Serial.println("------------STOP------------"); } current1.setTargetCurrent(current_target); current1.output(); Serial.print(" "); Serial.println(current1.getCurrent(), DEC); }
シリアルモニタからモータの入力電流を指定してそのときの測定値を出力させるプログラムです．実行した結果は，ハンチングが激しいですが電流値は一応制御できてるっぽいです．やったぜ．
（2016/09/19修正）ADCとDACを書き間違えていたので修正
倒立振子制作日記（そのいち）
はじめに
急に作りたくなったので．
私は一応，大学で機械工学を修め，ロボコンで実践的なロボット製作にも携わったことがある人間ですが，今ではすっかり情報系という半端モンでございます．研究でバーチャル空間にどっぷり漬かっている間にそろそろ機械工学の知識が風化し始めているので，ここいらでちょっと制御工学の復習をしとこうと思ってロボットの製作を思い立ったというわけであります．丁度夏休みだし．
そこで今回は，制御対象として比較的わかりやすい（であろう）部類であり，困ったときの先達がたーっくさんいる倒立振子を課題として，制御コードや回路から作って完成させようというプロジェクトを立ち上げてみました．勉強が目的なので，今後の制作では先達の情報は大いに参考にしつつも，できるだけ正解は見ないように（パクらないように）したいと思います．
部品集め
主にaitendo，秋月電子，千石電商で買いあさりました．それにしてもaitendoはセンサ・IC類がバカ安だな．．．とはいえ自分で全部揃えると結構かかりました．なるべく安く抑えようと思ったんだけどなー
（注）下表の他に，適当なジャンパ線，ケーブル類，はんだ，ピンヘッダ等が必要です．
| 品名・型番 | 購入場所 | 価格（円） | 用途 |
|---|---|---|---|
| びんぼうでいいの | aitendo | 999 | 制御マイコン |
| DCジャック | aitendo | 30 | 制御マイコン |
| Arduino用ピンヘッダ | aitendo | 95 | 制御マイコン |
| ATD6050 3軸加速度/ジャイロモジュール | aitendo | 450 | センサ |
| ATD712-M 高感度電流センサモジュール | aitendo | 750 * 2 | センサ |
| Arduino用ユニバーサル基板 | aitendo | 200 | 基板実装 |
| 電池ボックス単3x4本 | 秋月電子 | 60 | 電源 |
| 9V電池スナップ-DCプラグ変換アダプタ | 千石電商 | 100 | 電源 |
| L298N フルブリッジドライバ | 秋月電子 | 350 | モータードライバ |
| ブレッドボード | 千石電商 | 270 | 回路テスト |
| MYU-004 50mmタイヤセット(3mm六角シャフト用) | 千石電商 | 270 | 機体 |
| タミヤ 70097 ツインモーターギヤーボックス | 千石電商 | 790 | 機体 |
| タミヤ 70157 ユニバーサルプレート2枚セット | 千石電商 | 570 | 機体 |
| 合計 | 5,684 |
ちょっと解説
- 「びんぼうでいいの」とは，aitendo独自のArduino互換機です．ブートローダーが書き込まれていないので，Arduinoを持っていない場合は初めにこれを買うのは止した方がいいと思います．後述しますが結構書き込むの面倒だったので・・・
- 今回の倒立振子は，現在の角度のみを（ジャイロ＋加速度で）センシングして，その情報から駆動輪の出力を決定するフィードバック制御モデルを構成する予定です．
- ブラシ付DCモータの制御は電流制御入力の方がよいらしい(参考)ので，電流センサをわざわざ購入しました．結構高い・・・実は電流は2~5Ωくらいの抵抗があればオームの法則V=RIで直接計測できるんですが，どうしても抵抗で損失が生じてしまうせいで最大出力が落ちてしまうのが嫌だったので，購入に踏み切りました．実験的に抵抗でもやってみましたが，損失を少なくしようとするとどうしても出力電圧が小さくなります．まあオペアンプ買って電流センサ回路自作しても良かった気がするけど本質じゃないし・・・
- ギヤボックスのギヤ比は58.2:1を選びました．203.7:1はちょっと遅すぎた．
ハードウェア制作
テスト用回路
あとでアップロードします．
困った所
びんぼうでいいのへのブートローダー書き込み
Arduino IDEにはブートローダー書き込み機能がありますが，これを利用するにはブートローダー書き込み済みArduinoなどの書き込み装置が必要です．殆どの場合Arduino UNOを利用すると思います（ネットのTipsでもUNOばかり）が，生憎手元にはArduino Leonardoしかありませんでした．Leonardoの場合はちょっと特殊な操作が必要で，こちらの記事が参考になります．繋ぎ方はこちらの記事の通りにするとよいです．上記事を参考にする場合は，Arduino IDE 1.6.9では全くうまくいかなかったので，旧バージョンの1.0.6を使うと上手くいきました．参考にしてください．
L298Nのつなぎ方
実はデータシートの等価回路図をよく見ると分かるのですが，電流センサ用ピンのPin1とPin15はGNDに繋がないと動きません！ なぜかというと，これらはセンサ用と書いてありますが，実はHブリッジ回路のGNDに接続すべき部分になっているためです（センサと書いてあったので私は内臓センサのアナログ出力Pinだと誤解していました・・・）．つまりこれらのピンは各モータの電流の通り道そのものであって，センサを使わないからと言って浮かしてしまうとモータに電流が流れなくなってしまいます．従って，使わない時はGNDに接続し，電流センサに接続するときにも一方をGNDに接続する ことを忘れないようにしましょう．（このことに気づくまでに1日かかってしまったことは内緒）
以上，とりあえず文章だけ挙げておきます．余裕があったら追記していきます．
GLSLでOpenCV歪み補正を実装
OpenCVの歪み補正関数はcv::initUndistortRectifyMap()で予め作成したマッピング画像を使ってcv::remap()で歪んだカメラ画像を処理する，というものですが，このcv::remap()関数はOpenGLにやらせることもできます．それもプログラマブルシェーダGLSLを使った並列演算です．つまり，このテクニックを使えばGLSL上で並列化画像処理を実装することができるというわけです．（ここでは画像表示部をOpenGL/GLFWに任せています）
いやOpenCV使えよ
なぜOpenCVをすでに使っているのにわざわざそんなことをするかって？OK，これには理由が2つある．
まず一つ目は，画像処理を並列化することで実行速度の更なる向上が期待できるからだ．OpenCVの関数は既に強力な最適実装がされているものが多いが，それでもまだまだ最適化されていない実装もあるし，それに自分で作った画素アクセスを含む画像処理をもっと高速化するには並列化がおススメの場合が多い．せっかくOpenGLを使って画像表示しているのなら，プログラマブルシェーダを使ってGPUレベルで画像処理させるのが最も効率的なはず．だからGLSLを使うのだ．
もう一つの理由は，このような特殊なGLSLの使い方をすることでGPGPUによる並列演算テクニックの入り口に立つことができるからだ．まあつまり，GPU並列演算の初歩を学ぶという勉強目的ってこと．リアルタイム処理が必要なAR技術では今後必要になるであろう技術なので，勉強がてらOpenCVの機能を肩代わりするのがどういうことなのかを学んでおくことには価値がある．それにまだGLSLのこと良く知らないし．
というわけで，今回はcv::remap()をOpenGLレイヤーで実装することに挑戦しました．
実装
ソースコード
#pragma once #include <OpenGLHeader.h> #include <Shader.h> #include <opencv2\opencv.hpp> class Remapper { private: GLFWwindow *imgWindow; GLuint vao; GLuint vbo; GLuint imageObj, mapxObj, mapyObj; Shader s; int vertices; // バーテックスシェーダ const char *vertexSource = "#version 330 core \n" "layout(location = 0) in vec4 pv;\n" "void main(void)\n" "{\n" " gl_Position = pv;\n" "}\n"; // フラグメントシェーダ const char *fragmentSource = "#version 330 core \n" "uniform sampler2DRect image;\n" "uniform sampler2DRect map_x;\n" "uniform sampler2DRect map_y;\n" //"uniform vec2 resolution;\n" "layout(location = 0) out vec4 fc;\n" "void main(void)\n" "{\n" " vec4 mapx = texture(map_x, gl_FragCoord.xy);\n" " vec4 mapy = texture(map_y, gl_FragCoord.xy);\n" " fc = texture(image, vec2(mapx.x, mapy.x));\n" "}\n"; public: Remapper() { } Remapper(GLFWwindow *window) { init(window); } void init(GLFWwindow *window) { int w, h; glfwMakeContextCurrent(window); glfwGetWindowSize(window, &w, &h); imgWindow = window; // 頂点配列オブジェクト glGenVertexArrays(1, &vao); glBindVertexArray(vao); // 頂点バッファオブジェクト glGenBuffers(1, &vbo); glBindBuffer(GL_ARRAY_BUFFER, vbo); // [-1, 1] の正方形 static const GLfloat position[][2] = { { -1.0f, -1.0f }, { 1.0f, -1.0f }, { 1.0f, 1.0f }, { -1.0f, 1.0f } }; vertices = sizeof(position) / sizeof(position[0]); glBufferData(GL_ARRAY_BUFFER, sizeof(position), position, GL_STATIC_DRAW); glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, 0); glEnableVertexAttribArray(0); // テクスチャ glGenTextures(1, &imageObj); // source image について設定 glBindTexture(GL_TEXTURE_RECTANGLE, imageObj); glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_RGB8, w, h, 0, GL_BGR, GL_UNSIGNED_BYTE, NULL); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MAG_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MIN_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER); // map_x について設定 glGenTextures(1, &mapxObj); glBindTexture(GL_TEXTURE_RECTANGLE, mapxObj); glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_R32F, w, h, 0, GL_RED, GL_FLOAT, NULL); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MAG_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MIN_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER); // map_y について設定 glGenTextures(1, &mapyObj); glBindTexture(GL_TEXTURE_RECTANGLE, mapyObj); glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_R32F, w, h, 0, GL_RED, GL_FLOAT, NULL); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MAG_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_MIN_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER); glTexParameteri(GL_TEXTURE_RECTANGLE, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER); // シェーダのロード s.initInlineGLSL(vertexSource, fragmentSource); } void remap(cv::Mat src, cv::Mat &dst, cv::Mat map_x, cv::Mat map_y) { glfwMakeContextCurrent(imgWindow); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); s.enable();// シェーダプログラムの使用開始 GLint res[2] = { src.cols, src.rows }; // srcをテクスチャ0に転送する glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_RECTANGLE, imageObj); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_BGR, GL_UNSIGNED_BYTE, src.data); // map_xをテクスチャ1に転送する glActiveTexture(GL_TEXTURE1); glBindTexture(GL_TEXTURE_RECTANGLE, mapxObj); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_RED, GL_FLOAT, map_x.data); // map_yをテクスチャ2に転送する glActiveTexture(GL_TEXTURE2); glBindTexture(GL_TEXTURE_RECTANGLE, mapyObj); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_RED, GL_FLOAT, map_y.data); glUniform1i(glGetUniformLocation(s.program, "image"), 0); // uniform変数"image"にテクスチャユニット0(GL_TEXTURE0)を指定 glUniform1i(glGetUniformLocation(s.program, "map_x"), 1); glUniform1i(glGetUniformLocation(s.program, "map_y"), 2); //glUniform2fv(glGetUniformLocation(s.program, "resolution"), 1, res); glBindVertexArray(vao);// 描画に使う頂点配列オブジェクトの指定 glDrawArrays(GL_TRIANGLE_FAN, 0, vertices);// 図形の描画 glBindVertexArray(0);// 頂点配列オブジェクトの指定解除 s.disable();// シェーダプログラムの使用終了 glBindTexture(GL_TEXTURE_RECTANGLE, 0); // 読み込み dst = cv::Mat(src.size(), CV_8UC3); glReadBuffer(GL_BACK); glReadPixels(0, 0, dst.cols, dst.rows, GL_BGR, GL_UNSIGNED_BYTE, dst.data); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); } ~Remapper() { } };
例によって以前の記事で挙げたShaderクラスを使っています．OpenGLHeader.hにはOpenGLの他にGLFWとGLEW，GLMが含まれています．使い方は前回のGLImageクラスとほぼ同じです．Remapper::remap()のパラメータは順に，ソース画像(8bit,3ch)，出力画像(8bit,3ch)，マッピング画像(32bit,1ch)2枚，という構成で，ほぼcv::remap()と同様に使えます．
解説
フラグメントシェーダは，次のように書いてます．
const char *fragmentSource = "#version 330 core \n" "uniform sampler2DRect image;\n" "uniform sampler2DRect map_x;\n" "uniform sampler2DRect map_y;\n" //"uniform vec2 resolution;\n" "layout(location = 0) out vec4 fc;\n" "void main(void)\n" "{\n" " vec4 mapx = texture(map_x, gl_FragCoord.xy);\n" " vec4 mapy = texture(map_y, gl_FragCoord.xy);\n" " fc = texture(image, vec2(mapx.x, mapy.x));\n" "}\n";
やってることは簡単で，map_xとmap_yの画素値には浮動小数点で変形先の座標値が格納されているので，texture()で現在座標gl_FragCoord.xyの値（=移動先の画素値）を呼び出し，新たに得たリマップ先の座標値の画素値vec2(mapx.x, mapy.x)を本来の画像テクスチャimageからサンプルする，という実装になっています．ここのテクスチャサンプラsampler2DRect imageはサンプリング時に自動で線形補間するようにRemapper::init()内で設定していますので，これだけですべての画素について滑らかに変形できます．
ちなみに，わざわざ\nで改行するのは，コンパイルエラーが起きた場合に行数を把握しやすくするためです．\nを省略することもできますが，その場合ソースコードを1行で書いたとみなされてしまい，吐き出されるコンパイルエラーも行数が分からなくなってしまいます．
Remapper::remap()では，一度レンダリングしてからSwapせずに，直接Swap前のバッファにアクセスして画像データを読み込み，cv::Matに書き出しています．その後，バッファを全クリアして終了です．つまり，この一連の流れでcv::Mat型からcv::Mat型に変換する画像処理をシェーダレベルで実行したことになります．
なお，OpenGLとOpenCVでは画像テクスチャのY軸が逆ですが，ここでは反転処理を書いていません．なぜなら，今回の処理ではOpenGLに描かせた後にまたOpenCVに渡すので，Y軸が2回反転して結局元に戻るからです．もちろん，反転処理を入れればそれだけ遅くなります．GLSL上で反転させればそんなに速度変わらんかもしれんけど・・・
実装で苦労したところ
なぜかmap_xとmap_yテクスチャが転送できないという問題に陥りました．ほぼ前回のGLImageと同じだから，正しくソースコードが書けているはずなのに・・・
しかし，原因は（やっぱりって感じだけど）glTexImage2D()及びglTexSubImage2D()の設定にありました．本来ならば
glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_R32F, w, h, 0, GL_RED, GL_FLOAT, NULL); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_RED, GL_FLOAT, map_x.data);
と書くべきところを，
glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_R, w, h, 0, GL_R, GL_FLOAT, NULL); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_R, GL_FLOAT, map_x.data);
とか
glTexImage2D(GL_TEXTURE_RECTANGLE, 0, GL_LUMINANCE, w, h, 0, GL_LUMINANCE, GL_FLOAT, NULL); glTexSubImage2D(GL_TEXTURE_RECTANGLE, 0, 0, 0, res[0], res[1], GL_LUMINANCE, GL_FLOAT, map_x.data);
とか書いていたわけです．間違えどころはinternal format（正しいコードでGL_R32F）とformat（正しいコードでGL_RED）の値の組み合わせで，推奨される正しい組み合わせはOpenGLのリファレンスに書いてあります．GL_RとかGL_LUMINANCEとかはここには使えないくせに，こんなGL_REDとほぼ同じ意味のグローバル定数を別の値として用意しているOpenGLに俺は文句を言いたい．更に言うと，3chの画像でGL_RGBで2chの画像（某掲示板の事じゃないゾ！）でGL_RGって書くなら普通1chならGL_Rかなーとか思いこんじゃうでしょ…
くそが！！！！！
ふぅ，すっきりした．
結果
補正前
OpenCVで補正(cv::remap())
OpenGL(GLSL)で補正(Remapper::remap())
おお！ちゃんとできてる！これで並列化した高速画像処理演算ができるように…ん？
OpenCV: 0.014ms / 71.6fps
OpenGL: 0.016ms / 61.1fps
OpenCVの方が圧倒的に速いじゃんッッッ！！！
というわけで，まる2日かけて書いた今回の実装は完全に徒労に終わりました．でもきっと，今回の経験は将来の役に立つでしょう…
参考文献
床井研究室 - GLSL で画像処理 (1) 画像を取り込む
実は今回の記事は↑の丸パクリリファインです．床井先生いつも参考にさせてもらってます．ありがとうございます．
追記(2016/09/07)
この記事によるとどうやらglPixelRead()がめちゃ遅いらしく，PixelBufferObject(PBO)を活用するともっと高速化できそうです．PBOは4K映像も30FPSいけるくらい速いらしいし，今度試してみよう．