Marching squares

  描述

Marching squares是计算机图形学中一个用于矩阵网格点数据生成等值面的一个算法。

      
其主要思想是将每个2*2网格，其中4个顶点网格点根据等值线阈值处理为0，1标量。

   

• 1 代表数值大于阈值
  
• 0 代表数值小于阈值

这样处理之后有16中情况，2*2网格点根据左上角顺时针编号。

对于对角线这种情形，依据对角线大于两边处理。

算法主要流程：

      1. 遍历每个小网格，从16中固定情形中选择划线类别

      2.利用线性插值，结合网格点数值找寻交点

      3.  划线

代码

#coding:utf-8from pylab import *from matplotlib.ticker import MultipleLocator, FormatStrFormatterclass PlotDemo:def __init__(self,m,n):self._net = Noneself.m = mself.n = ndef _set_default_figure(self,m,n):ax = subplot(111) #注意:一般都在ax中设置,不再plot中设置#设置主刻度标签的位置,标签文本的格式xmajorLocator   = MultipleLocator(10) #将x主刻度标签设置为20的倍数xmajorFormatter = FormatStrFormatter('%1.1f') #设置x轴标签文本的格式xminorLocator   = MultipleLocator(1) #将x轴次刻度标签设置为5的倍数ymajorLocator   = MultipleLocator(10) #将y轴主刻度标签设置为0.5的倍数ymajorFormatter = FormatStrFormatter('%1.1f') #设置y轴标签文本的格式yminorLocator   = MultipleLocator(1) #将x轴次刻度标签设置为5的倍数ax.set_xlim(0,m)ax.set_ylim(0,n)ax.xaxis.set_major_locator(xmajorLocator)ax.xaxis.set_major_formatter(xmajorFormatter)ax.yaxis.set_major_locator(ymajorLocator)ax.yaxis.set_major_formatter(ymajorFormatter)#显示次刻度标签的位置,没有标签文本ax.xaxis.set_minor_locator(xminorLocator)ax.yaxis.set_minor_locator(yminorLocator)ax.xaxis.grid(True, which='minor') #x坐标轴的网格使用主刻度ax.yaxis.grid(True, which='minor') #y坐标轴的网格使用次刻度def show_source(self):self._set_default_figure(self.m, self.n)if self._net is None:raise  Exceptionnet = self._net.net_infoshape = net.shapexlen = shape[0]ylen = shape[1]for i in range(xlen):for j in range(ylen):if net[i][j] > 0:plot(i,j,'g+')show()def show_contour(self):self._set_default_figure(self.m, self.n)net = self._net.net_infoshape = net.shapexlen = shape[0]ylen = shape[1]for i in range(xlen):for j in range(ylen):if net[i][j] > 0:plot(i,j,'g+')net = self._netutils = MarchSquareUtlis(net)lines = utils.trancing_contours()width,height = net.net_info.shapearr = net.net_infoidx = 0for i in range(width-1):for j in range(height-1):x,y = i,jcount,v1,v2,v3,v4,v5,v6,v7,v8 = lines[idx]idx = idx + 1if count == 0:continueif count == 1:x1 = x + v1y1 = y + v2x2 = x + v3y2 = y + v4plot(x1,y1,x2,y2,'ro')if count == 2:x1 = x + v1y1 = y + v2x2 = x + v3y2 = y + v4plot(x1,y1,x2,y2,'ro')x1 = x + v5y1 = y + v6x2 = x + v7y2 = y + v8plot(x1,y1,x2,y2,'ro')show()def set_net_info(self, net_info):self._net = net_infoclass RandomGenNet(object):def __init__(self):self.arr = self._gen_random()def __init__(self,m,n):self.arr = self._gen_random(m,n)@propertydef net_info(self):    return self.arrdef _gen_random(self,m=100,n=100):return np.zeros((m,n),dtype='double')   def add_circle(self,center_x, center_y, radius,val):r = int(radius + 0.5)for x in range(center_x - r, center_x + r):rx = np.abs(center_x - x)ry = int( np.sqrt(radius * radius - rx * rx))for y in range(ry):self.arr[x][center_y - y] = valself.arr[x][center_y + y] = valreturndef add_retangle(self,lf_up_x, lf_up_y, width, height,val):for i in range(width):for j in range(height):self.arr[lf_up_x + i][lf_up_y + j] = val#clock-wise, top-left|top-right|bottom-right|bottom_leftdef get_retangle_bit(v1, v2, v3, v4):return v1 << 3 | v2 << 2 | v3 << 1 | v4#shift relative to top-leftdef get_retangle_shift(bitval):if bitval == 0 or bitval == 15:return (0,None,None,None,None,None,None,None,None)if bitval == 1 or 14:return (1,0,0.5,0.5,1,None,None,None,None)if bitval == 2 or bitval == 13:return (1,0.5,1,1,0.5,None,None,None,None)if bitval == 3 or bitval == 12:return (1,0,0.5,1,0.5,None,None,None,None)if bitval == 4 or bitval == 11:return (1,0,0.5,1,0.5,None,None,None,None)if bitval == 5:return (2,0,0.5,0.5,0,0.5,1,1,0.5)if bitval == 6 or bitval == 9:return (1,0.5,0,0.5,1,None,None,None,None)if bitval == 7 or bit == 8:return (1,0,0.5,0.5,0,None,None,None,None)if bitval == 10:return (2,0,0.5,0.5,1,0.5,0,1,0.5)class MarchSquareUtlis(object):def __init__(self, net):self.net = netdef trancing_contours(self):ret = []width,height = self.net.net_info.shapearr = self.net.net_infofor i in range(width-1):for j in range(height-1):v1 = int(arr[i][j])v2 = int(arr[i + 1][j])v3 = int(arr[i + 1][j + 1])v4 = int(arr[i][j + 1])bitv = get_retangle_bit(v1,v2,v3,v4)net_shift = get_retangle_shift(bitv)ret.append(net_shift)return retdemo = PlotDemo(100,100)netinfo = RandomGenNet(100,100)netinfo.add_circle(20,20,10,1)netinfo.add_circle(70,50,20,1)netinfo.add_retangle(20,45,40,20,1)netinfo.add_retangle(70,50,10,10,1)demo.set_net_info(netinfo)#demo.show_source()demo.show_contour()print get_retangle_bit(1,1,1,0)print range(3)

参考

Marching_squares_wiki