全象限bresenham algorithm python/c++实现及可视化

GoodluckTian

已于 2024-05-08 11:49:59 修改

阅读量927

点赞数 2

分类专栏：不知道放哪儿文章标签： python c++

于 2023-01-13 14:19:42 首次发布

本文为博主原创文章，未经博主允许禁止转载。如有问题，欢迎指正！

本文链接：https://blog.csdn.net/weixin_43152152/article/details/128673607

版权

不知道放哪儿专栏收录该内容

5 篇文章 0 订阅

订阅专栏

#! https://zhuanlan.zhihu.com/p/598780689

bresenham algorithm

全象限区域bresenham algorithm计算的python/c++实现

bresenham algorithm为计算机图形学中使用像素点显示直线的算法，算法使用整数运算，能大幅提升计算速度。最近概率栅格建图算法中涉及到直线绘制，故推导学习了下。

公式推导

主要是第一象限斜率小于1的区域的公式推导，其他区域可以转换到该区域计算。

斜率公式如下，其中d为上图中每移动一个像素，y轴上增加的数值，即斜率
$\frac{y_1 - y_0}{x_1 - x_0} = \frac{\delta y}{\delta x} = d$
此时判断每次x增长时，纵轴的增值，对第一个像素判断d>0.5，是则y+1，下一个像素则判断2d>1.5，否则y不变，下一个像素则判断3d>1.5。

因为上一个像素没变，对比的阈值依旧保持不变。如上图示例：

当 $x=x_1 = x_0+1$ 时，d>0.5，则 $y_1=y_0+1$ ，y向上移动一个像素。d-0.5>0 => 2*dy-dx > 0

当 $x=x_2 = x_0+2$ 时，2d<1.5，则 $y_2=y_1$ ，y保持不变，此时1.5保持不变 2d-1.5>0 => 2*dy-dx + 2*dy-2*dx<0

当 $x=x_2 = x_0+3$ 时，3d>1.5，则 $y_3=y_2+1$ ，y向上移动一个像素 3d-1.5>0 => 2*dy-dx +2*dy-2*dx + 2*dy>0

递推公式如下：
$D_k = \begin{cases} 2 * dy -dx &k = 0 \\ D_{k-1} + 2*dy & k>0且D_{k-1} < 0 \\ D_{k-1} + 2*dy-2*dx & k>0且D_{k-1} > 0 \\ \end{cases}$

伪代码如下：

plotLine(x0, y0, x1, y1)
    dx = x1 - x0
    dy = y1 - y0
    D = 2 * dy - dx
    y = y0

    for x from x0 to x1
        plot(x, y)
        if D > 0
            y = y + 1
            D = D - 2*dx
        end if
        D = D + 2*dy

再精简为如下：

plotLine(x0, y0, x1, y1)
    dx = x1 - x0
    dy = y1 - y0
    D = 0
    y = y0

    for x from x0 to x1
        plot(x, y)
        D += dy
        if 2 * D > 0
        	D = D - dx
            y = y + 1
        end if

所有情况讨论

以(x0,y0)为原点，可分为如下8种+平行xy轴的10种情况讨论

以下实现以(x0,y0)为原点实现，如果(x0,y0)不是原点，需要先把数据平移到原点，再平移回去，或者修改TransformQuadrant函数

python实现及可视化

import math

import matplotlib.pyplot as plt
import numpy as np


