辫状河流拓扑结构提取及分析

河流相关交流群

拓扑结构提取

卷积特征值

利用卷积值判断节点的特征 (10.1109/CIBIM.2011.5949223).
对于3x3卷积核来说，一个八位的二进制数就可以表达一种pattern，一共256种pattern（特征值256-511）。虽然比较笨，但是非常实用. python的cv2库可以直接算卷积

单点节点的特征值

• 单连接的端点：[257, 258, 260, 264, 272, 288, 320, 384]
• 二连接：2个4联通[296,266,386,416]，2个8联通[276,261,321,336]，一个4联通一个8联通[265,274,289,292,322,328,388,400]
• 三连接，边缘的点不相接：[277,293,297,298,325,329,330,337,338,340,394,402,404,418,420,424]
• 四连接，连接之间不互通：[426, 341]
• 四连接：[331, 333, 458, 421, 422, 299, 301, 425, 427, 428, 429, 430, 357, 361, 362, 363, 365, 490, 405, 406, 410, 339, 342, 343, 466, 468, 469, 470, 345, 346, 347, 349, 474, 309, 434, 436, 437, 438, 442, 373]
• 
• 四点节点的特征值
  • 即D8里面有三个0区域
  • [387, 391, 270, 271, 395, 398, 399, 399, 451, 455, 334, 335, 459, 462, 463, 463, 419, 423, 302, 303, 427, 430, 431, 431, 480, 481, 482, 483, 483, 484, 485, 486, 487, 487, 366, 367, 488, 489, 490, 491, 491, 492, 493, 494, 494, 495, 495, 495, 403, 407, 286, 287, 411, 414, 415, 415, 467, 471, 350, 351, 475, 478, 479, 479, 435, 439, 312, 313, 314, 315, 316, 317, 318, 318, 319, 319, 440, 441, 442, 443, 443, 444, 445, 446, 446, 447, 447, 447, 496, 497, 498, 499, 499, 500, 501, 502, 503, 503, 376, 377, 378, 379, 380, 381, 382, 382, 383, 383, 504, 504, 505, 505, 506, 506, 507, 507, 507, 508, 508, 509, 509, 510, 510, 510, 511, 511, 511, 511]
    • 
      在计算卷积之前，首先把多余的点去除，特征值[263, 387, 391, 270, 271, 399, 449, 451, 455, 463, 480, 481, 483, 487, 495, 284, 286, 287, 415, 479, 312, 316, 318, 319, 447, 368, 496, 497, 499, 503, 376, 380, 382, 383, 504, 505, 507, 508, 509, 510, 511]，即D8只有一类且点数>=3的情况。

def get_bnr(num):# 二进制数字列表转换
	mat = [int(k) for k in list(bin(num)[2:])]
	if len(mat) < 8: mat = [0] * (8 - len(mat)) + mat
	return mat

if input('测试卷积特征值'):#test
	fs=[]
	qs=[]
	end_poins=[257, 258, 260, 264, 272, 288, 320, 384]
	extra_pnts=[]
	for i in range(256):
		mat = get_bnr(num)
		mat.insert(4,1)
		mat=np.array(mat)

		f=np.multiply(np.array([16,32,64,8,256,128,4,2,1]),mat).sum()
		pattern = mat.reshape((3,3))
		f_label, num_features = ndimage.measurements.label( pattern-np.array([[0,0,0],[0,1,0],[0,0,0]]), structure=[[0,1,0],[1,0,1],[0,1,0]])

		selems=[np.array([[1,1,0],[1,1,0],[0,0,0]]),np.array([[0,1,1],[0,1,1],[0,0,0]]),np.array([[0,0,0],[1,1,0],[1,1,0]]),np.array([[0,0,0],[0,1,1],[0,1,1]])]
		for selem in selems:
			if np.sum(selem*pattern)==4:qs.append(f)

		if num_features>=3 and np.sum(f_label>0)>3:fs.append(f)
		if num_features==1 and np.sum(f_label>0)>=3:extra_pnts.append(f)

	print(extra_pnts)

RivMACN

同样概念的硬算卷积，matlab，RivMACN
bd1=data_binary(x-1,y)*32+data_binary(x-1,y+1)*64+data_binary(x,y+1)*128+data_binary(x+1,y+1)*1+data_binary(x+1,y)*2+data_binary(x+1,y-1)*4+data_binary(x,y-1)*8+data_binary(x-1,y-1)*16+data_binary(x,y)*256

Lowpath

这个工具根据算法，使用DEM来计算辫状河汊道
(https://github.com/tue-alga/ttga； https://zenodo.org/record/3518174#.Y4xQ8-gzbIV)

软件安装

1. 它是c语言和cmake写的，首先安装cmake，然后用cmake编译
   cmd中cd到ttga源代码主文件夹下

mkdir build
cd build
cmake ..
make

得到sln文件

2. sln文件要用VS2019打开，右键生成exe，在.\build\src\gui\Release或者.\build\src\gui\Debug下面

3. 运行的时候报错
   This application failed to start because it could not find or load the Qt platform plugin “windows” in “”.
   是因为没有装对应版本的qt，参考这篇文章
   ，完了再重新生成exe，还是会报同样的错，是因为缺少组件，按https://blog.csdn.net/qq_35488967/article/details/78504392解决，把bin里面所以带lib的dll都复制过去。

4. 输出的是txt，写脚本转换为npy

边界控制

'''lowpath的鸡肋功能
#boundary控制计算边界，很鸡肋的功能，还不能识别南北流向的河流边界，我写了方法判断流向，把南北向的影像旋转90°来适应这个边界控制的方式。
2
2
2
2
#source,C,R
1600
270
1600
271
#sink
1601
270
1601
271
#top
1600
270
1601
270
#bottm
1600
271
1601
271
'''

使用

使用软件算出来的links是txt文件，写代码把它转tif文件

import numpy as np
import pyris

def lines2mask(lines,size):
	#lines=[[(Row_idx),(Col_idx)],[],[],...]
	mask=np.zeros(size,dtype=int)
	for line in lines:
		Row_idx,Col_idx=line
		mask[Row_idx,Col_idx]+=1
	return mask>0

def read_from_links(txt_file):
	R_C=[]
	f=open(txt_file,'r')
	for idx,line in enumerate(f.readlines()):
		if idx==0:continue
		line=[int(i) for i in line.split(' ')[2:]]
		Row_idx,Col_idx=line[1::2],line[0::2]
		R_C.append([Row_idx,Col_idx])
	return R_C

R_C=read_from_links(r'C:\Users\23932\Desktop\TongTian\links.txt')
print(len(R_C))

mask=lines2mask(R_C[:2000],(777,2820))
proj,geotrans=pyris.load_Proj(r'C:\Users\23932\Desktop\TongTian\prj\2000_244.prj'),pyris.load_Trans(r'C:\Users\23932\Desktop\TongTian\geo\2000_244.p')
pyris.array2raster(mask[:,:,np.newaxis],r'C:\Users\23932\Desktop\TongTian\topo.tif',[0],proj=proj,GeoTransf=geotrans)
pyris.show_npy(mask)

RivMACNet

基于水体的Segementation提取的拓扑结构（10.1016/j.cageo.2022.105180）.
L789: 添加了一行 xx=double(uint64(xx));yy=double(uint64(yy));
if xx1 || yy1
continue
end