茂加部珈琲店

主にtech関連のメモ置き場です

有理Bスプライン曲線の非有理化

今回はNURBS曲線の小ネタです
非有理化と言っていますが、正しくは計算過程で一時的に次元を上げることができるよ、という話です
これを利用すれば、有理化を考慮していないアルゴリズム(Bスプライン用の計算など)に、NURBSを適用することができます

有理化されたn次Bスプライン曲線を、非有理(すべての重さが1)のn+1次Bスプライン曲線に変換する

n次の有理Bスプライン曲線は、簡単にn+1次の非有理曲線に変換することができます。   NURBS曲線の定義を思い出してみましょう

{ \displaystyle
C(u) = \frac{\sum_{i=0}^{k-1} N_{i,n}(u)w_iP_i } {\sum_{i=0}^{k-1} N_{i,n}(u)w_i}
}

注目するのは、{ w_iP_i}の部分です。例えば、{P_i}が三次元ベクトル {\begin{pmatrix}
 x_i \\
 y_i \\
 z_i
\end{pmatrix}} で、ウェイト{w_i}が設定されているとすると、
{w_iP_i=\begin{pmatrix} x_iw_i \\ y_iw_i \\ z_iw_i \end{pmatrix}}となります。
ここで、{P'_i = \begin{pmatrix} x_iw_i \\ y_iw_i \\ z_iw_i \\ w_i \end{pmatrix}, w'_i = 1} とすれば、{w'_iP'_i=\begin{pmatrix} x_iw_i \\ y_iw_i \\ z_iw_i \\ w_i \end{pmatrix}}です。 xyz成分が一致しているので、これでなんとかなりそうです。
問題は、重み{w}が変更されてしまうので、分母が変わってしまうことです。{w'_i}はすべて1ですので、正規化係数の分母は1になり消えます
よって、このときのNURBSは { \displaystyle
C(u) = \sum_{i=0}^{k-1} N_{i,n}(u)w'_iP'_i
= \begin{pmatrix}
\sum_{i=0}^{k-1} N_{i,n}(u)x_iw_i \\
\sum_{i=0}^{k-1} N_{i,n}(u)y_iw_i \\
\sum_{i=0}^{k-1} N_{i,n}(u)z_iw_i \\
\sum_{i=0}^{k-1} N_{i,n}(u)w_i
\end{pmatrix}
}

あら偶然、結果のw成分に、重みの変更で消えた分母が見えています。これで結果の全体を除算して、xyz成分だけを取り出すと、{\frac{\sum_{i=0}^{k-1} N_{i,n}(u)w_iP_i } {\sum_{i=0}^{k-1} N_{i,n}(u)w_i}}となり、正しい結果を得れることがわかります

まとめ

有理Bスプライン曲線を、Bスプライン曲線用のアルゴリズムに適用するときは、次のような手順で一時的にn+1次のBスプライン曲線に変換する

  • n+1次の新しいベクトルを用意し、先頭n次元にn次ベクトルの制御点座標をコピーする
  • コピーしたベクトルに、重さ{w}を掛ける
  • コピーしたベクトルのn+1次元要素に、重さ{w}をコピーする
  • (アルゴリズムを適用する)
  • 結果全体を、結果のn+1次元要素で割る
  • 結果の先頭n次元を取り出して終わり

多くの場合、重み付きn次元座標はn+1次元ベクトルとして保存されているので、コードにするとかなり単純になります。 ベクトルの計算やコピーは、最適化がきくかどうかなどで変わってくるので、一番早い方法を検証する必要がありそうです

type traits(型特性)クラスの利用

またまたc++勉強ネタ。

traitsクラス

traitsクラスは、templateを利用したテクニックの一つです。
標準ライブラリや、boostなどで多用される手法の一つでもあります。
今回は、traitsクラスの基本的な利用法を勉強したいと思います。

traitsクラス

type traits(型特性)は、型に関する特性を表現するために利用されるクラスです。
traitsには色々なバリエーションがありますが、まず例として、STLiterator_traitsクラステンプレートを見てみましょう.

traitsを利用した例 std::iterator_traits

std::iterator_traitsは、イテレータに関する型特性を表現するためのクラスで、以下のような構成をしています.

以下のコードはcpprefjp - C++日本語リファレンスから引用しました。

namespace std {
  template <class Iterator>
  struct iterator_traits {
    using difference_type   = typename Iterator::difference_type;
    using value_type        = typename Iterator::value_type;
    using pointer           = typename Iterator::pointer;
    using reference         = typename Iterator::reference;
    using iterator_category = typename Iterator::iterator_category;
  };

  // ポインタに対する特殊化
  template <class T>
  struct iterator_traits<T*> {
    using difference_type   = ptrdiff_t;
    using value_type        = T;
    using pointer           = T*;
    using reference         = T&;
    using iterator_category = random_access_iterator_tag;
  };

  template<class T>
  struct iterator_traits<const T*> {
    using difference_type   = ptrdiff_t;
    using value_type        = T;
    using pointer           = const T*;
    using reference         = const T&;
    using iterator_category = random_access_iterator_tag;
  };
}

コードを見ると、iterator_traitsは単にtypedefされた型情報を格納しているだけで、メンバを持っていません.
使用者は、適切にこのテンプレートを特殊化することによって、iterator用に作られたアルゴリズムを利用することができるようになります

traits活用: タグ・ディスパッチ

std::iterator_traitsには、タグ・ディスパッチのために利用する iterator_categoryという型を持っています
iterator_categoryには、タグ(tag)と呼ばれる空クラスを指定します
標準で用意されているクラスは以下

  • std::input_iterator_tag
  • std::output_iterator_tag
  • std::forward_iterator_tag
  • std::bidirectional_iterator_tag
  • std::random_access_iterator_tag

さて、このiterator_categoryをどう使うかというと、関数が引数の型によって選択できることを利用します

//こんな感じでディスパッチさせる
template <typename InputIterator>
void foo(InputIterator itr, std::input_iterator_tag){
}

template <typename OutputIterator>
void foo(OutputIterator itr, std::output_iterator_tag){
}

template <typename Iterator>
void foo(Iterator itr){
  foo(itr, 
      typename std::iterator_traits<Iterator>::iterator_category());
}

また、std::forward_iterator_tagなどは他のイテレータタグを継承しています
タグが変換可能であれば呼び出しを行えるようになるので、継承がうまく利用されていますね

以上、type traitsの基本でした。この他にもstatic関数と組み合わせたり、traitsには様々な活用法があるみたいです
いろんなライブラリをよんで勉強していきたいですね!

君はプライベートメンバにアクセスできるフレンズなんだね!

いまだc++を勉強中です。今回はちょっとマイナーな機能、friendについて

friend

friendを使えば、クラスのプライベートメンバにアクセスできる外部関数を指定できます。
privateを指定すると、基本的に外部からのアクセスは禁止されるのですが、friendを使えば特例を作ることができる機能です。
この機能を使わなければいけないことは少ないと思いますが、覚えておくと便利かもしれません。

ちなみに、friend指定はメンバ関数のように行いますが、アクセスが許可されるだけで、関数自体はメンバではありません。
つまりthisポインタは渡されませんので、参照やコピーを行ってオブジェクトを関数に渡す必要があります。
もちろん、staticメンバであれば、インスタンスは不要です。

使いみちを考えてみましたが、ポインタでリンクを構築する時とか?
トリッキーなので、一から設計をする際には、あまり選択したくない方法ではありますね

サンプルコードを以下に示します。

#include <iostream>
class A;
class B;
class C;

class A{
    private:
    friend bool connectAB( A& a, B& b); //フレンド指定は宣言のように行います
    B* ptrb = nullptr; 
};

class B{
    private:
    friend bool connectAB(A& a, B& b);
    A* ptra = nullptr;
};

class C{
    private:
    friend bool connectAB(A& a,B& b);
    static void foo(){
        std::cout << "すごーい!" << std::endl;
    }
};
bool connectAB(A& a,B& b){
  //privateにアクセスできます
    a.ptrb = &b;
    b.ptra = &a;
    C::foo(); //staticなら、インスタンスは不要
    return true;
}

int main(){
    A a;
    B b;
    connectAB(a,b); 
    return 0;

    //出力:
    //すごーい!
}

NURBS曲線の特徴

CGの分野では解像度に非依存な形状を保存するために、パラメトリック曲線を利用することがあります. 今回は、その中でも特に自由度の高いNURBS曲線の基本的な事項をおさらいします.

ベクタ形式ではBezier曲線が有名で、古くから使われていますが、Bezier曲線にはいくつかの欠点があります. その中でも代表的なものを纏めてみると

  • 次数が制御点数に影響される
  • 重み付けができない
  • 局所的な変更が難しい

これらの問題点を解決するために考案されたのが、BSpline曲線です. NURBS曲線(Non-Uniform Rational BSpline curve)とは、名前の通りBSpline曲線に便利な特性(重み付けと非一様ノット列)を追加した形のことを言います.
NURBS曲線はBezier曲線とBspline曲線の上位版にあたり、高度な表現力をもちながら扱いやすい表現です

NURBS曲線の定義

NURBS曲線は以下の式で表現されます.
{ \displaystyle
C(u) = \frac{\sum_{i=0}^{k-1} N_{i,n}(u)w_iP_i } {\sum_{i=0}^{k-1} N_{i,n}(u)w_i}
}
ただし、kは制御点の個数, nは次数です {N_{i,n}}(u)はBSpline基底関数と呼ばれ、以下のように定義されます

{n=0}のとき
 { \displaystyle
\begin{equation}N_{i,n}(u)= \left \{\begin{array}{l} 1 (u_i \leq u \lt u_{i+1} ) \\ 0 (otherwise) \end{array} \right.\end{equation}
}

{n\neq0}のとき

 { \displaystyle
N_{i,n}(u)= \frac{u-u_i}{u_{i+n}-u_i}N_{i,n-1}(u)+\frac{u_{i+n+1}-u}{u_{i+n+1}-u_{i+1}}N_{i+1,n-1}(u)
}

ここで、{n=0}のときに1を返す範囲は、

{  \displaystyle
(u_i \leq u \lt u_{i+1})
}

{  \displaystyle
(u_i \lt u \leq u_{i+1})
}

のどちらかを採用します.
区間になっているほうの端は範囲に含まれなくなるので、エラー処理を行う必要があることに注意してください. また、

{  \displaystyle
(u_i \leq u \leq u_{i+1})
}

を採用して、{n=1}で例外処理を行う方法もあります

以下に基底関数の素直な実装を載せておきます.(改良の余地はあります.エラー処理はお好みで追加してください)
今度はDeBoorのアルゴリズムとかをやりたいですね

/** @fun
 * @brief BSpline基底関数を計算します
 * @param t 独立変数t
 * @param i 制御点の番号i
 * @param n 次数n
 * @param knots ノットベクトル
 */
double NURBS::basis(double t, size_t i, size_t n,
                    const std::vector<double> &knots) {

  // N[i][0]
  if (n == 0)
    return (knots[i] <= t && t < knots[i + 1])
               ? 1
               : 0;

  double a1 = knots[i + n] - knots[i];
  double a2 = knots[i + n + 1] - knots[i + 1];

  double b1 = t - knots[i];
  double b2 = knots[i + n + 1] - t;

  // "ゼロで割られたらゼロ"という決まりです
  double r1 = (b1 == 0) ? 0 : basis(t, i, n - 1, knots) * b1 / a1;
  double r2 = (b2 == 0) ? 0 : basis(t, i + 1, n - 1, knots) * b2 / a2;

  return r1 + r2;
}


型情報の保存

あけましておめでとうございます.

突然ですが、基底クラスから元の型を復元したいと思うことはないでしょうか.C++ではそのような型情報の保存を簡単に行えるみたいです.
具体的には、仮想関数を利用して型情報を返す関数を定義します.

#include <iostream>
#include <typeinfo>

struct Base{
    //型情報を返す仮想関数
    virtual const std::type_info& type() const{
        return typeid(*this);
    }
};

struct ClassA : public Base{
    
};


int main(){
    //基底クラス
    Base base;
    //派生クラス
    ClassA classA;
    //基底クラスのポインタに格納された派生クラス
    Base* base_classA = &classA;

    std::cout << "type of base is " << base.type().name() << std::endl;
    std::cout << "type of classA is " << classA.type().name() << std::endl;
    std::cout << "type of base_classA is " << base_classA->type().name() << std::endl;

    return 0;
}

/* 実行結果:
type of base is 4Base
type of classA is 6ClassA
type of base_classA is 6ClassA
*/

このようにして、仮想関数を使うことで、基底クラスのポインタからも元の型名が取り出せます.
typeinfoが致すればキャストするようにすれば、比較的安全に型が復元できるはずですね.
このテクニックはboostライブラリのboost::anyでも使用されているようです.

virtualなデストラクタ

C++では、ポリモーフィズムを利用するために作成した基底クラスのデストラクタにはvirtualを付けることがあるようです.
virtualがないと、基底クラスのポインタで管理している場合は、基底クラスのデストラクタのみが呼ばれてしまうためです.
初心者は「そんなこと聞いてないよ」となりそうで、不親切だとは思いますが、様々な理由でこうなっているはずなので仕方ないでしょう.

以下のようなコードで確認してみましょう

#include <iostream>
#include <memory>

struct Base{
    Base(){
        std::cout << "  Base constructor" << std::endl;
    }
    ~Base(){
        std::cout << "  Base destructor" << std::endl;
    }
};

struct Child : public Base{
    Child(){
        std::cout << "  Child constructor" << std::endl;
    }
    ~Child(){
        std::cout << "  Child destructor" << std::endl;    
    }
};

int main(){
    {
    std::cout << "unique_ptr:" << std::endl;
    std::unique_ptr<Base> pBase = std::make_unique<Child>();
    }
    {
    std::cout << "shared_ptr" << std::endl;
    std::shared_ptr<Base> pBase = std::make_shared<Child>();
    }
    std::cout << "Raw pointer" << std::endl; 
    Base* pBaseRaw = new Child();
    delete(pBaseRaw);

    return 0;
} 

実行結果はこうなります

unique_ptr:
  Base constructor
  Child constructor
  Base destructor
shared_ptr
  Base constructor
  Child constructor
  Child destructor
  Base destructor
Raw pointer
  Base constructor
  Child constructor
  Base destructor

Rawポインタとunique_ptrでは、Childのデストラクタが呼ばれていないのが確認できます.
基底クラスのデストラクタにvirtualをつければ、全てChildのデストラクタを呼ぶようになります
また、shared_ptrはvirtualデストラクタがないにもかかわらず、Childのデストラクタが呼ばれているのも面白いですね.
これは、shared_ptrが作成時のクラスを利用してポインタを破棄してくれるからみたいです.
しかし、unique_ptrではやはりvirtualが必要となるので、あまり安心はできなさそうですね.

無名共用体、構造体とスマートポインタ

最近C++を勉強中の私です.
今日は無名構造体とスマートポインタを使ったコードが上手くいかなかったので、反省として残しておきます.

例えば、以下のようなコードを書いたとします.

template <typename T>
struct test{
    union{
        struct{std::unique_ptr<T> x,y,z;};
        //etc.
    };

    test(T _x, T _y,T _z){
        x = std::make_unique<T>(_x);
        y = std::make_unique<T>(_y);
        z = std::make_unique<T>(_z);
    }
  
};

こういった書き方は、ベクトルクラスを作ったりするときに便利ですよね
しかし、このような書き方は「行儀の良い」クラスにしかできないみたいです.
この構造体の初期化を行うコードをclangでコンパイルすると、こんなエラーが出ます

error: attempt to use a deleted function
    test<int> t = {1,2,3};
              ^
note: destructor of 'test<int>' is implicitly deleted because variant field '' has a non-trivial destructor
    struct{std::unique_ptr<T> x,y,z;};
    ^

どうやら、デストラクタを生成できないためスマートポインタは無名構造体には入れてはいけないようです.まだ良くわかってない気がしますが...
ところで、これをg++でコンパイルするとどうなるでしょうか?

segfaultで落ちます

ここでは何のエラーも表示されませんので、g++だけを使うと小一時間悩む羽目になります.
gdbはexceptionを検出してくれますが、ライブラリの奥深くに飛ばされて更に混乱するだけです.
エラーが良くわからない場合は別のコンパイラを試してみようね!という教訓でした. ..あときちんとC++を勉強しよう!

こういう構造をガンガン使ってそうなglm(OpenGL Mathematics)がどのように対処してるのか見ようと思いましたが、ソースコードを見て目が回ってしまいました...
勉強しながら読んでみたいと思います