Thanks for your answers,
the only problem I can see for the moment in the algorithm is that it=20
uses the Numpy package for matrix computation and for solving a linear=20
algebra system. ( with 2xn lines, where n is the initial number of=20
points ). Calculus are made with complex numbers. Maybe that could be a=20
problem if you wanted to put it inside PyX.
It was just for testing my understanding of the Hobby's article. In fact=20
I don't have it, I asked on the list fr.comp.text.tex (=20
http://groups.google.fr/group/fr.comp.text.tex/browse_thread/thread/b0963=
fc8e430a22b/?hl=3Dfr#=20
) if someone had a clue and I got the chance that J.C.Charpentier=20
studied it a little before. He then explained to me the basics, because=20
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=20
anything to calculate the argument ( sorry, I don't know the English=20
equivalent ) of a complex ? ( I've written my own function for this, but=20
I may miss something I believe). It's the same for the exponential form=20
of a complex "e^(j*theta) =3D cos(theta) + j*sin(theta)" .
 I don't know Numpy at all, I just started to use it yesterday. So=20
maybe the code looks ugly, the only reference I had is a table here :=20
http://www.scipy.org/NumPy_for_Matlab_Users
 I saw a bug once : when points are "badly" placed, the rendering is=20
rather strange. But J.C.Charpentier told me it can be like this in=20
MetaPost too, in fact I don't know. I'll test this week, but if some of=20
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=B2
Here's a little version for 3 points, don't care about the geoObject.=20
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 tonig=
ht.
# La fonction f de John Hobby:
def hobby(t, p):
sq2, sq5 =3D math.sqrt(2), math.sqrt(5)
sint, sinp =3D math.sin(t), math.sin(p)
cost, cosp =3D math.cos(t), math.cos(p)
k=3D 1.0/16.0
num =3D 2 + sq2*( sint  k*sinp)*(sinp  k*sint)*(cost  cosp)
den =3D 3*( 1 + 0.5*(sq51)*cost + 0.5*(3sq5)*cosp )
return float(num/den)
# Complexe de la forme z =3D e^(i*theta)
def e(theta):
return complex( math.cos(theta) , math.sin(theta) )
=20
# Calcule l'argument d'un complexe z donne
def arg(z):
# Computes the argument (in radians) of a given complex number
rez,imz =3D z.real, z.imag
if rez =3D=3D 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 =3D=3D 0.0:
return "Argument non defini pour un complexe nul \n"
=20
truc =3D 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.
=20
L'algorithme est de John Hobby dans le MetaFontBook.
"""
def __init__( self, p0, p1, p2):
import numpy
from numpy import linalg
q =3D pyx.path.path()
# z[] liste de complexes, on en compte n
# pour z0..z1..z2 (n=3D2)
z =3D [ 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 =3D [] # l[k] =3D distance entre z[k1] et z[k] (pour i dans [1=
,n])
=20
psi =3D [] # psi[k] =3D arg((z[k+1]z[k])/(z[k]z[k1])) (pour i=20
dans [1,n1], ici dans [1,2])
psi.append(0)
psi.append( arg((z[2]z[1])/(z[1]z[0])) )
# k varie de 0 =E0 n1 si n est le nombre de points donnes
l.append( 0 )
for k in range(2):
k +=3D1
l.append( abs( z[k]z[k1]) )
print "longueurs =3D",l[1],l[2]
theta =3D [] # n inconnues
phi =3D [] # n inconnues =20
a =3D=20
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=3D\n",a,"\n"
b =3D numpy.mat([ [psi[1]],[0.],[0.],[0.] ])
sol =3D 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] )
=20
# Control Points
u1 =3D z[0] + e(theta[0])*(z[1]z[0])*hobby( theta[0], phi[1] )
v1 =3D z[1]  e(phi[1])*(z[1]z[0])*hobby( phi[1], theta[0] )
u2 =3D z[1] + e(theta[1])*(z[2]z[1])*hobby( theta[1], phi[2] )
v2 =3D z[2]  e(phi[2])*(z[2]z[1])*hobby( phi[2], theta[1] )
=20
self.liste_controles =3D [ Point(u1.real, u1.imag), Point(v1.real=
,=20
v1.imag), Point(u2.real, u2.imag), Point(v2.real, v2.imag) ]
=20
q.append( pyx.path.moveto( z[0].real, z[0].imag ) )
q.append( pyx.path.curveto( u1.real, u1.imag, v1.real, v1.imag,=20
z[1].real, z[1].imag ) )
q.append( pyx.path.curveto( u2.real, u2.imag, v2.real, v2.imag,=20
z[2].real, z[2].imag ) )
q.normpath( epsilon=3DNone ).path()
pyx.path.path.__init__( self, *q.pathitems )
geoObject.__init__( self )
=20
def control_points( self ):
return self.liste_controles
=20
