Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.arcetri.astro.it/SWXCS/Files/sweepline.py Дата изменения: Mon Dec 17 02:20:48 2012 Дата индексирования: Sun Apr 10 01:41:12 2016 Кодировка: Поисковые слова: http astrokuban.info astrokuban

#!/usr/bin/env python
################################################################################
#sweepline.py
#Voronoi Tessellation with Fortune's sweep line algorithm
#A discrete Voronoi diagram can also be constructed.
#
#Copyright (C) 2011 Teng Liu
#
#This program is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#(at your option) any later version.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program. If not, see .
#
#Teng Liu's email: lewtonstein@gmail.com
################################################################################
import numpy as np
import pylab as pl
import pyfits as pf
import heapq
import itertools,pickle,cPickle
import time,sys,warnings,os
warnings.simplefilter('ignore')
color = lambda s,ncolor,nfont: "\033["+str(nfont)+";"+str(ncolor)+"m"+s+"\033[1;m"
#define:
#p #vertical,upward: y++
#horizontal,right: x++
debug = False
SLVround = 6
class HalfEdge(object):
def __init__(self,p0,p1,direct,base,twin=None,summit=None):
#if p0[1] # #the 2 points on 2 sides, p0 < p1
self.p0 = p0
self.p1 = p1
self.direct = direct # 0:upward -1:left 1:right
self.base = base
#history test
if self.direct == 0:
self.hist = (self.p0[0]+self.p1[0])/2.
else:
self.hist = [None,None]
self.twin = twin
self.summit = summit
def __str__(self):
return 'p0:'+str(self.p0)+' p1:'+str(self.p1)+' dir:'+str(self.direct)+' b:'+str(self.base)+' s:'+str(self.summit)
def copy(self):
return HalfEdge(self.p0,self.p1,self.direct,self.base,twin=self.twin,summit=self.summit)
def isRightOf(self,x,y):
#Intersection X > px
if self.p0 == (-1,-1) : return False
if self.p1 == (-2,-2): return True
if self.direct == 0: return x else:
if self.hist==[None,None] or y>self.hist[0]:
if self.p1[1]>self.p0[1]: u,v = self.p1,self.p0
else: v,u = self.p1,self.p0
insqrt = (u[1]*v[1]-u[1]*y-v[1]*y+y*y)*(u[0]*u[0]-2*u[0]*v[0]+u[1]*u[1]-2*u[1]*v[1]+v[0]*v[0]+v[1]*v[1])
#if insqrt<0: fuhao=-1
#else: fuhao=1
insqrt = round(insqrt,SLVround)
#this round is becaue there may be very small negative value
zx = (self.direct*np.sqrt(insqrt)-u[0]*v[1]+u[0]*y+u[1]*v[0]-v[0]*y)/(u[1]-v[1])
#if fuhao==-1: print insqrt,self,x,y,zx, (self.direct*(1)*np.sqrt(-insqrt)-u[0]*v[1]+u[0]*y+u[1]*v[0]-v[0]*y)/(u[1]-v[1]), (self.direct*(-1)*np.sqrt(-insqrt)-u[0]*v[1]+u[0]*y+u[1]*v[0]-v[0]*y)/(u[1]-v[1])
self.hist = [y,zx]
else:
zx = self.hist[1]
return x #the method from the paper doesn't work, so I use mine
#f=lambda u0,u1,v0,v1,y0,d: (d*np.sqrt((u1*v1-u1*y0-v1*y0+y0*y0)*(u0*u0-2*u0*v0+u1*u1-2*u1*v1+v0*v0+v1*v1))-u0*v1+u0*y0+u1*v0-v0*y0)/(u1-v1)
#by http://www.wolframalpha.com/
#solve[(x-u0)^2+(y-u1)^2=r^2,(x-v0)^2+(y-v1)^2=r^2,y0-y=r,x,y,r]
def Intersect(self,Next):
if self.p1!=Next.p0: sys.exit('ERROR: not consecutive')
if self.p0==(-1,-1): return None
if Next.p1==(-2,-2): return None
if self.direct<=0 and Next.direct>=0: return None #not Intersect
#1: a<0 b>0 2: a<0 b=0 3: a=0 b>0 4: a=0 b=0
ux,uy = self.p0
vx,vy = self.p1
wx,wy = Next.p1
a = (ux-vx)*(vy-wy)-(uy-vy)*(vx-wx)
if a<=0: return None
b1 = (ux-vx)*(ux+vx)+(uy-vy)*(uy+vy)
b2 = (vx-wx)*(vx+wx)+(vy-wy)*(vy+wy)
cx = (b1*(vy-wy)-b2*(uy-vy))/2./a
cy = (b2*(ux-vx)-b1*(vx-wx))/2./a
d = np.sqrt((cx-ux)**2+(cy-uy)**2)
return cx,cy+d
def complete(self,x,y):
xb,yb = self.base
if xb<=-0.5 and x<=-0.5 or xb>=Voronoi.ImgSizeX-0.5 and x>=Voronoi.ImgSizeX-0.5: self.summit = None
elif yb<=-0.5 and y<=-0.5 or yb>=Voronoi.ImgSizeY-0.5 and y>=Voronoi.ImgSizeY-0.5: self.summit = None
else:
x1,y1 = self.p0
x2,y2 = self.p1
if y>-0.5 and yb<-0.5: #crossing y==0
#for the edge with (x1,y1) and (x2,y2) on its two sides
#find the point (x,0) which is equi-distance to the two points
self.base = (x2**2 + (y2+0.5)**2 - x1**2 - (y1+0.5)**2) /2. /(x2-x1), -0.5
xb,yb = self.base
if ybVoronoi.ImgSizeY-0.5: #crossing y==max
x,y = (x2**2 + (y2+0.5-Voronoi.ImgSizeY)**2 - x1**2 - (y1+0.5-Voronoi.ImgSizeY)**2) /2. /(x2-x1), Voronoi.ImgSizeY-0.5
if y<-0.5 and yb>-0.5: #back crossing y==0
x,y = (x2**2 + (y2+0.5)**2 - x1**2 - (y1+0.5)**2) /2. /(x2-x1), -0.5
if x<-0.5 or x>Voronoi.ImgSizeX-0.5: #then check x
if self.direct>0: xx=Voronoi.ImgSizeX-0.5
elif self.direct<0: xx=-0.5
x,y = xx, (y-yb)/(x-xb)*(xx-xb) + yb
self.summit = x,y
#return self.summit
def line(self,x,y):
xb,yb = self.base
if yb<=-0.5 and y<=-0.5 or yb>=Voronoi.ImgSizeY-0.5 and y>=Voronoi.ImgSizeY-0.5:
return ""
if y>-0.5 and yb<-0.5:
#for the edge with (x1,y1) and (x2,y2) on its two sides
#find the point (x,0) which is equi-distance to the two points
x1,y1 = self.p0
x2,y2 = self.p1
xb = (x2**2 + (y2+0.5)**2 - x1**2 - (y1+0.5)**2) /2. /(x2-x1)
yb = -0.5
if ybVoronoi.ImgSizeY-0.5:
x1,y1 = self.p0
x2,y2 = self.p1
x = (x2**2 + (y2+0.5-Voronoi.ImgSizeY)**2 - x1**2 - (y1+0.5-Voronoi.ImgSizeY)**2) /2. /(x2-x1)
y = Voronoi.ImgSizeY-0.5
if y<-0.5 and yb>-0.5:
x1,y1 = self.p0
x2,y2 = self.p1
x = (x2**2 + (y2+0.5)**2 - x1**2 - (y1+0.5)**2) /2. /(x2-x1)
y = -0.5
return "line(%.3f,%.3f,%.3f,%.3f) # line=0 0 tag={%d,%d-%d,%d-%d} font=\"helvetica 1 normal roman\"" % (yb+1,xb+1,y+1,x+1,self.p0[0],self.p0[1],self.p1[0],self.p1[1],self.direct)
class SweepTable(dict):
#a dict of half edges
def __init__(self,Px,Py):
self.KL = [] #key list
self.ysweep=0
self.p = None

