Revisão | 59 (tree) |
---|---|
Hora | 2020-01-02 00:46:03 |
Autor | mocchi_2012 |
・最小二乗法で点列からNurbs曲線を作成する関数の追加
・ニュートン法によるNurbs曲線の最近点を求める関数の追加
@@ -1,6 +1,6 @@ | ||
1 | 1 | /* |
2 | 2 | * FitNurbsToPoints |
3 | - * Copylight (C) 2019 mocchi | |
3 | + * Copylight (C) 2019,2020 mocchi | |
4 | 4 | * mocchi_2003@yahoo.co.jp |
5 | 5 | * License: Boost ver.1 |
6 | 6 | */ |
@@ -8,6 +8,7 @@ | ||
8 | 8 | #include "ONGEO.h" |
9 | 9 | #include <cmath> |
10 | 10 | #include <numeric> |
11 | +#include <algorithm> | |
11 | 12 | #ifdef max |
12 | 13 | #undef max |
13 | 14 | #endif |
@@ -93,7 +94,7 @@ | ||
93 | 94 | } |
94 | 95 | }; |
95 | 96 | |
96 | -void InitializeSolver(int pt_cnt, int order, double *prm, double *knot, int knot_count, IEquationSolver *solver){ | |
97 | +void InitializeSolver_ExactlyPassThroughPoints(int pt_cnt, int order, double *prm, double *knot, int knot_count, IEquationSolver *solver){ | |
97 | 98 | int n = pt_cnt - 1; |
98 | 99 | int m = knot_count; |
99 | 100 | int p = order - 1; |
@@ -108,32 +109,91 @@ | ||
108 | 109 | } |
109 | 110 | knot[j+p-1] = s / static_cast<double>(p); |
110 | 111 | } |
111 | - ON_Matrix N(n+1, n+1); | |
112 | - for (int j = 0; j < n + 1; ++j){ | |
113 | - for (int i = 0; i < n + 1; ++i){ | |
114 | - N[j][i] = 0; | |
112 | + if (order > 2){ | |
113 | + ON_Matrix N(n+1, n+1); | |
114 | + for (int j = 0; j < n + 1; ++j){ | |
115 | + for (int i = 0; i < n + 1; ++i){ | |
116 | + N[j][i] = 0; | |
117 | + } | |
115 | 118 | } |
119 | + N[0][0] = 1; N[n][n] = 1; | |
120 | + | |
121 | + ON_SimpleArray<double> basis(order * order); basis.SetCount(basis.Capacity()); | |
122 | + for (int j = 1; j < n; ++j){ | |
123 | + double t = prm[j]; | |
124 | + int span_index = ON_NurbsSpanIndex(order, n+1, knot, t, 0, 0); | |
125 | + for (int i = 0; i < span_index; ++i){ | |
126 | + N[j][i] = 0; | |
127 | + } | |
128 | + ON_EvaluateNurbsBasis(order, knot + span_index, t, basis.First()); | |
129 | + for (int i = span_index; i < span_index + order; ++i){ | |
130 | + N[j][i] = basis[i-span_index]; | |
131 | + } | |
132 | + for (int i = span_index + order; i <= n; ++i){ | |
133 | + N[j][i] = 0; | |
134 | + } | |
135 | + } | |
136 | + | |
137 | + solver->Initialize(N); | |
116 | 138 | } |
117 | - N[0][0] = 1; N[n][n] = 1; | |
139 | +} | |
118 | 140 | |
141 | +// R のサイズは cpt_cnt - 2 | |
142 | +void InitializeSolver_PassThroughPoints_LSQ_EndPointsConstraint(int pt_cnt, int order, int cpt_cnt, double *prm, double *knot, int knot_count, const ON_3dPoint *pts, ON_3dPoint *R, IEquationSolver *solver){ | |
143 | + int n = cpt_cnt - 1; | |
144 | + int m = pt_cnt - 1; | |
145 | + int p = order - 1; | |
146 | + for (int i = 0; i < p; ++i){ | |
147 | + knot[i] = 0; | |
148 | + knot[knot_count-i-1] = 1; | |
149 | + } | |
150 | + | |
151 | + double d = static_cast<double>(m + 1) / static_cast<double>(n - p + 1); | |
152 | + for (double j = 1.0; j <= n - p; j += 1.0){ | |
153 | + double i = std::floor(j * d); | |
154 | + int ii = static_cast<int>(i); | |
155 | + double a = j * d - i; | |
156 | + knot[static_cast<int>(j)+p-1] = (1.0 - a) * prm[ii-1] + a * prm[ii]; | |
157 | + } | |
158 | + | |
159 | + ON_Matrix N(m-1, n-1); | |
119 | 160 | ON_SimpleArray<double> basis(order * order); basis.SetCount(basis.Capacity()); |
120 | - for (int j = 1; j < n; ++j){ | |
121 | - double t = prm[j]; | |
161 | + ON_ClassArray<ON_3dPoint> rpt(m-1); | |
162 | + | |
163 | + for (int k = 1; k <= m-1; ++k){ | |
164 | + int r = k; | |
165 | + double t = prm[k]; | |
166 | + | |
122 | 167 | int span_index = ON_NurbsSpanIndex(order, n+1, knot, t, 0, 0); |
123 | - for (int i = 0; i < span_index; ++i){ | |
124 | - N[j][i] = 0; | |
125 | - } | |
126 | 168 | ON_EvaluateNurbsBasis(order, knot + span_index, t, basis.First()); |
127 | - for (int i = span_index; i < span_index + order; ++i){ | |
128 | - N[j][i] = basis[i-span_index]; | |
169 | + double basis_0 = 0, basis_n = 0; | |
170 | + for (int c = 0; c <= n; ++c){ | |
171 | + double &dest = (c == 0) ? basis_0 : ((c == n) ? basis_n : N[r-1][c-1]); | |
172 | + int basis_idx = c - span_index; | |
173 | + if (basis_idx < 0 || basis_idx > order){ | |
174 | + dest = 0; | |
175 | + }else{ | |
176 | + dest = basis[basis_idx]; | |
177 | + } | |
129 | 178 | } |
130 | - for (int i = span_index + order; i <= n; ++i){ | |
131 | - N[j][i] = 0; | |
179 | + rpt.Append(pts[k] - pts[0] * basis_0 - pts[m] * basis_n); | |
180 | + } | |
181 | + | |
182 | + for (int c = 1; c <= n-1; ++c){ | |
183 | + ON_3dPoint Rc(0,0,0); | |
184 | + for (int r = 1; r <= m-1; ++r){ | |
185 | + Rc += rpt[r-1] * N[r-1][c-1]; | |
132 | 186 | } |
187 | + R[c-1] = Rc; | |
133 | 188 | } |
134 | 189 | |
135 | - solver->Initialize(N); | |
190 | + ON_Matrix NT = N; | |
191 | + NT.Transpose(); | |
192 | + ON_Matrix NN; | |
193 | + NN.Multiply(NT, N); | |
194 | + solver->Initialize(NN); | |
136 | 195 | } |
196 | + | |
137 | 197 | } |
138 | 198 | |
139 | 199 | bool ONGEO_FitNurbsCurveToPointArray(const ON_3dPoint *pts, int pt_cnt, int order, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, ONGEO_NI_Solver solver_id){ |
@@ -151,11 +211,14 @@ | ||
151 | 211 | ON_SimpleArray<double> prm(n+1); prm.SetCount(n+1); |
152 | 212 | CreateParameter(pts, pt_cnt, 1, pmethod, prm); |
153 | 213 | |
214 | + InitializeSolver_ExactlyPassThroughPoints(pt_cnt, order, prm.First(), nc.m_knot, nc.KnotCount(), solver); | |
154 | 215 | |
155 | - InitializeSolver(pt_cnt, order, prm.First(), nc.m_knot, nc.KnotCount(), solver); | |
156 | - | |
157 | 216 | ON_ClassArray<ON_3dPoint> Q(n+1); Q.SetCount(n+1); |
158 | - solver->Solve(n+1, pts, Q.First(), 1); | |
217 | + if (order == 2){ | |
218 | + for (int i = 0; i < pt_cnt; ++i) Q[i] = pts[i]; | |
219 | + }else{ | |
220 | + solver->Solve(n+1, pts, Q.First(), 1); | |
221 | + } | |
159 | 222 | for (int i = 0; i <= n; ++i){ |
160 | 223 | nc.SetCV(i, Q[i]); |
161 | 224 | } |
@@ -163,6 +226,217 @@ | ||
163 | 226 | return true; |
164 | 227 | } |
165 | 228 | |
229 | +double ONGEO_FitNurbsCurveToPointArray_LeastSquare(const ON_3dPoint *pts, int pt_cnt, int order, int cpt_cnt, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, double *prm_, double *err_, ONGEO_NI_Solver solver_id){ | |
230 | + Solver_ON_Matrix_Inverse solver_omi; | |
231 | + IEquationSolver *solver = 0; | |
232 | + if (solver_id == ONGEO_NI_Solver_ON_Matrix_Inverse){ | |
233 | + solver = &solver_omi; | |
234 | + } | |
235 | + | |
236 | + ON_NurbsCurve pol; | |
237 | + ONGEO_FitNurbsCurveToPointArray(pts, pt_cnt, 2, pmethod, pol, solver_id); | |
238 | + | |
239 | + int n = pt_cnt - 1; | |
240 | + ON_ClassArray<double> prm_buf; | |
241 | + double *prm = prm_; | |
242 | + if (!prm){ | |
243 | + prm_buf.SetCapacity(pt_cnt); | |
244 | + prm_buf.SetCount(pt_cnt); | |
245 | + prm = prm_buf.First(); | |
246 | + } | |
247 | + for (int i = 0; i < n+1; ++i){ | |
248 | + prm[i] = pol.Knot(i); | |
249 | + } | |
250 | + | |
251 | + nc.Create(3, false, order, cpt_cnt); | |
252 | + ON_ClassArray<ON_3dPoint> R(nc.CVCount()-2); R.SetCount(R.Capacity()); | |
253 | + InitializeSolver_PassThroughPoints_LSQ_EndPointsConstraint(pt_cnt, order, nc.CVCount(), prm, nc.m_knot, nc.KnotCount(), pts, R, solver); | |
254 | + ON_ClassArray<ON_3dPoint> P(nc.CVCount()); P.SetCount(P.Capacity()); | |
255 | + P[0] = pts[0]; | |
256 | + P[P.Count()-1] = pts[pt_cnt-1]; | |
257 | + solver->Solve(R.Count(), R.First(), P.First()+1, 1); | |
258 | + for (int i = 0; i < nc.CVCount(); ++i) nc.SetCV(i, P[i]); | |
259 | + | |
260 | + double err_max = 0; | |
261 | + for (int i = 1; i < pt_cnt - 1; ++i){ | |
262 | + double err = -1; | |
263 | + double t = prm[i]; | |
264 | + if (ONGEO_NearestPointNurbsCurve_Newton(nc, pts[i], t)){ | |
265 | + if (prm_) prm_[i] = t; | |
266 | + err = pts[i].DistanceTo(nc.PointAt(t)); | |
267 | + if (err_) err_[i] = err; | |
268 | + } | |
269 | + if (err_max < err) err_max = err; | |
270 | + } | |
271 | + return err_max; | |
272 | +} | |
273 | + | |
274 | +#if 0 | |
275 | +double CalculateBr(const ON_NurbsCurve &nc, int r, int s){ | |
276 | + int ord = nc.Order(); | |
277 | + int p = ord - 1; | |
278 | + int last = r - s; | |
279 | + int first = r - p; | |
280 | + int off = first - 1; | |
281 | + const double *knot = nc.Knot(); | |
282 | + | |
283 | + ON_ClassArray<ON_3dPoint> temp(last + 2 - off); | |
284 | + temp.SetCount(temp.Capacity()); | |
285 | + nc.GetCV(off, temp[0]); | |
286 | + nc.GetCV(last+1, temp[last + 1 - off]); | |
287 | + int i = first, j = last; | |
288 | + int ii = 1, jj = last - off; | |
289 | + while(j-i > 0){ | |
290 | + double alfi = (knot[r] - knot[i]) / (knot[i+ord] - knot[i]); | |
291 | + double alfj = (knot[r] - knot[j]) / (knot[j+ord] - knot[j]); | |
292 | + ON_3dPoint Pi, Pj; | |
293 | + nc.GetCV(i, Pi); nc.GetCV(j, Pj); | |
294 | + temp[ii] = (Pi - temp[ii-1] * (1.0 - alfi)) / alfi; | |
295 | + temp[jj] = (Pj - temp[jj+1] * alfj) / (1.0 - alfj); | |
296 | + i = i + 1, ii = ii + 1; | |
297 | + j = j - 1, jj = jj - 1; | |
298 | + } | |
299 | + if (j - 1 < 0){ | |
300 | + return temp[ii-1].DistanceTo(temp[jj+1]); | |
301 | + }else{ | |
302 | + double alfi = (knot[r] - knot[i]) / (knot[i+ord] - knot[i]); | |
303 | + ON_3dPoint Pi; | |
304 | + nc.GetCV(i, Pi); | |
305 | + return Pi.DistanceTo(temp[ii+1] * alfi + (1.0 - alfi) * temp[ii-1]); | |
306 | + } | |
307 | +} | |
308 | + | |
309 | +bool ONGEO_FitNurbsCurveToPointArray(const ON_3dPoint *pts, int pt_cnt, int order, double tolerance, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, ONGEO_NI_Solver solver_id){ | |
310 | + if (tolerance <= 0) return ONGEO_FitNurbsCurveToPointArray(pts, pt_cnt, order, pmethod, nc, solver_id); | |
311 | + | |
312 | + Solver_ON_Matrix_Inverse solver_omi; | |
313 | + IEquationSolver *solver = 0; | |
314 | + if (solver_id == ONGEO_NI_Solver_ON_Matrix_Inverse){ | |
315 | + solver = &solver_omi; | |
316 | + } | |
317 | + | |
318 | + ON_SimpleArray<double> basis(order * order); basis.SetCount(basis.Capacity()); | |
319 | + | |
320 | + ON_NurbsCurve nc_tmp; | |
321 | + ONGEO_FitNurbsCurveToPointArray(pts, pt_cnt, 2, pmethod, nc_tmp, solver_id); | |
322 | + | |
323 | + int n = pt_cnt - 1; | |
324 | + ON_SimpleArray<double> prm(n+1); prm.SetCount(n+1); | |
325 | + ON_SimpleArray<double> err(n+1); err.SetCount(n+1); | |
326 | + for (int i = 0; i < n+1; ++i){ | |
327 | + prm[i] = nc_tmp.Knot(i); | |
328 | + err[i] = 0; | |
329 | + } | |
330 | + | |
331 | + int p = order - 1; | |
332 | + // ノット除去 | |
333 | + for (int deg = 1; deg <= p; ++deg){ | |
334 | + ON_NurbsCurve nc_tmp2; | |
335 | + nc_tmp2.Create(3, false, deg+1, nc_tmp.CVCount()); | |
336 | + // ノット除去 | |
337 | + while(true){ | |
338 | + // 最小の Br 値と、そのBrに対応するノット番号を取得 | |
339 | + double Br_min = -1, err_min = -1; | |
340 | + const double *knot = nc_tmp.Knot(); | |
341 | + int s = nc_tmp.KnotMultiplicity(0); | |
342 | + int r_min = -1, r = s-1; | |
343 | + s = nc_tmp.KnotMultiplicity(r+1); | |
344 | + r += s; | |
345 | + while(s > 0){ | |
346 | + double Br = CalculateBr(nc_tmp, r, s); | |
347 | + | |
348 | + int p_s = deg + s; | |
349 | + s = nc_tmp.KnotMultiplicity(r+1); | |
350 | + if (r + s >= nc_tmp.KnotCount()-1) break; | |
351 | +#if 0 | |
352 | + nc_tmp2 = nc_tmp; | |
353 | + nc_tmp2.RemoveSpan(r-(p-1)); | |
354 | + double Br; | |
355 | + if ((p_s & 1) == 0){ | |
356 | + int k = p_s / 2; | |
357 | + double a = (knot[r] - knot[r-k]) / (knot[r-k+p+1] - knot[r-k]); | |
358 | + ON_3dPoint P, P1r, P1l; | |
359 | + nc_tmp.GetCV(r-k, P); | |
360 | + nc_tmp2.GetCV(r-k+1, P1r); | |
361 | + int rk_1 = r-k-1; | |
362 | + nc_tmp2.GetCV(rk_1 >= 0 ? rk_1 : 0, P1l); | |
363 | + Br = P.DistanceTo(P1r * a + P1l * (1.0-a)); | |
364 | + }else{ | |
365 | + int k = (p_s + 1) / 2; | |
366 | + ON_3dPoint P1r, P1l; | |
367 | + nc_tmp2.GetCV(r-k+1, P1r); | |
368 | + nc_tmp2.GetCV(r-k, P1l); | |
369 | + Br = P1r.DistanceTo(P1l); | |
370 | + } | |
371 | +#endif | |
372 | + if (Br_min < 0 || Br_min > Br) Br_min = Br, r_min = r; | |
373 | + r += s; | |
374 | + } | |
375 | + if (Br_min == -1) break; | |
376 | + { | |
377 | + int span_index = ON_NurbsSpanIndex(order, nc_tmp.CVCount(), knot, knot[r_min], 0, 0); | |
378 | + ON_EvaluateNurbsBasis(order, knot + span_index, knot[r_min], basis.First()); | |
379 | + int s_min = nc_tmp.KnotMultiplicity(r_min); | |
380 | + int p_s = deg + s_min; | |
381 | + if ((p_s & 1) == 0){ | |
382 | + int k = p_s / 2; | |
383 | + err_min = basis[r_min - k - span_index+1] * Br_min; | |
384 | + }else{ | |
385 | + int k = (p_s + 1) / 2; | |
386 | + double a = (knot[r_min] - knot[r_min-k+1])/(knot[r_min-k+deg+2] - knot[r_min-k+1]); | |
387 | + err_min = (1.0 - a) * basis[r_min - k + 1 - span_index+1] * Br_min; | |
388 | + } | |
389 | + } | |
390 | + if (err_min > tolerance) break; | |
391 | + int idx_s = ON_NurbsSpanIndex(order, nc_tmp.CVCount(), knot, knot[r_min], 0, 0); | |
392 | + int idx_n = idx_s + 2 * deg - 1; | |
393 | + | |
394 | + ptrdiff_t diff_s = std::lower_bound(prm.First(), prm.Last()+1, knot[idx_s]) - prm.First(); | |
395 | + ptrdiff_t diff_e = std::upper_bound(prm.First(), prm.Last()+1, knot[idx_n]) - prm.First(); | |
396 | + double *es = err.First() + diff_s, *ee = err.First() + diff_e; | |
397 | + | |
398 | + bool knot_removable = true; | |
399 | + for (double *e = es; e != ee; ++e){ | |
400 | + if (*e + err_min < tolerance) continue; | |
401 | + knot_removable = false; | |
402 | + break; | |
403 | + } | |
404 | + if (knot_removable){ | |
405 | + for (double *e = es; e != ee; ++e) *e += err_min; | |
406 | + nc_tmp.RemoveSpan(r_min-(order-2)); | |
407 | + }else break; | |
408 | + } | |
409 | + if (deg == p){ | |
410 | + nc = nc_tmp; | |
411 | + break; | |
412 | + } | |
413 | + | |
414 | + // 次数を1上げて最小二乗法で曲線を再生成 | |
415 | + nc_tmp.Create(3, false, deg+2, nc_tmp.CVCount()); | |
416 | + ON_ClassArray<ON_3dPoint> R(nc_tmp.CVCount()-2); R.SetCount(R.Capacity()); | |
417 | + InitializeSolver_PassThroughPoints_LSQ_EndPointsConstraint(pt_cnt, deg+2, nc_tmp.CVCount(), prm, nc_tmp.m_knot, nc_tmp.KnotCount(), pts, R, solver); | |
418 | + ON_ClassArray<ON_3dPoint> P(nc_tmp.CVCount()); P.SetCount(P.Capacity()); | |
419 | + P[0] = pts[0]; | |
420 | + P[P.Count()-1] = pts[pt_cnt-1]; | |
421 | + solver->Solve(R.Count(), R.First(), P.First()+1, 1); | |
422 | + for (int i = 0; i < nc_tmp.CVCount(); ++i) nc_tmp.SetCV(i, P[i]); | |
423 | + | |
424 | + double err_max = 0; | |
425 | + for (int i = 1; i < pt_cnt - 1; ++i){ | |
426 | + double t = prm[i]; | |
427 | + if (ONGEO_NearestPointNurbsCurve_Newton(nc_tmp, pts[i], t)){ | |
428 | + prm[i] = t; | |
429 | + err[i] = pts[i].DistanceTo(nc_tmp.PointAt(t)); | |
430 | + } | |
431 | + if (err_max < err[i]) err_max = err[i]; | |
432 | + } | |
433 | + } | |
434 | + | |
435 | + return true; | |
436 | +} | |
437 | +#endif | |
438 | + | |
439 | + | |
166 | 440 | bool ONGEO_FitNurbsSurfaceToPointGrid(const ON_3dPoint *pts, int u_count, int v_count, int u_order, int v_order, ONGEO_NI_ParameterMethod pmethod, ON_NurbsSurface &nf, ONGEO_NI_Solver solver_id){ |
167 | 441 | if (u_count < u_order || v_count < v_order || !pts) return false; |
168 | 442 | int n = u_count - 1, m = v_count - 1; |
@@ -187,7 +461,7 @@ | ||
187 | 461 | } |
188 | 462 | for (int i = 0; i <= n; ++i) prm_u[i] /= static_cast<double>(m+1); |
189 | 463 | |
190 | - InitializeSolver(n + 1, u_order, prm_u.First(), nf.m_knot[0], nf.KnotCount(0), solver); | |
464 | + InitializeSolver_ExactlyPassThroughPoints(n + 1, u_order, prm_u.First(), nf.m_knot[0], nf.KnotCount(0), solver); | |
191 | 465 | |
192 | 466 | ON_ClassArray<ON_3dPoint> R((n+1) * (m+1)); R.SetCount(R.Capacity()); |
193 | 467 | for (int j = 0; j <= m; ++j){ |
@@ -203,7 +477,7 @@ | ||
203 | 477 | } |
204 | 478 | for (int j = 0; j <= m; ++j) prm_v[j] /= static_cast<double>(n+1); |
205 | 479 | |
206 | - InitializeSolver(m + 1, v_order, prm_v.First(), nf.m_knot[1], nf.KnotCount(1), solver); | |
480 | + InitializeSolver_ExactlyPassThroughPoints(m + 1, v_order, prm_v.First(), nf.m_knot[1], nf.KnotCount(1), solver); | |
207 | 481 | |
208 | 482 | ON_ClassArray<ON_3dPoint> Q(m+1); Q.SetCount(m+1); |
209 | 483 | for (int i = 0; i <= n; ++i){ |
@@ -0,0 +1,54 @@ | ||
1 | +/* | |
2 | + * FitNurbsToPoints | |
3 | + * Copylight (C) 2020 mocchi | |
4 | + * mocchi_2003@yahoo.co.jp | |
5 | + * License: Boost ver.1 | |
6 | + */ | |
7 | + | |
8 | +#define NOMINMAX | |
9 | +#include "opennurbs.h" | |
10 | +#include "ONGEO.h" | |
11 | +#include <cmath> | |
12 | + | |
13 | +// Reference : Piegl, L. and Tiller, W. | |
14 | +// "The NURBS Book", Second edition. | |
15 | +// Springer-Verlag Berlin, Heidelberg, 1997. | |
16 | + | |
17 | +// pp.229 - 234 : Chapter 6.1, "Point Inversion and Projection for Curves and Surfaces" | |
18 | + | |
19 | +bool ONGEO_NearestPointNurbsCurve_Newton(const ON_NurbsCurve &nc, const ON_3dPoint &pt_query, double &t, ON_Interval *range_, int iteration_max){ | |
20 | + ON_3dPoint pt; | |
21 | + ON_3dVector der1, der2; | |
22 | + ON_Interval range; | |
23 | + if (range_) range = *range_; | |
24 | + else range = nc.Domain(); | |
25 | + for(int i = 0; i < iteration_max; ++i){ | |
26 | + nc.Ev2Der(t, pt, der1, der2); | |
27 | + double der1_len = der1.Length(); | |
28 | + ON_3dVector diff = pt - pt_query; | |
29 | + double diff_len = diff.Length(); | |
30 | + double numerator = ON_DotProduct(der1, diff); | |
31 | + | |
32 | + bool cond_1 = diff_len <= ON_ZERO_TOLERANCE; | |
33 | + bool cond_2 = std::abs(numerator / (der1_len * diff_len)) <= ON_ZERO_TOLERANCE; | |
34 | + bool cond_4 = false; | |
35 | + | |
36 | + if (!cond_1 || !cond_2){ | |
37 | + double tn = t - numerator / (ON_DotProduct(der2, diff) + der1_len * der1_len); | |
38 | + if (!nc.IsClosed()){ | |
39 | + if (tn < range.Min()) tn = range.Min(); | |
40 | + else if (tn > range.Max()) tn = range.Max(); | |
41 | + }else{ | |
42 | + if (tn < range.Min()) tn = range.Max() - (range.Min() - tn); | |
43 | + else if (tn > range.Max()) tn = range.Min() - (tn - range.Max()); | |
44 | + } | |
45 | + cond_4 = std::abs((tn - t) * der1_len) <= ON_ZERO_TOLERANCE; | |
46 | + t = tn; | |
47 | + } | |
48 | + if (cond_1 || cond_2 || cond_4){ | |
49 | + return true; | |
50 | + } | |
51 | + } | |
52 | + return false; | |
53 | +} | |
54 | + |
@@ -79,6 +79,15 @@ | ||
79 | 79 | /// @return クエリ点と最近点との距離 |
80 | 80 | ONGEO_DECL double ONGEO_NearestPointBezierCurve_BBClipping(const ON_BezierCurve *bc_begin, const ON_BezierCurve *bc_end, double tolerance, const ON_3dPoint &pt_query, const ON_BezierCurve *&bc_nearest, double &t, ON_3dPoint &pt_nearest); |
81 | 81 | |
82 | +/// 指定した点に最も近いNurbs曲線上のパラメータをニュートン法による収束計算で求める。 | |
83 | +/// @param [in] nc Nurbs曲線 | |
84 | +/// @param [in] pt_query クエリ点 | |
85 | +/// @param [in/out] t クエリ点に近いNurbs曲線上に位置を与えると、クエリ点に最も近いNurbs曲線上の位置を返す。 | |
86 | +/// @param [in] range_ 解が有りうる範囲を指定する。null_ptr を与えた場合は、曲線のパラメータ範囲全域を対象にする。 | |
87 | +/// @param [in] iteration_max 最も近い点の三次元座標 | |
88 | +/// @return true: 収束成功、 false: 収束失敗 | |
89 | +ONGEO_DECL bool ONGEO_NearestPointNurbsCurve_Newton(const ON_NurbsCurve &nc, const ON_3dPoint &pt_query, double &t, ON_Interval *range_ = 0, int iteration_max = 100000); | |
90 | + | |
82 | 91 | /// 有理ベジエ曲線の曲線長を再分割による収束計算で求める(制御点列長と始点-終点距離の差がトレランス以下になるまで)。 |
83 | 92 | /// @param [in] bc ベジエ曲線 |
84 | 93 | /// @param [in] tolerance トレランス |
@@ -452,11 +461,37 @@ | ||
452 | 461 | /// @param [in] pt_cnt 点列の個数 |
453 | 462 | /// @param [in] order Nurbs曲線の階数 |
454 | 463 | /// @param [in] pmethod 曲線パラメータの計算方法 |
455 | -/// @param [put] nc 作成した Nurbs曲線 | |
464 | +/// @param [out] nc 作成した Nurbs曲線 | |
456 | 465 | /// @param [in] solver 連立一次方程式をどの方法で解くか? |
457 | 466 | /// @return true:成功、false:失敗 |
458 | 467 | ONGEO_DECL bool ONGEO_FitNurbsCurveToPointArray(const ON_3dPoint *pts, int pt_cnt, int order, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, ONGEO_NI_Solver solver = ONGEO_NI_Solver_ON_Matrix_Inverse); |
459 | 468 | |
469 | +/// 与えられた点列から最小二乗法を用いて、指定した制御点数、階数の Nurbs曲線を作成する。 | |
470 | +/// @param [in] pts 点列 | |
471 | +/// @param [in] pt_cnt 点列の個数 | |
472 | +/// @param [in] order Nurbs曲線の階数 | |
473 | +/// @param [in] cpt_cnt 制御点数 | |
474 | +/// @param [in] pmethod 曲線パラメータの計算方法 | |
475 | +/// @param [out] nc 作成した Nurbs曲線 | |
476 | +/// @param [out] prm 各点を曲線に投影した際のパラメータ。呼び出し元で pt_cnt 分の領域を確保すること。null_ptr の時は格納しない。 | |
477 | +/// @param [out] err 各点と曲線との誤差。呼び出し元で pt_cnt 分の領域を確保すること。null_ptr の時は格納しない。 | |
478 | +/// @param [in] solver 連立一次方程式をどの方法で解くか? | |
479 | +/// @return 誤差の最大値 | |
480 | +ONGEO_DECL double ONGEO_FitNurbsCurveToPointArray_LeastSquare(const ON_3dPoint *pts, int pt_cnt, int order, int cpt_cnt, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, double *prm = 0, double *err = 0, ONGEO_NI_Solver solver_id = ONGEO_NI_Solver_ON_Matrix_Inverse); | |
481 | + | |
482 | +#if 0 | |
483 | +/// 与えられた点列から、許容誤差 tolerance 以内を通るNurbs曲線を作成する。 | |
484 | +/// @param [in] pts 点列 | |
485 | +/// @param [in] pt_cnt 点列の個数 | |
486 | +/// @param [in] order Nurbs曲線の階数 | |
487 | +/// @param [in] tolerance 許容誤差 | |
488 | +/// @param [in] pmethod 曲線パラメータの計算方法 | |
489 | +/// @param [out] nc 作成した Nurbs曲線 | |
490 | +/// @param [in] solver 連立一次方程式をどの方法で解くか? | |
491 | +/// @return true:成功、false:失敗 | |
492 | +ONGEO_DECL bool ONGEO_FitNurbsCurveToPointArray(const ON_3dPoint *pts, int pt_cnt, int order, double tolerance, ONGEO_NI_ParameterMethod pmethod, ON_NurbsCurve &nc, ONGEO_NI_Solver solver = ONGEO_NI_Solver_ON_Matrix_Inverse); | |
493 | +#endif | |
494 | + | |
460 | 495 | /// 与えられた点グリッドを通るNurbs曲面を作成する。 |
461 | 496 | /// @param [in] pts 点グリッド (P_u0v0、 P_u1v0、P_u2v0、...、P_u0v1、P_u1v1、... の順に並べる) |
462 | 497 | /// @param [in] u_count 点グリッドのU方向の個数 |
@@ -464,7 +499,7 @@ | ||
464 | 499 | /// @param [in] u_order Nurbs曲面のU方向階数 |
465 | 500 | /// @param [in] v_order Nurbs曲面のV方向階数 |
466 | 501 | /// @param [in] pmethod 曲線パラメータの計算方法 |
467 | -/// @param [put] nf 作成した Nurbs曲面 | |
502 | +/// @param [out] nf 作成した Nurbs曲面 | |
468 | 503 | /// @param [in] solver 連立一次方程式をどの方法で解くか? |
469 | 504 | /// @return true:成功、false:失敗 |
470 | 505 | ONGEO_DECL bool ONGEO_FitNurbsSurfaceToPointGrid(const ON_3dPoint *pts, int u_count, int v_count, int u_order, int v_order, ONGEO_NI_ParameterMethod pmethod, ON_NurbsSurface &nf, ONGEO_NI_Solver solver = ONGEO_NI_Solver_ON_Matrix_Inverse); |