如果CT机性能较差，图像中的金属伪影信息就无法消除，现在高级的CT机都带有金属伪影去除功能。
本文提出一种最简单的金属伪影的消除方法，对单金属图像进行伪影消除的效果还不错，简单的多金属图像进行伪影消除的效果还将就，复杂的多金属图像的伪影消除效果就不太好。伪影的消除依赖于线性插值方法，线性插值方法有需要金属的radon变换结果，多金属的randon变化就是多条正弦线的集合。
方法的原理参考文献：CT图像中金属伪影的快速校正（谷建伟）。
去除伪影的步骤如下：
1、对原始图像（图A）进行阈值分割，把金属图像分割出来，金属图像（图B），抠掉金属的图像（图C）。
%显示原始图像的灰度图srcImg=imread('image7.png');srcImg=rgb2gray(srcImg);srcImg=imresize(srcImg,[400,400]);figure;imshow(srcImg),title('原始图像');%选定阈值进行分割，选择阈值为220[srcImg_width,srcImg_height]=size(srcImg);metalImg = zeros(srcImg_width,srcImg_height);notMetalImg = zeros(srcImg_width,srcImg_height);thrshldVal = 220;for w=1:srcImg_widthfor h=1:srcImg_heightif(srcImg(w,h)<thrshldVal)metalImg(w,h) = 0;notMetalImg(w,h) = srcImg(w,h);else metalImg(w,h) = 255;notMetalImg(w,h) = 0;endendend% xmax = max(max(metalImg)); %求得InImg中的最大值% xmin = min(min(metalImg)); %求得InImg中的最小值% imageB = (metalImg-xmin)/(xmax-xmin); %归一化% figure;imshow(imageB),title('金属的阈值分割图像');% xmax = max(max(notMetalImg)); %求得InImg中的最大值% xmin = min(min(notMetalImg)); %求得InImg中的最小值% imageC = (notMetalImg-xmin)/(xmax-xmin); %归一化% figure;imshow(imageC),title('原始图像去除金属部分');
2、求图像A、图像B和图像C的Radon变换，得到三张正弦图，分别是图RA，图RB，图RC。
%Radon变换的角度范围radon_theta=0:179;% 金属的Radon变换[R_metal,X_metal]=radon(metalImg,radon_theta);figure,imagesc(radon_theta,X_metal,R_metal),title('金属的Radon变换');xlabel('\theta');ylabel('x\prime');colormap(gray);colorbar;%非金属部分的Radon变换[R_notMetal,X_notMetal]=radon(notMetalImg,radon_theta);figure,imagesc(radon_theta,X_notMetal,R_notMetal),title('非金属部分的Radon变换');xlabel('\theta');ylabel('x\prime');colormap(gray);colorbar;
3、由于图RB就是造成图A伪影的原因。在图RC中去除和图RB相交的部分，此时会留下间隙，间隙使用插值算法进行校正。此时得到图R_RES。
% 对金属图像的Radon变换结果进行插值校正[radon_width,radon_height]=size(R_metal);%radon变换的结果是横向坐标是距离，纵坐标是角度%metalPQBeta矩阵,存储金属的radon变换的边界metalPQBeta = zeros(radon_width,radon_height);%获取金属的Radon变换的边界for h = 1:radon_heightfor w = 2:radon_width-1preIndex = w -1;nextIndex = w + 1;currPix = R_metal(w,h);if(currPix > 0)if((R_metal(nextIndex,h) > 0 && R_metal(preIndex,h) == 0) || (R_metal(nextIndex,h) == 0 && R_metal(preIndex,h) > 0))metalPQBeta(w,h) = 100;endendendend% 开始插值，非金属的radon变换线性插值替换金属的radon变换for h = 1:radon_heightupIndex = 0;for w = 1:radon_widthif( metalPQBeta(w,h) == 100 && upIndex == 0)upIndex = w;elseif(metalPQBeta(w,h) == 100 && upIndex ~= 0)g_pBeta = R_notMetal(upIndex,h);g_qBeta = R_notMetal(w,h);midPixValue = R_metal(round((upIndex + w)*0.5),h); if(midPixValue > 0)for k = upIndex:wleft_val = g_pBeta*(w - k)/(w - upIndex);right_val = g_qBeta*(k - upIndex)/(w - upIndex);R_notMetal(k,h) = left_val + right_val;endupIndex = 0; elseupIndex = w;end endendend% figure,imagesc(radon_theta,X_notMetal, R_notMetal),title('对非金属部分radon变换结果进行插值校正');% xlabel('\theta');% ylabel('x\prime');% colormap(gray);% colorbar;
4、对图R_RES进行Radon的逆变换，得到图D，由于图D中此时缺少金属的信息，因此把图B中的金属部分加到图像D中，得到最终的金属伪影校正图F。
%对非金属的radon的变换的插值校正结果进行反变换irdnImg=iradon(R_notMetal,radon_theta,'linear','Hann',1,1);%金属信息补偿xmax = max(max(irdnImg)); %求得InImg中的最大值irdnImg=imresize(irdnImg,[400,400]);[outImg_width,outImg_height]=size(irdnImg);for w = 1:outImg_widthfor h = 1:outImg_heightif(metalImg(w,h) > 0)irdnImg(w,h) = xmax;endendend%图像的归一化处理xmin = min(min(irdnImg)); %求得InImg中的最小值resImg = (irdnImg-xmin)/(xmax-xmin); %归一化figure,imshow(resImg),title('金属伪影校正最终图像');
如果觉得《CT图像去除金属伪影-MATLAB实现》对你有帮助，请点赞、收藏，并留下你的观点哦！
