円形のグラデーションを描くアルゴリズム

円の縁が0で、円内部の指定した点Aに近づくに連れて値が上がり、点Aで1になるアルゴリズムを実装中なのですが、行き詰まったため質問です。

0を黒、1を白とすると上画像のような感じになるアルゴリズムを考えています(上画像は求めているアルゴリズム通りできていません)。
赤点が値が1の指定した点Aです。

点Aが円の中心と等しければ円内の各座標に対して中心からの距離を計算することで求めることができますが、
円の中心ではない箇所を1として、縁に近づくにつれ0に近くするというアルゴリズムがわかりません。

考え方だけでも構いませんのでお力添えいただけると幸いです。よろしくお願いいたします。

tiitoi

2021/01/12 12:01

「縁に近づくにつれ0にする」ならなぜ左側は縁からまだ距離があるのに途中で0になっているのでしょうか？

tiitoi

2021/01/12 12:05 編集

どのような経緯でそのような画像を作りたいのでしょうか？画像だけ見ると、3DCG でボールに右手前下の光源から光が当たっているような見え方ですね。 3DCGでやるライティングとかが参考になりそうかなと思いました。

mamyonya768

2021/01/12 12:17 編集

質問文の画像は求めているアルゴリズム通りできていません。文だけの説明では伝わりにくいと思い掲載しました。追記しておきます。値が段々と減るアルゴリズムが欲しいため実装しています。 3DCGのほうも調べようと思います。ありがとうございます。

行動規範の内容に同意します

回答5件

ベストアンサー

(x1, y1) のとき、t=1 で交点は t > 1 なので、

で [0, 255] の中間値を求めます。

直線と円の交点の求め方は高校数学の話なので解説は割愛しますが、以下に式の導出が載ってます。

Get location of vector/circle intersection? - Mathematics Stack Exchange

サンプルコード

Python で書きましたが、数式を愚直にコード化したので遅いです。
高速化したければ、C++ などで書き直してください。

import cv2
import numpy as np

# 真っ白の画像を作成
img = np.zeros((300, 300), np.uint8)

cx, cy = 150, 150  # 円の中心
r = 100  # 円の半径
x0, y0 = 200, 200  # 点 (質問画像赤の点)


def inside_circle(x, y):
    """円の内部かどうか
    """
    return np.sqrt((x - cx) ** 2 + (y - cy) ** 2) <= r


def get_intersections(x1, y1):
    """直線と円の交点を求める。
    """
    a = (x1 - x0) ** 2 + (y1 - y0) ** 2
    if a == 0:  # (x0, y0) の場合、t = 0
        return 0

    b = 2 * (x1 - x0) * (x0 - cx) + 2 * (y1 - y0) * (y0 - cy)
    c = (x0 - cx) ** 2 + (y0 - cy) ** 2 - r ** 2

    D = b ** 2 - 4 * a * c
    t1 = (-b + np.sqrt(D)) / (2 * a)
    t2 = (-b - np.sqrt(D)) / (2 * a)
    t = max(t1, t2)  # t > 0 の点が求める交点

    return t


# 素朴な実装
for y in range(img.shape[0]):
    for x in range(img.shape[1]):
        if inside_circle(x, y):
            # 円の内部の場合、t を計算する。
            t = get_intersections(x, y)
            # t に基づき、輝度値を決める。
            val = int((1 - 1 / t) * 255) if t > 0 else 255

            img[y, x] = val


cv2.imwrite("result.png", img)

追記

ppaul さんのコメントをうけて、他の補間方法での画像も作ってみました。
この関数を調整することで色合いがだいぶ変わってきますね
ガンマ補正に使われる関数を使って調整するといいかもしれません。

投稿2021/01/12 13:18

編集2021/01/13 16:05

tiitoi

総合スコア21956

mamyonya768

2021/01/12 13:49

迅速なご回答ありがとうございます。図もプログラムもわかりやすく、理解できました。ありがとうございました！

ppaul

2021/01/12 22:19

三角関数か二次関数をかぶせるともうすこし丸みを帯びた感じになると思います。また、numpyで配列演算を使うと、BLASという高速ライブラリを呼び出しますので、C/C++で自分で作るより速い場合が多いです。

tiitoi

