行列演算の応用 - mist-team/mist GitHub Wiki

MISTを用いると、本来は手間のかかる行列演算処理を、変数の演算と同じように記述することができます。

本チュートリアルでは、MISTの行列演算処理を利用して簡単な画像回転処理を記述するサンプルを紹介します。 行列演算に関しては​Document:行列演算と基礎編の行列演算も参考になります。

ヘッダの準備

基本的な行列演算（四則演算）を行うには、以下のファイルをインクルードします。

#include <mist/matrix.h>

また、逆行列の計算や固有値計算を行う場合は以下のファイルもインクルードする必要があります。 上の matrix.h とセットでインクルードしておくと良いでしょう。

#include <mist/numeric.h>

なお、 numeric.h を用いる場合は、行列計算ライブラリCLAPACKを準備しておかなければコンパイルできません。 もし準備がまだできていなければCLAPACK導入の方法を参照し、 numeric.h が使えるようにしておきましょう。パスを通すのと、コンパイル時のオプションを追加するのも忘れずに・・・ 行列の表現

行列の基本的な定義法は次のようになります。

mist::matrix< double > r( 3 , 3 );   // 3行3列の行列rを作成
mist::matrix< double > p( 3 , 1 );   // 3行1列の行列pを作成

行列へのアクセスは( 行番号, 列番号 )のように記述します。

p( 1 , 0 ) = 1.0;   // 行列pの2行1列目に1.0を代入

行列の応用演算

それではここで、 numeric.h を用いた行列の応用演算をいくらか紹介しましょう。

・正方行列 r の行列式を計算

double det_r = mist::det( r );

・行列 r の逆行列を計算

mist::matrix< double > r_inv = mist::inverse( r );

・係数行列 a と定数項ベクトル b から連立方程式を解く

// example
mist::matrix< double > a( 2 , 2 );
mist::matrix< double > b( 2 , 1 );
a( 0 , 0 ) = 2; a( 0 , 1 ) = 1;     b( 0 , 0 ) = 1;
a( 1 , 0 ) = 1; a( 1 , 1 ) = -1;    b( 1 , 0 ) = 5;

mist::matrix< double > x = mist::solve( a , b );

サンプルプログラム

逆行列計算を用いたサンプルプログラムを紹介します。 画像の回転を行う関数の中で、逆変換を求めるのに逆行列計算が用いられています。

#include <iostream>
#include <cmath>
#include <mist/mist.h>
#include <mist/io/bmp.h>
#include <mist/matrix.h>
#include <mist/numeric.h>

#define PI  3.14159265359

void rotation( const mist::array2< unsigned char > &image1 , mist::array2< unsigned char > &image2 , double arg );
unsigned char getpixel( double x , double y , const mist::array2< unsigned char > &img );
double min( const double a , const double b );
double max( const double a , const double b );


int main( int argc, char *argv[] )
{

    mist::array2< unsigned char > img1;
    mist::array2< unsigned char > img2;

    if( !mist::read_bmp( img1 , "input.bmp" ) ) // 入力画像
    {
        return( 1 );
    }

    rotation( img1 , img2 , 30.0 ); // 画像を30度時計回りに回転

    mist::write_bmp( img2 , "output.bmp" );

    return 0;
}

