若干计算机图形学算法实现

好久没更新了~

感觉研究生生活过的很快又很慢。快是因为一年时间转瞬即逝，基本上没有什么产出；慢是因为每天精力被一些琐碎的事情占据，一直在获得负反馈。

今天更新一下博客咯~（以下程序均是个人参考一些网络资料手动实现，经过若干千余行代码工程实践验证。由于此前网络资料参考并未记录，且多数质量欠佳，便不再附上。）

下面C++代码会涉及到一些基础数据结构的实现，在这里一并给出：

// ###################################################
// template point
// ###################################################
template <typename T>
struct TemplatePoint
{
    T x, y;
    TemplatePoint(T x_, T y_) : x(x_), y(y_) {}
    TemplatePoint() : x(0), y(0) {}
    friend std::ostream &operator<<(std::ostream &o, const TemplatePoint<T> &tp)
    {
        o << "(" << tp.x << ", " << tp.y << ")";
        return o;
    }
};
template <typename T>
bool operator==(const struct TemplatePoint<T> &l, const struct TemplatePoint<T> &r)
{
    return (l.x == r.x) && (l.y == r.y);
}
typedef struct TemplatePoint<double> dPoint;
typedef struct TemplatePoint<int> iPoint;

// scanline
typedef struct edge
{
    int ymax;
    double x;
    double dx;
    struct edge *next;
} Edge;
typedef Edge **EdgeTable;
typedef Edge *ActiveEdgeTable;

在Grid world中画一条直线：

// grid traversal
void gridTraversal(const dPoint &start, const dPoint &goal, const double resolution, std::vector<iPoint> &visited_grid)
{
    iPoint s_grid = {static_cast<int>(std::floor(start.x / resolution)), static_cast<int>(std::floor(start.y / resolution))};
    iPoint g_grid = {static_cast<int>(std::floor(goal.x / resolution)), static_cast<int>(std::floor(goal.y / resolution))};
    dPoint vector = {goal.x - start.x, goal.y - start.y};
    double stepX = (vector.x > 0) ? 1 : -1;
    double stepY = (vector.y > 0) ? 1 : -1;
    double next_grid_boundary_x = (s_grid.x + stepX) * resolution;
    double next_grid_boundary_y = (s_grid.y + stepY) * resolution;
    double tMaxX = (vector.x != 0) ? (next_grid_boundary_x - start.x) / vector.x : DBL_MAX;
    double tMaxY = (vector.y != 0) ? (next_grid_boundary_y - start.y) / vector.y : DBL_MAX;
    double tDeltaX = (vector.x != 0) ? resolution / vector.x * stepX : DBL_MAX;
    double tDeltaY = (vector.y != 0) ? resolution / vector.y * stepY : DBL_MAX;
    iPoint diff = {0, 0};
    iPoint c_grid = {s_grid.x, s_grid.y};
    visited_grid.push_back(c_grid);
    bool negative = false;
    if (s_grid.x != g_grid.x && vector.x < 0)
    {
        diff.x--, negative = true;
    }
    if (s_grid.y != g_grid.y && vector.y < 0)
    {
        diff.y--, negative = true;
    }
    if (negative)
    {
        c_grid.x += diff.x;
        c_grid.y += diff.y;
        visited_grid.push_back(c_grid);
    }
    double tx = tMaxX;
    double ty = tMaxY;
    while (!(c_grid == g_grid))
    {
        if (tx < ty)
        {
            c_grid.x += stepX;
            tx += tDeltaX;
        }
        else
        {
            c_grid.y += stepY;
            ty += tDeltaY;
        }
        visited_grid.push_back(c_grid);
    }
}

在Grid world里填充一个多边形（可凸可凹）：

