基于条件极值的四边形直角化方法研究

马鹏刚，王永红，周群强，赵宝军，张鹏岩

（1.国家测绘地理信息局第二地形测量队，陕西 西安 710054；2.河南大学 环境与规划学院，河南 开封 475004）

数字地形图是“数字城市”、“数字地球”的重要组成部分，在城市建设和规划中起着重要作用。在大比例尺地形图（1∶500、1∶1 000、1∶2 000）中，特定地物、建筑物、停车位、院落等均为重要要素，其属性的准确性直接影响用户的使用情况。航空摄影测量时，由于地形条件、投影变换、采集误差等因素，使得本来为矩形的地物要素变为一般的四边形，导致地图失真，为了能准确反映地物要素的实际情况，需要将四边形修正为矩形。然而，单纯人工修正的效率和准确度都很低，因此本文采用条件极值法，以北京地区停车位信息采集为例，利用GIS的二次开发技术编写了自动直角化软件。首先提取四边形的中点，以平行于长边的平行线为中线；再将四边形旋转，相邻两边互相垂直，中点不变，使其保持原来的中心位置。

1 四边形的产生和表现形式

一般情况下，野外测量或内业采集时得到的图形为矩形，但在进行预处理时，由于测量误差、投影转换、加密控制等因素，将导致矢量图形发生变形，甚至位移，使得矩形变成平行四边形[1-3]，如图1所示。

图1 四边形产生的原因与表现形式

2 直角化方法

伪矩形的存在，使得地图表达与实际情况存在差异，直角化时必须保证中点位置不变。目前运用最多的方法是最小二乘法[4]，其具有技术优化、方便快捷的优点，但只是将平行四边形简单直角化，车位位置却发生了偏移，不能满足采集北京地区停车位的位置和面积信息的要求。经过反复试验，将条件极值法与最小二乘法相结合，利用GIS的二次开发功能，以编程的形式开发了直角化软件。该方法确保了图形中点固定不变，长边中点以一定的角度旋转，直至四边相互垂直。

2.1 条件极值模型

在约束条件下，函数z=f(x1, x2, x3, …, xn)的阵列数为[4-5]：

令 F(x1, x2, x3, …, xn)=f(x1, x2, x3, ..., xm)+∑∮iβi(x1,x2, x3,…, xn)k∈n, ∮i∈m,分别对 x1, x2, x3, …, xn求偏导数。要使得函数z=f(x1, x2, x3,…, xn)存在极值，则其对各变量的偏导数为零。根据原始阵列和偏导数方程可求出z=f(x1, x2, x3, …, xn)极值时的x1, x2, x3, …, xn和∮1,∮2, …,∮m。

2.2 基于条件极值的中点固定模型

进行四边形直角化时，要确保其几何中心不变，反将其中线围绕中点进行旋转。这种算法类似一种近似值，即直角化后矩形的四角无限接近于90°，但不为90°，如图2a所示。首先，将原来的四边形（红色边线）沿长轴以A方向旋转[5-7]，使得长轴的方位角无限接近于零，则四边形的四角将无限接近于90°，由于条件极值的数学模型为阵列式，因此其永远不能到达设定值，即红色长轴与蓝色长轴近似重叠，但不重叠；然后，将红色长轴沿B方向旋转至蓝色长轴，使得长轴的方位角无限接近于零，则四边形的四角将无限接近于90°，由于条件极值的数学模型为阵列式，因此其永远不能到达设定值，即红色长轴与蓝色长轴近似重叠，但不重叠。两种方向所得结果一致，但旋转半径不同，因此利用数学基础，将A方向（即最短距离方向）设定为旋转长轴。

图2 旋转顺序与成果

对四边形（红色）进行直角化后，其相邻两边相互垂直，则直角对应的限制条件为：

当四边形存在非直角以及两个长轴没有重合时，进行旋转处理，将式（1）进行线性化处理可得到四角直角化后的参数。

2.3 基于GIS的程序实现过程

GIS的二次开发是面向对象的程序设计（OPP），本文以前期作业区北京地区停车位信息采集矢量数据（.mdb）为例，基于GIS的ArcObject，以Microsoft Visual Studio 2010为开发环境，以Visual Basic为开发语言进行程序开发。对象（Object）为停车位，类（Class）为四边形，继承（Inher）为可不旋转图形。

图3 四边形直角化实现过程

1）将四边形的坐标按顺时针方向排序，以中点为基点，四角分别为A、B、C、D，如图3所示。

2）为了避免单纯旋转得到如图2b所示的结果，即只进行旋转，而未进行垂直处理，需计算a角的度数，并将其与90°作减法，再按照得出的角度进行有目的的旋转，代码为：

