close all ; 
clear all ; 
clc; 
%%author - Praga
%%process time - 13 sec (without thresholding) 
%%             - 1 minute (with thresolding)
%%%%%%%%  Task 1  %%%%%%%%%%%%%%%%%%%%%%%
%%the conditions for R %%%%%%%%%%%%%%%%%%

%% 1) tan(theta/2) > R/Yp
%% 2) Yus*cos(theta/2) > 2R
%% 3) Xus*cos(theta/2) > 2R
%% 4) Yus - Yp > R
%% 5) Xus > 2R
%%

%%%%%%% Task 2 %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%defining all the variables%%%%%%%%

FOV=120;             %degree
theta=(pi*2)/180;    %radian 
R=2;                 %cm
Yp=6;                %cm
Xus=10;              %cm
Yus=10;              %cm
Sx=0.1;              %cm
Sy=0.1;              %cm
num_pixel=Xus/Sx;  %done only for X as both are equal

%%to get the main image over here
[im_main]=createmainimage(Xus,R,theta,num_pixel);

%%%%%%%%% Task 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%creation and visualization of the empty voxel%%%
VoxelMat=zeros(num_pixel,num_pixel,num_pixel);

for i=1:num_pixel
    for j=1:num_pixel
        for k=1:num_pixel
            if (i-num_pixel/2)^2+(j-num_pixel/2)^2+(k-num_pixel/2)^2<(R/Sx)^2
                VoxelMat(i,j,k)=255;
            end
            
        end
    end
end

%view of voxels
[vol_handle,FV]=VoxelPlotter(VoxelMat,1); 
view(3);
daspect([1,1,1]);
set(gca,'xlim',[0 num_pixel], 'ylim',[0 num_pixel], 'zlim',[0 num_pixel]);
title('actual sphere in Voxels without the reconstruction algorithem');
xlabel('X');ylabel('Y');zlabel('Z');

%push of taken image back into a 3D array as reconstruction 
%reconstruction view will be differnt as it makes coding easy
%rotating the reconstructed matrix by 90 degrees will give the original
%will be working with X and Z axis 

%%nearest neighbour voxel mapping

array_recon=zeros(num_pixel,num_pixel,num_pixel);
angle_1=-pi/3:theta:pi/3; %in radians

%special case handled first .......
%array_recon(:,:,num_pixel/2)=im_main(:,:,31);

%delete the angle zero from array 
%angle_1= angle_1([1:30, 32:end]);

for k=1:size(angle_1,2)
    %rotation matrix with respect to y axis 
    %Rot = [cos(angle(k)) 0 sin(angle(k));0 1 0;-sin(angle(k)) 0 cos(angle(k))];
    x=size(im_main,1);
    y=size(im_main,2);
    
    for i =1:1:x
        for j=1:1:y
            new_cord=rotx(angle_1(k)*180/pi)*[i,j,0]';
            new_cord(3)=new_cord(3)+50;
            new_cord=ceil(new_cord);
            if new_cord(3)>0
                array_recon(new_cord(1),new_cord(2),new_cord(3))=im_main(i,j,k);
                array_recon=uint8(array_recon);
            end
            
        end
    end    
end

%%print voxel after nearest neighbour map
figure;
[array_recon_1]=VoxelPlotter(array_recon,1); 
view(3);
daspect([1,1,1]);
set(gca,'xlim',[0 num_pixel], 'ylim',[0 num_pixel], 'zlim',[0 num_pixel]);
title('Voxels without the interpolation');
xlabel('X');ylabel('Y');zlabel('Z');

%%%%%bilinear interpolation %%%%%%%%%%%%
zoomed=zeros(size(array_recon));
%%%bilinear interpolation  
for j=1:size(array_recon,3)
    zoomed(:,:,j)=bilinear( array_recon(:,:,j) );
end

figure;
[zoomed_1]=VoxelPlotter(zoomed,1); 
view(3);
daspect([1,1,1]);
set(gca,'xlim',[0 num_pixel], 'ylim',[0 num_pixel], 'zlim',[0 num_pixel]);
title('Voxels after bilinear interpolation');
xlabel('X');ylabel('Y');zlabel('Z');

%%%%%%%safe images after reconstruction %%%%%%%%%
mkdir('final3Dimagesafterrecon');
cd('final3Dimagesafterrecon');
%save the image
image_index=1;

for i=1:100
    filename = 'p%d.bmp';
    filename = sprintf(filename,image_index);
    image_index = image_index + 1;
    imwrite(zoomed(:,:,i),filename);
end