#!/usr/bin/python2.6 #author: Bryan Bishop #date: 2011-04-26 from copy import copy, deepcopy from fractions import Fraction, gcd import inspect #for lineno() init_tol = [1, 1024] ###################### # NOTES # # A number of functions "return num_*". In many of these cases, they should # return whatever list they are manipulating, since we can find lengths easily. # ###################### bryan_message = "bryan hasn't implemented this yet" def lineno(): """returns the current line number in our program. author: Danny Yoo """ return str(inspect.currentframe().f_back.f_lineno) def sgn(number): if number > 0: return 1 elif number == 0: return 0 elif number < 0: return -1 raise Exception, "this should never happen" def den(thing): #return mpq_denref(thing) if type(thing) == float: return Fraction.from_float(thing).numerator return Fraction(str(thing)).denominator def num(thing): #return mpq_numref(thing) if type(thing) == float: return Fraction.from_float(thing).numerator return Fraction(str(thing)).numerator class K_SEGMENT: def __init__(self, *args): if len(args) == 2: first = args[0] second = args[1] if not isinstance(first, Point2D): first = Point2D(first) if not isinstance(second, Point2D): second = Point2D(second) self.start = first self.end = second def __repr__(self): return "K_SEGMENT(start=" + str(self.start) + ",end=" + str(self.end) + ")" def reverse(self): s = K_SEGMENT() s.start = self.end s.end = self.start return s def outer_box(self): b = self.start.bbox().merge(self.end.bbox()) i = 0 while i < 2: if b.low[i] == b.high[i]: b.low_open[i] = b.high_open[i] = 0 else: b.low_open[i] = b.high_open[i] = 1 i += 1 return b def inner_box(self): """ returns a pointer to the inner box for *this. returns False if the inner box for *this is not well-defined. """ b_proto = BOXCO2() s = self.start.bbox() e = self.end.bbox() well_defed = True i = 0 while well_defed and i < 2: if s.low[i] > e.high[i]: b_proto.low[i] = e.high[i] b_proto.high[i] = s.low[i] b_proto.low_open[i] = 1 - e.high_open[i] b_proto.high_open[i] = 1 - s.low_open[i] i += 1 elif s.high[i] < e.low[i]: b_proto.low[i] = s.high[i] b_proto.high[i] = e.low[i] b_proto.low_open[i] = 1 - s.high_open[i] b_proto.high_open[i] = 1 - e.low_open[i] i += 1 elif s.low[i] == e.high[i] and s.low_open[i] and e.high_open[i]: b_proto.low[i] = b_proto.high[i] = s.low[i] b_proto.low_open[i] = b_proto.high_open[i] = 0 i += 1 elif s.high[i] == e.low[i] and s.high_open[i] and e.low_open[i]: b_proto.low[i] = b_proto.high[i] = s.high[i] b_proto.low_open[i] = b_proto.high_open[i] = 0 i += 1 else: well_defed = False if well_defed: b = BOXCO2(b_proto) else: b = False return b class Curve: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 3: if isinstance(args[0], K_RATPOLY): self.init_ratpoly(args) elif isinstance(args[0], K_CURVE): self.init_curve(args) else: raise Exception, "argument not understood" else: raise Exception, "too many parameters" def init_noargs(self): self.poly = K_RATPOLY() self.segments = [] self.num_segments = 0 self.curve_in_other_dom = 0 self.dir_in_other_dom = 0 def init_ratpoly(self, args): #K_RATPOLY P, K_SEGMENT s[], n P = args[0] s = args[1] n = args[2] self.poly = K_RATPOLY(P) self.num_segments = n if self.num_segments > 0: self.segments = copy(s) else: self.segments = [] self.curve_in_other_dom = 0 self.dir_in_other_dom = 0 def assoc(self, curve, dir): assert curve self.curve_in_other_dom = curve curve.curve_in_other_dom = self self.dir_in_other_dom = curve.dir_in_other_dom = dir def bbox(self): self.num_segments = len(self.segments) assert self.num_segments > 0 b = self.segments[0].outer_box() i = 1 while i < self.num_segments: b = b.merge(self.segments[i].outer_box()) i += 1 return b def get_pts_proto_2(b_s, l_t, h_t, t_only, P1, P2, tol, count_edges): """ unsigned long get_pts_proto(const bigrational& b_s, const bigrational& l_t, const bigrational& h_t, const K_RATPOLY& t_only, const K_RATPOLY& P1, const K_RATPOLY& P2, K_POINT2D**& pts, const bigrational& tol, const int count_edges) computes the intersections "pts" of 2 bivariate polynomials P1 and P2 on/in the vertical line segment [(b_s, l_t), (b_s, h_t)]. by computing the roots of a univariate polynomial t_only on/in the interval [l_t, h_t]. returns the list of intersections if tol > 0 then a box for each intersection is no larger than tol in the t-direction. if tol = 0 then there is no limit on the size of boxes for intersections. if count_edges = 1 then the intersections at the endpoints of the line segment are counted. if count_edges = 0 then the intersections at the endpoints of the line segment are ignored. """ #kpoint2d.cc line 3515, 3684 print "Computing the intersection points of 2 bivariate polynomials:" print "\tP1 = ", str(P1) print "\tP2 = ", str(P2) print "on/in the vertical line segment [(" + str(b_s) + ", " + str(l_t) + "), (" + str(b_s) + ", " + str(h_t) + ")]." raise NotImplementedError, bryan_message def gen_curve_topo_proto(P, l_s, h_s, l_t, h_t, edge_l_s, edge_h_s, edge_l_t, edge_h_t, turn_s, turn_t): #kcurve.cc line 1074 to 3535 assert P.num_vars == 2 curves = [] raise NotImplementedError, bryan_message return curves def gen_curve_topo(P, l_s, h_s, l_t, h_t): """ resolves topology of a bivariate polynomial P in a region [l_s, h_s] x [l_t, h_t] to obtain monotone curves curves[0], curves[1], ..., curves[num_curves - 1]. returns curves. """ curves = [] assert P.num_vars == 2 assert l_s < h_s assert l_t < h_t P_l_s, P_h_s, P_l_t, P_h_t, L_s, H_s, L_t, H_t = K_RATPOLY(),K_RATPOLY(),K_RATPOLY(),K_RATPOLY(),K_RATPOLY(),K_RATPOLY(),K_RATPOLY(),K_RATPOLY() edge_l_s, edge_h_s, edge_l_t, edge_h_t = [], [], [], [] #or Point2D?? dP_ds, dP_dt = K_RATPOLY(), K_RATPOLY() turn_s, turn_t = [], [] #or Point2D?? print "kcurve: gen_curve_topo: -----------------------------" print "\tP = " + str(P) print "\t[" + str(l_s) + ", " + str(h_s) + "] x [" + str(l_t) + ", " + str(h_t) + "]" print "-----------------------------------------------------" P_l_s = P.subst_val(0, l_s) P_h_s = P.subst_val(0, h_s) P_l_t = P.subst_val(1, l_t) P_h_t = P.subst_val(1, h_t) L_s = K_RATPOLY(1, 0, l_s) H_s = K_RATPOLY(1, 0, h_s) L_t = K_RATPOLY(1, 0, l_t) H_t = K_RATPOLY(1, 0, h_t) print "\tP_l_s = ", P_l_s print "\tP_h_s = ", P_h_s print "\tP_l_t = ", P_l_t print "\tP_h_t = ", P_h_t init_tol = "unused" edge_l_s = get_pts_proto_2(l_s, l_t, h_t, P_l_s, P, L_s, init_tol, 1) edge_h_s = get_pts_proto_2(h_s, l_t, h_t, P_h_s, P, H_s, init_tol, 1) edge_l_t = get_pts_proto_2(l_s, h_s, P_l_t, l_t, P, L_t, init_tol, 0) edge_h_t = get_pts_proto_2(l_s, h_s, P_h_t, h_t, P, H_t, init_tol, 0) #lengths num_edge_l_s = len(edge_l_s) num_edge_h_s = len(edge_h_s) num_edge_l_t = len(edge_l_t) num_edge_h_t = len(edge_h_t) #kcurve.cc line 966 if num_edge_l_s > 1: edge_l_s = sort_t(edge_l_s) if num_edge_h_s > 1: edge_h_s = sort_t(edge_h_s) if num_edge_l_t > 1: edge_l_t = sort_s(edge_l_t) if num_edge_h_t > 1: edge_h_t = sort_s(edge_h_t) print "edge_l_s = ", edge_l_s print "edge_h_s = ", edge_h_s print "edge_l_t = ", edge_l_t print "edge_h_t = ", edge_h_t print " --------------------------------------------" dP_ds = P.derivative(0) dP_dt = P.derivative(1) print "\tdP_ds = ", dP_ds print "\tdP_dt = ", dP_dt print " --------------------------------------------" turn_s = get_pts(l_s, h_s, l_t, h_t, P, dP_ds, init_tol, 0) turn_t = get_pts(l_s, h_s, l_t, h_t, P, dP_dt, init_tol, 0) num_turn_s = len(turn_s) num_turn_t = len(turn_t) print "turn_s = ", turn_s print "turn_t = ", turn_t print " --------------------------------------------" curves = gen_curve_topo_proto(P, l_s, h_s, l_t, h_t, edge_l_s, edge_h_s, edge_l_t, edge_h_t, turn_s, turn_t) return curves class K_RATPOLY: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 1: self.init_args_from_ratpoly(args[0]) elif len(args) == 2: self.init_args_nv_d(args[0], args[1]) elif len(args) == 3: self.init_args_nvib(args[0], args[1], args[2]) elif len(args) == 4: self.init_args_nv_nc(args[0], args[1], args[2], args[3]) else: raise Exception, "wrong number of arguments" def init_noargs(self): self.num_vars = 0 self.degree = [] self.num_coeffs = 1 self.coeffs = [0] self.Sturm_seq = 0 def init_args_nv_d(self, nv, d): """constructs a polynomial of nv variables of degree d[0], d[1], ..., d[nv - 1] with anonymous coefficients.""" self.num_vars = nv if nv > 0: self.degree = [None] * nv else: self.degree = [] self.num_coeffs = 1 i = 0 while i < nv: self.degree[i] = d[i] self.num_coeffs = self.num_coeffs * (d[i] + 1) i += 1 #for i in range(0, nv): # self.degree[i] = d[i] # self.num_coeffs *= d[i] + 1 self.coeffs = [0] * self.num_coeffs self.Sturm_seq = 0 def init_args_nv_nc(self, nv, d, nc, c): """constructs a polynomial of nv variables of degree d[0], d[1], ..., d[nv - 1] with coefficients c[0], c[1], ..., c[nc - 1]""" self.num_vars = nv if nv > 0: self.degree = [None] * nv else: self.degree = [] self.num_coeffs = 1 for i in range(0, nv): self.degree[i] = d[i] self.num_coeffs *= d[i] + 1 assert nc == self.num_coeffs self.coeffs = [] * self.num_coeffs for i in range(0, self.num_coeffs): self.coeffs[i] = c[i] self.Sturm_seq = 0 self.reduce_deg() def init_args_nvib(self, nv, i, b): """constructs a polynomial X_i - b of nv variables""" assert i < nv # i < nv => nv > 0 self.num_vars = nv self.degree = [None] * self.num_vars for j in range(0, i): self.degree[j] = 0 self.degree[i] = 1 for j in range(i+1, self.num_vars): self.degree[j] = 0 self.num_coeffs = 2 self.coeffs = [1, -b] self.Sturm_seq = 0 def init_args_from_ratpoly(self, p): """the copy constructor please just use copy() or deepcopy()...""" if not isinstance(p, K_RATPOLY): raise Exception, "argument must be of type K_RATPOLY" self.num_vars = p.num_vars self.degree = copy(p.degree) self.num_coeffs = p.num_coeffs self.coeffs = copy(p.coeffs) self.Sturm_seq = p.Sturm_seq def __repr__(self): return "K_RATPOLY(coeffs=%s, degree=%s)" % (str(self.coeffs), str(self.degree)) def get_num_vars(self): return len(self.degree) def index_to_powers(self, index): """returns p = (p[0], p1[1], ..., p[num_vars - 1] s.t. coeffs[i] is the coefficient of the monomial X^p of *this.""" assert index < len(self.coeffs) if len(self.degree) > 0: j = index p = [None] * len(self.degree) k = len(self.degree) - 1 while k >= 0: p[k] = self.degree[k] - j % (self.degree[k] + 1) j /= self.degree[k] + 1 k -= 1 else: p = [] return p def index_to_deg(self, i): """returns s.t. coeffs[i] is the coefficient of the monomial of *this of degree t.""" j = i t = 0 k = len(self.coeffs) - 1 while k >= 0: t += self.degree[k] - j % (self.degree[k] + 1) j /= self.degree[k] + 1 k -= 1 return t def get_total_deg(self): """returns the total degree""" t = 0 i = 0 while i < len(self.coeffs): ti = self.index_to_deg(i) if sgn(self.coeffs[i]) and ti > t: t = ti return t def powers_to_index(self, p): """returns i s.t. coeffs[i] is the coefficient of teh monomial X^p of *this.""" i = 0 if p: j = 0 k = len(self.coeffs) while j < len(self.degree): assert p[j] <= self.degree[j] k /= self.degree[j] + 1 i += k * (self.degree[j] - p[j]) j += 1 assert i < len(self.coeffs) return i def get_coeff(self, p): """returns the coefficient of the monomial X^p of *this.""" return self.coeffs[self.powers_to_index(p)] def reduce_deg(self): """reduces deg[0], deg[1], ..., deg[num_vars - 1] to max. degrees in X_0, X_1, ..., X_{num_vars - 1}, resp. returns True if some deg[i] is reduced, and False otherwise.""" #locals: i, j, k, l, d, p, nc, c reduced = False if self.num_vars > 0: d = [0] * self.num_vars i = 0 while i < self.num_coeffs: p = self.index_to_powers(i) if sgn(self.coeffs[i]) and p: j = 0 while j < self.num_vars: if d[j] < p[j]: d[j] = p[j] j+= 1 i += 1 i = 0 while not reduced and i < self.num_vars: if d[i] < self.degree[i]: reduced = True else: #if d[i] >= self.deg[i] i += 1 if reduced: i = 0 nc = 1 while i < self.num_vars: nc *= d[i] + 1 i+= 1 c = [0] * nc i = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): p = self.index_to_powers(i) j = k = 0 l = nc while k < self.num_vars: assert p[k] <= d[k] l /= d[k] + 1 j += l * (d[k] - p[k]) k += 1 c[j] = self.coeffs[i] i+= 1 self.num_coeffs = nc self.coeffs = c self.deg = d return reduced def add_var(self, i): """returns a polynomial of num_vars + 1 variables of degree deg[0], deg[1], ..., deg[i - 1], 0, deg[i], ..., deg[num_vars - 1].""" assert i < (self.num_vars + 1) #locals: j, dp dp = [0] * (self.num_vars + 1) j = 0 while j < i: dp[j] = self.degree[j] j += 1 dp[i] = 0 j = i while j < self.num_vars: dp[j + 1] = self.degree[j] j += 1 P = K_RATPOLY(self.num_vars + 1, dp) j = 0 while j < self.num_coeffs: P.coeffs[j] = self.coeffs[j] j += 1 return P def remove_var(self, i): """PROVIDED self.degree[i] == 0, returns a polynomial of num_vars - 1 variables of degree deg[0], deg[1], ..., deg[i - 1], deg[i + 1], ..., deg[num_vars - 1].""" assert i < self.num_vars assert self.degree[i] == 0 raise NotImplementedError, bryan_message def add(self, P): """returns *this + P.""" if isinstance(P, K_RATPOLY): return self.add_ratpoly(P) raise NotImplementedError, bryan_message def add_ratpoly(self, P): """returns *this + P.""" assert self.num_vars == P.num_vars if self.num_vars > 0: dq = [] i = 0 while i < self.num_vars: #TODO: confirm that this is the correct form of dq[i] = deg[i] > P.deg[i] ? deg[i] : P.deg[i]; if self.degree[i] > P.degree[i]: dq.append(self.degree[i]) else: dq.append(P.degree[i]) i += 1 else: dq = [] Q = K_RATPOLY(self.num_vars, dq) i = 0 while i < self.num_coeffs: pq = self.index_to_powers(i) j = Q.powers_to_index(pq) Q.coeffs[j] = self.coeffs[i] i += 1 i = 0 while i < P.num_coeffs: pq = P.index_to_powers(i) j = Q.powers_to_index(pq) Q.coeffs[j] += P.coeffs[i] i += 1 Q.reduce_deg() return Q def sub(self, P): """returns *this - P.""" raise NotImplementedError, bryan_message def mul(self, P): """returns *this * P.""" assert self.num_vars == P.num_vars if self.num_vars > 0: dq = [] i = 0 while i < self.num_vars: dq.append(self.degree[i] + P.degree[i]) i += 1 else: #if self.num_vars == 0 dq = [] Q = K_RATPOLY(self.num_vars, dq) i = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): p = self.index_to_powers(i) j = 0 while j < P.num_coeffs: if sgn(P.coeffs[j]): pp = P.index_to_powers(j) if self.num_vars > 0: pq = [] k = 0 while k < self.num_vars: pq.append(p[k] + pp[k]) k += 1 else: pq = [] k = Q.powers_to_index(pq) Q.coeffs[k] += self.coeffs[i] * P.coeffs[j] j += 1 i += 1 Q.reduce_deg() return Q def neg(self): """returns - *this.""" P = K_RATPOLY(self.num_vars, self.degree) for i in range(0, self.num_coeffs): P.coeffs[i] = - self.coeffs[i] return P def exact_div(self, P): """returns *this / P PROVIDED *this is divisible by P.""" raise NotImplementedError, bryan_message def __add__(self, b): return self.add(b) def __sub__(self, b): return self.sub(b) def __mul__(self, b): return self.mul(b) def __neg__(self, b): return self.neg() def __div__(self, b): return self.div(b) def __eq__(self, b): #see cmp(const K_RATPOLY& P) raise NotImplementedError, bryan_message def derivative(self, i): """returns the partial derivative of *this w.r.t. X_i""" raise NotImplementedError, bryan_message def is_zero(self): """returns True if this is a zero polynomial, False otherwise.""" i = 0 while x < self.num_coeffs and not sgn(self.coeffs[i]): i += 1 if i < self.num_coeffs: return False else: return True #if i == num_coeffs def eq_upto_const(self, P): """returns True if this == P*b for some bigrational b, and False otherwise.""" raise NotImplementedError, bryan_message def subst_val_first_var_proto(self, b, reduced): """returns a polynomial of num_vars - 1 variables obtained by computing *this(b, X_1, ..., X_{num_vars - 1}), and renaming X_1, ..., X_{num_vars - 1} to X_0, X_1, ..., X_{num_vars - 2}. reduce_deg() is applied if reduced == 1.""" if self.num_vars > 0 and self.num_coeffs > 1: #FIXME: was just if self.num_vars > 0 if self.num_vars > 1: dp = [None] * (self.num_vars - 1) i = 1 while i < self.num_vars: dp[i - 1] = self.degree[i] i += 1 else: #if self.num_vars == 1 dp = [] #0? P = K_RATPOLY(self.num_vars - 1, dp) print "P.coeffs = ", P.coeffs print "self.coeffs = ", self.coeffs #should not be [0] P.num_coeffs = len(P.coeffs) i = 0 while i < P.num_coeffs: print "hi" P.coeffs[i] = self.coeffs[i] i += 1 i = 1 while i < self.degree[0]: j = 0 while j < P.num_coeffs: P.coeffs[j] *= b P.coeffs[j] += self.coeffs[j + i * P.num_coeffs] j += 1 i += 1 if reduced: P.reduce_deg() else: #if self.num_vars == 0 P = self return P def subst_val_first_var(self, b): """returns a polynomial of num_vars - 1 variables obtained by computing *this(b, X_1, ..., X_{num_vars - 1}), and renaming X_1, ..., X_{num_vars - 1} to X_0, X_1, ..., X_{num_vars - 2}.""" return self.subst_val_first_var_proto(b, 1) #There are multiple "evaluate" definitions. Make the one you need. #def evaluate(self, b): # """returns *this(b) PROVIDED *this is a univariate polynomial.""" # assert self.num_vars == 1 # return self.subst_val_first_var(b).coeffs[0] def fp_sgn_at(self, b): """PROVIDED *this is a univariate polynomial, returns 1 if *this(b) is certainly positive, - 1 if *this(b) is certainly negative, and 0 otherwise.""" raise NotImplementedError, bryan_message def sgn_at(self, b): """PROVIDED *this is a univariate polynomial, returns the sign of *this(b).""" assert self.num_vars == 1 if isinstance(b, list): if self.num_vars == 0: return self.sgn_at(b[0]) else: return sgn(self.evaluate(b)) else: fs = self.fp_sgn_at(b) if fs: s = fs else: s = sgn(self.evaluate(b)) return s def subst_val_proto(self, i, b, reduced): """returns a polynomial of num_vars - 1 variables obtained by computing *this(X_0, X_1, ..., X_{i - 1}, b, X_{i + 1}, ..., X_{num_vars - 1}) and, renaming X_{i + 1}, ..., X_{num_vars - 1} to X_i, ..., X_{num_vars - 2}. reduce_deg() is applied if reduced == 1.""" #kratpoly.cc line 1321 assert self.num_vars == 0 or i < self.num_vars if self.num_vars > 0: if i == 0: P = self.subst_val_first_var_proto(b, reduced) else: #if i > 0 if self.num_vars > 1: dp = [None] * (self.num_vars - 1) j = 0 while j < i: dp[j] = self.degree[j] j += 1 j = i + 1 while j < self.num_vars: dp[j - 1] = self.degree[j] j += 1 else: dp = [] #if num_vars == 1 P = K_RATPOLY(self.num_vars - 1, dp) bp = [None] * (self.degree[i] + 1) bp[0] = 1 if self.degree[1] > 0: bp[1] = b j = 1 while j <= self.degree[i]: bp[j] = b * bp[j - 1] j += 1 j = 0 while j < self.num_coeffs: if sgn(self.coeffs[j]): p = self.index_to_powers(j) if self.num_vars > 1: pp = [None] * (self.num_vars - 1) k = 0 while k < i: pp[k] = p[k] k += 1 k = i + 1 while k < self.num_vars: pp[k - 1] = p[k] k += 1 else: #if num_vars == 1 pp = [] #FIXME: 0? k = P.powers_to_index(pp) P.coeffs[k] += self.coeffs[j] * bp[p[i]] j += 1 if reduced: P.reduce_deg() else: P = self return P def subst_val(self, i, b): """returns a polynomial of num_vars - 1 variables obtained by computing *this(X_0, X_1, ..., X_{i - 1}, b, X_{i + 1}, ..., X_{num_vars - 1}) and, renaming X_{i + 1}, ..., X_{num_vars - 1} to X_i, ..., X_{num_vars - 2}.""" return self.subst_val_proto(i, b, 1) def subst_expr(self, i, P): """returns a polynomial of num_vars - 1 variables obtained by computing *this(X_o, X_1, ..., X_{i - 1}, P (X_0, X_1, ..., X_{i - 1}, X_{i + 1}, X_{num_vars - 1}), X_{i+1}, ..., X_{num_vars - 1}) and, renaming X_{i + 1}, ..., X_{num_vars - 1} to X_i, ..., X_{num_vars - 2}. """ assert i < self.num_vars assert self.num_vars == P.num_vars + 1 #locals: j, k, l, dq, PP, dpp0, p, p1, ppp, pq, Q if self.num_vars > 1: dq = [None] * (self.num_vars - 1) j = 0 while j < i: dq[j] = self.degree[j] + self.degree[i] * P.degree[j] j += 1 j = i + 1 while j < self.num_vars: dq[j - 1] = self.degree[j] + self.degree[i] * P.degree[j - 1] j += 1 else: dq = [] Q = K_RATPOLY(self.num_vars - 1, dq) PP = [K_RATPOLY()] * (self.degree[i] + 1) if P.num_vars > 0: dpp0 = [0] * P.num_vars else: dpp0 = [] PP[0] = K_RATPOLY(P.num_vars, dpp0) PP[0].coeffs[0] = 1 if self.degree[i] > 0: PP[1] = P j = 2 while j <= self.degree[i]: PP[j] = P * PP[j - 1] j += 1 j = 0 while j < self.num_coeffs: if sgn(self.coeffs[j]): p = self.index_to_powers(j) if self.num_vars > 1: p1 = [None] * (self.num_vars - 1) k = 0 while k < i: p1[k] = p[k] k += 1 k = i + 1 while k < self.num_vars: p1[k - 1] = p[k] k += 1 else: p1 = [] k = 0 while k < PP[p[i]].num_coeffs: if sgn(PP[p[i]].coeffs[k]): ppp = PP[p[i]].index_to_powers(k) if self.num_vars > 1: pq = [None] * (self.num_vars - 1) l = 0 while l < self.num_vars -1: pq[l] = p1[l] + ppp[l] l += 1 else: pq = [] l = Q.powers_to_index(pq) Q.coeffs[l] += self.coeffs[j] * PP[p[i]].coeffs[k] k += 1 j += 1 Q.reduce_deg() return Q #subst_expr(i, P, D) def subst_param_expr(self, X, Y, Z, W): """PROVIDED *this is a 3-variate polynomial and X, Y, Z, W are bivariate polynomials, returns *this(X/W, Y/W, Z/W).""" #see kratpoly.cc line 1744 assert self.num_vars == 3 assert X.num_vars == 2 assert Y.num_vars == 2 assert Z.num_vars == 2 assert W.num_vars == 2 d = [0, 0] #FIXME: this seems very, very broken I = K_RATPOLY(2, d) #<--- original, but results in I.coeffs = [0] ##I = K_RATPOLY(self.num_vars, self.degree) #I = K_RATPOLY(self) print "I = ", I i = 0 td = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): p = self.index_to_powers(i) t = p[0] + p[1] + p[2] if t > td: td = t i += 1 i = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): S = K_RATPOLY(2, d) S.coeffs[0] = self.coeffs[i] p = self.index_to_powers(i) dd = td - p[0] - p[1] - p[2] j = 0 while j < p[0]: S = S * X #FIXME: make sure K_RATPOLY::__mul__ is implemented j += 1 j = 0 while j < p[1]: S = S * Y j += 1 j = 0 while j < p[2]: S = S * Z j += 1 j = 0 while j < dd: S = S * W j += 1 I = I + S i += 1 I.reduce_deg() I.reduce_num_coeffs() return I #as_FLOATPOLY() #div(P, R) #div(P1, P2, R) def rem(self, P2): return self.div(P2) def set_Sturm_seq(self): """PROVIDED *this is a univariate polynomial, generates Sturm sequence of *this.""" raise NotImplementedError, bryan_message def num_Sturm_seq_perm(self, b): """PROVIDED *this is a univariate polynomial, returns the number of sign permanencies in the Sturm sequence of *this at b. """ raise NotImplementedError, bryan_message def gcd1(self, P): """PROVIDED *this and P are univariate polynomials, return their gcd computed by using Euclidean algorithm.""" raise NotImplementedError, bryan_message def Sylvester1(self, P): """PROVIDED *this and P are univariate polynomials, returns Sylvester resultant for *this and P.""" raise NotImplementedError, bryan_message def Sylvester2(self, P, i): """PROVIDED *this and P are bivariate polynomials, returns Sylvester resultant for *this and P w.r.t. X_i.""" raise NotImplementedError, bryan_message def GoodSylvester(self, P, i): """PROVIDED *this and P are bivariate polynomials, returns Sylvetser resultant for *this and P w.r.t. X_i.""" raise NotImplementedError, bryan_message def monic(self): """PROVIDED *this is a univariate polynomial, returns a monic polynomial obtained by dividing all the coeffs of *this by the leading coeff *this.""" assert self.num_vars == 1 P = K_RATPOLY(1, deg) lc = self.coeffs[0] if sgn(lc): for i in range(0, self.num_coeffs): if sgn(self.coeffs[i]): P.coeffs[i] = self.coeffs[i] / lc else: P.coeffs[0] = 0 #if this is a zero polynomial return P def monic_gcd1(self, P): """PROVIDED *this and P are univariate polynomials, returns monic gcd of *this and P.""" raise NotImplementedError, bryan_message def gcd2_pp(self, P): """PROVIDED *this and P are bivariate primitive polynomials, returns gcd(*this, P).""" raise NotImplementedError, bryan_message def gcd2(self, P): """PROVIDED *this and P are bivariate polynomials, returns gcd(*this, P).""" raise NotImplementedError, bryan_message def eval_range(self, low_in, high_in): """computes the range [low_out, high_out] which *this is evaluated at some value in [low_in, high_in] using affine arithmetic.""" i, j = 0, 0 d = [1] T = K_RATPOLY(self) err = 0 i = 0 while i < self.num_vars: N = K_RATPOLY(1, d) N.coeffs[0] = (high_in[i] - low_in[i]) / 2 N.coeffs[1] = (low_in[i] + high_in[i]) / 2 T = T.add_var(i) j = 0 while j < i: N = N.add_var(j) j += 1 j = i + 1 while j < self.num_vars: N = N.add_var(j) j += 1 T = T.subst_expr(i + 1, N) i += 1 err = 0 i = 0 while i < T.num_coeffs - 1: err += abs(T.coeffs[i]) i += 1 low_out = T.coeffs[T.num_coeffs - 1] - err high_out = T.coeffs[T.num_coeffs - 1] + err return low_out, high_out def eval_range_in_out(self, low_in, high_in, low_out, high_out): """computes the range [low_out, high_out] which *this is evaluated at some value in [low_in, high_in] using affine arithmetic.""" assert self.num_vars == len(low_in) assert self.num_vars == len(high_in) #local variables: i, j, err d = [1] N = K_RATPOLY() T = copy(self) i = 0 while i < self.num_vars: N = K_RATPOLY(1, d) N.coeffs[0] = (high_in[i] - low_in[i]) / 2 N.coeffs[1] = (low_in[i] + high_in[i]) / 2 T = T.add_var(i) j = 0 while j < i: N = N.add_var(j) j += 1 j = i + 1 while j < self.num_vars: N = N.add_var(j) j += 1 T = T.subst_expr(i + 1, N) i += 1 err = 0 i = 0 while i < T.num_coeffs - 1: err += abs(T.coeffs[i]) i += 1 low_out = T.coeffs[T.num_coeffs - 1] - err high_out = T.coeffs[T.num_coeffs - 1] + err return (low_out, high_out) def get_Mignotte_bd(self): """PROVIDED *this is a univariate polynomial, returns an upper bound on the size of the largest root of *this.""" raise NotImplementedError, bryan_message def reduce_num_coeffs(self): i, j, g = None, None, None d = 1 i = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): d *= den(self.coeffs[i]) #d *= Fraction(self.coeffs[i]).limit_denominator() #d *= den(self.coeffs[i]) i += 1 i = 0 while i < self.num_coeffs: if sgn(self.coeffs[i]): self.coeffs[i] *= d i += 1 i = 0 while i < self.num_coeffs and not sgn(self.coeffs[i]): i+= 1 g = 1 if i < self.num_coeffs: #g = Fraction(self.coeffs[i]).numerator g = num(self.coeffs[i]) j = i + 1 while j < self.num_coeffs and g != 1: if sgn(self.coeffs[j]): g = gcd(g, num(self.coeffs[j])) #might be self.gcd or self.num j += 1 if g != 1: j = 0 while j < self.num_coeffs: if sgn(self.coeffs[j]): self.coeffs[j] /= g j += 1 return d != 1 def reduce_coeffs(self): raise NotImplementedError, bryan_message def conv_to_Bernstein(self, deg_s, deg_t): raise NotImplementedError, bryan_message def transform_Impl(self, T): raise NotImplementedError, bryan_message def x_gcd(P1, P2): """PROVIDED P1 and P2 are uni or bivariate polynomials, returns their gcd.""" raise NotImplementedError, bryan_message def Sylvester1(P1, P2): return P1.Sylvester1(P2) def gen_VDM_row1(n): """returns the first row of Vandermonde matrix, i.e., a vector of n distinct rational numbers.""" R = [None] * n for i in range(0, n / 2): R[2 * i] = i + 1 R[2 * i + 1] = - i - 1 if n % 2: R[n - 1] = n / 2 + 1 return R def GoodSylvester(P1, P2, i): return P1.GoodSylvester(P2, i) class K_BOX3D: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 6: self.init_lows_and_highs(args) else: raise Exception, "too many arguments" def init_noargs(self): i = 3 self.low, self.high, self.low_infty, self.high_infty = [0] * i, [0] * i, [-1] * 3, [1] * 3 def init_lows_and_highs(self, args): l_s, h_s, l_t, h_t, l_u, h_u = args[0], args[1], args[2], args[3], args[4], args[5] self.low = [l_s, l_t, l_u] self.high = [h_s, h_t, h_u] self.low_infty = self.high_infty = [0] * 3 def overlap(self, b): o = 1 i = 0 while o and i < 3: if self.low_infty[i] < 0: if self.high_infty[i] < 0 and b.low_infty[i] >= 0: o = 0 elif not self.high_infty[i]: if b.low_infty[i] > 0: o = 0 elif not b.low_infty[i] and self.high[i] < b.low[i]: o = 0 elif not self.low_infty[i]: if b.high_infty[i] < 0: o = 0 elif not b.high_infty[i] and self.low[i] > b.high[i]: o = 0 elif not self.high_infty[i]: if b.low_infty[i] > 0: o = 0 elif not b.low_infty[i] and self.high[i] < b.low[i]: o = 0 else: #if self.low_infty[i] > 0 if b.high_infty[i] <= 0: o = 0 i += 1 return o #BOXCO2(self.PtBs, self.PtBs, self.PtBt, self.PtBt, 0, 0, 0, 0) class BOXCO2: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 4: self.init_medium(args) elif len(args) == 8: self.init_big(args) else: raise Exception, "number of arguments is wrong" def init_noargs(self): self.low = [0, 0] self.high = [0, 0] self.low_open = [1, 1] self.high_open = [1, 1] def init_big(self, args): l_s, h_s, l_t, h_t, l_open_s, h_open_s, l_open_t, h_open_t = args self.init_noargs() self.low = [l_s, l_t] self.high = [h_s, h_t] self.low_open = [l_open_s, l_open_t] self.high_open = [h_open_s, h_open_t] def init_medium(self, args): l, h, l_open, h_open = args self.low = l #copy(l) ?? self.high = h self.low_open = l_open self.high_open = h_open def merge(self, b): c = BOXCO2() i = 0 while i < 2: if self.low[i] < b.low[i]: c.low[i] = self.low[i] c.low_open[i] = self.low_open[i] elif self.low[i] > b.low[i]: c.low[i] = b.low[i] c.low_open[i] = b.low_open[i] else: #if self.low[i] == b.low[i] c.low[i] = self.low[i] c.low_open[i] = self.low_open[i] * b.low_open[i] if self.high[i] > b.high[i]: c.high[i] = self.high[i] c.high_open[i] = self.high_open[i] elif self.high[i] < b.high[i]: c.high[i] = b.high[i] c.high_open[i] = b.high_open[i] else: #if self.high[i] == b.high[i] c.high[i] = self.high[i] c.high_open[i] = self.high_open[i] * b.high_open[i] i += 1 return c class Point2D: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 2: self.init_xy(args[0], args[1]) else: raise Exception, "number of arguments is wrong" def init_noargs(self): self.type = 0 self.PtRs, self.PtRt, self.PtBs, self.PtBt = 0, 0, 0, 0 self.poly1, self.poly2 = K_RATPOLY(), K_RATPOLY() self.prev, self.next = self, self self.pt_in_other_dom = 0 def init_xy(self, x, y): self.init_noargs() self.type = 4 self.PtRs = 0 self.PtRt = 0 self.PtBs = x self.PtBt = y self.poly1 = K_RATPOLY(1, 0, x) self.poly2 = K_RATPOLY(1, 0, y) def __repr__(self): return "Point2D(type=%d, poly1=%s, poly2=%s)" % (self.type, str(self.poly1), str(self.poly2)) def get_low_s(self): assert (self.type) if self.type == 1 or self.type == 2: return self.PtRs.low else: #if self.type == 2 or self.type == 4 return self.PtBs def get_high_s(self): assert (self.type) if self.type == 1 or self.type == 2: return self.PtRs.high else: return self.PtBs def get_low_t(self): assert (self.type) if self.type == 1 or self.type == 3: return self.PtRt.low else: #if self.type == 2 or self.type == 4 return self.PtBt def get_high_t(self): assert (self.type) if self.type == 1 or self.type == 3: return self.PtRt.high else: #if self.type == 2 or self.type == 4 return self.PtBt def get_low(self): assert (self.type) if self.type == 1: h = [self.PtRs.low, self.PtRt.high] elif self.type == 2: h = [self.PtRs.low, self.PtBt] elif self.type == 3: h = [self.PtBs, self.PtRt.low] else: h = [self.PtBs, self.PtBt] #if self.type == 4 return h def get_high(self): assert(self.type) if self.type == 1: h = [self.PtRs.high, self.PtRt.high] elif self.type == 2: h = [self.PtRs.high, self.PtBt] elif self.type == 3: h = [self.PtBs, self.PtRt.high] else: h = [self.PtBs, self.PtBt] #if self.type == 4 return h def cut_s(self, b_s): #NOTE: cut_t and cut_s can be refactored to reduce code """ int K_POINT2D :: cut_s(const bigrational& b_s) const cut *this by the line s = b_s, i.e., refine the box for *this by setting get_low_s() or get_high_s() to be b_s. """ assert(self.type) if (self.type == 1 or self.type == 2) and self.PtRs.low < b_s < self.PtRs.high: assert(self.PtRs.num_roots == 1) if not self.PtRs.poly.sgn_at(b_s): if self.type == 1: self.type = 3 else: self.type = 4 #if self.type == 2 self.PtRs = 0 self.PtBs = b_s else: #if self.PtRt.poly.sgn_at(b_t) self.PtRs.cut(b_t) return 0 def cut_t(self, b_t): """ int K_POINT2D :: cut_t(const bigrational& b_t) const cut *this by the line t = b_t, i.e., refine the box for *this by setting get_low_t() or get_high_t() to be b_t. """ assert(self.type) if (self.type == 1 or self.type == 3) and self.PtRt.low < b_t < self.PtRt.high: assert(self.PtRt.num_roots == 1) if not self.PtRt.poly.sgn_at(b_t): if self.type == 1: self.type = 2 else: self.type = 4 #if self.type == 3 self.PtRt = 0 self.PtBt = b_t else: #if self.PtRt.poly.sgn_at(b_t) self.PtRt.cut(b_t) return 0 def halve_s(self): """ int K_POINT2D :: halve_s() const cut *this by the line s = half that halves the box for *this, i.e., refine the box for *this by setting get_low_s() or get_high_s() to be half. """ raise NotImplementedError, bryan_message def halve_t(self): """ int K_POINT2D :: halve_t() const cut *this by the line t = half that halves the box for *this, i.e., refine the box for *this by setting get_low_t() or get_high_t() to be half. """ raise NotImplementedError, bryan_message def reduce_s(self, num_bits): """ int K_POINT2D :: reduce_s(const unsigned long num_bits) const reduce *this s.t. [get_low_t(), get_high_t()] will contain all the numbers approximated by PtRs->float_est to num_bits precision. return 1 if some reduction occurs and 0 otherwise. """ raise NotImplementedError, bryan_message def reduce_t(self, num_bits): """ int K_POINT2D :: reduce_t(const unsigned long num_bits) const reduce *this s.t. [get_low_t(), get_high_t()] will contain all the numbers approximated by PtRt->float_est to num_bits precision. return 1 if some reduction occurs and 0 otherwise. """ raise NotImplementedError, bryan_message def reduce(self, num_bits_s, num_bits_t): """ int K_POINT2D :: reduce(const unsigned long num_bits_s, const unsigned long num_bits_t) const reduce *this s.t. [get_low_s(), get_high_s()] will contain all the numbers approximated by PtRs->float_est to num_bits_s precision and [get_low_t(), get_high_t()] will contain all the numbers approximated by PtRt->float_est to num_bits_t precision. return 1 if some reduction occurs and 0 otherwise. """ raise NotImplementedError, bryan_message def contract_s(self, tol): """ int K_POINT2D :: contract_s(const bigrational& tol) const contract *this s.t. [get_low_s(), get_high_s()] will be no larger than tol. """ raise NotImplementedError, bryan_message def contract_t(self, tol): """ int K_POINT2D :: contract_t(const bigrational& tol) const contract *this s.t. [get_low_t(), get_high_t()] will be no larger than tol. """ raise NotImplementedError, bryan_message def contract(self, tol_s, tol_t): """ int K_POINT2D :: contract(const bigrational& tol_s, const bigrational& tol_t) const contract *this s.t. [get_low_s(), get_high_s()] will be no larger than tol_s and [get_low_t(), get_high_t()] will be no larger than tol_t. """ raise NotImplementedError, bryan_message def shrink(self, fac_s, fac_t): """ int K_POINT2D :: shrink(const bigrational& fac_s, const bigrational& fac_t) const shrink *this s.t. [get_low_s(), get_high_s()] will become smaller by fac_s and [get_low_t(), get_high_t()] will become smaller by fac_t. """ raise NotImplementedError, bryan_message def compare_s(self, param): assert(self.type) if isinstance(param, Point2D): x = param assert(x.type) c = 0 if self.type == 3 or self.type == 4: if self.PtBs > x.PtBs: c = 1 elif self.PtBs < x.PtBs: c = -1 else: c = 0 #if self.PtBs == x.PtBs else: #if self.type == 1 or self.type == 2 if x.type == 3 or x.type == 4: if self.PtBs >= x.PtRs.high: c = 1 elif self.PtBs <= x.PtRs.low: c = -1 else: #if x.PtRs.low < self.PtBs < x.PtRs.high x.cut_s(PtBs) c = self.compare_s(x) else: #if x.type == 1 or x.type == 2 if self.PtRs.low > x.PtRs.low and self.PtRs.low < x.PtRs.high: x.cut_s(self.PtRs.low) c = self.compare_s(x) elif self.PtRs.high > x.PtRs.low and self.PtRs.high < x.PtRs.high: x.cut_s(self.PtRs.high) c = self.compare_s(x) elif self.PtRs.low < x.PtRs.low and self.PtRs.high > x.PtRs.low: self.cut_s(x.PtRs.low) c = self.compare_s(x) elif self.PtRs.low < x.PtRs.high and self.PtRs.high > x.PtRs.high: self.cut_s(x.PtRs.high) c = self.compare_s(x) elif self.PtRs.low >= x.PtRs.high: c = 1 elif self.PtRs.high <= x.PtRs.low: c = -1 else: c = 0 #if self.PtRs.low == x.PtRs.low and self.PtRs.high == x.PtRs.high return c else: b_s = param c = 0 if self.type == 3 or self.type == 4: if self.PtBs > b_s: c = 1 elif self.PtBs < b_s: c = -1 else: c = 0 #if self.PtBs == b_s else: #if self.type == 1 or self.type == 2 if self.PtRs.low >= b_s: c = 1 elif self.PtRs.high <= b_s: c = -1 else: #if self.PtRs.low < b_s < self.PtRs.high self.cut_s(b_s) c = self.compare_s(b_s) return c def compare_t(self, param): assert(self.type) if isinstance(param, Point2D): x = param assert(x.type) if self.type == 2 or self.type == 4: if x.type == 2 or x.type == 4: if self.PtBt > x.PtBt: c = 1 elif self.PtBt < x.PtBt: c = -1 else: c = 0 else: #if x.type == 1 or x.type == 3 if self.PtBt >= x.PtRt.high: c = 1 elif self.PtBt <= x.PtRt.low: c = -1 else: #if x.PtRt.low < PtBt < x.PtRt.high x.cut_t(PtBt) c = self.compare_t(x) else: #if self.type == 1 or self.type == 3 if x.type == 2 or x.type == 4: if self.PtRt.low >= x.PtBt: c = 1 elif self.PtRt.high <= x.Ptbt: c = -1 else: #if self.PtRt.low < x.PtBt < self.PtRt.high self.cut_t(x.PtBt) c = self.compare_t(x) else: #if x.type == 1 or x.type == 3 if self.PtRt.low > x.PtRt.low and PtRt.low < x.PtRt.high: x.cut_t(self.PtRt.low) c = self.compare_t(x) elif self.PtRt.high > x.PtRt.low and self.PtRt.high < x.PtRt.high: x.cut_t(self.PtRt.high) c = self.compare_t(x) elif self.PtRt.low < x.PtRt.low and self.PtRt.high > x.PtRt.low: self.cut_t(x.PtRt.low) c = self.compare_t(x) elif self.PtRt.low < x.PtRt.high and slef.PtRt.high > x.PtRt.high: self.cut_t(x.PtRt.high) c = self.compare_t(x) elif self.PtRt.low >= x.PtRt.high: c = 1 elif self.PtRt.high <= x.PtRt.low: c = -1 else: #if self.PtRt.low == x.PtRt.low and self.PtRt.high == x.PtRt.high c = 0 return c else: b_t = param c = 0 if self.type == 2 or self.type == 4: if self.PtBt > b_t: c = 1 elif self.PtBt < b_t: c = -1 else: c = 0 else: if self.PtRt.low >= b_t: c = 1 elif self.PtRt.high <= b_t: c = -1 else: self.cut_t(b_t) c = self.compare_t(b_t) return c #see kpoint1d.cc K_POINT1D::compare def compare(self, point): raise NotImplementedError, bryan_message def compare_(self, X, sort_type): if sort_type == "s": return self.compare_s(X) elif sort_type == "t": return self.compare_t(X) elif sort_type == None: return self.compare(X) else: raise Exception, "sort_type isn't s or t" def bbox(self): assert(self.type) if self.type == 1: b2 = BOXCO2(self.PtRs.low, self.PtRs.high, self.PtRt.low, self.PtRt.high, 1, 1, 1, 1) elif self.type == 2: b2 = BOXCO2(self.PtRs.low, self.PtRs.high, self.PtBt, self.PtBt, 1, 1, 0, 0) elif self.type == 3: b2 = BOXCO2(self.PtBs, self.PtBs, self.PtRt.low, self.PtRt.high, 0, 0, 1, 1) else: #if type == 4 b2 = BOXCO2(self.PtBs, self.PtBs, self.PtBt, self.PtBt, 0, 0, 0, 0) return b2 def overlap(self, point): """pretty sure this is right""" assert(self.type) assert(point.type) def o_setter(check): if check: return 1 return 0 ptrs_low_within_other_ptrs_range = x.PtRs.low <= Ptrs.low < x.PtRs.high ptrs_high_within_other_ptrs_range = x.PtRs.low < PtRs.high <= x.PtRs.high other_ptrs_low_within_ptrs_range = PtRs.low <= x.PtRs.low < PtRs.high other_ptrs_high_within_ptrs_range = PtRs.low < x.PtRs.high <= PtRs.high ptrs_interval_overlap = ptrs_low_within_other_ptrs_range \ or ptrs_high_within_other_ptrs_range \ or other_ptrs_low_within_ptrs_range \ or other_ptrs_high_within_ptrs_range ptrt_low_within_other_ptrt_range = x.PtRt.low <= PtRt.low < x.PtRt.high ptrt_high_within_other_ptrt_range = x.PtRt.low < PtRt.high <= x.PtRt.high other_ptrt_low_within_ptrt_range = PtRt.low <= x.PtRt.low < PtRt.high other_ptrt_high_within_ptrt_range = PtRt.low < x.PtRt.high <= PtRt.high ptrt_interval_overlap = ptrt_low_within_other_ptrt_range \ or ptrt_high_within_other_ptrt_range \ or other_ptrt_low_within_ptrt_range \ or other_ptrt_high_within_ptrt_range other_ptbt_within_ptrt_range = PtRt.low < x.PtBt < PtRt.high other_ptbs_within_ptrs_range = PtRs.low < x.PtBs < PtRs.high ptbs_within_other_ptrs_range = x.PtRs.low < PtBs < x.PtRs.high ptbt_within_other_ptrt_range = x.PtRt.low < PtBt < x.PtRt.high o = None if self.type == 1: if x.type == 1: o = o_setter( ptrs_interval_overlap and ptrt_interval_overlap ) elif x.type == 2: o = o_setter( ptrs_interval_overlap and other_ptbt_within_ptrt_range ) elif x.type == 3: o = o_setter( ptrt_interval_overlap and other_ptbs_within_ptrs_range ) elif x.type == 4: o = o_setter( other_ptbs_within_ptrs_range and other_ptbt_within_ptrt_range ) elif self.type == 2: if x.type == 1: o = o_setter( ptrs_interval_overlap and ptbt_within_other_ptrt_range ) elif x.type == 2: o = o_setter( ptrs_interval_overlap and ptbt_equal ) elif x.type == 3: o = o_setter( other_ptbs_within_ptrs_range and ptbt_within_other_ptrt_range ) elif x.type == 4: o = o_setter( other_ptbs_within_ptrs_range and ptbt_equal ) elif self.type == 3: if x.type == 1: o = o_setter( ptbs_within_other_ptrs_range and ptrt_interval_overlap ) elif x.type == 2: o = o_setter( ptbs_within_other_ptrs_range and other_ptbt_within_ptrt_range ) elif x.type == 3: o = o_setter( ptbs_equal and ptrt_interval_overlap ) elif x.type == 4: o = o_setter( ptbs_equal and other_ptbt_within_ptrt_range ) elif self.type == 4: if x.type == 1: o = o_setter( ptbs_within_other_ptrs_range and ptbt_within_other_ptrt_range ) elif x.type == 2: o = o_setter( ptbt_equal and ptbs_within_other_ptrs_range ) elif x.type == 3: o = o_setter( ptbs_equal and ptbt_within_other_ptrt_range ) elif x.type == 4: o = o_setter( ptbs_equal and ptbt_equal ) def equiv(self, point): return self is point def __eq__(self, param): if isinstance(param, Point2D): if self is param: return True else: self.cut_s(param.get_low_s()) self.cut_s(param.get_high_s()) self.cut_t(param.get_low_t()) self.cut_t(param.get_high_t()) param.cut_s(self.get_low_s()) param.cut_s(self.get_high_s()) param.cut_t(self.get_low_t()) param.cut_t(self.get_high_t()) if self.type != param.type: e = 0 elif not self.overlap(param): e = 0 else: if self.type == 1: #see if 2 univariate polynomials # that define the s-coordinate of self (and x) # have 1 and only 1 common root # in the interval for the s-coordinate of self (and x). num_pts_s = self.get_pts(self.PtRs.low, self.PtRs.high, gcd(self.PtRs.poly, param.PtRs.poly), pts_s, 0, 0) if num_pts_s == 1: num_pts_t = self.get_pts(self.PtRt.low, self.PtRt.high, gcd(self.PtRt.poly, param.PtRt.poly), pts_t, 0, 0) if num_pts_t == 1: prev.next = param.next param.next.prev = prev prev = param param.next = self e = 1 else: e = 0 else: e = 0 elif self.type == 2: num_pts_s = self.get_pts(self.PtRs.low, self.PtRs.high, gcd(self.PtRs.poly, param.PtRs.poly), pts_s, 0, 0) if num_pts_s == 1: prev.next = param.next param.next.prev = prev prev = param param.next = self e = 1 else: e = 0 elif self.type == 3: num_pts_t = self.get_pts(self.PtRt.low, self.PtRt.high, gcd(self.PtRt.poly, param.PtRt.poly), pts_t, 0, 0) if num_pts_t == 1: prev.next = param.next param.next.prev = prev prev = param param.next = self e = 1 else: e = 0 else: #self.type == 4 prev.next = param.next param.next.prev = prev prev = param param.next = self e = 1 return e return False def separate(self, point): assert(self.type) assert(point.type) if self.type == 4 and point.type == 3: point.cut_t(self.PtBt) elif self.type == 4 and point.type == 2: point.cut_s(self.PtBs) elif self.type == 4 and point.type == 1: point.cut_s(self.PtBs) point.cut_t(self.PtBt) elif self.type == 3 and point.type == 4: self.cut_t(point.PtBt) elif self.type == 3 and point.type == 3: self.cut_t(point.PtRt.low) self.cut_t(point.PtRt.high) point.cut_t(self.PtRt.low) point.cut_t(self.PtRt.high) elif self.type == 3 and point.type == 2: self.cut_t(point.PtBt) point.cut_s(self.PtBs) elif self.type == 3 and point.type == 1: self.cut_t(point.PtRt.low) self.cut_t(point.PtRt.high) point.cut_s(self.PtBs) point.cut_t(self.PtRt.low) point.cut_t(self.PtRt.high) elif self.type == 2 and point.type == 4: self.cut_s(point.PtBs) elif self.type == 2 and point.type == 3: self.cut_s(point.PtBs) point.cut_t(self.PtBt) elif self.type == 2 and point.type == 2: self.cut_s(point.PtRs.low) self.cut_s(point.PtRs.high) point.cut_s(self.PtRs.low) point.cut_s(self.PtRs.high) elif self.type == 2 and point.type == 1: self.cut_s(point.PtRs.low) self.cut_s(point.PtRs.high) point.cut_s(self.PtRs.low) point.cut_s(self.PtRs.high) point.cut_t(self.PtBt) elif self.type == 1 and point.type == 4: self.cut_s(point.PtBs) self.cut_t(point.PtBt) elif self.type == 1 and point.type == 3: self.cut_s(point.PtBs) self.cut_t(point.PtRt.low) self.cut_t(point.PtRt.high) point.cut_t(self.PtRt.low) point.cut_t(self.PtRt.high) elif self.type == 1 and point.type == 2: self.cut_s(point.PtRs.low) self.cut_s(point.PtRs.high) self.cut_t(point.PtBs) point.cut_s(self.PtRs.low) point.cut_s(self.PtRs.high) elif self.type == 1 and point.type == 1: self.cut_s(point.PtRs.low) self.cut_s(point.PtRs.high) self.cut_t(point.PtRt.low) self.cut_t(point.PtRt.high) point.cut_s(self.PtRs.low) point.cut_s(self.PtRs.high) point.cut_t(self.PtRt.low) point.cut_t(self.PtRt.high) while self.overlap(point): self.shrink(shrink_step, shrink_step) point.shrink(shrink_step, shrink_step) self.separate(x) return 0 def get_all_pts(self, ratpoly1, ratpoly2, points, tol): """ computes all the intersection points of 2 bivariate polynomials P1 and P2 returns the number of intersections if tol > 0 then a box for each intersection is no larger than tol in any direction. if tol = 0 then there is no limit on the size of boxes for intersections. """ P1, P2 = ratpoly1, ratpoly2 assert P1.get_num_vars() == 2 assert P2.get_num_vars() == 2 h_s = self.GoodSylvester(P1, P2, 1).get_Mignotte_bd() + s_epsilon l_s = - h_s - s_epsilon h_t = self.GoodSylvester(P1, P2, 0).get_Mignotte_bd() + t_epsilon l_t = - h_t - t_epsilon return self.get_pts(l_s, h_s, l_t, h_t, P1, P2, points, tol, 0) def get_fp_approx(self, a): if self.type == 1 or self.type == 2: if self.PtRs.ok_float: a[0] = self.PtRs.float_est else: a[0] = ((self.PtRs.low + self.PtRs.high) / 2).as_double() else: a[0] = self.PtBs.as_double() if self.type == 1 or self.type == 3: if self.PtRt.ok_float: a[1] = self.PtRt.float_est else: a[1] = ((self.PtRt.low + self.PtRt.high) / 2).as_double() else: a[1] = self.PtBt.as_double() return 0 def equal(self, point): return self == point def equal_s(self, point): """returns 1 if self and x are equal in s-coordinate, and 0 otherwise""" assert self.type assert point.type if self.equiv(point): e = 1 else: #if not self is point #cut self and point s.t. one is entirely contained in the other self.cut_s(point.get_low_s()) self.cut_s(point.get_high_s()) point.cut_s(self.get_low_s()) point.cut_s(self.get_high_s()) if self.compare(point): e = 0 elif (self.type == 1 or self.type == 2) and (point.type == 1 or point.type == 2): #see if 2 univariate polynomials # that define the s-coordinate of self (and x) # have 1 and only 1 common root # in the interval for the s-coordinate of self (and point). num_pts = self.get_pts(self.PtRs.low, self.PtRs.high, gcd(self.PtRs.poly, point.PtRs.poly), pts, 0, 0) if num_pts == 1: e = 1 else: e = 0 #if num_pts == 0 or num_pts > 1 elif (self.type == 3 or self.type == 4) and (point.type == 3 or point.type == 4): if self.PtBs == point.PtBs: e = 1 else: e = 0 #if self.PtBs != point.PtBs else: #if ((self.type == 1 or self.type == 2) and (point.type == 3 or point.type == 4) # or (self.type == 3 or self.type == 4) and (point.type == 1 or point.type == 2)) e = 0 return e def equal_t(self, point): """returns 1 if self and x are equal in the t-coordinate and 0 otherwise""" assert self.type assert point.type if self.equiv(point): e = 1 else: #if not self is point #cut self and point s.t. one is entirely contained in the other self.cut_t(point.get_low_t()) self.cut_t(point.get_high_t()) point.cut_t(self.get_low_t()) point.cut_t(self.get_high_t()) if self.compare(point): e = 0 elif (self.type == 1 or self.type == 3) and (point.type == 1 or point.type == 3): #see if 2 univariate polynomials # that define the s-coordinate of self (and x) # have 1 and only 1 common root # in the interval for the t-coordinate of self (and point). num_pts = self.get_pts(self.PtRt.low, self.PtRt.high, gcd(self.PtRt.poly, point.PtRt.poly), pts, 0, 0) if num_pts == 1: e = 1 else: e = 0 #if num_pts == 0 or num_pts > 1 elif (self.type == 2 or self.type == 4) and (point.type == 2 or point.type == 4): if self.PtBt == point.PtBt: e = 1 else: e = 0 #if self.PtBt != point.PtBt else: #if ((self.type == 1 or self.type == 3) and (point.type == 2 or point.type == 4) # or (self.type == 2 or self.type == 4) and (point.type == 1 or point.type == 3)) e = 0 return e def get_pts(self, l_s, h_s, l_t, h_t, P1, P2, pts, tol, count_edges): """ unsigned long get_pts(const bigrational& l_s, const bigrational& h_s, const bigrational& l_t, const bigrational& h_t, const K_RATPOLY& P1, const K_RATPOLY& P2, K_POINT2D**& pts, const bigrational& tol, const int count_edges) computes the intersections "pts" of 2 bivariate polynomials P1 and P2 on/in the region [l_s, h_s] x [l_t, h_t]. returns the number of intersections. if tol > 0 then a box for each intersection is no larger than tol in any direction. if tol = 0 then there is no limit on the size of boxes for intersections. if count_edges = 1 then the intersections on the boundary edges of the region are counted. if count_edges = 0 then the intersections on the boundary edges of the region are ignored. """ assert(P1.num_vars == 2) assert (P2.num_vars == 2) num_pts = 0 if not P1.deg[0] and not P1.deg[1] or not P2.deg[0] and not P2.deg[1]: #0.1 when P1 or P2 is a constant polynomial pts, num_pts = 0, 0 elif P1.eq_upto_const(P2): #0.2 when P1 and P2 are the same (up to constant multiple) pts, num_pts = 0, 0 elif l_s > h_s or l_t > h_t: #0.3 when the region [l_s, l_t] x [h_s, h_t] is not well-defined pts, num_pts = 0, 0 elif l_s == h_s and l_t == h_t: #1. when the region [l_s, l_t] x [h_s, h_t] is shrunken to the point (l_s, l_t) if count_edges: b = [] b.append(l_s) b.append(l_t) if not P1.sgn_at(b) and not P2.sgn_at(b): #if (l_s, l_t) is a common root of P1 and P2 then # pts consists of a single point, namely (l_s, l_t) num_pts = 1 pts = [None] * num_pts pts[0] = Point2D(l_s, h_s, P1, P2) else: #if P1.sgn_at(b) or P2.sgn_at(b) pts, num_pts = 0, 0 else: #if not count_edges pts, num_pts = 0, 0 elif l_s == h_s and l_t < h_t: #2. when the region [l_s, l_t] x [h_s, h_t] is # shrunken to the vertical line segment [(l_s, l_t), (l_s, h_t)] if count_edges: #see if 2 univariate polynomials P1 | s = l_s and P2 | s = l_s # have common roots in interval [l_t, h_t] num_pts = self.get_pts_proto(l_s, l_t, h_t, gcd(P1.subst_val(0, l_s), P2.subst_val(0, l_s)), P1, P2, pts, tol, 1) else: #if not count_edges pts, num_pts = 0, 0 elif l_s < h_s and l_t == h_t: #3. when the region [l_s, l_t] x [h_s, h_t] is # shrunken to the horizontal line segment [(l_s, l_t), (h_s, l_t)] if count_edges: self.get_pts_proto(l_s, h_s, gcd(P1.subst_val(1, l_t), P2.subst_val(1, l_t)), l_t, P1, P2, pts, tol, 1) else: pts, num_pts = 0, 0 elif P1.deg[0] == 1 and not P1.deg[1]: #4. when P1 = 0 is a vertical line s = b1_s b1_s = - P1.coeffs[1] / P1.coeffs[0] #if b1_s is in/on the interval [l_s, h_s] then # see if a univariate polynomial P2 | s = b1_s # has roots in/on the interval [l_t, h_t]. if b1_s > l_s and b1_s < h_s or count_edges and b1_s >= l_s and b1_s <= h_s: num_pts = self.get_pts_proto(b1_s, l_t, h_t, P2.subst_val(0, b1_s), P1, P2, pts, tol, count_edges) else: #if b1_s is not in/on the interval [l_s, h_s] pts, num_pts = 0, 0 elif not P1.deg[0] and P1.deg[1] == 1: #5. when P1 is a horizontal line t = b1_t b1_t = - P1.coeffs[1] / P1.coeffs[0] #if b1_t is in/on the interval [l_t, h_t] then # see if a univariate polynomial P2 | t = b1_t # has roots in/on the interval [l_s, h_s] if b1_t > l_t and b1_t < h_t or count_edges and b1_t >= l_t and b1_t <= h_t: num_pts = self.get_pts_proto(l_s, h_s, P2.subst_val(1, b1_t), b1_t, P1, P2, pts, tol, count_edges) else: #if b1_t is not in/on the interval [l_t, h_t] pts, num_pts = 0, 0 elif P2.deg[0] == 1 and not P2.deg[1]: #6. when P2 is a vertical line s = b2_s b2_s = - P2.coeffs[1] / P2.coeffs[0] #if b2_s is in/on the interval [l_s, h_s] then # see if a univariate polynomial P2 | s = b2_s # has roots in/on the interval [l_t, h_t] if b2_s > l_s and b2_s < h_s or count_edges and b2_s >= l_s and b2_s <= h_s: num_pts = self.get_pts_proto(b2_s, l_t, h_t, P1.subst_val(0, b2_s), P1, P2, pts, tol, count_edges) else: #if b2_s is not in/on the interval [l_s, h_s] pts, num_pts = 0, 0 elif not P2.deg[0] and P2.deg[1] == 1: #7. when P2 is a horizontal line t = b2_t b2_t = -P2.coeffs[1] / P2.coeffs[0] #if b2_t is in/on the interval [l_t, h_t] then # see if a univariate polynomial P2 | t = b2_t # has roots inon the interval [l_s, h_s] if b2_t > l_t and b2_t < h_t or count_edges and b2_t >= l_t and b2_t <= h_t: num_pts = self.get_pts_proto(l_s, h_s, P1.subst_val(1, b2_t), b2_t, P1, P2, pts, tol, count_edges) else: #if b2_t is not in/on the interval [l_t, h_t] pts, num_pts = 0, 0 elif P1.get_total_deg() == 1 and P2.get_total_deg() == 1: #8. When P1 and P2 are lienar, say, # P1 = c11 s + c12 t + c13 and P2 = c21 s + c22 t + c23, # use Cramer's formula to solve the system P1 = P2 = 0 and # see if it has a solution in/on the region [l_s, h_s] x [l_t, h_t] c11, c12, c13 = P1.coeffs[1], P1.coeffs[2], P1.coeffs[3] c21, c22, c23 = P2.coeffs[1], P2.coeffs[2], P2.coeffs[3] d = c11 * c22 - c12 * c21 if sgn(d): b_s = (- c13 * c22 + c12 * c23) / d b_t = (- c11 * c23 + c13 * c21) / d if b_s > l_s and b_s < h_s and b_t > l_t and b_t < h_t \ or count_edges and b_s >= l_s and b_s <= h_s and b_t >= l_t and b_t <= h_t: num_pts = 1 pts = [] pts[0] = Point2D(b_s, b_t, P1, P2) else: #if the solution (b_s, b_t) is not # in/on the region [l_s, h_s] x [l_t, h_t] pts, num_pts = 0, 0 else: #if d == 0, i.e. P1 = 0 and P2 = 0 are parallel pts, num_pts = 0, 0 else: #9. Otherwise, perform 2-D root finding. # We look the following objects: # # (l_s, h_t) edge h_t (h_s, h_t) # \ / # *----------* # | | # edge | | edge # l_s | | h_s # | | # | | # *----------* # / \ # (l_s, l_t) edge l_t (h_s, l_t) # # # CAUTION: Corners and boundary edges are looked at # only if count_edges == 1. if count_edges: #9.1 See if corners are intersections of P1 and P2. b = [] #9.1.1 See if (l_s, l_t) is an intersection of P1 and P2. b[0] = l_s b[1] = l_t if not P1.sgn_at(b) and not P2.sgn_at(b): num_c_l_l = 1 else: #if P1.sgn_at(b) or P2.sgn_at(b) num_c_l_l = 0 #9.1.2 See if (l_s, h_t) is an intersection of P1 and P2. b[0] = l_s b[1] = h_t if not P1.sgn_at(b) and not P2.sgn_at(b): num_c_l_h = 1 else: #if P1.sgn_at(b) or P2.sgn_at(b) num_c_l_h = 0 #9.1.3 See if (h_s, h_t) is an intersection of P1 and P2. b[0] = h_s b[1] = h_t if not P1.sgn_at(b) and not P2.sgn_at(b): num_c_h_h = 1 else: #if P1.sgn_at(b) or P2.sgn_at(b) num_c_h_h = 0 #9.1.4 See if (h_s, l_t) is an intersection of P1 and P2 b[0] = h_s b[1] = l_t if not P1.sgn_at(b) and not P2.sgn_at(b): num_c_h_l = 1 else: #if P1.sgn_at(b) or P2.sgn_at(b) num_c_h_l = 0 #9.2 See if there are intersections of P1 and P2 on boundary edges. #9.2.1 See if there are intersections of P1 and P2 on edge l_t. num_e_l_t = self.get_pts_proto(l_s, h_s, gcd(P1.subst_val(1, l_t), P2.subst_val(1, l_t)), l_t, P1, P2, e_l_t, tol, 0) #9.2.2 See if there are intersections of P1 and P2 on edge l_s. num_e_l_s = self.get_pts_proto(l_s, l_t, h_t, gcd(P1.subst_val(0, l_s), P2.subst_val(0, l_s)), P1, P2, e_l_s, tol, 0) #9.2.3 See if there are intersections of P1 and P2 on edge h_t. num_e_h_t = self.get_pts_proto(l_s, h_s, gcd(P1.subst_val(1, h_t), P2.subst_val(1, h_t)), h_t, P1, P2, e_h_t, tol, 0) #9.2.4 See if there are intersections of P1 and P2 on edge h_s. num_e_h_s = self.get_pts_proto(h_s, l_t, h_t, gcd(P1.subst_val(0, h_s), P2.subst_val(0, h_s)), P1, P2, e_h_s, tol, 0) else: #if not count_edges num_c_l_l = num_c_l_h = num_c_h_h = num_c_h_l = 0 num_e_l_t = num_e_l_s = num_e_h_t = num_e_h_s = 0 #9.3 See if there are intersections of P1 and P2 # in the open region (l_s, l_t) x (h_s, h_t) num_interior = self.get_pts_interior(l_s, l_t, h_s, h_t, P1, P2, interior, tol) #9.4 Gather all together. num_pts = num_c_l_l + num_c_l_h + num_c_h_h + num_c_h_l + num_e_l_t + num_e_l_s + num_e_h_t + num_e_h_s + num_interior if num_pts > 0: pts = [] i = 0 if num_c_l_l > 0: pts.append(Point2D(l_s, l_t, P1, P2)) i += 1 if num_c_l_h > 0: pts.append(Point2D(l_s, h_t, P1, P2)) i += 1 if num_c_h_h > 0: pts.append(Point2D(h_s, h_t, P1, P2)) i += 1 if num_c_h_l > 0: pts.append(Point2D(h_s, l_t, P1, P2)) i += 1 #set i to len(pts) ?? pts.extend(e_l_t[0:num_e_l_t]) i += num_e_l_t pts.extend(e_l_s[0:num_e_l_s]) i += num_e_l_s pts.extend(e_h_t[0:num_e_h_t]) i += num_e_h_t pts.extend(e_h_s[0:num_e_h_s]) i += num_e_h_s pts.extend(interior[0:num_interior]) i+= num_interior #eh? why are we doing this? #why not just set i = len(pts) assert i == num_pts else: #if num_pts == 0 pts = [] return pts def get_pts_interior(self, l_s, l_t, h_s, h_t, P1, P2, pts, tol): """ unsigned long get_pts_interior(const bigrational& l_s, const bigrational& l_t, const bigrational& h_s, const bigrational& h_t, const K_RATPOLY& P1, const K_RATPOLY& P2, K_POINT2D**& pts, const bigrational& tol) computes the intersections "pts" of 2 bivariate polynomials P1 and P2 in the open region (l_s, h_s) x (l_t, h_t) by using box-hit analysis. returns the number of intersections. if tol > 0 then a box for each intersection is no larger than tol in any direction. if tol = 0 then there is no limit on the size of boxes for intersections. """ s_only = self.GoodSylvester(P1, P2, 1) t_only = self.GoodSylvester(P1, P2, 0) #1. Compute all roots pts_s and pts_t of s_only and t_only # that are supersets of the s- and t- coordinate # of the intersections pts of P1 and P2. pts_s = self.get_pts(l_s, h_s, s_only, None, tol, 0) pts_t = self.get_pts(l_t, h_t, t_only, None, tol, 0) if not len(pts_s): pts, num_pts = 0, 0 elif not len(pts_t): pts, num_pts = 0, 0 else: #if len(pts_s) > 0 and len(pts_t) > 0 #2. Divide pts_s into pts_s1 and pts_s2 # that are the sets of type 1 and 2 points. # Also, divide pts_t into pts_t1 and pts_t2 # that are the sets of type 1 and 2 points. pts_s1, pts_s2 = [], [] #QUESTIONABLE for point in pts_s: if point.type == 1: pts_s1.append(point) else: #if point.type == 2 pts_s2.append(point) assert len(pts_s) == len(pts_s1) + len(pts_s2) #Make sure that the points in pts_s1 are separable. assert sort_(pts_s1) #line 2103 of kpoint2d.cc assert is_separable(pts_s1) pts_t1 = [] pts_t2 = [] for point in pts_t: if point.type == 1: pts_t1.append(point) else: pts_t2.append(point) assert len(pts_t) == len(pts_t1) + len(pts_t2) #Make sure that the points in pts_t1 are separable. assert sort_(pts_t1) assert is_separable(pts_t1) #3. See if a point in pts_s x pts_t # is actually an intersection of P1 and P2. # Let num_Pts_all to be # the number of possible intersections of P1 and P2. # Then, num_pts_all = # min { num_pts_s * num_pts_t # deg_s P1 * deg_t P2 + deg_t P1 * deg_s P2 }. num_pts_all = P1.deg[0] * P2.deg[1] + P1.deg[1] * P2.deg[0] num_pts_st = len(pts_s) * len(pts_t) if num_pts_all > num_pts_st: num_pts_all = num_pts_st pts_proto = [] num_pts = 0 #3.1 type 2 versus type 2 for point in pts_s2: b[0] = point.PtB for point2 in pts_t2: b[1] = point2.PtB if not P1.sgn_at(b) and not P2.sgn_at(b): pts_proto.append(Point2D(b[0], b[1], P1, P2)) #3.2 type 2 versus type 1 if len(pts_t1) > 0: for point in pts_s2: b_s = point.PtB G21 = gcd(P1.subst_val(0, b_s), P2.subst_val(0, b_s)) #(b_s, y) is an intersection of P1 and P2 <=> #y is a common root of P1 | s = b_s and P2 | s = b_s <=> #y is a root of G21 = gcd(P1 | s = b_s, P2 | s = b_s). if G21.deg[0] > 0: pts_21 = self.get_pts(pts_t1[0].get_low(), pts_t1[len(pts_t1) -1].get_high(), G21, None, tol, 0) if len(pts_21) > 0: #Make sure that the points in pts_21 are separable. assert sort_(pts_21) assert is_separable(pts_21) for point001 in pts_21: pts_proto.append(Point2D(b_s, point001, P1, P2)) #3.3 type 1 versus type 2 if len(pts_s1) > 0: for point in pts_t2: b_t = point.PtB G12 = gcd(P1.subst_val(1, b_t), P2.subst_val(1, b_t)) #(x, b_t) is an intersection of P1 and P2 <=> #x is a common root of P1 | t = b_t and P2 | t = b_t <=> #x is a root of G12 = gcd(P1 | t = b_t, P2 | t = b_t). if G12.deg[0] > 0: pts_12 = self.get_pts(pts_s1[0].get_low(), pts_s1[len(pts_s1) - 1].get_high(), G12, pts_12, tol, 0) if len(pts_12) > 0: #Make sure that the points in pts_12 are separable. assert sort(pts_12) assert is_separable(pts_12) for point2 in pts_12: pts_proto.append(Point2D(point2, b_t, P1, P2)) #3.4 type 1 versus type 1 # Perform box hit analysis. # For each side of the box for a point in pts_s1 x pts_t1. # count the number of intersections with P1 and P2. if len(pts_s1) > 0 and len(pts_t1) > 0: num_pts_s1 = len(pts_s1) s_l = [0] * num_pts_s1 s_h = [0] * num_pts_s1 t_l = [0] * num_pts_s1 t_h = [0] * num_pts_s1 i = 0 while i < num_pts_s1: s_l[i] = pts_s1[i].get_low() s_h[i] = pts_s1[i].get_high() i += 1 num_pts_t2 = len(pts_t2) i = 0 while i < num_pts_2: t_l[i] = pts_t1[i].get_low() t_h[i] = pts_t1[i].get_high() i += 1 num_pts_t1 = len(pts_t1) #FIXME: instantiate these arrays better #A = [ [None] * 3 for i in range(3) ] boxhits = [[[0] * MAX_NUM_PTS_HIT] * num_pts_t1] * num_pts_s1 boxhits1 = [[0] * num_pts_t1] * num_pts_s1 boxhits2 = [[0] * num_pts_t1] * num_pts_s1 #low_s: the left side of a box i = 0 while i < num_pts_s1: Q1 = P1.subst_val(0, s_l[i]) pts1 = self.get_pts(t_l[0], t_h[num_pts_t1 - 1], Q1, None, init_tol, 0) if len(pts1) > 0: poly_hit = pts1[0].poly sort_(pts1) match1 = match_intervals(t_l, t_h, num_pts_t1, pts1, None, num_pts1) else: #if len(pts1) == 0 poly_hit = 0 Q2 = P2.subst_val(0, s_l[i]) pts2 = self.get_pts(t_l[0], t_h[num_pts_t1 - 1], Q2, None, init_tol, 0) if len(pts2) > 0: sort_(pts2) match2 = match_intervals(t_l, t_h, num_pts_t1, pts2, None, num_pts2) count1 = count2 = 0 j = 0 while j < num_pts_t1: if not Q1.sgn_at(t_l[j]): b[0] = s_l[i] b[1] = t_l[j] sgn_dtds = - sgn(P1.derivative(0).evaluate(b)) * sgn(P1.derivative(1).evaluate(b)) if sgn_dtds > 0: boxhits[i][j][boxhits1[i][j] + boxhits2[i][j]] = 1 boxhits1[i][j] += 1 elif not sgn_dtds: boxhits1[i][j] += 5 if not Q2.sgn_at(t_l[j]): b[0] = s_l[i] b[1] = t_l[j] sgn_dtds = - sgn(P2.derivative(0).evaluate(b)) * sgn(P2.derivative(1).evaluate(b)) if sgn_dtds > 0: boxhits[i][j][boxhits1[i][j] + boxhits2[i][j]] = 2 boxhits2[i][j] += 1 elif not sgn_dtds: boxhits2[i][j] += 5 pts_hit = [0] * MAX_NUM_PTS_HIT #of Point1D num_pts_hit = 0 num_pts1 = len(pts1) while count1 < num_pts1 and match1[count] < j: count1 += 1 while count1 < num_pts1 and match1[count1] == j: boxhits1[i][j] += 1 count1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 num_pts2 = len(pts2) while count2 < num_pts2 and match2[count2] < j: count2 += 1 while count2 < num_pts2 and match2[count2] == j: boxhits2[i][j] += 1 count2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 num_pts_hit = len(pts_hit) if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) k = 0 while k < num_pts_hit: #FIXME: - l = boxhits1[i][j] + boxhits2[i][j] + num_pts_hit + k if poly_hit and pts_hit[k].poly is poly_hit: #might have to be == not is boxhits[i][j][l] = 1 else: boxhits[i][j][l] = 2 k += 1 j += 1 i += 1 #high_t: the top side of a box i = 0 while i < num_pts_t1: Q1 = P1.subst_val(1, t_h[i]) pts1 = self.get_pts(s_l[0], s_h[num_pts_s1 - 1], Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) match1 = match_intervals(s_l, s_h, num_pts_s1, pts1, None, num_pts1) else: #if num_pts1 == 0 poly_hit = 0 Q2 = P2.subst_val(1, t_h[i]) pts2 = self.get_pts(s_l[0], s_h[num_pts_s1 - 1], Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) match2 = match_intervals(s_l, s_h, num_pts_s1, pts2, None, num_pts2) count1 = count2 = 0 j = 0 while j < num_pts_s1: if not Q1.sgn_at(s_l[j]): b[0] = s_l[j] b[1] = t_h[i] sgn_dtds = - sgn(P1.derivative(0).evaluate(b)) * sgn(P1.derivative(1).evaluate(b)) if sgn_dtds < 0: boxhits[j][i][boxhits1[j][i] + boxhits2[j][i]] = 1 boxhits1[j][i] += 1 elif not sgn_dtds: boxhits1[j][i] += 5 if not Q2.sgn_at(s_l[j]): b[0] = s_l[j] b[1] = t_h[i] sgn_dtds = - sgn(P2.derivative(0).evaluate(b)) * sgn(P2.derivative(1).evaluate(b)) if sgn_dtds < 0: boxhits[j][i][boxhits1[j][i] + boxhits2[j][i]] = 2 boxhits2[j][i] += 1 elif not sgn_dtds: boxhits2[j][i] += 5 pts_hit = [] while (count1 < num_pts1 and match1[count1] < j): count1 += 1 while (count1 < num_pts1 and match1[count1] == j): boxhits1[j][i] += 1 count1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 while (count2 < num_pts2 and match2[count2] < j): count2 += 1 while (count2 < num_pts2 and match2[count2] == j): boxhits2[j][i] += 1 count2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) k = 0 while k < num_pts_hit: l = boxhits1[j][i] + boxhits2[j][i] - num_pts_hit + k if poly_hit and pts_hit[k].poly is poly_hit: boxhits[j][i][l] = 1 else: boxhits[j][i][l] = 2 k += 1 j += 1 i += 1 #high_s: the right side of a box: (counted in reverse order) i = 0 while i < num_pts_s1: Q1 = P1.subst_val(1, s_h[i]) pts1 = self.get_pts(t_l[0], t_h[num_pts_t1 - 1], Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) match1 = match_intervals(t_l, t_h, num_pts_t1, pts1, None, num_pts1) else: #if num_pts1 == 0 poly_hit = 0 Q2 = P2.subst_val(1, s_h[i]) pts2 = self.get_pts(t_l[0], t_h[num_pts_t1 - 1], Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) match2 = match_intervals(t_l, t_h, num_pts_t1, pts2, None, num_pts2) count1 = count2 = 0 j = 0 while j < num_pts_s1: if not Q1.sgn_at(t_h[j]): b[0] = s_h[i] b[1] = t_h[j] sgn_dtds = - sgn(P1.derivative(0).evaluate(b)) * sgn(P1.derivative(1).evaluate(b)) if sgn_dtds > 0: boxhits[i][j][boxhits1[i][j] + boxhits2[i][j]] = 1 boxhits1[i][j] += 1 elif not sgn_dtds: boxhits1[i][j] += 5 if not Q2.sgn_at(t_h[j]): b[0] = s_h[i] b[1] = t_h[j] sgn_dtds = - sgn(P2.derivative(0).evaluate(b)) * sgn(P2.derivative(1).evaluate(b)) if sgn_dtds > 0: boxhits[i][j][boxhits1[i][j] + boxhits2[i][j]] = 2 boxhits2[i][j] += 1 #line 2649 elif not sgn_dtds: boxhits2[i][j] += 5 pts_hit = [] while (count1 < num_pts1 and match1[count1] < j): count1 += 1 while (count1 < num_pts1 and match1[count1] == j): boxhits1[i][j] += 1 count1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 while (count2 < num_pts2 and match2[count2] < j): count2 += 1 while (count2 < num_pts2 and match2[count2] == j): boxhits2[i][j] += 1 count2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) k = 0 while k < num_pts_hit: l = boxhits1[i][j] + boxhits2[i][j] - k - 1 if poly_hit and pts_hit[k].poly is poly_hit: boxhits[i][j][l] = 1 else: boxhits[i][j][l] = 2 k += 1 j += 1 i += 1 #low_t: the bottom side of a box (counted in reverse order) i = 0 while i < num_pts_t1: Q1 = P1.subst_val(1, t_l[i]) pts1 = self.get_pts(s_l[0], s_h[num_pts_s1 - 1], Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) match1 = match_intervals(s_l, s_h, num_pts_s1, pts1, None, num_pts1) else: #if num_pts1 == 0 poly_hit = 0 Q2 = P2.subst_val(1, t_l[i]) pts2 = self.get_pts(s_l[0], s_h[num_pts_s1 - 1], Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) match2 = match_intervals(s_l, s_h, num_pts_s1, pts2, None, num_pts2) count1 = count2 = 0 j = 0 while j < num_pts_s1: if not Q1.sgn_at(s_l[j]): b[0] = s_h[j] b[1] = t_l[i] sgn_dtds = - sgn(P1.derivative(0).evaluate(b)) * sgn(P1.derivative(1).evaluate(b)) if sgn_dtds < 0: boxhits[j][i][boxhits1[j][i] + boxhits2[j][i]] = 1 boxhits1[j][i] += 1 elif not sgn_dtds: boxhits1[j][i] += 5 if not Q2.sgn_at(s_h[j]): b[0] = s_h[j] b[1] = t_l[i] sgn_dtds = - sgn(P2.derivative(0).evaluate(b)) * sgn(P2.derivative(1).evaluate(b)) if sgn_dtds < 0: boxhits[j][i][boxhits1[j][i] + boxhits2[j][i]] = 2 boxhits2[j][i] += 2 elif not sgn_dtds: boxhits2[j][i] += 5 pts_hit = [] while (count1 < num_pts1 and match1[count1] < j): count1 += 1 while (count1 < num_pts1 and match1[count1] == j): boxhits1[j][i] += 1 count1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 while (count2 < num_pts2 and match2[count2] < j): count2 += 1 while (count2 < num_pts2 and match2[count2] == j): boxhits2[j][i] += 1 count2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) #in reverse order k = num_pts_hit - 1 while k >= 0: l = boxhits1[j][i] + boxhits2[j][i] - k - 1 if poly_hit and pts_hit[k].poly is poly_hit: boxhits[j][i][l] = 1 else: boxhits[j][i][l] = 2 k -= 1 j += 1 i += 1 #See which of those boxes have roots in them. i = 0 while i < num_pts_s1: j = 0 while j < num_pts_t1: if boxhits1[i][j] >= 2 and boxhits2[i][j] >= 2: if boxhits1[i][j] + boxhits2[i][j] > 4: if self.refine_interior(pts_s1[i], pts_t1[j], P1, P2, pts_proto[num_pts], tol): num_pts += 1 elif boxhits[i][j][0] == 1 and boxhits[i][j][1] == 2 and \ boxhits[i][j][2] == 1 and boxhits[i][j][3] == 2 or \ boxhits[i][j][0] == 2 and boxhits[i][j][1] == 1 and \ boxhits[i][j][2] == 2 and boxhits[i][j][3] == 1: pts_proto[num_pts] = Point2D(pts_s1[i].PtR, pts_t1[j].PtR, P1, P2) num_pts += 1 j += 1 i += 1 if num_pts > 0: pts = [] pts.extend(pts_proto) else: pts = [] return pts def refine_interior(self, x_s, x_t, P1, P2, y, tol): """ int refine_interior(K_POINT1D* const x_s, K_POINT1D* const x_t, const K_RATPOLY& P1, const K_RATPOLY& P2, K_POINT2D*& y, const bigrational& tol) """ assert x_s.type assert x_t.type assert P1.num_vars == 2 assert P2.num_vars == 2 r = 0 if x_s.type == 1: x_s.shrink(shrink_step) if x_t.type == 1: x_t.shrink(shrink_step) if x_s.type == 2 and x_t.type == 2: y = Point2D(x_s.PtB, x_t.PtB, P1, P2) r = 1 elif x_s.type == 2 and x_t.type == 1: b_s = x_s.PtB G21 = gcd(P1.subst_val(0, b_s), P2.subst_val(0, b_s)) if G21.deg[0] > 0: pts_21 = self.get_pts(x_t.get_low(), x_t.get_high(), G21, None, tol, 0) num_pts_21 = len(pts_21) assert num_pts_21 <= 1 if num_pts_21 == 1: y = Point2D(b_s, pts_21[0], P1, P2) r = 1 elif x_s.type == 1 and x_t.type == 2: b_t = x_t.PtB G12 = gcd(P1.subst_val(1, b_t), P2.subst_val(1, b_t)) if G12.deg[0] > 0: pts_12 = self.get_pts(x_s.get_low(), x_s.get_high(), G12, pts_12, tol, 0) num_pts_12 = len(pts_12) assert num_pts_12 <= 1 if num_pts_12 == 1: y = Point2D(pts_12[0], b_t, P1, P2) r = 1 else: #if x_s.type == 1 and x_t.type == 1 boxhits = [0] * MAX_NUM_PTS_HIT #low_s Q1 = P1.subst_val(0, x_s.get_low()) pts1 = self.get_pts(x_t.get_low(), x_t.get_high(), Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) else: poly_hit = 0 Q2 = P2.subst_val(0, x_s.get_low()) pts2 = self.get_pts(x_t.get_low(), x_t.get_high(), Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) v = [None, None] if not Q1.sgn_at(x_t.get_low()): v[0] = x_s.get_low() v[1] = x_t.get_low() sgn_dtds = - sgn(P1.derivative(0).evaluate(v)) * sgn(P1.derivative(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 1 boxhits1 += 1 elif not sgn_dtds: boxhits1 += 5 if not Q2.sgn_at(x_t.get_low()): v[0] = x_s.get_low() v[1] = x_t.get_low() sgn_dtds = - sgn(P2.derivative(0).evaluate(v)) * sgn(P2.derivatives(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 2 boxhits2 += 1 elif not sgn_dtds: boxhits2 += 5 pts_hit = [] count1 = 0 while count1 < num_pts1: boxhits1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 count1 += 1 count2 = 0 while count2 < num_pts2: boxhits2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 count2 += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) i = 0 while i < num_pts_hit: j = boxhits1 + boxhits2 - num_pts_hit + i if poly_hit and pts_hit[i].poly is poly_hit: boxhits[j] = 1 else: boxhits[j] = 2 i += 1 #high_t Q1 = P1.subst_val(1, x_t.get_high()) pts1 = self.get_pts(x_s.get_low(), x_s.get_high(), Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) else: poly_hit = 0 Q2 = P2.subst_val(1, x_t.get_high()) pts2 = self.get_pts(x_s.get_low(), x_s.get_high(), Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) v = [None, None] if not Q1.sgn_at(x_s.get_low()): v[0] = x_s.get_low() v[1] = x_t.get_high() sgn_dtds = - sgn(P1.derivative(0).evaluate(v)) * sgn(P1.derivative(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 1 boxhits1 += 1 elif not sgn_dtds: boxhits1 += 5 if not Q2.sgn_at(x_s.get_low()): v[0] = x_s.get_low() v[1] = x_t.get_high() sgn_dtds = - sgn(P2.derivative(0).evaluate(v)) * sgn(P2.derivatives(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 2 boxhits2 += 1 elif not sgn_dtds: boxhits2 += 5 pts_hit = [] count1 = 0 while count1 < num_pts1: boxhits1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 count1 += 1 count2 = 0 while count2 < num_pts2: boxhits2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 count2 += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) i = 0 while i < num_pts_hit: j = boxhits1 + boxhits2 - num_pts_hit + i if poly_hit and pts_hit[i].poly is poly_hit: boxhits[j] = 1 else: boxhits[j] = 2 i += 1 #high_s (in reverse order) Q1 = P1.subst_val(0, x_s.get_high()) pts1 = self.get_pts(x_t.get_low(), x_t.get_high(), Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) else: poly_hit = 0 Q2 = P2.subst_val(0, x_s.get_high()) pts2 = self.get_pts(x_t.get_low(), x_t.get_high(), Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) v = [None, None] if not Q1.sgn_at(x_t.get_high()): v[0] = x_s.get_high() v[1] = x_t.get_high() sgn_dtds = - sgn(P1.derivative(0).evaluate(v)) * sgn(P1.derivative(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 1 boxhits1 += 1 elif not sgn_dtds: boxhits1 += 5 if not Q2.sgn_at(x_t.get_high()): v[0] = x_s.get_high() v[1] = x_t.get_high() sgn_dtds = - sgn(P2.derivative(0).evaluate(v)) * sgn(P2.derivatives(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 2 boxhits2 += 1 elif not sgn_dtds: boxhits2 += 5 pts_hit = [] count1 = 0 while count1 < num_pts1: boxhits1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 count1 += 1 count2 = 0 while count2 < num_pts2: boxhits2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 count2 += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) i = num_pts_hit - 1 while i >= 0: j = boxhits1 + boxhits2 - i - 1 if poly_hit and pts_hit[i].poly is poly_hit: boxhits[j] = 1 else: boxhits[j] = 2 i -= 1 #low_t (in reverse order) Q1 = P1.subst_val(1, x_t.get_low()) pts1 = self.get_pts(x_s.get_low(), x_s.get_high(), Q1, None, init_tol, 0) num_pts1 = len(pts1) if num_pts1 > 0: poly_hit = pts1[0].poly sort_(pts1) else: poly_hit = 0 Q2 = P2.subst_val(1, x_t.get_low()) pts2 = self.get_pts(x_s.get_low(), x_s.get_high(), Q2, None, init_tol, 0) num_pts2 = len(pts2) if num_pts2 > 0: sort_(pts2) v = [None, None] if not Q1.sgn_at(x_s.get_high()): v[0] = x_s.get_high() v[1] = x_t.get_low() sgn_dtds = - sgn(P1.derivative(0).evaluate(v)) * sgn(P1.derivative(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 1 boxhits1 += 1 elif not sgn_dtds: boxhits1 += 5 if not Q2.sgn_at(x_s.get_high()): v[0] = x_s.get_high() v[1] = x_t.get_low() sgn_dtds = - sgn(P2.derivative(0).evaluate(v)) * sgn(P2.derivatives(1).evaluate(v)) if sgn_dtds > 0: boxhits[boxhits1 + boxhits2] = 2 boxhits2 += 1 elif not sgn_dtds: boxhits2 += 5 pts_hit = [] count1 = 0 while count1 < num_pts1: boxhits1 += 1 pts_hit.append(pts1[count1]) num_pts_hit += 1 count1 += 1 count2 = 0 while count2 < num_pts2: boxhits2 += 1 pts_hit.append(pts2[count2]) num_pts_hit += 1 count2 += 1 if not sort(pts_hit): k = 0 while k < num_pts_hit - 1: l = k + 1 while l < num_pts_hit: pts_hit[k].separate(pts_hit[l]) l += 1 k += 1 assert sort(pts_hit) num_pts_hit = len(pts_hit) i = num_pts_hit - 1 while i >= 0: j = boxhits1 + boxhits2 - i - 1 if poly_hit and pts_hit[i].poly is poly_hit: boxhits[j] = 1 else: boxhits[j] = 2 i -= 1 #See which of those boxes have roots in them. if boxhits1 >= 2 and boxhits2 >= 2: if boxhits1 + boxhits2 > 4: r = refine_interior(x_s, x_t, P1, P2, y, tol) elif boxhits[0] == 1 and boxhits[1] == 2 and \ boxhits[2] == 1 and boxhits[3] == 2 \ or \ boxhits[0] == 2 and boxhits[1] == 1 and \ boxhits[2] == 2 and boxhits[3] == 1: y = Point2D(x_s.PtR, x_t.PtR, P1, P2) r = 1 return r def get_pts_proto_1(self): """ unsigned long get_pts_proto(const bigrational& l_s, const bigrational& h_s, const K_RATPOLY& s_only, const bigrational& b_t, const K_RATPOLY& P1, const K_RATPOLY& P2, K_POINT2D**& pts, const bigrational& tol, const int count_edges) computes the intersections "pts" of 2 bivariate polynomials P1 and P2 on/in the horizontal line segment [(l_s, b_t), (h_s, b_t)]. by computing the roots of a univariate polynomial s_only on/in the interval [l_s, h_s]. returns the number of intersections. if tol > 0 then a box for each intersection is no larger than tol in the s-direction. if tol = 0 then there is no limit on the size of boxes for intersections. if count_edges = 1 then the intersections at the endpoints of the line segment are counted. if count_edges = 0 then the intersections at the endpoints of the line segment are ignored. """ raise NotImplementedError, bryan_message def sort_(X, sort_type=None): """insertion sort of X in the (s or t)-coordinate""" i = 0 distinct = True while i < len(X): x = X[i] j = i - 1 c = x.compare_(X[j], sort_type) while j >= 0 and c < 0: X[j + 1] = X[j] j += 1 c = x.compare_(X[j], sort_type) X[j+1] = x if not c: distinct = False i += 1 return distinct def sort_s(X): return sort_(X, "s") def sort_t(X): return sort_(X, "t") class Surface: def __init__(self, *args): if len(args) == 1: if isinstance(args[0], K_RATPOLY): self.init_from_impl(args[0]) elif isinstance(args[0], Surface): self.init_from_surface(args[0]) else: raise Exception, "argument not understood" elif len(args) == 5: if isinstance(args[0], K_RATPOLY): self.init_from_rational_polynomials(args) else: raise Exception, "argument not understood" else: raise Exception, "too many arguments" def init_from_impl(self, impl): self.Impl = K_RATPOLY(impl) self.Impl_ok = True self.X, self.Y, self.Z, self.W = K_RATPOLY(), K_RATPOLY(), K_RATPOLY(), K_RATPOLY() self.mon_ok = True def init_from_surface(self, surface): self.Impl_ok = surface.Impl_ok if self.Impl_ok: self.Impl = surface.Impl else: self.Impl = K_RATPOLY() self.mon_ok = surface.mon_ok if self.mon_ok: self.X, self.Y, self.Z, self.W = surface.X, surface.Y, surface.Z, surface.W else: self.X, self.Y, self.Z, self.W = K_RATPOLY(), K_RATPOLY(), K_RATPOLY(), K_RATPOLY() def init_from_rational_polynomials(self, args): impl, x, y, z, w = args[0], args[1], args[2], args[3], args[4] self.Impl = impl if impl: self.Impl_ok = True else: self.Impl_ok = False if x and y and z and w: self.X, self.Y, self.Z, self.W = x, y, z, w self.mon_ok = True else: self.X, self.Y, self.Z, self.W = K_RATPOLY(), K_RATPOLY(), K_RATPOLY(), K_RATPOLY() self.mon_ok = False def get_range(self, l_s, h_s, l_t, h_t): assert self.mon_ok i = None b = K_BOX3D() v_l = [None] * 2 v_h = [None] * 2 X_l, X_h, Y_l, Y_h, Z_l, Z_h, W_l, W_h = 0, 0, 0, 0, 0, 0, 0, 0 v_l[0] = l_s v_l[1] = l_t v_h[0] = h_s v_h[1] = h_t self.X.eval_range_in_out(v_l, v_h, X_l, X_h) self.Y.eval_range_in_out(v_l, v_h, Y_l, Y_h) self.Z.eval_range_in_out(v_l, v_h, Z_l, Z_h) self.W.eval_range_in_out(v_l, v_h, W_l, W_h) if sgn(W_h) < 0: #if W_l <= W_h < 0 if sgn(X_l) > 0: b.low[0] = X_h / W_h b.high[0] = X_l / W_l elif sgn(X_h) < 0: b.low[0] = X_h / W_l b.high[0] = X_l / W_h else: b.low[0] = X_h / W_h b.high[0] = X_l / W_h if sgn(Y_l) > 0: b.low[1] = Y_h / W_h b.high[1] = Y_l / W_l elif sgn(Y_h) < 0: b.low[1] = Y_h / W_l b.high[1] = Y_l / W_h else: b.low[1] = Y_h / W_h b.high[1] = Y_l / W_h if sgn(Z_l) > 0: b.low[2] = Z_h / W_h b.high[2] = Z_l / W_l elif sgn(Z_h) < 0: b.low[2] = Z_h / W_l b.high[2] = Z_l / W_h else: b.low[2] = Z_h / W_h b.high[2] = Z_l / W_h i = 0 while i < 3: b.low_infty[i] = b.high_infty[i] = 0 i += 1 elif sgn(W_l) < 0 and not sgn(W_h): if sgn(X_l) > 0: b.low[0] = X_h / W_l b.low_infty[0] = 0 elif sgn(X_h) < 0: b.high[0] = X_l / W_l b.high_infty[0] = 0 if sgn(Y_l) > 0: b.low[1] = Y_h / W_l b.low_infty[1] = 0 elif sgn(Y_h) < 0: b.high[1] = Y_l / W_l b.high_infty[1] = 0 if sgn(Z_l) > 0: b.low[2] = Z_h / W_l b.low_infty[2] = 0 elif sgn(Z_h) < 0: b.high[2] = Z_l / W_l b.high_infty[2] = 0 elif not sgn(W_l) and not sgn(W_h): if sgn(X_l) > 0: b.low_infty[0] = 1 elif sgn(X_h) < 0: b.high_infty[0] = -1 if sgn(Y_l) > 0: b.low_infty[1] = 1 elif sgn(Y_h) < 0: b.high_infty[1] = -1 if sgn(Z_l) > 0: b.low_infty[2] = 1 elif sgn(Z_h) < 0: b.high_infty[2] = -1 elif not sgn(W_l) and sgn(W_h) > 0: if sgn(X_l) > 0: b.low[0] = X_l / W_h b.low_infty[0] = 0 elif sgn(X_h) < 0: b.high[0] = X_h / W_h b.high_infty[0] = 0 if sgn(Y_l) > 0: b.low[1] = Y_l / W_h b.low_infty[1] = 0 elif sgn(Y_h) < 0: b.high[1] = Y_h / W_h b.high_infty[1] = 0 if sgn(Z_l) > 0: b.low[2] = Z_l / W_h b.low_infty[2] = 0 elif sgn(Z_h) < 0: b.high[2] = Z_h / W_h b.high_infty[2] = 0 elif sgn(W_l) > 0: if sgn(X_l) > 0: b.low[0] = X_l / W_h b.high[0] = X_h / W_l elif sgn(X_h) < 0: b.low[0] = X_l / W_l b.high[0] = X_h / W_h else: b.low[0] = X_l / W_l b.high[0] = X_h / W_l if sgn(Y_l) > 0: b.low[1] = Y_l / W_h b.high[1] = Y_h / W_l elif sgn(Y_h) < 0: b.low[1] = Y_l / W_l b.high[1] = Y_h / W_h else: b.low[1] = Y_l / W_l b.high[1] = Y_h / W_l if sgn(Z_l) > 0: b.low[1] = Z_l / W_h b.high[1] = Z_h / W_l elif sgn(Z_h) < 0: b.low[1] = Z_l / W_l b.high[1] = Z_h / W_h else: b.low[1] = Z_l / W_l b.high[1] = Z_h / W_l i = 0 while i < 3: b.low_infty[i] = b.high_infty[i] = 0 i += 1 return b class Patch: def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 1: if isinstance(args[0], Partition): self.init_from_partition(args[0]) elif isinstacne(args[0], Patch): raise Exception, "please just use copy()" elif len(args) == 2: if isinstance(args[0], Surface) and isinstance(args[1], list): self.init_from_surface(args[0], args[1]) else: raise Exception, "argument not understood" else: raise Exception, "too many arguments" def init_noargs(self): self.surf = 0 self.is_head = True self.low_s = 0 self.high_s = 0 self.low_t = 0 self.high_t = 0 self.num_trim_curves = 0 self.trim_curves = [] self.adj_surfs = [] self.adj_patches = [] self.num_int_curves = 0 self.int_curves = [] self.adj_int_surfs = [] self.adj_int_patches = [] self.num_merged = 0 self.ic = 0 self.len_ic = 0 self.ic_closed = 0 def init_from_partition(self, p): self.init_noargs() self.surf = p.From.surf self.is_head = p.is_head self.num_trim_curves = p.num_trim_curves if self.num_trim_curves > 0: i = 0 while i < self.num_trim_curves: trim_curve = p.trim_curves[i] self.trim_curves.append(trim_curve) trim_curve.assoc(trim_curve.curve_in_other_dom, trim_curve.dir_in_other_dom) if p.rev_curves[i]: trim_curve.reverse() self.adj_surfs.append(0) self.adj_patches.append(0) i += 1 self.set_range() def init_from_surface(self, s, t, h=True): #surface = s; curve = t; self.init_noargs() self.surf = s self.is_head = h self.num_trim_curves = len(t) if self.num_trim_curves > 0: i = 0 while i < self.num_trim_curves: self.trim_curves.append(t[i]) self.adj_surfs.append(0) self.adj_patches.append(0) i += 1 self.set_range() def __repr__(self): return "Patch(trim_curves=" + str(len(self.trim_curves)) + ", adj_surfs=" + str(len(self.adj_surfs)) + ", adj_patches=" + str(len(self.adj_patches)) + ")" def set_range(self): """sets low_s, high_s, low_t and high_t""" #locals: i, b, diff_s, diff_t if self.num_trim_curves > 0: print self.trim_curves[0].__dict__ b = self.trim_curves[0].bbox() self.low_s = b.low[0] self.high_s = b.high[0] self.low_t = b.low[1] self.high_t = b.high[1] i = 0 while i < self.num_trim_curves: b = self.trim_curves[i].bbox() if b.low[0] < self.low_s: self.low_s = b.low[0] if b.high[0] > self.high_s: self.high_s = b.high[0] if b.low[1] < self.low_t: self.low_t = b.low[1] if b.high[1] > self.high_t: self.high_t = b.high[1] i += 1 diff_s = self.high_s - self.low_s assert sgn(diff_s) > 0 diff_t = self.high_t - self.low_t assert sgn(diff_t) > 0 diff_s /= 10 diff_t /= 10 self.low_s -= diff_s self.high_s += diff_s self.low_t -= diff_t self.high_t += diff_t return 0 def get_high_s(self): raise NotImplementedError, bryan_message def get_low_s(self): raise NotImplementedError, bryan_message def get_high_t(self): raise NotImplementedError, bryan_message def get_low_t(self): raise NotImplementedError, bryan_message def in_dom(self, point2d): raise NotImplementedError, bryan_message def contains(self, point2d): raise NotImplementedError, bryan_message def in_on_out(self, point2d): raise NotImplementedError, bryan_message def intersect(self, patch): #1. Compute intersection curves; # int_curve1 on this domain and int_curve2 on p's domain. int_curve1, int_curve2 = K_RATPOLY(), K_RATPOLY() v_l, v_h = [None, None] v_min, v_max = 0, 0 print "debug0: Patch.intersect" print "\tsurf.Impl = ", self.surf.Impl print "\tpatch.surf.Impl = ", patch.surf.Impl print "\t------------------------------------------" print "\t surf.X = ", self.surf.X print "\t surf.Y = ", self.surf.Y print "\t surf.Z = ", self.surf.Z print "\t surf.W = ", self.surf.W print "\t patch.surf.X = ", patch.surf.X print "\t patch.surf.Y = ", patch.surf.Y print "\t patch.surf.Z = ", patch.surf.Z print "\t patch.surf.W = ", patch.surf.W print "\t------------------------------------------" int_curve1 = patch.surf.Impl.subst_param_expr(self.surf.X, self.surf.Y, self.surf.Z, self.surf.W) print "int_curve1 = ", str(int_curve1) if not int_curve1.degree[0] and not int_curve1.degree[1] and not sgn(int_curve1.coeffs[0]): #if int_curve1 is 0 raise Exception, "Error: kpatch: intersect: Patches overlap and coplanar." v_l = [self.low_s, self.low_t] v_h = [self.high_s, self.high_t] v_min, v_max = int_curve1.eval_range(v_l, v_h) print "int_curve1 at [(" + str(v_l[0]) + ", " + str(v_l[1]) + ") x (" + str(v_h[0]) + ", " + str(v_h[1]) + ")] ranges (" + str(v_min) + ", " + str(v_max) + ")" if sgn(v_min) > 0 or sgn(v_max) < 0: raise Exception, "Error: kpatch: intersect: Patches do not intersect." int_curve2 = self.surf.Impl.subst_param_expr(patch.surf.X, patch.surf.Y, patch.surf.Z, patch.surf.W) print "int_curve2 = ", str(int_curve2) if not int_curve2.degree[0] and not int_curve2.degree[1] and not sgn(int_curve2.coeffs[0]): #if int_curve2 is 0 raise Exception, "Error: kpatch: intersect: Patches overlap and coplanar." v_l = [patch.low_s, patch.low_t] v_h = [patch.high_s, patch.high_t] v_min, v_max = int_curve2.eval_range(v_l, v_h) print "int_curve2 at [(" + str(v_l[0]) + ", " + str(v_l[1]) + ") x (" + str(v_h[0]) + ", " + str(v_h[1]) + ")] ranges (" + str(v_min) + ", " + str(v_max) + ")" if sgn(v_min) > 0 or sgn(v_max) < 0: raise Exception, "Error: kpatch: interesect: Patches do not intersect." print " -------------------------------------" #2. Resolve curve topologies. int_curves1 = gen_curve_topo(int_curve1, self.low_s, self.high_s, self.low_t, self.high_t) num_int_curves1 = len(int_curves1) print "\tnum_int_curves1 = ", num_int_curves1 #still a lot more to transfer over raise NotImplementedError, bryan_message def split_trim_curves(self, t, s): raise NotImplementedError, bryan_message def merge_curves(self): raise NotImplementedError, bryan_message def split_loops(self, new_patches): raise NotImplementedError, bryan_message def get_point_in(self): raise NotImplementedError, bryan_message def overlap(self, other): t_box = self.surf.get_range(self.low_s, self.high_s, self.low_t, self.high_t) s_box = other.surf.get_range(other.low_s, other.high_s, other.low_t, other.high_t) return t_box.overlap(s_box) class Graph: """ python conversion complete, minor simplifications not tested not refactored (but it should be..) """ def __init__(self, *args): """ constructor options; Graph() Graph(num_vertices) Graph(vertices) """ if len(args) == 0: self.init_noargs() elif len(args) == 1: if type(args[0]) == int: vertex_count = args[0] if vertex_count > 0: self.vertices = range(0, vertex_count) self.init_edges(count=vertex_count) else: self.init_noargs() elif type(args[0]) == list: vertices = args[0] if len(vertices) > 0: self.vertices = copy(vertices) self.init_edges() else: self.init_noargs() else: raise Exception, "argument type not understood" else: raise Exception, "too many arguments" def init_noargs(self): self.vertices = [] self.edges = [] def init_edges(self, count=None): if not count: count = len(self.vertices) self.edges = [[0] * count] * count def __repr__(self): return "Graph(vertex_count=" + len(self.vertices) + ")" def output(self): result = "" if len(self.vertices) > 0: result += "\n ------" v = 1 while v < len(self.vertices): result += "-------" v += 1 result += "-- \n" v = 0 while v < len(self.vertices): w = 0 while w < (len(self.vertices) - 1): #result.width(2) result += " | " + self.edges[v][w] w += 1 #result.width(2) result += " | " + self.edges[v][len(self.vertices) - 1] + " | \n" v += 1 result += "\n ------" v = 1 while v < len(self.vertices): result += "-------" v += 1 result += "-- \n" else: result = " NULL " return result def ID_to_index(self, id): assert len(self.vertices) > 0, "Graph must have vertices to use ID_to_index." return self.vertices.index(id) #v = 0 #while (v < len(self.vertices) and self.vertex[v] != id): v += 1 #return v def add_edge(self, x, y, ps): v = self.ID_to_index(x) w = self.ID_to_index(y) assert (v < len(self.vertices) and w < len(self.vertices)), "v and w must be in the Graph's vertices" if not ps: self.edges[v][w] = edges[w][v] = 1 else: edges[v][w] = edges[w][v] = -1 return 0 def del_edge(self, x, y): v = self.ID_to_index(x) w = self.ID_to_index(y) assert (v < len(self.vertices) and w < len(self.vertices)), "v and w must be in the Graph's ertices" edges[v][w] = edges[w][v] = 0 return 0 def gen_subgraph(self, n, V): v, w, v_o, w_o = 0, 0, 0, 0 S = Graph(n, V) while v < n: v_o = self.ID_to_index(V[v]) assert (v_o < len(self.vertices)) w = 0 while w < n: w_o = self.ID_to_index(V[w]) assert (w_o < len(self.vertices)) if self.edges[v_o][w_o]: S.edges[v][w] = S.edges[w][v] = 1 w += 1 v += 1 return S def DFS(self, v, visited, C, c): w = 0 visited[v] = True C[v] = c while w < len(self.vertices): if self.edges[v][w] > 0 and not visited[w]: self.DFS(w, visited, C, c) w += 1 return 0 def get_connected_components(self, C, G): v, w, visited, num_components = 0, 0, [], 0 if len(self.vertices) > 0: visited = [False] * len(self.vertices) C = [-1] * len(self.vertices) v = 0 num_components = 0 while v < len(self.vertices): self.DFS(v, visited, C, num_components) num_components += 1 v += 1 while v < len(self.vertices) and visited[v]: v += 1 G = Graph(num_components) v = 0 while v < len(self.vertices): w = v + 1 while w < len(self.vertices): if self.edges[v][w] < 0: G.add_edge(C[v], C[w], 0) w += 1 v += 1 else: #if not num_vertices C = [] G = 0 num_components = 0 #return num_components return C, G def DFS_for_PC(self, v, visited, C, cc): assert (cc == IN or cc == OUT) w, nc = 0, 0 visited[v] = True C[v] = cc if cc == IN: nc = OUT else: nc = IN w = 0 while w < len(self.vertices): if edges[v][w] > 0: if not visited[w]: self.DFS_for_PC(w, visited, C, nc) else: #if w has been visited assert(C[w] == nc) w += 1 return 0 def propagate_color(self, x, C, c): assert (len(self.vertices) > 0) v, w, visited, num_components = 0, 0, [], 0 v = self.ID_to_index(x) assert (v < len(self.vertices)) visited = [False] * len(self.vertices) C = [0] * len(self.vertices) self.DFS_for_PC(v, visited, C, c) return 0 class Partition: PARTITION_ID_COUNT = 0 def __init__(self, *args): if len(args) == 0: self.init_noargs() elif len(args) == 1: if isinstance(args[0], Patch): self.init_from_patch(args[0]) elif isinstance(args[0], Partition): self.init_from_partition(args[0]) else: raise Exception, "argument not understood" else: raise Exception, "too many arguments" def init_noargs(self): Partition.PARTITION_ID_COUNT += 1 self.id = Partition.PARTITION_ID_COUNT self.From = None #Patch self.is_head = True self.trim_curves = [] #Curve self.adj_curves = [] self.adj_surfs = [] self.adj_partitions = [] self.rev_curves = [] def init_from_patch(self, patch): self.init_noargs() self.From = patch if self.From: self.is_head = self.From.is_head self.trim_curves = self.From.trim_curves self.adj_curves = range(0, len(self.From.trim_curves)) self.adj_surfs = [0] * len(self.From.trim_curves) self.adj_partitions = [0] * len(self.From.trim_curves) self.rev_curves = [0] * len(self.From.trim_curves) else: pass #if not self.From; everything was done by init_noargs def init_from_partition(self, partition): self.init_noargs() self.From = partition.From self.is_head = partition.is_head self.trim_curves = partition.trim_curves if len(self.trim_curves) > 0: self.adj_curves = partition.adj_curves self.adj_surfs = partition.adj_surfs self.adj_partitions = partition.adj_partitions self.rev_curves = partition.rev_curves else: pass #init_noargs did enough def __repr__(self): return "Partition(id=" + self.id + ")" def gen_partitions(self, patch): """returns a list of partitions from a given patch""" assert patch, "gen_partitions needs a Patch" #### local variables #i, j, k = 0, 0, 0 #num_curves = 0 #int_end_pts = Point2D() #end_pt = Point2D() #partitions_fw, partitions_rev, list_trim, list_int, num_curves_on_partition, trim_used #i_s, i_c #num_partitions num_curves = patch.num_merged if num_curves > 0: partitions_fw = [-1] * num_curves partitions_rev = [-1] * num_curves int_end_pts = [None] * 2 * patch.num_int_curves i = 0 while i < num_curves: int_end_pts[2 * i] = patch.ic[i][0].start() int_end_pts[2 * i + 1] = patch.ic[i][patch.len_ic[i] - 1].end() i += 1 #list_trim = [[-1] * (num_curves + 1)] * (len(patch.trim_curves) + num_curves) #list_int = [[-1] * (num_curves + 1)] * (len(patch.trim_curves) + num_curves) list_trim = [[-1] * (len(patch.trim_curves) + num_curves)] * (num_curves + 1) list_int = [[-1] * (len(patch.trim_curves) + num_curves)] * (num_curves + 1) num_curves_on_partition = [0] * (num_curves + 1) trim_used = [False] * len(patch.trim_curves) num_partitions = 0 done = False while not done: i = 0 while i < len(patch.trim_curves) and trim_used[i]: i += 1 if i < len(patch.trim_curves): i_is = i go_on = True while go_on: list_trim[num_partitions][num_curves_on_partition[num_partitions]] = i trim_used[i] = True list_int[num_partitions][num_curves_on_partition[num_partitions]] = -1 num_curves_on_partition[num_partitions] += 1 j = 0 while ((j < 2 * num_curves) and not int_end_pts[j].equal(patch.trim_curves[i].end())): j += 1 if j < 2 * num_curves: #the next curve is an intersection curve list_trim[num_partitions][num_curves_on_partition[num_partitions]] = -1 list_int[num_partitions][num_curves_on_partition[num_partitions]] = j/2 num_curves_on_partition[num_partitions] += 1 if j%2: end_pt = patch.ic[j/2][0].start() partitions_rev[j/2] = num_partitions else: #if not j%2 end_pt = patch.ic[j/2][patch.len_ic[j/2]-1].end() partitinos_fw[j/2] = num_partitions #what? reusing the variable i in this loop? i = 0 while i < len(patch.trim_curves) and not end_pt.equal(patch.trim_curves[i].start()): i += 1 assert (i < len(patch.trim_curves)) else: #if (j == 2 * num_curves) #the next curve is a trimming curve i = (i + 1) % len(patch.trim_curves) if i == i_s: go_on = False num_partitions += 1 else: done = True #end while not done assert (num_partitions > 0) partitions = [None] * num_partitions #fill in partitions i = 0 while i < num_partitions: partition = Partition() partition.From = patch partition.is_head = patch.is_head partition.num_trim_curves = 0 j = 0 while j < num_curves_on_partition[i]: if list_trim[i][j] >= 0: partition.num_trim_curves += 1 else: partition.num_trim_curves += patch.len_ic[list_int[i][j]] j += 1 partition.trim_curves = [None] * partition.num_trim_curves partition.adj_curves = [None] * partition.num_trim_curves partition.adj_surfs = [None] * partition.num_trim_curves partition.adj_partitions = [None] * partition.num_trim_curves partition.rev_curves = [0] * partition.num_trim_curves partitions[i] = partition i += 1 #end while i < num_partitions i = 0 while i < num_partitions: i_c = 0 j = 0 while j < num_curves_on_partition[i]: if list_trim[i][j] >= 0: partitions[i].trim_curves[i_c] = patch.trim_curves[list_trim[i][j]] partitions[i].adj_curves[i_c] = list_trim[i][j] partitions[i].adj_surfs[i_c] = 0 partitions[i].adj_partitions[i_c] = 0 partitions[i].rev_curves[i_c] = 0 i_c += 1 else: #if (list_trim[i][j] < 0 and list_int[i][j] >= 0 if i == partitions_fw[list_int[i][j]]: k = 0 while k < patch.len_ic[list_int[i][j]]: partitions[i].trim_curves[i_c] = patch.ic[list_int[i][j]][k] partitions[i].adj_curves[i_c] = -1 partitions[i].adj_surfs[i_c] = 0 partitions[i].adj_partitions[i_c] = partitions[partitions_rev[list_int[i][j]]] partitions[i].rev_curves[i_c] = 0 i_c += 1 k += 1 else: #if i != partitions_fw[list_int[i][j]] k = patch.len_ic[list_int[i][j]] -1 while k >= 0: partitions[i].trim_curves[i_c] = patch.ic[list_int[i][j]][k] partitions[i].adj_curves[i_c] = -1 partitions[i].adj_surfs[i_c] = 0 partitions[i].adj_partitions[i_c] = partitions[partitions_fw[list_int[i][j]]] partitions[i].rev_curves[i_c] = 1 i_c += 1 k -= 1 j += 1 i += 1 #end while i < num_partitions #delete patch's old informations (omitted) patch.num_merged = 0 else: num_partitions = 1 partition = Partition(patch) partitions = [partition] return partitions def gen_adjacency(self, partitions, graph=Graph(), C=0, CG=Graph()): num_partitions = len(partitions) assert (num_partitions > 0) V = [] for partition in partitions: V.append(partition.id) G = Graph(V) for partition in partitions: f = partition.id j = 0 while j < len(partition.trim_curves): a = partition.adj_partitions[j] if a: t = a.id G.add_edge(f, t, 1) else: p = partition.From c = partition.adj_curves[j] a = 0 k = 0 while not a and k < num_partitions: if partitions[k].From == p.adj_patches[c]: l = 0 while l < len(partitions[k].trim_curves) and partitions[k].trim_curves[l] != p.trim_curves[c].curve_in_other_dom: l += 1 if l < len(partitions[k].trim_curves): a = partitions[k] k += 1 if a: t = a.id G.add_edge(f, t, 0) j += 1 return G.get_connected_components(C, CG) def get_point_in(self): assert (len(self.trim_curves) > 0) b = self.trim_curves[0].bbox() low_t = b.low[1] high_t = b.high[1] for trim_curve in self.trim_curves: b = trim_curve.bbox() if b.low[1] < low_t: low_t = b.low[1] if b.high[1] > high_t: high_t = b.high[1] cut_t = low_t + 3 * (high_t - low_t) / 7 poly_cut = K_RATPOLY(2, 1, cut_t) int_pts = [None] * MAX_NUM_INT_PTS #of Point2Ds num_int_pts = 0 for trim_curve in self.trim_curves: num_int_pts_proto = trim_curve.find_intersections(poly_cut, int_pts_proto, 1) if num_int_pts_proto > 0: j = 0 while j < num_int_pts_proto: assert num_int_pts < MAX_NUM_INT_PTS int_pts[num_int_pts] = int_pts_proto[j] j += 1 assert num_int_pts >= 2 sort_s(int_pts) #, num_int_pts) v_s = (int_pts[0].get_high_s() + int_pts[1].get_low_s()) / 2 return Point2D(v_s, cut_t) def select_relevant_partitions(self, partitions, G, components, colors, color, relevant_partitions, SG): pass class Solid: def __init__(self, patches=[]): self.patches = patches def __repr__(self): return "Solid(" + str(self.patches) + ")" def get_num_patches(self): return len(self.patches) num_patches = property(get_num_patches) def classify_point(self, x, y, z): """ returns IN if (x, y, z) lies inside *this and OUT if (x, y, z) lies outside *this. """ pass classify_pt=classify_point def classify_vector(self, vector): return self.classify_point(vector[0], vector[1], vector[2]) def gen_new_solid(self, partitions1, partitions2, operation): num_new_patches = partitions1.patches.length + partitions2.patches.length new_partitions = [Partition()] * num_new_patches new_partitions1 = [] #[Partition()] * partitions1.patches.length new_partitions2 = [] #[Partition()] * partitions2.patches.length new_patches = [Patch()] * num_new_patches new_partitions1.extend(deepcopy(partitions1)) new_partitions2.extend(deepcopy(partitions2)) self.process_new_partitions(new_partitions1, partitions1, partitions2, False, offset=0) self.process_new_partitions(new_partitions2, partitions2, partitions1, operation.lower()[0] == "d", offset=partitions1.length ) #combine the lists new_partitions = new_partitions1 + new_partitions2 for partition in new_partitions: new_patches.append(Patch(partition)) i = 0 while (i < new_patches.length): new_patch = new_patches[i] j = 0 while j < len(new_partitions[i].trim_curves): new_patch.adj_surfs[j] = new_partitions[i].adj_surfs[j] new_patch.adj_patches[j] = new_patches[new_partitions[i].adj_curves[j]] j += 1 i += 1 result = Solid(new_patches) return result def process_new_partitions(new_partitions, partitions1, partitions2, operation_flag=False, offset=0): #NOTE: you can do i+offset for the correct position if you use a global new_partitions for partition in partitions1: i = partitions1.index(partition) if operation_flag: new_partitions[i].is_head = not new_partitions[i].is_head for trim_curve in partition: j = partition.trim_curves.index(trim_curve) if trim_curve.has_adjacent_partitions: found = False k = 0 while (not found and k < partitions2.length): l = 0 while (not found and l < len(partitions2[k].trim_curves)): if partitions2[k].trim_curves[l] == trim_curve.curve_in_other_dom: found = True l += 1 if not found: k += 1 #assert k < partitions2.length new_partitions[i].adj_surfs[j] = partitions2[k].From.surf new_partitions[i].adj_curves[j] = k + partitions1.length elif not trim_curve.has_adjacent_partitions: found = False k = 0 while not found and k < partitions1.length: if partitions1[k].From == partition.From.adj_patches[partition.adj_curves[j]]: l = 0 while not found and l < len(partitions1[k].trim_curves): if (partitions1[k].trim_curves[l] == partition.From.trim_curves[partition.adj_curves[j]].curve_in_other_dom): found = True l += 1 if not found: k += 1 new_partitions[i].adj_surfs[j] = partition.From.adj_surfs[partition.adj_curves[j]] new_partitions[i].adj_curves[j] = k + offset def boolean(self, other_solid, operation): #intersect overlapping patches if len(self.patches) > 0 and len(other_solid.patches) > 0: for my_patch in self.patches: for other_patch in other_solid.patches: if my_patch.overlap(other_patch): print "Intersecting patch " + str(my_patch) + " and " + str(other_patch) + "." my_patch.intersect(other_patch) print "\tPatch " + str(my_patch) + " on solid 1 has " + str(my_patch.ic) + " intersection curves." print "\tPatch " + str(other_patch) + " on solid 2 has " + str(other_patch.ic) + " intersection curves." #split loops self.split_loops() other_solid.split_loops() #generate partitions t_partitions = self.generate_partitions() s_partitions = other_solid.generate_partitions() #generate graphs if t_partitions.length > 0 and s_partitions.length > 0: (t_partition_graph, t_components, t_component_graph) = self.gen_adjacency(t_partitions) (s_partition_graph, s_components, s_component_graph) = self.gen_adjacency(s_partitions) print "Shooting 3D rays." point = t_partitions[0].get_point_in() v = self.patches[0].surf.param_to_coord(point.get_low) t_classified = other_solid.classify_point(v) point = s_partitions[0].get_point_in() v = other_solid.patches[0].surf.param_to_coord(point.get_low()) s_classified = self.classify_point(v) t_colors = t_component_graph.propagate_color(t_components[0], t_classified) s_colors = s_component_graph.propagate_color(s_components[0], s_classified) print "Selecting partitions of the resulting solid." if operation.lower()[0] == "u": t_new_color, s_new_color = OUT, OUT elif operation.lower()[0] == "i": t_new_color, s_new_color = IN, IN elif operation.lower()[0] == "d": t_new_color, s_new_color = OUT, IN t_new_partitions, t_SG = select_relevant_partitions(t_partitions, t_partition_graph, t_components, t_colors, t_new_color) s_new_partitions, s_SG = select_relevant_partitions(s_partitions, s_partition_graph, s_components, s_colors, s_new_color) print t_new_partitions.length + " partitions are selected from solid 1." print s_new_partitions.length + " partitions are selected from solid 2." #build resulting solid print "Building the resulting solid." if t_partitions.length > 0 and s_partitions.length > 0: result = self.gen_new_solid(t_new_partitions, s_new_partitions, operation) elif t_partitions.length > 0: result = self elif s_partitions.length > 0: result = other_solid print "The resulting solid has " + result.patches.length + " patches." return result def generate_partitions(self): t_partitions = [] for patch in self.patches: partitions_proto = self.gen_partitions(patch) print "Patch " + patch + " forms " + partitions_proto.length + " partitions on solid " + self t_partitions.extend(partitions_proto) return t_partitions def split_loops(self): for patch in self.patches: num_closed_curves = patch.merge_curves() print "Patch " + str(patch) + " on solid 1 has " + str(num_closed_curves) + " loops." if num_closed_curves > 0: pass def weird_repetition(X, x, y, z, n): p = [0, 0] X.coeffs[X.powers_to_index(p)] = x[n] p[0] = 1 p[1] = 0 X.coeffs[X.powers_to_index(p)] = y[n] - x[n] p[0] = 0 p[1] = 1 X.coeffs[X.powers_to_index(p)] = z[n] - x[n] def get_param_plane(x, y, z): assert len(x) == 3 assert len(y) == 3 assert len(z) == 3 d = [1, 1] X = K_RATPOLY(2, d) Y = K_RATPOLY(2, d) Z = K_RATPOLY(2, d) W = K_RATPOLY(2, d) weird_repetition(X, x, y, z, 0) weird_repetition(Y, x, y, z, 1) weird_repetition(Z, x, y, z, 2) W.coeffs[W.powers_to_index([0, 0])] = 1 return (X, Y, Z, W) def weird_repetition2(X, points, n): X.coeffs[X.powers_to_index([0, 0])] = points[0][n] X.coeffs[X.powers_to_index([1, 0])] = points[1][n] - points[0][n] X.coeffs[X.powers_to_index([0, 1])] = points[3][n] - points[0][n] X.coeffs[X.powers_to_index([1, 1])] = points[2][n] + points[0][n] - points[3][n] - points[1][n] def get_param_bilin(points): d = [1, 1] X = K_RATPOLY(2, d) Y = K_RATPOLY(2, d) Z = K_RATPOLY(2, d) W = K_RATPOLY(2, d) weird_repetition2(X, points, 0) weird_repetition2(Y, points, 1) weird_repetition2(Z, points, 2) W.coeffs[W.powers_to_index([0, 0])] = 1 return (X, Y, Z, W) def get_impl_coeffs_bilin(param_coeffs): """ unsigned long get_impl_coeffs_bilin(const bigrational param_coeffs[12], bigrational*& impl_coeffs) computes the coefficients of the implicit formula for the bilinear surface X(s, t) = param_coeffs[0] + param_coeffs[1] s + param_coeffs[2] t + param_coeffs[3] st, Y(s, t) = param_coeffs[4] + param_coeffs[5] s + param_coeffs[6] t + param_coeffs[7] st, Z(s, t) = param_coeffs[8] + param_coeffs[9] s + param_coeffs[10] t + param_coeffs[11] st. returns the number of the coefficients of the implicit formula. """ #X(s, t) = x0 + x1 s + x2 t + x3 st, #Y(s, t) = y0 + y1 s + y2 t + y3 st, #Z(s, t) = z0 + z1 s + z2 t + z3 st, x0 = param_coeffs[0] x1 = param_coeffs[1] x2 = param_coeffs[2] x3 = param_coeffs[3] y0 = param_coeffs[4] y1 = param_coeffs[5] y2 = param_coeffs[6] y3 = param_coeffs[7] z0 = param_coeffs[8] z1 = param_coeffs[9] z2 = param_coeffs[10] z3 = param_coeffs[11] cx01 = x0 * x1 cx02 = x0 * x2 cx03 = x0 * x3 cx12 = x1 * x2 cx13 = x1 * x3 cx23 = x2 * x3 cy01 = y0 * y1 cy02 = y0 * y2 cy03 = y0 * y3 cy12 = y1 * y2 cy13 = y1 * y3 cy23 = y2 * y3 cz01 = z0 * z1 cz02 = z0 * z2 cz03 = z0 * z3 cz12 = z1 * z2 cz13 = z1 * z3 cz23 = z2 * z3 cx00 = x0 * x0 cx11 = x1 * x1 cx22 = x2 * x2 cx33 = x3 * x3 cy00 = y0 * y0 cy11 = y1 * y1 cy22 = y2 * y2 cy33 = y3 * y3 cz00 = z0 * z0 cz11 = z1 * z1 cz22 = z2 * z2 cz33 = z3 * z3 impl_coeffs = [] constant_term = \ cx00 * (cy12 * cz33 - cy13 * cz23 \ + cy33 * cz12 - cy23 * cz13) \ + cx11 * (cy02 * cz23 - cy03 * cz22 \ + cy23 * cz02 - cy22 * cz03) \ + cx22 * (cy01 * cz13 - cy03 * cz11 \ + cy13 * cz01 - cy11 * cz03) \ + cx33 * (cy00 * cz12 - cy01 * cz02 \ + cy12 * cz00 - cy02 * cz01) \ + cx01 * (- cy02 * cz33 + cy03 * cz23 - cy12 * cz23 + cy13 * cz22 \ - cy33 * cz02 + cy23 * cz03 - cy23 * cz12 + cy22 * cz13) \ + cx02 * (- cy01 * cz33 + cy03 * cz13 - cy12 * cz13 + cy23 * cz11 \ - cy33 * cz01 + cy13 * cz03 - cy13 * cz12 + cy11 * cz23) \ + cx13 * (- cy00 * cz23 + cy01 * cz22 - cy02 * cz12 + cy02 * cz03 \ - cy23 * cz00 + cy22 * cz01 - cy12 * cz02 + cy03 * cz02) \ + cx23 * (- cy00 * cz13 + cy02 * cz11 - cy01 * cz12 + cy01 * cz03 \ - cy13 * cz00 + cy11 * cz02 - cy12 * cz01 + cy03 * cz01) \ + 2 * cx03 * (- cy03 * cz12 - cy12 * cz03 + cy12 * cz12) \ + cx03 * (cy01 * cz23 + cy02 * cz13 - cy11 * cz22 \ + cy23 * cz01 + cy13 * cz02 - cy22 * cz11) \ + 2 * cx12 * (- cy03 * cz03 + cy03 * cz12 + cy12 * cz03) \ + cx12 * (cy00 * cz33 - cy01 * cz23 - cy02 * cz13 \ + cy33 * cz00 - cy23 * cz01 - cy13 * cz02) impl_coeffs.append(constant_term) x_coefficient = \ 2 * x0 * (- cy33 * cz12 + cy23 * cz13 \ - cy12 * cz33 + cy13 * cz23) \ + x1 * (cy33 * cz02 + cy23 * cz12 - cy23 * cz03 - cy22 * cz13 \ + cy02 * cz33 + cy12 * cz23 - cy03 * cz23 - cy13 * cz22) \ + x2 * (cy33 * cz01 + cy13 * cz12- cy13 * cz03 - cy11 * cz23 \ + cy01 * cz33 + cy12 * cz13 - cy03 * cz13 - cy23 * cz11) \ + 2 * x3 * (- cy12 * cz12 + cy12 * cz03 + cy03 * cz12) \ + x3 * (- cy23 * cz01 + cy22 * cz11 - cy13 * cz02 \ - cy01 * cz23 + cy11 * cz22 - cy02 * cz13) impl_coeffs.append(x_coefficient) y_coefficient = \ 2 * y0 * (- cx33 * cz12 + cx23 * cz13 \ - cx12 * cz33 + cx13 * cz23) \ + y1 * (cx33 * cz02 + cx23 * cz12 - cx23 * cz03 - cx22 * cz13 \ + cx02 * cz33 + cx12 * cz23 - cx03 * cz23 - cx13 * cz22) \ + y2 * (cx33 * cz01 + cx13 * cz12- cx13 * cz03 - cx11 * cz23 \ + cx01 * cz33 + cx12 * cz13 - cx03 * cz13 - cx23 * cz11) \ + 2 * y3 * (- cx12 * cz12 + cx12 * cz03 + cx03 * cz12) \ + y3 * (- cx23 * cz01 + cx22 * cz11 - cx13 * cz02 \ - cx01 * cz23 + cx11 * cz22 - cx02 * cz13) impl_coeffs.append(y_coefficient) z_coefficient = \ 2 * z0 * (- cx33 * cy12 + cx23 * cy13 \ - cx12 * cy33 + cx13 * cy23) \ + z1 * (cx33 * cy02 + cx23 * cy12 - cx23 * cy03 - cx22 * cy13 \ + cx02 * cy33 + cx12 * cy23 - cx03 * cy23 - cx13 * cy22) \ + z2 * (cx33 * cy01 + cx13 * cy12- cx13 * cy03 - cx11 * cy23 \ + cx01 * cy33 + cx12 * cy13 - cx03 * cy13 - cx23 * cy11) \ + 2 * z3 * (- cx12 * cy12 + cx12 * cy03 + cx03 * cy12) \ + z3 * (- cx23 * cy01 + cx22 * cy11 - cx13 * cy02 \ - cx01 * cy23 + cx11 * cy22 - cx02 * cy13) impl_coeffs.append(z_coefficient) xy_coefficient = \ - 2 * cz12 * x3 * y3 \ + cz13 * (x2 * y3 + x3 * y2) \ + cz23 * (x1 * y3 + x3 * y1) \ - cz33 * (x1 * y2 + x2 * y1) impl_coeffs.append(xy_coefficient) yz_coefficient = \ - 2 * cx12 * y3 * z3 \ + cx13 * (z2 * y3 + z3 * y2) \ + cx23 * (z1 * y3 + z3 * y1) \ - cx33 * (z1 * y2 + z2 * y1) impl_coeffs.append(yz_coefficient) xz_coefficient = \ - 2 * cy12 * x3 * z3 \ + cy13 * (x2 * z3 + x3 * z2) \ + cy23 * (x1 * z3 + x3 * z1) \ - cy33 * (x1 * z2 + x2 * z1) impl_coeffs.append(xz_coefficient) xsquared_coefficient = cy33 * cz12 - cy23 * cz13 + cy12 * cz33 - cy13 * cz23 impl_coeffs.append(xsquared_coefficient) ysquared_coefficient = cx33 * cz12 - cx23 * cz13 + cx12 * cz33 - cx13 * cz23 impl_coeffs.append(ysquared_coefficient) zsquared_coefficient = cx33 * cy12 - cx23 * cy13 + cx12 * cy33 - cx13 * cy23 impl_coeffs.append(zsquared_coefficient) return impl_coeffs def get_plane_coeffs(vector0, vector1, vector2): assert len(vector0) == 3 assert len(vector1) == 3 assert len(vector2) == 3 V1 = [None] * 3 V2 = [None] * 3 i = 0 while i < 3: V1[i] = vector1[i] - vector0[i] V2[i] = vector2[i] - vector0[i] i += 1 p_proto = [] p_proto.append(vector1[1] * vector2[2] - vector1[2] * vector2[1]) p_proto.append(vector1[2] * vector2[0] - vector1[0] * vector2[2]) p_proto.append(vector1[0] * vector2[1] - vector1[1] * vector2[0]) if sgn(p_proto[0]) or sgn(p_proto[1]) or sgn(p_proto[2]): p = [None] * 4 i = 0 while i < 3: p[i] = p_proto[i] i += 1 p[3] = -p[0] * vector0[0] - p[1] * vector0[1] - p[2] * vector0[2] else: #if vector0, vector1 and vector2 are collinear p = [0] * 4 return p def get_impl_plane_bilin(points): plane = get_plane_coeffs(points[0], points[1], points[2]) s = plane[0] * points[3][0] + plane[1] * points[3][1] + plane[2] * points[3][2] + plane[3] if not sgn(s): d = [1, 1, 1] impl = K_RATPOLY(3, d) impl.coeffs[impl.powers_to_index([0, 0, 0])] = plane[3] impl.coeffs[impl.powers_to_index([1, 0, 0])] = plane[0] impl.coeffs[impl.powers_to_index([0, 1, 0])] = plane[1] impl.coeffs[impl.powers_to_index([0, 0, 1])] = plane[2] impl.reduce_deg() impl.reduce_num_coeffs() is_plane = True else: param_coeffs = [ points[0][0], points[1][0] - points[0][0], points[3][0] - points[0][0], None, points[0][1], points[1][1] - points[0][1], points[3][1] - points[0][1], None, points[0][2], points[1][2] - points[0][2], points[3][2] - points[0][2], None ] param_coeffs[3] = points[2][0] - points[3][0] - param_coeffs[1] param_coeffs[7] = points[2][1] - points[3][1] - param_coeffs[5] param_coeffs[11] = points[2][2] - points[3][2] - param_coeffs[9] impl_coeffs = get_impl_coeffs_bilin(param_coeffs) assert len(impl_coeffs) == 10 d = [2, 2, 2] impl = K_RATPOLY(3, d) impl.coeffs[impl.powers_to_index([0, 0, 0])] = impl_coeffs[0] impl.coeffs[impl.powers_to_index([1, 0, 0])] = impl_coeffs[1] impl.coeffs[impl.powers_to_index([0, 1, 0])] = impl_coeffs[2] impl.coeffs[impl.powers_to_index([0, 0, 1])] = impl_coeffs[3] impl.coeffs[impl.powers_to_index([1, 1, 0])] = impl_coeffs[4] impl.coeffs[impl.powers_to_index([0, 1, 1])] = impl_coeffs[5] impl.coeffs[impl.powers_to_index([1, 0, 1])] = impl_coeffs[6] impl.coeffs[impl.powers_to_index([2, 0, 0])] = impl_coeffs[7] impl.coeffs[impl.powers_to_index([0, 2, 0])] = impl_coeffs[8] impl.coeffs[impl.powers_to_index([0, 0, 2])] = impl_coeffs[9] is_plane = False return is_plane, impl def get_patch4(points): is_plane, impl = get_impl_plane_bilin(points) if is_plane: (X, Y, Z, W) = get_param_plane(points[0], points[1], points[3]) else: (X, Y, Z, W) = get_param_bilin(points) surface = Surface(impl, X, Y, Z, W) trim_curves = [Curve()] * 4 X, Y, Z, W = surface.X, surface.Y, surface.Z, surface.W #locals: i, j, is_plane, t, #M = [[None, None, None], [None, None, None], [None, None, None]] M = [ [None] * 3 for temp1 in range(3) ] d = p = v = [None, None] P = K_RATPOLY() start = end = Point2D() segment = K_SEGMENT() num_trim_curves = len(trim_curves) if is_plane: #a planar patch #X p[0] = 1; p[1] = 0; M[0][0] = X.get_coeff(p) print "X.get_coeff([1,0]) == ", str(X.get_coeff([1,0])) print "M[0][0] is: ", str(M[0][0]) p[0] = 0; p[1] = 1; M[0][1] = X.get_coeff(p) p[0] = 0; p[1] = 0; M[0][2] = X.get_coeff(p) - points[2][0] #Y p[0] = 1; p[1] = 0; M[1][0] = Y.get_coeff(p) p[0] = 0; p[1] = 1; M[1][1] = Y.get_coeff(p) p[0] = 0; p[1] = 0; M[1][2] = Y.get_coeff(p) - points[2][1] #Z p[0] = 1; p[1] = 0; M[2][0] = Z.get_coeff(p) p[0] = 0; p[1] = 1; M[2][1] = Z.get_coeff(p) p[0] = 0; p[1] = 0; M[2][2] = Z.get_coeff(p) - points[2][2] if sgn(M[0][0]): i = 0 elif sgn(M[1][0]): i = 1 else: #if sgn(M[2][0]) i = 2 print lineno(), "M is: ", str(M) if i > 0: j = 0 while j < 3: t = M[0][j] M[0][j] = M[i][j] M[i][j] = t j += 1 i = 1 while i < 3: #TODO: check what M is in the esolid SWIG wrapper or C++ program #see genbox.cc line 694 print "M[i][0] is: ", str(M[i][0]) print "M[0][0] is: ", str(M[0][0]) print "M is: ", str(M) t = M[i][0] / M[0][0] #try: # t = M[i][0] / M[0][0] #except Exception, error: # t = 0 j = 0 while j < 3: M[i][j] -= t * M[0][j] j += 1 i += 1 if sgn(M[1][1]): v[1] = - M[1][2] / M[1][1] else: v[1] = - M[2][2] / M[2][1] v[0] = - (M[0][2] + M[0][1] * v[1]) / M[0][0] #set trim curves i = 0 while i < num_trim_curves: if i == 0: #trim_curves[0]: t = 0 from (0, 0) to (1, 0) d[0] = 0 d[1] = 1 P = K_RATPOLY(2, d) p[0] = 0 p[1] = 1 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(0, 0) end = Point2D(1, 0) segment = [K_SEGMENT(start, end)] elif i == 1: #trim_curves[1]: -1 + s + (1 - v[0]) / v[1] t = 0 from (1, 0) to v d[0] = d[1] = 1 P = K_RATPOLY(2, d) p[0] = p[1] = 0 P.coeffs[P.powers_to_index(p)] = -1 p[0] = 1; p[1] = 0 P.coeffs[P.powers_to_index(p)] = 1 p[0] = 0; p[1] = 1 P.coeffs[P.powers_to_index(p)] = (1 - v[0]) / v[1] P.reduce_deg() start = Point2D(1, 0) end = Point2D(v[0], v[1]) segment = [K_SEGMENT(start, end)] elif i == 2: #trim_curves[2]: - 1 + (1 - v[1]) / v[0] s + t from v to (0, 1) d[0] = d[1] = 1 P = K_RATPOLY(2, d) p[0] = p[1] = 0 P.coeffs[P.powers_to_index(p)] = -1 p[0] = 1; p[1] = 0 P.coeffs[P.powers_to_index(p)] = (1 - v[1]) / v[0] p[0] = 0; p[1] = 1 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(v[0], v[1]) end = Point2D(0, 1) segment = [K_SEGMENT(start, end)] else: #if i == 3 #trim_curves[3]: s = 0 from (0, 1) to (0, 0) d[0] = 1; d[1] = 0 P = K_RATPOLY(2, d) p[0] = 1; p[1] = 0 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(0, 1) end = Point2D(0, 0) segment = [K_SEGMENT(start, end)] trim_curves[i] = Curve(P, segment, 1) i += 1 else: #if not is_plane #A bilinear patch. Set trim_curves to be t=0, s=1, t=1 and s=0. i = 0 while i < num_trim_curves: if i == 0: #trim_curves[0]: t = 0 from (0, 0) to (1, 0) d[0] = 0 d[1] = 1 P = K_RATPOLY(2, d) p[0] = 0; p[1] = 1 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(0, 0) end = Point2D(1, 0) segment = [K_SEGMENT(start, end)] elif i == 2: #trim_curves[2]: -1 + t from (1, 1) to (0, 1) d[0] = 0 d[1] = 1 P = K_RATPOLY(2, d) p[0] = p[1] = 0 P.coeffs[P.powers_to_index(p)] = -1 p[0] = 0; p[1] = 1 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(1, 1) end = Point2D(0, 1) segment = [K_SEGMENT(start, end)] else: #if i == 3 --> trim_curves[3]: s = 0 from (0, 1) to (0, 0) d[0] = 1; d[1] = 0 P = K_RATPOLY(2, d) p[0] = 1; p[1] = 0 P.coeffs[P.powers_to_index(p)] = 1 P.reduce_deg() start = Point2D(0, 1) end = Point2D(0, 0) segment = [K_SEGMENT(start, end)] trim_curves[i] = Curve(P, segment, 1) i += 1 patch = Patch(surface, trim_curves) return patch def set_adj_surfs_and_adj_patches(patches, n, a, b, c, d): patches[n].adj_surfs = [] patches[n].adj_surfs.append(patches[a].surf) patches[n].adj_surfs.append(patches[b].surf) patches[n].adj_surfs.append(patches[c].surf) patches[n].adj_surfs.append(patches[d].surf) patches[n].adj_patches = [] patches[n].adj_patches.append(patches[a]) patches[n].adj_patches.append(patches[b]) patches[n].adj_patches.append(patches[c]) patches[n].adj_patches.append(patches[d]) def gen_box(points): """returns a box whose corners are the given points""" assert len(points) == 8 patches = [] # # 7-------------6 # /| /| # / | / | # / | / | # 4-------------5 | # | | | | # | | | | # | | | | # | 0---------|---1 # | / | / # | / | / # |/ |/ # 3-------------2 # # 1. Set pathces. # patches[0]: (0 1 2 3) patch = get_patch4([points[0], points[1], points[2], points[3]]) patches.append(patch) # patches[1]: (1 0 7 6) patch = get_patch4([points[1], points[0], points[7], points[6]]) patches.append(patch) # patches[2]: (2 1 6 5) patch = get_patch4([points[2], points[1], points[6], points[5]]) patches.append(patch) # patches[3]: (3 2 5 4) patch = get_patch4([points[3], points[2], points[5], points[4]]) patches.append(patch) # patches[4]: (0 3 4 7) patch = get_patch4([points[0], points[3], points[4], points[7]]) patches.append(patch) # patches[5]: (4 5 6 7) patch = get_patch4([points[4], points[5], points[6], points[7]]) patches.append(patch) # 2. For each patch, set its adj_surfs and adj_patches. set_adj_surfs_and_adj_patches(patches, 0, 1, 2, 3, 4) set_adj_surfs_and_adj_patches(patches, 1, 0, 4, 5, 2) set_adj_surfs_and_adj_patches(patches, 2, 0, 1, 5, 3) set_adj_surfs_and_adj_patches(patches, 3, 0, 2, 5, 4) set_adj_surfs_and_adj_patches(patches, 4, 0, 3, 5, 1) set_adj_surfs_and_adj_patches(patches, 5, 3, 2, 1, 4) # 3. Associate edges. patches[0].trim_curves[0].assoc(patches[1].trim_curves[0], -1) patches[0].trim_curves[1].assoc(patches[2].trim_curves[0], -1) patches[0].trim_curves[2].assoc(patches[3].trim_curves[0], -1) patches[0].trim_curves[3].assoc(patches[4].trim_curves[0], -1) patches[5].trim_curves[0].assoc(patches[3].trim_curves[2], -1) patches[5].trim_curves[1].assoc(patches[2].trim_curves[2], -1) patches[5].trim_curves[2].assoc(patches[1].trim_curves[2], -1) patches[5].trim_curves[3].assoc(patches[4].trim_curves[2], -1) patches[1].trim_curves[3].assoc(patches[2].trim_curves[1], -1) patches[2].trim_curves[3].assoc(patches[3].trim_curves[1], -1) patches[3].trim_curves[3].assoc(patches[4].trim_curves[1], -1) patches[4].trim_curves[3].assoc(patches[1].trim_curves[1], -1) #make some tasty cake for patch in patches: patch.surf.X.reduce_deg() patch.surf.Y.reduce_deg() patch.surf.Z.reduce_deg() patch.surf.W.reduce_deg() solid = Solid(patches) return solid