You are given two planar quadraticBézier curves, having control points and, respectively. They are parabolas intersecting at four points, as shown in the following diagram:
The task is to write a program that finds all four intersection points and prints their coordinates. You may use any algorithm you know of or can think of, including any of those that others have used.
This is a variant of the Bézier clipping method ofThomas W. Sederberg's lecture notes on "Computer Aided Geometric Design", using symmetric power (s-power) polynomials instead of the usual "control points". I do not know that it is any better than the method of the s-power polynomial flattening, which is done on this page in several languages. The flattening method is considerably simpler. Also, I have doubts about the numerical stability of the clipping method.
On the other hand, the code below introduces a new way of doing Bézier clipping, by transformation of coordinates. And it also hints at a geometric interpretation of s-power curves, in which all the constituent "points" undergo rotations and scalings, but only the endpoints are subject to translations.
Ada is an Algol-like language that is case-insensitive. I have used an all-lowercase style with manual indentation. Subprogram names are overloaded. Subprogram parameters may be passed either positionally or by name. Ada does have a few oddities, among which are that the term for a pointer is not "pointer" but "access".
pragmawide_character_encoding(utf8);pragmaada_2022;pragmaassertion_policy(check);---- The algorithm here is iterated "Bézier clipping" [1] of quadratic-- Bézier curves in a plane. However, the curves will be represented-- in the s-power basis [2,3], and there will be these other-- differences from the algorithms suggested by [1]:---- * Before clipping, the problem will be moved to a coordinate system-- in which the "skeleton" of the clipper curve—the line segment-- between its endpoints—extends from (0,0) and (1,0). In this-- coordinate system, the "distance function" is the y-coordinate.---- * The fat line will be the best possible parallel to the "skeleton"-- of the clipper curve. This is easily done in the s-power basis,-- once one has transformed coordinate systems as described above:-- the center coefficient of the (y,t)-polynomial is exactly four-- times the optimum signed fat line width.---- * Intersections with fat line boundaries will be found not by-- extending the sides of control polygons, but instead by the-- quadratic formula.---- * Mechanisms to ensure termination of the algorithm, and to improve-- the chances of finding intersections, may differ.---- Inputs and outputs to the algorithms will be in single-- precision. However, the internals will be in double precision. This-- decision helps resolve some minor numerical quandaries.---- References:---- [1] Sederberg, Thomas W., "Computer Aided Geometric Design"-- (2012). Faculty-- Publications. 1. https://scholarsarchive.byu.edu/facpub/1-- http://hdl.lib.byu.edu/1877/2822---- [2] J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial-- power basis’, ACM Transactions on Graphics, vol 16 no 3, July-- 1997, page 319.---- [3] J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis-- in geometry processing’, ACM Transactions on Graphics, vol 19-- no 1, January 2000, page 35.---- Suggestions for the future:---- * Combination of Bézier clipping and flatness testing. Remember-- from other Rosetta Code examples, that in the s-power basis it is-- simple to test the "flatness" of a parametric plane-- curve. Suppose two curves are running nearly parallel to each-- other, so they are difficult to reduce by clipping. Nevertheless,-- they can be made flatter and flatter, so eventually they become-- "segment-like" and can be tested by a line-intersection routine.--withada.text_io;useada.text_io;withada.float_text_io;useada.float_text_io;withada.numerics.elementary_functions;useada.numerics.elementary_functions;withada.numerics.long_elementary_functions;useada.numerics.long_elementary_functions;withada.unchecked_deallocation;procedurebezier_intersectionsissubtyperealisfloat;subtypelrealislong_float;functionbetween_0_and_1(x:real)returnbooleanisbeginreturn0.0<=xandx<=1.0;endbetween_0_and_1;functionbetween_0_and_1(x:lreal)returnbooleanisbeginreturn0.0<=xandx<=1.0;endbetween_0_and_1;procedurequadratic_formula(a,b,c:lreal;has_smaller:outboolean;smaller:outlreal;has_larger:outboolean;larger:outlreal)withpost=>has_smallerornothas_largeris---- This routine is based on QUADPL from the Naval Surface Warfare-- Center mathematical library nswc.f90 https://archive.ph/2IpWb---- We are interested only in real roots. Be aware that, if there-- are two roots, they might be a double root.--bh,e,d:lreal;-- Avogadro’s number, by definition (since 2019):arbitrary_value:constantlreal:=6.012214076e23;beginhas_smaller:=false;has_larger:=false;smaller:=arbitrary_value;larger:=arbitrary_value;ifa=0.0thenifb=0.0then-- The equation is a constant equation. I do not know what to-- do with that. There may be no solution, or there may be-- infinitely many solutions. Let us treat it as having no-- USEFUL solutions.null;else-- The equation is linear. There truly is only one root.has_smaller:=true;smaller:=-c/b;endif;elsifc=0.0then-- The equation is quadratic, but the constant term-- vanishes. One of the roots is trivially zero.has_smaller:=true;has_larger:=true;smaller:=0.0;larger:=-b/a;else-- Compute the discriminant, avoiding overflow.bh:=b/2.0;ifabs(bh)>=abs(c)thene:=1.0-(a/bh)*(c/bh);d:=sqrt(abs(e))*abs(bh);elsee:=(ifc<0.0then-aelsea);e:=bh*(bh/abs(c))-e;d:=sqrt(abs(e))*sqrt(abs(c));endif;ife<0.0then-- The solutions are complex conjugates, but we are interested-- only in real solutions.null;elsehas_smaller:=true;has_larger:=true;ifbh>=0.0thend:=-d;endif;larger:=(-bh+d)/a;iflarger=0.0thensmaller:=0.0;-- Obviously a double root.elsesmaller:=(c/larger)/a;endif;endif;endif;endquadratic_formula;typepointisrecordx,y:real;end record;typecontrol_point_curveisrecordpt0,pt1,pt2:point;end record;typespower_polyisrecordc0,c1,c2:real;end record;typespower_curveisrecordx,y:spower_poly;end record;functioneval(poly:spower_poly;t:real)returnrealisc0,c1,c2:real;beginc0:=poly.c0;c1:=poly.c1;c2:=poly.c2;return(ift<=0.5thenc0+t*((c2-c0)+((1.0-t)*c1))elsec2+(1.0-t)*((c0-c2)+(t*c1)));endeval;functioneval(curve:spower_curve;t:real)returnpointisbeginreturn(x=>eval(curve.x,t),y=>eval(curve.y,t));endeval;functionto_spower_curve(ctrl:control_point_curve)returnspower_curveiscurve:spower_curve;p0x,p0y:long_float;p1x,p1y:long_float;p2x,p2y:long_float;begincurve.x.c0:=ctrl.pt0.x;curve.y.c0:=ctrl.pt0.y;curve.x.c2:=ctrl.pt2.x;curve.y.c2:=ctrl.pt2.y;p0x:=long_float(ctrl.pt0.x);p0y:=long_float(ctrl.pt0.y);p1x:=long_float(ctrl.pt1.x);p1y:=long_float(ctrl.pt1.y);p2x:=long_float(ctrl.pt2.x);p2y:=long_float(ctrl.pt2.y);curve.x.c1:=float(p1x+p1x-p0x-p2x);curve.y.c1:=float(p1y+p1y-p0y-p2y);returncurve;endto_spower_curve;typeopointisrecord-- A point relative to the origin.x,y:lreal;end record;typedpointisrecord-- A point relative to some other point.dx,dy:lreal;end record;function"-"(pt1,pt2:opoint)returndpointisbeginreturn(dx=>pt1.x-pt2.x,dy=>pt1.y-pt2.y);end"-";function"+"(opt:opoint;dpt:dpoint)returnopointisbeginreturn(x=>opt.x+dpt.dx,y=>opt.y+dpt.dy);end"+";function"+"(pt1,pt2:dpoint)returndpointisbeginreturn(dx=>pt1.dx+pt2.dx,dy=>pt1.dy+pt2.dy);end"+";function"*"(scalar:lreal;pt:dpoint)returndpointisbeginreturn(dx=>scalar*pt.dx,dy=>scalar*pt.dy);end"*";typegeocurveisrecord---- A plane curve in s-power basis, represented geometrically.---- pt0 and pt2 form the curve’s ‘skeleton’. Their convex combination,-- ((1.0-t)*pt0)+(t*pt1), form the linear backbone of the curve.---- (1.0-t)*t*pt1 forms (x,t) and (y,t) parabolas that are-- convex-up/down, and which (because this is an s-power curve)-- vanish as the curve becomes linear.---- The actual geocurve is the locus of (1.0-t)*t*pt1 relative to-- the skeleton, respectively for each value of t in [0,1].---- If [t0,t1] is not [0.0,1.0], then this is a portion of some-- parent curve that is determined by context.---- (I feel as if all this can be crafted nicely in terms of a-- geometric algebra, with special points not only at infinity-- but also at the origin, but let us not worry about that for-- now. Likely the calculations would not change much, but the-- names of things might.)--pt0:opoint;pt1:dpoint;pt2:opoint;t0,t1:lreal;end record;functionto_geocurve(scurve:spower_curve)returngeocurveisgcurve:geocurve;begingcurve.pt0.x:=lreal(scurve.x.c0);gcurve.pt0.y:=lreal(scurve.y.c0);gcurve.pt1.dx:=lreal(scurve.x.c1);gcurve.pt1.dy:=lreal(scurve.y.c1);gcurve.pt2.x:=lreal(scurve.x.c2);gcurve.pt2.y:=lreal(scurve.y.c2);gcurve.t0:=0.0;gcurve.t1:=1.0;returngcurve;endto_geocurve;functioneval(gcurve:geocurve;t:lreal)returnopointispt0:opoint;pt1:dpoint;pt2:opoint;beginpt0:=gcurve.pt0;pt1:=gcurve.pt1;pt2:=gcurve.pt2;return(ift<=0.5thenpt0+t*((pt2-pt0)+((1.0-t)*pt1))elsept2+(1.0-t)*((pt0-pt2)+(t*pt1)));endeval;typeaffine_projectionisrecord---- An affine transformation used to project objects into the-- coordinate system of what we are calling a "skeleton".---- If it is an opoint, first add transpose[dx,dy]. If it is a-- dpoint, ignore dx,dy. Then, in either case, premultiply by-- [[a,b],[c,d]]. (One can view dx,dy as a dpoint.)--dx,dy:lreal;a,b:lreal;c,d:lreal;end record;functionproject(pt:opoint;trans:affine_projection)returnopointisx,y:lreal;begin-- Apply an affine transformation to an opoint, projecting it from-- one coordinate system to another.x:=pt.x+trans.dx;-- The origin may have moved.y:=pt.y+trans.dy;return(x=>trans.a*x+trans.b*y,-- Rotation, scaling, etc.y=>trans.c*x+trans.d*y);endproject;functionproject(pt:dpoint;trans:affine_projection)returndpointisx,y:lreal;begin-- Apply an affine transformation to a dpoint, projecting it from-- one coordinate system to another.x:=pt.dx;-- With dpoint, translation is not performed.y:=pt.dy;return(dx=>trans.a*x+trans.b*y,-- Rotation, scaling, &c.dy=>trans.c*x+trans.d*y);endproject;functionproject(gcurve:geocurve;trans:affine_projection)returngeocurveisbegin-- Apply an affine transformation to a geocurve, projecting it-- from one coordinate system to another.return(pt0=>project(gcurve.pt0,trans),pt1=>project(gcurve.pt1,trans),pt2=>project(gcurve.pt2,trans),t0=>0.0,t1=>1.0);endproject;functionsegment_projection_to_0_1(pt0,pt1:opoint)returnaffine_projectionisa,b:lreal;denom:lreal;begin-- Return the transformation that projects pt0 to (0,0) and pt1 to-- (1,0) without any distortions except translation, rotation, and-- scaling.a:=pt1.x-pt0.x;b:=pt1.y-pt0.y;denom:=(a*a)+(b*b);pragmaassert(denom/=0.0);a:=a/denom;b:=b/denom;return(dx=>-pt0.x,dy=>-pt0.y,a=>a,b=>b,c=>-b,d=>a);endsegment_projection_to_0_1;proceduretest_segment_projection_to_0_1istrans:affine_projection;pt:opoint;begin-- Some unit testing.trans:=segment_projection_to_0_1((1.0,2.0),(-3.0,5.0));pt:=project((1.0,2.0),trans);pragmaassert(abs(pt.x)<1.0e-10);pragmaassert(abs(pt.y)<1.0e-10);pt:=project((-3.0,5.0),trans);pragmaassert(abs(pt.x-1.0)<1.0e-10);pragmaassert(abs(pt.y)<1.0e-10);trans:=segment_projection_to_0_1((0.0,2.0),(0.0,-50.0));pt:=project((0.0,2.0),trans);pragmaassert(abs(pt.x)<1.0e-10);pragmaassert(abs(pt.y)<1.0e-10);pt:=project((0.0,-50.0),trans);pragmaassert(abs(pt.x-1.0)<1.0e-10);pragmaassert(abs(pt.y)<1.0e-10);endtest_segment_projection_to_0_1;functionskeleton_projection(clipper:geocurve)returnaffine_projectionisbegin-- Return the transformation that projects the "skeleton" of the-- clipper into a coordinate system where the skeleton goes from-- the origin (0,0) to the point (1,0).returnsegment_projection_to_0_1(clipper.pt0,clipper.pt2);endskeleton_projection;procedurefat_line_bounds(projected_clipper:geocurve;dmin,dmax:outlreal)isd:lreal;d1,d2:lreal;padding:constantlreal:=lreal(real'model_epsilon)/1000.0;begin---- (1-t)*t is bounded on [0,1] by 0 and 1/4. Thus we have our fat-- line boundaries. We leave some padding, however, for numerical-- safety. The padding is much smaller than our tolerance. (We are-- calculating in double precision, but to obtain results in-- single precision.)---- The fat line bounds will be y=dmin and y=dmax. If you ignore-- the padding, the fat line is the best possible among those-- parallel to the skeleton. It is likely to be much tighter than-- what you would infer from the Bernstein-basis control polygon.--d:=projected_clipper.pt1.dy/4.0;d1:=-padding*d;d2:=(1.0+padding)*d;dmin:=lreal'min(d1,d2);dmax:=lreal'max(d1,d2);endfat_line_bounds;typeclipping_placesis(clip_everywhere,clip_nowhere,clip_left_and_right);typeclipping_planisrecordwhere:clipping_places;t_lft:lreal;t_rgt:lreal;end record;functionmake_clipping_plan(projected_clippehend:geocurve;dmin,dmax:lreal)returnclipping_planisncrossings:integerrange0..4;tcrossings:array(1..4)oflreal;dcrossings:array(1..4)oflreal;has_smaller:boolean;has_larger:boolean;smaller:lreal;larger:lreal;c0,c1,c2:lreal;functionoutside_of_boundsreturnbooleanisdy,ymin,ymax:lreal;begin-- (1-t)*t is bounded on [0,1] by 0 and 1/4. Use that fact to-- make padding for a simple (although not especially tight)-- bounds on y.dy:=abs(projected_clippehend.pt1.dy)/4.0;ymin:=-dy+lreal'min(projected_clippehend.pt0.y,projected_clippehend.pt2.y);ymax:=dy+lreal'max(projected_clippehend.pt0.y,projected_clippehend.pt2.y);return(ymax<=dminordmax<=ymin);endoutside_of_bounds;procedurefind_crossings(yvalue:lreal)isbeginquadratic_formula(-c1,c1+c2-c0,c0-yvalue,has_smaller,smaller,has_larger,larger);ifhas_smallerthenifbetween_0_and_1(smaller)thenncrossings:=ncrossings+1;tcrossings(ncrossings):=smaller;dcrossings(ncrossings):=yvalue;endif;endif;ifhas_largerthenifbetween_0_and_1(larger)thenncrossings:=ncrossings+1;tcrossings(ncrossings):=larger;dcrossings(ncrossings):=yvalue;endif;endif;endfind_crossings;functioncreate_a_clipping_planreturnclipping_planisplan:clipping_plan;begin---- For simplicity, there is just one situation in which we will-- clip: when there is one crossing of dmin and one crossing of-- dmax.--ifncrossings=2and thendcrossings(1)/=dcrossings(2)thenplan.where:=clip_left_and_right;plan.t_lft:=lreal'min(tcrossings(1),tcrossings(2));plan.t_rgt:=lreal'max(tcrossings(1),tcrossings(2));elseplan.where:=clip_nowhere;-- Arbitrary values:plan.t_lft:=0.0;plan.t_rgt:=0.0;endif;returnplan;endcreate_a_clipping_plan;functionanalyze_and_planreturnclipping_planisbeginc0:=projected_clippehend.pt0.y;c1:=projected_clippehend.pt1.dy;c2:=projected_clippehend.pt2.y;ncrossings:=0;find_crossings(dmin);find_crossings(dmax);returncreate_a_clipping_plan;endanalyze_and_plan;beginreturn(ifoutside_of_boundsthen(where=>clip_everywhere,-- Arbitrary values:t_lft=>0.0,t_rgt=>0.0)elseanalyze_and_plan);endmake_clipping_plan;functionmake_portion(gcurve:geocurve;t0,t1:lreal)returngeocurveispt0:opoint;pt1:dpoint;pt2:opoint;beginpt0:=eval(gcurve,t0);pt1:=((t1-t0-t0)*t1+(t0*t0))*gcurve.pt1;pt2:=eval(gcurve,t1);return(pt0=>pt0,pt1=>pt1,pt2=>pt2,t0=>t0,t1=>t1);endmake_portion;proceduresegment_parameters(a0,a1:opoint;b0,b1:opoint;ta,tb:outlreal)isaxdiff:lreal;aydiff:lreal;bxdiff:lreal;bydiff:lreal;denom:lreal;ta1:lreal;tb1:lreal;begin-- Do the line segments (a0,a1) and (b0,b1) intersect? If so,-- return t-parameter values in ta,tb in [0,1] for the point of-- intersection, treating them as parametric splines of degree-- 1. Otherwise return ta=-1.0 and tb=-1.0.ta:=-1.0;tb:=-1.0;axdiff:=a1.x-a0.x;aydiff:=a1.y-a0.y;bxdiff:=b1.x-b0.x;bydiff:=b1.y-b0.y;denom:=(axdiff*bydiff)-(aydiff*bxdiff);ta1:=((bxdiff*a0.y)-(bydiff*a0.x)+(b0.x*b1.y)-(b1.x*b0.y));ta1:=ta1/denom;ifbetween_0_and_1(ta1)thentb1:=-((axdiff*b0.y)-(aydiff*b0.x)+(a0.x*a1.y)-(a1.x*a0.y));tb1:=tb1/denom;ifbetween_0_and_1(tb1)thenta:=ta1;tb:=tb1;endif;endif;endsegment_parameters;typet_pair;typet_pair_accessisaccesst_pair;typet_pairisrecordtp,tq:real;-- Single precision results.next:t_pair_access;end record;proceduret_pair_deleteisnewada.unchecked_deallocation(t_pair,t_pair_access);functionlength(params:t_pair_access)returnnaturalisn:natural;p:t_pair_access;beginn:=0;p:=params;whilep/=nullloopn:=n+1;p:=p.next;endloop;returnn;endlength;procedureinsert_t_pair(params:inoutt_pair_access;tp,tq:real)withpre=>between_0_and_1(tp)andbetween_0_and_1(tq)ispair:t_pair_access;begin-- Test for duplicates and also sort as we insert. This is a-- recursive implementation.ifparams=nullthenpair:=newt_pair;params:=pair;pair.next:=null;pair.tp:=tp;pair.tq:=tq;elsifabs(params.tp-tp)<=real'model_epsilonthennull;-- A duplicate, to within epsilon.elsifparams.tp<tptheninsert_t_pair(params.next,tp,tq);elsepair:=newt_pair;pair.next:=params.next;params.next:=pair;pair.tp:=params.tp;pair.tq:=params.tq;params.tp:=tp;params.tq:=tq;endif;endinsert_t_pair;procedureinsert_t_pair(params:inoutt_pair_access;tp,tq:lreal)withpre=>between_0_and_1(tp)andbetween_0_and_1(tq)isbegininsert_t_pair(params,real(tp),real(tq));endinsert_t_pair;typeclipping_rolesis(pportion_clips_qportion,qportion_clips_pportion);typeintersection_job;typeintersection_job_accessisaccessintersection_job;typeintersection_jobisrecordpportion:geocurve;qportion:geocurve;roles:clipping_roles;next:intersection_job_access;end record;procedureintersection_job_deleteisnewada.unchecked_deallocation(intersection_job,intersection_job_access);functionfind_intersections(p,q:geocurve)returnt_pair_accessisparams:t_pair_access;workload:intersection_job_access;pportion:geocurve;qportion:geocurve;roles:clipping_roles;tp,tq:lreal;tolerance:constantlreal:=0.4*lreal(real'model_epsilon);bisection_threshold:constantlreal:=0.2;procedureinsert_job(pportion:geocurve;qportion:geocurve;roles:clipping_roles)withpost=>workload/=nullisjob:intersection_job_access;beginjob:=newintersection_job;job.pportion:=pportion;job.qportion:=qportion;job.roles:=roles;job.next:=workload;workload:=job;endinsert_job;procedureinsert_job(pportion:geocurve;qportion:geocurve)withpost=>workload/=nullisroles:clipping_roles;beginifpportion.t1-pportion.t0<=qportion.t1-qportion.t0thenroles:=pportion_clips_qportion;elseroles:=qportion_clips_pportion;endif;insert_job(pportion,qportion,roles);endinsert_job;procedurechoose_jobwithpre=>workload/=nullisjob:intersection_job_access;beginjob:=workload;workload:=workload.next;pportion:=job.pportion;qportion:=job.qportion;roles:=job.roles;intersection_job_delete(job);endchoose_job;functiont_parameters_are_within_tolerancereturnbooleanisbeginreturn(pportion.t1-pportion.t0<=toleranceandqportion.t1-qportion.t0<=tolerance);endt_parameters_are_within_tolerance;procedureinsert_intersection_parametersisbeginsegment_parameters(pportion.pt0,pportion.pt2,qportion.pt0,qportion.pt2,tp,tq);ifbetween_0_and_1(tp)andbetween_0_and_1(tq)thentp:=(1.0-tp)*pportion.t0+tp*pportion.t1;tq:=(1.0-tq)*qportion.t0+tq*qportion.t1;insert_t_pair(params,tp,tq);endif;endinsert_intersection_parameters;procedureno_intersectionsisbegin-- Do nothing. The reason to have a procedure for this is so the-- code will document itself.null;endno_intersections;procedurebisect_themispportion1:geocurve;pportion2:geocurve;qportion1:geocurve;qportion2:geocurve;t0,t1,th:lreal;begin---- Sederberg employs heuristics to decide what to bisect and who-- should be the clipper. We will simply bisect both curves, and-- use some heuristic coded elsewhere to choose the-- clipper.---- (Bisection as it is normally done in the s-power basis is-- equivalent to change of homogeneous coordinates. It can be-- applied to p and q directly and, for quadratics, is a simple-- operation.)--t0:=pportion.t0;t1:=pportion.t1;th:=0.5*t0+0.5*t1;pportion1:=make_portion(p,t0,th);pportion2:=make_portion(p,th,t1);t0:=qportion.t0;t1:=qportion.t1;th:=0.5*t0+0.5*t1;qportion1:=make_portion(q,t0,th);qportion2:=make_portion(q,th,t1);insert_job(pportion1,qportion1);insert_job(pportion1,qportion2);insert_job(pportion2,qportion1);insert_job(pportion2,qportion2);endbisect_them;procedureclip_and_maybe_bisect(plan:clipping_plan)isclipping:geocurve;t_lft:lreal;t_rgt:lreal;t0,t1:lreal;begint_lft:=plan.t_lft;t_rgt:=plan.t_rgt;ifroles=pportion_clips_qportionthent0:=(1.0-t_lft)*qportion.t0+t_lft*qportion.t1;t1:=(1.0-t_rgt)*qportion.t0+t_rgt*qportion.t1;clipping:=make_portion(q,t0,t1);if1.0-(t_rgt-t_lft)<bisection_thresholdthenqportion:=clipping;bisect_them;elseinsert_job(pportion,clipping,qportion_clips_pportion);endif;elset0:=(1.0-t_lft)*pportion.t0+t_lft*pportion.t1;t1:=(1.0-t_rgt)*pportion.t0+t_rgt*pportion.t1;clipping:=make_portion(p,t0,t1);if1.0-(t_rgt-t_lft)<bisection_thresholdthenpportion:=clipping;bisect_them;elseinsert_job(clipping,qportion,pportion_clips_qportion);endif;endif;endclip_and_maybe_bisect;procedureperform_clippingisclippehend,clipper:geocurve;projected_clippehend:geocurve;projected_clipper:geocurve;trans:affine_projection;dmin,dmax:lreal;plan:clipping_plan;beginifroles=pportion_clips_qportionthenclipper:=pportion;clippehend:=qportion;elseclipper:=qportion;clippehend:=pportion;endif;trans:=skeleton_projection(clipper);projected_clipper:=project(clipper,trans);fat_line_bounds(projected_clipper,dmin,dmax);projected_clippehend:=project(clippehend,trans);plan:=make_clipping_plan(projected_clippehend,dmin,dmax);caseplan.whereiswhenclip_everywhere=>no_intersections;whenclip_nowhere=>bisect_them;whenclip_left_and_right=>clip_and_maybe_bisect(plan);endcase;endperform_clipping;beginparams:=null;workload:=null;insert_job(p,q);whileworkload/=nullloopchoose_job;ift_parameters_are_within_tolerancetheninsert_intersection_parameters;elseperform_clipping;endif;endloop;returnparams;endfind_intersections;functionfind_intersections(p,q:spower_curve)returnt_pair_accessisbeginreturnfind_intersections(to_geocurve(p),to_geocurve(q));endfind_intersections;pctrl,qctrl:control_point_curve;p,q:spower_curve;params,ptr:t_pair_access;begintest_segment_projection_to_0_1;pctrl:=((x=>-1.0,y=>0.0),(x=>0.0,y=>10.0),(x=>1.0,y=>0.0));qctrl:=((x=>2.0,y=>1.0),(x=>-8.0,y=>2.0),(x=>2.0,y=>3.0));p:=to_spower_curve(pctrl);q:=to_spower_curve(qctrl);params:=find_intersections(p,q);set_col(to=>9);put("convex up");set_col(to=>44);put_line("convex left");ptr:=params;whileptr/=nullloopset_col(to=>3);put(ptr.tp,fore=>1,aft=>6,exp=>0);set_col(to=>13);put("(");put(eval(p,ptr.tp).x,fore=>1,aft=>6,exp=>0);put(", ");put(eval(p,ptr.tp).y,fore=>1,aft=>6,exp=>0);put(")");set_col(to=>38);put(ptr.tq,fore=>1,aft=>6,exp=>0);set_col(to=>48);put("(");put(eval(q,ptr.tq).x,fore=>1,aft=>6,exp=>0);put(", ");put(eval(q,ptr.tq).y,fore=>1,aft=>6,exp=>0);put_line(")");ptr:=ptr.next;endloop;endbezier_intersections;
convex up convex left 0.072508 (-0.854983, 1.345017) 0.172508 (-0.854983, 1.345017) 0.159488 (-0.681025, 2.681026) 0.840512 (-0.681026, 2.681025) 0.827492 (0.654984, 2.854983) 0.927492 (0.654983, 2.854983) 0.940512 (0.881025, 1.118975) 0.059488 (0.881025, 1.118975)
This program flattens one of the curves (that is, converts it to a piecewise linear approximation) and finds intersections between the line segments and the other curve. This requires solving many quadratic equations, but that can be done by the quadratic formula.
(I do the flattening in part by using a representation of Bézier curves that probably is not widely known. For quadratic splines, the representation amounts to a line plus a quadratic term. When the quadratic term gets small enough, I simply remove it. Here, though, is an interesting blog post that talks about several other methods:https://raphlinus.github.io/graphics/curves/2019/12/23/flatten-quadbez.html )
(* In this program, one of the two curves is "flattened" (converted to a piecewise linear approximation). Then the problem is reduced to finding intersections of the other curve with line segments. I have never seen this method published in the literature, but somewhere saw it hinted at. Mainly to increase awareness of the representation, I flatten the one curve using the symmetric power polynomial basis. See J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, page 319. J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis in geometry processing’, ACM Transactions on Graphics, vol 19 no 1, January 2000, page 35. *)#include "share/atspre_staload.hats"%{^#include <math.h>%}(* One simple way to make a foreign function call. I want to use only the ATS prelude, but the prelude does not include support for the C math library. (The bundled libats/libc does, and separately available ats2-xprelude does.) *)extern fn sqrt : double -<> double = "mac#sqrt"macdef huge_val = $extval (double, "HUGE_VAL")#define NIL list_nil ()#define :: list_consfun eval_bernstein_degree2 (@(q0 : double, q1 : double, q2 : double), t : double) : double = let (* The de Casteljau algorithm. (The Schumaker-Volk algorithm also is good BTW and is faster. In this program it should make no noticeable difference, however.) *) val s = 1.0 - t val q01 = (s * q0) + (t * q1) val q12 = (s * q1) + (t * q2) val q012 = (s * q01) + (t * q12) in q012 end(* @(...) means an unboxed tuple. Also often can be written without the @, but then might be mistaken for argument parentheses. *)funbernstein2spower_degree2 (@(c0 : double, c1 : double, c2 : double)) : @(double, double, double) = (* Convert from Bernstein coefficients (control points) to symmetric power coefficients. *) @(c0, c1 + c1 - c0 - c2, c2)funspower_portion_degree2 (@(c0 : double, c1 : double, c2 : double), @(t0 : double, t1 : double)) : @(double, double, double) = (* Compose spower(c0, c1, c2) with spower(t0, t1). This will map the portion [t0,t1] onto [0,1]. (I got these expressions with Maxima, a while back.) *) let val t0_t0 = t0 * t0 and t0_t1 = t0 * t1 and t1_t1 = t1 * t1 and c2p1m0 = c2 + c1 - c0 val d0 = c0 + (c2p1m0 * t0) - (c1 * t0_t0) and d1 = (c1 * t1_t1) - ((c1 + c1) * t0_t1) + (c1 * t0_t0) and d2 = c0 + (c2p1m0 * t1) - (c1 * t1_t1) in @(d0, d1, d2) endfunsolve_linear_quadratic (@(px0 : double, px1 : double), @(py0 : double, py1 : double), @(qx0 : double, qx1 : double, qx2 : double), @(qy0 : double, qy1 : double, qy2 : double)) (* Returns the two real roots, or any numbers outside [0,1], if there are no real roots. *) : @(double, double) = let (* The coefficients of the quadratic equation can be found by the following Maxima commands, which implicitize the line segment and plug in the parametric equations of the parabola: /* The line. */ xp(t) := px0*(1-t) + px1*t$ yp(t) := py0*(1-t) + py1*t$ /* The quadratic (Bernstein basis). */ xq(t) := qx0*(1-t)**2 + 2*qx1*t*(1-t) + qx2*t**2$ yq(t) := qy0*(1-t)**2 + 2*qy1*t*(1-t) + qy2*t**2$ /* Implicitize and plug in. */ impl(t) := resultant(xq(t)-xp(tau), yq(t)-yp(tau), tau)$ impl(t); expand(impl(t)); Consequently you get a quadratic equation in t, which can be solved by the quadratic formula. Sometimes people solve this problem by projecting the line segment onto the x- or y-axis, and similarly projecting the parabola. However, the following is simpler to write, if you have Maxima to derive it for you. Whether it is better to use the expanded expression (as here) or not to, I do not know. *) val px0py1 = px0 * py1 and px1py0 = px1 * py0 and px0qy0 = px0 * qy0 and px0qy1 = px0 * qy1 and px0qy2 = px0 * qy2 and px1qy0 = px1 * qy0 and px1qy1 = px1 * qy1 and px1qy2 = px1 * qy2 and py0qx0 = py0 * qx0 and py0qx1 = py0 * qx1 and py0qx2 = py0 * qx2 and py1qx0 = py1 * qx0 and py1qx1 = py1 * qx1 and py1qx2 = py1 * qx2 val A = ~px1qy2 + px0qy2 - px1qy0 + py1qx0 + px0qy0 + py1qx2 - py0qx2 - py0qx0 + 2.0 * (px1qy1 - px0qy1 - py1qx1 + py0qx1) and B = 2.0 * (~px1qy1 + px0qy1 + px1qy0 - px0qy0 + py1qx1 - py0qx1 - py1qx0 + py0qx0) and C = ~px1qy0 + px0qy0 + py1qx0 - py0qx0 - px0py1 + px1py0 val discriminant = (B * B) - (4.0 * A * C) in if discriminant < g0i2f 0 then @(huge_val, huge_val) (* No real solutions. *) else let val sqrt_discr = sqrt (discriminant) val t1 = (~B - sqrt_discr) / (A + A) and t2 = (~B + sqrt_discr) / (A + A) fn check_t (t : double) : double = (* The parameter must lie in [0,1], and the intersection point must lie between (px0,py0) and (px1,py1). We will check only the x coordinate. *) if t < 0.0 || 1.0 < t then huge_val else let val x = eval_bernstein_degree2 (@(qx0, qx1, qx2), t) in if x < px0 || px1 < x then huge_val else t end in @(check_t t1, check_t t2) end endfunflat_enough (@(px0 : double, px1 : double, px2 : double), @(py0 : double, py1 : double, py2 : double), tol : double) : bool = (* The quadratic must be given in s-power coefficients. Its px1 and py1 terms are to be removed. Compare an error estimate to the segment length. *) let (* The symmetric power polynomials of degree 2 are 1-t t(1-t) t Conversion from quadratic to linear is effected by removal of the center term, with absolute error bounded by the value of the center coefficient, divided by 4 (because t(1-t) reaches a maximum of 1/4, at t=1/2). *) val error_squared = 0.125 * ((px1 * px1) + (py1 * py1)) and length_squared = (px2 - px0)**2 + (py2 - py0)**2 in error_squared / tol <= length_squared * tol end(* One might be curious why "t@ype" instead of "type". The answer is: the notation "type" is restricted to types that take up the same space as a C void-pointer, which includes ATS pointers, "boxed" types, etc. A "t@ype" can take up any amount of space, and so includes any type there is (except for linear types, which is a whole other subject). For instance, "int", "double", unboxed records, unboxed tuples, and so on. *)fun {a, b : t@ype} (* A polymorphic template function. *)list_any (pred : (a, b) -<cloref1> bool, obj : a, lst : List0 b) : bool = (* Does pred(obj, item) return true for any list item? Here the <cloref1> notation means that pred is a CLOSURE of the ordinary garbage-collected kind, such as functions tend implicitly to be in Lisps, MLs, Haskell, etc. *) case+ lst of | NIL => false | hd :: tl => if pred (obj, hd) then true else list_any (pred, obj, tl)funfind_intersection_parameters (px : @(double, double, double), py : @(double, double, double), qx : @(double, double, double), qy : @(double, double, double), tol : double, spacing : double) : List0 double = let val px = bernstein2spower_degree2 px and py = bernstein2spower_degree2 py fun within_spacing (t_candidate : double, t_in_list : double) :<cloref1> bool = abs (t_candidate - t_in_list) < spacing fun loop {n : nat} (params : list (double, n), n : int n, workload : List0 (@(double, double))) : List0 double = case+ workload of | NIL => params | hd :: tl => let val portionx = spower_portion_degree2 (px, hd) and portiony = spower_portion_degree2 (py, hd) in if flat_enough (portionx, portiony, tol) then let val @(portionx0, _, portionx2) = portionx and @(portiony0, _, portiony2) = portiony val @(root0, root1) = solve_linear_quadratic (@(portionx0, portionx2), @(portiony0, portiony2), qx, qy) in if 0.0 <= root0 && root0 <= 1.0 && ~list_any (within_spacing, root0, params) then begin if 0.0 <= root1 && root1 <= 1.0 && ~list_any (within_spacing, root1, params) then loop (root0 :: root1 :: params, n + 2, tl) else loop (root0 :: params, n + 1, tl) end else if 0.0 <= root1 && root1 <= 1.0 && ~list_any (within_spacing, root1, params) then loop (root1 :: params, n + 1, tl) else loop (params, n, tl) end else let val @(t0, t1) = hd val tmiddle = (0.5 * t0) + (0.5 * t1) val job1 = @(t0, tmiddle) and job2 = @(tmiddle, t1) in loop (params, n, job1 :: job2 :: tl) end end in loop (NIL, 0, @(0.0, 1.0) :: NIL) endimplementmain0 () = let val px = @(~1.0, 0.0, 1.0) val py = @(0.0, 10.0, 0.0) val qx = @(2.0, ~8.0, 2.0) val qy = @(1.0, 2.0, 3.0) val tol = 0.001 (* "Flatness ratio" *) val spacing = 0.0001 (* Min. spacing between parameters. *) val t_list = find_intersection_parameters (px, py, qx, qy, tol, spacing) (* For no particular reason, sort the intersections so they go from top to bottom. *) val t_list = list_vt2t (list_vt_reverse (list_mergesort t_list)) val () = println! ("From top to bottom:") fun loop {n : nat} .<n>. (t_list : list (double, n)) : void = case+ t_list of | NIL => () | t :: tl => begin println! ("(", eval_bernstein_degree2 (qx, t), ", ", eval_bernstein_degree2 (qy, t), ")"); loop tl end in loop t_list endFrom top to bottom:(0.654983, 2.854983)(-0.681024, 2.681025)(-0.854982, 1.345016)(0.881023, 1.118975)
/* The control points of a planar quadratic Bézier curve form a triangle--called the "control polygon"--that completely contains the curve. Furthermore, the rectangle formed by the minimum and maximum x and y values of the control polygon completely contain the polygon, and therefore also the curve. Thus a simple method for narrowing down where intersections might be is: subdivide both curves until you find "small enough" regions where these rectangles overlap.*/#include<stdio.h>#include<stdbool.h>#include<math.h>#include<assert.h>typedefstruct{doublex;doubley;}point;typedefstruct{doublec0;doublec1;doublec2;}quadSpline;// Non-parametric spline.typedefstruct{quadSplinex;quadSpliney;}quadCurve;// Planar parametric spline.// Subdivision by de Casteljau's algorithm.voidsubdivideQuadSpline(quadSplineq,doublet,quadSpline*u,quadSpline*v){doubles=1.0-t;doublec0=q.c0;doublec1=q.c1;doublec2=q.c2;u->c0=c0;v->c2=c2;u->c1=s*c0+t*c1;v->c1=s*c1+t*c2;u->c2=s*u->c1+t*v->c1;v->c0=u->c2;}voidsubdivideQuadCurve(quadCurveq,doublet,quadCurve*u,quadCurve*v){subdivideQuadSpline(q.x,t,&u->x,&v->x);subdivideQuadSpline(q.y,t,&u->y,&v->y);}// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.boolrectsOverlap(doublexa0,doubleya0,doublexa1,doubleya1,doublexb0,doubleyb0,doublexb1,doubleyb1){return(xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1);}doublemax3(doublex,doubley,doublez){returnfmax(fmax(x,y),z);}doublemin3(doublex,doubley,doublez){returnfmin(fmin(x,y),z);}// This accepts the point as an intersection if the boxes are small enough.voidtestIntersect(quadCurvep,quadCurveq,doubletol,bool*exclude,bool*accept,point*intersect){doublepxmin=min3(p.x.c0,p.x.c1,p.x.c2);doublepymin=min3(p.y.c0,p.y.c1,p.y.c2);doublepxmax=max3(p.x.c0,p.x.c1,p.x.c2);doublepymax=max3(p.y.c0,p.y.c1,p.y.c2);doubleqxmin=min3(q.x.c0,q.x.c1,q.x.c2);doubleqymin=min3(q.y.c0,q.y.c1,q.y.c2);doubleqxmax=max3(q.x.c0,q.x.c1,q.x.c2);doubleqymax=max3(q.y.c0,q.y.c1,q.y.c2);*exclude=true;*accept=false;if(rectsOverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)){*exclude=false;doublexmin=fmax(pxmin,qxmin);doublexmax=fmin(pxmax,qxmax);assert(xmax>=xmin);if(xmax-xmin<=tol){doubleymin=fmax(pymin,qymin);doubleymax=fmin(pymax,qymax);assert(ymax>=ymin);if(ymax-ymin<=tol){*accept=true;intersect->x=0.5*xmin+0.5*xmax;intersect->y=0.5*ymin+0.5*ymax;}}}}boolseemsToBeDuplicate(pointintersects[],inticount,pointxy,doublespacing){boolseemsToBeDup=false;inti=0;while(!seemsToBeDup&&i!=icount){pointpt=intersects[i];seemsToBeDup=fabs(pt.x-xy.x)<spacing&&fabs(pt.y-xy.y)<spacing;++i;}returnseemsToBeDup;}voidfindIntersects(quadCurvep,quadCurveq,doubletol,doublespacing,pointintersects[]){intnumIntersects=0;typedefstruct{quadCurvep;quadCurveq;}workset;worksetworkload[64];intnumWorksets=1;workload[0]=(workset){p,q};// Quit looking after having emptied the workload.while(numWorksets!=0){worksetwork=workload[numWorksets-1];--numWorksets;boolexclude,accept;pointintersect;testIntersect(work.p,work.q,tol,&exclude,&accept,&intersect);if(accept){// To avoid detecting the same intersection twice, require some// space between intersections.if(!seemsToBeDuplicate(intersects,numIntersects,intersect,spacing)){intersects[numIntersects++]=intersect;assert(numIntersects<=4);}}elseif(!exclude){quadCurvep0,p1,q0,q1;subdivideQuadCurve(work.p,0.5,&p0,&p1);subdivideQuadCurve(work.q,0.5,&q0,&q1);workload[numWorksets++]=(workset){p0,q0};workload[numWorksets++]=(workset){p0,q1};workload[numWorksets++]=(workset){p1,q0};workload[numWorksets++]=(workset){p1,q1};assert(numWorksets<=64);}}}intmain(){quadCurvep,q;p.x=(quadSpline){-1.0,0.0,1.0};p.y=(quadSpline){0.0,10.0,0.0};q.x=(quadSpline){2.0,-8.0,2.0};q.y=(quadSpline){1.0,2.0,3.0};doubletol=0.0000001;doublespacing=tol*10.0;pointintersects[4];findIntersects(p,q,tol,spacing,intersects);inti;for(i=0;i<4;++i){printf("(% f, %f)\n",intersects[i].x,intersects[i].y);}return0;}
( 0.654983, 2.854983)( 0.881025, 1.118975)(-0.681025, 2.681025)(-0.854983, 1.345017)
Unfortunately two of us were writing C implementations at the same time. Had I known this, I would have written the following in a different language.
This implementation uses a recursive function rather than a "workload container". (Another approach to adaptive bisection requires no extra storage at all. That is to work with integer values in a clever way.)
// If you are using GCC, compile with -std=gnu2x because there may be// C23-isms: [[attributes]], empty () instead of (void), etc./* In this program, both of the curves are adaptively "flattened": that is, converted to a piecewise linear approximation. Then the problem is reduced to finding intersections of line segments. How efficient or inefficient the method is I will not try to answer. (And I do sometimes compute things "too often", although a really good optimizer might fix that.) I will use the symmetric power basis that was introduced by J. Sánchez-Reyes: J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, page 319. J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis in geometry processing’, ACM Transactions on Graphics, vol 19 no 1, January 2000, page 35. Flattening a quadratic that is represented in this basis has a few advantages, which I will not go into here. */#include<stdio.h>#include<stdbool.h>#include<math.h>staticinlinevoiddo_nothing(){}structbernstein_spline{doubleb0;doubleb1;doubleb2;};structspower_spline{doublec0;doublec1;doublec2;};typedefstructbernstein_splinebernstein_spline;typedefstructspower_splinespower_spline;structspower_curve{spower_splinex;spower_spliney;};typedefstructspower_curvespower_curve;// Convert a non-parametric spline from Bernstein basis to s-power.[[gnu::const]]spower_splinebernstein_spline_to_spower(bernstein_splineS){spower_splineT={.c0=S.b0,.c1=(2*S.b1)-S.b0-S.b2,.c2=S.b2};returnT;}// Compose (c0, c1, c2) with (t0, t1). This will map the portion// [t0,t1] onto [0,1]. (To get these expressions, I did not use the// general-degree methods described by Sánchez-Reyes, but instead used// Maxima, some while ago.)//// This method is an alternative to de Casteljau subdivision, and can// be done with the coefficients in any basis. Instead of breaking the// spline into two pieces at a parameter value t, it gives you the// portion lying between two parameter values. In general that// requires two applications of de Casteljau subdivision. On the other// hand, subdivision requires two applications of the following.[[gnu::const]]inlinespower_splinespower_spline_portion(spower_splineS,doublet0,doublet1){doublet0_t0=t0*t0;doublet0_t1=t0*t1;doublet1_t1=t1*t1;doublec2p1m0=S.c2+S.c1-S.c0;spower_splineT={.c0=S.c0+(c2p1m0*t0)-(S.c1*t0_t0),.c1=(S.c1*t1_t1)-(2*S.c1*t0_t1)+(S.c1*t0_t0),.c2=S.c0+(c2p1m0*t1)-(S.c1*t1_t1)};returnT;}[[gnu::const]]inlinespower_curvespower_curve_portion(spower_curveC,doublet0,doublet1){spower_curveD={.x=spower_spline_portion(C.x,t0,t1),.y=spower_spline_portion(C.y,t0,t1)};returnD;}// Given a parametric curve, is it "flat enough" to have its quadratic// terms removed?[[gnu::const]]inlineboolflat_enough(spower_curveC,doubletol){// The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to// remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached// at t=1/2. That accounts for the 1/8=0.125 in the following:doublecx0=C.x.c0;doublecx1=C.x.c1;doublecx2=C.x.c2;doublecy0=C.y.c0;doublecy1=C.y.c1;doublecy2=C.y.c2;doubledx=cx2-cx0;doubledy=cy2-cy0;doubleerror_squared=0.125*((cx1*cx1)+(cy1*cy1));doublelength_squared=(dx*dx)+(dy*dy);return(error_squared<=length_squared*tol*tol);}// Given two line segments, do they intersect? One solution to this// problem is to use the implicitization method employed in the Maxima// example, except to do it with linear instead of quadratic// curves. That is what I do here, with the the roles of who gets// implicitized alternated. If both ways you get as answer a parameter// in [0,1], then the segments intersect.voidtest_line_segment_intersection(doubleax0,doubleax1,doubleay0,doubleay1,doublebx0,doublebx1,doubleby0,doubleby1,bool*they_intersect,double*x,double*y){doubleanumer=((bx1-bx0)*ay0-(by1-by0)*ax0+bx0*by1-bx1*by0);doublebnumer=-((ax1-ax0)*by0-(ay1-ay0)*bx0+ax0*ay1-ax1*ay0);doubledenom=((ax1-ax0)*(by1-by0)-(ay1-ay0)*(bx1-bx0));doubleta=anumer/denom;/* Parameter of segment a. */doubletb=bnumer/denom;/* Parameter of segment b. */*they_intersect=(0<=ta&&ta<=1&&0<=tb&&tb<=1);if(*they_intersect){*x=((1-ta)*ax0)+(ta*ax1);*y=((1-ta)*ay0)+(ta*ay1);}}booltoo_close(doublex,doubley,doublexs[],doubleys[],size_tnum_points,doublespacing){booltoo_close=false;size_ti=0;while(!too_close&&i!=num_points){too_close=(fabs(x-xs[i])<spacing&&fabs(y-ys[i])<spacing);i+=1;}returntoo_close;}voidrecursion(doubletp0,doubletp1,doubletq0,doubletq1,spower_curveP,spower_curveQ,doubletol,doublespacing,size_tmax_points,doublexs[max_points],doubleys[max_points],size_t*num_points){if(*num_points==max_points)do_nothing();elseif(!flat_enough(spower_curve_portion(P,tp0,tp1),tol)){doubletp_half=(0.5*tp0)+(0.5*tp1);if(!(flat_enough(spower_curve_portion(Q,tq0,tq1),tol))){doubletq_half=(0.5*tq0)+(0.5*tq1);recursion(tp0,tp_half,tq0,tq_half,P,Q,tol,spacing,max_points,xs,ys,num_points);recursion(tp0,tp_half,tq_half,tq1,P,Q,tol,spacing,max_points,xs,ys,num_points);recursion(tp_half,tp1,tq0,tq_half,P,Q,tol,spacing,max_points,xs,ys,num_points);recursion(tp_half,tp1,tq_half,tq1,P,Q,tol,spacing,max_points,xs,ys,num_points);}else{recursion(tp0,tp_half,tq0,tq1,P,Q,tol,spacing,max_points,xs,ys,num_points);recursion(tp_half,tp1,tq0,tq1,P,Q,tol,spacing,max_points,xs,ys,num_points);}}elseif(!(flat_enough(spower_curve_portion(Q,tq0,tq1),tol))){doubletq_half=(0.5*tq0)+(0.5*tq1);recursion(tp0,tp1,tq0,tq_half,P,Q,tol,spacing,max_points,xs,ys,num_points);recursion(tp0,tp1,tq_half,tq1,P,Q,tol,spacing,max_points,xs,ys,num_points);}else{spower_curveP1=spower_curve_portion(P,tp0,tp1);spower_curveQ1=spower_curve_portion(Q,tq0,tq1);boolthey_intersect;doublex,y;test_line_segment_intersection(P1.x.c0,P1.x.c2,P1.y.c0,P1.y.c2,Q1.x.c0,Q1.x.c2,Q1.y.c0,Q1.y.c2,&they_intersect,&x,&y);if(they_intersect&&!too_close(x,y,xs,ys,*num_points,spacing)){xs[*num_points]=x;ys[*num_points]=y;*num_points+=1;}}}voidfind_intersections(spower_curveP,spower_curveQ,doubleflatness_tolerance,doublepoint_spacing,size_tmax_points,doublexs[max_points],doubleys[max_points],size_t*num_points){*num_points=0;recursion(0,1,0,1,P,Q,flatness_tolerance,point_spacing,max_points,xs,ys,num_points);}intmain(){bernstein_splinebPx={.b0=-1,.b1=0,.b2=1};bernstein_splinebPy={.b0=0,.b1=10,.b2=0};bernstein_splinebQx={.b0=2,.b1=-8,.b2=2};bernstein_splinebQy={.b0=1,.b1=2,.b2=3};spower_splinePx=bernstein_spline_to_spower(bPx);spower_splinePy=bernstein_spline_to_spower(bPy);spower_splineQx=bernstein_spline_to_spower(bQx);spower_splineQy=bernstein_spline_to_spower(bQy);spower_curveP={.x=Px,.y=Py};spower_curveQ={.x=Qx,.y=Qy};doubleflatness_tolerance=0.001;doublepoint_spacing=0.000001;/* Max norm minimum spacing. */constsize_tmax_points=10;doublexs[max_points];doubleys[max_points];size_tnum_points;find_intersections(P,Q,flatness_tolerance,point_spacing,max_points,xs,ys,&num_points);for(size_ti=0;i!=num_points;i+=1)printf("(%f, %f)\n",xs[i],ys[i]);return0;}
(-0.854982, 1.345017)(-0.681024, 2.681024)(0.881023, 1.118977)(0.654983, 2.854982)
usingSystem;usingSystem.Collections.Generic;publicclassBezierCurveIntersection{publicstaticvoidMain(string[]args){QuadCurvevertical=newQuadCurve(newQuadSpline(-1.0,0.0,1.0),newQuadSpline(0.0,10.0,0.0));// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)QuadCurvehorizontal=newQuadCurve(newQuadSpline(2.0,-8.0,2.0),newQuadSpline(1.0,2.0,3.0));// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)Console.WriteLine("The points of intersection are:");List<Point>intersects=FindIntersects(vertical,horizontal);foreach(Pointintersectinintersects){Console.WriteLine($"( {intersect.X,9:0.000000}, {intersect.Y,9:0.000000} )");}}privatestaticList<Point>FindIntersects(QuadCurvep,QuadCurveq){List<Point>result=newList<Point>();Stack<QuadCurve>stack=newStack<QuadCurve>();stack.Push(p);stack.Push(q);while(stack.Count>0){QuadCurvepp=stack.Pop();QuadCurveqq=stack.Pop();List<object>objects=TestIntersection(pp,qq);boolaccepted=(bool)objects[0];boolexcluded=(bool)objects[1];Pointintersect=(Point)objects[2];if(accepted){if(!SeemsToBeDuplicate(result,intersect)){result.Add(intersect);}}elseif(!excluded){QuadCurvep0=newQuadCurve();QuadCurveq0=newQuadCurve();QuadCurvep1=newQuadCurve();QuadCurveq1=newQuadCurve();SubdivideQuadCurve(pp,0.5,p0,p1);SubdivideQuadCurve(qq,0.5,q0,q1);stack.Push(p0);stack.Push(q0);stack.Push(p0);stack.Push(q1);stack.Push(p1);stack.Push(q0);stack.Push(p1);stack.Push(q1);}}returnresult;}privatestaticboolSeemsToBeDuplicate(List<Point>intersects,Pointpoint){foreach(Pointintersectinintersects){if(Math.Abs(intersect.X-point.X)<Spacing&&Math.Abs(intersect.Y-point.Y)<Spacing){returntrue;}}returnfalse;}privatestaticList<object>TestIntersection(QuadCurvep,QuadCurveq){doublepxMin=Math.Min(Math.Min(p.X.C0,p.X.C1),p.X.C2);doublepyMin=Math.Min(Math.Min(p.Y.C0,p.Y.C1),p.Y.C2);doublepxMax=Math.Max(Math.Max(p.X.C0,p.X.C1),p.X.C2);doublepyMax=Math.Max(Math.Max(p.Y.C0,p.Y.C1),p.Y.C2);doubleqxMin=Math.Min(Math.Min(q.X.C0,q.X.C1),q.X.C2);doubleqyMin=Math.Min(Math.Min(q.Y.C0,q.Y.C1),q.Y.C2);doubleqxMax=Math.Max(Math.Max(q.X.C0,q.X.C1),q.X.C2);doubleqyMax=Math.Max(Math.Max(q.Y.C0,q.Y.C1),q.Y.C2);boolaccepted=false;boolexcluded=true;Pointintersect=newPoint(0.0,0.0);if(RectanglesOverlap(pxMin,pyMin,pxMax,pyMax,qxMin,qyMin,qxMax,qyMax)){excluded=false;doublexMin=Math.Max(pxMin,qxMin);doublexMax=Math.Min(pxMax,pxMax);if(xMax-xMin<=Tolerance){doubleyMin=Math.Max(pyMin,qyMin);doubleyMax=Math.Min(pyMax,qyMax);if(yMax-yMin<=Tolerance){accepted=true;intersect=newPoint(0.5*(xMin+xMax),0.5*(yMin+yMax));}}}returnnewList<object>{accepted,excluded,intersect};}privatestaticboolRectanglesOverlap(doublexa0,doubleya0,doublexa1,doubleya1,doublexb0,doubleyb0,doublexb1,doubleyb1){returnxb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1;}privatestaticvoidSubdivideQuadCurve(QuadCurveq,doublet,QuadCurveu,QuadCurvev){SubdivideQuadSpline(q.X,t,u.X,v.X);SubdivideQuadSpline(q.Y,t,u.Y,v.Y);}// de Casteljau's algorithmprivatestaticvoidSubdivideQuadSpline(QuadSplineq,doublet,QuadSplineu,QuadSplinev){doubles=1.0-t;u.C0=q.C0;v.C2=q.C2;u.C1=s*q.C0+t*q.C1;v.C1=s*q.C1+t*q.C2;u.C2=s*u.C1+t*v.C1;v.C0=u.C2;}publicstructPoint{publicdoubleX{get;}publicdoubleY{get;}publicPoint(doublex,doubley){X=x;Y=y;}}publicclassQuadSpline{publicdoubleC0{get;set;}publicdoubleC1{get;set;}publicdoubleC2{get;set;}publicQuadSpline(doublec0,doublec1,doublec2){C0=c0;C1=c1;C2=c2;}publicQuadSpline():this(0.0,0.0,0.0){}}publicclassQuadCurve{publicQuadSplineX{get;set;}publicQuadSplineY{get;set;}publicQuadCurve(QuadSplinex,QuadSpliney){X=x;Y=y;}publicQuadCurve():this(newQuadSpline(),newQuadSpline()){}}privateconstdoubleTolerance=0.000_000_1;privateconstdoubleSpacing=10*Tolerance;}
The points of intersection are:( 0.654983, 2.854983 )( -0.681025, 2.681025 )( 0.881025, 1.118975 )( -0.854983, 1.345017 )
#include<cmath>#include<iomanip>#include<iostream>#include<stack>#include<tuple>#include<vector>constdoubleTOLERANCE=0.000'000'1;constdoubleSPACING=10*TOLERANCE;typedefstd::pair<double,double>point;classquad_spline{public:quad_spline(doublec0,doublec1,doublec2):c0(c0),c1(c1),c2(c2){};quad_spline():c0(0.0),c1(0.0),c2(0.0){};doublec0,c1,c2;};classquad_curve{public:quad_curve(quad_splinex,quad_spliney):x(x),y(y){};quad_curve():x(quad_spline()),y(quad_spline()){};quad_splinex,y;};// de Casteljau's algorithmvoidsubdivide_quad_spline(constquad_spline&q,constdouble&t,quad_spline&u,quad_spline&v){constdoubles=1.0-t;u.c0=q.c0;v.c2=q.c2;u.c1=s*q.c0+t*q.c1;v.c1=s*q.c1+t*q.c2;u.c2=s*u.c1+t*v.c1;v.c0=u.c2;}voidsubdivide_quad_curve(constquad_curve&q,constdoublet,quad_curve&u,quad_curve&v){subdivide_quad_spline(q.x,t,u.x,v.x);subdivide_quad_spline(q.y,t,u.y,v.y);}boolrectangles_overlap(constdouble&xa0,constdouble&ya0,constdouble&xa1,constdouble&ya1,constdouble&xb0,constdouble&yb0,constdouble&xb1,constdouble&yb1){returnxb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1;}std::tuple<bool,bool,point>test_intersection(constquad_curve&p,constquad_curve&q){constdoublepx_min=std::min(std::min(p.x.c0,p.x.c1),p.x.c2);constdoublepy_min=std::min(std::min(p.y.c0,p.y.c1),p.y.c2);constdoublepx_max=std::max(std::max(p.x.c0,p.x.c1),p.x.c2);constdoublepy_max=std::max(std::max(p.y.c0,p.y.c1),p.y.c2);constdoubleqx_min=std::min(std::min(q.x.c0,q.x.c1),q.x.c2);constdoubleqy_min=std::min(std::min(q.y.c0,q.y.c1),q.y.c2);constdoubleqx_max=std::max(std::max(q.x.c0,q.x.c1),q.x.c2);constdoubleqy_max=std::max(std::max(q.y.c0,q.y.c1),q.y.c2);boolaccepted=false;boolexcluded=true;pointintersect=std::make_pair(0.0,0.0);if(rectangles_overlap(px_min,py_min,px_max,py_max,qx_min,qy_min,qx_max,qy_max)){excluded=false;constdoublex_min=std::max(px_min,qx_min);constdoublex_max=std::min(px_max,px_max);if(x_max-x_min<=TOLERANCE){constdoubley_min=std::max(py_min,qy_min);constdoubley_max=std::min(py_max,qy_max);if(y_max-y_min<=TOLERANCE){accepted=true;intersect=std::make_pair(0.5*(x_min+x_max),0.5*(y_min+y_max));}}}returnstd::make_tuple(accepted,excluded,intersect);}boolseems_to_be_duplicate(conststd::vector<point>&intersects,constpoint&pt){for(pointintersect:intersects){if(fabs(intersect.first-pt.first)<SPACING&&fabs(intersect.second-pt.second)<SPACING){returntrue;}}returnfalse;}std::vector<point>find_intersects(constquad_curve&p,constquad_curve&q){std::vector<point>result;std::stack<quad_curve>stack;stack.push(p);stack.push(q);while(!stack.empty()){quad_curvepp=stack.top();stack.pop();quad_curveqq=stack.top();stack.pop();std::tuple<bool,bool,point>objects=test_intersection(pp,qq);boolaccepted,excluded;pointintersect;std::tie(accepted,excluded,intersect)=objects;if(accepted){if(!seems_to_be_duplicate(result,intersect)){result.push_back(intersect);}}elseif(!excluded){quad_curvep0{},q0{},p1{},q1{};subdivide_quad_curve(pp,0.5,p0,p1);subdivide_quad_curve(qq,0.5,q0,q1);for(quad_curveitem:{p0,q0,p0,q1,p1,q0,p1,q1}){stack.push(item);}}}returnresult;}intmain(){quad_curvevertical(quad_spline(-1.0,0.0,1.0),quad_spline(0.0,10.0,0.0));// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)quad_curvehorizontal(quad_spline(2.0,-8.0,2.0),quad_spline(1.0,2.0,3.0));// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)std::cout<<"The points of intersection are:"<<std::endl;std::vector<point>intersects=find_intersects(vertical,horizontal);for(constpoint&intersect:intersects){std::cout<<"( "<<std::setw(9)<<std::setprecision(6)<<intersect.first<<", "<<intersect.second<<" )"<<std::endl;}}
The points of intersection are:( 0.654983, 2.85498 )( -0.681025, 2.68102 )( 0.881025, 1.11898 )( -0.854983, 1.34502 )
This program subdivides both curves by de Casteljau's algorithm, until only very small subdivisions with overlapping control polygons remain. (You could use recursion instead of the workload container. With the container it is easier to terminate early, and also the program then uses only constant stack space.)
tol, so that the check sometimes comes out positive.// The control points of a planar quadratic Bézier curve form a// triangle--called the "control polygon"--that completely contains// the curve. Furthermore, the rectangle formed by the minimum and// maximum x and y values of the control polygon completely contain// the polygon, and therefore also the curve.//// Thus a simple method for narrowing down where intersections might// be is: subdivide both curves until you find "small enough" regions// where these rectangles overlap.importstd.algorithm;importstd.container.slist;importstd.math;importstd.range;importstd.stdio;structpoint{doublex,y;}structquadratic_spline// Non-parametric spline.{doublec0,c1,c2;}structquadratic_curve// Planar parametric spline.{quadratic_splinex,y;}voidsubdivide_quadratic_spline(quadratic_splineq,doublet,refquadratic_splineu,refquadratic_splinev){// Subdivision by de Casteljau's algorithm.immutables=1-t;immutablec0=q.c0;immutablec1=q.c1;immutablec2=q.c2;u.c0=c0;v.c2=c2;u.c1=(s*c0)+(t*c1);v.c1=(s*c1)+(t*c2);u.c2=(s*u.c1)+(t*v.c1);v.c0=u.c2;}voidsubdivide_quadratic_curve(quadratic_curveq,doublet,refquadratic_curveu,refquadratic_curvev){subdivide_quadratic_spline(q.x,t,u.x,v.x);subdivide_quadratic_spline(q.y,t,u.y,v.y);}boolrectangles_overlap(doublexa0,doubleya0,doublexa1,doubleya1,doublexb0,doubleyb0,doublexb1,doubleyb1){// It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1.return(xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1);}voidtest_intersection(quadratic_curvep,quadratic_curveq,doubletol,refboolexclude,refboolaccept,refpointintersection){// I will not do a lot of checking for intersections, as one might// wish to do in a particular application. If the boxes are small// enough, I will accept the point as an intersection.immutablepxmin=min(p.x.c0,p.x.c1,p.x.c2);immutablepymin=min(p.y.c0,p.y.c1,p.y.c2);immutablepxmax=max(p.x.c0,p.x.c1,p.x.c2);immutablepymax=max(p.y.c0,p.y.c1,p.y.c2);immutableqxmin=min(q.x.c0,q.x.c1,q.x.c2);immutableqymin=min(q.y.c0,q.y.c1,q.y.c2);immutableqxmax=max(q.x.c0,q.x.c1,q.x.c2);immutableqymax=max(q.y.c0,q.y.c1,q.y.c2);exclude=true;accept=false;if(rectangles_overlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)){exclude=false;immutablexmin=max(pxmin,qxmin);immutablexmax=min(pxmax,qxmax);assert(xmax>=xmin);if(xmax-xmin<=tol){immutableymin=max(pymin,qymin);immutableymax=min(pymax,qymax);assert(ymax>=ymin);if(ymax-ymin<=tol){accept=true;intersection=point((0.5*xmin)+(0.5*xmax),(0.5*ymin)+(0.5*ymax));}}}}boolseems_to_be_a_duplicate(point[]intersections,pointxy,doublespacing){boolseems_to_be_dup=false;inti=0;while(!seems_to_be_dup&&i!=intersections.length){immutablept=intersections[i];seems_to_be_dup=fabs(pt.x-xy.x)<spacing&&fabs(pt.y-xy.y)<spacing;i+=1;}returnseems_to_be_dup;}point[]find_intersections(quadratic_curvep,quadratic_curveq,doubletol,doublespacing){point[]intersections;intnum_intersections=0;structworkset{quadratic_curvep,q;}SList!worksetworkload;// Initial workload.workload.insertFront(workset(p,q));// Quit looking after having /*found four intersections*/ or emptied// the workload.while(/*num_intersections != 4 &&*/!workload.empty){autowork=workload.front;workload.removeFront();boolexclude;boolaccept;pointintersection;test_intersection(work.p,work.q,tol,exclude,accept,intersection);if(accept){// This is a crude way to avoid detecting the same// intersection twice: require some space between// intersections. For, say, font design work, this method// should be adequate.if(!seems_to_be_a_duplicate(intersections,intersection,spacing)){intersections.length=num_intersections+1;intersections[num_intersections]=intersection;num_intersections+=1;}}elseif(!exclude){quadratic_curvep0,p1,q0,q1;subdivide_quadratic_curve(work.p,0.5,p0,p1);subdivide_quadratic_curve(work.q,0.5,q0,q1);workload.insertFront(workset(p0,q0));workload.insertFront(workset(p0,q1));workload.insertFront(workset(p1,q0));workload.insertFront(workset(p1,q1));}}returnintersections;}intmain(){quadratic_curvep,q;p.x.c0=-1.0;p.x.c1=0.0;p.x.c2=1.0;p.y.c0=0.0;p.y.c1=10.0;p.y.c2=0.0;q.x.c0=2.0;q.x.c1=-8.0;q.x.c2=2.0;q.y.c0=1.0;q.y.c1=2.0;q.y.c2=3.0;immutabletol=0.0000001;immutablespacing=10*tol;autointersections=find_intersections(p,q,tol,spacing);for(inti=0;i!=intersections.length;i+=1)printf("(%f, %f)\n",intersections[i].x,intersections[i].y);return0;}
(0.654983, 2.854983)(0.881025, 1.118975)(-0.681025, 2.681025)(-0.854983, 1.345017)
programbezier_intersections! This program does not do any subdivision, but instead takes! advantage of monotonicity.!! It is possible for points accidentally to be counted twice, for! instance if they lie right on an interval boundary. We will avoid! that by the crude (but likely satisfactory) mechanism of requiring! a minimum max norm between intersections.implicit none! Constantsinteger,parameter::max_intersections=4integer,parameter::max_start_interv=4integer,parameter::max_workload=1000! Maximum size for workload stackdouble precision,parameter::px0=-1.0,px1=0.0,px2=1.0double precision,parameter::py0=0.0,py1=10.0,py2=0.0double precision,parameter::qx0=2.0,qx1=-8.0,qx2=2.0double precision,parameter::qy0=1.0,qy1=2.0,qy2=3.0double precision,parameter::tol=0.0000001double precision,parameter::spacing=100.0*tol! Variableslogical::px_has_extreme_pt,py_has_extreme_ptlogical::qx_has_extreme_pt,qy_has_extreme_ptdouble precision::px_extreme_pt,py_extreme_ptdouble precision::qx_extreme,qy_extremeinteger::p_num_start_interv,q_num_start_intervdouble precision::p_start_interv(max_start_interv)double precision::q_start_interv(max_start_interv)double precision::workload(4,max_workload)! tp0, tp1, tq0, tq1integer::workload_sizeinteger::num_intersectionsdouble precision::intersections_x(max_intersections)double precision::intersections_y(max_intersections)integer::i,j,kdouble precision::tp0,tp1,tq0,tq1double precision::xp0,xp1,xq0,xq1double precision::yp0,yp1,yq0,yq1logical::exclude,acceptdouble precision::x,ydouble precision::tp_middle,tq_middle! Main program logic! Find monotonic sections of the curves, and use those as the! starting jobs.p_num_start_interv=2p_start_interv(1)=0.0p_start_interv(2)=1.0callpossibly_insert_extreme_point(px0,px1,px2,p_num_start_interv,p_start_interv)callpossibly_insert_extreme_point(py0,py1,py2,p_num_start_interv,p_start_interv)q_num_start_interv=2q_start_interv(1)=0.0q_start_interv(2)=1.0callpossibly_insert_extreme_point(qx0,qx1,qx2,q_num_start_interv,q_start_interv)callpossibly_insert_extreme_point(qy0,qy1,qy2,q_num_start_interv,q_start_interv)workload_size=0doi=2,p_num_start_intervdoj=2,q_num_start_intervcalldefer_work(p_start_interv(i-1),p_start_interv(i),&q_start_interv(j-1),q_start_interv(j))end do end do! Go through the workload, deferring work as necessary.num_intersections=0do while(.not.work_is_done())! The following code recomputes values of the splines! sometimes. You may wish to store such values in the work pile,! to avoid recomputing them.calldo_some_work(tp0,tp1,tq0,tq1)xp0=schumaker_volk(px0,px1,px2,tp0)yp0=schumaker_volk(py0,py1,py2,tp0)xp1=schumaker_volk(px0,px1,px2,tp1)yp1=schumaker_volk(py0,py1,py2,tp1)xq0=schumaker_volk(qx0,qx1,qx2,tq0)yq0=schumaker_volk(qy0,qy1,qy2,tq0)xq1=schumaker_volk(qx0,qx1,qx2,tq1)yq1=schumaker_volk(qy0,qy1,qy2,tq1)calltest_intersection(xp0,xp1,yp0,yp1,xq0,xq1,yq0,yq1,tol,exclude,accept,x,y)if(accept)then callmaybe_add_intersection(x,y,spacing)else if(.not.exclude)thentp_middle=0.5*(tp0+tp1)tq_middle=0.5*(tq0+tq1)calldefer_work(tp0,tp_middle,tq0,tq_middle)calldefer_work(tp0,tp_middle,tq_middle,tq1)calldefer_work(tp_middle,tp1,tq0,tq_middle)calldefer_work(tp_middle,tp1,tq_middle,tq1)end if end do if(num_intersections==0)then print*,'no intersections'else dok=1,num_intersectionswrite(*,'(A, F12.8, A, F12.8, A)')'(',intersections_x(k),', ',intersections_y(k),')'end do end ifcontains! Schumaker's and Volk's algorithm for evaluating a Bezier spline in! Bernstein basis. This is faster than de Casteljau, though not quite! as numerical stable.double precisionfunctionschumaker_volk(c0,c1,c2,t)double precision,intent(IN)::c0,c1,c2,tdouble precision::s,u,vs=1.0-tif(t<=0.5)then! Horner form in the variable u = t/s, taking into account the! binomial coefficients = 1,2,1.u=t/sv=c0+(u*(c1+c1+(u*c2)))! Multiply by s raised to the degree of the spline.v=v*s*selse! Horner form in the variable u = s/t, taking into account the! binomial coefficients = 1,2,1.u=s/tv=c2+(u*(c1+c1+(u*c0)))! Multiply by t raised to the degree of the spline.v=v*t*tend ifschumaker_volk=vend functionschumaker_volk! Find extreme point of a quadratic Bezier splinesubroutinefind_extreme_point(c0,c1,c2,lies_inside_01,extreme_point)double precision,intent(IN)::c0,c1,c2logical,intent(OUT)::lies_inside_01double precision,intent(OUT)::extreme_pointdouble precision::numer,denom! If the spline has c0=c2 but not c0=c1=c2, then treat it as having! an extreme point at 0.5.if(c0==c2.and.c0/=c1)thenlies_inside_01=.true.extreme_point=0.5else! Find the root of the derivative of the spline.lies_inside_01=.false.numer=c0-c1denom=c2-c1-c1+c0if(denom/=0.0.and.numer*denom>=0.0.and.numer<=denom)thenlies_inside_01=.true.extreme_point=numer/denomend if end if end subroutinefind_extreme_point! Insert extreme point into start intervals if it lies in (0,1)subroutinepossibly_insert_extreme_point(c0,c1,c2,num_start_interv,start_interv)double precision,intent(IN)::c0,c1,c2integer,intent(INOUT)::num_start_intervdouble precision,intent(INOUT)::start_interv(max_start_interv)logical::lies_inside_01double precision::extreme_ptcallfind_extreme_point(c0,c1,c2,lies_inside_01,extreme_pt)if(lies_inside_01.and.extreme_pt>0.0.and.extreme_pt<1.0)then if(num_start_interv==2)thenstart_interv(3)=1.0start_interv(2)=extreme_ptnum_start_interv=3else if(extreme_pt<start_interv(2))thenstart_interv(4)=1.0start_interv(3)=start_interv(2)start_interv(2)=extreme_ptnum_start_interv=4else if(extreme_pt>start_interv(2))thenstart_interv(4)=1.0start_interv(3)=extreme_ptnum_start_interv=4end if end if end subroutinepossibly_insert_extreme_point! Minimum of two valuesdouble precisionfunctionminimum2(x,y)double precision,intent(IN)::x,yminimum2=min(x,y)end functionminimum2! Maximum of two valuesdouble precisionfunctionmaximum2(x,y)double precision,intent(IN)::x,ymaximum2=max(x,y)end functionmaximum2! Check if two rectangles overlaplogicalfunctionrectangles_overlap(xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1)double precision,intent(IN)::xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1! It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1.rectangles_overlap=(xb0<=xa1.and.xa0<=xb1.and.yb0<=ya1.and.ya0<=yb1)end functionrectangles_overlap! Test for intersection between two line segmentssubroutinetest_intersection(xp0,xp1,yp0,yp1,xq0,xq1,yq0,yq1,tol,exclude,accept,x,y)double precision,intent(IN)::xp0,xp1,yp0,yp1,xq0,xq1,yq0,yq1,tollogical,intent(OUT)::exclude,acceptdouble precision,intent(OUT)::x,ydouble precision::xpmin,ypmin,xpmax,ypmaxdouble precision::xqmin,yqmin,xqmax,yqmaxdouble precision::xmin,xmax,ymin,ymaxxpmin=minimum2(xp0,xp1)ypmin=minimum2(yp0,yp1)xpmax=maximum2(xp0,xp1)ypmax=maximum2(yp0,yp1)xqmin=minimum2(xq0,xq1)yqmin=minimum2(yq0,yq1)xqmax=maximum2(xq0,xq1)yqmax=maximum2(yq0,yq1)exclude=.true.accept=.false.if(rectangles_overlap(xpmin,ypmin,xpmax,ypmax,xqmin,yqmin,xqmax,yqmax))thenexclude=.false.xmin=maximum2(xpmin,xqmin)xmax=minimum2(xpmax,xqmax)if(xmax<xmin)stop'Assertion failed: xmax >= xmin'if(xmax-xmin<=tol)thenymin=maximum2(ypmin,yqmin)ymax=minimum2(ypmax,yqmax)if(ymax<ymin)stop'Assertion failed: ymax >= ymin'if(ymax-ymin<=tol)then accept=.true.x=0.5*(xmin+xmax)y=0.5*(ymin+ymax)end if end if end if end subroutinetest_intersection! Check if workload is emptylogicalfunctionwork_is_done()work_is_done=(workload_size==0)end functionwork_is_done! Add work to the workload stacksubroutinedefer_work(tp0,tp1,tq0,tq1)double precision,intent(IN)::tp0,tp1,tq0,tq1if(workload_size>=max_workload)stop'Error: Workload stack overflow'workload_size=workload_size+1workload(1,workload_size)=tp0workload(2,workload_size)=tp1workload(3,workload_size)=tq0workload(4,workload_size)=tq1end subroutinedefer_work! Remove and return work from the workload stacksubroutinedo_some_work(tp0,tp1,tq0,tq1)double precision,intent(OUT)::tp0,tp1,tq0,tq1if(work_is_done())stop'Assertion failed: Workload is empty'tp0=workload(1,workload_size)tp1=workload(2,workload_size)tq0=workload(3,workload_size)tq1=workload(4,workload_size)workload_size=workload_size-1end subroutinedo_some_work! Add intersection point if it's not too close to existing onessubroutinemaybe_add_intersection(x,y,spacing)double precision,intent(IN)::x,y,spacinginteger::ilogical::too_closeif(num_intersections==0)thenintersections_x(1)=xintersections_y(1)=ynum_intersections=1elsetoo_close=.false.doi=1,num_intersectionsif(abs(x-intersections_x(i))<spacing.and.&abs(y-intersections_y(i))<spacing)thentoo_close=.true.exit end if end do if(.not.too_close)then if(num_intersections>=max_intersections)stop'Too many intersections'num_intersections=num_intersections+1intersections_x(num_intersections)=xintersections_y(num_intersections)=yend if end if end subroutinemaybe_add_intersectionend programbezier_intersections
( 0.65498343, 2.85498342)( 0.88102499, 1.11897501)( -0.68102493, 2.68102500)( -0.85498342, 1.34501657)
' The control points of a planar quadratic Bézier curve form a' triangle--called the "control polygon"--that completely contains' the curve. Furthermore, the rectangle formed by the minimum and' maximum x and y values of the control polygon completely contain' the polygon, and therefore also the curve.'' Thus a simple method for narrowing down where intersections might' be is: subdivide both curves until you find "small enough" regions' where these rectangles overlap.#defineMin(a,b)iif((a)<(b),(a),(b))#defineMax(a,b)iif((a)>(b),(a),(b))' Note that these are all mutable as we need to pass by reference.TypePuntoxAsDoubleyAsDoubleEndTypeTypeQuadSplinec0AsDoublec1AsDoublec2AsDouble' Non-parametric splineEndTypeTypeQuadCurvexAsQuadSplineyAsQuadSpline' Planar parametric splineEndType' Subdivision by de Casteljau's algorithmSubsubdivideQuadSpline(qAsQuadSpline,tAsDouble,uAsQuadSpline,vAsQuadSpline)DimAsDoubles=1.0-tDimAsDoublec0=q.c0DimAsDoublec1=q.c1DimAsDoublec2=q.c2u.c0=c0v.c2=c2u.c1=s*c0+t*c1v.c1=s*c1+t*c2u.c2=s*u.c1+t*v.c1v.c0=u.c2EndSubSubsubdivideQuadCurve(qAsQuadCurve,tAsDouble,uAsQuadCurve,vAsQuadCurve)subdivideQuadSpline(q.x,t,u.x,v.x)subdivideQuadSpline(q.y,t,u.y,v.y)EndSub' It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.FunctionrectsOverlap(xa0AsDouble,ya0AsDouble,xa1AsDouble,ya1AsDouble,_xb0AsDouble,yb0AsDouble,xb1AsDouble,yb1AsDouble)AsBooleanReturn(xb0<=xa1Andxa0<=xb1Andyb0<=ya1Andya0<=yb1)EndFunctionFunctionmax3(xAsDouble,yAsDouble,zAsDouble)AsDoubleReturnMax(Max(x,y),z)EndFunctionFunctionmin3(xAsDouble,yAsDouble,zAsDouble)AsDoubleReturnMin(Min(x,y),z)EndFunction' This accepts the point as an intersection if the boxes are small enough.SubtestIntersect(pAsQuadCurve,qAsQuadCurve,tolAsDouble,_ByrefexcludeAsBoolean,ByrefacceptAsBoolean,ByrefintersectAsPunto)DimAsDoublepxmin=min3(p.x.c0,p.x.c1,p.x.c2)DimAsDoublepymin=min3(p.y.c0,p.y.c1,p.y.c2)DimAsDoublepxmax=max3(p.x.c0,p.x.c1,p.x.c2)DimAsDoublepymax=max3(p.y.c0,p.y.c1,p.y.c2)DimAsDoubleqxmin=min3(q.x.c0,q.x.c1,q.x.c2)DimAsDoubleqymin=min3(q.y.c0,q.y.c1,q.y.c2)DimAsDoubleqxmax=max3(q.x.c0,q.x.c1,q.x.c2)DimAsDoubleqymax=max3(q.y.c0,q.y.c1,q.y.c2)exclude=Trueaccept=FalseIfrectsOverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)Thenexclude=FalseDimAsDoublexmin=Max(pxmin,qxmin)DimAsDoublexmax=Min(pxmax,qxmax)Assert(xmax>=xmin)Ifxmax-xmin<=tolThenDimAsDoubleymin=Max(pymin,qymin)DimAsDoubleymax=Min(pymax,qymax)Assert(ymax>=ymin)Ifymax-ymin<=tolThenaccept=Trueintersect.x=0.5*xmin+0.5*xmaxintersect.y=0.5*ymin+0.5*ymaxEndIfEndIfEndIfEndSubFunctionseemsToBeDuplicate(intersects()AsPunto,icountAsInteger,_xyAsPunto,spacingAsDouble)AsBooleanDimAsBooleanseemsToBeDup=FalseDimAsIntegeri=0WhileNotseemsToBeDupAndi<>icountDimAsPuntopt=intersects(i)seemsToBeDup=Abs(pt.x-xy.x)<spacingAndAbs(pt.y-xy.y)<spacingi+=1WendReturnseemsToBeDupEndFunctionSubfindIntersects(pAsQuadCurve,qAsQuadCurve,tolAsDouble,_spacingAsDouble,intersects()AsPunto)DimAsIntegernumIntersects=0TypeworksetpAsQuadCurveqAsQuadCurveEndTypeDimAsworksetworkload(64)DimAsIntegernumWorksets=1workload(0)=Type<workset>(p,q)' Quit looking after having emptied the workload.WhilenumWorksets<>0DimAsworksetwork=workload(numWorksets-1)numWorksets-=1DimAsBooleanexclude,acceptDimAsPuntointersecttestIntersect(work.p,work.q,tol,exclude,accept,intersect)IfacceptThen' To avoid detecting the same intersection twice, require some' space between intersections.IfNotseemsToBeDuplicate(intersects(),numIntersects,intersect,spacing)Thenintersects(numIntersects)=intersectnumIntersects+=1EndIfElseifNotexcludeThenDimAsQuadCurvep0,p1,q0,q1subdivideQuadCurve(work.p,0.5,p0,p1)subdivideQuadCurve(work.q,0.5,q0,q1)workload(numWorksets)=Type<workset>(p0,q0)numWorksets+=1workload(numWorksets)=Type<workset>(p0,q1)numWorksets+=1workload(numWorksets)=Type<workset>(p1,q0)numWorksets+=1workload(numWorksets)=Type<workset>(p1,q1)numWorksets+=1EndIfWendEndSubDimAsQuadCurvep,qp.x=Type<QuadSpline>(-1.0,0.0,1.0)p.y=Type<QuadSpline>(0.0,10.0,0.0)q.x=Type<QuadSpline>(2.0,-8.0,2.0)q.y=Type<QuadSpline>(1.0,2.0,3.0)DimAsDoubletol=0.0000001DimAsDoublespacing=tol*10.0DimAsPuntointersects(4)findIntersects(p,q,tol,spacing,intersects())ForiAsInteger=0To3Print"(";intersects(i).x;", ";intersects(i).y;")"NextiSleep
( 0.6549834311008453, 2.854983419179916)( 0.8810249865055084, 1.118975013494492)(-0.6810249347860431, 2.681024998426437)(-0.8549834191799164, 1.345016568899155)
' In this program, both of the curves are adaptively "flattened":' that is, converted to a piecewise linear approximation. Then the' problem is reduced to finding intersections of line segments.'' How efficient or inefficient the method is I will not try to answer.' (And I do sometimes compute things "too often", although a really good' optimizer might fix that.)'' I will use the symmetric power basis that was introduced by J. Sánchez-Reyes:'' J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power' basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, page 319.'' J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis' in geometry processing’, ACM Transactions on Graphics, vol 19' no 1, January 2000, page 35.'' Flattening a quadratic that is represented in this basis has a few' advantages, which I will not go into here. */Typebernstein_splineb0AsDoubleb1AsDoubleb2AsDoubleEndTypeTypespower_splinec0AsDoublec1AsDoublec2AsDoubleEndTypeTypespower_curvexAsspower_splineyAsspower_splineEndType' Convert a non-parametric spline from Bernstein basis to s-power.Functionbernstein_spline_to_spower(SAsbernstein_spline)Asspower_splineDimAsspower_splineTT.c0=S.b0T.c1=(2*S.b1)-S.b0-S.b2T.c2=S.b2ReturnTEndFunction' Compose (c0, c1, c2) with (t0, t1). This will map the portion [t0,t1] onto' [0,1]. (To get these expressions, I did not use the general-degree methods' described by Sánchez-Reyes, but instead used Maxima, some while ago.)'' This method is an alternative to de Casteljau subdivision, and can be done' with the coefficients in any basis. Instead of breaking the spline into two' pieces at a parameter value t, it gives you the portion lying between two' parameter values. In general that requires two applications of de Casteljau' subdivision. On the other hand, subdivision requires two applications of the' following.Functionspower_spline_portion(SAsspower_spline,t0AsDouble,t1AsDouble)Asspower_splineDimAsDoublet0_t0=t0*t0DimAsDoublet0_t1=t0*t1DimAsDoublet1_t1=t1*t1DimAsDoublec2p1m0=S.c2+S.c1-S.c0DimAsspower_splineTT.c0=S.c0+(c2p1m0*t0)-(S.c1*t0_t0)T.c1=(S.c1*t1_t1)-(2*S.c1*t0_t1)+(S.c1*t0_t0)T.c2=S.c0+(c2p1m0*t1)-(S.c1*t1_t1)ReturnTEndFunctionFunctionspower_curve_portion(CAsspower_curve,t0AsDouble,t1AsDouble)Asspower_curveDimAsspower_curveDD.x=spower_spline_portion(C.x,t0,t1)D.y=spower_spline_portion(C.y,t0,t1)ReturnDEndFunction' Given a parametric curve, is it "flat enough" to have its quadratic' terms removed?Functionflat_enough(CAsspower_curve,tolAsDouble)AsBoolean' The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to' remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached' at t=1/2. That accounts for the 1/8=0.125 in the following:DimAsDoublecx0=C.x.c0DimAsDoublecx1=C.x.c1DimAsDoublecx2=C.x.c2DimAsDoublecy0=C.y.c0DimAsDoublecy1=C.y.c1DimAsDoublecy2=C.y.c2DimAsDoubledx=cx2-cx0DimAsDoubledy=cy2-cy0DimAsDoubleerror_squared=0.125*((cx1*cx1)+(cy1*cy1))DimAsDoublelength_squared=(dx*dx)+(dy*dy)Return(error_squared<=length_squared*tol*tol)EndFunction' Given two line segments, do they intersect? One solution to this problem is' to use the implicitization method employed in the Maxima example, except to' do it with linear instead of quadratic curves. That is what I do here, with' the the roles of who gets implicitized alternated. If both ways you get as' answer a parameter in [0,1], then the segments intersect.Subtest_line_segment_intersection(ax0AsDouble,ax1AsDouble,_ay0AsDouble,ay1AsDouble,bx0AsDouble,bx1AsDouble,_by0AsDouble,by1AsDouble,Byrefthey_intersectAsBoolean,_ByrefxAsDouble,ByrefyAsDouble)DimAsDoubleanumer=((bx1-bx0)*ay0-(by1-by0)*ax0_+bx0*by1-bx1*by0)DimAsDoublebnumer=-((ax1-ax0)*by0-(ay1-ay0)*bx0_+ax0*ay1-ax1*ay0)DimAsDoubledenom=((ax1-ax0)*(by1-by0)_-(ay1-ay0)*(bx1-bx0))DimAsDoubleta=anumer/denom' Parameter of segment a.DimAsDoubletb=bnumer/denom' Parameter of segment b.they_intersect=(0<=taAndta<=1And0<=tbAndtb<=1)Ifthey_intersectThenx=((1-ta)*ax0)+(ta*ax1)y=((1-ta)*ay0)+(ta*ay1)EndIfEndSubFunctiontoo_close(xAsDouble,yAsDouble,xs()AsDouble,ys()AsDouble,_num_pointsAsInteger,spacingAsDouble)AsBooleanDimAsBooleantooclose=FalseDimAsIntegeri=0WhileNottoocloseAndi<>num_pointstooclose=(Abs(x-xs(i))<spacingAndAbs(y-ys(i))<spacing)i+=1WendReturntoocloseEndFunctionSubrecursion(tp0AsDouble,tp1AsDouble,tq0AsDouble,tq1AsDouble,_PAsspower_curve,QAsspower_curve,tolAsDouble,spacingAsDouble,_max_pointsAsInteger,xs()AsDouble,ys()AsDouble,Byrefnum_pointsAsInteger)Ifnum_points=max_pointsThen' do nothingElseifNotflat_enough(spower_curve_portion(P,tp0,tp1),tol)ThenDimAsDoubletp_half=(0.5*tp0)+(0.5*tp1)IfNotflat_enough(spower_curve_portion(Q,tq0,tq1),tol)ThenDimAsDoubletq_half=(0.5*tq0)+(0.5*tq1)recursion(tp0,tp_half,tq0,tq_half,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)recursion(tp0,tp_half,tq_half,tq1,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)recursion(tp_half,tp1,tq0,tq_half,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)recursion(tp_half,tp1,tq_half,tq1,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)Elserecursion(tp0,tp_half,tq0,tq1,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)recursion(tp_half,tp1,tq0,tq1,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)EndIfElseifNotflat_enough(spower_curve_portion(Q,tq0,tq1),tol)ThenDimAsDoubletq_half=(0.5*tq0)+(0.5*tq1)recursion(tp0,tp1,tq0,tq_half,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)recursion(tp0,tp1,tq_half,tq1,P,Q,tol,_spacing,max_points,xs(),ys(),num_points)ElseDimAsspower_curveP1=spower_curve_portion(P,tp0,tp1)DimAsspower_curveQ1=spower_curve_portion(Q,tq0,tq1)DimAsBooleanthey_intersectDimAsDoublex,ytest_line_segment_intersection(P1.x.c0,P1.x.c2,_P1.y.c0,P1.y.c2,_Q1.x.c0,Q1.x.c2,_Q1.y.c0,Q1.y.c2,_they_intersect,x,y)Ifthey_intersectAndNottoo_close(x,y,xs(),ys(),num_points,spacing)Thenxs(num_points)=xys(num_points)=ynum_points+=1EndIfEndIfEndSubSubfind_intersections(PAsspower_curve,QAsspower_curve,_flatness_toleranceAsDouble,point_spacingAsDouble,_max_pointsAsInteger,xs()AsDouble,ys()AsDouble,_Byrefnum_pointsAsInteger)num_points=0recursion(0,1,0,1,P,Q,flatness_tolerance,point_spacing,_max_points,xs(),ys(),num_points)EndSubDimAsbernstein_splinebPx=Type<bernstein_spline>(-1,0,1)DimAsbernstein_splinebPy=Type<bernstein_spline>(0,10,0)DimAsbernstein_splinebQx=Type<bernstein_spline>(2,-8,2)DimAsbernstein_splinebQy=Type<bernstein_spline>(1,2,3)DimAsspower_splinePx=bernstein_spline_to_spower(bPx)DimAsspower_splinePy=bernstein_spline_to_spower(bPy)DimAsspower_splineQx=bernstein_spline_to_spower(bQx)DimAsspower_splineQy=bernstein_spline_to_spower(bQy)DimAsspower_curveP=Type<spower_curve>(Px,Py)DimAsspower_curveQ=Type<spower_curve>(Qx,Qy)DimAsDoubleflatness_tolerance=0.001DimAsDoublepoint_spacing=0.000001' Max norm minimum spacing.Constmax_pointsAsInteger=10DimAsDoublexs(max_points)DimAsDoubleys(max_points)DimAsIntegernum_pointsfind_intersections(P,Q,flatness_tolerance,point_spacing,_max_points,xs(),ys(),num_points)ForiAsInteger=0Tonum_points-1Print"(";xs(i);", ";ys(i);")"NextiSleep
/* The control points of a planar quadratic Bézier curve form a triangle--called the "control polygon"--that completely contains the curve. Furthermore, the rectangle formed by the minimum and maximum x and y values of the control polygon completely contain the polygon, and therefore also the curve. Thus a simple method for narrowing down where intersections might be is: subdivide both curves until you find "small enough" regions where these rectangles overlap.*/packagemainimport("fmt""log""math")typepointstruct{x,yfloat64}typequadSplinestruct{// Non-parametric spline.c0,c1,c2float64}typequadCurvestruct{// Planar parametric spline.x,yquadSpline}// Subdivision by de Casteljau's algorithm.funcsubdivideQuadSpline(qquadSpline,tfloat64,u,v*quadSpline){s:=1.0-tc0:=q.c0c1:=q.c1c2:=q.c2u.c0=c0v.c2=c2u.c1=s*c0+t*c1v.c1=s*c1+t*c2u.c2=s*u.c1+t*v.c1v.c0=u.c2}funcsubdivideQuadCurve(qquadCurve,tfloat64,u,v*quadCurve){subdivideQuadSpline(q.x,t,&u.x,&v.x)subdivideQuadSpline(q.y,t,&u.y,&v.y)}// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.funcrectsOverlap(xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1float64)bool{return(xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1)}funcmax3(x,y,zfloat64)float64{returnmath.Max(math.Max(x,y),z)}funcmin3(x,y,zfloat64)float64{returnmath.Min(math.Min(x,y),z)}// This accepts the point as an intersection if the boxes are small enough.functestIntersect(p,qquadCurve,tolfloat64,exclude,accept*bool,intersect*point){pxmin:=min3(p.x.c0,p.x.c1,p.x.c2)pymin:=min3(p.y.c0,p.y.c1,p.y.c2)pxmax:=max3(p.x.c0,p.x.c1,p.x.c2)pymax:=max3(p.y.c0,p.y.c1,p.y.c2)qxmin:=min3(q.x.c0,q.x.c1,q.x.c2)qymin:=min3(q.y.c0,q.y.c1,q.y.c2)qxmax:=max3(q.x.c0,q.x.c1,q.x.c2)qymax:=max3(q.y.c0,q.y.c1,q.y.c2)*exclude=true*accept=falseifrectsOverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax){*exclude=falsexmin:=math.Max(pxmin,qxmin)xmax:=math.Min(pxmax,pxmax)ifxmax<xmin{log.Fatalf("Assertion failure: %f < %f\n",xmax,xmin)}ifxmax-xmin<=tol{ymin:=math.Max(pymin,qymin)ymax:=math.Min(pymax,qymax)ifymax<ymin{log.Fatalf("Assertion failure: %f < %f\n",ymax,ymin)}ifymax-ymin<=tol{*accept=trueintersect.x=0.5*xmin+0.5*xmaxintersect.y=0.5*ymin+0.5*ymax}}}}funcseemsToBeDuplicate(intersects[]point,xypoint,spacingfloat64)bool{seemsToBeDup:=falsei:=0for!seemsToBeDup&&i!=len(intersects){pt:=intersects[i]seemsToBeDup=math.Abs(pt.x-xy.x)<spacing&&math.Abs(pt.y-xy.y)<spacingi++}returnseemsToBeDup}funcfindIntersects(p,qquadCurve,tol,spacingfloat64)[]point{varintersects[]pointtypeworksetstruct{p,qquadCurve}workload:=[]workset{workset{p,q}}// Quit looking after having emptied the workload.forlen(workload)>0{l:=len(workload)work:=workload[l-1]workload=workload[0:l-1]varexclude,acceptboolintersect:=point{0,0}testIntersect(work.p,work.q,tol,&exclude,&accept,&intersect)ifaccept{// To avoid detecting the same intersection twice, require some// space between intersections.if!seemsToBeDuplicate(intersects,intersect,spacing){intersects=append(intersects,intersect)}}elseif!exclude{varp0,p1,q0,q1quadCurvesubdivideQuadCurve(work.p,0.5,&p0,&p1)subdivideQuadCurve(work.q,0.5,&q0,&q1)workload=append(workload,workset{p0,q0})workload=append(workload,workset{p0,q1})workload=append(workload,workset{p1,q0})workload=append(workload,workset{p1,q1})}}returnintersects}funcmain(){varp,qquadCurvep.x=quadSpline{-1.0,0.0,1.0}p.y=quadSpline{0.0,10.0,0.0}q.x=quadSpline{2.0,-8.0,2.0}q.y=quadSpline{1.0,2.0,3.0}tol:=0.0000001spacing:=tol*10intersects:=findIntersects(p,q,tol,spacing)for_,intersect:=rangeintersects{fmt.Printf("(% f, %f)\n",intersect.x,intersect.y)}}
( 0.654983, 2.854983)( 0.881025, 1.118975)(-0.681025, 2.681025)(-0.854983, 1.345017)
importjava.util.ArrayList;importjava.util.List;importjava.util.Stack;publicfinalclassBezierCurveIntersection{publicstaticvoidmain(String[]aArgs){QuadCurvevertical=newQuadCurve(newQuadSpline(-1.0,0.0,1.0),newQuadSpline(0.0,10.0,0.0));// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)QuadCurvehorizontal=newQuadCurve(newQuadSpline(2.0,-8.0,2.0),newQuadSpline(1.0,2.0,3.0));// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)System.out.println("The points of intersection are:");List<Point>intersects=findIntersects(vertical,horizontal);for(Pointintersect:intersects){System.out.println(String.format("%s%9.6f%s%9.6f%s","( ",intersect.aX,", ",intersect.aY," )"));}}privatestaticList<Point>findIntersects(QuadCurveaP,QuadCurveaQ){List<Point>result=newArrayList<Point>();Stack<QuadCurve>stack=newStack<QuadCurve>();stack.push(aP);stack.push(aQ);while(!stack.isEmpty()){QuadCurvepp=stack.pop();QuadCurveqq=stack.pop();List<Object>objects=testIntersection(pp,qq);finalbooleanaccepted=(boolean)objects.get(0);finalbooleanexcluded=(boolean)objects.get(1);Pointintersect=(Point)objects.get(2);if(accepted){if(!seemsToBeDuplicate(result,intersect)){result.add(intersect);}}elseif(!excluded){QuadCurvep0=newQuadCurve();QuadCurveq0=newQuadCurve();QuadCurvep1=newQuadCurve();QuadCurveq1=newQuadCurve();subdivideQuadCurve(pp,0.5,p0,p1);subdivideQuadCurve(qq,0.5,q0,q1);stack.addAll(List.of(p0,q0,p0,q1,p1,q0,p1,q1));}}returnresult;}privatestaticbooleanseemsToBeDuplicate(List<Point>aIntersects,PointaPoint){for(Pointintersect:aIntersects){if(Math.abs(intersect.aX-aPoint.aX)<SPACING&&Math.abs(intersect.aY-aPoint.aY)<SPACING){returntrue;}}returnfalse;}privatestaticList<Object>testIntersection(QuadCurveaP,QuadCurveaQ){finaldoublepxMin=Math.min(Math.min(aP.x.c0,aP.x.c1),aP.x.c2);finaldoublepyMin=Math.min(Math.min(aP.y.c0,aP.y.c1),aP.y.c2);finaldoublepxMax=Math.max(Math.max(aP.x.c0,aP.x.c1),aP.x.c2);finaldoublepyMax=Math.max(Math.max(aP.y.c0,aP.y.c1),aP.y.c2);finaldoubleqxMin=Math.min(Math.min(aQ.x.c0,aQ.x.c1),aQ.x.c2);finaldoubleqyMin=Math.min(Math.min(aQ.y.c0,aQ.y.c1),aQ.y.c2);finaldoubleqxMax=Math.max(Math.max(aQ.x.c0,aQ.x.c1),aQ.x.c2);finaldoubleqyMax=Math.max(Math.max(aQ.y.c0,aQ.y.c1),aQ.y.c2);booleanaccepted=false;booleanexcluded=true;Pointintersect=newPoint(0.0,0.0);if(rectanglesOverlap(pxMin,pyMin,pxMax,pyMax,qxMin,qyMin,qxMax,qyMax)){excluded=false;finaldoublexMin=Math.max(pxMin,qxMin);finaldoublexMax=Math.min(pxMax,pxMax);if(xMax-xMin<=TOLERANCE){finaldoubleyMin=Math.max(pyMin,qyMin);finaldoubleyMax=Math.min(pyMax,qyMax);if(yMax-yMin<=TOLERANCE){accepted=true;intersect=newPoint(0.5*(xMin+xMax),0.5*(yMin+yMax));}}}returnList.of(accepted,excluded,intersect);}privatestaticbooleanrectanglesOverlap(doubleaXa0,doubleaYa0,doubleaXa1,doubleaYa1,doubleaXb0,doubleaYb0,doubleaXb1,doubleaYb1){returnaXb0<=aXa1&&aXa0<=aXb1&&aYb0<=aYa1&&aYa0<=aYb1;}privatestaticvoidsubdivideQuadCurve(QuadCurveaQ,doubleaT,QuadCurveaU,QuadCurveaV){subdivideQuadSpline(aQ.x,aT,aU.x,aV.x);subdivideQuadSpline(aQ.y,aT,aU.y,aV.y);}// de Casteljau's algorithmprivatestaticvoidsubdivideQuadSpline(QuadSplineaQ,doubleaT,QuadSplineaU,QuadSplineaV){finaldoubles=1.0-aT;aU.c0=aQ.c0;aV.c2=aQ.c2;aU.c1=s*aQ.c0+aT*aQ.c1;aV.c1=s*aQ.c1+aT*aQ.c2;aU.c2=s*aU.c1+aT*aV.c1;aV.c0=aU.c2;}privatestaticrecordPoint(doubleaX,doubleaY){}privatestaticclassQuadSpline{publicQuadSpline(doubleaC0,doubleaC1,doubleaC2){c0=aC0;c1=aC1;c2=aC2;}publicQuadSpline(){this(0.0,0.0,0.0);}privatedoublec0,c1,c2;}privatestaticclassQuadCurve{publicQuadCurve(QuadSplineaX,QuadSplineaY){x=aX;y=aY;}publicQuadCurve(){this(newQuadSpline(),newQuadSpline());}privateQuadSplinex,y;}privatestaticfinaldoubleTOLERANCE=0.000_000_1;privatestaticfinaldoubleSPACING=10*TOLERANCE;}
The points of intersection are:( 0.654983, 2.854983 )( -0.681025, 2.681025 )( 0.881025, 1.118975 )( -0.854983, 1.345017 )
classPoint{constructor(x,y){this.x=x;this.y=y;}}classQuadSpline{constructor(c0,c1,c2){this.c0=c0;this.c1=c1;this.c2=c2;}}classQuadCurve{constructor(x,y){this.x=x;this.y=y;}}// Subdivision by de Casteljau's algorithm.functionsubdivideQuadSpline(q,t,u,v){consts=1.0-t;constc0=q.c0;constc1=q.c1;constc2=q.c2;u.c0=c0;v.c2=c2;u.c1=s*c0+t*c1;v.c1=s*c1+t*c2;u.c2=s*u.c1+t*v.c1;v.c0=u.c2;}functionsubdivideQuadCurve(q,t,u,v){subdivideQuadSpline(q.x,t,u.x,v.x);subdivideQuadSpline(q.y,t,u.y,v.y);}// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.functionrectsOverlap(xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1){return(xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1);}functionmax3(x,y,z){returnMath.max(Math.max(x,y),z);}functionmin3(x,y,z){returnMath.min(Math.min(x,y),z);}// This accepts the point as an intersection if the boxes are small enough.functiontestIntersect(p,q,tol){constpxmin=min3(p.x.c0,p.x.c1,p.x.c2);constpymin=min3(p.y.c0,p.y.c1,p.y.c2);constpxmax=max3(p.x.c0,p.x.c1,p.x.c2);constpymax=max3(p.y.c0,p.y.c1,p.y.c2);constqxmin=min3(q.x.c0,q.x.c1,q.x.c2);constqymin=min3(q.y.c0,q.y.c1,q.y.c2);constqxmax=max3(q.x.c0,q.x.c1,q.x.c2);constqymax=max3(q.y.c0,q.y.c1,q.y.c2);letexclude=true;letaccept=false;letintersect=newPoint(0,0);if(rectsOverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)){exclude=false;constxmin=Math.max(pxmin,qxmin);constxmax=Math.min(pxmax,qxmax);if(xmax<xmin){thrownewError(`Assertion failure:${xmax} <${xmin}`);}if(xmax-xmin<=tol){constymin=Math.max(pymin,qymin);constymax=Math.min(pymax,qymax);if(ymax<ymin){thrownewError(`Assertion failure:${ymax} <${ymin}`);}if(ymax-ymin<=tol){accept=true;intersect.x=0.5*xmin+0.5*xmax;intersect.y=0.5*ymin+0.5*ymax;}}}return{exclude,accept,intersect};}functionseemsToBeDuplicate(intersects,xy,spacing){letseemsToBeDup=false;leti=0;while(!seemsToBeDup&&i!==intersects.length){constpt=intersects[i];seemsToBeDup=Math.abs(pt.x-xy.x)<spacing&&Math.abs(pt.y-xy.y)<spacing;i++;}returnseemsToBeDup;}functionfindIntersects(p,q,tol,spacing){constintersects=[];constworkload=[{p,q}];// Quit looking after having emptied the workload.while(workload.length>0){constwork=workload.pop();const{exclude,accept,intersect}=testIntersect(work.p,work.q,tol);if(accept){// To avoid detecting the same intersection twice, require some// space between intersections.if(!seemsToBeDuplicate(intersects,intersect,spacing)){intersects.push(intersect);}}elseif(!exclude){constp0=newQuadCurve(newQuadSpline(0,0,0),newQuadSpline(0,0,0));constp1=newQuadCurve(newQuadSpline(0,0,0),newQuadSpline(0,0,0));constq0=newQuadCurve(newQuadSpline(0,0,0),newQuadSpline(0,0,0));constq1=newQuadCurve(newQuadSpline(0,0,0),newQuadSpline(0,0,0));subdivideQuadCurve(work.p,0.5,p0,p1);subdivideQuadCurve(work.q,0.5,q0,q1);workload.push({p:p1,q:q1});workload.push({p:p1,q:q0});workload.push({p:p0,q:q1});workload.push({p:p0,q:q0});}}returnintersects;}functionmain(){constp=newQuadCurve(newQuadSpline(-1.0,0.0,1.0),newQuadSpline(0.0,10.0,0.0));constq=newQuadCurve(newQuadSpline(2.0,-8.0,2.0),newQuadSpline(1.0,2.0,3.0));consttol=0.0000001;constspacing=tol*10;constintersects=findIntersects(p,q,tol,spacing);for(constintersectofintersects){console.log(`(${intersect.x},${intersect.y})`);}}main();
(-0.8549834191799164, 1.3450165688991547)(-0.6810249984264374, 2.681024934786043)(0.8810248545815753, 1.1189751454184247)(0.6549834311008453, 2.8549834191799164)
This implementation combines rectangle overlap (as in the Modula-2) and curve flattening (as in one of the C implementations), and tries to arrange the calculations efficiently. I think it may be a new algorithm worth serious consideration. Curve flattening alone requires too many line-segment-intersection checks, but rectangle overlap "prunes the tree".
Furthermore, the algorithm returns-parameter pairs, which is ideal.
Icon is a very good language in which to express this algorithm,if the audience can read Icon. But Icon is not designed to be readable by the programming public. See instead thePascal. Pascal is meant to be readable. There also instead of a "workload" set there is a recursive procedure, which probably more programmers will recognize as "how to do adaptive bisection". (But both methods are common.)
See also:
# This program combines the methods of the 2nd C implementation (which# by itself is inefficient) with those of the Modula-2 implementation,# and then rearranges the computations to try to achieve greater# efficiency.## The algorithm actually returns t-parameters for two curves, as a# pair for each intersection point. This is exactly what one might# want: for instance, to break a font glyph, made from two or more# other glyphs, into pieces at the points of intersection of all the# outlines.## The code below is written to illustrate the algorithm rather than to# squeeze performance out of Icon. For instance, I use a "set" to# store the workload, and, when choosing the next workitem-pair to# work on, do so by random selection. It would be faster, certainly,# to use an Icon "list", as either a stack or a queue or both.## It is also possible, of course, to cast the algorithm as a recursive# procedure.### References on the s-power basis:## J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power# basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997,# page 319.## J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis# in geometry processing’, ACM Transactions on Graphics, vol 19# no 1, January 2000, page 35.#recordpoint(x,y)recordspower(c0,c1,c2)recordcurve(x,y)recordworkitem(P,t0,t1,pt0,pt1)$defineP_controls[point(-1,0),point(0,10),point(1,0)]$defineQ_controls[point(2,1),point(-8,2),point(2,3)]$defineDEFAULT_NUMPIECES2# Bisection.# Tolerance of the ratio of a bound on the non-linear component to the# length of the segment. I use a max norm but you can use your# favorite norm.$defineDEFAULT_FLATNESS_TOLERANCE0.0001# For demonstration, I choose a minimum spacing between intersection# points equal to several times single precision machine epsilon. I# measure distance using a max norm, but you can use your favorite# norm.$defineDEFAULT_MINIMUM_SPACING1e-6proceduremain()localP,Qlocalintersections,xsectP:=controls_to_curve!P_controlsQ:=controls_to_curve!Q_controlsintersections:=find_intersections(P,Q)write()write_tabbed_line(" convex up\t"||" convex left\t")everyxsect:=!intersectionsdowrite_tabbed_line(xsect[1]||"\t("||xsect[2].x||", "||xsect[2].y||")\t"||xsect[3]||"\t("||xsect[4].x||", "||xsect[4].y||")")write()endprocedurewrite_tabbed_line(line)write(detab(line,18,56,74))endprocedurefind_intersections(P,Q,tol,spacing)# Return a list of [tP, pointP, tQ, pointQ] for the intersections,# sorted by tP values.localworkload,tslocaltP,ptPlocaltQ,ptQlocaltbl,intersections/tol:=DEFAULT_FLATNESS_TOLERANCE/spacing:=DEFAULT_MINIMUM_SPACINGworkload:=initial_workload(P,Q)tbl:=table()everyts:=process_workload(tol,workload)do{tP:=ts[1];ptP:=curve_eval(P,tP)tQ:=ts[2];ptQ:=curve_eval(Q,tQ)not(max(abs((!tbl)[2].x-ptP.x),abs((!tbl)[2].y-ptP.y))<spacing)¬(max(abs((!tbl)[4].x-ptQ.x),abs((!tbl)[4].y-ptQ.y))<spacing)&tbl[tP]:=[tP,ptP,tQ,ptQ]}tbl:=sort(tbl,1)everyput(intersections:=[],(!tbl)[2])returnintersectionsendprocedureprocess_workload(tol,workload)# Generate pairs of t-parameters.localpair,tswhile*workload~=0do{pair:=?workloaddelete(workload,pair)ifrectangles_overlap(pair[1].pt0,pair[1].pt1,pair[2].pt0,pair[2].pt1)then{ifflat_enough(tol,pair[1])then{ifflat_enough(tol,pair[2])then{ifts:=segment_parameters(pair[1].pt0,pair[1].pt1,pair[2].pt0,pair[2].pt1)thensuspend[(1-ts[1])*pair[1].t0+ts[1]*pair[1].t1,(1-ts[2])*pair[2].t0+ts[2]*pair[2].t1]}elseeveryinsert(workload,[pair[1],split_workitem(pair[2])])}else{ifflat_enough(tol,pair[2])theneveryinsert(workload,[split_workitem(pair[1]),pair[2]])elseeveryinsert(workload,[split_workitem(pair[1]),split_workitem(pair[2])])}}}endproceduresplit_workitem(W,num_pieces)# Split a workitem in pieces and generate the pieces.localfraction,t1_t0,ts,pts,i/num_pieces:=DEFAULT_NUMPIECESfraction:=1.0/num_piecest1_t0:=W.t1-W.t0everyput(ts:=[],W.t0+(1tonum_pieces-1)*fraction*t1_t0)everyput(pts:=[],curve_eval(W.P,!ts))ts:=[W.t0]|||ts|||[W.t1]pts:=[W.pt0]|||pts|||[W.pt1]everyi:=1to*pts-1dosuspend(workitem(W.P,ts[i],ts[i+1],pts[i],pts[i+1]))endprocedureinitial_workload(P,Q)# Create an initial workload set, by breaking P and Q at any# critical points.localworkloadeveryinsert(workload:=set(),[break_at_critical_points(P),break_at_critical_points(Q)])returnworkloadendprocedurebreak_at_critical_points(P)# Generate workitems for the curve P, after breaking it at any# critical points.localts,pts,its:=[0]|||sort(curve_critical_points(P))|||[1]everyput(pts:=[],curve_eval(P,!ts))everyi:=1to*pts-1dosuspend(workitem(P,ts[i],ts[i+1],pts[i],pts[i+1]))endprocedureflat_enough(tol,P,t0,t1,pt0,pt1)# Is the [t0,t1] portion of the curve P flat enough to be treated as# a straight line between pt0 and pt1, where pt0 and pt1 are the# endpoints of the portion?localerror,length# Let flat_enough be called this way, where W is a workitem:# flat_enough(tol,W)if/t0then{pt1:=P.pt1pt0:=P.pt0t1:=P.t1t0:=P.t0P:=P.P}# pt0 and pt1 probably have been computed before and saved, but if# necessary they could be computed now:/pt0:=curve_eval(P,t0)/pt1:=curve_eval(P,t1)# The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to# remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached# at t=1/2. That accounts for the 1/4=0.25 in the following, which# uses "max norm" length measurements. (Substitute your favorite# norm.)error:=0.25*max(abs(spower_center_coef(P.x,t0,t1)),abs(spower_center_coef(P.y,t0,t1)))length:=max(abs(pt1.x-pt0.x),abs(pt1.y-pt0.y))((error<=length*tol)&return)|failendprocedurecurve_eval(P,t)# Return the point that lies on the curve P at parameter value t.returnpoint(spower_eval(P.x,t),spower_eval(P.y,t))endprocedurecurve_critical_points(P)# Return a set containing parameter values (values of t) for the# critical points of curve P.localtsts:=set()insert(ts,spower_critical_point(P.x))insert(ts,spower_critical_point(P.y))returntsendprocedurespower_eval(p,t)# Evaluate the s-power spline p at t.return(p.c0+(p.c1*t))*(1-t)+(p.c2*t)endprocedurespower_center_coef(p,t0,t1)# Return the center coefficient for the [t0,t1] portion of the# s-power spline p.if/t1then{t1:=t0[2];t0:=t0[1]}# Allow a list as t0.returnp.c1*((t1-t0-t0)*t1+(t0*t0))endprocedurespower_critical_point(p)# Return t in (0,1) where p is at a critical point, else fail.localtp.c1=0&fail# The spline is linearp.c0=p.c2&return0.5# The spline is "pulse-like".t:=(0.5*(p.c2+p.c1-p.c0))/p.c1# Root of the derivative.0<t<1&returntfailendprocedurerectangles_overlap(ptA0,ptA1,ptB0,ptB1)# Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1)# overlap?localax0,ay0,ax1,ay1localbx0,by0,bx1,by1ax0:=ptA0.x;ax1:=ptA1.xbx0:=ptB0.x;bx1:=ptB1.xifax1<ax0thenax0:=:ax1ifbx1<bx0thenbx0:=:bx1bx1<ax0&failax1<bx0&failay0:=ptA0.y;ay1:=ptA1.yby0:=ptB0.y;by1:=ptB1.yifay1<ay0thenay0:=:ay1ifby1<by0thenby0:=:by1by1<ay0&failay1<by0&failreturnendproceduresegment_parameters(ptA0,ptA1,ptB0,ptB1)# Return the respective [0,1] parameters of line segments# (ptA0,ptA1) and (ptB0,ptB1), for their intersection point. Fail if# there are not such parameters.localax0,ax1,ay0,ay1localbx0,bx1,by0,by1localax1_ax0,ay1_ay0localbx1_bx0,by1_by0localanumer,bnumer,denomlocaltA,tBax0:=ptA0.x;ax1:=ptA1.xay0:=ptA0.y;ay1:=ptA1.ybx0:=ptB0.x;bx1:=ptB1.xby0:=ptB0.y;by1:=ptB1.yax1_ax0:=ax1-ax0ay1_ay0:=ay1-ay0bx1_bx0:=bx1-bx0by1_by0:=by1-by0denom:=(ax1_ax0*by1_by0)-(ay1_ay0*bx1_bx0)anumer:=(bx1_bx0*ay0)-(by1_by0*ax0)+(bx0*by1)-(bx1*by0)tA:=anumer/denom0<=tA<=1|failbnumer:=-((ax1_ax0*by0)-(ay1_ay0*bx0)+(ax0*ay1)-(ax1*ay0))tB:=bnumer/denom0<=tB<=1|failreturn[tA,tB]endprocedurecontrols_to_curve(p0,p1,p2)# Convert control points to a curve in s-power basis.returncurve(spower(p0.x,(2*p1.x)-p0.x-p2.x,p2.x),spower(p0.y,(2*p1.y)-p0.y-p2.y,p2.y))endprocedureabs(x)return(ifx<0then-xelsex)endproceduremax(x,y)return(ifx<ythenyelsex)end
For each estimated intersection point, the program prints out the-parameters and the corresponding values of the curves.
convex up convex left 0.07250828117 (-0.8549834377, 1.345016607) 0.1725082997 (-0.8549837251, 1.345016599)0.1594875309 (-0.6810249382, 2.681025168) 0.8405124691 (-0.681025168, 2.681024938)0.8274917003 (0.6549834005, 2.854983725) 0.9274917188 (0.6549833933, 2.854983438)0.9405124667 (0.8810249334, 1.118975334) 0.05948753331 (0.8810246662, 1.118975067)
#= The control points of a planar quadratic Bézier curve form a triangle--called the "control polygon"--that completely contains the curve. Furthermore, the rectangle formed by the minimum and maximum x and y values of the control polygon completely contain the polygon, and therefore also the curve. Thus a simple method for narrowing down where intersections might be is: subdivide both curves until you find "small enough" regions where these rectangles overlap.=#mutablestructPointx::Float64y::Float64Point()=new(0,0)endmutablestructQuadSpline# Non-parametric spline.c0::Float64c1::Float64c2::Float64QuadSpline(a=0.0,b=0.0,c=0.0)=new(a,b,c)endmutablestructQuadCurve# Planar parametric spline.x::QuadSpliney::QuadSplineQuadCurve(x=QuadSpline(),y=QuadSpline())=new(x,y)endconstWorkset=Tuple{QuadCurve,QuadCurve}""" Subdivision by de Casteljau's algorithm. """functionsubdivide!(q::QuadSpline,t::Float64,u::QuadSpline,v::QuadSpline)s=1.0-tc0=q.c0c1=q.c1c2=q.c2u.c0=c0v.c2=c2u.c1=s*c0+t*c1v.c1=s*c1+t*c2u.c2=s*u.c1+t*v.c1v.c0=u.c2endfunctionsubdivide!(q::QuadCurve,t::Float64,u::QuadCurve,v::QuadCurve)subdivide!(q.x,t,u.x,v.x)subdivide!(q.y,t,u.y,v.y)end""" It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1. """functionrectsoverlap(xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1)returnxb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1end""" This accepts the point as an intersection if the boxes are small enough. """functiontestintersect(p::QuadCurve,q::QuadCurve,tol::Float64,intersect::Point)pxmin=min(p.x.c0,p.x.c1,p.x.c2)pymin=min(p.y.c0,p.y.c1,p.y.c2)pxmax=max(p.x.c0,p.x.c1,p.x.c2)pymax=max(p.y.c0,p.y.c1,p.y.c2)qxmin=min(q.x.c0,q.x.c1,q.x.c2)qymin=min(q.y.c0,q.y.c1,q.y.c2)qxmax=max(q.x.c0,q.x.c1,q.x.c2)qymax=max(q.y.c0,q.y.c1,q.y.c2)exclude=trueaccept=falseifrectsoverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)exclude=falsexmin=max(pxmin,qxmin)xmax=min(pxmax,pxmax)@assertxmin<xmax"Assertion failure:$xmax <$xmin"ifxmax-xmin<=tolymin=max(pymin,qymin)ymax=min(pymax,qymax)@assertymin<ymax"Assertion failure:$ymax <$ymin"ifymax-ymin<=tolaccept=trueintersect.x=0.5*xmin+0.5*xmaxintersect.y=0.5*ymin+0.5*ymaxendendendreturnexclude,acceptend""" return true if there seems to be a duplicate of xy in the intersects `isects` """seemsduplicate(isects,xy,sp)=any(p->abs(p.x-xy.x)<sp&&abs(p.y-xy.y)<sp,isects)""" find intersections of p and q """functionfindintersects(p,q,tol,spacing)intersects=Point[]workload=Workset[(p,q)]while!isempty(workload)work=popfirst!(workload)intersect=Point()exclude,accept=testintersect(first(work),last(work),tol,intersect)ifaccept# To avoid detecting the same intersection twice, require some# space between intersections.if!seemsduplicate(intersects,intersect,spacing)pushfirst!(intersects,intersect)endelseif!excludep0,p1,q0,q1=QuadCurve(),QuadCurve(),QuadCurve(),QuadCurve()subdivide!(first(work),0.5,p0,p1)subdivide!(last(work),0.5,q0,q1)pushfirst!(workload,(p0,q0),(p0,q1),(p1,q0),(p1,q1))endendreturnintersectsend""" test the methods """functiontestintersections()p,q=QuadCurve(),QuadCurve()p.x=QuadSpline(-1.0,0.0,1.0)p.y=QuadSpline(0.0,10.0,0.0)q.x=QuadSpline(2.0,-8.0,2.0)q.y=QuadSpline(1.0,2.0,3.0)tol=0.0000001spacing=tol*10forintersectinfindintersects(p,q,tol,spacing)println("($(intersect.x)$(intersect.y))")endendtestintersections()
(0.6549834311008453 2.8549834191799164)(0.8810249865055084 1.1189750134944916)(-0.6810249984264374 2.681024934786043)(-0.8549834191799164 1.3450165688991547)
ClearAll[BezierFunc]BezierFunc[pts_List,\[Alpha]_]:=Nest[BlockMap[Apply[(1-\[Alpha])#1+\[Alpha]#2&],#,2,1]&,pts,Length[pts]-1][[1]]b1=BezierFunc[{{-1,0},{0,10},{1,0}},s];b2=BezierFunc[{{2,1},{-8,2},{2,3}},t];eqs=Thread[b1==b2];{s,t,b1}/.Solve[%,{s,t}];Grid[N[%]]
0.07250830.172508{-0.854983,1.34502}0.8274920.927492{0.654983,2.85498}0.1594880.840512{-0.681025,2.68102}0.9405120.0594875{0.881025,1.11898}This Maxima batch script finds an implicit equation for one of the curves, plugs the parametric equations of the other curve into the implicit equation, and then solves the resulting quartic equation.
In theory, doing just the above could find an intersection that lies outside the parameter range of the implicitized curve. But we know all four points lie in that range. One could doublecheck by reversing the roles of the two curves.
/*The method of implicitization:1. Find an implicit equation for one of the curves, in x and y.2. Plug the parametric equations of the other curve into the implicit equation.3. Find the roots of the resulting polynomial equation in t.4. Plug those roots into the parametric equations of step (2).*//* The Bernstein basis of degree 2. See https://en.wikipedia.org/w/index.php?title=Bernstein_polynomial&oldid=1144565695 */b02(t):=1-2*t+t**2$b12(t):=2*t-2*t**2$b22(t):=t**2$/* The convex-up parabola, with its control points as coefficients of the Bernstein basis. */xu(t):=-b02(t)+b22(t)$yu(t):=10*b12(t)$/* The convex-left parabola, with its control points as coefficients of the Bernstein basis. */xl(t):=2*b02(t)-8*b12(t)+2*b22(t)$yl(t):=b02(t)+2*b12(t)+3*b22(t)$/* One can implicitize the convex-up Bézier curve by computing the resultant of x - xu and y - yu. The method is mentioned at https://en.wikipedia.org/w/index.php?title=Gr%C3%B6bner_basis&oldid=1152603392#Implicitization_of_a_rational_curve although they are describing a more general method that I do not know how to do. Here I combine forming the resultant with plugging in xl(t) and yl(t). */quartic_poly:resultant(xl(t)-xu(tau),yl(t)-yu(tau),tau)$/* Find all the roots of the quartic equation that lie in [0,1]. */roots:ev(realroots(quartic_poly=0),float)$roots:sublist(roots,lambda([item],0<=rhs(t)andrhs(t)<=1))$/* Plug them into xl(t) and yl(t). */fori:1thrulength(roots)doblock(display(expand(xl(roots[i]))),display(expand(yl(roots[i]))))$/* As an afterword, I would like to mention some drawbacks of implicitization. * It cannot find self-intersections. This is a major problem for curves of degree 3 or greater. * It gives you the t-parameter values for only one of the two curves. If you just need t-parameter values for both curves (such as to break them up at intersection points), then you could perform implicitization both ways. But, if you need to know which t corresponds to which, you need more than just implicitization. (A method for finding t from given (x,y), for instance.) * It requires first constructing a polynomial of degree 4, 9, 16, etc., and then finding its roots in [0,1]. There are serious difficulties associated with both of those tasks. */
(%i2) b02(t):=1-2*t+t^2(%i3) b12(t):=2*t-2*t^2(%i4) b22(t):=t^2(%i5) xu(t):=-b02(t)+b22(t)(%i6) yu(t):=10*b12(t)(%i7) xl(t):=2*b02(t)-8*b12(t)+2*b22(t)(%i8) yl(t):=b02(t)+2*b12(t)+3*b22(t)(%i9) quartic_poly:resultant(xl(t)-xu(tau),yl(t)-yu(tau),tau)(%i10) roots:ev(realroots(quartic_poly = 0),float)(%i11) roots:sublist(roots,lambda([item],0 <= rhs(t) and rhs(t) <= 1))(%i12) for i thru length(roots) do block(display(expand(xl(roots[i]))),display(expand(yl(roots[i])))) 2 20 t - 20 t + 2 = 0.8810253968010677 2 t + 1 = 1.1189749836921692 2 20 t - 20 t + 2 = - 0.8549833297165428 2 t + 1 = 1.3450165390968323 2 20 t - 20 t + 2 = - 0.681024960552616 2 t + 1 = 2.681024968624115 2 20 t - 20 t + 2 = 0.6549829805579925 2 t + 1 = 2.854983389377594
Compile with the "-fiso" flag.
This program is similar to theD but, instead of using control points to form rectangles, uses values of the curves to form the rectangles. Instead of subdivision, there is function evaluation.
(* This program does not do any subdivision, but instead takes advantage of monotonicity. It is possible for points accidentally to be counted twice, for instance if they lie right on an interval boundary. We will avoid that by the crude (but likely satisfactory) mechanism of requiring a minimum max norm between intersections. *)MODULEbezierIntersectionsInModula2;(* ISO Modula-2 libraries. *)FROMStorageIMPORTALLOCATE,DEALLOCATE;FROMSYSTEMIMPORTTSIZE;IMPORTSLongIO;IMPORTSTextIO;(* GNU Modula-2 gm2-libs *)FROMAssertionIMPORTAssert;(* Schumaker's and Volk's algorithm for evaluating a Bézier spline in Bernstein basis. This is faster than de Casteljau, though not quite as numerical stable. *)PROCEDURESchumakerVolk(c0,c1,c2,t:LONGREAL):LONGREAL;VARs,u,v:LONGREAL;BEGINs:=1.0-t;IFt<=0.5THEN(* Horner form in the variable u = t/s, taking into account the binomial coefficients = 1,2,1. *)u:=t/s;v:=c0+(u*(c1+c1+(u*c2)));(* Multiply by s raised to the degree of the spline. *)v:=v*s*s;ELSE(* Horner form in the variable u = s/t, taking into account the binomial coefficients = 1,2,1. *)u:=s/t;v:=c2+(u*(c1+c1+(u*c0)));(* Multiply by t raised to the degree of the spline. *)v:=v*t*t;END;RETURNv;ENDSchumakerVolk;PROCEDUREFindExtremePoint(c0,c1,c2:LONGREAL;VARLiesInside01:BOOLEAN;VARExtremePoint:LONGREAL);VARnumer,denom:LONGREAL;BEGIN(* If the spline has c0=c2 but not c0=c1=c2, then treat it as having an extreme point at 0.5. *)IF(c0=c2)AND(c0<>c1)THENLiesInside01:=TRUE;ExtremePoint:=0.5ELSE(* Find the root of the derivative of the spline. *)LiesInside01:=FALSE;numer:=c0-c1;denom:=c2-c1-c1+c0;IF(denom<>0.0)AND(numer*denom>=0.0)AND(numer<=denom)THENLiesInside01:=TRUE;ExtremePoint:=numer/denomENDENDENDFindExtremePoint;TYPEStartIntervCount=[2..4];StartIntervArray=ARRAY[1..4]OFLONGREAL;PROCEDUREPossiblyInsertExtremePoint(c0,c1,c2:LONGREAL;VARnumStartInterv:StartIntervCount;VARstartInterv:StartIntervArray);VARliesInside01:BOOLEAN;extremePt:LONGREAL;BEGINFindExtremePoint(c0,c1,c2,liesInside01,extremePt);IFliesInside01AND(0.0<extremePt)AND(extremePt<1.0)THENIFnumStartInterv=2THENstartInterv[3]:=1.0;startInterv[2]:=extremePt;numStartInterv:=3ELSIFextremePt<startInterv[2]THENstartInterv[4]:=1.0;startInterv[3]:=startInterv[2];startInterv[2]:=extremePt;numStartInterv:=4ELSIFextremePt>startInterv[2]THENstartInterv[4]:=1.0;startInterv[3]:=extremePt;numStartInterv:=4ENDENDENDPossiblyInsertExtremePoint;PROCEDUREminimum2(x,y:LONGREAL):LONGREAL;VARw:LONGREAL;BEGINIFx<=yTHENw:=xELSEw:=yEND;RETURNw;ENDminimum2;PROCEDUREmaximum2(x,y:LONGREAL):LONGREAL;VARw:LONGREAL;BEGINIFx>=yTHENw:=xELSEw:=yEND;RETURNw;ENDmaximum2;PROCEDURERectanglesOverlap(xa0,ya0,xa1,ya1:LONGREAL;xb0,yb0,xb1,yb1:LONGREAL):BOOLEAN;BEGIN(* It is assumed that xa0<=xa1, ya0<=ya1, xb0<=xb1, and yb0<=yb1. *)RETURN((xb0<=xa1)AND(xa0<=xb1)AND(yb0<=ya1)AND(ya0<=yb1))ENDRectanglesOverlap;PROCEDURETestIntersection(xp0,xp1:LONGREAL;yp0,yp1:LONGREAL;xq0,xq1:LONGREAL;yq0,yq1:LONGREAL;tol:LONGREAL;VARexclude,accept:BOOLEAN;VARx,y:LONGREAL);VARxpmin,ypmin,xpmax,ypmax:LONGREAL;xqmin,yqmin,xqmax,yqmax:LONGREAL;xmin,xmax,ymin,ymax:LONGREAL;BEGINxpmin:=minimum2(xp0,xp1);ypmin:=minimum2(yp0,yp1);xpmax:=maximum2(xp0,xp1);ypmax:=maximum2(yp0,yp1);xqmin:=minimum2(xq0,xq1);yqmin:=minimum2(yq0,yq1);xqmax:=maximum2(xq0,xq1);yqmax:=maximum2(yq0,yq1);exclude:=TRUE;accept:=FALSE;IFRectanglesOverlap(xpmin,ypmin,xpmax,ypmax,xqmin,yqmin,xqmax,yqmax)THENexclude:=FALSE;xmin:=maximum2(xpmin,xqmin);xmax:=minimum2(xpmax,xqmax);Assert(xmax>=xmin);IFxmax-xmin<=tolTHENymin:=maximum2(ypmin,yqmin);ymax:=minimum2(ypmax,yqmax);Assert(ymax>=ymin);IFymax-ymin<=tolTHENaccept:=TRUE;x:=(0.5*xmin)+(0.5*xmax);y:=(0.5*ymin)+(0.5*ymax);END;END;END;ENDTestIntersection;TYPEWorkPile=POINTERTOWorkTask;WorkTask=RECORDtp0,tp1:LONGREAL;tq0,tq1:LONGREAL;next:WorkPileEND;PROCEDUREWorkIsDone(workload:WorkPile):BOOLEAN;BEGINRETURNworkload=NILENDWorkIsDone;PROCEDUREDeferWork(VARworkload:WorkPile;tp0,tp1:LONGREAL;tq0,tq1:LONGREAL);VARwork:WorkPile;BEGINALLOCATE(work,TSIZE(WorkTask));work^.tp0:=tp0;work^.tp1:=tp1;work^.tq0:=tq0;work^.tq1:=tq1;work^.next:=workload;workload:=workENDDeferWork;PROCEDUREDoSomeWork(VARworkload:WorkPile;VARtp0,tp1:LONGREAL;VARtq0,tq1:LONGREAL);VARwork:WorkPile;BEGINAssert(NOTWorkIsDone(workload));work:=workload;tp0:=work^.tp0;tp1:=work^.tp1;tq0:=work^.tq0;tq1:=work^.tq1;workload:=work^.next;DEALLOCATE(work,TSIZE(WorkTask));ENDDoSomeWork;CONSTpx0=-1.0;px1=0.0;px2=1.0;py0=0.0;py1=10.0;py2=0.0;qx0=2.0;qx1=-8.0;qx2=2.0;qy0=1.0;qy1=2.0;qy2=3.0;tol=0.0000001;spacing=100.0*tol;TYPEIntersectionCount=[0..4];IntersectionRange=[1..4];VARpxHasExtremePt,pyHasExtremePt:BOOLEAN;qxHasExtremePt,qyHasExtremePt:BOOLEAN;pxExtremePt,pyExtremePt:LONGREAL;qxExtremePt,qyExtremePt:LONGREAL;pNumStartInterv,qNumStartInterv:StartIntervCount;pStartInterv,qStartInterv:StartIntervArray;workload:WorkPile;i,j:StartIntervCount;numIntersections,k:IntersectionCount;intersectionsX:ARRAYIntersectionRangeOFLONGREAL;intersectionsY:ARRAYIntersectionRangeOFLONGREAL;tp0,tp1,tq0,tq1:LONGREAL;xp0,xp1,xq0,xq1:LONGREAL;yp0,yp1,yq0,yq1:LONGREAL;exclude,accept:BOOLEAN;x,y:LONGREAL;tpMiddle,tqMiddle:LONGREAL;PROCEDUREMaybeAddIntersection(x,y:LONGREAL;spacing:LONGREAL);VARi:IntersectionRange;VARTooClose:BOOLEAN;BEGINIFnumIntersections=0THENintersectionsX[1]:=x;intersectionsY[1]:=y;numIntersections:=1;ELSETooClose:=FALSE;FORi:=1TOnumIntersectionsDOIF(ABS(x-intersectionsX[i])<spacing)AND(ABS(y-intersectionsY[i])<spacing)THENTooClose:=TRUEENDEND;IFNOTTooCloseTHENnumIntersections:=numIntersections+1;intersectionsX[numIntersections]:=x;intersectionsY[numIntersections]:=yENDENDENDMaybeAddIntersection;BEGIN(* Find monotonic sections of the curves, and use those as the starting jobs. *)pNumStartInterv:=2;pStartInterv[1]:=0.0;pStartInterv[2]:=1.0;PossiblyInsertExtremePoint(px0,px1,px2,pNumStartInterv,pStartInterv);PossiblyInsertExtremePoint(py0,py1,py2,pNumStartInterv,pStartInterv);qNumStartInterv:=2;qStartInterv[1]:=0.0;qStartInterv[2]:=1.0;PossiblyInsertExtremePoint(qx0,qx1,qx2,qNumStartInterv,qStartInterv);PossiblyInsertExtremePoint(qy0,qy1,qy2,qNumStartInterv,qStartInterv);workload:=NIL;FORi:=2TOpNumStartIntervDOFORj:=2TOqNumStartIntervDODeferWork(workload,pStartInterv[i-1],pStartInterv[i],qStartInterv[j-1],qStartInterv[j])END;END;(* Go through the workload, deferring work as necessary. *)numIntersections:=0;WHILENOTWorkIsDone(workload)DO(* The following code recomputes values of the splines sometimes. You may wish to store such values in the work pile, to avoid recomputing them. *)DoSomeWork(workload,tp0,tp1,tq0,tq1);xp0:=SchumakerVolk(px0,px1,px2,tp0);yp0:=SchumakerVolk(py0,py1,py2,tp0);xp1:=SchumakerVolk(px0,px1,px2,tp1);yp1:=SchumakerVolk(py0,py1,py2,tp1);xq0:=SchumakerVolk(qx0,qx1,qx2,tq0);yq0:=SchumakerVolk(qy0,qy1,qy2,tq0);xq1:=SchumakerVolk(qx0,qx1,qx2,tq1);yq1:=SchumakerVolk(qy0,qy1,qy2,tq1);TestIntersection(xp0,xp1,yp0,yp1,xq0,xq1,yq0,yq1,tol,exclude,accept,x,y);IFacceptTHENMaybeAddIntersection(x,y,spacing)ELSIFNOTexcludeTHENtpMiddle:=(0.5*tp0)+(0.5*tp1);tqMiddle:=(0.5*tq0)+(0.5*tq1);DeferWork(workload,tp0,tpMiddle,tq0,tqMiddle);DeferWork(workload,tp0,tpMiddle,tqMiddle,tq1);DeferWork(workload,tpMiddle,tp1,tq0,tqMiddle);DeferWork(workload,tpMiddle,tp1,tqMiddle,tq1);ENDEND;IFnumIntersections=0THENSTextIO.WriteString("no intersections");STextIO.WriteLn;ELSEFORk:=1TOnumIntersectionsDOSTextIO.WriteString("(");SLongIO.WriteReal(intersectionsX[k],10);STextIO.WriteString(", ");SLongIO.WriteReal(intersectionsY[k],10);STextIO.WriteString(")");STextIO.WriteLn;ENDENDENDbezierIntersectionsInModula2.
(0.65498343, 2.85498342)(0.88102499, 1.11897501)(-0.6810249, 2.68102500)(-0.8549834, 1.34501657)
importstd/strformattypePoint=tuple[x,y:float]QuadSpline=tuple[c0,c1,c2:float]# Non-parametric spline.QuadCurve=tuple[x,y:QuadSpline]# Planar parametric spline.funcsubdivideQuadSpline(q:QuadSpline;t:float;u,v:varQuadSpline)=## Subdivision by de Casteljau's algorithm.lets=1.0-tu.c0=q.c0v.c2=q.c2u.c1=s*q.c0+t*q.c1v.c1=s*q.c1+t*q.c2u.c2=s*u.c1+t*v.c1v.c0=u.c2funcsubdivideQuadCurve(q:QuadCurve;t:float;u,v:varQuadCurve)=subdivideQuadSpline(q.x,t,u.x,v.x)subdivideQuadSpline(q.y,t,u.y,v.y)# It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.funcrectsOverlap(xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1:float):bool=xb0<=xa1andxa0<=xb1andyb0<=ya1andya0<=yb1funcmax(x,y,z:float):float=max(max(x,y),z)funcmin(x,y,z:float):float=min(min(x,y),z)# This accepts the point as an intersection if the boxes are small enough.functestIntersect(p,q:QuadCurve;tol:float;exclude,accept:varbool;intersect:varPoint)=letpxmin=min(p.x.c0,p.x.c1,p.x.c2)pymin=min(p.y.c0,p.y.c1,p.y.c2)pxmax=max(p.x.c0,p.x.c1,p.x.c2)pymax=max(p.y.c0,p.y.c1,p.y.c2)qxmin=min(q.x.c0,q.x.c1,q.x.c2)qymin=min(q.y.c0,q.y.c1,q.y.c2)qxmax=max(q.x.c0,q.x.c1,q.x.c2)qymax=max(q.y.c0,q.y.c1,q.y.c2)exclude=trueaccept=falseifrectsOverlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax):exclude=falseletxmin=max(pxmin,qxmin)letxmax=min(pxmax,pxmax)assertxmin<=xmax,&"Assertion failure: {xmin} <= {xmax}"ifxmax-xmin<=tol:letymin=max(pymin,qymin)letymax=min(pymax,qymax)assertymin<=ymax,&"Assertion failure: {ymin} <= {ymax}"ifymax-ymin<=tol:accept=trueintersect=(0.5*xmin+0.5*xmax,0.5*ymin+0.5*ymax)funcseemsToBeDuplicate(intersects:openArray[Point];xy:Point;spacing:float):bool=vari=0whilenotresultandi!=intersects.len:letpt=intersects[i]result=abs(pt.x-xy.x)<spacingandabs(pt.y-xy.y)<spacingincifuncfindIntersects(p,q:QuadCurve;tol,spacing:float):seq[Point]=varworkload=@[(p:p,q:q)]# Quit looking after having emptied the workload.whileworkload.len>0:varexclude,accept:boolvarintersect:Pointletwork=workload.pop()testIntersect(work.p,work.q,tol,exclude,accept,intersect)ifaccept:# To avoid detecting the same intersection twice, require some# space between intersections.ifnotseemsToBeDuplicate(result,intersect,spacing):result.addintersectelifnotexclude:varp0,p1,q0,q1:QuadCurvesubdivideQuadCurve(work.p,0.5,p0,p1)subdivideQuadCurve(work.q,0.5,q0,q1)workload.add@[(p0,q0),(p0,q1),(p1,q0),(p1,q1)]varp,q:QuadCurvep.x=(-1.0,0.0,1.0)p.y=(0.0,10.0,0.0)q.x=(2.0,-8.0,2.0)q.y=(1.0,2.0,3.0)constTol=0.0000001constSpacing=Tol*10letintersects=findIntersects(p,q,Tol,Spacing)forintersectinintersects:echo&"({intersect.x: .6f}, {intersect.y: .6f})"
( 0.654983, 2.854983)( 0.881025, 1.118975)(-0.681025, 2.681025)(-0.854983, 1.345017)
This program is mostly a translation of the Python.
Note that that the Icon program could itself be run as an Object Icon program, after a few necessary modifications. ("import io" would have to be added, the implementations of "abs" and "max" would have to be removed, and maybe a few other minor changes would be needed.)
See also:
# This is the algorithm of the Icon and Python implementations.# Primarily it is a translation of the Python. Snippets are taken from# the Icon.## References on the symmetric power basis:## J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power# basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997,# page 319.## J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis# in geometry processing’, ACM Transactions on Graphics, vol 19# no 1, January 2000, page 35.#import ipl.printf(printf)procedure main () local flatness_tolerance local minimum_spacing local p, q local tbl, params, tp, tq, ppoint, qpoint local keys, k flatness_tolerance := 0.0001 minimum_spacing := 0.000001 p := Curve.from_controls (Point (-1.0, 0.0), Point (0.0, 10.0), Point (1.0, 0.0)) q := Curve.from_controls (Point (2.0, 1.0), Point (-8.0, 2.0), Point ( 2.0, 3.0)) tbl := table () every params := find_intersections (p, q, flatness_tolerance) do { tp := params[1]; ppoint := p.eval(tp) tq := params[2]; qpoint := q.eval(tq) not (length ((!tbl)[2].x - ppoint.x, (!tbl)[2].y - ppoint.y) < minimum_spacing) & not (length ((!tbl)[4].x - qpoint.x, (!tbl)[4].y - qpoint.y) < minimum_spacing) & tbl[tp] := [tp, ppoint, tq, qpoint] } every put (keys := [], key (tbl)) keys := sort (keys) printf ("\n") printf (" convex up " || " convex left\n") every k := !keys do { tp := tbl[k][1] ppoint := tbl[k][2] tq := tbl[k][3] qpoint := tbl[k][4] printf (" %11.8r (%11.8r, %11.8r) " || "%11.8r (%11.8r, %11.8r)\n", tp, ppoint.x, ppoint.y, tq, qpoint.x, qpoint.y) } printf ("\n")end# Generate t-parameter pairs detected as corresponding to intersection# points of p and q. There may be duplicate detections. It is assumed# those will be weeded out by later processing. The tol parameter# specifies the "flatness tolerance" and is a relative tolerance.procedure find_intersections (p, q, tol) local i, j local tp, tq, pportion, qportion local workload, workpair, params # The initial workload is the cartesian product of the monotonic # portions of p and q, respectively. tp := [0.0] ||| sort (p.critical_points()) ||| [1.0] tq := [0.0] ||| sort (q.critical_points()) ||| [1.0] workload := set () every i := 1 to *tp - 1 do every j := 1 to *tq - 1 do { pportion := Portion (p, tp[i], tp[i + 1], p.eval(tp[i]), p.eval(tp[i + 1])) qportion := Portion (q, tq[j], tq[j + 1], q.eval(tq[j]), q.eval(tq[j + 1])) insert (workload, [pportion, qportion]) } while *workload ~= 0 do { workpair := ?workload delete (workload, workpair) pportion := workpair[1] qportion := workpair[2] if rectangles_overlap (pportion.endpt0, pportion.endpt1, qportion.endpt0, qportion.endpt1) then { if pportion.flat_enough(tol) then { if qportion.flat_enough(tol) then { if params := segment_parameters(pportion.endpt0, pportion.endpt1, qportion.endpt0, qportion.endpt1) then { tp := params[1] tq := params[2] tp := (1 - tp) * pportion.t0 + tp * pportion.t1 tq := (1 - tq) * qportion.t0 + tq * qportion.t1 suspend [tp, tq] } } else every insert (workload, [pportion, qportion.split()]) } else { if qportion.flat_enough(tol) then every insert (workload, [pportion.split(), qportion]) else every insert (workload, [pportion.split(), qportion.split()]) } } }endclass Point () public const x, y public new (x, y) self.x := x self.y := y return endendclass SPower () # Non-parametric spline in s-power basis. public const c0, c1, c2 private critpoints public new (c0, c1, c2) self.c0 := c0 self.c1 := c1 self.c2 := c2 return end # Evaluate at t. public eval (t) return (c0 + (c1 * t)) * (1.0 - t) + (c2 * t) end # Return the center coefficient for the [t0,t1] portion. (The other # coefficients can be found with the eval method.) public center_coef (t0, t1) return c1 * ((t1 - t0 - t0) * t1 + (t0 * t0)) end # Return a set of independent variable values for the critical # points that lie in (0,1). (This is a memoizing implementation.) public critical_points () local t if /critpoints then { critpoints := set () if c1 ~= 0.0 then # If c1 is zero then the spline is linear. { if c0 = c2 then insert (critpoints, 0.5) # The spline is "pulse-like". else { # The root of the derivative is the critical point. t = (0.5 * (c2 + c1 - c0)) / c1 0 < t < 1 | fail insert (critpoints, t) } } } return critpoints endendclass Curve () # Parametric spline in s-power basis. public const x, y public new (x, y) self.x := x self.y := y return end public static from_controls (ctl0, ctl1, ctl2) local c0x, c0y, c1x, c1y, c2x, c2y c0x := ctl0.x c0y := ctl0.y c1x := (2.0 * ctl1.x) - ctl0.x - ctl2.x c1y := (2.0 * ctl1.y) - ctl0.y - ctl2.y c2x := ctl2.x c2y := ctl2.y return Curve (SPower (c0x, c1x, c2x), SPower (c0y, c1y, c2y)) end # Evaluate at t. public eval (t) return Point (x.eval(t), y.eval(t)) end # Return a set of t-parameter values for the critical points that # lie in (0,1). public critical_points () return (x.critical_points() ++ y.critical_points()) endendclass Portion () # The [t0,t1]-portion of a parametric spline. public const curve, t0, t1, endpt0, endpt1 public static default_num_pieces private static init () default_num_pieces := 2 end public new (curve, t0, t1, endpt0, endpt1) self.curve := curve self.t0 := t0 self.t1 := t1 self.endpt0 := endpt0 self.endpt1 := endpt1 return end # Is the Portion close enough to linear to be treated as a line # segment? public flat_enough (tol) local xcentercoef, ycentercoef local xlen, ylen # The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to # remove the terms in t(1-t). The maximum of t(1-t) is 1/4, # reached at t=1/2. That accounts for the 1/4=0.25 in the # following. xcentercoef := curve.x.center_coef(t0, t1) ycentercoef := curve.y.center_coef(t0, t1) xlen := endpt1.x - endpt0.x ylen := endpt1.y - endpt0.y return compare_lengths (0.25 * xcentercoef, 0.25 * ycentercoef, tol * xlen, tol * ylen) <= 0 end # Generate num_pieces (or default_num_pieces) sections of the # Portion. public split (num_pieces) local k1, k, ts, vals, i /num_pieces := default_num_pieces k1 := 1.0 / num_pieces every put (ts := [], (k := k1 * (1 to num_pieces - 1) & (1.0 - k) * t0 + k * t1)) every put (vals := [], curve.eval(!ts)) ts := [t0] ||| ts ||| [t1] vals := [endpt0] ||| vals ||| [endpt1] every i := 1 to *ts - 1 do suspend Portion (curve, ts[i], ts[i + 1], vals[i], vals[i + 1]) endend# Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1)# overlap?procedure rectangles_overlap (ptA0, ptA1, ptB0, ptB1) local ax0, ay0, ax1, ay1 local bx0, by0, bx1, by1 ax0 := ptA0.x; ax1 := ptA1.x bx0 := ptB0.x; bx1 := ptB1.x if ax1 < ax0 then ax0 :=: ax1 if bx1 < bx0 then bx0 :=: bx1 bx1 < ax0 & fail ax1 < bx0 & fail ay0 := ptA0.y; ay1 := ptA1.y by0 := ptB0.y; by1 := ptB1.y if ay1 < ay0 then ay0 :=: ay1 if by1 < by0 then by0 :=: by1 by1 < ay0 & fail ay1 < by0 & fail returnend# Return the respective [0,1] parameters of line segments# (ptA0,ptA1) and (ptB0,ptB1), for their intersection point. Fail if# there are not such parameters.procedure segment_parameters (ptA0, ptA1, ptB0, ptB1) local ax0, ax1, ay0, ay1 local bx0, bx1, by0, by1 local axdiff, aydiff local bxdiff, bydiff local anumer, bnumer, denom local tA, tB ax0 := ptA0.x; ax1 := ptA1.x ay0 := ptA0.y; ay1 := ptA1.y bx0 := ptB0.x; bx1 := ptB1.x by0 := ptB0.y; by1 := ptB1.y axdiff := ax1 - ax0 aydiff := ay1 - ay0 bxdiff := bx1 - bx0 bydiff := by1 - by0 denom := (axdiff * bydiff) - (aydiff * bxdiff) anumer := (bxdiff * ay0) - (bydiff * ax0) + (bx0 * by1) - (bx1 * by0) tA := anumer / denom 0 <= tA <= 1 | fail bnumer := -((axdiff * by0) - (aydiff * bx0) + (ax0 * ay1) - (ax1 * ay0)) tB := bnumer / denom 0 <= tB <= 1 | fail return [tA, tB]end# Length according to some norm, where (ax,ay) is a "measuring stick"# vector. Here I use the max norm.procedure length (a_x, a_y) return max (abs (a_x), abs (a_y))end# Having a compare_lengths function lets one compare lengths in the# euclidean metric by comparing the squares of the lengths, and thus# avoiding the square root. The following, however, is a general# implementation.procedure compare_lengths (a_x, a_y, b_x, b_y) local len_a, len_b, cmpval len_a := length (a_x, a_y) len_b := length (b_x, b_y) if len_a < len_b then cmpval := -1 else if len_a > len_b then cmpval := 1 else cmpval := 0 return cmpvalendconvex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
This is the algorithm of theIcon example, recast as a recursive procedure.
See also:
{$mode ISO}{ Tell Free Pascal Compiler to use "ISO 7185" mode. }{This is the algorithm of the Icon example, recast as a recursiveprocedure.The "var" notation in formal parameter lists means pass byreference. All other parameters are implicitly passed by value.Pascal is case-insensitive.In the old days, when Pascal was printed as a means to expressalgorithms, it was usually in a fashion similar to Algol 60 referencelanguage. It was printed mostly in lowercase and did not haveunderscores. Reserved words were in boldface and variables, etc., werein italics. The effect was like that of Algol 60 reference language.Code entry practices for Pascal were another matter. It may have beenall uppercase, with ML-style comment braces instead of squigglybraces. It may have had uppercase reserved words and "Pascal case"variables, etc., as one sees also in Modula-2 and Oberon-2 code.Here I have deliberately adopted an all-lowercase style.References on the s-power basis: J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, page 319. J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis in geometry processing’, ACM Transactions on Graphics, vol 19 no 1, January 2000, page 35.}programbezierintersections;constflatnesstolerance=0.0001;minimumspacing=0.000001;maxintersections=10;typepoint=recordx,y:realend;spower={ non-parametric spline in s-power basis }recordc0,c1,c2:realend;curve={ parametric spline in s-power basis }recordx,y:spowerend;portion={ portion of a parametric spline in [t0,t1] }recordcurv:curve;t0,t1:real;endpt0,endpt1:point{ pre-computed for efficiency }end;intersectionscount=0..maxintersections;intersectionsrange=1..maxintersections;tparamsarray=array[intersectionsrange]ofreal;coordsarray=array[intersectionsrange]ofpoint;varnumintersections:intersectionscount;tparamsp:tparamsarray;coordsp:coordsarray;tparamsq:tparamsarray;coordsq:coordsarray;pglobal,qglobal:curve;i:integer;{ Minimum of two real. }functionrmin(x,y:real):real;beginifx<ythenrmin:=xelsermin:=yend;{ Maximum of two real. }functionrmax(x,y:real):real;beginifx<ythenrmax:=yelsermax:=xend;{ Insertion sort of an array of real. }procedurerealsort(n:integer;vara:arrayofreal);vari,j:integer;x:real;done:boolean;begini:=low(a)+1;whilei<ndobeginx:=a[i];j:=i-1;done:=false;whilenotdonedobeginifj+1=low(a)thendone:=trueelseifa[j]<=xthendone:=trueelsebegina[j+1]:=a[j];j:=j-1endend;a[j+1]:=x;i:=i+1endend;{ "Length" according to some definition. Here I use a max norm. The "distance" between two points is the "length" of the differences of the corresponding coordinates. (The sign of the difference should be immaterial.) }functionlength(ax,ay:real):real;beginlength:=rmax(abs(ax),abs(ay))end;{ Having a "comparelengths" function makes it possible to use a euclidean norm for "length", and yet avoid square roots. One compares the squares of lengths, instead of the lengths themselves. However, here I use a more general implementation. }functioncomparelengths(ax,ay,bx,by:real):integer;varlena,lenb:real;beginlena:=length(ax,ay);lenb:=length(bx,by);iflena<lenbthencomparelengths:=-1elseiflena>lenbthencomparelengths:=1elsecomparelengths:=0end;functionmakepoint(x,y:real):point;beginmakepoint.x:=x;makepoint.y:=yend;functionmakeportion(curv:curve;t0,t1:real;endpt0,endpt1:point):portion;beginmakeportion.curv:=curv;makeportion.t0:=t0;makeportion.t1:=t1;makeportion.endpt0:=endpt0;makeportion.endpt1:=endpt1;end;{ Convert from control points (that is, Bernstein basis) to the symmetric power basis. }functioncontrolstospower(ctl0,ctl1,ctl2:point):curve;begincontrolstospower.x.c0:=ctl0.x;controlstospower.y.c0:=ctl0.y;controlstospower.x.c1:=(2.0*ctl1.x)-ctl0.x-ctl2.x;controlstospower.y.c1:=(2.0*ctl1.y)-ctl0.y-ctl2.y;controlstospower.x.c2:=ctl2.x;controlstospower.y.c2:=ctl2.yend;{ Evaluate an s-power spline at t. }functionspowereval(spow:spower;t:real):real;beginspowereval:=(spow.c0+(spow.c1*t))*(1.0-t)+(spow.c2*t)end;{ Evaluate a curve at t. }functioncurveeval(curv:curve;t:real):point;begincurveeval.x:=spowereval(curv.x,t);curveeval.y:=spowereval(curv.y,t)end;{ Return the center coefficient for the [t0,t1] portion of an s-power spline. (The endpoint coefficients can be found with spowereval.) }functionspowercentercoef(spow:spower;t0,t1:real):real;beginspowercentercoef:=spow.c1*((t1-t0-t0)*t1+(t0*t0))end;{ Return t in (0,1) where spow is at a critical point, else return -1.0. }functionspowercriticalpt(spow:spower):real;vart:real;beginspowercriticalpt:=-1.0;ifspow.c1<>0.0then{ If c1 is zero, then the spline is linear. }beginifspow.c1=spow.c2thenspowercriticalpt:=0.5{ The spline is "pulse-like". }elsebegin{ t = root of the derivative }t:=(spow.c2+spow.c1-spow.c0)/(spow.c1+spow.c1);if(0.0<t)and(t<1.0)thenspowercriticalpt:=tendendend;{ Bisect a portion and pre-compute the new shared endpoint. }procedurebisectportion(port:portion;varport1,port2:portion);beginport1.curv:=port.curv;port2.curv:=port.curv;port1.t0:=port.t0;port1.t1:=0.5*(port.t0+port.t1);port2.t0:=port1.t1;port2.t1:=port.t1;port1.endpt0:=port.endpt0;port1.endpt1:=curveeval(port.curv,port1.t1);port2.endpt0:=port1.endpt1;port2.endpt1:=port.endpt1;end;{ Do the rectangles with corners at (a0,a1) and (b0,b1) overlap at all? }functionrectanglesoverlap(a0,a1,b0,b1:point):boolean;beginrectanglesoverlap:=((rmin(a0.x,a1.x)<=rmax(b0.x,b1.x))and(rmin(b0.x,b1.x)<=rmax(a0.x,a1.x))and(rmin(a0.y,a1.y)<=rmax(b0.y,b1.y))and(rmin(b0.y,b1.y)<=rmax(a0.y,a1.y)))end;{ Set the respective [0,1] parameters of line segments (a0,a1) and (b0,b1), for their intersection point. If there are not two such parameters, set both values to -1.0. }proceduresegmentparameters(a0,a1,b0,b1:point;varta,tb:real);varanumer,bnumer,denom:real;axdiff,aydiff,bxdiff,bydiff:real;beginaxdiff:=a1.x-a0.x;aydiff:=a1.y-a0.y;bxdiff:=b1.x-b0.x;bydiff:=b1.y-b0.y;denom:=(axdiff*bydiff)-(aydiff*bxdiff);anumer:=((bxdiff*a0.y)-(bydiff*a0.x)+(b0.x*b1.y)-(b1.x*b0.y));ta:=anumer/denom;if(ta<0.0)or(1.0<ta)thenbeginta:=-1.0;tb:=-1.0endelsebeginbnumer:=-((axdiff*b0.y)-(aydiff*b0.x)+(a0.x*a1.y)-(a1.x*a0.y));tb:=bnumer/denom;if(tb<0.0)or(1.0<tb)thenbeginta:=-1.0;tb:=-1.0endendend;{ Is a curve portion flat enough to be treated as a line segment between its endpoints? }functionflatenough(port:portion;tol:real):boolean;varxcentercoef,ycentercoef:real;begin{ The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached at t=1/2. That accounts for the 1/4=0.25 in the following. }{ The "with" construct here is a shorthand to implicitly use fields of the "port" record. Thus "curv.x" means "port.curv.x", etc. }withportdobeginxcentercoef:=spowercentercoef(curv.x,t0,t1);ycentercoef:=spowercentercoef(curv.y,t0,t1);flatenough:=comparelengths(0.25*xcentercoef,0.25*ycentercoef,tol*(endpt1.x-endpt0.x),tol*(endpt1.y-endpt0.y))<=0endend;{ If the intersection point corresponding to tp and tq is not already listed, insert it into the arrays, sorted by the value of tp. }procedureinsertintersection(p:curve;tp:real;q:curve;tq:real);varppoint,qpoint:point;lenp,lenq:real;i:intersectionscount;insertionpoint:intersectionscount;beginifnumintersections<>maxintersectionsthenbeginppoint:=curveeval(p,tp);qpoint:=curveeval(q,tq);insertionpoint:=numintersections+1;{ Insert at end. }i:=0;while(0<insertionpoint)and(i<>numintersections)dobegini:=i+1;lenp:=length(coordsp[i].x-ppoint.x,coordsp[i].y-ppoint.y);lenq:=length(coordsq[i].x-qpoint.x,coordsq[i].y-qpoint.y);if(lenp<minimumspacing)and(lenq<minimumspacing)theninsertionpoint:=0{ The point is already listed. }elseiftp<tparamsp[i]thenbegininsertionpoint:=i;{ Insert here instead of at end. }i:=numintersectionsendend;ifinsertionpoint<>numintersections+1thenfori:=numintersections+1downtoinsertionpoint+1dobegintparamsp[i]:=tparamsp[i-1];coordsp[i]:=coordsp[i-1];tparamsq[i]:=tparamsq[i-1];coordsq[i]:=coordsq[i-1]end;tparamsp[insertionpoint]:=tp;coordsp[insertionpoint]:=ppoint;tparamsq[insertionpoint]:=tq;coordsq[insertionpoint]:=qpoint;numintersections:=numintersections+1endend;{ Find intersections between portions of two curves. }procedurefindportionintersections(pportion,qportion:portion);vartp,tq:real;pport1,pport2:portion;qport1,qport2:portion;beginifrectanglesoverlap(pportion.endpt0,pportion.endpt1,qportion.endpt0,qportion.endpt1)thenbeginifflatenough(pportion,flatnesstolerance)thenbeginifflatenough(qportion,flatnesstolerance)thenbeginsegmentparameters(pportion.endpt0,pportion.endpt1,qportion.endpt0,qportion.endpt1,tp,tq);if0.0<=tpthenbegintp:=(1.0-tp)*pportion.t0+tp*pportion.t1;tq:=(1.0-tq)*qportion.t0+tq*qportion.t1;insertintersection(pportion.curv,tp,qportion.curv,tq)endendelsebeginbisectportion(qportion,qport1,qport2);findportionintersections(pportion,qport1);findportionintersections(pportion,qport2)endendelsebeginbisectportion(pportion,pport1,pport2);ifflatenough(qportion,flatnesstolerance)thenbeginfindportionintersections(pport1,qportion);findportionintersections(pport2,qportion)endelsebeginbisectportion(qportion,qport1,qport2);findportionintersections(pport1,qport1);findportionintersections(pport1,qport2);findportionintersections(pport2,qport1);findportionintersections(pport2,qport2)endendendend;{ Find intersections in [0,1]. }procedurefindintersections(p,q:curve);vartpx,tpy,tqx,tqy:real;tp,tq:array[1..4]ofreal;ppoints,qpoints:array[1..4]ofpoint;np,nq,i,j:integer;pportion,qportion:portion;procedurepfindcriticalpts;vari:integer;begintp[1]:=0.0;tp[2]:=1.0;np:=2;tpx:=spowercriticalpt(p.x);tpy:=spowercriticalpt(p.y);if(0.0<tpx)and(tpx<1.0)thenbeginnp:=np+1;tp[np]:=tpxend;if(0.0<tpy)and(tpy<1.0)and(tpy<>tpx)thenbeginnp:=np+1;tp[np]:=tpyend;realsort(np,tp);fori:=1tonpdoppoints[i]:=curveeval(p,tp[i])end;procedureqfindcriticalpts;vari:integer;begintq[1]:=0.0;tq[2]:=1.0;nq:=2;tqx:=spowercriticalpt(q.x);tqy:=spowercriticalpt(q.y);if(0.0<tqx)and(tqx<1.0)thenbeginnq:=nq+1;tq[nq]:=tqxend;if(0.0<tqy)and(tqy<1.0)and(tqy<>tqx)thenbeginnq:=nq+1;tq[nq]:=tqyend;realsort(nq,tq);fori:=1tonqdoqpoints[i]:=curveeval(q,tq[i])end;begin{ Break the curves at critical points, so one can assume the portion between two endpoints is monotonic along both axes. }pfindcriticalpts;qfindcriticalpts;{ Find intersections in the cartesian product of portions of the two curves. (If you would like to compare with the Icon code: In the Icon, goal-directed evaluation is inserting such cartesian products into the "workload" set. However, to do this requires only one "every" construct instead of two, and there is no need for loop/counter variables.) }fori:=1tonp-1doforj:=1tonq-1dobeginpportion:=makeportion(p,tp[i],tp[i+1],ppoints[i],ppoints[i+1]);qportion:=makeportion(q,tq[j],tq[j+1],qpoints[j],qpoints[j+1]);findportionintersections(pportion,qportion);endend;beginpglobal:=controlstospower(makepoint(-1.0,0.0),makepoint(0.0,10.0),makepoint(1.0,0.0));qglobal:=controlstospower(makepoint(2.0,1.0),makepoint(-8.0,2.0),makepoint(2.0,3.0));numintersections:=0;findintersections(pglobal,qglobal);writeln;writeln(' convex up ',' convex left');fori:=1tonumintersectionsdowriteln(' ',tparamsp[i]:11:8,' (',coordsp[i].x:11:8,', ',coordsp[i].y:11:8,') ',tparamsq[i]:11:8,' (',coordsq[i].x:11:8,', ',coordsq[i].y:11:8,')');writelnend.
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
Aside: at long last found my first ever real-world use of sq_atom()... and o/c it had a silly bug!
enumX,Ytypequadratic_spline(sequence/*c*/)-- return apply(c,sq_atom)={1,1,1} -- oops, requires 1.0.3...returntrue-- this will do instead for 1.0.2 and earlierendtypetypequadratic_curve(sequence/*c*/)-- return apply(c,sq_atom)={{1,1,1},{1,1,1}} -- dittoreturntrueendtypefunctionsubdivide_quadratic_spline(quadratic_splineq,atomt)// Subdivision by de Casteljau's algorithm.atom{c1,c2,c3}=q,s=1-t,u1=(s*c1)+(t*c2),v1=(s*c2)+(t*c3),m=(s*u1)+(t*v1)return{{c1,u1,m},{m,v1,c3}}endfunctionfunctionsubdivide_quadratic_curve(quadratic_curveq,atomt)sequence{px,qx}=subdivide_quadratic_spline(q[X],t),{py,qy}=subdivide_quadratic_spline(q[Y],t)return{{px,py},{qx,qy}}endfunctionfunctionrectangles_overlap(atomxa1,ya1,xa2,ya2,xb1,yb1,xb2,yb2)assert(xa1<=xa2andya1<=ya2andxb1<=xb2andyb1<=yb2)returnxb1<=xa2andxa1<=xb2andyb1<=ya2andya1<=yb2endfunctionfunctiontest_intersection(quadratic_curvep,q,atomtolerance)atompxmin=min(p[X]),pymin=min(p[Y]),pxmax=max(p[X]),pymax=max(p[Y]),qxmin=min(q[X]),qymin=min(q[Y]),qxmax=max(q[X]),qymax=max(q[Y])ifrectangles_overlap(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)thenatomxmin=max(pxmin,qxmin),xmax=min(pxmax,qxmax),ymin=max(pymin,qymin),ymax=min(pymax,qymax)assert(xmax>=xminandymax>=ymin)ifxmax-xmin<=toleranceandymax-ymin<=tolerancethen-- we found a suitable intersection!return{(xmin+xmax)/2,(ymin+ymax)/2}endifreturntrue-- accept/further subdivideendifreturnfalse-- excludeendfunctionfunctionseems_to_be_a_duplicate(sequenceintersections,xy,atomspacing)forptinintersectionsdoifabs(pt[X]-xy[X])<spacingandabs(pt[Y]-xy[Y])<spacingthenreturntrueendifendforreturnfalseendfunctionfunctionfind_intersections(quadratic_curvep,q,atomtolerance)sequenceinsects={},todo={{p,q}}whilelength(todo)do{{p,q},todo}={todo[1],todo[2..$]}objectinsect=test_intersection(p,q,tolerance)ifsequence(insect)thenifnotseems_to_be_a_duplicate(insects,insect,tolerance*10)theninsects&={insect}endifelsifinsectthensequence{p1,p2}=subdivide_quadratic_curve(p,0.5),{q1,q2}=subdivide_quadratic_curve(q,0.5)todo&={{p1,q1},{p1,q2},{p2,q1},{p2,q2}}endifendwhileinsects=sort_columns(insects,{-Y,X})returninsectsendfunctionquadratic_curvep={{-1,0,1},{0,10,0}},q={{2,-8,2},{1,2,3}}sequenceintersections=find_intersections(p,q,0.000001)printf(1,"Intersections from top to bottom:\n")pp(intersections,{pp_Nest,1,pp_FltFmt,"%9.6f"})
Intersections from top to bottom:{{ 0.654983, 2.854984}, {-0.681025, 2.681025}, {-0.854984, 1.345017}, { 0.881025, 1.118975}}This is the recursive version of the algorithm that appeared first (in non-recursive form) in theIcon example.
-- -*- mode: indented-text; tab-width: 2; -*-enumX,YenumC0,C1,C2enumCURV,T0,T1,ENDPT0,ENDPT1functionnormlength(atoma_x,a_y)-- A "length" according to our chosen norm (which happens to be the-- max norm). Here (a_x,a_y) is a vector used as measuring rod.returnmax(abs(a_x),abs(a_y))endfunctionfunctioncompare_normlengths(atoma_x,a_y,b_x,b_y)-- Here is a general implementation for comparison of two-- "normlength". For the euclidean norm, you might wish to skip-- taking square roots.atomlen_a=normlength(a_x,a_y),len_b=normlength(b_x,b_y),cmpval=0iflen_a<len_bthencmpval=-1elsiflen_a>len_bthencmpval=1endifreturncmpvalendfunctionfunctioncontrols_to_spower(sequencecontrols)-- Convert from control points (that is, Bernstein basis) to the-- symmetric power basis.sequencept0=controls[1],pt1=controls[2],pt2=controls[3]return{{pt0[X],(2*pt1[X])-pt0[X]-pt2[X],pt2[X]},{pt0[Y],(2*pt1[Y])-pt0[Y]-pt2[Y],pt2[Y]}}endfunctionfunctionspower_eval(sequencespow,atomt)-- Evaluate an s-power spline at t.return(spow[C0]+(spow[C1]*t))*(1-t)+(spow[C2]*t)endfunctionfunctioncurve_eval(sequencecurv,atomt)-- Evaluate a curve at t.return{spower_eval(curv[X],t),spower_eval(curv[Y],t)}endfunctionfunctionspower_center_coef(sequencespow,atomt0,t1)-- Return the center coefficient for the [t0,t1] portion of an-- s-power spline. (The endpoint coefficients can be found with-- spower_eval.) }returnspow[C1]*((t1-t0-t0)*t1+(t0*t0))endfunctionfunctionspower_critical_pt(sequencespow)-- Return t in (0,1) where p is at a critical point, else return -1.atomcrit_pt=-1ifspow[C1]!=0then-- If c1 is zero, then the spline is linear.ifspow[C1]=spow[C2]thencrit_pt=0.5-- The spline is "pulse-like".else-- The critical point is the root of the derivative.atomt=(0.5*(spow[C2]+spow[C1]-spow[C0]))/spow[C1]if(0<t)and(t<1)thencrit_pt=tendifendifendifreturncrit_ptendfunctionfunctionbisect_portion(sequenceport)-- Bisect a portion and pre-compute the new shared endpoint.atomt_mid=0.5*(port[T0]+port[T1])sequencept_mid=curve_eval(port[CURV],t_mid)return{{port[CURV],port[T0],t_mid,port[ENDPT0],pt_mid},{port[CURV],t_mid,port[T1],pt_mid,port[ENDPT1]}}endfunctionfunctionrectangles_overlap(sequencea0,a1,b0,b1)-- Do the rectangles with corners at (a0,a1) and (b0,b1) overlap at-- all?return((min(a0[X],a1[X])<=max(b0[X],b1[X]))and(min(b0[X],b1[X])<=max(a0[X],a1[X]))and(min(a0[Y],a1[Y])<=max(b0[Y],b1[Y]))and(min(b0[Y],b1[Y])<=max(a0[Y],a1[Y])))endfunctionfunctionsegment_parameters(sequencea0,a1,b0,b1)-- Return the respective [0,1] parameters of line segments (a0,a1)-- and (b0,b1), for their intersection point. If there are not two-- such parameters, return -1 for both values.atomaxdiff=a1[X]-a0[X],aydiff=a1[Y]-a0[Y],bxdiff=b1[X]-b0[X],bydiff=b1[Y]-b0[Y],denom=(axdiff*bydiff)-(aydiff*bxdiff),anumer=((bxdiff*a0[Y])-(bydiff*a0[X])+(b0[X]*b1[Y])-(b1[X]*b0[Y])),ta=anumer/denom,tb=-1if(ta<0.0)or(1.0<ta)thenta=-1elseatombnumer=-((axdiff*b0[Y])-(aydiff*b0[X])+(a0[X]*a1[Y])-(a1[X]*a0[Y]))tb=bnumer/denomif(tb<0.0)or(1.0<tb)thenta=-1tb=-1endifendifreturn{ta,tb}endfunctionfunctionflat_enough(sequenceport,atomtol)-- Is a curve portion flat enough to be treated as a line segment-- between its endpoints?-- The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to-- remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached-- at t=1/2. That accounts for the 1/4=0.25 in the following.atomxcentercoef=spower_center_coef(port[CURV][X],port[T0],port[T1]),ycentercoef=spower_center_coef(port[CURV][Y],port[T0],port[T1]),xlen=port[ENDPT1][X]-port[ENDPT0][X],ylen=port[ENDPT1][Y]-port[ENDPT0][Y]return(compare_normlengths(0.25*xcentercoef,0.25*ycentercoef,tol*xlen,tol*ylen)<=0)endfunctionfunctionfind_portion_intersections(sequencepportion,qportion,atomtol)-- Find intersections between portions of two curves. Return pairs-- of t-parameters. There may be duplicates and (due to truncation-- error) near-intersections.sequenceintersections={}ifrectangles_overlap(pportion[ENDPT0],pportion[ENDPT1],qportion[ENDPT0],qportion[ENDPT1])thenifflat_enough(pportion,tol)thenifflat_enough(qportion,tol)thenatomtp,tq{tp,tq}=segment_parameters(pportion[ENDPT0],pportion[ENDPT1],qportion[ENDPT0],qportion[ENDPT1])if0<=tpthentp=(1-tp)*pportion[T0]+tp*pportion[T1]tq=(1-tq)*qportion[T0]+tq*qportion[T1]intersections&={{tp,tq}}endifelsesequenceqport1,qport2{qport1,qport2}=bisect_portion(qportion)intersections&=find_portion_intersections(pportion,qport1,tol)intersections&=find_portion_intersections(pportion,qport2,tol)endifelseifflat_enough(qportion,tol)thensequencepport1,pport2{pport1,pport2}=bisect_portion(pportion)intersections&=find_portion_intersections(pport1,qportion,tol)intersections&=find_portion_intersections(pport2,qportion,tol)elsesequencepport1,pport2sequenceqport1,qport2{pport1,pport2}=bisect_portion(pportion){qport1,qport2}=bisect_portion(qportion)intersections&=find_portion_intersections(pport1,qport1,tol)intersections&=find_portion_intersections(pport1,qport2,tol)intersections&=find_portion_intersections(pport2,qport1,tol)intersections&=find_portion_intersections(pport2,qport2,tol)endifendifendifreturnintersectionsendfunctionfunctionfind_intersections(sequencep,q,atomtol)-- Find intersections in [0,1]. Return pairs of t-parameters.-- There may be duplicates and (due to truncation error)-- near-intersections.-- Break the curves at critical points, so one can assume the-- portion between two endpoints is monotonic along both axes.sequencetp={0,1}atomtpx=spower_critical_pt(p[X]),tpy=spower_critical_pt(p[Y])if0<tpxthentp&={tpx}endifif0<tpyandtpy!=tpxthentp&={tpy}endiftp=sort(tp)sequencetq={0,1}sequencepvals={}fortintpdopvals&={curve_eval(p,t)}endforatomtqx=spower_critical_pt(q[X]),tqy=spower_critical_pt(q[Y])if0<tqxthentq&={tqx}endifif0<tqyandtqy!=tqxthentq&={tqy}endiftq=sort(tq)sequenceqvals={}fortintqdoqvals&={curve_eval(q,t)}endfor-- Find intersections in the cartesian product of the monotonic-- portions of the two curves.sequenceintersections={}fori=1tolength(tp)-1doforj=1tolength(tq)-1dosequencepportion={p,tp[i],tp[i+1],pvals[i],pvals[i+1]},qportion={q,tq[j],tq[j+1],qvals[j],qvals[j+1]}intersections&=find_portion_intersections(pportion,qportion,tol)endforendforreturnintersectionsendfunctionsequencepcontrols={{-1,0},{0,10},{1,0}},qcontrols={{2,1},{-8,2},{2,3}},p=controls_to_spower(pcontrols),q=controls_to_spower(qcontrols)atomflatness_tolerance=0.0001---- With this algorithm and its extension to higher degrees:---- I suspect (albeit VERY, VERY, VERY weakly) merely removing exact-- duplicates from the returned pairs of t-parameters will suffice to-- eliminate repeated detections, because (aside from intersections-- with multiplicities) these SHOULD occur only at endpoints, which-- adjacent portions share, and where the t-parameters are exact zeros-- and ones.---- In any case, comparing t-parameters is certainly an alternative to-- comparing point distances, especially if you want to let-- intersections have multiplicity (as can easily happen with-- cubics). Scheme's SRFI-1 has "delete-duplicates", which lets one-- specify an equivalence predicate other than the default "equal?"---- the predicate can be defined as a closure to test closeness to-- within some tolerance. That is how the code below SHOULD be-- written.---- But I do not know how to do the same thing so simply in Phix, and-- thus will merely say "unique" here and let others update the code-- if they wish. :)--sequencet_pairs=unique(find_intersections(p,q,flatness_tolerance))printf(1,"\n")printf(1," convex up ")printf(1," convex left\n")forttint_pairsdoatomtp,tq{tp,tq}=ttsequenceppoint=curve_eval(p,tp),qpoint=curve_eval(q,tq)printf(1," %11.8f (%11.8f, %11.8f) %11.8f (%11.8f, %11.8f)\n",{tp,ppoint[X],ppoint[Y],tq,qpoint[X],qpoint[Y]})endforprintf(1,"\n")
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
This is mostly based on the Icon, but with aspects of the Pascal. In particular, using a set instead of recursion is like the Icon, but having subprograms for computing and comparing lengths is borrowed from the Pascal. (Perhaps that enhancement will be retrofitted into the Icon at some time.)
See also:
#!/bin/env python3## * * *## This is the algorithm that was introduced with the Icon example, and# perhaps is new (at least in its details). It works by putting both# curves into the symmetric power basis, then first breaking them at# their critical points, then doing an adaptive flattening process# until the problem is reduced to the intersection of two# lines. Possible lines of inquiry are pruned by looking for overlap# of the rectangles formed by the endpoints of curve portions.## Unlike Icon, Python does not have goal-directed evaluation# (GDE). What Python DOES have are "iterables" and# "comprehensions". Where you see "yield" and comprehensions in the# Python you will likely see "suspend" and "every" in the Icon.## To put it another way: In Python, there are objects that "can be# iterated over". In Icon, there are objects that "can produce values# more than once". In either case, the objects are equivalent to a set# (albeit an ordered set), and really what THIS algorithm deals with# is (unordered) sets.## Another thing about Icon to be aware of, when comparing this# algorithm's implementations, is that Icon does not have boolean# expressions. It has "succeed" and "fail". An Icon expression either# "succeeds" and has a value or it "fails" and has no value. An "if"# construct tests whether an expression succeeded, not what the# expression's value is. (Booleans are easily "faked", of course, if# you want them. The language variant Object Icon explicitly# introduces &yes and &no as boolean values.)## * * *## References on the symmetric power basis:## J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power# basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997,# page 319.## J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis# in geometry processing’, ACM Transactions on Graphics, vol 19# no 1, January 2000, page 35.#deflength(ax,ay):'''Length according to some norm, where (ax,ay) is a "measuring stick" vector. Here I use the max norm.'''assertisinstance(ax,float)assertisinstance(ay,float)returnmax(abs(ax),abs(ay))defcompare_lengths(ax,ay,bx,by):'''Having a compare_lengths function lets one compare lengths in the euclidean metric by comparing the squares of the lengths, and thus avoiding the square root. The following, however, is a general implementation.'''assertisinstance(ax,float)assertisinstance(ay,float)assertisinstance(bx,float)assertisinstance(by,float)len_a=length(ax,ay)len_b=length(bx,by)iflen_a<len_b:cmpval=-1eliflen_a>len_b:cmpval=1else:cmpval=0returncmpvaldefrectangles_overlap(a0,a1,b0,b1):'''Do the rectangles with corners at (a0,a1) and (b0,b1) overlap at all?'''assertisinstance(a0,Point)assertisinstance(a1,Point)assertisinstance(b0,Point)assertisinstance(b1,Point)return((min(a0.x,a1.x)<=max(b0.x,b1.x))and(min(b0.x,b1.x)<=max(a0.x,a1.x))and(min(a0.y,a1.y)<=max(b0.y,b1.y))and(min(b0.y,b1.y)<=max(a0.y,a1.y)))defsegment_parameters(a0,a1,b0,b1):'''Do the line segments (a0,a1) and (b0,b1) intersect? If so, return a tuple of their t-parameter values for the point of intersection, treating them as parametric splines of degree 1. Otherwise return None.'''assertisinstance(a0,Point)assertisinstance(a1,Point)assertisinstance(b0,Point)assertisinstance(b1,Point)retval=Noneaxdiff=a1.x-a0.xaydiff=a1.y-a0.ybxdiff=b1.x-b0.xbydiff=b1.y-b0.ydenom=(axdiff*bydiff)-(aydiff*bxdiff)anumer=((bxdiff*a0.y)-(bydiff*a0.x)+(b0.x*b1.y)-(b1.x*b0.y))ta=anumer/denomif0.0<=taandta<=1.0:bnumer=-((axdiff*b0.y)-(aydiff*b0.x)+(a0.x*a1.y)-(a1.x*a0.y))tb=bnumer/denomif0.0<=tbandtb<=1.0:retval=(ta,tb)returnretvalclassPoint:def__init__(self,x,y):assertisinstance(x,float)assertisinstance(y,float)self.x=xself.y=yclassSPower:'''Non-parametric spline in s-power basis.'''def__init__(self,c0,c1,c2):assertisinstance(c0,float)assertisinstance(c1,float)assertisinstance(c2,float)self.c0=c0self.c1=c1self.c2=c2defval(self,t):'''Evaluate at t.'''assertisinstance(t,float)return(self.c0+(self.c1*t))*(1.0-t)+(self.c2*t)defcenter_coef(self,t0,t1):'''Return the center coefficient for the [t0,t1] portion. (The other coefficients can be found with the val method.)'''assertisinstance(t0,float)assertisinstance(t1,float)returnself.c1*((t1-t0-t0)*t1+(t0*t0))defcritical_points(self):'''Return a set of independent variable values for the critical points that lie in (0,1).'''critpoints=set()ifself.c1!=0:# If c1 is zero then the spline is linear.ifself.c0==self.c2:critpoints={0.5}# The spline is "pulse-like".else:# The root of the derivative is the critical point.t=(0.5*(self.c2+self.c1-self.c0))/self.c1if0.0<tandt<1.0:critpoints={t}returncritpointsclassCurve:'''Parametric spline in s-power basis.'''def__init__(self,x,y):assertisinstance(x,SPower)assertisinstance(y,SPower)self.x=xself.y=y@staticmethoddeffrom_controls(ctl0,ctl1,ctl2):assertisinstance(ctl0,Point)assertisinstance(ctl1,Point)assertisinstance(ctl2,Point)c0x=ctl0.xc0y=ctl0.yc1x=(2.0*ctl1.x)-ctl0.x-ctl2.xc1y=(2.0*ctl1.y)-ctl0.y-ctl2.yc2x=ctl2.xc2y=ctl2.yreturnCurve(SPower(c0x,c1x,c2x),SPower(c0y,c1y,c2y))defval(self,t):'''Evaluate at t.'''assertisinstance(t,float)returnPoint(self.x.val(t),self.y.val(t))defcritical_points(self):'''Return a set of t-parameter values for the critical points that lie in (0,1).'''return(self.x.critical_points()|self.y.critical_points())classPortion:'''Portion of a parametric spline in [t0,t1].'''default_num_pieces=2def__init__(self,curve,t0,t1,endpt0,endpt1):assertisinstance(curve,Curve)assertisinstance(t0,float)assertisinstance(t1,float)assertisinstance(endpt0,Point)assertisinstance(endpt1,Point)self.curve=curveself.t0=t0self.t1=t1self.endpt0=endpt0self.endpt1=endpt1defflat_enough(self,tol):'''Is the Portion close enough to linear to be treated as a line segment?'''# The degree-2 s-power polynomials are 1-t, t(1-t), t. We want# to remove the terms in t(1-t). The maximum of t(1-t) is 1/4,# reached at t=1/2. That accounts for the 1/4=0.25 in the# following.xcentercoef=self.curve.x.center_coef(self.t0,self.t1)ycentercoef=self.curve.y.center_coef(self.t0,self.t1)xlen=self.endpt1.x-self.endpt0.xylen=self.endpt1.y-self.endpt0.yreturncompare_lengths(0.25*xcentercoef,0.25*ycentercoef,tol*xlen,tol*ylen)<=0defsplit(self,num_pieces=default_num_pieces):'''Generate num_pieces sections of the Portion.'''assertisinstance(num_pieces,int)assert2<=num_piecesk=1.0/num_piecests=[(1.0-(k*i))*self.t0+(k*i)*self.t1foriinrange(1,num_pieces)]vals=[self.curve.val(t)fortints]ts=[self.t0]+ts+[self.t1]vals=[self.endpt0]+vals+[self.endpt1]foriinrange(len(ts)-1):yieldPortion(self.curve,ts[i],ts[i+1],vals[i],vals[i+1])deffind_intersections(p,q,tol):'''Generate t-parameter pairs detected as corresponding to intersection points of p and q. There may be duplicate detections. It is assumed those will be weeded out by later processing. The tol parameter specifies the "flatness tolerance" and is a relative tolerance.'''assertisinstance(p,Curve)assertisinstance(q,Curve)# The initial workload is the cartesian product of the monotonic# portions of p and q, respectively.tp=[0.0]+sorted(p.critical_points())+[1.0]tq=[0.0]+sorted(q.critical_points())+[1.0]workload={(Portion(p,tp[i],tp[i+1],p.val(tp[i]),p.val(tp[i+1])),Portion(q,tq[j],tq[j+1],q.val(tq[j]),q.val(tq[j+1])))foriinrange(len(tp)-1)forjinrange(len(tq)-1)}whilelen(workload)!=0:(pportion,qportion)=workload.pop()ifrectangles_overlap(pportion.endpt0,pportion.endpt1,qportion.endpt0,qportion.endpt1):ifpportion.flat_enough(tol):ifqportion.flat_enough(tol):params=segment_parameters(pportion.endpt0,pportion.endpt1,qportion.endpt0,qportion.endpt1)ifparamsisnotNone:(tp,tq)=paramstp=(1-tp)*pportion.t0+tp*pportion.t1tq=(1-tq)*qportion.t0+tq*qportion.t1yield(tp,tq)else:workload|={(pportion,qport)forqportinqportion.split()}else:ifqportion.flat_enough(tol):workload|={(pport,qportion)forpportinpportion.split()}else:workload|={(pport,qport)forpportinpportion.split()forqportinqportion.split()}if__name__=='__main__':flatness_tolerance=0.0001minimum_spacing=0.000001p=Curve.from_controls(Point(-1.0,0.0),Point(0.0,10.0),Point(1.0,0.0))q=Curve.from_controls(Point(2.0,1.0),Point(-8.0,2.0),Point(2.0,3.0))intersections=dict()for(tp,tq)infind_intersections(p,q,flatness_tolerance):pval=p.val(tp)qval=q.val(tq)ifall([(minimum_spacing<=length(pval.x-intersections[t][1].x,pval.y-intersections[t][1].y))and(minimum_spacing<=length(qval.x-intersections[t][3].x,qval.y-intersections[t][3].y))fortinintersections]):intersections[tp]=(tp,pval,tq,qval)print()print(' convex up ',' convex left');forkinsorted(intersections):(tp,pointp,tq,pointq)=intersections[k]print(("%11.8f (%11.8f,%11.8f) "+"%11.8f (%11.8f,%11.8f)")%(tp,pointp.x,pointp.y,tq,pointq.x,pointq.y))print()
convex up convex left 0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) 0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) 0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) 0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507)
# 20231025 Raku programming solutionclassPoint {has ($.x,$.y)isrwisdefault(0) }classQuadSpline {has ($.c0,$.c1,$.c2)isrwisdefault(0) }classQuadCurve {hasQuadSpline ($.x,$.y)isrw }classWorkset {hasQuadCurve ($.p,$.q) }subsubdivideQuadSpline($q,$t) {my$s =1.0 -$t;my ($c0,$c1,$c2) =dogiven$q { .c0, .c1, .c2 };my$u_c1 =$s*$c0 +$t*$c1;my$v_c1 =$s*$c1 +$t*$c2;my$u_c2 =$s*$u_c1 +$t*$v_c1;return ($c0,$u_c1,$u_c2), ($u_c2,$v_c1,$c2)}subsubdivideQuadCurve($q,$t,$uisrw,$visrw) {with (subdivideQuadSpline($q.x,$t),subdivideQuadSpline($q.y,$t))».List.flat {$u=QuadCurve.new(x =>QuadSpline.new(c0 => .[0],c1 => .[1],c2 => .[2]),y =>QuadSpline.new(c0 => .[6],c1 => .[7],c2 => .[8]));$v=QuadCurve.new(x =>QuadSpline.new(c0 => .[3],c1 => .[4],c2 => .[5]),y =>QuadSpline.new(c0 => .[9],c1 => .[10],c2 => .[11])) }}subseemsToBeDuplicate(@intersects,$xy,$spacing) {my$seemsToBeDup =False;for@intersects {$seemsToBeDup =abs(.x -$xy.x) <$spacing &&abs(.y -$xy.y) <$spacing;lastif$seemsToBeDup; }return$seemsToBeDup;}subrectsOverlap($xa0,$ya0,$xa1,$ya1,$xb0,$yb0,$xb1,$yb1) {return$xb0 <=$xa1 &&$xa0 <=$xb1 &&$yb0 <=$ya1 &&$ya0 <=$yb1}subtestIntersect($p,$q,$tol,$excludeisrw,$acceptisrw,$intersectisrw) {my$pxmin =min($p.x.c0,$p.x.c1,$p.x.c2);my$pymin =min($p.y.c0,$p.y.c1,$p.y.c2);my$pxmax =max($p.x.c0,$p.x.c1,$p.x.c2);my$pymax =max($p.y.c0,$p.y.c1,$p.y.c2);my$qxmin =min($q.x.c0,$q.x.c1,$q.x.c2);my$qymin =min($q.y.c0,$q.y.c1,$q.y.c2);my$qxmax =max($q.x.c0,$q.x.c1,$q.x.c2);my$qymax =max($q.y.c0,$q.y.c1,$q.y.c2); ($exclude,$accept) =True,False;ifrectsOverlap($pxmin,$pymin,$pxmax,$pymax,$qxmin,$qymin,$qxmax,$qymax) {$exclude =False;my ($xmin,$xmax) =max($pxmin,$qxmin),min($pxmax,$pxmax);if ($xmax <$xmin) {die"Assertion failure: $xmax < $xmin\n" }my ($ymin,$ymax) =max($pymin,$qymin),min($pymax,$qymax);if ($ymax <$ymin) {die"Assertion failure: $ymax < $ymin\n" }if$xmax -$xmin <=$toland$ymax -$ymin <=$tol {$accept =True;given$intersect { (.x, .y) =0.5X*$xmin+$xmax,$ymin+$ymax } } }}subfind-intersects($p,$q,$tol,$spacing) {myPoint@intersects;my@workload =Workset.new(p =>$p,q => $q); while my $work =@workload.pop {my ($intersect,$exclude,$accept) =Point.new,False,False;testIntersect($work.p,$work.q, $tol,$exclude,$accept,$intersect);if$accept {unlessseemsToBeDuplicate(@intersects,$intersect,$spacing) {@intersects.push:$intersect; } }elsifnot$exclude {myQuadCurve ($p0,$p1,$q0,$q1);subdivideQuadCurve($work.p,0.5,$p0,$p1);subdivideQuadCurve($work.q, 0.5,$q0,$q1);for$p0,$p1X$q0,$q1 {@workload.push:Workset.new(p => .[0],q => .[1]) } } } return @intersects;}my $p =QuadCurve.new(x =>QuadSpline.new(c0 => -1.0,c1 =>0.0,c2 =>1.0),y =>QuadSpline.new(c0 =>0.0,c1 =>10.0,c2 =>0.0));my$q =QuadCurve.new(x =>QuadSpline.new(c0 =>2.0,c1 => -8.0,c2 =>2.0),y =>QuadSpline.new(c0 =>1.0,c1 =>2.0,c2 =>3.0));my$spacing = (my$tol =0.0000001 ) *10;.sayforfind-intersects($p,$q,$tol,$spacing);
You mayAttempt This Online!
usestd::f64;constTOLERANCE:f64=0.000_000_1;constSPACING:f64=10.0*TOLERANCE;//Replacesstd::pair<double,double>#[derive(Debug,Clone,Copy,PartialEq,Default)]structPoint{x:f64,y:f64,}implPoint{fnnew(x:f64,y:f64)->Self{Point{x,y}}}//Replacesquad_splineclass#[derive(Debug,Clone,Copy,Default)]structQuadSpline{c0:f64,c1:f64,c2:f64,}implQuadSpline{fnnew(c0:f64,c1:f64,c2:f64)->Self{QuadSpline{c0,c1,c2}}///deCasteljau'salgorithmfor1Dsplinesubdivision///Returnstwonewsplines,representingthecurvefrom0->tandt->1fnsubdivide(&self,t:f64)->(QuadSpline,QuadSpline){lets=1.0-t;letu_c0=self.c0;letv_c2=self.c2;letu_c1=s*self.c0+t*self.c1;letv_c1=s*self.c1+t*self.c2;letu_c2=s*u_c1+t*v_c1;letv_c0=u_c2;letu=QuadSpline::new(u_c0,u_c1,u_c2);letv=QuadSpline::new(v_c0,v_c1,v_c2);(u,v)}//Helpertogetmin/maxcontrolpointsforboundingboxfnmin_max(&self)->(f64,f64){letmin_val=self.c0.min(self.c1).min(self.c2);letmax_val=self.c0.max(self.c1).max(self.c2);(min_val,max_val)}}//Replacesquad_curveclass#[derive(Debug,Clone,Copy,Default)]structQuadCurve{x:QuadSpline,y:QuadSpline,}implQuadCurve{fnnew(x:QuadSpline,y:QuadSpline)->Self{QuadCurve{x,y}}///deCasteljau'salgorithmfor2Dcurvesubdivision///Returnstwonewcurves,representingthecurvefrom0->tandt->1fnsubdivide(&self,t:f64)->(QuadCurve,QuadCurve){let(ux,vx)=self.x.subdivide(t);let(uy,vy)=self.y.subdivide(t);(QuadCurve::new(ux,uy),QuadCurve::new(vx,vy))}///Calculatetheaxis-alignedboundingbox(AABB)fnbounding_box(&self)->(f64,f64,f64,f64){let(px_min,px_max)=self.x.min_max();let(py_min,py_max)=self.y.min_max();(px_min,py_min,px_max,py_max)}}///Checksiftwoaxis-alignedrectanglesoverlapfnrectangles_overlap(xa0:f64,ya0:f64,xa1:f64,ya1:f64,xb0:f64,yb0:f64,xb1:f64,yb1:f64,)->bool{//Ensurecoordinatesareorderedmin->maxifnecessary,althoughbounding_boxshouldhandlethis.//let(xa0,xa1)=ifxa0>xa1{(xa1,xa0)}else{(xa0,xa1)};//let(ya0,ya1)=ifya0>ya1{(ya1,ya0)}else{(ya0,ya1)};//let(xb0,xb1)=ifxb0>xb1{(xb1,xb0)}else{(xb0,xb1)};//let(yb0,yb1)=ifyb0>yb1{(yb1,yb0)}else{(yb0,yb1)};xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1}///Testintersectioncandidacybasedonboundingboxoverlapandsize.///Returns:(accepted,excluded,potential_intersect_point)fntest_intersection(p:&QuadCurve,q:&QuadCurve)->(bool,bool,Point){let(px_min,py_min,px_max,py_max)=p.bounding_box();let(qx_min,qy_min,qx_max,qy_max)=q.bounding_box();letmutaccepted=false;letmutexcluded=true;letmutintersect=Point::default(); // Default point (0.0, 0.0)ifrectangles_overlap(px_min,py_min,px_max,py_max,qx_min,qy_min,qx_max,qy_max){excluded=false;//Calculateoverlapregionletx_min=px_min.max(qx_min);letx_max=px_max.min(qx_max); // Corrected from C++ possible typo px_max.min(px_max)ifx_max-x_min<=TOLERANCE{lety_min=py_min.max(qy_min);lety_max=py_max.min(qy_max);ify_max-y_min<=TOLERANCE{accepted=true;//Usemidpointofthetinyoverlapboxastheintersectionpointintersect=Point::new(0.5*(x_min+x_max),0.5*(y_min+y_max));}}}(accepted,excluded,intersect)}///Checkifapointistooclosetoexistingintersectionpointsfnseems_to_be_duplicate(intersects:&[Point],pt:Point)->bool{for&intersectinintersects{if(intersect.x-pt.x).abs()<SPACING&&(intersect.y-pt.y).abs()<SPACING{returntrue;}}false}///FindintersectionpointsbetweentwoquadraticBeziercurvesusingrecursivesubdivision.fnfind_intersects(p:&QuadCurve,q:&QuadCurve)->Vec<Point>{letmutresult:Vec<Point>=Vec::new();//UseaVecasastack,storingpairsofcurvestocheck//Pushing/poppingindividualcurvestomimicC++stackbehaviorexactlyletmutstack:Vec<QuadCurve>=Vec::new();stack.push(*q); // Push q first (will be popped as qq)stack.push(*p); // Push p second (will be popped as pp)whilestack.len()>=2{//PopintheorderthatmatchestheC++version'sstackprocessingletpp=stack.pop().unwrap(); // Pop p-curve segmentletqq=stack.pop().unwrap(); // Pop q-curve segmentlet(accepted,excluded,intersect)=test_intersection(&pp,&qq);ifaccepted{if!seems_to_be_duplicate(&result,intersect){result.push(intersect);}}elseif!excluded{//Boundingboxesoverlapsignificantly,subdividebothcurveslet(p0,p1)=pp.subdivide(0.5);let(q0,q1)=qq.subdivide(0.5);//Pushthe4pairsofsub-curvesontothestackforfurthertesting.//OrdermatchestheC++version'spushorder{p0,q0,p0,q1,p1,q0,p1,q1}//Processedpairswillbe(p0,q0),(p0,q1),(p1,q0),(p1,q1)eventually.//Pushinreverseorderofdesiredprocessing:stack.push(q1); stack.push(p1); // For pair (p1, q1)stack.push(q0); stack.push(p1); // For pair (p1, q0)stack.push(q1); stack.push(p0); // For pair (p0, q1)stack.push(q0); stack.push(p0); // For pair (p0, q0)}//Ifexcluded,donothing,thispairofsegmentscannotintersect.}result}fnmain(){//QuadCurveverticalrepresentstheBeziercurvehavingcontrolpoints(-1,0),(0,10)and(1,0)letvertical=QuadCurve::new(QuadSpline::new(-1.0,0.0,1.0),QuadSpline::new(0.0,10.0,0.0),);//QuadCurvehorizontalrepresentstheBeziercurvehavingcontrolpoints(2,1),(-8,2)and(2,3)lethorizontal=QuadCurve::new(QuadSpline::new(2.0,-8.0,2.0),QuadSpline::new(1.0,2.0,3.0),);println!("The points of intersection are:");letintersects=find_intersects(&vertical,&horizontal);ifintersects.is_empty(){println!("No intersections found.");}else{forintersectinintersects{//FormatoutputsimilartoC++versionprintln!("( {:9.6}, {:9.6} )",intersect.x,intersect.y);}}}
The points of intersection are:( -0.854983, 1.345017 )( -0.681025, 2.681025 )( 0.881025, 1.118975 )( 0.654983, 2.854983 )
For R7RS-small, plus SRFI-132 and SRFI-144, both adopted as part of R7RS-large. SRFI-144 is used only to get the value of fl-epsilon.
This is the algorithm of theIcon. The implementation here has support for three different norms and for both relative and absolute tolerance.
;;;; For R7RS-small plus very limited R7RS-large. The implementation;; must support libraries, Unicode identifiers, IEEE arithmetic, etc.;;;; These will work:;;;; gosh -r7 file.scm;; csc -X r7rs -R r7rs file.scm;;(define-library(spowerquadratic)(exportspowerspower?spower-c0spower-c1spower-c2plane-curveplane-curve?plane-curve-xplane-curve-ycontrol-points->plane-curvespower-evalplane-curve-evalcenter-coefcritical-points)(import(schemebase)(schemecase-lambda)(srfi132);; R7RS-large name: (scheme sort))(begin(define-record-type<spower>(spowerc0c1c2)spower?(c0spower-c0)(c1spower-c1)(c2spower-c2))(define-record-type<plane-curve>(plane-curvexy)plane-curve?(xplane-curve-x)(yplane-curve-y))(define(point->valuespoint);; A bit of playfulness on my part: a few ways to write an;; (x,y)-point: (cons x y), (list x y), (vector x y), and their;; equivalents.(cond((and(pair?point)(real?(carpoint)))(cond((real?(cdrpoint))(values(carpoint)(cdrpoint)))((and(real?(cadrpoint))(null?(cddrpoint)))(values(carpoint)(cadrpoint)))(else(error"not a plane point"point))))((and(vector?point)(=(vector-lengthpoint)2))(values(vector-refpoint0)(vector-refpoint1)))(else(error"not a plane point"point))))(definecontrol-points->plane-curve;; A bit more playfulness: the control points can be given;; separately, as a list, or as a vector.(case-lambda((pt0pt1pt2)(let-values(((cx0cy0)(point->valuespt0))((cx1cy1)(point->valuespt1))((cx2cy2)(point->valuespt2)))(plane-curve(spowercx0(-(+cx1cx1)cx0cx2)cx2)(spowercy0(-(+cy1cy1)cy0cy2)cy2))))((pt-list-or-vector)(applycontrol-points->plane-curve(if(vector?pt-list-or-vector)(vector->listpt-list-or-vector)pt-list-or-vector)))))(define(spower-evalpolyt);; Evaluate the polynomial at t.(let((c0(spower-c0poly))(c1(spower-c1poly))(c2(spower-c2poly)))(+(*(+c0(*c1t))(-1t))(*c2t))))(define(plane-curve-evalcurvet);; Evaluate the plane curve at t, returning x and y as multiple;; values.(values(spower-eval(plane-curve-xcurve)t)(spower-eval(plane-curve-ycurve)t)))(define(center-coefpolyt0t1);; Return the center coefficient for the [t0,t1] portion. (The;; other coefficients can be found with the spower-eval;; procedured.)(let((c1(spower-c1poly)))(*c1(+(*(-t1t0t0)t1)(*t0t0)))))(define(critical-pointspoly-or-curve);; Return a sorted list of independent variable values for the;; critical points that lie in (0,1). Duplicates are removed.(cond((plane-curve?poly-or-curve)(let((X(plane-curve-xpoly-or-curve))(Y(plane-curve-ypoly-or-curve)))(list-delete-neighbor-dups=(list-sort<(append(critical-pointsX)(critical-pointsY))))))((spower?poly-or-curve)(let((c0(spower-c0poly-or-curve))(c1(spower-c1poly-or-curve))(c2(spower-c2poly-or-curve)))(cond((zero?c1)'()); The spline is linear.((=c0c2)'(1/2)); The spline is "pulse-like".(else;; The root of the derivative is the critical;; point.(let((t(/(-(+c2c1)c0)(+c1c1))))(if(and(positive?t)(<t1))(listt)'()))))))(else(error"not an spower polynomial or plane-curve"poly-or-curve))))));; end library (spower quadratic)(define-library(spowerquadraticintersections);; Parameters. (That is, variables whose value settings have;; "dynamic" rather than lexical scope.)(export*tolerance-norm**flatness-tolerance**absolute-tolerance*)(exportfind-intersection-parameters)(import(schemebase)(schemecase-lambda)(schemeinexact)(spowerquadratic)(srfi144);; R7RS-large name: (scheme flonum))(begin(define-record-type<portion>(make-portioncurvet0t1endpt0endpt1)portion?(curvecurve@)(t0t0@)(t1t1@)(endpt0endpt0@)(endpt1endpt1@))(define(curve-evalcurvet)(call-with-values(lambda()(plane-curve-evalcurvet))cons))(define(0<=x<=1x)(and(<=0x)(<=x1)))(define(lerptab); "linear interpolation"(+(*(-1t)a)(*tb)))(define(bisect-portionportion);; Bisect portion and return the two new portions as multiple;; values.(let((curve(curve@portion))(t0(t0@portion))(t1(t1@portion))(pt0(endpt0@portion))(pt1(endpt1@portion)))(let*((t½(*1/2(+t0t1)))(pt½(curve-evalcurvet½)))(values(make-portioncurvet0t½pt0pt½)(make-portioncurvet½t1pt½pt1)))))(define(check-normx)(cond((and(positive?x)(infinite?x))x)((=x1)(exactx))((=x2)(exactx))(else(error"not a norm we can handle"x))))(define(check-flatness-tolerancex)(cond((zero?x)x)((positive?x)x)(else(error"not a flatness (relative) tolerance we can handle"x))))(define(check-absolute-tolerancex)(cond((zero?x)x)((positive?x)x)(else(error"not an absolute tolerance we can handle"x))))(define*tolerance-norm*;; To be fancier, we could take strings such as "taxicab",;; "euclidean", "max", etc., as arguments to the;; constructor. But here we are taking only the value of p in a;; p-norm (= pth root of the sum of |x| raised p and |y| raised;; p), and then only for p = 1, 2, +inf.0.(make-parameter+inf.0check-norm));; Default is the max norm.;; A relative tolerance. This setting seems to me rather strict;; for a lot of applications. But you can override it.(define*flatness-tolerance*(make-parameter0.0001check-flatness-tolerance))(define*absolute-tolerance*(make-parameter(*50fl-epsilon)check-absolute-tolerance))(define(compare-lengthsnormaxaybxby)(define(compare-lengths-1-norm)(let((‖a‖(+(absax)(absay)))(‖b‖(+(absbx)(absby))))(cond((<‖a‖‖b‖)-1)((>‖a‖‖b‖)1)(else0))))(define(compare-lengths-2-norm)(let((‖a‖²(*axay))(‖b‖²(*bxby)))(cond((<‖a‖²‖b‖²)-1)((>‖a‖²‖b‖²)1)(else0))))(define(compare-lengths-inf-norm)(let((‖a‖(max(absax)(absay)))(‖b‖(max(absbx)(absby))))(cond((<‖a‖‖b‖)-1)((>‖a‖‖b‖)1)(else0))))(cond((=norm1)(compare-lengths-1-norm))((=norm2)(compare-lengths-2-norm))(else(compare-lengths-inf-norm))))(define(compare-to-tolnormaxaytol)(define(compare-to-tol-1-norm)(let((‖a‖(+(absax)(absay))))(cond((<‖a‖tol)-1)((>‖a‖tol)1)(else0))))(define(compare-to-tol-2-norm)(let((‖a‖²(*axay))(tol²(*toltol)))(cond((<‖a‖²tol²)-1)((>‖a‖²tol²)1)(else0))))(define(compare-to-tol-inf-norm)(let((‖a‖(max(absax)(absay))))(cond((<‖a‖tol)-1)((>‖a‖tol)1)(else0))))(cond((=norm1)(compare-to-tol-1-norm))((=norm2)(compare-to-tol-2-norm))(else(compare-to-tol-inf-norm))))(define(flat-enough?portionnormrtolatol);; Is the portion flat enough or small enough to be treated as;; if it were a line segment?;; The degree-2 s-power polynomials are 1-t, t(1-t), t. We;; want to remove the terms in t(1-t). The maximum of t(1-t);; is 1/4, reached at t=1/2. That accounts for the 1/4 in the;; following.(let((curve(curve@portion))(t0(t0@portion))(t1(t1@portion))(pt0(endpt0@portion))(pt1(endpt1@portion)))(let((X(plane-curve-xcurve))(Y(plane-curve-ycurve)))(let((errx(*1/4(center-coefXt0t1)))(erry(*1/4(center-coefYt0t1)))(testx(*rtol(-(carpt1)(carpt0))))(testy(*rtol(-(cdrpt1)(cdrpt0)))))(or(<=(compare-lengthsnormerrxerrytestxtesty)0)(<=(compare-to-tolnormerrxerryatol)0))))))(define(rectangles-overlapa0a1b0b1);;;; Do the rectangles with corners at (a0,a1) and (b0,b1) have;; overlapping interiors or edges? The corner points must be;; represented as cons-pairs, with x as the car and y as the;; cdr.;;;; (A side note: I had thought for a few hours that checking;; only for overlapping interiors offered some advantages, but I;; changed my mind.);;(let((a0x(min(cara0)(cara1)))(a1x(max(cara0)(cara1)))(b0x(min(carb0)(carb1)))(b1x(max(carb0)(carb1))))(cond((<b1xa0x)#f)((<a1xb0x)#f)(else(let((a0y(min(cdra0)(cdra1)))(a1y(max(cdra0)(cdra1)))(b0y(min(cdrb0)(cdrb1)))(b1y(max(cdrb0)(cdrb1))))(cond((<b1ya0y)#f)((<a1yb0y)#f)(else#t)))))))(define(segment-parametersa0a1b0b1);; Return (as multiple values) the respective [0,1] parameters;; of line segments (a0,a1) and (b0,b1), for their intersection;; point. Return #f and #f if there is no intersection. The;; endpoints must be represented as cons-pairs, with x as the;; car and y as the cdr.(let((a0x(cara0))(a0y(cdra0))(a1x(cara1))(a1y(cdra1))(b0x(carb0))(b0y(cdrb0))(b1x(carb1))(b1y(cdrb1)))(let((axdiff(-a1xa0x))(aydiff(-a1ya0y))(bxdiff(-b1xb0x))(bydiff(-b1yb0y)))(let*((denom(-(*axdiffbydiff)(*aydiffbxdiff)))(anumer(+(*bxdiffa0y)(-(*bydiffa0x))(*b0xb1y)(-(*b1xb0y))))(ta(/anumerdenom)))(if(not(0<=x<=1ta))(values#f#f)(let*((bnumer(-(+(*axdiffb0y)(-(*aydiffb0x))(*a0xa1y)(-(*a1xa0y)))))(tb(/bnumerdenom)))(if(not(0<=x<=1tb))(values#f#f)(valuestatb))))))))(define(initial-workloadPQ);; Generate an initial workload from plane curves P and Q. One;; does this by splitting the curves at their critical points;; and constructing the Cartesian product of the resulting;; portions. The workload is a list of cons-pairs of;; portions. (The list represents a set, but the obvious place,;; from which to get an arbitrary pair to work on next, is the;; head of the list. Thus the list effectively is a stack.)(define(t->pointcurve)(lambda(t)(curve-evalcurvet)))(let*((pparams0`(0,@(critical-pointsP)1))(pvalues0(map(t->pointP)pparams0))(qparams0`(0,@(critical-pointsQ)1))(qvalues0(map(t->pointQ)qparams0)))(letloop((pparamspparams0)(pvaluespvalues0)(qparamsqparams0)(qvaluesqvalues0)(workload'()))(cond((null?(cdrpparams))workload)((null?(cdrqparams))(loop(cdrpparams)(cdrpvalues)qparams0qvalues0workload))(else(let((pportion(make-portionP(carpparams)(cadrpparams)(carpvalues)(cadrpvalues)))(qportion(make-portionQ(carqparams)(cadrqparams)(carqvalues)(cadrqvalues))))(looppparamspvalues(cdrqparams)(cdrqvalues)(cons(conspportionqportion)workload))))))))(define(params≅?ab)(and(≅?(cara)(carb))(≅?(cdra)(cdrb))))(define(≅?xy);; I do not know this test is any better, for THIS algorithm,;; than an exact equality test. But it is no worse.(<=(abs(-xy))(*0.5fl-epsilon(max(absx)(absy)))))(define(include-paramstptqlst)(let((params(constptq)))(if(not(memberparamslstparams≅?))(consparamslst)lst)))(definefind-intersection-parameters(case-lambda((PQ)(find-intersection-parametersPQ#f#f#f))((PQnorm)(find-intersection-parametersPQnorm#f#f))((PQnormrtol)(find-intersection-parametersPQnormrtol#f))((PQnormrtolatol)(let((norm(check-norm(ornorm(*tolerance-norm*))))(rtol(check-flatness-tolerance(orrtol(*flatness-tolerance*))))(atol(check-absolute-tolerance(oratol(*absolute-tolerance*)))))(%%find-intersection-parametersPQnormrtolatol)))))(define(%%find-intersection-parametersPQnormrtolatol)(letloop((workload(initial-workloadPQ))(params'()))(if(null?workload)params(let((pportion(caarworkload))(qportion(cdarworkload))(workload(cdrworkload)))(if(not(rectangles-overlap(endpt0@pportion)(endpt1@pportion)(endpt0@qportion)(endpt1@qportion)))(loopworkloadparams)(if(flat-enough?pportionnormrtolatol)(if(flat-enough?qportionnormrtolatol)(let-values(((tptq)(segment-parameters(endpt0@pportion)(endpt1@pportion)(endpt0@qportion)(endpt1@qportion))))(iftp(let*((tp0(t0@pportion))(tp1(t1@pportion))(tq0(t0@qportion))(tq1(t1@qportion))(tp(lerptptp0tp1))(tq(lerptqtq0tq1)))(loopworkload(include-paramstptqparams)))(loopworkloadparams)))(let-values(((qport1qport2)(bisect-portionqportion)))(loop`(,(conspportionqport1),(conspportionqport2).,workload)params)))(if(flat-enough?qportionnormrtolatol)(let-values(((pport1pport2)(bisect-portionpportion)))(loop`(,(conspport1qportion),(conspport2qportion).,workload)params))(let-values(((pport1pport2)(bisect-portionpportion))((qport1qport2)(bisect-portionqportion)))(loop`(,(conspport1qport1),(conspport1qport2),(conspport2qport1),(conspport2qport2).,workload)params)))))))))));; end library (spower quadratic intersections)(import(schemebase)(schemewrite)(spowerquadratic)(spowerquadraticintersections)(srfi132);; R7RS-large name: (scheme sort))(defineP(control-points->plane-curve'((-10)(010)(10))))(defineQ(control-points->plane-curve'((21)(-82)(23))));; Sort the results by the parameters for P. Thus the intersection;; points will be sorted left to right.(defineparams(list-sort(lambda(tpair1tpair2)(<(cartpair1)(cartpair2)))(find-intersection-parametersPQ)))(for-each(lambda(pair)(let((tp(carpair))(tq(cdrpair)))(let-values(((axay)(plane-curve-evalPtp))((bxby)(plane-curve-evalQtq)))(display(inexacttp))(display"\t(")(display(inexactax))(display", ")(display(inexactay))(display")\t")(display(inexacttq))(display"\t(")(display(inexactbx))(display", ")(display(inexactby))(display")\n"))))params)
There is no standard for formatted output in R7RS Scheme, so in the following I do not bother. (This does not mean there are not excellentunstandardized ways to do formatted output in Schemes.) This is the output from Gauche Scheme.
0.07250828117222911 (-0.8549834376555417, 1.3450166066735616) 0.17250829973466453 (-0.8549837251463935, 1.345016599469329)0.15948753092173137 (-0.6810249381565372, 2.6810251680444233) 0.8405124690782686 (-0.6810251680444231, 2.681024938156537)0.8274917002653355 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542)0.9405124666929737 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525)
It is possible to reduce the number of bisections by solving a quadratic equation as soon as just one of the curve portions is flat enough to be treated as if it were a line segment. That refinement is implemented here.
;;;; For R7RS-small plus very limited R7RS-large. The implementation;; must support libraries, Unicode identifiers, IEEE arithmetic, etc.;;;; These will work:;;;; gosh -r7 file.scm;; csc -X r7rs -R r7rs file.scm;;(define-library(spowerquadratic)(exportspowerspower?spower-c0spower-c1spower-c2plane-curveplane-curve?plane-curve-xplane-curve-ycontrol-points->plane-curvespower-evalplane-curve-evalcenter-coefcritical-points)(import(schemebase)(schemecase-lambda)(srfi132);; R7RS-large name: (scheme sort))(begin(define-record-type<spower>(spowerc0c1c2)spower?(c0spower-c0)(c1spower-c1)(c2spower-c2))(define-record-type<plane-curve>(plane-curvexy)plane-curve?(xplane-curve-x)(yplane-curve-y))(define(point->valuespoint);; A bit of playfulness on my part: a few ways to write an;; (x,y)-point: (cons x y), (list x y), (vector x y), and their;; equivalents.(cond((and(pair?point)(real?(carpoint)))(cond((real?(cdrpoint))(values(carpoint)(cdrpoint)))((and(real?(cadrpoint))(null?(cddrpoint)))(values(carpoint)(cadrpoint)))(else(error"not a plane point"point))))((and(vector?point)(=(vector-lengthpoint)2))(values(vector-refpoint0)(vector-refpoint1)))(else(error"not a plane point"point))))(definecontrol-points->plane-curve;; A bit more playfulness: the control points can be given;; separately, as a list, or as a vector.(case-lambda((pt0pt1pt2)(let-values(((cx0cy0)(point->valuespt0))((cx1cy1)(point->valuespt1))((cx2cy2)(point->valuespt2)))(plane-curve(spowercx0(-(+cx1cx1)cx0cx2)cx2)(spowercy0(-(+cy1cy1)cy0cy2)cy2))))((pt-list-or-vector)(applycontrol-points->plane-curve(if(vector?pt-list-or-vector)(vector->listpt-list-or-vector)pt-list-or-vector)))))(define(spower-evalpolyt);; Evaluate the polynomial at t.(let((c0(spower-c0poly))(c1(spower-c1poly))(c2(spower-c2poly)))(+(*(+c0(*c1t))(-1t))(*c2t))))(define(plane-curve-evalcurvet);; Evaluate the plane curve at t, returning x and y as multiple;; values.(values(spower-eval(plane-curve-xcurve)t)(spower-eval(plane-curve-ycurve)t)))(define(center-coefpolyt0t1);; Return the center coefficient for the [t0,t1] portion. (The;; other coefficients can be found with the spower-eval;; procedured.)(let((c1(spower-c1poly)))(*c1(+(*(-t1t0t0)t1)(*t0t0)))))(define(critical-pointspoly-or-curve);; Return a sorted list of independent variable values for the;; critical points that lie in (0,1). Duplicates are removed.(cond((plane-curve?poly-or-curve)(let((X(plane-curve-xpoly-or-curve))(Y(plane-curve-ypoly-or-curve)))(list-delete-neighbor-dups=(list-sort<(append(critical-pointsX)(critical-pointsY))))))((spower?poly-or-curve)(let((c0(spower-c0poly-or-curve))(c1(spower-c1poly-or-curve))(c2(spower-c2poly-or-curve)))(cond((zero?c1)'()); The spline is linear.((=c0c2)'(1/2)); The spline is "pulse-like".(else;; The root of the derivative is the critical;; point.(let((t(/(-(+c2c1)c0)(+c1c1))))(if(and(positive?t)(<t1))(listt)'()))))))(else(error"not an spower polynomial or plane-curve"poly-or-curve))))));; end library (spower quadratic)(define-library(spowerquadraticintersections);; Parameters. (That is, variables whose value settings have;; "dynamic" rather than lexical scope.)(export*tolerance-norm**flatness-tolerance**absolute-tolerance*)(exportfind-intersection-parameters)(import(schemebase)(schemecase-lambda)(schemeinexact)(spowerquadratic)(srfi144);; R7RS-large name: (scheme flonum));; REMOVE THIS FOR A PRACTICAL APPLICATION.(import(schemewrite))(begin(define-record-type<portion>(make-portioncurvet0t1endpt0endpt1)portion?(curvecurve@)(t0t0@)(t1t1@)(endpt0endpt0@)(endpt1endpt1@))(define(curve-evalcurvet)(call-with-values(lambda()(plane-curve-evalcurvet))cons))(define(0<=x<=1x)(and(<=0x)(<=x1)))(define(lerptab); "linear interpolation"(+(*(-1t)a)(*tb)))(define(bisect-portionportion);; Bisect portion and return the two new portions as multiple;; values.(let((curve(curve@portion))(t0(t0@portion))(t1(t1@portion))(pt0(endpt0@portion))(pt1(endpt1@portion)))(let*((t½(*1/2(+t0t1)))(pt½(curve-evalcurvet½)))(values(make-portioncurvet0t½pt0pt½)(make-portioncurvet½t1pt½pt1)))))(define(check-normx)(cond((and(positive?x)(infinite?x))x)((=x1)(exactx))((=x2)(exactx))(else(error"not a norm we can handle"x))))(define(check-flatness-tolerancex)(cond((zero?x)x)((positive?x)x)(else(error"not a flatness (relative) tolerance we can handle"x))))(define(check-absolute-tolerancex)(cond((zero?x)x)((positive?x)x)(else(error"not an absolute tolerance we can handle"x))))(define*tolerance-norm*;; To be fancier, we could take strings such as "taxicab",;; "euclidean", "max", etc., as arguments to the;; constructor. But here we are taking only the value of p in a;; p-norm (= pth root of the sum of |x| raised p and |y| raised;; p), and then only for p = 1, 2, +inf.0.(make-parameter+inf.0check-norm));; Default is the max norm.;; A relative tolerance. This setting seems to me rather strict;; for a lot of applications. But you can override it.(define*flatness-tolerance*(make-parameter0.0001check-flatness-tolerance))(define*absolute-tolerance*(make-parameter(*50fl-epsilon)check-absolute-tolerance))(define(compare-lengthsnormaxaybxby)(define(compare-lengths-1-norm)(let((‖a‖(+(absax)(absay)))(‖b‖(+(absbx)(absby))))(cond((<‖a‖‖b‖)-1)((>‖a‖‖b‖)1)(else0))))(define(compare-lengths-2-norm)(let((‖a‖²(*axay))(‖b‖²(*bxby)))(cond((<‖a‖²‖b‖²)-1)((>‖a‖²‖b‖²)1)(else0))))(define(compare-lengths-inf-norm)(let((‖a‖(max(absax)(absay)))(‖b‖(max(absbx)(absby))))(cond((<‖a‖‖b‖)-1)((>‖a‖‖b‖)1)(else0))))(cond((=norm1)(compare-lengths-1-norm))((=norm2)(compare-lengths-2-norm))(else(compare-lengths-inf-norm))))(define(compare-to-tolnormaxaytol)(define(compare-to-tol-1-norm)(let((‖a‖(+(absax)(absay))))(cond((<‖a‖tol)-1)((>‖a‖tol)1)(else0))))(define(compare-to-tol-2-norm)(let((‖a‖²(*axay))(tol²(*toltol)))(cond((<‖a‖²tol²)-1)((>‖a‖²tol²)1)(else0))))(define(compare-to-tol-inf-norm)(let((‖a‖(max(absax)(absay))))(cond((<‖a‖tol)-1)((>‖a‖tol)1)(else0))))(cond((=norm1)(compare-to-tol-1-norm))((=norm2)(compare-to-tol-2-norm))(else(compare-to-tol-inf-norm))))(define(flat-enough?portionnormrtolatol);; Is the portion flat enough or small enough to be treated as;; if it were a line segment?;; The degree-2 s-power polynomials are 1-t, t(1-t), t. We;; want to remove the terms in t(1-t). The maximum of t(1-t);; is 1/4, reached at t=1/2. That accounts for the 1/4 in the;; following.(let((curve(curve@portion))(t0(t0@portion))(t1(t1@portion))(pt0(endpt0@portion))(pt1(endpt1@portion)))(let((X(plane-curve-xcurve))(Y(plane-curve-ycurve)))(let((errx(*1/4(center-coefXt0t1)))(erry(*1/4(center-coefYt0t1)))(testx(*rtol(-(carpt1)(carpt0))))(testy(*rtol(-(cdrpt1)(cdrpt0)))))(or(<=(compare-lengthsnormerrxerrytestxtesty)0)(<=(compare-to-tolnormerrxerryatol)0))))))(define(rectangles-overlapa0a1b0b1);;;; Do the rectangles with corners at (a0,a1) and (b0,b1) have;; overlapping interiors or edges? The corner points must be;; represented as cons-pairs, with x as the car and y as the;; cdr.;;;; (A side note: I had thought for a few hours that checking;; only for overlapping interiors offered some advantages, but I;; changed my mind.);;(let((a0x(min(cara0)(cara1)))(a1x(max(cara0)(cara1)))(b0x(min(carb0)(carb1)))(b1x(max(carb0)(carb1))))(cond((<b1xa0x)#f)((<a1xb0x)#f)(else(let((a0y(min(cdra0)(cdra1)))(a1y(max(cdra0)(cdra1)))(b0y(min(cdrb0)(cdrb1)))(b1y(max(cdrb0)(cdrb1))))(cond((<b1ya0y)#f)((<a1yb0y)#f)(else#t)))))))(define(segment-parametersa0a1b0b1);; Return (as multiple values) the respective [0,1] parameters;; of line segments (a0,a1) and (b0,b1), for their intersection;; point. Return #f and #f if there is no intersection. The;; endpoints must be represented as cons-pairs, with x as the;; car and y as the cdr.(let((a0x(cara0))(a0y(cdra0))(a1x(cara1))(a1y(cdra1))(b0x(carb0))(b0y(cdrb0))(b1x(carb1))(b1y(cdrb1)))(let((axdiff(-a1xa0x))(aydiff(-a1ya0y))(bxdiff(-b1xb0x))(bydiff(-b1yb0y)))(let*((denom(-(*axdiffbydiff)(*aydiffbxdiff)))(anumer(+(*bxdiffa0y)(-(*bydiffa0x))(*b0xb1y)(-(*b1xb0y))))(ta(/anumerdenom)))(if(not(0<=x<=1ta))(values#f#f)(let*((bnumer(-(+(*axdiffb0y)(-(*aydiffb0x))(*a0xa1y)(-(*a1xa0y)))))(tb(/bnumerdenom)))(if(not(0<=x<=1tb))(values#f#f)(valuestatb))))))))(define(initial-workloadPQ);; Generate an initial workload from plane curves P and Q. One;; does this by splitting the curves at their critical points;; and constructing the Cartesian product of the resulting;; portions. The workload is a list of cons-pairs of;; portions. (The list represents a set, but the obvious place,;; from which to get an arbitrary pair to work on next, is the;; head of the list. Thus the list effectively is a stack.)(define(t->pointcurve)(lambda(t)(curve-evalcurvet)));;;; There are endless ways to write the loop or recursion to;; compute the Cartesian product. The original Scheme did it one;; way. Here is a completely different way. It involves having;; two pointers go through a list at the same time, spaced one;; node apart. This is done on more than one list at a;; time. Also, this time the algorithm is done "procedurally";; rather than "functionally".;;(let*((pparams0`(0,@(critical-pointsP)1))(pvalues0(map(t->pointP)pparams0))(qparams0`(0,@(critical-pointsQ)1))(qvalues0(map(t->pointQ)qparams0))(workload'()))(do((ppipparams0(cdrppi))(ppj(cdrpparams0)(cdrppj))(pvipvalues0(cdrpvi))(pvj(cdrpvalues0)(cdrpvj)))((null?ppj))(do((qpiqparams0(cdrqpi))(qpj(cdrqparams0)(cdrqpj))(qviqvalues0(cdrqvi))(qvj(cdrqvalues0)(cdrqvj)))((null?qpj))(let((pportion(make-portionP(carppi)(carppj)(carpvi)(carpvj)))(qportion(make-portionQ(carqpi)(carqpj)(carqvi)(carqvj))))(set!workload`((,pportion.,qportion).,workload)))))workload))(define(params≅?ab)(and(≅?(cara)(carb))(≅?(cdra)(cdrb))))(define(≅?xy);; PERHAPS IT WOULD BE BETTER TO HAVE THIS DEFINITION BE A;; PARAMETER.(<=(abs(-xy))(*0.5fl-epsilon(max(absx)(absy)))))(define(include-paramstptqlst)(let((params(constptq)))(if(not(memberparamslstparams≅?))(consparamslst)lst)))(definefind-intersection-parameters(case-lambda((PQ)(find-intersection-parametersPQ#f#f#f))((PQnorm)(find-intersection-parametersPQnorm#f#f))((PQnormrtol)(find-intersection-parametersPQnormrtol#f))((PQnormrtolatol)(let((norm(check-norm(ornorm(*tolerance-norm*))))(rtol(check-flatness-tolerance(orrtol(*flatness-tolerance*))))(atol(check-absolute-tolerance(oratol(*absolute-tolerance*)))))(%%find-intersection-parametersPQnormrtolatol)))));; REMOVE THIS FOR A PRACTICAL APPLICATION.(defineNUM-ITERATIONS0)(define(%%find-intersection-parametersPQnormrtolatol);;;; Among other things that this version of the program;; demonstrates: you can safely break up a standard Scheme loop;; into a mutual tail recursion. There will be no "stack;; blow-up". (At least if you do not count as "stack blow-up";; the very strange way CHICKEN Scheme works.);;;; (It is interesting, by the way, to know in what languages one;; can do this sort of thing, to what degree. In standard;; Scheme, it is possible without any restrictions. In ATS, one;; can do it safely as long as an "fnx" construct is possible,;; because this is precisely what "fnx" is for. But I have tried;; a very, very simple mutual recursion in Standard ML, and had;; it work fine in Poly/ML but blow up the stack in MLton, even;; though MLton is overall the more aggressively optimizing;; compiler.);;(define(start-hereworkloadparams);; REMOVE THIS FOR A PRACTICAL APPLICATION.(set!NUM-ITERATIONS(+NUM-ITERATIONS1))(if(null?workload)(begin;; REMOVE THIS FOR A PRACTICAL APPLICATION.(displayNUM-ITERATIONS)(display"\n")params)(let((pportion(caarworkload))(qportion(cdarworkload))(workload(cdrworkload)))(if(not(rectangles-overlap(endpt0@pportion)(endpt1@pportion)(endpt0@qportion)(endpt1@qportion)))(start-hereworkloadparams)(if(flat-enough?pportionnormrtolatol)(if(flat-enough?qportionnormrtolatol)(linear-linearpportionqportionworkloadparams)(linear-quadraticpportionqportionworkloadparams#f))(if(flat-enough?qportionnormrtolatol)(linear-quadraticqportionpportionworkloadparams#t)(bisect-bothpportionqportionworkloadparams)))))))(define(linear-linearpportionqportionworkloadparams);; Both portions are flat enough to treat as linear. Treat;; them as line segments and see if the segments intersect.;; Include the intersection if they do.;; REMOVE THIS FOR A PRACTICAL APPLICATION.(display"linear-linear\n")(let-values(((tptq)(segment-parameters(endpt0@pportion)(endpt1@pportion)(endpt0@qportion)(endpt1@qportion))))(iftp(let((tp(lerptp(t0@pportion)(t1@pportion)))(tq(lerptq(t0@qportion)(t1@qportion))))(start-hereworkload(include-paramstptqparams)))(start-hereworkloadparams))))(define(linear-quadraticpportionqportionworkloadparamsreversed-roles?);; Only pportion is flat enough to treat as linear. Find its;; intersections with the quadratic qportion, and include;; either or both if they are within the correct;; intervals. (Use the qportion argument instead of Q;; directly, so the roles can be reversed for;; quadratic-linear.);;;; The following Maxima commands will supply formulas for;; finding values of t for a quadratic in s-power basis;; plugged into a line in s-power (or Bernstein) basis:;;;; /* The line. */;; xp(t) := px0*(1-t) + px1*t$;; yp(t) := py0*(1-t) + py1*t$;;;; /* The quadratic (s-power basis). */;; xq(t) := qx0*(1-t) + qx1*t*(1-t) + qx2*t$;; yq(t) := qy0*(1-t) + qy1*t*(1-t) + qy2*t$;;;; /* Implicitize and plug in. */;; impl(t) := resultant(xq(t)-xp(tau), yq(t)-yp(tau), tau)$;; expand(impl(t));;;;; REMOVE THIS FOR A PRACTICAL APPLICATION.(ifreversed-roles?(display"quadratic-linear\n")(display"linear-quadratic\n"))(let*((p0(endpt0@pportion))(p1(endpt1@pportion))(QX(plane-curve-x(curve@qportion)))(QY(plane-curve-y(curve@qportion)))(px0(carp0))(py0(cdrp0))(px1(carp1))(py1(cdrp1))(qx0(spower-c0QX))(qy0(spower-c0QY))(qx1(spower-c1QX))(qy1(spower-c1QY))(qx2(spower-c2QX))(qy2(spower-c2QY))(px0×py1(*px0py1))(px1×py0(*px1py0))(px0×qy0(*px0qy0))(px0×qy1(*px0qy1))(px0×qy2(*px0qy2))(px1×qy0(*px1qy0))(px1×qy1(*px1qy1))(px1×qy2(*px1qy2))(py0×qx0(*py0qx0))(py0×qx1(*py0qx1))(py0×qx2(*py0qx2))(py1×qx0(*py1qx0))(py1×qx1(*py1qx1))(py1×qx2(*py1qx2));; Construct the equation A×t² + B×t + C = 0 and solve;; it by the quadratic formula.(A(+px1×qy1(-px0×qy1)(-py1×qx1)py0×qx1))(B(+(-px1×qy2)px0×qy2(-px1×qy1)px0×qy1px1×qy0(-px0×qy0)py1×qx2(-py0×qx2)py1×qx1(-py0×qx1)(-py1×qx0)py0×qx0))(C(+(-px1×qy0)px0×qy0py1×qx0(-py0×qx0)(-px0×py1)px1×py0))(discr(-(*BB)(*4AC))))(define(invert-paramtq);;;; If one of the t-parameter solutions to the quadratic is;; in [0,1], invert the corresponding (x,y)-coordinates,;; to find the corresponding t-parameter solution for the;; linear. If it is in [0,1], then the (x,y)-coordinates;; are an intersection point. Return that corresponding;; t-parameter. Otherwise return #f.;;(and(0<=x<=1tq)(let((dx(-px1px0))(dy(-py1py0)))(if(>=(absdx)(absdy))(let*((x(spower-evalQXtq))(tp(/(-xpx0)dx)))(and(0<=x<=1tp)(lerptp(t0@pportion)(t1@pportion))))(let*((y(spower-evalQYtq))(tp(/(-ypy0)dy)))(and(0<=x<=1tp)(lerptp(t0@pportion)(t1@pportion))))))))(unless(negative?discr)(let*((rootd(sqrtdiscr))(tq1(/(+(-B)rootd)(+AA)))(tq2(/(-(-B)rootd)(+AA)))(tp1(invert-paramtq1))(tp2(invert-paramtq2)))(whentp1(set!params(ifreversed-roles?(include-paramstq1tp1params)(include-paramstp1tq1params))))(whentp2(set!params(ifreversed-roles?(include-paramstq2tp2params)(include-paramstp2tq2params))))))(start-hereworkloadparams)))(define(bisect-bothpportionqportionworkloadparams);; Neither portion is flat enough to treat as linear. Bisect;; them and add the four new portion-pairs to the workload.(let-values(((pport1pport2)(bisect-portionpportion))((qport1qport2)(bisect-portionqportion)))(start-here`((,pport1.,qport1)(,pport1.,qport2)(,pport2.,qport1)(,pport2.,qport2).,workload)params)))(start-here(initial-workloadPQ)'()))));; end library (spower quadratic intersections)(import(schemebase)(schemewrite)(spowerquadratic)(spowerquadraticintersections)(srfi132);; R7RS-large name: (scheme sort))(defineP(control-points->plane-curve'((-10)(010)(10))))(defineP1(control-points->plane-curve'((-10)(-1/25)(3350))))(defineP2(control-points->plane-curve'((05)(1/25)(10))))(defineQ(control-points->plane-curve'((21)(-82)(23))));; Sort the results by the parameters for P. Thus the intersection;; points will be sorted left to right.(defineparams(list-sort(lambda(tpair1tpair2)(<(cartpair1)(cartpair2)))(find-intersection-parametersPQ)))(for-each(lambda(pair)(let((tp(carpair))(tq(cdrpair)))(let-values(((axay)(plane-curve-evalPtp))((bxby)(plane-curve-evalQtq)))(display(inexacttp))(display"\t(")(display(inexactax))(display", ")(display(inexactay))(display")\t")(display(inexacttq))(display"\t(")(display(inexactbx))(display", ")(display(inexactby))(display")\n"))))params)(newline)(defineparams(list-sort(lambda(tpair1tpair2)(<(cartpair1)(cartpair2)))(find-intersection-parametersP1Q)))(for-each(lambda(pair)(let((tp(carpair))(tq(cdrpair)))(let-values(((axay)(plane-curve-evalP1tp))((bxby)(plane-curve-evalQtq)))(display(inexacttp))(display"\t(")(display(inexactax))(display", ")(display(inexactay))(display")\t")(display(inexacttq))(display"\t(")(display(inexactbx))(display", ")(display(inexactby))(display")\n"))))params)(newline)(defineparams(list-sort(lambda(tpair1tpair2)(<(cartpair1)(cartpair2)))(find-intersection-parametersQP1)))(for-each(lambda(pair)(let((tp(carpair))(tq(cdrpair)))(let-values(((axay)(plane-curve-evalQtp))((bxby)(plane-curve-evalP1tq)))(display(inexacttp))(display"\t(")(display(inexactax))(display", ")(display(inexactay))(display")\t")(display(inexacttq))(display"\t(")(display(inexactbx))(display", ")(display(inexactby))(display")\n"))))params)(newline)(defineparams(list-sort(lambda(tpair1tpair2)(<(cartpair1)(cartpair2)))(find-intersection-parametersP2Q)))(for-each(lambda(pair)(let((tp(carpair))(tq(cdrpair)))(let-values(((axay)(plane-curve-evalP2tp))((bxby)(plane-curve-evalQtq)))(display(inexacttp))(display"\t(")(display(inexactax))(display", ")(display(inexactay))(display")\t")(display(inexacttq))(display"\t(")(display(inexactbx))(display", ")(display(inexactby))(display")\n"))))params)
It turns out that the only bisection branch being executed, for the posed problem, was the one where both portions are bisected. Thus, in the new implementation here, the quadratic formula is never employed. However, fiddling with the numbers turns up problems for which this is not the case.
Also I have added a fourth problem, in which the first the convex-up parabola is cut in half. Interestingly, it takes more iterations to solvethat problem than the original! Note, though, that I have not categorized iterations according to their kinds.
$ gosh -r7 bezier_intersections_2.scmlinear-linearlinear-linearlinear-linearlinear-linearlinear-linearlinear-linear2170.07250828117222911 (-0.8549834376555417, 1.3450166066735616) 0.17250829973466453 (-0.8549837251463935, 1.345016599469329)0.15948753092173137 (-0.6810249381565372, 2.6810251680444233) 0.8405124690782686 (-0.6810251680444231, 2.681024938156537)0.8274917002653355 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542)0.9405124666929737 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525)quadratic-linearquadratic-linearquadratic-linearquadratic-linearquadratic-linearquadratic-linear4040.09485068481613748 (-0.6082597856508847, 1.3083729445649852) 0.15418647228249246 (-0.6082600809514531, 1.3083729445649848)0.1670105772077877 (0.0874641628839754, 2.7858070880490136) 0.892903544024507 (0.08746389814035438, 2.7858070880490144)linear-quadraticlinear-quadraticlinear-quadraticlinear-quadraticlinear-quadraticlinear-quadratic5910.15418647228249246 (-0.6082600809514531, 1.3083729445649848) 0.09485068481613748 (-0.6082597856508847, 1.3083729445649852)0.892903544024507 (0.08746389814035438, 2.7858070880490144) 0.1670105772077877 (0.0874641628839754, 2.7858070880490136)linear-linearlinear-linearlinear-linear7020.654983400530671 (0.654983400530671, 2.8549837251463934) 0.9274917188277709 (0.6549833933264384, 2.854983437655542)0.8810249333859473 (0.8810249333859473, 1.118975333761435) 0.059487533307026316 (0.881024666238565, 1.1189750666140525)
/* The control points of a planar quadratic Bézier curve form a triangle--called the "control polygon"--that completely contains the curve. Furthermore, the rectangle formed by the minimum and maximum x and y values of the control polygon completely contain the polygon, and therefore also the curve. Thus a simple method for narrowing down where intersections might be is: subdivide both curves until you find "small enough" regions where these rectangles overlap.*/import"./dynamic"forStructimport"./trait"forByRefimport"./math"forMath,Numsimport"./assert"forAssertimport"./seq"forStackimport"./fmt"forFmt// Note that these are all mutable as we need to pass by reference.varPoint=Struct.create("Point",["x","y"])varQuadSpline=Struct.create("QuadSpline",["c0","c1","c2"])// non-parametricvarQuadCurve=Struct.create("QuadCurve",["x","y"])// planar parametricvarWorkset=Struct.create("Workset",["p","q"])// Subdivision by de Casteljau's algorithmvarsubdivideQuadSpline=Fn.new{|q,t,u,v|vars=1-tvarc0=q.c0varc1=q.c1varc2=q.c2u.c0=c0v.c2=c2u.c1=s*c0+t*c1v.c1=s*c1+t*c2u.c2=s*u.c1+t*v.c1v.c0=u.c2}varsubdivideQuadCurve=Fn.new{|q,t,u,v|subdivideQuadSpline.call(q.x,t,u.x,v.x)subdivideQuadSpline.call(q.y,t,u.y,v.y)}// It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.varrectsOverlap=Fn.new{|xa0,ya0,xa1,ya1,xb0,yb0,xb1,yb1|return(xb0<=xa1&&xa0<=xb1&&yb0<=ya1&&ya0<=yb1)}// This accepts the point as an intersection if the boxes are small enough.vartestIntersect=Fn.new{|p,q,tol,exclude,accept,intersect|varpxmin=Nums.min([p.x.c0,p.x.c1,p.x.c2])varpymin=Nums.min([p.y.c0,p.y.c1,p.y.c2])varpxmax=Nums.max([p.x.c0,p.x.c1,p.x.c2])varpymax=Nums.max([p.y.c0,p.y.c1,p.y.c2])varqxmin=Nums.min([q.x.c0,q.x.c1,q.x.c2])varqymin=Nums.min([q.y.c0,q.y.c1,q.y.c2])varqxmax=Nums.max([q.x.c0,q.x.c1,q.x.c2])varqymax=Nums.max([q.y.c0,q.y.c1,q.y.c2])exclude.value=trueaccept.value=falseif(rectsOverlap.call(pxmin,pymin,pxmax,pymax,qxmin,qymin,qxmax,qymax)){exclude.value=falsevarxmin=Math.max(pxmin,qxmin)varxmax=Math.min(pxmax,qxmax)Assert.ok(xmax>=xmin)if(xmax-xmin<=tol){varymin=Math.max(pymin,qymin)varymax=Math.min(pymax,qymax)Assert.ok(ymax>=ymin)if(ymax-ymin<=tol){accept.value=trueintersect.x=0.5*xmin+0.5*xmaxintersect.y=0.5*ymin+0.5*ymax}}}}varseemsToBeDuplicate=Fn.new{|intersects,xy,spacing|varseemsToBeDup=falsevari=0while(!seemsToBeDup&&i!=intersects.count){varpt=intersects[i]seemsToBeDup=(pt.x-xy.x).abs<spacing&&(pt.y-xy.y).abs<spacingi=i+1}returnseemsToBeDup}varfindIntersects=Fn.new{|p,q,tol,spacing|varintersects=[]varworkload=Stack.new()workload.push(Workset.new(p,q))// Quit looking after having emptied the workload.while(!workload.isEmpty){varwork=workload.peek()workload.pop()varexclude=ByRef.new(false)varaccept=ByRef.new(false)varintersect=Point.new(0,0)testIntersect.call(work.p,work.q,tol,exclude,accept,intersect)if(accept.value){// To avoid detecting the same intersection twice, require some// space between intersections.if(!seemsToBeDuplicate.call(intersects,intersect,spacing)){intersects.add(intersect)}}elseif(!exclude.value){varp0=QuadCurve.new(QuadSpline.new(0,0,0),QuadSpline.new(0,0,0))varp1=QuadCurve.new(QuadSpline.new(0,0,0),QuadSpline.new(0,0,0))varq0=QuadCurve.new(QuadSpline.new(0,0,0),QuadSpline.new(0,0,0))varq1=QuadCurve.new(QuadSpline.new(0,0,0),QuadSpline.new(0,0,0))subdivideQuadCurve.call(work.p,0.5,p0,p1)subdivideQuadCurve.call(work.q,0.5,q0,q1)workload.push(Workset.new(p0,q0))workload.push(Workset.new(p0,q1))workload.push(Workset.new(p1,q0))workload.push(Workset.new(p1,q1))}}returnintersects}varp=QuadCurve.new(QuadSpline.new(-1,0,1),QuadSpline.new(0,10,0))varq=QuadCurve.new(QuadSpline.new(2,-8,2),QuadSpline.new(1,2,3))vartol=0.0000001varspacing=10*tolvarintersects=findIntersects.call(p,q,tol,spacing)for(intersectinintersects)Fmt.print("($ f, $f)",intersect.x,intersect.y)
( 0.654983, 2.854983)( 0.881025, 1.118975)(-0.681025, 2.681025)(-0.854983, 1.345017)