#include /* Last edited on 2008-02-16 00:15:47 by stolfi */ #include #include #include #include #include #include #include #include #include struct vsignal_t { vsignal_frame_count_t nt; /* Number of sampling times (data frames). */ int nc; /* Number of channels per data frame. */ bool_t per; /* TRUE if the signal is periodic. */ vsignal_cont_t c; /* Continuity order, including {-1} for piecewise continuous. */ vsignal_kind_t kind; /* Reconstruction filter kind. */ double *ts; /* {ts} is NULL, or {ts[i]} is the time of sample frame {i}. */ double *vs; /* {vs[i*nch + j]} is the sample in channel {j} of data frame {i}. */ } vsignal_t; /* If {ts} is NULL, we assume {ts[i]} to be {i}. ??? Should we generalize this convention for uniform (but not unit) steps ??? The reconstruction filter actually defines two functions {f.v(z)} and {f.t(z)}, where {z} is an auxiliary parameter roughly corresponding to the frame index. We hope !!! that {f.t(z)} increases monotonically with {z}. The externally visible function {f(t)} of time is defined as {f(t) = f.w(f.t^{-1}(t))}, where {f.t^{-1}} is the functional inverse of {f.t}. For the reconstruction, the the {z}-range of {f.t} and {f.v} are divided into consecutive unit-length intervals (called /cells/). Cell number {i} stretches from {zmin[i] == i-r} to {zmax[i] == i+r)}; where {r == } .??? Note that the cell endpoints are half-integers when {cont} is odd, and integers when {cont} is even. {cont+1} consecutive sampling times and the corresponding frames to define the value of {V(z)} and {T(z)} within each cell the {t[i]} is the time currently associated with {p[i]}, then between times {t[i-1]} and {t[i]} the curve is assumed to have constant speed {h[i] == s[i]/(t[i]-t[i-1])}. Therefore the arc length is supposed to vary linearly with time between consecutive samples. The curve's position at a time {t'} is estimated from the arclength {s'} by Hermite interpolation between {0 -> (p[i-1], h[i]*d[i-1])} and {s[i+1] -> (p[i],h[i]*d[i])} or even by linear interpolation between {p[i-1]} and {p[i]}, depending on the context. Thus the position estimates are not fully consistent with the constant-speed assumtion nor with the arclengths {s[i]}. Hopefully the discrepancies will have negligible effect if the samples are sufficiently dense. The directions {d[i]} are automatically computed so as to minimize the velocity variation along the whole curve. The label function {c.r(t)} is defined by periodic linear interpolation of the values {c.r[i]}. !!! Revise !!! The vectors {c.d} and {c.s} are used internally to speed up some operations. The vector {d[i]} is the estimated (unit-length) direction of the curve at point {p[i]}, and {s[i]} is the estimated arc length between {p[i-1]} and {p[i]}, for {i IN [0..m-1]}. */ int vsignal_channel_count(vsignal_t *f) { return f->nc; } bool_t vsignal_is_periodic(vsignal_t *f) { return f->per; } vsignal_cont_t vsignal_continuity(vsignal_t *f) { return f->c; } vsignal_kind_t vsignal_kind(vsignal_t *f) { return f->kind; } bool_t vsignal_is_empty(vsignal_t *f) { return (f->nt == 0); } vsignal_frame_count_t vsignal_frame_count(vsignal_t *f) { return f->nt; } double vsignal_time_period(vsignal_t *f) { if (f->per) { assert(f->nt >= 2); return vsignal_get_sample_time(f, f->nt-1) - vsignal_get_sample_time(f, 0); } else { return 0; } } void vsignal_data_period(vsignal_t *f, double v[]) { int j; if (f->per) { assert(f->nt >= 2); for (j = 0; j < f->nc; f++) { v[j] = vsignal_get_data_sample(f, f->nt-1, j) - vsignal_get_data_sample(f, 0, j); } } else { for (j = 0; j < f->nc; f++) { v[j] = 0; } } } void vsignal_set_sample_time(vsignal_t *f, vsignal_frame_index_t i, double t) { if (f->per) { /* Reduce {i} and {t} modulo the fundamental period: */ double T = vsignal_time_period(vsignal_t *f); assert(f->nt >= 2); vsignal_frame_index_t ir = imod(i, f->nt); int64_t q = (i - ir)/f->nt; t += q*T; } else { ??? } } double vsignal_get_sample_time(vsignal_t *f, vsignal_frame_index_t i) void vsignal_set_sample(vsignal_t *f, vsignal_frame_index_t i, int j, double v) double vsignal_get_sample(vsignal_t *f, vsignal_frame_index_t i, int j) void vsignal_set_frame(vsignal_t *f, vsignal_frame_index_t i, double v[]) void vsignal_get_frame(vsignal_t *f, vsignal_frame_index_t i, double v[]) void vsignal_set_sampling_times(vsignal_t *f, int ts[]) interval_t vsignal_fundamental_domain(vsignal_t *f) { } void vsignal_eval(vsignal_t *f, double t, double v[]) { } void vsignal_eval_diff(vsignal_t *f, double t, double v[], double dv[]) { } vsignal_t vsignal_new(int nc, vsignal_frame_count_t nt, bool_t per, vsignal_kind_t k, vsignal_cont_t c) void vsignal_set_sample_counts(vsignal_t *f, int nc, vsignal_frame_count_t nt) void vsignal_make_periodic(vsignal_t *f, bool_t per) void vsignal_set_kind(vsignal_t *f, vsignal_kind_t k) void vsignal_set_cont(vsignal_t *f, vsignal_cont_t c) vsignal_t *vsignal_copy(vsignal_t *f) void vsignal_copy_contents(vsignal_t *d, vsignal_t *f) void vsignal_uniformize(vsignal_t *f, double ) void vsignal_filter(vsignal_t *d, double sigma, vsignal_t *f) void vsignal_diff(vsignal_t *d, vsignal_t *f) T New( unsigned m ) { { /* with */ ??? c = NEW(T; m := m; p := pz_r3_chain_new(m); t := pz_double_chain_new(m); r := pz_double_chain_new(m); d := pz_r3_chain_new(m); s := double_vec_new(m); a := double_vec_new(m); tPeriod := 1.0e0; rPeriod := 1.0e0 ); ??? p = c.p^; ??? t = c.t^; ??? r = c.r^; ??? d = c.d^; ??? s = c.s^; /* do */ for (i = 0; i <= m-1 ; i++){ p[i] = r3_t{0.0e0, 0.0e0, 0.0e0}; t[i] = FLOAT(i,double)/FLOAT(m,double); r[i] = t[i]; d[i] = r3_t{0.0e0, 0.0e0, 0.0e0}; s[i] = 0.0e0;; }; return c; } } /* New */ T Copy( T c ) { { /* with */ ??? m = c.m; ??? d = NEW(T; m := m; p := pz_r3_chain_new(m); t := double_vec_new(m); r := pz_double_chain_new(m); d := pz_r3_chain_new(m); s := double_vec_new(m); a := double_vec_new(m); tPeriod := c.tPeriod; rPeriod := c.rPeriod ); /* do */ d.p^ = c.p^; d.t^ = c.t^; d.r^ = c.r^; d.d^ = c.d^; d.s^ = c.s^; return d; } } /* Copy */ void ComputeArcParameters( r3_t *pa, r3_t *da, r3_t *pb, r3_t *db, LR2.T *ctr, LR2.T *hab, LR2.T *wa, LR2.T *wb ) /* Given the endpoints {pa,pb} of a spline arc, and the unit tangent vectors {da,db} at its endpoints, computes the arc parameters {ctr, hab, wa, wb} needed by {EvalArcPos} and {EvalArcVel}. */ { affirm(pa[2] == 0.0e0 ); affirm(pb[2] == 0.0e0 ); ctr = LR2.T{0.5e0*(pa[0] + pb[0]), 0.5e0*(pa[1] + pb[1])}; hab = LR2.T{0.5e0*(pb[0] - pa[0]), 0.5e0*(pb[1] - pa[1])}; if (( hab[0] == 0.0e0 ) ANDAND ( hab[1] == 0.0e0 )){ wa = LR2.T{0.0e0, 0.0e0}; wb = LR2.T{0.0e0, 0.0e0}; }else{ { /* with */ ??? u = LR2.Dir(hab); /* Directions {dra,drb} are {da,db} relative to {u} */ ??? dra = LR2.T{da[0]*u[0] + da[1]*u[1], da[1]*u[0] - da[0]*u[1]}; ??? drb = LR2.T{db[0]*u[0] + db[1]*u[1], db[1]*u[0] - db[0]*u[1]}; ??? sumx = dra[0] + drb[0]; ??? difx = dra[0] - drb[0]; ??? sumy = dra[1] + drb[1]; ??? dify = dra[1] - drb[1]; ??? ux2 = (sumx*sumx + sumy*sumy)/4.0e0; ??? uy2 = (difx*difx + dify*dify)/4.0e0; ??? mag = (2.0e0*ux2 + 3.0e0*uy2)/(ux2+uy2 + 1.0e-50); /* do */ wa = LR2.Scale(mag, LR2.T{dra[0], dra[1]}); wb = LR2.Scale(mag, LR2.T{drb[0], drb[1]}); }; } } /* ComputeArcParameters */ r3_t EvalArcPos( LR2.T *ctr, LR2.T *hab, LR2.T *wa, LR2.T *wb, double s ) /* Given the parameters {ctr, hab, wa,wb} of an arc parametrized in {[-1..+1]}, computes its position at parameter value {s}. */ { { /* with */ ??? sumx = wa[0]+wb[0]; ??? difx = wb[0]-wa[0]; ??? kx = 0.125e0 * ((((sumx - 4.0e0)*s + difx)*s - sumx + 12.0e0)*s - difx); ??? sumy = wa[1]+wb[1]; ??? dify = wb[1]-wa[1]; ??? ky = 0.125e0 * (((sumy * s + dify)*s - sumy)*s - dify); ??? xt = ctr[0] + kx*hab[0] - ky*hab[1]; ??? yt = ctr[1] + kx*hab[1] + ky*hab[0]; /* do */ return r3_t{xt, yt, 0.0e0}; } } /* EvalArcPos */ r3_t EvalArcVel( <*UNUSED*> READONLY ctr: LR2.T, LR2.T *hab, LR2.T *wa, LR2.T *wb, double s ) /* Given the parameters {ctr, hab, wa,wb} of an arc parametrized in {[-1..+1]}, computes the derivative {dp/ds} of its position {p(s)} rel to {s}, at parameter value {s}. */ { { /* with */ ??? sumx = wa[0]+wb[0]; ??? difx = wb[0]-wa[0]; ??? dkx = 0.125e0 * (((3.0e0*sumx - 12.0e0)*s + 2.0e0*difx)*s - sumx + 12.0e0); ??? sumy = wa[1]+wb[1]; ??? dify = wb[1]-wa[1]; ??? dky = 0.125e0 * (((3.0e0*sumy)*s + 2.0e0*dify)*s - sumy); ??? dxs = dkx*hab[0] - dky*hab[1]; ??? dys = dkx*hab[1] + dky*hab[0]; /* do */ return r3_t{dxs, dys, 0.0e0}; } } /* EvalArcVel */ <*UNUSED*> r3_t EvalMidArcVel( <*UNUSED*> READONLY ctr: LR2.T, LR2.T *hab, LR2.T *wa, LR2.T *wb ) /* Same as {EvalArcVel(ctr, hab, wa,wb, 0)} but a bit faster. */ { { /* with */ ??? sumx = wa[0]+wb[0]; ??? dkx = 0.125e0 * (12.0e0 - sumx); ??? sumy = wa[1]+wb[1]; ??? dky = 0.125e0 * (- sumy); ??? dxt = dkx*hab[0] - dky*hab[1]; ??? dyt = dkx*hab[1] + dky*hab[0]; /* do */ return r3_t{dxt, dyt, 0.0e0}; } } /* EvalMidArcVel */ double ComputeArcLength( <*UNUSED*> READONLY ctr: LR2.T, LR2.T *hab, LR2.T *wa, LR2.T *wb ) /* Estimates the length of the Hermite curve with the given paramteres. */ /* Estimates the length of the arc with the given Arc parameters. */ { { /* with */ /* Compute the derivative of the Hermite interpolant. */ /* Let the Hermite interpolant be {h(u) = h3*u^3 + ... + h0}; where {u} ranges from -1 to +1. Note that {dh/du(-1) = wa/2}, and {dh/du(1) = wb/2}. The polynomial {h}, relative to the {hab} frame, is then | h(u) = | | + (-1,0)*(+u^3 - 3*u + 2)/4 | + (+1,0)*(-u^3 + 3*u + 2)/4 | + wa*(u^3 - u^2 - u + 1)/8 | + wb*(u^3 + u^2 - u - 1)/8 | = | + (+ 1/8*(wa+wb) - (1/2,0))*u^3 | + (+ 1/8*(wb-wa) )*u^2 | + (- 1/8*(wa+wb) + (3/2,0))*u | + (- 1/8*(wb-wa) ) Its derivative with respect to {u} is then | dh/du(u) = | + (+ 3/8*(wa+wb) - (3/2,0)) * u^2 | + (+ 1/4*(wb-wa) ) * u | + (- 1/8*(wa+wb) + (3/2,0)) */ ??? dxdu_2 = -1.5e0 + 0.375e0*(wa[0]+wb[0]); ??? dydu_2 = + 0.375e0*(wa[1]+wb[1]); ??? dxdu_1 = 0.250e0*(wb[0] - wa[0]); ??? dydu_1 = 0.250e0*(wb[1] - wa[1]); ??? dxdu_0 = +1.5e0 - 0.125e0*(wa[0]+wb[0]); ??? dydu_0 = - 0.125e0*(wa[1]+wb[1]); /* Computes the speed squared {|dh/du|^2}, as a function of {u}. */ /* We only need the even powers: */ ??? s_4 = dxdu_2*dxdu_2 + dydu_2*dydu_2; s_2 = 2.0e0 * (dxdu_2*dxdu_0 + dydu_2*dydu_0) + (dxdu_1*dxdu_1 + dydu_1*dydu_1); ??? s_0 = dxdu_0*dxdu_0 + dydu_0*dydu_0; /* Computes the approximate square root of the polynomial. */ /* Assumes the constant term {s0} dominates, so */ /* {sqrt(U + s_0) ~ sqrt(s_0) + U/sqrt(s_0)/2} */ ??? r_0 = sqrt(s_0 + 1.0e-50); ??? r_2 = 0.5e0 * s_2/r_0; ??? r_4 = 0.5e0 * s_4/r_0; /* Integrates the speed from -1 to +1 to get the length: */ ??? length = 2.0e0*(r_4/5.0e0 + r_2/3.0e0 + r_0); /* do */ return LR2.Norm(hab)*length; } } /* ComputeArcLength */ T Diff( T c ) { { /* with */ ??? m = c.m; ??? ct = c.t^, cp = c.p^, cr = c.r^; ??? d = New(m); ??? dt = d.t^, dp = d.p^, dr = d.r^; /* do */ if (( m <= 2 )){ for (i = 0; i <= m-1 ; i++){ dp[i] = r3_t{0.0e0, 0.0e0, 0.0e0} ;}; }else{ for (io = 0; io <= m-1 ; io++){ { /* with */ ??? im = (io-1) MOD m; ??? tm = ct[im] + FLOAT((io-1) DIV m, double) * c.tPeriod; ??? to = ct[io]; ??? ip = (io+1) MOD m; ??? tp = ct[ip] + FLOAT((io+1) DIV m, double) * c.tPeriod; ??? vel = pz_geo.EstimateVelocityC(tm, cp[im], to, cp[io], tp, cp[ip]); /* do */ dp[io] = r3_t{vel[0], vel[1], 0.0e0}; affirm(dp[io][0] + 1.0e0 > dp[io][0] );; }; };; }; for (i = 0; i <= m-1 ; i++){ dt[i] = ct[i]; dr[i] = cr[i] ;}; d.tPeriod = c.tPeriod; d.rPeriod = c.rPeriod; RecomputeDirections(d); RecomputeLengths(d); return d; } } /* Diff */ /* ----------------------------------- Old Diff T Diff( T c ) VAR ctr, hab, wa, wb: LR2.T; { { /* with */ ??? ct = c.t^, cp = c.p^, cr = c.r^, cd = c.d^; ??? d = New(c.m); ??? dt = d.t^, dp = d.p^, dr = d.r^; /* do */ d.tPeriod = c.tPeriod; d.rPeriod = c.rPeriod; for (i = 0; i < (cp.nel ) ; i++){ { /* with */ ??? ia = (i - 1) MOD c.m; ??? qa = ((double)ia DIV c.m); ??? pa = cp[ia]; ??? da = cd[ia]; ??? ta = ct[ia] + c.tPeriod * qa; ??? ra = cr[ia] + c.rPeriod * qa; ??? ib = i; ??? pb = cp[ib]; ??? db = cd[ib]; ??? tb = ct[ib]; ??? rb = cr[ib]; /* do */ ComputeArcParameters(pa, da, pb, db, ctr, hab, wa, wb); dt[i] = 0.5e0 * (ta + tb); dp[i] = EvalMidArcVel(ctr, hab, wa, wb); dr[i] = 0.5e0 * (ra + rb);; }; }; RecomputeDirections(d); RecomputeLengths(d); return d; } } /* Diff */ -----------------------------------------------------------*/ void SetSamples( T c, pz_r3_chain_t *p ) { affirm((p.nel) == c.m ); c.p^ = p; RecomputeDirections(c); RecomputeLengths(c); } /* SetSamples */ void SetTimes( T c, pz_double_chain_t *t, double tPeriod ) { affirm((t.nel) == c.m ); c.t^ = t; c.tPeriod = tPeriod; } /* SetTimes */ void SetLabels( T c, pz_double_chain_t *r, double rPeriod ) { affirm((r.nel) == c.m ); c.r^ = r; c.rPeriod = rPeriod; } /* SetLabels */ void RecomputeDirections( T c ) /* Recomputes the directions {c.d[i]} from the samples {c.p[i]}. Does not use the times {c.t}, the lengths {c.s}, or the labels {c.r}. */ VAR im, ip: unsigned; { { /* with */ ??? m = c.m; ??? p = c.p^; ??? d = c.d^; /* do */ if (( m <= 2 )){ for (i = 0; i <= m-1 ; i++){ d[i] = r3_t{0.0e0, 0.0e0, 0.0e0} ;}; }else{ im = m-1; for (io = 0; io <= m-1 ; io++){ ip = (io+1) MOD m; { /* with */ ??? tm = r3_dist(p[im], p[io]) + 1.0e-50; ??? tp = r3_dist(p[io], p[ip]) + 1.0e-50; ??? vel = pz_geo.EstimateVelocityC(-tm, p[im], 0.0e0, p[io], +tp, p[ip]); /* do */ if (( vel[0] == 0.0e0 ) ANDAND ( vel[1] == 0.0e0 )){ d[io] = r3_t{0.0e0, 0.0e0, 0.0e0} }else{ d[io] = r3_dir(r3_t{vel[0], vel[1], 0.0e0});; }; affirm(d[io][0] + 1.0e0 > d[io][0] );; }; im = io;; };; }; } } /* RecomputeDirections */ void RecomputeLengths( T c ) /* Recomputes the estimated arclength {s[i]} between points {p[i-1]} and {p[i]}. Uses the points {c.p} and the estimated direction {c.s}; does not use the times {c.t} or the labels {c.r}. */ VAR ia: unsigned; VAR ctr, hab, wa, wb: LR2.T; { { /* with */ ??? m = c.m; ??? p = c.p^; ??? d = c.d^; ??? s = c.s^; /* do */ if (( m <= 1 )){ for (i = 0; i <= m-1 ; i++){ s[0] = 0.0e0 ;}; }else{ ia = m-1; for (ib = 0; ib <= m-1 ; ib++){ ComputeArcParameters(p[ia], d[ia], p[ib], d[ib], ctr, hab, wa, wb); { /* with */ ??? si = ComputeArcLength(ctr, hab, wa, wb); /* do */ s[ib] = si; ia = ib; }; }; }; } } /* RecomputeLengths */ void Uniformize( T c ) VAR totLen: double; { { /* with */ ??? m = c.m; ??? s = c.s^; ??? t = c.t^; /* do */ totLen = 0.0e0; for (i = 0; i <= m-1 ; i++){ t[i] = totLen; totLen = totLen + s[i];; }; c.tPeriod = totLen /* No need to recompute lengths since they don't depend on times. */; } } /* Uniformize */ void SetToFiltered( T c, T d, double sigma ) { { /* with */ ??? ct = c.t^; ??? m = c.m; ??? T = d.tPeriod; ??? arg = c.a^; /* do */ /* Temporary hack: filters the *linear* interpolation of the points {p[i]}. */ pz_smooth_polygonal(sigma, d.p^, d.t^, d.tPeriod, 0.0e0, c.p^); /* {pz_smooth} will resample at {m} uniform times in {[0..T]}: */ for (i = 0; i <= m-1 ; i++){ ct[i] = T * ((double)i)/((double)m); }; c.tPeriod = d.tPeriod; /* Recompute the velocities and lengths: */ RecomputeDirections(c); RecomputeLengths(c); /* Compute the labels of the sampled points: */ pz_double_chain_args_from_values(d.t^, d.tPeriod, ct, arg); pz_double_chain_values_from_args(d.r^, d.rPeriod, arg, c.r^); c.rPeriod = d.rPeriod;; }; } /* SetToFiltered */ void Eval( T c, double t, r3_t *p, r3_t *v, double *r ) VAR ia, ib: unsigned; k: int; VAR ctr, hab, wa, wb: LR2.T; /* Computed for arc {c.p[ih-1]} .. {c.p[ih]} */ { { /* with */ ??? cm = c.m; ??? cp = c.p^; ??? cd = c.d^; ??? ct = c.t^; ??? cr = c.r^; ??? tLo = ct[0]; ??? tHi = tLo + c.tPeriod; /* do */ affirm((ct.nel) == cm ); /* Find consecutive times {ct[ia]}, {ct[ib]} that bracket {t} */ ia = 0; ib = 1; /* Invariant: {ia IN [0..cm-1]}, {ib == ia+1} */ pz_proc.ReduceToPeriod(t, tLo, tHi, k); /* Should do binary search ! */ while ( t < ct[ia] ){ ib = ia; ia-- ;}; while ( ib < cm ) ANDAND ( t > ct[ib] ){ ia = ib; ib++ ;}; { /* with */ ??? ta = ct[ia]; ??? tb = ct[ib MOD cm] + c.tPeriod * ((double)ib DIV cm); ??? ra = cr[ia]; ??? rb = cr[ib MOD cm] + c.rPeriod * ((double)ib DIV cm); /* do */ affirm(ta <= t ) ANDAND ( t <= tb ); ComputeArcParameters( cp[ia], cd[ia], cp[ib MOD cm], cd[ib MOD cm], ctr, hab, wa, wb ); { /* with */ ??? sh = (t - ta)/(tb - ta); ??? s = 2.0e0*sh - 1.0e0; /* do */ p = EvalArcPos(ctr, hab, wa, wb, s); v = EvalArcVel(ctr, hab, wa, wb, s); r = ra + sh * (rb - ra) + FLOAT(k,double) * c.rPeriod; };; }; ; } } /* Eval */ void Sample( T c, pz_double_chain_t *t, pz_r3_chain_t *p, pz_r3_chain_t *d, pz_double_chain_t *r ) VAR ia, ib, ih: unsigned; k: int; tau: double; VAR ctr, hab, wa, wb: LR2.T; /* Computed for arc {c.p[ih-1]} .. {c.p[ih]} */ { { /* with */ ??? cm = c.m; ??? cp = c.p^; ??? cd = c.d^; ??? ct = c.t^; ??? cr = c.r^; ??? n = (p.nel); ??? tLo = ct[0]; ??? tHi = tLo + c.tPeriod; /* do */ affirm((ct.nel) == cm ); affirm((t.nel) >= n ); ia = 0; ib = 1; ih = (unsigned.nel - 1); for (j = 0; j <= n-1 ; j++){ /* Invariant: {ia IN [0..cm-1]}, {ib == ia+1} */ tau = t[j]; pz_proc.ReduceToPeriod(tau, tLo, tHi, k); while ( tau < ct[ia] ){ ib = ia; ia-- ;}; while ( ib < cm ) ANDAND ( tau > ct[ib] ){ ia = ib; ib++ ;}; { /* with */ ??? ta = ct[ia]; ??? tb = ct[ib MOD cm] + c.tPeriod * ((double)ib DIV cm); ??? ra = cr[ia]; ??? rb = cr[ib MOD cm] + c.rPeriod * ((double)ib DIV cm); /* do */ affirm(ta <= tau ) ANDAND ( tau <= tb ); if (( ib != ih )){ ComputeArcParameters( cp[ia], cd[ia], cp[ib MOD cm], cd[ib MOD cm], ctr, hab, wa, wb ); ih = ib; }; { /* with */ ??? sh = (tau - ta)/(tb - ta); ??? s = 2.0e0*sh - 1.0e0; /* do */ p[j] = EvalArcPos(ctr, hab, wa, wb, s); d[j] = EvalArcVel(ctr, hab, wa, wb, s); r[j] = ra + sh * (rb - ra) + FLOAT(k,double) * c.rPeriod; };; }; ; };; } } /* Sample */ { /* Copyright © 2001 Universidade Estadual de Campinas (UNICAMP). Authors: Helena C. G. Leitão and Jorge Stolfi. This file can be freely distributed, used, and modified, provided that this copyright and authorship notice is preserved, and that any modified versions are clearly marked as such. This software has NO WARRANTY of correctness or applicability for any purpose. Neither the authors nor their employers chall be held responsible for any losses or damages that may result from its use. */