ymin = Py>Py.min()
self.Q = EventQueue(Px[ymin],Py[ymin])

#initialize self dict
ymin = np.where(Py==Py.min())[0]
x0,y0 = Px[ymin[0]],Py[ymin[0]]
self[-1] = HalfEdge((-1,-1),(x0,y0),-1,(-1,-1)) #0: left boundary of image
self.KL.append(-1)
self.counter = itertools.count(1) #to make sure each key is different
for n in np.arange(1,ymin.size):
x = Px[ymin[n]]
y = Py[ymin[n]]
count = self.counter.next()
self[count] = HalfEdge((x0,y0),(x,y),0,((x0+x)/2.,-0.5))
self.KL.append(count)
x0,y0 = x,y
self[-2] = HalfEdge((x0,y0),(-2,-2),1,(-2,-2)) #1: right boundary of image
self.KL.append(-2)

def __str__(self):
s = ""
for k in self.KL:
s += "\t"+str(k)+' : '+self[k].__str__()+'\n'
return s

def locate(self,p):
#find the key of first edge behind which p locates
#X(key) <= px < X(next key)
Nlow = 0
Nhigh = len(self.KL)-1
#print Nlow,'--',Nhigh
while True:
if Nlow+1 == Nhigh:
break
Nmid = (Nlow+Nhigh)/2
#print 'mid:',Nmid
if self[self.KL[Nmid]].isRightOf(p[0],p[1]):
Nhigh = Nmid
else:
Nlow = Nmid
return Nlow

def insert(self,n,x,y):
#input: self[KL[n]].X <= p < self[KL[n+1]].X
#p: y,x,...
x0,y0 = self[self.KL[n]].p1
base = x,(y**2-y0**2-(x-x0)**2)/2./(y-y0)
count = self.counter.next()
self.KL.insert(n+1,count)
self[count] = HalfEdge((x0,y0),(x,y),-1,base,count+1) #new minus half edge
#if debug: print 'new: ',count,self[count]
count = self.counter.next()
self.KL.insert(n+2,count)
self[count] = HalfEdge((x,y),(x0,y0),1,base,count-1) #new plus half edge
#if debug: print 'new: ',count,self[count]

