(**************************************************************************) (* Last modified 6/11/98 *) (**************************************************************************) (* A rational version of the de Casteljau algorithm. *) (* 11/25/2002 *) (* Control polygon cpoly = (b0,...,bm) *) (* Each control point is either a weighed point {x1,...xn,w}, w <> 0, * or a control vector {x1,...xn,0}. * * * To handle poly curves, set w = 1 * * * * Affine frame on real line: (r, s) * pdecas[cpoly, r, s, t]: Computation of the current point F(t) * To display pdecas, use showpic2D, or showpic3D * gsubdiv[cpoly, r, s, l, m, n]: perform n steps of subdivision. * starting with t = l, and ending with t = m. First, recomputes new control * polygon, and then use interval [l, m], and subdivide with $t = (l + m)/2. * To display gsubdiv, use subdiv2D, or subdiv3D * subdiv[cpoly, r, s, n]: perform n steps of subdivision (using t = (r + 2)/2) * To display subdiv, use showsub2D, or showsub3D * To display subdiv with flagpol, use showsub2, or showsub3 * ggdecas[cpoly, tt, r, s]: computes the polar value f(t1,...,tm), * where tt = {t1,...,tm}, using the control polygon cpoly with m edges * To display ggdecas, use showc2D, or showc3D * showclo2D[cpoly, r, s, lx, mx, ly, my, n] will render * a closed rational curve, performing n subdivision steps * This program first computes a new control polygon wrt [-1, 1], * and then computes the opposite control polygon (In the hat space) * where oppoly[[i]] = (-1)^(m - i) poly[[i]] * showclo3D[cpoly, r, s, lx, mx, ly, my, lz, mz, n] is the 3D version * rendclo2D[cpoly, r, s, lx, mx, ly, my, flagpol, n] will render * a closed rational curve, performing n subdivision steps * first computes new control poly wrt [-1, 1] * If flagpol =!= 1, do not display control polygon * If flagpol === -1, show original arc in black, and dual arc in red * rendlo3D[cpoly, r, s, lx, mx, ly, my, lz, mz, flagpol, n] is the 3D version * rendcloz2D[cpoly, r, s, lx, mx, ly, my, flagpol, n] and * rendcloz3D[cpoly, r, s, lx, mx, ly, my, flagpol, n] work like * rendclo2D and rendclo3D, but with to the original control poly wrt [0, 1] ***************************************************************************) (* Performs general affine interpolation between two points p, q *) (* w.r.t. affine basis [r, s], and interpolating value t *) lerp[p_List,q_List,r_,s_,t_] := (s - t)/(s - r) p + (t - r)/(s - r) q; olerp[{p1__},{p2__},r_,s_,t_] := Block[ {p = {p1}, q = {p2}, res}, res = (s - t)/(s - r) p + (t - r)/(s - r) q; res ]; (* computes the point F(t) and the line segments involved in computing F(t) using the de Casteljau algorithm *) decas[{cpoly__}, r_, s_, t_] := Block[ {bb = {cpoly}, b = {}, m, i, j, lseg = {}, res}, (m = Length[bb] - 1; Do[ Do[ b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t]]; If[i > 1, lseg = Append[lseg, {b[[i - 1]], b[[i]]}]] , {i, 1, m - j + 1} ]; bb = b; b = {}, {j, 1, m} ]; res := Append[lseg, bb[[1]]]; res ) ]; (* computes a point F(t) on a curve using the de Casteljau algorithm *) (* this is the simplest version of de Casteljau *) (* the auxiliary points involved in the algorithm are not computed *) badecas[{cpoly__}, r_, s_, t_] := Block[ {bb = {cpoly}, b = {}, m, i, j}, (m = Length[bb] - 1; Do[ Do[ b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t]], {i, 1, m - j + 1} ]; bb = b; b = {}, {j, 1, m} ]; bb[[1]] ) ]; (* computes a polar value f(t1, ..., tm) *) (* using the de Casteljau algorithm *) (* the input is a list tt of m numbers *) (* the auxiliary points involved in the algorithm are not computed *) bdecas[{cpoly__}, {tt__}, r_, s_] := Block[ {bb = {cpoly}, b = {}, t = {tt}, m, i, j, lseg = {}}, (m = Length[bb] - 1; Do[ Do[ b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t[[j]]]], {i, 1, m - j + 1} ]; bb = b; b = {}, {j, 1, m} ]; bb[[1]] ) ]; (* computes a polar value using the de Casteljau algorithm *) (* and the line segments involved in computing f(t1, ..., tm) *) (* the input is a list tt of m numbers *) gdecas[{cpoly__}, {tt__}, r_, s_] := Block[ {bb = {cpoly}, b = {}, t = {tt}, m, i, j, lseg = {}, res}, (m = Length[bb] - 1; Do[ Do[ b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t[[j]]]]; If[i > 1, lseg = Append[lseg, {b[[i - 1]], b[[i]]}]] , {i, 1, m - j + 1} ]; bb = b; b = {}, {j, 1, m} ]; res := Append[lseg, bb[[1]]]; res ) ]; (* Maps a net to the hat space by multiplying all coords except *) (* the weight, by the weight (for each point) *) prepare[{poly__}] := Block[ {lpoly = {poly}, newp = {}, pt, h, w, i, l1}, l1 = Length[lpoly]; Do[ pt = lpoly[[i]]; w = Last[pt]; h = Drop[pt, -1]; If[w =!= 0, h = w * h; pt = Append[h, w] ]; newp = Append[newp,pt], {i, 1, l1} ]; newp ]; (* To project a point in the hat space back onto the affine space *) aproj[{poly__}] := Block[ {pt = {poly}, res, h, w}, w = Last[pt]; h = Drop[pt, -1]; If[Abs[w] > 10^(-20), If[Abs[w] <= 10^(-16), Print["Warning: Point with very small weight: ", pt]]; h = h/w, h = 10^(20) * h; Print["*** Warning: Point at infinity!: ", pt, " ***"] ]; h ]; (* To project a list of points in the hat space back onto the affine space *) (* points in the hat space which are almost the zero vector are discarded *) proj[{poly__}] := Block[ {sl = {poly}, np = {}, pt, h, j, l2, flag}, l2 = Length[sl]; Do[ pt = sl[[j]]; flag = nearzero[pt]; (* If pt is almost the zero vector, discard *) If[flag =!= 1, h = aproj[pt]; np = Append[np,h]], {j, 1, l2} ]; np ]; (* To project a list of control polygons in the hat space *) (* back onto the affine space *) projlis[{netlis__}] := Block[ {slis = {netlis}, newlis = {}, anet, newnet, j, l2}, l2 = Length[slis]; Do[ anet = slis[[j]]; newnet = proj[anet]; newlis = Append[newlis,newnet], {j, 1, l2} ]; newlis ]; (* this function calls decas, projects down to the affine space *) (* and computes the line segments in Mathematica, with colors *) pdecas[{cpoly__}, r_, s_, t_] := Block[ {bb = {cpoly}, pt, ll, res, i, l1, edge, lt, rt}, res = decas[bb, r, s, t]; pt = Last[res]; res = Drop[res, -1]; l1 = Length[res]; ll = {}; Do[ edge = res[[i]]; lt = edge[[1]]; rt = edge[[2]]; lt = aproj[lt]; rt = aproj[rt]; edge = {lt, rt}; ll = Append[ll, Line[edge]], {i, 1, l1} ]; pt = aproj[pt]; res = Append[ll, {RGBColor[1,0,0], PointSize[0.01], Point[pt]}]; res ]; (* Computes a polar value, and the line segments involved in the computation, as Mathematica Line[] elements Also projects down from the hat space *) ggdecas[{cpoly__}, {tt__}, r_, s_] := Block[ {bb, pt, ll, res, i, l1, edge, lt, rt}, bb = {cpoly}; res= gdecas[bb, {tt}, r, s]; pt = Last[res]; res = Drop[res, -1]; l1 = Length[res]; ll = {}; Do[ edge = res[[i]]; lt = edge[[1]]; rt = edge[[2]]; lt = aproj[lt]; rt = aproj[rt]; edge = {lt, rt}; ll = Append[ll, Line[edge]], {i, 1, l1} ]; pt = aproj[pt]; res = Append[ll, {PointSize[0.01], Point[pt]}]; res ]; (* Computes new control net by computing polar values *) (* f(l,..,l,m,...,m) *) (* This version is inefficient *) oldnewcpoly[{cpoly__}, r_, s_, l_, m_] := Block[ {bb, npoly = {}, pt, i, j, l1, tt}, bb = {cpoly}; l1 = Length[bb]; Do[ tt = {}; Do[tt = Append[tt, l], {i, 1, l1 - j} ]; Do[tt = Append[tt, m], {i, l1 - j + 1, l1 - 1} ]; pt = bdecas[bb, tt, r, s]; npoly = Append[npoly, pt], {j, 1, l1} ]; npoly ]; (* computes the control polygon for display. *) (* Control vectors cause extra trouble *) controlfig2D[{poly__}] := Block[ {p = {poly}, polylin, cpoly, cvec, cpt, conpt, res, totnum, nodenum, ptnum, vecnum, w, i, n, pt, orig = 1}, n := Length[p]; polylin = {}; cpoly = {}; cvec = {}; cpt = {}; conpt = {}; totnum = 0; nodenum = 0; ptnum = 0; vecnum = 0; Do[w := Last[p[[i]]]; If[w === 0, pt = Drop[p[[i]],-1] + orig; cvec = Append[cvec, Line[{{orig, orig}, pt}]]; pt = pt + {0.1,0}; cvec = Append[cvec, Text["v", pt, {-1,0}]]; vecnum = vecnum + 1; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}] , pt = Drop[p[[i]],-1]; cpt = Append[cpt, Point[pt]]; ptnum = ptnum + 1; cpoly = Append[cpoly, pt]; nodenum = nodenum + 1 ], {i, 1, n} ]; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}]; If[cvec =!= {}, cvec = Append[cvec, Point[{orig,orig}]]]; conpt = Prepend[cpt, PointSize[0.01]]; res = Join[conpt,cvec,polylin]; res = Prepend[res, RGBColor[0,1,0]]; res ]; (* rfig2 works as controlfig2D, except that it keeps track of the list *) (* of weights in wll, and appends it as second argument *) rfig2D[{poly__}] := Block[ {p = {poly}, polylin, cpoly, cvec, cpt, conpt, res, totnum, nodenum, ptnum, vecnum, w, i, n, pt, wll, orig = 1}, n := Length[p]; polylin = {}; cpoly = {}; cvec = {}; cpt = {}; conpt = {}; wll = {}; totnum = 0; nodenum = 0; ptnum = 0; vecnum = 0; cww = {}; Do[w := Last[p[[i]]]; wll = Append[wll, w]; If[w === 0, pt = Drop[p[[i]],-1] + orig; cvec = Append[cvec, Line[{{orig, orig}, pt}]]; pt = pt + {0.1,0}; cvec = Append[cvec, Text["v", pt, {-1,0}]]; vecnum = vecnum + 1; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}] , pt = Drop[p[[i]],-1]; cpt = Append[cpt, Point[pt]]; ptnum = ptnum + 1; cww = Append[cww, Text[w, pt + {0.7, 0}, {1,0}]]; cpoly = Append[cpoly, pt]; nodenum = nodenum + 1 ], {i, 1, n} ]; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}]; If[cvec =!= {}, cvec = Append[cvec, Point[{orig,orig}]]]; conpt = Prepend[cpt, PointSize[0.01]]; res = Join[conpt,cww,cvec,polylin]; res = Prepend[res, RGBColor[1,0,0]]; res = Append[res, wll]; res ]; (* 3D version of controlfig2D *) controlfig3D[{poly__}] := Block[ {p = {poly}, polylin, cpoly, cvec, cpt, conpt, res, totnum, nodenum, ptnum, vecnum, w, i, n, pt, orig = 1}, n := Length[p]; polylin = {}; cpoly = {}; cvec = {}; cpt = {}; conpt = {}; totnum = 0; nodenum = 0; ptnum = 0; vecnum = 0; Do[w := Last[p[[i]]]; If[w === 0, pt = Drop[p[[i]],-1] + orig; cvec = Append[cvec, Line[{{orig, orig, orig}, pt}]]; pt = pt + {0.1,0,0}; cvec = Append[cvec, Text["v", pt, {-1,0}]]; vecnum = vecnum + 1; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}] , pt = Drop[p[[i]],-1]; cpt = Append[cpt, Point[pt]]; ptnum = ptnum + 1; cpoly = Append[cpoly, pt]; nodenum = nodenum + 1 ], {i, 1, n} ]; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}]; If[cvec =!= {}, cvec = Append[cvec, Point[{orig, orig, orig}]]]; conpt = Prepend[cpt, PointSize[0.01]]; res = Join[conpt,cvec,polylin]; res = Prepend[res, RGBColor[0,1,0]]; res ]; (* 3D version of rfig2D *) rfig3D[{poly__}] := Block[ {p = {poly}, polylin, cpoly, cvec, cpt, conpt, res, totnum, nodenum, ptnum, vecnum, w, i, n, pt, wll, orig = 1}, n := Length[p]; polylin = {}; cpoly = {}; cvec = {}; cpt = {}; conpt = {}; wll = {}; totnum = 0; nodenum = 0; ptnum = 0; vecnum = 0; Do[w := Last[p[[i]]]; wll = Append[wll, w]; If[w === 0, pt = Drop[p[[i]],-1] + orig; cvec = Append[cvec, Line[{{orig, orig, orig}, pt}]]; pt = pt + {0.1,0,0}; cvec = Append[cvec, Text["v", pt, {-1,0}]]; vecnum = vecnum + 1; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}] , pt = Drop[p[[i]],-1]; cpt = Append[cpt, Point[pt]]; ptnum = ptnum + 1; cpoly = Append[cpoly, pt]; nodenum = nodenum + 1 ], {i, 1, n} ]; If[nodenum > 1, polylin = Append[polylin, Line[cpoly]]; totnum = totnum + 1; nodenum = 0; cpoly = {} , nodenum = 0; cpoly = {}]; If[cvec =!= {}, cvec = Append[cvec, Point[{orig, orig, orig}]]]; conpt = Prepend[cpt, PointSize[0.01]]; res = Join[conpt,cvec,polylin]; res = Prepend[res, RGBColor[1,0,0]]; res = Append[res,wll]; res ]; (* projects control polygon in hat space back onto affine space *) (* keeps the weights *) newpoly[{cpoly__}] := Block[ {res = {cpoly}, l1, i, npt, pt, newpt, w}, l1 = Length[res]; newp = {}; Do[ w = Last[res[[i]]]; pt = Drop[res[[i]],-1]; If[w =!= 0, pt = pt/w; npt = Append[pt, w], npt = res[[i]]]; newp = Append[newp, npt], {i, 1, l1} ]; newp ]; (* To display the construction of a point on the curve *) (* and the line segments involved in its construction + *) (* control polygon *) showpic2D[{poly__},r_,s_,lx_,mx_,ly_,my_,t_] := Block[ {pol1}, pol1 = prepare[{poly}]; Show[ Graphics[{controlfig2D[{poly}], {Thickness[0.001], pdecas[pol1, r, s, t]}}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction] ]; (* To display the construction of a polar value *) (* and the line segments involved in its construction + *) (* control polygon *) showc2D[{poly__},{tt__},r_,s_,lx_,mx_,ly_,my_] := Block[ {pol1}, pol1 = prepare[{poly}]; Show[ Graphics[{controlfig2D[{poly}], {Thickness[0.001], ggdecas[pol1, {tt}, r, s]}}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction] ]; (* 3D version of showpic2D *) showpic3D[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,t_] := Block[ {pol1}, pol1 = prepare[{poly}]; Show[ Graphics3D[{controlfig3D[{poly}], {Thickness[0.001], pdecas[pol1, r, s, t]}}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* 3D version of showc2D *) showc3D[{poly__}, {tt__}, r_,s_,lx_,mx_,ly_,my_,lz_,mz_] := Block[ {pol1}, pol1 = prepare[{poly}]; Show[ Graphics3D[{controlfig3D[{poly}], {Thickness[0.001], ggdecas[pol1, {tt}, r, s]}}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* Performs a single subdivision step using the de Casteljau algorithm *) (* Returns the control poly (f(r,...,r, t, ..., t)) and *) (* (f(t, ..., t, s, ..., s)) *) subdecas[{cpoly__}, r_, s_, t_] := Block[ {bb = {cpoly}, b = {}, ud = {}, ld = {}, m, i, j, lseg = {}, res}, (m = Length[bb] - 1; ud = {bb[[1]]}; ld = {bb[[m + 1]]}; Do[ Do[ b = Append[b, lerp[bb[[i]], bb[[i+1]], r, s, t]], {i, 1, m - j + 1} ]; ud = Append[ud, b[[1]]]; ld = Prepend[ld, b[[m - j + 1]]]; bb = b; b = {}, {j, 1, m} ]; res := Join[{ud},{ld}]; res ) ]; (* Computes the control polygon wrt new affine frame (a, b) *) (* Assumes that a = (1 - lambda) r + lambda t and *) (* that b = (1 - mu) r + mu t, wrt original frame (s, t) *) (* Returns control poly (f(a, ..., a, b, ..., b)) *) (* This version is much faster than oldnewcpoly *) newcpoly[{cpoly__}, r_, s_, a_, b_] := Block[ {poly = {cpoly}, m, i, pol1, pola, pol2, npoly, pt}, ( If[a =!= r, pol1 = subdecas[poly, r, s, a]; pola = pol1[[1]]; pol2 = {}; m = Length[pola]; Do[ pt = pola[[i]]; pol2 = Prepend[pol2, pt], {i, 1, m} ]; npoly = subdecas[pol2, a, r, b], (* Print[" npoly: ", npoly] *) npoly= subdecas[poly, r, s, b] ]; npoly[[1]] ) ]; (* subdivides each control polygon in a list of control polygons *) (* using subdecas. Uses t = (r + s)/2 *) subdivstep[{poly__}, r_, s_] := Block[ {cpoly = {poly}, lpoly = {}, t, l, i}, (l = Length[cpoly]; t = (r + s)/2; Do[ lpoly = Join[lpoly, subdecas[cpoly[[i]], r, s, t]] , {i, 1, l} ]; lpoly ) ]; (* To create a list of line segments from a list of control polygons *) makedgelis[{poly__}] := Block[ {res, sl, newsl = {poly}, i, j, l1, l2}, (l1 = Length[newsl]; res = {}; Do[ sl = newsl[[i]]; l2 = Length[sl]; Do[If[j > 1, res = Append[res, Line[{sl[[j-1]], sl[[j]]}]]], {j, 1, l2} ], {i, 1, l1} ]; res ) ]; (* calls subdivstep n times *) (* old version also projecting down to the hat space *) oldsubdiv[{poly__}, r_, s_, n_] := Block[ {res = {}, pol1 = {poly}, newp = {}, newsl = {}, i, j}, ( newp = {pol1}; Do[ newp = subdivstep[newp, r, s], {i, 1, n} ]; newsl = projlis[newp]; res = makedgelis[newsl]; res ) ]; (* calls subdivstep n times *) subdiv[{poly__}, r_, s_, n_] := Block[ {pol1 = {poly}, newp = {}, i}, ( newp = {pol1}; Do[ newp = subdivstep[newp, r, s], {i, 1, n} ]; newp ) ]; (* This old version of gsubdiv also projects down from the hat space *) oldgsubdiv[{poly__}, r_, s_, l_, m_, n_] := Block[ {res = {}, pol1 = {poly}, newp = {}, newsl = {}, i}, ( If[(l =!= r) || (m =!= s), newp = {newcpoly[pol1, r, s, l, m]} , newp = {pol1} ]; Do[ newp = subdivstep[newp, l, m], {i, 1, n} ]; newsl = projlis[newp]; res = makedgelis[newsl]; res ) ]; (* calls subdivstep n times, but also recomputes control net *) (* wrt the affine base (m, n) *) gsubdiv[{poly__}, r_, s_, l_, m_, n_] := Block[ {pol1 = {poly}, newp = {}, newsl = {}, i}, ( If[(l =!= r) || (m =!= s), newp = {newcpoly[pol1, r, s, l, m]} , newp = {pol1} ]; Do[ newp = subdivstep[newp, l, m], {i, 1, n} ]; newp ) ]; (* Projects down the list of control polygons from the hat space *) (* onto the affine space, and creates the polyline. *) tolineseg[{poly__}] := Block[ {pol = {poly}, res, newsl}, ( newsl = projlis[pol]; res = makedgelis[newsl]; res ) ]; (* To display the result of n steps of subdivision *) (* curve segment over [r, s] *) showsub2D[{poly__},r_,s_,lx_,mx_,ly_,my_,n_] := Block[ {pol1, curv}, pol1 = prepare[{poly}]; curv = subdiv[pol1, r, s, n]; curv = tolineseg[curv]; Show[ Graphics[{controlfig2D[{poly}], curv}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction] ]; (* To display the result of n steps of subdivision *) (* curve segment over [r, s] *) (* If flagpol =!= 1, do not display control polygon *) showsub2[{poly__},r_,s_,lx_,mx_,ly_,my_,flagpol_,n_] := Block[ {pol1, curv, image}, pol1 = prepare[{poly}]; curv = subdiv[pol1, r, s, n]; curv = tolineseg[curv]; If[flagpol =!= 1, image = {curv}, image = {controlfig2D[{poly}], curv} ]; Show[ Graphics[image], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction] ]; (* 3D version of showsub2D *) showsub3D[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,n_] := Block[ {pol1, curv}, pol1 = prepare[{poly}]; curv = subdiv[pol1, r, s, n]; curv = tolineseg[curv]; Show[ Graphics3D[{controlfig3D[{poly}], curv}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* 3D version of showsubD (with flagpol *) showsub3[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,flagpol_,n_] := Block[ {pol1, curv, image}, pol1 = prepare[{poly}]; curv = subdiv[pol1, r, s, n]; curv = tolineseg[curv]; If[flagpol =!= 1, image = {curv}, image = {controlfig3D[{poly}], curv} ]; Show[ Graphics3D[image], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction] ]; (* To display the result of n steps of subdivision *) (* Also recomputes control net to display curve segment over [l, m] *) subdiv2D[{poly__},r_,s_,l_,m_,lx_,mx_,ly_,my_,n_] := Block[ {cpoly = {poly}, np, wl, curv, pol1, pt, ipt, lpt, ttt, ppt}, ttt = {}; pol1 = prepare[cpoly]; np = newcpoly[pol1, r, s, l, m]; np = newpoly[np]; ipt = First[np]; ipt = aproj[ipt]; lpt = Last[np]; lpt = aproj[lpt]; np = rfig2D[np]; wl = Last[np]; np = Drop[np,-1]; ttt = Append[ttt, Text["F(" <> ToString[l] <> ")", ipt + {1, -0.5}, {1,0}]]; ttt = Append[ttt, Text["F(" <> ToString[m] <> ")", lpt + {1, -0.5}, {1,0}]]; Print["Weights: ", wl]; pt = badecas[pol1, l, m, (l + m)/2]; pt = aproj[pt]; (* ppt = InputForm[(l+m)/2]; *) ppt = oo; ttt = Append[ttt, Text["F(" <> ToString[ppt] <> ")", pt + {1, -0.5}, {1,0}]]; curv = gsubdiv[pol1, r, s, l, m, n]; curv = tolineseg[curv]; Print["Ready to display! "]; Show[ Graphics[{controlfig2D[{poly}], Point[pt], ttt, np, curv}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; ]; (* To display the result of n steps of subdivision *) (* Also recomputes control net to display curve segment over [l, m] *) (* version for paper, shows F((1 + m)/2) and F(infty) *) subdiv2Db[{poly__},r_,s_,l_,m_,lx_,mx_,ly_,my_,n_,ab_] := Block[ {cpoly = {poly}, np, wl, curv, pol1, pt, ipt, lpt, ttt, ppt}, ttt = {}; pol1 = prepare[cpoly]; np = newcpoly[pol1, r, s, l, m]; np = newpoly[np]; ipt = First[np]; ipt = aproj[ipt]; lpt = Last[np]; lpt = aproj[lpt]; np = rfig2D[np]; wl = Last[np]; np = Drop[np,-1]; ttt = Append[ttt, Text["F(" <> ToString[l] <> ")", ipt + {1, -0.5}, {1,0}]]; ttt = Append[ttt, Text["F(" <> ToString[m] <> ")", lpt + {1, -0.5}, {1,0}]]; Print["Weights: ", wl]; pt = badecas[pol1, l, m, (l + m)/2]; pt = aproj[pt]; If[ab === 1, ppt = InputForm[(l+m)/2] , ppt = oo ]; ttt = Append[ttt, Text["F(" <> ToString[ppt] <> ")", pt + {1, -0.5}, {1,0}]]; curv = gsubdiv[pol1, r, s, l, m, n]; curv = tolineseg[curv]; Print["Ready to display! "]; Show[ Graphics[{controlfig2D[{poly}], Point[pt], ttt, np, curv}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; ]; (* 3D version of subdiv2D *) subdiv3D[{poly__},r_,s_,l_,m_,lx_,mx_,ly_,my_,lz_,mz_,n_] := Block[ {cpoly = {poly}, np, wl, curv, pol1}, pol1 = prepare[cpoly]; np = newcpoly[pol1, r, s, l, m]; np = newpoly[np]; np = rfig3D[np]; wl = Last[np]; np = Drop[np,-1]; Print["Weights: ", wl]; curv = gsubdiv[pol1, r, s, l, m, n]; curv = tolineseg[curv]; Print["Ready to display! "]; Show[ Graphics3D[{controlfig3D[{poly}], np, curv}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* computes the weights of the complementary part of a closed *) (* curve, i.e. the part corresponding to t outside ]-1, 1[ *) negpoly[{poly__}] := Block[ {cpoly = {poly}, l1, i, pol, pt, w, h}, pol = {}; l1 = Length[cpoly]; Do[ pt = cpoly[[i]]; If[OddQ[l1 - i], pt = -pt]; pol = Prepend[pol, pt], {i, 1, l1} ]; pol ]; (* to display a closed rational 2D curve. *) showclo2D[{poly__},r_,s_,lx_,mx_,ly_,my_,n_] := Block[ {cpoly = {poly}, np1, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2}, l = -1; m = 1; (* maps control poly to hat space *) pol1 = prepare[cpoly]; (* computes control polygon wrt [-1, 1] *) np1 = newcpoly[pol1, r, s, l, m]; (* projects down to affine space *) npol1 = newpoly[np1]; npol1 = rfig2D[npol1]; wl1 = Last[npol1]; npol1 = Drop[npol1,-1]; Print["Weights: ", wl1]; Print["New Cont. Poly.: ", np1]; (* computes dual control polygon *) pol2 = negpoly[np1]; np2 = newpoly[pol2]; npol2 = rfig2D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = gsubdiv[pol1, r, s, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = gsubdiv[pol2, l, m, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; Show[ Graphics[{controlfig2D[{poly}], npol1, npol2, curv1, curv2}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; ]; (* to display a closed rational 2D curve. First computes control poly wrt [-1 , 1] *) (* If flagpol =!= 1, do not display control polygon *) (* If flagpol === -1, show original arc in black, and dual arc in red *) rendclo2D[{poly__},r_,s_,lx_,mx_,ly_,my_,flagpol_,n_] := Block[ {cpoly = {poly}, np1, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2, image}, l = -1; m = 1; (* maps control poly to hat space *) pol1 = prepare[cpoly]; (* computes control polygon wrt [-1, 1] *) np1 = newcpoly[pol1, r, s, l, m]; (* projects down to affine space *) npol1 = newpoly[np1]; npol1 = rfig2D[npol1]; wl1 = Last[npol1]; npol1 = Drop[npol1,-1]; Print["Weights: ", wl1]; Print["New Cont. Poly.: ", np1]; (* computes dual control polygon *) pol2 = negpoly[np1]; np2 = newpoly[pol2]; npol2 = rfig2D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[np1, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; If[flagpol === 1, image = {controlfig2D[{poly}], npol1, curv1, {RGBColor[1,0,0], curv2}}, If[flagpol === -1, image = {curv1, {RGBColor[1,0,0], curv2}}, image = {curv1, curv2} ] ]; Show[Graphics[{image}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; ]; (* to display a closed rational 2D curve. Use original control polygon wrt [0, 1] *) (* If flagpol =!= 1, do not display control polygon *) (* If flagpol === -1, show original arc in black, and dual arc in red *) rendcloz2D[{poly__},r_,s_,lx_,mx_,ly_,my_,flagpol_,n_] := Block[ {cpoly = {poly}, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2, image}, l = -1; m = 1; (* maps control poly to hat space *) pol1 = prepare[cpoly]; npol1 = rfig2D[pol1]; wl1 = Last[npol1]; Print["Weights: ", wl1]; (* computes dual control polygon *) pol2 = negpoly[pol1]; np2 = newpoly[pol2]; npol2 = rfig2D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[pol1, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; If[flagpol === 1, image = {controlfig2D[cpoly], curv1, {RGBColor[1,0,0], curv2}}, If[flagpol === -1, image = {curv1, {RGBColor[1,0,0], curv2}}, image = {curv1, curv2} ] ]; Show[Graphics[{image}], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; ]; (* 3D version of showclo2D *) showclo3D[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,n_] := Block[ {cpoly = {poly}, np1, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2}, l = -1; m = 1; pol1 = prepare[cpoly]; np1 = newcpoly[pol1, r, s, l, m]; npol1 = newpoly[np1]; npol1 = rfig3D[npol1]; wl1 = Last[npol1]; npol1 = Drop[npol1,-1]; Print["Weights: ", wl1]; Print["New Cont. Poly.: ", np1]; pol2 = negpoly[np1]; np2 = newpoly[pol2]; npol2 = rfig3D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[np1, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; Show[ Graphics3D[{controlfig3D[{poly}], npol1, curv1, curv2}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* 3D version of rendclo2D *) (* to display a closed rational 3D curve. *) (* If flagpol =!= 1, do not display control polygon *) rendclo3D[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,flagpol_,n_] := Block[ {cpoly = {poly}, np1, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2, image}, l = -1; m = 1; pol1 = prepare[cpoly]; np1 = newcpoly[pol1, r, s, l, m]; npol1 = newpoly[np1]; npol1 = rfig3D[npol1]; wl1 = Last[npol1]; npol1 = Drop[npol1,-1]; Print["Weights: ", wl1]; Print["New Cont. Poly.: ", np1]; pol2 = negpoly[np1]; np2 = newpoly[pol2]; npol2 = rfig3D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[np1, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; If[flagpol === 1, image = {controlfig3D[{poly}], npol1, curv1, {RGBColor[1,0,0], curv2}}, If[flagpol === -1, image = {curv1, {RGBColor[1,0,0], curv2}}, image = {curv1, curv2} ] ]; Show[Graphics3D[{image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (* 3D version of rendcloz2D *) (* to display a closed rational 3D curve. *) (* Use original control polygon wrt [0, 1] *) (* If flagpol =!= 1, do not display control polygon *) rendcloz3D[{poly__},r_,s_,lx_,mx_,ly_,my_,lz_,mz_,flagpol_,n_] := Block[ {cpoly = {poly}, np2, npol1, npol2, wl1, wl2, l, m, pol, curv1, curv2, pol1, pol2, image}, l = -1; m = 1; (* maps control poly to hat space *) pol1 = prepare[cpoly]; npol1 = rfig3D[pol1]; wl1 = Last[npol1]; Print["Weights: ", wl1]; (* computes dual control polygon *) pol2 = negpoly[pol1]; np2 = newpoly[pol2]; npol2 = rfig3D[np2]; wl2 = Last[npol2]; npol2 = Drop[npol2,-1]; Print["Dual Weights: ", wl2]; Print["Dual Cont. Poly.: ", pol2]; curv1 = subdiv[pol1, l, m, n]; curv1 = tolineseg[curv1]; Print["First curve segment done! "]; curv2 = subdiv[pol2, l, m, n]; curv2 = tolineseg[curv2]; Print["Second curve segment done! "]; Print["Ready to display! "]; If[flagpol === 1, image = {controlfig3D[{poly}], curv1, {RGBColor[1,0,0], curv2}}, If[flagpol === -1, image = {curv1, {RGBColor[1,0,0], curv2}}, image = {curv1, curv2} ] ]; Show[Graphics3D[{image}], PlotRange -> {{lx, mx}, {ly, my}, {lz, mz}}, Axes -> True, AxesLabel -> {"x", "y", "z"}, DisplayFunction -> $DisplayFunction]; ]; (*** parametric quartic ***) xx1[t_] := t^4 - 2 * t^3 - t^2 + 2 * t; yy1[t_] := t^3; xx2[t_] := t^3 - 2 * t^2 - t^ + 2; yy2[t_] := t * (t^3 - 2 * t^2 - t^ + 2); xx3[t_] := (t - 1) * (t^2 + t + 1); yy3[t_] := t * (t - 1) * (t^2 + t + 1); xx4[t_] := t * (t - 1) * (t + 1) * (t - 2); yy4[t_] := (t - 0.5) * (t - 1.5) * (t + 0.5); xx5[t_] := t^4 - t^3; yy5[t_] := t^4 + t; curv2D[ll_,mm_] := ParametricPlot[{xx5[t], yy5[t]}, {t, ll, mm}, DisplayFunction -> Identity]; showpic2D[ll_,mm_,lx_,mx_,ly_,my_] := Show[curv2D[ll,mm], AspectRatio -> Automatic, PlotRange -> {{lx, mx}, {ly, my}}, DisplayFunction -> $DisplayFunction]; (*** star-shape control polygons for curves with double nodes ***) conspoly[{dd__},{ww__},m_] := Block[ {poly = {}, w = {ww}, dist = {dd}, x, y, i, p, step, inc, n}, (n = m + 1; If[OddQ[m], inc = n/2 -1]; Do[ If[OddQ[n], p = (n - 1)/2; step = Mod[i * p, n]; x = N[dist[[i+1]] * Cos[(step * 2 * Pi)/n]]; y = N[dist[[i+1]] * Sin[(step * 2 * Pi)/n]] , p = n/2; If[OddQ[i], inc = inc + p, inc = Mod[inc + p + 1, n]]; x = N[dist[[i+1]] * Cos[(inc * 2 * Pi)/n]]; y = N[dist[[i+1]] * Sin[(inc * 2 * Pi)/n]] ]; poly = Append[poly, {x, y, w[[i+1]]}], {i, 0, m} ]; res = poly; res ) ];