%% 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.

%% File information

mainFold = 'G:\PHOTO\3D SCANS\180305orange2\'; % file location of images

imStart=206; % number of first image saved in scan (e.g. DSC_0206.JPG)
N=200; % number of images in scan
imArr =imStart:(imStart+N); % array of image numbers

%% Load test image to determine crop

picNum = num2str(imStart);
switch length(picNum)
    
    case 1
        
        picNum = ['000',picNum];
        
    case 2
        picNum = ['00',picNum];
        
    case 3
        picNum = ['0',picNum];
        
    case 4
        picNum = picNum;
        
end

I=imread([mainFold,'DSC_',picNum,'.JPG']); % load first image

imagesc(I); hold on; axis image; % display image
title('Select top-left and bottom-right for crop');

[Xcrop1, Ycrop1]=ginput(1); % select top-left corner of image
plot(Xcrop1, Ycrop1,'or');

[Xcrop2, Ycrop2]=ginput(1); % select bottom-right corner of image

Xcrop1 =round(Xcrop1);
Xcrop2 =round(Xcrop2);
Ycrop1 =round(Ycrop1);
Ycrop2 =round(Ycrop2);

clf;
imagesc(I(Ycrop1:Ycrop2,Xcrop1:Xcrop2,:)); axis image; % image after cropping

%% Load images and store them in an array

clear Z

for c =1:N

    picNum = num2str(imArr(c)); % picture number
    lenFrame = length(picNum);
    
    switch lenFrame
        
        case 1
            
            picNum = ['000',picNum];
            
        case 2
            picNum = ['00',picNum];
            
        case 3
            picNum = ['0',picNum];
            
        case 4
            picNum = picNum;
            
    end
    
    I=imread([mainFold,'DSC_',picNum,'.JPG']); % read in image  
    I = I(Ycrop1:Ycrop2,Xcrop1:Xcrop2,:); % crop image
    I = rgb2gray(I); % convert to grayscale
    I = imcomplement(I); % invert image intensities. 0 now corresponds to no light attenuation and 255 corresponds to light completed absorped by object
    I = imresize(I,0.3); % downsample image
    
    Z(:,:,c)=I; % store image in array
    
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