2021/01/13 01:30 編集

> ppaul さんコメントありがとうございます。 numpy でループ除去して vectorization が効く書き方をすれば高速化できるのはわかるのですが、今回は言語の指定がとくになかったため、アルゴリズムをそのままコード化したものを載せました。 > 三角関数か二次関数をかぶせるこれは輝度値の値を計算するときに線形補間でなくすということでしょうか？

ppaul

2021/01/13 14:18

現在の描画は、xy平面上のx^2+y^2=1という円を底面とし、点A(p,q,1)を頂点とする斜円錐を考え、点(x,y,0)に対して斜円錐状の点(x,y,h)を取って、255×ｈで色(輝度)を付けているのですね。これだと、円錐なのでとがった感じがします。二次関数をかぶせるというのは、255×(-(h-1)^2+1)で色(輝度)を付けることです。三角関数をかぶせるというのは、255×sin(h×2/π)で色(輝度)を付けることです。すると、頂点の付近で微分が0になるので見た目も丸くなります。

tiitoi

2021/01/13 16:07 編集

なるほど。理解できました。試してみましたが、滑らかな見た目になりますね。補足ありがとうございます。

行動規範の内容に同意します

点Aと外周のみが指定されていて，「その中間がどのように変化していくか」は指定されていないから，こんなのでも良いのではないだろうか．

他の回答では，「画素位置と点Aとを結ぶ直線と円周との交点までの距離（計算がちょっとだけ面倒）」を使っているが，上述のように別にそうせねばならない理由も無いであろうから，「画素位置から円周までの距離（すなわち，画素位置と円の中心とを結ぶ直線と円周との交点までの距離．計算が簡単）」を使っても良いだろう，と．

C++
1const int IMG_SIZE = 240;   //画像サイズ(x,yで共用)
2const int CENTER = IMG_SIZE/2;  //てきとーに決めた円の中心(x,yで共用)
3const int RADIUS = 100; //円の半径
4
5int main()
6{
7    //点Aの座標 = (Ax,Ay)
8    const int Ax = CENTER + 45;
9    const int Ay = CENTER + 25;
10
11    //画像バッファ準備
12    cv::Mat Img( IMG_SIZE, IMG_SIZE, CV_8UC1 ); //IMG_SIZE*IMG_SIZEの画像バッファを用意
13    Img = 100;  //画像全体を背景の輝度値で埋めておく．値はてきとー．
14    
15    //全画素について走査し，その画素の輝度値を決めていく
16    for( int y=0; y<IMG_SIZE; ++y )
17    {
18        unsigned char *pI = Img.ptr<unsigned char>( y );
19        for( int x=0; x<IMG_SIZE; ++x, ++pI )
20        {
21            //(x,y)の外周からの距離
22            double DistFromCircum;
23            {
24                const int y_ = (y - CENTER);
25                const int x_ = (x - CENTER);
26                const double R = sqrt( double(x_*x_ + y_*y_) );
27                if( R > RADIUS )continue;
28                 DistFromCircum = RADIUS - R;
29            }
30            //(x,y)の点Aからの距離
31            double DistFromPointA;
32            {
33                int dx = x-Ax;
34                int dy = y-Ay;
35                DistFromPointA = sqrt( (double)( dx*dx + dy*dy ) );
36            }
37            //てきとーに２つの距離の比率で．
38            *pI = cvRound( 255 * DistFromCircum / (DistFromPointA + DistFromCircum) );
39        }
40    }
41    //結果画像の表示
42    cv::imshow( "Img", Img );
43    if( cv::waitKey() == 's' ){ cv::imwrite( "Result.png", Img );   }
44    return 0;
45}

投稿2021/01/13 02:44

編集2021/01/13 02:52

fana

総合スコア11708

mamyonya768

2021/01/13 17:16

変化方法に関してはあまり気にしていなかったのでご回答頂いた手法でも問題なさそうですね。自分では思いつかなかった考え方でした。ありがとうございました！

行動規範の内容に同意します

点Aが円の中心と一致する場合の，円の中心からの距離と明るさの関係のグラフを下図の左側だとしよう．
中央は明るさが1であり，円周では明るさが0である．

点Aの場所がずれた場合，単純に，下図の真ん中のグラフのようになるのだと考える．

