-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain_v2.m
78 lines (63 loc) · 2.15 KB
/
main_v2.m
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
Nx = 259; % number of grid points in x
Ny = 100; % number of grid points in y
Nr = 279; % number of detector elements
Nphi = 180; % number of projection angles
muest = 200;
xLen = Nx * muest; % side of reconstructed rectangle
yLen = Ny * muest; % side of reconstructed rectangle
phiLen = pi; % angular interval
rLen = Nr * muest; % length of detector
phiVec = (0 : Nphi - 1) * phiLen / Nphi;
rVec = (0.5 : Nr - 0.5) * rLen / Nr - rLen / 2;
xVec = (0.5 : Nx - 0.5) * xLen / Nx - xLen / 2;
yVec = (0.5 : Ny - 0.5) * yLen / Ny - yLen / 2;
figure(1); colormap gray;
Nslices = 496;
Nproj = 180;
m = zeros(496,290,180);
mouse3D = zeros(Ny,Nx,Nslices);
p = zeros(290,180);
k = 1;
% Loading multiimage into memory
for i = 1:Nproj
mouse = imread('./data/mouse.tif',i);
m(:,:,i) = mouse(:,:);
end
for i = 201:201
% Resampling to create Sinograms
p(:,:) = m(i,:,:);
% Filtered Backprojecting Sinograms
q = rampfilter(p, rVec, 'signal');
f = backproject_linear(q, rVec, phiVec, xVec, yVec, 'linear');
subplot(1,2,1);imagesc(p); title(['Sinogram (slice ', int2str(i),')']);pause(0);
subplot(1,2,2);imagesc(yVec, xVec, f);title(['Reconstruction (Sinogram ', int2str(i),')'])
pause(0);
% Stacking 2Ds to build the 3D reconstruction
mouse3D(:,:,k) = f;
% Saving subplot(1,2,2) for video generation
M(k) = getframe;
k = k+1;
end
% Visualizing 3D reconstruction
%load mri;
%img = squeeze(mouse3D);
%figure(2);
%imshow3D(img, [], 1);
%pause(0);
% Saving video
v = VideoWriter('mouse_3D_reconstruction.avi');
v.FrameRate = 60;
open(v);
writeVideo(v,M);
close(v);
% Re-normalization and saving in a .tif
mn=min(min(min(mouse3D)))
mx=max(max(max(mouse3D)))
B_P_P = uint8(((mouse3D)-mn).*(255/(mx-mn)));
for k = 1:Nslices
if(k==1)
imwrite(B_P_P(:,:,k),'R.tif','WriteMode','overwrite');
else
imwrite(B_P_P(:,:,k),'R.tif','WriteMode','append');
end
end