sino图的两种投影方案

sino图的两种投影方案

1.旋转坐标
2.旋转图像

准备阶段

close all; clear; clc;
%% Parameters
N = 512;
views = 360;
bins = 512;
dv = 2 * pi / views;
x1 = 0; y1 = 0; r1 = 200; mu1 = 0.5;
x2 = 80; y2 = 80; r2 = 30; mu2 = 1;
%% Generate model
vec = -N/2 + 0.5: N/2 - 0.5;
[xx, yy] = meshgrid(vec, -vec);
img = zeros(N, N);
img((xx - x1).^2 + (yy - y1).^2 <= r1^2) = mu1;
img((xx - x2).^2 + (yy - y2).^2 <= r2^2) = mu2;
figure('name','cirle in circle'); imshow(img, []); title('Original image');

旋转坐标

%% projection1
dn = 0.5;
inte_vec = -N/2 + dn/2: dn: N/2 - dn/2;
[ox, oy] = meshgrid(vec, -inte_vec);
alpha = dv/2: dv: 2*pi - dv/2;
proj = zeros(views, bins);
for i = 1: views
    theta = alpha(i);
    X = ox .* cos(theta) - oy .* sin(theta);
    Y = ox .* sin(theta) + oy .* cos(theta);
    temp = interp2(xx, yy, img, X, Y, 'linear', 0);
    proj(i, :) = sum(temp);
end
figure('name','projection1');
imshow(proj, []);
title('Projection image by rotate axis');

旋转图像

我们实现把图像保存为’cirInCir.bmp’

%% projection2
img = imread('cirInCir.bmp');%读图
gImg = rgb2gray(img); %这里用的是灰度图
B = zeros(360,1);
for i = 1:360
    rImg = imrotate(gImg,-(i-1));
    [w,l] = size(rImg);
    B(i) = length(round(sum(rImg)));
end
[j,k] = max(B);
A = zeros(360,j);
for m = 1:360
    rImg = imrotate(gImg,-(m-1));
    [w,l] = size(rImg);
    C = round(sum(rImg)/l);
    k = round((j-B(m))/2);
    for n = 1:B(m)
        A(m,k+n) = C(n);        
    end
end
figure('name','prjection2');
imshow(A,[])
title('Projection image by rotate picture');