Code covered by the BSD License  

Highlights from
Mth Order Piecewise Spline Interpolation

image thumbnail

Mth Order Piecewise Spline Interpolation

by

 

07 Feb 2009 (Updated )

A GUI that plots and gives equations for variable order intermediate splines.

gsp()
function gsp()
% GSP General Piecewise Spline Interpolation
%
% GSP is a graphical user interface which takes x and y points as inputs, along with
% the order m, then outputs mth order splines between each of the x points. In most 
% cases additional constraints are to be entered as well.
% 
% Example:
% 
% Input:
% 
% x points : [1 2 3 4]
% y points : [2 7 -1 5]
% Order    : 2
% *One Additional Constraint is Required*
% f''(1) = 0
% 
% Output:
% Polynomials:
%   0             5            -3
% -13            57           -55
%  27          -183           305
% 
% In the Range:
% 1  2
% 2  3
% 3  4
% 
% This means the following:
% 
% for 1 < x < 2   y = 0*x^2 + 5*x - 3
% for 2 < x < 3   y = -13*x^2 + 57*x - 55
% for 3 < x < 4   y = 27*x^2 - 183*x + 305
% 
% You can obtain the mentioned forms using poly2sym (requires Symbolic Math Toolbox).
% 
% If you want to delete constraints, put its values as '*'. For more help and examples,
% see the pictures accompanied in the zip file. An optional feature can be
% accesessed if the function MQUAKE is present. mquake.m can be found here:
% <http://www.mathworks.com/matlabcentral/fileexchange/22816>
% 
% To launch the GUI, type gsp in the command window with this file in the current directory.
% Alternatively, you can choose Debug -> Run from this editor window, or press F5.
% 
% 
% Husam Aldahiyat, 2009
% numandina@gmail.com
%
global aA

%% figure and uicontrols
figure('units','normalized','position',[.2 .2 .65 .65],'menubar','none','numbertitle','off','color','w','name','General Splines')
axes('position',[.35 .1 .5 .5])
ed1=uicontrol('style','edit','units','normalized','position',[.025 .875 .1 .05],'backgroundcolor','w');
uicontrol('style','text','units','normalized','position',[.025 .94 .1 .025],'backgroundcolor','w','string','x Points');
ed2=uicontrol('style','edit','units','normalized','position',[.025 .75 .1 .05],'backgroundcolor','w');
uicontrol('style','text','units','normalized','position',[.025 .81 .1 .025],'backgroundcolor','w','string','Y Points');
px=uicontrol('style','listbox','units','normalized','position',[.2 .75 .1 .2],'max',2,'backgroundcolor','w',...
	'string',{'Example 1';'Example 2';'Example 3';'Example 4';'Example 5';'Example 6'},'callback',@pn);
ed3=uicontrol('style','edit','units','normalized','position',[.045 .63 .05 .05],'backgroundcolor','w','callback',@cond);
uicontrol('style','text','units','normalized','position',[.025 .69 .1 .025],'backgroundcolor','w','string','Spline Order');
err=uicontrol('style','text','foregroundcolor','r','units','normalized','position',[.025 .55 .22 .05],'backgroundcolor','w',...
	'horizontalalignment','left');
uicontrol('style','pushbutton','units','normalized','position',[.15 .63 .05 .05],'backgroundcolor','w','callback',@go,'string','Solve');
pv=uicontrol('style','pushbutton','units','normalized','position',[.22 .63 .075 .05],'backgroundcolor','w','callback',@dv,...
	'string','Continuity','visible','off');
thih=uicontrol('style','edit','units','normalized','position',[.025 .025 .125 .525],'backgroundcolor','w','max',2,'string','',...
	'horizontalalignment','right','fontsize',13,'fontname','calibri');
edc=uicontrol('style','edit','units','normalized','position',[.15 .025 .075 .525],'backgroundcolor','w','max',2,...
	'horizontalalignment','left','fontsize',13,'fontname','calibri','string','');
