WiscPlan_v2/WiscPlanPhotonkV125/matlab_frontend/nrrd2geometry.m

function [Geometry patient_dir] = nrrd2geometry
    
    %% -----===== Get CT files =====-----
    [IMG_in, IMG_path, filterIdx] = uigetfile([{'*.nrrd'; '*.am'}], 'Select CT image', 'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\');
    
%     IMG_in = 'test_ct.nrrd';
%     IMG_path = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\';
    
    total_num_steps = 4;
    hWaitbar = waitbar(0);
    
    % 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 - self explanatory
    % Geometry.start - self explanatory
    
    % 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
    
    switch filterIdx
        case 0
            warning('No file selected!')
        case 1
            [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
%             warning('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 = meta.spacedirections;
            Geometry.start=meta.spaceorigin;
        case 2
            imgIn=am2mat([IMG_path, IMG_in]);
            Geometry.data = permute(imgIn.data, [2,1,3]);
%             warning('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
        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', 'C:\Box Sync\Optimal Stopping\Data\CDP\CDP001\');
    
    switch filterIdx2
        case 0
            warning('No file selected!')
        case 1
            [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
            warning('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]);
            warning('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
            % something selected
        otherwise
%             error('You should never see this.')
    end
    
    
        %% -----===== Save geometry files =====-----
    patient_dir = uifile('getdir', 'Save the patient data to directory');
%     patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\data\PatientData\';
    
    % make directories
    mkdir(patient_dir);
    mkdir(fullfile(patient_dir, 'beamlet_batch_files'));
    mkdir(fullfile(patient_dir, 'geometry_files'));
    mkdir(fullfile(patient_dir, 'matlab_files'));
    mkdir(fullfile(patient_dir, 'opt_input'));
    mkdir(fullfile(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;
%         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(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(patient_dir, 'geometry_files', 'rhomw.bin'), 'w');
    fwrite(fid, Geometry.rhomw, 'single');
    fclose(fid);
    
    fid = fopen(fullfile(patient_dir, 'geometry_files', 'Smw.bin'), 'w');
    fwrite(fid, Geometry.Smw, 'single');
    fclose(fid);

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

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

    % hWaitbar = waitbar(1/total_num_steps, 'Step 4, Create optimization geometry');

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

end