这是写在作业里的，所以描述偏正式= =放在博客里主要是因为这个矢量化的过程都是自己从零一步步构想，并用代码实现出来的，值得纪念一下。或许相关算法已有很多，并且比这个要优秀，但是百度了一圈实在找不到不破坏拓扑关系的线程库函数（cv2.FindContours点名批评）。
为了在一个栅格图像中提取像元值相同，图形不能有自相交的矢量要素，并且不破坏原有图像中的拓扑关系。我们小组进行了一次算法上的尝试。主要思路：首先，使用区域生长法识别出每一个图形；其次，对图形进行轮廓提取；最后，根据图形轮廓的点序列生成相应的几何要素。
对图形的识别
这里，我们使用了最简单易懂的区域生长法。由于图形不能自相交，所以我们只使用四邻域的生长方式。若是使用八邻域生长方式，可能会出现对角生长的情况，最终使得一个图形被一个交点分成多个部分。
区域生长法的基本思想就是从任意给定的像元点出发，检索四周像元的值。若符合条件，则归入与中心像元相同的一类。然后再从被归入的像元开始，继续向四周检索，直到周围再无符合条件的像元则停止。
算法实现上的思路也十分简单，首先需要定义一个寻找四周像元坐标的函数。这里的坐标，以行、列数减一表示。例如，图像左上角的像元，坐落于第一行第一列，那么其坐标为(0,0)。那么其四周点的坐标也就是分别在x坐标和y坐标上加减一了，但这里还需要注意图像边界的情况，做一个边界判断即可。这一块用getaround函数来实现，详见文末。
在获取像元四周的坐标后，需要做两个判断。一是四周的像元是否符合条件，二是四周的像元是否已被归入某个类别。第一项十分容易判断。而第二项，可以通过新建一个与栅格图像等大小的矩阵来存贮每个像元是否被拾取的状态，这里以0表示未被归类，1表示已被归类。而被归类的像元又是生长的起始像元，所以这里需要用到递归。我们使用变量shape用来存储目前类别中的像元。
那“任意给定的点”该如何选择呢？由于栅格图像中，所有的像元都需要被归类，因此，只需要不断遍历未被归类的像元坐标即可，遍历终止的条件即是所有像元均被归类。我们使用shapes来存储所有的图形与对应的像元值。
对轮廓的提取
现在，我们得到同类别的图形了，需要提取出图形的边界点。为了不破坏数据的空间关系，必须严格按照栅格像元的边界来提取，不可有位移、放缩或者简化的情况，否则会导致距离关系与拓扑关系的不正确。
首先，需要确立一个坐标系。这里以栅格图像左上角的顶点为(0,0)，水平为x轴，指向右侧，垂直为y轴，指向下侧，构建一个平面直角坐标系。之后我们所有的图形边界点的坐标都是在这个坐标系下的。
而对于边界的提取，采用围绕法。即记从图形的一点出发，绕着图形走一圈得到的路径为图形的边界。这里默认为顺时针的走法。而每步行走过程中方向需要根据当前格点四周四块像元是否在图像内作判断，分为三种情况。
1.当前格点四周栅格仅有一个为同类时，则向右拐。
2.当前格点四周栅格有三个为同类，则向左拐。
- 当前格点四周栅格有两个为同类，若同类分布为对角，则向右拐，否则走直线。
应当不存在其他的情况，若存在，则说明之前的判断出现了错误。具体的代码详见材料的其他附件。
几何要素的生成
这里使用ogr库来进行矢量图形的创建与一系列操作。对于之前得到的点序列，通过坐标转换后，直接写入ogr.wkbLinearRing再转为ogr.wkbPolygon即可完成矢量图形的创建。
代码：findedge.py
import numpy as np
import sys
sys.setrecursionlimit(4000)
#00110
#00000
#01100
#00011
#00111
#随手写的，不要吐槽变量名
def getaround(p,boardx,boardy):
res=[]
if(-1<p[0]-1<boardx):
res.append((p[0]-1,p[1]))
if(-1<p[0]+1<boardx):
res.append((p[0]+1,p[1]))
if(-1<p[1]-1<boardy):
res.append((p[0],p[1]-1))
if(-1<p[1]+1<boardy):
res.append((p[0],p[1]+1))
return res
def insp(sp,ep,gs):
sr=[(ep[1],ep[0]),(ep[1]-1,ep[0]),(ep[1]-1,ep[0]-1),(ep[1],ep[0]-1)]
tmpsr=[]
case=0
for i in sr:
if i in gs:
case+=1
tmpsr.append(1)
else:
tmpsr.append(0)
case2=0
se=(ep[0]-sp[0],ep[1]-sp[1])
if(se[0]==0):
case2=1
if(case==1):
if(case2==1):
er=(se[1],se[0])
else:
er=(-se[1],-se[0])
p_res=(ep[0]+er[0],ep[1]+er[1])
if(case==2):
if(tmpsr==[0,1,0,1] or tmpsr==[1,0,1,0]):
if(case2==1):
er=(se[1],se[0])
else:
er=(-se[1],-se[0])
p_res=(ep[0]+er[0],ep[1]+er[1])
else:
p_res=(ep[0]+se[0],ep[1]+se[1])
if(case==3):
if(case2==1):
er=(-se[1],-se[0])
else:
er=(se[1],se[0])
p_res=(ep[0]+er[0],ep[1]+er[1])
#print(case,sp,ep,sr)
return p_res
class FindEdge:
def __init__(self,ls):
self.ls=ls
self.bdx=len(ls)
self.bdy=len(ls[0])
#self.uniqV=np.unique(ls).tolist()
self.lsTa=np.zeros((self.bdx,self.bdy))
self.shapes=[]
self.frms=[]
def f(self,p,dm,shape,boardx,boardy):
zhouwei=getaround(p,boardx,boardy)
for i in zhouwei:
if (self.lsTa[i[0]][i[1]]==0 and self.ls[i[0]][i[1]]==dm):
self.lsTa[i[0]][i[1]]=1
shape.append(i)
self.f(i,dm,shape,boardx,boardy)
return
def calshape(self):
while((self.lsTa != np.ones((self.bdy,self.bdx)))) :
for i in range(self.bdx):
for j in range(self.bdy):
if(self.lsTa[i][j]==0):
shape=[]
self.f((i,j),self.ls[i][j],shape,self.bdx,self.bdy)
self.shapes.append((shape,self.ls[i][j]))
#print(self.shapes)
def caledge(self):
for i in self.shapes:
i[0].sort()
frm=[]
startp=(i[0][0][1],i[0][0][0])
secp=(i[0][0][1],i[0][0][0]+1)
frm.append(startp)
frm.append(secp)
while secp != startp:
secp=insp(frm[-2],frm[-1],i[0])
frm.append(secp)
self.frms.append((frm,i[1]))
return self.frms
if __name__=="__main__":
ls=[[0,0,1,1,0],[0,0,0,0,0],[0,1,1,0,0],[0,0,0,1,1],[0,0,1,1,1]]
ls=np.ones([62,65], dtype = int)
_fe=FindEdge(ls)
_fe.calshape()
edges=_fe.caledge()
for i in edges:
print(i)