def shrink(self):
x,y,l,r = self.p
n = self.KL.index(l)
if r!=self.KL[n+1]: sys.exit('l-not-r')
else:
L = self.KL[n-1]
R = self.KL[n+2]
self.Q.deleteifhas(L)
self.Q.deleteifhas(r)
#(x,y) is the top point of the circle which goes through the three points (l.p0,l.p1,r.p1), l.p1==r.p0
#the center of the circle is (x,y-r)
#let x0,y0 be anyone of the three points
#(x-x0)^2 + (y-r-y0)^2 = r^2
if self[l].p1!=self[r].p0: sys.exit('no')
if self[l].p0[1]==self[r].p1[1]:
direct = 0
else:
if self[l].p0[1]>self[r].p1[1]:
x0 = self[l].p0[0]
if x==x0: direct = 1
else:
x0 = self[r].p1[0]
if x==x0: direct = -1
if x < x0: direct = -1
elif x > x0: direct = 1
x0,y0 = self[l].p0
xb,yb = x,(y**2-y0**2-(x-x0)**2)/2./(y-y0)
#if debug: print 'shrink:',self[l],self[r],x,y,(xb,yb)
count = self.counter.next()
self[count] = HalfEdge(self[l].p0,self[r].p1,direct,(xb,yb))
self.KL[n] = count
self.KL.pop(n+1)
self[l].complete(xb,yb)
self[r].complete(xb,yb)
i = self[L].Intersect(self[count])
if i is not None:
#if debug: print '->V',i,self[L],self[count]
self.Q.insert(i,L,count)
i = self[count].Intersect(self[R])
if i is not None:
#if debug: print '->V',i,self[count],self[R]
self.Q.insert(i,count,R)
return l,r
def renewEnds(self):
y = self.p[1]
todo = []
if y <= self.ysweep or len(self)<=2: return todo
if self[self.KL[-2]].isRightOf(Voronoi.ImgSizeX-0.5,y):
#if debug: print 'EdgeR:',self[-2],self[self.KL[-2]]
if self[-2].p0!=self[self.KL[-2]].p1: sys.exit('should not')
u0,u1 = self[self.KL[-2]].p0
v0,v1 = self[self.KL[-2]].p1
if u1==v1: sys.exit('cannot')
yb = (u0**2-2*u0*(Voronoi.ImgSizeX-0.5)+u1**2-v0**2+2*v0*(Voronoi.ImgSizeX-0.5)-v1**2)/2./(u1-v1)
self[self.KL[-2]].complete(Voronoi.ImgSizeX-0.5,yb)
todo.append([self.KL[-2],self[self.KL[-2]].copy()])
if self.KL[-2] in self.Q:
print 'Please check'
print self[self.KL[-2]]
print self[self.KL[-1]]
exit()
#maybe the VertexEvt to delete is just the p which has just been deleted from Q as the minimum.
if self.p is not None and len(self.p)==4 and self.KL[-3]==self.p[2]: self.p = None
else: self.Q.deleteifhas(self.KL[-3])
self[-2] = self.pop(self.KL[-2])
self.KL.pop(-2)
self[-2].p1 = (-2,-2)
self[-2].direct = 1
#if debug: print '-->',self[-2]
if not self[self.KL[1]].isRightOf(-0.5,y):
#if debug: print 'EdgeL:',self[-1],self[self.KL[1]]
if self[-1].p1!=self[self.KL[1]].p0: sys.exit('should not')
u0,u1 = self[self.KL[1]].p0
v0,v1 = self[self.KL[1]].p1
if u1==v1: sys.exit('cannot')
yb = (u0**2-2*u0*(-0.5)+u1**2-v0**2+2*v0*(-0.5)-v1**2)/2./(u1-v1)
self[self.KL[1]].complete(-0.5,yb)
todo.append([self.KL[1],self[self.KL[1]].copy()])
if self.p is not None and len(self.p)==4 and self.KL[1]==self.p[2]: self.p = None
else: self.Q.deleteifhas(self.KL[1])
self[-1] = self.pop(self.KL[1])
self.KL.pop(1)
self[-1].p0 = (-1,-1)
self[-1].direct = -1
#if debug: print '-->:',self[-1]
self.ysweep = y
return todo
class EventQueue(dict):
def __init__(self,Px,Py):
if Px.shape==Py.shape==(Px.size,):
self.S = np.zeros((Px.size,2))
self.S[:,0] = Py
self.S[:,1] = Px
self.SMin = 0
self.counter = itertools.count(1)
self.VerEvt = []
else:
sys.exit('ERROR: shape error')
def __str__(self):
return self.__repr__()+'\n'+str(self.VerEvt)
def delMin(self):
while len(self.VerEvt)>0 and self.VerEvt[0][2] == 0: heapq.heappop(self.VerEvt)
#The interesting property of a heap is that its smallest element is always the root, heap[0].
if self.SMin self.SMin += 1
return self.S[self.SMin-1][[1,0]]
elif len(self.VerEvt)>0 and (self.SMin>=self.S.shape[0] or self.S[self.SMin].tolist() > self.VerEvt[0][:2]):
y,x,count,l,r = heapq.heappop(self.VerEvt)
self.pop(l)
return [x,y,l,r]
else:
return None
def insert(self,p,l,r,count=None):
if l in self: sys.exit(str(l)+'already in')
if count is None: count=self.counter.next()
evt = [p[1],p[0],count,l,r]
self[l] = evt
heapq.heappush(self.VerEvt,evt)
def deleteifhas(self,l):
#if debug: print '><:',l
evt = self.pop(l,None)
if evt is not None: evt[2] = 0
class Voronoi(object):
ImgSizeX = 0
ImgSizeY = 0
def __init__(self,img,**kwargs):
self.FileName = kwargs.pop('FileName','FileName')
print color('Voronoi Construction '+self.FileName,34,1)
Voronoi.ImgSizeX,Voronoi.ImgSizeY = img.shape
self.pmap = np.zeros(img.shape,dtype='int32') #dvd
self.Py,self.Px = np.where(img.T>0) #coordinates (y,x)
self.pmap[self.Px,self.Py] = np.arange(1,self.Px.size+1)
self.Edges = {}
self.construct(**kwargs)
self.Amap = None
self.PPdis = None
self.Adj = None
self.EdgePoint = None
if kwargs.pop('calarea',None): self.calArea()
if kwargs.pop('caldvd',None): self.calDVD()
if kwargs.pop('caladj',None): self.calAdj()
#if __name__ != '__main__': pickle.dump(self, file(self.FileName+'_vor.dat','w'))
#END_OF_init

