Get orthogonal points on other spline
7 次查看(过去 30 天)
显示 更早的评论
Hi There,
my problem is, that I have 2 splines. The first one is devided in 20 equally sized parts. Now I have to get the 20 points on the other Spline, that are orthogonal on der blue spline. I tried it by just deviding the red one as well into 20 parts, but I realised, that the points are not orthogonal then.
I hope somebody can help me with my problem.
Thanks in advance.
(You got to have the "interparc" function by John D'Errico)
n = 3;
n1 = n-1;
a = 20; % Breite des Krummers entlang x
b = 10; % Höhe des Krümmers entlang y
P = [0 b;0 0;a 0];
x = P(:,1); % x-Werte für scatter
y = P(:,2); % y-Werte für scatter
T = 20; % Anazhl an Teilpunkten für den Plot
H = 5;
R = 1;
t = linspace(0,1);
B = bernsteinMatrix(n1,t);
bezierCurve = B*P;
px = bezierCurve(:,1);
py = bezierCurve(:,2);
pt = interparc(T,px,py,'spline'); % x,y-Werte!!
plot(px,py,pt(:,1),pt(:,2),'b-o')
axis equal
axis tight
grid on
hold on
syms t
B = bernsteinMatrix(n1,t);
bezierCurve = B*P;
normalToCurve = diff(bezierCurve)*[0 1;-1 0];
normalToCurve = normalToCurve/norm(normalToCurve);
newNormalToCurve = [double(subs(normalToCurve(1),t,linspace(0,1)))' double(subs(normalToCurve(2),t,linspace(0,1)))'];
f = linspace(0,1);
B = bernsteinMatrix(n1,f);
bezierCurve = B*P;
px = bezierCurve(:,1);
py = bezierCurve(:,2);
assume(t>=0 & t<=1)
B = bernsteinMatrix(n1,t);
bezierCurve = simplify(B*P);
s(t)=int(norm(diff(bezierCurve)),0,t);
for i = 1:100
S = double(s((i-1)/99));
d = R*(1-S/double(s(1)))+(H/2)*S/double(s(1));
delta(i) = d;
end
delta = delta';
pxnew1 = px+newNormalToCurve(:,1).*delta;
pynew1 = py+newNormalToCurve(:,2).*delta;
plot(pxnew1, pynew1,'r')
That's how I tried it, but did not work as I wanted
ptnew1 = interparc(T,pxnew1,pynew1,'spline');
plot(pxnew1,pynew1,ptnew1(:,1),ptnew1(:,2),'r-o')
function [pt,dudt,fofthandle] = interparc(t,px,py,varargin)
% interparc: interpolate points along a curve in 2 or more dimensions
% usage: pt = interparc(t,px,py) % a 2-d curve
% usage: pt = interparc(t,px,py,pz) % a 3-d curve
% usage: pt = interparc(t,px,py,pz,pw,...) % a 4-d or higher dimensional curve
% usage: pt = interparc(t,px,py,method) % a 2-d curve, method is specified
% usage: [pt,dudt,fofthandle] = interparc(t,px,py,...) % also returns derivatives, and a function handle
%
% Interpolates new points at any fractional point along
% the curve defined by a list of points in 2 or more
% dimensions. The curve may be defined by any sequence
% of non-replicated points.
%
% arguments: (input)
% t - vector of numbers, 0 <= t <= 1, that define
% the fractional distance along the curve to
% interpolate the curve at. t = 0 will generate
% the very first point in the point list, and
% t = 1 yields the last point in that list.
% Similarly, t = 0.5 will yield the mid-point
% on the curve in terms of arc length as the
% curve is interpolated by a parametric spline.
%
% If t is a scalar integer, at least 2, then
% it specifies the number of equally spaced
% points in arclength to be generated along
% the curve.
%
% px, py, pz, ... - vectors of length n, defining
% points along the curve. n must be at least 2.
% Exact Replicate points should not be present
% in the curve, although there is no constraint
% that the curve has replicate independent
% variables.
%
% method - (OPTIONAL) string flag - denotes the method
% used to compute the points along the curve.
%
% method may be any of 'linear', 'spline', or 'pchip',
% or any simple contraction thereof, such as 'lin',
% 'sp', or even 'p'.
%
% method == 'linear' --> Uses a linear chordal
% approximation to interpolate the curve.
% This method is the most efficient.
%
% method == 'pchip' --> Uses a parametric pchip
% approximation for the interpolation
% in arc length.
%
% method == 'spline' --> Uses a parametric spline
% approximation for the interpolation in
% arc length. Generally for a smooth curve,
% this method may be most accurate.
%
% method = 'csape' --> if available, this tool will
% allow a periodic spline fit for closed curves.
% ONLY use this method if your points should
% represent a closed curve.
%
% If the last point is NOT the same as the
% first point on the curve, then the curve
% will be forced to be periodic by this option.
% That is, the first point will be replicated
% onto the end.
%
% If csape is not present in your matlab release,
% then an error will result.
%
% DEFAULT: 'spline'
%
%
% arguments: (output)
% pt - Interpolated points at the specified fractional
% distance (in arc length) along the curve.
%
% dudt - when a second return argument is required,
% interparc will return the parametric derivatives
% (dx/dt, dy/dt, dz/dt, ...) as an array.
%
% fofthandle - a function handle, taking numbers in the interval [0,1]
% and evaluating the function at those points.
%
% Extrapolation will not be permitted by this call.
% Any values of t that lie outside of the interval [0,1]
% will be clipped to the endpoints of the curve.
%
% Example:
% % Interpolate a set of unequally spaced points around
% % the perimeter of a unit circle, generating equally
% % spaced points around the perimeter.
% theta = sort(rand(15,1))*2*pi;
% theta(end+1) = theta(1);
% px = cos(theta);
% py = sin(theta);
%
% % interpolate using parametric splines
% pt = interparc(100,px,py,'spline');
%
% % Plot the result
% plot(px,py,'r*',pt(:,1),pt(:,2),'b-o')
% axis([-1.1 1.1 -1.1 1.1])
% axis equal
% grid on
% xlabel X
% ylabel Y
% title 'Points in blue are uniform in arclength around the circle'
%
%
% Example:
% % For the previous set of points, generate exactly 6
% % points around the parametric splines, verifying
% % the uniformity of the arc length interpolant.
% pt = interparc(6,px,py,'spline');
%
% % Convert back to polar form. See that the radius
% % is indeed 1, quite accurately.
% [TH,R] = cart2pol(pt(:,1),pt(:,2))
% % TH =
% % 0.86005
% % 2.1141
% % -2.9117
% % -1.654
% % -0.39649
% % 0.86005
% % R =
% % 1
% % 0.9997
% % 0.9998
% % 0.99999
% % 1.0001
% % 1
%
% % Unwrap the polar angles, and difference them.
% diff(unwrap(TH))
% % ans =
% % 1.2541
% % 1.2573
% % 1.2577
% % 1.2575
% % 1.2565
%
% % Six points around the circle should be separated by
% % 2*pi/5 radians, if they were perfectly uniform. The
% % slight differences are due to the imperfect accuracy
% % of the parametric splines.
% 2*pi/5
% % ans =
% % 1.2566
%
%
% See also: arclength, spline, pchip, interp1
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 3/15/2010
% unpack the arguments and check for errors
if nargin < 3
error('ARCLENGTH:insufficientarguments', ...
'at least t, px, and py must be supplied')
end
t = t(:);
if (numel(t) == 1) && (t > 1) && (rem(t,1) == 0)
% t specifies the number of points to be generated
% equally spaced in arclength
t = linspace(0,1,t)';
elseif any(t < 0) || any(t > 1)
error('ARCLENGTH:impropert', ...
'All elements of t must be 0 <= t <= 1')
end
% how many points will be interpolated?
nt = numel(t);
% the number of points on the curve itself
px = px(:);
py = py(:);
n = numel(px);
% are px and py both vectors of the same length?
if ~isvector(px) || ~isvector(py) || (length(py) ~= n)
error('ARCLENGTH:improperpxorpy', ...
'px and py must be vectors of the same length')
elseif n < 2
error('ARCLENGTH:improperpxorpy', ...
'px and py must be vectors of length at least 2')
end
% compose px and py into a single array. this way,
% if more dimensions are provided, the extension
% is trivial.
pxy = [px,py];
ndim = 2;
% the default method is 'linear'
method = 'spline';
% are there any other arguments?
if nargin > 3
% there are. check the last argument. Is it a string?
if ischar(varargin{end})
method = varargin{end};
varargin(end) = [];
% method may be any of {'linear', 'pchip', 'spline', 'csape'.}
% any any simple contraction thereof.
valid = {'linear', 'pchip', 'spline', 'csape'};
[method,errstr] = validstring(method,valid);
if ~isempty(errstr)
error('INTERPARC:incorrectmethod',errstr)
end
end
% anything that remains in varargin must add
% an additional dimension on the curve/polygon
for i = 1:numel(varargin)
pz = varargin{i};
pz = pz(:);
if numel(pz) ~= n
error('ARCLENGTH:improperpxorpy', ...
'pz must be of the same size as px and py')
end
pxy = [pxy,pz]; %#ok
end
% the final number of dimensions provided
ndim = size(pxy,2);
end
% if csape, then make sure the first point is replicated at the end.
% also test to see if csape is available
if method(1) == 'c'
if exist('csape','file') == 0
error('CSAPE was requested, but you lack the necessary toolbox.')
end
p1 = pxy(1,:);
pend = pxy(end,:);
% get a tolerance on whether the first point is replicated.
if norm(p1 - pend) > 10*eps(norm(max(abs(pxy),[],1)))
% the two end points were not identical, so wrap the curve
pxy(end+1,:) = p1;
nt = nt + 1;
end
end
% preallocate the result, pt
pt = NaN(nt,ndim);
% Compute the chordal (linear) arclength
% of each segment. This will be needed for
% any of the methods.
chordlen = sqrt(sum(diff(pxy,[],1).^2,2));
% Normalize the arclengths to a unit total
chordlen = chordlen/sum(chordlen);
% cumulative arclength
cumarc = [0;cumsum(chordlen)];
% The linear interpolant is trivial. do it as a special case
if method(1) == 'l'
% The linear method.
% which interval did each point fall in, in
% terms of t?
[junk,tbins] = histc(t,cumarc); %#ok
% catch any problems at the ends
tbins((tbins <= 0) | (t <= 0)) = 1;
tbins((tbins >= n) | (t >= 1)) = n - 1;
% interpolate
s = (t - cumarc(tbins))./chordlen(tbins);
% be nice, and allow the code to work on older releases
% that don't have bsxfun
pt = pxy(tbins,:) + (pxy(tbins+1,:) - pxy(tbins,:)).*repmat(s,1,ndim);
% do we need to compute derivatives here?
if nargout > 1
dudt = (pxy(tbins+1,:) - pxy(tbins,:))./repmat(chordlen(tbins),1,ndim);
end
% do we need to create the spline as a piecewise linear function?
if nargout > 2
spl = cell(1,ndim);
for i = 1:ndim
coefs = [diff(pxy(:,i))./diff(cumarc),pxy(1:(end-1),i)];
spl{i} = mkpp(cumarc.',coefs);
end
%create a function handle for evaluation, passing in the splines
fofthandle = @(t) foft(t,spl);
end
% we are done at this point
return
end
% If we drop down to here, we have either a spline
% or csape or pchip interpolant to work with.
% compute parametric splines
spl = cell(1,ndim);
spld = spl;
diffarray = [3 0 0;0 2 0;0 0 1;0 0 0];
for i = 1:ndim
switch method
case 'pchip'
spl{i} = pchip(cumarc,pxy(:,i));
case 'spline'
spl{i} = spline(cumarc,pxy(:,i));
nc = numel(spl{i}.coefs);
if nc < 4
% just pretend it has cubic segments
spl{i}.coefs = [zeros(1,4-nc),spl{i}.coefs];
spl{i}.order = 4;
end
case 'csape'
% csape was specified, so the curve is presumed closed,
% therefore periodic
spl{i} = csape(cumarc,pxy(:,i),'periodic');
nc = numel(spl{i}.coefs);
if nc < 4
% just pretend it has cubic segments
spl{i}.coefs = [zeros(1,4-nc),spl{i}.coefs];
spl{i}.order = 4;
end
end
% and now differentiate them
xp = spl{i};
xp.coefs = xp.coefs*diffarray;
xp.order = 3;
spld{i} = xp;
end
% catch the case where there were exactly three points
% in the curve, and spline was used to generate the
% interpolant. In this case, spline creates a curve with
% only one piece, not two.
if (numel(cumarc) == 3) && (method(1) == 's')
cumarc = spl{1}.breaks;
n = numel(cumarc);
chordlen = sum(chordlen);
end
% Generate the total arclength along the curve
% by integrating each segment and summing the
% results. The integration scheme does its job
% using an ode solver.
% polyarray here contains the derivative polynomials
% for each spline in a given segment
polyarray = zeros(ndim,3);
seglen = zeros(n-1,1);
% options for ode45
opts = odeset('reltol',1.e-9);
for i = 1:spl{1}.pieces
% extract polynomials for the derivatives
for j = 1:ndim
polyarray(j,:) = spld{j}.coefs(i,:);
end
% integrate the arclength for the i'th segment
% using ode45 for the integral. I could have
% done this part with quad too, but then it
% would not have been perfectly (numerically)
% consistent with the next operation in this tool.
[tout,yout] = ode45(@(t,y) segkernel(t,y),[0,chordlen(i)],0,opts); %#ok
seglen(i) = yout(end);
end
% and normalize the segments to have unit total length
totalsplinelength = sum(seglen);
cumseglen = [0;cumsum(seglen)];
% which interval did each point fall into, in
% terms of t, but relative to the cumulative
% arc lengths along the parametric spline?
[junk,tbins] = histc(t*totalsplinelength,cumseglen); %#ok
% catch any problems at the ends
tbins((tbins <= 0) | (t <= 0)) = 1;
tbins((tbins >= n) | (t >= 1)) = n - 1;
% Do the fractional integration within each segment
% for the interpolated points. t is the parameter
% used to define the splines. It is defined in terms
% of a linear chordal arclength. This works nicely when
% a linear piecewise interpolant was used. However,
% what is asked for is an arclength interpolation
% in terms of arclength of the spline itself. Call s
% the arclength traveled along the spline.
s = totalsplinelength*t;
% the ode45 options will now include an events property
% so we can catch zero crossings.
opts = odeset('reltol',1.e-9,'events',@ode_events);
ti = t;
for i = 1:nt
% si is the piece of arc length that we will look
% for in this spline segment.
si = s(i) - cumseglen(tbins(i));
% extract polynomials for the derivatives
% in the interval the point lies in
for j = 1:ndim
polyarray(j,:) = spld{j}.coefs(tbins(i),:);
end
% we need to integrate in t, until the integral
% crosses the specified value of si. Because we
% have defined totalsplinelength, the lengths will
% be normalized at this point to a unit length.
%
% Start the ode solver at -si, so we will just
% look for an event where y crosses zero.
[tout,yout,te,ye] = ode45(@(t,y) segkernel(t,y),[0,chordlen(tbins(i))],-si,opts); %#ok
% we only need that point where a zero crossing occurred
% if no crossing was found, then we can look at each end.
if ~isempty(te)
ti(i) = te(1) + cumarc(tbins(i));
else
% a crossing must have happened at the very
% beginning or the end, and the ode solver
% missed it, not trapping that event.
if abs(yout(1)) < abs(yout(end))
% the event must have been at the start.
ti(i) = tout(1) + cumarc(tbins(i));
else
% the event must have been at the end.
ti(i) = tout(end) + cumarc(tbins(i));
end
end
end
% Interpolate the parametric splines at ti to get
% our interpolated value.
for L = 1:ndim
pt(:,L) = ppval(spl{L},ti);
end
% do we need to compute first derivatives here at each point?
if nargout > 1
dudt = zeros(nt,ndim);
for L = 1:ndim
dudt(:,L) = ppval(spld{L},ti);
end
end
% create a function handle for evaluation, passing in the splines
if nargout > 2
fofthandle = @(t) foft(t,spl);
end
% ===============================================
% nested function for the integration kernel
% ===============================================
function val = segkernel(t,y) %#ok
% sqrt((dx/dt)^2 + (dy/dt)^2 + ...)
val = zeros(size(t));
for k = 1:ndim
val = val + polyval(polyarray(k,:),t).^2;
end
val = sqrt(val);
end % function segkernel
% ===============================================
% nested function for ode45 integration events
% ===============================================
function [value,isterminal,direction] = ode_events(t,y) %#ok
% ode event trap, looking for zero crossings of y.
value = y;
isterminal = ones(size(y));
direction = ones(size(y));
end % function ode_events
end % mainline - interparc
% ===============================================
% end mainline - interparc
% ===============================================
% begin subfunctions
% ===============================================
% ===============================================
% subfunction for evaluation at any point externally
% ===============================================
function f_t = foft(t,spl)
% tool allowing the user to evaluate the interpolant at any given point for any values t in [0,1]
pdim = numel(spl);
f_t = zeros(numel(t),pdim);
% convert t to a column vector, clipping it to [0,1] as we do.
t = max(0,min(1,t(:)));
% just loop over the splines in the cell array of splines
for i = 1:pdim
f_t(:,i) = ppval(spl{i},t);
end
end % function foft
function [str,errorclass] = validstring(arg,valid)
% validstring: compares a string against a set of valid options
% usage: [str,errorclass] = validstring(arg,valid)
%
% If a direct hit, or any unambiguous shortening is found, that
% string is returned. Capitalization is ignored.
%
% arguments: (input)
% arg - character string, to be tested against a list
% of valid choices. Capitalization is ignored.
%
% valid - cellstring array of alternative choices
%
% Arguments: (output)
% str - string - resulting choice resolved from the
% list of valid arguments. If no unambiguous
% choice can be resolved, then str will be empty.
%
% errorclass - string - A string argument that explains
% the error. It will be one of the following
% possibilities:
%
% '' --> No error. An unambiguous match for arg
% was found among the choices.
%
% 'No match found' --> No match was found among
% the choices provided in valid.
%
% 'Ambiguous argument' --> At least two ambiguous
% matches were found among those provided
% in valid.
%
%
% Example:
% valid = {'off' 'on' 'The sky is falling'}
%
%
% See also: parse_pv_pairs, strmatch, strcmpi
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 3/25/2010
ind = find(strncmpi(lower(arg),valid,numel(arg)));
if isempty(ind)
% No hit found
errorclass = 'No match found';
str = '';
elseif (length(ind) > 1)
% Ambiguous arg, hitting more than one of the valid options
errorclass = 'Ambiguous argument';
str = '';
return
else
errorclass = '';
str = valid{ind};
end
end % function validstring
3 个评论
Bruno Luong
2022-10-13
"I need a matrix of 20 points on curve 2, that are cloesest to the 20 points on curve 1.
I think in this context the closest points also mean to be orthogonal to the curve 1"
If I read you correctly this is wrong.
The line that link 2 points that are the closest points (P) on a curve C2 from a given point (Q, on C1) is orthogonal to C2.
So in your case not necessary orthogonal to C1.
The way you describe the problem is so ambigeous that I don't know what you try to do.
回答(1 个)
John D'Errico
2022-10-9
编辑:John D'Errico
2022-10-9
I did NOT say it was impossible to solve. I did say that what you are asking does NOT mean the closest poiints on curve 2 will happen to lie on lines that are orthogonal to BOTH curves. That is provably wrong, in fact.
What you are asking is in fact pretty easy, since I've already posted a tool that does exactly what you asked on the file exchange. (That is not to say the solution was trivial. Not difficult though.) The function distance2curve does exaclty what you are asking to do. That is, given a curve C, and a set of points in space, find the closest point on C to each point in the list. As long as that closest point does not lie at an end point of the curve, it will be on a line that is orthogonal to the curve. So if the curve is a smooth and closed curve, the minimal distance point will always lie on a normal vector.
Find distance2curve here:
However...
As you will see using the curve on a test example, suppose you have two curves. I'll call them C1 and C2. Now, choose an arbitrary set of points on C1. Find the closest point on C2 to each of them. Next, start with that same set of points on C2, and find the closest points on C1. You will NOT in general recover the first set of points for two completely general curves! (It is easy to give a counter-example here that shows my claim to be true. Think about it!)
Again, points of MUTUAL orthogonality on both curves will only lie at points that are at a locally minimal distance between the two curves. There will only be a finite set of such points, and you cannot choose them arbitrarily, unless the curves are quite simple. For example, two parallel lines or two concentric circles would be the simple cases where you have an infinite number of points that have closest points that also lie always on joint normal lines. But those a special cases, and fairly rare.
Anyway, distance2curve solves your problem as asked. I wrote it not long after I wrote interparc. I did consider writing a code that would compute the points of minimal overall distance between two space curves (it would be a little similar in nature), but then I got busy trying to solve world peace, and that just wore me out.
7 个评论
Torsten
2022-10-10
If the three curves were parts of concentric circles, you could expect the three points to lie on a straight line. But why should they in the above constellation ? You are still in the illusion that the lines connecting the three points are orthogonal to all three curves, aren't you ?
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Splines 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!