# -*- coding: iso-8859-7 -*-
from math import hypot
from p_ggen import xfrange
from p_gmath import linint
def prnone(t): pass
class ThanYpyka:
"Computes the isocurves from a triangulation."
def __init__(self, ls, saveis, prt=prnone):
"Inintialize the object with the links of a triangulation."
self.ls = ls #Triangulation links
self.saveis = saveis #User supplied function to accept isocurves as they are computed
self.prt = prt #User supplied function to print messages
self.cis = [] #The coordinates of an isocurve
self.seen = set() #Visited edges
self.dhl = 1.0 #The iso-dimension
self.apOrioMax = 50.0 #Longer edges than this are considered nonexistant
def ypyka(self, dhl=None, apmax=None):
"Compute all isocurves."
if dhl != None: self.dhl = dhl
if apmax != None: self.apOrioMax = apmax
assert len(self.ls) > 2, " 3 "
#----- a a
for k in self.ls: break
hmin = hmax = k[2]
for k in self.ls:
h1 = k[2]
if h1 > hmax: hmax = h1
if h1 < hmin: hmin = h1
#-----
h1 = int(hmin / dhl) * dhl
for his in xfrange(h1, hmax, dhl):
pr = 1.0
if his > (hmin+hmax)*0.5: pr = -1.0
self.ypyka1(his, pr)
self.seen.clear() # Free memory
del self.cis[:] # Free memory
#=====================================================================
def ypyka1(self, his, pr):
"Compute all isocurves of elevation his."
self.prt(' %.2f' % his)
#----- a' a his
dhis = 0.01*self.dhl
seen = self.seen
seen.clear()
for karx, linksarx in self.ls.iteritems():
harx = karx[2]
if harx == his: harx -= dhis
if harx*pr > his*pr: continue
for iarx, larx in enumerate(linksarx):
if frozenset((karx, larx)) in seen: continue
hl = larx[2]
if hl == his: hl -= dhis
if hl*pr < his*pr: continue
if self.orioper(karx, larx): continue
icod = self.ypyka2(his, karx, iarx)
self.saveis(icod, self.cis) #
#==========================================================================
def orioper(self, k, l):
"Tests if distance is out of bounds."
d2 = hypot(l[0]-k[0], l[1]-k[1])
return d2 > self.apOrioMax
#==========================================================================
def ypyka2(self, his, karx, iarx):
"Compute 1 isocurve."
#-----
del self.cis[:]
icod = self.mishyp(his, karx, iarx, 1)
#-----
if icod < 0:
if self.cis[0] != self.cis[-1]:
assert 0, " : ."
icod = 0
#----
else:
self.cis.reverse()
del self.cis[-1] # ..
self.seen.remove(frozenset((karx, self.ls[karx][iarx]))) # ..
icod = self.mishyp(his, karx, iarx, -1)
if icod == -1:
assert 0, " : ."
icod = 0
return icod
#=====================================================================
def mishyp (self, his, karx, iarx, ibhm):
"Compute half of a isocurve."
seen = self.seen
dhis = 0.01*self.dhl
k = karx
hk = k[2]
if hk == his: hk -= dhis
linksk = self.ls[k]
i = iarx
l = linksk[i]
hl = l[2]
if hl == his: hl -= dhis
edge = frozenset((k, l))
#----- - a k
while True:
self.tomis(k, hk, l, hl, his)
seen.add(edge)
lp = l
#--------- ksp(k)
i = (i + ibhm) % len(linksk)
#--------- ,
l = linksk[i]
hl = l[2]
if hl == his: hl -= dhis
if lp not in self.ls[l]: return 0 #
edge = frozenset((k, l))
if edge in seen: # ();
assert (hk-his) * (his-hl) >= 0.0, " . !"
self.tomis(k, hk, l, hl, his)
return -1
if (hk-his) * (his-hl) >= 0.0:
if self.orioper(k, l): return 0 #
continue
#---------aa
k = l
hk = k[2]
if hk == his: hk -= dhis
linksk = self.ls[k]
for i, l in enumerate(linksk):
if l == lp: break
else:
assert 0, " : !"
hl = l[2]
if hl == his: hl -= dhis
edge = frozenset((k, l))
if edge in seen: # ();
self.tomis(k, hk, l, hl, his)
return -1
assert (hk-his) * (his-hl) >= 0.0, " !"
if self.orioper(k, l): return 0 #
#==========================================================================
def tomis(self, k, hk, l, hl, his):
"Find intersection and save it."
#-----
dh = hl - hk
assert dh != 0.0, " 0.01*dhl, ;"
ca = list(k)
ca[2] = hk
cb = list(l)
cb[2] = hl
if dh == 0.0:
ct = [(a+b)*0.5 for a,b in zip(ca, cb)]
else:
ct = [linint(hk, a, hl, b, his) for a,b in zip(ca, cb)]
#-----
self.cis.append(ct)
def test():
"Test the computation of contour lines."
from p_ggen import prg
from p_gchart import ThanChart,vis
from tt import ThanTri
ch = ThanChart()
def showis(icod, cis):
xx = [c[0] for c in cis]
yy = [c[1] for c in cis]
ch.curveAdd(xx, yy)
print "Reading tri.."
tri = ThanTri()
tri.readtri(open("013.tri"))
print "Producing contours.."
p = ThanYpyka(tri.ls, showis, prg)
p.ypyka(1, 200)
del p
vis(ch)
if __name__ == "__main__": test()
|