%% Optical CT processing 

% Created by Jonathan Bumstead
% Date: 15 March 2018
% Full instructable: https://www.instructables.com/id/Desktop-CT-and-3D-Scanner-With-Arduino/
% Contact: jonathan.bumstead at gmail.com

% Code for processing images collected with desktop CT scanner. Images are
% uploaded, cropped, and downsampled. Then using the inverse Radon
% transform, the volume of transparent objects is constructed. Finally,
% the program displays volumetric data in several ways.

%% Load cropped and downsampled datasets

N=200; % number of images in scan

foldIn = 'C:\Users\Jon\Downloads\tomato\';
clear Z
for p=1:N
    filename = [foldIn,'I',num2str(p),'.jpg'];
    I = imread(filename);
    I = rgb2gray(I);
    I = imcomplement(I);
    Z(:,:,p)= I;
end

Z = double(Z);
[nPy,nPx,~]=size(Z); 

%% Watch video of the scan (optional)
close all;

Ma = max(Z(:));
Mi = min(Z(:));

figure;
for d=1:N
    
    imagesc((Z(:,:,d)),[Mi Ma]); axis image;
    pause(0.01);
    clf
    
end

%% Watch sinogram (optional)

for z=1:nPy
    
    S=squeeze(Z(z,:,:)); % reshape data into sinogram
    imagesc(S,[Mi Ma]); axis image; % display sinogram
    pause(0.01);
    clf
    
end

%% Inverse Radon transform with filtered backprojection

theta = linspace(0,360,N); % array of angles from scan

clear FB
for z=1:nPy
    
    S=squeeze(Z(z,:,:)); % reshape data into sinogram
    FB(:,:,z) = iradon(S,theta); % filtered backprojection
    
end

%% Display results 1: cross sections

M=max(FB(:));
Mi=min(FB(:));

figure;
for c =1:nPy
    
    imagesc(FB(:,:,c),[0 M*0.1]); axis image off;
    pause(0.01);
    clf;
    
end

%% Display results 2: 3D rendering

[nFy,nFx,nFz]=size(FB);
FBt =FB(:,:,:);
factorDown = 2;
ny=round(nFy/factorDown); nx=round(nFy/factorDown); nz=round(nFz/factorDown); %% desired output dimensions

[y, x, z]=ndgrid(linspace(1,size(FBt,1),ny),linspace(1,size(FBt,2),nx),linspace(1,size(FBt,3),nz));
FBdown=interp3(FBt,x,y,z);
FBdown(FBdown<0.3)=0;
FBdown=FBdown(:,:,nz:-1:1);
FBdown = FBdown./max(FBdown(:))*255;

figure;
vol3d('cdata', FBdown, 'xdata', [0 1], 'ydata', [0 1], 'zdata', [0 1]);
alphamap([0 linspace(0.1, 0, 256)]);
colormap('jet');
axis equal off
set(gcf, 'color', 'w');
view(3);

%% Display results 3: 3D rendering 2

intensityThreshold = 15;

figure;
p = patch(isosurface(x,y,z,FBdown,intensityThreshold));
isonormals(x,y,z,FBdown,p)
p.EdgeColor = 'none';
p.FaceColor = 'red';
p.Marker = 'none';
p.FaceAlpha = 0.2;

daspect([1 1 1])
view(3); 

camlight 
lighting gouraud

set(gcf, 'color', 'w');
set(gca,'CameraViewAngleMode','Manual','XTick',0:100:600,'YTick',0:100:600,'ZTick',0:100:600,'XTickLabel',' ','YTickLabel',' ','ZTickLabel',' ')
% xlim([-0.5 600]);
% ylim([-0.5 600]);
% zlim([-0.5 600]);
axis tight equal 
grid on;

figure;
[faces,verts,colors] = isosurface(x,y,z,FBdown,intensityThreshold,y-max(y(:))/2); 
patch('Vertices',verts,'Faces',faces,'FaceVertexCData',colors,...
    'FaceColor','interp','EdgeColor','interp')
view(30,-15)
axis vis3d
colormap jet