OSDN > Pesquisar Software > ongeo > SCM > SVN > Commit

ongeo

R/O
SSH
HTTPS

ongeo: Commit

Commit MetaInfo

Revisão	52 (tree)
Hora	2019-08-11 14:39:09
Autor	mocchi_2012

Mensagem de Log

曲線長計算の追加

Mudança Sumário

added: trunk/LengthBezierCurve_NumericalIntegration.cpp (diff)
added: trunk/LengthBezierCurve_Subdivision.cpp (diff)
added: trunk/LengthParamNurbsCurve_AdaptiveQuadrate.cpp (diff)
added: trunk/NumericalIntegration.cpp (diff)
added: trunk/ONGEO_CurveSubdivision.hpp (diff)
added: trunk/ONGEO_NumericalIntegration.hpp (diff)
added: trunk/ONGEO_QuasiInterpolating.hpp (diff)
added: trunk/ONGEO_RootFinding.hpp (diff)
added: trunk/TessellateCurve.cpp (diff)

Diff

--- trunk/LengthBezierCurve_NumericalIntegration.cpp (nonexistent)

+++ trunk/LengthBezierCurve_NumericalIntegration.cpp (revision 52)

		@@ -0,0 +1,60 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+#include "ONGEO.h"
	4	+#include <cmath>
	5	+#include <cstdio>
	6	+
	7	+#include "ONGEO_NumericalIntegration.hpp"
	8	+
	9	+// Reference : B. Guenter, R. Parent,
	10	+// "Computing the arc length of parametric curves".
	11	+// IEEE Computer Graphics and Applications, Vol 10, No.3, May. 1990, pp.72-78.
	12	+
	13	+namespace{
	14	+ struct IntegralFunc{
	15	+ ON_BezierCurve bc;
	16	+ double operator()(double t) const{
	17	+ ON_3dPoint pt;
	18	+ ON_3dVector dt;
	19	+ bc.Ev1Der(t, pt, dt);
	20	+ return dt.Length();
	21	+ }
	22	+ };
	23	+}
	24	+
	25	+double ONGEO_LengthBezierCurve_NumericalIntegration(const ON_BezierCurve &bc, const double subdivided_prms, int cnt, double estimated_deviation){
	26	+ ONGEO_NumericalIntegration_GaussKronrod<5> integ(0, 0);
	27	+
	28	+ IntegralFunc func;
	29	+ func.bc = bc;
	30	+
	31	+ double t_prv = *subdivided_prms;
	32	+ double sum_l = 0, sum_err = 0;
	33	+ for (int i = 1; i < cnt; ++i){
	34	+ double t_cur = subdivided_prms[i];
	35	+ integ.SetRange(t_prv, t_cur);
	36	+ double err, l = integ.Calc(func, err);
	37	+ sum_l += l;
	38	+ sum_err += err;
	39	+ t_prv = t_cur;
	40	+ }
	41	+ if (estimated_deviation) *estimated_deviation = sum_err;
	42	+ return sum_l;
	43	+}
	44	+
	45	+double ONGEO_LengthBezierCurve_AdaptiveQuadrate(const ON_BezierCurve &bc, const double subdivided_prms, int cnt, double tolerance, double estimated_deviation){
	46	+ IntegralFunc func;
	47	+ func.bc = bc;
	48	+
	49	+ double t_prv = *subdivided_prms;
	50	+ double sum_l = 0, sum_err = 0;
	51	+ for (int i = 1; i < cnt; ++i){
	52	+ double t_cur = subdivided_prms[i];
	53	+ double err, l = ONGEO_AdaptiveQuadrature<ONGEO_NumericalIntegration_GaussKronrod<5>, IntegralFunc>(func, t_prv, t_cur, tolerance * (t_cur - t_prv), err);
	54	+ sum_l += l;
	55	+ sum_err += err;
	56	+ t_prv = t_cur;
	57	+ }
	58	+ if (estimated_deviation) *estimated_deviation = sum_err;
	59	+ return sum_l;
	60	+}

--- trunk/LengthBezierCurve_Subdivision.cpp (nonexistent)

+++ trunk/LengthBezierCurve_Subdivision.cpp (revision 52)

		@@ -0,0 +1,52 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+#include "ONGEO.h"
	4	+#include <cmath>
	5	+#include <cstdio>
	6	+
	7	+#include "ONGEO_CurveSubdivision.hpp"
	8	+
	9	+// Reference : Maharavo Randrianarivony,
	10	+// "Length Estimation of Rational Bezier Curves and Application to CAD Parametrization".
	11	+
	12	+inline void GetLengthBounds(const ON_BezierCurve &bc_cur, double &lb, double &ub){
	13	+ ON_3dPoint pt_start, pt_end;
	14	+ bc_cur.GetCV(0, pt_start);
	15	+ bc_cur.GetCV(bc_cur.CVCount()-1, pt_end);
	16	+ lb = pt_end.DistanceTo(pt_start);
	17	+ ub = bc_cur.ControlPolygonLength();
	18	+}
	19	+
	20	+struct LengthOperator{
	21	+ // input
	22	+ double tolerance;
	23	+
	24	+ // output
	25	+ double sum_length;
	26	+ LengthOperator(double tolerance_) : tolerance(tolerance_), sum_length(0){
	27	+ }
	28	+ bool TestConverge(const ON_BezierCurve &bc_cur, const ON_Interval &t_int){
	29	+ double lb, ub;
	30	+ GetLengthBounds(bc_cur, lb, ub);
	31	+ bool converged = ((ub - lb) < tolerance * t_int.Length());
	32	+ if (converged){
	33	+ sum_length += 0.5 * (ub + lb);
	34	+ }
	35	+ return converged;
	36	+ }
	37	+};
	38	+
	39	+double ONGEO_LengthBezierCurve_SimpleSubdivision(const ON_BezierCurve &bc, double tolerance){
	40	+ double lb_0, ub_0;
	41	+ GetLengthBounds(bc, lb_0, ub_0);
	42	+
	43	+ // ほぼ直線
	44	+ if ((ub_0 - lb_0) < tolerance){
	45	+ return 0.5 * (ub_0 + lb_0);
	46	+ }
	47	+
	48	+ LengthOperator lo(tolerance);
	49	+ ONGEO_CurveSubdivision(lo, bc);
	50	+ return lo.sum_length;
	51	+}
	52	+

