23 if (! globalEulerSpiralLookupTable)
51 double sangle, eangle;
53 if (start_angle>
M_PI) sangle = start_angle-2*
M_PI;
54 else sangle = start_angle;
56 if (end_angle>M_PI) eangle = end_angle-2*
M_PI;
57 else eangle = end_angle;
60 int ilow, ihigh, jlow, jhigh;
62 ilow = (
int)floor((sangle+M_PI)/
_dt);
63 ihigh = (
int)ceil((sangle+M_PI)/
_dt);
65 jlow = (
int)floor((eangle+M_PI)/
_dt);
66 jhigh = (
int)ceil((eangle+M_PI)/
_dt);
68 double slow =
_theta[ilow];
69 double elow =
_theta[jlow];
71 double a = (sangle - slow)/
_dt;
72 double b = (eangle - elow)/
_dt;
74 double k0 = (1-a)*(1-b)*
ES_k0[ilow][jlow] + a*(1-b)*
ES_k0[ihigh][jlow] +
75 (1-a)*b*
ES_k0[ilow][jhigh] + a*b*
ES_k0[ihigh][jhigh];
85 double sangle, eangle;
87 if (start_angle>
M_PI) sangle = start_angle-2*
M_PI;
88 else sangle = start_angle;
90 if (end_angle>M_PI) eangle = end_angle-2*
M_PI;
91 else eangle = end_angle;
94 int ilow, ihigh, jlow, jhigh;
96 ilow = (
int)floor((sangle+M_PI)/
_dt);
97 ihigh = (
int)ceil((sangle+M_PI)/
_dt);
99 jlow = (
int)floor((eangle+M_PI)/
_dt);
100 jhigh = (
int)ceil((eangle+M_PI)/
_dt);
102 double slow =
_theta[ilow];
103 double elow =
_theta[jlow];
105 double a = (sangle - slow)/
_dt;
106 double b = (eangle - elow)/
_dt;
108 double k1 = (1-a)*(1-b)*
ES_k1[ilow][jlow] + a*(1-b)*
ES_k1[ihigh][jlow] +
109 (1-a)*b*
ES_k1[ilow][jhigh] + a*b*
ES_k1[ihigh][jhigh];
119 double sangle, eangle;
121 if (start_angle>
M_PI) sangle = start_angle-2*
M_PI;
122 else sangle = start_angle;
124 if (end_angle>M_PI) eangle = end_angle-2*
M_PI;
125 else eangle = end_angle;
128 int ilow, ihigh, jlow, jhigh;
130 ilow = (
int)floor((sangle+M_PI)/
_dt);
131 ihigh = (
int)ceil((sangle+M_PI)/
_dt);
133 jlow = (
int)floor((eangle+M_PI)/
_dt);
134 jhigh = (
int)ceil((eangle+M_PI)/
_dt);
136 double slow =
_theta[ilow];
137 double elow =
_theta[jlow];
139 double a = (sangle - slow)/
_dt;
140 double b = (eangle - elow)/
_dt;
153 double sangle, eangle;
155 if (start_angle>
M_PI) sangle = start_angle-2*
M_PI;
156 else sangle = start_angle;
158 if (end_angle>M_PI) eangle = end_angle-2*
M_PI;
159 else eangle = end_angle;
162 int ilow, ihigh, jlow, jhigh;
164 ilow = (
int)floor((sangle+M_PI)/
_dt);
165 ihigh = (
int)ceil((sangle+M_PI)/
_dt);
167 jlow = (
int)floor((eangle+M_PI)/
_dt);
168 jhigh = (
int)ceil((eangle+M_PI)/
_dt);
170 double slow =
_theta[ilow];
171 double elow =
_theta[jlow];
173 double a = (sangle - slow)/
_dt;
174 double b = (eangle - elow)/
_dt;
176 double L = (1-a)*(1-b)*
ES_L[ilow][jlow] + a*(1-b)*
ES_L[ihigh][jlow] +
177 (1-a)*b*
ES_L[ilow][jhigh] + a*b*
ES_L[ihigh][jhigh];
220 double prev_error = error;
222 double k0 = k0_init_est;
223 double L = L_init_est;
225 double e1, e2, e3, e4 = 0;
239 if (error>prev_error)
245 if (error==e1) k0 = k0 + dstep;
246 else if (error==e2) k0 = k0 - dstep;
247 else if (error==e3) L = L + dstep;
248 else if (error==e4) L = L - dstep;
264 if (ds==0 && NPts==0){
276 for (
int i=1; i<NPts; i++,s+=ds){
278 spiral.push_back(cur_pt);
322 end_pt.
setX(start_pt.
getX()+L*cos(theta));
323 end_pt.
setY(start_pt.
getY()+L*sin(theta));
328 double const_term = 1.0/k0;
329 end_pt.
setX(start_pt.
getX()+const_term*(sin(k0*L+theta)-sin(theta)));
330 end_pt.
setY(start_pt.
getY()-const_term*(cos(k0*L+theta)-cos(theta)));
335 double const1 = sqrt(
M_PI*fabs(gamma));
336 double const2 = sqrt(
M_PI/fabs(gamma));
341 double C = (fresnel1.
getX() - fresnel2.
getX());
345 double S = fresnel1.
getY() - fresnel2.
getY();
347 double cos_term = cos(theta-((k0*k0)/(2.0*gamma)));
348 double sin_term = sin(theta-((k0*k0)/(2.0*gamma)));
350 end_pt.
setX(start_pt.
getX() + const2*(C*cos_term - S*sin_term));
351 end_pt.
setY(start_pt.
getY() + const2*(C*sin_term + S*cos_term));
376 #define FPMIN 1.0e-30
384 double a,ax,fact,pix2,
sign,sum,sumc,sums,term,test;
386 std::complex<double> b,cc,d,h,del,cs;
390 if (ax < sqrt(
FPMIN))
403 fact=(
M_PI/2.0)*ax*ax;
408 for (k=1;k<=
MAXIT;k++)
424 if (term < test)
break;
430 std::cout <<
"series failed in fresnel" << std::endl;
438 b = std::complex<double>(1.0,-pix2);
439 cc = std::complex<double>(1.0/
FPMIN,0.0);
440 d=h = std::complex<double>(1.0,0.0)/b;
443 for (k=2;k<=
MAXIT;k++)
447 b= b+std::complex<double>(4.0,0.0);
448 d=(std::complex<double>(1.0,0.0)/((a*d)+b));
451 cc=(b+(std::complex<double>(a,0.0)/cc));
455 if ((fabs(del.real()-1.0)+fabs(del.imag())) <
EPS)
459 std::cout <<
"cf failed in frenel" << std::endl;
461 h=std::complex<double>(ax,-ax)*h;
462 cs=std::complex<double>(0.5,0.5)*(std::complex<double>(1.0,0.0) - std::complex<double>(cos(0.5*pix2),sin(0.5*pix2))*h );
464 result.
setX(cs.real());
465 result.
setY(cs.imag());
double L(double start_angle, double end_angle)
static EulerSpiralLookupTable * globalEulerSpiralLookupTable
void set_end_params(Point2D< double > end_pt, double end_angle)
double compute_error(double k0, double L)
void set_start_params(Point2D< double > start_pt, double start_angle)
~EulerSpiralLookupTable()
static EulerSpiralLookupTable * get_globalEulerSpiralLookupTable()
double angle0To2Pi(double angle)
double euc_distance(Point2D< point_type1 > pt1, Point2D< point_type2 > pt2)
void compute_biarc_params()
double CCW(double reference, double angle1)
void computeSpiral(std::vector< Point2D< double > > &spiral, double ds=0, int NPts=0)
double k0(double start_angle, double end_angle)
Point2D< double > get_fresnel_integral(double value)
Point2D< double > compute_es_point(EulerSpiralParams &es_params, double arclength, bool bNormalized=false)
Point2D< double > compute_end_pt(double arclength, bool bNormalized=false)
double gamma(double start_angle, double end_angle)
Point2D< double > start_pt
double k1(double start_angle, double end_angle)
#define MAX_NUM_ITERATIONS