python实现可视域算法

#-*-coding:utf8-*-
“““
written: 2014.04.02; by guoning; vision: 0.1
function： 比较三种算法的效率
“““

import gdal
import time
import math
import numpy


#边界扫描算法（基于bresenham近似）
def Los_sweep（arearadius）:
    blockvisib = [[0 for i in range（radius+1）]for i in range（radius+1）]  #可视性存储矩阵（r+1）*（r+1）
    blockvisib[0][0] = 1
    #area = band.ReadAsArray（int（x0）int（y0）int（radius+1）int（radius+1））
    area = numpy.transpose（area）  #转置，使x，y符合使用习惯
    viewheight = area[0][0]
    mslope = -65535
    x = 0
    y = 0
    dx = radius
    for i in range（radius+1）:
        dy = i
        dp = 2*dy-dx
        while（x < dx）:
            x+=1
            if（dp < 0）:
                dp+=2*dy
            else:
			    y+=1
			    dp+=2*（dy-dx）
			
            tslope = （area[x][y] - viewheight）/math.sqrt（x*x+y*y）
            if（tslope > mslope）:
		        blockvisib[y][x] = 1
		        mslope = tslope

    #将两矩阵转置进行相同的运算	
	area = numpy.transpose（area）
	blockvisib = numpy.transpose（blockvisib）
	mslope = -65535
	x = 0
	y = 0
    dx = radius
    for i in range（radius+1）:
        dy = i
        dp = 2*dy-dx
        while（x < dx）:
            x+=1
            if（dp < 0）:
                dp+=2*dy
            else:
			    y+=1
			    dp+=2*（dy-dx）
			    
            tslope = （area[x][y] - viewheight）/math.sqrt（x*x+y*y）
            if（tslope > mslope）:
		        blockvisib[y][x] = 1
		        mslope = tslope
	#得到完整区块可视性	        
	blockvisib = numpy.transpose（blockvisib）
	return blockvisib


#xdraw算法
def XDraw（arearadius）:
    blockvisib = [[0 for i in range（radius+1）]for i in range（radius+1）]  #可视性存储矩阵（r+1）*（r+1）
    blockvisib[0][0] = 1
    
    #area = band.ReadAsArray（int（x0）int（y0）int（radius+1）int（radius+1））
    area = numpy.transpose（area）  #转置，使x，y符合使用习惯
    lowestvisibalh = [[0 for i in range（radius+1）]for i in range（radius+1）]  #初始化最低可视海拔存储矩阵为原海拔矩阵
    lowestvisibalh[0][0] = area[0][0]
    viewheight = area[0][0]  #视点高度
    
    #判断x轴方向可视性并存入可视矩阵，更新最低可视海拔矩阵
    heightlist1 = [area[i+1][0] for i in range（radius）]
    mslope = -65535
    for i in range（radius）:
        tslope = float（heightlist1[i] - viewheight）/（i+1）
        if（tslope > mslope）:
            blockvisib[i+1][0] = 1
            mslope = tslope
            lowestvisibalh[i+1][0] = area[i+1][0]
        else:
            lowestvisibalh[i+1][0] = lowestvisibalh[i][0]+mslope
            
    #判断y轴方向可视性并存入可视矩阵，更新最低可视海拔矩阵
    heightlist2 = [area[0][i+1] for i in range（radius）]
    mslope = -65535
    for i in range（radius）:
        tslope = float（heightlist2[i] - viewheight）/（i+1）
        if（tslope > mslope）:
            blockvisib[0][i+1] = 1
            mslope = tslope
            lowestvisibalh[0][i+1] = area[0][i+1]
        else:
            lowestvisibalh[0][i+1] = lowestvisibalh[0][i]+mslope
    
    #判断对角线方向可视性并存入可视矩阵，更新最低可视海拔矩阵
    heightlist3 = [area[i+1][i+1] for i in range（radius）]
    sqrt2 = math.sqrt（2）
    mslope = -65535
    for i in range（radius）:
        tslope = float（heightlist3[i] - viewheight