// 角度を入力とし画像を回転させる
void rotation( const mist::array2< unsigned char > &image1 , mist::array2< unsigned char > &image2 , double arg )
{
        
    // 回転行列
    mist::matrix< double > r( 3 , 3 );
    r( 0 , 0 ) = std::cos( arg * PI / 180.0 );
    r( 0 , 1 ) = -std::sin( arg * PI / 180.0 );
    r( 0 , 2 ) = 0.0;
    r( 1 , 0 ) = std::sin( arg * PI / 180.0 );
    r( 1 , 1 ) = std::cos( arg * PI / 180.0 );
    r( 1 , 2 ) = 0.0;
    r( 2 , 0 ) = 0.0;
    r( 2 , 1 ) = 0.0;
    r( 2 , 2 ) = 1.0;

    // 画像の左上座標
    mist::matrix< double > lt;
    lt = mist::matrix< double >::_31( 0.0 , 0.0 , 1.0 );

    // 画像の右上座標
    mist::matrix< double > rt;
    rt = mist::matrix< double >::_31( image1.width() , 0.0 , 1.0 );

    // 画像の左下座標
    mist::matrix< double > lb;
    lb = mist::matrix< double >::_31( 0.0 , image1.height() , 1.0 );

    // 画像の右下座標
    mist::matrix< double > rb;
       rb = mist::matrix< double >::_31( image1.width() , image1.height() , 1.0 );

    // 各座標の回転後の座標を求める
    lt = r * lt;
    rt = r * rt;
    lb = r * lb;
    rb = r * rb;

    // 回転後の平行移動行列
    mist::matrix< double > t( 3 , 3 );
    t( 0 , 0 ) = 1.0;
    t( 0 , 1 ) = 0.0;
    t( 0 , 2 ) = -min( min( rb( 0 , 0 ) , lt( 0 , 0 ) ) , min( lb( 0 , 0 ) , rt( 0 , 0 ) ) );
    t( 1 , 0 ) = 0.0;
    t( 1 , 1 ) = 1.0;
    t( 1 , 2 ) = -min( min( rb( 1 , 0 ) , lt( 1 , 0 ) ) , min( lb( 1 , 0 ) , rt( 1 , 0 ) ) );
    t( 2 , 0 ) = 0.0;
    t( 2 , 1 ) = 0.0;
    t( 2 , 2 ) = 1.0;

    // 逆変換行列を求める
    const mist::matrix< double >  inv_tr = mist::inverse( t * r );

    //回転後の画像の縦横幅を求める
    const unsigned int w 
        = static_cast< unsigned int >( max( fabs( rb( 0 , 0 ) - lt( 0 , 0 ) ) , fabs( lb( 0 , 0 ) - rt( 0 , 0 ) ) ) + 0.5 );
    const unsigned int h 
        = static_cast< unsigned int >( max( fabs( rb( 1 , 0 ) - lt( 1 , 0 ) ) , fabs( lb( 1 , 0 ) - rt( 1 , 0 ) ) ) + 0.5 );

    image2.resize( w , h );
    mist::matrix< double > p;
    unsigned int  x, y;
    for( x = 0 ; x < w ; x++ )
    {
        for( y = 0 ; y < h ; y++ )
        {
            // 回転処理後の画像中のピクセルp
            p = mist::matrix< double >::_31( x , y , 1.0 );

            // 回転処理前のピクセル位置を求める
            p = inv_tr * p;
            image2( x , y ) = getpixel( p( 0 , 0 ) , p( 1 , 0 ) , image1 );
        }
    }

    return;
}

// 線形補間によるサブピクセルの輝度の取得
unsigned char getpixel( double x , double y , const mist::array2< unsigned char > &img )
{
    if( x >= 0 && x < img.width( ) - 1 && y >= 0 && y < img.height( ) - 1 )
    {
        const unsigned int uint_x = static_cast< unsigned int >( x );
        const unsigned int uint_y = static_cast< unsigned int >( y );

        return static_cast< unsigned char > ( ( 1.0 - ( x - uint_x ) ) * ( 1.0 - ( y - uint_y ) ) * img( uint_x    , uint_y     )
                                            + ( x - uint_x )           * ( 1.0 - ( y - uint_y ) ) * img( uint_x + 1, uint_y     )
                                            + ( 1.0 - ( x - uint_x ) ) * ( y - uint_y )           * img( uint_x    , uint_y + 1 )
                                            + ( x - uint_x )           * ( y - uint_y )           * img( uint_x + 1, uint_y + 1 )
                                            + 0.5 );
    }
    else
    {
        return 0;
    }
}

double max( const double a , const double b )
{
    return ( a > b ) ? a : b;
}

double min( const double a , const double b )
{
    return ( a < b ) ? a : b;
}
⚠️ **GitHub.com Fallback** ⚠️