Line Segment Detectorで遊びたかった

はじめに

某学会で道路の白線検出にLine Segment Detectorを使っていると紹介があり、気になったので調べてみました。

なにこれ

エッジを検出して、直線部分を見つけ出して線を引きます。白線の形が出れば道路の状況が分かり、車をどの方向に動かせば良いかを求めることができます。

実装

論文見ながらさくっと作ってみました。ただ、後半のノイズ除去部分がよく分からなかったので検索してみたら本家サイトでソースコードを配布していたので、自分で作るのをやめてしまいました(--;)。
制作途中ですがソースコードおいときます。

f:id:wildpie:20141013212157p:plain

C++
public ref class Pixel
{
public:
	Pixel(int x, int y, double angle, double magnitude) 
            : x(x), y(y), angle(angle), magnitude(magnitude) {}
	int x;
	int y;
	double angle;
	double magnitude;
	bool isUsed;
};

void LSD::Run()
{
	// USB cameraで撮影(OpenCV)
	auto color_image = camera_->GrabImage();
	// ぼかす
	cv::GaussianBlur(color_image, color_image, cv::Size(7, 7), 0.6/0.8);

	pixels_ = gcnew array<Pixel^, 2>(color_image.rows, color_image.cols);

	// 1. Image Scaling
	// Scale factor
	const double s = 0.8; 
	cv::resize(color_image, color_image, cv::Size(), s, s, cv::INTER_CUBIC);

	cv::Mat_<uchar> image;
	cv::cvtColor(color_image, image, CV_BGR2GRAY);

	// 2. Gradient Computation
	data->Clear();
	cv::Mat_<cv::Vec3b> magnitude(image.size());
	ComputeGradient(image, magnitude);

        // System::Drawing::Bitmap^に変換
	camera_image_ = GetBitmap(color_image, 0);
	magnitude_image_ = GetBitmap(magnitude, 1);
}

void LSD::ComputeGradient(cv::Mat_<uchar>& image, cv::Mat_<cv::Vec3b>& m)
{
	for (int y = 0; y < image.rows - 1; y++)
	{
		for (int x = 0; x < image.cols - 1; x++)
		{
			double gx = (image(y, x + 1) + image(y + 1, x + 1) - image(y, x) - image(y + 1, x)) / 2.0;
			double gy = (image(y + 1, x) + image(y + 1, x + 1) - image(y, x) - image(y, x + 1)) / 2.0;

			double angle = atan2(gx, -gy);
			double magnitude = sqrt(gx*gx + gy*gy);

			data->Add(gcnew Pixel(x, y, angle, magnitude));
			pixels_[y, x] = gcnew Pixel(x, y, angle, magnitude);

			m(y, x)[0] = magnitude;
			m(y, x)[1] = magnitude;
			m(y, x)[2] = magnitude;
		}
	}
}
C#

こっちはきれいにしていないので、人に見せるコードじゃないですが・・。

private async Task CaptureCamera()
{
    var rnd = new Random();

    await Task.Run(() =>
        {
            while (true)
            {
                _lsd.Run();

                List<Pixel> data = _lsd.data;

                var rho = (double) numericUpDown1.Value;

                List<Pixel> orderedList = data.Where(v => v.magnitude > rho)
                                              .OrderByDescending(v => v.magnitude)
                                              .ToList();

                foreach (Pixel pixel in orderedList)
                {
                    _lsd.pixels_[pixel.y, pixel.x].isUsed = true;
                }

                Bitmap magnitude = _lsd.magnitude_image_;
                var mi = (Bitmap) _lsd.magnitude_image_.Clone();

                foreach (Pixel p in orderedList)
                {
                    IEnumerable<Pixel> region = RegionGrow(p);

                    double cx = 0;
                    double cy = 0;
                    double angle = CalcRectangle(region, ref cx, ref cy);

                    Color c = Color.FromArgb(rnd.Next(256), rnd.Next(256), rnd.Next(256));
                    foreach (Pixel pp in region)
                        magnitude.SetPixel(pp.x, pp.y, c);

                    //if (!_lsd.pixels_[p.y, p.x].isUsed)
                    //mi.SetPixel((int) cx, (int) cy, c);

                    double minX = region.Select(v => v.x).Min();
                    double maxX = region.Select(v => v.x).Max();
                    double minY = region.Select(v => v.y).Min();
                    double maxY = region.Select(v => v.y).Max();

                    double length = Math.Sqrt((maxY - minY)*(maxY - minY) + (maxX - minX)*(maxX - minX))/2;

                    Graphics grfx = Graphics.FromImage(mi);
                    grfx.DrawLine(Pens.Magenta,
                                  (float) (cx - length*Math.Cos(angle - Math.PI/2.0)),
                                  (float) (cy - length*Math.Sin(angle - Math.PI/2.0)),
                                  (float) (cx + length*Math.Cos(angle - Math.PI/2.0)),
                                  (float) (cy + length*Math.Sin(angle - Math.PI/2.0)));
                }

                pictureBox1.Image = _lsd.camera_image_;
                pictureBox2.Image = magnitude;
                pictureBox3.Image = mi;
            }
        });
}

