Savitsky-Golay Filter Problem - Smoothing 3D line

I wanted to smooth a 3d line using the Savitzky-Golay filter, but for this example is seems to not work properly. Any ideas why, and how to fix it?

% Savitzky–Golay filter (sgolayfilt) - smoothing individual axes
windowWidth = 27;    %Standard example values
polynomialOrder = 3;
xsg=sgolayfilt(points(:,1),polynomialOrder, windowWidth);
ysg=sgolayfilt(points(:,2),polynomialOrder, windowWidth);
zsg=sgolayfilt(points(:,3),polynomialOrder, windowWidth);
xyzsg = [xsg,ysg,zsg];

clf()
hold on
plot3(points(:,1),points(:,2),points(:,3),'bo')
plot3(xyzsg(:,1),xyzsg(:,2),xyzsg(:,3),'gx')
hold off

Green crosses show the smoothed version - blue circles is the original data

Cheers

Answer 1

[Image Analyst]

How many elements are in your array? 27 looks like an incredibly wide window for your data. I would have used a width of 5 to 9 elements. Try that and see how it works out.

[William]

167 elements. 5-9 makes a random looking distribution of data points.

I've attached the data if that helps

[William]

[William]

Could it be because the spiral changes z-direction i.e. goes up and then down and then up again? Not too sure...

[Image Analyst]

William, the problem is your data is not parameterized - not sorted. Look what happens if you draw lines between adjacent initial elements when you plot:

You need to fix that. You can fix the bad data by searching for the closest point, or better is to just create the arrays correctly in the first place. Can you do that?

[William]

ah! No, I'm not sure how to do that.

http://uk.mathworks.com/help/matlab/ref/dsearchn.html?

[William]

I've had a go but not quite there yet [UPDATE]

sort_points=zeros(size(points));
sort_points(1,:)=points(1,:);
store_sorted_indices=zeros(size(points,1),1);
for i = 1 : size(points,1)-2;
  X = points(i,1);
  Y = points(i,2);
  Z = points(i,3);
  distances = sqrt((X-points(i+1:end,1)).^2+(Y-points(i+1:end,2)).^2+(Z-points(i+1:end,3)).^2);
  [sortedDistances, sortedIndices] = sort(distances, 'ascend');
  % store sorted indices
  for d = 1 : size(sortedIndices,1)
      % loop through stores finding closest unstored value
      if 0 < find(store_sorted_indices == sortedIndices(d) )
          % if stored index is found -> skip 
          continue
      else
          % store distance index as closest
          store_sorted_indices(i) = sortedIndices(d);
      end
  sort_points(i,:)=points(d,:);
  end
end

Any tips?

[Image Analyst]

No, don't use sort. Just use min(). And there are some other problems with it. You need to build up the list as you go and remove "taken/used" points from the list of points you're checking distances against. I'll see what I can do.....stay tuned....

[Image Analyst]

This works beautifully. The red curve moves smoothly through the noisy original blue data points.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clc;
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 2;

% Savitzky–Golay filter (sgolayfilt) - smoothing individual axes
storedStucture = load('dataline.mat');
points = storedStucture.points;
originalx = points(:, 1);
originaly = points(:, 2);
originalz = points(:, 3);
numPoints = length(originalx);
% The points are not in order.  We need to sort them.
% I'll assume that element 1 is at one end of the curve, then
% I'll take the closest of the remaining points to be the next one.
% And so on until there are no points remaining.
% Create a list of
output = zeros(1, numPoints);
% Start with the first element.
xSoFar = originalx(1);
ySoFar = originaly(1);
zSoFar = originalz(1);
remainingx = originalx(2:end); % Initialize remaining points.
remainingy = originaly(2:end); % Initialize remaining points.
remainingz = originalz(2:end); % Initialize remaining points.
for k = 1 : numPoints-1
  index = output(end);
  % Get this coordinate.
  thisX = xSoFar(k);
  thisY = ySoFar(k);
  thisZ = zSoFar(k);
  % Find distances in 3D space to all of the remaining points.
  distances = sqrt(...
    (thisX - remainingx).^2 + ...
    (thisY - remainingy).^2 + ...
    (thisZ - remainingz).^2);
  % Find the index in the remaining points of the
  % point that is closest to this one.
  [minDistance, indexOfMin] = min(distances);
  % Assign the closest coordinates to the next point
  xSoFar(k+1) = remainingx(indexOfMin);
  ySoFar(k+1) = remainingy(indexOfMin);
  zSoFar(k+1) = remainingz(indexOfMin);
  % Remove that point from the list of remaining points.
  remainingx(indexOfMin) = [];
  remainingy(indexOfMin) = [];
  remainingz(indexOfMin) = [];
end
% Now they're arranged properly and we can proceed.
% Rename the arrays for convenience
x = xSoFar;
y = ySoFar;
z = zSoFar;
% Plot them
plot3(x, y, z, 'b.-', 'LineWidth', 2)
grid on;
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
hold on;

% Do a Savitzky-Golay Filter.
windowWidth = 9;    %Standard example values
polynomialOrder = 3;
% Filter each coordinate list in turn.
xsg=sgolayfilt(x, polynomialOrder, windowWidth);
ysg=sgolayfilt(y, polynomialOrder, windowWidth);
zsg=sgolayfilt(z, polynomialOrder, windowWidth);
% Plot the smoothed curve.
plot3(xsg, ysg, zsg,'r*-', 'LineWidth', 2, 'MarkerSize', 3)
hold off

[William]

nice!

[William]

Any idea how to avoid this? Seems to of missed 2 points and then jumps back to them from the top. -> perhaps just ignore last few points?

[Image Analyst]

That did not happen with the data I used. Is that a new/different set of data?

[William]

yes - i've attached it to the previous post