Thanks for your answers, the only problem I can see for the moment in the algorithm is that it uses the Numpy package for matrix computation and for solving a linear algebra system. ( with 2xn lines, where n is the initial number of points ). Calculus are made with complex numbers. Maybe that could be a problem if you wanted to put it inside PyX.
It was just for testing my understanding of the Hobby's article. In fact I don't have it, I asked on the list fr.comp.text.tex ( http://groups.google.fr/group/fr.comp.text.tex/browse_thread/thread/b0963fc8e430a22b/?hl=fr# ) if someone had a clue and I got the chance that J.C.Charpentier studied it a little before. He then explained to me the basics, because he doesn't remind all the formulas with the tensions and other options. Now, the bad part : - I found the Python's complex module rather light : isn't there anything to calculate the argument ( sorry, I don't know the English equivalent ) of a complex ? ( I've written my own function for this, but I may miss something I believe). It's the same for the exponential form of a complex "e^(j*theta) = cos(theta) + j*sin(theta)" . - I don't know Numpy at all, I just started to use it yesterday. So maybe the code looks ugly, the only reference I had is a table here : http://www.scipy.org/NumPy_for_Matlab_Users - I saw a bug once : when points are "badly" placed, the rendering is rather strange. But J.C.Charpentier told me it can be like this in MetaPost too, in fact I don't know. I'll test this week, but if some of you use MetaPost from time to time, maybe we can compare our results ? I'll try to do my best to post all this stuff tomorrow, see you : Kib² Here's a little version for 3 points, don't care about the geoObject. p0, p1, p2 are just Points instance with x and y attributes. Sorry, some comments are still in French but I don't have much time tonight. # La fonction f de John Hobby: def hobby(t, p): sq2, sq5 = math.sqrt(2), math.sqrt(5) sint, sinp = math.sin(t), math.sin(p) cost, cosp = math.cos(t), math.cos(p) k= 1.0/16.0 num = 2 + sq2*( sint - k*sinp)*(sinp - k*sint)*(cost - cosp) den = 3*( 1 + 0.5*(sq5-1)*cost + 0.5*(3-sq5)*cosp ) return float(num/den) # Complexe de la forme z = e^(i*theta) def e(theta): return complex( math.cos(theta) , math.sin(theta) ) # Calcule l'argument d'un complexe z donne def arg(z): # Computes the argument (in radians) of a given complex number rez,imz = z.real, z.imag if rez == 0.0: # on est dur l'axe des ordonnees if imz >0: return math.pi/2.0 if imz <0: return -math.pi/2.0 if imz == 0.0: return "Argument non defini pour un complexe nul \n" truc = math.atan(imz/rez) if rez > 0.0: return truc if rez < 0: return truc + math.pi class HobbyBezier( pyx.path.path, geoObject ): """ Trace une courbe de Bezier passant par 3 Points donnes. L'algorithme est de John Hobby dans le MetaFontBook. """ def __init__( self, p0, p1, p2): import numpy from numpy import linalg q = pyx.path.path() # z[] liste de complexes, on en compte n # pour z0..z1..z2 (n=2) z = [ complex(p0.x,p0.y) , complex(p1.x,p1.y) , complex(p2.x,p2.y) ] # On decoupe l'intervalle de temps t en n morceaux l = [] # l[k] = distance entre z[k-1] et z[k] (pour i dans [1,n]) psi = [] # psi[k] = arg((z[k+1]-z[k])/(z[k]-z[k-1])) (pour i dans [1,n-1], ici dans [1,2]) psi.append(0) psi.append( arg((z[2]-z[1])/(z[1]-z[0])) ) # k varie de 0 à n-1 si n est le nombre de points donnes l.append( 0 ) for k in range(2): k +=1 l.append( abs( z[k]-z[k-1]) ) print "longueurs =",l[1],l[2] theta = [] # n inconnues phi = [] # n inconnues a = numpy.mat([[0.,1.,1.,0.],[l[2],2.*l[1],-2.*l[2],-l[1]],[1.,0.,-1.,0.],[0.,1.,0.,-1.]]) print "a=\n",a,"\n" b = numpy.mat([ [-psi[1]],[0.],[0.],[0.] ]) sol = linalg.solve(a,b) phi.append(0) # phi[0] does not count for i in range(4): if i < 2: theta.append( sol[i,0] ) else: phi.append( sol[i,0] ) # Control Points u1 = z[0] + e(theta[0])*(z[1]-z[0])*hobby( theta[0], phi[1] ) v1 = z[1] - e(-phi[1])*(z[1]-z[0])*hobby( phi[1], theta[0] ) u2 = z[1] + e(theta[1])*(z[2]-z[1])*hobby( theta[1], phi[2] ) v2 = z[2] - e(-phi[2])*(z[2]-z[1])*hobby( phi[2], theta[1] ) self.liste_controles = [ Point(u1.real, u1.imag), Point(v1.real, v1.imag), Point(u2.real, u2.imag), Point(v2.real, v2.imag) ] q.append( pyx.path.moveto( z[0].real, z[0].imag ) ) q.append( pyx.path.curveto( u1.real, u1.imag, v1.real, v1.imag, z[1].real, z[1].imag ) ) q.append( pyx.path.curveto( u2.real, u2.imag, v2.real, v2.imag, z[2].real, z[2].imag ) ) q.normpath( epsilon=None ).path() pyx.path.path.__init__( self, *q.pathitems ) geoObject.__init__( self ) def control_points( self ): return self.liste_controles ------------------------------------------------------------------------- Take Surveys. Earn Cash. Influence the Future of IT Join SourceForge.net's Techsay panel and you'll get the chance to share your opinions on IT & business topics through brief surveys - and earn cash http://www.techsay.com/default.php?page=join.php&p=sourceforge&CID=DEVDEV _______________________________________________ PyX-user mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pyx-user
