%
% This script can be executed from shell with the command
% matlab -nodisplay -r "howToPrepareData_AT;quit"
%
%% How to call autoSegment function
% init
 clear all
 close all
% clc

%% parameters of the segmentation. These parameters are used and you can use them as they are:
% opt.isLeft = true;
opt.isLeft = false;
opt.doSmoothing = true;
opt.ParallelProcess = true;
opt.numTemplates = 5;
opt.threshold = 190;
opt.volThr = 35000;
opt.ransacMaxIter = 5.0e6;
opt.ransacDist    = 6; %mm
opt.ransacNo      = 4; % number of samples for ransac!
opt.searchWindowSize = [65, 65, 45];

% The folder that contains images to segment
% folderIn = './sampleData/samplesLeft/';
% folderIn = './sampleData/samplesRight/';
% folderIn = '../../../../data/N/0/3/N34-85688/CT-N34-85688-1/'; % should check it exist
% folderIn = '../../../../data/P/2/2/P221-657235/CT-P221-657235-1/';
% folderIn = '../../../../data/N/2/0/N209-149341/CT-N209-149341-1/';
%folderIn = '../../../../data/P/0/0/P421-430069/CT-P421-430069-1/';
folderIn = '../../../../data/P/0/0/P1-490313/CT-P1-490313-1/';
folderOut = [folderIn 'matlab']; % set output in a matlab folder within folderIn
if ~exist(folderOut, 'dir') % create folder if not exist 
    mkdir(folderOut);
end
    
modelPath = '../Model/';

% This line adds sift feature extraction toolbox to the path
 run('../Libraries/SIFT3D-master/build/lib/wrappers/matlab/setupSift3D.m')

% all nenessary functions are stored in the folder fonctions. We add these
% functions to the path
addpath(genpath('./functions'));

%% load model and extract the sift features on the meanShape
display('Load SS models (right & left)...')
tic;
% model = loadModel(modelPath, opt.isLeft);
modelLeft  = loadModel(modelPath, true);  % isLeft = true;
modelRight = loadModel(modelPath, false); % isLeft = false;
toc;

if(folderIn(end)~=filesep)
    folderIn = [folderIn filesep];
end

if(folderOut(end)~=filesep)
    folderOut = [folderOut filesep];
end
%% load samples one by one and perform segmentation and then save the result
% This example assumes that all images are listed in a single folder. or in
% folders inside one folder specified in folderIn. For example, let's say
% folderIn is structured as follows:
% folderIn
%           Img1.am
%           Img2
%                  Img2.mhd
%                  Img2.raw
%           Img3
%               Img3_0001.dcm
%               Img3_0002.dcm
%               ...
%               Img3_0250.dcm
%           Img4.am
%           ...

% If the structure of your folders are different, then you should adapt
% this folder


% list of all images in the folder
bones_in = dir([folderIn '*']);

%Exclude "." and ".." from the list. Put dicom images in a folder!
bones_in = bones_in(3:end);

% Restrict to folder named "dicom", for use in /home/shoulder/data/
bones_in_idx = find(strcmp({bones_in.name},'dicom') == 1);
bones_in = bones_in(bones_in_idx);
% we should check it exists

disp(['folderIn: ' folderIn]);
disp(['folderOut: ' folderOut]);
disp(bones_in);



