円の高速描画についてのページ

Counter : 79280

とある経緯から、簡単に円を描画する必要が出てきまして、その時にいろいろといじってみたのでここにまとめます。
　最初にググるのはもう普通ですね。円については数学の世界ではかなり昔から考察されていたようで、実数計算を使わない高速な描画方法も「ブレゼンハム」「ミッチェナー」などの名称でいくつか存在するようです。個別の手法の詳細な説明については他の方にお任せするとして、今回使ったのは英語版Wikiの中にあった「Midpoint circle algorithm」の項目。最初は上半分にある普通のアルゴリズムで見ていたのですが、下半分にある「Optimization」のプログラムを見て驚愕。ほとんど計算してない！でも円が描ける！ってのに感動しちゃって、これまでにいろいろ試していたのを全て放棄しちゃって組み込んでみました(汗) そうして出来たプログラムソースが以下の通り。

typedef struct {
  int xc,yc,r,flag;		// 入力円データ
  int p[8],x[8],y[8];           // 描画する座標
  int vx,vy;                    // 相対ベクトル
  int e;                        // ワーク用
} sg_circle_data;

sg_circle_data *init_sg_circle_data(sg_circle_data *scd,
                                    int x,int y,int r,int flag)
{
  int i;
  if(scd == NULL) scd = xalloc(sg_circle_data,1);

  scd->xc = x;
  scd->yc = y;
  scd->r = r;
  scd->flag = flag;

  for(i=0;i<8;i++)
    {
      scd->p[i] = FALSE;
      scd->x[i] = 0;
      scd->y[i] = 0;
    }
  scd->vx = 0;
  scd->vy = r;

  scd->e = - r;
  return(scd);
}

int end_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return(FALSE);
  return(scd->vx <= scd->vy);
}

void step_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return;
  if(scd->flag)
    {
      scd->p[0] = TRUE;
      scd->x[0] = scd->xc+scd->vx+1; scd->y[0] = scd->yc+scd->vy+1;
      scd->p[1] = TRUE;
      scd->x[1] = scd->xc-scd->vx+0; scd->y[1] = scd->yc+scd->vy+1;
      scd->p[2] = TRUE;
      scd->x[2] = scd->xc+scd->vx+1; scd->y[2] = scd->yc-scd->vy+0;
      scd->p[3] = TRUE;
      scd->x[3] = scd->xc-scd->vx+0; scd->y[3] = scd->yc-scd->vy+0;
      scd->p[4] = (scd->vx != scd->vy);
      scd->x[4] = scd->xc+scd->vy+1; scd->y[4] = scd->yc+scd->vx+1;
      scd->p[5] = (scd->vx != scd->vy);
      scd->x[5] = scd->xc-scd->vy+0; scd->y[5] = scd->yc+scd->vx+1;
      scd->p[6] = (scd->vx != scd->vy);
      scd->x[6] = scd->xc+scd->vy+1; scd->y[6] = scd->yc-scd->vx+0;
      scd->p[7] = (scd->vx != scd->vy);
      scd->x[7] = scd->xc-scd->vy+0; scd->y[7] = scd->yc-scd->vx+0;
    }
  else
    {
      scd->p[0] = TRUE;
      scd->x[0] = scd->xc+scd->vx; scd->y[0] = scd->yc+scd->vy;
      scd->p[1] = (scd->vx != 0);
      scd->x[1] = scd->xc-scd->vx; scd->y[1] = scd->yc+scd->vy;
      scd->p[2] = (scd->vy != 0);
      scd->x[2] = scd->xc+scd->vx; scd->y[2] = scd->yc-scd->vy;
      scd->p[3] = (scd->vx != 0) && (scd->vy != 0);
      scd->x[3] = scd->xc-scd->vx; scd->y[3] = scd->yc-scd->vy;
      scd->p[4] = (scd->vx != scd->vy);
      scd->x[4] = scd->xc+scd->vy; scd->y[4] = scd->yc+scd->vx;
      scd->p[5] = (scd->vx != scd->vy) && (scd->vy != 0);
      scd->x[5] = scd->xc-scd->vy; scd->y[5] = scd->yc+scd->vx;
      scd->p[6] = (scd->vx != scd->vy) && (scd->vx != 0);
      scd->x[6] = scd->xc+scd->vy; scd->y[6] = scd->yc-scd->vx;
      scd->p[7] = (scd->vx != scd->vy) && (scd->vx != 0) && (scd->vy != 0);
      scd->x[7] = scd->xc-scd->vy; scd->y[7] = scd->yc-scd->vx;
    }
}

