基于matlab的医学成像技术滤波反投影仿真,包括直接反投影,S-L滤波,R-L滤波,Lewitt滤波

147 阅读1分钟

1.算法描述

       医学成像技术滤波反投影  含R-L滤波,  R-S滤波,Lewitt滤波  重建后图像清晰。对人脑体模仿真数据与实际的CT图像数据进行滤波反投影重建实验,比较直接反投影与R——L和S——L滤波方法的重建结果。

 

2.仿真效果预览

matlab2022a仿真结果如下:

 

1.png

 

3.MATLAB部分代码预览 `function produce(filter)

        for m=1:180/delta;

            creat_nProj(m);%补零,用于卷积

            

            c_h1=conv(nProj,h1);%求卷积

            c_h2=conv(nProj,h2);

            c_h3=conv(nProj,h3);

            c=0.5*(N-1)(1-cos(mdelta)-sin(m*delta));

            for i=1:N;

                for j=1:N;

                    L(i,j)=c+(i-1)cos(mdelta)+(j-1)sin(mdelta); %内插

                    n = fix(L(i,j));

                    cL=L(i,j)-n;

                    if strcmp(filter,'ori')

                       if (n>0)&(n<255)

                            rProj1(i,j)=rProj1(i,j)+(1-cL)nProj(n+N)+cLnProj(n+1+N);

                        elseif n==255

                            rProj1(i,j) =rProj1(i,j)+nProj(n+N);

                        elseif n==0

                            rProj1(i,j) =rProj1(i,j)+nProj(n+N+1);

                       end

                    elseif strcmp(filter,'sl')

                        if (n>0)&(n<255)

                            rProj2(i,j)=rProj2(i,j)+(1 - cL)c_h1(n+N+N-1)+cLc_h1(n+N+N);

                        elseif n==255

                            rProj2(i,j) =rProj2(i,j)+c_h1(n+N+N-1);

                        elseif n==0

                            rProj2(i,j)=rProj2(i,j)+c_h1(n+N+N);

                       end

                    elseif strcmp(filter,'rl')

                       if (n>0)&(n<255)

                            rProj3(i,j)=rProj3(i,j)+(1 - cL)c_h2(n+N+N-1)+cLc_h2(n+N+N);

                        elseif n==255

                            rProj3(i,j)=rProj3(i,j)+c_h2(n+N+N-1);

                        elseif n==0

                            rProj3(i,j) =rProj3(i,j)+c_h2(n+N+N);

                       end

                    elseif strcmp(filter,'lew')

                        if (n>0)&(n<255)

                            rProj4(i,j)=rProj4(i,j)+(1 - cL)c_h3(n+N+N-1)+cLc_h3(n+N+N);

                        elseif n==255

                            rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N-1);

                        elseif n==0

                            rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N);

                        end

                    end

                end

            end

        end

    end

A_021`