因みにOpenCV HoughCircleで円を検出
HoughCircleを使うなら以下
var circles = Cv2.HoughCircles(blur, HoughMethods.Gradient, dp, minDist, param1, param2, minRadius, maxRadius);
パラメータの意味についてはこちらが分かりやすかったです。
しかし、私の検出したかった円がうまく検出できません。
パラメータを調整しながら様子を見ましたが、検出しすぎるか、足りないか、と言った結果になりました。
検出数を多くしてフィルタリングすれば・・・！と思ったのですが、処理に結構な時間がかかったり、期待した円と違う結果になったりします。
そこで自前で円検出を作って目的に近いものが出来るか、やってみました。
自前で円の検出に挑戦
コード中のベタ書きの数字は対象画像によってかなり変わります。
今回の対象は
- 元画像がほぼグレースケール
- 検出対象はほぼ真っ黒
- 画像中の2割くらいの面積は検出対象
- 検出対象の周りはほぼ白
- 複雑な線はほとんど存在しない
といった画像です。
普通の風景写真などから検出しようとする場合、この数値のままでは厳しいと思います。
まだ実験用の為見苦しい点が多いですがご容赦ください。
円の計算式(円弧上の座標から円の中心と半径を計算する)
円弧から半径、中心点を計算する為に、まず計算式を調べます。
最小二乗法で近似円計算すれば良いんでしょ？位の知識から調べました。
ここで丁寧に説明されてます。
上の参考のページに一部間違いがある様です。
下の方のコード中でmatrix2としている部分の符号がマイナスでなくプラスの様です。
誤
\left(
\begin{matrix}
-\sum_{i=1}^n (x_i^3 - x_iy_i^2) \\
-\sum_{i=1}^n (x_i^2y_i - y_i^3) \\
-\sum_{i=1}^n (x_i^2 - y_i^2)
\end{matrix}
\right)
正
\left(
\begin{matrix}
-\sum_{i=1}^n (x_i^3 + x_iy_i^2) \\
-\sum_{i=1}^n (x_i^2y_i + y_i^3) \\
-\sum_{i=1}^n (x_i^2 + y_i^2)
\end{matrix}
\right)
・・・行列の反転、積を計算しなければならないんですね。
そういえば今までC#で行列演算をしてこなかったような。
.NetにはDirectX向けか何かのMatrixクラスは有りましたが、どうも今回の目的には使えない模様。
Math.Net Numericsを使うのが楽そうですが、
訳有ってこれが使いづらい状況ですので自作してみることにします。
円弧計算のプログラム
計算結果保持用クラス
OpenCVにはCircleSegmentと言う構造体が有り、円の情報を持たせるにはこれが良さそうです。
しかし、
- パラメータの数値をfloatで持つため内部の桁が大きくなる箇所があり変換が必要。 (円弧のパラメータ計算時に大きな桁になる恐れがあります)
- デバッグ用に輪郭の座標を保持したい
と思い自作クラスを作りました。
このあたりは状況次第でしょうか。
/// <summary>
/// 円弧の計算結果保持クラス
/// OpenCvSharp.CircleSegmentでは核パラメータがfloatであった。
/// doubleで持ちたかったので自作
/// </summary>
public class Arc
{
/// <summary>
/// 中心点
/// </summary>
public Point2d Center { get; protected set; }
/// <summary>
/// 円弧の中心位置 X
/// 返した先でOpenCvを参照しなくてもいいように残します
/// </summary>
public double X { get { return Center == null ? double.NaN : Center.X; } }
// <summary>
/// 円弧の中心位置 Y
/// 返した先でOpenCvを参照しなくてもいいように残します
/// </summary>
public double Y { get { return Center == null ? double.NaN : Center.Y; } }
/// <summary>
/// 円弧の半径
/// </summary>
public double Radius { get; set; }
/// <summary>
/// 検出した輪郭の座標リスト
/// デバッグ用の意味合いが強い
/// </summary>
public List<Point2d> Contour { get; protected set; }
/// <summary>
/// コンストラクタ
/// </summary>
/// <param name="center"></param>
/// <param name="radius"></param>
/// <param name="contour"></param>
public Arc(Point2d center, double radius, IEnumerable<Point2d> contour)
{
Center = center;
Radius = radius;
Contour = contour.ToList();
}
public Arc(ValueType x, ValueType y, double radius, IEnumerable<Point2d> contour) : this(new Point2d((double)x, (double)y), radius, contour) { }
}
計算の主体部分
結構な行数になってしまいました。
最初1000点超えのPointを計算したらLINQのSumでOverflowExceptionが発生したり、時間がかかったりと紆余曲折有りました。
Overflowは久々に直面しましたね。コレクションの要素の2乗3乗の総和があるのでこの辺りが結構桁が大きくなります。
これに渡す点数は最低3点。
上手く取れていれば5点でかなり近い円が取れました。
現在10点でやってますが、
Overflow対策で点数を絞るなり、要素の数値が余りにも大きい場合は一旦小さくするなどの工夫が必要かもしれません。
using OpenCvSharp;
~~
/// <summary>
/// 円弧計算のクラス
/// 線上の座標から最小二乗法に則って円の中心点と半径を計算します。
/// </summary>
public static class CalcArc
{
/// <summary>
/// 線上の座標から最小二乗法に則って円の中心点と半径を計算します。
/// </summary>
/// <param name="points">線上の座標</param>
/// <returns>計算結果を保持する円弧クラス</returns>
public static Arc Calc(IEnumerable<Point2d> points)
{
Func<IEnumerable<Point2d>, double> Func_sigmaXsquare = (points_local) =>
{
return points_local.Select(p => p.X * p.X).Sum();
};
Func<IEnumerable<Point2d>, double> Func_sigmaYsquare = (points_local) =>
{
return points_local.Select(p => p.Y * p.Y).Sum();
};
Func<IEnumerable<Point2d>, double> Func_sigmaXY = (points_local) =>
{
return points_local.Select(p => p.X * p.Y).Sum();
};
Func<IEnumerable<Point2d>, double> Func_sigmaX = (points_local) =>
{
return points_local.Select(p => p.X).Sum();
};
Func<IEnumerable<Point2d>, double> Func_sigmaY = (points_local) =>
{
return points_local.Select(p => p.Y).Sum();
};
//Func<IEnumerable<double>, double> LocalFunc_Sum = (list) =>
// {
// double rtn = 0;
// foreach (var element in list)
// {
// rtn += element;
// }
// return rtn;
// };
double sigmaXSquare = Func_sigmaXsquare(points);
double sigmaYSquare = Func_sigmaYsquare(points);
double sigmaXY = Func_sigmaXY(points);
double sigmaX = Func_sigmaX(points);
double sigmaY = Func_sigmaY(points);
double sigma1 = (double)points.Count();
var matrix1 = new double[,]
{
{ sigmaXSquare, sigmaXY, sigmaX }
,{ sigmaXY, sigmaYSquare, sigmaY }
,{ sigmaX, sigmaY, sigma1 }
};
matrix1 = InvertMatrix(matrix1);
var matrix2 = new double[,]
{
//System.Drawin.Point または OpenCvSharp.Point
//を多く詰め込んだpointsを渡されると、int型で演算されここでOverflowExceptionが発生します。
//doubleで扱う様にするとかなり大きな数値を扱わない限りOverflowしません
//が、しかし、画像の画素数が多いと行列演算中のx^3、y^3等で結構な桁までいってしまいます。
//よってここで渡すpointsはそもそも数を絞ってやるのが安全です。
//今回はそこまで書いてません。
{ -1* points.Select(p => p.X * p.X * p.X + p.X * p.Y * p.Y).Sum() }
,{ -1* points.Select(p => p.X * p.X * p.Y + p.Y * p.Y * p.Y).Sum() }
,{ -1* points.Select(p => p.X * p.X + p.Y * p.Y).Sum() }
//foreachで計算するのとSumでは特に差が無いのでSum()にします↑
//{ -1* LocalFunc_Sum(points.Select(p => p.X * p.X * p.X + p.X * p.Y * p.Y)) }
//,{ -1* LocalFunc_Sum(points.Select(p => p.X * p.X * p.Y + p.Y * p.Y * p.Y)) }
//,{ -1* LocalFunc_Sum(points.Select(p => p.X * p.X + p.Y * p.Y)) }
};
var abcMat = MultipleMatrix(matrix1, matrix2);
var a = (-1 * abcMat[0, 0]) / 2;
var b = (-1 * abcMat[1, 0]) / 2;
var rad = Math.Sqrt((-1 * abcMat[2, 0]) + a * a + b * b);
return new Arc(a,b, rad, points);
}
/// <summary>
/// 3*3行列と1*3行列の積を計算します。
/// </summary>
/// <param name="m1">3*3行列</param>
/// <param name="m2">1*3行列</param>
/// <returns>演算結果の行列(1*3)</returns>
private static double[,] MultipleMatrix(double[,] m1, double[,] m2)
{
if (m1.GetLongLength(0) != 3 || m1.GetLongLength(1) != 3 || m2.GetLength(0) != 3 || m2.GetLength(1) != 1)
{
throw new ArgumentOutOfRangeException("行列の大きさが対応外です");
}
double a = m1[0, 0], b = m1[0, 1], c = m1[0, 2];
double d = m1[1, 0], e = m1[1, 1], f = m1[1, 2];
double g = m1[2, 0], h = m1[2, 1], i = m1[2, 2];
double j = m2[0, 0], k = m2[1, 0], l = m2[2, 0];
return new double[,]
{
{ c * l + b * k + a * j}
,{ f * l + e * k + d * j}
,{ i * l + h * k + g * j}
};
}
/// <summary>
/// 3*3行列を反転します。(逆数)
/// </summary>
/// <param name="m">3*3行列</param>
/// <returns>反転した行列(3*3)</returns>
private static double[,] InvertMatrix(double[,] m)
{
var denom = CalcDenominator(m);
if (m.GetLength(0) != 3
|| m.GetLength(1) != 3)
{
throw new ArgumentOutOfRangeException("行列の長さは3物にしか対応しません");
}
double a = m[0, 0], b = m[0, 1], c = m[0, 2];
double d = m[1, 0], e = m[1, 1], f = m[1, 2];
double g = m[2, 0], h = m[2, 1], i = m[2, 2];
return new double[,]
{
{(e * i - f * h) * denom,
-1*((b * i - c * h) * denom),
(b * f - c * e) * denom,
}
,{-1*((d * i - f * g) * denom),
(a * i - c * g) * denom,
-1*((a * f - c * d) * denom),
}
,{(d * h - e * g) * denom,
-1*((a * h - b * g) * denom),
(a * e - b * d) * denom,
}
};
}
/// <summary>
/// 行列反転演算の係数を計算します。
/// </summary>
/// <param name="inputMatrix">反転する行列</param>
/// <returns>係数</returns>
private static double CalcDenominator(double[,] inputMatrix)
{
if (inputMatrix.GetLength(0) != inputMatrix.GetLength(1))
{
throw new ArgumentOutOfRangeException("行列の長さが揃ったものにしか対応できません");
}
var vals1 = new double[inputMatrix.GetLength(0)];
for (int i = 0; i < inputMatrix.GetLength(0); i++)
{
vals1[i] = inputMatrix[0, i];
for (int j = 1; j < inputMatrix.GetLength(0); j++)
{
int index = i + j;
if (index >= inputMatrix.GetLength(0))
{
index -= inputMatrix.GetLength(0);
}
vals1[i] *= inputMatrix[j, index];
}
}
var vals2 = new double[inputMatrix.GetLength(0)];
for (int i = 0; i < inputMatrix.GetLength(0); i++)
{
vals2[i] = inputMatrix[0, inputMatrix.GetLength(1) - 1 - i];
for (int j = 1; j < inputMatrix.GetLength(0); j++)
{
int index = inputMatrix.GetLength(1) - 1 - i - j;
if (index < 0)
{
index += inputMatrix.GetLength(0);
}
vals2[i] *= inputMatrix[j, index];
}
}
return 1 / (vals1.Sum() - vals2.Sum());
}
}
円の輪郭検出
OpenCvSharpで輪郭検出しました。
ここは自前より数段いいですね。
色々と調整項目が有りました。
領域を限定してしまえば2値化してしまうのが楽そうです。
が、それがちょっとやりづらかったので輪郭検出の肝はCannyで行いました。
対象画像では割と思った通りの輪郭が取れました。
using OpenCvSharp;
using OpenCvSharp.Extensions;
~~
/// <summary>
/// 円弧検出クラス
/// </summary>
public static class DetectCircle
{
/// <summary>
/// 円弧を計算して返します
/// </summary>
/// <param name="bmp">円弧が描かれているBitmap画像</param>
/// <param name="outLine">計算対象エリア(この範囲でトリミングして計算します。 )</param>
/// <param name="maskArea">画像中のマスク領域</param>
/// <param name="maxRadius">最大半径 : この値を超える半径の円弧は返り値から除外します。</param>
/// <param name="minRadius">最小半径 : この値を下回る半径の円弧は返り値から除外します。</param>
/// <returns>計算した円弧の結果を詰めた配列</returns>
public static Arc[] GetArcs(Bitmap bmp, Rectangle outLine, IEnumerable<Rectangle> maskArea = null, double maxRadius = double.MaxValue, double minRadius = 0)
{
//円弧計算結果の値チェック 数値か？上限以下？下限以下？
Func<double, double, double, bool> checkValueRange = (input, max, min) =>
{
return !double.IsInfinity(input)
&& !double.IsNaN(input)
&& max >= input && min <= input;
};
//直線らしき座標を間引く為、XY座標の変化量がある点のみに絞る
Func<OpenCvSharp.Point[], int, int> findNextIndex = (inputContour, startIndex) =>
{
int j = startIndex;
do
{
j++;
if (inputContour.Length <= j)
{
return -1;
}
} while (Math.Abs(inputContour[startIndex].X - inputContour[j].X) <= 1
|| Math.Abs(inputContour[startIndex].Y - inputContour[j].Y) <= 1);
return j;
};
//直線らしき座標を間引く為、変化点の傾きの違う点を抽出する (直線は傾きが均一)
//contourは概ね近い連続点となっているので順番にアクセスすれば輪郭の端から順番にアクセスできる。
Func<OpenCvSharp.Point[], IEnumerable<Point2d>> searchNotStraightContour = (inputContour) =>
{
var notStraightContour = new List<Point2d>();
for (int i = 0; i < inputContour.Length;)
{
//iの次の座標変化がある点を見つける
int j = findNextIndex(inputContour, i);
if (j < 0)
{
break;
}
//jの次の座標変化がある点を見つける
int k = findNextIndex(inputContour, j);
if (k < 0)
{
break;
}
//それぞれの点の差を取り
var x01 = inputContour[j].X - inputContour[i].X;
var y01 = inputContour[j].Y - inputContour[i].Y;
var x02 = inputContour[k].X - inputContour[i].X;
var y02 = inputContour[k].Y - inputContour[i].Y;
//1点目と2点目の傾きを算出
double xy01 = (double)y01 / (double)x01;
//1点目と3点面の傾きを算出
double xy02 = (double)y02 / (double)x02;
if (!double.IsNaN(xy01) && !double.IsInfinity(xy01)
&& !double.IsNaN(xy02) && !double.IsInfinity(xy02)
&& Math.Abs(x01) < maxRadius * 2 //最大直径以上の距離が離れた点なら判定しない
&& Math.Abs(y01) < maxRadius * 2 //最大直径以上の距離が離れた点なら判定しない
&& Math.Abs(x02) < maxRadius * 2 //最大直径以上の距離が離れた点なら判定しない
&& Math.Abs(y02) < maxRadius * 2 //最大直径以上の距離が離れた点なら判定しない
&& Math.Abs(xy01 - xy02) > 0.025//傾きの差が小さすぎる点は直線であろう。<--------要調整
)
{
notStraightContour.Add(inputContour.ElementAt(j));
}
i = k + 1;
}
return notStraightContour;
};
//画像中の無視する範囲の輪郭を除去します。
Func<OpenCvSharp.Point[], OpenCvSharp.Point[]> getMaskContuour = (inputPoints) =>
{
List<OpenCvSharp.Point> maskedContour = new List<OpenCvSharp.Point>();
if (maskArea == null)
{
//画像中の無視領域は指定されていないのでそのまま返す。
maskedContour.AddRange(inputPoints);
}
else
{
foreach (var con in inputPoints)
{
if (maskArea
.Where(a => a.Left < con.X && a.Right > con.X
&& a.Top < con.Y && a.Bottom > con.Y).Any())
{
continue;
}
maskedContour.Add(con);
}
}
return maskedContour.ToArray();
};
var inputMat = BitmapConverter.ToMat(bmp).CvtColor(ColorConversionCodes.RGB2BGR);
var trimMat = inputMat[new Rect(outLine.Left, outLine.Top, outLine.Width, outLine.Height)];
var gray = trimMat.CvtColor(ColorConversionCodes.BGR2GRAY);
//binaryにすると円以外の輪郭が強く出てしまい、これを除去する必要がある。今回は不採用
//var bin_img = gray.Threshold(12, 255, ThresholdTypes.Binary);//.Otsu);
//ガウシアンフィルタのカーネルサイズは画像(画素数)によって要調整。画素数が多いほど大きくした方が良さそうである。
var blur = gray.GaussianBlur(new OpenCvSharp.Size(9, 9), 0);
//画像次第でここの数値は要調整 なるべく数値を上げて輪郭数を少なくした方が速いが、輪郭が短くブツ切れになると円の検出が厳しい。
var canny_img = blur.Canny(100, 250);
OpenCvSharp.Point[][] contours;
HierarchyIndex[] hierarchyIndexes;
canny_img.FindContours(out contours, out hierarchyIndexes, RetrievalModes.CComp, ContourApproximationModes.ApproxNone);
//検出した円弧を詰めるList
var rtnList = new List<Arc>();
//見つけた輪郭を走査
foreach (var contour in contours)
{
//輪郭中の計算点数
int numOfPoints = 10;
//画像中の処理対象外にある輪郭を除去する
var maskedContour = getMaskContuour(contour);
//完全な直線と疑われるものは除く
var notStraightContour = searchNotStraightContour(maskedContour);
//対象が無ければ処理しない
if (!notStraightContour.Any()
|| notStraightContour.Count() < numOfPoints)
{
continue;
}
//間引きしてdouble型のPointコレクションを取得
var subContour = notStraightContour
//間引きする
.Where((a, index) => index % (notStraightContour.Count() / numOfPoints) == 0)
//外形トリミング分の座標補正
.Select(a => new Point2d(a.X + outLine.X, a.Y + outLine.Y));
//円弧を計算
var arc = CalcArc.Calc(subContour);
if (checkValueRange(arc.X, double.MaxValue, 0)
&& checkValueRange(arc.Y, double.MaxValue, 0)
&& checkValueRange(arc.Radius, maxRadius, minRadius)
)
{
rtnList.Add(arc);
}
}
return rtnList.ToArray();
}
}
所感
- 自前計算はそこそこ速度も出たので満足
- HoughCicle結構使いづらい
- でもOpenCV凄い