Perpendicular Points on Splines
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi there,
I have a problem by trying to get specific points on a spline.
What I tried ws to get the points on the red curves, that are perpendicular to the blue curves and corresponding to the points. Which means, that they have the least distance betwee both curves, measured at the blue points.
What I ended up with are points that are perpendicular to the red line and have the least distance to the blue line.
So it is like the perpendicularity is measured on the red line instaed of being measured on the blue line.
This ends up in the fact, that the red points of the other generated splines, are not part of the same flat surface like the vlue points.
I hope this is understandbale in any way.
n = 3;
n1 = n-1;
a = 25;
b = 10;
P = [0 b;0 0;a 0];
T = 15;
H = 7;
R = 1.5;
t = linspace(0,1);
B = bernsteinMatrix(n1,t);
bezierCurve = B*P;
px = bezierCurve(:,1);
pz = bezierCurve(:,2);
py = zeros(100,1);
pty = zeros(T,1);
pt = interparc(T,px,pz,'spline'); % x,y-Werte!!
plot3(px,py,pz,pt(:,1),pty,pt(:,2),'b-o',"MarkerSize",4)
view([0 0])
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);
%Obere Linien
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;
pznew1 = py+newNormalToCurve(:,2).*delta;
d2c = distance2curve([pxnew1 pznew1],pt,'spline');
for i = 1:2
d2ynew = ones(T,1)*R*(-1)^i;
pynew1n = ones(100,1)*R*(-1)^i;
plot3(pxnew1,pynew1n,pznew1,'r')
plot3(d2c(:,1),d2ynew,d2c(:,2),"r-o","MarkerSize",4)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms t
B = bernsteinMatrix(n1,t);
bezierCurve = B*P;
normalToCurve = diff(bezierCurve)*[0 1;-1 0]; %x-Werte sind Abstand der y-Werte zueinander und umgekehrt
normalToCurve = normalToCurve/norm(normalToCurve); %unit normal vector
newNormalToCurve2 = [double(subs(normalToCurve(1),t,linspace(0,1,T)))' double(subs(normalToCurve(2),t,linspace(0,1,T)))'];
f = linspace(0,1,T);
B = bernsteinMatrix(n1,f);
bezierCurve = B*P;
px2 = bezierCurve(:,1);
py2 = bezierCurve(:,2);
% Variable s entlang der Mittellinie definiert
assume(t>=0 & t<=1)
B = bernsteinMatrix(n1,t);
bezierCurve = simplify(B*P);
s(t)=int(norm(diff(bezierCurve)),0,t);
for i = 1:T
S = double(s((i-1)/(T-1)));
v = R*(1-S/double(s(1)))+(H/2)*S/double(s(1));
delta2(i) = v;
end
delta2 = delta2';
pxn2 = pt(:,1)+newNormalToCurve2(:,1).*delta2;
pzn2 = pt(:,2)+newNormalToCurve2(:,2).*delta2;
d2c2 = distance2curve([pxn2 pzn2],pt,'spline');
for i = 1:2
d2ynew = ones(T,1)*R*(-1)^i;
pynew2n = ones(T,1)*R*(-1)^i;
plot3(pxn2,pynew2n,pzn2,'g')
plot3(d2c2(:,1),d2ynew,d2c2(:,2),"g-o","MarkerSize",4)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Untere Linien
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
pxnew2 = px+newNormalToCurve(:,1).*delta;
pznew2 = py+newNormalToCurve(:,2).*delta;
d2c2 = distance2curve([pxnew2 pznew2],pt,'spline');
for i = 1:2
d2ynew = ones(T,1)*R*(-1)^i;
pynew2n = ones(100,1)*R*(-1)^i;
plot3(pxnew2,pynew2n,pznew2,'r')
plot3(d2c2(:,1),d2ynew,d2c2(:,2),"r-o","MarkerSize",4)
end
hold off
%clear all
function [pt,dudt,fofthandle] = interparc(t,px,py,varargin)
% 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
% 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
if (numel(cumarc) == 3) && (method(1) == 's')
cumarc = spl{1}.breaks;
n = numel(cumarc);
chordlen = sum(chordlen);
end
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
[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;
s = totalsplinelength*t;
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
[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
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
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)
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
function [xy,distance,t_a] = distance2curve(curvexy,mapxy,interpmethod)
% check for errors, defaults, etc...
if (nargin < 2)
error('DISTANCE2CURVE:insufficientarguments', ...
'at least curvexy and mapxy must be supplied')
elseif nargin > 3
error('DISTANCE2CURVE:abundantarguments', ...
'Too many arguments were supplied')
end
% get the dimension of the space our points live in
[n,p] = size(curvexy);
if isempty(curvexy) || isempty(mapxy)
% empty begets empty. you might say this was a pointless exercise.
xy = zeros(0,p);
distance = zeros(0,p);
t_a = zeros(0,p);
return
end
% do curvexy and mapxy live in the same space?
if size(mapxy,2) ~= p
error('DISTANCE2CURVE:improperpxorpy', ...
'curvexy and mapxy do not appear to live in the same dimension spaces')
end
% do the points live in at least 2 dimensions?
if p < 2
error('DISTANCE2CURVE:improperpxorpy', ...
'The points MUST live in at least 2 dimensions for any curve to be defined.')
end
% how many points to be mapped to the curve?
m = size(mapxy,1);
% make sure that curvexy and mapxy are doubles, as uint8, etc
% would cause problems down the line.
curvexy = double(curvexy);
mapxy = double(mapxy);
% test for complex inputs
if ~isreal(curvexy) || ~isreal(mapxy)
error('DISTANCE2CURVE:complexinputs','curvexy and mapxy may not be complex')
end
% default for interpmethod
if (nargin < 3) || isempty(interpmethod)
interpmethod = 'linear';
elseif ~ischar(interpmethod)
error('DISTANCE2CURVE:invalidinterpmethod', ...
'Invalid method indicated. Only ''linear'',''pchip'',''spline'' allowed')
else
validmethods = {'linear' 'pchip' 'spline'};
ind = strmatch(lower(interpmethod),validmethods);
if isempty(ind) || (length(ind) > 1)
error('DISTANCE2CURVE:invalidinterpmethod', ...
'Invalid method indicated. Only ''linear'',''pchip'',''spline'' allowed')
end
interpmethod = validmethods{ind};
end
% if the curve is a single point, stop here
if n == 1
% return the appropriate parameters
xy = repmat(curvexy,m,1);
t_a = zeros(m,1);
% 2 norm distance, or sqrt of sum of squares of differences
distance = sqrt(sum(bsxfun(@minus,curvexy,mapxy).^2,2));
% we can drop out here
return
end
% compute the chordal linear arclengths, and scale to [0,1].
seglen = sqrt(sum(diff(curvexy,[],1).^2,2));
t0 = [0;cumsum(seglen)/sum(seglen)];
% We need to build some parametric splines.
% compute the splines, storing the polynomials in one 3-d array
ppsegs = cell(1,p);
% the breaks for the splines will be t0, unless spline got fancy
% on us here.
breaks = t0;
for i = 1:p
switch interpmethod
case 'linear'
dt = diff(t0);
ind = 1:(n-1);
ppsegs{i} = [(curvexy(ind + 1,i) - curvexy(ind,i))./dt,curvexy(ind,i)];
case 'pchip'
spl = pchip(t0,curvexy(:,i));
ppsegs{i} = spl.coefs;
case 'spline'
spl = spline(t0,curvexy(:,i));
breaks = spl.breaks';
nc = numel(spl.coefs);
if nc < 4
% just pretend it has cubic segments
spl.coefs = [zeros(1,4-nc),spl{i}.coefs];
spl.order = 4;
end
ppsegs{i} = spl.coefs;
end
end
nbr = numel(breaks);
% for each point in mapxy, find the closest point to those
% in curvexy. This part we can do in a vectorized form.
pointdistances = ipdm(mapxy,curvexy,'metric',2, ...
'result','structure','subset','nearestneighbor');
% initialize the return variables, using the closest point
% found in the set curvexy.
xy = curvexy(pointdistances.columnindex,:);
distance = pointdistances.distance;
t = t0(pointdistances.columnindex);
% we must now do at least some looping, still vectorized where possible.
% the piecewise linear case is simpler though, so do it separately.
if strcmp(interpmethod,'linear');
% loop over the individual points, vectorizing in the number of
% segments, when there are many segments, but not many points to map.
if n >= (5*m)
% many segments, so loop over the points in mapxy
for i = 1:m
% the i'th point in mapxy
xyi = mapxy(i,:);
% Compute the location (in t) of the minimal distance
% point to xyi, for all lines.
tnum = zeros(nbr - 1,1);
tden = tnum;
for j = 1:p
ppj = ppsegs{j};
tden = tden + ppj(:,1).^2;
tnum = tnum + ppj(:,1).*(xyi(j) - ppj(:,2));
end
tmin = tnum./tden;
% toss out any element of tmin that is less than or equal to
% zero, or or is greater than dt for that segment.
tmin((tmin <= 0) | (tmin >= diff(t0))) = NaN;
% for any segments with a valid minimum distance inside the
% segment itself, compute that distance.
dmin = zeros(nbr - 1,1);
for j = 1:p
ppi = ppsegs{j};
dmin = dmin + (ppi(:,1).*tmin + ppi(:,2) - xyi(j)).^2;
end
dmin = sqrt(dmin);
% what is the minimum distance among these segments?
[mindist,minind] = min(dmin);
if ~isnan(mindist) && (distance(i) > mindist)
% there is a best segment, better than the
% closest point from curvexy.
distance(i) = mindist;
t(i) = tmin(minind) + t0(minind);
for j = 1:p
ppj = ppsegs{j};
xy(i,j) = ppj(minind,1).*tmin(minind) + ppj(minind,2);
end
end
end
else
for i = 1:(n-1)
% the i'th segment of the curve
t1 = t0(i);
t2 = t0(i+1);
% Compute the location (in t) of the minimal distance
% point to mapxy, for all points.
tnum = zeros(m,1);
tden = 0;
for j = 1:p
ppj = ppsegs{j};
tden = tden + ppj(i,1).^2;
tnum = tnum + ppj(i,1).*(mapxy(:,j) - ppj(i,2));
end
tmin = tnum./tden;
% We only care about those points for this segment where there
% is a minimal distance to the segment that is internal to the
% segment.
k = find((tmin > 0) & (tmin < (t2-t1)));
nk = numel(k);
if nk > 0
% for any points with a valid minimum distance inside the
% segment itself, compute that distance.
dmin = zeros(nk,1);
xymin = zeros(nk,p);
for j = 1:p
ppj = ppsegs{j};
xymin(:,j) = ppj(i,1).*tmin(k) + ppj(i,2);
dmin = dmin + (xymin(:,j) - mapxy(k,j)).^2;
end
dmin = sqrt(dmin);
L = dmin < distance(k);
% this segment has a closer point
% closest point from curvexy.
if any(L)
distance(k(L)) = dmin(L);
t(k(L)) = tmin(k(L)) + t0(i);
xy(k(L),:) = xymin(L,:);
end
end
end
end
% for the linear case, t is identical to the fractional arc length
% along the curve.
t_a = t;
else
% cubic segments. here it is simplest to loop over the
% distinct curve segments. We need not test the endpoints
% of the segments, since the call to ipdm did that part.
xytrans = zeros(1,p);
polydiff = @(dp) dp(1:6).*[6 5 4 3 2 1];
for j = 1:(n-1)
% the j'th curve segment
t1 = t0(j);
t2 = t0(j+1);
distpoly0 = zeros(1,7);
for i = 1:p
ppi = ppsegs{i};
% this will allow us to translate each poly to pass through
% (0,0) (i.e., at t = 0)
xytrans(i) = ppi(j,4);
distpoly0(1:2) = distpoly0(1:2) + ppi(j,1).*[ppi(j,1),2*ppi(j,2)];
distpoly0(3) = distpoly0(3) + 2.*ppi(j,1).*ppi(j,3) + ppi(j,2).^2;
distpoly0(4) = distpoly0(4) + 2.*ppi(j,2).*ppi(j,3);
distpoly0(5) = distpoly0(5) + ppi(j,3).^2;
end
for i = 1:m
xyi = mapxy(i,:) - xytrans;
% update the poly for this particular point
% (-2*a1*x0)*t^3
% (-2*a2*x0)*t^2
% (-2*a3*x0)*t
% x0^2
distpoly = distpoly0;
for k = 1:p
ppk = ppsegs{k};
distpoly(4:6) = distpoly(4:6) - 2.*ppk(j,1:3).*xyi(k);
distpoly(7) = distpoly(7) + xyi(k).^2;
end
diffpoly = polydiff(distpoly);
tstationary = roots(diffpoly);
% discard any with an imaginary part, those that are less
% than 0, or greater than t2-t1.
k = (imag(tstationary) ~= 0) | ...
(real(tstationary) <= 0) | ...
(real(tstationary) >= (t2 - t1));
tstationary(k) = [];
% for any solutions that remain, compute the distance.
if ~isempty(tstationary)
mindist = zeros(size(tstationary));
xyij = zeros(numel(tstationary),p);
for k = 1:p
xyij(:,k) = polyval(ppsegs{k}(j,:),tstationary);
mindist = mindist + (mapxy(i,k) - xyij(:,k)).^2;
end
mindist = sqrt(mindist);
% just in case there is more than one stationary point
[mindist,ind] = min(mindist);
if mindist < distance(i)
% we found a point on this segment that is better
% than the endpoint values for that segment.
distance(i) = mindist;
xy(i,:) = xyij(ind,:);
t(i) = tstationary(ind) + t0(j);
end
end % if ~isempty(tstationary)
end % for i = 1:n
end % for j = 1:(n-1)
if nargout >= 2
% build new piecewise polynomials for each segment that
% represent (dx/dt)^2 + (dy/dt)^2 + ...
%
% Since each poly must be cubic at this point, the result will be
% a 4th degree piecewise polynomial.
kernelcoefs = zeros(nbr-1,5);
for i = 1:p
ppi = ppsegs{i};
kernelcoefs = kernelcoefs + [9*ppi(:,1).^2, ...
12*ppi(:,1).*ppi(:,2), ...
4*ppi(:,2).^2 + 6*ppi(:,1).*ppi(:,3), ...
4*ppi(:,2).*ppi(:,3), ppi(:,3).^2];
end
% get the arc length for each segment. quadgk will suffice here
% since we need to integrate the sqrt of each poly
arclengths = zeros(nbr-1,1);
for i = 1:(nbr - 1)
lengthfun = @(t) sqrt(polyval(kernelcoefs(i,:),t));
arclengths(i) = quadgk(lengthfun,0,t0(i+1) - t0(i));
end
% get the cumulative arclengths, then scale by the sum
% this gives us fractional arc lengths.
arclengths = cumsum(arclengths);
totallength = arclengths(end);
arclengths = [0;arclengths/totallength];
% where does each point fall in terms of fractional cumulative
% chordal arclength? (i.e., t0?)
[tbin,tbin] = histc(t,t0);
tbin(tbin < 1) = 1; % being careful at the bottom end
tbin(tbin >= nbr) = nbr - 1; % if the point fell at the very top...
% the total length below the segment in question
t_a = arclengths(tbin);
% now get the piece in the tbin segment
for i = 1:m
lengthfun = @(t) sqrt(polyval(kernelcoefs(tbin(i),:),t));
t_a(i) = t_a(i) + quadgk(lengthfun,0,t(i) - t0(tbin(i)))/totallength;
end
end
end % if strcmp(interpmethod,'linear');
end
% ==========================================================
function d = ipdm(data1,varargin)
% Default property values
params.Metric = 2;
params.Result = 'array';
params.Subset = 'all';
params.Limit = [];
params.ChunkSize = 2^25;
% untangle the arguments
if nargin<1
% if called with no arguments, then the user probably
% needs help. Give it to them.
help ipdm
return
end
% were two sets of data provided?
pvpairs = {};
if nargin==1
% only 1 set of data provided
dataflag = 1;
data2 = [];
else
if ischar(varargin{1})
dataflag = 1;
data2 = [];
pvpairs = varargin;
else
dataflag = 2;
data2 = varargin{1};
if nargin>2
pvpairs = varargin(2:end);
end
end
end
% get data sizes for later
[n1,dim] = size(data1);
if dataflag == 2
n2 = size(data2,1);
end
% Test the class of the input variables
if ~(isa(data1,'double') || isa(data1,'single')) || ...
((dataflag == 2) && ~(isa(data2,'double') || isa(data2,'single')))
error('data points must be either single or double precision variables.')
end
% do we need to process any property/value pairs?
if nargin>2
params = parse_pv_pairs(params,pvpairs);
% check for problems in the properties
% was a legal Subset provided?
if ~isempty(params.Subset) && ~ischar(params.Subset)
error('If provided, ''Subset'' must be character')
elseif isempty(params.Subset)
params.Subset = 'all';
end
valid = {'all','maximum','minimum','largestfew','smallestfew', ...
'nearestneighbor','farthestneighbor'};
ind = find(strncmpi(params.Subset,valid,length(params.Subset)));
if (length(ind)==1)
params.Subset = valid{ind};
else
error(['Invalid Subset: ',params.Subset])
end
% was a limit provided?
if ~ismember(params.Subset,{'all','nearestneighbor','farthestneighbor'}) && ...
isempty(params.Limit)
error('No limit provided, but a Subset that requires a limit value was specified')
end
% check the limit values for validity
if length(params.Limit)>1
error('Limit must be scalar or empty')
end
switch params.Subset
case {'largestfew', 'smallestfew'}
% must be at least 1, and an integer
if (params.Limit<1) || (round(params.Limit)~=params.Limit)
error('Limit must be a positive integer for LargestFew or NearestFew')
end
end
% was a legal Result provided?
if isempty(params.Result)
params.result = 'Array';
elseif ~ischar(params.Result)
error('If provided, ''Result'' must be character or empty')
end
valid = {'array','structure'};
ind = find(strncmpi(params.Result,valid,length(params.Result)));
if (length(ind)==1)
params.Result = valid{ind};
else
error(['Invalid Result: ',params.Subset])
end
% check for the metric
if isempty(params.Metric)
params.Metric = 2;
elseif (length(params.Metric)~=1) || ~ismember(params.Metric,[0 1 2 inf])
error('If supplied, ''Metric'' must be a scalar, and one of [0 1 2 inf]')
end
end % if nargin>2
if (dim == 1) && (params.Metric == 2)
params.Metric = 1;
end
params.usebsxfun = (5==exist('bsxfun','builtin'));
% check for dimension mismatch if 2 sets
if (dataflag==2) && (size(data2,2)~=dim)
error('If 2 point sets provided, then both must have the same number of columns')
end
% Total number of distances to compute, in case I must do it in batches
if dataflag==1
n2 = n1;
end
ntotal = n1*n2;
% FINALLY!!! Compute inter-point distances
switch params.Subset
case 'all'
% One set or two?
if dataflag == 1
d = distcomp(data1,data1,params);
else
d = distcomp(data1,data2,params);
end
% Must we return it as a struct?
if params.Result(1) == 's'
[rind,cind] = ndgrid(1:size(d,1),1:size(d,2));
ds.rowindex = rind(:);
ds.columnindex = cind(:);
ds.distance = d(:);
d = ds;
end
case {'minimum' 'maximum'}
if ((ntotal*8)<=params.ChunkSize) || (params.Result(1) == 'a')
% its small enough to do it all at once
% One set or two?
if dataflag == 1
d = distcomp(data1,data1,params);
else
d = distcomp(data1,data2,params);
end
% Must we return it as a struct?
if params.Result(1) == 'a'
% its an array, fill the unwanted distances with +/- inf
if params.Subset(2) == 'i'
% minimum
d(d<=params.Limit) = -inf;
else
% maximum
d(d>=params.Limit) = +inf;
end
else
% a struct will be returned
if params.Subset(2) == 'i'
% minimum
[dist.rowindex,dist.columnindex] = find(d>=params.Limit);
else
% maximum
[dist.rowindex,dist.columnindex] = find(d<=params.Limit);
end
dist.distance = d(dist.rowindex + n1*(dist.columnindex-1));
d = dist;
end
else
% we need to break this into chunks. This branch
% will always return a struct.
% this is the number of rows of data1 that we will
% process at a time.
bs = floor(params.ChunkSize/(8*n2));
bs = min(n1,max(1,bs));
accum = cell(0,1);
% now loop over the chunks
batch = 1:bs;
while ~isempty(batch)
% One set or two?
if dataflag == 1
dist = distcomp(data1(batch,:),data1,params);
else
dist = distcomp(data1(batch,:),data2,params);
end
% big or small as requested
if ('i'==params.Subset(2))
% minimum value specified
[I,J,V] = find(dist>=params.Limit);
else
% maximum limit
[I,J] = find(dist<=params.Limit);
I = I(:);
J = J(:);
V = dist(I + (J-1)*length(batch));
I = I + (batch(1)-1);
end
% and stuff them into the cell structure
if ~isempty(V)
accum{end+1,1} = [I,J,V(:)]; %#ok
end
% increment the batch
batch = batch + bs;
if batch(end)>n1
batch(batch>n1) = [];
end
end
% convert the cells into one flat array
accum = cell2mat(accum);
if isempty(accum)
d.rowindex = [];
d.columnindex = [];
d.distance = [];
else
% we found something
% sort on the second column, to put them in a reasonable order
accum = sortrows(accum,[2 1]);
d.rowindex = accum(:,1);
d.columnindex = accum(:,2);
d.distance = accum(:,3);
end
end
case {'smallestfew' 'largestfew'}
% find the k smallest/largest distances. k is
% given by params.Limit
% if only 1 set, params.Limit must be less than n*(n-1)/2
if dataflag == 1
params.Limit = min(params.Limit,n1*(n1-1)/2);
end
% is this a large problem?
if ((ntotal*8) <= params.ChunkSize)
% small potatoes
% One set or two?
if dataflag == 1
dist = distcomp(data1,data1,params);
% if only one data set, set the diagonal and
% below that to +/- inf so we don't find it.
temp = find(tril(ones(n1,n1),0));
if params.Subset(1) == 's'
dist(temp) = inf;
else
dist(temp) = -inf;
end
else
dist = distcomp(data1,data2,params);
end
% sort the distances to find those we need
if ('s'==params.Subset(1))
% smallestfew
[val,tags] = sort(dist(:),'ascend');
else
% largestfew
[val,tags] = sort(dist(:),'descend');
end
val = val(1:params.Limit);
tags = tags(1:params.Limit);
% recover the row and column index from the linear
% index returned by sort in tags.
[d.rowindex,d.columnindex] = ind2sub([n1,size(dist,2)],tags);
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse(d.rowindex,d.columnindex,val,n1,size(dist,2));
else
% a structure
d.distance = val;
end
else
% chunks
% this is the number of rows of data1 that we will
% process at a time.
bs = floor(params.ChunkSize/(8*n2));
bs = min(n1,max(1,bs));
if (8*params.Limit*n1/bs) <= params.ChunkSize
accum = cell(0,1);
% now loop over the chunks
batch = (1:bs)';
while ~isempty(batch)
% One set or two?
if dataflag == 1
dist = distcomp(data1(batch,:),data1,params);
k = find(tril(ones(length(batch),n2),batch(1)-1));
if ('s'==params.Subset(1))
dist(k) = inf;
else
dist(k) = -inf;
end
else
dist = distcomp(data1(batch,:),data2,params);
end
% big or small as requested, keeping only the best
% params.Limit number of elements
if ('s'==params.Subset(1))
% minimum value specified
[tags,tags] = sort(dist(:),1,'ascend'); %#ok
tags = tags(1:bs);
[I,J] = ndgrid(batch,1:n2);
ijv = [I(tags),J(tags),dist(tags)];
else
% maximum limit
[tags,tags] = sort(dist(:),1,'descend'); %#ok
tags = tags(1:bs);
[I,J] = ndgrid(batch,1:n2);
ijv = [I(tags),J(tags),dist(tags)];
end
% and stuff them into the cell structure
accum{end+1,1} = ijv; %#ok
% increment the batch
batch = batch + bs;
if batch(end)>n1
batch(batch>n1) = [];
end
end
% convert the cells into one flat array
accum = cell2mat(accum);
% keep only the params.Limit best of those singled out
accum = sortrows(accum,3);
if ('s'==params.Subset(1))
% minimum value specified
accum = accum(1:params.Limit,:);
else
% minimum value specified
accum = accum(end + 1 - (1:params.Limit),:);
end
d.rowindex = accum(:,1);
d.columnindex = accum(:,2);
d.distance = accum(:,3);
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse(d.rowindex,d.columnindex,d.distance,n1,size(dist,2));
end
else
ijv = zeros(0,3);
% loop over the chunks
batch = (1:bs)';
while ~isempty(batch)
% One set or two?
if dataflag == 1
dist = distcomp(data1(batch,:),data1,params);
k = find(tril(ones(length(batch),n2),batch(1)-1));
if ('s'==params.Subset(1))
dist(k) = inf;
else
dist(k) = -inf;
end
else
dist = distcomp(data1(batch,:),data2,params);
end
[I,J] = ndgrid(batch,1:n2);
ijv = [ijv;[I(:),J(:),dist(:)]]; %#ok
% big or small as requested, keeping only the best
% params.Limit number of elements
if size(ijv,1) > params.Limit
if ('s'==params.Subset(1))
% minimum value specified
[tags,tags] = sort(ijv(:,3),1,'ascend'); %#ok
else
[tags,tags] = sort(ijv(:,3),1,'ascend'); %#ok
end
ijv = ijv(tags(1:params.Limit),:);
end
% increment the batch
batch = batch + bs;
if batch(end)>n1
batch(batch>n1) = [];
end
end
% They are fully trimmed down. stuff a structure
d.rowindex = ijv(:,1);
d.columnindex = ijv(:,2);
d.distance = ijv(:,3);
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse(d.rowindex,d.columnindex,d.distance,n1,size(dist,2));
end
end
end
case {'nearestneighbor' 'farthestneighbor'}
% find the closest/farthest neighbor for every point
% is this a large problem? Or a 1-d problem?
if dim == 1
% first split it into the farthest versus nearest cases.
if params.Subset(1) == 'f'
% farthest away
% One set or two?
if dataflag == 1
[d2min,minind] = min(data1);
[d2max,maxind] = max(data1);
else
[d2min,minind] = min(data2);
[d2max,maxind] = max(data2);
end
d.rowindex = (1:n1)';
d.columnindex = repmat(maxind,n1,1);
d.distance = repmat(d2max,n1,1);
% which endpoint was further away?
k = abs((data1 - d2min)) >= abs((data1 - d2max));
if any(k)
d.columnindex(k) = minind;
d.distance(k) = d2min;
end
else
% nearest. this is mainly a sort and some fussing around.
d.rowindex = (1:n1)';
d.columnindex = ones(n1,1);
d.distance = zeros(n1,1);
% One set or two?
if dataflag == 1
% if only one data point, then we are done
if n1 == 2
% if exactly two data points, its trivial
d.columnindex = [2 1];
d.distance = repmat(abs(diff(data1)),2,1);
elseif n1>2
% at least three points. do a sort.
[sorted_data,tags] = sort(data1);
% handle the first and last points separately
d.columnindex(tags(1)) = tags(2);
d.distance(tags(1)) = sorted_data(2) - sorted_data(1);
d.columnindex(tags(end)) = tags(end-1);
d.distance(tags(end)) = sorted_data(end) - sorted_data(end-1);
ind = (2:(n1-1))';
d1 = sorted_data(ind) - sorted_data(ind-1);
d2 = sorted_data(ind+1) - sorted_data(ind);
k = d1 < d2;
d.distance(tags(ind(k))) = d1(k);
d.columnindex(tags(ind(k))) = tags(ind(k)-1);
k = ~k;
d.distance(tags(ind(k))) = d2(k);
d.columnindex(tags(ind(k))) = tags(ind(k)+1);
end % if n1 == 2
else
% Two sets of data. still really a sort and some fuss.
if n2 == 1
% there is only one point in data2
d.distance = abs(data1 - data2);
% d.columnindex is already set correctly
else
ind12 = [1:n1,1:n2]';
bool12 = [zeros(n1,1);ones(n2,1)];
[sorted_data,tags] = sort([data1;data2]);
ind12 = ind12(tags);
bool12 = bool12(tags);
% where did each point end up after the sort?
loc1 = find(~bool12);
loc2 = find(bool12);
% for each point in data1, what is the (sorted) data2
% element which appears most nearly to the left of it?
cs = cumsum(bool12);
leftelement = cs(loc1);
% any points which fell below the minimum element in data2
% will have a zero for the index of the element on their
% left. fix this.
leftelement = max(1,leftelement);
% likewise, any point greater than the max in data2 will
% have an n2 in left element. this too will be a problem
% later, so fix it.
leftelement = min(n2-1,leftelement);
% distance to the left hand element
dleft = abs(sorted_data(loc1) - sorted_data(loc2(leftelement)));
dright = abs(sorted_data(loc1) - sorted_data(loc2(leftelement+1)));
% find the points which are closer to the left element in data2
k = (dleft < dright);
d.distance(ind12(loc1(k))) = dleft(k);
d.columnindex(ind12(loc1(k))) = ind12(loc2(leftelement(k)));
k = ~k;
d.distance(ind12(loc1(k))) = dright(k);
d.columnindex(ind12(loc1(k))) = ind12(loc2(leftelement(k)+1));
end % if n2 == 1
end % if dataflag == 1
end % if params.Subset(1) == 'f'
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
end
elseif (ntotal>1000) && (((params.Metric == 0) && (params.Subset(1) == 'n')) || ...
((params.Metric == inf) && (params.Subset(1) == 'f')))
% do the first dimension
if dataflag == 1
d = ipdm(data1(:,1),data1(:,1),'subset',params.Subset,'metric',params.Metric,'result','struct');
else
d = ipdm(data1(:,1),data2(:,1),'subset',params.Subset,'metric',params.Metric,'result','struct');
end
% its slightly different for nearest versus farthest here
% now, loop over dimensions
for i = 2:dim
if dataflag == 1
di = ipdm(data1(:,i),data1(:,i),'subset',params.Subset,'metric',params.Metric,'result','struct');
else
di = ipdm(data1(:,i),data2(:,i),'subset',params.Subset,'metric',params.Metric,'result','struct');
end
% did any of the distances change?
if params.Metric == 0
% the 0 norm, with nearest neighbour, so take the
% smallest distance in any dimension.
k = d.distance > di.distance;
else
% inf norm. so take the largest distance across dimensions
k = d.distance < di.distance;
end
if any(k)
d.distance(k) = di.distance(k);
d.columnindex(k) = di.columnindex(k);
end
end
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
end
elseif ((ntotal*8) <= params.ChunkSize)
% None of the other special cases apply, so do it using brute
% force for the small potatoes problem.
% One set or two?
if dataflag == 1
dist = distcomp(data1,data1,params);
else
dist = distcomp(data1,data2,params);
end
% if only one data set and if a nearest neighbor
% problem, set the diagonal to +inf so we don't find it.
if (dataflag==1) && (n1>1) && ('n'==params.Subset(1))
diagind = (1:n1) + (0:n1:(n1^2-1));
dist(diagind) = +inf;
end
if ('n'==params.Subset(1))
% nearest
[val,j] = min(dist,[],2);
else
% farthest
[val,j] = max(dist,[],2);
end
% create the matrix as a sparse one or a struct?
if params.Result(1)=='a'
% its an array, so make the array sparse.
d = sparse((1:n1)',j,val,n1,size(dist,2));
else
% a structure
d.rowindex = (1:n1)';
d.columnindex = j;
d.distance = val;
end
else
% break it into chunks
bs = floor(params.ChunkSize/(8*n2));
bs = min(n1,max(1,bs));
% pre-allocate the result
d.rowindex = (1:n1)';
d.columnindex = zeros(n1,1);
d.distance = zeros(n1,1);
% now loop over the chunks
batch = 1:bs;
while ~isempty(batch)
% One set or two?
if dataflag == 1
dist = distcomp(data1(batch,:),data1,params);
else
dist = distcomp(data1(batch,:),data2,params);
end
% if only one data set and if a nearest neighbor
% problem, set the diagonal to +inf so we don't find it.
if (dataflag==1) && (n1>1) && ('n'==params.Subset(1))
diagind = 1:length(batch);
diagind = diagind + (diagind-2+batch(1))*length(batch);
dist(diagind) = +inf;
end
% big or small as requested
if ('n'==params.Subset(1))
% nearest
[val,j] = min(dist,[],2);
else
% farthest
[val,j] = max(dist,[],2);
end
% and stuff them into the result structure
d.columnindex(batch) = j;
d.distance(batch) = val;
% increment the batch
batch = batch + bs;
if batch(end)>n1
batch(batch>n1) = [];
end
end
% did we need to return a struct or an array?
if params.Result(1) == 'a'
% an array. make it a sparse one
d = sparse(d.rowindex,d.columnindex,d.distance,n1,n2);
end
end % if dim == 1
end % switch params.Subset
end
% End of mainline
% ======================================================
% begin subfunctions
% ======================================================
function d = distcomp(set1,set2,params)
% Subfunction to compute all distances between two sets of points
dim = size(set1,2);
if params.usebsxfun
% its a recent enough version of matlab that we can
% use bsxfun at all.
n1 = size(set1,1);
n2 = size(set2,1);
if (dim>1) && ((n1*n2*dim)<=params.ChunkSize)
% its a small enough problem that we might gain by full
% use of bsxfun
switch params.Metric
case 2
d = sum(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim])).^2,3);
case 1
d = sum(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),3);
case inf
d = max(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),[],3);
case 0
d = min(abs(bsxfun(@minus,reshape(set1,[n1,1,dim]),reshape(set2,[1,n2,dim]))),[],3);
end
else
% too big, so that the ChunkSize will have been exceeded, or just 1-d
if params.Metric == 2
d = bsxfun(@minus,set1(:,1),set2(:,1)').^2;
else
d = abs(bsxfun(@minus,set1(:,1),set2(:,1)'));
end
for i=2:dim
switch params.Metric
case 2
d = d + bsxfun(@minus,set1(:,i),set2(:,i)').^2;
case 1
d = d + abs(bsxfun(@minus,set1(:,i),set2(:,i)'));
case inf
d = max(d,abs(bsxfun(@minus,set1(:,i),set2(:,i)')));
case 0
d = min(d,abs(bsxfun(@minus,set1(:,i),set2(:,i)')));
end
end
end
else
% Cannot use bsxfun. Sigh. Do things the hard (and slower) way.
n1 = size(set1,1);
n2 = size(set2,1);
if params.Metric == 2
d = (repmat(set1(:,1),1,n2) - repmat(set2(:,1)',n1,1)).^2;
else
d = abs(repmat(set1(:,1),1,n2) - repmat(set2(:,1)',n1,1));
end
for i=2:dim
switch params.Metric
case 2
d = d + (repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)).^2;
case 1
d = d + abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1));
case inf
d = max(d,abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)));
case 0
d = min(d,abs(repmat(set1(:,i),1,n2) - repmat(set2(:,i)',n1,1)));
end
end
end
% if 2 norm, then we must sqrt at the end
if params.Metric==2
d = sqrt(d);
end
end
% ==============================================================
% end main ipdm
% begin included function - parse_pv_pairs
% ==============================================================
function params=parse_pv_pairs(params,pv_pairs)
npv = length(pv_pairs);
n = npv/2;
if n~=floor(n)
error 'Property/value pairs must come in PAIRS.'
end
if n<=0
% just return the defaults
return
end
if ~isstruct(params)
error 'No structure for defaults was supplied'
end
% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
p_i = lower(pv_pairs{2*i-1});
v_i = pv_pairs{2*i};
ind = strmatch(p_i,lpropnames,'exact');
if isempty(ind)
ind = find(strncmp(p_i,lpropnames,length(p_i)));
if isempty(ind)
error(['No matching property found for: ',pv_pairs{2*i-1}])
elseif length(ind)>1
error(['Ambiguous property name: ',pv_pairs{2*i-1}])
end
end
p_i = propnames{ind};
% override the corresponding default in params.
% Use setfield for comptability issues with older releases.
params = setfield(params,p_i,v_i); %#ok
end
end
25 Kommentare
Antworten (1)
Matt J
am 12 Okt. 2022
I hope this is understandbale in any way.
A figure might help. Regardless, though, you should probably use pdist2 to find nearest points. It should be much simpler than the code you've posted.
7 Kommentare
Matt J
am 13 Okt. 2022
Bearbeitet: Matt J
am 13 Okt. 2022
Then i tried to get the coordinates of 15 points that are on the red line and perpendicular to the blue curve in the blue line.
I'm not sure I understand what it means for a point to be perpendicular to a line. Only lines can be perpendicular to lines.
Perhaps you are saying you want to extend a line from a point on the blue curve, so that the line is perpendicular to the blue curve, and find where this line intersects the red curve. There are ways to do so, but nothing can gaurantee that the red intersection point that you find will minimize the distance from the blue point to the red curve, which is what you said you wanted in your original post: "What I tried ws to get the points on the red curves, that are perpendicular to the blue curves and ... have the least distance betwee both curves"
Below, in any case, is what I imagined you doing with pdist2 to find the minimum distance points on the red curve. The example is for simpler curves that I just made up.
fBlue=@(t) exp(-2*t); %example curve - blue
fRed=@(t) exp(-2*t)+0.5; %example curve - red
xblue=linspace(0,2,10)';
yblue=fBlue(xblue);
xred=linspace(0,30,1e4)';
yred=fRed(xred);
[~,I]=pdist2([xred,yred] ,[xblue,yblue],'euc','Smallest',1);
xred=xred(I); yred=yred(I);
plot(xblue,yblue,'o-b', xred,yred,'o-r'); axis equal
Siehe auch
Kategorien
Mehr zu Spline Construction finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!