/* Copyright (C) 2012-2020, Dirk Krause SPDX-License-Identifier: BSD-3-Clause */ /* WARNING: This file was generated by the dkct program (see http://dktools.sourceforge.net/ for details). Changes you make here will be lost if dkct is run again! You should modify the original source and run dkct on it. Original source: dk3xsp.ctr */ /** @file dk3xsp.c The dk3xsp module. */ #line 285 "dk3xsp.ctr" #include "dk3all.h" #include "dk3fig.h" #include "dk3xsp.h" #line 291 "dk3xsp.ctr" /* ************************************************************************ */ /* * * */ /* * Basic X-spline functions. * */ /* * * */ /* ************************************************************************ */ #line 298 "dk3xsp.ctr" /** Initialize spline point (set all components to 0.0). @param ppt Point to initialize. */ static void dk3xsp_init_spline_point(dk3_fig_spline_point_t *ppt) { ppt->x = 0.0; ppt->y = 0.0; ppt->s = 0.0; } /** Copy spline point or initialize a destination point. @param pd Destination pointer. @param ps Source pointer, may be NULL. */ static void dk3xsp_copy_spline_point( dk3_fig_spline_point_t *pd, dk3_fig_spline_point_t const *ps ) { if(pd) { if(ps) { pd->x = ps->x; pd->y = ps->y; pd->s = ps->s; } else { dk3xsp_init_spline_point(pd); } } } void dk3xsp_reset(dk3_xspline_segment_t *pseg) { if(pseg) { dk3mem_res((void *)pseg, sizeof(dk3_xspline_segment_t)); dk3xsp_init_spline_point(&(pseg->a)); dk3xsp_init_spline_point(&(pseg->b)); dk3xsp_init_spline_point(&(pseg->c)); dk3xsp_init_spline_point(&(pseg->d)); pseg->pa = pseg->pb = pseg->pc = pseg->pd = 2.0; pseg->qa = pseg->qb = pseg->qc = pseg->qd = 0.0; pseg->x = pseg->y = pseg->dxdt = pseg->dydt = 0.0; pseg->duadt = pseg->dubdt = pseg->ducdt = pseg->duddt = 0.0; pseg->cb = 0; } } void dk3xsp_set_cb(dk3_xspline_segment_t *pseg, int v) { if(pseg) { pseg->cb = ((v) ? 1 : 0); } } /** The f(u) function. @param u The u argument. @param p The p parameter. @return Function result. */ static double dk3xsp_f(double u, double p) { double back; back = u * u * u *(10.0 - p + u * (2.0 * p - 15.0 + u * (6.0 - p))); #line 378 "dk3xsp.ctr" return back; } /** The f(u) function derived. @param u The u argument. @param p The p parameter. @return Function result. */ static double dk3xsp_dfdu(double u, double p) { double back; back = u * u * (30.0 - 3.0*p + u * (8.0*p - 60.0 + u * (30.0 - 5.0*p))); #line 395 "dk3xsp.ctr" return back; } /** The g(u) function. @param u Argument u. @param p Parameter p. @param q Parameter q. @return Function result. */ static double dk3xsp_g(double u, double p, double q) { double back; back = u*(q+u*(2.0*q+u*(10.0-12.0*q-p+u*(2.0*p+14.0*q-15+u*(6.0-5.0*q-p))))); #line 413 "dk3xsp.ctr" return back; } /** The g(u) function derived. @param u Argument u. @param p Parameter p. @param q Parameter q. @return Function result. */ static double dk3xsp_dgdu(double u, double p, double q) { double back; back = q+u*(4.0*q+u*(30.0-3.0*p-36.0*q+u*(56.0*q+8.0*p-60.0+u*(30.0-5.0*p-25.0*q)))); #line 432 "dk3xsp.ctr" return back; } /** The h(u) function. @param u Argument u. @param q Parameter q. @return Function result. */ static double dk3xsp_h(double u, double q) { double back; back = u*(q+u*(2.0*q-u*u*(q*(u+2.0)))); #line 449 "dk3xsp.ctr" return back; } /** The h(u) function derived. @param u Argument u. @param q Parameter q. @return Function result. */ static double dk3xsp_dhdu(double u, double q) { double back; back = q+u*(4.0*q-u*u*(q*(8.0+5.0*u))); #line 466 "dk3xsp.ctr" return back; } double dk3xsp_get_x(dk3_xspline_segment_t const *pseg) { double back = 0.0; if(pseg) { back = pseg->x; } return back; } double dk3xsp_get_y(dk3_xspline_segment_t const *pseg) { double back = 0.0; if(pseg) { back = pseg->y; } return back; } double dk3xsp_get_dxdt(dk3_xspline_segment_t const *pseg) { double back = 0.0; if(pseg) { back = pseg->dxdt; } return back; } double dk3xsp_get_dydt(dk3_xspline_segment_t const *pseg) { double back = 0.0; if(pseg) { back = pseg->dydt; } return back; } void dk3xsp_set( dk3_xspline_segment_t *pseg, dk3_fig_spline_point_t const *a, dk3_fig_spline_point_t const *b, dk3_fig_spline_point_t const *c, dk3_fig_spline_point_t const *d ) { #line 529 "dk3xsp.ctr" if(pseg) { pseg->ha = ((a) ? 1 : 0); pseg->hd = ((d) ? 1 : 0); dk3xsp_copy_spline_point(&(pseg->a), a); dk3xsp_copy_spline_point(&(pseg->b), b); dk3xsp_copy_spline_point(&(pseg->c), c); dk3xsp_copy_spline_point(&(pseg->d), d); pseg->pa = pseg->pb = pseg->pc = pseg->pd = 0.0; pseg->qa = pseg->qb = pseg->qc = pseg->qd = 0.0; pseg->duadt = pseg->dubdt = pseg->ducdt = pseg->duddt = 0.0; if((pseg->b).s < 0.0) { if(pseg->cb) { pseg->qa = pseg->qc = 0.0 - ((pseg->b).s); } else { pseg->qa = pseg->qc = -0.5 * ((pseg->b).s); } pseg->duadt = -1.0; pseg->ducdt = 1.0; } else { pseg->pa = pseg->pc = 2.0 * (1.0 + (pseg->b).s) * (1.0 + (pseg->b).s); pseg->ducdt = 1.0 / (1.0 + (pseg->b).s); pseg->duadt = -1.0 * pseg->ducdt; } if((pseg->c).s < 0.0) { if(pseg->cb) { pseg->qb = pseg->qd = 0.0 - ((pseg->c).s); } else { pseg->qb = pseg->qd = -0.5 * ((pseg->c).s); } pseg->dubdt = -1.0; pseg->duddt = 1.0; } else { pseg->pb = pseg->pd = 2.0 * (1.0 + (pseg->c).s) * (1.0 + (pseg->c).s); pseg->duddt = 1.0 / (1.0 + (pseg->c).s); pseg->dubdt = -1.0 * pseg->duddt; } #line 564 "dk3xsp.ctr" #line 565 "dk3xsp.ctr" #line 566 "dk3xsp.ctr" #line 567 "dk3xsp.ctr" } #line 568 "dk3xsp.ctr" } int dk3xsp_calculate( dk3_xspline_segment_t *pseg, double t, int fder ) { double ua = 0.0; /* Argument u for blending function from point A. */ double ub = 0.0; /* = B. */ double uc = 0.0; /* = C. */ double ud = 0.0; /* = D. */ double fa = 0.0; /* Result from blending function from point A. */ double fb = 0.0; /* = B. */ double fc = 0.0; /* = C. */ double fd = 0.0; /* = D. */ double zx = 0.0; /* Counter in fraction to calculate x. */ double zy = 0.0; /* Counter in fraction to calculate y. */ double n = 0.0; /* Denominator in both fractions. */ double dfadt = 0.0; /* Derivative result for point A blending function. */ double dfbdt = 0.0; /* = B. */ double dfcdt = 0.0; /* = C. */ double dfddt = 0.0; /* = D. */ double nder = 0.0; /* Derivative of fraction denominator. */ double zxder = 0.0; /* Derivative of x fraction counter. */ double zyder = 0.0; /* Derivative of y fraction counter. */ double nsq = 0.0; /* Denominator square value. */ int ec = 0; /* Error code for mathematical operations. */ int back = 0; #line 601 "dk3xsp.ctr" if(pseg) { back = 1; /* Initialize variables. */ #if VERSION_BEFORE_20150821 ua = ub = uc = ud = fa = fb = fc = fd = zx = zy = n = 0.0; #endif /* Set ua, ub, uc, ud. */ if((pseg->b).s < 0.0) { ua = -1.0 * t; uc = t; } else { ua = ((pseg->b).s - t) / (1.0 + (pseg->b).s); uc = ((pseg->b).s + t) / (1.0 + (pseg->b).s); } if((pseg->c).s < 0.0) { ub = 1.0 - t; ud = t - 1.0; } else { ub = (1.0 + (pseg->c).s - t) / (1.0 + (pseg->c).s); ud = (t + (pseg->c).s - 1.0) / (1.0 + (pseg->c).s); } #line 627 "dk3xsp.ctr" #line 628 "dk3xsp.ctr" #line 629 "dk3xsp.ctr" #line 630 "dk3xsp.ctr" /* Calculate fa, fb, fc, fd. */ if((pseg->b).s < 0.0) { fa = dk3xsp_h(ua, pseg->qa); fc = dk3xsp_g(uc, pseg->pc, pseg->qc); } else { if(t < (pseg->b).s) { fa = dk3xsp_f(ua, pseg->pa); } fc = dk3xsp_f(uc, pseg->pc); } if((pseg->c).s < 0.0) { fb = dk3xsp_g(ub, pseg->pb, pseg->qb); fd = dk3xsp_h(ud, pseg->qd); } else { fb = dk3xsp_f(ub, pseg->pb); if(t > (1.0 - (pseg->c).s)) { fd = dk3xsp_f(ud, pseg->pd); } } #line 652 "dk3xsp.ctr" #line 653 "dk3xsp.ctr" #line 654 "dk3xsp.ctr" #line 655 "dk3xsp.ctr" /* Calculate counters and denominator. */ zx = dk3ma_d_add_ok( dk3ma_d_add_ok( ((pseg->ha) ? (fa * (pseg->a).x) : 0.0), (fb * (pseg->b).x), &ec ), dk3ma_d_add_ok( (fc * (pseg->c).x), ((pseg->hd) ? (fd * (pseg->d).x) : 0.0), &ec ), &ec ); zy = dk3ma_d_add_ok( dk3ma_d_add_ok( ((pseg->ha) ? (fa * (pseg->a).y) : 0.0), (fb * (pseg->b).y), &ec ), dk3ma_d_add_ok( (fc * (pseg->c).y), ((pseg->hd) ? (fd * (pseg->d).y) : 0.0), &ec ), &ec ); n = dk3ma_d_add_ok( dk3ma_d_add_ok(((pseg->ha) ? fa : 0.0), fb, &ec), dk3ma_d_add_ok(fc, ((pseg->hd) ? fd : 0.0), &ec), &ec ); #line 690 "dk3xsp.ctr" #line 691 "dk3xsp.ctr" #line 692 "dk3xsp.ctr" /* Calculate x and y value. */ pseg->x = dk3ma_d_div_ok(zx, n, &ec); pseg->y = dk3ma_d_div_ok(zy, n, &ec); #line 698 "dk3xsp.ctr" #line 699 "dk3xsp.ctr" /* Calculate derivatives if required. */ if(fder) { /* Initialize variables. */ dfadt = dfbdt = dfcdt = dfddt = nder = zxder = zyder = 0.0; /* Calculate dfadt, dfbdt, dfcdt, dfddt. */ if((pseg->b).s < 0.0) { dfadt = dk3ma_d_mul_ok(dk3xsp_dhdu(ua, pseg->qa), pseg->duadt, &ec); dfcdt = dk3ma_d_mul_ok( dk3xsp_dgdu(uc, pseg->pc, pseg->qc), pseg->ducdt, &ec ); } else { if(t < (pseg->b).s) { dfadt = dk3ma_d_mul_ok(dk3xsp_dfdu(ua, pseg->pa), pseg->duadt, &ec); } dfcdt = dk3ma_d_mul_ok(dk3xsp_dfdu(uc, pseg->pc), pseg->ducdt, &ec); } if((pseg->c).s < 0.0) { dfbdt = dk3ma_d_mul_ok( dk3xsp_dgdu(ub, pseg->pb, pseg->qb), pseg->dubdt, &ec ); dfddt = dk3ma_d_mul_ok(dk3xsp_dhdu(ud, pseg->qd), pseg->duddt, &ec); } else { dfbdt = dk3ma_d_mul_ok(dk3xsp_dfdu(ub, pseg->pb), pseg->dubdt, &ec); if(t > (1.0 - (pseg->c).s)) { dfddt = dk3ma_d_mul_ok(dk3xsp_dfdu(ud, pseg->pd), pseg->duddt, &ec); } } #line 737 "dk3xsp.ctr" #line 738 "dk3xsp.ctr" #line 739 "dk3xsp.ctr" #line 740 "dk3xsp.ctr" /* Calculate counters and denominator. */ zxder = dk3ma_d_add_ok( dk3ma_d_add_ok( ((pseg->ha) ? dk3ma_d_mul_ok(dfadt, (pseg->a).x, &ec) : 0.0), dk3ma_d_mul_ok(dfbdt, (pseg->b).x, &ec), &ec ), dk3ma_d_add_ok( dk3ma_d_mul_ok(dfcdt, (pseg->c).x, &ec), ((pseg->hd) ? dk3ma_d_mul_ok(dfddt, (pseg->d).x, &ec) : 0.0), &ec ), &ec ); zyder = dk3ma_d_add_ok( dk3ma_d_add_ok( ((pseg->ha) ? dk3ma_d_mul_ok(dfadt, (pseg->a).y, &ec) : 0.0), dk3ma_d_mul_ok(dfbdt, (pseg->b).y, &ec), &ec ), dk3ma_d_add_ok( dk3ma_d_mul_ok(dfcdt, (pseg->c).y, &ec), ((pseg->hd) ? dk3ma_d_mul_ok(dfddt, (pseg->d).y, &ec) : 0.0), &ec ), &ec ); nder = dk3ma_d_add_ok( dk3ma_d_add_ok(((pseg->ha) ? dfadt : 0.0), dfbdt, &ec), dk3ma_d_add_ok(dfcdt, ((pseg->hd) ? dfddt : 0.0), &ec), &ec ); nsq = n * n; #line 776 "dk3xsp.ctr" #line 777 "dk3xsp.ctr" #line 778 "dk3xsp.ctr" #line 779 "dk3xsp.ctr" pseg->dxdt = dk3ma_d_div_ok( dk3ma_d_sub_ok( dk3ma_d_mul_ok(n, zxder, &ec), dk3ma_d_mul_ok(zx, nder, &ec), &ec ), nsq, &ec ); pseg->dydt = dk3ma_d_div_ok( dk3ma_d_sub_ok( dk3ma_d_mul_ok(n, zyder, &ec), dk3ma_d_mul_ok(zy, nder, &ec), &ec ), nsq, &ec ); #line 798 "dk3xsp.ctr" #line 799 "dk3xsp.ctr" } if(ec) { back = 0; } } #line 804 "dk3xsp.ctr" return back; } /* ************************************************************************ */ /* * * */ /* * Length calculations for single segments. * */ /* * * */ /* ************************************************************************ */ #line 812 "dk3xsp.ctr" int dk3xsp_is_null(double s) { int back = 0; if(fabs(s) < 1.0e-6) { back = 1; } return back; } /** Convert long to double. @param l Long value to convert. @return Conversion result. */ static double dk3xsp_l_to_d(long l) { return ((double)l); } double dk3xsp_segment_length( dk3_xspline_segment_t *pseg, double prec, int *ec ) { double back = -1.0; double lastpassres; /* Result of previous pass. */ double passres; /* Result of current pass. */ double deltax; /* Delta x between points. */ double deltay; /* Delta y between points. */ double lastx; /* Previous x. */ double lasty; /* Previous y. */ double x = 0.0; /* Current x. */ double y = 0.0; /* Current y. */ double t; /* Current t. */ double dsegs; /* Double value of segs. */ long segs = 64L; /* Number of segments. */ long i; /* Current segment. */ int mec = 0; /* Mathematical error code. */ int mc = 1; /* Flag: Must continue. */ int cc = 1; /* Flag: Can continue. */ #line 862 "dk3xsp.ctr" if((pseg) && (prec > 0.0)) { if(dk3xsp_is_null((pseg->b).s) && dk3xsp_is_null((pseg->c).s)) { /* Simple line. */ deltax = dk3ma_d_sub_ok((pseg->c).x, (pseg->b).x, &mec); deltay = dk3ma_d_sub_ok((pseg->c).y, (pseg->b).y, &mec); back = sqrt( dk3ma_d_add_ok( dk3ma_d_mul_ok(deltax, deltax, &mec), dk3ma_d_mul_ok(deltay, deltay, &mec), &mec ) ); #line 876 "dk3xsp.ctr" } else { #line 877 "dk3xsp.ctr" /* Spline calculation really necessary. */ lastpassres = 0.0 - 2.0 * prec; if(lastpassres > -1.0) { lastpassres = -1.0; } while((cc) && (mc) && (!(mec))) { #line 883 "dk3xsp.ctr" passres = 0.0; if(dk3xsp_calculate(pseg, 0.0, 0)) { lastx = dk3xsp_get_x(pseg); lasty = dk3xsp_get_y(pseg); dsegs = dk3xsp_l_to_d(segs); for(i = 1L; i <= segs; i++) { if(i == segs) { t = 1.0; } else { t = dk3xsp_l_to_d(i) / dsegs; } if(dk3xsp_calculate(pseg, t, 0)) { x = dk3xsp_get_x(pseg); y = dk3xsp_get_y(pseg); } else { cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } if(cc) { deltax = dk3ma_d_sub_ok(x, lastx, &mec); deltay = dk3ma_d_sub_ok(y, lasty, &mec); passres = dk3ma_d_add_ok( passres, sqrt( dk3ma_d_add_ok( dk3ma_d_mul_ok(deltax, deltax, &mec), dk3ma_d_mul_ok(deltay, deltay, &mec), &mec ) ), &mec ); lastx = x; lasty = y; } } #line 918 "dk3xsp.ctr" if(fabs(passres - lastpassres) < prec) { #line 919 "dk3xsp.ctr" mc = 0; /* Finished successfully. */ back = passres; #line 922 "dk3xsp.ctr" } else { #line 923 "dk3xsp.ctr" lastpassres = passres; segs = segs * 2L; if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) { #line 926 "dk3xsp.ctr" cc = 0; if(ec) { *ec = DK3_ERROR_ITERATION; } } } } else { cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } } } if(mec) { back = -1.0; if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 940 "dk3xsp.ctr" return back; } double dk3xsp_segment_partial_length( dk3_xspline_segment_t *pseg, double tstart, int dir, double prec, int *ec ) { double back = -1.0; double lastpassres; /* Result of previous pass. */ double passres; /* Result of current pass. */ double deltax; /* Delta x between points. */ double deltay; /* Delta y between points. */ double lastx; /* Previous x. */ double lasty; /* Previous y. */ double x = 0.0; /* Current x. */ double y = 0.0; /* Current y. */ double t; /* Current t. */ double dsegs; /* Double value of segs. */ double tdiff; /* 1.0 - tstart. */ long segs = 64L; /* Number of segments. */ long i; /* Current segment. */ int mec = 0; /* Mathematical error code. */ int mc = 1; /* Flag: Must continue. */ int cc = 1; /* Flag: Can continue. */ #line 972 "dk3xsp.ctr" if((pseg) && (prec > 0.0) && (tstart >= 0.0) && (tstart<= 1.0)) { #line 973 "dk3xsp.ctr" if(dk3xsp_is_null((pseg->b).s) && dk3xsp_is_null((pseg->c).s)) { /* Simple line. */ #line 978 "dk3xsp.ctr" if(dk3xsp_calculate(pseg, tstart, 0)) { x = dk3xsp_get_x(pseg); y = dk3xsp_get_y(pseg); } else { #if VERSION_BEFORE_20150821 cc = 0; #endif mec = DK3_ERROR_MATH_OVERFLOW; } if(dir) { deltax = dk3ma_d_sub_ok((pseg->b).x, x, &mec); deltay = dk3ma_d_sub_ok((pseg->b).y, y, &mec); } else { deltax = dk3ma_d_sub_ok((pseg->c).x, x, &mec); deltay = dk3ma_d_sub_ok((pseg->c).y, y, &mec); } back = sqrt( dk3ma_d_add_ok( dk3ma_d_mul_ok(deltax, deltax, &mec), dk3ma_d_mul_ok(deltay, deltay, &mec), &mec ) ); } else { /* Spline calculations really necessary. */ #line 1006 "dk3xsp.ctr" if(dir) { #line 1007 "dk3xsp.ctr" lastpassres = 0.0 - 2.0 * prec; if(lastpassres > -1.0) { lastpassres = -1.0; } while((cc) && (mc) && (!(mec))) { passres = 0.0; dsegs = dk3xsp_l_to_d(segs); if(dk3xsp_calculate(pseg, 0.0, 0)) { lastx = dk3xsp_get_x(pseg); lasty = dk3xsp_get_y(pseg); for(i = 1L; i <= segs; i++) { if(i == segs) { t = tstart; } else { t = (tstart * dk3xsp_l_to_d(i)) / dsegs; } #line 1021 "dk3xsp.ctr" if(dk3xsp_calculate(pseg, t, 0)) { #line 1023 "dk3xsp.ctr" #line 1024 "dk3xsp.ctr" x = dk3xsp_get_x(pseg); #line 1025 "dk3xsp.ctr" y = dk3xsp_get_y(pseg); #line 1026 "dk3xsp.ctr" deltax = dk3ma_d_sub_ok(x, lastx, &mec); #line 1027 "dk3xsp.ctr" deltay = dk3ma_d_sub_ok(y, lasty, &mec); #line 1028 "dk3xsp.ctr" #line 1029 "dk3xsp.ctr" passres = dk3ma_d_add_ok( passres, sqrt( dk3ma_d_add_ok( dk3ma_d_mul_ok(deltax, deltax, &mec), dk3ma_d_mul_ok(deltay, deltay, &mec), &mec ) ), &mec ); #line 1040 "dk3xsp.ctr" lastx = x; lasty = y; } else { cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } } if(fabs(passres - lastpassres) < prec) { mc = 0; back = passres; } else { lastpassres = passres; segs = segs * 2L; if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) { cc = 0; if(ec) { *ec = DK3_ERROR_ITERATION; } } } } else { cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } } } else { #line 1063 "dk3xsp.ctr" tdiff = 1.0 - tstart; lastpassres = 0.0 - 2.0 * prec; if(lastpassres > -1.0) { lastpassres = -1.0; } while((cc) && (mc) && (!(mec))) { passres = 0.0; dsegs = dk3xsp_l_to_d(segs); if(dk3xsp_calculate(pseg, tstart, 0)) { lastx = dk3xsp_get_x(pseg); lasty = dk3xsp_get_y(pseg); for(i = 1L; i <= segs; i++) { if(i == segs) { #line 1074 "dk3xsp.ctr" t = 1.0; } else { t = tstart + ((tdiff * dk3xsp_l_to_d(i)) / dsegs); } #line 1079 "dk3xsp.ctr" if(dk3xsp_calculate(pseg, t, 0)) { x = dk3xsp_get_x(pseg); y = dk3xsp_get_y(pseg); } else { #line 1083 "dk3xsp.ctr" cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } if(cc) { #line 1088 "dk3xsp.ctr" #line 1089 "dk3xsp.ctr" #line 1090 "dk3xsp.ctr" #line 1091 "dk3xsp.ctr" deltax = dk3ma_d_sub_ok(x, lastx, &mec); deltay = dk3ma_d_sub_ok(y, lasty, &mec); #line 1094 "dk3xsp.ctr" #line 1095 "dk3xsp.ctr" #line 1096 "dk3xsp.ctr" passres = dk3ma_d_add_ok( passres, sqrt( dk3ma_d_add_ok( dk3ma_d_mul_ok(deltax, deltax, &mec), dk3ma_d_mul_ok(deltay, deltay, &mec), &mec ) ), &mec ); #line 1107 "dk3xsp.ctr" lastx = x; lasty = y; } } } else { cc = 0; mec = DK3_ERROR_MATH_OVERFLOW; } if(fabs(passres - lastpassres) < prec) { mc = 0; back = passres; } else { lastpassres = passres; segs = segs * 2L; if(segs > DK3XSP_MAX_ITERATION_SEGMENTS) { cc = 0; if(ec) { *ec = DK3_ERROR_ITERATION; } } } } } } if(mec) { back = -1.0; if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1133 "dk3xsp.ctr" return back; } /* ************************************************************************ */ /* * * */ /* * Length calculations for entire splines. * */ /* * * */ /* ************************************************************************ */ #line 1141 "dk3xsp.ctr" double dk3xsp_open_spline_length( dk3_fig_spline_point_t const *points, size_t np, double prec, double *fslc, int cb, int *ec ) { dk3_xspline_segment_t seg; double val; double back = 0.0; size_t i; int mec = 0; #line 1160 "dk3xsp.ctr" if((points) && (1 < np) && (prec > 0.0)) { dk3xsp_reset(&seg); if(cb) { dk3xsp_set_cb(&seg, 1); } for(i = 0; ((i < (np - 1)) && (back > -1.0)); i++) { #line 1165 "dk3xsp.ctr" val = -1.0; if(fslc) { val = fslc[i]; } if(val < -0.5) { dk3xsp_set( &seg, ((i > 1) ? (&(points[i-1])) : NULL), &(points[i]), &(points[i+1]), ((i < (np - 2)) ? &(points[i+2]) : NULL) ); val = dk3xsp_segment_length(&seg, (prec / (double)(np - 1)), ec); if(val >= 0.0) { if(fslc) { fslc[i] = val; } } } if(val < 0.0) { back = -1.0; i = np - 1; } else { back = dk3ma_d_add_ok(back, val, &mec); if(mec) { back = -1.0; i = np - 1; if(ec) { *ec = mec; } } } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1197 "dk3xsp.ctr" return back; } /** Find Index of segment start point. @param tstart Parameter t. @param ec Pointer to error code variable, required. @return Start segment index. */ static size_t dk3xsp_find_segment_start(double tstart, int *ec) { #if VERSION_BEFORE_20140809 unsigned long ul; size_t back = 0; int mec = 0; #line 1216 "dk3xsp.ctr" ul = dk3ma_d_to_ul_ok(floor(tstart), &mec); if(mec) { *ec = mec; } else { back = (size_t)ul; if((unsigned long)back != ul) { *ec = DK3_ERROR_MATH_OVERFLOW; } } #line 1225 "dk3xsp.ctr" return back; #else size_t back = 0; int mec = 0; #line 1230 "dk3xsp.ctr" back = dk3ma_d_to_sz_ok(floor(tstart), &mec); if (mec) { if (ec) { *ec = mec; } } #line 1233 "dk3xsp.ctr" return back; #endif } int dk3xsp_calculate_position( dk3_xspline_segment_t *pseg, dk3_fig_spline_point_t const *points, size_t np, double t, int cb, int der, int *ec ) { size_t i = 0; int back = 0; int mec = 0; if((pseg) && (points) && (np > 1) && ((t >= 0.0) && (t <= (double)(np - 1)))) { i = dk3xsp_find_segment_start(t, &mec); if(!(mec)) { dk3xsp_reset(pseg); dk3xsp_set_cb(pseg, cb); if(i <= (np - 2)) { /* Not final point */ dk3xsp_set( pseg, ((i > 0) ? &(points[i - 1]) : NULL), &(points[i]), &(points[i+1]), ((i < (np - 2)) ? &(points[i + 2]) : NULL) ); back = dk3xsp_calculate(pseg, (t - (double)i), der); } else { /* In final point */ dk3xsp_set( pseg, ((np > 2) ? &(points[np - 3]): NULL), &(points[np - 2]), &(points[np - 1]), NULL ); back = dk3xsp_calculate(pseg, 1.0, der); } } else { if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS;} } return back; } double dk3xsp_open_spline_partial_length( dk3_fig_spline_point_t const *points, size_t np, double prec, double tstart, double *fslc, int dir, int cb, int *ec ) { dk3_xspline_segment_t seg; double back = -1.0; double val; size_t startseg; size_t i; int mec = 0; #line 1312 "dk3xsp.ctr" if((points) && (1 < np) && (prec > 0.0)) { if((tstart >= 0.0) && (tstart <= ((double)(np - 1)))) { dk3xsp_reset(&seg); if(cb) { dk3xsp_set_cb(&seg, 1); } startseg = dk3xsp_find_segment_start(tstart, &mec); if(mec) { #line 1318 "dk3xsp.ctr" if(ec) { *ec = mec; } } else { #line 1320 "dk3xsp.ctr" if(startseg == (np - 1)) { #line 1321 "dk3xsp.ctr" /* tstart points to the final point. So we have 0 for forward length or the full spline length for backward length. */ if(dir) { #line 1327 "dk3xsp.ctr" back = dk3xsp_open_spline_length(points, np, prec, fslc, cb, ec); } else { #line 1329 "dk3xsp.ctr" back = 0.0; } } else { #line 1332 "dk3xsp.ctr" if(tstart > 0.0) { #line 1333 "dk3xsp.ctr" /* tstart points to somewhere in the spline. We have to calculate one segment partially, all other segments completely. */ dk3xsp_set( &seg, ((startseg > 0) ? &(points[startseg-1]) : NULL), &(points[startseg]), &(points[startseg+1]), ((startseg < (np - 2)) ? &(points[startseg+2]) : NULL) ); val = dk3xsp_segment_partial_length( &seg,(tstart - (double)startseg),dir,(prec / (double)(np - 1)),ec ); if(val >= 0.0) { #line 1349 "dk3xsp.ctr" back = val; /* Partial segment calculated successfully. Now calculate other segments. */ if(dir) { #line 1355 "dk3xsp.ctr" /* Calculate all segments before the partial one. */ for(i = 0; ((i < startseg) && (back > -1.0)); i++) { val = -1.0; #line 1360 "dk3xsp.ctr" if(fslc) { val = fslc[i]; } if(val < -0.5) { dk3xsp_set( &seg, ((i > 0) ? &(points[i-1]) : NULL), &(points[i]), &(points[i+1]), ((i < (np - 2)) ? &(points[i+2]) : NULL) ); val = dk3xsp_segment_length( &seg, (prec / (double)(np - 1)), ec ); if(val >= 0.0) { if(fslc) { fslc[i] = val; } } } if(val >= 0.0) { #line 1381 "dk3xsp.ctr" back = dk3ma_d_add_ok(back, val, &mec); if(mec) { back = -1.0; i = np; if(ec) { *ec = mec; } } } else { back = -1.0; i = np; } } } else { #line 1391 "dk3xsp.ctr" /* Calculate all segments after the partial one. */ for(i = (startseg + 1); ((i < (np - 1)) && (back > -1)); i++) { val = -1.0; #line 1396 "dk3xsp.ctr" if(fslc) { val = fslc[i]; } if(val < -0.5) { dk3xsp_set( &seg, ((i > 0) ? &(points[i-1]) : NULL), &(points[i]), &(points[i+1]), ((i < (np - 2)) ? &(points[i+2]) : NULL) ); val = dk3xsp_segment_length( &seg, (prec / (double)(np - 1)), ec ); if(val >= 0.0) { if(fslc) { fslc[i] = val; } } } if(val >= 0.0) { #line 1417 "dk3xsp.ctr" back = dk3ma_d_add_ok(back, val, &mec); if(mec) { back = -1.0; i = np; if(ec) { *ec = mec; } } } else { #line 1423 "dk3xsp.ctr" back = -1.0; i = np; } } } #line 1427 "dk3xsp.ctr" } else { #line 1428 "dk3xsp.ctr" } } else { #line 1430 "dk3xsp.ctr" /* tstart points to the start point. So we have the full spline length for forward length, 0 for backward length. */ if(0 == dir) { #line 1436 "dk3xsp.ctr" back = dk3xsp_open_spline_length(points, np, prec, fslc, cb, ec); } else { #line 1438 "dk3xsp.ctr" back = 0.0; } } } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1449 "dk3xsp.ctr" return back; } double dk3xsp_get_t_for_distance_bisection( dk3_fig_spline_point_t const *points, size_t np, double dist, double prec, double *fslc, int dir, int cb, int *ec ) { double back = -1.0; double xl; /* Left x. */ double xr; /* Right x. */ double xc; /* Center x. */ double yl; /* Left y. */ double yr; /* Right y. */ double yc; /* Center y. */ double l; /* Spline length. */ double fprec; /* factor * precision. */ double ayl = 0.0; /* abs yl. */ double ayr = 0.0; /* abs yr. */ unsigned long passno = 0UL; /* Current pass. */ int mec = 0; /* Error code variable. */ int cc = 1; /* Flag: Can continue. */ int mc = 1; /* Flag: Must continue. */ #line 1482 "dk3xsp.ctr" fprec = 0.2 * prec; if((points) && (1 < np) && (prec > 0.0) && (dist >= 0.0)) { l = dk3xsp_open_spline_length(points, np, fprec, fslc, cb, &mec); #line 1486 "dk3xsp.ctr" if(l > dist) { xl = 0.0; xr = (double)(np - 1); if(dir) { yl = 0.0 - dist; yr = dk3ma_d_sub_ok(l, dist, &mec); #line 1492 "dk3xsp.ctr" } else { yl = dk3ma_d_sub_ok(l, dist, &mec); #line 1494 "dk3xsp.ctr" yr = 0.0 - dist; } while((cc) && (mc) && (!(mec))) { xc = 0.5 * dk3ma_d_add_ok(xl, xr, &mec); #line 1499 "dk3xsp.ctr" yc = dk3ma_d_sub_ok( dk3xsp_open_spline_partial_length( points, np, fprec, xc, fslc, dir, cb, ec ), dist, &mec ); #line 1506 "dk3xsp.ctr" #line 1507 "dk3xsp.ctr" #line 1508 "dk3xsp.ctr" #line 1509 "dk3xsp.ctr" if(dir) { if(yc > 0.0) { xr = xc; yr = yc; } else { xl = xc; yl = yc; } } else { if(yc > 0.0) { xl = xc; yl = yc; } else { xr = xc; yr = yc; } } ayl = fabs(yl); ayr = fabs(yr); if((ayl < prec) && (ayr < prec)) { mc = 0; } else { passno++; #line 1528 "dk3xsp.ctr" if(passno > 1000UL) { #line 1529 "dk3xsp.ctr" mec = DK3_ERROR_ITERATION; } else { } } } if((cc) && (!(mec))) { #line 1535 "dk3xsp.ctr" back = xl + ((xr - xl) * ayl)/(ayl + ayr); } if(mec) { #line 1538 "dk3xsp.ctr" back = -1.0; if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1547 "dk3xsp.ctr" return back; } #if DK3_HAVE_ISFINITE #define MY_ISFINITE(x) (isfinite(x)) #else #if DK3_HAVE__FINITE #define MY_ISFINITE(x) (_finite(x)) #else #define MY_ISFINITE(x) (1) #endif #endif #if DK3_HAVE_ISLESS #define MY_ISLESS(a,b) (isless((a),(b))) #else #define MY_ISLESS(a,b) (((MY_ISFINITE(a)) && (MY_ISFINITE(b))) && ((a) < (b))) #endif static double dk3xsp_gamma(double fg, double fx) { double back; back = 1.0 - fx / fg; if (!(MY_ISFINITE(back))) { back = 0.5; } else { if (0.0 >= back) { back = 0.5; } else { if (1.0 < back) { back = 1.0; } } } return back; } double dk3xsp_get_t_for_distance( dk3_fig_spline_point_t const *points, size_t np, double dist, double prec, double *fslc, int dir, int cb, int *ec ) { double back = -1.0; double x = 0.0; double xo; double a; double b; double ya; double yb; double yx; double fprec; double l; double gamma; unsigned long passno = 0UL; int ls = 0; int mec = 0; int cc = 1; #line 1619 "dk3xsp.ctr" fprec = prec / 8.0; if((points) && (1 < np) && (prec > 0.0) && (dist >= 0.0)) { l = dk3xsp_open_spline_length(points, np, fprec, fslc, cb, &mec); if ((l > dist) && (0 == mec)) { a = 0.0; b = (double)(np - 1); if(dir) { ya = 0.0 - dist; yb = dk3ma_d_sub_ok(l, dist, &mec); } else { ya = dk3ma_d_sub_ok(l, dist, &mec); yb = 0.0 - dist; } if (0 == mec) { while (1 == cc) { #line 1635 "dk3xsp.ctr" #line 1636 "dk3xsp.ctr" #line 1637 "dk3xsp.ctr" xo = x; x = (a * yb - b * ya) / (yb - ya); #line 1640 "dk3xsp.ctr" if (MY_ISFINITE(x)) { yx = dk3ma_d_sub_ok( dk3xsp_open_spline_partial_length( points, np, fprec, x, fslc, dir, cb, ec ), dist, &mec ); #line 1647 "dk3xsp.ctr" if (MY_ISFINITE(yx)) { if ((passno > 0UL) && (MY_ISLESS((x-xo),prec)) && (fabs(yx) < prec)) { cc = 0; #line 1651 "dk3xsp.ctr" back = x; } else { if (passno >= 1000UL) { #line 1656 "dk3xsp.ctr" cc = -1; /* ##### ERROR: TOO MANY PASSES */ } else { if (dir) { if (yx > 0.0) { #line 1661 "dk3xsp.ctr" if (1 == ls) { #line 1662 "dk3xsp.ctr" gamma = dk3xsp_gamma(yb, yx); ya *= gamma; } b = x; yb = yx; ls = 1; } else { #line 1668 "dk3xsp.ctr" if (-1 == ls) { #line 1669 "dk3xsp.ctr" gamma = dk3xsp_gamma(ya, yx); yb *= gamma; } a = x; ya = yx; ls = -1; } } else { if (yx > 0.0) { #line 1677 "dk3xsp.ctr" if (-1 == ls) { #line 1678 "dk3xsp.ctr" gamma = dk3xsp_gamma(ya, yx); yb *= gamma; } a = x; ya = yx; ls = -1; } else { #line 1684 "dk3xsp.ctr" if (1 == ls) { #line 1685 "dk3xsp.ctr" gamma = dk3xsp_gamma(yb, yx); ya *= gamma; } b = x; yb = yx; ls = 1; } } } } } else { #line 1695 "dk3xsp.ctr" cc = -1; /* ##### MATH ERROR */ } } else { #line 1699 "dk3xsp.ctr" cc = -1; /* ##### MATH ERROR */ } if (ULONG_MAX > passno) { passno++; } } } else { #line 1705 "dk3xsp.ctr" /* ##### ERROR: Math error */ } } else { #line 1708 "dk3xsp.ctr" if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } } else { #line 1711 "dk3xsp.ctr" if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1714 "dk3xsp.ctr" return back; } #line 1719 "dk3xsp.ctr" double dk3xsp_direction( dk3_fig_spline_point_t const *points, size_t np, double t, int cb, int *ec ) { dk3_xspline_segment_t seg; double back = -10.0 * M_PI; double dxdt; double dydt; double delta; double tleft; double tright; double xleft; double xright; double yleft; double yright; int mec = 0; #line 1743 "dk3xsp.ctr" if((points) && (1 < np)) { if((t >= 0.0) && (t <= ((double)(np - 1)))) { dk3xsp_reset(&seg); dk3xsp_set_cb(&seg, ((cb) ? 1 : 0)); if(dk3xsp_calculate_position(&seg, points, np, t, cb, 1, &mec)) { dxdt = dk3xsp_get_dxdt(&seg); dydt = dk3xsp_get_dydt(&seg); if((fabs(dxdt) > 1.0e-6) || (fabs(dydt) > 1.0e-6)) { back = dk3ma_d_atan2(dydt, dxdt); } else { #line 1753 "dk3xsp.ctr" delta = 0.001; tleft = tright = t; if(t > delta) { tleft = t - delta; } if(t < (((double)(np - 1)) - delta)) { tright = t + delta; } if(dk3xsp_calculate_position(&seg, points, np, tleft, cb, 0, &mec)) { xleft = dk3xsp_get_x(&seg); yleft = dk3xsp_get_y(&seg); if(dk3xsp_calculate_position(&seg, points, np, tright, cb, 0, &mec)) { xright = dk3xsp_get_x(&seg); yright = dk3xsp_get_y(&seg); dxdt = dk3ma_d_sub_ok(xright, xleft, &mec); dydt = dk3ma_d_sub_ok(yright, yleft, &mec); back = dk3ma_d_atan2(dydt, dxdt); } else { #line 1767 "dk3xsp.ctr" } } else { #line 1769 "dk3xsp.ctr" } } } else { #line 1772 "dk3xsp.ctr" } if(mec) { if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1780 "dk3xsp.ctr" return back; } double dk3xsp_direction_and_position( double *xpos, double *ypos, dk3_fig_spline_point_t const *points, size_t np, double t, int cb, int *ec ) { dk3_xspline_segment_t seg; double back = -10.0 * M_PI; double dxdt; double dydt; double delta; double tleft; double tright; double xleft; double xright; double yleft; double yright; int mec = 0; #line 1809 "dk3xsp.ctr" if((points) && (1 < np) && (xpos) && (ypos)) { if((t >= 0.0) && (t <= ((double)(np - 1)))) { dk3xsp_reset(&seg); dk3xsp_set_cb(&seg, ((cb) ? 1 : 0)); if(dk3xsp_calculate_position(&seg, points, np, t, cb, 1, &mec)) { *xpos = dk3xsp_get_x(&seg); *ypos = dk3xsp_get_y(&seg); dxdt = dk3xsp_get_dxdt(&seg); dydt = dk3xsp_get_dydt(&seg); #line 1818 "dk3xsp.ctr" if((fabs(dxdt) > 1.0e-6) || (fabs(dydt) > 1.0e-6)) { back = dk3ma_d_atan2(dydt, dxdt); #line 1820 "dk3xsp.ctr" } else { #line 1821 "dk3xsp.ctr" delta = 0.001; tleft = tright = t; if(t > delta) { tleft = t - delta; } if(t < (((double)(np - 1)) - delta)) { tright = t + delta; } if(dk3xsp_calculate_position(&seg, points, np, tleft, cb, 0, &mec)) { xleft = dk3xsp_get_x(&seg); yleft = dk3xsp_get_y(&seg); if(dk3xsp_calculate_position(&seg, points, np, tright, cb, 0, &mec)) { xright = dk3xsp_get_x(&seg); yright = dk3xsp_get_y(&seg); dxdt = dk3ma_d_sub_ok(xright, xleft, &mec); dydt = dk3ma_d_sub_ok(yright, yleft, &mec); back = dk3ma_d_atan2(dydt, dxdt); } else { #line 1835 "dk3xsp.ctr" } } else { #line 1837 "dk3xsp.ctr" } } } else { #line 1840 "dk3xsp.ctr" } if(mec) { if(ec) { *ec = mec; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } } else { if(ec) { *ec = DK3_ERROR_INVALID_ARGS; } } #line 1848 "dk3xsp.ctr" return back; }