Labfans是一个针对大学生、工程师和科研工作者的技术社区。 | 论坛首页 | 联系我们(Contact Us) |
![]() |
![]() |
#1 |
初级会员
注册日期: 2013-10-08
帖子: 2
声望力: 0 ![]() |
![]()
第一次做,找不到头绪!
计算出来的值全是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); |
![]() |
![]() |
![]() |
#2 |
初级会员
注册日期: 2013-12-23
住址: 南昌
年龄: 35
帖子: 7
声望力: 0 ![]() |
![]()
你的spline呢?你根本就没进行样条插值呀!!!
|
![]() |
![]() |