文章目录

• 引言
• 1算法
• 1.1矢量点的栅格化
• 1.1.1思想
• 1.1.2算法实现

• 1.2.1八方向栅格化
• 思想
• 算法流程
• 算法实现

• 思想
• 算法流程
• 算法实现

• 2.1八方向栅格化结果
• 2.2全路径栅格化结果

引言

矢量数据结构和栅格数据结构的相互转换，是地理信息系统的常用功能。其中，矢量线转栅格的两种算法较为基础且常用，八方向栅格化和全方向栅格化。本文介绍这两种算法的思想、流程和使用C#语言实现。

1算法

1.1矢量点的栅格化

在介绍线的栅格化之前，先看点，因为点的栅格化是线的栅格化的基础。

1.1.1思想

如下图1所示，x-o-y轴是矢量点所在的平面坐标。O’（x0,y0）是栅格像元行列号计算的起点，从O‘开始向X正方向增加列号，向Y负方向行号增加。因此我们可以根据这个关系建立其栅格行列号与平面坐标系的数学关系，公式如1.1所示。
{I=1+[y0−yDy]J=1+[x−x0Dx]\\left\\{\\begin{array}{l}I=1+\\left[\\frac{y_{0}-y}{D_{y}}\\right] \\\\J=1+\\left[\\frac{x-x_{0}}{D_{x}}\\right]\\end{array}\\right.⎩⎨⎧​I=1+[Dy​y0​−y​]J=1+[Dx​x−x0​​]​
式中，Dx Dy为栅格的长、框，x y为矢量点的二维坐标，x0 y0为栅格的起点坐标；[] 为向下取整。

​ 图1

1.1.2算法实现

实现函数

static RowColumnPair PointToRaster(Point2d P0, Point2d P, double d){int j = (int)(1 + Math.Floor(Math.Abs(P.X - P0.X) / d));int i = (int)(1 + Math.Floor(Math.Abs(P0.Y - P.Y) / d));return new RowColumnPair(i, j);}

其中, RowColumnPair是自定义类型，用于存储一个点的行号和列号的数据结构。其实也可以用长度为2的int[]数组替代，int[0]代表行/列，int[1]代表列/行。取决于个人编程习惯，我比较喜欢面向对象的方法自定义类型，因为这样拓展性比较好。另外还有自定义类型Point2d，用于存储双精度浮点型的二维点坐标。

RowColumnPair类

/// <summary>/// 行、列数据类型/// </summary>public class RowColumnPair{/// <summary>/// 行/// </summary>public int Row { set; get; }/// <summary>/// 列/// </summary>public int Column { set; get; }public RowColumnPair() { }public RowColumnPair(int r, int c){this.Row = r;this.Column = c;}}

Point2d类

/// <summary>/// 点类/// </summary>public class Point2d{public double X { set; get; }public double Y { set; get; }public Point2d() { }public Point2d(double x, double y){this.X = x;this.Y = y;}public Point2d(Point2d p){this.X = p.X;this.Y = p.Y;C#}}

1.2矢量线的栅格化

1.2.1八方向栅格化

思想

根据矢量线的倾角情况，栅格的每行或每列只要一个像元被填充。

算法流程

预定义：线的端点为A、B。为方便行、列的遍历，我们使行、列的其中一个使按正方向遍历。在这里我选择列。即使A点确保在B点的左边，若不是，则交换一下。

步骤

（1）判断A、B两点位置，若A.X>B.X，则A在B的右边，需要交换。

（2）先将端点坐标栅格化，计算出行、列。存入栅格点集合中。

（3）根据两对行号、列号，判断 :

​ a.线是向上或向下【这个很重要。如果是向下，垂直于Y轴的直线就要沿着负方向动，计算其于矢量线的交线；如果是向上，垂直于Y轴的直线就要沿着正方向动，计算其于矢量线的交线。当然，如果是列差更大的情况就没有这个考虑，因为列差大是按垂直于X轴的直线运动的。】。向上则记录符号 1，反之记录符号2.

​ b.行差、和列差哪个大。

（4）若行差大于列差，则逐行计算每行的中心线，求出本行的中心线于矢量线的交点；
y=y中心线x=(y−y1)b+x1\\begin{array}{l}y=y _{ 中心线 } \\\\x=\\left(y-y_{1}\\right) b+x_{1}\\end{array}y=y中心线​x=(y−y1​)b+x1​​
​ 若列差大于行差，则逐行计算每列的中心线，求出本列的中心线于矢量线的交点。
x˙=x中心线y=(x−x1)b′+y1\\begin{array}{l}\\dot{x}=x_{\\text {中心线 }} \\\\y=\\left(x-x_{1}\\right) b^{\\prime}+y_{1}\\end{array}x˙=x中心线​y=(x−x1​)b′+y1​​
​ 该问题，我们可以简化为，已知y求已知直线上x；已知x求已知直线上y的过程。

（5）上一步计算的二维点坐标，使用点转栅格方法进行栅格化存入集合中。

（6）结束

算法实现

public   static List<RowColumnPair> VectorToRaster_EightDirection(Point2d p0, Line line, double d){List<RowColumnPair> lists = new List<RowColumnPair>();//为方便计算，让Line的A点一定在B的左边if (line.A.X > line.B.X){Point2d temp = new Point2d(line.B);line.B = new Point2d(line.A);line.A = new Point2d(temp);}//端点的栅格化var p1 = PointToRaster(p0, line.A, d);var p2 = PointToRaster(p0, line.B, d);lists.Add(p1);lists.Add(p2);//线段穿越列数int columns = p2.Column - p1.Column;columns += 1;//线段穿越行数int rows = p1.Row - p2.Row;//记录A走向B时，行（Y）的递增方向有：正方向、负方向。int fh = 1;if (rows < 0){rows = -rows;fh = -1;}rows += 1;//如果穿越的列大于穿越的行，则计算线段经过的每个栅格中线的X方向与线段的交点。if (columns > rows){double x = (p1.Column - 0.5) * d + p0.X;for (int i = 0; i < columns - 2; i++){x += d;double y = line.GetYByX(x);var row_column = PointToRaster(p0, new Point2d(x, y), d);lists.Add(row_column);}}else{double y = -(p1.Row - 0.5) * d + p0.Y;for (int i = 0; i < rows - 2; i++){y += d * fh;double x = line.GetXByY(y);var row_column = PointToRaster(p0, new Point2d(x, y), d);lists.Add(row_column);}}return lists;}

​ 其中，自定义数据类型Line ，用于存储两个端点和相关辅助计算方法。

Line类型

/// <summary>/// 线段/// </summary>public class Line{public Point2d A { set; get; }public Point2d B { set; get; }public Line() { }public Line(Point2d a, Point2d b){A = new Point2d(a);B = new Point2d(b);}/// <summary>/// 线段上任意X对应的Y/// </summary>/// <param name=\"x\"></param>/// <returns></returns>public double GetYByX(double x){double y = (B.Y - A.Y) / (B.X - A.X) * (x - A.X) + A.Y;return y;}/// <summary>/// 线段上任意Y对应的X/// </summary>/// <param name=\"y\"></param>/// <returns></returns>public double GetXByY(double y){double x = (B.X - A.X) * (y - A.Y) / (B.Y - A.Y) + A.X;return x;}}

1.2.2 全路径矢量化

思想

实质是一种“分带法”，将栅格按一行或一列分带，找出每一条带上于矢量线的相交的像元，全部填充。

算法流程

（1）将端点矢量化。

（2)基于矢量线的始末点和倾角大小，计算每一行/列要填充的的首末列号/行号。同样，根据列差还是行差大也有不同的计算方法。

​ 第一种情况：|X2-X1|≥|Y2-Y1|时，计算矢量倾角K=（Y2-Y1）/(X2-X1)，为方便计算，确保A–>B的行号是增加的，否则交换A B。

​ 计算起始列号:
ja=[(y0−(i−1)m−y1k+x1−x0)/m]+1j_{a}=\\left[\\left(\\frac{y_{0}-(i-1) m-y_{1}}{k}+x_{1}-x_{0}\\right) / m\\right]+1ja​=[(ky0​−(i−1)m−y1​​+x1​−x0​)/m]+1
​ 计算终止列号：
je=[(y0−im−y1k+x1−x0)/m]+1j_{e}=\\left[\\left(\\frac{y_{0}-i m-y_{1}}{k}+x_{1}-x_{0}\\right) / m\\right]+1je​=[(ky0​−im−y1​​+x1​−x0​)/m]+1
​ i为当前行号，i从A.Row开始取，到B.Row结束；

​ m为像元宽度;(x0,y0)为栅格的起点坐标（左上角）；

​ (x1 ,y1)为A点坐标；

​ []为取整。

​ **只有第一行需要计算ja,之后的每一行的ja都为上一行的je。**每一行的ja 到je之间的所有像元全部涂黑。

​ 第二种情况：|X2-X1|<|Y2-Y1|时,计算矢量倾角K=（Y2-Y1）/(X2-X1)，为方便计算，确保A–>B的列号是增加的，否则交换A B。

​ 计算起始行号：
ia=[((x1−x0−(j−1)∗m)∗tan⁡α+y0−y1)/m]+1i_a=[((x_1-x_0-(j-1) * m) * \\tan \\alpha+y_0-y_1 )/ m]+1ia​=[((x1​−x0​−(j−1)∗m)∗tanα+y0​−y1​)/m]+1
​ 计算终止列号：
ie=[((x1−x0−j∗m)∗tan⁡α+y0−y1)/m]+1i_e=[((x_1-x_0-j * m) * \\tan \\alpha+y_0-y_1 )/ m]+1ie​=[((x1​−x0​−j∗m)∗tanα+y0​−y1​)/m]+1
​ 参数解释同第一种情况。略。

​ **只有第一行需要计算ia,之后的每一行的ia都为上一行的je。**每一行的ia 到ie之间的所有像元全部涂黑。

算法实现

/// <summary>/// 第五题：5、编程实现全路径矢量线栅格化/// </summary>/// <param name=\"p0\"></param>/// <param name=\"line\"></param>/// <param name=\"d\"></param>/// <returns></returns>public  static List<RowColumnPair> VectorToRaster_AllDirection(Point2d p0,Line line ,double d){List<RowColumnPair> lists = new List<RowColumnPair>();//端点的栅格化var p1 = PointToRaster(p0, line.A, d);var p2 = PointToRaster(p0, line.B, d);lists.Add(p1);lists.Add(p2);if ((Math.Abs(line.B.X-line.A.X)>Math.Abs(line.B.Y-line.A.Y))){//确保A-->B的行号是增加的if (line.A.Y<line.B.Y){Point2d temp = new Point2d(line.B);line.B = new Point2d(line.A);line.A = new Point2d(temp);}var startRaster = PointToRaster(p0, line.A, d);var endRaster = PointToRaster(p0, line.B, d);;                double k = (line.B.Y - line.A.Y) / (line.B.X - line.A.X);int ja = (int)Math.Floor(((p0.Y - (startRaster.Row - 1) * d - line.A.Y) / k + line.A.X - p0.X) / d )+ 1;int startJ = ja;int lastJ = 0;//上一次的je，防止因je与startJ的交换破坏了je的逻辑。for (int i = startRaster.Row; i <= endRaster.Row; i++){int je = (int)Math.Floor(((p0.Y - i * d - line.A.Y) / k + line.A.X - p0.X) / d) + 1;lastJ = je;if (je<startJ){var temp = je;je = startJ;startJ = temp;}//每行的起始列要比终止列小for (int j = startJ; j <=je; j++){RowColumnPair pair = new RowColumnPair(i, j);lists.Add(pair);}startJ = lastJ;}}else{//确保A-->B的列号是增加的if (line.A.X>line.B.X){Point2d temp = new Point2d(line.B);line.B = new Point2d(line.A);line.A = new Point2d(temp);}var startRaster = PointToRaster(p0, line.A, d);var endRaster = PointToRaster(p0, line.B, d);double k = (line.B.Y - line.A.Y) / (line.B.X - line.A.X);int ia = (int)Math.Floor(((line.A.X - p0.X - (startRaster.Column - 1) * d) * k + p0.Y - line.A.Y) / d) + 1;int startI = ia;int lastI = 0;//上一次的je，防止因ie与startI的交换破坏了ie的逻辑。for (int j = startRaster.Column; j <=endRaster.Column; j++){int ie= (int)Math.Floor(((line.A.X - p0.X - j  * d) * k + p0.Y - line.A.Y) / d) + 1;lastI = ie;if (startI>ie){var temp = startI;startI = ie;ie = temp;}//每列的起始行要比终止行小for (int i = startI; i <=ie; i++){RowColumnPair pair = new RowColumnPair(i, j);lists.Add(pair);}startI = ie;}}return lists;}

2实验结果

本来想用GDI+做可视化渲染看出结果，觉得太麻烦了。毕竟复现算法的目的是为了熟悉算法而已。于是用绘图软件花了个简单的示意图。凑活着看吧。

2.1八方向栅格化结果

测试数据1：栅格起点80, 580 栅格边长50 长宽数量10X10 线段(500, 500) (200, 300)

测试数据2：栅格起点80, 580 栅格边长50 长宽数量10X10 线段(400, 300),(200, 500)

2.2全路径栅格化结果

测试数据1：栅格起点80, 580 栅格边长50 长宽数量10X10 线段(500, 500) (200, 300)

测试数据2：栅格起点80, 580 栅格边长50 长宽数量10X10 线段(400, 300),(200, 500)

小结

如果是垂直和水平的线，可直接计算出栅格行列，加上这个判断也很简单。代码中如果加入对直线形状的判断可能会提高时间复杂度。而且本文是复现算法，算法能实现所有的情况即可，对于水平和垂直线的特殊情况自己心里知道就好。
另外，可以观察到全路径栅格化方法，在线段的两端位置会比实际线段长，导致栅格失真，在实际项目中我们需要做边界处理。