WiscPlan_v2/WiscPlanPhotonkV125/matlab_frontend/load2geometry.m

function [Geometry] = load2geometry
    
    %% -----===== Get CT files =====-----
    
    load('WiscPlan_preferences.mat')
    if ~isfield(WiscPlan_preferences,'inDataPath')
        Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
    elseif WiscPlan_preferences.inDataPath ~= 0
        Geometry.data_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the directory with patient data');
    else
        Geometry.data_dir = uigetdir('C:', 'Select the directory with patient data');
    end
    WiscPlan_preferences.inDataPath = Geometry.data_dir;
    thisDir = mfilename('fullpath');
    idcs   = strfind(thisDir,'\');
    prefsdir = thisDir(1:idcs(end)-1);
    save([prefsdir '\WiscPlan_preferences.mat'], 'WiscPlan_preferences')
    
%     Geometry.data_dir = uigetdir('C:\Box Sync\Optimal Stopping\Data\CDP', 'Select folder with patient data');
%     Geometry.data_dir = uigetdir('\\Mpufs5\data_wnx1\_Data\QTBI-clinical\Medivation\0391811\T2', 'Select folder with patient data');
    

    [IMG_in, IMG_path, filterIdx] = uigetfile([{'*.nrrd'; '*.am'}], 'Select CT image', Geometry.data_dir);
%     'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\'
    
    total_num_steps = 4;
    hWaitbar = waitbar(0);
    
    % --  geometry explained --
    % geometry_file important parts:
    % Geometry.ROIS - actually only to get the binary mask of the target
    % Geometry.data - CT image in HU
    % Geometry.voxel_size - size of voxels in cm
    % Geometry.start - image coordinate origin
    
    % these are derived from above:
    % +Geometry.rhomw - CT image with values in relative water density 
    % +Geometry.Smw - very similair to rhomw
    % +Geometry.Fmw2
    % +Geometry.BTV 
    % +Geometry.ring - a ring around PTV used to pull down IMRT optimization
    % --------------------------
    
    switch filterIdx
        case 0
            warning('No file selected, aborting!')
            return 
        case 1
            [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
            
            % offset the CT, because that is what WiscPlan wants
            Geometry.data = Geometry.data + 1000;
            
            disp('NRRD file loaded!')
            Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
            Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
            [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
            
            disp(['current voxel size: ' num2str(meta.spacedirections(1)) ' ' ...
                num2str(meta.spacedirections(2)) ' ' num2str(meta.spacedirections(3)) ' cm.'])
            
            x = input('Is this correct? (y/n)', 's')
            switch x
                case 'n'
                    rescale_f = input('specify resampling factor')
                case 'y'
                    rescale_f = 1;
            end
            
            Geometry.voxel_size = meta.spacedirections * rescale_f;
            Geometry.start=meta.spaceorigin;
        case 2
            imgIn=am2mat([IMG_path, IMG_in]);
            Geometry.data = permute(imgIn.data, [2,1,3]);
            % offset the CT, because that is what WiscPlan wants
            Geometry.data = Geometry.data + 1000;
            
            disp('NRRD file loaded!')
            Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
            Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
            [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
            
            Geometry.voxel_size = imgIn.voxel_size;
            Geometry.start=imgIn.start;
        case 3
            % something selected
            error('please do not select *all* in file selection')
        otherwise
%             error('You should never see this.')
    end
    
    
    % Geometry.ROIS
    
    %% -----===== Get mask/contour files =====-----
    [TRGT_in, TRGT_path, filterIdx2] = uigetfile([{'*.nrrd'; '*.am'}], 'Select file with Target', Geometry.data_dir);
    
    switch filterIdx2
        case 0
            warning('No file selected!')
        case 1
            [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
            disp('NRRD file loaded!')
            Geometry.ROIS{1}.ind = find(TRGT_img>0);
            Geometry.ROIS{1}.name= 'Target';
        case 2
            TRGT_in=am2mat([TRGT_path, TRGT_in]);
            disp('AM file loaded!')
            TRGT_img = permute(TRGT_in.data, [2,1,3]);
            Geometry.ROIS{1}.ind = find(TRGT_img>0);
            Geometry.ROIS{1}.name= 'Target';
        case 3
            error('please do not select *all* in file selection')
            % something selected
        otherwise
%             error('You should never see this.')
    end
    
    
        %% -----===== Save geometry files =====-----
    patient_dir = uifile('getdir', 'Save the patient data to directory');
    Geometry.patient_dir = patient_dir;
    
%     patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
    
    % make directories
    mkdir(Geometry.patient_dir);
    mkdir(fullfile(Geometry.patient_dir, 'beamlet_batch_files'));
    mkdir(fullfile(Geometry.patient_dir, 'geometry_files'));
    mkdir(fullfile(Geometry.patient_dir, 'matlab_files'));
    mkdir(fullfile(Geometry.patient_dir, 'opt_input'));
    mkdir(fullfile(Geometry.patient_dir, 'opt_output'));

    %% -----===== BTV and Ring creation =====-----
    waitbar(1/total_num_steps, hWaitbar, 'Step 2: Assign target and BTV margin');
    ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);

    [target_idx okay] = listdlg('ListString', ROI_names, ...
        'SelectionMode', 'single', 'Name', 'Target Selection', ...
        'PromptString', 'Please select the target ROI. ');
    if okay ~= 1
        msgbox('Plan creation aborted');
        delete(hWaitbar);
        return;
    end

    [BTV_margin_answer] = inputdlg({sprintf('Please enter the BTV margin (cm):\n(default 0.6 cm or 1 sigma, enter 0 to skip)')}, ...
        'BTV margin specification', 1, {'0.6'});
    if isempty(BTV_margin_answer)
        BTV_margin = 0.6;
    else
        BTV_margin = str2double(BTV_margin_answer{1});
    end

    % target_idx and BTV_margin are set. Expand PTV to BTV
    PTVmask = false(size(Geometry.rhomw));
    % for target_idx = target_indices
            PTVmask(Geometry.ROIS{target_idx}.ind) = 1;
    % end
    if BTV_margin > 0
        if exist('BTV_margin', 'var') && BTV_margin >= min(Geometry.voxel_size)
            bwD = bwdistsc(PTVmask, Geometry.voxel_size);
            Geometry.BTV = bwD <= BTV_margin;
        elseif BTV_margin < min(Geometry.voxel_size)
            warning('BTV margin too small!')
            Geometry.BTV = zeros(size(Geometry.data));
        end
    end

    % Create btv
    Geometry.ROIS{end+1} = Geometry.ROIS{end};
    Geometry.ROIS{end}.name = 'BTV';
    Geometry.ROIS{end}.num_curves = 0;
    Geometry.ROIS{end}.curves = {};
    Geometry.ROIS{end}.ind = find(Geometry.BTV);
    Geometry.ROIS{end}.visible = false;

    % Create ring
    [ring_margin_answer] = inputdlg({sprintf('Please enter the ring margin (cm):')}, ...
        'Ring margin specification', 1, {'1'});
    if isempty(ring_margin_answer)
        ring_margin = 1;
    else
        ring_margin = str2double(ring_margin_answer{1});
    end
    bwD = bwdistsc(PTVmask, Geometry.voxel_size);
    Geometry.Ring = bwD <= ring_margin; % default ring radius 3 cm
    Geometry.Ring = xor(Geometry.Ring, PTVmask);
    Geometry.ROIS{end+1} = Geometry.ROIS{end};
    Geometry.ROIS{end}.name = 'Ring';
    Geometry.ROIS{end}.num_curves = 0;
    Geometry.ROIS{end}.curves = {};
    Geometry.ROIS{end}.ind = find(Geometry.Ring);
    Geometry.ROIS{end}.visible = false;

    
    %% -----===== Save the matlab files =====-----
    waitbar(2.33/total_num_steps, hWaitbar, 'Step 3: Save matlab geometry files');
    % Save matlab geometry file
    save(fullfile(Geometry.patient_dir, 'matlab_files', 'Geometry.mat'), 'Geometry');

    waitbar(2.66/total_num_steps, hWaitbar, 'Step 3: Save raw geometry files');
    % Write binary geometry files
    fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'rhomw.bin'), 'w');
    fwrite(fid, Geometry.rhomw, 'single');
    fclose(fid);
    
    fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'Smw.bin'), 'w');
    fwrite(fid, Geometry.Smw, 'single');
    fclose(fid);

    fid = fopen(fullfile(Geometry.patient_dir, 'geometry_files', 'Fmw2.bin'), 'w');
    fwrite(fid, Geometry.Fmw2, 'single');
    fclose(fid);

    delete(hWaitbar);
    waitfor(msgbox(['Plan geometry created successfully in ' '"' Geometry.patient_dir '"']));
    

end