How to use smoothing spline to fit a closed curve?

조회 수: 29 (최근 30일)
Wesley
Wesley 2021년 1월 5일
댓글: Sanchit Sinha 2023년 9월 27일
How to use smoothing spline to fit a closed curve?
  댓글 수: 2
KSSV
KSSV 2021년 1월 5일
What data you have?
Wesley
Wesley 2021년 1월 5일
Thank you very much for helping me answer this question, smoothingspline is very effective for non-closed curve fitting. I would like to try to use smoothingspline for closed curve fitting.

댓글을 달려면 로그인하십시오.

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 1월 5일
hello
this will do the trick
all the best ,
load('exampledata.mat')
x = data1(:,1);
y = data1(:,2);
centroid_x = mean(x);
centroid_y = mean(y);
[theta,r] = cart2pol(x-centroid_x,y-centroid_y);
% sort theta in ascending order
[theta_sorted,ind] = sort(theta);
r_sorted = r(ind);
% remove duplicates before interpolation
[theta_unique,IA,IC] = unique(theta_sorted);
r_unique = r_sorted(IA);
% sliding average (smoothing)
N = 50; % adjust smoothing factor
rr = myslidingavg(r_unique, N);
[u,v] = pol2cart(theta_unique,rr);
u = u + centroid_x;
v = v + centroid_y;
plot(x,y,'.b',u,v,'.r');grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = myslidingavg(in, N)
% OUTPUT_ARRAY = MYSLIDINGAVG(INPUT_ARRAY, N)
%
% The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
% the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
% vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
% of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
% * The input argument INPUT_ARRAY is the numerical data array to be processed.
% * The input argument N is the number of neighboring data points to average over for each point of IN.
%
% * The output argument OUTPUT_ARRAY is the output data array.
if (isempty(in)) | (N<=0) % If the input array is empty or N is non-positive,
disp(sprintf('SlidingAvg: (Error) empty input data or N null.')); % an error is reported to the standard output and the
return; % execution of the routine is stopped.
end % if
if (N==1) % If the number of neighbouring points over which the sliding
out = in; % average will be performed is '1', then no average actually occur and
return; % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if % is stopped.
nx = length(in); % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.
if (N>=(2*(nx-1))) % If the number of neighbouring points over which the sliding
out = mean(in)*ones(size(in)); % average will be performed is large enough, then the average actually covers all the points
return; % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if % The execution of the routine is stopped.
out = zeros(size(in)); % In all the other situations, the initialization of the output data structure is performed.
if rem(N,2)~=1 % When N is even, then we proceed in taking the half of it:
m = N/2; % m = N / 2.
else % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
m = (N-1)/2; % m = (N-1) / 2.
end % if
for i=1:nx, % For each element (i-th) contained in the input numerical array, a check must be performed:
dist2start = i-1; % index distance from current index to start index (1)
dist2end = nx-i; % index distance from current index to end index (nx)
if dist2start<m || dist2end<m % if we are close to start / end of data, reduce the mean calculation on centered data vector reduced to available samples
dd = min(dist2start,dist2end); % min of the two distance (start or end)
else
dd = m;
end % if
out(i) = mean(in(i-dd:i+dd)); % mean of centered data , reduced to available samples at both ends of the data vector
end % for i
end
  댓글 수: 2
Bjorn Gustavsson
Bjorn Gustavsson 2021년 1월 5일
Use consolidator to handle multiple data-points at the same independent coordinate, otherwise one "random" point out of each non-unique group will be selected.
Modify myslidingavg to handle the periodic case - i.e. for the last few points include data-points from the start instead of reducing the filtering window etc.
Sanchit Sinha
Sanchit Sinha 2023년 9월 27일
Thankyou your code is wonderful

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Bjorn Gustavsson
Bjorn Gustavsson 2021년 1월 5일
You'll get some progress moving to radial coordinates (it will not be the same properties of a spline in radial coordinates as in cartesian, but...):
subplot(1,2,1)
plot(data1(:,1),data1(:,2),'.')
r0 = mean(data1);
hold on
plot(r0(1),r0(2),'c.','markersize',23)
theta = atan2(data1(:,2)-r0(2),data1(:,1)-r0(1)); % yeah there's a function for this,
r = ((data1(:,2)-r0(2)).^2+(data1(:,1)-r0(1)).^2).^.5;% cart2shpere?
[thetaS,idxS] = sort(theta);
subplot(1,2,2)
plot(thetaS([1:end]),r(idxS([1:end])))
% make a spline-fit of r of theta, make sure to have overlap to get periodic behaviour
knots = -(pi*(1+1/64)):(pi/64):(pi*(1+1/64)); % Check this for rounding problems
sp = spap2(knots([1 1 1:end end end end]),4,...
[theta(end-323:end)-2*pi;thetaS([1:end]);2*pi+thetaS(1:323)],...
[r(idxS(end-323:end));r(idxS([1:end]));r(idxS(1:323))]);
rS = fnval(sp,thetaS);
hold on
plot(thetaS,rS)
subplot(1,2,1)
plot(r0(1)+rS.*cos(thetaS),r0(2)+rS.*sin(thetaS),'r-')
Also have a look at the SLM - shape language modeling toolbox on the file exchange.
HTH

카테고리

Help CenterFile Exchange에서 Get Started with Curve Fitting Toolbox에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by