WiscPlan_v2/WiscPlanPhotonkV125/matlab_frontend/resample_geometry.m

function resample_geometry(Geometry, Nx, Ny, Nz)
    Geometry_old = Geometry;
    str = inputdlg({'Enter X voxel dimension (cm):', 'Enter Y voxel dimension (cm):', ...
        'Enter Z voxel dimension (cm):'}, 'input', [1,8], {'0.2', '0.2', '0.2'});
    % dimensions in mm
    X_dim = str2double(str{1});
    Y_dim = str2double(str{2});
    Z_dim = str2double(str{3});
    disp(['Desired voxel dimensions: ' num2str(X_dim) ' ' num2str(Y_dim) ' ' num2str(Z_dim) ' mm'])
    
    % round to appropriate number of voxels
    Nx = round(size(Geometry_old.data,1)*Geometry_old.voxel_size(1)/X_dim);
    Ny = round(size(Geometry_old.data,2)*Geometry_old.voxel_size(2)/Y_dim);
    Nz = round(size(Geometry_old.data,3)*Geometry_old.voxel_size(3)/Z_dim);
    
    % find the correct scale factor
    fx = Nx/size(Geometry_old.data,1);
    fy = Ny/size(Geometry_old.data,2);
    fz = Nz/size(Geometry_old.data,3);
    
    %% start converting the Geometry file
    Geometry.voxel_size = Geometry_old.voxel_size./[fx, fy, fz];
    disp(['New voxel sizes: ' num2str(Geometry.voxel_size(1)) ' ' num2str(Geometry.voxel_size(2)) ...
        ' ' num2str(Geometry.voxel_size(3)) ' mm'])
    
    load('WiscPlan_preferences.mat')
    if ~isfield(WiscPlan_preferences,'inDataPath')
        Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
    elseif WiscPlan_preferences.inDataPath ~= 0
        Geometry.patient_dir = uigetdir(WiscPlan_preferences.inDataPath, 'Select the new patient directory');
    else
        Geometry.patient_dir = uigetdir('C:', 'Select the new patient directory');
    end
    
    Geometry.data = imresize3(Geometry_old.data, [Nx, Ny, Nz]);
    Geometry.rhomw = imresize3(Geometry_old.rhomw, [Nx, Ny, Nz]);
    Geometry.Smw = imresize3(Geometry_old.Smw, [Nx, Ny, Nz]);
    Geometry.Fmw2 = imresize3(Geometry_old.Fmw2, [Nx, Ny, Nz]);
    Geometry.BTV = logical(round(imresize3(double(Geometry_old.BTV), [Nx, Ny, Nz])));
    Geometry.Ring = logical(round(imresize3(double(Geometry_old.Ring), [Nx, Ny, Nz])));

    Geometry.x = Geometry.start(1) +linspace(0, Nx-1, Nx)*Geometry.voxel_size(1);
    Geometry.y = Geometry.start(2) +linspace(0, Ny-1, Ny)*Geometry.voxel_size(2);
    Geometry.z = Geometry.start(3) +linspace(0, Nz-1, Nz)*Geometry.voxel_size(3);
    
    Geometry.ROIS = {};
    for ROI_i = 1:numel(Geometry_old.ROIS)
       Geometry.ROIS{ROI_i}.name = Geometry_old.ROIS{ROI_i}.name;
       tabula = zeros(size(Geometry_old.data));
       tabula(Geometry_old.ROIS{ROI_i}.ind) = 1;
       tabula = imresize3(tabula, [Nx, Ny, Nz]);
       tabula = imgaussfilt3(tabula, 0.5);
       Geometry.ROIS{ROI_i}.ind = find(tabula>0.3);
    end
    
    %% create export folder and start populating with files
    mkdir([Geometry.patient_dir '\matlab_files']);
    save([Geometry.patient_dir '\matlab_files\Geometry.mat'], 'Geometry');
    
    % resample beamlets
    if false
%     if isfile([Geometry_old.patient_dir '\batch_dose_backup.bin'])
        disp('calculated beamlet file found - resampling')
        
        % resample the beamlets file
        f = waitbar(0, 'Resampling beamlets');
        
        beamlets_in = read_ryan_beamlets([Geometry_old.patient_dir '\batch_dose_backup.bin'], 'ryan');
        for beam_i = 1:numel(beamlets_in)
            
            waitbar(beam_i/numel(beamlets_in),f, ['Correcting beamlet shift: ' num2str(beam_i)])
            tabula = zeros(beamlets_in{beam_i}.x_count, beamlets_in{beam_i}.y_count, beamlets_in{beam_i}.z_count);
            
            beamlets_out{beam_i}.num = beamlets_in{beam_i}.num-1; %reading beamlets adds +1 to number, this fixes it
            beamlets_out{beam_i}.x_count = Nx;
            beamlets_out{beam_i}.y_count = Ny;
            beamlets_out{beam_i}.z_count = Nz;
            
            tabula(beamlets_in{beam_i}.non_zero_indices) = beamlets_in{beam_i}.non_zero_values;
            tabula = imresize3(tabula, [Nx, Ny, Nz]);
            
            beamlets_out{beam_i}.non_zero_indices = find(tabula>0.005*max(tabula(:)));
            beamlets_out{beam_i}.non_zero_values = tabula( beamlets_out{beam_i}.non_zero_indices);

        end
        
        write_ryan_beamlets([Geometry.patient_dir '\batch_dose.bin'],beamlets_out);
        write_ryan_beamlets([Geometry.patient_dir '\batch_dose_backup.bin'],beamlets_out);
        
        close(f)
        
    else
        disp('no beamlet file found')
    end
    disp('done!')

end