bezier_

 function [bezcurve, intcurveyy] = bezier_(points, numofpbc, intcurvexx, fig)
 
 Creates Bezier curve (output 'bezcurve') from 'points' (1st input argument) and the number of points (2nd input argument) 
   and creates from it another interpolated curve (whose x-coordinates are in the input 'intcurvexx' and whose y-coordinates 
   are in the output 'intcurveyy').
 
 INPUTS:
   points:   matrix ((n+1) x 2) with the original points in xy plane
   numofpbc: number of points in the Bezier curve (by default 100)
   intcurvexx: vector with x-coordinates of the interpolation curve. If this argument does not exist or is empty, 
               the program generates Bezier curve, but no interpolation curve
   fig:      any value if you want a figure of points and curve (otherwise, do not enter 4th argument)
                You can enter here the representing symbol for the points (for instance, 'ks' for black squares).
 
 OUTPUTS:
   bezcurve: the Bezier curve, not interpolated, in the format [x y], i.e. a (numofpbc x 2) matrix.
   intcurveyy: vector with y-coordinates (by non-parametric interpolation from intcurvexx) of the interpolation curve;
             it has sense only if intcurvexx elements are monotonically increasing.
 
 Example:    t = (1:100)';
             x = 0.2*randn(size(t)) - sin(pi*t/100) + 0.5*t/100;
             points = [t x];
             bezier_(points, 500, [], 1);  % Creates only the Bezier curve and represents it together with the original points
 Or:         bc = bezier_(points);         % You want the Bezier curve (with 100 points), but no graph nor interpolated curve
 Or:         [bc, intcyy] = bezier_(points, 500, (1:0.1:20)', 1); % You want all the Bezier curve, the interpolation of part of it 
               and the graph of all.

0001 function [bezcurve, intcurveyy] = bezier_(points, numofpbc, intcurvexx, fig)
0002 % function [bezcurve, intcurveyy] = bezier_(points, numofpbc, intcurvexx, fig)
0003 %
0004 % Creates Bezier curve (output 'bezcurve') from 'points' (1st input argument) and the number of points (2nd input argument)
0005 %   and creates from it another interpolated curve (whose x-coordinates are in the input 'intcurvexx' and whose y-coordinates
0006 %   are in the output 'intcurveyy').
0007 %
0008 % INPUTS:
0009 %   points:   matrix ((n+1) x 2) with the original points in xy plane
0010 %   numofpbc: number of points in the Bezier curve (by default 100)
0011 %   intcurvexx: vector with x-coordinates of the interpolation curve. If this argument does not exist or is empty,
0012 %               the program generates Bezier curve, but no interpolation curve
0013 %   fig:      any value if you want a figure of points and curve (otherwise, do not enter 4th argument)
0014 %                You can enter here the representing symbol for the points (for instance, 'ks' for black squares).
0015 %
0016 % OUTPUTS:
0017 %   bezcurve: the Bezier curve, not interpolated, in the format [x y], i.e. a (numofpbc x 2) matrix.
0018 %   intcurveyy: vector with y-coordinates (by non-parametric interpolation from intcurvexx) of the interpolation curve;
0019 %             it has sense only if intcurvexx elements are monotonically increasing.
0020 %
0021 % Example:    t = (1:100)';
0022 %             x = 0.2*randn(size(t)) - sin(pi*t/100) + 0.5*t/100;
0023 %             points = [t x];
0024 %             bezier_(points, 500, [], 1);  % Creates only the Bezier curve and represents it together with the original points
0025 % Or:         bc = bezier_(points);         % You want the Bezier curve (with 100 points), but no graph nor interpolated curve
0026 % Or:         [bc, intcyy] = bezier_(points, 500, (1:0.1:20)', 1); % You want all the Bezier curve, the interpolation of part of it
0027 %               and the graph of all.
0028 
0029 bezcurve = [];
0030 intcurveyy = [];
0031 
0032 if nargin < 1, error('Function bezier_() needs at least 1 argument (the points)'); end
0033 if size(points, 2)~=2, error('Points (1st argument) must be in column format (as a 2 column-matrix [x y]).'); end
0034 n = size(points, 1) - 1;  % number of points is n + 1  (their index goes from 0 to n)
0035 if n > 1000, error('Bézier curve for more than 1000 points is not allowed'); end
0036 
0037 
0038 if nargin < 2, numofpbc = 100; end
0039 if isempty(numofpbc), numofpbc = 100; end
0040 
0041 if nargin < 3, intcurvexx = []; end   % points(:, 1); end
0042 %if isempty(intcurvexx), intcurvexx = points(:, 1); end
0043 
0044 vectort = 0:1/(numofpbc-1):1;
0045 bezcurve = zeros(length(vectort), 2);
0046 
0047 
0048 % In this block we compute each point in the Bezier curve
0049 numpoint = 1;
0050 for t = vectort
0051   suma = [0 0];
0052   for i=0:n
0053     % In next line we add the next point multiplied by the corresponding Bernstein polynomial
0054     suma = suma + points(i+1, :)*nchoosekJH(n, i)*(t^i)*((1 - t)^(n - i));  
0055   end
0056   bezcurve(numpoint, :) = suma;
0057   numpoint = numpoint + 1;
0058 end
0059 
0060 
0061 if ~isempty(intcurvexx)    % The user wants interpolation curve
0062   intcurveyy = interp1(bezcurve(:, 1), bezcurve(:, 2), intcurvexx, 'cubic');  % Makes cubic interpolation from the Bezier curve. Other options: 'linear' or 'spline'
0063 end
0064 
0065 if nargin < 4, return; end     % The user does not want any graphics
0066 
0067 
0068 figure
0069 if ischar(fig)  % If 'fig' is of char type (for instance 'o' or '.'), it is used as plot symbol; otherwise, the program decides
0070   symbol = fig; % the user has decided the symbol
0071 elseif n < 30  % Few points
0072   symbol = 'o';
0073 else  % Many points
0074   symbol = '.';
0075 end  
0076 for j = 1:n
0077   hPoints = plot(points(:,1), points(:,2), symbol);   % plots the points with 'symbol'
0078   hold on;
0079 end
0080 hold on, hbc = plot(bezcurve(:, 1), bezcurve(:, 2), 'LineWidth', 3, 'Color', 'r');  % plots the Bezier curve
0081 if ~isempty(intcurvexx)
0082   hold on, hic = plot(intcurvexx, intcurveyy, 'LineWidth', 2, 'Color', 'g');  % plots the interpolated curve
0083   %h = legend('Bezier curve', 'interpolated curve');
0084   %set(h,'Interpreter','none')
0085   set(get(get(hPoints,'Annotation'),'LegendInformation'),...
0086     'IconDisplayStyle','off'); % Exclude line from legend
0087   legend([hbc hic], {'Bezier curve', 'interp. curve'}); 
0088   %legend([hPoints hbc hic],{'points', 'Bezier curve', 'interp. curve'});
0089   title(sprintf('# pts: %d; # pts in Bezier curve: %d; # pts in interp.curve: %d', n+1, numofpbc, length(intcurvexx)))
0090 else  % No interpolation asked for
0091   %h = legend('Bezier curve');
0092   %set(h,'Interpreter','none')
0093   set(get(get(hPoints, 'Annotation'), 'LegendInformation'), 'IconDisplayStyle', 'off'); % Exclude line from legend
0094   legend([hbc], {'Bezier curve'})
0095   %legend([hPoints hbc], {'points', 'Bezier curve'})
0096   title(sprintf('# pts: %d; # pts in Bezier curve: %d', n+1, numofpbc))
0097 end
0098 
0099 
0100 
0101 function res = nchoosekJH(n ,k)   % (faster and less problematic) alternative to nchoosek(n, k)
0102   if n > 1000, error('Binomial oefficients for n greater than 1000 are not allowed'); end
0103   if k==0, res = 1; return; end
0104   if n==0, res = 0; return; end
0105 res = 1;
0106 for i=1:k
0107   res = res*((n - k + i)/i);
0108 end