Radon变换（Matlb）_iradon变换和radon变换

clear
clc
%只是尽量的大概的模拟出一个radon变换
%I = imread('Mark.bmp');
%I = imread('gantrycrane.png');
 I = imread('2Line.png');
I = rgb2gray(I);
% I(find(I==254))=0;  %2Line.png这个图用的
[m,n] = size(I);
x0 =(floor(m/2)+1);y0 =(floor(n/2)+1);%坐标原点 放中间了 所以p有正负两个值
xptemp = floor((floor(sqrt(m^2+n^2))+1)/2)+1;%p值得大小 就是图像对角线的长度的一半
xp = [-xptemp:1:xptemp];%为了画图用
xp = xp';
P = zeros(2*xptemp+1,180);%留着盛放变换结果
tic;
for t =1:180  %1:180 或者 0:179 为了方便 选择了1:180  
    ptemp = zeros(m,n);%为了刻画狄拉克δ函数
    ctemp = floor([(1-x0):(m-x0)].*cos((pi/180)*(t+90)))+1;%cos x 列 n  盛放了所有x*cos(t)
    stemp = floor([(1-y0):(n-y0)].*sin((pi/180)*(90+t)))+1;%sin y 行 m  盛放了所有x*sin(t)
        for i =1:n
        ptemp(:,i) = ctemp(:)';  %将对应的所有的x*cos(t)放入到ptemp中
        end
        for j =1:m
            ptemp(j,:) = ptemp(j,:)+stemp(:)';
            %将对应的所有的y*sin(t)与已经对应放入的x*cos(t)相加放入到ptemp中
        end
        
        %从这应该可以用积分,但是由于I是个二维矩阵,求积分时会报矩阵维度不一致,暂未找到解决办法。
%  附（不能用）：   syms x y t p;
%         f = @(x,y,t,p) I.*(1.*(x*cos(t)+y*sin(t)-p~=0)+0.*(x*cos(t)+y*sin(t)-p==0));
%         g = @(t,p) integral2(@(x,y)f(x,y,t,p), , , ,);
%  最后代数，在f位置会报错 把I去掉就不会，不知道怎么改，也不确定式子是否有误。
        for p = -xptemp:1:xptemp             
            dtemp = dirac(p - ptemp);
            idx = dtemp ==Inf;dtemp(idx)=1;   
            % 将狄拉克δ函数刻画好后就是一条白色直线(二值图像)
            %因为t角度确定,p距离确定,dtemp和I是一样大的所以对应x,y,即x,y确定
            
            P(p+xptemp+1,t) = sum(sum(dtemp.*double(I)));
            %这用求和应该算是在算积分,最后会变成一个点,和自带radon变换对比貌似应该是
            %如果不求和，选出dtemp.*double(I)中出现次数最多的数(除0外)
            %达到一定阀值应该就是条直线,应该算是变成了用线模板进行线检测了
            
            
            %dtemp和I的点成的和应该是I在dtemp图像中这条线的所有像素值的和
            %radon变换的原理是将线通过线积分变换为点,直线检测是点的累计次数
            %一条线就是一个点,点越多线越清晰,
            %但是百度查到直线检测检测的是二维矩阵的峰值
            %使用findpeaks函数，可以通过限定阈值控制检测峰值的个数
            %得到峰值点,二维矩阵中有x,y,p,t所以能画出直线
            %但是注意坐标原点位置和角度变换，而且画出的为直线不是线段
            %百度有各种优化的和不同计算方式的radon变换,但是看不懂。
        end
end
toc;
t =1:180;
subplot(1,2,1),imagesc(t,xp,P),title('编的');colormap(hot);colorbar
tic;[R,XP] = radon(I);toc;
subplot(1,2,2),imagesc(t,XP,R),title('自带的');colormap(hot);colorbar
%由于没有使用积分，最后的结果全是整数,并且成像不太好。并且程序速度极低,太久了，
%这个程序算是只为了和自带的radon变换的到相似的结果
%radon变换将所有角度加了90度 ，没有想通为什么
%最后最高值比自带的radon变换要大
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

结果**

速度差很多，效果也不太好，图像越大，速度差的越多