//Console.WriteLine("Stand:{0},L1:{1},L2:{2},L3:{3},L4:{4},LS-L1:{5},LS-L2:{6},LS-L3:{7},LS-L4:{8}", lStand.Angle, l1.Angle, l2.Angle, l3.Angle, l4.Angle,

//lStand.Angle-l1.Angle, lStand.Angle-l2.Angle,lStand.Angle-l3.Angle, lStand.Angle-l4.Angle);

double AngleDiff1= Math.Abs(lStand.Angle-l1.Angle);

double AngleDiff2= Math.Abs(lStand.Angle-l2.Angle);

double AngleDiff3= Math.Abs(lStand.Angle-l3.Angle);

double AngleDiff4= Math.Abs(lStand.Angle-l4.Angle);

double ad1= Math.Min(AngleDiff1, AngleDiff3);

double ad2= Math.Min(AngleDiff2, AngleDiff4);

‘和轴线平行

double ParallelLineLength= 0;

‘和轴线垂直

double VerticalLine= 0;

if (ad1 ＜ ad2)//1-2 或 3-4

{

ParallelLineLength= (l1.Length + l3.Length) / 2;

VerticalLine= (l2.Length + l4.Length) / 2;

else//2-3或4-1

{ParallelLineLength= (l2.Length + l4.Length) / 2;

VerticalLine= (l1.Length + l3.Length) / 2;}

double Angle= lStand.Angle;

double HalfPi= Math.PI / 2;

PointClass cp= new PointClass();

‘沿轴线前进一半

cp.ConstructAngleDistance(pi.Centroid, Angle,ParallelLineLength / 2);

PointClass cp1= new PointClass();

PointClass cp2= new PointClass();

PointClass cp3= new PointClass();

PointClass cp4= new PointClass();

Angle += HalfPi;

cp1.ConstructAngleDistance(cp, Angle, VerticalLine / 2);

Angle += HalfPi;

cp2.ConstructAngleDistance(cp1, Angle, ParallelLineLength);

Angle += HalfPi;

cp3.ConstructAngleDistance(cp2, Angle, VerticalLine);

Angle += HalfPi;

cp4.ConstructAngleDistance(cp3, Angle, ParallelLineLength);

/Console.WriteLine("({0},{1})({2},{3})({4},{5})({6},{7})", cp1.X,cp1.Y,cp2.X,cp2.Y,cp3.X,cp3.Y,cp4.X,cp4.Y);

‘找到每组中心点距离最大的两个点连成直线

‘考虑只有一个图形的极端情况

List＜ParallelogramInfo＞ lstGroupParallelogramInfo=new List＜ParallelogramInfo＞();

List＜int＞ lstOID= new List＜int＞();

do

{

‘由于存在跳跃，每次都全部遍历

for (int i= 0; i ＜ lstParallelogramInfo.Count; i++)

{

ParallelogramInfo pi= lstParallelogramInfo[i];

if (lstOID.Contains(pi.OID))

continue;

Console.WriteLine(string.Format("{0}/{1}", lstOID.Count, lstParallelogramInfo.Count));

‘只有一个的情况

if (lstGroupParallelogramInfo.Count== 1)

‘找到每组中心点距离最大的两个点连成直线

‘考虑只有一个图形的极端情况

// List＜ParallelogramInfo＞ lstGroupParallelogramInfo= new List＜ParallelogramInfo＞();

// List＜int＞ lstOID= new List＜int＞();

‘由于存在跳跃，每次都全部遍历

/ for (int i= 0; i ＜ lstParallelogramInfo.Count; i++)

3）将原始四边形的中点与旋转后的四边形（此时由于还未作垂直处理，不能称其为矩形）定位重合。

4）已计算得到a角的值，根据角度值作垂直处理，并设计界面，代码为：

＜Window x:Class="Revise2Rect.FrmMain"

xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"

xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"

Title="图形校正" Height="800" Width="800" Wi ndowStartupLocation="CenterScreen" Icon="default.ico"Closing="Window_Closing"＞

＜Window.Background＞

＜LinearGradientBrush EndPoint="1,1"StartPoint="0,0"＞

＜GradientStop Color="LightBlue" Offset="0" /＞

＜GradientStop Color="LightGreen" Offset="0.5" /＞

＜GradientStop Color="#FFBE5A5A" Offset="1" /＞

＜/LinearGradientBrush＞

＜/Window.Background＞

输出结果如图4所示。

图4 输出结果

3 结 语

实验证明，本文所提出的方法对于四边形（多边形还未实现）进行直角化的效果明显，最终成果已无限接近矩形，避免了部分软件直角化时产生的四边形失真或中心位移问题。对于四边形四边较小的图形，其边长直接影响直角化效果，四边越接近，直角化效果越明显。本文设定的条件极值模型能快速、简单、精准地进行直角化处理，但目前还不适用于其他多边形的直角化。四边形以上的多边形直角化算法相对复杂一些，后续可根据项目要求继续探讨。