void fillPolygon(const std::vector<Point> &vertices, double resolution, std::vector<intPoint> &polygArea)
{
    if (vertices.size() < 3)
        return;
    std::vector<intPoint> intVertices;
    for (int i = 0; i < vertices.size(); i++)
    {
        intPoint intVertex = {(int)(vertices[i].x / resolution), (int)(vertices[i].y / resolution)};
        intVertices.push_back(intVertex);
    }
    std::vector<int> tempY;
    int minY, maxY, scanLineNums;
    for (int i = 0; i < intVertices.size(); i++)
    {
        tempY.push_back(intVertices[i].y);
    }
    minY = *std::min_element(tempY.begin(), tempY.end());
    maxY = *std::max_element(tempY.begin(), tempY.end());
    scanLineNums = maxY - minY + 1;
    /*建立边表*/
    EdgeTable ET = new Edge *[scanLineNums];
    for (int y = 0; y < scanLineNums; y++)
    {
        ET[y] = new Edge;
        ET[y]->next = NULL;
    }
    intPoint predP, currP1, currP2, nextP;
    /*初始化边表*/
    for (int i = 0; i < intVertices.size(); i++)
    {
        currP1.y = intVertices[i].y;
        currP2.y = intVertices[(i + 1) % intVertices.size()].y;
        if (currP1.y == currP2.y) // 舍弃平行X轴的边
            continue;
        currP1.x = intVertices[i].x;
        currP2.x = intVertices[(i + 1) % intVertices.size()].x;
        predP.x = intVertices[(i - 1 + intVertices.size()) % intVertices.size()].x;
        predP.y = intVertices[(i - 1 + intVertices.size()) % intVertices.size()].y;
        nextP.x = intVertices[(i + 2) % intVertices.size()].x;
        nextP.y = intVertices[(i + 2) % intVertices.size()].y;
        int ymin = std::min(currP1.y, currP2.y);
        int ymax = std::max(currP1.y, currP2.y);
        double x = currP1.y > currP2.y ? currP2.x : currP1.x;
        double dx = (double)(currP1.x - currP2.x) / (double)(currP1.y - currP2.y);
        if (((currP2.y >= currP1.y) && (currP1.y >= predP.y)) || ((currP1.y >= currP2.y) && (currP2.y >= nextP.y)))
        {
            ymin++;
            x += dx;
        }
        Edge *tempE = new Edge;
        tempE->ymax = ymax;
        tempE->x = x;
        tempE->dx = dx;
        tempE->next = ET[ymin - minY]->next;
        ET[ymin - minY]->next = tempE;
    }
    /*建立活动边表*/
    ActiveEdgeTable AET = new Edge;
    AET->next = NULL;
    /*扫描线扫描*/
    for (int y = minY; y < maxY + 1; y++)
    {
        /*取出ET中当前扫描行的所有边并按x的递增顺序（若x相等则按dx的递增顺序）插入AET*/
        while (ET[y - minY]->next)
        {
            Edge *tempE = ET[y - minY]->next;
            Edge *tempAET = AET;
            while (tempAET->next)
            {
                if ((tempE->x > tempAET->next->x) || ((tempE->x == tempAET->next->x) && (tempE->dx > tempAET->next->dx)))
                {
                    tempAET = tempAET->next;
                    continue;
                }
                break;
            }
            ET[y - minY]->next = tempE->next;
            tempE->next = tempAET->next;
            tempAET->next = tempE;
        }
        /*将AET中的边两两配对并将中间点添加到polygArea中*/
        Edge *tempAET = AET;
        while (tempAET->next && tempAET->next->next)
        {
            for (int x = tempAET->next->x; x < tempAET->next->next->x; x++)
            {
                intPoint occ = {x, y};
                polygArea.push_back(occ);
            }
            tempAET = tempAET->next->next;
        }
        /*删除AET中满足y=ymax的边*/
        tempAET = AET;
        while (tempAET->next)
        {
            if (tempAET->next->ymax == y)
            {
                Edge *tempE = tempAET->next;
                tempAET->next = tempE->next;
                tempE->next = NULL;
                delete tempE;
            }
            else
            {
                tempAET = tempAET->next;
            }
        }
        /*更新AET中边的x值，进入下一循环*/
        tempAET = AET;
        while (tempAET->next)
        {
            tempAET->next->x += tempAET->next->dx;
            tempAET = tempAET->next;
        }
    }
    /*释放边表内存*/
    for (int y = 0; y < scanLineNums; y++)
    {
        while (ET[y]->next)
        {
            Edge *next = ET[y]->next->next;
            delete ET[y]->next;
            ET[y]->next = next;
        }
        delete ET[y];
    }
    delete[] ET;
    /*释放活动边表内存*/
    while (AET)
    {
        Edge *next = AET->next;
        delete AET;
        AET = next;
    }
}

 这一部分代码由Python实现，且较为简单，便不附上C++实现，需要时可以自己参考手动实现。

在Grid world 中画一个圆：

def circle_kernel_generation(
    circle_radius: float, grid_resolution: float, occ_value: int = 1
):
    # 中点圆算法
    r = np.ceil(circle_radius / grid_resolution).astype(int)  # 向上取整
    kernel = np.zeros((2 * r + 1, 2 * r + 1)).astype(int)  # 初始化kernel
    d, x, y = 1 - r, 0, r
    while x <= y:
        kernel[r + x, r + y] = occ_value
        kernel[r + y, r + x] = occ_value
        kernel[r + x, r - y] = occ_value
        kernel[r + y, r - x] = occ_value
        kernel[r - x, r + y] = occ_value
        kernel[r - y, r + x] = occ_value
        kernel[r - x, r - y] = occ_value
        kernel[r - y, r - x] = occ_value
        if d < 0:
            d += 2 * x + 3
        else:
            d += 2 * (x - y) + 5
            y -= 1
        x += 1
    return kernel

在Grid world中填充一个圆：

def disk_kernel_generation(
    disk_radius: float, grid_resolution: float, occ_value: int = 1
):
    # 这部分代码参考论文：Fast collision checking for intelligent vehicle motion planning
    r = np.ceil(disk_radius / grid_resolution).astype(int)  # 向上取整
    kernel = np.zeros((2 * r + 1, 2 * r + 1)).astype(int)  # 初始化kernel
    f, x, y = 1 - r, r, 0
    delta_up_left, delta_up = -2 * r, 1
    while x >= y:  # 原始论文错误
        if f > 0:
            kernel[r - x : r + x + 1, r - y : r + y + 1] = occ_value
            kernel[r - y : r + y + 1, r - x : r + x + 1] = occ_value
            x -= 1
            delta_up_left += 2
            f += delta_up_left
        y += 1
        delta_up += 2
        f += delta_up
    return kernel

以上就是本篇博客全部内容咯~

希望自己能早日做上真正的科研工作，有自己的成果产出；也希望一年后申博、套磁顺利呀~ 

Shaw Shaw加油！往前看，别回头！