2次元の等高線を作成する関数テンプレート

2次元の等高線を作成する関数テンプレートを作成したよ。よかったら使ってみてくれ。縦横の格子の各頂点に高さ、温度、密度などの「値」がセットされていたとして、あとは閾値さえ与えてやれば、

  • 等高線の輪郭線形状
  • 等高線の輪郭線形状で囲まれた図形領域のテセレーション情報

を返してくれるよ。Vec2は四則演算をサポートする普通の空間ベクトルクラスであれば多分何でも行けると思う。

  template< class CellMatrix, class Vec2 > 
  void create_2d_contour
    ( 
    CellMatrix& cmat, 
    const float& threshold,
    vector< pair< Vec2, Vec2 > >* edges,
    vector< Vec2 >* tesselation
    )
    // policy concepts: ( use for your own concept template.. )
    //class yourCellMatrix{
    //public:
    //	float operator()( const int i, const int j, Vec2* pos  ){ /* return val(i,j) and pos(i,j), here.. */ } 
    //	const int NumRow( void ){ /* return number of rows, here.. */ }
    //	const int NumCol( void ){ /* return number of columns, here.. */ }
    //};
  {
    typedef pair< Vec2, Vec2 > Vec2Pair;
    typedef vector< Vec2Pair > Vec2PairVec;
    typedef vector< Vec2 > Vec2Vec;
    // 本関数は、等高線の構成要素である各セル内に発生させられた、
    // LINEのセットを、edges の中に、登録します。edges の中に登録された
    // 全てのエッジ要素は、閾値よりも値の大きい閉じた領域をCCWで囲む
    // 向きに定義されます。
    class RenderContourEdge
    {
    public:
      RenderContourEdge( Vec2PairVec* l, Vec2Vec* ts, const float thre )
        :l_(l), ts_(ts), thre_(thre){}
      void AddTesselation( const Vec2& a, const Vec2& b, const Vec2& c ){
        ts_->push_back( a );
        ts_->push_back( b );
        ts_->push_back( c );
        return;
      }
      // d--c  3--2
      // |  |  |  |
      // a--b  0--1
      void RenderCorner( const int o, bool rev ){
        // こちらは出力の輪郭線のCW,CCWを正しく決めるための大事な処理。
        Vec2 a2 = this->med_p( 0, 3, o );
        Vec2 c2 = this->med_p( 3, 2, o );
        if( !rev ){
          l_->push_back( make_pair( a2, c2 ) );
          this->AddTesselation( c2, v_[ (3+o)%4 ].second, a2 );
        }else{
          l_->push_back( make_pair( c2, a2 ) );
          this->AddTesselation( v_[ (1+o)%4 ].second, a2, v_[ (0+o)%4 ].second );
          this->AddTesselation( v_[ (1+o)%4 ].second, c2, a2 );
          this->AddTesselation( v_[ (1+o)%4 ].second, v_[ (2+o)%4 ].second, c2 );
        }
      }
      void RenderDiagonal( bool rev ){ 
        Vec2 ab = this->med_p( 0, 1 );
        Vec2 bc = this->med_p( 1, 2 );
        Vec2 dc = this->med_p( 2, 3 );
        Vec2 ad = this->med_p( 3, 0 );
        if( !rev ){
          l_->push_back( make_pair( bc, ab ) );
          l_->push_back( make_pair( ad, dc ) );
          this->AddTesselation( bc, ab, v_[1].second );
          this->AddTesselation( ad, dc, v_[3].second );
        }else{
          l_->push_back( make_pair( ab, bc ) );
          l_->push_back( make_pair( dc, ad ) );
          this->AddTesselation( ab, ad, v_[0].second );
          this->AddTesselation( ad, ab, dc );
          this->AddTesselation( dc, ab, bc );
          this->AddTesselation( dc, bc, v_[2].second );
        }
        return;
      }
      void RenderRectangle( const bool rev ){
        if( !rev ){
          this->AddTesselation( v_[0].second, v_[1].second, v_[2].second );
          this->AddTesselation( v_[0].second, v_[2].second, v_[3].second );
        }
        return;
      }
      void RenderSplit( const int o, bool rev ){
        // こちらは出力の輪郭線のCW,CCWを正しく決めるための大事な処理。
        Vec2 ad = this->med_p( 0, 3, o );
        Vec2 bc = this->med_p( 1, 2, o );
        if( !rev ){
          l_->push_back( make_pair( ad, bc ) );
          this->AddTesselation( v_[ (2+o)%4 ].second, v_[ (3+o)%4 ].second, ad );
          this->AddTesselation( ad, bc, v_[ (2+o)%4 ].second );
        }else{
          l_->push_back( make_pair( bc, ad ) );
          this->AddTesselation( ad, v_[ (0+o)%4 ].second, bc );
          this->AddTesselation( bc, v_[ (0+o)%4 ].second, v_[ (1+o)%4 ].second );
        }		
      }
    private:
      Vec2 med_p( const int c0, const int d0, const int o=0 )
      {
        int c = ( c0 + o )%4;
        int d = ( d0 + o )%4;
        {
          // 微小な計算誤差ゆえにセル間の要素が繋がらなくなるのを防ぐ
          // ためには、辺の計算は必ず。各軸の値の小さい方がc,大きいほうが
          // dと、統一させた上で、補間計算をしなくてはならない。
          int e = c*10+d;
          if( e==23 || e==30 || e==10 || e==21 ) swap( c, d);
        }
        return
          v_[c].second +
          ( thre_-v_[c].first ) / ( v_[d].first - v_[c].first ) *
          ( v_[d].second - v_[c].second );
      }
    public:
      pair< float, Vec2 > v_[4];
    private:
      Vec2PairVec* l_;
      Vec2Vec*	ts_;
      float thre_;
    };
    RenderContourEdge rce( edges, tesselation, threshold );;
    for( int i=0; i<cmat.NumRow()-1; ++i ){
      for( int j=0; j<cmat.NumCol()-1; ++j ){
        // d--c  3--2
        // |  |  |  |
        // a--b  0--1
        const int di[] = { 0,0,1,1,0 };
        bool rvsd=0;
        unsigned char bits=0;//  bits: abcd
        {
          for( int k=0; k<4; ++k ){
            pair< float, Vec2 >& ipair = rce.v_[k];
            ipair.first = cmat( i+di[k], j+di[k+1], &ipair.second );
            bits <<=1;
            bits |= ( ipair.first >= threshold ); 
          }// k
          if( bits & 0x08 ){// ビット対象性を利用 
            bits = ~bits;
            bits &= 0x0f;
            rvsd=!0;
          }
        }
        // 2の4乗=16通りのうち、一番目のビットは常に0という
        // ビット対象性を考慮して8通りは面倒見ましょう。
        switch( bits ){
        case 0x0:	rce.RenderRectangle( !rvsd ); break;// --
        case 0x1:	rce.RenderCorner( 0, rvsd ); break;// ad-dc
        case 0x2:	rce.RenderCorner( 3, rvsd ); break;// bc-dc
        case 0x3:	rce.RenderSplit( 0, rvsd ); break;// ad-bc
        case 0x4:	rce.RenderCorner( 2, rvsd ); break;// ab-bc
        case 0x5:	rce.RenderDiagonal( rvsd ); break;// bc-ab, ad-dc
        case 0x6:	rce.RenderSplit( 3, rvsd ); break;// ab-dc
        case 0x7:	rce.RenderCorner( 1, !rvsd ); break;// ab-ad
        default: break;
        }// switch
      }// j
    }// i
    return;
  }

使い方はこんな感じだ。

  
  // CellMatrix ポリシーを実装(関数内クラスでOK)
  class myCellMatrix{
  public:
	myCellMatrix( myGrid* g ):g_(g){};
	float operator()( const int i, const int j, Vec3* pos ){ 
		g_->CalcPos( i, j, pos );
		return g_->Matf()->operator()( i, j ); 
	}
	const int NumRow( void ){ return g_->Matf()->NumRow(); }
	const int NumCol( void ){ return g_->Matf()->NumCol(); }
  private:
	myGrid* g_;
  };

 // 等高線計算実行
  {
    vector< pair< Vec2, Vec2 > > edges;
    Vec2Vec tess;
    pocGeom::create_2d_contour( myCellMatrix( this ), threshold, &edges, &tess );
  }
 // 等高線計算終了