--- trunk/LengthParamNurbsCurve_AdaptiveQuadrate.cpp (nonexistent)

+++ trunk/LengthParamNurbsCurve_AdaptiveQuadrate.cpp (revision 52)

		@@ -0,0 +1,150 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+#include "ONGEO.h"
	4	+#include <cmath>
	5	+#include <cstdio>
	6	+#include <algorithm>
	7	+
	8	+#include "ONGEO_NumericalIntegration.hpp"
	9	+#include "ONGEO_RootFinding.hpp"
	10	+
	11	+// Reference : B. Guenter, R. Parent,
	12	+// "Computing the arc length of parametric curves".
	13	+// IEEE Computer Graphics and Applications, Vol 10, No.3, May. 1990, pp.72-78.
	14	+
	15	+namespace{
	16	+ struct IntegralFunc{
	17	+ ON_NurbsCurve nc;
	18	+ double operator()(double t) const{
	19	+ ON_3dPoint pt;
	20	+ ON_3dVector dt;
	21	+ nc.Ev1Der(t, pt, dt);
	22	+ return dt.Length();
	23	+ }
	24	+ };
	25	+}
	26	+
	27	+ONGEO_LengthParamNurbsCurve_AdaptiveQuadrate::ONGEO_LengthParamNurbsCurve_AdaptiveQuadrate(const ON_NurbsCurve &nc, double tolerance){
	28	+ IntegralFunc func;
	29	+ func.nc = nc;
	30	+ this->nc = nc;
	31	+ this->tolerance = tolerance;
	32	+
	33	+ if (!nc.IsValid()) return;
	34	+
	35	+ // Nurbs 曲線を区間分割
	36	+ {
	37	+ ON_SimpleArray<double> span_vect(nc.SpanCount()+1);
	38	+ if (span_vect.Capacity() == 1) return;
	39	+
	40	+ span_vect.SetCount(span_vect.Capacity());
	41	+ nc.GetSpanVector(span_vect.First());
	42	+
	43	+ double cpl = nc.ControlPolygonLength();
	44	+
	45	+ nc_prms.AppendNew();
	46	+ nc.GetDomain(&nc_prms[0], 0);
	47	+
	48	+ for (int j = 0; j <= nc.CVCount() - nc.Order(); ++j){
	49	+ ON_BezierCurve bc;
	50	+ if (!nc.ConvertSpanToBezier(j, bc)) continue;
	51	+ ON_Interval p_int(*nc_prms.Last(), nc.Knot(j+nc.Order()-1));
	52	+
	53	+ ON_SimpleArray<double> tess_prm;
	54	+ ONGEO_TessellateBezierCurve_QuasiInterpolating(bc, cpl / 1000.0, tess_prm);
	55	+
	56	+ for (int i = 1; i < tess_prm.Count(); ++i){
	57	+ nc_prms.Append(p_int.ParameterAt(tess_prm[i]));
	58	+ }
	59	+ }
	60	+ }
	61	+
	62	+ double prm_range = nc_prms.Last() - nc_prms.First();
	63	+
	64	+ // 適応求積法で区間ごとの曲線長を算出し、区間ごとの距離の積算値を nc_lengths に保持
	65	+ double t_prv = nc_prms[0];
	66	+ double sum_l = 0, sum_err = 0;
	67	+ nc_lengths.Append(0);
	68	+ for (int i = 1; i < nc_prms.Count(); ++i){
	69	+ double t_cur = nc_prms[i];
	70	+ double err, l = ONGEO_AdaptiveQuadrature<ONGEO_NumericalIntegration_GaussKronrod<5>, IntegralFunc>(func, t_prv, t_cur, tolerance * (t_cur - t_prv) / prm_range, err);
	71	+ sum_l += l;
	72	+ sum_err += err;
	73	+ t_prv = t_cur;
	74	+
	75	+ nc_lengths.Append(sum_l);
	76	+ }
	77	+ estimated_deviation = sum_err;
	78	+}
	79	+
	80	+double ONGEO_LengthParamNurbsCurve_AdaptiveQuadrate::Length() const{
	81	+ if (nc_lengths.Count() == 0) return 0;
	82	+ return *nc_lengths.Last();
	83	+}
	84	+
	85	+double ONGEO_LengthParamNurbsCurve_AdaptiveQuadrate::ParamToLength(double nc_prm) const{
	86	+ IntegralFunc func;
	87	+ func.nc = nc;
	88	+ double prm_range = nc_prms.Last() - nc_prms.First();
	89	+
	90	+ // 指定されたパラメータ値を超えない最も大きなパラメータ値を保持した配列から探し、
	91	+ // そのインデックスを取得する。
	92	+ const double *p = std::upper_bound(nc_prms.First(), nc_prms.Last()+1, nc_prm);
	93	+ int p_idx = p - nc_prms.First();
	94	+ if (p_idx == 0) return -1;
	95	+
	96	+ // 指定されたパラメータ値と同じ値が配列に入っている場合は、対応する距離そのままを返す。
	97	+ if (*(p-1) == nc_prm) return nc_lengths[p_idx-1];
	98	+
	99	+ // 配列から探したパラメータ値と指定されたパラメータ値で決まる区間の曲線長を算出し、足して返す。
	100	+ double err;
	101	+ double l = ONGEO_AdaptiveQuadrature<ONGEO_NumericalIntegration_GaussKronrod<5>, IntegralFunc>(func, (p-1), nc_prm, tolerance (nc_prm - *(p-1)) / prm_range, err);
	102	+
	103	+ return nc_lengths[p_idx-1] + l;
	104	+}
	105	+double ONGEO_LengthParamNurbsCurve_AdaptiveQuadrate::LengthToParam(double nc_length) const{
	106	+ IntegralFunc func;
	107	+ func.nc = nc;
	108	+ double prm_range = nc_prms.Last() - nc_prms.First();
	109	+
	110	+ const double *l = std::upper_bound(nc_lengths.First(), nc_lengths.Last()+1, nc_length);
	111	+ int l_idx = l - nc_lengths.First();
	112	+ if (l_idx == 0) return -1;
	113	+
	114	+ // 指定された距離と同じ値が配列に入っている場合は、対応するパラメータそのままを返す。
	115	+ if (*(l-1) == nc_length) return nc_prms[l_idx-1];
	116	+
	117	+
	118	+ ON_Interval p_range(nc_prms[l_idx-1], nc_prms[l_idx]);
	119	+ ON_Interval l_range(nc_lengths[l_idx-1], nc_lengths[l_idx]);
	120	+ double diff_l = nc_length - *(l-1);
	121	+#if 0
	122	+ // 配列から見つけた距離から区間を特定し、区間内の距離からパラメータをはさみうち法で算出し、返す。
	123	+ for(;;){
	124	+ double prm_estimate = p_range.ParameterAt(l_range.NormalizedParameterAt(nc_length));
	125	+ double err;
	126	+ double itg_l = ONGEO_AdaptiveQuadrature<ONGEO_NumericalIntegration_GaussKronrod<5>, IntegralFunc>(func, p_range.Min(), prm_estimate, tolerance * (prm_estimate - p_range.Min()) / prm_range, err);
	127	+ if (diff_l >= itg_l){
	128	+ if (diff_l - itg_l <= tolerance) return prm_estimate;
	129	+ l_range.m_t[0] += itg_l;
	130	+ p_range.m_t[0] = prm_estimate;
	131	+ diff_l = nc_length - l_range.m_t[0];
	132	+ }else{
	133	+ if (itg_l - diff_l < tolerance) return prm_estimate;
	134	+ l_range.m_t[1] = l_range.m_t[0] + itg_l;
	135	+ p_range.m_t[1] = prm_estimate;
	136	+ }
	137	+ }
	138	+#else
	139	+ return ONGEO_FindRootByBrentMethod(p_range.m_t[0], p_range.m_t[1],
	140	+ [&] (double t) {
	141	+ double err, itg_l = 0;
	142	+ if (std::abs(t - p_range.m_t[0]) > tolerance){
	143	+ itg_l = ONGEO_AdaptiveQuadrature<ONGEO_NumericalIntegration_GaussKronrod<5>, IntegralFunc>(func, p_range.m_t[0], t, tolerance * (t - p_range.m_t[0]) / prm_range, err);
	144	+ }
	145	+ return l_range.m_t[0] +itg_l - nc_length;
	146	+ }
	147	+ , tolerance);
	148	+#endif
	149	+}
	150	+

