function RtomoAll = tomoBeamOnOffSinogram(GeomROI,xpmin,xpmax,Mxp,Nphi,Nrot)

% Figures out which beams should be on or off for tomotherapy treatment.


ptv = [];
ptv.x_count = GeomROI.dim(1);
ptv.y_count = GeomROI.dim(2);
ptv.z_count = GeomROI.dim(3);
ptv.non_zero_indices = GeomROI.ind;
ptv.non_zero_values = single(ones(length(GeomROI.ind),1));

ptvFull = full3D(ptv);

tumorMask = [];
tumorMask.start = GeomROI.start;
tumorMask.voxel_size = GeomROI.pixdim;
tumorMask.data = ptvFull;

[M,N,Q] = size(tumorMask.data);

% create axes
x = [0:M-1]*tumorMask.voxel_size(1) + tumorMask.start(1);
y = [0:N-1]*tumorMask.voxel_size(2) + tumorMask.start(2);
z = [0:Q-1]*tumorMask.voxel_size(3) + tumorMask.start(3);

% Calculate the maximum tumor cross section
bigSlice = double(~~sum(tumorMask.data,3));

% grow the tumor by 5 voxels
BW = bwdist(bigSlice);
bigSlice(BW < 5) = 1;

bigSlice = rot90(bigSlice);

xpTomo = [-(Mxp-1)/2:(Mxp-1)/2]*(xpmax - xpmin)/Mxp;

NphiTomo = 51;
phiTomo = [0:NphiTomo-1]*360/NphiTomo;

phi = [0:359];
% calculate the radon transform of the mask
[R,xp] = radon(bigSlice,phi);
xp = xp*tumorMask.voxel_size(1);

% create downsampling matrices
[PHI,XP] = meshgrid(phi,xp);
[PHITOMO,XPTOMO] = meshgrid(phiTomo,xpTomo);

% resample the sinogram
Rinterp = interp2(PHI,XP,R,PHITOMO,XPTOMO);
Rinterp(isnan(Rinterp) == 1) = 0;
Rtomo = ~~Rinterp;

% make copies of sinogram for 3D dose calculation

RtomoDouble = double(Rtomo);

% convert 2D sinograms to 3D for dose calc

RtomoDouble = reshape(RtomoDouble,Mxp*NphiTomo,1);
rowVector = ones(1,Nrot);
RtomoAll = RtomoDouble*rowVector;
RtomoAll = reshape(RtomoAll,Mxp,NphiTomo*Nrot);

fprintf('Unused beams turned off.\n')