%% options copyright owner = Dirk Krause copyright year = 2018-xxxx SPDX-License-Identifier: BSD-3-Clause %% header /** @file wxd2lxs.h X-spline handling in wxd2lat. */ #ifndef DK4ERROR_H_INCLUDED #include "dk4error.h" #endif #ifndef WXDTYPES_H_INCLUDED #include "wxdtypes.h" #endif #ifdef __cplusplus extern "C" { #endif /** Set up X-spline segment structure from X-spline data object. @param xs X-spline segment structure to set up. @param pobj X-spline object to obtain data from. @param ib Segment number. */ void wxd2lxs_set_xspline( dk4_xsp_2d_t *xs, wxd_object_t *pobj, uint16_t ib ); /** Check spline object, if necessary calculate the segment lengths. @param job Job structure. @param drw Drawing structure. @param pobj Object to process. @paam backptr Address of success variable to reset on error. */ void wxd2lxs_calculate_spline_segments( wxd2lat_job_t *job, wxd_object_t *pobj, int *backptr ); /** Calculate coordinates and optionally direction for a given t value. This function can only be used for open splines with at least one arrowhead attached. @param rp Address of result variable to set on success. If d=0 a 2 elements array is required to store x and y. If d!=0 a 3 elements array is required to store x, y and alpha (direction angle). @param pobj X-spline object. @param t Parameter t. @param d Flag: Calculate derivative values. @param erp Error report, may be NULL. @return 1 on success, 0 on error (calculation failed). */ int wxd2lxs_calculate( double *rp, wxd_object_t *pobj, double t, int d, dk4_er_t *erp ); /** Calculate partial X-spline length from 0 to a given t value. This function can only be used for open splines with at least one arrowhead attached. @param rp Address of result variable to set on success. @param pobj X-spline object. @param t Parameter t. @param eps Required precision. @param erp Error report, may be NULL. @return 1 on success, 0 on error (calculation failed). */ int wxd2lxs_length_to_parameter( double *rp, wxd_object_t *pobj, double t, double eps, dk4_er_t *erp ); /** Find t parameter for a given partial X-spline length l. @param rp Address of result variable for t value. @param pobj Open X-spline with arrowhead(s). @param l Length to search to for. @param xslp Precision for length calculation. @param xspp Iteration precision for t. @param erp Error report, may be NULL. @return 1 on success, 0 on error. */ int wxd2lxs_parameter_for_length( double *rp, wxd_object_t *pobj, double l, double xslp, double xspp, dk4_er_t *erp ); #ifdef __cplusplus } #endif /* vim: set ai sw=4 ts=4 : */ %% module #ifndef WXD2LAT_H_INCLUDED #include "wxd2lat.h" #endif $!trace-include void wxd2lxs_set_xspline( dk4_xsp_2d_t *xs, wxd_object_t *pobj, uint16_t ib ) { wxd_spline_point_t *p; uint16_t n; uint16_t ia; uint16_t ic; uint16_t id; if ((NULL != xs) && (NULL != pobj)) { if ( ((WXD_OT_C_SPLINE == pobj->ot) && (ib < (pobj->det).s.n)) || ((WXD_OT_O_SPLINE) && (ib < ((pobj->det).s.n) - (uint16_t)1U)) ) { p = (pobj->det).s.p; n = (pobj->det).s.n; ia = ib - 1; ic = ib + 1; id = ib + 2; if ((uint16_t)0U == ib) { ia = n - (uint16_t)1U; } else { if (n - (uint16_t)1U == ib) { ic = (uint16_t)0U; id = (uint16_t)1U; } else { if (n - (uint16_t)2U == ib) { id = (uint16_t)0U; } } } dk4xsp2d_reset(xs); #if 0 dk4xsp2d_reset_points(xs); #endif dk4xsp2d_set_point(xs, p[ib].x, p[ib].y, p[ib].s, 1); dk4xsp2d_set_point(xs, p[ic].x, p[ic].y, p[ic].s, 2); if (((uint16_t)0U < ib) || (WXD_OT_C_SPLINE == pobj->ot)) { dk4xsp2d_set_point(xs, p[ia].x, p[ia].y, p[ia].s, 0); } if ((n - (uint16_t)2U > ib) || (WXD_OT_C_SPLINE == pobj->ot)) { dk4xsp2d_set_point(xs, p[id].x, p[id].y, p[id].s, 3); } } } } int wxd2lxs_calculate( double *rp, wxd_object_t *pobj, double t, int d, dk4_er_t *erp ) { dk4_xsp_2d_t xs; double val[4]; double t1; double t2; size_t ib; #if VERSION_BEFORE_20190110 size_t ia; size_t ic; size_t id; #endif int back = 0; int res; $? "+ wxd2lxs_calculate" /* Check arguments */ if ((NULL == rp) || (NULL == pobj) || (WXD_OT_O_SPLINE != pobj->ot)) { dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if ((0.0 > t) || ((double)((pobj->det).s.n) < t)) { dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if (NULL == (pobj->det).s.sl) { dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } /* Calculate indices */ ib = (size_t)floor(t); t = t - (double)ib; #if VERSION_BEFORE_20190110 ic = ib + 1; id = ib + 2; ia = ib - 1; /* Set up X-spline */ dk4xsp2d_reset(&xs); dk4xsp2d_reset_points(&xs); if ((size_t)0U < ib) { dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ia].x, ((pobj->det).s.p)[ia].y, ((pobj->det).s.p)[ia].s, 0 ); } dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ib].x, ((pobj->det).s.p)[ib].y, ((pobj->det).s.p)[ib].s, 1 ); dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ic].x, ((pobj->det).s.p)[ic].y, ((pobj->det).s.p)[ic].s, 2 ); if ((pobj->det).s.n - 2 > ib) { dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[id].x, ((pobj->det).s.p)[id].y, ((pobj->det).s.p)[id].s, 3 ); } #else wxd2lxs_set_xspline(&xs, pobj, ib); #endif /* Retrieve values */ if (0 != d) { res = dk4xsp2d_calculate_value_derivative(val, &xs, t, erp); } else { res = dk4xsp2d_calculate_value(val, &xs, t, erp); } if (0 == res) { goto finished; } /* Transfer x and y */ if (0 != d) { rp[0] = val[0]; rp[1] = val[2]; if ((1.0e-8 < fabs(val[1])) || (1.0e-8 < fabs(val[3]))) { rp[2] = atan2(val[3], val[1]); back = 1; if (0 == dk4ma_is_finite(rp[2])) { t1 = t2 = t; if ((1.0/64.0) < t) { t1 = t - (1.0/64.0); } if ((1.0 - 1.0/64.0) > t) { t2 = t + (1.0/64.0); } if (0 != dk4xsp2d_calculate_value(val, &xs, t1, erp)) { if (0 != dk4xsp2d_calculate_value(&(val[2]),&xs,t2,erp)) { rp[2] = atan2((val[3] - val[1]), (val[2] - val[0])); back = 1; } } } } else { t1 = t2 = t; if ((1.0/64.0) < t) { t1 = t - (1.0/64.0); } if ((1.0 - 1.0/64.0) > t) { t2 = t + (1.0/64.0); } if (0 != dk4xsp2d_calculate_value(val, &xs, t1, erp)) { if (0 != dk4xsp2d_calculate_value(&(val[2]), &xs, t2, erp)) { rp[2] = atan2((val[3] - val[1]), (val[2] - val[0])); back = 1; } } } } else { rp[0] = val[0]; rp[1] = val[1]; back = 1; } finished: $? "- wxd2lxs_calculate %d", back return back; } int wxd2lxs_length_to_parameter( double *rp, wxd_object_t *pobj, double t, double eps, dk4_er_t *erp ) { dk4_xsp_2d_t xs; double v = 0.0; double pv = 0.0; size_t i; #if VERSION_BEFORE_20190110 size_t ia; size_t ic; size_t id; #endif size_t ib; int back = 0; $? "+ wxd2lxs_length_to_parameter %g", t /* Check arguments */ if ((NULL == rp) || (NULL == pobj) || (WXD_OT_O_SPLINE != pobj->ot)) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if ((0.0 > t) || ((double)((pobj->det).s.n) < t)) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if (NULL == (pobj->det).s.sl) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } /* Sum up complete segments. */ back = 1; ib = (size_t)floor(t); for (i = (size_t)0U; i < ib; i++) { v += ((pobj->det).s.sl)[i]; } $? ". sum of %u segments: %g", (unsigned)ib, v /* Calculate partial segment. */ t -= (double)ib; if (0.0 < t) { #if VERSION_BEFORE_20190110 ia = ib - (size_t)1U; ic = ib + (size_t)1U; id = ib + (size_t)2U; dk4xsp2d_reset(&xs); dk4xsp2d_reset_points(&xs); dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ib].x, ((pobj->det).s.p)[ib].y, ((pobj->det).s.p)[ib].s, 1 ); dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ic].x, ((pobj->det).s.p)[ic].y, ((pobj->det).s.p)[ic].s, 2 ); if ((size_t)0U < ib) { dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[ia].x, ((pobj->det).s.p)[ia].y, ((pobj->det).s.p)[ia].s, 0 ); } if ((size_t)((pobj->det).s.n) - (size_t)2U > ib) { dk4xsp2d_set_point( &xs, ((pobj->det).s.p)[id].x, ((pobj->det).s.p)[id].y, ((pobj->det).s.p)[id].s, 3 ); } #else wxd2lxs_set_xspline(&xs, pobj, ib); #endif back = dk4xsp2d_calculate_partial_length(&pv, &xs, t, eps, erp); if (0 != back) { v += pv; $? ". partial segment length %g", pv } #if TRACE_DEBUG else { $? "! calculation" } #endif } if (0 != back) { if (0 != dk4ma_is_finite(v)) { *rp = v; $? ". result = %g", v } else { back = 0; dk4error_set_simple_error_code(erp, DK4_E_MATH_OVERFLOW); } } finished: $? "- wxd2lxs_length_to_parameter %d", back return back; } /** Function used in iteration. @param d Address of destination variable for function value at x. @param x X position (the t parameter). @param ps A wxd2lxs_iter_t containing X-spline object, length and precision. @return 1 on success, 0 on error. */ static int iter_fct(double *d, double x, const void *ps) { const wxd2lxs_iter_t *piter; double v; int back = 0; piter = (const wxd2lxs_iter_t *)ps; back = wxd2lxs_length_to_parameter( &v, piter->pobj, x, piter->xslp, piter->erp ); if (0 != back) { *d = v - piter->l; } return back; } int wxd2lxs_parameter_for_length( double *rp, wxd_object_t *pobj, double l, double xslp, double xspp, dk4_er_t *erp ) { wxd2lxs_iter_t iter; dk4_iter_ctx_t ctx; unsigned long pp = 0UL; int back = 0; int res; $? "+ wxd2lxs_parameter_for_length %g", l /* Check arguments. */ if ((NULL == rp) || (NULL == pobj)) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if ((0.0 >= xslp) || (0.0 >= xspp)) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } if ((0.0 > l) || ((pobj->det).s.suml < l)) { $? "! args" dk4error_set_simple_error_code(erp, DK4_E_INVALID_ARGUMENTS); goto finished; } /* Set up iteration structures */ iter.pobj = pobj; iter.l = l; iter.xslp = xslp; iter.erp = erp; dk4iter_ctx_init(&ctx); dk4iter_ctx_set_eps_x(&ctx, xspp); dk4iter_ctx_set_eps_y(&ctx, xspp); dk4iter_ctx_set_maxpass(&ctx, 1024UL); dk4iter_ctx_set_exact(&ctx, 0); dk4iter_ctx_set_algorithm(&ctx, DK4_ITER_ALG_RF_ANDERSON_BJOERCK); dk4iter_ctx_set_min(&ctx, 0.0); dk4iter_ctx_set_max(&ctx, (double)((pobj->det).s.n - (uint16_t)1U)); /* Run iteration */ $? ". run iteration" res = dk4iter_interval( rp, &pp, iter_fct, (const void *)(&iter), 0.0, (double)((pobj->det).s.n - (uint16_t)1U), &ctx ); /* Process iteration result */ if (DK4_ITER_RESULT_SUCCESS == res) { back = 1; $? ". success after %lu passes", pp } #if TRACE_DEBUG else { $? "! iteration failed %d", res } #endif finished: $? "- wxd2lxs_parameter_for_length %d", back return back; } void wxd2lxs_calculate_spline_segments( wxd2lat_job_t *job, wxd_object_t *pobj, int *backptr ) { dk4_xsp_2d_t xsp; dk4_er_t er; size_t nelem; uint16_t i; int res; int failed; $? "+ wxd2lxs_calculate_spline_segments" failed = 0; if (((uint16_t)1U < (pobj->det).s.n) && (NULL == (pobj->det).s.sl)) { nelem = (size_t)((pobj->det).s.n - (uint16_t)1U); (pobj->det).s.sl = dk4mem_new(double,nelem,NULL); if (NULL != (pobj->det).s.sl) { dk4xsp2d_reset(&xsp); (pobj->det).s.suml = 0.0; for (i = 0; i < nelem; i++) { ((pobj->det).s.sl)[i] = 0.0; /* Calculate segment length */ #if VERSION_BEFORE_20190110 dk4xsp2d_reset_points(&xsp); dk4xsp2d_set_point( &xsp, ((pobj->det).s.p)[i].x, ((pobj->det).s.p)[i].y, ((pobj->det).s.p)[i].s, 1 ); dk4xsp2d_set_point( &xsp, ((pobj->det).s.p)[i+1].x, ((pobj->det).s.p)[i+1].y, ((pobj->det).s.p)[i+1].s, 2 ); if (0 < i) { $? ". not left outer segment" dk4xsp2d_set_point( &xsp, ((pobj->det).s.p)[i-1].x,((pobj->det).s.p)[i-1].y, ((pobj->det).s.p)[i-1].s, 0 ); } if (nelem - 1 > i) { $? ". not right outer segment" dk4xsp2d_set_point( &xsp, ((pobj->det).s.p)[i+2].x,((pobj->det).s.p)[i+2].y, ((pobj->det).s.p)[i+2].s, 3 ); } #else wxd2lxs_set_xspline(&xsp, pobj, i); #endif res = dk4xsp2d_calculate_length( &(((pobj->det).s.sl)[i]), &xsp, (job->grco).xslp, NULL ); if (0 != res) { $? ". seg %u length %g", i, ((pobj->det).s.sl)[i] (pobj->det).s.suml += ((pobj->det).s.sl)[i]; } else { $? "! segment length calculation failed" *backptr = 0; failed = 1; /* ERROR: Segment length calculation failed */ dk4app_log_1( job->app, job->msg, job->sz_msg, DK4_LL_ERROR, 79 ); } } $? ". spline length = %g", (pobj->det).s.suml if (0 == failed) { (pobj->det).s.ts = 0.0; (pobj->det).s.te = (pobj->det).s.suml; if (0 != (pobj->ab).type) { dk4error_init(&er); res = wxd2lxs_parameter_for_length( &((pobj->det).s.ts), pobj, (pobj->ab).cut, (job->grco).xslp, (job->grco).xspp, &er ); if (0 == res) { $? "! cut failed" (pobj->det).s.ts = -1.0; *backptr = 0; /* ERROR: Cut calculation failed */ dk4app_log_1( job->app, job->msg, job->sz_msg, DK4_LL_ERROR, 80 ); } else { $? ". ts = %g", (pobj->det).s.ts } } if (0 != (pobj->af).type) { dk4error_init(&er); res = wxd2lxs_parameter_for_length( &((pobj->det).s.te), pobj, ((pobj->det).s.suml - (pobj->af).cut), (job->grco).xslp, (job->grco).xspp, &er ); if (0 == res) { $? "! cut failed" (pobj->det).s.te = -1.0; *backptr = 0; /* ERROR: Cut calculation failed */ dk4app_log_1( job->app, job->msg, job->sz_msg, DK4_LL_ERROR, 80 ); } else { $? ". te = %g", (pobj->det).s.te } } } } else { $? "! memory" /* ERROR: Memory allocation failed */ dk4app_log_1( job->app, job->msg, job->sz_msg, DK4_LL_ERROR, 81 ); *backptr = 0; } } #if TRACE_DEBUG else { $? "! args" } #endif $? "- wxd2lxs_calculate_spline_segments" } /* vim: set ai sw=4 ts=4 : */