--- trunk/NumericalIntegration.cpp (nonexistent)

+++ trunk/NumericalIntegration.cpp (revision 52)

		@@ -0,0 +1,33 @@
	1	+#include "ONGEO.h"
	2	+#include "ONGEO_NumericalIntegration.hpp"
	3	+
	4	+// 数値積分用の係数列
	5	+
	6	+// Gauss-Kronrod 公式用の係数
	7	+// 参考
	8	+// http://www.ktech.biz/jp/archives/1009
	9	+// http://www.ktech.biz/jp/archives/1045
	10	+
	11	+template<>
	12	+ONGEO_DECL void ONGEO_GaussKronrod_Coefs<3>(const double &t_i_, const double &wg_i_, const double *&wk_i_){
	13	+ static const double t_i[] = { 0, 0.434243749346802, 0.774596669241483, 0.96049126870802 };
	14	+ static const double wg_i[] = { 0.555555555555556, 0.888888888888889, 0.555555555555556 };
	15	+ static const double wk_i[] = { 0.450916538658474, 0.401397414775962, 0.268488089868333, 0.104656226026467 };
	16	+ t_i_ = t_i, wg_i_ = wg_i, wk_i_ = wk_i;
	17	+}
	18	+
	19	+template<>
	20	+ONGEO_DECL void ONGEO_GaussKronrod_Coefs<5>(const double &t_i_, const double &wg_i_, const double *&wk_i_){
	21	+ static const double t_i[] = { 0, 0.279630413161783, 0.538469310105683, 0.754166726570849, 0.906179845938664, 0.984085360094842 };
	22	+ static const double wg_i[] = { 0.236926885056189, 0.478628670499367, 0.568888888888889, 0.478628670499367, 0.236926885056189 };
	23	+ static const double wk_i[] = { 0.282987417857491, 0.272849801912559, 0.241040339228648, 0.186800796556493, 0.115233316622474, 0.0425820367510818 };
	24	+ t_i_ = t_i, wg_i_ = wg_i, wk_i_ = wk_i;
	25	+}
	26	+
	27	+template<>
	28	+ONGEO_DECL void ONGEO_GaussKronrod_Coefs<7>(const double &t_i_, const double &wg_i_, const double *&wk_i_){
	29	+ static const double t_i[] = { 0, 0.207784955007898, 0.405845151377397, 0.586087235467691, 0.741531185599394, 0.864864423359769, 0.949107912342758, 0.991455371120813 };
	30	+ static const double wg_i[] = { 0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277, 0.129484966168870 };
	31	+ static const double wk_i[] = { 0.209482141084728, 0.204432940075299, 0.190350578064785, 0.169004726639268, 0.140653259715526, 0.104790010322250, 0.0630920926299790, 0.0229353220105292 };
	32	+ t_i_ = t_i, wg_i_ = wg_i, wk_i_ = wk_i;
	33	+}

