function Geometry = dicomrt2geometry(dicomrt_dir)
%DICOMRT2GEOMETRY Import DicomRT dataset to Ryan's "Geometry" struct
% Usage:
%   [ Geometry ] = dicomrt2geometry ( [input_dir] )
% Input:
%   input_dir = directory contains DicomRT data
%   [] = prompt to select input_dir
% Output:
%   Geometry = Ryan's Geometry struct used with RDX TPS
%
% This function imports DicomRT data set and create RDX TPS compatible Geometry 
% struct.
%
% TODO: Switch from filename based Dicom identification to StudyID based.
%
% See also readPinnGeometry.m
%
% Author: Xiaohu Mo

Geometry = [];

if nargin < 1
    dicomrt_dir = uifile('getdir', 'Select the directory contains the DICOM-RT images');
end

if exist(dicomrt_dir, 'dir') ~= 7
    msgbox('Cannot find input directory');
    return;
end

%% -----===== Start Importing =====-----
% set dicom dictionary before any dicominfo or dicomread
dicomdict('factory');

% change path due to the limitation of DicomRT toolbox
% will cd back later
previousPath = pwd;
cd(dicomrt_dir);

% get the list of filenames
dirlist = dir(dicomrt_dir);
filenames = cell(0);
% remove directories
for i = 1:numel(dirlist)
    [tilde, tilde, ext] = fileparts(dirlist(i).name);
    if ~dirlist(i).isdir && strcmpi(ext, '.dcm')
        filenames{end+1} = dirlist(i).name;
    end
end

% get all dicominfo to a cell array the same size as filelist
for i = 1:numel(filenames)
    dcminfo{i} = dicominfo(filenames{i});
end

% TODO integrity check
% UID should match
% should be only one RTSTRUCT file present

%% -----===== Import Dicom CT data set =====-----
% Create DicomRT ct list file
fid = fopen('ctlist.txt', 'w');
for i = 1:numel(dcminfo)
    if strcmpi(dcminfo{i}.Modality, 'CT')
        fprintf(fid, '%s\n', filenames{i});
    end
end
fclose(fid);

[cellCt, ct_xmesh, ct_ymesh, ct_zmesh] = dicomrt_loadct('ctlist.txt');

%% -----===== Set Downsampling flag =====-----
button = questdlg('Do you want to downsample the CT dataset?', 'Downsampling', 'Yes', 'No', 'Yes');
if strcmpi(button, 'yes')
    downsample_flag = true;
else
    downsample_flag = false;
end

%% -----===== Import DicomRT structure =====-----
% File DicomRT structure file
for i = 1:numel(dcminfo)
    if strcmpi(dcminfo{i}.Modality, 'RTSTRUCT')
        rsfilename = fullfile(dicomrt_dir, filenames{i});
    end
end

% prompt if can not find RTSTRUCT file in this directory
if isempty(rsfilename)
    [temp_FileName temp_PathName] = uifile('get', '*.dcm', 'Please select the DicomRT-Struct file');
    rsfilename = fullfile(temp_PathName, temp_FileName);
end

% Import the struct
cellVoi = dicomrt_loadvoi(rsfilename);

% Downsampling
if downsample_flag == true
    if numel(ct_xmesh) < 128
        downsample_factor = 0.5;
        dx = ct_xmesh(2)-ct_xmesh(1);
        dy = ct_ymesh(2)-ct_ymesh(1);
        ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);
        ct_ymesh = linspace(ct_ymesh(1)+dy/2, ct_ymesh(end)-dy/2, numel(ct_ymesh)*downsample_factor);
    else
        downsample_factor = 0.125;
        dx = ct_xmesh(2)-ct_xmesh(1);
        dy = ct_ymesh(2)-ct_ymesh(1);
        ct_xmesh = linspace(ct_xmesh(1)+dx/2, ct_xmesh(end)-dx/2, numel(ct_xmesh)*downsample_factor);
        ct_ymesh = linspace(ct_ymesh(1)+dy/2, ct_ymesh(end)-dy/2, numel(ct_ymesh)*downsample_factor);        
    end
end
ROIS = dicomrt_rasterizeVoi(cellVoi, cellCt, ct_xmesh, ct_ymesh, ct_zmesh);

%% -----===== Done Importing =====-----
% change the path back
cd(previousPath);

%% -----===== Create Geometry struct =====-----
Geometry.ROIS = ROIS;
Geometry.data = cellCt{2};
Geometry.voxel_size = [ct_xmesh(2)-ct_xmesh(1) ...
    ct_ymesh(2)-ct_ymesh(1) ...
    ct_zmesh(2)-ct_zmesh(1)];
Geometry.start = [ct_xmesh(1) ct_ymesh(1) ct_zmesh(1)];

if downsample_flag == true
    Geometry.data = resamplegrid(Geometry.data, [downsample_factor downsample_factor 1]);
    % Geometry.data = resamplegrid(Geometry.data, 0.5);
    % note that ct_xmesh, ct_ymesh, ct_zmesh is already downsampled
    Geometry.start(1:2) = Geometry.start(1:2) + Geometry.voxel_size(1:2) / 4;
end

Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
[Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);