これを，下図右側の絵のように，たくさんの円盤（緑横線）を位置をずらしながら重ねたような形だと考えれば，
各円盤の径と中心位置は，明るさに対して線形に変化する形だ．
話を簡単にするために中心が原点の単位円で考えれば，明るさをt(0～1)としたとき，明るさtの円盤は

半径(t) = 1-t
中心(t) = t * A, ここでAは点Aの座標．

である．

上記より，ある画素位置(x,y)に対する明るさを求めることは，その画素が無数の円盤のうちのどの円盤の周上にあるのか(すなわちtの値)を求めることだ．
中心(t)から(x,y)までの距離 = 半径(t) という式を t に関して解けばいい．

C++
1const int IMG_SIZE = 240;   //画像サイズ(x,yで共用)
2const int CENTER = IMG_SIZE/2;  //てきとーに決めた円の中心(x,yで共用)
3const int RADIUS = 100; //円の半径
4
5int main()
6{
7    //点Aの座標 = (Ax,Ay)
8    const int Ax = CENTER + 45;
9    const int Ay = CENTER + 25;
10
11    const double Ax_ = double(Ax-CENTER)/RADIUS;
12    const double Ay_ = double(Ay-CENTER)/RADIUS;
13
14    //画像バッファ準備
15    cv::Mat Img( IMG_SIZE, IMG_SIZE, CV_8UC1 ); //IMG_SIZE*IMG_SIZEの画像バッファを用意
16    Img = 100;  //画像全体を背景の輝度値で埋めておく．値はてきとー．
17    
18    //全画素について走査し，その画素の輝度値を決めていく
19    for( int y=0; y<IMG_SIZE; ++y )
20    {
21        unsigned char *pI = Img.ptr<unsigned char>( y );
22        double y_ = double(y-CENTER) / RADIUS;
23        double y2 = y_*y_;
24        for( int x=0; x<IMG_SIZE; ++x, ++pI )
25        {
26            double x_ = double(x-CENTER) / RADIUS;
27            double x2 = x_*x_;
28            if( x2+y2 > 1.0 )continue;
29
30            double a = 1.0 - Ax_*Ax_ - Ay_*Ay_;
31            double b = 2 * ( -1.0 + Ax_*x_ + Ay_*y_ );
32            double c = 1.0 - x2 - y2;
33
34            double sqrt_part = sqrt( b*b - 4*a*c );
35            double nume = -b - sqrt_part;
36            double denom = 2 * a;
37
38            double t = nume / denom;
39            *pI = cvRound( t * 255 );
40        }
41    }
42    //結果画像の表示
43    cv::imshow( "Img", Img );
44    if( cv::waitKey() == 's' ){ cv::imwrite( "Result.png", Img );   }
45    return 0;
46}

円形の　グラデーションを　描画　することが目的（結果の絵が欲しいだけ）なのであれば，
（コメントにも書いたように）単にその話をそのまま実施する，すなわち，明るさを変えながら円を複数回描画するのが最も素直と言える．（結果画像は割愛）

C++
1const int IMG_SIZE = 240;   //画像サイズ(x,yで共用)
2const int CENTER = IMG_SIZE/2;  //てきとーに決めた円の中心(x,yで共用)
3const int RADIUS = 100; //円の半径
4
5int main()
6{
7    //点Aの座標 = (Ax,Ay)
8    const int Ax = CENTER + 45;
9    const int Ay = CENTER + 25;
10
11    //MEMO : 正規化した世界で考えている変数は，変数名末尾にアンダーバー'_'を付与しています．
12    const double Ax_ = double(Ax-CENTER)/RADIUS;
13    const double Ay_ = double(Ay-CENTER)/RADIUS;
14
15    //画像バッファ準備
16    cv::Mat Img( IMG_SIZE, IMG_SIZE, CV_8UC1 ); //IMG_SIZE*IMG_SIZEの画像バッファを用意
17    Img = 100;  //画像全体を背景の輝度値で埋めておく．値はてきとー．
18
19    //塗りつぶし円を256回描画する
20    for( int t=0; t<=255; ++t )
21    {
22        double t_ = t / 255.0;  //0～1
23        double cx_ = Ax_ * t_;  //円の中心
24        double cy_ = Ay_ * t_;  //円の中心
25        double r_ = 1.0 - t_;   //円の半径
26
27        int cx = CENTER + cvRound( cx_ * RADIUS );
28        int cy = CENTER + cvRound( cy_ * RADIUS );
29        int r = cvRound( r_ * RADIUS );
30        cv::circle( Img, cv::Point(cx,cy), r, cv::Scalar(t), -1 );
31    }
32    
33    //結果画像の表示
34    cv::imshow( "Img", Img );
35    if( cv::waitKey() == 's' ){ cv::imwrite( "Result.png", Img );   }
36    return 0;
37}

