-
Notifications
You must be signed in to change notification settings - Fork 0
/
foo
241 lines (205 loc) · 7.47 KB
/
foo
1
2
3
4
5
6
7
8
9
berny.m:function val = berny(n, i, x)% This function is a non-recursive formula for cubic Bernstein% Polynomials which form the basis for the cubic Bezier curves% which are used in the supporting programs. The inputs are the% degree of the polynomial, n; the particular curve that is ass-% igned a value of zero up to and including the degree, ii and% the points between [0,1) to be evaluated, x. The output is% points on the curve. It was written by M. R. Holmes.ni = [1 3 3 1]; m = size(x);if n < i val = zeros (m) ;elseif i < 0 val = zeros (m) ;elseif ((n == 0) & (i == 0)) val = 1;else val = ni(i+1) * (x.^i) .* ((ones(m) - x) .^(n-i));end
distEJL.m:% function dt = dist(P). This function computes the initial% distances from the knot points to their adjacent control% points for the initial guess curve. It returns the vector% of distances to iguess.m. The function was written by% E. J. Lane.
distEJL.m:d2 = sqrt(sum(d1.^2))/3; % Computes the initial distances.dt = [d2;d2]; % Assembles the vector of distances.
iguess.m:function [IG, k] = iguess(Q)% function [IG;,k]=iguess(l). This routine takes a set of data point, !% and picks out subset of the data point for knot points, P; computes the% position of the knot points,k; computes the initial distances, dt, to% place the interior control points, C, which are also computed; compute% the angles, ang, of the unit tangent vectors at each knot point; and% assembles the vector IG of paramaters P, ang, and dt, for the curve.% The routine returns the "vector" of parmeters and plots the curve,% its polygon, and the data points in Q. It was written by M. R. Holmes% and revised by E. J. Lane.global dpkpc;[r ,m] = size(Q);disp('Give the number of knotpoints.');n =input(' ');disp('Type "1" for default knot position or "2" to input your own.');h = input (' ');if h == 1 k = defk(m,n); % Calls for default knot position.elseif h == 2 disp('Input initial knot sequence as follows "[1 4 8 ... n]".') k =input(' ');elseif h ~= 1 | h ~= 2 disp('Error! Start over and choose "1" or "2".') ,pause(2) iguessenddpkpc = k; % Position of knot points passed globally.P = knots(Q,k); % call to compute the knotpoints.dt = distEJL(P); % Call to compute the distance between% successive ve knot points.ang = tang(Q,k); % Call to compute the angles for% the unit tangent vectors.c = ctpts(P,ang,dt); % Call to compute the control points;% for the curve.pltC(C,Q,P); % Call to plot the initial guess curve,% its control polygon, and points in Q.% Assemble the composite vector of the initial guess curve parametersIG = [P(1,:) P(2,:) ang dt(1,:) dt(2,:)];
newk.m:function nk = newk(Q,P)% function nk = newk(Q,P). This function takes data points, Q,and% knot points, P, as input. The function finds the closest data% point of the cubic segment that is associated with that knot% point, and returns a new k-array, nk. It ensures the data points% of Q are properly associated with the proper segment of the curve.% "dpkpc" is a global variable that is initially equal to the old% k. This was written by M. R. Holmes and revised by E. J. Lane.global dpkpc[r,m] = size(Q);[s,n] = size(P);nk(1) = 1; nk(n) = m; % Ensures the knot sequence starts and % ends with the 1st and last points in Q.for i = 2:n-1 js = dpkpc(i-1); je = dpkpc(i+l); % variables to pick out jm = dpkpc(i); % interior knot positions. z = je - js + l; mm = jm - js + 1; R = Q(:,js:je) - P(:,i) * ones(l,z); % Finds differences % between data points and % knot point being checked. for jj = l:z D(jj) = R(:,jj)'*R(:,jj); % Ensures differences are % positive for comparison. end if mm < z sd = sign(D(mm) - D(mm+1)) ; % Compares for smallest % difference to find elseif mm > 1 % new dividing points. sd = sign(D(mm-1) - D(mm)); else sd = 0; end while D(mm) - D(mm+sd) > 0 if (mm == 2) && (sd < 0) break; end if (mm == m-1) && (sd > 0) break; end mm = mm + sd; end nk(i)= mm+ js - 1; % knot positions.enddpkpc = nk; % knot positions or sequence.
objf2.m:function se = objf2(x,Q,t,k)% function se = objf2(x,Q,t,k). This is objective function that % will be minimized by "fmins”. The input arguments are the vector x of% parameters for the curve, data points Q, a toggle t if a new knot has% been inserted or one removed, and the knot sequence, k. The output is% the sum from the function “sod” plus the the distance square from the% first and last data points to the first and last knot points, respec0-% tively. This was written by M. R. Holmes and revised by E. J. Lane.global dpkpc% LLoop to change dpkpc if a knot was inserted or removed.if t == 1 dpkpc = k; t = n; global dpkpcendif t == 0 global dpkpc [r,s] = size(Q); [P,ang,dt] = ktangdt(x); % Call to separate x into its subcomponents. C = ctpts(P,ang,dt); % Call to compute control points. dpkpc = newk(Q, P); % Calls function computes the % new dividing point positions . m = lenqth(x); n = round(m/5); fp = P(:,1) - Q (:, 1); % Compute the distance squared lp = P(:,n) - Q(:,s); % from the first and last data points to % the first and last knot points, respectively. % Calls the function that computes the sums of the square % of the distances from the data points to the nearest point % on the cubic segment. se = sod(C,Q,dpkpc) + fp'*fp + lp'*lp;endend
opdist.m:function se2 = opdist(id,Q,P,ang)% function se2 = opdist(id,Q,P,ang). This function is the obj% ective function to be minimized during segment optimization.% It receives subcomponents from a curve's composite vector, a% segment at a time. It returns the sum error for a segment.% It was written by E. J. Lane.n = length(Q);C = ctpts(P,ang,id); % Call to compute the segments% control points.se2 = 0;% Loop to find distance error in a segment.for j = 1:n np = NearestPoint(C', Q(:,j)'); if (j==1) && (np==C(:,1)) d = zeros(1,2); elseif (j==n) && (np==C(:,4)') d = zeros(1,2); else d = (Q(:,j)' - np); end se2 = se2 + d*d';end
tang.m:% function ang = tang(Q,k). This function computes the angles% for the unit tangent vectors at the knot points.u = unitv(Q,k); % Call to compute unit tangents.ang = atan2(u(2,:), u(1,:)); % Converts tangents to angles.
unitv.m:function uv = unitv(Q,k)% function uv = unitv(Q,k). This function takes data points,% Q, and the position of knotpoints, k, as input variables.% It uses chord length parameterization to fit a parametric% quadratic curve to five data points. The unit tangent vec-% tors are approximated by the unit tangent vectors for these% quadratic functions. It returns the set of unit tangent% vectors in the direction of the knot points. It was written% by M. R. Holmes[r, m] = size (Q);n = length(k);for j = 1:n % Loop to index knot positions. if j == 1 k(j) = 1; kt = 1; elseif j == n k(j) = m-4; kt = 5; else k(j) = k(j)-2; kt = 3; end % Extracting the knot point and four adjacent points. x = Q(1, k(j):k(j)+4)'; y = Q(2, k(j):k(j)+4)'; xd = diff(x); yd= diff(y); % Get chord length. d = sqrt(xd.*xd + yd.*yd); t(1) = 0; t(2) = d(1); t(3) = t(2) + d(2); t(4) = t(3) + d(3); t(5) = t(4) + d(4); c = [ones(5,1) t' (t.*t)'] \ [x y]; u = c(2, :) + 2*c(3,:)*t(kt); u = u/norm(u); uv(:,j) = u'; % Approximation of unit tangents % by unit tangent to quadratic.end