[MATLAB图像处理] 哪位大侠能帮小弟改一下这个有关三次样条 插值多项式的程序

[a57649314]

w=imread('11.jpg');
[x,y]=size(w);
n = size(x,2); % 插值点个数, 注意: 不是书上的 n
for j = 1 : n-1
h(j) = x(j+1) - x(j); % 步长
f(j) = (y(j+1) - y(j)) / h(j); % 一阶差商
end
for j = 2 : n-1
lambda(j) = h(j)/(h(j-1) + h(j));
mu(j) = h(j-1)/(h(j-1) + h(j));
g(j) = 3 * (lambda(j) * f(j-1) + mu(j) * f(j));
end

A = 2 * eye(n,n);
for i = 2 : n-1
A(i,i+1) = mu(i);
A(i,i-1) = lambda(i);
end

b = zeros(n,1);
for i = 2 : n-1
b(i) = g(i);
end

% 2 阶导数边界条件
%b(1) = 3*f(1);
%A(1,2) = 1;
%b(n) = 3*f(n-1);
%A(n,n-1) = 1;

m = A \ b; % 求解线性方程组, 得到各节点的一阶导数值
S = [];
for j = 1 : n-1
s = (y(j)/(h(j)^3))*conv([1 -2*x(j+1) x(j+1)^2],[2 h(j)-2*x(j)]) + ...
(y(j+1)/(h(j)^3))*conv([1 -2*x(j) x(j)^2],[-2 h(j)+2*x(j+1)]) + ...
(m(j)/(h(j)^2))*conv([1 -2*x(j+1) x(j+1)^2],[1 -x(j)]) + ...
(m(j+1)/(h(j)^2))*conv([1 -2*x(j) x(j)^2],[1 -x(j+1)]);
S = [S;s];
end