马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
计算出来的值全是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); [m1,m2] = 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);
|