点群の凸包を求める

1. 問題
   
   　与えられた2次元の点群を包み込む多角形を求める問題です。下の□が与えられた点群で、多角形がそれを覆う凸包です。
   
   
   
   
2. 手法
   
   
   1. GrahamScan
      
      ここでの手法（プログラムは）　セジビック著　「アルゴリズムC++」 を参考にしました。図も著書を引用させていただきます。ここでは、概略しか説明しません。詳細は、本を参照してください。最初に以下の点が与えられたとします。まず、一番下でかつ右にある基礎となる点を求めます。ここでは、B点になります。
      
      セジビック著　「アルゴリズムC++」より
      
      　B点から各点への角度を求め、この角度で点を整列します。
      　B,M,J,L,N,P,K,...
      となります。
      
      
      　　　セジビック著　「アルゴリズムC++」より
      
      BとMは凸包の点です。次にJを加えます。この時点では、J も凸包の点とします。しかし、次にLを加えますとBJLは凸でなくなります。そこで、Jを凸包から除外しLを加えます。ここまでで、BMLが4点の凸包になります。どうように、LNPK と進みます。ここで、P点が除外されます。
      
      
      　　　　セジビック著　「アルゴリズムC++」より
      
      つぎにFIとつづき、E点に達すると、IFK点が除外されます。このように、除外走査は、複数遡って起こる場合があります。
      
      
      
      
   2. 処理時間
      　この方法では最初に、角度による点群のソートが必要ですが、このソート時間がこの手法の時間依存性を定めます。ソートが終われば、凸包 の点を求める時間は点の数に比例します。点を進めていく段階で、後戻りが起こりますが、この戻りは各点に対し1回のみです。
      　したがって、GrahamScanはソートを除けば、点の数に比例する時間しか必要としません。
      
      
3. プログラム
   
   
   1. プログラム構成
      　VC6のダイアログベース型MFCで作成します。画面構成は、左にnext ボタンを配置するのみです。
      
      
      
      
   2. クラス構成
      
      　ダイアログを管理する　CClosureDlg クラス以外に、点集合を扱う　CPoints　クラス、線分を扱う　CLine クラス、点を扱う　CGpoint クラス　を作成します。
      
      void CClosureDlg::OnPaint() でダイアログへの表示を行います。、点集合クラスのインスタンスを作成し、ランダムに点を生成し、表示します（ps.Mark()）。ps.GrahamScan(ps.Max); で凸包を求めるグラハムスキャンを行い、凸包上の点の数をnumに取得します。num個の点を折れ線で接続し、表示します（ps.Show()）。
      
      

           CClientDC *pDC;
         　pDC = new CClientDC(this);
      
           CPoints ps;
           ps.Mark(pDC,ps.Max);
           int num=ps.GrahamScan(ps.Max);
            ps.Show(pDC,num);

      
      nextボタンを押すとイベント処理で、void CClosureDlg::OnButtonNext() を実行します。ここでは、
      　Invalidate(TRUE);
      で画面の再描画を行います。点集合を再生成するため、別の点集合に対する凸包処理を行います。
      
      
   3. クラスの詳細
      CGPointは MFC のCPointを継承します。したがって、座標（ｘ、ｙ）や 点の代入演算　= はCPointクラスのそれを利用します。
      
      class CGpoint : public CPoint
      
      CGpoint クラスでは、点位置を表示する関数を用意します。
      　CLineクラスでは、二つの点から線分を定義するコンストラクタと、線分を表示する関数　show()、線分の水平線に対する傾きを計算する Theta() から構成されます。Theta() は線分の傾きを比から計算します。詳細は参考文献を参照してください。
      
      

      float CLine::Theta()
      {
        int dx,dy,ax,ay;
          float t;
          dx=p2.x-p1.x;
          ax=abs(dx);
          dy=p2.y-p1.y;
          ay=abs(dy);
          t= (ax+ay ==0) ? 0: (float)dy/(ax+ay);
          if (dx<0) t=2-t;
          else if(dy<0) t=4+t;
      
          return t*90.0f;
      }

      CPointsクラスでは配列points[]で点群を記録します。点の最大数は29です。コンストラクタでランダムな点を生成します。ここでは、配列の 1..Max を利用します。0 番目の要素は　「センチネル」 として利用します。
      

      CPoints::CPoints()
      {
          //200*200 の範囲のランダムな点を1..Maxに生成
          //points[0]はsentinel、points[MAX]も利用
          Max=29;
          for(int i=1;i<=Max;i++){
          points[i].x=(int)((rand()/32000.0)*200)+20;
          points[i].y=(int)((rand()/32000.0)*200)+20;
          }
      }

      
      次の GrahamScan がこの手法の中核です。右下点（表示では右上の点になります）を min に求め、minを1番目の数と交換します（swap(1,min)）。
      次に、SortbyTheta(N);で、この1番の点との角度で各点をソートします。ソート順位回りながら、凹んだ点をバックして削除しながら、凸包の点を求めていきます。Mが現在までの凸包の点の数です。
      
      

      int CPoints::GrahamScan(int N)
      {
          //Find Base
          int i,min,M;
          //右下の点をp[1]にする
         for(min=1,i=2;i<=N;i++)
           if(points[i].y < points[min].y) min=i;
           for(i=1;i<=N;i++)
               if(points[i].y == points[min].y)
                   if(points[i].x > points[min].x) min =i;
              //TRACE("base %d %d\n",min,points[min].y);
          swap(1,min);//baseの点を1番にする
      
          //走査開始
          SortbyTheta(N);
          points[0]=points[N];//Sentinel
          for(M=3,i=4;i<=N;i++){
                //M > M-1＞ iの点は時計回りでなくてはいけない
                while(points[M].Ccw(points[M],points[M-1],points[i]) >= 0) M--;
                M++;
                swap(i,M);
           }
          return M;
      
      }

      
      最後に、SortbyTheta(N); の例を示します。簡単のため選択法でソートしています。これでは、点の数が多くなると、時間がかかります。ヒープなど高速な手法に交換してください。
      
      

      void CPoints::SortbyTheta(int N)
      {
         int i,j,min;
           float thj,thmin=360.0;
           CLine L;
        
           for(i=1;i<=N;i++){
               min=i;
               thmin=360.0;
               L.p1=points[1];
               for(j=i;j<=N;j++){
                   L.p2=points[j];
                   thj=L.Theta();
                   if(thj < thmin){
                       thmin=thj;
                       min=j;
                   }
               }
               swap(i,min);
               //TRACE("sort:%d %f\n",min,thmin);
           }
      }

      
      
   4. ダウンロード
      このプロジェクトをダウンロードできます。次の行をクリックして、*.lzhファイルを適当なフォルダに保存します。
      ダウンロード開始
      　このファイルを解凍すれば、VC6.0 のプロジェクトファイルになります。.NET　でも実行可能です。