如果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_width
    for h=1:srcImg_height
        if(srcImg(w,h)<thrshldVal)
          metalImg(w,h) = 0;
          notMetalImg(w,h) = srcImg(w,h);
        else 
           metalImg(w,h) = 255;
           notMetalImg(w,h) = 0;
        end
    end
end

% 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_height
    for w = 2:radon_width-1
          preIndex = 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;
            end
         end
    end
end

% 开始插值，非金属的radon变换线性插值替换金属的radon变换
for h = 1:radon_height
    upIndex = 0;
    for w = 1:radon_width
        if( 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:w
                 left_val = g_pBeta*(w - k)/(w - upIndex);
                 right_val = g_qBeta*(k - upIndex)/(w - upIndex);
                 R_notMetal(k,h) = left_val + right_val;
                end
                 upIndex = 0; 
            else
                upIndex = w;
            end        
        end
    end
end


% 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_width
    for h = 1:outImg_height
        if(metalImg(w,h) > 0)
           irdnImg(w,h) = xmax;
        end
    end
end

%图像的归一化处理
xmin = min(min(irdnImg)); %求得InImg中的最小值
resImg = (irdnImg-xmin)/(xmax-xmin); %归一化
figure,imshow(resImg),title('金属伪影校正最终图像');