@staticmethod
def calculate(img,**kwargs):
try:
FileName = kwargs.get('FileName','not exist')
v = pickle.load(file(FileName+'_vor.dat','r'))
except:
v = Voronoi(img,**kwargs)
else:
print "<< "+FileName+'_vor.dat'
if kwargs.pop('caldvd',None): #in case loaded vor.dat has no DVD calculated
if v.pmap.max()==0: v.calDVD()
if kwargs.pop('caldst',None):
st=time.time()
print '>> '+v.FileName+'_dst.fits'
dmap = np.zeros(img.shape,dtype='float64')
Pvalue = img[v.Px,v.Py] / v.Amap[v.Px,v.Py]
dmap[v.Px,v.Py] = Pvalue
pf.writeto(v.FileName+'_dst.fits',dmap,clobber=True)
dmap[v.pmap>0] = Pvalue[v.pmap[v.pmap>0]-1]
pf.writeto(v.FileName+'_Dst.fits',dmap,clobber=True)
if debug: print time.time()-st

def construct(self,**kwargs):
StartTime=time.time()
T = SweepTable(self.Px,self.Py)
Q = T.Q
T.p = Q.delMin()
while T.p is not None:
for k,e in T.renewEnds():
self.Edges[k] = e
if T.p is None: #p is deleted during renewing ending half edges.
pass
elif len(T.p)==2: #Site Event
#if debug: print 'SiteEvent:',T.p
n = T.locate(T.p)
L = T.KL[n]
R = T.KL[n+1]
T.insert(n,T.p[0],T.p[1])
l = T.KL[n+1]
r = T.KL[n+2]
if R!= T.KL[n+3]: sys.exit('no')
Q.deleteifhas(L)
i1 = T[L].Intersect(T[l])
i2 = T[r].Intersect(T[R])
if i1 is not None:
#if debug: print '->V',i1,T[L],T[l]
Q.insert(i1,L,l)
if i2 is not None:
#if debug: print '->V',i2,T[r],T[R]
Q.insert(i2,r,R)

elif len(T.p)==4: #Vertex Event
if T.p[2] not in T or T.p[3] not in T:
print 'Useless VertexEvent:',T.p
else:
#if debug: print 'VertexEvent:',T.p[0],T.p[1]
l,r = T.shrink()
self.Edges[l] = T.pop(l)
self.Edges[r] = T.pop(r)
else: sys.exit('error2')
#if debug: print 'T',T
#if debug: print 'Q',Q
T.p = Q.delMin()
#END_OF_while p is not None:
T.pop(-1)
T.pop(-2)
for l in T.keys():
e = T[l]
if e.direct==0: e.complete(e.base[0],Voronoi.ImgSizeY-0.5)
else: e.complete(e.base[0],Voronoi.ImgSizeY) #out of region, summit will be recalculated
self.Edges[l] = T.pop(l)

for e in self.Edges.values():
if e.summit is not None:
if round(e.base[0],SLVround)==round(e.summit[0],SLVround) and round(e.base[1],SLVround)==round(e.summit[1],SLVround):
e.summit = None
continue
if e.twin is not None:
e2 = self.Edges[e.twin]
if e2.summit is not None:
assert e.base==e2.base
assert self.Edges[e2.twin]==e
e.base = e2.summit
e2.summit = None
assert e.p0==e2.p1 and e.p1==e2.p0
#Delete nouse edges
for k in self.Edges.keys():
if self.Edges[k].summit is None: self.Edges.pop(k)