Ap=uicontrol('style','listbox','units','normalized','position',[.35 .675 .5 .25],'backgroundcolor','w','max',2,'fontname','courier',...
	'callback',@chv);
Ac=uicontrol('style','listbox','units','normalized','position',[.85 .675 .125 .25],'backgroundcolor','w','max',2,'fontname','courier');
uicontrol('style','text','units','normalized','position',[.35 .935 .5 .025],'backgroundcolor','w','max',2,'string','Spline Polynomials');
uicontrol('style','text','units','normalized','position',[.85 .935 .12 .025],'backgroundcolor','w','max',2,'string','For x in the range');

	function cond(varargin)		
		m=str2double(get(ed3,'string'));
		if m>1
			x=str2num(get(ed1,'string')); %#ok
			if isempty(x)
				set(err,'string','Please give values for x')
				return
			end
			set(err,'string',sprintf('%d Additional Condition(s) Needed',m-1))
			trip=cell(m*2,1);
			for k1=1:m
				koi=repmat('''',1,k1);
				trip{k1}=['f ',koi,'(',num2str(x(1)),') = '];
				trip{k1+m}=['f ',koi,'(',num2str(x(end)),') = '];
			end
			set(thih,'string',trip);
		else
			set(err,'string','No Additional Constraints Required');
			set(edc,'string','')
			set(thih,'string','')
		end
		if m>1
			set(edc,'string',repmat('*',m*2,1))
		end
	end

	function pn(varargin)
		switch get(px,'value')
			case 1
				 set(ed1,'string','0 10')
				 set(ed2,'string','0 0')
				 set(ed3,'string','2')
				 cond
				 set(err,'string','')
				 set(edc,'string',['5';'*';'*';'*'])
			case 2
				 set(ed1,'string','0 1 2')
				 set(ed2,'string','0 1 0')
				 set(ed3,'string','3')
				 cond
				 set(err,'string','')
				 set(edc,'string',['0';'*';'*';'0';'*';'*'])	
			case 3
				 set(ed1,'string','0 1 2')
				 set(ed2,'string','0 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',['0';'0';'*';'*';'*';'0';'0';'*';'*';'*'])				
			case 4
				 set(ed1,'string','1 2 3 4 7')
				 set(ed2,'string','2 -4 0 3 0')
				 set(ed3,'string','3')
				 cond
				 set(err,'string','')
				 set(edc,'string',['*';'0';'*';'*';'0';'*'])				
			case 5
				 set(ed1,'string','1 2 4 6 8 10')
				 set(ed2,'string','2 4 -2 3 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',strvcat('5','0','*','*','*','-4','3','*','*','*')) %#ok			
			case 6			
				 set(ed1,'string','1 2 4 6 8 10')
				 set(ed2,'string','2 4 -2 3 1 0')
				 set(ed3,'string','5')
				 cond
				 set(err,'string','')
				 set(edc,'string',strvcat('-5','0','*','*','*','4','3','*','*','*')) %#ok
		end		
	end

%% Solution

	function go(varargin)
		set(pv,'visible','off')
		x=str2num(get(ed1,'string')); %#ok
		y=str2num(get(ed2,'string')); %#ok
		if length(y)~=length(x)
			set(err,'string','Points are not of equal lengths')
			return
		end
		if any(find(abs(sort(x)-x)>1e-10))
			set(err,'string','x Points need to be ascending')
			return
		end
		if length(unique(x))<length(x)
			set(err,'string','Repeated x points found')
			return
		end
		m=str2double(get(ed3,'string'));
		n=length(x)-1; % number of splines
		g=fliplr(0:m);
		DAT=zeros((m+1)*n);
		B=zeros(n*(1+m),1);
		B(1:n+1)=y;
		S=zeros(n,m+1);
		for k=1:n
			DAT(k,k+m*k-m:k+m*k)=x(k).^g;
		end
		DAT(k+1,k+m*k-m:k+m*k)=x(k+1).^g;
		k=k+2;
		poln=ones(1,m+1);
		for k2=1:m
			pol=[poln,zeros(1,(k2-1))];
			for p=1:n-1
				DAT(k+p+k2*n-k2-n,p*m+p-m:p*m+p+m+1)=...
					repmat(pol,1,2).*repmat(x(p+1),1,(m+1)*2).^repmat(g-(k2-1),1,2).*[ones(1,m+1),-ones(1,m+1)];
			end			
			poln=polyder(poln);
		end
		aa=(get(edc,'string'));	
		if m>1
			con=find(~ismember(1:m*2,findstr(aa(:,1)','*'))==1);
			if length(con)>m-1
				set(err,'string','Too many constraints entered')
				return
			end
			if isempty(con)
				set(err,'string','Too few constraints entered')
				return
			end
			for lp=1:m-1
				poln=ones(1,m+1);
				wo=con(1);
				if wo>m
					con(1)=con(1)-m;
				end
				for h1=1:con(1)
					poln=polyder(poln);
				end
				if wo>m					
					DAT(end-(m-1-lp),end-m:end)=repmat(x(end),1,m+1).^(g-con(1)).*[poln,zeros(1,m+1-length(poln))];
					B(end-(m-1-lp))=str2double(aa(wo,:));
				end
				if wo<=m
					DAT(end-(m-1-lp),1:m+1)=[poln,zeros(1,m+1-length(poln))].*repmat(x(1),1,m+1).^(g-con(1));
					B(end-(m-1-lp))=str2double(aa(con(1),:));
				end
				con(1)=[];
			end
		end
		DAT(isnan(DAT))=0;		
		s=DAT\B;
		C=zeros(n,2);
		for k=1:n
			S(k,:)=s(k*(m+1)-m:k*(m+1));
			C(k,:)=x(k:k+1);
		end
		set(err,'string','')

%% Displaying Results

		set(Ap,'string',num2str(S))
		set(Ac,'string',num2str(C))
		cla reset
		hold on
		st=[1 0 0;0 1 0];
		for k=1:n
			t=linspace(C(k,1),C(k,2),100);
			f=polyval(S(k,:),t);
			aA(k)=plot(t,f);
			set(findobj('color','b'),'color',st(mod(k,2)+1,:),'linewidth',2)			
		end
		plot(x,y,'.')
		set(findobj('marker','.'),'markersize',10,'color','k')
		xlabel('x')
		ylabel('y')
		xlim([x(1)-1 x(end)+1])
		h=axis;
		if h(3)>min(y)-1
			h(3)=min(y)-1;
		end
		if h(4)<max(y)+1
			h(4)=max(y)+1;
		end
		axis(h)
		set(pv,'visible','on')
	end

	function dv(varargin)
		S=str2num(get(Ap,'string')); %#ok
		C=str2num(get(Ac,'string')); %#ok
		st=[1 0 0;0 1 0];
		m=str2double(get(ed3,'string'));
		n=length(str2num(get(ed1,'string')))-1; %#ok
		figure('color','w','menubar','none','numbertitle','off','name','Continuity')
		for k2=1:m
			for g=1:n
				S(g,:)=[zeros(m+1-length([polyder(S(g,1:end-k2+1)),zeros(1,k2)])),polyder(S(g,1:end-k2+1)),zeros(1,k2)];
			end
			subplot(ceil(m/2),ceil(m/ceil(m/2)),k2)
			hold on
			for k=1:n
				t=linspace(C(k,1),C(k,2),100);
				f=polyval(S(k,1:end-k2),t);
				plot(t,f);
				set(findobj('color','b'),'color',st(mod(k,2)+1,:),'linewidth',2)
			end
			title(sprintf('Derivate Order %d',k2))
		end
	end

	function chv(varargin)		
		h=axis;
		try
			mquake(aA(get(Ap,'value')),(h(2)+h(4))/4,.5,1)
		catch
		end
	end

end

Contact us