--- trunk/ONGEO_CurveSubdivision.hpp (nonexistent)

+++ trunk/ONGEO_CurveSubdivision.hpp (revision 52)

		@@ -0,0 +1,27 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+
	4	+template <typename Operator>
	5	+void ONGEO_CurveSubdivision(Operator &op, const ON_BezierCurve &bc){
	6	+ ON_SimpleArray<double> stack;
	7	+
	8	+ stack.Append(0);
	9	+
	10	+ while(stack.Count()){
	11	+ double t_s = *stack.Last();
	12	+ stack.SetCount(stack.Count()-1);
	13	+
	14	+ ON_Interval t_int;
	15	+ t_int.m_t[0] = t_s;
	16	+ t_int.m_t[1] = (stack.Count()) ? *stack.Last() : 1.0;
	17	+
	18	+ ON_BezierCurve bc_cur = bc;
	19	+ bc_cur.Trim(t_int);
	20	+ if (!op.TestConverge(bc_cur, t_int)){
	21	+ double t_m = t_int.Mid();
	22	+ stack.Append(t_m);
	23	+ stack.Append(t_s);
	24	+ }
	25	+ }
	26	+}
	27	+

--- trunk/ONGEO_NumericalIntegration.hpp (nonexistent)

+++ trunk/ONGEO_NumericalIntegration.hpp (revision 52)

		@@ -0,0 +1,73 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+
	4	+template <int N> void ONGEO_GaussKronrod_Coefs(const double &t_i_, const double &wg_i_, const double *&wk_i_);
	5	+
	6	+template <int N> struct ONGEO_NumericalIntegration_GaussKronrod{
	7	+ double t[N2+1], wg[N], wk[N2+1];
	8	+ ONGEO_NumericalIntegration_GaussKronrod(double a, double b){
	9	+ SetRange(a, b);
	10	+ }
	11	+ void SetRange(double a, double b){
	12	+ if (a == b) return;
	13	+ const double t_i, wg_i, *wk_i;
	14	+ ONGEO_GaussKronrod_Coefs<N>(t_i, wg_i, wk_i);
	15	+
	16	+ // 区間変更
	17	+ double ba_diff = (b - a) * 0.5, ba_mean = (a + b) * 0.5;
	18	+ for (int i = 1; i <= N; ++i) t[N+i] = ba_diff * t_i[i] + ba_mean, t[N-i] = - ba_diff * t_i[i] + ba_mean;
	19	+ t[N] = ba_diff * t_i[0] + ba_mean;
	20	+ for (int i = 0; i < N; ++i) wg[i] = ba_diff * wg_i[i];
	21	+ for (int i = 1; i <= N; ++i) wk[N+i] = wk[N-i] = ba_diff * wk_i[i];
	22	+ wk[N] = ba_diff * wk_i[0];
	23	+ }
	24	+
	25	+ template <typename Func>
	26	+ double Calc(const Func &f, double &err){
	27	+ double fv[N*2+1];
	28	+ for (int i = 0; i < N*2+1; ++i) fv[i] = f(t[i]);
	29	+ double g = 0;
	30	+ for (int i = 0; i < N; ++i){
	31	+ g += fv[i2+1] wg[i];
	32	+ }
	33	+ double k = 0;
	34	+ for (int i = 0; i < N2+1; ++i) k += fv[i] wk[i];
	35	+ err = std::pow(200.0 * std::abs(k - g), 1.5);
	36	+ return k;
	37	+ }
	38	+};
	39	+
	40	+template <typename Integrator, typename Func> double ONGEO_AdaptiveQuadrature(const Func &f, double a, double b, double tolerance, double &err_result){
	41	+ double range = b - a;
	42	+ Integrator integ(0, 0);
	43	+ ON_SimpleArray<double> stack;
	44	+
	45	+ double sum = 0;
	46	+ err_result = 0;
	47	+
	48	+ stack.Append(b);
	49	+ stack.Append(a);
	50	+
	51	+ while(stack.Count() > 1){
	52	+ double t_s = *stack.Last();
	53	+ stack.SetCount(stack.Count()-1);
	54	+
	55	+ double s_a = t_s, s_b = *stack.Last();
	56	+
	57	+ integ.SetRange(s_a, s_b);
	58	+
	59	+ double sect_tolerance = tolerance * ((s_b - s_a) / range);
	60	+
	61	+ double err;
	62	+ double r = integ.Calc(f, err);
	63	+ if (err < sect_tolerance){
	64	+ sum += r;
	65	+ err_result += err;
	66	+ }else{
	67	+ double t_m = (s_a + s_b) * 0.5;
	68	+ stack.Append(t_m);
	69	+ stack.Append(t_s);
	70	+ }
	71	+ }
	72	+ return sum;
	73	+};

