#目的
何がしたいかというと，画像処理において細線化（Thinning, Skeletonization）は２値画像に変換された画像を，その中心を通る１ピクセル幅の線画像に変換する方法ですが，そこからさらに各終点および分岐点を端点とした複数のセグメントに分割するという話です．
要するに
こういう画像を...
こういう感じで個々のセグメントに分割して，それぞれを点配列で取得したいというわけです．
意外と手こずった上に恐らく自分以外に需要が無さそうなので敢えてシェアしたい衝動に駆られた次第ですが，もっとセーフティで効率のいい方法があれば教えて下さい．
#条件
- セグメント＝２つの端点からなる１本の線とする
- 入力画像は１ピクセル幅の線画像とし，各ピクセルの階調は背景を0，前景を255（0以上）とする
- 出力データは j 個の点配列 p を１つのセグメントとする２重配列 P とする
P_{i, j} = (p_0, p_1, ... p_{j-1} ) \\
p_j = (x, y)
- セグメントは線画像の分岐点および終点を両端とする
- 接合点となるピクセルはそこから得られる全てのセグメントの端点としてそれらに含める
- 線がループ状になっている場合はそれを１つのセグメントとする，ただし，ループ上の線上に２つ以上の分岐点がある場合はそれらの分岐点を端点とするセグメントとして個別に取得する
- 入力画像の線の連結パターンは４近傍および８近傍の両方に対応する
- セグメントの取得順序や，取得されたセグメントの端点が終点であるか分岐点であるかまでは考えなくて良い（とりあえず）
#考え方
いろいろ悩んだ末に以下の手順でまとまりました．
- 左上から１画素ずつ走査し，取得済みでない0以上の画素値のもの（有効画素）があれば2の処理へ，無ければ処理を終了する．
- 最初の端点を探す，すなわち近傍の有効画素を伝わって行き，終点か分岐点に到達した時点で3へ進む．
- 端点の近傍にある有効画素を全て取得し，それらの方向と座標を開始点スタックに貯める．
- 開始点スタックから１つ取り出し，保持していた座標と方向を元に近傍の有効画素を伝わっていく，その際に通った点は出力データとして保持しておく．
- 終点へたどり着いたら場合は通ってきた点を１つのセグメントとして出力データに保持する，分岐点にたどり着いた場合はセグメント取得と同時に，まだ捜査してない近傍の有効画素を開始点スタックに貯める
- 開始点スタックが無くなるまで４〜５の処理を続け，無くなったら７へ
- セグメントとして取得した点はすべて「取得済み」とし，最初の点に戻って再び１の処理へ
#画像の細線化処理
画像の細線化処理についてはこの前Hilditchの方法を書いたので，今回はその計算結果の画像を使います．
#実装
#include <iostream>
#include <vector>
#include <array>
/// Output data type
typedef std::vector<std::array<int, 2> > Point2dArray;
/**
* @brief return point from index
*
* @param i Pixel index
* @param width Image width
*
* @return Point (x, y)
*/
static std::array<int, 2> i2p(int i, int width)
{
int x = i % width;
int y = i / width;
return {x, y};
}
/**
* @brief Divide branched chain.
*
* @param pix Input image pixels (unsigned char, single channel)
* @param width Input image width
* @param height Input image height
* @param dstPts Output points array (vector<Point2dArray>)
* @param bNeighborType4 Input image type, true=4-neighbor, false=8-neighbor.
*
* @return true=Process was succeed
*/
bool divideBranchedChain(const unsigned char* pix, const int width, const int height, std::vector<Point2dArray> *dstPts, bool bNeighborType4)
{
std::vector<std::vector<std::array<int, 2> > > dstBuf;
// offset: [x, y, direction], direction=0...7, 0=right, 1=upper right, ...7=lower right.
std::vector<std::array<int, 3> > offset = {{1,0,0}, {1,-1,1}, {0,-1,2}, {-1,-1,3}, {-1,0,4}, {-1,1,5}, {0,1,6}, {1,1,7}};
// segment: [current focused pixel index, neighbor pixel index, neighbor pixel direction]
std::vector<std::array<int, 3> > neighbors;
std::vector<std::array<int, 3> > nextSearchPixels;
std::unique_ptr<char[]> tPix(new char[width * height]); // searched pix
for (int i = 0; i < width * height; ++i) tPix[i] = 0;
/// (function) find neighbor pixels, and update neighbors.
auto findNeighbor = [&](int focusIdx, int ignoreIdx, int startDir)
{
neighbors.clear();
for (int i = 0; i < offset.size(); ++i)
{
int ti = (i + startDir) % 8;
if (bNeighborType4 && ti % 2 != 0) continue;
int x = offset[ti][0] + (focusIdx % width);
int y = offset[ti][1] + (focusIdx / width);
int testIdx = y * width + x;
// tests:
if (testIdx != ignoreIdx && // is not ignore?
(x >= 0 && x < width && y >= 0 && y < height) && // is not outside?
tPix[testIdx] == 0 && // is not already searched pixel?
pix[testIdx] > 0) // is white pixel
{
neighbors.push_back( {focusIdx, testIdx, offset[ti][2]} );
}
}
};
/// (function) get direction
auto getDirection = [&offset, &width](int focusIdx, int targetIdx)
{
for (const auto& e : offset)
{
int x = e[0] + (focusIdx % width);
int y = e[1] + (focusIdx / width);
int testIdx = y * width + x;
if (testIdx == targetIdx)
return e[2];
}
throw "can't found direction."; // error
};
// scan
for (int i = 0; i < width * height; ++i)
{
if (pix[i] > 0 && tPix[i] == 0)
{
int firstIdx = i;
int focusIdx = firstIdx;
int lastFocus = -1;
int beginDir = 0;
bool bFirst = true;
int count = 0;
try {
// find begin point(s) on a segment
while (count < width * height)
{
findNeighbor(focusIdx, lastFocus, beginDir);
// tests
if (neighbors.empty())
{
if (bFirst)
{
// single pixel
dstBuf.push_back(Point2dArray(1, i2p(focusIdx, width)));
tPix[focusIdx] = 1;
break;
}
// end point
nextSearchPixels.push_back({focusIdx, lastFocus, getDirection(focusIdx, lastFocus)});
break;
}
else {
if (bFirst)
{
if (neighbors.size() == 1)
{
// first pixel is end point
nextSearchPixels.swap(neighbors);
break;
}
}
else {
if (neighbors.size() >= 2 || focusIdx == firstIdx)
{
// branched point, or repeated point
neighbors.push_back({focusIdx, lastFocus, getDirection(focusIdx, lastFocus)});
nextSearchPixels.swap(neighbors);
break;
}
}
}
// continue
lastFocus = focusIdx;
focusIdx = neighbors[0][1];
beginDir = neighbors[0][2];
bFirst = false;
++count;
}
if (count == width * height - 1)
{
throw "endless loop exception";
}
// pick up points on the chains.
while (nextSearchPixels.empty() == false)
{
bFirst = true;
lastFocus = -1;
auto it = nextSearchPixels.begin();
firstIdx = (*it)[0];
focusIdx = firstIdx;
beginDir = (*it)[2];
nextSearchPixels.erase(it);
count = 0;
Point2dArray chainPts;
do {
findNeighbor(focusIdx, lastFocus, beginDir);
// tests
if (neighbors.empty() || (bFirst == false && focusIdx == firstIdx))
{
break;
}
else if (bFirst == false && neighbors.size() >= 2)
{
std::copy(begin(neighbors), end(neighbors), back_inserter(nextSearchPixels));
break;
}
// continue
if (bFirst)
{
chainPts.push_back(i2p(focusIdx, width));
tPix[focusIdx] = 1;
}
lastFocus = focusIdx;
focusIdx = neighbors[0][1];
beginDir = neighbors[0][2];
bFirst = false;
chainPts.push_back(i2p(focusIdx, width));
tPix[focusIdx] = 1;
}
while (count < width * height);
if (count == width * height - 1)
{
throw "endless loop exception";
}
if (chainPts.empty() == false) dstBuf.push_back(chainPts);
}
}
catch (const char* e)
{
std::cout << "[exception] " << e << std::endl;
return false;
}
}
}
dstPts->swap(dstBuf);
return true;
}
for文の中にwhile(true)が２つも登場するというハラハラドキドキなコードですね→ループ回数上限を付けてちょっとだけセーフティにしました．入力画像はあまり多くのパターン試してないですが，無限ループに陥る可能性が無いとはいえないのでいつかもっと整理したいです．
##補足
関数 divideBranchedChain は第１引数に入力画像のポインタ（unsigned char*, シングルチャンネル），第２，３引数にその横幅と高さ，第４引数（出力）に出力データのポインタ，第５引数はtrueだと入力の線画像を４近傍連結，falseだと8近傍連結として処理します．
たぶんC++14じゃないとコンパイルできないかもです．
ところでvector<array<T, N> >っていうの便利ですね．出力データはint型の２次元座標（x, y）なので，本来ならint変数を２個もつ構造体を作ったりするべきだと思いますが，そこは横着してstd::array<int, 2>をvector配列に持たせたりしてます．このvector<array<T, N> > っていうの最近よく使ってます．
あと，とくに意味もなくサブルーチンをラムダ式で書いてたりしますが別の関数として外に出しても大丈夫なはずです．
#計算結果
入力画像
細線化（Hilditchの方法）
セグメント分割
各セグメント毎に色を「赤→青→緑→赤→青...」と分けています．
クリックすると多分もっとよく見えます．
一応ループ状の線や，１ピクセルしかないセグメントもちゃんと取れています．