#d = np.array([[e.base[1],e.base[0],e.summit[1],e.summit[0]] for e in self.Edges.values() if e.summit is not None])+1
d = np.array([[e.base[1],e.base[0],e.summit[1],e.summit[0],e.p0[1],e.p0[0],e.p1[1],e.p1[0]] for e in self.Edges.values() if e.summit is not None])+1
file(self.FileName+'.reg','w').write('image\n')
np.savetxt(self.FileName+'.reg',d,fmt="line(%.3f,%.3f,%.3f,%.3f) # tag={%d,%d, %d,%d}")
if debug: print 'time:',time.time()-StartTime
#END_OF_construct()

def drawLine(self,e):
if np.isinf(e.summit[0]): print e
if e.direct == 0:
if e.base[0]!=e.summit[0]: sys.exit('ERROR')
if e.base[1] y1 = e.base[1]
y2 = e.summit[1]
else:
y2 = e.base[1]
y1 = e.summit[1]
x = round(e.base[0],SLVround)
if x<0: x=0
if x == int(x):
y1 = round(y1,SLVround)
if y1<0: y1=0
elif y1==int(y1): y1 = int(y1)
else: y1 = int(y1)+1
y2 = int(y2)+1
self.pmap[e.base[0],y1:y2] = -1
return
elif e.direct < 0:
if e.p0[1]>e.p1[1]: sys.exit('error')
N1 = self.pmap[e.p0]
N2 = self.pmap[e.p1]
else: #e.direct > 0
if e.p1[1]>e.p0[1]: sys.exit('error')
N1 = self.pmap[e.p1]
N2 = self.pmap[e.p0]
if e.base[0] < e.summit[0]: #x1 x1,y1 = e.base
x2,y2 = e.summit
else:
x1,y1 = e.summit
x2,y2 = e.base
if x2==x1 and debug:
print e.base,e.summit
print 'Zero-Length Edges not eliminated clearly'
a = 1.*(y2-y1)/(x2-x1)
b = y1-a*x1
x1 = round(x1,SLVround)
if x1<0: x1 = 0
elif x1==int(x1): x1 = int(x1)
else: x1 = int(x1)+1
x2=int(x2)+1
if x2>Voronoi.ImgSizeX: x2 = Voronoi.ImgSizeX
X = np.arange(x1,x2)
Y = a*X+b
for n in np.arange(X.size):
Y[n] = round(Y[n],SLVround)
if Y[n]<0: Y[n]=0
nYn = int(Y[n])
if Y[n] == nYn:
self.pmap[X[n],Y[n]] = -1
y1 = nYn-1
y2 = nYn+1
else:
y1 = nYn
y2 = nYn+1
if y1>=0:
N = self.pmap[X[n],y1]
if N==0 or N>N1: self.pmap[X[n],y1] = N1
if y2 N = self.pmap[X[n],y2]
if N>=0 and N return

def calDVD(self):
for k in self.Edges.keys():
#print self.Edges[k]
self.drawLine(self.Edges[k])
if debug: pf.writeto(self.FileName+'_line.fits',self.pmap,clobber=True)

#Fill self.pmap using edge lines as marks
for x in np.arange(Voronoi.ImgSizeX):
if self.pmap[x,0]==0 or self.pmap[x,1]==0:
y0 = 0
N = self.pmap[x,0]
setNew = False
else:
setNew = True
for y in np.arange(1,Voronoi.ImgSizeY):
if self.pmap[x,y] == -1:
if not setNew:
if N==0: print "Can't do anything"
self.pmap[x,y0:y] = N
setNew = True
elif self.pmap[x,y] > 0:
if not setNew:
if y+1 if N==0: N = self.pmap[x,y]
self.pmap[x,y0:y] = N
setNew = True
else:
if y+1 y0 = y
N = self.pmap[x,y]
setNew = False
if not setNew:
self.pmap[x,y0:y+1] = N

if debug: pf.writeto(self.FileName+'_line1.fits',self.pmap,clobber=True)
#include multiple pixels which are equidistance and marked as -1
if (self.pmap==0).any(): sys.exit('ERROR: sth left')
VeryRareCase=[]
for x,y in np.argwhere(self.pmap==-1):
Ns = np.array([[fx,fy] for [fx,fy] in [[x+1,y],[x,y+1],[x-1,y],[x,y-1]] \
if fx>=0 and fy>=0 and fx0])
if Ns.size==0:
Ns = np.array([[fx,fy] for [fx,fy] in [[x+1,y+1],[x+1,y-1],[x-1,y+1],[x-1,y-1]] \
if fx>=0 and fy>=0 and fx0])
if Ns.size==0:
#very rare case, exist at a very low probability
VeryRareCase.append([x,y])
self.pmap[x,y] = 0
continue
Ns = self.pmap[Ns[:,0],Ns[:,1]]
self.pmap[x,y] = 0
for N in Ns:
if (Ns==N).sum()>1: self.pmap[x,y] = -N
if self.pmap[x,y]==0: self.pmap[x,y] = -Ns[np.random.randint(0,Ns.size)]
self.pmap[self.pmap<0] *= -1
#Finished
for x,y in VeryRareCase:
Ns = np.array([[fx,fy] for [fx,fy] in [[x+1,y],[x,y+1],[x-1,y],[x,y-1],[x+1,y+1],[x+1,y-1],[x-1,y+1],[x-1,y-1]] \
if fx>=0 and fy>=0 and fx0])
if Ns.size==0: sys.exit('ERROR: pixel %d,%d has no neighbor' %(y+1,x+1))
Ns = self.pmap[Ns[:,0],Ns[:,1]]
self.pmap[x,y] = 0
for N in Ns:
if (Ns==N).sum()>1: self.pmap[x,y] = N
if self.pmap[x,y]==0: self.pmap[x,y] = Ns[np.random.randint(0,Ns.size)]
if (self.pmap==0).any(): sys.exit('ERROR: calDVD not completed')
if __name__ == '__main__': pf.writeto(self.FileName+'_dvd.fits',self.pmap,clobber=True)
#END_OF_def calDVD(self)

