// math about cairo lines // This lib does not draw. It has computational aids. //manefest constants, tc style: SIZE_INT[return 4] INFINITY[return 2147483647] // NEW/MODIFIED LIB STUFF dot1 int x(1)[dot x(0),x(1)] // define pto=ptb offset by x,y pixels offset int pto(1), ptb(1),x,y [ pto(0)=ptb(0)+x; pto(1)=ptb(1)+y ] // define x given path a to b. rl is + for rt turn, - for left. // mag(rl_len)>1 implies trim to that length, else length is // equal to alen(a,b). corner90pt int ptx(1), pta(1), ptb(1), rl_len [ int rl rl=1; if(rl_len<0)rl= -1 _corner90 ptx, pta(0), pta(1), ptb(0), ptb(1), rl if(mag(rl_len)>1)[ int t(3), len setsegment t, ptb(0), ptb(1), ptx(0), ptx(1) len = mag(rl_len) trimto ptx, t, len ] ] xersectpt int ix(1), saf(1),sat(1), sbf(1), sbt(1) [ int res res = _xersect ix, saf(0),saf(1),sat(0),sat(1),sbf(0),sbf(1),sbt(0),sbt(1) if(res)[ pl"xersectpt: Lines are parallel" ] return res ] // draw a gentle arc from cp to tp with bow b at middle. bow is the // distance from straight line to the arc. arcbow int cp(1), tp(1), bow [ int m(1), a(1), cpa(1), tpa(1), tpac(1), cpac(1), ctr(1), rad // points: mid, arc, center, radius midpoint m,cp,tp // define midpoint m corner90pt a, cp,m, bow // define arcpoint a, rt turn distance bow midpoint cpa, cp, a // define cp->a midpoint, cpa midpoint tpa, a, tp // define a->tp midpoint, tpa corner90pt tpac,tp,tpa,+1 // points towards center, rt turns, corner90pt cpac,a,cpa,+1 // arbitrary distance // define center, intersection of extended tpac,cpac xersectpt ctr, tpa,tpac, cpa,cpac rad = alen ctr(0),ctr(1), a(0),a(1) next int s,e,dx,dy dx = tp(0)-ctr(0) dy = tp(1)-ctr(1) e = aarctan2 dy,dx dx = cp(0)-ctr(0) dy = cp(1)-ctr(1) s = aarctan2 dy,dx if(bow>0)arcneg ctr(0), ctr(1), rad, s,e else arc ctr(0), ctr(1), rad, s,e // showtext" <<== END, bow = ";shownum bow ] // end NEW/MODIFIED LIB STUFF // Turn 90 degree corners, +1 for right, -1 for left corner90 int ix(1), rl [ int x0, y0, x1, y1 x0=lastmoveto(0); y0=lastmoveto(1) x1=lastlineto(0); y1=lastlineto(1) return _corner90 ix, x0, y0, x1, y1, rl ] _corner90 int ix(1), x0, y0, x1, y1, rl [ int dbg;dbg=0 if(dbg)[pl"";pn x0;pn y0;pn x1;pn y1;pn rl;] int dx, dy, sdir, sdir, ndir, dtab(3), sdtab, temp // given a seg get two deltas, dx dy. dx=x1-x0; dy=y1-y0 if(dbg)[pl"dx,dy";pn dx;pn dy] // compute signdir from their respective signs sdir = 0; if(dx>0)sdir=2; if(dy>0)sdir=sdir+1 if(dbg)[pl"sdir";pn sdir] // from sign map to ndir, ==>> Ar*/notes.txt ~9558 setsegment dtab, 2,1,3,0 // FIX, this is a table, not a seg if(dbg)[pl"dtab";psegment dtab] // from the rl arg compute new ndir (ndir +- 1)%4 ndir = (dtab(sdir)+rl+4)%4 // new direction if(dbg)[pl"ndir";pn ndir] // map back to signdir setsegment dtab, 3,2,0,1 //2nd FIX sdir = dtab(ndir) if(dbg)[pl"sdir";pn sdir] // flip the MAGNITUDEs of dx and dy if(dx>0)temp = dx; else temp=-dx if(dy>0)dx = dy; else dx = -dy dy = temp if(dbg)[pl"flipped mags of dx,dy";pn dx;pn dy] // apply the signs to dx and dy if(sdir%2==0)dx = -dx if(sdir/2==0)dy = -dy if(dbg)[pl"signed dx,dy";pn dx;pn dy] // set ix to point 1 plus new delta ix(0) = dx+x1 ix(1) = dy+y1 if(dbg)[pl"ix";ppoint(ix)] // return length return alen(x0,y0,x1,y1) ] // convenience routines rt90 int ix(1)[ return corner90 ix, +1 ] lf90 int ix(1)[ return corner90 ix, -1 ] // Extend a line segment // computes a point ix that would extend a line segment // to length len. Positive len extends from the lines end, // negative from the lines beginning. seg is not modified. extend int ix(1), seg(3), len [ int dbg;dbg=0 if(dbg)[ pl" ~150 in extend: " psegment seg;ps" len";pn len ] int slen, fac slen = seglen(seg) if(len>0)fac = 1+len/slen else fac = -1+len/slen if(dbg)[ pl"~159 len,fac";pn len;pn fac ] _extend ix,seg,fac if(fac>0)[ seg(2)=ix(0) seg(3)=ix(1) if(dbg)[pl"~165 new last point"; ppoint seg+2*SIZE_INT] ] else if(fac<0)[ seg(0)=ix(0) seg(1)=ix(1) if(dbg)[pl"~159 new first point"; ppoint seg] ] slen = seglen seg if(dbg)[ pl"new seg";psegment seg;ps" slen";pn slen pl"";pl"~172 calling trim"; ] if(dbg)[pl"~178: len slen";pn len;pn slen] if(len>0)[ trim ix,seg,slen-len seg(2)=ix(0) seg(3)=ix(1) if(dbg)[pl"lineLib~99 new last point"; ppoint seg+2*SIZE_INT] ] else if(len<0)[ trim ix,seg,-slen-len seg(0)=ix(0) seg(1)=ix(1) if(dbg)[pl"lineLib~104 new first point"; ppoint seg] ] ] // compute point ix that would extend seg's length by a factor, // fac. Positive fac extends the lines end, negative the // beginning. seg is not modified. _extend int ix(1), seg(3), fac [ int dbg;dbg=0 if(dbg)[pl" ~7(ex): fac";pn fac] int mfac, dx, dy if(fac>0)mfac=fac; else mfac=-fac if(dbg)[pl" ~9(ex): mfac";pn mfac] dx = (seg(2)-seg(0))*mfac dy = (seg(3)-seg(1))*mfac if(dbg)[pl" ~11: dx,dy";pn dx;pn dy] if(fac>0)[ ix(0) = seg(0)+dx ix(1) = seg(1)+dy ] else [ ix(0) = seg(2)-dx ix(1) = seg(3)-dy ] if(dbg)[pl" ~21: ix ";ppoint ix] if(dbg)[pl" _extend done";pl""] ] // if length of seg>mag(len) return in ix a new point that would // trim seg's length by mag(len). Trim from seg(2..3) if len>0, // else seg(0..1). seg is not modified. trim int ix(1),seg(3),len [ int dbg;dbg=0 if(dbg)[pl"~137 trim seg len "; psegment seg; pn len ] int slen, amt, dx, dy dx = seg(2)-seg(0) dy = seg(3)-seg(1) if(dbg)[pl" ~32(tr): dx,dy";pn dx;pn dy] slen = seglen(seg) if(len>0)[ amt = (1000*len+500)/slen if(dbg)[pl" ~36: len,slen,amt"; pn len;pn slen;pn amt ] ix(0) = seg(2) - amt*dx/1000 ix(1) = seg(3) - amt*dy/1000 ] else[ amt = (-1000*len-500)/slen if(dbg)[pl" len,slen,amt";pn len;pn slen;pn amt] ix(0) = seg(0) + amt*dx/1000 ix(1) = seg(1) + amt*dy/1000 ] ] // trims the seg TO lenght len trimto int ix(1),seg(3),to [ int amt amt = seglen(seg) - to trim ix,seg,amt ] // sets a segment array, seg, with two points setsegment int seg(3), ax,ay,bx,by [ seg(0)=ax; seg(1)=ay; seg(2)=bx; seg(3)=by; ] // sets a point array, p, with x,y setpoint int p(1), ax,ay[ p(0)=ax; p(1)=ay ] // prints a segment in brackets psegment int seg(3) [ ps"["; ppoint seg; ppoint seg+2*SIZE_INT; ps"]" ] // prints a point in brackets ppoint int pt(1) [ ps"["; pn pt(0); pn pt(1); ps"]" ] // returns the length of a segment seglen int seg(3) [ return alen seg(0),seg(1),seg(2),seg(3) ] // intersection of two lines // args: answer, line a, line b // returns: -1 for parallel lines, 0 for OK xersectseg int ix(1), sa(3), sb(3) [ return _xersect ix, sa(0),sa(1),sa(2),sa(3),sb(0),sb(1),sb(2),sb(3) ] _xersect int ix(1), ax0,ay0,ax1,ay1, bx0,by0,bx1,by1 [ int xdb;xdb=0 if(xdb)[pl"line a ="; pn ax0;pn ay0;pn ax1;pn ay1] if(xdb)[pl"line b ="; pn bx0;pn by0;pn bx1;pn by1] int day,dax,dby,dbx,f,e,h,g,y,x,avert,bvert dax = ax1-ax0 day = ay1-ay0 dbx = bx1-bx0 dby = by1-by0 if(xdb)[pl"dax,day,dbx,dby";pn dax;pn day;pn dbx;pn dby] if(dax==0) [ // line a is vertical if(dbx==0) return -1 // parallel x = ax0 avert=1 if(xdb)[pl"a is vertical at x ="; pn x] f=INFINITY // a slope hopefully never encountered ] else [ f = day*1000/dax //slope e = ay0 - f*ax0/1000 // y=0 intersept if(xdb)[pl"a: y=fx+e f,e="; pn f;pn e] ] if(dbx==0) [ // line b is vertical x = bx0 bvert=1 if(xdb)[pl"b is vertical at x ="; pn x] ] else [ h = dby*1000/dbx g = by0 - h*bx0/1000 if(xdb)[pl"b: y=hx+g ="; pn h;pn g] if(h==f)return -1; // parallel ] if( 0==(avert+bvert) ) x = ((e-g)*1000)/(h-f) if(xdb)[pl"x ="; pn x] if( 0==bvert ) y = g + h*x/1000 else y = e + f*x/1000 if(xdb)[pl"bvert y ="; pn bvert;pn y] ix(0) = x ix(1) = y if(xdb)[pl"intersect is at "; ppoint ix] return 0 ] bisectseg int bsect(1), icntr(1), sega(3),segb(3) [ int xret, segin(3),segout(3), ptout(1),ptin(1) xret = _xersect icntr,sega(0),sega(1),sega(2),sega(3),segb(0),segb(1),segb(2),segb(3) if(xret)return xret // parallel, no angle bisect bsect, sega(0),sega(1),icntr(0),icntr(1),segb(2),segb(3) ] // see SamplePrograms/testBisect.tc for comments and embedded debug code bisect int bsect(1), b0,b1,m0,m1,e0,e1 [ int segin(3), segout(3), ptout(1),ptin(1),inoutlen setsegment segin, m0,m1,b0,b1 //both segs middle outwards setsegment segout, m0,m1,e0,e1 if(parallel(segin,segout)) return -1 extend ptin, segin, 500 extend ptout, segout, 500 inoutlen = alen(ptin(0),ptin(1),ptout(0),ptout(1)) if(inoutlen<=800)[ bsect(0) = (ptout(0)+ptin(0))/2 bsect(1) = (ptout(1)+ptin(1))/2 return 0 ] else [ reverseseg segin; reverseseg segout int in90(1), out90(1) _corner90 in90, segin(0), segin(1), segin(2), segin(3), -1 _corner90 out90, segout(0), segout(1), segout(2), segout(3), 1 return bisect(bsect,in90(0),in90(1),m0,m1,out90(0),out90(1)) ] ] // assign p1=p2 assignpt int p1(1), p2(1)[ p1(0)=p2(0) p1(1)=p2(1) ] // swap two points swappt int p1(1), p2(1)[ int temp(1) assignpt temp,p1 assignpt p1,p2 assignpt p2,temp ] // swap first and last points of s reverseseg int s(3)[ swappt s, s+2*SIZE_INT ] // set n=nearest and f=furthest from ii points of seg // return 0 if ii is between n and f, else +1 for nearest // is segs end or -1 if its beginning nearfar int n(1), f(1), ii(1), seg(3)[ int len1, len2 len1 = alen ii(0),ii(1),seg(0),seg(1) len2 = alen ii(0),ii(1),seg(2),seg(3) if(len1>len2)[ assign f,seg assign n,seg+2*SIZE_INT which = +1 ] else [ assign n,seg assign f,seg+2*SIZE_INT which = -1 ] if(len1+len2==(seglen(seg)) return 0 // NOT QUITE ROBUST else return which ] // test that three pts are roughly colinear. Return 0 if not, // else 1,2,3 to indicate which is in the middle. Return -1 if // points are tightly clustered. NEEDS TESTING. colinear int p1(1),p2(1),p3(1)[ int d12, d13, d23, dlargest d12=distance p1,p2 d13=distance p1,p3 d23=distance p2,p3 dlargest=d12; middle=3; dsum=d13+d23 if(d13>dlargest)[ dlargest = d13; middle = 2; dsum=d12+d23 ] if(d23>dlargest)[ dlargest = d23; middle = 1; dsum=d12+d13 ] if(dlargest<2) return -1 // tightly clustered ddiff = dsum-dlargest if((-1