# 全象限
def BresenhamAlgorithm(x0, y0, x1, y1):
    # 1.process parallel situation
    if x0 == x1 and y0 == y1:
        return [x0, y0]
    elif x0 == x1:
        if y0 < y1:
            y_min = y0
            y_max = y1
        else:
            y_min = y1 + 1
            y_max = y0 + 1

        result = []
        for y in range(y_min, y_max):
            result.append([x0, y])
        return result
    elif y0 == y1:
        if x0 < x1:
            x_min = x0
            x_max = x1
        else:
            x_min = x1 + 1
            x_max = x0 + 1
        result = []
        for x in range(x_min, x_max):
            result.append([x, y0])
        return result

    origin_x = x0
    origin_y = y0
    x0 -= origin_x
    y0 -= origin_y
    x1 -= origin_x
    y1 -= origin_y

    situation = 0
    if x1 > x0 and y1 > y0:
        if (y1 - y0) < (x1 - x0):
            situation = 11
        else:
            situation = 12
    elif x1 < x0 and y1 > y0:
        if (y1 - y0) < (x0 - x1):
            situation = 24
        else:
            situation = 23
    elif x1 < x0 and y1 < y0:
        if (y0 - y1) < (x0 - x1):
            situation = 35
        else:
            situation = 36
    elif x1 > x0 and y1 < y0:
        if (y0 - y1) < (x1 - x0):
            situation = 48
        else:
            situation = 47

    # transform quadrant-2/3/4 to quadrant-1, or transform back
    def Swap(xt, yt):
        tmp = xt
        xt = yt
        yt = tmp
        return xt, yt

    def TransformQuadrant(xt, yt, pose, back=False):
        if pose == 12:
            xt, yt = Swap(xt, yt)
        elif pose == 23:
            xt = -xt
            xt, yt = Swap(xt, yt)
        elif pose == 24:
            xt = -xt
        elif pose == 35:
            xt = -xt
            yt = -yt
        elif pose == 36:
            xt = -xt
            yt = -yt
            xt, yt = Swap(xt, yt)
        elif pose == 47:
            yt = -yt
            xt, yt = Swap(xt, yt)
        elif pose == 48:
            yt = -yt

        if back:
            if pose == 23 or pose == 47:
                xt = -xt
                yt = -yt
        return xt, yt

    # 3. transform to quadrant-1_1
    # print(f"before {x1 - x0}, {y1 - y0}")
    x1, y1 = TransformQuadrant(x1, y1, situation)
    # print(f"after {x1 - x0}, {y1 - y0}")

    # 4. compute grid in line
    delta_x = x1 - x0
    delta_y = y1 - y0
    error = 0
    y = y0
    result = []
    for x in range(x0, x1):
        result.append([x, y])
        error += delta_y
        if 2 * error > delta_x:
            error -= delta_x
            y += 1

    # 5. transform back to original quadrant
    for pos in result:
        pos[0], pos[1] = TransformQuadrant(pos[0], pos[1], situation, True)

    for pos in result:
        pos[0] += origin_x
        pos[1] += origin_y

    return result


# grid map range
range_max = 70
grid_map = [[0 for i in range(range_max)] for j in range(range_max)]

# # 设置x0,y0
xx, yy = 10, 10
# # xx, yy = 0, 0
# test_pts = [
#     [xx - 10, yy], [xx + 10, yy], [xx, yy - 10], [xx, yy + 10],  # 平行轴
#     [xx + 15, yy + 10], [xx + 10, yy + 15],  # 第一象限
#     [xx - 15, yy + 10], [xx - 10, yy + 15],  # 第二象限
#     [xx - 15, yy - 10], [xx - 10, yy - 15],  # 第三象限
#     [xx + 15, yy - 10], [xx + 10, yy - 15]  # 第四象限
# ]

test_pts = []
dist = 20
for i in range(0, 360, 10):
    x = int(dist * math.cos(i / 180 * math.pi))
    y = int(dist * math.sin(i / 180 * math.pi))
    test_pts.append([x + xx, y + yy])


grid_idx = []
move_dis = 20
for pt in test_pts:
    plt.plot([xx + move_dis, pt[0] + move_dis], [yy + move_dis, pt[1] + move_dis])
    res = BresenhamAlgorithm(xx, yy, pt[0], pt[1])
    for a in res:
        # a[0] += xx
        # a[1] += yy
        grid_map[a[1] + move_dis][a[0] + move_dis] = 5
        grid_idx.append(a)

print(grid_idx)

plt.imshow(grid_map, cmap=plt.cm.hot, interpolation='nearest', vmin=0, vmax=10)
plt.colorbar()
plt.xlim(0, 20)
plt.ylim(0, 20)
my_x_ticks = np.arange(0, range_max, 1)
my_y_ticks = np.arange(0, range_max, 1)
plt.xticks(my_x_ticks)
plt.yticks(my_y_ticks)
plt.grid(True)
plt.show()

C++实现

/**
 * @details implement of bresenham algorithm from https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
 * @param x0
 * @param y0
 * @param x1
 * @param y1
 * @return intersection grid indexes
 */
