光線伝達のお勉強

はじめに

光学の分類に光線光学があります。これは光の電磁波的特性を無視したもので、波長が小さいと仮定しています。こうすることで線形の式で光の進む向きを計算できます。

光線伝達マトリックス

光の入力をv=(y1, θ1)、出力をw=(y2, θ2)とするとw=Mvにより入出力の関係が記述できます。

平面境界

スネルの公式n1sinθ1=n2sinθ2を使います。

下の例だとn1=1.51、n2=1.0です。入力側の方が大きい場合、大きい側に近い方に光が曲がります。

Vector2d RefractionNative(double y1, double theta1)
{
  double n1 = 1.51;
  double n2 = 1.0;

  Vector2d v;
  v << y1, theta1;

  Matrix2d M;
  M << 1.0, 0, 0, n1/n2;

  Vector2d w = M * v;

  return w;
}

球面境界

薄肉レンズ

おなじみのレンズはこんな感じです。

おわり

光学は前から勉強したいと思っているんですが、数式むずかしいです(・∪・)
シミュレーターあれば良いんですけどね。
参考文献：基本 光工学１

ExcelからOpenCVで画像を開く

はじめに

Excelにはプログラミングのできる環境としてVBAが用意されています。ただあまり使いやすいとはいえないので、別の方法があるか調べてみました。
COMというインターフェースを使うとVisual Studioで作成したDLLを呼べるみたいなので試してみます。

実装

ここではC++/CLIのOpenCVで取得したカメラ画像をExcelで表示するプログラムを作ってみます。VBAから見えてほしいインターフェースは、画像の座標を入れるとRGBの値が返ってくるものです。これを実装すると以下のようになります。
参考
Extend your VBA code with C#, VB.Net or C++/CLI | Pragmateek

#pragma once
using namespace System::Runtime::InteropServices;

#include <opencv2/highgui/highgui.hpp>

[Guid("FCB93926-3C39-447E-81D5-A1FB109A6EF3")]
public interface class ICamera
{
  array<int>^ GetPixelColor(int r, int c);
};

[Guid("A1AFD711-99A3-42C2-B6CD-1A2B2B1B3CB8")]
[ClassInterface(ClassInterfaceType::None)]
public ref class Camera : ICamera
{
public:
  Camera();
  ~Camera();

  virtual array<int>^ GetPixelColor(int r, int c);
private:
  cv::VideoCapture *capture_;
  cv::Mat_<cv::Vec3b> *image_;
};

VBAから呼びたい関数はインターフェースクラスで宣言します。ここだと座標rとcを入れると輝度値の配列が返ってきます。
GUIDは他とかぶらない値でVisual Studioの 　ツール->GUIDの作成　からコピペできます。cv::VideoCaptureとcv::Matをポインタにしているのは、C++/CLIがアンマネージドはポインタじゃないと動かないからです。これ結構つらい制約なんですよね。

#include "Camera.h"

#include <iostream>

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

Camera::Camera()
{
  capture_ = new cv::VideoCapture(CV_CAP_DSHOW);

  image_ = new cv::Mat_<cv::Vec3b>();
  *capture_ >> *image_;

  // 小さめにしてドット絵ぽく
  cv::resize(*image_, *image_, cv::Size(640 / 10, 480 / 10), 0, 0, cv::INTER_NEAREST);
}

Camera::~Camera()
{
  delete capture_;
  delete image_;
}

array<int>^ Camera::GetPixelColor(int r, int c)
{
  array<int>^ pixel = gcnew array<int>(3);
  pixel[0] = (*image_)(r, c)[0];
  pixel[1] = (*image_)(r, c)[1];
  pixel[2] = (*image_)(r, c)[2];

  return pixel;
}

作ったDLLはregasm.exeで登録してVBAから見えるようにします。Visual StudioのTool Commandから管理者権限で
RegAsm /codebase /tlb path-to-dll
を実行します。そうするとVBAのツール->参照設定から登録したDLLが見えます。

Option Explicit

Sub Display()
    Dim camera As ExcelCamera2.camera
    Set camera = New ExcelCamera2.camera
    
    Dim dataSheet As Worksheet
    Set dataSheet = Sheets("Sheet1")
    
    Dim r As Long
    Dim c As Long
    For r = 1 To 480 / 10
        For c = 1 To 640 / 10
            Dim pixel As Variant
            pixel = camera.GetPixelColor(r - 1, c - 1)
            dataSheet.Cells(r, c).Interior.Color = RGB(pixel(2), pixel(1), pixel(0))
        Next
    Next
End Sub

F5キー押して実行するとこんな感じになります。めでたし。

（続き書きました：パワポのスライド上でプログラムを動かす - wildpieの日記）

射影変換のお勉強

はじめに

射影変換はある平面を別の平面に射影することができる変換です。斜めから見たものを、もし正面から見たらどうなるかを計算できます。

変換式

変換前の座標(x,y)を(x',y')に変換するための行列Hを求めることが目的です。
（参考：CiNii 論文 -  反復を要しない射影変換の高精度解法）

x'とHxは平行になるので、x' × Hx=0と書けます。

とすると、内積

に変形できます。

最小二乗法

