WiscPlan_v2/nrrd2geometry.m

function [Geometry patient_dir] = nrrd2geometry
    
    %% -----===== Get CT files =====-----
%     [IMG_in, IMG_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\PhantomData\', '*.nrrd']);
    
    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
    
    if IMG_in==0 % no file selected, no change to CT data
        warning('No file selected!')
    else 
        [Geometry.data, meta]=nrrdread([IMG_path, IMG_in]);
    end
    
    % Geometry.rhomw
    Geometry.rhomw = CT2dens(Geometry.data, 'pinn');
    Geometry.rhomw(Geometry.rhomw < 0.0012) = 0.0012;
    [Geometry.Smw Geometry.Fmw2] = dens2mstp(Geometry.rhomw);
    
    % Geometry.voxel_size
    VoxSizeString = meta.spacedirections;
    C=textscan(VoxSizeString(2:end-1), '%f', 9,'Delimiter',{',',') ('});
    C= C{1};
    Geometry.voxel_size = [C(1), C(5), C(9)];
    
    % Geometry.start
    StartString = meta.spaceorigin;
    C=textscan(StartString(2:end-1), '%f', 3,'Delimiter',',');
    Geometry.start=C{1};
    Geometry.start=Geometry.start';
    
    % Geometry.ROIS
    
    %% -----===== Get mask/contour files =====-----
%     [TRGT_in, TRGT_path, ~] = uigetfile(['C:\010-work\003_localGit\WiscPlan_v2\PhantomData\','*.nrrd']);
    TRGT_in = 'test_seg.nrrd';
    TRGT_path = 'C:\010-work\003_localGit\WiscPlan_v2\PhantomData\';
    
    if TRGT_in==0 % no file selected, no change to CT data
        warning('No file selected!')
    else 
        [TRGT_img, meta]=nrrdread([TRGT_path, TRGT_in]);
    end
    Geometry.ROIS{1}.ind = find(TRGT_img>0);
    Geometry.ROIS{1}.name= 'Target';

    
        %% -----===== Save geometry files =====-----
%     patient_dir = uifile('getdir', 'Save the patient data to directory');
    patient_dir = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\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;

    
    
    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