std::vector<std::pair<int, int>> BresenhamAlgorithm(int x0, int y0, int x1, int y1)
{
    using GridIndex = std::pair<int, int>;

    // 1. process parallel situation
    if (x1 == x0 && y1 == y0)
    {
        return {GridIndex(x0, y0)};
    }
    else if (x1 == x0)
    {
        GridIndex              tmp_index;
        std::vector<GridIndex> result;
        result.reserve(static_cast<size_t>(std::abs(y1 - y0)));
        int y_min = 0;
        int y_max = 0;
        if (y0 < y1)
        {
            y_min = y0;
            y_max = y1;
        }
        else
        {
            y_min = y1 + 1;
            y_max = y0 + 1;
        }
        for (int y = y_min; y < y_max; ++y)
        {
            tmp_index.first  = x1;
            tmp_index.second = y;
            result.emplace_back(tmp_index);
        }
        return result;
    }
    else if (y1 == y0)
    {
        GridIndex              tmp_index;
        std::vector<GridIndex> result;
        result.reserve(static_cast<size_t>(std::abs(x1 - y0)));
        int x_min = 0;
        int x_max = 0;
        if (x0 < x1)
        {
            x_min = x0;
            x_max = x1;
        }
        else
        {
            x_min = x1 + 1;
            x_max = x0 + 1;
        }
        for (int x = x_min; x < x_max; ++x)
        {
            tmp_index.first  = x;
            tmp_index.second = y1;
            result.emplace_back(tmp_index);
        }
        return result;
    }

    // situation include eight parts of quadrant(counterclockwise)
    enum class Situation : int8_t
    {
        kQuadrant1_1 = 0,
        kQuadrant1_2 = 2,
        kQuadrant2_3 = 3,
        kQuadrant2_4 = 4,
        kQuadrant3_5 = 5,
        kQuadrant3_6 = 6,
        kQuadrant4_7 = 7,
        kQuadrant4_8 = 8
    };

    // 2. get situation from grid position
    auto origin_x = x0;
    auto origin_y = y0;
    x0 -= origin_x;
    y0 -= origin_y;
    x1 -= origin_x;
    y1 -= origin_y;

    Situation situation = Situation::kQuadrant1_1;
    if (x1 > x0 && y1 > y0)
    {
        situation = ((y1 - y0) < (x1 - x0)) ? Situation::kQuadrant1_1 : Situation::kQuadrant1_2;
    }
    else if (x1 < x0 && y1 > y0)
    {
        situation = ((y1 - y0) < (x0 - x1)) ? Situation::kQuadrant2_4 : Situation::kQuadrant2_3;
    }
    else if (x1 < x0 && y1 < y0)
    {
        situation = ((y0 - y1) < (x0 - x1)) ? Situation::kQuadrant3_5 : Situation::kQuadrant3_6;
    }
    else if (x1 > x0 && y1 < y0)
    {
        situation = ((y0 - y1) < (x1 - x0)) ? Situation::kQuadrant4_8 : Situation::kQuadrant4_7;
    }

    // transform quadrant-2/3/4 to quadrant-1, or transform back
    auto TransformQuadrant = [&](int& x1, int& y1, Situation situation, bool back = false) {
        switch (situation)
        {
            case Situation::kQuadrant1_2:
                std::swap(x1, y1);
                break;
            case Situation::kQuadrant2_3:
                x1 = -x1;
                std::swap(x1, y1);
                break;
            case Situation::kQuadrant2_4:
                x1 = -x1;
                break;
            case Situation::kQuadrant3_5:
                x1 = -x1;
                y1 = -y1;
                break;
            case Situation::kQuadrant3_6:
                x1 = -x1;
                y1 = -y1;
                std::swap(x1, y1);
                break;
            case Situation::kQuadrant4_7:
                y1 = -y1;
                std::swap(x1, y1);
                break;
            case Situation::kQuadrant4_8:
                y1 = -y1;
                break;
            case Situation::kQuadrant1_1:
            default:
                break;
        }
        if (back)
        {
            if (situation == Situation::kQuadrant2_3 || situation == Situation::kQuadrant4_7)
            {
                x1 = -x1;
                y1 = -y1;
            }
        }
    };

    // 3. transform to quadrant-1_1
    TransformQuadrant(x1, y1, situation);
    int delta_x = x1 - x0;
    int delta_y = y1 - y0;
    int error   = 0;
    int y       = y0;

    // 4. compute grid in line
    GridIndex              tmp_index;
    std::vector<GridIndex> result;
    for (int x = x0; x < x1; ++x)
    {
        tmp_index.first  = x;
        tmp_index.second = y;
        result.emplace_back(tmp_index);

        error += delta_y;
        if (2 * error > delta_x)
        {
            error -= delta_x;
            y += 1;
        }
    }

    // 5. transform back to original quadrant
    for (auto& grid_index : result)
    {
        TransformQuadrant(grid_index.first, grid_index.second, situation, true);
    }
    for (auto& grid_index : result)
    {
        grid_index.first += origin_x;
        grid_index.second += origin_y;
    }
    return result;
}