void next_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return;
  scd->e += scd->vx;
  (scd->vx)++;
  scd->e += scd->vx;
  if(scd->e >= 0)
    {
      scd->e -= scd->vy;
      (scd->vy)--;
      scd->e -= scd->vy;
    }
}

実際に円を描画する部分は以下のとおり。

// 円の描画関数
void draw_circle(int x,int y,int r,int flag)
{
  int i;
  sg_circle_data *scd;

  for(scd=init_sg_circle_data(NULL,x,y,r,flag);
      end_sg_circle_data(scd);next_sg_circle_data(scd))
    {
      step_sg_circle_data(scd);
      for(i=0;i<8;i++)
        {
          if(scd->p[i])
            {
              pset_pixel(scd->x[i],scd->y[i]);  // 実際に点を描画する関数にしてね
            }
        }
    }
  free(scd);
}

// 中身の詰まった円の描画関数
void draw_fill_circle(int x,int y,int r,int flag)
{
  int i,px;
  sg_circle_data *scd;

  for(scd=init_sg_circle_data(NULL,x,y,r,flag);
      end_sg_circle_data(scd);next_sg_circle_data(scd))
    {
      step_sg_circle_data(scd);
      for(i=0;i<8;i+=2)
        {
          for(px=scd->x[i+1];px<=scd->x[i];px++)
            {
              pset_pixel(px,scd->y[i]);  // 実際に点を描画する関数にしてね
            }
        }
    }
  free(scd);
}

ちょっと無駄な処理に見える部分もあるかもしれないですけど、おそらく本当に最速なんじゃないかと思えるくらいのアルゴリズムです。最近のCPUは整数演算より浮動小数演算の方が速いよ、って点はおいとくとして(汗)
ちなみに、flagは余計かもしれないけど、FALSEの時は普通に指定画素の中心を原点にした、直径が「奇数画素」の描画をするもの。 TRUEは指定画素の右下(正確には(x+0.5,y+0.5)の点)を原点にした、直径が「偶数画素」の描画をするもの。 TRUE時は現時点では単に1画素ずらして描画しているだけなので、正確な円で無いことは確定なんだけど、これを正確にする計算方法が分かんないのと、理論上1画素以上はズレていないハズなのとでとりあえず保留。もし正確な計算方法が分かればそれはそれで、 flagによって計算式を切り替えれば済む話だし。
で、とりあえずflagはFALSEで固定して、早速半径0～10の円をそれぞれさらっと描画してみる。 (分かりやすくするために2倍に拡大してますが、それでも見えない人はブラウザで拡大して見てください(汗))

一番左が半径0で1点のみの描画。次が半径1、半径2…と続いて一番右が半径10の円。ん！？半径4の時の描画がなんかおかしい。他の円は8近傍なんだけど、半径4の時だけ一部4近傍っぽい描画になっている。これは是非修正したいって思っていろいろ調べてみたところ、vx++した時に vyを変化させない(vy--しない)状態でvx=vyになった場合でも描画してしまい、結果的に4近傍で繋がってるように描画する点が問題になることが判明した。そこで必要最小限の変更で修正を試みた結果、 next_sg_circle_dataを以下のように変更。

void next_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return;
  scd->e += scd->vx;
  (scd->vx)++;
  scd->e += scd->vx;
  if(scd->e >= 0)
    {
      scd->e -= scd->vy;
      (scd->vy)--;
      scd->e -= scd->vy;
    }
  else
    {
      // vyそのままの場合、vx==vyになるのなら描画しない
      if(scd->vx == scd->vy)
        {
          (scd->vy)--;
        }
    }
}

ちなみに半径4と半径25の時にこの条件にハマるようである。それ以上の半径の時は発生しなさそうなんだけど、確証が無い。
で、修正したサブルーチンで改めて描画してみる。

