fourier_fit

 FOURIER_FIT: use fourier series to interpolate onto a boundary

 [pp] = fourier_fit(points) fits a Fourier series
    points - [x y] contour to be fitted
 [pp] = fourier_fit(points,N) fits a Fourier series and downsamples
    N is the number of Fourier components to fit to.

 [xy] = fourier_fit(pp,  linear_frac, start) returns points at linear_frac
 distance along the contour
    pp          - piecewise polynomial structure
    linear_frac - vector of length fractions (0 to 1) to calculate points
    start       - (optional) a seed for the first point
    xy          - cartesian coordinates of the points


 xy = fourier_fit(points,'fit_spacing',spacing);
 example
  xy = fourier_fit(point,'fit_spacing',0:.1:999);
  th = atan2(xy(:,2)-centroid(2),xy(:,1)-centroid(1));

0001 function [C,th] = fourier_fit(points,N,start);
0002 % FOURIER_FIT: use fourier series to interpolate onto a boundary
0003 %
0004 % [pp] = fourier_fit(points) fits a Fourier series
0005 %    points - [x y] contour to be fitted
0006 % [pp] = fourier_fit(points,N) fits a Fourier series and downsamples
0007 %    N is the number of Fourier components to fit to.
0008 %
0009 % [xy] = fourier_fit(pp,  linear_frac, start) returns points at linear_frac
0010 % distance along the contour
0011 %    pp          - piecewise polynomial structure
0012 %    linear_frac - vector of length fractions (0 to 1) to calculate points
0013 %    start       - (optional) a seed for the first point
0014 %    xy          - cartesian coordinates of the points
0015 %
0016 %
0017 % xy = fourier_fit(points,'fit_spacing',spacing);
0018 % example
0019 %  xy = fourier_fit(point,'fit_spacing',0:.1:999);
0020 %  th = atan2(xy(:,2)-centroid(2),xy(:,1)-centroid(1));
0021 %
0022 
0023 % (C) Andy Adler, 2010. Licenced under GPL v2 or v3
0024 % $Id: fourier_fit.m 6041 2019-12-31 03:55:43Z aadler $
0025 
0026 if ischar(points) && strcmp(points,'UNIT_TEST'); do_unit_test; return ; end
0027 
0028 if nargin<2; N= size(points,1); end
0029 if ischar(N) && strcmp(N,'fit_spacing')
0030    pp = fourier_fit(points, size(points,1),points(1,:));
0031    C = fourier_fit(pp, start);
0032    return
0033 end
0034 
0035 
0036 if size(points,2)==2 % calling analysis function
0037    C = fit_to_fourier(points,N);
0038 elseif size(points,2)==1 % calling synthesis function
0039    if nargin<3; start = []; end
0040    C = fit_from_fourier(points,N,start);
0041 else
0042    error('size of first parameter to fourier_fit not undersood');
0043 end
0044 
0045 % this will crash if N>length(points)
0046 function C = fit_to_fourier(points,N);
0047    z = points*[1;1i];
0048    Z = fft(z,max(N,size(points,1)));
0049    if rem(N,2)==0 % Even
0050      N2 = N/2;
0051      Zlp = Z([1,2:N2,1,end-(N2-2:-1:0)]);
0052      Zlp(N2+1) = 0.5*(Z(N2+1) + Z(end-N2+1)); %centre point
0053    else 
0054      N2 = floor(N/2);
0055      Zlp = Z([1,2:N2+1,end-(N2-1:-1:0)]);
0056    end
0057    C = length(Zlp)/length(Z)*Zlp; % Amplitude scaling
0058     
0059 function xy = fit_from_fourier(C,linear_frac,start);
0060    % Step 1: oversample
0061    N = length(C);
0062    n2 = ceil(N/2);
0063 
0064    pad = zeros(10000,1);
0065    if rem(N,2)==0 % even
0066       Zos = [C(1:n2); C(n2+1)/2; pad; C(n2+1)/2; C(n2+2:end)];
0067    else
0068       Zos = [C(1:n2); pad; C(n2+1:end)];
0069    end
0070    Zos = length(Zos)/length(C)*Zos;
0071    zos = ifft(Zos);
0072    % rearrange the points such that they start as close as possible to the
0073    % seed
0074    if ~isempty(start)
0075        if size(start,2) == 1; start = start'; end
0076        start = start*[1; 1i];
0077        dist = abs(zos-start);
0078        [jnk,p] = min(dist);
0079        zos = circshift(zos,-p+1);
0080    end
0081    
0082    % Step 2:
0083    zos(end+1) = zos(1); % make sure the loop is closed
0084    dpath= abs(diff(zos));
0085    pathlen = [0;cumsum(dpath)];
0086 
0087    % interpolate points onto path
0088    npath = pathlen/max(pathlen);
0089    linear_frac= mod(linear_frac,1);
0090    z_int = interp1(npath, zos, linear_frac);
0091 
0092    xy= [real(z_int(:)), imag(z_int(:))];
0093    
0094 function do_unit_test
0095 
0096 subplot(211);
0097 a=[-0.8981  0.1404;-0.7492  0.5146;-0.2146  0.3504;
0098     0.3162  0.5069; 0.7935  0.2702; 0.9615 -0.2339;
0099     0.6751 -0.8677; 0.0565 -0.6997;-0.3635 -0.8563;
0100    -0.9745 -0.4668];
0101 
0102    C= fourier_fit(a,10);
0103    xy = fourier_fit(C, linspace(.05,1.04,100));
0104    xy2= fourier_fit(C, [.3,.4,.5,.6]);
0105    plot(a(:,1),a(:,2),'r',xy(:,1),xy(:,2),'b',xy2(:,1),xy2(:,2),'m*');
0106 
0107    a(5,:)= [];
0108    C= fourier_fit(a);
0109    xy = fourier_fit(C, linspace(.05,1.04,100));
0110    xy2= fourier_fit(C, [.3,.4,.5,.6]);
0111    plot(a(:,1),a(:,2),'r',xy(:,1),xy(:,2),'b',xy2(:,1),xy2(:,2),'m*');
0112    eidors_msg('VIEW GRAPH TO VERIFY',0);
0113 
0114 subplot(212);
0115  n_elecs = 8; centroid = mean(a);
0116  p = linspace(0,1,n_elecs+1)'; p(end) = [];
0117  xy = fourier_fit(a,'fit_spacing',p);
0118  xyC= fourier_fit(C, linspace(.05,1.04,100));
0119  plot(a(:,1),a(:,2),'r',xy(:,1),xy(:,2),'b*',xyC(:,1),xyC(:,2),'g');
0120 
0121   centroid = mean(a);
0122  th = atan2(xy(:,2)-centroid(2), xy(:,1)-centroid(1));
0123  hold on;
0124  aa = [xy(:,1),centroid(1)+0*xy(:,1)]';
0125  bb = [xy(:,2),centroid(1)+0*xy(:,2)]';
0126  plot(aa,bb,'k-');
0127   unit_test_cmp('fit th',diff(th(1:4)), [ -0.671324440246111;
0128   -1.011260720375461; -0.791975282338955 ],1e-10);
0129 
0130