TEST(test_bresenham, test_bresenham)
{
	std::vector<std::pair<int, int>> ground_truth
		{{ -9, 0 }, { -8, 0 }, { -7, 0 }, { -6, 0 }, { -5, 0 }, { -4, 0 }, { -3, 0 }, { -2, 0 }, { -1, 0 }, { 0, 0 },
		 { 0, 0 }, { 1, 0 }, { 2, 0 }, { 3, 0 }, { 4, 0 }, { 5, 0 }, { 6, 0 }, { 7, 0 }, { 8, 0 }, { 9, 0 }, { 0, -9 },
		 { 0, -8 }, { 0, -7 }, { 0, -6 }, { 0, -5 }, { 0, -4 }, { 0, -3 }, { 0, -2 }, { 0, -1 }, { 0, 0 }, { 0, 0 },
		 { 0, 1 }, { 0, 2 }, { 0, 3 }, { 0, 4 }, { 0, 5 }, { 0, 6 }, { 0, 7 }, { 0, 8 }, { 0, 9 }, { 0, 0 }, { 1, 1 },
		 { 2, 1 }, { 3, 2 }, { 4, 3 }, { 5, 3 }, { 6, 4 }, { 7, 5 }, { 8, 5 }, { 9, 6 }, { 10, 7 }, { 11, 7 },
		 { 12, 8 }, { 13, 9 }, { 14, 9 }, { 0, 0 }, { 1, 1 }, { 1, 2 }, { 2, 3 }, { 3, 4 }, { 3, 5 }, { 4, 6 },
		 { 5, 7 }, { 5, 8 }, { 6, 9 }, { 7, 10 }, { 7, 11 }, { 8, 12 }, { 9, 13 }, { 9, 14 }, { 0, 0 }, { -1, 1 },
		 { -2, 1 }, { -3, 2 }, { -4, 3 }, { -5, 3 }, { -6, 4 }, { -7, 5 }, { -8, 5 }, { -9, 6 }, { -10, 7 }, { -11, 7 },
		 { -12, 8 }, { -13, 9 }, { -14, 9 }, { 0, 0 }, { -1, 1 }, { -1, 2 }, { -2, 3 }, { -3, 4 }, { -3, 5 }, { -4, 6 },
		 { -5, 7 }, { -5, 8 }, { -6, 9 }, { -7, 10 }, { -7, 11 }, { -8, 12 }, { -9, 13 }, { -9, 14 }, { 0, 0 },
		 { -1, -1 }, { -2, -1 }, { -3, -2 }, { -4, -3 }, { -5, -3 }, { -6, -4 }, { -7, -5 }, { -8, -5 }, { -9, -6 },
		 { -10, -7 }, { -11, -7 }, { -12, -8 }, { -13, -9 }, { -14, -9 }, { 0, 0 }, { -1, -1 }, { -1, -2 }, { -2, -3 },
		 { -3, -4 }, { -3, -5 }, { -4, -6 }, { -5, -7 }, { -5, -8 }, { -6, -9 }, { -7, -10 }, { -7, -11 }, { -8, -12 },
		 { -9, -13 }, { -9, -14 }, { 0, 0 }, { 1, -1 }, { 2, -1 }, { 3, -2 }, { 4, -3 }, { 5, -3 }, { 6, -4 },
		 { 7, -5 }, { 8, -5 }, { 9, -6 }, { 10, -7 }, { 11, -7 }, { 12, -8 }, { 13, -9 }, { 14, -9 }, { 0, 0 },
		 { 1, -1 }, { 1, -2 }, { 2, -3 }, { 3, -4 }, { 3, -5 }, { 4, -6 }, { 5, -7 }, { 5, -8 }, { 6, -9 }, { 7, -10 },
		 { 7, -11 }, { 8, -12 }, { 9, -13 }, { 9, -14 }};

	std::vector<std::pair<int, int>> test_pts{
		{ -10, 0 }, { 10, 0 }, { 0, -10 }, { 0, 10 },
		{ +15, +10 }, { +10, +15 },
		{ -15, +10 }, { -10, +15 },
		{ -15, -10 }, { -10, -15 },
		{ +15, -10 }, { +10, -15 }
	};

	std::vector<std::pair<int, int>> results;
	for (auto& pt : test_pts)
	{
		auto res = BresenhamAlgorithm(0, 0, pt.first, pt.second);
		results.insert(results.end(), res.begin(), res.end());
	}
	EXPECT_EQ(results, ground_truth);
	// for (auto& pt : results)
	// {
	// 	printf("[%d, %d], ", pt.first, pt.second);
	// }
	// printf("\n");
}
int main(int argc, char** argv)
{
	::testing::InitGoogleTest(&argc, argv);
	return RUN_ALL_TESTS();

	return 0;
}

cmake_minimum_required(VERSION 3.21)
project(test_bresenham)

set(CMAKE_CXX_STANDARD 14)

add_executable(test_bresenham main.cpp)
target_link_libraries(test_bresenham gtest pthread)

参考

https://www.bilibili.com/video/BV1364y1d7Lo/?vd_source=9408e8cab54b943547dbc522a9112342
https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
https://www.cs.montana.edu/courses/spring2009/425/dslectures/Bresenham.pdf