diff --git a/pNbody/Movie.py b/pNbody/Movie.py index 8e322cd..09f29af 100644 --- a/pNbody/Movie.py +++ b/pNbody/Movie.py @@ -1,453 +1,453 @@ import string from numpy import * from PIL import Image from PIL import ImageDraw from PIL import ImageFont from PIL import ImagePalette from PIL import ImageFont class Movie: ''' a Movie class ''' ########################## def __init__(self,name,mode=None): ########################## self.name = name self.mode = mode self.starttime = 0.0 self.stoptime = 0.0 self.dt = 0.0 self.npic = 0 self.numByte = 0 self.numLine = 0 self.headerlength = 0 self.blocksize = 0 self.filmsize = 0 self.shape = (0,0) self.current_npic = 0 # self.film ########################## def open(self,mode='r',readall=0): ########################## """ open a file """ numByte = 256 # open the file self.film = open(self.name, mode) ############################################### # read the header and find numByte and numLine ############################################### self.header = self.film.read(256) if len(self.header) != 256: raise "film length < 256, ,stop" self.film.seek(240) str1 = self.film.read(8) str2 = self.film.read(8) if (str1 !=" " and str2 !=" "): self.numByte = string.atoi(str1) self.numLine = string.atoi(str2) else: self.numByte = 256 self.numLine = 192 if self.mode == 'RGB': self.numLine = self.numLine*3 if self.numByte < 256: self.headerlength = 256 else: self.headerlength = self.numByte ############################################### # read the header and find numByte and numLine ############################################### self.blocksize = self.numLine*self.numByte + self.headerlength self.shape = (self.numByte,self.numLine) ############################# # check film length ############################# self.film.seek(0,2) self.filmsize = self.film.tell() self.film.seek(0) fmod(self.filmsize,self.blocksize) self.npic = self.filmsize/self.blocksize ############################# # read times ############################# self.film.seek(0) self.starttime = float(self.film.read(self.headerlength)[:8]) self.film.seek(0) self.moveto(self.npic) self.stoptime = float(self.film.read(self.headerlength)[:8]) self.film.seek(0) if self.npic > 1: self.dt = (self.stoptime-self.starttime)/self.npic else: self.dt = 0 self.current_npic = 0 ''' ############################# # go to the top of the file ############################# self.film.seek(0) if readall: # reading next label time = self.film.read(self.headerlength) if len(time) == self.headerlength: # record starttime self.starttime = string.atof(time[:8]) # reading next data data = self.film.read(self.numByte*self.numLine) if len(data) == self.numByte*self.numLine: # add one pict self.npic = self.npic + 1 ####################### # loop over all images ####################### ok = 1 while ok: ok = 0 # reading next label time = self.film.read(self.headerlength) if len(time) == self.headerlength: # reading next data data = self.film.read(self.numByte*self.numLine) if len(data) == self.numByte*self.numLine: ok = 1 self.npic = self.npic + 1 self.stoptime = string.atof(time[:8]) # go to the top of the file self.film.seek(0) if self.npic > 1: self.dt = (self.stoptime-self.starttime)/(self.npic-1) else: self.dt = 0. ''' ########################## def info(self): ########################## """ give info """ print ''' -- %s -- initial time : %8.3f final time : %8.3f dt : %8.3f number of images : %d number of collumns : %d number of lines : %d length of header : %d length of block : %d length of film : %d current_npic : %d ''' % (self.name,self.starttime,self.stoptime,self.dt,self.npic,self.numByte,self.numLine,self.headerlength,self.blocksize,self.filmsize,self.current_npic) ########################## def new(self,numByte,numLine): ########################## self.numByte = numByte self.numLine = numLine if self.numByte < 256: self.headerlength = 256 else: self.headerlength = self.numByte self.film = open(self.name, "w") ########################## def close(self): ########################## self.film.close() ########################## def read_one(self,mode=None): ########################## # reading next label time = self.film.read(self.headerlength) if len(time) == self.headerlength: # record starttime self.current_time = string.atof(time[:8]) # reading next data data = self.film.read(self.numByte*self.numLine) if len(data) == self.numByte*self.numLine: self.current_image = data if mode == "array": return reshape(fromstring(data,uint8),(self.numLine,self.numByte)) elif mode == "image": - return Image.fromstring("P",shape,(self.numByte,self.numLine)) + return Image.frombytes("P",shape,(self.numByte,self.numLine)) else: return data ########################## def read_one_with_time(self,mode=None): ########################## # reading next label time = self.film.read(self.headerlength) if len(time) == self.headerlength: # record starttime self.current_time = string.atof(time[:8]) shape = (self.numByte,self.numLine) if mode == "array": data = self.film.read(self.numByte*self.numLine) return time,reshape(fromstring(data,'b'),(self.numLine,self.numByte)) elif mode == "image": data = self.film.read(self.numByte*self.numLine) - return time,Image.fromstring("P",shape,data) + return time,Image.frombytes("P",shape,data) elif mode == "image_rgb": dataR = self.film.read(self.numByte*self.numLine) self.film.read(self.headerlength) dataG = self.film.read(self.numByte*self.numLine) self.film.read(self.headerlength) dataB = self.film.read(self.numByte*self.numLine) - imageR = Image.fromstring("L",shape,dataR) - imageG = Image.fromstring("L",shape,dataG) - imageB = Image.fromstring("L",shape,dataB) + imageR = Image.frombytes("L",shape,dataR) + imageG = Image.frombytes("L",shape,dataG) + imageB = Image.frombytes("L",shape,dataB) return time,Image.merge('RGB',(imageR,imageG,imageB)) else: data = self.film.read(self.numByte*self.numLine) return time,data else: return None,None ########################## def read(self,skip=0,mode='array'): ########################## ''' skip = 0 : read image at the current position skip = 1 : skip an image skip = -1 : read the image before (go back) skip = -2 : skip an image before (go back) ''' # move relative to the current position try: self.film.seek(skip*self.blocksize,1) except IOError: self.moveto(0) return -1,0,-1 self.current_npic = self.film.tell()/self.blocksize if self.current_npic>self.npic-1: self.moveto(self.npic) return -2,self.current_npic,-2 # reading next label time = self.film.read(self.headerlength) # record starttime self.current_time = string.atof(time[:8]) if mode == "array": data = self.film.read(self.numByte*self.numLine) return self.current_time,self.current_npic,transpose(reshape(fromstring(data,'b'),(self.numLine,self.numByte))) elif mode == "image": data = self.film.read(self.numByte*self.numLine) - return self.current_time,self.current_npic,Image.fromstring("P",self.shape,data) + return self.current_time,self.current_npic,Image.frombytes("P",self.shape,data) elif mode == "image_rgb": dataR = self.film.read(self.numByte*self.numLine) self.film.read(self.headerlength) dataG = self.film.read(self.numByte*self.numLine) self.film.read(self.headerlength) dataB = self.film.read(self.numByte*self.numLine) - imageR = Image.fromstring("L",self.shape,dataR) - imageG = Image.fromstring("L",self.shape,dataG) - imageB = Image.fromstring("L",self.shape,dataB) + imageR = Image.frombytes("L",self.shape,dataR) + imageG = Image.frombytes("L",self.shape,dataG) + imageB = Image.frombytes("L",self.shape,dataB) return self.current_time,self.current_npic,Image.merge('RGB',(imageR,imageG,imageB)) else: data = self.film.read(self.numByte*self.numLine) return self.current_time,self.current_npic,data ########################## def moveto(self,npic): ########################## npic = min(npic,self.npic-1) npic = max(npic,0) current_npic = self.film.tell()/self.blocksize dnpic = npic-current_npic db = dnpic*self.blocksize self.film.seek(db,1) self.current_npic = self.film.tell()/self.blocksize ########################### def write_pic(self,time,data): ########################### import string recsize = self.numByte record = '%8.3f' % time s1=string.ljust(record,240) s2=string.ljust(`self.numByte`,8) if self.mode == 'RGB': s3=string.ljust(`self.numLine/3`,8) else: s3=string.ljust(`self.numLine`,8) record = s1+s2+s3 record= string.ljust(record,recsize) self.film.write(record) self.film.write(data) ########################### def write_pic_rgb(self,time,dataR,dataG,dataB): ########################### import string recsize = self.numByte record = '%8.3f' % time s1=string.ljust(record,240) s2=string.ljust(`self.numByte`,8) s3=string.ljust(`self.numLine`,8) record = s1+s2+s3 record= string.ljust(record,recsize) self.film.write(record) self.film.write(dataR) self.film.write(dataG) self.film.write(dataB) ########################### def get_img(self,data): ########################### ''' can be replaced by read_one with option "image" ''' shape = (self.numByte,self.numLine) - image = Image.fromstring("P",shape,data) + image = Image.frombytes("P",shape,data) return image ################################### def append_h(numByte,numLine,datas): ################################### newdata = '' # loop over the lines for j in range(numLine): j1 = j*numByte # set cursors j2 = j1 + numByte # loop over the images for data in datas: newdata = newdata + data[j1:j2] # simply sum return newdata diff --git a/pNbody/libutil.py b/pNbody/libutil.py index 342f2b9..7c35b14 100644 --- a/pNbody/libutil.py +++ b/pNbody/libutil.py @@ -1,1897 +1,1897 @@ # -*- coding: iso-8859-1 -*- import types from copy import deepcopy from numpy import * from nbodymodule import * from mapping import * from palette import * from myNumeric import * try: import Tkinter as tk from PIL import ImageTk is_tk = True except ImportError: is_tk = False from PIL import Image from PIL import ImageDraw from PIL import ImageFont from PIL import ImagePalette import mpi import param import thread from mapping import * try: from scipy import ndimage from scipy import convolve as conv is_scipy = True except ImportError: is_scipy = False try: import libqt is_qt = True except: is_qt = False def get_fake_fct(nx=256,ny=256): # create a function mx = 2*pi my = 2*pi x,y = indices((nx,ny)) x = mx*(x-nx/2)/(nx/2) y = my*(y-ny/2)/(ny/2) R = sqrt(x**2+y**2) R1 = sqrt((x-pi)**2+(y-pi)**2) mat = cos(R)/(1+R) + 0.5*cos(R1)/(1+R1) return mat def get_n_per_task(ntot,NTask): n_per_task = zeros(NTask,int) for Task in range(NTask-1,-1,-1): n_per_task[Task] = ntot/NTask + ntot%NTask*(Task==0) return n_per_task def get_npart_all(npart,NTask): npart_all = zeros((NTask,len(npart)),int) i=0 for n in npart: npart_all[:,i] = get_n_per_task(n,mpi.mpi_NTask()) i = i + 1 return npart_all def old_get_n_per_task(ntot,NTask): nleft = ntot n_per_task = zeros(NTask,int) for Task in range(NTask-1,-1,-1): if nleft < 2*ntot/NTask: n_per_task[Task] = nleft else: n_per_task[Task] = ntot/NTask nleft = nleft - n_per_task[Task] if nleft!=0: raise "get_n_per_task : nleft != 0 (%d)"%(nleft) return n_per_task def cross_product(x,y): if x.shape!=y.shape: raise "cross_product error","x and y do not have the same shape" n = len(x) p = zeros((n,3),x.dtype) p[:,0] = x[:,1]*y[:,2] - x[:,2]*y[:,1] p[:,1] = x[:,2]*y[:,0] - x[:,0]*y[:,2] p[:,2] = x[:,0]*y[:,1] - x[:,1]*y[:,0] return p def myhistogram(a,bins): ''' Return the histogram (n x 1 float array) of the n x 1 array "a". "bins" (m x 1 array) specify the bins of the histogram. ''' n = searchsorted(sort(a),bins) n = concatenate([n,[len(a)]]) return n[1:]-n[:-1] def compress_from_lst(x,num,lst,reject=False): ''' Return the compression of x ''' # 1) sort the list ys = sort(lst) # 2) sort index in current file xs = sort(num) zs = take(arange(len(x)),argsort(num)) # sort 0,1,2,n following xs (or self.num) # 3) apply mask on sorted lists (here, getmask need xs and ys to be sorted) m = getmask(xs.astype(int),ys.astype(int)) if reject: m = logical_not(m) # 4) revert mask, following zs inverse transformation c = take(m,argsort(zs)) return compress(c,x,axis=0) def tranfert_functions(rmin,rmax,g=None,gm=None): ''' This function computes the normalized tranfer fonction from g and gm It is very usefull to tranform a linear vetor in a non linear one example of g: g = lambda r:log(r/rc+1) gm = lambda r:rc*(exp(r)-1) ''' rmin = float(rmin) rmax = float(rmax) f = lambda r:r # by default, identity fm = lambda r:r # by default, identity if g != None and gm != None: f = lambda r: ( g(r)-g(rmin) )/( g(rmax)-g(rmin) ) *(rmax-rmin) + rmin fm = lambda f: gm( (f-rmin)*( g(rmax)-g(rmin) )/( rmax-rmin ) + g(rmin) ) return f,fm def geter(n,rmin,rmax,g,gm): ''' Generate a one dimentional non linear array of r ''' n = int(n) dr = (rmax-rmin)/float(n-1) f,fm = tranfert_functions(rmin,rmax,g,gm) Rs = arange(0,rmax+dr,dr) Rs= fm(Rs) return Rs def geter2(n,rmin,rmax,g,gm): ''' Generate a one dimentional non linear array of r ''' n = int(n) dr = (rmax-rmin)/float(n-1) f,fm = tranfert_functions(rmin,rmax,g,gm) Rs = arange(0,rmax+dr,dr) Rs= f(Rs) return Rs def getr(nr=31,nt=32,rm=100.): ''' Return a sequence of number (n x 1 array), where n=nr+1 defined by: Pfenniger & Friedli (1994) ''' j = arange(nr+1) ct1 = 0.5+(nt/(pi+pi)) ct2 = (exp((nr)/ct1)-1.) r = rm*(exp(j/ct1)-1)/ct2 return r def invgetr(r,nr=31,nt=32,rm=100.): ''' From r, return the corresponding indexes. Inverse of getr function. ''' ct1 = 0.5+(nt/(pi+pi)) ct2 = (exp((nr)/ct1)-1.) i = ct1*log( ct2 * r/rm + 1) return i def get_eyes(x0,xp,alpha,dr): ''' Return the position of two eyes. x0 : position of the head xp : looking point theta : rotation of the head dr : distance of the eyes ''' x0 = array(x0,float) xp = array(xp,float) # distance between the observer and the looking point Robs = sqrt((x0[0]-xp[0])**2 + (x0[1]-xp[1])**2 + (x0[2]-xp[2])**2) # init x = array([[dr,-Robs,0.],[-dr,-Robs,0.]],float32) m = Nbody(pos=x,status='new') # first : put xp in the origin x0 = x0 - xp # angle in azimuth phi = arctan2(x0[1],x0[0]) + pi/2 # angle in nutation R = sqrt(x0[0]**2 + x0[1]**2) if R==0: if x0[2] >= 0: theta = pi else: theta = -pi else: theta = arctan(x0[2]/R) # rotate # rotations alpha if alpha != 0.: m.rotate(alpha,axis='y',mode='p') # rotation in nutation if theta!=0.: m.rotate(-theta,axis='x',mode='p') # rotation in azimuth if phi != 0.: m.rotate(phi,axis='z',mode='p') # translate m.translate(xp) x1 = m.pos[0,:] x2 = m.pos[1,:] return x1,x2 def apply_filter(mat,name=None,opt=None): ''' Apply a filter to an image ''' if name == "convol": nx = opt[0] ny = opt[1] sx = float(opt[2]) sy = float(opt[3]) nxd = (nx-1)/2 nyd = (ny-1)/2 # the kernel b = fromfunction(lambda j,i: exp(-(i-nxd)**2/sx**2 + -(j-nyd)**2/sy**2),(nx,ny)) s = sum(sum(b)) b = b/s # normalisation # conversion: b = b.astype(float) mat=mat.astype(float) return convol(mat,b) elif name == "convolve": nx = opt[0] ny = opt[1] sx = float(opt[2]) sy = float(opt[3]) nxd = (nx-1)/2 nyd = (ny-1)/2 # the kernel b = fromfunction(lambda j,i: exp(-(i-nxd)**2/sx**2 + -(j-nyd)**2/sy**2),(nx,ny)) s = sum(sum(b)) b = b/s # normalisation # conversion: b = b.astype(float) mat=mat.astype(float) return conv.convolve2d(mat,b,output=None,fft=0) if name == "boxcar": nx = opt[0] ny = opt[1] boxshape = (nx,ny) return conv.boxcar(mat,boxshape,mode='reflect') if name == "gaussian": sigma = float(opt[0]) return ndimage.gaussian_filter(mat,sigma,mode='nearest') if name == "uniform": sigma = float(opt[0]) return ndimage.uniform_filter(mat,sigma,mode='nearest') elif name == None: pass else: print "unknown filter type" return mat def set_ranges(mat,scale='log',mn=None,mx=None,cd=None): ''' Transform an n x m float array into an n x m int array that will be used to create an image. The float values are rescaled and cutted in order to range between 0 and 255. mat : the matrice scale : lin or log mn : lower value for the cutoff mx : higer value for the cutoff cd : parameter ''' rm = ravel(mat) if mn == None: mn = min(rm) if mx == None: mx = max(rm) if mn == mx: mn = min(rm) mx = max(rm) mat = clip(mat,mn,mx) if scale == 'log': if cd == None or cd == 0: cd = rm.mean() try: mat = 255.*log(1.+(mat-mn)/(cd)) / log(1.+(mx-mn)/(cd)) except: mat = mat*0. elif scale == 'lin': mat = 255.*(mat-mn)/(mx-mn) cd = 0 return mat.astype(int),mn,mx,cd def contours(m,matint,nl,mn,mx,kx,ky,color,crush='no'): ''' Compute iso-contours on a n x m float array. If "l_min" equal "l_max", levels are automatically between the minimum and maximum values of the matrix "mat". m = matrice (real values) matint = matrice (interger values) kx = num of sub-boxes ky = num of sub-boxes nl = # of levels mn = min mx = max color = color of contours ''' # create an empty matrice newm = zeros(m.shape,float32) if color != 0: # transform color color = array(color,int8) if mx == mn: rm = ravel(m) mn = min(rm) mx = max(rm) levels = arange(mn+(mx-mn)/nl,mx,(mx-mn)/(nl+1)) else: levels = arange(mn,mx+(mx-mn)/(nl-1),(mx-mn)/(nl-1))[:nl] #print levels m = transpose(m) # !!! X = [(),(),(),(),()] rect = [(1,2),(2,3),(3,4),(4,1)] # number of sub-boxes per axis nx = (m.shape[0]-1)/(kx-1) ny = (m.shape[1]-1)/(ky-1) for i in range(0,nx): for j in range(0,ny): ix = (kx-1)*i # here, we could add an offset iy = (ky-1)*j X[1]=(ix, iy, m[iy ][ix]) X[2]=(ix+(kx-1),iy, m[iy ][ix+(kx-1)]) X[3]=(ix+(kx-1),iy+(ky-1),m[iy+(ky-1)][ix+(kx-1)]) X[4]=(ix, iy+(ky-1),m[iy+(ky-1)][ix]) # find the center X[0]=(ix+0.5*kx,iy+0.5*ky,0.25*(X[1][2]+X[2][2]+X[3][2]+X[4][2])) # loop over the levels for l in levels: # loop over the trianles for r in rect: z1 = X[r[0]][2] z2 = X[r[1]][2] z3 = X[0][2] # find the maximum if z1 > z2: if z1 > z3: if z2 > z3: c = X[r[0]] b = X[r[1]] a = X[0] else: c = X[r[0]] b = X[0] a = X[r[1]] else: c = X[0] b = X[r[0]] a = X[r[1]] else: if z2 > z3: if z1 > z3: c = X[r[1]] b = X[r[0]] a = X[0] else: c = X[r[1]] b = X[0] a = X[r[0]] else: c = X[0] b = X[r[1]] a = X[r[0]] # a,b,c are the tree triangle points #create_line(newm,a[0],a[1],b[0],b[1],color) #create_line(newm,b[0],b[1],c[0],c[1],color) #create_line(newm,c[0],c[1],a[0],a[1],color) if l >= a[2] and l <= c[2] and a[2] != c[2] : # the level cut the triangle lamb = (l-a[2])/(c[2]-a[2]) xx1 = int(lamb*(c[0]-a[0]) + a[0]) yy1 = int(lamb*(c[1]-a[1]) + a[1]) if l >= b[2] and l <= c[2] and b[2] != c[2]: lamb = (l-b[2])/(c[2]-b[2]) xx2 = int(lamb*(c[0]-b[0]) + b[0]) yy2 = int(lamb*(c[1]-b[1]) + b[1]) elif b[2] != a[2] : lamb = (l-a[2])/(b[2]-a[2]) xx2 = int(lamb*(b[0]-a[0]) + a[0]) yy2 = int(lamb*(b[1]-a[1]) + a[1]) create_line(newm,xx1,yy1,xx2,yy2,color) matc = newm.astype(int8) if crush == 'yes': matint = ones(matc.shape) matint = matint*255 matint = where(matc,matc,matint) else: matint = where(matc,matc,matint) return matint def add_box(matint,shape=(256,256),size=(30.,30.),center=None,box_opts=(1,None,None,255)): ''' add a box on the frame ''' lweight = box_opts[0] xticks = box_opts[1] yticks = box_opts[2] color = box_opts[3] box = sbox(shape,size,lweight=lweight,xticks=xticks,yticks=yticks,color=color) matint = where(box!=0,box,matint) return matint def mplot(mat,palette='light',save=None,marker=None,header=None): ''' plot a 2d array ''' from pNbody import io if mpi.mpi_IsMaster(): if save!=None: if os.path.splitext(save)[1] == ".fits": io.WriteFits(transpose(mat).astype(float32),save, header) return # if the matrix is real if mat.dtype != int32: matint,mn_opt,mx_opt,cd_opt = set_ranges(mat,scale='lin',cd=0,mn=0,mx=0) mat = matint # add marker if marker!=None: if marker=='cross': nx,ny = mat.shape ix,iy = indices(mat.shape) mat = where((ix==nx/2)+(ix==nx/2-1)+(iy==ny/2)+(iy==ny/2-1),255,mat) if marker=='circle': nx,ny = mat.shape ix,iy = indices(mat.shape) r = 100 dr = 1 Rs = sqrt((ix-nx/2)**2 + (iy-ny/2)**2) c = (Rsr-dr) mat = where(c,255,mat) image = get_image(mat,name=None,palette_name=palette) if save==None: display(image) else: image.save(save) def get_image(mat,name=None,palette_name='light',mode='RGB'): ''' Return an image (PIL object). data : numpy 2x2 object name : name of the output palette_name : name of a palette ''' # modif 09.03.05 mat = transpose(mat) mat = mat.astype(int8) #mat = mat.astype(int) # ok #mat = transpose(mat) # la transposee fait changer aussi autre chose ??? c'est peut-etre la sortie de Py_BuildValue("(Of)",mat,cdopt); pas bon... - image = Image.fromstring("P",(mat.shape[1],mat.shape[0]),mat.tostring()) + image = Image.frombytes("P",(mat.shape[1],mat.shape[0]),mat.tostring()) # include the palette palette = Palette(palette_name) image.putpalette(palette.palette) # set mode if mode=='RGB': # to convert with ImageQt, need to be in RGB image = image.convert('RGB') # save it if name != None: image.save(name) return image def display(image): if mpi.mpi_IsMaster(): if is_qt: libqt.display(image) elif is_tk: #root = tk.Tk() root = tk.Toplevel() canvas = tk.Canvas(root) canvas.pack() imagetk = ImageTk.PhotoImage(image) canvas.config(width=image.size[0],height=image.size[1]) canvas.create_image(0.,0.,anchor=tk.NW,image=imagetk) root.mainloop() else: raise "tk or qt not present" def sbox(shape,size,lweight=1,xticks=None,yticks=None,color=255): ''' simple box return a matrix of integer, containing a box with labels xticks = (m0,d0,h0,m1,d1,h1) 0 = big 1 = small m0,m1 = dist between ticks d0,d1 = first tick h0,h1 = height of the ticks ''' center = None # parameters nn = 4. bticks = array([1.,2.,5.]) color = int(color) # create matrix from scratch matint = zeros(shape,int) # add the outer box for l in range(lweight): # in x matint[:,l]=(zeros((shape[0],))+color).astype(int) matint[:,shape[1]-1-l]=(zeros((shape[0],))+color).astype(int) # may be commented # in y matint[l,:]=(zeros((shape[1],))+color).astype(int) matint[shape[0]-1-l,:]=(zeros((shape[1],))+color).astype(int) # may be commented # add the ticks in x if xticks == None: #cx = center[0] #cy = center[1] cx = 0 cy = 0 hx = shape[0] hy = shape[1] dx = size[0] dy = size[1] rat = (2.*dx)/nn # dist between ticks ratn = rat / 10**int(log10(rat)) d0 = bticks[argmin(fabs(bticks-ratn))] * 10**int(log10(rat)) n0 = int(fmod(2.*dx/d0,2)+(2.*dx/d0)) # number of ticks m0 = cx - (n0*d0)/2. # first tick h0 = hy/20. # height of the ticks d1 = d0/4. # dist between ticks n1 = n0*4 m1 = m0 # first tick h1 = h0/2. # height of the ticks else: m0 = xticks[0] d0 = xticks[1] h0 = xticks[2] m1 = xticks[3] d1 = xticks[4] h1 = xticks[5] # big ticks x matint = drawxticks(matint,m0,d0,n0,h0,shape,size,center,color) # small ticks x matint = drawxticks(matint,m1,d1,n1,h1,shape,size,center,color) #print "dx = %8.3f"%(d0) # add the ticks in y if xticks == None: #cx = center[0] #cy = center[1] cx = 0 cy = 0 hx = shape[0] hy = shape[1] dx = size[0] dy = size[1] rat = (2.*dx)/nn # dist between ticks ratn = rat / 10**int(log10(rat)) d0 = bticks[argmin(fabs(bticks-ratn))] * 10**int(log10(rat)) n0 = int(fmod(2.*dx/d0,2)+(2.*dx/d0)) # number of ticks m0 = cy - (n0*d0)/2. # first tick h0 = hx/20. # height of the ticks d1 = d0/4. # dist between ticks n1 = n0*4 m1 = m0 # first tick h1 = h0/2. # height of the ticks else: m0 = xticks[0] d0 = xticks[1] h0 = xticks[2] m1 = xticks[3] d1 = xticks[4] h1 = xticks[5] # big ticks y matint = drawyticks(matint,m0,d0,n0,h0,shape,size,center,color) # small ticks y matint = drawyticks(matint,m1,d1,n1,h1,shape,size,center,color) #print "dy = %8.3f"%(d0) return matint def drawxticks(matint,m0,d0,n0,h0,shape,size,center,color): ''' draw x ticks in a matrix ''' x = m0 for nx in range(n0): ix,iy = phys2img(shape,size,center,x,0) try: # !! version dependant !! ix = int(ix[0]) except: ix = int(ix) for iy in range(int(h0)): matint[ix,iy]=array([color],int)[0] matint[ix,shape[0]-1-iy]=array([color],int)[0] x = x + d0 return matint def drawyticks(matint,m0,d0,n0,h0,shape,size,center,color): ''' draw x ticks in a matrix ''' x = m0 for nx in range(n0): ix,iy = phys2img(shape,size,center,x,0) try: # !! version dependant !! ix = int(ix[0]) except: ix = int(ix) for iy in range(int(h0)): matint[iy,ix]=array([color],int)[0] matint[shape[1]-1-iy,ix]=array([color],int)[0] x = x + d0 return matint def draw_line(m,x0,x1,y0,y1,c,z0=None,z1=None): dd = 0.5 # half a pixel imax = m.shape[0]-1 jmax = m.shape[1]-1 if abs(x0-x1) > abs(y0-y1) : ixmin = int(round(min(x0,x1))) ixmax = int(round(max(x0,x1))) a = (y1-y0)/(x1-x0) #az = (z1-z0)/(x1-x0) for i in xrange(ixmin,ixmax+1): j = round((i - x0)*a + y0) #jz = (i - x0)*az + z0 if (i<=imax) and (i>=0) and (j<=jmax) and (j>=0): m[i,j] = c elif abs(x0-x1) <= abs(y0-y1) and abs(y0-y1) > 0: iymin = int(round(min(y0,y1))) iymax = int(round(max(y0,y1))) a = (x1-x0)/(y1-y0) #az = (z1-z0)/(y1-y0) for j in xrange(iymin,iymax+1): i = round((j - y0)*a + x0) #iz = (j - y0)*az + z0 if (i<=imax) and (i>=0) and (j<=jmax) and (j>=0): m[i,j] = c else: i = round(x0) j = round(y0) if (i<=imax) and (i>=0) and (j<=jmax) and (j>=0): m[i,j] = c return m def draw_points(m,x,y,c): n = len(x) for i in xrange(n): m = draw_line(m,x[i],x[i],y[i],y[i],c) return m def draw_lines(m,x,y,c): n = len(x) for i in xrange(n-1): m = draw_line(m,x[i],x[i+1],y[i],y[i+1],c) return m def draw_polygon(m,x,y,c): n = len(x) for i in xrange(n-1): m = draw_line(m,x[i],x[i+1],y[i],y[i+1],c) m = draw_line(m,x[n-1],x[0],y[n-1],y[0],c) return m def draw_segments(m,x,y,c,zp=None): n = len(x) p = int(n/2) for j in xrange(p): i = 2*j if zp!=None: m = draw_line(m,x[i],x[i+1],y[i],y[i+1],c,zp[i],zp[i+1]) else: m = draw_line(m,x[i],x[i+1],y[i],y[i+1],c) return m def draw_polygonN(m,x,y,c,N): n = len(x) p = int(n/N) for j in xrange(p): i = N*j for k in xrange(N-1): m = draw_line(m,x[i+k],x[i+k+1],y[i+k],y[i+k+1],c) m = draw_line(m,x[i+N-1],x[i],y[i+N-1],y[i],c) return m def phys2img(shape,size,center,x,y): ''' convert physical position into the image pixel ''' #cx = center[0] #cy = center[1] cx = 0 cy = 0 hx = shape[0] hy = shape[1] dx = size[0] dy = size[1] ax = hx/(2.*dx) bx = (hx/2.)*(1.-cx/dx) ay = hy/(2.*dy) by = (hy/2.)*(1.-cy/dy) ix = int(ax*x+bx) iy = int(ay*y+by) ix = clip(ix,0,hx-1) iy = clip(iy,0,hy-1) return ix,iy ################################# def can_to_carth(l,scale_fact,x_can,y_can): ################################# x = float(x_can - l[0]/2.)/scale_fact[0] y = float(-y_can + l[1]/2.)/scale_fact[1] return x,y ################################# def carth_to_can(l,scale_fact,x,y): ################################# x_can = int( x*scale_fact[0] + l[0]/2.) y_can = int(-y*scale_fact[1]+ l[1]/2.) return x_can,y_can ################################# def getval(nb,mode='m',obs=None): ################################# """ For each point, return a specific value linked to this point modes : ----- 0 : moment 0 m : moment 0 x : first moment in x y : first moment in y z : first moment in z y2 : second moment in x y2 : second moment in y z2 : second moment in z vx : first velocity moment in x vy : first velocity moment in y vz : first velocity moment in z vx2 : second velocity moment in x vy2 : second velocity moment in y vz2 : second velocity moment in z Lx : kinetic momemtum in x Ly : kinetic momemtum in y Lz : kinetic momemtum in z lx : specific kinetic momemtum in x ly : specific kinetic momemtum in y lz : specific kinetic momemtum in z u : specific energy rho : density T : temperature A : entropy P : pressure Tcool : cooling time Lum : luminosity Ne : local electro density # depends on projection r : first momemtum of radial distance r2 : second momemtum of radial distance vr : first momemtum of radial velocity vr2 : second momemtum of radial velocity vxyr : first momemtum of radial velocity in the plane vxyr2 : second momemtum of radial velocity in the plane vtr : first momemtum of tangential velocity in the plane vtr2 : second momemtum of tangential velocity in the plane 3 types of values ----------------- (1) scalar values pos,m,tem,u,rho,ne vx,vx2,vy,vy2,vz,vz2,vxyr,vxyr2,vxyt,vxyt2 (2) values computed with respect to xp z,z2 Lx,Ly,Lz,lx,ly,ly (3) values computed with respect to absolute position in space empty for the moment (do not have the initial positions...) """ ###################################### # values independent of angle of view ###################################### ########### # scalar ########### # moment 0 if mode =='0': val = ones(len(nb.pos)) # moment 0 elif mode =='m': val = ones(len(nb.pos)) # u elif mode=='u': val = nb.U() # rho elif mode=='rho': val = nb.Rho() # T elif mode=='T': val = nb.T() # A elif mode=='A': val = nb.A() # P elif mode=='P': val = nb.P() # Tcool elif mode=='Tcool': val = nb.Tcool() # Ne elif mode=='Ne': val = nb.Ne() # Hsml elif mode=='Hsml': val = nb.Hsml() # Lum elif mode=='Lum': val = nb.Lum() ########### # vectors ########### # x elif mode=='x': val = nb.pos[:,0] # y elif mode=='y': val = nb.pos[:,1] # z elif mode=='z': val = nb.pos[:,2] # x2 elif mode=='x2': val = nb.pos[:,0]**2 # y2 elif mode=='y2': val = nb.pos[:,1]**2 # z2 elif mode=='z2': val = nb.pos[:,2]**2 # vx elif mode=='vx': val = nb.vel[:,0] # vy elif mode=='vy': val = nb.vel[:,1] # vz elif mode=='vz': val = nb.vel[:,2] # vx2 elif mode=='vx2': val = nb.vel[:,0]**2 # vy2 elif mode=='vy2': val = nb.vel[:,1]**2 # vz2 elif mode=='vz2': val = nb.vel[:,2]**2 # kinetic momentum in x elif mode=='Lx': val = nb.mass*(nb.pos[:,1]*nb.vel[:,2] - nb.pos[:,2]*nb.vel[:,1]) # kinetic momentum in y elif mode=='Ly': val = nb.mass*(nb.pos[:,2]*nb.vel[:,0] - nb.pos[:,0]*nb.vel[:,2]) # kinetic momentum in z elif mode=='Lz': val = nb.mass*(nb.pos[:,0]*nb.vel[:,1] - nb.pos[:,1]*nb.vel[:,0]) # specific kinetic momentum in x elif mode=='lx': val = (nb.pos[:,1]*nb.vel[:,2] - nb.pos[:,2]*nb.vel[:,1]) # specific kinetic momentum in y elif mode=='ly': val = (nb.pos[:,2]*nb.vel[:,0] - nb.pos[:,0]*nb.vel[:,2]) # specific kinetic momentum in z elif mode=='lz': val = (nb.pos[:,0]*nb.vel[:,1] - nb.pos[:,1]*nb.vel[:,0]) # luminosity elif mode=='lum': val = nb.luminosity_spec() ######################################## # values dependent on the angle of view ######################################## # first moment in z elif mode=='r': # center xp nb.translate(obs[0]-obs[1]) val = nb.pos[:,2] # second moment in z elif mode=='r2': # center xp nb.translate(obs[0]-obs[1]) val = nb.pos[:,2]**2 # first moment in z elif mode=='vr': val = nb.vel[:,2] # second moment in z elif mode=='vr2': # center xp val = nb.vel[:,2]**2 # first velocity moment in radial velocity in the plane elif mode=='vxyr': val = (nb.pos[:,0]*nb.vel[:,0] + nb.pos[:,1]*nb.vel[:,1]) / sqrt(nb.pos[:,0]**2 + nb.pos[:,1]**2) # second velocity moment in radial velocity in the plane elif mode=='vxyr2': val = ((nb.pos[:,0]*nb.vel[:,0] + nb.pos[:,1]*nb.vel[:,1]) / sqrt(nb.pos[:,0]**2 + nb.pos[:,1]**2))**2 # first moment in tangential velocity in the plane elif mode=='vtr': val = (nb.pos[:,0]*nb.vel[:,1] - nb.pos[:,1]*nb.vel[:,0] )/ sqrt(nb.pos[:,0]**2 + nb.pos[:,1]**2) # second moment in tangential velocity in the plane elif mode=='vtr2': val = ((nb.pos[:,0]*nb.vel[:,1] - nb.pos[:,1]*nb.vel[:,0] )/ sqrt(nb.pos[:,0]**2 + nb.pos[:,1]**2))**2 # other mode else: #print "getval : unknown mode %s"%(mode) cmd = "val = %s"%(mode) #print "trying command %s"%(cmd) exec(cmd) del nb return val.astype(float32) ################################# def getvaltype(mode='m'): ################################# """ list values that depends on projection """ if mode=='r': valtype = 'in projection' elif mode=='r2': valtype = 'in projection' elif mode=='vr': valtype = 'in projection' elif mode=='vr2': valtype = 'in projection' elif mode=='vxyr': valtype = 'in projection' elif mode=='vxyr2': valtype = 'in projection' elif mode=='vtr': valtype = 'in projection' elif mode=='vtr2': valtype = 'in projection' # other mode else: valtype = 'normal' return valtype ################################# def extract_parameters(arg,kw,defaultparams): ################################# """ this function extract parameters given to a function it returns a dictionary of parameters with respective value defaultparams : dictionary of default parameters """ params = {} if len(kw) == 0 and len(arg) >= 1: if type(arg[0]) == types.DictionaryType: params = arg[0] elif len(kw) >=1 and len(arg) >= 1: if type(arg[0]) == types.DictionaryType: params = arg[0] # add other keywords for key in kw.keys(): params[key]=kw[key] elif len(kw) >= 1 and len(arg) == 0: params = kw newparams = deepcopy(defaultparams) for key in params.keys(): if defaultparams.has_key(key): newparams[key] = params[key] return newparams ########################################################## # functions relative to mapping ########################################################## ################################# def GetNumberMap(pos,shape): ################################# ''' ''' val = ones(pos.shape).astype(float32) if len(shape) == 1: # compute zero momentum m0 = mkmap1d(pos,val,val,shape) elif len(shape) == 2: # compute zero momentum m0 = mkmap2d(pos,val,val,shape) elif len(shape) == 3: # compute zero momentum m0 = mkmap3d(pos,val,val,shape) else: return return m0 ################################# def GetMassMap(pos,mass,shape): ################################# ''' ''' val = ones(pos.shape).astype(float32) if len(shape) == 1: # compute zero momentum m0 = mkmap1d(pos,mass,val,shape) elif len(shape) == 2: # compute zero momentum m0 = mkmap2d(pos,mass,val,shape) elif len(shape) == 3: # compute zero momentum m0 = mkmap3d(pos,mass,val,shape) else: return return m0 ################################# def GetMeanValMap(pos,mass,val,shape): ################################# ''' ''' if len(shape) == 1: # compute zero momentum m0 = mkmap1d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap1d(pos,mass,val,shape) elif len(shape) == 2: # compute zero momentum m0 = mkmap2d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap2d(pos,mass,val,shape) elif len(shape) == 3: # compute zero momentum m0 = mkmap3d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap3d(pos,mass,val,shape) else: return # combi map return GetMeanMap(m0,m1) ################################# def GetSigmaValMap(pos,mass,val,shape): ################################# ''' ''' if len(shape) == 1: # compute zero momentum m0 = mkmap1d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap1d(pos,mass,val,shape) # compute second momentum m2 = mkmap1d(pos,mass,val*val,shape) elif len(shape) == 2: # compute zero momentum m0 = mkmap2d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap2d(pos,mass,val,shape) # compute second momentum m2 = mkmap2d(pos,mass,val*val,shape) elif len(shape) == 3: # compute zero momentum m0 = mkmap3d(pos,mass,ones(len(pos),float32),shape) # compute first momentum m1 = mkmap3d(pos,mass,val,shape) # compute second momentum m2 = mkmap3d(pos,mass,val*val,shape) else: return # combi map return GetSigmaMap(m0,m1,m2) ################################# def GetMeanMap(m0,m1): ################################# ''' Return a MeanMap using the 0 and 1 momentum m0 : zero momentum m1 : first momentum ''' m1 = where(m0==0,0,m1) m0 = where(m0==0,1,m0) mat = m1/m0 return mat ################################# def GetSigmaMap(m0,m1,m2): ################################# ''' Return a MeanMap using the 0 and 1 and 2 momentum m0 : zero momentum m1 : first momentum m2 : second momentum ''' m1 = where(m0==0,0,m1) m2 = where(m0==0,0,m2) m0 = where(m0==0,1,m0) mat = m2/m0 - (m1/m0)**2 mat = where((mat>0),sqrt(mat),0) return mat ########################################################## # extract from map functions ########################################################## def Extract1dMeanFrom2dMap(x,y,mass,val,kx,ky,nmin,momentum=0): ''' Extract the mean along one axis, from a 2d mean or sigma matrix x : pos in first dim. (beetween 0 and 1) y : pos in sec dim. (beetween 0 and 1) mass : mass of particles val : value to compute kx : number of bins in x ky : number of bins in y nmin : min number of particles needed to compute value momentum : 0,1,2 (-1=number) ''' shape = (kx,ky) # normalize values pos = transpose(array([x,y,y],float32)) mat_num = GetNumberMap(pos,shape) if momentum == -1: mat_val = GetNumberMap(pos,shape) if momentum == 0: mat_val = GetMassMap(pos,mass,shape) elif momentum == 1: mat_val = GetMeanValMap(pos,mass,val,shape) elif momentum == 2: mat_val = GetSigmaValMap(pos,mass,val,shape) ################ # 2d mean ################ mat_num = transpose(mat_num) mat_val = transpose(mat_val) m1 = sum(mat_num,0) m0 = sum(ones(mat_val.shape),0) vec_num = where((m0!=0),m1/m0,0) c = (mat_num>nmin) # *(mat_val!=0) m1 = sum(mat_val*c,0) m0 = sum(ones(mat_val.shape)*c,0) vec_sigma = where((m0!=0),m1/m0,0) return vec_sigma def get1dMeanFrom2dMap(mat_val,mat_num,nmin=32,axis=0): m1 = sum(mat_num,axis) m0 = sum(ones(mat_val.shape),axis) vec_num = where((m0!=0),m1/m0,axis) c = (mat_num>nmin) m1 = sum(mat_val*c,axis) m0 = sum(ones(mat_val.shape)*c,axis) vec_sigma = where((m0!=0),m1/m0,axis) return vec_sigma ########################################################## # filter ########################################################## ################################# def log_filter(x,xmin,xmax,xc,kx=1.0): ################################# ''' map a value between 0 and kx ''' if xc==0: return kx * (x-xmin)/(xmax-xmin) else: return kx * log(1+(x-xmin)/xc) / log(1+(xmax-xmin)/xc) ################################# def log_filter_inv(k,xmin,xmax,xc,kx=1.0): ################################# ''' map a value betwen xmin and xmax ''' if xc==0: return (k/kx*(xmax-xmin)) + xmin else: A = log(1+(xmax-xmin)/xc) return xc*(exp(A*k/kx)-1.0) + xmin ########################################################## # change of coordinate ########################################################## ###################### # cylindrical coord ###################### def vel_cyl2cart(pos=None,vel=None): """ Transform velocities in cylindrical coordinates vr,vt,vz into carthesian coodinates vx,vy,vz. Pos is the position of particles in cart. coord. Vel is the velocity in cylindrical coord. Return a 3xn float array. """ x = pos[:,0] y = pos[:,1] z = pos[:,2] vr = vel[:,0] vt = vel[:,1] vz = vel[:,2] r = sqrt(x**2+y**2) vx = where(r>0,(vr*x - vt*y)/r,0) vy = where(r>0,(vr*y + vt*x)/r,0) vz = vz return transpose(array([vx,vy,vz])).astype(float32) def vel_cart2cyl(pos,vel): """ Transform velocities in carthesian coordinates vx,vy,vz into cylindrical coodinates vr,vz,vz. Pos is the position of particles in cart. coord. Vel is the velocity in cart. coord. Return a 3xn float array. """ x = pos[:,0] y = pos[:,1] z = pos[:,2] vx = vel[:,0] vy = vel[:,1] vz = vel[:,2] r = sqrt(x**2+y**2+z**2) vr = where(r>0,(x*vx + y*vy)/r,0) vt = where(r>0,(x*vy - y*vx)/r,0) vz = vz return transpose(array([vr,vt,vz])).astype(float32) ########################################################## # gl2nbody ########################################################## def RotateAround(angle,axis,point,ObsM): ''' this should be C ''' # this work with OpenGL #Q = ones(16,float) #glLoadIdentity(); #glTranslated(point[0],point[1],point[2]); #glRotated(angle,axis[0],axis[1],axis[2]); #glTranslated(-point[0],-point[1],-point[2]); ##Q = glGetDoublev(GL_MODELVIEW_MATRIX); #glMultMatrixd(ObsM); #ObsM = glGetDoublev(GL_MODELVIEW_MATRIX); #return ravel(ObsM) angle = angle*pi/180 point = concatenate((point,array([0]))) M = ObsM M.shape = (4,4) # center point M = M-point # construction of the rotation matrix norm = sqrt(axis[0]**2 + axis[1]**2 + axis[2]**2) if norm == 0: return x sn = sin(-angle/2.) e0 = cos(-angle/2.) e1 = axis[0]*sn/norm e2 = axis[1]*sn/norm e3 = axis[2]*sn/norm a = zeros((4,4),float) a[0,0] = e0**2 + e1**2 - e2**2 - e3**2 a[1,0] = 2.*(e1*e2 + e0*e3) a[2,0] = 2.*(e1*e3 - e0*e2) a[3,0] = 0. a[0,1] = 2.*(e1*e2 - e0*e3) a[1,1] = e0**2 - e1**2 + e2**2 - e3**2 a[2,1] = 2.*(e2*e3 + e0*e1) a[3,1] = 0. a[0,2] = 2.*(e1*e3 + e0*e2) a[1,2] = 2.*(e2*e3 - e0*e1) a[2,2] = e0**2 - e1**2 - e2**2 + e3**2 a[3,2] = 0. a[0,3] = 0. a[1,3] = 0. a[2,3] = 0. a[3,3] = 1. a = a.astype(float) # multiply x and v M = dot(M,a) # decenter point M = M+point return ravel(M) def gl2pNbody(glparam,nbparam=None): EYE = 0 PTS = 4 HEA = 8 ARM = 12 gwinShapeX = 512; # to change... gwinShapeY = 512; ProjectionMode=glparam["ProjectionMode"] if (ProjectionMode): # frustum gwinPerspectiveTop = glparam["PerspectiveNear"]*tan(glparam["PerspectiveFov"]*0.5*pi/180.); gwinPerspectiveRight = gwinPerspectiveTop*float(gwinShapeX)/float(gwinShapeY); else: # ortho gwinPerspectiveLeft = -5*glparam["PerspectiveNear"]; gwinPerspectiveRight = -gwinPerspectiveLeft; gwinPerspectiveTop = float(gwinShapeY)/float(gwinShapeX)*gwinPerspectiveRight; gwinClip1 = glparam["PerspectiveNear"] gwinClip2 = glparam["PerspectiveFar"] #/*********************/ #/* observer position */ #/*********************/ # /* copy the observer matrix */ #ObsObject.ComputeEyes(); M = zeros(16,float) axis = zeros(3,float) point = zeros(3,float) M[0] = glparam["M0"] M[1] = glparam["M1"] M[2] = glparam["M2"] M[3] = glparam["M3"] M[4] = glparam["M4"] M[5] = glparam["M5"] M[6] = glparam["M6"] M[7] = glparam["M7"] M[8] = glparam["M8"] M[9] = glparam["M9"] M[10] = glparam["M10"] M[11] = glparam["M11"] M[12] = glparam["M12"] M[13] = glparam["M13"] M[14] = glparam["M14"] M[15] = glparam["M15"] # rotate axis[0] = M[EYE+0]-M[PTS+0]; axis[1] = M[EYE+1]-M[PTS+1]; axis[2] = M[EYE+2]-M[PTS+2]; point[0] = M[EYE+0]; point[1] = M[EYE+1]; point[2] = M[EYE+2]; M = RotateAround(90,axis,point,M); # this is also bad !!! # this is not correct, dist must come from projec. params */ dist = sqrt( pow(M[0]-M[8],2) + pow(M[1]-M[9],2) + pow(M[2]-M[10],2)) # head obs1 = M[0]; obs2 = M[1]; obs3 = M[2]; # lookat point */ norm = sqrt( pow(M[0]-M[4],2) + pow(M[1]-M[5],2) + pow(M[2]-M[6],2)); obs4 = obs1 + (M[4] -obs1)/norm*dist; obs5 = obs2 + (M[5] -obs2)/norm*dist; obs6 = obs3 + (M[6] -obs3)/norm*dist; # arm norm = sqrt( pow(M[0]-M[8],2) + pow(M[1]-M[9],2) + pow(M[2]-M[10],2)); obs7 = obs1 + (M[8] - obs1)/norm; obs8 = obs2 + (M[9] - obs2)/norm; obs9 = obs3 + (M[10]- obs3)/norm; # head norm = sqrt( pow(M[0]-M[12],2) + pow(M[1]-M[13],2) + pow(M[2]-M[14],2)); obs10 = obs1 + (M[12] - obs1)/norm; obs11 = obs2 + (M[13] - obs2)/norm; obs12 = obs3 + (M[14] - obs3)/norm; obs = array([obs1,obs2,obs3,obs4,obs5,obs6,obs7,obs8,obs9,obs10,obs11,obs12],float) obs.shape = (4,3) if nbparam == None: nbparam = param.Params(PARAMETERFILE,None) nbparam.set('obs',obs) nbparam.set('X0','None') nbparam.set('xp','None') nbparam.set('alpha','None') nbparam.set('view','xy') nbparam.set('r_obs',dist) nbparam.set('clip',(gwinClip1,gwinClip2)) nbparam.set('eye','None') nbparam.set('dist_eye',0.1) nbparam.set('foc',0.1) if ProjectionMode: nbparam.set('persp','on') nbparam.set('cut','yes') else: nbparam.set('persp','off') nbparam.set('cut','no') nbparam.set('shape',(gwinShapeX,gwinShapeY)) nbparam.set('size',(gwinPerspectiveRight,gwinPerspectiveTop)) return nbparam diff --git a/pNbody/palette.py b/pNbody/palette.py index 4e2830d..7d200e3 100644 --- a/pNbody/palette.py +++ b/pNbody/palette.py @@ -1,652 +1,652 @@ ''' this module is used to deal with color palettes. ''' try: import Tkinter as tk from PIL import ImageTk is_tk = True except ImportError: is_tk = False from numpy import * import os, string import glob # import parameters from parameters import * # !!! scipy import automatiquement Image qui entre en conflit avec Image # le mieux est donc de remplacer explicitement Image par Image de scipy (qui est le meme d'ailleurs) try: from scipy.pilutil import Image from scipy.interpolate import splrep,splev is_scipy = True except ImportError: from PIL import Image is_scipy = False #################################################################################### # read lut (palette) #################################################################################### def readlut(filename=os.path.join(PALETTEDIR,DEFAULTPALETTE)): ''' Read a lut file. ''' numByte = 256 file = open(filename, "r") data = file.read(numByte) data = file.read(numByte) # read red data = file.read(numByte) red = fromstring(data,'b') # read green data = file.read(numByte) green = fromstring(data,'b') # read blue data = file.read(numByte) blue = fromstring(data,'b') file.close() # combine the colors pal = chr(red[0]) + chr(green[0]) + chr(blue[0]) for i in range (1,256): pal = pal + chr(red[i]) + chr(green[i]) + chr(blue[i]) return pal def read_gimp_palette(name): f = open(name) f.readline() pal = '' for i in range(255): line = f.readline() line = string.split(line) r = int(line[0]) g = int(line[1]) b = int(line[2]) pal = pal + chr(r) + chr(g) + chr(b) return pal #################################################################################### # Create a palette #################################################################################### def getpalette(palettename): try: pal = read_gimp_palette(palettename) except: try: pal = readlut(palettename) except: pal = readlut() return(pal) ########################## class EditParams: ########################## def __init__(self,master, c,w,d,r): self.master = master self.root = Toplevel() self.frame = Frame(self.root) self.frame.grid(column=0,row=0) # lablel c self.labelc = Label(self.frame,text='valeur c = ') self.labelc.grid(column=0,row=0) # entry c self.entryc = Entry(self.frame) self.entryc.grid(column=1,row=0) self.entryc.insert(INSERT,"%s"%c) # lablel w self.labelw = Label(self.frame,text='valeur w = ') self.labelw.grid(column=0,row=1) # entry w self.entryw = Entry(self.frame) self.entryw.grid(column=1,row=1) self.entryw.insert(INSERT,"%s"%w) # lablel d self.labeld = Label(self.frame,text='valeur d = ') self.labeld.grid(column=0,row=2) # entry d self.entryd = Entry(self.frame) self.entryd.grid(column=1,row=2) self.entryd.insert(INSERT,"%s"%d) # lablel r self.labelr = Label(self.frame,text='valeur r = ') self.labelr.grid(column=0,row=3) # entry r self.entryr = Entry(self.frame) self.entryr.grid(column=1,row=3) self.entryr.insert(INSERT,"%s"%r) # buttons self.sendbutton = Button(self.frame,text='send',command=self.send) self.sendbutton.grid(column=0,row=4) self.okbutton = Button(self.frame,text='ok',command=self.ok) self.okbutton.grid(column=1,row=4) self.cancelbutton = Button(self.frame,text='cancel',command=self.cancel) self.cancelbutton.grid(column=2,row=4) def ok(self): try: c = int(self.entryc.get()) except: print "invalid value for c" try: w = int(self.entryw.get()) except: print "invalid value for w" try: d = float(self.entryd.get()) except: print "invalid value for d" try: r = int(self.entryr.get()) except: print "invalid value for r" # send values to main self.master.c = c self.master.w = w self.master.d = d self.master.invert = r self.master.get(c=c,w=w,d=d,invert=r) self.master.draw() self.root.destroy() def send(self): try: c = int(self.entryc.get()) except: print "invalid value for c" try: w = int(self.entryw.get()) except: print "invalid value for w" try: d = float(self.entryd.get()) except: print "invalid value for d" try: r = int(self.entryr.get()) except: print "invalid value for r" # send values to main self.master.c = c self.master.w = w self.master.d = d self.master.invert = r self.master.get(c=c,w=w,d=d,invert=r) self.master.draw() def cancel(self): self.root.destroy() #################################################################################### # Class Palette #################################################################################### class Palette: def __init__(self,name='light'): self.name = self.check_palette_name(name) self.Canvas = None self.tables = glob.glob(os.path.join(PALETTEDIR,'*')) self.c = 128 self.w = 256 self.d = 256. self.invert = 0 self.read(self.name) ########################### def check_palette_name(self,name): ########################### if os.path.isfile(name): pass else: name = os.path.join(PALETTEDIR,name) if not os.path.exists(name): print name,"do not exists, using %s instead"%(DEFAULTPALETTE) name = os.path.join(PALETTEDIR,DEFAULTPALETTE) return name ################################ def read(self,name): ################################ t = [] r = [] g = [] b = [] f = open(name) f.readline() pal = '' for i in range(256): line = f.readline() line = string.split(line) t.append(float(i)) r.append(float(line[0])) g.append(float(line[1])) b.append(float(line[2])) self.t = array(t,float) self.r = array(r,float) self.g = array(g,float) self.b = array(b,float) # store initial values self.t0 = self.t self.r0 = self.r self.g0 = self.g self.b0 = self.b self.mkspline() # spline should not be done here self.get(c=self.c,w=self.w,d=self.d) ################################ def write(self,name): ################################ f = open(name,'w') f.write("# table rgb (%s)\n"%name) for i in range(256): f.write("%3d %3d %3d\n"%(self.r[i],self.g[i],self.b[i])) f.close() ################################ def mkspline(self): ################################ if is_scipy: # splines self.ar = splrep(self.t,self.r,s=0) self.ag = splrep(self.t,self.g,s=0) self.ab = splrep(self.t,self.b,s=0) ################################ def append(self,p): ################################ ''' add a new palette ''' c = self.r + self.g + self.b self.r = where(c==0,p.r,self.r) self.g = where(c==0,p.g,self.g) self.b = where(c==0,p.b,self.b) ################################ def getr(self,vals): ################################ ''' return r value ''' vals = clip(vals,0,255) return self.r[vals] ################################ def getg(self,vals): ################################ ''' return g value ''' vals = clip(vals,0,255) return self.g[vals] ################################ def getb(self,vals): ################################ ''' return b value ''' vals = clip(vals,0,255) return self.b[vals] ################################ def get(self,c=0,w=0,d=256.,invert=0): ################################ ''' return the palette w : width of the interval [0,255] c : center of the interval [0,255] ''' if w==0 or c == 0: xmin = 0 xmax = 256 else: xmin = c - w/2 xmax = c + w/2 self.t = arange(xmin,xmax,(xmax-xmin)/256.) # fonction logarithmique pour exemple if d<256.: self.t = (log((self.t-xmin)/d+1) / log((xmax-xmin)/d+1) )*(xmax-xmin) + xmin if is_scipy : self.r = splev(self.t,self.ar) self.g = splev(self.t,self.ag) self.b = splev(self.t,self.ab) self.r = clip(self.r,0.,255.) self.g = clip(self.g,0.,255.) self.b = clip(self.b,0.,255.) # inversion if invert: self.r = 255. - self.r self.g = 255. - self.g self.b = 255. - self.b self.palette = '' for i in range(256): r = int(self.r[i]) g = int(self.g[i]) b = int(self.b[i]) self.palette = self.palette + chr(r) + chr(g) + chr(b) ################################ def setrange(self,mn,mx): ################################ ''' clip the palette between mn and mx NB: this function should be added to get mn = minimum mx = maximum ''' self.t = arange(0,255,255./(mx-mn)) if is_scipy : self.r = splev(self.t,self.ar) self.g = splev(self.t,self.ag) self.b = splev(self.t,self.ab) self.r = clip(self.r,0.,255.) self.g = clip(self.g,0.,255.) self.b = clip(self.b,0.,255.) # reconstuct the palette z = zeros(256,float) z[mn:mn+len(self.r)]=self.r self.r = z z = zeros(256,float) z[mn:mn+len(self.g)]=self.g self.g = z z = zeros(256,float) z[mn:mn+len(self.b)]=self.b self.b = z self.palette = '' for i in range(256): r = int(self.r[i]) g = int(self.g[i]) b = int(self.b[i]) self.palette = self.palette + chr(r) + chr(g) + chr(b) ########################### def addFrame(self,frame): ########################### if is_tk: ########################### # frame ########################### self.Frame = frame self.Frame.bind('',self.revert) # marche pas... self.Frame.bind('',self.edit) # marche pas... ########################### # canvas ########################### self.Canvas = tk.Canvas(frame,height=48,width=512) self.Canvas.grid(column=0,row=0) self.Canvas.bind('',self.motion) self.Canvas.bind('',self.move1) self.Canvas.bind('',self.move2) self.Canvas.bind('',self.move3) self.Canvas.bind('',self.mouse1) self.Canvas.bind('',self.edit) self.Canvas.bind('',self.mouse3) self.draw() else: print "addFrame : tk is not present" ########################### def change(self,name): ########################### self.name = self.check_palette_name(name) self.read(self.name) self.draw() ########################### def draw(self): ########################### if is_tk: palette_shape = (512,48) data = chr(0) for i in range (0,48): for j in range (0,256): data = data + chr(j) data = data + chr(j) data = data[1:] # converting data - image = Image.fromstring("P",palette_shape,data) + image = Image.frombytes("P",palette_shape,data) # include the palette image.putpalette(self.palette) # create a Tk photo self.palette_pic = ImageTk.PhotoImage(image) # insert photo in the pannel if self.Canvas != None: self.Canvas.create_image(0.,0.,anchor=tk.NW,image=self.palette_pic) else: print "draw : tk is not present" ########################### def edit(self,event): ########################### EditParams(self,self.c,self.w,self.d,self.invert) ########################### def mouse1(self,event): ########################### # find index in tables i = self.tables.index(self.name)+1 if i=0: self.change(self.tables[i]) ########################### def move1(self,event): ########################### self.c = self.c - (event.x - self.oevent.x) self.get(c=self.c,w=self.w,d=self.d) self.draw() self.oevent = event ########################### def move2(self,event): ########################### self.w = self.w + (event.x - self.oevent.x) self.get(c=self.c,w=self.w,d=self.d) self.draw() self.oevent = event ########################### def move3(self,event): ########################### self.d = self.d + (event.x - self.oevent.x)*(0.01*self.d) if self.d < 0. : self.d = 1e-3 if self.d > 256. : self.d = 256. self.get(c=self.c,w=self.w,d=self.d) self.draw() self.oevent = event ########################### def mouse3(self,event): ########################### self.c = 128 self.w = 256 self.d = 256. self.invert = 0 self.get(c=self.c,w=self.w,d=self.d) self.draw() self.oevent = event ########################### def motion(self,event): ########################### self.oevent = event ########################### def revert(self,event): ########################### print "revert"