for b = 1 : size(bones_in, 1)
disp('Load images...');
tic;
    boneName = [folderIn bones_in(b).name];
    
    %if the folder contains DICOMDIR or each image is in one directory!
    if (isdir(boneName)) 
        imageFile = dir(boneName);
        imageFile = imageFile(3:end);
        for file_counter = 1:size(imageFile, 1)
            if(strcmp(imageFile(file_counter).name(end-2:end), 'mhd') || strcmp(imageFile(file_counter).name(end-2:end), 'dcm'))
                boneName = [boneName, filesep imageFile(file_counter).name];
                break;
            end
        end
    end
    
    fname = boneName;
    ind = strfind(fname, '.');
    if(isempty(ind))
        continue;
    end
    ind = ind(end);
    fileFormat = fname(ind+1:end);
    % right now I have code for ".am", ".mhd" and ".dcm" file formats. You
    % can add more cases here to read data.
    switch fileFormat 
        case 'am'
            try
                [imgInfo, Img] = LoadData_Amira(boneName);
            catch
                display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']);
            end
        case 'mhd'
            try
                [imgMHD, info]=read_mhd(boneName);
            catch
                display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']);
            end
            eval(['Img = ' class(imgMHD.data) '(zeros(size(imgMHD.data, 2), size(imgMHD.data, 1), size(imgMHD.data, 3)));']);
            for ii = 1:size(imgMHD.data, 3)
                Img(:, :, ii) = imgMHD.data(:, :, ii)';
            end
            imgInfo.BoundingBox.xMin = imgMHD.origin(1);
            imgInfo.BoundingBox.yMin = imgMHD.origin(2);
            imgInfo.BoundingBox.zMin = imgMHD.origin(3);
            imgInfo.PixelSpacing(1)  = imgMHD.spacing(1);
            imgInfo.PixelSpacing(2)  = imgMHD.spacing(2);
            imgInfo.SliceThickness   = imgMHD.spacing(3);
            imgInfo.xSize = imgMHD.size(1);
            imgInfo.ySize = imgMHD.size(2);
            imgInfo.zSize = imgMHD.size(3);
            imgInfo.BoundingBox.xMax = imgMHD.origin(1) + (imgMHD.size(1)-1) * imgInfo.PixelSpacing(1);
            imgInfo.BoundingBox.yMax = imgMHD.origin(2) + (imgMHD.size(2)-1) * imgInfo.PixelSpacing(2);
            imgInfo.BoundingBox.zMax = imgMHD.origin(3) + (imgMHD.size(3)-1) * imgInfo.SliceThickness;
            
            if imgInfo.BoundingBox.zMax - imgInfo.BoundingBox.zMin > 700
                display ([ bones_in(b).name ': skipped processing, due to the large size of scan, please crop this image and then try again (' num2str(size(imgMHD.data)) ')']);
                continue;
            end
            clear imgMHD;
        case 'dcm'
            dicomDir = strsplit(boneName, filesep);
            dicomDir = strjoin(dicomDir(1, 1:end-1), filesep);
            [Img, imgInfo] = readDicomStack(dicomDir);
            if (isempty(Img))
                display ([ bones_in(b).name ': Could not read the dicom images.']);
                continue;
            end
        otherwise   
            continue;
    end
    fname =  bones_in(b).name;
    ind = strfind(fname, '.');
    if(isempty(ind))
        ind = length(fname)+1;
    else
        ind = ind(end);
    end
toc;

% Identify scapula side
disp('Identify scapula side');
tic;
[Left, Right] = IdentifyRightAndLeftScapula(Img, imgInfo, modelLeft, modelRight); % Img and imgInfo is similar to autoSegment function.
toc;

% Segment scapula
disp('Segment CT...');
tic;

if(isempty(Right) && ~isempty(Left))
    disp('Image contains only the left scapula.');
    opt.isLeft = true;
    finalSurface = autoSegment(Img, imgInfo, modelLeft, opt);
 end
 if(isempty(Left) && ~isempty(Right))
    disp('Image contains only the right scapula.');
    opt.isLeft = false;
    finalSurface = autoSegment(Img, imgInfo, modelRight, opt);
 end
 if (~isempty(Right) && ~isempty(Left))
    disp('Image contains both left and right scapulae.');
    opt.isLeft = false;
    finalSurface_Right = autoSegment(Right.img, Right.imgInfo, modelRight, opt);
    opt.isLeft = true;
    finalSurface_Left = autoSegment(Left.img, Left.imgInfo, modelLeft, opt);
 end
 if (isempty(Right) && isempty(Left))
    disp('Could not recognize, which scapulae it is.');
 end

toc;

% Write surface file
disp('Write surface file...');
tic;
    simplePlyWrite(finalSurface, [folderOut fname(1:ind-1) '_auto.ply']);
toc;

disp('Measure anatomy...');
tic;
    measurements = doMeasurements(finalSurface);
%    save( [folderOut fname(1:ind-1) '.mat'],'measurements');
    save( [folderOut 'anatData.mat'],'measurements');
toc;
    
    disp(['Processing of sample ' fname(1:ind-1) ' is finished.'])
end