eeg_topo_surfacespline_script - explore spline toolbox for topographic mapping
0001 % eeg_topo_surfacespline_script - explore spline toolbox for topographic mapping 0002 0003 0004 % If you want to interpolate at sites other than the breaks and/or by 0005 % splines other than cubic splines with simple knots, then you use 0006 % the spapi command. In its simplest form, you would say: 0007 % 0008 % sp = spapi(k,x,y); 0009 % 0010 % in which the first argument, k, specifies the order of the interpolating 0011 % spline; this is the number of coefficients in each polynomial piece, 0012 % i.e., 1 more than the nominal degree of its polynomial pieces. For 0013 % example, the next figure shows a linear, a quadratic, and a quartic 0014 % spline interpolant to our data, as obtained by the statements: 0015 % 0016 % x = (4*pi)*[0 1 rand(1,20)]; y = sin(x); 0017 % figure 0018 % sp2 = spapi(2,x,y); fnplt(sp2), hold on 0019 % sp3 = spapi(3,x,y); fnplt(sp3,'--') 0020 % sp5 = spapi(5,x,y); fnplt(sp5,'-.'), plot(x,y,'o') 0021 % legend('linear','quadratic','quartic','data'), hold off 0022 0023 0024 % What if the data are noisy? For example, suppose that the given values are: 0025 % 0026 % noisy = y + .1*(rand(size(x))-.5); 0027 % 0028 % Then you might prefer to approximate instead. For example, you might try 0029 % the cubic smoothing spline, obtained by the command: 0030 % 0031 % scs = csaps(x,noisy); 0032 % fnplt(scs), hold on, plot(x,noisy,'o') 0033 % legend('smoothing spline','noisy data'), hold off 0034 0035 0036 % If you don't like the level of smoothing done by csaps(x,y), you can change 0037 % it by specifying the smoothing parameter, p, as an optional third argument. 0038 % Choose this number anywhere between 0 and 1. As p changes from 0 to 1, the 0039 % smoothing spline changes, correspondingly, from one extreme, the least-squares 0040 % straight-line approximation to the data, to the other extreme, the "natural" 0041 % cubic spline interpolant to the data. Since csaps returns the smoothing 0042 % parameter actually used as an optional second output, you could now experiment, 0043 % as follows: 0044 % 0045 % [scs,p] = csaps(x,noisy); fnplt(scs), hold on 0046 % fnplt(csaps(x,noisy,p/2),'--') 0047 % fnplt(csaps(x,noisy,(1+p)/2),':'), plot(x,noisy,'o') 0048 % legend('smoothing spline','more smoothed','less smoothed','noisy data'), hold off 0049 0050 0051 % At times, you might prefer simply to get the smoothest cubic spline sp that is 0052 % within a specified tolerance tol of the given data in the sense that: 0053 % 0054 % norm(noisy - fnval(sp,x))^2 <= tol 0055 % 0056 % This spline is provided by the command: 0057 % 0058 % sp = spaps(x,noisy,tol); 0059 0060 0061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0062 % If f is one of the cs, ch, or sp splines then, it can be displayed by the statement 0063 % fnplt(f) 0064 % Its value at a is given by the statement: 0065 % fnval(f,a); 0066 % Its second derivative is constructed by the statement: 0067 % DDf = fnder(fnder(f)); 0068 % or by the statement: 0069 % DDf = fnder(f,2); 0070 % Its definite integral over the interval [a..b] is supplied by the statement: 0071 % [1 -1]*fnval(fnint(f),[b;a]); 0072 % and the difference between the spline in cs and the one in ch can be computed as 0073 % fncmb(cs,'-',sp); 0074 % 0075 % 0076 % The toolbox supports vector-valued splines. For example, if you want a spline curve 0077 % through given planar points (x(i), y(i)), i = 1, ..., n, then the statements 0078 % 0079 % xy = [x;y]; df = diff(xy.').'; 0080 % t = cumsum([0, sqrt([1 1]*(df.*df))]); 0081 % cv = csapi(t,xy); 0082 % 0083 % provide such a spline curve, using chord-length parametrization and cubic spline 0084 % interpolation with the not-a-knot end condition, as can be verified by the 0085 % statements 0086 % 0087 % fnplt(cv), hold on, plot(x,y,'o'), hold off 0088 % 0089 % 0090 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0091 % 0092 % Vector-valued splines are also used in the approximation to gridded data, in 0093 % any number of variables, using tensor-product splines. The same spline-construction 0094 % commands are used, only the form of the input differs. For example, 0095 % if x is an m-vector, y is an n-vector, and z is an array of size [m,n], then 0096 % 0097 % cs = csapi({x,y},z); 0098 % 0099 % describes a bicubic spline f satisfying f( x(i), y(j) ) = z(i,j) for 0100 % i = 1:m, j = 1:n. Such a multivariate spline can be vector-valued. For example, 0101 % 0102 % x = 0:4; y=-2:2; s2 = 1/sqrt(2); 0103 % z(3,:,:) = [0 1 s2 0 -s2 -1 0].'*[1 1 1 1 1]; 0104 % z(2,:,:) = [1 0 s2 1 s2 0 -1].'*[0 1 0 -1 0]; 0105 % z(1,:,:) = [1 0 s2 1 s2 0 -1].'*[1 0 -1 0 1]; 0106 % sph = csape({x,y},z,{'clamped','periodic'}); 0107 % fnplt(sph), axis equal, axis off 0108 % 0109 % gives a perfectly acceptable sphere. Its projection onto the (x,z)-plane is plotted by 0110 % 0111 % fnplt(fncmb(sph,[1 0 0; 0 0 1])) 0112 % 0113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0114 0115 0116 0117 clear; close all; 0118 0119 x = .0001 + [-4:.2:4]; 0120 y = -3:.2:3; 0121 0122 [yy,xx] = meshgrid(y,x); 0123 0124 r = pi * sqrt(xx.^2 + yy.^2); % sinc function 0125 z = sin(r)./r; 0126 0127 bcs = csapi( {x,y}, z); % cubic spline (cs) approximate (ap) with interpolation (i) 0128 0129 surf(xx,yy,z) 0130 figure 0131 fnplt(bcs) % plot cubic spline interpolation 0132 0133 0134 0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0136 0137 a = 0:4; 0138 b = -2:2; 0139 0140 R = 4; 0141 r = 2; 0142 0143 V(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].'*[1 1 1 1 1]; 0144 V(2,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[0 1 0 -1 0]; 0145 V(1,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[1 0 -1 0 1]; 0146 0147 dough0 = csape( {a,b}, V, 'periodic'); 0148 0149 figure; fnplt(dough0) 0150 0151 % normals to surface 0152 nx = 43; 0153 xy = [ones(1,nx); linspace(2,-2,nx)]; 0154 points = fnval(dough0,xy)'; 0155 0156 ders = fnval(fndir(dough0, eye(2)), xy); 0157 0158 normals = cross(ders(4:6,:), ders(1:3,:)); 0159 normals = (normals ./ repmat(sqrt(sum(normals .* normals)), 3, 1))'; 0160 pn = [points; points + normals]; 0161 hold on 0162 for j = 1:nx 0163 plot3( pn( [j,j+nx], 1), pn( [j,j+nx], 2), pn( [j,j+nx], 3) ), 0164 end 0165 hold off 0166 0167 % 2D projection 0168 figure; fnplt( fncmb( dough0, [1 0 0; 0 1 0] )) 0169