これでOKかなぁと思いつつ、今度は中心点を動かさないで0,2,4,6…と、2ずつ半径を増やして描画してみる。ついでに1,3,5,7…も隣に描いてみる。

むむむ、なんかスキマが妙に目立つ・・・よくみると、上で修正した半径4の円と半径6の円の間が際立っている。半径4の円をある意味小さくしたから、その分半径6の円と離れすぎてしまったんだよね。これまたどうにかならないかいろいろいじってみたんだけど、どうやら判定に使っている変数eの値が-1の時にvy--を実行して、概念的に円を内側に寄せるとうまくいきそう。これを、またまた必要最小限の変更で済ませようと思い、最初から変数eを大きめにすることで対応することにして、 init_sg_circle_dataを以下のように変更。

sg_circle_data *init_sg_circle_data(sg_circle_data *scd,
                                    int x,int y,int r,int flag)
{
  int i;
  if(scd == NULL) scd = xalloc(sg_circle_data,1);

  scd->xc = x;
  scd->yc = y;
  scd->r = r;
  scd->flag = flag;

  for(i=0;i<8;i++)
    {
      scd->p[i] = FALSE;
      scd->x[i] = 0;
      scd->y[i] = 0;
    }
  scd->vx = 0;
  scd->vy = r;

  scd->e = - r;
  if(r == 6) (scd->e)++;  // ピンポイント修正
  return(scd);
}

ただ、このままでは半径3～半径6の円が全て八角形になってしまい、並んだ円を見てるとあまり面白くない。ここまで半径4と半径6をいじってたんだから、その間の半径5もいじっちゃえということで、さらに以下のように変更。

sg_circle_data *init_sg_circle_data(sg_circle_data *scd,
                                    int x,int y,int r,int flag)
{
  int i;
  if(scd == NULL) scd = xalloc(sg_circle_data,1);

  scd->xc = x;
  scd->yc = y;
  scd->r = r;
  scd->flag = flag;

  for(i=0;i<8;i++)
    {
      scd->p[i] = FALSE;
      scd->x[i] = 0;
      scd->y[i] = 0;
    }
  scd->vx = 0;
  scd->vy = r;

  scd->e = - r;
  if((r == 5)||(r == 6)) (scd->e)++;  // ピンポイント修正
  return(scd);
}

以上の修正を入れて改めて描画してみる。

うん、目立っていたスキマは消えてるし、半径3～半径6の並びもそんなに悪くない・・・・よね？ということで(今回の)最終版まとめ。

typedef struct {
  int xc,yc,r,flag;		// 入力円データ
  int p[8],x[8],y[8];		// 描画する座標
  int vx,vy;			// 相対ベクトル
  int e;			// ワーク用
} sg_circle_data;

sg_circle_data *init_sg_circle_data(sg_circle_data *scd,
				    int x,int y,int r,int flag)
{
  int i;
  if(scd == NULL) scd = xalloc(sg_circle_data,1);
  
  scd->xc = x;
  scd->yc = y;
  scd->r = r;
  scd->flag = flag;
  
  for(i=0;i<8;i++)
    {
      scd->p[i] = FALSE;
      scd->x[i] = 0;
      scd->y[i] = 0;
    }
  scd->vx = 0;
  scd->vy = r;
  
  scd->e = - r;
  if((r == 5)|| (r == 6)) (scd->e)++;  // ピンポイント修正
  return(scd);
}

int end_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return(FALSE);
  return(scd->vx <= scd->vy);
}

