#-*-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?????????????x+=1
????????????if（dp?????????????????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?????????????x+=1
????????????if（dp?????????????????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