[PyX-user] MetaPost like Bezier curves

[kib2 <ki...@fr...>]

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*(sq5-1)*cost + 0.5*(3-sq5)*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[k-1] et z[k] (pour i dans [1=
,n])
      =20
        psi =3D [] # psi[k] =3D arg((z[k+1]-z[k])/(z[k]-z[k-1])) (pour i=20
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 =E0 n-1 si n est le nombre de points donnes
        l.append( 0 )
        for k in range(2):
            k +=3D1
            l.append( abs( z[k]-z[k-1]) )
        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