--- trunk/ONGEO_QuasiInterpolating.hpp (nonexistent)

+++ trunk/ONGEO_QuasiInterpolating.hpp (revision 52)

		@@ -0,0 +1,188 @@
	1	+#define NOMINMAX
	2	+#include "ONGEO.h"
	3	+#include <cmath>
	4	+#include <stack>
	5	+#include <map>
	6	+#include <algorithm>
	7	+#include <vector>
	8	+#include <string>
	9	+#include <cstdio>
	10	+
	11	+
	12	+// Reference : Yohan D.Fougerolle, Sandrine Lanquetin, Marc Neveu, and Thierry Lauthelier,
	13	+// "A Geometric Algorithm for Ray/Bezier Surfaces Intersection using Quasi-interpolating Control Net".
	14	+// Signal Image Technology and Internet Based Systems, 2008.
	15	+
	16	+
	17	+// ======
	18	+// === 曲線、曲面に共通な処理 ===
	19	+// ======
	20	+inline double ONGEO_QI_GetZInf(int dim){
	21	+ if (dim == 1) return 0;
	22	+ if (dim == 2) return 0.0625; // 1/16
	23	+ return static_cast<double>(((dim + 1) / 2) * (dim / 2))/(2.0*static_cast<double>(dim)) - 0.25;
	24	+}
	25	+
	26	+inline void ONGEO_QI_GetQuasiInterpCP(const ON_3dPoint &pp, const ON_3dPoint &pc, const ON_3dPoint &pn, ON_3dPoint &po){
	27	+ po.x = (pp.x + pc.x * 2.0 + pn.x) * 0.25;
	28	+ po.y = (pp.y + pc.y * 2.0 + pn.y) * 0.25;
	29	+ po.z = (pp.z + pc.z * 2.0 + pn.z) * 0.25;
	30	+}
	31	+
	32	+inline void ONGEO_QI_GetQuasiInterpCP(const ON_4dPoint &pp, const ON_4dPoint &pc, const ON_4dPoint &pn, ON_4dPoint &po){
	33	+ po.x = (pp.x + pc.x * 2.0 + pn.x) * 0.25;
	34	+ po.y = (pp.y + pc.y * 2.0 + pn.y) * 0.25;
	35	+ po.z = (pp.z + pc.z * 2.0 + pn.z) * 0.25;
	36	+ po.w = (pp.w + pc.w * 2.0 + pn.w) * 0.25;
	37	+}
	38	+
	39	+inline double ONGEO_QI_GetDelta2LengthSquared(const ON_3dPoint &pp, const ON_3dPoint &pc, const ON_3dPoint &pn){
	40	+ double d, dd = 0;
	41	+ d = pp.x - pc.x * 2.0 + pn.x, d *= d, dd += d;
	42	+ d = pp.y - pc.y * 2.0 + pn.y, d *= d, dd += d;
	43	+ d = pp.z - pc.z * 2.0 + pn.z, d *= d, dd += d;
	44	+ return dd;
	45	+// return ON_3dVector(pp - pc * 2.0 + pn).LengthSquared();
	46	+}
	47	+
	48	+inline double ONGEO_QI_GetDelta2LengthSquared(const ON_4dPoint &pp, const ON_4dPoint &pc, const ON_4dPoint &pn){
	49	+ double d, dd = 0;
	50	+ d = pp.x - pc.x * 2.0 + pn.x, d *= d, dd += d;
	51	+ d = pp.y - pc.y * 2.0 + pn.y, d *= d, dd += d;
	52	+ d = pp.z - pc.z * 2.0 + pn.z, d *= d, dd += d;
	53	+ d = pp.w - pc.w * 2.0 + pn.w, d *= d, dd += d;
	54	+ return dd;
	55	+}
	56	+
	57	+template <typename T, typename S> struct ONGEO_QI_PointDimTraits{
	58	+};
	59	+// ======
	60	+// === 曲線用の処理 ===
	61	+// ======
	62	+template <> struct ONGEO_QI_PointDimTraits<ON_3dPoint, ON_BezierCurve>{
	63	+ static const int DIM = 3;
	64	+ static double GetErrFromErr2(const ON_BezierCurve &bez, double errmax2){
	65	+ return ONGEO_QI_GetZInf(bez.m_order-1) * std::sqrt(errmax2);
	66	+ }
	67	+};
	68	+
	69	+template <> struct ONGEO_QI_PointDimTraits<ON_4dPoint, ON_BezierCurve>{
	70	+ static const int DIM = 4;
	71	+ static double GetErrFromErr2(const ON_BezierCurve &bez, double errmax2){
	72	+ double wmin = bez.Weight(0);
	73	+ for (int i = 0; i < bez.m_order; ++i){
	74	+ if (wmin > bez.Weight(i)) wmin = bez.Weight(i);
	75	+ }
	76	+ return ONGEO_QI_PointDimTraits<ON_3dPoint, ON_BezierCurve>::GetErrFromErr2(bez, errmax2) / wmin;
	77	+ }
	78	+};
	79	+
	80	+template <typename T> void ONGEO_QI_GetQuasiInterpControlPoints(const ON_BezierCurve &bez, ON_SimpleArray<T> &qicp){
	81	+ ON_SimpleArray<T> work_cp;
	82	+ work_cp.SetCapacity(bez.m_order);
	83	+ work_cp.SetCount(work_cp.Capacity());
	84	+
	85	+ qicp.SetCapacity(work_cp.Count());
	86	+ qicp.SetCount(qicp.Capacity());
	87	+ for (int i = 0; i < bez.m_order; ++i){
	88	+ bez.GetCV(i, j, work_cp[i]);
	89	+ }
	90	+ for (int i = 1; i < bez.m_order-1; ++i){
	91	+ ONGEO_QI_GetQuasiInterpCP(work_cp[i-1], work_cp[i], work_cp[i+1], qicp[i]);
	92	+ }
	93	+ qicp[0] = work_cp[0];
	94	+ qicp[bez.m_order-1] = work_cp[bez.m_order-1];
	95	+}
	96	+
	97	+template <typename T> double ONGEO_QI_GetError(const ON_BezierCurve &bez){
	98	+ double errmax2 = 0;
	99	+ const int ord = bez.m_order;
	100	+ ON_SimpleArray<T> cp(ord);
	101	+ for (int i = 0; i < ord; ++i){
	102	+ ONGEO_GetCVHomogeneous(bez, i, cp[i]);
	103	+ }
	104	+ for (int i = 1; i < ord-1; ++i){
	105	+ double d = ONGEO_QI_GetDelta2LengthSquared(cp[i-1], cp[i], cp[i+1]);
	106	+ if (errmax2 < d) errmax2 = d;
	107	+ }
	108	+
	109	+ return ONGEO_QI_PointDimTraits<T, ON_BezierCurve>::GetErrFromErr2(bez, errmax2);
	110	+}
	111	+
	112	+// ======
	113	+// === 曲面用の処理 ===
	114	+// ======
	115	+template <> struct ONGEO_QI_PointDimTraits<ON_3dPoint, ON_BezierSurface>{
	116	+ static const int DIM = 3;
	117	+ static double GetErrFromErr2(const ON_BezierSurface &bez, double errmax2u, double errmax2v){
	118	+ return ONGEO_QI_GetZInf(bez.m_order[0]-1) * std::sqrt(errmax2u) + ONGEO_QI_GetZInf(bez.m_order[1]-1) * std::sqrt(errmax2v);
	119	+ }
	120	+};
	121	+
	122	+template <> struct ONGEO_QI_PointDimTraits<ON_4dPoint, ON_BezierSurface>{
	123	+ static const int DIM = 4;
	124	+ static double GetErrFromErr2(const ON_BezierSurface &bez, double errmax2u, double errmax2v){
	125	+ double wmin = bez.Weight(0, 0);
	126	+ for (int j = 0; j < bez.m_order[1]; ++j){
	127	+ for (int i = 0; i < bez.m_order[0]; ++i){
	128	+ if (wmin > bez.Weight(i, j)) wmin = bez.Weight(i, j);
	129	+ }
	130	+ }
	131	+ return ONGEO_QI_PointDimTraits<ON_3dPoint, ON_BezierSurface>::GetErrFromErr2(bez, errmax2u, errmax2v) / wmin;
	132	+ }
	133	+};
	134	+
	135	+template <typename T> void ONGEO_QI_GetQuasiInterpControlPoints(const ON_BezierSurface &bez, ON_SimpleArray<T> &qicp){
	136	+ ON_SimpleArray<T> work_cp;
	137	+ work_cp.SetCapacity(bez.m_order[0] * bez.m_order[1]);
	138	+ work_cp.SetCount(work_cp.Capacity());
	139	+ // U方向
	140	+ ON_SimpleArray<T> work_u(bez.m_order[0]);
	141	+ for (int j = 0, ji = 0; j < bez.m_order[1]; ++j, ji += bez.m_order[0]){
	142	+ for (int i = 0; i < bez.m_order[0]; ++i){
	143	+ bez.GetCV(i, j, work_u[i]);
	144	+ }
	145	+ for (int i = 1; i < bez.m_order[0]-1; ++i){
	146	+ ONGEO_QI_GetQuasiInterpCP(work_u[i-1], work_u[i], work_u[i+1], work_cp[ji+i]);
	147	+ }
	148	+ work_cp[ji] = work_u[0], work_cp[ji+bez.m_order[0]-1] = work_u[bez.m_order[0]-1];
	149	+ }
	150	+// ON_SimpleArray<T> work_v(bez.m_order[1]);
	151	+
	152	+ qicp.SetCapacity(work_cp.Count());
	153	+ qicp.SetCount(qicp.Capacity());
	154	+ // V方向
	155	+ for (int i = 0; i < bez.m_order[0]; ++i){
	156	+ int j, ji;
	157	+ for (j = 1, ji = bez.m_order[0]; j < bez.m_order[1]-1; ++j, ji += bez.m_order[0]){
	158	+ ONGEO_QI_GetQuasiInterpCP(work_cp[ji+i-bez.m_order[0]], work_cp[ji+i], work_cp[ji+i+bez.m_order[0]], qicp[ji+i]);
	159	+ }
	160	+ qicp[i] = work_cp[i], qicp[ji+i] = work_cp[ji+i];
	161	+ }
	162	+}
	163	+
	164	+template <typename T> double ONGEO_QI_GetError(const ON_BezierSurface &bez){
	165	+ // U方向
	166	+ double errmax2u = 0;
	167	+ const int ord_u = bez.m_order[0], ord_v = bez.m_order[1];
	168	+ ON_SimpleArray<T> cp(ord_u * ord_v);
	169	+ for (int j = 0, ji = 0; j < ord_v; ++j, ji += ord_u){
	170	+ for (int i = 0; i < ord_u; ++i){
	171	+ ONGEO_GetCVHomogeneous(bez, i, j, cp[ji+i]);
	172	+ }
	173	+ for (int i = 1; i < ord_u-1; ++i){
	174	+ double d = ONGEO_QI_GetDelta2LengthSquared(cp[ji+i-1], cp[ji+i], cp[ji+i+1]);
	175	+ if (errmax2u < d) errmax2u = d;
	176	+ }
	177	+ }
	178	+ // V方向
	179	+ double errmax2v = 0;
	180	+ for (int i = 0; i < ord_u; ++i){
	181	+ for (int j = 1, ji = ord_u; j < ord_v-1; ++j, ji += ord_u){
	182	+ double d = ONGEO_QI_GetDelta2LengthSquared(cp[ji-ord_u], cp[ji], cp[ji+ord_u]);
	183	+ if (errmax2v < d) errmax2v = d;
	184	+ }
	185	+ }
	186	+
	187	+ return ONGEO_QI_PointDimTraits<T, ON_BezierSurface>::GetErrFromErr2(bez, errmax2u, errmax2v);
	188	+}

