Savitzky Golay for cyclic data

I have 2D angular bodies of which I want to smooth the boundaries using a savitzky golay filter. The boundaries of the bodies are defined by matrices of x- and y-coordinates. The first row and last row of each matrix is identical since the bodies are cyclic. How can I adapt the savitzky Golay filter so it includes the last rows for smoothing the first data point and so that it includes the first coordinates for smoothing the last data point? I want the filter thus to understand that it is a cylcic structure.

When aplying the filter as it is, I do not get a closed circular object as output.

Answer 1

[John D'Errico]

So the first and last rows form a cyclic boundary? Easy peasy. :)

First, make up some garbage data.

R1 = rand(1,10);
A = [R1;rand(10);R1];

A is a 12x10 array, with the desired property.

Now you want to generate a Savitsky-Golay filter for this problem. Say a window of length 5. I'll arbitrarily pick a polynomial order as 1, so a linear model.

The filter is easy enough to generate, although I know what it will look like in advance.

M = pinv([ones(5,1),(-2:2)']);
F = M(1,:);

It should be [1 1 1 1 1]/5. A quick check verifies that fact:

F
F =
          0.2          0.2          0.2          0.2          0.2

So we could now use conv2, with the convolution kernel of F.' to do the Savitsky-Golay. In order to force it to be cyclic though, just one more trick. Copy a couple of rows at the top and bottom, the number of which is based on the window length of the filter.

Ahat = [A(10:11,:);A;A(2:3,:)];
Afilt = conv2(Ahat,F.','valid');
Afilt
Afilt =
      0.52952      0.60457      0.47513      0.50322      0.50705       0.7234      0.49297      0.50984      0.35622      0.69495
      0.48703      0.75185       0.5751      0.57635      0.50425      0.71976      0.59578       0.4411      0.45476      0.72483
      0.64224      0.57291      0.48437      0.55766      0.46489      0.67685       0.5031      0.32243      0.54791      0.63962
      0.54804      0.41666       0.5809      0.67033      0.46353        0.678      0.54297      0.33712      0.60857      0.50357
      0.51663      0.37418      0.47428      0.72159       0.3828      0.52293      0.70295      0.36182      0.70249      0.46339
      0.37867      0.27461      0.40356      0.55024      0.40237       0.4292      0.75475      0.46202      0.76706      0.53382
      0.40356      0.18139      0.27915      0.43136      0.47726      0.55353      0.60434      0.48566      0.68725      0.52565
      0.30876      0.21009      0.26532      0.49641      0.58333      0.42446      0.55892      0.56039      0.56067      0.51497
      0.25662      0.40364      0.29951      0.33418      0.57106       0.4526      0.64386      0.61905      0.48827      0.62547
      0.30555       0.5171       0.3188      0.22171        0.626       0.5964      0.53163      0.64753      0.41491       0.6742
       0.4321      0.59299      0.41492      0.34121      0.56364      0.67361      0.36325      0.58987       0.3041       0.7128
      0.52952      0.60457      0.47513      0.50322      0.50705       0.7234      0.49297      0.50984      0.35622      0.69495

No problem. Afilt is a 12x10 array,

Afilt(1,:) == Afilt(end,:)
ans =
  1×10 logical array
   1   1   1   1   1   1   1   1   1   1

See there was no need to "adapt" the filter kernel. The trick was to properly expand the matrix itself.