void step_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return;
  if(scd->flag)
    {
      scd->p[0] = TRUE;
      scd->x[0] = scd->xc+scd->vx+1; scd->y[0] = scd->yc+scd->vy+1;
      scd->p[1] = TRUE;
      scd->x[1] = scd->xc-scd->vx+0; scd->y[1] = scd->yc+scd->vy+1;
      scd->p[2] = TRUE;
      scd->x[2] = scd->xc+scd->vx+1; scd->y[2] = scd->yc-scd->vy+0;
      scd->p[3] = TRUE;
      scd->x[3] = scd->xc-scd->vx+0; scd->y[3] = scd->yc-scd->vy+0;
      scd->p[4] = (scd->vx != scd->vy);
      scd->x[4] = scd->xc+scd->vy+1; scd->y[4] = scd->yc+scd->vx+1;
      scd->p[5] = (scd->vx != scd->vy);
      scd->x[5] = scd->xc-scd->vy+0; scd->y[5] = scd->yc+scd->vx+1;
      scd->p[6] = (scd->vx != scd->vy);
      scd->x[6] = scd->xc+scd->vy+1; scd->y[6] = scd->yc-scd->vx+0;
      scd->p[7] = (scd->vx != scd->vy);
      scd->x[7] = scd->xc-scd->vy+0; scd->y[7] = scd->yc-scd->vx+0;
    }
  else
    {
      scd->p[0] = TRUE;
      scd->x[0] = scd->xc+scd->vx; scd->y[0] = scd->yc+scd->vy;
      scd->p[1] = (scd->vx != 0);
      scd->x[1] = scd->xc-scd->vx; scd->y[1] = scd->yc+scd->vy;
      scd->p[2] = (scd->vy != 0);
      scd->x[2] = scd->xc+scd->vx; scd->y[2] = scd->yc-scd->vy;
      scd->p[3] = (scd->vx != 0) && (scd->vy != 0);
      scd->x[3] = scd->xc-scd->vx; scd->y[3] = scd->yc-scd->vy;
      scd->p[4] = (scd->vx != scd->vy);
      scd->x[4] = scd->xc+scd->vy; scd->y[4] = scd->yc+scd->vx;
      scd->p[5] = (scd->vx != scd->vy) && (scd->vy != 0);
      scd->x[5] = scd->xc-scd->vy; scd->y[5] = scd->yc+scd->vx;
      scd->p[6] = (scd->vx != scd->vy) && (scd->vx != 0);
      scd->x[6] = scd->xc+scd->vy; scd->y[6] = scd->yc-scd->vx;
      scd->p[7] = (scd->vx != scd->vy) && (scd->vx != 0) && (scd->vy != 0);
      scd->x[7] = scd->xc-scd->vy; scd->y[7] = scd->yc-scd->vx;
    }
}

void next_sg_circle_data(sg_circle_data *scd)
{
  if(scd == NULL) return;
  scd->e += scd->vx;
  (scd->vx)++;
  scd->e += scd->vx;
  if(scd->e >= 0)
    {
      scd->e -= scd->vy;
      (scd->vy)--;
      scd->e -= scd->vy;
    }
  else
    {
      // vyそのままの場合、vx==vyになるのなら描画しない
      // scd-rが小さい時にのみ発生？r=4,25で起きるのは確認したが、
      // それ以上では絶対起こらないという確証が無い。
      if(scd->vx == scd->vy)
	{
	  (scd->vy)--;
	}
    }
}

// 円の描画関数
void draw_circle(x,y,r)
{
  int i;
  sg_circle_data *scd;
  
  for(scd=init_sg_circle_data(NULL,x,y,r,FALSE);
      end_sg_circle_data(scd);next_sg_circle_data(scd))
    {
      step_sg_circle_data(scd);
      for(i=0;i<8;i++)
	{
	  if(scd->p[i])
	    {
	      pset_pixel(scd->x[i],scd->y[i]);  // 実際に描画する関数にしてね
	    }
	}
    }
  free(scd);
}

// 中身の詰まった円の描画関数
void draw_fill_circle(int x,int y,int r)
{
  int i,px;
  sg_circle_data *scd;
  
  for(scd=init_sg_circle_data(NULL,x,y,r,FALSE);
      end_sg_circle_data(scd);next_sg_circle_data(scd))
    {
      step_sg_circle_data(scd);
      for(i=0;i<8;i+=2)
        {
          for(px=scd->x[i+1];px<=scd->x[i];px++)
            {
              pset_pixel(px,scd->y[i]);  // 実際に描画する関数にしてね
            }
        }
    }
  free(scd);
}

これで完璧、のハズ・・・
ただ、数学者ではないので、今回の修正が果たして妥当なのかどうか証明できません。また個人的には、(半径3は仕方ないとしても)半径4,半径6の時の描画が八角形になってしまう部分を正直どうにかしたいんだけど、そこを直そうとすると他の描画結果がどうにもこうにも美しくない。
でもまぁ、高速なのは間違いないところなので今回はこれで終了、ということで。

E-mail address : saibara@big.or.jp