interpolation

引用:原帖由 fitcat07 於 2014-3-9 09:28 AM 发表
系 spline interpolation？定用 Bezier curves 做？

复制内容到剪贴板代码:template <typename type>
type PointToBezier(const multiarray<vector<type>> &bezier, const vector<type> &point)
{
multiarray<vector<type>> _bezier = BezierDerivative(bezier), _bezier2 = BezierDerivative(_bezier);
auto getInterval = [](const type x) -> type
{
if(x < 0)
return 0;
if(x > 1)
return 1;
return x;
};
auto _NewtonsMethod = [&](type x)
{
x = getInterval(x);
//**********Newton's method
vector<type> d = BezierPoint<type>(bezier, x) - point, _d = BezierPoint<type>(_bezier, x);
x -= Dot(d, _d) / (Dot(d, BezierPoint<type>(_bezier2, x)) + Dot(_d, _d));

return x;
};
auto _Correction = [&](type x) -> type
{
x = getInterval(x);
//**********Correction
type head = Norm(point - bezier.front()),
middle = Norm(point - BezierPoint<type>(bezier, x)),
end = Norm(point - bezier.back());
if(middle > head && head < end)
x = 0;
if(middle > end && head > end)
x = 1;

return x;
};

return _Correction(_NewtonsMethod(_NewtonsMethod(0.5)));
}
template <typename type>
multiarray<vector<type>> &BezierIntersection(const multiarray<vector<type>> &bezier1,
const multiarray<vector<type>> &bezier2, multiarray<vector<type>> &result = multiarray<vector<type>>())
{
multiarray<vector<type>> _bezier1 = BezierDerivative(bezier1), _bezier2 = BezierDerivative(bezier2);
auto getInterval = [](const type x) -> type
{
if (x < 0)
return 0;
if (x > 1)
return 1;
return x;
};
auto _NewtonsMethod = [&](type &x, type &y)
{
x = getInterval(x);
y = getInterval(y);
//**********Newton's method
vector<type> d = BezierPoint<type>(bezier1, x) - BezierPoint<type>(bezier2, y),
_dx = BezierPoint<type>(_bezier1, x), _dy = -BezierPoint<type>(_bezier2, y);
type det = _dx[0] * _dy[1] - _dy[0] * _dx[1];
x -= (_dy[1] * d[0] - _dy[0] * d[1]) / det;
y -= (_dx[0] * d[1] - _dx[1] * d[0]) / det;
};

result.clear();
for (type i = 0; i <= 1; i += (type)1 / ((bezier1.size() - 1) * (bezier2.size() - 1)))
{
type x = PointToBezier(bezier1, BezierPoint(bezier2, i)), y = i;
_NewtonsMethod(x, y);
_NewtonsMethod(x, y);

if (x >= 0 && x <= 1 && y >= 0 && y <= 1
&& Norm(BezierPoint(bezier1, x) - BezierPoint(bezier2, y)) < 1
)
result.push_back(BezierPoint(bezier1, x));
}
return result;
}

附件 未命名.png (150.16 KB)

2014-3-9 08:13 PM

未命名2.png (159.68 KB)

2014-3-9 08:13 PM

未命名.png (150.16 KB)

2014-3-9 08:13 PM

未命名2.png (159.68 KB)

2014-3-9 08:13 PM