def calAdj(self):
if debug: print color("\nCalculating Adjacency List",32,0)
assert (self.pmap[self.Px,self.Py] == np.arange(1,self.Px.size+1)).all()
self.Adj = [[] for _ in itertools.repeat(None,self.Px.size+1)]
for e in self.Edges.values():
self.Adj[self.pmap[tuple(e.p0)]].append(self.pmap[tuple(e.p1)])
self.Adj[self.pmap[tuple(e.p1)]].append(self.pmap[tuple(e.p0)])
#for _ in self.Adj[1:]: assert _
#END_OF_def calAdj(self):

def calArea(self):
#for e in self.Edges.values():
# assert e.summit is not None
if debug: print color("\nCalculating Voronoi Cell Areas",32,0)
self.Amap = np.zeros(self.pmap.shape,dtype='float64')
P0 = np.int32([e.p0 for e in self.Edges.values()])
P1 = np.int32([e.p1 for e in self.Edges.values()])
E0 = np.array([e.base for e in self.Edges.values()]).round(SLVround)
E1 = np.array([e.summit for e in self.Edges.values()]).round(SLVround)

DisPPList = lambda P0,P1: np.sqrt((P0[:,0]-P1[:,0])**2+(P0[:,1]-P1[:,1])**2)
self.PPdis = DisPPList(P0,P1)

#First calculate "Edge area list"
#"Edge area": area of the triangle made of each edge and one of its according points
Earea = DisPPList(E0,E1) * self.PPdis / 4.

#Then calculate Parea
#"Point area": area of each point (cell)
for n in np.arange(P0.shape[0]):
self.Amap[tuple(P0[n])] += Earea[n]
self.Amap[tuple(P1[n])] += Earea[n]

#Modify open cells at the image edge
#Ignore the special cases of the cells at the four corners
TriArea = lambda p1,p2,p3: abs(p2[0]*p3[1]+p3[0]*p1[1]+p1[0]*p2[1]-p1[0]*p3[1]-p2[0]*p1[1]-p3[0]*p2[1])/2.
EE = np.vstack((E0,E1))

