怎样用样条插值进行层间插值？
计算出来的值全是0。求指点！x0=dicomread('C:\Users\Felix\Desktop\t2_haste_tra_p2_mbh_320_14\IM-0012-0007.dcm');
y0=dicomread('C:\Users\Felix\Desktop\t2_haste_tra_p2_mbh_320_14\IM-0012-0008.dcm');
x=dicomread('C:\Users\Felix\Desktop\IM-0012-1007.dcm');
%x是所求图像
n = length(x0); = size(x);
y = zeros(m1, m2);
p(n) = 0;q(1) = 0;d(1) = 0;d(n) = 0;
for k = 2:n-1,
h(1) = x0(k) - x0(k-1);h(2) = x0(k+1) - x0(k);
p(k) = h(1) / (h(1) + h(2));
q(k) = h(2) / (h(1) + h(2));
d(k) = 6*((y0(k+1) - y0(k))/h(2) - (y0(k) - y0(k-1)) / h(1)) / (h(1) + h (2));
end
b(1) = q(1) / 2;
for k = 2:n-1,
b(k) = q(k) / (2 - p(k)*b(k-1));
end
temp(1) = d(1) / 2;
for k = 2:n,
temp(k) = (d(k) - p(k)*temp(k-1)) / (2 - p(k)*b(k-1));
end
M(n) = temp(n);
for k = n-1:-1:1,
M(k) = temp(k) - b(k)*M(k+1);
end
% 三次自然边界样条插值
for i = 1:m1*m2
for k = 1:n-1,
if (x(i) >= x0(k)) && (x(i) <= x0(k+1)),
h = x0(k+1) - x0(k);
y(i) = M(k)*(x0(k+1) - x(i))^3 / (6*h)+ M(k+1)*(x(i) - x0(k))^3 / (6*h)+ (y0(k) - M(k)*h^2 / 6)*((x0(k+1) - x(i))/h)+ (y0(k+1) - M(k+1)*h^2 / 6)*((x(i) - x0(k))/h);
end
end
end
dicomwrite(y,'C:\Users\Felix\Desktop\IM-0012-1007.dcm');
T=imadjust(x0);
T1=imadjust(y0);
F=imadjust(y);
subplot(131),imshow(T);
subplot(132),imshow(F);
subplot(133),imshow(T1);
页:
[1]