投稿2021/01/13 05:08

編集2021/01/19 01:57

fana

総合スコア11708

fana

2021/01/13 05:14

サンプルコードでは，解くべき式の両辺を２乗して，tを２次関数の解の公式を用いて求めている．公式には±があるので，「+と-のどっちを使うべきなの？」という問題があったが，それを決める理屈がよくわからんので，単に２パターン動かしてみてうまく動いた側にしてある．（double nume = -b - sqrt_part;　の箇所）

fana

2021/01/13 07:01

「単色で塗りつぶした円を描く手段」があるなら，それを用いて諧調数（e.g. 8bitなら256）と同じ回数だけ円を描いてやれば同じような絵が描ける．

mamyonya768

2021/01/13 17:04

まだ理解できていませんが、ご回答ありがとうございます！

行動規範の内容に同意します

円を，球を正射影した像だと考え，画像上の画素位置(x,y)に対して，（３次元空間での）円の中心から(x,y)に対応する球面上の点までのベクトルを考えることができます．
各画素の輝度値は，点Aから求めたベクトルと，各画素位置から求めたベクトルとの内積に基づき決定するのが簡単でしょう．

↓のサンプルコードでは絵を作るのに OpenCV を使ってますが，計算内容はコードからわかるかと．

C++
1const int IMG_SIZE = 240;   //画像サイズ(x,yで共用)
2const int CENTER = IMG_SIZE/2;  //てきとーに決めた円の中心(x,yで共用)
3const int RADIUS = 100; //円の半径
4
5//作業関数
6inline bool Unit3DVec( int x, int y, double *pX, double *pY, double *pZ )
7{
8    const double y_ = double(y - CENTER)/RADIUS;
9    const double sqy = y_*y_;
10    if( sqy > 1.0 )return false;
11
12    const double x_ = double(x - CENTER)/RADIUS;
13    const double sqx = x_*x_;
14    if(  sqx + sqy > 1.0 )return false;
15
16    *pX = x_;
17    *pY = y_;
18    *pZ = sqrt( 1.0 - sqx - sqy );
19    return true;
20}
21
22int main()
23{
24    const double POW = 2.0; //明暗変化具合調整パラメタ．
25
26    //点Aの座標 = (Ax,Ay)
27    const int Ax = CENTER + 45;
28    const int Ay = CENTER + 25;
29    double AX,AY,AZ;
30    if( !Unit3DVec( Ax,Ay, &AX,&AY,&AZ ) )return 0;
31
32    //
33    cv::Mat Img( IMG_SIZE, IMG_SIZE, CV_8UC1 ); //IMG_SIZE*IMG_SIZEの画像バッファを用意
34    Img = 100;  //画像全体を背景の輝度値で埋めておく．値はてきとー．
35    
36    //全画素について走査し，その画素の輝度値を決めていく
37    for( int y=0; y<IMG_SIZE; ++y )
38    {
39        unsigned char *pI = Img.ptr<unsigned char>( y );
40
41        for( int x=0; x<IMG_SIZE; ++x, ++pI )
42        {
43            double X,Y,Z;
44            if( !Unit3DVec( x,y, &X,&Y,&Z ) )continue;
45
46            const double InnerProd = X*AX + Y*AY + Z*AZ;
47            if( InnerProd > 0 )
48            {   *pI = cvRound( pow(InnerProd,POW)*255 );    }
49            else
50            {   *pI = 0;    }
51        }
52    }
53    //結果画像の表示
54    cv::imshow( "Img", Img );
55    if( cv::waitKey() == 's' ){ cv::imwrite( "Result.png", Img );   }
56    return 0;
57}