Matlab实现血管曲率指数 (VTI).zip

function [arclen,seglen] = arc_length(px,py,varargin) % arclength: compute arc length of a space curve, or any curve represented as a sequence of points % usage: [arclen,seglen] = arclength(px,py) % a 2-d curve % usage: [arclen,seglen] = arclength(px,py,pz) % a 3-d space curve % usage: [arclen,seglen] = arclength(px,py,method) % specifies the method used % % Computes the arc length of a function or any % general 2-d, 3-d or higher dimensional space % curve using various methods. % % arguments: (input) % px, py, pz, ... - vectors of length n, defining points % along the curve. n must be at least 2. Replicate % points should not be present in the curve. % % method - (OPTIONAL) string flag - denotes the method % used to compute the arc length of the curve. % % method may be any of 'linear', 'spline', or 'pchip', % or any simple contraction thereof, such as 'lin', % 'sp', or even 'p'. % % method == 'linear' --> Uses a linear chordal % approximation to compute the arc length. % This method is the most efficient. % % method == 'pchip' --> Fits a parametric pchip % approximation, then integrates the % segments numerically. % % method == 'spline' --> Uses a parametric spline % approximation to fit the curves, then % integrates the segments numerically. % Generally for a smooth curve, this % method may be most accurate. % % DEFAULT: 'linear' % % % arguments: (output) % arclen - scalar total arclength of all curve segments % % seglen - arclength of each independent curve segment % there will be n-1 segments for which the % arc length will be computed. % % % Example: % % Compute the length of the perimeter of a unit circle % theta = linspace(0,2*pi,10); % x = cos(theta); % y = sin(theta); % % % The exact value is % 2*pi % % ans = % % 6.28318530717959 % % % linear chord lengths % arclen = arclength(x,y,'l') % % arclen = % % 6.1564 % % % Integrated pchip curve fit % arclen = arclength(x,y,'p') % % arclen = % % 6.2782 % % % Integrated spline fit % arclen = arclength(x,y,'s') % % arclen = % % 6.2856 % % Example: % % A (linear) space curve in 5 dimensions % x = 0:.25:1; % y = x; % z = x; % u = x; % v = x; % % % The length of this curve is simply sqrt(5) % % since the "curve" is merely the diagonal of a % % unit 5 dimensional hyper-cube. % [arclen,seglen] = arclength(x,y,z,u,v,'l') % % arclen = % % 2.23606797749979 % % seglen = % % 0.559016994374947 % % 0.559016994374947 % % 0.559016994374947 % % 0.559016994374947 % % % See also: interparc, spline, pchip, interp1 % % Author: John D'Errico % e-mail: woodchips@rochester.rr.com % Release: 1.0 % Release date: 3/10/2010 % unpack the arguments and check for errors if nargin < 2 error('ARCLENGTH:insufficientarguments', ... 'at least px and py must be supplied') end n = length(px); % are px and py both vectors of the same length? if ~isvector(px) || ~isvector(py) || (length(py) ~= n) error('ARCLENGTH:improperpxorpy', ... 'px and py must be vectors of the same length') elseif n < 2 error('ARCLENGTH:improperpxorpy', ... 'px and py must be vectors of length at least 2') end % compile the curve into one array data = [px(:),py(:)]; % defaults for method and tol method = 'linear'; % which other arguments are included in varargin? if numel(varargin) > 0 % at least one other argument was supplied for i = 1:numel(varargin) arg = varargin{i}; if ischar(arg) % it must be the method validmethods = {'linear' 'pchip' 'spline'}; ind = strmatch(lower(arg),validmethods); if isempty(ind) || (length(ind) > 1) error('ARCLENGTH:invalidmethod', ... 'Invalid method indicated. Only ''linear'',''pchip'',''spline'' allowed.') end method = validmethods{ind}; else % it must be pz, defining a space curve in higher dimensions if numel(arg) ~= n error('ARCLENGTH:inconsistentpz', ... 'pz was supplied, but is inconsistent in size with px and py') end % expand the data array to be a 3-d space curve data = [data,arg(:)]; %#ok end end end % what dimension do we live in? nd = size(data,2); % compute the chordal linear arclengths seglen = sqrt(sum(diff(data,[],1).^2,2)); arclen = sum(seglen); % we can quit if the method was 'linear'. if strcmpi(method,'linear') % we are now done. just exit return end % 'spline' or 'pchip' must have been indicated, % so we will be doing an integration. Save the % linear chord lengths for later use. chordlen = seglen; % compute the splines spl = cell(1,nd); spld = spl; diffarray = [3 0 0;0 2 0;0 0 1;0 0 0]; for i = 1:nd switch method case 'pchip' spl{i} = pchip([0;cumsum(chordlen)],data(:,i)); case 'spline' spl{i} = spline([0;cumsum(chordlen)],data(:,i)); nc = numel(spl{i}.coefs); if nc < 4 % just pretend it has cubic segments spl{i}.coefs = [zeros(1,4-nc),spl{i}.coefs]; spl{i}.order = 4; end end % and now differentiate them xp = spl{i}; xp.coefs = xp.coefs*diffarray; xp.order = 3; spld{i} = xp; end % numerical integration along the curve polyarray = zeros(nd,3); for i = 1:spl{1}.pieces % extract polynomials for the derivatives for j = 1:nd polyarray(j,:) = spld{j}.coefs(i,:); end % integrate the arclength for the i'th segment % using quadgk for the integral. I could have % done this part with an ode solver too. seglen(i) = quadgk(@(t) segkernel(t),0,chordlen(i)); end % and sum the segments arclen = sum(seglen); % ========================== % end main function % ========================== % begin nested functions % ========================== function val = segkernel(t) % sqrt((dx/dt)^2 + (dy/dt)^2) val = zeros(size(t)); for k = 1:nd val = val + polyval(polyarray(k,:),t).^2; end val = sqrt(val); end % function segkernel end % function arclength