How to use smoothing spline to fit a closed curve?

How to use smoothing spline to fit a closed curve?

[KSSV]

What data you have?

[Wesley]

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.

Answer 1

[Mathieu NOE]

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

[Bjorn Gustavsson]

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]

Thankyou your code is wonderful

Answer 2

[Bjorn Gustavsson]

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