--- trunk/ONGEO_RootFinding.hpp (nonexistent)

+++ trunk/ONGEO_RootFinding.hpp (revision 52)

		@@ -0,0 +1,66 @@
	1	+#define NOMINMAX
	2	+#include "ONGEO.h"
	3	+#include <cmath>
	4	+
	5	+// Brent法を用いて根を求める。区間 ti1 - ti2 の中で単調増加、または単調減少であること。
	6	+// @param ti1 [in] 狭めたい区間の下限
	7	+// @param ti2 [in] 狭めたい区間の上限 ti1 < ti2であること。
	8	+// @return 根となる値
	9	+// Reference : R.P.Brent,
	10	+// "An algorithm with guaranteed convergence for finding a zero of a function".
	11	+// The Computer Journal, 14(1971), pp.422-425.
	12	+template <typename EvalFunc>
	13	+double ONGEO_FindRootByBrentMethod(double ti1, double ti2, EvalFunc f, double tolerance) {
	14	+ double a = ti1, b = ti2;
	15	+ double fa = f(a);
	16	+ double fb = f(b);
	17	+ double c = a, fc = fa, d, dprev;
	18	+ d = dprev = b - a;
	19	+ // d : 収束処理でbを動かす量
	20	+ for(;;){
	21	+ // 根に近い方をb,fb、遠い方をa=c,fa=fcとする。
	22	+ if (std::abs(fc) < std::abs(fb)){
	23	+ a = b, b = c, c = a;
	24	+ fa = fb, fb = fc, fc = fa;
	25	+ }
	26	+ double tol = 2.0 * DBL_EPSILON * std::abs(b) + tolerance;
	27	+ double c_b = c - b;
	28	+ double m = 0.5 * c_b;
	29	+ if (std::abs(m) <= tol \|\| fb == 0) break;
	30	+
	31	+ // bisectionの場合のd、dprev
	32	+ d = dprev = m;
	33	+ if (std::abs(dprev) >= tol && std::abs(fa) > std::abs(fb)){
	34	+ double r3 = fb / fa;
	35	+ double p, q;
	36	+ if (a == c){
	37	+ // linear interpolation
	38	+ p = c_b * r3;
	39	+ q = 1.0 - r3;
	40	+ }else{
	41	+ // inverse quadratic interpolation
	42	+ double ifc = 1.0 / fc;
	43	+ double r1 = fa * ifc, r2 = fb * ifc;
	44	+ p = r3 * (c_b * r1 * (r1 - r2) - (b - a) * (r2 - 1));
	45	+ q = (r1 - 1) * (r2 - 1) * (r3 - 1);
	46	+ }
	47	+ if (p > 0) q = -q;
	48	+ else p = -p;
	49	+ if (2.0 * p < 1.5 * c_b * q - std::abs(tol * q) && p < std::abs(0.5 * dprev * q)){
	50	+ // 上記条件に合う場合は、dとdprevをlinear interpolation、
	51	+ // またはinverse quadratic interpolation由来のものに書き換える。
	52	+ dprev = d;
	53	+ d = p / q;
	54	+ }
	55	+ }
	56	+ // a, fa, b, fbを更新
	57	+ a = b, fa = fb;
	58	+ b += (std::abs(d) > tol) ? d : ((m > 0) ? tol : -tol);
	59	+ fb = f(b);
	60	+ if (fb * fc > 0){
	61	+ c = a, fc = fa;
	62	+ dprev = b - a;
	63	+ }
	64	+ }
	65	+ return b;
	66	+}

