piece_poly_fit

 PIECE_POLY_FIT: piecewise polynomial fitting toolset

 [pp m] = pp_fit(points, fstr) fits a piecewise polynomial to a contour
    points - [x y] contour to be fitted
    fstr   - function to use: 'spline' or 'pchip' (default)
    m      - returns the mean of the original points which was subtracted
             for fitting

 [th xy] = pp_fit(pp, start_th, linear_frac) returns points at linear_frac
 distance along the contour using start_th as the starting angle.
    pp          - piecewise polynomial structure
    start_th    - starting angle
    linear_frac - vector of length fractions (0 to 1) to calculate points
    th          - angles of the desired points (-pi to pi) (use ppval to
                  obtain the radii.
    xy          - cartesian coordinates of the points

0001 function [pp, m] = piece_poly_fit(points, fstr, linear_frac)
0002 % PIECE_POLY_FIT: piecewise polynomial fitting toolset
0003 %
0004 % [pp m] = pp_fit(points, fstr) fits a piecewise polynomial to a contour
0005 %    points - [x y] contour to be fitted
0006 %    fstr   - function to use: 'spline' or 'pchip' (default)
0007 %    m      - returns the mean of the original points which was subtracted
0008 %             for fitting
0009 %
0010 % [th xy] = pp_fit(pp, start_th, linear_frac) returns points at linear_frac
0011 % distance along the contour using start_th as the starting angle.
0012 %    pp          - piecewise polynomial structure
0013 %    start_th    - starting angle
0014 %    linear_frac - vector of length fractions (0 to 1) to calculate points
0015 %    th          - angles of the desired points (-pi to pi) (use ppval to
0016 %                  obtain the radii.
0017 %    xy          - cartesian coordinates of the points
0018 
0019 % (C) Bartlomiej Grychtol, 2010. Licenced under GPL v2 or v3
0020 % $Id: piece_poly_fit.m 5112 2015-06-14 13:00:41Z aadler $
0021 
0022 if ischar(points) && strcmp(points,'UNIT_TEST'); do_unit_test; return ; end
0023 
0024 if nargin < 2
0025     fstr = 'pchip';
0026 end
0027 m = [];
0028 if isstruct(points)
0029     pts = max(500,length(linear_frac)*3);
0030     start_th = fstr;
0031     [pp m] = path_len(points, pts, start_th, linear_frac );
0032 else
0033     [pp m] = fit_to_pp(points, fstr);
0034 end
0035 
0036 %%
0037 function [pp centroid] = fit_to_pp(points, fstr)
0038 
0039 % 0. Subtract the centroid and convert to polar coords
0040 centroid = mean(points);
0041 n_points = size(points,1);
0042 points = points - repmat(centroid, [n_points,1]);
0043 [ppoints(:,1), ppoints(:,2)] = cart2pol(points(:,1), points(:,2));
0044 % 1. close the loop:
0045 ppoints = sortrows(ppoints,1);
0046 r = ppoints(:,2);rho = ppoints(:,1);
0047 % add points at +/- pi (if none present)
0048 if rho(1) == -pi
0049     r = [r; r(1)];
0050     rho = [rho; pi];
0051 elseif rho(end) == pi
0052     r = [r(end); r];
0053     rho = [-pi; rho];
0054 else
0055     dist = 2*pi - rho(end) + rho(1);
0056     m = r(end)*(pi+rho(1))/dist + r(1)*(pi-rho(end))/dist;
0057     r = [m; r; m];
0058     rho = [-pi; ppoints(:,1); pi  ];
0059 end
0060 % 2. fit
0061 switch fstr
0062     case 'pchip'
0063         pp=pchip(rho,r);
0064     case 'spline'
0065         df = (r(2) - r(end-1)) / ( rho(end-1) + rho(2));
0066         pp=spline(rho, [df;r;df]);
0067 %         pp=spline(rho, r);
0068 end
0069 
0070 function r = fit_from_pp(pp,rho)
0071 r = ppval(pp,rho);
0072 
0073 % start_th is starting angle for interpolation
0074 % linear_frac is length fraction at which to find the theta => [0.1, 0.5,]
0075 function [th_frac xy]  = path_len( pp, pts, start_th, linear_frac )
0076    th = linspace(start_th, start_th+2*pi,pts+1)';
0077    [x,y] = pol2cart(th,ones(size(th)));
0078    [th,jnk] = cart2pol(x,y); % 2nd param for octave bug (3.6.2)
0079    th(end) = [];
0080    th = sortrows(th);
0081    dth= diff(th);
0082    rad = fit_from_pp(pp, th);
0083    drad= diff(rad);
0084    rad_= 0.5*(rad(1:end-1) + rad(2:end));
0085    dlen= sqrt( (rad_ .* dth).^2 + drad.^2 );
0086    pathlen = [0;cumsum(dlen)];
0087    
0088    npath = pathlen/max(pathlen);
0089    linear_frac = linear_frac + 0.5 + start_th / (2*pi) ;
0090    idx = find(linear_frac > 1);
0091    linear_frac(idx) = linear_frac(idx) - 1;
0092    th_frac = interp1(npath, th, linear_frac);
0093 
0094    lrad = fit_from_pp(pp,th_frac);
0095    [xe,ye]= pol2cart( th_frac, lrad);
0096    xy = [xe ye];
0097    
0098 %%
0099 function do_unit_test
0100 %     a = [
0101 %    -0.8981   -0.7492   -0.2146    0.3162    0.7935    0.9615    0.6751    0.0565   -0.3635   -0.9745
0102 %     0.1404    0.5146    0.3504    0.5069    0.2702   -0.2339   -0.8677   -0.6997   -0.8563   -0.4668 ]';
0103 %
0104 % centroid = mean(a);
0105 % n_points = size(a,1);
0106 % a = a - repmat(centroid, [n_points,1]);
0107 
0108 th=linspace(0,2*pi,33)'; th(end)=[];
0109 a=[sin(th)*0.3,cos(th)];
0110 
0111 
0112   pp= piece_poly_fit(a,'pchip');
0113    
0114    th2 = linspace(-pi,pi,33)';th2(end)=[];
0115    r2 = ppval(pp,th2);
0116    [xf,yf] = pol2cart(th2,r2);
0117 
0118 
0119    
0120    [lfrac xy] = path_len( pp, 100, 0, linspace(.5,.8,6)' );
0121    
0122    clf
0123    subplot(121)
0124    plot(a(:,1),a(:,2),'*',xf,yf,'b-', xy(:,1),xy(:,2),'r+');
0125    axis equal
0126 
0127    pp= piece_poly_fit(a,'spline');
0128    th2 = linspace(-pi,pi,33)';th2(end)=[];
0129    r2 = ppval(pp,th2);
0130    [x,y] = pol2cart(th2,r2);
0131    subplot(122)
0132    plot(a(:,1),a(:,2),'r+'); hold on
0133    plot(x,y,'.');
0134    axis equal