From adba95c339f20334865d904a92aa749b570f377b Mon Sep 17 00:00:00 2001 From: =?utf8?q?Andr=C3=A9=20Wobst?= Date: Mon, 6 Feb 2006 19:22:08 +0000 Subject: [PATCH] fix the intersection bug and add some linear algebra functions to helper.py git-svn-id: https://pyx.svn.sourceforge.net/svnroot/pyx/trunk/pyx@2541 069f4177-920e-0410-937b-c2a4a81bcd90 --- CHANGES | 9 +- pyx/helper.py | 496 +++++++++++++++++++++++++++++++++++++++++++++----------- pyx/normpath.py | 113 +++++++++---- 3 files changed, 496 insertions(+), 122 deletions(-) rewrite pyx/helper.py (64%) diff --git a/CHANGES b/CHANGES index b76c820b..d8f7166b 100644 --- a/CHANGES +++ b/CHANGES @@ -2,7 +2,6 @@ PyX TODO: - normpath module: - - fix intersection (see test/functional) - style module: - make dash length available - text module: @@ -24,6 +23,7 @@ TODO: this depends on the plattform and is not only an issue within the graph system) - graph.style module: - add styles using curves instead of lines + - histogram doesn't stroke lines to the baseline for steps=0 when two subsequent values are equal - canvas module: - contructor should only take a texrunner keyword argument - set method??? @@ -69,6 +69,11 @@ TODO: graph data +TODO for 0.9: + - text module: + - add encoding; check whether encoding=ascii works well for the common existing cases + + 0.9 (2005/??/??): - most important changes (to be included in the release notes): - mathtree removal (warning about integer division) @@ -116,6 +121,7 @@ TODO: - correct closepath (do not invalidate currentpoint but set it to the beginning of the current subpath); structural rework of pathitems - calculate real bboxes for Bezier curves + - fix intersection due to non-linear parametrization of bezier curves - add rotate methods to path, normpath, normsubpath, and normsubpathitems - add flushskippedline to normsubpath - add arclentoparam to normsubpath and normsubpathitems @@ -142,6 +148,7 @@ TODO: certain PostScript files (reported by Gert Ingold) - helper module: - remove nodefault marker in favour of separate _marker classes + - add some linear algebra functions - dvifile module: - silently ignore TrueType fonts in font mapping files (reported by Gabriel Vasseur) diff --git a/pyx/helper.py b/pyx/helper.py dissimilarity index 64% index d500aa01..1fcdd2b3 100644 --- a/pyx/helper.py +++ b/pyx/helper.py @@ -1,93 +1,403 @@ -#!/usr/bin/env python -# -*- coding: ISO-8859-1 -*- -# -# -# Copyright (C) 2002-2004 Jörg Lehmann -# Copyright (C) 2002-2005 André Wobst -# Copyright (C) 2005 Michael Schindler -# -# This file is part of PyX (http://pyx.sourceforge.net/). -# -# PyX 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 2 of the License, or -# (at your option) any later version. -# -# PyX 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 PyX; if not, write to the Free Software -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA - - -def sign(x): - return (x >= 0) and 1 or -1 - -# XXX fallback for Numeric (eigenvalue computation) to be implemented along -# know algorithms (like from numerical recipes) - -import Numeric, LinearAlgebra - -def realpolyroots(coeffs, epsilon=1e-5, polish=0): - - """returns the roots of a polynom with given coefficients - - This helper routine uses the package Numeric to find the roots - of the polynomial with coefficients given in coeffs: - 0 = \sum_{i=0}^N x^i coeffs[i] - The solution is found via an equivalent eigenvalue problem - """ - - if not coeffs: - return [] - - try: - 1.0 / coeffs[-1] - except ZeroDivisionError: - roots = realpolyroots(coeffs[:-1], epsilon=epsilon) - else: - # divide the coefficients by the maximal order - # this is to be done in the matrix, anyhow and - # makes checking more stable - coeffs = [coeff/coeffs[-1] for coeff in coeffs] - - N = len(coeffs) - 1 - # build the Matrix of the polynomial problem - mat = Numeric.zeros((N, N), Numeric.Float) - for i in range(N-1): - mat[i+1][i] = 1 - for i in range(N): - mat[0][i] = -coeffs[N-1-i] - # find the eigenvalues of the matrix (== the roots of the polynomial) - roots = [complex(root) for root in LinearAlgebra.eigenvalues(mat)] - # take only the real roots - roots = [root.real for root in roots if -epsilon < root.imag < epsilon] - - # polish the roots with Newton-Raphson - if polish: - def polish(root, epsilon): - polynom = 2*epsilon - while abs(polynom) > epsilon: - polynom = coeffs[N]*root + coeffs[N-1] - poprime = coeffs[N]*N - for i in range(N-2,-1,-1): - polynom = polynom*root + coeffs[i] - poprime = poprime*root + coeffs[i+1]*(i+1) - root -= polynom / poprime - return root - - roots = [polish(root, epsilon) for root in roots] - - # # check if the roots are really roots! - # for root in roots: - # polynom = coeffs[N]*root + coeffs[N-1] - # for i in range(N-2,-1,-1): - # polynom = polynom*root + coeffs[i] - # if abs(polynom) > epsilon: - # raise Exception("value %f instead of 0" % polynom) - - return roots - +#!/usr/bin/env python +# -*- coding: ISO-8859-1 -*- +# +# +# Copyright (C) 2002-2004 Jörg Lehmann +# Copyright (C) 2002-2006 André Wobst +# Copyright (C) 2005 Michael Schindler +# +# This file is part of PyX (http://pyx.sourceforge.net/). +# +# PyX 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 2 of the License, or +# (at your option) any later version. +# +# PyX 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 PyX; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA + +import math, types + +def sign(x): + return (x >= 0) and 1 or -1 + +try: + import Numeric, LinearAlgebra + _has_numeric = 1 +except: + _has_numeric = 0 + +def evrealtrisym(d, e): + """returns eigenvalues of a tridiagonal real symmetric matrix + + d are the diagonal elements and e are the subdiagonal elements. + len(d) == len(e) + 1; d and e are destroyed + """ + n = len(d) + assert len(d) == len(e) + 1 + e = e + [0.0] + for l in range(n): + for iter in range(30): + for m in range(l, n-1): + dd = abs(d[m]) + abs(d[m+1]) + if abs(e[m]) + dd == dd: + break + else: + m = n-1 + if m == l: + break + g = (d[l+1]-d[l]) / (2.0*e[l]) + r = math.hypot(g, 1.0) + if g >= 0: + g = d[m]-d[l]+e[l]/(g+abs(r)) + else: + g = d[m]-d[l]+e[l]/(g-abs(r)) + s = c = 1.0 + p = 0.0 + for i in range(m, l, -1): + f = s*e[i-1] + b = c*e[i-1] + e[i] = r = math.hypot(f, g) + if r == 0.0: + d[i] -= p + e[m] = 0.0 + break + s = f/r + c = g/r + g = d[i]-p + r = (d[i-1]-g)*s + 2.0*c*b + p = s*r + d[i] = g+p + g = c*r-b + else: + d[l] -= p + e[l] = g + e[m] = 0.0 + else: + raise RuntimeError("Too many iterations in evrealtrisym") + return d + +def realsymtotrisym(a): + """creates a real tridiagonal matrix out of a real symmetric matrix + + Eigenvalues are not altered by the transformation; a is the matrix, i.e. a + list of lists but only the diagonal and off-diagonal elements left of the + diagonal are accessed (the other values don't need to be provided and the + lists might be shortend); a is destroyed + """ + n = len(a) + d = [None] * n + e = [None] * (n-1) + x = [None] * (n-1) + for i in range(n-1, 0, -1): + l = i + h = scale = 0.0 + if l > 1: + for k in range(l): + scale += abs(a[i][k]) + if scale == 0.0: + e[i-1] = a[i][l-1] + else: + for k in range(l): + a[i][k] /= scale + h += a[i][k]*a[i][k] + f = a[i][l-1] + g = math.sqrt(h) + if f >= 0: + g = -g + e[i-1] = scale*g + h -= f*g + a[i][l-1] = f-g + f = 0.0 + for j in range(l): + g = 0.0 + for k in range(j+1): + g += a[j][k]*a[i][k] + for k in range(j+1, l): + g += a[k][j]*a[i][k] + x[j] = g/h + f += x[j]*a[i][j] + hh = 0.5*f/h + for j in range(l): + f = a[i][j] + x[j] = g = x[j]-hh*f + for k in range(j+1): + a[j][k] -= f*x[k] + g*a[i][k] + else: + e[i-1] = a[i][l-1] + d[i] = h + for i in range(n): + d[i] = a[i][i] + return d, e + +def evrealsym(a): + """returns eigenvalues of a real symmetric matrix + + a is a real symmetric matrix, i.e. a list of lists but only the diagonal + and off-diagonal elements below the diagonal are accessed (other values are + not accessed, i.e. the lists might be shortend). + """ + # make a copy and ensure floats + a = [[float(x) for x in row] for row in a] + return evrealtrisym(*realsymtotrisym(a)) + +def balancerealsquare(a, radix=2): + """modify a general real square matrix to become balanced + + This routine replaces the general real square matrix a by a balanced matrix + with identical eigenvalues. A symmetric matrix is already balanced and is + unaffected by this procedure. The parameter radix should be the machine's + floating-point radix to prevent any numerical inaccuracies to be introduced + by this method. + """ + n = len(a) + radix2 = radix*radix + last = 0 + while not last: + last = 1 + for i in range(n): + r = c = 0.0 + for j in range(n): + if j+1 != i+1: + c += abs(a[j][i]) + r += abs(a[i][j]) + if c != 0 and r != 0: + g = r/radix + f = 1.0 + s = c+r + while c < g: + f *= radix + c *= radix2 + g = r*radix + while c > g: + f /= radix + c /= radix2 + if (c+r)/f < 0.95*s: + last = 0 + g = 1.0/f + for j in range(n): + a[i][j] *= g + for j in range(n): + a[j][i] *= f + +def hessrealsquare(a): + """modify a general real square matrix to the upper hessenberg form + + This routine replaces the general real square matrix a by a hessenberg + formed matrix with identical eigenvalues. + """ + n = len(a) + for m in range(1, n-1): + x = 0.0 + i = m + for j in range(m, n): + if (abs(a[j][m-1]) > abs(x)): + x = a[j][m-1] + i = j + if i != m: + for j in range(m-1, n): + a[i][j], a[m][j] = a[m][j], a[i][j] + for j in range(n): + a[j][i], a[j][m] = a[j][m], a[j][i] + if x: + for i in range(m+1, n): + y = a[i][m-1] + if y: + y /= x + a[i][m-1] = y + for j in range(m, n): + a[i][j] -= y*a[m][j] + for j in range(n): + a[j][m] += y*a[j][i] + for i in range(2, n): + for j in range(i-1): + a[i][j] = 0 + +def evhessrealsquare(a, realonly=0): + """returns eigenvalues of a real square matrix in upper hessenberg form + + a is a real square matrix in upper hessenberg form. a is destroyed. + The return value is a list containing a mixture of floats and complex + values. + """ + n = len(a) + e = [] + + anorm = 0.0 + for i in range(n): + for j in range(n): + anorm += abs(a[i][j]) + t = 0.0 + while n > 0: + its = 0 + first = 1 + while first or l < n-2: + first = 0 + for l in range(n-1, 0, -1): + s = abs(a[l-1][l-1]) + abs(a[l][l]) + if s: + s = anorm + if abs(a[l][l-1]) + s == s: + break + else: + l = 0 + x = a[n-1][n-1] + if l == n-1: + e.append(x+t) + n -= 1 + else: + y = a[n-2][n-2] + w = a[n-1][n-2]*a[n-2][n-1] + if l == n-2: + p = 0.5*(y-x) + q = p*p+w + z = math.sqrt(abs(q)) + x += t + if q >= 0.0: + if p >= 0: + z = p + abs(z) + else: + z = p - abs(z) + e.append(x+z) + if z: + e.append(x-w/z) + else: + e.append(x+z) + elif not realonly: + e.append(x+p - z*1J) + e.append(x+p + z*1J) + n -= 2 + else: + if its == 30: + raise RuntimeError("Too many iterations in evhessrealsquare") + if its == 10 or its == 20: + t += x + for i in range(n): + a[i][i] -= x + s = abs(a[n-1][n-2]) + abs(a[n-2][n-3]) + y = x = 0.75*s + w = -0.4375*s*s + its += 1 + for m in range(n-3, l-1, -1): + z = a[m][m] + r = x-z + s = y-z + p = (r*s-w)/a[m+1][m] + a[m][m+1] + q = a[m+1][m+1]-z-r-s + r = a[m+2][m+1] + s = abs(p)+abs(q)+abs(r) + p /= s + q /= s + r /= s + if m == l: + break + u = abs(a[m][m-1])*(abs(q)+abs(r)) + v = abs(p) * (abs(a[m-1][m-1]) + abs(z) + abs(a[m+1][m+1])) + if u + v == v: + break + for i in range(m+2, n): + a[i][i-2] = 0.0 + if i != m+2: + a[i][i-3]=0.0 + for k in range(m, n-1): + if k != m: + p = a[k][k-1] + q = a[k+1][k-1] + r = 0.0 + if k != n-2: + r = a[k+2][k-1] + x = abs(p) + abs(q) + abs(r) + if x: + p /= x + q /= x + r /= x + if p >= 0: + s = math.sqrt(p*p+q*q+r*r) + else: + s = -math.sqrt(p*p+q*q+r*r) + if s: + if k == m: + if l != m: + a[k][k-1] = -a[k][k-1] + else: + a[k][k-1] = -s*x + p += s + x = p/s + y = q/s + z = r/s + q /= p + r /= p + for j in range(k, n): + p = a[k][j]+q*a[k+1][j] + if k != n-2: + p += r*a[k+2][j] + a[k+2][j] -= p*z + a[k+1][j] -= p*y + a[k][j] -= p*x + mmin = min(n-1, k+3) + for i in range(l, mmin+1): + p = x*a[i][k] + y*a[i][k+1] + if k != n-2: + p += z*a[i][k+2] + a[i][k+2] -= p*r + a[i][k+1] -= p*q + a[i][k] -= p + else: + i = mmin + 1 + return e + +def evrealsquare(a, copy=1, realonly=0): + """returns eigenvalues of a real square matrix + + a is a real square matrix, i.e. a list of lists. The return value is a list + containing a mixture of floats and complex values. + """ + if copy: + # make a copy and ensure floats + a = [[float(x) for x in row] for row in a] + balancerealsquare(a) + hessrealsquare(a) + return evhessrealsquare(a, realonly=realonly) + +def realpolyroots(coeffs, epsilon=1e-5): + """returns the roots of a polynom with given coefficients + + polynomial with coefficients given in coeffs: + 0 = \sum_i x^i coeffs[i] + The solution is found via an equivalent eigenvalue problem + """ + if not coeffs: + return [] + try: + f = -1.0 / coeffs[-1] + except ZeroDivisionError: + return realpolyroots(coeffs[:-1]) + else: + n = len(coeffs) - 1 + if _has_numeric: + a = Numeric.zeros((n, n), Numeric.Float) + for i in range(n-1): + a[i+1][i] = 1 + for i in range(n): + a[0][i] = f*coeffs[n-1-i] + roots = [] + for root in LinearAlgebra.eigenvalues(a): + if type(root) == types.ComplexType: + if not root.imag: + roots.append(root.real) + else: + roots.append(root) + return roots + else: + a = [[0.0 for i in range(n)] for j in range(n)] + for i in range(n-1): + a[i+1][i] = 1.0 + for i in range(n): + a[0][i] = f*coeffs[n-1-i] + return evrealsquare(a, copy=0, realonly=1) diff --git a/pyx/normpath.py b/pyx/normpath.py index a204c356..1f32fc1a 100644 --- a/pyx/normpath.py +++ b/pyx/normpath.py @@ -4,7 +4,7 @@ # # Copyright (C) 2002-2005 Jörg Lehmann # Copyright (C) 2003-2004 Michael Schindler -# Copyright (C) 2002-2005 André Wobst +# Copyright (C) 2002-2006 André Wobst # # This file is part of PyX (http://pyx.sourceforge.net/). # @@ -32,7 +32,7 @@ except ImportError: def radians(x): return x*math.pi/180 def degrees(x): return x*180/math.pi -import bbox, canvas, path, trafo, unit +import bbox, canvas, helper, path, trafo, unit try: sum([]) @@ -406,31 +406,31 @@ class normcurve_pt(normsubpathitem): # bound for the length, and the length of the straight line # between start and end point of the normcurve (i.e. |P3-P1|), # which represents a lower bound. - upperlen1 = (math.hypot(x01_pt - self.x0_pt, y01_pt - self.y0_pt) + - math.hypot(x01_12_pt - x01_pt, y01_12_pt - y01_pt) + - math.hypot(xmidpoint_pt - x01_12_pt, ymidpoint_pt - y01_12_pt)) - lowerlen1 = math.hypot(xmidpoint_pt - self.x0_pt, ymidpoint_pt - self.y0_pt) - if upperlen1-lowerlen1 < epsilon: - c1 = normline_pt(self.x0_pt, self.y0_pt, xmidpoint_pt, ymidpoint_pt) + l0_pt = math.hypot(xmidpoint_pt - self.x0_pt, ymidpoint_pt - self.y0_pt) + l1_pt = math.hypot(x01_pt - self.x0_pt, y01_pt - self.y0_pt) + l2_pt = math.hypot(x01_12_pt - x01_pt, y01_12_pt - y01_pt) + l3_pt = math.hypot(xmidpoint_pt - x01_12_pt, ymidpoint_pt - y01_12_pt) + if l1_pt+l2_pt+l3_pt-l0_pt < epsilon: + a = leftnormline_pt(self.x0_pt, self.y0_pt, xmidpoint_pt, ymidpoint_pt, l1_pt, l2_pt, l3_pt) else: - c1 = normcurve_pt(self.x0_pt, self.y0_pt, - x01_pt, y01_pt, - x01_12_pt, y01_12_pt, - xmidpoint_pt, ymidpoint_pt) - - upperlen2 = (math.hypot(x12_23_pt - xmidpoint_pt, y12_23_pt - ymidpoint_pt) + - math.hypot(x23_pt - x12_23_pt, y23_pt - y12_23_pt) + - math.hypot(self.x3_pt - x23_pt, self.y3_pt - y23_pt)) - lowerlen2 = math.hypot(self.x3_pt - xmidpoint_pt, self.y3_pt - ymidpoint_pt) - if upperlen2-lowerlen2 < epsilon: - c2 = normline_pt(xmidpoint_pt, ymidpoint_pt, self.x3_pt, self.y3_pt) + a = leftnormcurve_pt(self.x0_pt, self.y0_pt, + x01_pt, y01_pt, + x01_12_pt, y01_12_pt, + xmidpoint_pt, ymidpoint_pt) + + l0_pt = math.hypot(self.x3_pt - xmidpoint_pt, self.y3_pt - ymidpoint_pt) + l1_pt = math.hypot(x12_23_pt - xmidpoint_pt, y12_23_pt - ymidpoint_pt) + l2_pt = math.hypot(x23_pt - x12_23_pt, y23_pt - y12_23_pt) + l3_pt = math.hypot(self.x3_pt - x23_pt, self.y3_pt - y23_pt) + if l1_pt+l2_pt+l3_pt-l0_pt < epsilon: + b = rightnormline_pt(xmidpoint_pt, ymidpoint_pt, self.x3_pt, self.y3_pt, l1_pt, l2_pt, l3_pt) else: - c2 = normcurve_pt(xmidpoint_pt, ymidpoint_pt, - x12_23_pt, y12_23_pt, - x23_pt, y23_pt, - self.x3_pt, self.y3_pt) + b = rightnormcurve_pt(xmidpoint_pt, ymidpoint_pt, + x12_23_pt, y12_23_pt, + x23_pt, y23_pt, + self.x3_pt, self.y3_pt) - return c1, c2 + return a, b def _arclentoparam_pt(self, lengths_pt, epsilon): a, b = self._midpointsplit(epsilon) @@ -439,9 +439,9 @@ class normcurve_pt(normsubpathitem): params = [] for param_a, param_b, length_pt in zip(params_a, params_b, lengths_pt): if length_pt > arclen_a_pt: - params.append(0.5+0.5*param_b) + params.append(b.convertparam(param_b)) else: - params.append(0.5*param_a) + params.append(a.convertparam(param_a)) return params, arclen_a_pt + arclen_b_pt def arclentoparam_pt(self, lengths_pt, epsilon): @@ -540,8 +540,8 @@ class normcurve_pt(normsubpathitem): a, b = self._midpointsplit(epsilon) # To improve the performance in the general case we alternate the # splitting process between the two normsubpathitems - return ( [( 0.5*a_t, o_t) for o_t, a_t in other.intersect(a, epsilon)] + - [(0.5+0.5*b_t, o_t) for o_t, b_t in other.intersect(b, epsilon)] ) + return ( [(a.convertparam(a_t), o_t) for o_t, a_t in other.intersect(a, epsilon)] + + [(b.convertparam(b_t), o_t) for o_t, b_t in other.intersect(b, epsilon)] ) def modifiedbegin_pt(self, x_pt, y_pt): return normcurve_pt(x_pt, y_pt, @@ -702,6 +702,63 @@ class normcurve_pt(normsubpathitem): return 6*self.y3_pt-18*self.y2_pt+18*self.y1_pt-6*self.y0_pt +# curve replacements used by midpointsplit: +# The replacements are normline_pt and normcurve_pt instances with an +# additional convertparam function taking care of the proper conversion +# the parametrization. Note that we need only one parametric conversion +# (when a parameter gets calculated), since in the other direction no +# midpointsplits are needed at all + +class leftnormline_pt(normline_pt): + + __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" + + def __init__(self, x0_pt, y0_pt, x1_pt, y1_pt, l1_pt, l2_pt, l3_pt): + normline_pt.__init__(self, x0_pt, y0_pt, x1_pt, y1_pt) + self.l1_pt = l1_pt + self.l2_pt = l2_pt + self.l3_pt = l3_pt + + def convertparam(self, param): + if 0 <= param <= 1: + params = helper.realpolyroots([-param*(self.l1_pt+self.l2_pt+self.l3_pt), + 3*self.l1_pt, + -3*self.l1_pt+3*self.l2_pt, + self.l1_pt-2*self.l2_pt+self.l3_pt]) + # we might get several solutions and choose the one closest to 0.5 + params.sort(lambda t1, t2: cmp(abs(t1-0.5), abs(t2-0.5))) + return 0.5*params[0] + else: + # when we are outside the proper parameter range, we skip the non-linear + # transformation, since it becomes slow and it might even start to be + # numerically instable + return 0.5*param + + +class rightnormline_pt(leftnormline_pt): + + __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" + + def convertparam(self, param): + return 0.5+leftnormline_pt.convertparam(self, param) + + +class leftnormcurve_pt(normcurve_pt): + + __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "x2_pt", "y2_pt", "x3_pt", "y3_pt" + + def convertparam(self, param): + return 0.5*param + + +class rightnormcurve_pt(normcurve_pt): + + __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "x2_pt", "y2_pt", "x3_pt", "y3_pt" + + def convertparam(self, param): + return 0.5+0.5*param + + ################################################################################ # normsubpath ################################################################################ -- 2.11.4.GIT