[PyX-checkins] pyx/pyx mathutils.py,NONE,1.1 connector.py,1.33,1.34 deformer.py,1.32,1.33 normpath.py,1.12,1.13 unit.py,1.53,1.54 helper.py,1.22,NONE From: André Wobst - 2006-02-10 16:48:31 Update of /cvsroot/pyx/pyx/pyx In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv14952/pyx Modified Files: connector.py deformer.py normpath.py unit.py Added Files: mathutils.py Removed Files: helper.py Log Message: move helper to mathutils; implement cubic and quartic polynom root solver Index: connector.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/connector.py,v retrieving revision 1.33 retrieving revision 1.34 diff -C2 -d -r1.33 -r1.34 *** connector.py 27 Sep 2005 08:36:57 -0000 1.33 --- connector.py 10 Feb 2006 16:48:21 -0000 1.34 *************** *** 24,28 **** import math from math import pi, sin, cos, atan2, tan, hypot, acos, sqrt ! import path, unit, helper, normpath try: from math import radians, degrees --- 24,28 ---- import math from math import pi, sin, cos, atan2, tan, hypot, acos, sqrt ! import path, unit, mathutils, normpath try: from math import radians, degrees *************** *** 137,141 **** else: radius = abs(0.5 * (bulge + 0.25 * distance**2 / bulge)) ! centerdist = helper.sign(bulge) * (radius - abs(bulge)) center = (0.5 * (self.box1.center[0] + self.box2.center[0]) - tangent[1]*centerdist, 0.5 * (self.box1.center[1] + self.box2.center[1]) + tangent[0]*centerdist) --- 137,141 ---- else: radius = abs(0.5 * (bulge + 0.25 * distance**2 / bulge)) ! centerdist = mathutils.sign(bulge) * (radius - abs(bulge)) center = (0.5 * (self.box1.center[0] + self.box2.center[0]) - tangent[1]*centerdist, 0.5 * (self.box1.center[1] + self.box2.center[1]) + tangent[0]*centerdist) Index: deformer.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/deformer.py,v retrieving revision 1.32 retrieving revision 1.33 diff -C2 -d -r1.32 -r1.33 *** deformer.py 12 Oct 2005 08:02:15 -0000 1.32 --- deformer.py 10 Feb 2006 16:48:21 -0000 1.33 *************** *** 24,28 **** import math, warnings ! import attr, helper, path, normpath, unit, color from path import degrees --- 24,28 ---- import math, warnings ! import attr, mathutils, path, normpath, unit, color from path import degrees *************** *** 108,112 **** # extreme case: curvA geometrically too big if fallback_threshold * abs(curvA*AB) > 1: ! a = math.sqrt(abs(D / (1.5 * curvA))) * helper.sign(D*curvA) b = (D - 1.5*curvA*a*abs(a)) / T return [(a, b)] --- 108,112 ---- # extreme case: curvA geometrically too big if fallback_threshold * abs(curvA*AB) > 1: ! a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA) b = (D - 1.5*curvA*a*abs(a)) / T return [(a, b)] *************** *** 114,118 **** # extreme case: curvB geometrically too big if fallback_threshold * abs(curvB*AB) > 1: ! b = math.sqrt(abs(E / (1.5 * curvB))) * helper.sign(E*curvB) a = (E - 1.5*curvB*b*abs(b)) / T return [(a, b)] --- 114,118 ---- # extreme case: curvB geometrically too big if fallback_threshold * abs(curvB*AB) > 1: ! b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB) a = (E - 1.5*curvB*b*abs(b)) / T return [(a, b)] *************** *** 140,144 **** # extreme case: T much smaller than both curvatures try: ! a = math.sqrt(abs(D / (1.5 * curvA))) * helper.sign(D*curvA) except ZeroDivisionError: # we have tried the case with small curvA already --- 140,144 ---- # extreme case: T much smaller than both curvatures try: ! a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA) except ZeroDivisionError: # we have tried the case with small curvA already *************** *** 147,151 **** else: try: ! b = math.sqrt(abs(E / (1.5 * curvB))) * helper.sign(E*curvB) except ZeroDivisionError: # we have tried the case with small curvB already --- 147,151 ---- else: try: ! b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB) except ZeroDivisionError: # we have tried the case with small curvB already *************** *** 170,177 **** for sign_a in [+1, -1]: for sign_b in [+1, -1]: ! coeffs_a = (sign_b*1.5*curvB*D*D - T*T*E, T**3, -sign_b*sign_a*4.5*curvA*curvB*D, 0.0, sign_b*3.375*curvA*curvA*curvB) ! coeffs_b = (sign_a*1.5*curvA*E*E - T*T*D, T**3, -sign_a*sign_b*4.5*curvA*curvB*E, 0.0, sign_a*3.375*curvA*curvB*curvB) ! candidates_a += [root for root in helper.realpolyroots(coeffs_a) if sign_a*root >= 0] ! candidates_b += [root for root in helper.realpolyroots(coeffs_b) if sign_b*root >= 0] # check the combined equations --- 170,177 ---- for sign_a in [+1, -1]: for sign_b in [+1, -1]: ! coeffs_a = (sign_b*3.375*curvA*curvA*curvB, 0.0, -sign_b*sign_a*4.5*curvA*curvB*D, T**3, sign_b*1.5*curvB*D*D - T*T*E) ! coeffs_b = (sign_a*3.375*curvA*curvB*curvB, 0.0, -sign_a*sign_b*4.5*curvA*curvB*E, T**3, sign_a*1.5*curvA*E*E - T*T*D) ! candidates_a += [root for root in mathutils.realpolyroots(*coeffs_a) if sign_a*root >= 0] ! candidates_b += [root for root in mathutils.realpolyroots(*coeffs_b) if sign_b*root >= 0] # check the combined equations *************** *** 301,305 **** radius = abs(unit.topt(self.radius)) turnangle = degrees(self.turnangle) ! sign = helper.sign(self.sign) cosTurn = math.cos(turnangle) --- 301,305 ---- radius = abs(unit.topt(self.radius)) turnangle = degrees(self.turnangle) ! sign = mathutils.sign(self.sign) cosTurn = math.cos(turnangle) *************** *** 514,519 **** # this results in a discontinuous curvature # (but the absolute value is still continuous) ! s1 = helper.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0])) ! s2 = helper.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0])) c1 = s1 * abs(c1) c2 = s2 * abs(c2) --- 514,519 ---- # this results in a discontinuous curvature # (but the absolute value is still continuous) ! s1 = mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0])) ! s2 = mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0])) c1 = s1 * abs(c1) c2 = s2 * abs(c2) *************** *** 977,981 **** dist = self.dist_pt # do not look closer than epsilon: ! dist_relerr = helper.sign(dist) * max(abs(self.relerr*dist), epsilon) checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams] --- 977,981 ---- dist = self.dist_pt # do not look closer than epsilon: ! dist_relerr = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon) checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams] Index: normpath.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/normpath.py,v retrieving revision 1.12 retrieving revision 1.13 diff -C2 -d -r1.12 -r1.13 *** normpath.py 6 Feb 2006 19:22:08 -0000 1.12 --- normpath.py 10 Feb 2006 16:48:21 -0000 1.13 *************** *** 33,37 **** def degrees(x): return x*180/math.pi ! import bbox, canvas, helper, path, trafo, unit try: --- 33,37 ---- def degrees(x): return x*180/math.pi ! import bbox, canvas, mathutils, path, trafo, unit try: *************** *** 412,421 **** 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: ! 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) --- 412,421 ---- 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: ! 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) *************** *** 424,433 **** 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: ! b = rightnormcurve_pt(xmidpoint_pt, ymidpoint_pt, ! x12_23_pt, y12_23_pt, ! x23_pt, y23_pt, ! self.x3_pt, self.y3_pt) return a, b --- 424,433 ---- 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: ! b = _rightnormcurve_pt(xmidpoint_pt, ymidpoint_pt, ! x12_23_pt, y12_23_pt, ! x23_pt, y23_pt, ! self.x3_pt, self.y3_pt) return a, b *************** *** 440,446 **** for param_a, param_b, length_pt in zip(params_a, params_b, lengths_pt): if length_pt > arclen_a_pt: ! params.append(b.convertparam(param_b)) else: ! params.append(a.convertparam(param_a)) return params, arclen_a_pt + arclen_b_pt --- 440,446 ---- for param_a, param_b, length_pt in zip(params_a, params_b, lengths_pt): if length_pt > arclen_a_pt: ! params.append(b.subparamtoparam(param_b)) else: ! params.append(a.subparamtoparam(param_a)) return params, arclen_a_pt + arclen_b_pt *************** *** 541,546 **** # To improve the performance in the general case we alternate the # splitting process between the two normsubpathitems ! 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): --- 541,546 ---- # To improve the performance in the general case we alternate the # splitting process between the two normsubpathitems ! return ( [(a.subparamtoparam(a_t), o_t) for o_t, a_t in other.intersect(a, epsilon)] + ! [(b.subparamtoparam(b_t), o_t) for o_t, b_t in other.intersect(b, epsilon)] ) def modifiedbegin_pt(self, x_pt, y_pt): *************** *** 705,714 **** # 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" --- 705,714 ---- # curve replacements used by midpointsplit: # The replacements are normline_pt and normcurve_pt instances with an ! # additional subparamtoparam function for proper conversion of the ! # parametrization. Note that we only one direction (when a parameter ! # gets calculated), since the other way around direction midpointsplit ! # is not needed at all ! class _leftnormline_pt(normline_pt): __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" *************** *** 720,730 **** 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] --- 720,733 ---- self.l3_pt = l3_pt ! def subparamtoparam(self, param): if 0 <= param <= 1: ! params = mathutils.realpolyroots(self.l1_pt-2*self.l2_pt+self.l3_pt, ! -3*self.l1_pt+3*self.l2_pt, ! 3*self.l1_pt, ! -param*(self.l1_pt+self.l2_pt+self.l3_pt)) # we might get several solutions and choose the one closest to 0.5 + # (we want the solution to be in the range 0 <= param <= 1; in case + # we get several solutions in this range, they all will be close to + # each other since l1_pt+l2_pt+l3_pt-l0_pt < epsilon) params.sort(lambda t1, t2: cmp(abs(t1-0.5), abs(t2-0.5))) return 0.5*params[0] *************** *** 736,760 **** ! 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 --- 739,763 ---- ! class _rightnormline_pt(_leftnormline_pt): __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" ! def subparamtoparam(self, param): ! return 0.5+_leftnormline_pt.subparamtoparam(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 subparamtoparam(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 subparamtoparam(self, param): return 0.5+0.5*param Index: unit.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/unit.py,v retrieving revision 1.53 retrieving revision 1.54 diff -C2 -d -r1.53 -r1.54 *** unit.py 8 Jul 2005 13:03:48 -0000 1.53 --- unit.py 10 Feb 2006 16:48:21 -0000 1.54 *************** *** 23,27 **** import types - import helper scale = { 't':1, 'u':1, 'v':1, 'w':1, 'x':1 } --- 23,26 ---- --- NEW FILE: mathutils.py --- #!/usr/bin/env python # -*- coding: ISO-8859-1 -*- # # # Copyright (C) 2005-2006 Michael Schindler # Copyright (C) 2006 André Wobst # # 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 try: import Numeric, LinearAlgebra _has_numeric = 1 except: _has_numeric = 0 def sign(x): """sign of x, i.e. +1 or -1; returns 1 for x == 0""" if x >= 0: return 1 return -1 def asinh(x): """Return the arc hyperbolic sine of x.""" return math.log(x+math.sqrt(x*x+1)) def acosh(x): """Return the arc hyperbolic cosine of x.""" return math.log(x+math.sqrt(x*x-1)) def _realroots_quadratic(a1, a0): """gives the real roots of x**2 + a1 * x + a0 = 0""" D = a1*a1 - 4*a0 if D < 0: return [] SD = math.sqrt(D) return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)] def _realroots_cubic(a2, a1, a0): """gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0""" # see http://mathworld.wolfram.com/CubicFormula.html for details Q = (3*a1 - a2*a2) / 9.0 R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0 D = Q*Q*Q + R*R if D > 0: # one real and two complex roots SD = math.sqrt(D) if R + SD >= 0: S = (R + SD)**(1/3.0) else: S = -(-R - SD)**(1/3.0) if R - SD >= 0: T = (R - SD)**(1/3.0) else: T = -(SD - R)**(1/3.0) return [S + T - a2/3.0] elif D == 0: if Q == 0: # one real root (R==0) return [-a2/3.0] else: # two real roots (R>0, Q<0) S = -math.sqrt(-Q) return [2*S - a2/3.0, -S - a2/3.0] else: # three real roots (Q<0) SQ = math.sqrt(-Q) theta = math.acos(R/(SQ**3)) return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)] def _realroots_quartic(a3, a2, a1, a0): """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0""" # see http://mathworld.wolfram.com/QuarticEquation.html for details ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3) ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and y*y-4*a0 >= 0] if not ys: return [] y1 = min(ys) if a3*y1-2*a1 < 0: return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) + _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0)))) else: return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) + _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0)))) def realpolyroots(*cs): """returns the roots of a polynom with given coefficients polynomial with coefficients given in cs: 0 = \sum_i cs[i] * x^(len(cs)-i-1) """ if not cs: return [0] try: f = 1.0/cs[0] cs = [f*coeff for coeff in cs[1:]] except ArithmeticError: return realpolyroots(*cs[1:]) else: n = len(cs) if n == 0: return [] elif n == 1: return [-cs[0]] elif n == 2: return _realroots_quadratic(*cs) elif n == 3: return _realroots_cubic(*cs) elif n == 4: return _realroots_quartic(*cs) else: raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4") def realpolyroots_eigenvalue(*cs): # as realpolyroots but using an equivalent eigenvalue problem # (this code is currently used for functional tests only) if not _has_numeric: raise RuntimeError("realpolyroots_eigenvalue depends on Numeric") if not cs: return [0] try: f = 1.0/cs[0] cs = [f*coeff for coeff in cs[1:]] except ArithmeticError: return realpolyroots_eigenvalue(*cs[1:]) else: if not cs: return [] n = len(cs) 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] = -cs[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 --- helper.py DELETED --- 
 [PyX-checkins] pyx/pyx mathutils.py,NONE,1.1 connector.py,1.33,1.34 deformer.py,1.32,1.33 normpath.py,1.12,1.13 unit.py,1.53,1.54 helper.py,1.22,NONE From: André Wobst - 2006-02-10 16:48:31 Update of /cvsroot/pyx/pyx/pyx In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv14952/pyx Modified Files: connector.py deformer.py normpath.py unit.py Added Files: mathutils.py Removed Files: helper.py Log Message: move helper to mathutils; implement cubic and quartic polynom root solver Index: connector.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/connector.py,v retrieving revision 1.33 retrieving revision 1.34 diff -C2 -d -r1.33 -r1.34 *** connector.py 27 Sep 2005 08:36:57 -0000 1.33 --- connector.py 10 Feb 2006 16:48:21 -0000 1.34 *************** *** 24,28 **** import math from math import pi, sin, cos, atan2, tan, hypot, acos, sqrt ! import path, unit, helper, normpath try: from math import radians, degrees --- 24,28 ---- import math from math import pi, sin, cos, atan2, tan, hypot, acos, sqrt ! import path, unit, mathutils, normpath try: from math import radians, degrees *************** *** 137,141 **** else: radius = abs(0.5 * (bulge + 0.25 * distance**2 / bulge)) ! centerdist = helper.sign(bulge) * (radius - abs(bulge)) center = (0.5 * (self.box1.center[0] + self.box2.center[0]) - tangent[1]*centerdist, 0.5 * (self.box1.center[1] + self.box2.center[1]) + tangent[0]*centerdist) --- 137,141 ---- else: radius = abs(0.5 * (bulge + 0.25 * distance**2 / bulge)) ! centerdist = mathutils.sign(bulge) * (radius - abs(bulge)) center = (0.5 * (self.box1.center[0] + self.box2.center[0]) - tangent[1]*centerdist, 0.5 * (self.box1.center[1] + self.box2.center[1]) + tangent[0]*centerdist) Index: deformer.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/deformer.py,v retrieving revision 1.32 retrieving revision 1.33 diff -C2 -d -r1.32 -r1.33 *** deformer.py 12 Oct 2005 08:02:15 -0000 1.32 --- deformer.py 10 Feb 2006 16:48:21 -0000 1.33 *************** *** 24,28 **** import math, warnings ! import attr, helper, path, normpath, unit, color from path import degrees --- 24,28 ---- import math, warnings ! import attr, mathutils, path, normpath, unit, color from path import degrees *************** *** 108,112 **** # extreme case: curvA geometrically too big if fallback_threshold * abs(curvA*AB) > 1: ! a = math.sqrt(abs(D / (1.5 * curvA))) * helper.sign(D*curvA) b = (D - 1.5*curvA*a*abs(a)) / T return [(a, b)] --- 108,112 ---- # extreme case: curvA geometrically too big if fallback_threshold * abs(curvA*AB) > 1: ! a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA) b = (D - 1.5*curvA*a*abs(a)) / T return [(a, b)] *************** *** 114,118 **** # extreme case: curvB geometrically too big if fallback_threshold * abs(curvB*AB) > 1: ! b = math.sqrt(abs(E / (1.5 * curvB))) * helper.sign(E*curvB) a = (E - 1.5*curvB*b*abs(b)) / T return [(a, b)] --- 114,118 ---- # extreme case: curvB geometrically too big if fallback_threshold * abs(curvB*AB) > 1: ! b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB) a = (E - 1.5*curvB*b*abs(b)) / T return [(a, b)] *************** *** 140,144 **** # extreme case: T much smaller than both curvatures try: ! a = math.sqrt(abs(D / (1.5 * curvA))) * helper.sign(D*curvA) except ZeroDivisionError: # we have tried the case with small curvA already --- 140,144 ---- # extreme case: T much smaller than both curvatures try: ! a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA) except ZeroDivisionError: # we have tried the case with small curvA already *************** *** 147,151 **** else: try: ! b = math.sqrt(abs(E / (1.5 * curvB))) * helper.sign(E*curvB) except ZeroDivisionError: # we have tried the case with small curvB already --- 147,151 ---- else: try: ! b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB) except ZeroDivisionError: # we have tried the case with small curvB already *************** *** 170,177 **** for sign_a in [+1, -1]: for sign_b in [+1, -1]: ! coeffs_a = (sign_b*1.5*curvB*D*D - T*T*E, T**3, -sign_b*sign_a*4.5*curvA*curvB*D, 0.0, sign_b*3.375*curvA*curvA*curvB) ! coeffs_b = (sign_a*1.5*curvA*E*E - T*T*D, T**3, -sign_a*sign_b*4.5*curvA*curvB*E, 0.0, sign_a*3.375*curvA*curvB*curvB) ! candidates_a += [root for root in helper.realpolyroots(coeffs_a) if sign_a*root >= 0] ! candidates_b += [root for root in helper.realpolyroots(coeffs_b) if sign_b*root >= 0] # check the combined equations --- 170,177 ---- for sign_a in [+1, -1]: for sign_b in [+1, -1]: ! coeffs_a = (sign_b*3.375*curvA*curvA*curvB, 0.0, -sign_b*sign_a*4.5*curvA*curvB*D, T**3, sign_b*1.5*curvB*D*D - T*T*E) ! coeffs_b = (sign_a*3.375*curvA*curvB*curvB, 0.0, -sign_a*sign_b*4.5*curvA*curvB*E, T**3, sign_a*1.5*curvA*E*E - T*T*D) ! candidates_a += [root for root in mathutils.realpolyroots(*coeffs_a) if sign_a*root >= 0] ! candidates_b += [root for root in mathutils.realpolyroots(*coeffs_b) if sign_b*root >= 0] # check the combined equations *************** *** 301,305 **** radius = abs(unit.topt(self.radius)) turnangle = degrees(self.turnangle) ! sign = helper.sign(self.sign) cosTurn = math.cos(turnangle) --- 301,305 ---- radius = abs(unit.topt(self.radius)) turnangle = degrees(self.turnangle) ! sign = mathutils.sign(self.sign) cosTurn = math.cos(turnangle) *************** *** 514,519 **** # this results in a discontinuous curvature # (but the absolute value is still continuous) ! s1 = helper.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0])) ! s2 = helper.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0])) c1 = s1 * abs(c1) c2 = s2 * abs(c2) --- 514,519 ---- # this results in a discontinuous curvature # (but the absolute value is still continuous) ! s1 = mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0])) ! s2 = mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0])) c1 = s1 * abs(c1) c2 = s2 * abs(c2) *************** *** 977,981 **** dist = self.dist_pt # do not look closer than epsilon: ! dist_relerr = helper.sign(dist) * max(abs(self.relerr*dist), epsilon) checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams] --- 977,981 ---- dist = self.dist_pt # do not look closer than epsilon: ! dist_relerr = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon) checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams] Index: normpath.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/normpath.py,v retrieving revision 1.12 retrieving revision 1.13 diff -C2 -d -r1.12 -r1.13 *** normpath.py 6 Feb 2006 19:22:08 -0000 1.12 --- normpath.py 10 Feb 2006 16:48:21 -0000 1.13 *************** *** 33,37 **** def degrees(x): return x*180/math.pi ! import bbox, canvas, helper, path, trafo, unit try: --- 33,37 ---- def degrees(x): return x*180/math.pi ! import bbox, canvas, mathutils, path, trafo, unit try: *************** *** 412,421 **** 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: ! 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) --- 412,421 ---- 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: ! 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) *************** *** 424,433 **** 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: ! b = rightnormcurve_pt(xmidpoint_pt, ymidpoint_pt, ! x12_23_pt, y12_23_pt, ! x23_pt, y23_pt, ! self.x3_pt, self.y3_pt) return a, b --- 424,433 ---- 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: ! b = _rightnormcurve_pt(xmidpoint_pt, ymidpoint_pt, ! x12_23_pt, y12_23_pt, ! x23_pt, y23_pt, ! self.x3_pt, self.y3_pt) return a, b *************** *** 440,446 **** for param_a, param_b, length_pt in zip(params_a, params_b, lengths_pt): if length_pt > arclen_a_pt: ! params.append(b.convertparam(param_b)) else: ! params.append(a.convertparam(param_a)) return params, arclen_a_pt + arclen_b_pt --- 440,446 ---- for param_a, param_b, length_pt in zip(params_a, params_b, lengths_pt): if length_pt > arclen_a_pt: ! params.append(b.subparamtoparam(param_b)) else: ! params.append(a.subparamtoparam(param_a)) return params, arclen_a_pt + arclen_b_pt *************** *** 541,546 **** # To improve the performance in the general case we alternate the # splitting process between the two normsubpathitems ! 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): --- 541,546 ---- # To improve the performance in the general case we alternate the # splitting process between the two normsubpathitems ! return ( [(a.subparamtoparam(a_t), o_t) for o_t, a_t in other.intersect(a, epsilon)] + ! [(b.subparamtoparam(b_t), o_t) for o_t, b_t in other.intersect(b, epsilon)] ) def modifiedbegin_pt(self, x_pt, y_pt): *************** *** 705,714 **** # 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" --- 705,714 ---- # curve replacements used by midpointsplit: # The replacements are normline_pt and normcurve_pt instances with an ! # additional subparamtoparam function for proper conversion of the ! # parametrization. Note that we only one direction (when a parameter ! # gets calculated), since the other way around direction midpointsplit ! # is not needed at all ! class _leftnormline_pt(normline_pt): __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" *************** *** 720,730 **** 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] --- 720,733 ---- self.l3_pt = l3_pt ! def subparamtoparam(self, param): if 0 <= param <= 1: ! params = mathutils.realpolyroots(self.l1_pt-2*self.l2_pt+self.l3_pt, ! -3*self.l1_pt+3*self.l2_pt, ! 3*self.l1_pt, ! -param*(self.l1_pt+self.l2_pt+self.l3_pt)) # we might get several solutions and choose the one closest to 0.5 + # (we want the solution to be in the range 0 <= param <= 1; in case + # we get several solutions in this range, they all will be close to + # each other since l1_pt+l2_pt+l3_pt-l0_pt < epsilon) params.sort(lambda t1, t2: cmp(abs(t1-0.5), abs(t2-0.5))) return 0.5*params[0] *************** *** 736,760 **** ! 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 --- 739,763 ---- ! class _rightnormline_pt(_leftnormline_pt): __slots__ = "x0_pt", "y0_pt", "x1_pt", "y1_pt", "l1_pt", "l2_pt", "l3_pt" ! def subparamtoparam(self, param): ! return 0.5+_leftnormline_pt.subparamtoparam(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 subparamtoparam(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 subparamtoparam(self, param): return 0.5+0.5*param Index: unit.py =================================================================== RCS file: /cvsroot/pyx/pyx/pyx/unit.py,v retrieving revision 1.53 retrieving revision 1.54 diff -C2 -d -r1.53 -r1.54 *** unit.py 8 Jul 2005 13:03:48 -0000 1.53 --- unit.py 10 Feb 2006 16:48:21 -0000 1.54 *************** *** 23,27 **** import types - import helper scale = { 't':1, 'u':1, 'v':1, 'w':1, 'x':1 } --- 23,26 ---- --- NEW FILE: mathutils.py --- #!/usr/bin/env python # -*- coding: ISO-8859-1 -*- # # # Copyright (C) 2005-2006 Michael Schindler # Copyright (C) 2006 André Wobst # # 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 try: import Numeric, LinearAlgebra _has_numeric = 1 except: _has_numeric = 0 def sign(x): """sign of x, i.e. +1 or -1; returns 1 for x == 0""" if x >= 0: return 1 return -1 def asinh(x): """Return the arc hyperbolic sine of x.""" return math.log(x+math.sqrt(x*x+1)) def acosh(x): """Return the arc hyperbolic cosine of x.""" return math.log(x+math.sqrt(x*x-1)) def _realroots_quadratic(a1, a0): """gives the real roots of x**2 + a1 * x + a0 = 0""" D = a1*a1 - 4*a0 if D < 0: return [] SD = math.sqrt(D) return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)] def _realroots_cubic(a2, a1, a0): """gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0""" # see http://mathworld.wolfram.com/CubicFormula.html for details Q = (3*a1 - a2*a2) / 9.0 R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0 D = Q*Q*Q + R*R if D > 0: # one real and two complex roots SD = math.sqrt(D) if R + SD >= 0: S = (R + SD)**(1/3.0) else: S = -(-R - SD)**(1/3.0) if R - SD >= 0: T = (R - SD)**(1/3.0) else: T = -(SD - R)**(1/3.0) return [S + T - a2/3.0] elif D == 0: if Q == 0: # one real root (R==0) return [-a2/3.0] else: # two real roots (R>0, Q<0) S = -math.sqrt(-Q) return [2*S - a2/3.0, -S - a2/3.0] else: # three real roots (Q<0) SQ = math.sqrt(-Q) theta = math.acos(R/(SQ**3)) return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)] def _realroots_quartic(a3, a2, a1, a0): """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0""" # see http://mathworld.wolfram.com/QuarticEquation.html for details ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3) ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and y*y-4*a0 >= 0] if not ys: return [] y1 = min(ys) if a3*y1-2*a1 < 0: return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) + _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0)))) else: return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) + _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0)))) def realpolyroots(*cs): """returns the roots of a polynom with given coefficients polynomial with coefficients given in cs: 0 = \sum_i cs[i] * x^(len(cs)-i-1) """ if not cs: return [0] try: f = 1.0/cs[0] cs = [f*coeff for coeff in cs[1:]] except ArithmeticError: return realpolyroots(*cs[1:]) else: n = len(cs) if n == 0: return [] elif n == 1: return [-cs[0]] elif n == 2: return _realroots_quadratic(*cs) elif n == 3: return _realroots_cubic(*cs) elif n == 4: return _realroots_quartic(*cs) else: raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4") def realpolyroots_eigenvalue(*cs): # as realpolyroots but using an equivalent eigenvalue problem # (this code is currently used for functional tests only) if not _has_numeric: raise RuntimeError("realpolyroots_eigenvalue depends on Numeric") if not cs: return [0] try: f = 1.0/cs[0] cs = [f*coeff for coeff in cs[1:]] except ArithmeticError: return realpolyroots_eigenvalue(*cs[1:]) else: if not cs: return [] n = len(cs) 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] = -cs[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 --- helper.py DELETED ---