asaのブログ

OpenGL、WebGL、画像処理についてまとめています

凸包を求める(Grahamスキャン)

凸包を求めます。

1. 凸包について
凸包 - Wikipedia

凸包とはある点の集合が与えられれた時、それを包含する最小の凸多角形あるいは凸多面体のことを言います。

2. 凸包を求めるアルゴリズム

凸包を求めるアルゴリズムにはいくつかの方法があり
(1) ギフト包装法
(2) 分割統治法
(3) Grahamスキャン
などがあります。

今回はGrahamスキャンを使って凸包を求めました。

3. Grahamスキャン

Animation depicting the Monotone algorithm.gif
By Maonus - Own work, CC BY-SA 4.0, Link

上側凸包、下側凸包を分けて求めていきます。上のgifが非常にわかり良いです。各凸包を求める際一番外側の点を探しているのがわかるかと思います。この点は外積を使用して求めていきます。今の点より外側にあるか内側にあるか、あるいは直線上にあるかを求めるのに外積が役に立ちます。

このアルゴリズムを使用するには、まず先に点列をx座標順にソートしておく必要があります。

C++でこれを実装するに当たって、標準ライブラリのvectorとsort関数を使用しました。自作のVectorクラスはsortに対応していなかったので、比較関数を次のように定義しました。

bool Vector::compareX(const Vector& vec1, const Vector& vec2) {
    return vec1.x < vec2.x || (vec1.x == vec2.x && vec1.y < vec2.y);
}

ソートが完了した後、上のgifで見たように順番に凸包を求めていきます。凸包を構成する点列が取得できれば、あとはそれをCanvasあるいはWebGLなどで描画すると確認できるでしょう(下のコードはAizu Onlineジャッジで動作を確認しました)。

// 凸包を取得する
// https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#C++
vector<Vector> Polygon::getConvexHull(vector<Vector>& src) {
    sort(src.begin(), src.end(), Vector::compareX);
    int n = (int)src.size();
    int k = 0;
    vector<Vector> ver(2 * n);
    // 下側凸包を求める
    for (int i = 0; i < n; ++i) {
        while(k >= 2 && (ver[k - 1] - ver[k - 2]).cross(src[i] - ver[k - 2]) < 0.0) { // 3点の外積を求める
            k--; // 反時計周りのとき、src[i]は内側の点なので除く
        }
        ver[k++] = src[i];
    }
    // 上側凸包を求める
    for (int i = n - 1, t = k + 1; i > 0; --i) {
        while(k >= t && (ver[k - 1] - ver[k - 2]).cross(src[i - 1] - ver[k - 2]) < 0.0) { // 3点の外積を求める
            k--; // 反時計周りのとき、src[i - 1]は内側の点なので除く
        }
        ver[k++] = src[i - 1];
    }
    ver.resize(k - 1);
    return ver;
}

多角形の面積

多角形の面積を求める。

1. 公式について
多角形の座標が分かっているときには、面積を求める公式が存在する
n個の頂点を起き、i番目の頂点を vi とする。2つのベクトルの外積をv1. cross (v2)のように表すとする。
このとき多角形の面積は

1 / 2 Σ vi. cross(v(i+1))

で求めることができる。 ただしiがnに達したときi+1は最初の要素を指すものとする。