--- trunk/TessellateCurve.cpp (nonexistent)

+++ trunk/TessellateCurve.cpp (revision 52)

		@@ -0,0 +1,108 @@
	1	+#define NOMINMAX
	2	+#include "opennurbs.h"
	3	+#include "ONGEO.h"
	4	+#include <cmath>
	5	+#include <cstdio>
	6	+
	7	+#include "ONGEO_CurveSubdivision.hpp"
	8	+#include "ONGEO_QuasiInterpolating.hpp"
	9	+
	10	+inline double GetRoughSag(const ON_BezierCurve &bc_cur){
	11	+ ON_Line l;
	12	+ bc_cur.GetCV(0, l.from);
	13	+ bc_cur.GetCV(bc_cur.CVCount()-1, l.to);
	14	+ double sag_max = 0;
	15	+ for (int i = 1; i < bc_cur.CVCount() - 1; ++i){
	16	+ ON_3dPoint pt;
	17	+ bc_cur.GetCV(i, pt);
	18	+ double dist = l.MinimumDistanceTo(pt);
	19	+ if (sag_max < dist) sag_max = dist;
	20	+ }
	21	+ return sag_max;
	22	+}
	23	+
	24	+struct SimpleTessellateOperator{
	25	+ // input
	26	+ double tolerance;
	27	+
	28	+ // output
	29	+ ON_SimpleArray<double> prms;
	30	+ SimpleTessellateOperator(double tolerance_) : tolerance(tolerance_){}
	31	+ void Initialize(double tolerance_){
	32	+ tolerance = tolerance_;
	33	+ prms.Empty();
	34	+ }
	35	+ bool TestConverge(const ON_BezierCurve &bc_cur, const ON_Interval &t_int){
	36	+ double rsag = GetRoughSag(bc_cur);
	37	+ bool converged = rsag < tolerance;
	38	+ if (converged){
	39	+ prms.Append(t_int.m_t[0]);
	40	+ }
	41	+ return converged;
	42	+ }
	43	+ void EndConverge(){
	44	+ prms.Append(1.0);
	45	+ }
	46	+};
	47	+
	48	+void ONGEO_TessellateBezierCurve_Simple(const ON_BezierCurve &bc, double tolerance, ON_SimpleArray<double> &prms){
	49	+ SimpleTessellateOperator to(tolerance);
	50	+ ONGEO_CurveSubdivision(to, bc);
	51	+ to.EndConverge();
	52	+
	53	+ ONGEO_Swap(prms, to.prms);
	54	+}
	55	+
	56	+template<typename T> struct QuasiInterpolatingTessellateOperator{
	57	+ // input
	58	+ double tolerance;
	59	+
	60	+ // output
	61	+ ON_SimpleArray<double> prms;
	62	+ QuasiInterpolatingTessellateOperator(double tolerance_) : tolerance(tolerance_){}
	63	+ void Initialize(double tolerance_){
	64	+ tolerance = tolerance_;
	65	+ prms.Empty();
	66	+ }
	67	+ bool TestConverge(const ON_BezierCurve &bc_cur, const ON_Interval &t_int){
	68	+ double qi_sag = ONGEO_QI_GetError<T>(bc_cur);
	69	+ bool qi_converged = qi_sag < tolerance;
	70	+ if (qi_converged){
	71	+#if 1
	72	+ int idx_c = bc_cur.CVCount() / 2;
	73	+ double t0 = static_cast<double>(idx_c) / (bc_cur.CVCount()-1);
	74	+ double t1 = static_cast<double>(idx_c+1) / (bc_cur.CVCount()-1);
	75	+ ON_3dPoint ptcrv = bc_cur.PointAt((t0+t1)*0.5);
	76	+ ON_3dPoint ptlin = (bc_cur.PointAt(t0) + bc_cur.PointAt(t1)) * 0.5;
	77	+ if (ptcrv.DistanceTo(ptlin) > tolerance) return false;
	78	+#endif
	79	+
	80	+ for (int i = 0; i < bc_cur.CVCount() - 1; ++i){
	81	+ double t = static_cast<double>(i) / (bc_cur.CVCount() - 1);
	82	+ prms.Append(t_int.ParameterAt(t));
	83	+ }
	84	+ return true;
	85	+ }
	86	+ return false;
	87	+ }
	88	+ void EndConverge(){
	89	+ prms.Append(1.0);
	90	+ }
	91	+};
	92	+
	93	+void ONGEO_TessellateBezierCurve_QuasiInterpolating(const ON_BezierCurve &bc, double tolerance, ON_SimpleArray<double> &prms){
	94	+ if (bc.IsRational()){
	95	+ QuasiInterpolatingTessellateOperator<ON_4dPoint> qo(tolerance);
	96	+ ONGEO_CurveSubdivision(qo, bc);
	97	+ qo.EndConverge();
	98	+
	99	+ ONGEO_Swap(prms, qo.prms);
	100	+ }else{
	101	+ QuasiInterpolatingTessellateOperator<ON_3dPoint> qo(tolerance);
	102	+ ONGEO_CurveSubdivision(qo, bc);
	103	+ qo.EndConverge();
	104	+
	105	+ ONGEO_Swap(prms, qo.prms);
	106	+ }
	107	+}
	108	+

Show on old repository browser