########################################################################
#L R B T [min,max]
VyL = [Voronoi.ImgSizeY, -1]
VxB = [Voronoi.ImgSizeX, -1]
VyR = [Voronoi.ImgSizeY, -1]
VxT = [Voronoi.ImgSizeX, -1]
EpL = {}
EpB = {}
EpR = {}
EpT = {}
EdgePoint=[]
################################################################################
n = np.where((EE[:,0]==-0.5)&(EE[:,1]==-0.5))[0]
if n:
n = n[0]
x1,y1 = P0[n-P0.shape[0]]
x2,y2 = P1[n-P0.shape[0]]
if x1 < x2:
assert y1 > y2
EpL[self.pmap[x1,y1]] = [-0.5]
EpB[self.pmap[x2,y2]] = [-0.5]
else:
assert y1 < y2
EpL[self.pmap[x2,y2]] = [-0.5]
EpB[self.pmap[x1,y1]] = [-0.5]
EE[n] = [-1,-1]
n = np.where((EE[:,0]==-0.5)&(EE[:,1]==Voronoi.ImgSizeY-0.5))[0]
if n:
n = n[0]
x1,y1 = P0[n-P0.shape[0]]
x2,y2 = P1[n-P0.shape[0]]
if x1 < x2:
assert y1 < y2
EpL[self.pmap[x1,y1]] = [Voronoi.ImgSizeY-0.5]
EpT[self.pmap[x2,y2]] = [-0.5]
else:
assert y1 > y2
EpL[self.pmap[x2,y2]] = [Voronoi.ImgSizeY-0.5]
EpT[self.pmap[x1,y1]] = [-0.5]
EE[n] = [-1,Voronoi.ImgSizeY]
n = np.where((EE[:,0]==Voronoi.ImgSizeX-0.5)&(EE[:,1]==-0.5))[0]
if n:
n = n[0]
x1,y1 = P0[n-P0.shape[0]]
x2,y2 = P1[n-P0.shape[0]]
if x1 < x2:
assert y1 < y2
EpB[self.pmap[x1,y1]] = [Voronoi.ImgSizeX-0.5]
EpR[self.pmap[x2,y2]] = [-0.5]
else:
assert y1 > y2
EpB[self.pmap[x2,y2]] = [Voronoi.ImgSizeX-0.5]
EpR[self.pmap[x1,y1]] = [-0.5]
EE[n] = [Voronoi.ImgSizeX,-1]
n = np.where((EE[:,0]==Voronoi.ImgSizeX-0.5)&(EE[:,1]==Voronoi.ImgSizeY-0.5))[0]
if n:
n = n[0]
x1,y1 = P0[n-P0.shape[0]]
x2,y2 = P1[n-P0.shape[0]]
if x1 < x2:
assert y1 > y2
EpT[self.pmap[x1,y1]] = [Voronoi.ImgSizeX-0.5]
EpR[self.pmap[x2,y2]] = [Voronoi.ImgSizeY-0.5]
else:
assert y1 < y2
EpT[self.pmap[x2,y2]] = [Voronoi.ImgSizeX-0.5]
EpR[self.pmap[x1,y1]] = [Voronoi.ImgSizeY-0.5]
EE[n] = [Voronoi.ImgSizeX,Voronoi.ImgSizeY]
################################################################################
for n in np.where(EE[:,0]==-0.5)[0]:
N = self.pmap[tuple(P0[n-P0.shape[0]])]
EdgePoint.append(N)
EpL[N] = EpL.get(N,[])+[EE[n,1]]
N = self.pmap[tuple(P1[n-P0.shape[0]])]
EdgePoint.append(N)
EpL[N] = EpL.get(N,[])+[EE[n,1]]
VyL[0] = min(VyL[0],EE[n,1])
VyL[1] = max(VyL[1],EE[n,1])
for n in np.where(EE[:,1]==-0.5)[0]:
N = self.pmap[tuple(P0[n-P0.shape[0]])]
EdgePoint.append(N)
EpB[N] = EpB.get(N,[])+[EE[n,0]]
N = self.pmap[tuple(P1[n-P0.shape[0]])]
EdgePoint.append(N)
EpB[N] = EpB.get(N,[])+[EE[n,0]]
VxB[0] = min(VxB[0],EE[n,0])
VxB[1] = max(VxB[1],EE[n,0])
for n in np.where(EE[:,0]==Voronoi.ImgSizeX-0.5)[0]:
N = self.pmap[tuple(P0[n-P0.shape[0]])]
EdgePoint.append(N)
EpR[N] = EpR.get(N,[])+[EE[n,1]]
N = self.pmap[tuple(P1[n-P0.shape[0]])]
EdgePoint.append(N)
EpR[N] = EpR.get(N,[])+[EE[n,1]]
VyR[0] = min(VyR[0],EE[n,1])
VyR[1] = max(VyR[1],EE[n,1])
for n in np.where(EE[:,1]==Voronoi.ImgSizeY-0.5)[0]:
N = self.pmap[tuple(P0[n-P0.shape[0]])]
EdgePoint.append(N)
EpT[N] = EpT.get(N,[])+[EE[n,0]]
N = self.pmap[tuple(P1[n-P0.shape[0]])]
EdgePoint.append(N)
EpT[N] = EpT.get(N,[])+[EE[n,0]]
VxT[0] = min(VxT[0],EE[n,0])
VxT[1] = max(VxT[1],EE[n,0])
self.EdgePoint = np.unique(EdgePoint)
if debug:
EPfile = file('EdgePoint.reg','w')
for N in self.EdgePoint:
print >>EPfile,"circle(%d,%d,2) # color=red" % (self.Py[N-1]+1,self.Px[N-1]+1)
EPfile.close()
########################################################################
# For such case where ther is no crossing point in one edge of the image
if len(EpL)==0:
assert VyL == [Voronoi.ImgSizeY, -1]
n = np.argmin(self.Px)
N = self.pmap[self.Px[n],self.Py[n]]
assert len(EpB[N])==1 and len(EpT[N])==1 and EpB[N][0]==VxB[0] and EpT[N][0]==VxT[0]
EpL[N] = [-0.5, Voronoi.ImgSizeX-0.5]
EpB[N].append(-0.5)
EpT[N].append(-0.5)
if len(EpB)==0:
assert VxB == [Voronoi.ImgSizeX, -1]
n = np.argmin(self.Py)
N = self.pmap[self.Px[n],self.Py[n]]
assert len(EpL[N])==1 and len(EpR[N])==1 and EpL[N][0]==VyL[0] and EpR[N][0]==VyR[0]
EpB[N] = [-0.5, Voronoi.ImgSizeY-0.5]
EpL[N].append(-0.5)
EpR[N].append(-0.5)
if len(EpR)==0:
assert VyR == [Voronoi.ImgSizeY, -1]
n = np.argmax(self.Px)
N = self.pmap[self.Px[n],self.Py[n]]
assert len(EpB[N])==1 and len(EpT[N])==1 and EpB[N][0]==VxB[1] and EpT[N][0]==VxT[1]
EpR[N] = [-0.5, Voronoi.ImgSizeX-0.5]
EpB[N].append(Voronoi.ImgSizeX-0.5)
EpT[N].append(Voronoi.ImgSizeX-0.5)
if len(EpT)==0:
assert VxT == [Voronoi.ImgSizeX, -1]
n = np.argmax(self.Py)
N = self.pmap[self.Px[n],self.Py[n]]
assert len(EpL[N])==1 and len(EpR[N])==1 and EpL[N][0]==VyL[1] and EpR[N][0]==VyR[1]
EpT[N] = [-0.5, Voronoi.ImgSizeY-0.5]
EpL[N].append(Voronoi.ImgSizeY-0.5)
EpR[N].append(Voronoi.ImgSizeY-0.5)
########################################################################
for N in EpL.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpL[N])==1:
if EpL[N][0]==VyL[0] and N in EpB.keys():
#N refers to the left bottom pixel
assert len(EpB[N])==1 and EpB[N][0]==VxB[0]
EpL[N].append(-0.5)
EpB[N].append(-0.5)
elif EpL[N][0]==VyL[1] and N in EpT.keys():
#Only one cross point in this edge. N refers to the left top pixel
assert len(EpT[N])==1 and EpT[N][0]==VxT[0]
EpL[N].append(Voronoi.ImgSizeX-0.5)
EpT[N].append(-0.5)
else: raise RuntimeError('ERROR:noleftnoright')
for N in EpR.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpR[N])==1:
if EpR[N][0]==VyR[0] and N in EpB.keys():
assert len(EpB[N])==1 and EpB[N][0]==VxB[1]
EpR[N].append(-0.5)
EpB[N].append(Voronoi.ImgSizeX-0.5)
elif EpR[N][0]==VyR[1] and N in EpT.keys():
assert len(EpT[N])==1 and EpT[N][0]==VxT[1]
EpR[N].append(Voronoi.ImgSizeX-0.5)
EpT[N].append(Voronoi.ImgSizeX-0.5)
else: raise RuntimeError('ERROR:noleftnoright')
########################################################################
for N in EpL.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpL[N])!=2: print y+1,x+1,EpL[N]
assert len(EpL[N])==2
self.Amap[x,y] += (x+0.5)*abs(EpL[N][0]-EpL[N][1])/2.
#print "circle(%.1f,%.1f,1.3) # text={%.3f}" %(y+1,x+1,(x+0.5)*abs(EpL[N][0]-EpL[N][1])/2.)
for N in EpB.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpB[N])!=2: print y+1,x+1,EpB[N]
assert len(EpB[N])==2
self.Amap[x,y] += (y+0.5)*abs(EpB[N][0]-EpB[N][1])/2.
#print "circle(%.1f,%.1f,1.3) # text={%.3f}" %(y+1,x+1,(y+0.5)*abs(EpB[N][0]-EpB[N][1])/2.)
for N in EpR.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpR[N])!=2: print y+1,x+1,EpR[N]
assert len(EpR[N])==2
self.Amap[x,y] += (Voronoi.ImgSizeX-0.5-x)*abs(EpR[N][0]-EpR[N][1])/2.
#print "circle(%.1f,%.1f,1.3) # text={%.3f}" %(y+1,x+1,(Voronoi.ImgSizeX-0.5-x)*abs(EpR[N][0]-EpR[N][1])/2.)
for N in EpT.keys():
x,y = self.Px[N-1],self.Py[N-1]
if len(EpT[N])!=2: print y+1,x+1,EpT[N],N
assert len(EpT[N])==2
self.Amap[x,y] += (Voronoi.ImgSizeY-0.5-y)*abs(EpT[N][0]-EpT[N][1])/2.
#print "circle(%.1f,%.1f,1.3) # text={%.3f}" %(y+1,x+1,(Voronoi.ImgSizeY-0.5-y)*abs(EpT[N][0]-EpT[N][1])/2.)
########################################################################

if __name__ == '__main__': pf.writeto(self.FileName+'_area.fits',self.Amap,clobber=True)
#END_OF_def calArea(self):

#1.Useless check is not necessary anylonger, it's useless
#2.locate maybe improved to be faster
##By profile.run(), I found that __str__() for large data structure is very slow, sometimes printing it can be the most time-consuming function.
if __name__ == '__main__':
def usage():
print """sweepline.py imagefile.fits [options]
OPTIONS
--caldvd Calculate Discrete Voronoi Diagram
--caldst Calculate Density image
--calarea Calculate Area Map
"""
exit()
ImgFile=""
caldvd=False
caldst=False
calarea=False
for arg in sys.argv[1:]:
if arg == '--caldvd': caldvd=True
elif arg == '--calarea': calarea=True
elif arg == '--caldst':
caldvd=True
calarea=True
caldst=True
elif os.path.isfile(arg): ImgFile=arg
elif arg == '-d': debug = True
elif arg == '-h' or arg == '--help': usage()
if ImgFile!="": Voronoi.calculate(pf.getdata(ImgFile),FileName=ImgFile.rsplit('.',1)[0],caldvd=caldvd,caldst=caldst,calarea=calarea)
else: usage()