コードは下のような感じになります。gitに幾何ライブラリをまとめていますので、Polygonクラスの詳細等についてはこちらから確認お願いいたします(https://github.com/asaren08/GLib)。

// 面積を計算する
double Polygon::calculateArea() {
    double sum = 0;
    for (int i = 0; i < vertex.size(); i++) {
        Vector vec1, vec2;
        vec1 = vertex[i];
        if (i == vertex.size() - 1) {
            vec2 = vertex[0];
        } else {
            vec2 = vertex[i + 1];
        }
        sum += vec1.cross(vec2);
    }
    return sum / 2.0;
}


X. 参考
http://www004.upp.so-net.ne.jp/s_honma/urawaza/area2.htm

線分間の距離について

線分間の距離を求める。

1. 場合分け

線分が交差している時距離は0とする。また線分同士が接している時も同様とする。
線分同士の最短距離は、線分上にある場合とない場合に分けて考えることができる。線分Bが線分Aの直上にあるとき、線分間の距離は線分Bの端点で線分Aに近いものが、線分Aに落とした垂線の距離でとることができる。また直上にない時は線分Aのいずれかの端点と線分Bのいずれかの端点との距離が最短距離となる。

2. 線分が自分からみてどこにいるか?

LineSegmentという始点と終点のベクトルを持つクラスをおく。最初に接しているかどうかを判定し、接していなければ距離を計算していく。
 計算のための場合分けはまず始点よりも外にあるか内にあるかで考える。始点よりも外にあれば、始点との距離を求めるだけでよい。内にあれば、さらに場合分けし、直上にあるのか、終点より外側にあるのかを判定する。直上にあれば、始点との正射影ベクトルを求め、垂線の距離を計算し、外にあれば、終点との距離を計算する。

// 任意の点が線分の始点より外にあるか
bool LineSegment::isOutOfStart(const Vector& vec) const {
    Vector p = end - start;
    Vector q = vec - start;
    double d = p.dot(q);
    if (d < 0) { // 内積が負のとき
        return true;
    }
    return false;
}
// 任意の点との距離を求める
// http://gihyo.jp/dev/serial/01/as3/0053?page=2
double LineSegment::getDistance(const Vector& vec) const {
    if (isPointInLine(vec)) { // 線分内に点を含むとき
        return 0.0;
    }
    if (!isOutOfStart(vec)) { // 線分の始点より内側に点があるとき
        Vector p = end - start;
        Vector q = vec - start;
        Vector proj = Vector::getProjection(q, p);
        if (p.norm() > proj.norm()) { // 正射影したベクトルが線分より短いとき
            return (q - proj).norm();
        } else { // 正射影したベクトルが線分と等しいか、長いとき
            return (vec - end).norm();
        }
    } else { // 線分の始点より外にあるとき
        return (vec - start).norm();
    }
}

上の計算をそれぞれの端の点、正味4点について計算する。

// 任意の線分との距離を求める
double LineSegment::getDistance(const LineSegment& l) const {
    if (isCross(l)) { // 交差しているとき
        return 0.0;
    }
    
    LineSegment ls;
    ls.setStartEnd(start.x, start.y, end.x, end.y);
    double dis1 = getDistance(l.getStart());
    double dis2 = getDistance(l.getEnd());
    double dis3 = l.getDistance(ls.getStart());
    double dis4 = l.getDistance(ls.getEnd());
    
    double min1 = fmin(dis1, dis2);
    double min2 = fmin(dis3, dis4);
    return fmin(min1, min2);
}

線分の交差判定について

線分Aと線分Bがあったとき、AとBが交差しているか求める。

1. ベクトルを使う

線分は始点と終点を持ったベクトルで表すことができる。
LineSegmentというクラスをおいて、メンバに始点ベクトルと終点ベクトルを持たせる。

2. 外積を使って交差判定する
あるベクトルに対して別のベクトルが右にあるか、左にあるかは外積を使って求めることができる。

外積は簡単に求めることができて、次の計算で求める。

    // 外積を求める
    float cross(const Vector& vec) const {
        return x * vec.y - y * vec.x;
    }

線分Aが線分Bと交差しているかを判定するには、線分Aから見て線分Bの端の点が自分の右と左にあるかを調べる。また線分Bから見て、線分Aの端の点が右と左にそれぞれあるかを見る。

平行の時は、線分Aと線分Bが重なっているかを見て、交差を判定する。

    // 指定の線分と交差しているか
    // @return 0(交差していない), 1(交差している)
    int isCross(const LineSegment& l) {
        Vector l1 = end - start;
        float c1 = l1.cross(l.getStart() - start);
        float c2 = l1.cross(l.getEnd() - start);
        
        Vector l2 = l.getEnd() - l.getStart();
        float c3 = l2.cross(start - l.getStart());
        float c4 = l2.cross(end - l.getStart());
        
        if (c1 * c2 == 0.0 && c3 * c4 == 0.0) { // 平行のとき
            LineSegment ls;
            ls.setStartEnd(start.x, start.y, end.x, end.y);
            if (isInLine(l) || l.isInLine(ls)) {
                return 1;
            }
            return 0;
        } else if ((c1 * c2 < 0.0 && c3 * c4 <= 0.0)
                   || (c1 * c2 <= 0.0 && c3 * c4 < 0.0)
                   || (c1 * c2 < 0.0 && c3 * c4 < 0.0)) { // 交差するとき
            return 1;
        } else { // 交差しないとき
            return 0;
        }
    }

平行のとき、線分がある線分の中に含まれるかどうかは、それぞれの端の点が線分の中に含まれるかを見る。

線分の中に点が含まれるときは、線分ベクトルの単位ベクトルと始点から点までのベクトルの単位ベクトルを求めて一致しているか、そして線分ベクトルの長さの中に始点から点までのベクトルの長さがおさまっていることを確認する。

    // 指定の線分を含むか
    // このメソッドは線分同士が平行の時しか使用しない
    bool isInLine(const LineSegment& l) const {
        bool b1 = isPointInLine(l.getStart());
        bool b2 = isPointInLine(l.getEnd());
        if (b1 || b2) {
            return true;
        } else {
            return false;
        }
    }
    
    // 指定の点を線分の中に含むか
    // このメソッドは平行の時しか使用しない
    bool isPointInLine(const Vector& p) const {
        if (start == p) {
            return true;
        }
        
        Vector l1 = end - start;
        Vector l2 = p - start;
        Vector u1 = l1.getUnit();
        Vector u2 = l2.getUnit();

        if (u1 == u2 && l2.norm() <= l1.norm()) {
            return true;
        } else {
            return false;
        }
    }

線対称なベクトルを求める

正射影が求まると線対称なベクトルも簡単に求まる。
求めた正射影ベクトルへの垂線となるベクトルの根元を正射影ベクトルまで持って来れば良い。


gist8d799076706d3bf23c97480a4aa0a9b6

正射影ベクトルを計算する

▪️正射影を求める
たまに思い出すとこうしてあれしてと手順が浮かんでも、うっかり間違う可能性があるのでメモ。

VectorAとVectorBがあったとき、VectorAをVectorBに正射影する。
VectorBの単位ベクトルを e とおいて、2つのベクトルの余弦を cos 、VectorAの長さを norm とすると
norm *e * cos
で求まる。VectorAを最後に足してやるのをたまに忘れるので注意。

▪️コード


gist7044b79909905c22a2ca8ba2fb255347

【WebGL】鮮鋭化

■鮮鋭化後 
f:id:asa_r:20180225123530p:plain

■元画像
f:id:asa_r:20180225123628j:plain

カーネルの変更
 鮮鋭化はカーネルの変更だけでいけます。

gist031eaed25b867f5739548db196c41feb


 なお、一緒に載っていたエンボスのカーネルは下の通り。

gist07f80c810336822ea865f397c7eed5da

 使って見るとこんな感じ。
f:id:asa_r:20180225124857p:plain

■参考
8.2. Convolution Matrix