以下の式のJを最小とするhを求めることが最小二乗法になります。

これは２次形式なので、行列Mの最小固有値に対する固有ベクトルがhとなります。
実行結果が以下の図です。４点マウスで選択して左の画像を変換すると下の画像になります。

void Homography::Process(List<Point^>^ list1, List<Point^>^ list2)
{
  if (list1->Count != list2->Count)
    throw gcnew Exception("Data count unmatched");

  const double f0 = 1.0;
  if (is_first)
  {
     h = CalcH(list1, list2, f0); // 最小二乗法でHを計算

     h.resize(3, 3);
     MatrixXd temp = h.transpose();
     h = temp;	
     is_first = false;
  }

  cv::Mat_<cv::Vec3b> camera_image = usbcamra_->GrabImage();

  cv::Mat_<double> cvh;	
  cv::eigen2cv(h, cvh);

  cv::Mat_<cv::Vec3b> image;
  cv::warpPerspective(camera_image, image, cvh, image.size());

  cv::imshow("Result", image);
}

MatrixXd Homography::CalcH(List<Point^>^ list1, List<Point^>^ list2, double f0)
{
  const int N = list1->Count;
  MatrixXd mls = MatrixXd::Zero(9, 9);
  for (int i = 0; i < N; i++)
  {
    const double x1 = list1[i]->X;
    const double y1 = list1[i]->Y;
    const double x2 = list2[i]->X;
    const double y2 = list2[i]->Y;

    VectorXd xi[3] = {VectorXd::Zero(9), VectorXd::Zero(9), VectorXd::Zero(9)};
    xi[0] << 0, 0, 0, -f0*x1, -f0*y1, -f0*f0, x1*y2, y1*y2, f0*y2;
    xi[1] << f0*x1, f0*y1, f0*f0, 0, 0, 0, -x1*x2, -y1*x2, -f0*x2;
    xi[2] << -x1*y2, -y1*y2, -f0*y2, x1*x2, y1*x2, f0*x2, 0, 0, 0; 

    for (int j = 0; j < 3; j++)
      mls += xi[j] * xi[j].transpose();
  }
  mls /= N;

  // 最小固有ベクトルを求めて、最小２乗法をとく
  int min_eigen_index = 0;
  EigenSolver<MatrixXd> es(mls);
  es.eigenvalues().real().minCoeff(&min_eigen_index);
  return es.eigenvectors().col(min_eigen_index).real();
}

FFTでウェーブレット変換する

はじめに

フーリエ変換のように周波数の特性を得ることができる連続ウェーブレット変換を試しました。

連続ウェーブレット変換

母関数によっていろいろな種類がありますが、一番よく目にするのはGabor のウェーブレットです。こちらのサイトによると、（ガボールウェーブレット変換(VC++)）

を評価するみたいです。場所とスケールを変えた関数との類似度を計算していて、場所に応じた周波数の特性を求めることができます。フーリエ変換ですと、この場所が周波数高くてここは低いなどの議論はできませんが、ウェーブレット変換ならオッケーです。

実装

前の記事（EigenでFFTする - wildpieの日記）でFFTを使えるようになったので、コンボリューションは周波数領域で行います。周波数領域だとかけ算になるやつです。
たとえばガウシアンフィルタだとこんな感じです（たぶん）。

void CalcGaussKernel(VectorXd& kernel)
{
    const double sigma = 1.0;
    for (int t = 0; t < kernel.size(); t++)
    {
        kernel[t] = exp(-t*t/(2.0 * sigma))  / (sqrt(2.0*M_PI) * sigma);
    }
}

void GaussianFilter(const VectorXd& data, VectorXd& result)
{
    VectorXd kernel(data.size());
    CalcGaussKernel(kernel);

    VectorXcd data_freq = fft.fwd(data);
    VectorXcd kernel_freq = fft.fwd(kernel);
    VectorXcd product = data_freq.array() * kernel_freq.array();
    result = fft.inv(product);
}

これのKernelをGaborに変えます。

void CalcGaborKernel(VectorXcd& kernel, double a)
{
    const double sigma = 1.0;
    for (int i = 0; i < kernel.size(); i++)
    {
        double t = i / a / 8000.0;
        double omega = 2.0 * M_PI;
        kernel[i] = exp(-t*t/(sigma * sigma)) / (2.0 * sqrt(M_PI) * sigma)
                          * exp(std::complex<double>(0.0, omega * t)); // 共役?
    }
}

VectorXcd data_freq = fft.fwd(data);
MatrixXcd map(16, data.size());
for (int m = 1; m <= map.rows(); m++)
{
    double f = 2.0 * m; // 周波数, 0はどうしよう
    double a = 1.0 / f;
    VectorXcd kernel(data.size());
    CalcGaborKernel(kernel, a);

    VectorXcd kernel_freq = fft.fwd(kernel);

    VectorXcd product = data_freq.array() * kernel_freq.array();
    VectorXcd result = fft.inv(product);
    map.row(m - 1) = result / sqrt(a);
}

mapを表示するGUIを作りました。マイク音声をリアルタイムに処理しています
（c#でマイク音声をFFTする - wildpieの日記）。
合っているかあやしいですが、周波数の高い場所を検出できています。