private IEnumerable<Pixel> RegionGrow(Pixel pixel)
{
    var region = new List<Pixel> {pixel};

    double thetaRegion = pixel.angle;
    double sx = Math.Cos(thetaRegion);
    double sy = Math.Sin(thetaRegion);

    const double tau = Math.PI / 8.0; // 22.5;

    for (int i = 0; i < region.Count(); i++)
    {
        Pixel p = region[i];
        var temp = new List<Pixel>();

        foreach (Pixel q in GetNeighborhood(p))
        {
            if (q == null || _lsd.pixels_[q.y, q.x].isUsed)
                continue;

            if (Math.Abs(thetaRegion - q.angle) < tau)
            {
                temp.Add(q);
                _lsd.pixels_[q.y, q.x].isUsed = true;

                sx += Math.Cos(q.angle);
                sy += Math.Sin(q.angle);
                thetaRegion = Math.Atan2(sy, sx);
            }
        }

        if (temp.Any())
            region.AddRange(temp);
    }

    return region;
}

private IEnumerable<Pixel> GetNeighborhood(Pixel p)
{
    var list = new List<Pixel>();

    int x = p.x;
    int y = p.y;
    if (0 < x && x < _lsd.pixels_.GetLength(1) - 1
        && 0 < y && y < _lsd.pixels_.GetLength(0) - 1)
    {
        list.Add(_lsd.pixels_[y - 1, x - 1]);
        list.Add(_lsd.pixels_[y - 1, x]);
        list.Add(_lsd.pixels_[y - 1, x + 1]);
        list.Add(_lsd.pixels_[y, x - 1]);
        list.Add(_lsd.pixels_[y, x + 1]);
        list.Add(_lsd.pixels_[y + 1, x - 1]);
        list.Add(_lsd.pixels_[y + 1, x]);
        list.Add(_lsd.pixels_[y + 1, x + 1]);
    }

    return list;
}

private double CalcRectangle(IEnumerable<Pixel> region, ref double ccx, ref double ccy)
{
    double sumG = region.Sum(v => v.magnitude);
    ccx = region.Select(v => v.magnitude*v.x).Sum()/sumG;
    ccy = region.Select(v => v.magnitude*v.y).Sum()/sumG;

    double cx = ccx;
    double cy = ccy;
    double mxx = region.Select(v => v.magnitude*Math.Pow(v.x - cx, 2.0)).Sum()/sumG;
    double myy = region.Select(v => v.magnitude*Math.Pow(v.y - cy, 2.0)).Sum()/sumG;
    double mxy = region.Select(v => v.magnitude*(v.x - cx)*(v.y - cy)).Sum()/sumG;

    Func<double, double, double, double> getEigenMax
        = (a, b, c) => ((a + c) + Math.Sqrt((a + c)*(a + c) - 4.0*(a*c - b*b)))/2.0;

    double eigen = getEigenMax(mxx, mxy, myy);

    return Math.Atan2(mxy, mxx - eigen);
}